Evaluating inverse distance weighting and kriging methods in estimation of some physical and chemical properties of soil in Qazvin Plain

327 Evaluating inverse distance weighting and kriging methods in estimation of some physical and chemical properties of soil in Qazvin Plain Sayed Roholla Mousavi a,*, Feraidon Sarmadian a, Somayeh Dehghani b, Mahmood Reza Sadikhani c , Abass Taati b a Laboratory of Remote Sensing and GIS, Department of Science and Soil Engineering, University of Tehran, Karaj, Iran b Soil Science Department, College of Agriculture, Shahrkord University, Chaharmahal Bakhtiari, Iran c Soil Science Department, College of Agriculture, Lorestan University, Lorestan, Iran


Introduction
Information about variability of the soil properties is very important in ecological modeling, environmental predictions, accurate farming and natural resources management (Lin et al., 2004).Spatial variability of soil input data can be highly effective on results of deductive, experimental and theoretical soil models (Wilding et al., 1994).Geostatistics is one of the very useful methods in estimation of spatial distribution of data.Spatial variability of soil properties or non-uniformity resulting from spatial differences in the observed soil properties includes two structural or non-structural components (Pang et al., 2011;Mohammadi, 2006).These variations result from both inherent process (soil constituent factors) and managerial process (such as fertilizer consumption, crop rotation and type of cultivation) in any spatial and temporal scale (Castrignanò et al., 2000).Estimation accuracy of each variable depends on two factors: 1-Selection of desirable interpolation method for obtaining soil properties in the non-sampled points, 2-suitable application of the interpolation methods considering nature and properties of data.The most common interpolation methods in agricultural researches are inverse distance weighing and kriging methods (Kravchenko and Bullock, 1999).Many researchers have compared IDW and kriging.In some cases, the performance of kriging was generally better than IDW (Kravchenko and Bullock, 1999;Kravchenko, 2003).Robinson and Metternicht (2005) used three inverse distance weighing and kriging and spline methods for estimation of electric conductivity and organic carbon.Gallichand et al. (1992) and Heisel et al. (1992) studied and evaluated different interpolation methods to prepare surface soil clays and found that weighted moving average (WMA), kriging and co-kriging methods are similar to each other.Reza et al. (2010) Evaluated and comprised of some soil chemical parameters by ordinary kriging and inverse distance weighting methods and they reported spatial variability for different soil parameters (except available nitrogen) may be better understood by ordinary kriging than by IDW method.Moustafa and Yomota (1998) measured soil salinity and prepared contour-maps for drainage projects with geostatistical methods such as kriging.Result of their studies showed that kriging method showed more accurate estimation than other methods.Laslett et al. (1987) conducted researches on study and evaluation of some Interpolation methods for estimation of surface soil pH.They used inverse distance weighing and kriging methods and used data obtained from Digital Elevation Model and climate for estimation of pH and finally recognized the kriging method as the most suitable method.Shakoori Katigari et al. (2011) reported the best semivariogram model for the fitted model on bulk densityas spherical model and range of influence of bulk density as 500m.Mohammadzamani et al. (2007) in a research in wheat farm managed by the local farmers in Sorkhan kalateh located in 25 km of Guilan Province reported the fitted models to bulk density, sand , silt , clay , pH, ESP, lime as spherical and reported electrical conductivity, soil saturation moisture and cation exchange capacity as Gaussian model.Importance of the spatial variations in soil properties is evident.Rosemary et al. (2017) using geostatistical spatial variability of some soil properties (EC, pH, OC, CEC and textures) in different applications examined in Sri Lanka.The researchers reported that the spatial dependence of different soil properties and spatial dependence class for EC and pH stronger won.Although factors of variations are different in various points, understanding sources of variations helps us do better management and considering that no studies have been conducted on spatial variability of soil properties in this region, the present research was conducted to study spatial variations of some physical and chemical properties of soil in Qazvin plain and two inverse distance weighing and kriging methods were used for zoning of points and two mentioned methods were evaluated for estimation of these properties.

