MODELING SOIL WATER REDISTRIBUTION UNDER SURFACE DRIP IRRIGATION

ABSTRACT

The main objective of this study was to develop a numerical model to simulate water distribution, and the shape of the wetted soil volume, resulting from a point source irrigation at the ground surface (dripper), given soil hydraulic properties and irrigation system parameters. The mathematical model was developed in Fortran 90 at the Department of Biosystems Engineering, at the Luiz de Queiroz College of Agriculture, ESALQ/USP. The computer program was structured to allow the user to input information such as: a) flow and transport related to soil properties, b) information on irrigation type, c) boundary conditions, d) simulation time, and e) water application time through irrigation. Data from another experiment carried out at the Department of Biosystems Engineering were used to validate the model. Model performance was evaluated based on Willmott concordance index, coefficient of efficiency, root mean square error, mean error, and maximum absolute error. Based on comparisons of the model deviations from the measured data and with other results reported in the literature, it was clear that the predictions of the proposed numerical model were very satisfactory.

KEYWORDS
Wet Bulb; Finite Differences; Trickle Irrigation

INTRODUCTION

The scarcity of water in many regions of the world, combined with the large amount of water used in agriculture, revels the need to use more efficient irrigation practices. Micro irrigation allows the application of small water level in a localized way implying a greater application efficiency. Drip irrigation is one of the most efficient systems of transporting water and nutrients to the root zone of the plant, with the objective of providing water to a limited volume of soil in the region where the largest water extraction by plants takes place (Naglič et al., 2014; Kandelous & Šimůnek et al., 2010). Once the form of distribution of moisture within a volume of wet soil is known, also called of wet bulb, the emitter or emitters can be arranged in the camp design so that the plant can consume water and nutrients more efficiently.

Many researches were carried out in order to determine the distribution of water in the soil by drip irrigation using mathematical methods, from the physical properties of the soil and the irrigation system with satisfactory results (Siyal & Skaggs, 2009; Kandelous et al., 2011; Samadianfard et al., 2012; Arbat et al., 2013; Subbauah & Mashru, 2013).

The physical approximation that describes the movement of water in unsaturated soils can be represented by Richards equation (van Dam & Feddes, 2000; Li et al., 2015). There are different analytical solutions for this equation, which can be applied to drip irrigation, such as, for example, the one developed by Abid et al. (2012). However, because Richards is a nonlinear second order partial differential equation, the analytical solutions are obtained mainly from simplifications of the hydraulic properties of the soil. For more complex problems, and with irregular geometry, it is necessary to use the numerical techniques (van Dam & Feddes, 2000; Yao et al., 2011; Abid, 2014; Li et al. 2015; Šimůnek et al., 2016). Therefore, the objective of the present study was to develop a numerical model capable of simulating the water distribution and the shape of the wetted soil from the irrigation by a point source in the soil surface (dripper), using the hydraulic properties of the soil and the soil irrigation system as input variables.

MATERIAL AND METHODS

The mathematical model was developed in the Department of Biosystems Engineering at the Luiz de Queiroz College of Agriculture, ESALQ/USP, using the Fortran 90 programming language. The computational program was structured to allow the user to enter information such as: a) soil profile data, with regard to their physical-hydric properties, b) information about the irrigation system, c) boundary conditions, d) simulation time and e) water application time through irrigation. The model presents some simplifications in the solution of Richards' equation, such as, not considering the environment as isotropic and isothermal, the flow in macropores and the flow of steam in the soil. These simplifications were performed with the intention of leaving the model with a smaller number of input parameters.

The Richards' equation which describes the movement of water in an isothermal porous medium, two-dimensional, with the positive vertical coordinate downwards, and under unsaturated conditions, can be described by [eq. (1)], expressing the hydraulic conductivity as a function of the matric potential:

(1) $$C\left(h\right)\frac{\partial h}{\partial t}=\frac{\partial }{\partial x}\left[K_x\left(h\right)\frac{\partial h}{\partial x}\right]+\frac{\partial }{\partial z}\left[K_z\left(h\right)\frac{\partial \left(h\right)}{\partial z}\right]+\frac{\partial K_z\left(h\right)}{\partial z}$$

Where,

h is the matric potential [L];

K (h) is the unsaturated hydraulic conductivity [L T⁻¹], due to the matric potential,

C (h) is the specific water capacity function [L⁻¹], which is equal to the slope of the soil water retention curve.

In this expression, partial derivatives of h appear with respect to space and time, which can be replaced by finite differences. The coefficients are function of the dependent variable h, having their values estimated for the different situations of time and space. The finite differences approximation implies that the calculation domain and the time are discretized. Thus, the calculation domain (Ω) is represented by a set of points that occupy the “knots” of a rectangular mesh. The mesh divides Ω on a regular basis (Figure 1).

FIGURE 1
Discretization scheme for domain calculation (space).

Among several possible forms of resolution by finite differences, an implicit scheme was chosen. This means that, between instants t and t + 1, the partial derivatives in order of space were evaluated based on the values of h at time t + 1 (implicit scheme).

Thus, at the time step Δt, between the instants t and t + 1, at point (i, j), the [eq. (1)] was discretized as follows (equation 2):

(2) $$\begin{array}{l}{C}_{i,j}^t\frac{\left(h_{i,j}^{t+1}-h_{i,j}^t\right)}{\Delta t^t}=\frac{K_{i+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,j}^t\left(\frac{h_{i+1,j}^{t+1}-h_{i,j}^{t+1}}{\Delta x}\right)-K_{i-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,j}^t\left(\frac{h_{i,j}^{t+1}-h_{i-1,j}^{t+1}}{\Delta x}\right)}{\Delta x}\hfill \\ +\frac{K_{i,j+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}^t\left(\frac{h_{i,j+1}^{t+1}-h_{i,j}^{t+1}}{\Delta z_l}+1\right)-K_{i,j-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}^t\left(\frac{h_{i,j}^{t+1}-h_{i,j-1}^{t+1}}{\Delta z_u}+1\right)}{\Delta z_{i,j}}\hfill \end{array}$$

In which the indices indicate i the column of the mesh, j the line, and t the time, $\Delta Z_l=Z_{i,j+1}-Z_{i,j}$ and $\Delta Z_u=Z_{i,j}-Z_{i,j-1}$

Due to the high nonlinearity of the specific water capacity C (h), we have for each step of time, mass balance errors, when highly transient conditions are simulated (van Dam & Feddes, 2000). In the present model, we used the modification in the solution of [eq. (2)], proposed by van Dam & Feddes (2000), generating [eq. (3)].

(3) $${\theta }_{i,j}^{t+1}-{\theta }_{i,j}^t=C_{i,j}^{t+1,p-1}\left(h_{i,j}^{t+1,p}-h_{i,j}^{t+1,p-1}\right)+{\theta }_{i,j}^{t+1,p-1}-{\theta }_{i,j}^t$$

Where,

p is the number of iterations.

Substituting [eq. (3)] into [eq. (2)], we have [eq. (4)]:

(4) $$\begin{array}{l}{C}_{i,j}^{t+1,p-1}\left(h_{i,j}^{t+1,p}-h_{i,j}^{t+1,p-1}\right)+{\theta }_{i,j}^{t+1,p-1}-{\theta }_{i,j}^t=\\ \Delta t^t\left[\frac{K_{i+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,j}^t\left(\frac{h_{i+1,j}^{t+1,p}-h_{i,j}^{t+1,p}}{\Delta x}\right)-K_{i-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,j}^t\left(\frac{h_{i,j}^{t+1,p}-h_{i-1,j}^{t+1,p}}{\Delta x}\right)}{\Delta x}\right]\\ +\Delta t^t\left[\frac{K_{i,j+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}^t\left(\frac{h_{i,j+1}^{t+1,p}-h_{i,j}^{t+1,p}}{\Delta z_l}+1\right)-K_{i,j-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}^t\left(\frac{h_{i,j}^{t+1,p}-h_{i,j-1}^{t+1,p}}{\Delta z_u}+1\right)}{\Delta z_{i,j}}\right]\end{array}$$

