- Email: [email protected]
Drip-Irriwater: Computer software to simulate soil wetting patterns under surface drip irrigation
Computers and Electronics in Agriculture 98 (2013) 183–192
Contents lists available at ScienceDirect
Computers and Electronics in Agriculture journal homepage: www.elsevier.com/locate/compag
Drip-Irriwater: Computer software to simulate soil wetting patterns under surface drip irrigation G. Arbat a,⇑, J. Puig-Bargués a, M. Duran-Ros a, J. Barragán b, F. Ramírez de Cartagena a a b
Department of Chemical and Agricultural Engineering and Technology, University of Girona, C. de Maria Aurèlia Capmany, 61, 17071 Girona, Spain Department of Agricultural and Forestry Engineering, University of Lleida, Av. de l’Alcalde Rovira Roure, 191, 25198 Lleida, Spain
a r t i c l e
i n f o
Article history: Received 3 January 2013 Received in revised form 5 August 2013 Accepted 11 August 2013
Keywords: Soil water modeling Drip irrigation design Graphical user interface
a b s t r a c t Drip irrigation provides greater efficiency in terms of water usage and energy. These factors are very important in light of the current competition for water resources between the various users, especially in the Mediterranean region due to water scarcity. The shape and dimensions of the volume of wet soil below the emitter are some of the most influential variables in the optimal design and management of drip irrigation systems. This paper presents Drip-Irriwater, a code for determining soil wetting patterns under a single emitter in order to suggest planning factors for a drip irrigation system. The code solves Richards’ equation using the finite difference method subject to appropriate boundary conditions for drip irrigation. The boundary conditions on the surface of the soil allow us to consider forming a pond under the emitter. The user interface enables users to simply enter the input parameters (discharge rate, soil type and horizons, irrigation time, initial water content and total simulation time). The code then displays the soil water content and pressure heads, enabling a visual distinction between the wetted radius and depth, the parameters required for subsequent drip irrigation design. The results obtained with Drip-Irriwater were compared with those obtained with HYDRUS and with results from field tests carried out on three different soil types. The soil water distribution, as well as the wetted radius and depth, calculated from the Drip-Irriwater and HYDRUS code, were very similar. Field tests conducted to verify the results of the code confirmed that Drip-Irriwater gave a good approximation of the wetted radius and depth, and could therefore be used to design a drip irrigation system. Ó 2013 Elsevier B.V. All rights reserved.
1. Introduction Water scarcity in many regions of the world, combined with the large amount of water used in agriculture, is promoting the adoption of more efficient irrigation practices. Micro-irrigation, which allows irrigation at a low rate, achieves high efficiency in terms of water usage and is easily automated, is being used more and more. For instance, in Spain in 1999, micro-irrigation was adopted on 17% of the total irrigated land, while in 2011, this increased to 48% (MAGRAMA, 2011). Nevertheless, the use of micro-irrigation alone does not guarantee high efficiency. It requires adequate design and conscientious management (Phene, 1995). The agronomic design must be performed prior to designing the hydraulic network. This initial stage of the design process consists of establishing the discharge of the emitters, the distance between them and the emitter volume in order to wet a specific soil volume according to the extension of the
⇑ Corresponding author. Tel.: +34 972 418 459; fax: +34 972 418 399. E-mail address: [email protected] (G. Arbat). 0168-1699/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.compag.2013.08.009
rooting system and the needs for irrigation water. This set of variables is critical in order to achieve optimal application efficiency. Soil type plays a major role in adequate irrigation planning. Carrying out field tests which consist of excavating a soil pit to measure the extension and depth of the wetted soil volume under an emitter, is probably the most reliable method for selecting the adequate emitter rate and separation between emitters (Keller and Bliesner, 1990; Battam et al., 2003). Nevertheless, field tests are costly, and it is sometimes difficult to visually distinguish the wetted soil from the dry soil. For this reason, they are rarely conducted. The application of empirical or physical models could be an alternative. Schwartzmass and Zur (1986) equations, obtained empirically, are easy to use and allow us to obtain the width and depth of the wetted soil volume. However, due to the great variation in the hydraulic properties of the soils, their predictions are not reliable (Keller and Bliesner, 1990). Amin and Ekhmaj (2006) and Malek and Peters (2011) developed similar equations using several different experimental datasets. Kandelous et al. (2008) developed several equations to characterise the wetting front from the subsurface emitter, applying dimensional analysis using experimental data from a single soil type and emitter flow.
184
G. Arbat et al. / Computers and Electronics in Agriculture 98 (2013) 183–192
The physical approach describing water movement in unsaturated soils can be adopted using Richards’ equation (Richards, 1931). There are different analytical solutions to this equation, which can be applied to drip irrigation, including the ones developed by Raats (1971), Warrick (1974), Philip (1984) or Revol et al. (1997). In order to analytically solve the Richards equation, certain hypotheses must be assumed concerning the shape of the source, the particular form of the soil hydraulic functions, or in some cases, the steady state conditions. Using analytical instead of numerical solutions is preferred due to their speed (Cook et al., 2003, 2006). Nevertheless, nowadays personal computers are powerful enough and so this is not an issue. On the other hand, numerical solutions enable us to treat flow domain, boundary conditions and soil properties more realistically than analytical solutions (Warrick, 2003). Numerous models to analyze water flow under drip irrigation have been presented in scientific papers over the last 40 years (Brandt et al., 1971; Van der Ploeg and Benecke, 1974; Armstrong and Wilson, 1983; Taghavi et al., 1984; Ghali, 1986; Lafolie et al., 1989; Ramírez de Cartagena, 1994; Annandale et al., 2003), but their lack of a user-friendly interface could explain their underutilisation in drip irrigation design. HYDRUS 2D/3D (Šimu˚nek et al., 2006), a general purpose code allowing the creation of adequate flow domains and boundary conditions to simulate soil water distribution under micro-irrigation, has been used extensively in scientific papers (e.g. Cote et al., 2003; Cook et al., 2003; Skaggs et al., 2004; Arbat et al., 2008; Arbat et al., 2013; Kandelous and Šimu˚nek, 2010), but again their usage in drip irrigation design is uncommon. It is probable that the expertise required to use a general purpose model such as HYDRUS 2D/3D, which requires defining the flow domain and boundary conditions, is not found in the majority of technicians and consultancy services who design irrigation systems. The lack of models to solve specific problems was mentioned by Bastiaanssen et al. (2007) as one of the reasons why soil water models are not used by most technicians. At present, there are two different software products that can provide the planning factors needed for micro-irrigation design in a user-friendly environment: (1) WetUp (Cook et al., 2003), which displays wetting patterns from a database for specific soil types and drip surface and subsurface rates. The results are based on the approximate method of Thorburn et al. (2003) for calculating the radial and vertical distances and they assume that the wetted perimeter can be approximated by an ellipse, in accordance with Hachum et al. (1976). They also compared the calculation method for wetted distances with Philip’s (1984) approximate analytical equations. (2) Neuro-Drip (Hinnell et al., 2010), which displays soil wetting patterns for subsurface drip irrigation, is based on an artificial neural network which was trained with numerical simulations carried out by HYDRUS-2D (Šimu˚nek et al., 2006). The present paper describes Drip-Irriwater, a new, user-friendly software package that calculates the soil wetting pattern in homogeneous or layered soils under surface drip irrigation. The code is based on the Richards’ equation numerical solution using the finite difference method. The code validation is also presented, comparing their results with the ones obtained with the widely used HYDRUS 2D/3D and with field tests conducted under different soil conditions. Finally, an illustrative example of its use is shown. The flexibility the code offers the user for the precise definition of the soil type, the possibilities provided to estimate the hydraulic functions, the ponded area considered under the emitter during irrigation, and the output format of the results are the main features, contributions and advantages of this newly developed code.
2. Material and methods 2.1. Theory Soil water infiltration and redistribution under a drip source can be studied using Richards’ equation. Assuming symmetry around the vertical axis passing through the emitter, the equation can be expressed as:
@h @ @h @ @H þ K K ¼ @t @x @x @z @z
ð1Þ
where h is the volumetric soil water content, t is time, x and z are the respective horizontal and vertical co-ordinates, H is the hydraulic head, h represents the pressure head and K the hydraulic conductivity function, depends on h. The Richards’ equation is applied to a cylinder of soil under the emitter using a similar scheme to that described in the cylindrical model of Brandt et al. (1971). The numerical procedure to solve Eq. (1) was based on the scheme applied previously by Van der Ploeg and Benecke (1974), Armstrong and Wilson (1983) and Ramírez de Cartagena (1994). The numerical procedure applied is characterised by its mathematical simplicity, which enabled the presence of different soil layers and the ponded region on the soil surface to be taken into account. It has also been demonstrated to be numerically stable. 2.2. Flow domain and boundary conditions The origin of the coordinates (x = 0 and z = 0) is placed at the centre of the source, the negative z coordinate being downward (Fig. 1). It is assumed that there is no flow beyond the region of interest, which is fixed at a radial distance of 300 cm and a vertical distance of 300 cm, assuming that the influence of the drip source at this distance is negligible. A ‘‘no flow’’ boundary condition at 300 cm depth represents an impediment to the water to infiltrate deeper into the soil profile, but the effect of this impermeable layer will be inconsequential in the region of interest near the emitter. The soil cylinder is divided into concentric rings, 5 cm wide (Di) and 5 cm deep (Dj), centred on the vertical axis passing through the emitter. As a result, this provides an equally spaced mesh of 300 300 cm divided into cells of 5 5 cm each (Fig. 1). The boundary condition at the soil surface was established in order to take into account the variable extension in time of the ponded region formed under the emitter. According to experimental observations, an approximately circular saturated ponded area develops in the vicinity of the emitter. At the start of irrigation, this area is very small and grows over time until it reaches a maximum size (Brandt et al., 1971; Arbat, 2006). In the developed model, this area is the zone where the water enters into the soil, which is assumed to be perfectly circular, and a function of time. Initially, when irrigation starts, the size of the boundary condition is the same as the size of a cell (5 cm), and it will grow when more cells reach saturation. 2.3. Numerical procedure The differential Eq. (1) is transformed into a set of algebraic functions. The variables for a certain position and time take values associated with their position within the flow domain and boundary conditions described in previous section. The procedure applied to solve Eq. (1) is described below: The flow rate, Q (L3 T1) between two consecutive cells (Fig. 1) is calculated as follows:
Q ¼ qA
ð2Þ
185
G. Arbat et al. / Computers and Electronics in Agriculture 98 (2013) 183–192
X = j max
ρ (t)
Radius (j)
Δj
Boundaryconditions
Δi
No flow Evaporation Water t iinfiltration
1,jj ii-1 i,j-1 j
Z = i max
i,j
i,j+1
i+1,j Flow in a generic cell
Depth (i) Fig. 1. Chart of the flow domain, mesh, boundary conditions and water transfers in a Drip-Irriwater generic cell.
where q is the water flow density and A the surface area of the water, flowing from one cell to a consecutive one. The flow rate (Q) is calculated using Buckingham–Darcy’s law:
The change in volumetric soil water content in the ring corresponding to the generic cell (i, j) for a time step (Dt) is:
DH Q ¼ KðhÞ A DL
Dhi;j ¼ ð3Þ
where K(h) is the hydraulic conductivity corresponding to a water content value, DH the difference of hydraulic potential between two consecutive cells (H is defined as negative when the soil is unsaturated) and DL the distance between two cells. As described in the flow domain and boundary conditions section, the model assumes symmetrical flow with respect to the vertical axis passing through the emitter. Therefore, each cell can be associated to a ring with a depth equal to Di and a width equal to Dj. Considering the ring associated with a generic cell within the flow domain (i, j) and its four neighbouring cells (i, j 1), (i 1, j), (i, j + 1),(i + 1, j) (Fig. 1) and applying Eq. (3), the flow rates throughout the contours of the ring associated with the cell (i, j) are computed as: i;j Q 1 ¼ Q i;j i1;j ¼ K i1;j
i Hi;j Hi1;j h 2 p Rjþ1 R2j Di
ð4Þ
i;j Q 2 ¼ Q i;j i;jþ1 ¼ K i;jþ1
Hi;j Hi;jþ1 ð2pRjþ1 DiÞ Dj
ð5Þ
i Hi;j Hiþ1;j h 2 p Rjþ1 R2j Di
ð6Þ
Hi;j Hi;j1 ð2pRj DiÞ Dj
ð7Þ
Q3 ¼
Q i;j iþ1;j
¼
K i;j i1;j
i;j Q 4 ¼ Q i;j i;j1 ¼ K i;j1
ð9Þ
where Vi,j is the volume of the ring associated with cell i, j. Therefore, the soil water change (Dhi,j) is calculated as:
Dhi;j ¼
DQ i;j ðR2jþ1
p
R2j ÞDi
Dt
ð10Þ
The volumetric soil water content at ring i, j for the next time step (t + Dt) is:
htþDt ¼ ht þ Dhi;j
ð11Þ
All of these calculations are made for all of the cells in the flow domain, starting with the cell (1,1), closest to the emitter, and ending with the cell (imax, jmax). The water contents are updated for every cell. The same procedure must be repeated for every time step Dt until the final time specified in the simulation is reached. The boundary conditions of a certain cell depend on its position. Therefore, the flow domain was divided into nine different regions (Fig. 2) in order to calculate the net flow rate in the cells of each region:
For region 1 : i ¼ 1; 1 6 j 6 NðtÞ
where Q1 and Q4 are the vertical water flows through the ring associated with the cell i, j and Q2 and Q3 the horizontal flows through the same ring. K super-script is the mean hydraulic conductivity of the sub-script cells indicated in the subscript and superscript used for the water flow computations. Hsubscript is the hydraulic potential corresponding to the cell indicated in the subscript. Rj and Rj+1 are the interior and exterior radii of the ring associated with the generic cell (i, j). When the sign of the flow rate (Q) computed with Eqs. (4)–(7) is negative, the water flow leaves generic cell i, j. When it is positive, the water flow enters cell i, j. The net flow rate (DQ) in the generic cell (i, j) is given as:
DQ ¼ Q 1 þ Q 2 þ Q 3 þ Q 4
DQ i;j Dt V i;j
ð8Þ
X = jmax
ρ (t) 4
7
2
5
8
3
6
9
1.a
1b
Radius (j)
Z = i max Depth (i) Fig. 2. Regions of the flow domain to define the particular boundary conditions in each group of cells.
186
G. Arbat et al. / Computers and Electronics in Agriculture 98 (2013) 183–192
where N(t) is the number of saturated cells. At the beginning of the simulation, this value is equal to 1. The net flow rate in region 1 is calculated as:
DQ i;j ¼
Qe i;j þ Q i;j i;jþ1 þ Q iþ1;j NðtÞ
ð12Þ
For region 2: 1 < i < imax, j = 1 i;j i;j DQ i;j ¼ Q i;j i1;j þ Q i;jþ1 þ Q iþ1;j
ð13Þ
For region 3: i = imax, j = 1 i;j DQ i;j ¼ Q i;j i1;j þ Q i;jþ1
ð14Þ
For region 4: i = 1, 1 < j < jmax i;j i;j DQ i;j ¼ Q i;j i;jþ1 þ Q iþ1;j þ Q i;j1
ð15Þ
For region 5: 1 < i < imax, 1 < j < jmax i;j i;j i;j DQ i;j ¼ Q i;j i1;j þ Q i;jþ1 þ Q iþ1;j þ Q i;j1
ð16Þ
For region 6: 1 = imax, 1 < j < jmax i;j i;j DQ i;j ¼ Q i;j i1;j þ Q i;jþ1 þ Q i;j1
ð17Þ
For region 7: i = 1, j = jmax i;j DQ i;j ¼ Q i;j iþ1;j þ Q i;j1
ð18Þ
For region 8: 1 < i < imax, j = jmax i;j i;j DQ i;j ¼ Q i;j i1;j þ Q iþ1;j þ Q i;j1
ð19Þ
For region 9: i = imax, j = jmax i;j DQ i;j ¼ Q i;j i1;j þ Q i;j1
ð20Þ
2.4. The Drip-Irriwater code The Drip-Irriwater code is composed of two separate software programs: the graphical user interface (GUI) and the program to perform the numerical procedure described in the numerical procedure section.
The graphical user interface (GUI) was programmed in C#, using the Microsoft.NET framework. It was designed to introduce the input variables and to display the results (soil water and hydraulic potential distribution). The program to perform the numerical procedure was programmed in FORTRAN, using Compaq Digital Visual Fortran 6.0. The code was registered in 2010 with the number 02/2010/ 3036 on the General Registry of Intellectual Property in Spain. Drip-Irriwater software can be obtained from the authors upon request. 2.5. Input and output parameters In the GUI, the input parameters are classified in three groups (Fig. 3): (a) Soil input parameters: the number of soil layers, the model to account for the soil hydraulic functions and the initial water content in each soil layer. A maximum of 8 different soil layers can be defined by the user. The depth of each layer must be a multiple of 5 (the depth of a cell defined in the code) and can have different associated hydraulic functions and initial soil water contents. There are two different options to account for the soil hydraulic functions (h(h) and K(h)) in the Drip-Irriwater code. These are the Saxton model (Saxton et al., 1986) and the van Genuchten– Mualem model (van Genuchten, 1980). When the Saxton model is selected to define the soil hydraulic functions, the user must define the percentage of sand, clay and loam for each particular layer. On the other hand, when the van Genuchten–Mualem model is selected, the user must select the textural class and then the soil hydraulic parameters will be assigned according to the mean values presented by Carsel and Parrish (1988). The user will then have the option to modify some or all of the proposed values. (b) Irrigation parameters: the emitter discharge and irrigation time. (c) Simulation parameters: the simulation time and the time step in the numerical solution (1 s advised).
Fig. 3. Example of the data input interface in Drip-Irriwater.
187
G. Arbat et al. / Computers and Electronics in Agriculture 98 (2013) 183–192
Once the input parameters are specified using the GUI, they must be exported to be written in files. After that, the program to perform the numerical procedure must be executed from the GUI and, when the calculations are finished, the results (distribution of soil water content and hydraulic heads) can be imported and displayed using the GUI (Fig. 4). These output results are displayed numerically and graphically in a window. The cells where these variables have increased after irrigation are automatically coloured in different intensities of blue, according to the value reached by these variables, in order to help identify the soil volume wetted by the emitter more easily (Fig. 4). When multi-layered soil is simulated, high water content has been observed in the cells at the interface between the two horizons, due to the abrupt change in soil properties. In this case and in this position, the cell coloration is less useful in detecting the wetted volume.
characteristics of the different horizons (texture, apparent density and water content of the soil at 33 and 1500 kPa) for each soil profile (Table 1). The bulk density was measured from three undisturbed soil samples for each horizon of each profile and the soil hydraulic characterisation was performed with an equal number of samples and measured in the laboratory with a Richard’s pressure plate apparatus. The soils were classified as Typic Calcixerepts, Typic Haploxeralfs and Oxyaquic Xerofluvents (SSS, 2010). The hydraulic functions of the soils follow the van Genuchten– Mualem model (VG–M). The parameters of these equations were determined from pedotransfer functions using the ROSETTA Lite (Schaap et al., 2001) computer program. The input parameters were the percentages of clay, silt and sand, the apparent density, and the water content at 33 and 1500 kPa (Table 1). The parameters of the VG–M equations are presented in Table 2. 3.2. Irrigation trials
3. Field studies 3.1. Soil characterisation and hydraulic properties In order to verify the model in the field, a series of tests were carried out during the summers of 2001 and 2002 on three representative agricultural soil types. Two of them were located in Monells and the other in Canet de la Tallada (all located in Girona, Spain). A soil trench (3 m long, 0.60 m wide, 2 m deep) was opened at each of the different field sites to allow a detailed pedological description of the soil profile. Determinations were made of the
In each of the three soils, the irrigation tests were carried out with a discharge of 2, 4, 8.5 and 25 L h1, respectively, using an isolated emitter, in each case applying a total volume of 25 L. In order to identify the water input zones below the emitters, the sizes of the pools were measured throughout the testing period. Each test determined the water content before, during and after irrigation with a portable TDR probe (IMKO TRIME-T) and 1.6 m long access tubes placed at horizontal distances of 5, 20, 35, 50, and 75 cm in two perpendicular directions away from the emitter and at 100 and 150 cm in only one of the directions. The position of the tubes
Fig. 4. Example of output results (volumetric soil water content) obtained with Drip-Irriwater.
Table 1 Physical characteristics of different horizons of soils used in the experiments. Soil type and field site
Horizon
Depth (cm)
Clay (%)
Silt (%)
Sand (%)
Bulk density (g cm3)
Volumetric water content at 33 kPa (%)
Volumetric water content at 1500 kPa (%)
Typic Calcixerepts, Monells
Ap Bkn C
0–35 35–80 >80
20.3 20.6 17.0
45.4 45.8 66.5
34.3 33.6 16.5
1.38 1.40 1.42
36.57 23.84 41.60
14.31 13.41 12.80
Typic Haploxeralfs, Monells
Ap Bw1 Bt2
0–20 20–50 >50
23.1 25.3 34.8
28.4 28.3 27.6
48.5 46.4 37.6
1.65 1.59 1.63
29.70 25.44 35.86
14.85 12.72 19.56
Oxyaquic Xerofluvents, Canet de la Tallada
Ap
0–35
11.54
38.87
49.59
1.65
26.83
11.73
Bw1 Bw2 Bw3
35–75 75–120 >120
12.64 4.94 7.01
48.84 26.00 17.34
38.52 69.06 75.65
1.53 1.44 1.68
28.21 30.79 –
12.07 13.95 –
188
G. Arbat et al. / Computers and Electronics in Agriculture 98 (2013) 183–192
Table 2 Parameters of the van Genuchten–Mualem equations for the three different soil types and corresponding horizons. Soil type and field site
Horizon
Depth (cm)
hr (cm3 cm3)
hs (cm3 cm3)
a (cm1)
n
Ks (cm h1)
Typic Calcixerepts, Monells
Ap Bkn C
0–35 35–80 >80
0.0632 0.0553 0.0847
0.4311 0.4035 0.4550
0.0023 0.0273 0.0010
1.7285 1.3377 2.3304
0.52 0.97 0.41
Typic Haploxeralfs, Monells
Ap and Bw1 Bt2
0–50 >50
0.0467 0.0569
0.3691 0.3968
0.098 0.0039
1.3476 1.3779
0.28 0.07
Oxyaquic Xerofluvents, Canet de la Tallada
Ap and Bw1 Bw2 Bw3
0–75 75–120 >120
0.0399 0.0424 0.0384
0.3593 0.4013 0.3354
0.0071 0.0071 0.0460
1.4200 1.4188 1.5295
0.47 1.12 1.73
is shown in Fig. 5. Measurements were taken at every 15 cm from 20 to 140 cm soil depth. The probe readings were corrected using a specific calibration for each of the horizons.
4. Model validation 4.1. Comparison of the results obtained with Drip-Irriwater and HYDRUS HYDRUS (Šimu˚nek et al., 2006), a well-known general purpose commercial software for simulating soil water movement and solute transport in variable-saturated media, was used to assess the soil water distribution under an isolated emitter with the Drip-Irriwater code. In order to compare the results obtained with both codes, the soil water content was simulated on a soil of sandy–loam texture, which was irrigated with an emitter discharge of 4 L h1 over 3 h and 17 min. The water was then left to be redistributed for 12 h from the start of irrigation. The initial soil water content was considered to be 0.184 cm3 cm3 throughout the entire flow domain. The hydraulic properties of the soil were estimated using van Genuchten–Mualem functions, and the parameters were selected from a textural class database incorporated in HYDRUS and DripIrriwater (Carsel and Parrish, 1988), corresponding to a sandy– loam soil. The soil water infiltration was considered to be axisymmetric, and the size of the computational flow domain was 100 cm width in the horizontal and 150 cm deep. This was considered to be large enough to prevent interferences in the water flow from the emitter and not so large as to increase the simulation time unnecessarily. The rectangular flow region (100 cm wide 150 cm depth) was discretised into a structured triangular finite element mesh. The number of nodes on the horizontal and vertical sides of the rectangular region was 50 and 70 respectively. The element sizes were then distributed in the rectangular domain resulting in
triangular finite elements of approximately 2 cm. Nevertheless the size of the elements at the top of the domain was defined to be slightly smaller than the ones at the bottom because more active flow was expected. Initially the size of the elements in HYDRUS was selected to be 5 cm, to match the mesh size used by the finite differences method used in Drip-Irriwater, but these large elements produced problems during the numerical solution. Therefore, a smaller size of the elements was defined. No-flow boundary conditions were imposed on the right, left and bottom sides of the computational flow domain in the same way as the boundary conditions used by Drip-Irriwater. These boundary conditions have been reported in previous studies to simulate water infiltration from an isolated emitter using HYDRUS (Kandelous and Šimu˚nek, 2010). As the standard version of HYDRUS-2D (Šimu˚nek et al., 2006) does not have a dynamic boundary condition implemented that is suitable for simulating variable ponding at the soil surface beneath the drip emitter during irrigation, the water entry at the soil surface was simulated as a constant water flux within an area of 30 cm from the upper right corner of the flow domain, according to the observed saturated radius in similar textural soil characteristics with an emitter discharge rate of 4 L h1. The soil water contents computed with Drip-Irriwater at 5, 15, 25, 35, 45 and 55 cm horizontal distances from the emitter and at 5, 15, 25, 35, 45, 55 and 60 cm depths were visually compared with those at the same positions obtained with HYDRUS and represented in a two dimensional plot of iso-water contents interpolating the values between these points using the Kriging method. The comparison was also performed between the soil water content vertically at 30 cm distance from the dripper and laterally at 30 cm depth from the soil surface. In order to determine the agreement of both models the root mean square error (RMSE) and the coefficient of determination (r2) were computed in the same grid of points as follows:
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Pn 0 2 i¼1 ðP i P i Þ RMSE ¼ n
75 50
ð21Þ
35 20
12 0 0 P P P P i C B i i¼1 C rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r2 ¼ B @rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 A Pn 0 0 2 Pn i¼1 P i P i¼1 P i P 0
Pn
ð22Þ
Emitter 20
35
50
100
75
150
Y
5
X
Fig. 5. Chart showing the position of the TDR access tubes for the measurement of the soil water content in the field tests.
where Pi is soil water content predictedby Drip-Irriwater, P0i the soil water content obtained with HYDRUS, P the average soil water con tent predicted with Drip-Irriwater, P 0 the average soil water content obtained with HYDRUS, and n the number of observations. The units of the RMSE are cm3 cm3.
189
G. Arbat et al. / Computers and Electronics in Agriculture 98 (2013) 183–192
4.2. Comparison between simulation and experimental results
drip irrigation experiments (Skaggs et al., 2004; Arbat et al., 2008; Kandelous and Šimu˚nek, 2010).
Simulated and measured soil water contents, and the maximum radius (R) and depth (D) of the wetted soil volume were compared at the end of irrigation and at 24 h after the end of irrigation. R and D were computed using the criteria that soil water content was 5% higher than initially measured. As the computational mesh size was 5 5 cm, R and D were always multiples of 5. Simulated and measured soil water contents were compared at all the measured points by calculating the RMSE and r2 in each of the tests at the end of irrigation and after 24 h. The calculations were made with Eqs. (21) and (22), in this case O indicates the measured water contents, and P the water contents with predicted Drip-Irriwater. The measured soil water content from 0 to 75 cm was the average soil water content measured in two perpendicular directions (Fig. 5). r2 indicates the degree of error variance between two values. Higher r2 value means less error variance between the observed and predicted values. The RMSE has the advantage of expressing the error with the same units as the variable. Both statistics were selected because they had previously been used to compare predicted and measured soil water contents in different
5. Results and discussion 5.1. Comparison between Drip-Irriwater estimations and HYDRUS The comparison of soil water distribution obtained from DripIrriwater and HYDRUS showed a similar moisture distribution pattern (Fig. 6). When comparing the simulated water contents within the soil wetted region (55 cm depth and 55 cm from the emitter source) obtained with both models, the RMSE was 0.018 cm3 cm3 and r2 was 0.86. Fig. 7 shows a comparison between the soil water content vertically at selected positions of 30 cm distance from the emitter and laterally at 30 cm depth from the soil surface. The RMSEs comparing the water contents were 0.007 cm3 cm3 for the vertical line at 30 cm from the emitter, and 0.013 cm3 cm3 for the horizontal line 30 cm deep, while r2 for the same positions was 0.931 and 0.934, respectively. The statistics comparing the simulated water contents within the soil wetted region and the ones obtained in selected positions indicate small differences in
Fig. 6. Soil water content distribution at 12 h after starting irrigation with a drip flow of 4 L h1 and irrigation time of 3 h 17 min, with HYDRUS (a) and Drip-Irriwater (b).
Vertical sección at 30cm from the emitter 0.00 0
Depth (cm) D
10 20
0.05
0.10
HYDRUS Drip-Irriwater
30 40 50 60 70
RMSE = 0.007 (cm3/cm3) R2 = 0.931
0.15
0.20
0.25
0.30
Vo olume etric s soil wa ater c conten nt (cm m 3 cm m -3)
Volumetric soil water content (cm 3 cm -3)
Horitzontal section at 30 cm depth
0.35
HYDRUS
0.30
Drip-Irriwater
0.25 0. 20 0 15 0.
RMSE = 0.013 0 013 (cm3/cm3) R 2 = 0.934
0.10 0.05 0.00 0
5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90
Horizontal distance (cm) Fig. 7. Comparison of soil water content between HYDRUS (red dots) and Drip-Irriwater (blue dots): (i) vertically at 30 cm away from the emitter and (ii) horizontally at 30 cm depth from the soil surface corresponding a drip flow of 4 L h1 for an irrigation time of 3 h and 17 min 12 h after stopping irrigation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
190
G. Arbat et al. / Computers and Electronics in Agriculture 98 (2013) 183–192
the soil content at certain points, which resulted in part from the different dimensions of the computational mesh used by the two codes. While in Drip-Irriwater the mesh was composed of square cells of 5 5 cm, in HYDRUS the triangular finite elements had vertical and horizontal sides ranging from 1.5 to 2 cm, allowing greater resolution in the resulting soil water distribution. The different boundary conditions, used to simulate the soil volume where the water infiltrated under the emitter, also explains the slightly different results obtained in the two computer programs. With an accuracy of 5 cm, equal to the length of the cells used in Drip-Irriwater, the radius of the wetted soil volume at a depth of 30 cm was approximately 40–45 cm and the maximum wetted depth was 55–60 cm, irrespective of the software used (Fig. 6). 5.2. Comparison of the field tests and Drip-Irriwater results The comparison between the water contents simulated with the Drip-Irriwater code and the ones measured at the end of irrigation (Table 3) shows that the average RMSE was 0.028 cm3 cm3 (with a minimum value of 0.019 and a maximum of 0.035) and an average r2 of 0.54 (with a minimum value of 0.27 and a maximum of 0.83). Twenty-four hours after completion of irrigation (Table 4) the average RMSE was 0.026 cm3 cm3 (with a minimum value of 0.015 and maximum of 0.040) and the average r2 was 0.56 (with a minimum value of 0.45 and maximum of 0.75). The differences between the water contents in the simulations and in the tests could be partly due to the errors inherent in the method of measurement. In the calibrations performed on the same sites where tests were performed, comparing soil water contents measured with TDR and soil samples, the RMSE was between 0.030 cm3 cm3
and 0.050 cm3 cm3, and the coefficient of determination (r2) between 0.57 and 0.89 (Arbat, 2006). This shows that sometimes the measurement method could affect the determination of the wetted radius and depth, which explains the significant differences found between measured and simulated values in some comparisons (Table 3). The Drip-Irriwater code simulates a slightly greater wetted radius than the wetted radius determined from TDR measurements, with a mean difference of 7.3 cm at the end of irrigation and 4.5 cm after 24 h of redistribution (Tables 3 and 4). Similarly the wetted depth determined with the Drip-Irriwater code was smaller than the wetted depth determined from measurements with TDR, with a mean difference of 10 cm at the end of irrigation and 2.7 cm after 24 h of redistribution. The differences between the values of the wetted radius and depth in the simulated and those determined from measurements with TDR were lower 24 h after the end of irrigation. It is considered that the differences between the values of the wetted radius and depth were smaller overall and the Drip-Irriwater code allowed us to estimate these parameters appropriately for the agronomic design of a drip irrigation system. Kandelous and Šimu˚nek (2010) compared the extent of wetted soil volume, using analytical and empirical models, with the one obtained with the numerical model HYDRUS, concluding that the predictions of the latter are more realistic. However, in surface drip irrigation, HYDRUS may underestimate the wetted horizontal length and overestimate the wetted depth when pounded area is formed on the soil surface. This is because the static boundary condition defined in HYDRUS does not reflect the extension of the pounded area. In this regard, the
Table 3 Statistics of comparison between measured and simulated soil water contents, and values of measured and simulated wetted radii and depth with an irrigation volume of 25 L at the end of the irrigation event. Soil type
Flow rate (L h1)
RMSE (cm3 cm3)
r2
Measured radii (cm)
Simulated radii (cm)
Measured depth (cm)
Simulated depth (cm)
Typic Calcixerepts
2.0 4.0 8.5 25.0
0.035 0.021 0.028 0.033
0.36 0.83 0.67 0.49
25 40 40 40
45 45 40 40
50 40 45 40
35 35 35 30
Typic Haploxeralfs
2.0 4.0 8.5 25.0
0.034 0.034 0.028 0.019
0.67 0.59 0.61 0.43
45 30 30 20
40 40 40 40
40 30 25 30
35 25 20 15
Oxyaquic Xerofluvents
4.0 8.5 25.0
0.029 0.023 0.027
0.46 0.58 0.27
35 40 35
45 45 40
45 40 40
35 30 20
0.028
0.54
Average 2
RMSE: Root mean square error, r : coefficient of determination. Table 4 Statistics of comparison between the measured and simulated soil water contents, and values of measured and simulated wetted radii and depth with an applied irrigation volume of 25 L, 24 h after the end of irrigation. Soil type
Flow rate (L h1)
RMSE (cm3 cm3)
r2
Measured radii (cm)
Simulated radii (cm)
Measured depth (cm)
Simulated depth (cm)
Typic Calcixerepts
2.0 4.0 8.5 25.0
0.026 0.015 0.018 0.018
0.47 0.63 0.50 0.51
35 45 40 40
50 50 45 50
50 40 45 35
35 35 35 35
Typic Haploxeralfs
2.0 4.0 8.5 25.0
0.027 0.035 0.040 0.034
0.74 0.67 0.59 0.67
45 45 40 40
45 40 45 45
45 45 55 50
50 40 45 45
Oxyaquic Xerofluvents
4.0 8.5 25.0
0.026 0.024 0.018
0.43 0.45 0.52
60 40 40
50 55 45
50 45 40
60 50 40
0.026
0.56
Average
RMSE: Root mean square error, r2: coefficient of determination.
G. Arbat et al. / Computers and Electronics in Agriculture 98 (2013) 183–192
191
Fig. 8. Input and output GUI of Drip-Irriwater for 2 soils (loamy sand and sandy) irrigated with a drip discharge of 4 L h1 for 2 h.
Drip-Irriwater code adds the effect of the pond as a dynamic boundary condition to the Richards equation. As described in the Numerical Procedure section, the horizontal extension of the ponded area is recalculated for each time step, increasing the extension of the infiltration zone when a new cell is saturated (region 1 in Fig. 2). Therefore the Drip-Irriwater code handles a dynamic boundary condition for the infiltration zone under the emitter, while in the standard version of HYDRUS this condition is fixed by the user and kept constant throughout the simulation. The use of this dynamic boundary condition in the Drip-Irriwater code is more realistic when the application rate is greater than the infiltration capacity of the soil. 5.3. Illustrative examples Two examples are given in this section in order to guide the user. Similar amount of irrigation water (8 L) was applied in both cases with a 4 L h1 flow rate emitter and the water content were only displayed after 2 h of irrigation in two different soil textures. In the first case, the soil was of loamy–sand texture and in the second it was of sandy texture. The soil hydraulic parameters inferred from the texture were used in the simulations. The initial soil water content was set up as 15.5% in both cases. Fig. 8 shows the input data for the loamy–sand texture, and the output soil water contents for both soils. It is worth noting that the soil texture plays an important role in the shape and size of the wetting patterns. In the sandy soil, the maximum wetted depth after 2 h was approximately 65 cm, clearly deeper than that of 50 cm, observed in the loamy–sand soil. 6. Conclusions The Drip-Irriwater code allows the simulation of soil water distribution under surface drip irrigation, considering the main
parameters involved in the shape and dimensions of the volume of wet soil such as: soil type, emitter flow and irrigation time, the initial soil water content and redistribution time after irrigation. The soil water content distribution obtained with Drip-Irriwater was compared with those obtained with the widely recognised HYDRUS code and the comparison statistics (RMSE, r2) indicate good agreement between both. The small differences observed may be partially due to the fact that Drip-Irriwater incorporates a dynamic boundary condition, which takes the ponded surface into account. This condition is not available in the standard version of HYDRUS. In particular, the design variables obtained with both codes (wetted radius and depth) were identical for the simulated case. The field verification of the model revealed an acceptable agreement between measured and simulated soil water content for surface drip irrigation design. Overall, the differences between the wetted radius and depth were low and could be partly explained by the errors inherent to the measurement method. It can be concluded that the Drip-Irriwater software is accessible to design engineers, irrigation technicians and practitioners, as the input parameters can be easily known or estimated. The model could also be used to estimate the soil planning factors (wetted radius and wetted depth) with sufficient approximation to efficiently design surface drip irrigation systems. Acknowledgments The authors would like to express their gratitude to the Institut de Recerca i Tecnologia Agroalimentàries (IRTA), specially to the centres of IRTA-Mas Badia in Canet de la Tallada, and IRTA-Finca Camps i Armet in Monells for their help in carrying out this investigation.
192
G. Arbat et al. / Computers and Electronics in Agriculture 98 (2013) 183–192
References Amin, M., Ekhmaj, A., 2006. DIPAC-drip irrigation water distribution calculator. In: 7th International Micro Irrigation Congress, 10–16 September, PWTC, Kuala Lumpur, Malaysia. Annandale, J.G., Jovanovic, N.Z., Campbell, G.S., Du Sautoy, N., Benadé, N., 2003. A two-dimensional water balance model for micro-irrigated hedgerow tree crops. Irrig. Sci. 22 (3–4), 157–170. Arbat, G., 2006. Desarrollo y validación de la dinámica del agua en el suelo. Aplicación al diseño agronómico y al manejo en riego localizado. Tesis Doctoral, Universitat de Lleida, Lleida, Spain. Arbat, G., Puig-Bargués, J., Barragán, J., Bonany, J., Ramírez de Cartagena, F., 2008. Monitoring soil water status for micro-irrigation management versus modelling approach. Biosyst. Eng. 100 (2), 286–296. Arbat, G., Roselló, A., Domingo Olivé, F., Puig-Bargués, J., González Llinàs, E., DuranRos, M., Pujol, J., Ramírez de Cartagena, F., 2013. Soil water and nitrate distribution under drip irrigated corn receiving pig slurry. Agric. Water Manage. 120, 11–22. Armstrong, C.F., Wilson, T.V., 1983. Computer model for moisture distribution in stratified soils under a trickle irrigation source. Trans. ASABE 26 (6), 1704–1709. Bastiaanssen, W.G.M., Allen, R.G., Droogers, P., D’urso, G., Steduto, P., 2007. Twentyfive years modeling irrigated and drained soils: State of the Art. Agric. Water Manage. 92 (3), 111–125. Battam, M.A., Boughton, D.G., Sutton, B.G., 2003. Soil pits as a simple design aid for subsurface drip irrigation systems. Irrig. Sci. 22 (3), 135–142. Brandt, A., Bresler, E., Diner, N., Ben Asher, J., Heller, J., Goldberg, D., 1971. Soil Sci. Soc. Am. Proc. 35 (5), 675–682. Carsel, R., Parrish, R., 1988. Developing joint probability of soil water retention characteristics. Water Resour. Res. 24 (5), 755–769. Cook, F.J., Thorburn, P.J., Fitch, P., Bristow, K.L., 2003. WetUp: a software tool to display approximate wetting patterns from drippers. Irrig. Sci. 22 (3–4), 129– 134. Cook, F.J., Thorburn, P.J., Fitch, P., Charlesworth, P.B., Bristow, K.L., 2006. Modelling trickle irrigation: comparison of analytical and numerical models for estimation of wetting front position with time. Environ. Model Soft. 21 (9), 1353–1359. Cote, C.M., Bristow, K.L., Charlesworth, P.B., Cook, F.J., Thornburn, P.J., 2003. Analysis of soil wetting and solute transport in subsurface trickle irrigation. Irrig. Sci. 22 (3–4), 143–156. Ghali, G.S., 1986. Mathematical Modelling of Soil Moisture Dynamics in Trickle Irrigated Fields. Ph.D. Thesis, University of Southampton. Hachum, A.Y., Alfaro, J.F., Willardson, L.S., 1976. Water movement in soil from a trickle source. J. Irrig. Drain. Div. Proc. ASCE 102 (IR2), 179–192. Hinnell, A.C., Lazarovitch, N., Furman, A., Poulton, M., Warrick, A.W., 2010. NeuroDrip: estimation of subsurface wetting patterns for drip irrigation using neural networks. Irrig. Sci. 28 (6), 535–544. Kandelous, M., Liagahat, A., Abbasi, F., 2008. Estimation of soil moisture pattern in subsurface drip irrigation using dimensional analysis method. J. Agric. Sci. 39 (2), 371–378. Kandelous, M., Šimu˚nek, J., 2010. Comparison of numerical, analytical, and empirical models to estimate wetting patterns for surface and subsurface drip irrigation. Irrig. Sci. 28 (5), 435–444. Keller, J., Bliesner, R.D., 1990. Trickle irrigation planning factors. In: Van Nostrand Reinhold (Ed.), Sprinkle and Trickle Irrigation. Avi Book, New York, pp. 453–477.
Lafolie, F., Guennelon, R., van Genuchten, M.T., 1989. Analysis of water flow under trickle irrigation: I. Theory and numerical solution. Soil Sci. Soc. Am. J. 53 (5), 1310–1318. MAGRAMA (Ministerio de Agricultura, Alimentación y Medio Ambiente de España), 2011. Encuesta sobre superficies y rendimientos de cultivos. Anàlisis de los regadíos españoles. Secretaría General Técnica, Madrid. Malek, K., Peters, R.T., 2011. Wetting pattern models for drip irrigation: new empirical model. J. Irrig. Drain. Eng. 137 (8), 530–536. Phene, C.J., 1995. Research trends in microirrigation. Microirrigation for a changing world: conserving resources/preserving the environment. In: Proceedings of the Fifth Microirrigation Congress, 2–6 April, Orlando, Florida. Philip, J.R., 1984. Travel times from buried and surface infiltration point sources. Water Resour. Res. 20 (7), 990–994. Raats, P.A.C., 1971. Steady infiltration from point sources, cavities and basins. Soil Sci. Soc. Am. Proc. 35, 689–694. Ramírez de Cartagena, F., 1994. Simulación numérica de la dinámica del agua en el suelo. Aplicación al diseño de sistemas de riego LAF. Tesis Doctoral, Universitat de Lleida, Lleida, Spain. Revol, P., Clothier, B.E., Mailhol, J.C., Vachaud, G., Vauclin, M., 1997. Infiltration from a surface point source and drip irrigation. 2. An approximate time-dependent solution for a wet-front position. Water Resour. Res. 33 (8), 1869–1874. Richards, L.A., 1931. Capillary conduction of liquids through porous mediums. Physics 1, 318–333. SSS (Soil Survey Staff), 2010. Keys to Soil Taxonomy, 11th ed. USDA, Natural Resources Conservation Service. Saxton, K.E., Rawls, W.J., Romberger, J.S., Papendick, R.I., 1986. Estimating generalized soil–water characteristics from texture. Soil Sci. Soc. Am. J. 50 (4), 1031–1036. Schaap, M.G., Leij, F.J., van Genuchten, M.T., 2001. ROSETTA: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol. 251 (3–4), 163–176. Schwartzmass, M., Zur, B., 1986. Emitter spacing and geometry of wetted soil volume. J. Irrig. Drain. Eng. 112 (3), 242–253. Šimu˚nek, J., van Genuchten, M.T., Sejna, M., 2006. HYDRUS, software package for simulating the two- and three-dimensional movement of water, heat and solute transport in variably-saturated media. Technical Manual v. 1. PC Progress, Prague, Czech Republic. Skaggs, T.H., Trou, T., Šimu˚nek, J., Shouse, P.J., 2004. Comparison of HYDRUS-2D simulations of drip irrigation with experimental observations. J. Irrig. Drain. Eng. 130 (4), 304–310. Taghavi, S.A., Mariño, M.A., Rolston, D.E., 1984. Infiltration from trickle irrigation source. J. Irrig. Drain. Eng. 110 (4), 331–341. Thorburn, P.J., Cook, F.J., Bristow, K.L., 2003. Soil dependent wetting from trickle emitters: implication for design and management. Irrig. Sci. 22 (3–4), 121–127. Van der Ploeg, R.R., Benecke, P., 1974. Unsteady, unsaturated, n-dimensional moisture flow in soil: a computer simulation program. Soil Sci. Soc. Am. Proc. 38 (6), 881–885. van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892–898. Warrick, A.W., 1974. Time-dependent linearized infiltration: I. Point sources. Soil Sci. Soc. Am. Proc. 38 (3), 383–386. Warrick, A.W., 2003. Soil Water Dynamics. Oxford University Press, New York.
Contents lists available at ScienceDirect
Computers and Electronics in Agriculture journal homepage: www.elsevier.com/locate/compag
Drip-Irriwater: Computer software to simulate soil wetting patterns under surface drip irrigation G. Arbat a,⇑, J. Puig-Bargués a, M. Duran-Ros a, J. Barragán b, F. Ramírez de Cartagena a a b
Department of Chemical and Agricultural Engineering and Technology, University of Girona, C. de Maria Aurèlia Capmany, 61, 17071 Girona, Spain Department of Agricultural and Forestry Engineering, University of Lleida, Av. de l’Alcalde Rovira Roure, 191, 25198 Lleida, Spain
a r t i c l e
i n f o
Article history: Received 3 January 2013 Received in revised form 5 August 2013 Accepted 11 August 2013
Keywords: Soil water modeling Drip irrigation design Graphical user interface
a b s t r a c t Drip irrigation provides greater efficiency in terms of water usage and energy. These factors are very important in light of the current competition for water resources between the various users, especially in the Mediterranean region due to water scarcity. The shape and dimensions of the volume of wet soil below the emitter are some of the most influential variables in the optimal design and management of drip irrigation systems. This paper presents Drip-Irriwater, a code for determining soil wetting patterns under a single emitter in order to suggest planning factors for a drip irrigation system. The code solves Richards’ equation using the finite difference method subject to appropriate boundary conditions for drip irrigation. The boundary conditions on the surface of the soil allow us to consider forming a pond under the emitter. The user interface enables users to simply enter the input parameters (discharge rate, soil type and horizons, irrigation time, initial water content and total simulation time). The code then displays the soil water content and pressure heads, enabling a visual distinction between the wetted radius and depth, the parameters required for subsequent drip irrigation design. The results obtained with Drip-Irriwater were compared with those obtained with HYDRUS and with results from field tests carried out on three different soil types. The soil water distribution, as well as the wetted radius and depth, calculated from the Drip-Irriwater and HYDRUS code, were very similar. Field tests conducted to verify the results of the code confirmed that Drip-Irriwater gave a good approximation of the wetted radius and depth, and could therefore be used to design a drip irrigation system. Ó 2013 Elsevier B.V. All rights reserved.
1. Introduction Water scarcity in many regions of the world, combined with the large amount of water used in agriculture, is promoting the adoption of more efficient irrigation practices. Micro-irrigation, which allows irrigation at a low rate, achieves high efficiency in terms of water usage and is easily automated, is being used more and more. For instance, in Spain in 1999, micro-irrigation was adopted on 17% of the total irrigated land, while in 2011, this increased to 48% (MAGRAMA, 2011). Nevertheless, the use of micro-irrigation alone does not guarantee high efficiency. It requires adequate design and conscientious management (Phene, 1995). The agronomic design must be performed prior to designing the hydraulic network. This initial stage of the design process consists of establishing the discharge of the emitters, the distance between them and the emitter volume in order to wet a specific soil volume according to the extension of the
⇑ Corresponding author. Tel.: +34 972 418 459; fax: +34 972 418 399. E-mail address: [email protected] (G. Arbat). 0168-1699/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.compag.2013.08.009
rooting system and the needs for irrigation water. This set of variables is critical in order to achieve optimal application efficiency. Soil type plays a major role in adequate irrigation planning. Carrying out field tests which consist of excavating a soil pit to measure the extension and depth of the wetted soil volume under an emitter, is probably the most reliable method for selecting the adequate emitter rate and separation between emitters (Keller and Bliesner, 1990; Battam et al., 2003). Nevertheless, field tests are costly, and it is sometimes difficult to visually distinguish the wetted soil from the dry soil. For this reason, they are rarely conducted. The application of empirical or physical models could be an alternative. Schwartzmass and Zur (1986) equations, obtained empirically, are easy to use and allow us to obtain the width and depth of the wetted soil volume. However, due to the great variation in the hydraulic properties of the soils, their predictions are not reliable (Keller and Bliesner, 1990). Amin and Ekhmaj (2006) and Malek and Peters (2011) developed similar equations using several different experimental datasets. Kandelous et al. (2008) developed several equations to characterise the wetting front from the subsurface emitter, applying dimensional analysis using experimental data from a single soil type and emitter flow.
184
G. Arbat et al. / Computers and Electronics in Agriculture 98 (2013) 183–192
The physical approach describing water movement in unsaturated soils can be adopted using Richards’ equation (Richards, 1931). There are different analytical solutions to this equation, which can be applied to drip irrigation, including the ones developed by Raats (1971), Warrick (1974), Philip (1984) or Revol et al. (1997). In order to analytically solve the Richards equation, certain hypotheses must be assumed concerning the shape of the source, the particular form of the soil hydraulic functions, or in some cases, the steady state conditions. Using analytical instead of numerical solutions is preferred due to their speed (Cook et al., 2003, 2006). Nevertheless, nowadays personal computers are powerful enough and so this is not an issue. On the other hand, numerical solutions enable us to treat flow domain, boundary conditions and soil properties more realistically than analytical solutions (Warrick, 2003). Numerous models to analyze water flow under drip irrigation have been presented in scientific papers over the last 40 years (Brandt et al., 1971; Van der Ploeg and Benecke, 1974; Armstrong and Wilson, 1983; Taghavi et al., 1984; Ghali, 1986; Lafolie et al., 1989; Ramírez de Cartagena, 1994; Annandale et al., 2003), but their lack of a user-friendly interface could explain their underutilisation in drip irrigation design. HYDRUS 2D/3D (Šimu˚nek et al., 2006), a general purpose code allowing the creation of adequate flow domains and boundary conditions to simulate soil water distribution under micro-irrigation, has been used extensively in scientific papers (e.g. Cote et al., 2003; Cook et al., 2003; Skaggs et al., 2004; Arbat et al., 2008; Arbat et al., 2013; Kandelous and Šimu˚nek, 2010), but again their usage in drip irrigation design is uncommon. It is probable that the expertise required to use a general purpose model such as HYDRUS 2D/3D, which requires defining the flow domain and boundary conditions, is not found in the majority of technicians and consultancy services who design irrigation systems. The lack of models to solve specific problems was mentioned by Bastiaanssen et al. (2007) as one of the reasons why soil water models are not used by most technicians. At present, there are two different software products that can provide the planning factors needed for micro-irrigation design in a user-friendly environment: (1) WetUp (Cook et al., 2003), which displays wetting patterns from a database for specific soil types and drip surface and subsurface rates. The results are based on the approximate method of Thorburn et al. (2003) for calculating the radial and vertical distances and they assume that the wetted perimeter can be approximated by an ellipse, in accordance with Hachum et al. (1976). They also compared the calculation method for wetted distances with Philip’s (1984) approximate analytical equations. (2) Neuro-Drip (Hinnell et al., 2010), which displays soil wetting patterns for subsurface drip irrigation, is based on an artificial neural network which was trained with numerical simulations carried out by HYDRUS-2D (Šimu˚nek et al., 2006). The present paper describes Drip-Irriwater, a new, user-friendly software package that calculates the soil wetting pattern in homogeneous or layered soils under surface drip irrigation. The code is based on the Richards’ equation numerical solution using the finite difference method. The code validation is also presented, comparing their results with the ones obtained with the widely used HYDRUS 2D/3D and with field tests conducted under different soil conditions. Finally, an illustrative example of its use is shown. The flexibility the code offers the user for the precise definition of the soil type, the possibilities provided to estimate the hydraulic functions, the ponded area considered under the emitter during irrigation, and the output format of the results are the main features, contributions and advantages of this newly developed code.
2. Material and methods 2.1. Theory Soil water infiltration and redistribution under a drip source can be studied using Richards’ equation. Assuming symmetry around the vertical axis passing through the emitter, the equation can be expressed as:
@h @ @h @ @H þ K K ¼ @t @x @x @z @z
ð1Þ
where h is the volumetric soil water content, t is time, x and z are the respective horizontal and vertical co-ordinates, H is the hydraulic head, h represents the pressure head and K the hydraulic conductivity function, depends on h. The Richards’ equation is applied to a cylinder of soil under the emitter using a similar scheme to that described in the cylindrical model of Brandt et al. (1971). The numerical procedure to solve Eq. (1) was based on the scheme applied previously by Van der Ploeg and Benecke (1974), Armstrong and Wilson (1983) and Ramírez de Cartagena (1994). The numerical procedure applied is characterised by its mathematical simplicity, which enabled the presence of different soil layers and the ponded region on the soil surface to be taken into account. It has also been demonstrated to be numerically stable. 2.2. Flow domain and boundary conditions The origin of the coordinates (x = 0 and z = 0) is placed at the centre of the source, the negative z coordinate being downward (Fig. 1). It is assumed that there is no flow beyond the region of interest, which is fixed at a radial distance of 300 cm and a vertical distance of 300 cm, assuming that the influence of the drip source at this distance is negligible. A ‘‘no flow’’ boundary condition at 300 cm depth represents an impediment to the water to infiltrate deeper into the soil profile, but the effect of this impermeable layer will be inconsequential in the region of interest near the emitter. The soil cylinder is divided into concentric rings, 5 cm wide (Di) and 5 cm deep (Dj), centred on the vertical axis passing through the emitter. As a result, this provides an equally spaced mesh of 300 300 cm divided into cells of 5 5 cm each (Fig. 1). The boundary condition at the soil surface was established in order to take into account the variable extension in time of the ponded region formed under the emitter. According to experimental observations, an approximately circular saturated ponded area develops in the vicinity of the emitter. At the start of irrigation, this area is very small and grows over time until it reaches a maximum size (Brandt et al., 1971; Arbat, 2006). In the developed model, this area is the zone where the water enters into the soil, which is assumed to be perfectly circular, and a function of time. Initially, when irrigation starts, the size of the boundary condition is the same as the size of a cell (5 cm), and it will grow when more cells reach saturation. 2.3. Numerical procedure The differential Eq. (1) is transformed into a set of algebraic functions. The variables for a certain position and time take values associated with their position within the flow domain and boundary conditions described in previous section. The procedure applied to solve Eq. (1) is described below: The flow rate, Q (L3 T1) between two consecutive cells (Fig. 1) is calculated as follows:
Q ¼ qA
ð2Þ
185
G. Arbat et al. / Computers and Electronics in Agriculture 98 (2013) 183–192
X = j max
ρ (t)
Radius (j)
Δj
Boundaryconditions
Δi
No flow Evaporation Water t iinfiltration
1,jj ii-1 i,j-1 j
Z = i max
i,j
i,j+1
i+1,j Flow in a generic cell
Depth (i) Fig. 1. Chart of the flow domain, mesh, boundary conditions and water transfers in a Drip-Irriwater generic cell.
where q is the water flow density and A the surface area of the water, flowing from one cell to a consecutive one. The flow rate (Q) is calculated using Buckingham–Darcy’s law:
The change in volumetric soil water content in the ring corresponding to the generic cell (i, j) for a time step (Dt) is:
DH Q ¼ KðhÞ A DL
Dhi;j ¼ ð3Þ
where K(h) is the hydraulic conductivity corresponding to a water content value, DH the difference of hydraulic potential between two consecutive cells (H is defined as negative when the soil is unsaturated) and DL the distance between two cells. As described in the flow domain and boundary conditions section, the model assumes symmetrical flow with respect to the vertical axis passing through the emitter. Therefore, each cell can be associated to a ring with a depth equal to Di and a width equal to Dj. Considering the ring associated with a generic cell within the flow domain (i, j) and its four neighbouring cells (i, j 1), (i 1, j), (i, j + 1),(i + 1, j) (Fig. 1) and applying Eq. (3), the flow rates throughout the contours of the ring associated with the cell (i, j) are computed as: i;j Q 1 ¼ Q i;j i1;j ¼ K i1;j
i Hi;j Hi1;j h 2 p Rjþ1 R2j Di
ð4Þ
i;j Q 2 ¼ Q i;j i;jþ1 ¼ K i;jþ1
Hi;j Hi;jþ1 ð2pRjþ1 DiÞ Dj
ð5Þ
i Hi;j Hiþ1;j h 2 p Rjþ1 R2j Di
ð6Þ
Hi;j Hi;j1 ð2pRj DiÞ Dj
ð7Þ
Q3 ¼
Q i;j iþ1;j
¼
K i;j i1;j
i;j Q 4 ¼ Q i;j i;j1 ¼ K i;j1
ð9Þ
where Vi,j is the volume of the ring associated with cell i, j. Therefore, the soil water change (Dhi,j) is calculated as:
Dhi;j ¼
DQ i;j ðR2jþ1
p
R2j ÞDi
Dt
ð10Þ
The volumetric soil water content at ring i, j for the next time step (t + Dt) is:
htþDt ¼ ht þ Dhi;j
ð11Þ
All of these calculations are made for all of the cells in the flow domain, starting with the cell (1,1), closest to the emitter, and ending with the cell (imax, jmax). The water contents are updated for every cell. The same procedure must be repeated for every time step Dt until the final time specified in the simulation is reached. The boundary conditions of a certain cell depend on its position. Therefore, the flow domain was divided into nine different regions (Fig. 2) in order to calculate the net flow rate in the cells of each region:
For region 1 : i ¼ 1; 1 6 j 6 NðtÞ
where Q1 and Q4 are the vertical water flows through the ring associated with the cell i, j and Q2 and Q3 the horizontal flows through the same ring. K super-script is the mean hydraulic conductivity of the sub-script cells indicated in the subscript and superscript used for the water flow computations. Hsubscript is the hydraulic potential corresponding to the cell indicated in the subscript. Rj and Rj+1 are the interior and exterior radii of the ring associated with the generic cell (i, j). When the sign of the flow rate (Q) computed with Eqs. (4)–(7) is negative, the water flow leaves generic cell i, j. When it is positive, the water flow enters cell i, j. The net flow rate (DQ) in the generic cell (i, j) is given as:
DQ ¼ Q 1 þ Q 2 þ Q 3 þ Q 4
DQ i;j Dt V i;j
ð8Þ
X = jmax
ρ (t) 4
7
2
5
8
3
6
9
1.a
1b
Radius (j)
Z = i max Depth (i) Fig. 2. Regions of the flow domain to define the particular boundary conditions in each group of cells.
186
G. Arbat et al. / Computers and Electronics in Agriculture 98 (2013) 183–192
where N(t) is the number of saturated cells. At the beginning of the simulation, this value is equal to 1. The net flow rate in region 1 is calculated as:
DQ i;j ¼
Qe i;j þ Q i;j i;jþ1 þ Q iþ1;j NðtÞ
ð12Þ
For region 2: 1 < i < imax, j = 1 i;j i;j DQ i;j ¼ Q i;j i1;j þ Q i;jþ1 þ Q iþ1;j
ð13Þ
For region 3: i = imax, j = 1 i;j DQ i;j ¼ Q i;j i1;j þ Q i;jþ1
ð14Þ
For region 4: i = 1, 1 < j < jmax i;j i;j DQ i;j ¼ Q i;j i;jþ1 þ Q iþ1;j þ Q i;j1
ð15Þ
For region 5: 1 < i < imax, 1 < j < jmax i;j i;j i;j DQ i;j ¼ Q i;j i1;j þ Q i;jþ1 þ Q iþ1;j þ Q i;j1
ð16Þ
For region 6: 1 = imax, 1 < j < jmax i;j i;j DQ i;j ¼ Q i;j i1;j þ Q i;jþ1 þ Q i;j1
ð17Þ
For region 7: i = 1, j = jmax i;j DQ i;j ¼ Q i;j iþ1;j þ Q i;j1
ð18Þ
For region 8: 1 < i < imax, j = jmax i;j i;j DQ i;j ¼ Q i;j i1;j þ Q iþ1;j þ Q i;j1
ð19Þ
For region 9: i = imax, j = jmax i;j DQ i;j ¼ Q i;j i1;j þ Q i;j1
ð20Þ
2.4. The Drip-Irriwater code The Drip-Irriwater code is composed of two separate software programs: the graphical user interface (GUI) and the program to perform the numerical procedure described in the numerical procedure section.
The graphical user interface (GUI) was programmed in C#, using the Microsoft.NET framework. It was designed to introduce the input variables and to display the results (soil water and hydraulic potential distribution). The program to perform the numerical procedure was programmed in FORTRAN, using Compaq Digital Visual Fortran 6.0. The code was registered in 2010 with the number 02/2010/ 3036 on the General Registry of Intellectual Property in Spain. Drip-Irriwater software can be obtained from the authors upon request. 2.5. Input and output parameters In the GUI, the input parameters are classified in three groups (Fig. 3): (a) Soil input parameters: the number of soil layers, the model to account for the soil hydraulic functions and the initial water content in each soil layer. A maximum of 8 different soil layers can be defined by the user. The depth of each layer must be a multiple of 5 (the depth of a cell defined in the code) and can have different associated hydraulic functions and initial soil water contents. There are two different options to account for the soil hydraulic functions (h(h) and K(h)) in the Drip-Irriwater code. These are the Saxton model (Saxton et al., 1986) and the van Genuchten– Mualem model (van Genuchten, 1980). When the Saxton model is selected to define the soil hydraulic functions, the user must define the percentage of sand, clay and loam for each particular layer. On the other hand, when the van Genuchten–Mualem model is selected, the user must select the textural class and then the soil hydraulic parameters will be assigned according to the mean values presented by Carsel and Parrish (1988). The user will then have the option to modify some or all of the proposed values. (b) Irrigation parameters: the emitter discharge and irrigation time. (c) Simulation parameters: the simulation time and the time step in the numerical solution (1 s advised).
Fig. 3. Example of the data input interface in Drip-Irriwater.
187
G. Arbat et al. / Computers and Electronics in Agriculture 98 (2013) 183–192
Once the input parameters are specified using the GUI, they must be exported to be written in files. After that, the program to perform the numerical procedure must be executed from the GUI and, when the calculations are finished, the results (distribution of soil water content and hydraulic heads) can be imported and displayed using the GUI (Fig. 4). These output results are displayed numerically and graphically in a window. The cells where these variables have increased after irrigation are automatically coloured in different intensities of blue, according to the value reached by these variables, in order to help identify the soil volume wetted by the emitter more easily (Fig. 4). When multi-layered soil is simulated, high water content has been observed in the cells at the interface between the two horizons, due to the abrupt change in soil properties. In this case and in this position, the cell coloration is less useful in detecting the wetted volume.
characteristics of the different horizons (texture, apparent density and water content of the soil at 33 and 1500 kPa) for each soil profile (Table 1). The bulk density was measured from three undisturbed soil samples for each horizon of each profile and the soil hydraulic characterisation was performed with an equal number of samples and measured in the laboratory with a Richard’s pressure plate apparatus. The soils were classified as Typic Calcixerepts, Typic Haploxeralfs and Oxyaquic Xerofluvents (SSS, 2010). The hydraulic functions of the soils follow the van Genuchten– Mualem model (VG–M). The parameters of these equations were determined from pedotransfer functions using the ROSETTA Lite (Schaap et al., 2001) computer program. The input parameters were the percentages of clay, silt and sand, the apparent density, and the water content at 33 and 1500 kPa (Table 1). The parameters of the VG–M equations are presented in Table 2. 3.2. Irrigation trials
3. Field studies 3.1. Soil characterisation and hydraulic properties In order to verify the model in the field, a series of tests were carried out during the summers of 2001 and 2002 on three representative agricultural soil types. Two of them were located in Monells and the other in Canet de la Tallada (all located in Girona, Spain). A soil trench (3 m long, 0.60 m wide, 2 m deep) was opened at each of the different field sites to allow a detailed pedological description of the soil profile. Determinations were made of the
In each of the three soils, the irrigation tests were carried out with a discharge of 2, 4, 8.5 and 25 L h1, respectively, using an isolated emitter, in each case applying a total volume of 25 L. In order to identify the water input zones below the emitters, the sizes of the pools were measured throughout the testing period. Each test determined the water content before, during and after irrigation with a portable TDR probe (IMKO TRIME-T) and 1.6 m long access tubes placed at horizontal distances of 5, 20, 35, 50, and 75 cm in two perpendicular directions away from the emitter and at 100 and 150 cm in only one of the directions. The position of the tubes
Fig. 4. Example of output results (volumetric soil water content) obtained with Drip-Irriwater.
Table 1 Physical characteristics of different horizons of soils used in the experiments. Soil type and field site
Horizon
Depth (cm)
Clay (%)
Silt (%)
Sand (%)
Bulk density (g cm3)
Volumetric water content at 33 kPa (%)
Volumetric water content at 1500 kPa (%)
Typic Calcixerepts, Monells
Ap Bkn C
0–35 35–80 >80
20.3 20.6 17.0
45.4 45.8 66.5
34.3 33.6 16.5
1.38 1.40 1.42
36.57 23.84 41.60
14.31 13.41 12.80
Typic Haploxeralfs, Monells
Ap Bw1 Bt2
0–20 20–50 >50
23.1 25.3 34.8
28.4 28.3 27.6
48.5 46.4 37.6
1.65 1.59 1.63
29.70 25.44 35.86
14.85 12.72 19.56
Oxyaquic Xerofluvents, Canet de la Tallada
Ap
0–35
11.54
38.87
49.59
1.65
26.83
11.73
Bw1 Bw2 Bw3
35–75 75–120 >120
12.64 4.94 7.01
48.84 26.00 17.34
38.52 69.06 75.65
1.53 1.44 1.68
28.21 30.79 –
12.07 13.95 –
188
G. Arbat et al. / Computers and Electronics in Agriculture 98 (2013) 183–192
Table 2 Parameters of the van Genuchten–Mualem equations for the three different soil types and corresponding horizons. Soil type and field site
Horizon
Depth (cm)
hr (cm3 cm3)
hs (cm3 cm3)
a (cm1)
n
Ks (cm h1)
Typic Calcixerepts, Monells
Ap Bkn C
0–35 35–80 >80
0.0632 0.0553 0.0847
0.4311 0.4035 0.4550
0.0023 0.0273 0.0010
1.7285 1.3377 2.3304
0.52 0.97 0.41
Typic Haploxeralfs, Monells
Ap and Bw1 Bt2
0–50 >50
0.0467 0.0569
0.3691 0.3968
0.098 0.0039
1.3476 1.3779
0.28 0.07
Oxyaquic Xerofluvents, Canet de la Tallada
Ap and Bw1 Bw2 Bw3
0–75 75–120 >120
0.0399 0.0424 0.0384
0.3593 0.4013 0.3354
0.0071 0.0071 0.0460
1.4200 1.4188 1.5295
0.47 1.12 1.73
is shown in Fig. 5. Measurements were taken at every 15 cm from 20 to 140 cm soil depth. The probe readings were corrected using a specific calibration for each of the horizons.
4. Model validation 4.1. Comparison of the results obtained with Drip-Irriwater and HYDRUS HYDRUS (Šimu˚nek et al., 2006), a well-known general purpose commercial software for simulating soil water movement and solute transport in variable-saturated media, was used to assess the soil water distribution under an isolated emitter with the Drip-Irriwater code. In order to compare the results obtained with both codes, the soil water content was simulated on a soil of sandy–loam texture, which was irrigated with an emitter discharge of 4 L h1 over 3 h and 17 min. The water was then left to be redistributed for 12 h from the start of irrigation. The initial soil water content was considered to be 0.184 cm3 cm3 throughout the entire flow domain. The hydraulic properties of the soil were estimated using van Genuchten–Mualem functions, and the parameters were selected from a textural class database incorporated in HYDRUS and DripIrriwater (Carsel and Parrish, 1988), corresponding to a sandy– loam soil. The soil water infiltration was considered to be axisymmetric, and the size of the computational flow domain was 100 cm width in the horizontal and 150 cm deep. This was considered to be large enough to prevent interferences in the water flow from the emitter and not so large as to increase the simulation time unnecessarily. The rectangular flow region (100 cm wide 150 cm depth) was discretised into a structured triangular finite element mesh. The number of nodes on the horizontal and vertical sides of the rectangular region was 50 and 70 respectively. The element sizes were then distributed in the rectangular domain resulting in
triangular finite elements of approximately 2 cm. Nevertheless the size of the elements at the top of the domain was defined to be slightly smaller than the ones at the bottom because more active flow was expected. Initially the size of the elements in HYDRUS was selected to be 5 cm, to match the mesh size used by the finite differences method used in Drip-Irriwater, but these large elements produced problems during the numerical solution. Therefore, a smaller size of the elements was defined. No-flow boundary conditions were imposed on the right, left and bottom sides of the computational flow domain in the same way as the boundary conditions used by Drip-Irriwater. These boundary conditions have been reported in previous studies to simulate water infiltration from an isolated emitter using HYDRUS (Kandelous and Šimu˚nek, 2010). As the standard version of HYDRUS-2D (Šimu˚nek et al., 2006) does not have a dynamic boundary condition implemented that is suitable for simulating variable ponding at the soil surface beneath the drip emitter during irrigation, the water entry at the soil surface was simulated as a constant water flux within an area of 30 cm from the upper right corner of the flow domain, according to the observed saturated radius in similar textural soil characteristics with an emitter discharge rate of 4 L h1. The soil water contents computed with Drip-Irriwater at 5, 15, 25, 35, 45 and 55 cm horizontal distances from the emitter and at 5, 15, 25, 35, 45, 55 and 60 cm depths were visually compared with those at the same positions obtained with HYDRUS and represented in a two dimensional plot of iso-water contents interpolating the values between these points using the Kriging method. The comparison was also performed between the soil water content vertically at 30 cm distance from the dripper and laterally at 30 cm depth from the soil surface. In order to determine the agreement of both models the root mean square error (RMSE) and the coefficient of determination (r2) were computed in the same grid of points as follows:
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Pn 0 2 i¼1 ðP i P i Þ RMSE ¼ n
75 50
ð21Þ
35 20
12 0 0 P P P P i C B i i¼1 C rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r2 ¼ B @rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 A Pn 0 0 2 Pn i¼1 P i P i¼1 P i P 0
Pn
ð22Þ
Emitter 20
35
50
100
75
150
Y
5
X
Fig. 5. Chart showing the position of the TDR access tubes for the measurement of the soil water content in the field tests.
where Pi is soil water content predictedby Drip-Irriwater, P0i the soil water content obtained with HYDRUS, P the average soil water con tent predicted with Drip-Irriwater, P 0 the average soil water content obtained with HYDRUS, and n the number of observations. The units of the RMSE are cm3 cm3.
189
G. Arbat et al. / Computers and Electronics in Agriculture 98 (2013) 183–192
4.2. Comparison between simulation and experimental results
drip irrigation experiments (Skaggs et al., 2004; Arbat et al., 2008; Kandelous and Šimu˚nek, 2010).
Simulated and measured soil water contents, and the maximum radius (R) and depth (D) of the wetted soil volume were compared at the end of irrigation and at 24 h after the end of irrigation. R and D were computed using the criteria that soil water content was 5% higher than initially measured. As the computational mesh size was 5 5 cm, R and D were always multiples of 5. Simulated and measured soil water contents were compared at all the measured points by calculating the RMSE and r2 in each of the tests at the end of irrigation and after 24 h. The calculations were made with Eqs. (21) and (22), in this case O indicates the measured water contents, and P the water contents with predicted Drip-Irriwater. The measured soil water content from 0 to 75 cm was the average soil water content measured in two perpendicular directions (Fig. 5). r2 indicates the degree of error variance between two values. Higher r2 value means less error variance between the observed and predicted values. The RMSE has the advantage of expressing the error with the same units as the variable. Both statistics were selected because they had previously been used to compare predicted and measured soil water contents in different
5. Results and discussion 5.1. Comparison between Drip-Irriwater estimations and HYDRUS The comparison of soil water distribution obtained from DripIrriwater and HYDRUS showed a similar moisture distribution pattern (Fig. 6). When comparing the simulated water contents within the soil wetted region (55 cm depth and 55 cm from the emitter source) obtained with both models, the RMSE was 0.018 cm3 cm3 and r2 was 0.86. Fig. 7 shows a comparison between the soil water content vertically at selected positions of 30 cm distance from the emitter and laterally at 30 cm depth from the soil surface. The RMSEs comparing the water contents were 0.007 cm3 cm3 for the vertical line at 30 cm from the emitter, and 0.013 cm3 cm3 for the horizontal line 30 cm deep, while r2 for the same positions was 0.931 and 0.934, respectively. The statistics comparing the simulated water contents within the soil wetted region and the ones obtained in selected positions indicate small differences in
Fig. 6. Soil water content distribution at 12 h after starting irrigation with a drip flow of 4 L h1 and irrigation time of 3 h 17 min, with HYDRUS (a) and Drip-Irriwater (b).
Vertical sección at 30cm from the emitter 0.00 0
Depth (cm) D
10 20
0.05
0.10
HYDRUS Drip-Irriwater
30 40 50 60 70
RMSE = 0.007 (cm3/cm3) R2 = 0.931
0.15
0.20
0.25
0.30
Vo olume etric s soil wa ater c conten nt (cm m 3 cm m -3)
Volumetric soil water content (cm 3 cm -3)
Horitzontal section at 30 cm depth
0.35
HYDRUS
0.30
Drip-Irriwater
0.25 0. 20 0 15 0.
RMSE = 0.013 0 013 (cm3/cm3) R 2 = 0.934
0.10 0.05 0.00 0
5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90
Horizontal distance (cm) Fig. 7. Comparison of soil water content between HYDRUS (red dots) and Drip-Irriwater (blue dots): (i) vertically at 30 cm away from the emitter and (ii) horizontally at 30 cm depth from the soil surface corresponding a drip flow of 4 L h1 for an irrigation time of 3 h and 17 min 12 h after stopping irrigation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
190
G. Arbat et al. / Computers and Electronics in Agriculture 98 (2013) 183–192
the soil content at certain points, which resulted in part from the different dimensions of the computational mesh used by the two codes. While in Drip-Irriwater the mesh was composed of square cells of 5 5 cm, in HYDRUS the triangular finite elements had vertical and horizontal sides ranging from 1.5 to 2 cm, allowing greater resolution in the resulting soil water distribution. The different boundary conditions, used to simulate the soil volume where the water infiltrated under the emitter, also explains the slightly different results obtained in the two computer programs. With an accuracy of 5 cm, equal to the length of the cells used in Drip-Irriwater, the radius of the wetted soil volume at a depth of 30 cm was approximately 40–45 cm and the maximum wetted depth was 55–60 cm, irrespective of the software used (Fig. 6). 5.2. Comparison of the field tests and Drip-Irriwater results The comparison between the water contents simulated with the Drip-Irriwater code and the ones measured at the end of irrigation (Table 3) shows that the average RMSE was 0.028 cm3 cm3 (with a minimum value of 0.019 and a maximum of 0.035) and an average r2 of 0.54 (with a minimum value of 0.27 and a maximum of 0.83). Twenty-four hours after completion of irrigation (Table 4) the average RMSE was 0.026 cm3 cm3 (with a minimum value of 0.015 and maximum of 0.040) and the average r2 was 0.56 (with a minimum value of 0.45 and maximum of 0.75). The differences between the water contents in the simulations and in the tests could be partly due to the errors inherent in the method of measurement. In the calibrations performed on the same sites where tests were performed, comparing soil water contents measured with TDR and soil samples, the RMSE was between 0.030 cm3 cm3
and 0.050 cm3 cm3, and the coefficient of determination (r2) between 0.57 and 0.89 (Arbat, 2006). This shows that sometimes the measurement method could affect the determination of the wetted radius and depth, which explains the significant differences found between measured and simulated values in some comparisons (Table 3). The Drip-Irriwater code simulates a slightly greater wetted radius than the wetted radius determined from TDR measurements, with a mean difference of 7.3 cm at the end of irrigation and 4.5 cm after 24 h of redistribution (Tables 3 and 4). Similarly the wetted depth determined with the Drip-Irriwater code was smaller than the wetted depth determined from measurements with TDR, with a mean difference of 10 cm at the end of irrigation and 2.7 cm after 24 h of redistribution. The differences between the values of the wetted radius and depth in the simulated and those determined from measurements with TDR were lower 24 h after the end of irrigation. It is considered that the differences between the values of the wetted radius and depth were smaller overall and the Drip-Irriwater code allowed us to estimate these parameters appropriately for the agronomic design of a drip irrigation system. Kandelous and Šimu˚nek (2010) compared the extent of wetted soil volume, using analytical and empirical models, with the one obtained with the numerical model HYDRUS, concluding that the predictions of the latter are more realistic. However, in surface drip irrigation, HYDRUS may underestimate the wetted horizontal length and overestimate the wetted depth when pounded area is formed on the soil surface. This is because the static boundary condition defined in HYDRUS does not reflect the extension of the pounded area. In this regard, the
Table 3 Statistics of comparison between measured and simulated soil water contents, and values of measured and simulated wetted radii and depth with an irrigation volume of 25 L at the end of the irrigation event. Soil type
Flow rate (L h1)
RMSE (cm3 cm3)
r2
Measured radii (cm)
Simulated radii (cm)
Measured depth (cm)
Simulated depth (cm)
Typic Calcixerepts
2.0 4.0 8.5 25.0
0.035 0.021 0.028 0.033
0.36 0.83 0.67 0.49
25 40 40 40
45 45 40 40
50 40 45 40
35 35 35 30
Typic Haploxeralfs
2.0 4.0 8.5 25.0
0.034 0.034 0.028 0.019
0.67 0.59 0.61 0.43
45 30 30 20
40 40 40 40
40 30 25 30
35 25 20 15
Oxyaquic Xerofluvents
4.0 8.5 25.0
0.029 0.023 0.027
0.46 0.58 0.27
35 40 35
45 45 40
45 40 40
35 30 20
0.028
0.54
Average 2
RMSE: Root mean square error, r : coefficient of determination. Table 4 Statistics of comparison between the measured and simulated soil water contents, and values of measured and simulated wetted radii and depth with an applied irrigation volume of 25 L, 24 h after the end of irrigation. Soil type
Flow rate (L h1)
RMSE (cm3 cm3)
r2
Measured radii (cm)
Simulated radii (cm)
Measured depth (cm)
Simulated depth (cm)
Typic Calcixerepts
2.0 4.0 8.5 25.0
0.026 0.015 0.018 0.018
0.47 0.63 0.50 0.51
35 45 40 40
50 50 45 50
50 40 45 35
35 35 35 35
Typic Haploxeralfs
2.0 4.0 8.5 25.0
0.027 0.035 0.040 0.034
0.74 0.67 0.59 0.67
45 45 40 40
45 40 45 45
45 45 55 50
50 40 45 45
Oxyaquic Xerofluvents
4.0 8.5 25.0
0.026 0.024 0.018
0.43 0.45 0.52
60 40 40
50 55 45
50 45 40
60 50 40
0.026
0.56
Average
RMSE: Root mean square error, r2: coefficient of determination.
G. Arbat et al. / Computers and Electronics in Agriculture 98 (2013) 183–192
191
Fig. 8. Input and output GUI of Drip-Irriwater for 2 soils (loamy sand and sandy) irrigated with a drip discharge of 4 L h1 for 2 h.
Drip-Irriwater code adds the effect of the pond as a dynamic boundary condition to the Richards equation. As described in the Numerical Procedure section, the horizontal extension of the ponded area is recalculated for each time step, increasing the extension of the infiltration zone when a new cell is saturated (region 1 in Fig. 2). Therefore the Drip-Irriwater code handles a dynamic boundary condition for the infiltration zone under the emitter, while in the standard version of HYDRUS this condition is fixed by the user and kept constant throughout the simulation. The use of this dynamic boundary condition in the Drip-Irriwater code is more realistic when the application rate is greater than the infiltration capacity of the soil. 5.3. Illustrative examples Two examples are given in this section in order to guide the user. Similar amount of irrigation water (8 L) was applied in both cases with a 4 L h1 flow rate emitter and the water content were only displayed after 2 h of irrigation in two different soil textures. In the first case, the soil was of loamy–sand texture and in the second it was of sandy texture. The soil hydraulic parameters inferred from the texture were used in the simulations. The initial soil water content was set up as 15.5% in both cases. Fig. 8 shows the input data for the loamy–sand texture, and the output soil water contents for both soils. It is worth noting that the soil texture plays an important role in the shape and size of the wetting patterns. In the sandy soil, the maximum wetted depth after 2 h was approximately 65 cm, clearly deeper than that of 50 cm, observed in the loamy–sand soil. 6. Conclusions The Drip-Irriwater code allows the simulation of soil water distribution under surface drip irrigation, considering the main
parameters involved in the shape and dimensions of the volume of wet soil such as: soil type, emitter flow and irrigation time, the initial soil water content and redistribution time after irrigation. The soil water content distribution obtained with Drip-Irriwater was compared with those obtained with the widely recognised HYDRUS code and the comparison statistics (RMSE, r2) indicate good agreement between both. The small differences observed may be partially due to the fact that Drip-Irriwater incorporates a dynamic boundary condition, which takes the ponded surface into account. This condition is not available in the standard version of HYDRUS. In particular, the design variables obtained with both codes (wetted radius and depth) were identical for the simulated case. The field verification of the model revealed an acceptable agreement between measured and simulated soil water content for surface drip irrigation design. Overall, the differences between the wetted radius and depth were low and could be partly explained by the errors inherent to the measurement method. It can be concluded that the Drip-Irriwater software is accessible to design engineers, irrigation technicians and practitioners, as the input parameters can be easily known or estimated. The model could also be used to estimate the soil planning factors (wetted radius and wetted depth) with sufficient approximation to efficiently design surface drip irrigation systems. Acknowledgments The authors would like to express their gratitude to the Institut de Recerca i Tecnologia Agroalimentàries (IRTA), specially to the centres of IRTA-Mas Badia in Canet de la Tallada, and IRTA-Finca Camps i Armet in Monells for their help in carrying out this investigation.
192
G. Arbat et al. / Computers and Electronics in Agriculture 98 (2013) 183–192
References Amin, M., Ekhmaj, A., 2006. DIPAC-drip irrigation water distribution calculator. In: 7th International Micro Irrigation Congress, 10–16 September, PWTC, Kuala Lumpur, Malaysia. Annandale, J.G., Jovanovic, N.Z., Campbell, G.S., Du Sautoy, N., Benadé, N., 2003. A two-dimensional water balance model for micro-irrigated hedgerow tree crops. Irrig. Sci. 22 (3–4), 157–170. Arbat, G., 2006. Desarrollo y validación de la dinámica del agua en el suelo. Aplicación al diseño agronómico y al manejo en riego localizado. Tesis Doctoral, Universitat de Lleida, Lleida, Spain. Arbat, G., Puig-Bargués, J., Barragán, J., Bonany, J., Ramírez de Cartagena, F., 2008. Monitoring soil water status for micro-irrigation management versus modelling approach. Biosyst. Eng. 100 (2), 286–296. Arbat, G., Roselló, A., Domingo Olivé, F., Puig-Bargués, J., González Llinàs, E., DuranRos, M., Pujol, J., Ramírez de Cartagena, F., 2013. Soil water and nitrate distribution under drip irrigated corn receiving pig slurry. Agric. Water Manage. 120, 11–22. Armstrong, C.F., Wilson, T.V., 1983. Computer model for moisture distribution in stratified soils under a trickle irrigation source. Trans. ASABE 26 (6), 1704–1709. Bastiaanssen, W.G.M., Allen, R.G., Droogers, P., D’urso, G., Steduto, P., 2007. Twentyfive years modeling irrigated and drained soils: State of the Art. Agric. Water Manage. 92 (3), 111–125. Battam, M.A., Boughton, D.G., Sutton, B.G., 2003. Soil pits as a simple design aid for subsurface drip irrigation systems. Irrig. Sci. 22 (3), 135–142. Brandt, A., Bresler, E., Diner, N., Ben Asher, J., Heller, J., Goldberg, D., 1971. Soil Sci. Soc. Am. Proc. 35 (5), 675–682. Carsel, R., Parrish, R., 1988. Developing joint probability of soil water retention characteristics. Water Resour. Res. 24 (5), 755–769. Cook, F.J., Thorburn, P.J., Fitch, P., Bristow, K.L., 2003. WetUp: a software tool to display approximate wetting patterns from drippers. Irrig. Sci. 22 (3–4), 129– 134. Cook, F.J., Thorburn, P.J., Fitch, P., Charlesworth, P.B., Bristow, K.L., 2006. Modelling trickle irrigation: comparison of analytical and numerical models for estimation of wetting front position with time. Environ. Model Soft. 21 (9), 1353–1359. Cote, C.M., Bristow, K.L., Charlesworth, P.B., Cook, F.J., Thornburn, P.J., 2003. Analysis of soil wetting and solute transport in subsurface trickle irrigation. Irrig. Sci. 22 (3–4), 143–156. Ghali, G.S., 1986. Mathematical Modelling of Soil Moisture Dynamics in Trickle Irrigated Fields. Ph.D. Thesis, University of Southampton. Hachum, A.Y., Alfaro, J.F., Willardson, L.S., 1976. Water movement in soil from a trickle source. J. Irrig. Drain. Div. Proc. ASCE 102 (IR2), 179–192. Hinnell, A.C., Lazarovitch, N., Furman, A., Poulton, M., Warrick, A.W., 2010. NeuroDrip: estimation of subsurface wetting patterns for drip irrigation using neural networks. Irrig. Sci. 28 (6), 535–544. Kandelous, M., Liagahat, A., Abbasi, F., 2008. Estimation of soil moisture pattern in subsurface drip irrigation using dimensional analysis method. J. Agric. Sci. 39 (2), 371–378. Kandelous, M., Šimu˚nek, J., 2010. Comparison of numerical, analytical, and empirical models to estimate wetting patterns for surface and subsurface drip irrigation. Irrig. Sci. 28 (5), 435–444. Keller, J., Bliesner, R.D., 1990. Trickle irrigation planning factors. In: Van Nostrand Reinhold (Ed.), Sprinkle and Trickle Irrigation. Avi Book, New York, pp. 453–477.
Lafolie, F., Guennelon, R., van Genuchten, M.T., 1989. Analysis of water flow under trickle irrigation: I. Theory and numerical solution. Soil Sci. Soc. Am. J. 53 (5), 1310–1318. MAGRAMA (Ministerio de Agricultura, Alimentación y Medio Ambiente de España), 2011. Encuesta sobre superficies y rendimientos de cultivos. Anàlisis de los regadíos españoles. Secretaría General Técnica, Madrid. Malek, K., Peters, R.T., 2011. Wetting pattern models for drip irrigation: new empirical model. J. Irrig. Drain. Eng. 137 (8), 530–536. Phene, C.J., 1995. Research trends in microirrigation. Microirrigation for a changing world: conserving resources/preserving the environment. In: Proceedings of the Fifth Microirrigation Congress, 2–6 April, Orlando, Florida. Philip, J.R., 1984. Travel times from buried and surface infiltration point sources. Water Resour. Res. 20 (7), 990–994. Raats, P.A.C., 1971. Steady infiltration from point sources, cavities and basins. Soil Sci. Soc. Am. Proc. 35, 689–694. Ramírez de Cartagena, F., 1994. Simulación numérica de la dinámica del agua en el suelo. Aplicación al diseño de sistemas de riego LAF. Tesis Doctoral, Universitat de Lleida, Lleida, Spain. Revol, P., Clothier, B.E., Mailhol, J.C., Vachaud, G., Vauclin, M., 1997. Infiltration from a surface point source and drip irrigation. 2. An approximate time-dependent solution for a wet-front position. Water Resour. Res. 33 (8), 1869–1874. Richards, L.A., 1931. Capillary conduction of liquids through porous mediums. Physics 1, 318–333. SSS (Soil Survey Staff), 2010. Keys to Soil Taxonomy, 11th ed. USDA, Natural Resources Conservation Service. Saxton, K.E., Rawls, W.J., Romberger, J.S., Papendick, R.I., 1986. Estimating generalized soil–water characteristics from texture. Soil Sci. Soc. Am. J. 50 (4), 1031–1036. Schaap, M.G., Leij, F.J., van Genuchten, M.T., 2001. ROSETTA: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol. 251 (3–4), 163–176. Schwartzmass, M., Zur, B., 1986. Emitter spacing and geometry of wetted soil volume. J. Irrig. Drain. Eng. 112 (3), 242–253. Šimu˚nek, J., van Genuchten, M.T., Sejna, M., 2006. HYDRUS, software package for simulating the two- and three-dimensional movement of water, heat and solute transport in variably-saturated media. Technical Manual v. 1. PC Progress, Prague, Czech Republic. Skaggs, T.H., Trou, T., Šimu˚nek, J., Shouse, P.J., 2004. Comparison of HYDRUS-2D simulations of drip irrigation with experimental observations. J. Irrig. Drain. Eng. 130 (4), 304–310. Taghavi, S.A., Mariño, M.A., Rolston, D.E., 1984. Infiltration from trickle irrigation source. J. Irrig. Drain. Eng. 110 (4), 331–341. Thorburn, P.J., Cook, F.J., Bristow, K.L., 2003. Soil dependent wetting from trickle emitters: implication for design and management. Irrig. Sci. 22 (3–4), 121–127. Van der Ploeg, R.R., Benecke, P., 1974. Unsteady, unsaturated, n-dimensional moisture flow in soil: a computer simulation program. Soil Sci. Soc. Am. Proc. 38 (6), 881–885. van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892–898. Warrick, A.W., 1974. Time-dependent linearized infiltration: I. Point sources. Soil Sci. Soc. Am. Proc. 38 (3), 383–386. Warrick, A.W., 2003. Soil Water Dynamics. Oxford University Press, New York.