Spatial Variability Analysis of Soil Field Capacity at Lalgudi Block in Tamil Nadu using Geostatistical Models

Soil properties, by nature vary extensively over space and time (Greenholtz  et al., 1988). Spatial variability of soil properties has a great effect on the site-specific management of soil properties. Field capacity is amount of water held in soil after the gravitational water had drained from the saturated soil. It is a soil parameter that is widely used in soil hydrology, land management and irrigation and drainage engineering. The spatial variation of field capacity throughout the irrigation regions must be considered for proper scheduling of irrigation. Geostatistics accounts for identifying the best interpolation method to acquire the spatial map of soil parameters.
       
Geostatistical tools were used to study the spatial variation of soil moisture, texture and yield of two rainfed crop namely vetch-oat and durum wheat (Munoz-Pardo  et al., 1990). Brooker  et al. (1995) presented a geostatistical analysis on depth of topsoil and root zone readily available water at a Vineyard near Waikerie in South Australia and prepared spatial maps which gave better understanding of the soil in that region. Santra  et al. (2008) highlighted that spatial variation of soil properties like field capacity (FC) and permanent wilting point (PWP) is important in precision farming and environmental modelling. Anane  et al. (2012) established a GIS model integrated with fuzzy approach to map and rank suitable sites for irrigation with treated wastewater. De Benedetto  et al. (2013) applied geostatistical techniques to estimate the soil water content and explored the capability of these methods. De Paz  et al. (2015) estimated the spatial variation of maximum irrigation rate using the geostatistical techniques and presented the spatial maps that could serve irrigation managers for optimizing the irrigation rates for each field. Li  et al. (2017) examined the temporal Irrigation Water Productivity trend of cereal crops over the Hexi Corridor in Northwest China by employing descriptive analysis, trend analysis and change-point analysis and showed that area supported by unit of irrigation water use, fertilization and agricultural film had dominant impacts during the whole period. Shit  et al. (2016) measured the soil properties like pH, electrical conductivity, phosphorus, potassium and organic carbon. Random samples were taken from Medinippur Sadar block of Paschim Medinipur district in West Bengal (India). Laekemariam  et al. (2018) assessed the spatial variability of soil fertility and mapped application of fertilizer types using geostatistics. Singha and Swain (2020) used exponential model generated by Ordinary Kriging to estimate the lentil land use suitability. Thakur  et al. (2021) used GPS and GIS Techniques to assess the soil fertility status to propose effective strategy on fertilizer use and cropping pattern. Geostatistics was used for generating soil textural maps for precision agriculture and land management in the research region (Rajalakshimi  et al. 2023). Spatial variability maps of various soil physical properties will be useful for site-specific farming, for example, variable rate irrigation (Rajalakshimi et al., 2023; Ramachandran et al., 2024).
               
With this contextual enlightenment, the main objective of this study is to identify the best geostatistical interpolation technique for spatial estimation of field capacity of soils in Lalgudi block of Tiruchirapalli district, Tamil Nadu, India.

Study area
 
The study was conducted in Lalgudi block of Tiruchirapalli district, Tamil Nadu, India (Fig 1). Lalgudi block is situated at 10^o52’27” N and 78^o48’57” E geo-coordinate and located 70 m above mean sea level. The total geographical area of the block is 20558 hectares. Lalgudi block has semi-arid climate with an average rainfall of 877 mm.


 
Soil sample collection and laboratory analysis
 
Soil samples were collected from different parts of Lalgudi block at a depth of 30 cm. Around 20 samples were collected within Lalgudi block at random locations as shown in Fig 2. The locations include paddy field, banana field, sugarcane field and forest area. Using Global Positioning System (GPS), the latitude and longitude was recorded for the places where soil samples were taken. Field capacity was determined using Pressure Plate Apparatus for soil samples collected at different locations (L1 to L20). The experiment was done twice to get concurrent results.

​
 
Geostatistical modelling
 
Exploratory spatial data analysis for soil parameters
 