Developing the equation (4) and grouping the terms, we obtain:

(5) $$A_{i,j}{h}_{i,j-1}^{t+1,p}+B_{i,j}{h}_{i,j+1}^{t+1,p}+D_{i,j}{h}_{i,j}^{t+1,p}+E_{i,j}{h}_{i-1,j}^{t+1,p}+F_{i,j}{h}_{i+1,j}^{t+1,p}=G_{i,j}$$

Where,

(6) $$A_{i,j}=-\frac{\Delta t^t}{\Delta z_{i,j}\Delta z_u}{K}_{i,j-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}^t$$

(7) $$B_{i,j}=-\frac{\Delta t^t}{\Delta z_{i,j}\Delta z_l}{K}_{i,j+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}^t$$

(8) $$D_{i,j}=\left[\frac{\Delta t^t}{\Delta x^2}\left(K_{i+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,j}^t+K_{i-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,j}^t\right)+\frac{\Delta t^t}{\Delta z_{i,j}\Delta z_l}{K}_{i,j+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}^t+\frac{\Delta t^t}{\Delta z_{i,j}\Delta z_u}{K}_{i,j-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}^t+C_{i,j}^{t+1,p-1}\right]$$

(9) $$E_{i,j}=-\frac{\Delta t}{\Delta x^2}{K}_{i-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,j}^t$$

(10) $$F_{i,j}=-\frac{\Delta t}{\Delta x^2}{K}_{i+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,j}^t$$

(11) $$G_{i,j}=C_{i,j}^{t+1,p-1}{h}_{i,j}^{t+1,p-1}-{\theta }_{i,j}^{t+1,p-1}+{\theta }_{i,j}^t-\frac{\Delta t^t}{\Delta z_{i,j}}\left(K_{{}^{i,j+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}^t-K_{i,j+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}^t\right)$$

For the convergence analysis of the iterative solution of the [eq. (4)] the recommendation of van Dam & Feddes (2000) was used, in which the criterion based on θ was implemented as convergence analysis. van Dam & Feddes (2000) report that using the criterion based on θ, the simulations are performed in less time, without sacrificing precision in the mass balance. The value of the function K (h) is explicitly linearized in [eq. (4)]. In this model the arithmetic mean was used. The model used to describe the water retention curve in the soil was the van Genuchten model. For the calculation of the hydraulic conductivity, the model developed by Mualem was used (further details of the van Genuchten and Mualem equation can be checked in Li et al., 2015).

Contrary to the great majority of the models of water movement simulation in two dimensions, we opted for a more realistic way, by the hypothesis of a non-uniform initial water profile. For the first irrigation, it was assumed that the initial matric potential h_o depended only on the depth, while at the beginning of the subsequent irrigations, a variation along the horizontal axis was also considered.

Therefore:

For the first irrigation:

(12) $$h(x,z,0)=h_0(z)0\le x\le X_{\max },0\le z\le Z_{\max }$$

For the following irrigations we have:

(13) $$h(x,z,0)=h_0(x,z)0\le x\le X_{\max },0\le z\le Z_{\max }$$

The computation domain Ω is a rectangle (Figure 1) whose four sides constitute its contours (boundaries). In this way, a Cartesian coordinate system was considered in which flow directions X and Z were established. The following boundary conditions were adopted:

a) AC Boundary (Figure 1)

Because it is a Boundary of the domain where its neighboring cells in the negative direction of the X axis belong to one of the quadrants of the total soil volume, given the symmetry, we are faced with a situation equivalent to a null flow (or Neumann condition), which is:

