Soils are inherently heterogeneous in nature, diverse and dynamic system  and its properties change in time and space continuously  . Heterogeneity in soil properties with depth and across landscapes can be accounted for by several interacting factors that operate with different intensities and at different scales and acting simultaneously  . Estimating spatial variability of soil properties is important for evaluating environment  and provides the factors and processes controlling potential in agriculture production  . It is also an important determinant of efficiency of farm inputs and yield as well as crop management and design and effectiveness of field research trials.
The availability of soil nutrients for plant growth and yield production is а function of different parameters, including soil pH, soil organic matter and texture, and soil biological activities  . Hence, determination of such parameters is important for evaluating nutrient behavior in the soil and for suggesting appropriate methods of enhancing nutrient availability to plant.
The important way to gather knowledge in this respect is to prepare maps through spatial interpolation of point based measurements of soil properties using geostatistics. There have been growing interests in the study of spatial variation of soil properties using geostatistics since 1970s, as geostatistics techniques were well developed and successful in characterizing the spatial variations of soil properties  . While many studies have been carried out at a small-scale  , relatively few have been done at large-scale  .
The current study was undertaken in the Research farm of SKUAST-K, Shalimar for analyzing the spatial variability of soil properties and for identification of nutrient deficiency zones for site specific nutrient management.
2. Materials and Methods
2.1. Study Area
The present investigation was carried out in a Research farm of SKUAST-K, Shalimar (34˚8'42" and 34˚9'3"N latitudes and 74˚39'5" and 74˚53'5.6"E longitude) Srinagar (Figure 1). It has 1615 m average altitude above sea level and covers an area of 23.8 ha. The climate is temperate and characterized by mild summers and chilling winters having normal annual maximum temperature of 19.53˚C and minimum of 6.80˚C with normal annual rainfall of 786.2 mm. Dominant vegetation are cereals (wheat, rice, maize, oats), vegetables, fruits and floriculture.
2.2. Soil Sampling and Analysis
A total of seventy seven (77) samples were selected in a systematic grid design using Arc GIS. Each grid was specified at a fixed distance of 50 × 50 m2 grid from 0 - 22.5 cm depth. Samples were thoroughly mixed and ground to pass through 2 mm sieve, then stored in plastic bags prior to chemical analysis. Soil pH was determined in 1:2.5 soil: water suspension with digital glass electrode pH
Figure 1. Georeferenced sampling sites of research farm of SKUAST-K, Shalimar.
meter  . The organic matter (OM) was determined using the K2Cr2O7 titration method  . The Available nitrogen (AvN) was determined by the Alkaline Potassium Permanganate method  . The Olsen method was used to determine available phosphorus (AvP) using a molybdate reaction for colorimetric detection  . The neutral 1N ammonium acetate extraction method was used to determine exchangeable/available potassium (AvK)  . Available sulphur was determined by following the turbidimetric method of  .
3. Statistical Analysis
3.1. Exploratory Statistical Analysis
Statistical parameters which are generally accepted as indicators of the central tendency and spread of the data, were analyzed. These include description of mean, minimum and maximum values, standard deviation and coefficient of variation. To decide whether or not the data followed the normal frequency distribution, the coefficients of skewness and kurtosis were examined  . These statistical parameters were calculated using SPSS 20.0 release software  . The coefficient of variation (CV) was mainly used to assess the variability of the different data sets. Exploratory data analyses for normality tests were conducted. Normality tests were conducted using Q-Q plots   Non-normal data were transformed to stabilize the variance. Then normality tests were recalculated using the transformed data, as asymmetry in the distribution of data has an important effect on the geostatistical analysis  .
3.2. Geostatistical Analysis
Spatial analysis was carried out by the use of geostatistical method (Arc GIS) and mapping software (Surfer). Spatial variability in soil fertility parameters were calculated for 0 to 25 cm depth. Firstly variograms were applied to measure the spatial variability of sampled locations, which also provides the parameters that are necessary for interpolation of unsampled areas
Variograms and kriging interpolation were performed in ArcGIS 10.2. The formula applied to the variogram  is:
where γ(h) is the experimental semivariogram value at a distance interval h, m(h) is number of sample value pairs within the distance interval h, Z(Xi), Z(Xi + h) are sample values at two points separated by the distance h. Several semivariogram functions were evaluated to choose the best fit with the data. Spherical, Exponential or Gaussian models were fitted to the empirical semivariograms. The stationary models, i.e., Gaussian (Equation (2)), Exponential (Equation (3)) and Spherical model (Equation (4)) that fitted to experimental semivariograms were defined in the following equations  :
where C0 is the nugget, C1 is the partial sill, and a is the range of spatial dependence to reach the sill (C0 + C1). The semivariance generally increases with sample separation distance before reaching an asymptote а (the range value). Samples separated by distances greater than range value are considered to be spatially independent where as, within the range, samples show greater similarity when they are nearer to each other  . Variance that exists at а scale smaller than the field sampling is found at zero lag distance and is known as the nugget variance (C0)  .
The sill represented the amount of variation defined by the spatial correlation structure and it is the value of the semi-variogram at which the model first levels out (given as partial sill plus the nugget  . Partial sill (C) is the lag distance between measurements at which one value for а variable does not influence neighbouring values. The partial sill (C) is the variance caused by factors such parent material variability, and vegetation and topographic differences. The nugget-to sill ratio (SH) designates the degree of spatial heterogeneity arising from random components to that of the total spatial heterogeneity.
The nugget/sill ratio was used as a criterion for classifying the spatial dependence of soil properties. The variable has strong spatial dependence if the ratio is less than 25%; between 25% and 75%, the variable has moderate spatial dependence; and the variable shows only weak spatial dependence if the ratio is greater than 75%  . A value close to 0% indicates that the variable has strong spatial auto-relationship while that close to 100% indicates spatial heterogeneity is dominated by randomness, or nugget effect  .
As the semivariogram models of the soil data were evaluated, they were used in the development of maps by ordinary kriging interpolation  . The cross validation is applied to evaluate and compare the performance of different interpolation methods through mean square error (MSE), average standard error (ASE), root-mean-square error (RMSE) and the standardized root mean square error (RMSSE)  . For best fitted model, there must be minimum error   .
4. Results and Discussion
4.1. Descriptive Statistics
The descriptive statistics of the soil fertility parameters are given in Table 1. The coefficient of variation values (CV) was used to interpret the variability in soil properties. The criteria proposed by Gomes and Garcia  proposed was used to classify the parameters into low (<10%), medium (10% - 20%); high (20% - 30%) and very high (>30%) variabilities. Accordingly, present study indicates low to very high variability of soil fertility parameters within the fields. The greatest and the least CVs for soil parameters were obtained for soil available phosphorus (AvP) (56.87%) and pH (7.06%), respectively, indicating very high variability for AvP relative to the other soil parameters. A high CV is the first indicator of data heterogeneity  . Generally, pH and OC are considered to be
Table 1. Univariate statistical analysis for soil parameters.
SD = standard deviation, CV = coefficient of variation, OM = organic matter (%), N = Available nitrogen (kg∙ha−1), P = Available phosphorus (kg∙ha−1), K = Available potassium (kg∙ha−1), S = Available sulphur (kg∙ha−1).
stable soil parameters  . However, at the study area high variability (OC) observed for OC could be ascribed to pedogenic processes influenced by the micro-topographical variations operating over different periods of time   . Similar results were also figured out by Aishah et al.  who found CV value of 4% for pH and Tagore et al.  reported least variability (CV = 6.37%) for pH among all the analyzed soil parameters.
The normal distribution of data was examined by Quantile-Quantile (QQ) plot. The quantile-quantile plot (QQ plot) is a simple graphical method for comparing two sets of sample quantiles  . In this study, the normal Q-Q plots for the raw data were produced (Figure 2) and it was found that only OC and S followed a straight diagonal line The underlying reason for soil elements being distributed normally or non-normally may be associated with differences in management practices, land use, vegetation cover, and topographic, and topographic effects  .
Descriptive statistics in this study indicates moderate to high skewness. The value of skewness varies from −0.01 to 1.42 depicting moderate to high skewness (Table 1). Highly skewed parameters indicate that these properties have a local distribution; that is, high values were found for these properties at some points, but most values were low  . These high soil test values may not always be an outlier but a form of natural or management induced variation  .
4.2. Analysis of Spatial Dependence of Soil Fertility Parameters
For geostatistical analyses of soil parameters, an appropriate model was chosen. Best suited models for various parameters are presented in Table 2. Exponential models were fitted to the experimental semivariograms for the OC, pH, AvP, AvK and S while as only N was best suited to the Gaussian model. Tagore et al.  reported that Exponential model fits well with experimental semi-variogram of pH, EC, OC, available N, P, K, S and Zn. However, Reza et al.  while working on alluvial soils of India describes Spherical model to be the best fit for N, P and Zn contents. Moreover the findings are consistent with the researches of Some’e et al.  and Mahmoudabadi et al.  .
When the distribution of soil properties is moderately or strongly spatially correlated, the average extent of these patches is given by the range of the semivаriogram  . Comparing range values, longer range value was observed for N (763.20 m) and lowest for OM (99.60 m) (Table 2 and Figure 3). A larger range indicates that observed values of the soil variable are influenced by other values of this variable over greater distances than soil variables which have smaller ranges  . The soil sampling distance in the range of 50 × 50 m2 in this study was close with models range value of all the parameters. So, the collected samples can be valid and applicable in a two- or three-fold larger area. According to Kerry and Oliver  and Fu et al.  , the sampling interval should be less than half the semivariogram range. Thus it can be an effective criterion for the evaluation of sampling design and the mapping of soil properties.
(a) (b) (c) (d)
(e) (f) (g) (h)
Figure 2. Normal P-P plot for different soil parameters.
Variation at microscales smaller than the sampling distances will appear as a part of the nugget effect   . The value of nugget varied widely for soil properties. It was highest for sulphur (S) and the lowest for soil pH. Lower values of nugget effect (C0) indicate low errors in measurements. The high nugget
Table 2. Values of model parameters used to find the best semivariogram.
C0 = nugget effect; C1 = partial sill; C0 + C1 = sill; degree of spatial dependence (DSD) = C0 /(C0 + C1) DSD; strong DSD (<25%); moderate DSD (>25 to <75%); weak DSD (>75%). SD: Spatial dependence; MSE: Mean square error; ASE: Average standard error; RMS: Root-mean-square error; MSE: Mean standard error; RMSSE: Root-mean-square standardized error. alog transformed.
(a) (b) (c) (d)
Figure 3. Experimental and fitted semivariogram models for soil parameters: (a) OM; (b) N; (c) P; (d) K; (e) S; (f) pH.
of macronutrients was probably because of high soil heterogeneity resulting in large spatial variability of these nutrients.
The spatial dependence can indicate the level of similarity or disturbance of the soil condition  . The spatial dependency of the data was assessed from the ratio of nugget and sill (Figure 3). Cаmbаrdellа et al. (24) defined this ratio of <25 for strong, 25 to 75 for moderate, and >75 as weak spatial dependence. According to this classification, OM, K, S and pH showed а strong spatial dependence; and N and P exhibited weak degree of spatial dependence (Table 2).
(a) (b) (c) (d)
Figure 4. Spatial variability map of different soil properties (a) OM (b) N (c) P (d) K (e) S (f) pH.
Based on the results of the present study we may conclude that moderate and weak spatial dependence of soil fertility parameters can be usually attributed to soil and crop management practices  . These results are in line with the observations reported by Vasu et al.  .
Cross-validation was used to estimate which of the semivariogram models could give the most accurate predictions of the unknown values of the study area. It was shown that the error terms ME and MSE were close to zero. Subsequently, with implementing these best fit theoretical models and corresponding semivariogram parameters, spatial variability maps of soil properties were created using the ordinary Krigging (Figure 4). These maps will be potentially helpful for researchers for precision agriculture and site-specific nutrient management.
This study demonstrated that the classical statistics of the soil elements indicated a coefficient of variation up to 56.87%, which could not support to identify the sources of variability. This indicates that the classical statistical techniques were utilized to identify an overall variability of soil elements but lacked the necessary techniques to identify the kind of systematic spatial variability at farm scale. However, the geostatistical techniques offer alternative methods over the classical statistics for estimating the parameters spatial dependence and variability in the farm. According to the results, the semivariogram analyses show the presence of a moderate to weak spatial dependence of the selected soil properties within the study area. The Krigged maps of soil parameters can help the researchers to become familiar with the characteristics related to the analyzed soil properties and accordingly can plan appropriate agricultural strategies, including fertilization. Such analyses can save time and expenses, while being statistically of great precision and usability.