Firstly, exploratory spatial data analysis aids in examining the data in different ways. It was done with Explore Data option in Geostatistical Analyst tool of ArcGIS. The exploration of data was done via. 1) Examine the distribution of the dataset; 2) identify the trends in the dataset; 3) Understand the spatial autocorrelation and directional influence of the dataset. The histogram tool was used to display the frequency of distribution of the dataset and summary statistics. The normal QQ plot tool was used to get the quantile-quantile (QQ) plots. It is a graph on which quantiles from two distributions are plotted relative to each other. The trend analysis helps in identifying the global trend in the input dataset. A unique feature of the Trend Analysis tool is that the values are then projected onto the x,z plane and the y,z plane as scatterplots. Polynomials were fitted to the scatter plots on the projected planes. If the curve through the projected points is flat, then no trend exists. Semivariogram modeling is a key step between spatial description and spatial prediction. The empirical semivariogram provides information on the spatial autocorrelation of datasets.
 
Kriging methods
 
The most common interpolation technique used is Kriging. It is an advanced geostatistical procedure that generates an estimated surface from a scattered set of points with z-values. In this study, ordinary kriging and disjunctive kriging was used for the interpolation of dataset. Ordinary kriging is the most wide spread method resulting in a smoothed surface and inexact interpolation. When the dataset follows a normal distribution, ordinary kriging is used. The standard model for ordinary kriging is given by:

 

                                                                                                                                                                                                                               
Where,
Z(s)= Variable of interest, decomposed into a deterministic trend.
µ(s)= A random, autocorrelated errors form e(s).
s = Location i.e. the spatial x- (longitude) and y- (latitude)  coordinates.

Disjunctive kriging is a non-linear type in which the original dataset is transformed using a series of addictive functions, typically Hermite polynomials. Therefore the standard model is altered and given as:
                                                                                                                                                                                                                                 

                                                                                                                                                      

 
Where,
F()= An arbitrary function of Z(s).
       
In general, disjunctive kriging tries to do more than ordinary kriging.
 
Experimental semivariogram
 
Let Z be an intrinsic random variable function and let Z(X), for i = 1, 2 …N, be a sampling of size N. Then the semivariogram of the random function is given by:
 

 

 
Where,
n(h) = Number of pairs at the lag distance h apart.
        This is called as experimental semivariogram.
 
Characteristics of semivariogram
 
A semivariogram is a graphical representation of semivariance and lag distance. The distance at which the semivariogram first flattens out is known as the range (a). Sample locations separated by distances closer than the range are spatially autocorrelated, whereas locations farther apart than the range are not. The value of semivariogram at the range on the y-axis is called the sill (C). The partial sill is the sill minus the nugget. Theoretically, at zero separation distance (lag = 0), the semivariogram value is 0. However, at an infinitesimally small separation distance, the semivariogram often exhibits a nugget effect, which is a value greater than 0. All the parameters were estimated using the Geostatistical Tool in the ArcGIS.
 
Dependency index (DI)
 
The spatial correlation of the parameters was estimated by calculating the Dependency Index. It is the ratio between Nugget and Sill. It is given by the following formula:

                                                                                                                                                               
  When DI is less than 0.25, then there is a strong correlation; DI is between 0.25 and 0.75, there is a moderate correlation and DI is greater than 0.75, there exists a weak correlation of dataset (Cambardella et al., 1994; Bo et al., 2003).
 
Semivariogram models
 
The semivariogram depicts the spatial autocorrelation of the measured sample points. Once each pair of locations was plotted, a model is fit through them. Geostatistical Analyst provides several semivariogram/covariance functions to model the empirical semivariogram among which the following six models were used in this study.
       
The circular semivariogram was estimated by the following model:

 

 
 
 A spherical semivariogram is line in shape near the origin. It was estimated by the following model:
 

 

 
 
The model is said to be transitive because it reaches a finite sill at a fine range. It is given as default model in some packages.
       