(14) $$K(h)\frac{\partial h}{\partial x}=0t\ge 0,0\le z\le Z_{\max },\text{x}=0$$

b) BD Boundary (Figure 1)

Once the hypothesis of insulation of the bulbs or wetting front has been established, and as this boundary is defined so that the wetting front does not reach the boundary, there is also a null flow condition (Neumann condition).

(15) $$K\left(h\right)\frac{\partial h}{\partial x}=0t\ge 0,0\le z\le Z_{\max },x=X_{\max }$$

c) CD Boundary (Figure 1)

The lower boundary was displaced in such a way that the influence of the irrigation water in this zone was null (Dirichlet condition).

(16) $$K\left(h\right)\frac{\partial h}{\partial x}=0t\ge 0,0\le x\le X_{\max },z=Z_{\max }$$

d) AB Boundary (Figure 1)

As in the surface drip irrigation only one “knot” receives all the flow of the dripper, that is, source point, and considering from it the dampness front spreads through the domain Ω, there are two distinct zones: a region of the domain that receives the flow of the dripper (equation 17) and another region that undergoes evaporation (if the user wishes to simulate) or the absence of flow, eqs (18) and (19), respectively.

(17) $$-K\left(h\right)\frac{\partial h}{\partial z}+K\left(h\right)=q_{i}t\le t_{p}x=0,\text{z}=0$$

(18) $$-K\left(h\right)\frac{\partial h}{\partial z}+K\left(h\right)=Evt>t_{p}0<x\le X_{\max },\text{z}=0$$

(19) $$-K\left(h\right)\frac{\partial h}{\partial z}+K\left(h\right)=0t>t_{p}0<x\le X_{\max },\text{z}=0$$

Where,

q_i is the water flow.

For validation of the model, the data obtained by Rivera (2004) were used in an experiment conducted in the Department of Biosystems Engineering of the Luiz de Queiroz College of Agriculture. The soil material used originated from a profile classified as Red Latosol, sandy phase, collected inside the ESALQ/USP. The collection was done from a uniform layer that “extended” from the surface to a depth of 30 centimeters. Table 1 shows the physical-water characteristics of the soil and Table 2 shows the parameters of the retention curve.

TABLE 1
Soil physical and hydraulic properties of the experimental soil.

Texture			Density		Porosity	Ko
Sand (%)	Silt (%)	Clay (%)	Soil (g cm−3)	Particles (g cm−3)	(%)	(cm h−1)
67	6	27	1.369	2.575	46.835	10.221

Rivera (2004)

TABLE 2
Retention curve parameters, according to van Genuchten model.

θr (cm3 cm−3)	θs (cm3 cm−3)	α (cm−1)	N	M
0.113	0.482	0.029428	1.828069	0.452975

Rivera (2004)

In order to simulate a dripper, Rivera (2004) used a serum calibrator previously calibrated for a flow rate of 3 L h⁻¹. This dispenser was located in the center of the polyethylene carton containing the soil and was coupled to a 10-liter capacity Mariotte flask by means of capillary tube, keeping the hydraulic charge constant in the flask.

The application time was two hours, and a solution volume of 6 liters was therefore applied. Soil moisture, after the test, was determined using the gravimetric method. The sampling points were located along a mesh, taking as the central axis the point where the emitter was located; from that point it was shown every 10 cm in the horizontal direction and 10 cm in the vertical along two rays, so that every schematized ring was sampled twice. The total sampled rays were six (two replicates for each time), arranged to form on the surface of the soil angles of 60 degrees, that is, the bulb was divided into six slices of equal sizes. In both radial and vertical directions, 5 samples were taken, totaling 25 samples per radius. The sampling times were: before irrigation; 24; 48; and 72 hours after the end of irrigation.

