Evaluation of Geostatistical Techniques for Mapping Spatial Distribution of Soil PH , Salinity and Plant Cover Affected by Environmental Factors in Southern Iran

The study presented in this paper attempts to evaluate some interpolation techniques for mapping spatial distribution of soil pH, salinity and plant cover in Hormozgan province, Iran. The relationships among environmental factors and distribution of vegetation types were also investigated. Plot sampling was applied in the study area. Landform parameters of each plot were recorded and canopy cover percentages of each species were measured while stoniness and browsing damage were estimated. Results indicated that there was a significant difference in vegetation cover for high and low slope steepness. Also, vegetation cover was greater than other cases in the mountains with calcareous lithology. In general, there were no significant relationships among vegetation cover and soil properties such as pH, EC, and texture. Other soil properties, such as soil depth and gravel percentage were significantly affected by vegetation cover. Moreover, the geostatistical results showed that kriging and cokriging methods were better than inverse distance weighting (IDW) method for prediction of the spatial distribution of soil properties. Also, the results indicated that all the concerned soil and plant parameters were better determined by means of a cokriging method. Land elevation, which was highly correlated with studied parameters, was used as an auxiliary parameter.


Introduction
Soil quality preservation is one of the most important issues in sustainable management of ecosystems.In addition, identification of the spatial distribution of soil characteristics has a vital role in most of the bioenvironmental systems.Variability is one of the intrinsic characteristics of soil quality and within the same ecosystem soil properties may show significant spatial variations (Robinson and Metternicht, 2006).These variations are mainly arising from factors and processes of pedogenesis and land use (Ersahin, 2003).Given that classic statistics consider soil datasets as independent, their implementations lead very often to unrealistic results (Hasany-Pak, 1998).For this reason, we opted to apply geostatistical methods in order to better understand spatial variations of soil characteristics.
Under natural conditions, i.e. without anthropogenic intervention, vegetation cover is determined only environmental factors, such as the local climate, the geologic and soil background, and topography, e.g.aspect, slopesteepness, elevation and terrain position.Soil properties can strongly influence plant growth (Hoveizeh, 1997;Jafari et al., 2004;Taghizadeh Mehrjardi et al., 2009).Jenny (1980) explained that soil parent materials and time of soil development are the two main factors powerfully affecting the plant composition.In addition, soil chemistry has been reported to affect plant species composition through levels of salinity (Sharma and Shankar, 1991;Abbadi and El Sheikh, 2002), pH, calcium and organic carbon (Abd El-Ghani et al., 2002).Low species richness has been recorded where levels of salinity and CaCO 3 are high (Abd El-Ghani et al., 2002), while species composition and richness patterns have also been linked to trends in community development through time (succession) as soil conditions for plant growth improve in association with increasing soil depth, organic matter and water holding capacity, and decreasing pH and CaCO 3 (Shaukat et al., 1981).Elevation, slope and aspects are three environmental factors affecting soil and climate.Therefore the role of elevation, slope and aspects in soil properties such as its depth, moisture and organic matter content is highly complex ( Jafari et al., 2004;Ashcroft, 2006).
Today different geostatistical techniques are used widely for prediction of spatial variations of the soil properties.Hosseini et al. (1994) investigated the spatial variations of soil salinity and alkalinity in Alberta region.The results showed that the kriging and cokriging methods were suit-for specific soil and crop properties in Hormozgan province of Iran.
b) Mapping spatial distribution of soil pH, salinity and plant cover in Hormozgan province of Iran

Study area
The study area (BandarAbbas region in Hormozgan province) is located in northern costal zones of the Persian Gulf in south of Iran (Fig. 1).The extent of the study area was about 18,000 ha and was situated in mountain and hill area having elevation ranging from 860 m to 3,081 m.Annual rainfall is 214 mm, mean temperature is 24.33°C, average maximum temperature is 31.25°C,and average minimum is 17.35°C.
The soils of the study area are mostly shallow and therefore are only suitable as rangeland for herd grazing.According to the American system of soil taxonomy (USDA, 2003), the soils of the region are classified in as Entisols.
The area was divided into tree main units based on geomorphology.A terrain map unit (TMU) was created by stratification the study area into relatively homogenous areas that are Terrain Units.Then each main unit subdivided into sub units based on lithology and morhpology.The final TMU shows the distribution of 15 land unit in the study area (Tab.1).Stratification was achieved with stereoscopic aerial photo interpretation prior to fieldwork.The preliminary boundaries were then corrected during the fieldwork.