Study area and experimental design
This research was conducted in Abyek City located in east of Qazvin province.Range of variations in elevation of the region varied from 1150 to 1450 m above open sea level (Figure 1).Range of variations of slope also varied from 1 to 15%.The region with approximate area of 16630 hectares is located in latitude of 36´1´´ to 36´9´´ and longitude of 50´21´´ to 50´14´´.It has mean annual precipitation of 257 mm and average annual temperature of 14 ˚C and the coldest month of the year is January and the hottest month of the year is July.For this purpose, geographical position of the sampling points was recorded with global positioning device at the beginning of sampling.Then, the soil samples were extracted from 62 observational points which are located in an irregular network with mean distance of 800 m from each other.After transferring the samples to the laboratory, they were dried under air.For the next tests, it was passed through 2-mm sieve.Schematic Position of the sampling points was presented in the studied region (Figure 2).

Physical and chemical properties measurements
Percent of sand, silt and clay particles was calculated with hydrometric method (Bouyoucos, 1962), bulk density of soil was calculated with sampling cylinder (Klute and Dirksen, 1986), and total porosity of soil samples was calculated with apparent and real specific mass of soil as a follow (Eq.1): Where ρb is bulk density (g.cm -3 ) and ρs is soil real specific mass (g.cm -3 ) and F is total porosity of soil (%).The pH in saturated extract was measured with pH meter (Thomas, 1996), electrical conductivity was measured with electric conductivity meter in saturated extract (Richards, 1954), Percent of calcium carbonate was measured with neutralization method with normal Hydrochloric acid 1 Normal and titration with sodium hydroxide (Nelson, 1982).

Geostatistical Analysis
Geostatistical Analysis was performed with the measured specifications in 62 points of the studied region.
To study frequency distribution and determine descriptive statistics of each variable, maximum and minimum statistics, mean, standard deviation, skewness and Kurtosis were determined.To study normality of data in confidence level of 95%, Kolmogorov-Smirnov test was used and to study normality of data, logarithmic transform method was used.To study spatial variations of the studied properties, experimental semivariogram was calculated for all variables and spatial structure of the data was studied in the studied zone.Variogram is a mathematical model and a vector quantity which shows spatial relationship between values of the measured variable in terms of squared value difference of two points and considering their distance and direction.Considering that it is not possible to calculate the whole studied population, semivariogram is estimated in a specified separation distance with the follow (Eq.2): is the measured value of variable in position (xi) and N(h) is the number of paired comparisons in distance of (h) in the studied zone.In structural analysis process, the standard models have been obtained considering values of (R 2 ) and Residual sum of squares (RSS).The model which had the maximum R 2 and minimum RSS was selected as the best semivariogram model.Then, the experimental semivariogram related to each one of the properties was calculated and different theoretical models including spherical and exponential models were fitted to them.Definition of the fitted models is as following equations (Eq. 3 and 4) (Mohammadi, 2006): (3) Where semivariogram include variables such as C0, nugget effect (the minimum value of semivariogram) which has been calculated and shows discontinuity of semivariogram curve adjacent to the origin coordinates reflecting variance of sampling errors and spatial variance in shorter distances from the sampling distance (Li and Heap, 2008).C1 is Variogram threshold, total nugget effect (Co) and structured part of semivariogram (C) which is equivalent to total variance of the studied variable and relatively constant value with random variations is called threshold), a is radius of influence or range (the distance in which semivariogram is fixed and reaches sill), h is separation distance.To determine different classes of spatial dependence of the studied variables, ratio of structured effect to total variance was used.if ratio of structured effect (C) to sill (C+C0) is less than 25%, between 25 and 75% and more than 75%, that specification will be in range of weak, medium and strong spatial dependence (Cambardella et al., 1994).To evaluate the estimator methods, several evaluation indices can be used concurrently.In this research, statistical indices of coefficient of determination (R 2 ), mean bias error (MBE), mean absolute error (MAE) were used (Mishra et al., 2010).MBE indicates deviation of the method and MAE indicates accuracy of each method.Positive values of MBE indicate overestimation of the model and its negative values indicate its underestimation.The manner of calculating statistics of (mean bias error (MBE) and mean absolute error (MAE)) is as follows equations (Eq. 5 and 6) (Wakernagel, 2003): where Z* (xi) is the estimated value in point X, Z (xi) is the observed value in point Xi and n is the number of points.In the next stage, the desired variables were estimated with kriging interpolation methods and inverse distance weighing method.Kriging is a weighted moving average and is defined as follows: (Eq.7) (Pang et al., 2011): And inverse distance weighing method is defined as following (Equation 8) Where Di is the i -th observed point distance to the estimated point, (α) is the inverse distance weighing power and (n) is the number of neighborhood points.The contour maps were drawn for each one of the physical and chemical properties measured in the studied region with results of suitable estimation method.All geostatistical operations were performed with GS + 5.1 software and properties zoning in Arc GIS9.3 software.