Based on the temporal and spatial distribution of volumetric moisture (θ), the model was evaluated by comparing the values of θ obtained experimentally by Rivera (2004) with those obtained by the proposed model. For the comparisons between observed and simulated data the following statistical indices were used, as suggested by Legates & Mccabe (1999): Willmott concordance index (id); efficiency coefficient (E), square root of the mean error (RMSE), mean error (ME) and absolute maximum error (AMAXE). These indices are defined by eqs. (20-24):

(20) $$Id=1-\frac{{\displaystyle \sum _{i=1}^n{\left({\text{P}}_{\text{i}}-{\text{O}}_{\text{i}}\right)}^2}}{{\displaystyle \sum _{i=1}^{\text{n}}{\left(\left|{\text{P}}_{\text{i}}-\text{O}\right|+\left|{\text{O}}_{\text{i}}-\text{O}\right|\right)}^2}}$$

(21) $$E=1-\frac{{\displaystyle \sum _{i=1}^n{\left({\text{P}}_{\text{i}}-{\text{O}}_{\text{i}}\right)}^2}}{{\displaystyle \sum _{i=1}^{\text{n}}{\left({\text{O}}_{\text{i}}-\text{O}\right)}^2}}$$

(22) $$SRME=\sqrt{{\displaystyle \sum _{i=1}^N\frac{{\left({\text{P}}_{\text{i}}-{\text{O}}_{\text{i}}\right)}^2}{N}}}$$

(23) $$EM=\frac{{\displaystyle \sum _{i=1}^n\left({\text{P}}_{\text{i}}-{\text{O}}_{\text{i}}\right)}}{N}$$

(24) $$AMAXE=\max {\left\{\left|\text{Pi}-\text{Oi}\right|\right\}}_{i=1}^N$$

Where,

O_i are the observed data, obtained in the experimental tests;

P_i are the data simulated by the model;

N is the number of observations, and

O is the mean of the observed values.

RESULTS AND DISCUSSION

Figure 2 shows a comparison between the measured and simulated results by the finite difference method for the distribution of soil moisture (θ) after 24 hours of redistribution at different horizontal and vertical distances. Simulated results showed agreement with the observed data, the largest difference being observed for the distance of 5 cm and depth of 35 cm (Figure 2A), in which the model presented an underestimation of the value of θ.

FIGURE 2
Observed and simulated soil water content (θ) at different depths and distances, 24 hours after the end of irrigation: (A) 5 cm; (B) 15 cm; (C) 35 cm; (D) 55 cm.

For the time of 24 hours after the simulation beginning, the model estimated the θ values more precisely for the farthest points of the emitter (Figure 2C and Figure 2D); this result can be attributed to the fact that at the end of 24 hours the water had not yet been distributed throughout the domain. Therefore, for the points located farthest from the emitter, the soil moisture values would still be close to the initial simulation condition.

Figure 3 shows the distribution profiles shown by means of θ isolines for 24 hours after irrigation, generated by the implicit finite difference methods and observed values. It can be observed that soil moisture, both observed as simulated, ranged about 0.12 to 0.20 cm³ cm⁻³. According to Figure 3, it can be seen that the soil moisture values obtained by Rivera (2004) varied in the range of 0.17 to 0.20 cm³ cm⁻³, between the distance of 33 cm and depth of 35 cm compared to the emitter, decreasing both horizontally and vertically, as it moves away from the emitter, the bulb acquires a hemispherical shape. This same behavior was simulated with great precision by the proposed model.

FIGURE 3
Observed and simulated soil water content 24 hours after the end of irrigation.

From the observation of the isolines shown in Figure 3, in general, the water distribution in the simulated soil showed good agreement when compared to that observed in both depth and width of the wet zone.