Sampling
In this study the quadrate plots were used for measuring vegetation parameters.For each TMU unit 4 to 19 plots were selected using a stratified sampling method.The observation sites were established with the help of topography, TMU and vegetation type maps.The sizes of samples were chosen based on a minimum area that was 100 m 2 (10×10 m).Within each plot 3 plots of 1 m 2 were selected.For determining the number of samples, a TMU map was crossed with vegetation type map and three plots (10×10 m) taken in each unit.Totally, 104 plots with 100 m 2 size were taken.At each plot (100 m 2 ) a relieve data sheet was filled in.Then, a visual estimate of browsing damage, percentages for grass cover, forbs cover and shrub cover was recorded.Three 1 m 2 plots were randomly taken inside the main plots (100 m 2 ).Percentage cover of each species of vegetation type was measured in each plot (1×1 m).Parameters of slope such as slope length, slope shape, slope percent and aspect were measured by compass and altimeter in plots of 100 m 2 .One soil sample was taken for measurement of EC, pH, texture and colour of soil.
Soil samples to laboratory the samples were air-dried at room temperature and then squashed.The dry samples were ground to pass through a 2 mm sieve.Although, in some of the experiments such as the determination of able for estimation of soil salinity (EC) and sodium adsorption ratio (SAR), respectively.Mohammadi (2000) determined some of the topsoil properties including salinity, water at saturation percentage, sodium adsorption ratio and percentage of lime using geostatistical predictors and Landsat TM sensing-numerical information as a secondary variable.Results showed that geostatistical predictors had relative superiority to the equations having linear correlations.Sokoti et al. (2006) used different geostatistical methods to predict the soil salinity distribution in Urmia plain of Iran.Results indicated that the kriging method having a gaussian model had the highest accuracy among other geostatistical techniques for estimating of salinity level in areas without any information.Similarly, Mc Bratney et al. (2003) developed some comprehensive maps for physical, chemical and biological soil properties in big areas of Australia by means of geostatistic methods, Geographical Information System (GIS) and Remote sensing techniques.Meul and Van Meirvenne (2003) used four techniques including the ordinary kriging, comprehensive kriging, simple kriging, and cokriging methods for estimation of the silt content in Belgium.They also used Digital Elevation Model (DEM) as a secondary variable.Results showed that the comprehensive kriging method had the lowest estimating error between all studied interpolation techniques.Ersahin (2003) used kriging and cokriging methods for investigation of spatial variations of infiltration rate in North West of Tookat of Turkey.The soil bulk density was used as an auxiliary variable in cokriging method in this study.Results illustrated that the cokriging method was a suitable technique for estimation of infiltration rate.Robinson and Metternicht (2006) used three different techniques including IDW, cokriging, and spilain methods for predicting the level of soil salinity, acidity and organic matter in South West of Australia.Results showed that the cokriging and spilain methods were the best techniques for estimation of the soil salinity level and organic matter content.Results also showed that the IDW method was suitable for prediction of soil acidity level.
A large part of Iran is located in the arid and semi-arid region, where low and erratic rainfall is the most outstanding characteristic of the land.Because of lack of accuracy or reliable information about vegetation, most of the watershed management plans cannot be applied successfully.Therefore, knowledge of the vegetation state is very important for watershed management.Considering all the above, this research was conducted in the framework of improving the management in rangelands towards soil conservation and prevention of land degradation in Iran.The main aim of the study was to determine significant relationships among vegetation parameters and environmental factors.This overall aim was attempted through the following specific objectives of the study: (a) Evaluation of the performance of different interpolation techniques including kriging, cokriging and IDW percentage of organic matter the soil samples were passed through a 70 mesh sieve (having diameter less than 0.5 mm).The majority of laboratory studies were performed according to Methods of Soil Analysis presented by Sparks et al. (1996).Organic carbon was determined by Walkley-Black method (Nelson and Sommers, 1982); carbonate and bicarbonate were determined by Bernard's calcimetric method; pH in aqueous suspensions and 1 M KCl (Peech, 1965); CEC by the method of Chapman (1965); EC of the saturated extract by the Bower and Wilcox (1965) method.Also, soluble (meq l -1 ) and exchangeable basic cations, CaCO 3 , gypsum, SAR and ESP were measured with respect to standard methods (Page, 1986).Particlesize distribution was determined after dissolution of Ca-CO 3 with 2 M HCl and decomposition of organic matter with 30% H 2 O 2 .After repeated washing to remove salts, samples were dispersed using sodium hexametaphosphate for determination of sand, silt and clay fractions by the pipette method (Day, 1965).Field data were entered into Excel software.Exploratory data was analyzed using box plot, correlation, regression, and Kruskal-Wallis test for finding the relationship between environmental variables and vegetation.Spatial prediction methods Geostatistical prediction includes two stages: (a) Identification and modeling of spatial structure where continuity, homogeneity and spatial structure of a given variable is studied using a variogram.
(b) Geostatistical estimation using kriging technique which depends on the properties of the fitted variogram which affects all stages of the process.

