Modelling infiltration and geostatistical analysis of spatial variability of sorptivity and transmissivity in a flood spreading area

Knowledge of infiltration characteristics is useful in hydrological studies of agricultural soils. Soil hydraulic parameters such as steady infiltration rate, sorptivity, and transmissivity can exhibit appreciable spatial variability. The main objectives of this study were to examine several mathematical models of infiltration and to analyze the spatial variability of observed final infiltration rate, estimated sorptivity and estimated transmissivity in flood spreading and control areas in Ilam province, Iran. The suitability of geostatistics to describe such spatial variability was assessed using data from 30 infiltration measurements sampled along three lines. The Horton model provided the most accurate simulation of infiltration considering all measurements and the Philip’s two-term model provided less accurate simulation. A comparison of the measured values and the estimated final infiltration rates showed that the Kostiakov-Lewis, Kostiakov, and SCS models could not estimate the final infiltration rate as well as Horton model. Estimated sorptivity and transmissivity parameters of the Philip’s two-term model and final infiltration rate had spatial structure, and were considered to be structural variables over the transect pattern. The Gaussian model provided the best-fit theoretical variogram for these three parameters. Variogram values ranged from 99 and 88 m for sorptivity and final infiltration rate to 686 (spherical) and 384 m (Gaussian) for transmissivity. Sorptivity, transmissivity and final infiltration attributes showed a high degree of spatial dependence, being 0.99, 0.81 and 1, respectively. Results showed that kriging could be used to predict the studied parameters in the study area.


278
F. Haghighi-Fashi et al. / Span J Agric Res (2014) 12 (1): 277-288 2005). Infiltration rate and soil sorptivity are relevant to water movement in the soil medium and transmissivity provides information about relevant hydraulic and structural soil properties. Sorptivity which is a function of soil suction is important for understanding hydraulic properties. It is defined as a physical quantity that shows the capacity of porous media for capillary uptake and release of water in the soil (Philip, 1957). This parameter is also def ined as the relationship between soil hydraulic properties and other more easily measured properties such as porosity (Bouma, 1989). Haghighi et al. (2010) found that sorptivity exhibits spatial variability under field conditions. They found that the spatial variation of sorptivity showed no clear pattern in the study region, which was in accordance with the f indings of Machiwal et al. (2006). Transmissivity in the Philip's model may approximate to saturated hydraulic conductivity (K s ) over long periods (Swartzendruber & Youngs, 1974).
An important issue in infiltration modeling is how to describe spatial variability in relevant soil properties. Previous studies have shown large spatial variability in infiltration behavior of agricultural land (Nielsen et al., 1973;Carvallo et al., 1976) and have presented advances in the use of spatial statistics related to soil variability (Burrough, 1993;Machiwal et al., 2006). The coefficient of variation (Jury et al., 1991;Wollenhaupt et al., 1997) is one approach employed in the study of spatial variability. Other technique associated with the spatial variability includes soil mapping unit classification (Mulla & McBratney, 2002). Soil hydraulic properties can vary appreciably even within similar soil mapping units (Comegna & Vitale, 1993). Scaling is one approach that has been used to describe the spatial variability of soil properties (Peck et al., 1977;Machiwal et al., 2006). The geostatistical technique was applied to analyze the spatial variability of soil sorptivity (Sepaskhah et al., 2005) and soil texture (Camarinha et al., 2011). In geostatistic analysis, a variogram or semi-variogram is employed to determine spatial variability of a variable, and provide the input parameters for the spatial interpolation of kriging (Webster & Oliver, 2007). The concept of a geostatistic method has been described by Isaaks & Srivastava (1989) and Kitanidis (1997). The semi-variogram model provides an estimate of the model required to characterize the spatial pattern of the variables. The goal of variogram analysis is to fit a model to the points that form the empirical semi-variogram. The kriging approach is suitable to investigate the spatial variability of soil water properties that can be estimated with a low sampling cost to make management decisions (Vieira et al., 1981).
The spatial variability of infiltration and related parameters (sorptivity and transmissivity) can be evaluated and predicted using geostatistical techniques (Cressie, 1993;Orjuela-Matta et al., 2012). The objectives of the present study were to model the infiltration rate of flood spreading and control areas and evaluate common models to estimate the f inal infiltration rate. In addition, geostatistical analysis was performed on steady-state infiltration rates, estimated sorptivity, and estimated transmissivity parameters of a Philip's two-term model in the flood spreading and control areas of the province of Ilam in Iran.