The exponential semivariogram was estimated by the following model:
 

 

                                                                                                                                                                                                                              

The sill is approached asymptotically. It is useful when there is a larger nugget and a slow rise to the sill.

The Gaussian semivariogram was estimated by the following model:

 

                                                                                                                                                                                                                          

 The model approaches the sill asymptotically. A graph model has a parabolic form near the origin. It provides more S-shaped curve.    
       
The pentaspherical semivariogram was estimated by the following model:
 

 
The Sine Hole effect semivariogram was estimated by the following model:
 

The model is not monotonic. It reaches a global maximum and then it continues a damped oscillation around the sill. Each model is designed to fit different types of phenomena more accurately.
 
Cross validation statistics
 
In order to study the interpolator performance, statistical formulas like mean error (ME), root-mean-square error (RMSE), average standard error (ASE), mean standard error (MSE) and root-mean-square standardized error (RMSSE) were used. The formulas are given below:

 

 
Generally, the best prediction model is the one that has the standardized mean nearest to zero, the smallest root-mean-square prediction error, the average standard error nearest to the root-mean-square prediction error and the standardized root mean square prediction error nearest to one. If RMSE is less than ASE and RMSSE is less than 1, then the model overestimates the given dataset.

Exploratory data analysis
 
The histogram of field capacity along with the mean, maximum, minimum, standard deviation, skewness, kurtosis and median of soil properties is shown in Fig 3. It was found that field capacity of twenty soil samples were non-normally distributed. It was confirmed by the large difference between mean and median of each dataset. For example, the mean of field capacity was 23.77 and the median was 22.77 percentage. Hence ordinary kriging could not be applied for the dataset without applying data transformation.

​
       
Therefore, log transformation was made to the given dataset and the histograms were estimated and shown in Fig 4. The mean and median of the log transformed dataset was close to each other. For example, mean of field capacity was 3.10 and median was 3.12. This indicated normal distribution of log transformed dataset of field capacity.

​
       
The Normal QQ plots without transformation of dataset and with transformation of dataset were obtained for field capacity. Fig 5a and 5b represents the QQ plots of field capacity without transformation and with log transformation respectively. The points were closer to the normal straight line in case of log-transformed dataset indicating the closeness of the dataset to be normally distributed when data transformation was done.

​
       
The trend exhibited by field capacity is shown in Fig 6 which helps in identifying the global trend in the dataset. The polynomial curve (blue line on the xz plane) in all the dataset indicates that there exists a trend exhibited by field capacity. Field capacity demonstrated a strong u-shape curve and it indicates that field capacity had strongest influence from the center of Lalgudi block towards all its boundaries.  The curve through projected points in the yz plane is almost flat in all the three dataset which indicated no trend in the yz plane for all dataset.

​
       
The semivariogram cloud provides information on spatial autocorrelation of dataset and look for outliers. The semivariogram cloud for field capacity is shown in Fig 7. Each red dot indicates empirical semivarigram value of two soil samples plotted against the distance of separation. When closely examined, as the distance increases, the pair of points of field capacity are closer together and it indicates that the points are alike.

​
       
The experimental semivariogram fitted with six models using ordinary and disjunctive kriging for field capacity is shown in Fig 8 and 9 respectively. From those semivariograms, the characteristics of semivariogram like nugget, sill and range and dependency index were estimated and presented in Table 1 and 2 respectively. The results showed that the field capacity dataset had moderate spatial dependency in case of both kriging methods and with six models (Table 1 and 2) except Sine-Hole effect model in disjunctive kriging. Sine-Hole effect model in disjunctive kriging exhibited a weak spatial dependency for field capacity dataset. The results indicated that theoretical models satisfactorily represent the spatial variability of field capacity in the study area.

​

​





       
For field capacity dataset, a maximum range of 10709 m was observed in Exponential and Sine-Hole Effect model in ordinary kriging. The higher range value indicates more continuity and smoother spatial variability of soil property.  However, the range was lesser in case of disjunctive kriging when compared with ordinary kriging.
       