Kriging
The presence of a spatial structure where observations close to each other are more alike than those that are far apart (spatial autocorrelation) is a prerequisite to the application of geostatistics (Goovaerts, 1999;Robinson and Metternicht, 2006).The experimental variogram measures the average degree of dissimilarity between unsampled values and a nearby data value (Deutsch and Journel, 1998), and thus can depict autocorrelation at various distances.The value of the experimental variogram for a separation distance of h (referred to as the lag) is half the average squared difference between the value at Z(xi) and the value at Z(xi+h) (Lark, 2000;Robinson and Metternicht, 2006) (Eq.1): Where: n(h) is the number of data pairs within a given class of distance and direction.If the values at Z(xi) and Z(xi+h) are auto correlated the result of Eq. (1) will be small, relative to an uncorrelated pair of points.From analysis of the experimental variogram, a suitable model (e.g.spherical, exponential) is then fitted, usually by weighted least squares, and the parameters (e.g.range, nugget and sill) are then used in the kriging procedure.

Inverse distance weighting (IDW)
In interpolation with IDW method, a weight is attributed to the point to be measured.The amount of this weight is dependent on the distance of the point to another unknown point.These weights are controlled on the bases of power of ten.With increase of power of ten, the effect of the points that are farther apart diminishes.Lesser power distributes the weights more uniformly between neighboring points.We should keep in mind that in this method the distance between the points count, so the points of equal distance have equal weights (Burrough and McDonnell, 1998).In this method the weight factor was calculated with (Eq.2): Where: λ i is the weight of point, D i denotes the distance between point i and the unknown point and α is the power ten of weight.

Cokriging
The "co-regionalization" (expressed as correlation) between two variables, can be exploited to advantage for estimation purposes by the cokriging technique.In this sense, the advantages of cokriging were realized through reductions in costs or sampling effort.The cross-semivariogram was used to quantify cross-spatial auto-covariance between the original variable and the covariate (Webster and Oliver, 2001;Stefanoni and Hernandez, 2006).The cross-semivariance was computed by (Eq.3): Where: λ uv (h) is cross-semivariance between u and v variable, Z u (x) is primary variable and Z v (x) is the secondary variable.

Comparison between the different methods
Finally, the criterion of root mean square error (RSME) was applied to evaluate model performances in cross-validation mode.The smallest RMSE indicates the most accurate predictions.The RMSE was derived according to (Eq.4): Where: Z(xi) is observed value at point xi, Z * (xi) is predicted value at is predicted value at point xi and N is number of samples.