To better quantitatively evaluate the accuracy of the model predictions, the following performance indices were calculated: the Willmott concordance index (id); the coefficient of efficiency (E), the root mean square error (RMSE), the mean error (ME) and the absolute maximum error (AMAXE), involving the values of θ observed and those estimated by the model. These indices can be checked by Table 3; the value of the id was 0.9869, while the coefficient of efficiency (E), also known as the Nash and Sutcliffe coefficient, was 0.9522. According to Santos (2011), the coefficient of efficiency (E) represents the ratio between the mean error square of the estimate and the observed data variance, subtracted from the unit; this coefficient varies from -∞ to 1. Values of E equal to zero indicate that the mean of the observed data is an estimate as good as the values predicted by the model; when E ≤ 0, the mean of the observed values is a better estimate than the values estimated by the model. This index, therefore, presents a superiority in relation to the Willmott index in terms of interpretation. Therefore, the values simulated by the model according to the id and E indices indicated that the values of θ were accurately estimated.

TABLE 3
Statistical comparison of simulated and observed soil moisture content at different redistribution times.

Time (hours)	Id	E	RMSE	ME	AMAXE
24	0.9869	0.9522	0.0044	0.0002	0.0173
48	0.9384	0.7931	0.0072	0.0018	0.0236
72	0.9851	0.9444	0.0033	0.0005	0.0151

The simulation of water redistribution in the soil after 24 hours from the beginning by the finite difference method showed RMSE values equal to 0.0044, the MS equal to 0.0002 and the AMAXE equal to 0.0173, thus showing the satisfactory results of the model proposed that uses the technique of the finite differences when compared to the sampling (Table 3). There is a range of studies reporting good model efficiency in the wet bulb forecast, comparing simulation results with experimental data (Siyal & Skaggs, 2009; Kandelous et al., 2011; Samadianfard et al., 2012; Arbat et al., 2013; Subbauah & Mashru, 2013; Abid, 2014).

From Figure 4, it is possible to observe a comparison between the soil water content vertically, at different distance positions from the emitter, for the redistribution time of 48 hours after irrigation. A significant agreement can be noted between the values simulated by the proposed model and the observed values. The largest differences between the observed and simulated values were for the distance of 35 cm (Figure 4C) and for the distance of 55 cm (Figure 4D). It can be observed from Figure 4C that the model underestimated the values at the lowest depths and overestimated them to greater depths. In other words, the model was not able to predict lower percolation at the most superficial depths and most percolation in the deeper layers of the soil. At the distance of 55 cm from the emitter (Figure 4D), the model overestimated the values of θ at almost all depths. For the distance of 55 cm the model estimated a greater percolation than it was observed by the data.

FIGURE 4
Observed and simulated soil water content (θ) at different depths and distances, 48 hours after the end of irrigation: (A) 5 cm; (B) 15 cm; (C) 35 cm; (D) 55 cm.

The comparison of the water distribution in the soil obtained by the experimental (Rivera, 2004) and simulated values, evidenced a similar water distribution pattern (Figure 5). It was observed that the dimensions of the bulb remained statistically constant when compared with the time of 24 hours of redistribution. However, there was a decrease in soil bulk values, especially in the cells near the emitter point, for both the observed and simulated data. At the distances closer to the emitter there was a decrease in humidity when compared to the 24 hours time. However, the decrease in the number of upper cells resulted in an increase in humidity in the adjacent lateral and inferior cells, a product of the redistribution of water in the soil.

FIGURE 5
Observed and simulated soil water content 48 hours after the end of irrigation.

Comparing the performance of the model in the simulation for water redistribution after 48 hours, it was observed that the coefficient of efficiency (E) was 0.7931 (Table 3). As the value of E> 0 (Table 3), it was concluded that the values simulated by the methodology tested were better than the mean values observed.

A positive or negative value of ME indicates that the model overestimates or underestimates the experimental results, respectively. According to Table 3, the models overestimated the experimental results by 0.0018 cm³ cm⁻³. Siyal & Skaggs (2009) evaluated the efficacy of the traditional Hydrus 2/3D model compared to experimental data and found the mean value of MS equal to −0.006 cm³ cm^-3.

