Estimating soil wetting patterns for drip irrigation using genetic programming

Drip irrigation is considered as one of the most efficient irrigation systems. Knowledge of the soil wetted perimeter arising from infiltration of water from drippers is important in the design and management of efficient irrigation systems. To this aim, numerical models can represent a powerful tool to analyze the evolution of the wetting pattern during irrigation, in order to explore drip irrigation management strategies, to set up the duration of irrigation, and finally to optimize water use efficiency. This paper examines the potential of genetic programming (GP) in simulating wetting patterns of drip irrigation. First by considering 12 different soil textures of USDA-SCS soil texture triangle, different emitter discharge and duration of irrigation, soil wetting patterns have been simulated by using HYDRUS 2D software. Then using the calculated values of depth and radius of wetting pattern as target outputs, two different GP models have been considered. Finally, the capability of GP for simulating wetting patterns was analyzed using some values of data set that were not used in training. Results showed that the GP method had good agreement with results of HYDRUS 2D software in the case of considering full set of operators with R2 of 0.99 and 0.99 and root mean squared error of 2.88 and 4.94 in estimation of radius and depth of wetting patterns, respectively. Also, field experimental results in a sandy loam soil with emitter discharge of 4 L h–1 showed reasonable agreement with GP results. As a conclusion, the results of the study demonstrate the usefulness of the GP method for estimating wetting patterns of drip irrigation. Additional key words: genetic programming; HYDRUS 2D; infiltration; numerical models; soil texture triangle.

the measurement of soil infiltration parameters, as well as many of the complexities and challenges in applying current understanding to irrigation situations.
Analytical techniques have been proposed for the study of infiltration from a surface point-source (Wooding, 1968;Warrick, 1974;Bresler, 1978;Ben Asher et al., 1986), but these are all limited by one or more simplifying assumptions.Numerical methods also have been developed to simulate this phenomenon (Brandt et al., 1971;Taghavi et al., 1984;Healy, 1987).For instance, HYDRUS 2D is a model based on finite-element numerical solutions of the flow equations (Simunek et al., 2006) allowing simulations of threedimensional axially symmetric water flow, solute transport and root water and nutrient uptake.Modeling studies by Assouline (2002) and Abbasi et al. (2003a,b) showed that HYDRUS 2D simulations of soil water content and solute distributions were reasonably close to measured values.Cote et al. (2003) used the HYDRUS 2D model to simulate soil water and solute transport under subsurface drip irrigation.Skaggs et al. (2004) concluded that HYDRUS 2D had predicted soil water distribution for drip tape irrigation that agreed well with experimental observations.Bufon et al. (2011) experimentally validated HYDRUS 2D simulations in an Amarillo soil for cotton subsurface drip irrigation in Texas high plains.Kandelous & Simunek