Relationship beteween environmental factores and vegetation cover
There were differences in vegetation cover in different geopedologic map units.The highest vegetation cover was observed in Hill with glacis relief-type, moderate slope percent, dominant class altitude of 1200-1400 m and deep soil.Also, the lowest cover percentage of vegetation occurred in Valley.There was a significant difference in percentage canopy cover of vegetation between glacis (Hi211) and other landforms except some parts of Mountain with soil cover (Mo221).Also significant differences were found in the vegetation cover of Valley, Mountain and Hill (Tab.2).There was a low positive correlation between altitude and vegetation cover (r=0.22,p=0.02).There was a significant difference in percentage canopy cover of vegetation between where altitude was less than 1200 m and where altitude was between 1200-1400 m (Kruskal-Wallis test P<0.01).There was a significant difference in vegetation cover among slopes more than 65° and slope 0-10° and 10-35° (Tab.3).There was no significant difference in vegetation cover between different aspects.As high tem-tween soil map unit A, D, E and C (α=0.05, p<0.03).Soil map unit influences the distribution of plant species and each species assessed had a different canopy cover in a different soil map unit.Highest canopy cover of Cymbopogon and Platychaete occurred in the E soil map unit, whereas highest canopy cover of Asteragalus and Artemisia occurred in the D soil map unit.
There was no significant difference between EC and percentage cover of vegetation.There was no difference for vegetation cover between different soil textures.Browsing damage is a factor influencing percentage of vegetation cover therefore the relationship between browsing damage and vegetation cover was assessed in different geopolitical map units, slope classes and altitude.There was a negative relationship between percentage cover and browsing damage in Hill, Piedmont and Valley landscapes.Lowest vegetation cover occurred in Valley where it had the highest browsing damage.There was no clear relationship between vegetation cover and browsing damage in different altitude classes.There was a negative relationship between browsing and cover on different slopes meaning when slope steepness increased vegetation cover increased and browsing damage decreased.
Generally, canopy cover percentage in the Mountain landscape is higher than other landscapes because the study area is located in an arid zone, where the important factor influencing vegetation is moisture.In the "Mountain" with calcareous lithology, moisture is available more than other landscapes.This result was the same as those found under desert conditions, where cliffs and rocky outcrops often enjoy conditions more favorable to plant life than other stony habitats (Zohary, 1973).This was because in rock crevices fine soil material was accumulated and most rain water that runs into these clefts is well preserved and protected against evaporation.In addition, rocks are favorable habitat for a number of shade demanding plants.Rocks are also inhabited by true lithophytes, whose roots are able to break the solid rock into pieces.Many studies confirm that there is significant relationship between altitude and vegetation because of the high correlation between precipitation and elevation which has been observed in semiarid peratures and low soil depth on the hillsides has caused low top soil moisture.
The properties of soil vary from place to place.Natural soil bodies are the result of climate and living organisms acting on parent soil material, with topography or local relief exerting a modifying influence over the time required for soil-forming processes to act.For the most part, soils were the same wherever all elements of the five factors were the same (Soil survey manual).Based on soil analyses for texture, EC, pH, there were no significant variation in the area, except in one of the mountainous landforms where the parent material was salt dominated.Also there were some variations in the soil depth with changing landscape, in mountains and hills the soil depth was shallow while in piedmonts it was deeper.There were differences in vegetation cover percentage in different soil map units and some of them were significant (Tab.4).Lowest canopy cover occurs in B soil map unit that consist of Typic Xerofluvents, Loamy Skeletal and Typic Xerorthents, Loamy Skeletal.There was significant difference in canopy cover between soil map unit A, B, D and E (α=0.01, p<0.003) and be- The first step for making use of kriging and cokriging methods in the present study was to investigate the presence of spatial structure among the available data by means of variogram analysis.This was achieved using those data which were normalized.Subsequently, variograms for the kriging method were calculated (Fig. 2).After that, the best model for fitting to the experimental variogram was selected based on the smallest Residual Sums of Squares (RSS) values (Tab.6).Therefore, the spherical model was selected as the most suitable model for estimation of percentage plant cover, while the gaussian model was the most suitable for evaluating of soil pH and EC.
The ratio of nugget variance to sill expressed in percentages (C O /C O +C) was regarded as a criterion for classifying the spatial dependence of soil parameters.Since this ratio for all of the soil properties were less than 25 percent, these parameters showed strong spatial dependence.Also, the range effect for plant cover, EC and pH were about 34, 67 and 51 km, respectively (Tab.7).
The first step in cokriging method was to compute of cross-variograms.The cross-variograms were modeled the same as variograms so that the same restricted set of functions were available for them.Cross-variograms were modeled to predict the spatial relations between two variables in the cokriging method.Typically the major aim was estimating of just one variable in this technique.However, more variables which were regarded as auxiliary variable were estimated in the present study.In general, cokriging reduces the estimation of variance but, this reduction strongly depends on the number of samples taken.Moreover, in the present study the factor of land elevation from a Digital Elevation Model (DEM) was used as an auxiliary variable to develop the cross-variograms (Fig. 3).
After modeling the variograms, three different techniques: kriging, cokriging and IDW were used to predict the spatial distribution of soil characteristics and percent plant cover.The RMSE was used to evaluate the three geostatistical techniques.Results showed that the kriging and cokriging methods were expected to be superior to IDW method for estimation of soil properties (Tab.8).Results also indicated that for all of studied soil and plant parameters cokriging method had the greatest accuracy.Finally, maps of spatial distribution of studied soil properties and landscapes (Smith et al., 1990).As slope increases vegetation cover percentage increases because human activity, such as cutting and browsing damage, is greater on shallower than steeper slopes.

