MODELING SOIL WATER REDISTRIBUTION UNDER SURFACE DRIP IRRIGATION

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.


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, twodimensional, 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: Where, h is the matric potential [L]; K (h) is the unsaturated hydraulic conductivity [L T -1 ], due to the matric potential, C (h) is the specific water capacity function [L -1 ], 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): Engenharia Agrícola, Jaboticabal, v.39, n.1, p.55-64, jan./feb.2019 In which the indices indicate i the column of the mesh, j the line, and t the time, 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)].
  Where, p is the number of iterations.
Substituting [eq.( 3)] into [eq.( 2)], we have [eq.( 4)]: Developing the equation ( 4) and grouping the terms, we obtain: Where, Engenharia Agrícola, Jaboticabal, v.39, n.1, p.55-64, jan./feb.2019 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 ho 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: For the following irrigations we have: 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: 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).
The lower boundary was displaced in such a way that the influence of the irrigation water in this zone was null (Dirichlet condition).

 
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.
Where, qi 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 physicalwater characteristics of the soil and Table 2 shows the parameters of the retention curve.In order to simulate a dripper, Rivera (2004) used a serum calibrator previously calibrated for a flow rate of 3 L h -1 .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 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): Where, Oi are the observed data, obtained in the experimental tests; Pi 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 θ.
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.
Engenharia Agrícola, Jaboticabal, v.39, n.1, p.55-64, jan./feb.2019 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 3 cm -3 .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 3 cm -3 , 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.
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.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.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.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)

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

FIGURE 5 .
FIGURE 5. Observed and simulated soil water content 48 hours after the end of irrigation.
, 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 -3 .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.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.

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

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

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