Introduction
Drip irrigation has been regarded as a potentially efficient method of irrigation.The potential benefits of using drip irrigation include higher yields (Camp, 1998), improved trafficability (Steele et al.., 1996) and lower water use (Camp, 1998).Despite these potential advantages, poor design and/or poor management can result in water losses from drip irrigation comparable with those from more traditional irrigation systems.Given the high installation costs often required for drip irrigation (Darusman et al., 1997), it is crucial that systems are designed and managed correctly if the benefits of using drip irrigation are to be fully exploited.Many of the design and management decisions require understanding the wetted zone pattern around the emitter (Bresler, 1978;Lubana & Narda, 2001) and its relation to the root system.Mathematical models have proven very useful in predicting water and nutrient movement in the soil (Mmolawa & Or, 2000a) enabling improved drip irrigation system design and management.Or (1995) investigated the effects of mild spatial variation in soil hydraulic properties on wetting patterns and the consequences on soil water sensor placement and interpretation.Mmolawa & Or (2000b) discussed aspects of soil water and solute dynamics as affected by the irrigation and fertigation methods, in the presence of active plant uptake of water and solutes.Vrugt et al. (2001) tested the suitability of a three dimensional root water uptake model for simultaneous simulation of transient soil water flow around an almond tree and compared performance and results of mentioned model with oneand two-dimensional root water uptake models.Smith & Warrick (2007) presented basic relations of soil water and soil water flow which are important in irrigation design and outlined methods to measure soil water content, pressure head and conductivity.They also discussed the calculation of infiltration rates and Abbreviations used: GA (genetic algorithm); GEP (gene expression programming); GP (genetic programming); GP 4basic (genetic programming with four basic operators); GP fullset (genetic programming with full set of operators); MLR (multiple linear regression), RMSE (root mean squared error).Symbols: d m (distance from emitter computed by HYDRUS 2D, cm); d s (distance from emitter computed by different methods, cm); E(ij) (error of an individual program); F (set of functions); f i (fitness); h (pressure head in soil-hydraulic function, m); h (length of the head in genetic programming; K (unsaturated hydraulic conductivity function, cm day -1 ); K ij A (components of a dimensionless anisotropy tensor K A , cm day -1 ); K r (relative hydraulic conductivity, cm day -1 ); K s (saturated hydraulic conductivity, cm day -1 ); l (pore connectivity parameter); n (shape parameter in soil-hydraulic function); n (number of values in evaluation parameters); p (precision); P (ij) (predicted value by the individual program); Q (emitter discharge, L h -1 ); r (radius of wetting patterns, cm); R (selection range); RMSE (root mean squared error); S (sink term, s -1 ); S e (degree of saturation); SS A (between-group some of squares); SS E (within-group or error sum of squares); SS T (total sum of squares); T (set of terminals); t (time, h); T j (target value for fitness case); x i (i(1,2) (spatial coordinates, m); Y j .
(mean of each group); Y .. (grand mean); Y ij (single score); z (depth of wetting patterns, cm); α (shape parameter); θ (volumetric water content); θ i (initial water content); θ r (residual water content); θ s (saturated water content).Drip wetting patterns with genetic programming putationally intensive reducing the application of numerical models for irrigation system design and management decision-making (Hinnell et al., 2010).Moradi-Jalal et al. (2004) presented a new management model, WAPIRRA Scheduler, for the optimal design and operation of water distribution systems.Their model used genetic algorithm (GA) optimization to automatically determine annually the least cost of pumping stations while satisfying target hydraulic performance requirements.Application of their model to a study case showed considerable savings in cost and energy.Kumar et al. (2006) developed a model based on genetic algorithm (GA) for obtaining an optimal operating policy and optimal crop water allocations from an irrigation reservoir.They reported that the model can be used for optimal utilization of the available water resources of any reservoir system to obtain maximum benefits.Reca & Martinez (2006) developed a computer model called Genetic Algorithm Pipe Network Optimization Model (GENOME) with the aim of optimizing the design of new looped irrigation water distribution networks and they optimized a real complex irrigation network to evaluate the potential of the genetic algorithm for the optimal design of large-scale networks.They reported that although the mentioned model produced satisfactory results, some adjustments would be desirable to improve the performance of genetic algorithms when the complexity of the network requires it.Genetic programming has been implemented in different fields of water engineering such as rainfall-runoff modeling (Aytac & Alp, 2008), filling up gaps in wave data (Ustoorikar & Deo, 2008) and evapotranspiration modeling (Kisi & Guven, 2010) but, to the best of our knowledge, the input-output mapping capability of GP technique in estimating soil wetting patterns for drip irrigation has not been studied.The goal of this study is to propose an alternative way which produces precise results in comparison with numerical models for computation of the spatial and temporal wetting patterns during infiltration from surface drip emitters using genetic programming.

Governing flow equation
Consider two-and/or three-dimensional isothermal uniform Darcian flow of water in a variably saturated rigid porous medium and assume that the air phase plays an insignificant role in the liquid flow process.The governing flow equation for these conditions is given by the following modified form of the Richards' equation (Simunek et al., 2006): where ) ( ) [2]   where K r = relative hydraulic conductivity and K s = saturated hydraulic conductivity [LT -1 ].The anisotropy tensor [1] is used to account for an anisotropic medium.The diagonal entries of K ij A equal one and the offdiagonal entries equals zero for an isotropic medium.Also soil-hydraulic functions of van Genuchten (1980) who used the statistical pore-size distribution model of Mualem (1976) to obtain a predictive equation for the unsaturated hydraulic conductivity function in terms of soil water retention parameters are given by where The above equations contain six independent parameters: θ r , θ s , α, n, K s , and l.The pore connectivity parameter l in the hydraulic conductivity function was estimated (Mualem, 1976) to be about 0.5 as an average for many soils.

General overview of genetic programming
In this section, a brief overview of the genetic programming (GP) and gene expression programming (GEP) is given.GP is a generalization of genetic algorithms (GAs) (Goldberg, 1989).Detailed explanations of GP and GEP are provided by Koza (1992) and Ferreira (2006), respectively.The fundamental difference between GA, GP and GEP is due to the nature of the individuals.In the GA, the individuals are linear strings of fixed length (chromosomes).In the GP, the individuals are nonlinear entities of different sizes and shapes (parse trees), and in GEP the individuals are encoded as linear strings of fixed length (the genome or chromosomes), which are afterwards expressed as nonlinear entities of different sizes and shapes (Ferreira, 2001a,b).GP is a search technique that allows the solution of problems by automatically generating algorithms and expressions.These expressions are coded or represented as a tree structure with its terminals (leaves) and nodes (functions).GP applies GAs to a "population" of programs, i.e., typically encoded as tree-structures.Trial programs are evaluated against a "fitness function" and the best solutions selected for modification and reevaluation.This modification-evaluation cycle is repeated until a "correct" program is produced.
There are five preliminary steps for solving a problem by using GP.These are the determination of (i) the set of terminals, (ii) the set of functions, (iii) the fitness measure, (iv) the values of the numerical parameters and qualitative variables for controlling the run, and (v) the criterion for designating a result and terminating a run (Koza, 1992).
A GEP flowchart improved by Ferreira (2001b) is presented in Suppl.Fig. 1

(pdf).
The automatic program generation is carried out by means of a process derived from Darwin's evolution theory, in which, after subsequent generations, new trees (individuals) are produced from old ones via crossover, copy, and mutation (Fuchs, 1998;Luke & Spector, 1998).Based on natural selection, the best trees will have more chances of being chosen to become part of the next generation.Thus, a stochastic process is established where, after successive generations, a well-adapted tree is obtained.
There are five steps in preparing to use GEP of which the first is to choose the fitness function.The fitness of an individual program i for fitness case j is evaluated by Ferreira (2006) using: where p is the precision and E(ij) is the error of an individual program i for fitness case j.For the absolute error, this is expressed by: Again for the absolute error, the fitness f i of an individual program i is expressed by: where R is the selection range, P (ij) is the value predicted by the individual program i for fitness case j (out of n fitness cases) and T j is the target value for fitness case j.The second step consists of choosing the set of terminals T and the set of functions F to create the chromosomes.In this problem, the terminal set obviously consists of the independent variables.The choice of the appropriate function set is not so obvious.However, a good guess can always be helpful in order to include all of the necessary functions.In this study, four basic arithmetic operators, i.e. ( +, -, ×, /) and some basic mathematical functions, i.e. ( +, -, ×, /, Log, Ln, Power, Sin, Cosine, Arctangent) were utilized.The third step is to choose the chromosomal architecture, i.e., the length of the head and the number of genes.Values of the length of the head, h = 8, and three genes per chromosome were employed.The fourth step is to choose the linking function.In this study, the subprograms were linked by addition.Finally, the last step is to choose the set of genetic operators that cause variation and their rates.A combination of all 1159 Drip wetting patterns with genetic programming genetic operators, i.e., mutation, transposition and recombination, was used for this purpose.The parameters of the training of the GP are given in Table 1.

Numerical simulation of wetting pattern for drip irrigation
HYDRUS-2D which uses the Galerkin finiteelement method to solve Eqs.[1] to [5] was applied to simulate the three dimensional axis-symmetric water flow.Simulations were carried out considering 150 cm deep and 120 cm wide soil profile, where an emitter was placed on the soil surface.The computational flow domain was made large enough to ensure that the right and bottom boundaries did not affect the simulations.Absence of flux was considered along the surface and the lateral boundaries and free drainage along the bottom boundary of the soil profile.A constant flux density corresponding to the emitter discharge was assumed along the emitter boundary surface and free drainage was considered in the bottom of the domain.Also initial water content in whole domain was assumed as initial condition.An unstructured mesh was automatically generated to discretize the flow domain into triangles.A total of 7,320 nodes were used to represent the entire simulation domain.Fig. 1 shows the scheme of the grid used for the numerical simulations.

Field experiment
For evaluating the accuracy of the proposed genetic programming method, an experiment of water infiltration under drip irrigation was conducted on a sandy loam soil at the Khalatpooshan Agricultural Sciences Center (37°03' N, 46°37' E and 1567.3 m asl), located in Tabriz, northwest of Iran.A polyethylene drip pipeline with a length of 15 m was installed on the soil surface which has a 16 mm outside diameter, a wall thickness of 2 mm, an emitter spacing of 2 m and emitter discharge of 4 L h -1 .During irrigation, soil was excavated around each emitter and wetted pattern dimensions were measured.A coordinate system was established on the profile with the origin at the soil surface directly above the drip line.The observed soil wetting had a high degree of horizontal symmetry.
The bulk and particle density of the soil was 1.62 and 2.61 g cm -3 , respectively.Saturated hydraulic conductivity, measured by double rings apparatus, was 33 cm day -1 .Soil water retention curve was determined during a drying process by using pressure membrane extractor (model 1000) which was manufactured by Soil Moisture Equipment Corporation.This model of pressure membrane extractor incorporates disposable cellulose membranes in the extraction of water from soil samples over a pressure range of 0 to 15 bars.

Evaluation parameters
Several parameters can be considered for the evaluation of radius and depth of wetting patterns estimates.In this study the following statistic criteria were used: root mean squared error (RMSE) and correlation coefficient (R 2 ).
where d m = distance from emitter computed by HYD-RUS 2D (cm), d s = distance from emitter computed by different methods (cm) and n = number of values.

Simulation of wetting patterns using genetic programming
In this research the HYDRUS 2D software was used to compute depth and radius of 1,500 simulations using 12 different soil textures of USDA-SCS soil texture triangle.Emitter discharge ranged from 2 to 8 L h -1 and irrigation duration varied from 0.5 to 8 h.Then 1,125 sets of these data were used for training of GP.Finally capability of GP for simulating wetting patterns of drip irrigation was analyzed using some dataset values that were not used for training.
One of the advantages of GP in comparison with other tools is producing analytical formula for determination of output parameters.After processing, Eqs.[11] and [12] for determination of radius and depth of wetting patterns with full set of operators ( +, -, ×, /, Log, Ln, Power, Sin, Cosine, Arctangent), respectively, and Eqs.[13] and [14] where r = radius of wetting patterns [cm], z = depth of wetting patterns [cm], θ r = residual water content [-], Due to the fact that calculation of radius and depth of wetting patterns of drip irrigation with these equations is somehow time consuming, a code with programming language of Wolfram Mathematica 7.0 was written.This code, which a picture of that is shown in Fig. 2, can dynamically calculate radius and depth of wetting patterns.In addition, the developed code has ability to animate wetting patterns during irrigation.

Simulation of wetting patterns using multiple linear regression
Furthermore, the multiple linear regression method (MLR) has been used to simulate wetting patterns of drip irrigation.For this purpose, Minitab 15 software has been used for calculation of MLR parameters.Resulted formulas for computation of radius and depth of wetting patterns are presented in Eqs.[15] and [16], respectively.where r = radius of wetting patterns [cm], z = depth of wetting patterns [cm], θ r = residual water content [-],

Comparison of the results of different methods
The results of two different GP models including GP fullset model with the functions +, -, ×, /, Log, Ln, Power, Sin, Cosine, Arctangent and GP 4basic model with the functions +, -, ×, / and MLR in estimation of radial and vertical depth of wetting patterns of drip irrigation are shown in Table 2.It can be concluded from Table 2 that the GP fullset has more precision in validation stage in comparison with GP 4basic and MLR models.In addition, the analysis showed that Drip wetting patterns with genetic programming there are not constant trends in the estimation of radial and vertical depth wetting patterns.For instance, when estimating radial depth wetting pattern from 375 points in validation stage, 131 estimated values by GP fullset have absolute differences less than 1 cm from corresponded HYDRUS data.On the other hand, regarding the estimation of vertical depth, 139 estimated values have absolute differences less than 1 cm.However, it should be noted that radial depths were slightly better estimated than vertical depths.This is because the variances between radial depths was less than the corresponded value related to vertical depths.In fact, the RMSE values in a case of radial depths is less than the corresponded RMSE value in a case of vertical depths.
Scatter plots of GP fullset, GP 4basic and MLR in validation stage for estimation of radial depth of wetting patterns are shown in Fig. 3.The estimates of the both GP models are closer to the corresponding observed values than those of the MLR model (Fig. 3).In other words, the MLR model has much more scattered estimates than those of the GP models.In Fig. 3, the fit line equations (assuming that the equation is y = ax + b) given in scatterplots indicates that a and b coefficients of the GP fullest and GP 4basic models are closer to the 1 and 0 than those of the MLR model.A comparison of two GP models reveals that the estimates of the GP fullset model seems to be closer to the exact (1:1) line than the GP 4basic model especially for radial distances from emitter values higher than 60 cm.The GP fullset model underestimates the observed high radial distance values while the GP 4basic model overestimates them.
On the other hand, in Fig. 4 it is illustrated the scatter plots of vertical distances from emitter (z) esti-mated by the GP fullest, GP 4basic and MLR models versus HYDRUS 2D in validation stage.Similarly to what reported for the radial distance from the emitter, the GP model estimates of the vertical distance are closer to the observed values than those of the MLR MLR models for irrigation durations of 0.5 and 1.0 h (Fig. 5).The estimated wetting patterns with GP fullset from beginning of irrigation to 1 h after it have not very good agreement with the observed curves.GP fullset and MLR models give similar estimates for the time duration of 1.5 h.In the case of 2 h after beginning of irrigation, the MLR model was closer to the observed wetting pattern than the other models for the x-y range of 20-30 cm.It can be concluded that the MLR model is the tool that better predicts the observed trends for irrigation durations longer than 2 h.These results suggest that in those situations, genetic programming might not give any significant advantage over the MLR procedure which can be used more easily.
Genetic programming combining either full set or four basic operators, or multiple linear regression method were used to simulate wetting patterns of drip irrigation systems in comparison with the more commonly employed HYDRUS software.The results are satisfactory and allow the users estimating wetting patterns dimensions for any given time, emitter discharge and soil hydraulic properties without having to perform a detailed numerical simulation using the HYDRUS software.GP fullset model performed better than the GP 4basic and MLR models for irrigation durations of 0.5, 1.0, 1.5 and 2.0 h.The MLR model was found to be better than the GP models for irrigation durations longer than 2 h.The comparison results in field experiments revealed that the GP method could be employed successfully in modeling wetting patterns of drip irrigation.

Figure 1 .
Figure 1.Scheme of the finite element grid used in the numerical simulations.

Figure 2 .
Figure 2. View of the developed program for computation of wetting patterns of drip irrigation by using Wolfram Mathematica 7.0.

Table 1 .
Parameters used for the genetic programming (GP) analysis of the soil wetted patterns

Table 2 .
Calibration and validation results of the GP fullest, GP 4basic and multiple linear regression (MLR) models in estimation of radial and vertical distance from emitter RMSE: root mean square error, R 2 : correlation coefficient.