Spatial prediction and mapping
Statistical data which were related to studied parameters were calculated using raw and normalized values (Tab.5).Normalization was applied according to the level of skewness for pH and percentage plant cover, using the square root method.The data for soil EC which had a skew value greater than one were normalized using the logarithm method.When distribution of data is complicated and fitting them by common statistical distributions is difficult, inter-percent plant cover were prepared using the best interpolation method in GIS environment (Fig. 4).
Many variables exhibit a non-normal distribution of measured values and therefore don't initially satisfy the basic assumption in geostatistics of statistical normality.This restriction is eliminated by applying a data transform to the sample values that make them more amenable to anal- These properties have a large and prominent role for increasing the estimation accuracy.
According to the geostatistical results the kriging and cokriging methods were better than the IDW method for predicting the spatial distribution of soil and plant properties.These findings are similar to the works done by Mo-polation by kriging method is the best option.After drawing a variogram and fitting the most suitable model, different parameters related to them was extracted.Having a stability and firmness in spatial structure for all the studied soil properties shows strong spatial connection and high accuracy among the models extracted from these fittings.rameters.Furthermore, the results obtained from evaluation of different methods showed that the cokriging method increased prediction accuracy in estimation of soil hammadi (2000) and Sokoti et al. (2006).They similarly reported that the geostatistical methods had considerably more accuracy than the IDW method for all studied pa-  types of lands.Therefore by means of the most developed techniques such as geostatistical methods the maps of soil and crop properties can be provide using the lowest number of data.Hence, besides using interpolation methods, the data derived from satellite information can be used for cokriging method as an auxiliary variable in next studies.In addition, most likely this methodology could be used for mapping more permanent physical soil properties, even when aerial colour photographs were available from an archive of standard surveys (Lopez-Granados et al., 2005).geostatistical methods for prediction of soil properties.Moreover, Sokoti et al. (2006) reported that the cokriging method was a suitable method for estimating the percentage of sand, clay and Saturation Percentage (SP).In addition, the great efficiency of the geostatistical methods is supported by other researchers.Wei et al. (2006) evaluated the soil organic matter distribution in north east of China.They found that the kriging method could predict organic matter distribution with a high accuracy.Xiaopeng and Lingqing (2008) used the ordinary kriging method for estimating the spatial distribution of mercury pollution in the city of Baoji in Chaine.They similarly confirmed the high accuracy of the kriging method.Shao et al. (2006) also used some geostatistical techniques and kriging method to determine the spatial distribution of soil nutrients in a vegetable production area of Hebei province in China.The results revealed that all the geostatistical techniques increased the accuracy of predictions.Shi et al. (2005) provided the salinity distribution maps in a coastal saline field in China using geostatistical methods.They introduced the cokriging method with auxiliary parameters of Na + and Ca +2 + Mg +2 , respectively as the most suitable method for predicting the level of SAR and EC.In the same way, Hosseini et al. (1994) found that the cokriging method had the most accuracy for estimation of SAR level.
Maps from these kriged estimates showed that a combination of geostatistical techniques and digital data in a GIS environment were accurate enough to improve the identification of soil management zones, which is the first step for site-specific management.Moreover for improving predictions of sparse information from soil surveys it could be very useful to use more intensive, densely sampled but less expensive ancillary data, even if they are from an archive of standard surveys with coordinates not exactly at the same locations (Lopez-Granados et al., 2005).

Conclusions
Based on this study geopedological unit was a good indicator for determining vegetation type.So that the final result showed that these boundaries were very adequate to determine the vegetation type map.This shows how soil and vegetation can be integrated.The soils of the study area were very young without subsurface diagnostic horizons, so all the soils were classified as Entisols.These soils were located on the hill, mountains foot slopes, and young alluvial culluvial fans with very high gravel (more than 50 percent).This reduces the ability of the soil as a growing environment for dense vegetation.Hence, there was no significant relationship between vegetation cover and vegetation index in the study area.The sustainable agriculture will take place just with classification of the lands according to their abilities and proportions for different types of land use.For achievement to this issue, the map of the soil and crop properties in this study can be one of the basic resources for giving valuable information about various

Fig. 1 .
Fig. 1.Location of study area in southern Iran and sampling points ysis and estimation.The most useful data transform is the log-transform.Since natural log values can be back transformed to real values, we can use a semivariogram model derived from the transformed sample values to predict the spatial variation of logarithmic values.

Fig. 2 .
Fig. 2. Variograms related to different factors in kriging method

Fig. 4 .
Fig. 4. Maps of soil and plant properties