According to Figure 6, a good agreement can be observed between the values simulated by the proposed model using the finite difference technique and the observed values of the soil water content vertically at all different distance positions from the emitter for the redistribution time of 72 hours after irrigation.

FIGURE 6
Observed and simulated soil water content (θ) at different depths and distances, 72 hours after the end of irrigation: (A) 5 cm; (B) 15 cm; (C) 35 cm; (D) 55 cm.

For the redistribution time of 72 hours the model did not show any tendency to overestimate or underestimate the values of θ (Figure 6), showing that it can be used to simulate the distribution and redistribution of water in the soil for problems of drip irrigation.

From Figure 7 we can observe the distribution profiles shown by means of θ isolines for 72 hours after irrigation, generated by the finite difference and observed methods. In both conditions (simulated and observed data) the cells near the emitter continued to “lose” water to adjacent cells, but with low intensity; however, this did not contribute significantly to changes in bulb dimensions. Both moistures (simulated and observed) had an approximate behavior regarding the distribution of the water content in the bulb (Figure 7).

FIGURE 7
Observed and simulated soil water content 72 hours after the end of irrigation.

The model presented again a good estimate of the distribution of water in the soil, considering that the θ isolines obtained by the model coincide as the isolines obtained from the values observed in the field.

Rivera et al. (2006) commented that in this type of soil it can be said that, after 24 hours, the redistribution of the water inside the wet bulb practically ends, being the major alterations observed in the cells near the emitter, that is, the water flow that occurred after this time was minimal, so that the dimensions of the bulb began to remain practically unchanged. This same behavior was successfully simulated by the proposed model, thus showing the potentiality of the model technique to be used for studies of wet bulb estimation, which is an important parameter in the design and management of the microirrigation system (Siyal & Skaggs, 2009; Subbaiah, 2011; Li et al., 2015).

From the values observed in Table 3, it can be seen that the value of id was higher than 0.9, so it can be said that the model estimated with agreement the redistribution of water in the soil after 72 hours from the end of irrigation. Also, according to Table 3, we have the value of E> 0, meaning that the values simulated by the methodology tested were better than the mean values observed.

It is important to note that, for the “ME” and “AMAXE” indices, there is the difference of humidity as analyzed variable (cm³ cm⁻³). Thus, the value of AMAXE for the finite difference methods represented a difference of 0.0151 cm³ cm⁻³ (1.51 %) concerning to the observed value (Table 3). The positive value of “ME”, occurred for the simulation, indicates that the models overestimated the values obtained experimentally (Table 3).

Analyzing the values obtained for the RMSE index (Table 3), it was observed that for the methodology tested for the solution of water redistribution in the soil, the model presented low RMSE value, 0.0033. Kandelous & Šimůnek (2010) found RMSE values ranging up to a maximum of 0.045 cm³ cm⁻³, in simulations of a subsurface drip system as Hydrus 2/3D model.

Analyzing the performance of the proposed model used for the three times of water redistribution in the soil (Table 3), it can be said that the model achieved a good performance according to the data observed; such an assertion can be justified by the low RMSE values (close to zero), low ME values, and higher values of “E” indices as well as id values always close to unity.

The small differences observed between the θ estimated and observed values by the model are possibly due to the limitations of the model. However, comparing the deviations of the model with respect to the measured data, and comparing the deviations of this model with other results presented in the literature, it can be affirmed that the predictions of the proposed numerical model were very satisfactory.

CONCLUSIONS

The developed model allowed the distribution of water in the soil under superficial drip irrigation, considering the main parameters involved in the shape and dimensions of the volume of the wet soil. According to the statistical parameters of the model evaluation, it can be concluded that a very satisfactory agreement was obtained in a water distribution test with the respective parameters in the experiments.

ACKNOWLEDGMENTS