Statistical analysis
Statistical description of physical and chemical properties of the measured soil has been summarized in 62 points of the studied region (Table 1).Results of Kolmogorov-Smirnov test also showed that the properties have relatively normal distribution and results of skewness coefficient in Table 1 confirm the said results.Study of the physical and chemical properties of soil showed that mean percent of all three particles have relatively equal means and generally, the soil tissue of the region is clay, sandyclay -loam and sandy loam considering mean size distribution of the soil particles.Mean of the lime in the studied region was regarded as the lime soils based on the criterion provided by Pansu and Gautheyrou (2006) who states that the soils of which calcium carbonate is above or equal to 6.As Soil and Water Research Institute (1989) reported that the soils which have electrical conductivity of lower than 4 dS/m are in unrestricted range.Considering that insignificant area of the southern part of the region has electrical conductivity of above 4dS and with high restriction, other soils of the region will lack salinity restriction, also electrical conductivity has the highest CV, more than 35% (Table 1) this result shows that there was high variation between Electronical conductivity range in northern and southern part of study area.Based on classification by Pansu and Gautheyrou (2006), the studied soil is classified as a part of soils with low bulk density.Variability class of the variations coefficient has been obtained based on criterion presented Wilding (1985) such that if variations coefficient is below 15%, it will be in low variability class, if variations coefficient is between 15% and 35%, it will be in medium variability class and if variations coefficient is above 35%, it will be in high variability class.Considering the obtained results, bulk density, porosity and pH were included in low variability class and silt percent was included in medium class and sand, clay, electrical conductivity and lime were included in high variability class.Results also showed that the maximum and minimum variance related to percent of sand and Total porosity respectively.

Spatial correlation analysis
Results of fitting standard models to experimental variograms showed that all of these variables had spatial structure and the spherical and exponential model has had the best fitting with variations of the studied properties (Table 2).Considering results of this research, other properties than silt percent and total porosity (for which the best fitted model to semivariogram was exponential were spherical).Mohammadzamani et al. (2007) reported spherical model as the best fitted model to bulk density, sand and clay.The best fitted model for pH was aspherical (Yang et al., 2011).Variogram study of the studied properties in the studied region showed that their variability in the region was independent of the special geographical direction.Results show that the minimum nugget effect of the chemical and physical properties of the region relates to bulk density and total porosity of soil and this result shows that relative variance and sampling size are suitable for clarifying spatial structures and the highest nugget effect relates to sand percent indicating strong random variance in short intervals resulting from error, measurement and short range variations and the studied property in shorter distance from the shortest sampling distance.The smallest sill of the studied properties relates to total porosity and the highest hill relates to percent of sand indicating random or unstructured variations or total variance of all samples which have been applied in calculation of semivariogram.Oji et al. (2012) in their studies, reported that the maximum sill related to percent of sand (159.3%).
The maximum effective range among the studied properties is 21100 m related to bulk density, total porosity and electrical conductivity and pH indicating extensive influence of the mentioned variables on their adjacent points.With increasing effective range, sampling distances have increased and the numbers of necessary samples for sampling and necessary costs for sampling are reduced.The minimum radius of influence is 10740 m related to silt percent.Oji et al. (2012) reported that clay percent as the maximum radius of influence (2296 m).Results showed that the highest ratio of structured effect to semivariogram related to lime (0.88) indicating that 88% had spatial structure and only 12% was random and the minimum related to silt percent, bulk density and total porosity of soil (0.5) indicating that about 50% of total variations relating to the desired properties of the soil had spatial structure and about 50% of these variations was random and lacked spatial structure.Spatial dependence class was medium in all properties except lime and percent of clay (strong class) and other properties had medium spatial dependence.Strong spatial structure of the studied properties means that the use of geostatistical methods can be useful in analysis of the variability models of the studied variables.Strong spatial dependence of the soil properties can depend on inherent properties of soil (formation of soil) while medium spatial dependence is majorly attributed to external factors (soil management methods).Cambardella et al. (1994) concluded that there is a medium dependence for bulk density.López-Granados et al. (2002) also reported the strong dependence for silt in two depths of 0 to 10 and 25 to 35 cm.
As shown in the Table 3, the maximum coefficient of determination and the minimum MBE and MAE relate to total porosity of soil.Results show that kriging method shows better estimation of the studied properties than the inverse distance weighing method and the intended maps were zoned in the region with kriging estimator (Figure 3).Table3.Suitable method for estimation of physical properties of soil and the calculated criteria for evaluation and accuracy of the performed estimations in the studied region ** R 2 , MBE, and MAE are coefficient of determination, mean bias error, and mean absolute error, respectively.
The maximum and minimum quantity of sand are found in northeast and southwest of the region and the highest quantity of silt is found from west to center of the region and the lowest quantity is found in southwest and southeast.
It is assumed that soils with higher clay rate have been formed in physiographic units with low slope (alluvial piedmont plain and lowlands) and there has been enough chance of depositing the carried particles due to erosive processes.Physiographical unit of the alluvial piedmont plain has been due to accumulation of fine deposits with higher clay quantity to which water of upstream lands cause them to be transferred.The highest bulk density is found in southeast and we also find minimum total porosity in southeast.Considering high clay quantity in south of the region, bulk density of soil is low in this part.considering increase of electrical conductivity in southeast of the region, we find the highest pH in this part and also soil texture of the region has become finer in this part and the main reason for such results is formation of soils in this region in physiographical unit of the lowlands and one of the major properties of soils in these regions is observation of the mentioned characteristics.The minimum electrical conductivity and pH are observed in northeast of the region.Soils of these regions are alternatively affected by flood flows of the region due to location in hills with high slope (10-16%).For this reason, there has been shorter time for accumulation of salts due to soil erosion.Equivalent calcium carbonate map in the studied zone indicates calcification of the soils in the studied region.The highest quantity of the equivalent calcium is found in west to south and the minimum quantity is found in east of the region.This can be related to role of geogenetic factor of parent material considering that soils with parent material of calcareous marl were found in these regions.The properties interpolation results showed that the kriging estimator better can show variability of the studied properties in the region in comparison to IDW method (Smith and Halvorson, 2011).Also in another study, similar results were achieved by Zare-Mehrjardi et al. ( 2010) who reported that the kriging method for estimating spatial distribution of soil properties better than the inverse distance weighting method.Thereby, kriging method for all properties has the highest value for R 2 and has lowest value for MAE and MBE.At the end, considering the best interpolation method, spatial variability map of each of the properties was prepared in ArcGIS 9.3 software.