The exploratory data analysis resulted in similar results as reported by to De Paz  et al. (2015) where, a non-normal distribution for field capacity was obtained and those variables were log-transformed to ensure normal distribution. The range indicates the distance in a field where measured properties are no longer spatially correlated. Measured properties of the samples at a distance less than the range become more alike with decreasing distances between them (Tabi and Ogunkunle, 2007). De Paz  et al. (2015) also reported that the range was between 3000 and 5200 m for field capacity. When comparing results from this study with De Paz  et al. (2015), it must be pointed out that all soil properties dataset exhibited higher range indicating a smoother spatial variability and also the number of samples collected was sufficient for interpolation within the block.
 
Cross validation statistics of semivariograms
 
The cross validation statistics is useful in predicting the best performing model for interpolation of data. The results of cross validation parameters like ME, RMSE, ASE, MSE and RMSSE is presented in Table 3 and 4. From the results, it was inferred that field capacity dataset was overestimated by the six models in both kriging methods. The best prediction method should have lesser RMSE value and hence the circular model in disjunctive kriging gave lesser RMSE value. The nugget effect, which represents random variation caused mainly by the undetectable experimental error and field variation within the minimum sampling space (Cerri  et al., 2004; Askin and Kizilkaya, 2006) is close to zero for the field capacity in this study. Gülser et al. (2016) also indicated that nugget values close to zero for the soil physical properties which revealed that all variances of the soil properties were reasonably well explained at the sampling distance used in that study by the lag distance.




 
Spatial interpolation of field capacity
 
The spatial distribution maps of the field capacity is shown in Fig 10 and 11. The spatial interpolation of field capacity using Ordinary Kriging gave field capacity which varied from 13.0 % to 28 % for different models and disjunctive kriging gave field capacity which varied from 14.0% to 30.0%. Overall, in all the spatial distribution maps, the field capacity is low in the northern part of the block and low at the central and western part. The eastern part of the Lalgudi block is having high field capacity as indicated in the spatial distribution maps.

​

​
       
The results showed that the field capacity had moderate spatial dependency in case of both kriging methods and indicated that theoretical models satisfactorily represent the spatial variability of field capacity in the study area. But Gülser  et al. (2016) reported strong spatial dependency in the soil physical properties.
               
Santra  et al. (2008) showed that spatial prediction of basic soil properties using semivariogram parameters is better than assuming mean of observed value for any unsampled location when cross validation of the kriged map was done. Kriged maps illustrated positional similarity between the field capacity along the Lagudi block. Santra  et al. (2008) also repored that evaluation of spatial maps of field capacity showed reasonable accuracy for farm-level or regional-scale application.

Geostatistics was applied for identifying the best interpolation method to acquire the spatial map of field capacity. Initial data exploration resulted in non-normal distribution of field capacity and hence log-transformation was done before kriging was used. The cross validation statistics showed that circular model with disjunctive kriging of field capacity was the best model for interpolation with lesser root mean square error and also had moderate spatial dependency. It overestimates the field capacity when interpolation was done. Consequently, the final maps and calculated results can also be used in the decision processes for land and water managements and soil conservation practices by authorities, as well as by farmers for irrigated fields in this semi-arid areas. These maps play key roles in crop selection for different blocks of a farm and in scheduling irrigation of crops in a field.

Disclaimers
 
The views and conclusions expressed in this article are solely those of the authors and do not necessarily represent the views of their affiliated institutions. The authors are responsible for the accuracy and completeness of the information provided, but do not accept any liability for any direct or indirect losses resulting from the use of this content.

Informed consent
 
All animal procedures for experiments were approved by the Committee of Experimental Animal care and handling techniques were approved by the University of Animal Care Committee.

The authors declare that there are no conflicts of interest regarding the publication of this article. No funding or sponsorship influenced the design of the study, data collection, analysis, decision to publish or preparation of the manuscript.