This work was supported by Coordination for the Improvement of Higher Education Personnel (CAPES). The authors would like to thank the team of the Laboratory of Soil Physics of the Department of Biosystems Engineering (Luiz de Queiroz College of Agriculture (ESALQ/USP)). We are also grateful for all important comments from the anonymous reviewers, which contributed to a significant upgrade of the manuscript.

REFERENCES

• Abid MB (2014) Numerical Simulation of Two-Dimensional Unsaturated Flow from a Trickle Irrigation Source Using the Finite-Volume Method. Journal of Irrigation and Drainage Engineering 141(4):06014005.
• Abid MB, Hamad SN, Hussain SM (2012) Analytical solution of unsaturated soil water flow from a point source. Journal of Engineering 18(1):95-106.
• Arbat G, Puig-Bargués J, Duran-Ros M, Barragán J, Ramírez de Cartagena F (2013) Drip-irriwater: computer software to simulate soil wetting patterns under surface drip irrigation. Computers and Electronics in Agriculture 98:183-192.
• Kandelous MM, Š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.
• Kandelous MM, Šimůnek J, Van Genuchten MT, Malek K (2011) Soil water content distributions between two emitters of a subsurface drip irrigation system. Soil Science Society of America Journal 75(2):488.
• Legates DR, Mccabe GJ (1999) Evaluating the use of “goodness-of-fit” measures in hydrologic and hydroclimatic model validation. Water Resources Research 35(1):233-241.
• Li X, Shi H, Šimůnek J, Gong X, Peng Z (2015) Modeling soil water dynamics in a drip-irrigated intercropping field under plastic mulch. Irrigation Science 33(4):289-302.
• Naglič B, Kechavarzi C, Coulon F, Pintar M (2014) Numerical investigation of the influence of texture, surface drip emitter discharge rate and initial soil moisture condition on wetting pattern size. Irrigation Science 32(6):421-436.
• Rivera RNC (2004) Modelagem da dinâmica da água e do potássio na irrigação por gotejamento superficial. Tese Doutorado, Piracicaba, Escola Superior de Agricultura “Luiz de Queiroz” Universidade de São Paulo.
• Rivera RNC, Duarte SN, Miranda JH, Botrel TA (2006) Modelagem da dinâmica do potássio no solo sob irrigação por gotejamento: validação do modelo. Engenharia Agrícola 26(2):388-394.
• Samadianfard S, Sadraddini AA, Nazemi AH, Provenzano G, Kisi O (2012) Estimating soil wetting patterns for drip irrigation using genetic programming. Spanish Journal of Agricultural Research 10(4):1155.
• Santos MA (2011) Extração de água do solo por plantas de soja: modelagem hidrofísica e empírica. Dissertação Mestrado, Piracicaba, Escola Superior de Agricultura “Luiz de Queiroz”, Universidade de São Paulo.
• Šimůnek J, Bristow KL, Helalia SA (2016) The effect of different fertigation strategies and furrow surface treatments on plant water and nitrogen use. Irrigation Science 34(1):53-69,
• Siyal AA, Skaggs TH (2009) Measured and simulated soil wetting patterns under porous clay pipe sub-surface irrigation. Agricultural Water Management 96(6):893-904.
• Subbaiah R (2011) A review of models for predicting soil water dynamics during trickle irrigation. Irrigation Science 31(3):225-258.
• Subbauah R, Mashru HH (2013) Modeling for predicting soil wetting radius under point source surface trickle irrigation. Agricultural Engineering International: CIGR Journal 15(3):1-10.
• van Dam JC, Feddes RA (2000) Numerical simulation of infiltration, evaporation and shallow groundwater levels with the Richards equation. Journal of Hydrology 233(1):72-85.
• Yao WW, Ma XY, Li J, Parkes M (2011) Simulation of point source wetting pattern of subsurface drip irrigation. Irrigation Science 29:331-339.

Publication Dates

• Publication in this collection
  Jan-Feb 2019

History

• Received
  31 Mar 2016
• Accepted
  25 Aug 2017