Conclusion
The present research showed that Geographical Information System and Geostatistical System had high ability to spatial distribution and zoning of some physical and chemical properties of soil.Results of fitting standard models to experimental variograms showed that all of these variables had spatial structure and the spherical and exponential model had the best fitting to variations of the studied properties.The highest coefficient of determination and the minimum MBE and MAE related to total soil porosity.Considering the applied error indices, ordinary kriging has acted as a suitable estimator in zoning of the studied properties in the studied region.Spatial structure and range of influence of variables were highly affected by non-inherent variability and managerial factors.Accuracy of the resulting maps relating to soil properties showed suitable agreement with observational results.Therefore, Variogram and its related properties can be applied as an efficient tool for design of sampling networks and identification of managerial zones in accurate farming.For this reason, soil can be managed better for saving agricultural inputs and environmental protection with geostatistical technique and kriging technique and zoning of farms and creation of isolated zones.It is necessary to note that results of this research can be used in the studied region and cannot be generalized to other regions and spatial scattering of soil variables should be studied for each region.

Table 1 .
Descriptive statistics for soil Physical and chemicals properties.Descriptive statistics BD, dry bulk density; F, Total porosity of soil; EC, soil saturated paste electrical conductivity; CCE, calcium carbonate equivalent.*Min., minimum value; Max., maximum value; CV, coefficient of variation.

Table 2 .
Cambardella et al. (1994)fitted model to semivariogram of some physical and chemical properties of soil measured in 62 points of the studied region EC, F, BD, clay, silt, sand * are percent of sand, silt, clay, bulk density (gr.cm -3 ), total porosity percent, electrical conductivity (dS m -1 ), pH, equivalent calcium carbonate (percent).**ififratio of structured effect (C)to sill (C+C0) is less than 25%, between 25 and 75% and more than 75%, that specification will be included in weak, medium and strong spatial dependenceCambardella et al. (1994).RSS is Residual sums of squares and R 2 is coefficient of determination (R 2 ).