Characteristics of the study flood spreading area
The data used in this study are from a study by Soleimani et al. (2011) on the Dehloran flood spreading area, Ilam province, Iran (32°27' to 32°35' N latitude and 47°25' to 47°42' E longitude, elevation of 200 m) (Fig. 1). The study area encompassed 5200 ha. The flood spreading system is one of the suitable methods for aquifers recharge, flood control, and flood water conservation containing sediment load. This flood spreading area is classified as Entisol. The average annual precipitation is approximately 235 mm and the average annual temperature is 26°C. The climate is arid with a hot summer. Average annual evaporation is 4300 mm. The topography of the area is moderately flat. The soils have low fertility and poor capacity for moisture retention. The texture of the soil is sandy based on USDA soil classification system. Current vegetation types in the area include Gypsophila stipa, Stipa peteropyron, Stipa astragalus, and Cornalaca alhaji.

Experimental design and infiltration tests
A total of 30 measurement sites (27 from floodspreading areas and three from control areas (which were areas that did not receive flooding), were selected over the study area to model inf iltration and for geostatistical analysis. The inf iltration data were measured on transect (3 lines) pattern (1500 m in size) over the study area (Fig. 1). Three replicate infiltration measurements were made at each of the 30 measurement sites. The first, second and third lines of the flood spreading system and the control line were selected. A total of nine measurements were made along each line of the flood spreading system. The measurement sites were considered to be four triangles as in Fig 1 and that in one triangle, there was a 4 m circle and those three points 90 degrees to one another were selected as the three sampling points. A base center and circular radius of 4 m was considered to select the f irst, second, and third inf iltration measurement points (i 1 , i 2 , i 3 ). The mean infiltration rates measured on the transect (line) provided an opportunity for analysis of the spatial variability of the measured inf iltration and estimated sorptivity and transmissivity of the Philip's two-term model (Philip, 1957).
Inf iltration was measured using a double-ring infiltrometer (Bouwer et al., 1999). The diameter of the inner ring was 30 cm and that of the outer ring was 60 cm with a height of 50 cm. The decreasing rate of water was measured in the inner ring at time intervals of 1, 2, 3, 5, 10, 15, 20, 25, and 30 min. The initial soil water content at the sites ranged from 2.4% to 14.9%. The optimal model-estimated final infiltration rate, sorptivity and transmissivity of the Philip's two-term model, and other model parameters were determined with the least square method using MATLAB TM (The MathWorks Inc, Natick, MA, USA).

Estimation of infiltration model parameters
In the present study, the data from each of the 30 sites across the study area were fitted to commonlyused inf iltration models being the Kostiakov (Kostiakov, 1932), Kostiakov-Lewis (Mezencev, 1948), Horton (Horton, 1940), Soil Conservation Service (SCS) (NRCS, 1974), and Philip's two-term (Philip, 1957) models, and the accuracy to estimate infiltration was evaluated. A summary of the infiltration models applied in the current study follows.
-Kostiakov model -Kostiakov-Lewis model. The modified model of Kostiakov for long time is defined as: where i(t) is the infiltration rate (LT -1 ) as a function of time, a and b are the equation ' s parameters (a > 0 and 0 < b < 1) and i c is the steady infiltration rate (LT -1 ).
-Horton's model. The Horton ' s infiltration model (Horton, 1940) is expressed as follows: where i c is the steady infiltration rate (LT -1 ), i 0 is the initial infiltration rate (LT -1 ), and t is time (T). β is the inf iltration decay factor. In this study, the Horton model has been applied in the form of Equation [4].
where I is the cumulative infiltration rate (L) and c, m, and a are the constant parameters and t is time (T).
-Philip's two-parameter model. The Philip's twoterm model (Philip, 1957) is defined as: where S is sorptivity (LT -0.5 ), that is function of soil matric potential, A is transmissivity factor (LT -1 ) as a function of soil properties and water content, and t is time (T).
-SCS model. Soil Conservation Service (SCS) is expressed as follows: where a and b are the equation ' s constants. I is the cumulative infiltration rate (L).

Best-fit infiltration model
The coefficient of determination (R 2 ) and root mean squared error (RMSE) were used to evaluate the bestfit infiltration model for the study area. The observed and estimated final infiltration rates and the ability of the models to estimate the model parameters were evaluated.

Geostatistical analysis
Analysis of the spatial characteristics of the variables (measured steady infiltration rate, estimated sorptivity, and estimated transmissivity) was performed by kriging to assess the viability of spatial interpolation of the parameters. The variogram (Á(h)) is a major tool in geostatistics to describe the spatial relationship between neighboring observations. The variogram is defined as follows: [7] where N(h) is the number of pairs of points that are separated by distance h (Yates & Warrick, 1999), z is the variable or parameter to be interpolated, and x i and x i + h are the values of each pair of points.
For geostatistical analysis, a variogram value can be interpolated for any sampling interval. The commonly-used models are linear, logarithmic, spherical, exponential, Gaussian, and pure nugget effect (Isaaks & Srivastava, 1989). Webster & Oliver (2007) reviewed the characteristics and conditions of these models. They have common parameters, including the nugget (C 0 ), range (a), and sill (C 0 + C). The nugget represents discontinuity of the variogram near the origin and non-spatial variation. The range is a measure of the maximum distance at which the constant correlation is reached. The sill is the value of the semi-variance where the model is stabilized and has a constant value.
After assessing the spatial relationship between neighboring observations, the kriging approach was employed for spatial prediction at sites not sampled and to build contour maps (Sepaskhah et al., 2005). The kriging equations are defined as follows: where S*(x p ) is the krigged value at location x p ; S(x i ) is the known value at location x i ; λ i is the weight related to the data; μ is the Lagrange multiplier; and γ(x i , x j ) is the value of the variogram corresponding to a vector with its origin in x i and extremity in x j . GS+ (Gamma Design Software, Plainwell, MI, USA) and ArcGIS (Redlands, CA, USA) geostatistical analyst were used for geostatistical analysis in this research. The best model was selected based on the highest R 2 and the lowest RSS. The degree of spatial dependence was assessed, based on the relationship between the nugget and the sill (C / C 0 + C). Spatial dependence is classif ied as weak when it is <0.25, moderate when it is 0.25 to 0.75, and strong when it is >0.75 (Cambardella et al., 1994).

Selection of best-fit infiltration model
Parameters values for the selected inf iltration models and goodness-of-fit parameters are shown in Tables 1 and 2, respectively. The mean final infiltration rate for the 30 infiltration measurements was 0.15 cm min -1 with a minimum of 0.08 cm min -1 and maximum of 0.27 cm min -1 in the sandy soils of study areas ( Table 1). The lowest RMSE was obtained for the Horton model at nearly all sites, demonstrating that the Horton model was the best-fit infiltration model for the entire study area.
The values of the final steady infiltration and RMSE showed that the estimated inf iltration rates by the Horton infiltration model were closer to the measured values at the selected flood spreading and control areas (Table 1 and 2). The results showed that the Kostiakov-Lewis model inaccurately estimated the final steady inf iltration rate at most sites. A comparison of the measured values and the estimated final infiltration rates showed that the Philip's two-term, Kostiakov, and SCS models were less accurate than the Horton model; which was supported by the higher RMSE values compared to the Horton model ( Table 2). Table 1 showed that the range of variation of estimated transmissivity was large, which demonstrated considerable spatial variability in the study area. Philip's model parameter A (transmissivity) at 10 sites was found to be 0 (Table 1). These f indings demonstrate that the Philip's two-term model did not fit the observed infiltration data well at these 10 sites. A at these 10 sites was negative. The lower limit for the transmissivity estimation was considered to be 0 by default. Despite the negative A, estimated sorptivity (S) ranged from 2.80 to 6.20 cm min -0.5 and A from 0.010 to 0.137 cm min -1 at the different sites.

Properties of sorptivity and transmissivity
Estimated transmissivity values from the Philip's two-term infiltration model from 20 sites were used to determine the spatial variability because the values Haghighi-Fashi et al. / Span J Agric Res (2014) 12 (1): 277-288 from the other 10 sites were negative. Overall, the variation in estimated sorptivity and transmissivity values over the study area showed spatial structure. The relationships between estimated sorptivity and estimated transmissivity were not signif icant (R 2 = 0.31) and there was no trend observed for an increase in sorptivity as transmissivity increased (Fig. 2).

Geostatistical analysis
Data analysis showed that all data had a symmetric distribution. All attributes indicated a normal distribution and the mean and median values were similar (Table 3). The skewness was between -1 and +1 and the kurtosis coeff icient was close to 3, demonstrating a behavior that approaches normal distribution (Grego & Vieira, 2005) (Table 3). Fig. 3 shows the contour maps for attributes obtained using interpolation. The contour lines are more intensive for the final infiltration rate and the variability of its values is greater than for the estimated sorptivity and transmissivity values (Fig. 3). The spatial correlation of final infiltration rate, estimated sorptivity and estimated transmissivity was determined using the appropriate semi-variogram. Fig. 4 shows the best-fitted variograms for the data series. These figures show that estimated transmissivity revealed a spatial structure for both Gaussian and spherical models; estimated sorptivity and measured final infiltration rate revealed spatial structures for the Gaussian model, indicating that the spatial structure was determined for the experiments.
The properties of the fitted semi-variograms are shown in Table 4. The number of lags and the variogram model type were determined using trial and error. The range of variograms was 98.4 m and 87.9 m for sorptivity and the final infiltration rate and 384 m (Guassian) for transmissivity. Based on the Table 4, all three variograms showed a spatial structure and the range values for these variograms are the same those indicated in their corresponding figures. The predominant variogram model was Gaussian, followed by spherical, exponential, and linear models.   The R 2 values for the present study for estimation of transmissivity, sorptivity, and final infiltration rate were 0.82, 0.68 and 0.85, respectively. The parameters showing the least and highest coeff icients of determination were for sorptivity and final infiltration rate, with values of 0.68 and 0.85, respectively. The R 2 for final infiltration rate demonstrates a good fit with the data of the semi-variogram.

Best-fit infiltration model
The results show that the Horton model was the most accurate model for the flood spreading and control areas. The suitability of the infiltration models was tested under different conditions, because final infiltration rate is generally a soil-dependent parameter (Singh, 1992;Mishra et al., 2003). The applicability of the Kostiakov and Philip's models to estimate the infiltration rate is site-specific (Mbagwu, 1993). The results from the present and from previously mentioned studies indicate that not all models are applicable in all soils. The different RMSE for each model can be attributed to regional soil variation at the test sites, such as in texture class because of the dependency of infiltration rate on soil texture (Haghighi et al., 2010).
The complexity of field conditions can appreciably influence the relationships between observed and estimated values. Spatial variation and preferential  C / (C 0 + C): degree of spatial dependence; R 2 : coefficient of determination; RSS: residual sum of squares.
flow can affect the observed values and the model parameters. The poor performance of the Kostiakov-Lewis model in estimating the final infiltration rate may be because it takes longer to obtain a steady infiltration rate (Navar & Synnott, 2000;Haghighi et al., 2010). In addition, it seems that a period of more than 2 h may be required for application of the Philip's two-term model (Philip, 1969;Mbagwu, 1993). In general, spatial soil variation is a very important consideration in infiltration modelling.

Sorptivity and transmissivity
Negative values for estimated transmissivity mean that, after a certain point of time, the infiltration rates estimated by the Philip's equation became negative, demonstrating that extra filtration occurred or that the ground exuded water (Jaynes & Gifford, 1981;Machiwal et al., 2006). It is not possible to physically interpret the negative value for transmissivity; the 10 test sites, where these negative values were determined, were assessed for possible reasons for the negative model parameter. It was expected that sites having negative transmissivity would show the lowest basic infiltration rate (Machiwal et al., 2006); however, the results did not prove this assumption.
The parameters of sorptivity and transmissivity accounted for the effects of factors such as macropores, initial soil water content, soil texture, and soil structure that affect infiltration (Sharma et al., 1980). The variability of these soil factors highly influenced A and S. The results showed that the relationships between estimated sorptivity and estimated transmissivity were not significant (R 2 = 0.31) (Fig. 2). This is in agreement with Haghighi et al. (2010), who found a poor relationship between A and S (R 2 = 0.127). Talsma (1969) also found poor agreement between saturated hydraulic conductivity and sorptivity, which agrees with the present study.
In summer, A and S are good estimations of infiltration rates and are useful for planning irrigation systems and for management of soil and water resources necessary for precision agriculture. The findings of the present study can be used for planning and management of flood spreading areas in Iran.
The relationship between sorptivity and transmissivity indicated no distinct pattern in the study region. This finding is in agreement with results reported by Machiwal et al. (2006) and Haghighi et al. (2010). The variation in soil effective porosity, pore size distribution and soil matric potential are important factors in the spatial variability of Philip's model A and S (Sharma et al., 1980) and should be considered for spatial variability analysis of these parameters. The relationship between A and S and porosity is a topic for future study. In addition, non-uniform initial soil moisture within the test sites can cause appreciable variation in model parameters.

Geostatistical analysis
The R 2 values for estimation of transmissivity, sorptivity, and final infiltration rate were acceptable. This result agrees with the findings of Sepaskhah et al. (2005) who found that R 2 = 0.81 is acceptable only for estimation of sorptivity. Russo & Bresler (1981), Gupta et al. (1994) and Sepaskhah et al. (2005) reported variogram ranges of 35, 55 and 3.42 m for sorptivity, respectively, which proves that soil properties vary significantly at the different study locations. Warrick & Nielsen (1980) and Ramírez-López et al. (2008) found that sorptivity showed medium variability, while the cumulative infiltration indicated high variability. The soil of their study area was an oxisol and the region experiences high temperatures and rainfall in excess of 2000 mm yr -1 , which contrasts with the present study. These parameters have demonstrated high variability in other studies (Rodríguez-Vásquez et al., 2008;Cucunubá-Melo et al., 2011).
The semi-variogram model for the three properties was Gaussian and spherical. The estimated sorptivity indicated a nugget effect (0.14), demonstrating discontinuity of the variogram close to the origin. Where the contour lines are intensive, the variability of values is more severe. The variability of the final infiltration rate was higher than for sorptivity and transmissivity (Fig. 3). The random behavior of this parameter can be attributed to variability of soil properties, although statistical analysis showed a nearly normal behavior for this attribute.
The spatial variability of soil physical properties (Russo & Bresler, 1981;Camarinha et al., 2011), especially soil pore size distribution, are major factors in spatial variation of the final infiltration rate, Philip sorptivity and transmissivity (Sharma et al., 1980;Machiwal et al., 2006;Haghighi et al., 2010) in the flood plain and control areas of the studied region. Hallett et al. (2004) stated that macro-pores and repellency is the likely cause of severe spatial variability.
The estimated sorptivity and transmissivity and the measured final infiltration attributes showed a high degree of spatial dependence at 0.99, 0.81 and 1, respectively. The closeness of spatial dependence to one demonstrated the reliability of the fit of data to theoretical semi-variograms, which is necessary for the interpolation of variables in the contour maps (Orjuela-Matta et al., 2012).
In conclusion, contour maps allow identification of soil sectors where differences in values of estimated water sorptivity and transmissivity and measured final inf iltration rate exist. This identif ication helps to optimize irrigation management and to reduce surface runoff (Orjuela-Matta et al., 2012). The results showed that the kriging estimator was able to estimate the sorptivity, transmissivity and final infiltration rate if the variables had spatial variability.

Selection of best-fit infiltration model and geostatistical analysis for water management
Selection of best inf iltration model is a major component of irrigation efficiency. The present study showed that the Horton model can be used to adequately estimate the final infiltration rate in the study area. The Philip's sorptivity was site-dependent and varied widely between sampling points. Usually, in irrigation projects, a large amount of water is applied with low irrigation efficiency because of the spatial variability of soil hydraulic parameters such as transmissivity, sorptivity and final infiltration rate that influence the infiltration model. Dividing the area into sub-areas, each with a unique infiltration model, may result in better water management (Sepaskhah et al., 2005). As shown, the contour lines for the studied parameters were intensive for the final infiltration rates compared to estimated transmissivity and sorptivity. In conclusion, dividing the areas into sub-areas with relatively similar final infiltration values and two other parameters may help to manage water properly and design an efficient irrigation system and flood spreading area plan. The use of kriging in large areas to extrapolate the values related to inf iltration by generating contour maps can be a suitable tool for managing water resources and to partitioning the areas.
In the present study, the goodness of the fit of five models and their ability to estimate the final infil-tration rate of soils was assessed for a flood spreading area in Ilam, Iran, using root mean squared error (RMSE). The results showed that the Horton model can be applied to estimate the final infiltration rate in the study area. The Kostiakov-Lewis model did not estimate the final infiltration rate at most sites. The spatial variability of the infiltration process was also analyzed. The spatial variability of the final infiltration rate, estimated sorptivity, and estimated transmissivity were evaluated for soil in the flood plain and in nonflooded, control conditions. The attributes were analyzed for spatial variability and contour and prediction maps were obtained by kriging for the parameters fitted to the theoretical models of the semivariograms. It was found that final infiltration rate, estimated sorptivity and estimated transmissivity are variable parameters under field conditions.
The results showed that the final infiltration rate, estimated sorptivity, and estimated transmissivity had spatial distributions. The variograms ranged from 99 m and 88 m for estimated sorptivity and final infiltration rate and 384 m (Gaussian) for estimated transmissivity, respectively. The f inal inf iltration rate, estimated sorptivity and estimated transmissivity showed high degrees of spatial dependence and 0.99, 0.81, and 1, respectively. Geostatistical analysis helped to identify the variability of f inal inf iltration rate, estimated sorptivity and estimated transmissivity, and allowed prediction of them in non-sampled zones using kriging as an interpolation algorithm. It is suggested that the geostatistical approach obtains practical information that can aid in the management of water and soil resources in agricultural production and the operation of irrigation systems.