Digital soil mapping in the Bara district of Nepal using kriging tool in ArcGIS

Abstract

Digital soil mapping has been widely used to develop statistical models of the relationships between environmental variables and soil attributes. This study aimed at determining and mapping the spatial distribution of the variability in soil chemical properties of the agricultural floodplain lands of the Bara district in Nepal. The study was carried out in 23 Village Development Committees with 12,516 ha total area, in the southern part of the Bara district. A total of 109 surface soil samples (0 to 15 cm depth) were collected and analyzed for pH, organic matter (OM), nitrogen (N), phosphorus (P, expressed as P2O5), potassium (K, expressed as K2O), zinc (Zn), and boron (B) status. Descriptive statistics showed that most of the measured soil chemical variables (other than pH and P2O5) were skewed and non-normally distributed and logarithmic transformation was then applied. A geostatistical tool, kriging, was used in ArcGIS to interpolate measured values for those variables and several digital map layers were developed based on each soil chemical property. Geostatistical interpolation identified a moderate spatial variability for pH, OM, N, P2O5, and a weak spatial variability for K2O, Zn, and B, depending upon the use of amendments, fertilizing methods, and tillage, along with the inherent characteristics of each variable. Exponential (pH, OM, N, and Zn), Spherical (K2O and B), and Gaussian (P2O5) models were fitted to the semivariograms of the soil variables. These maps allow farmers to assess existing farm soils, thus allowing them to make easier and more efficient management decisions and maintain the sustainability of productivity.

About the Authors:

Dinesh Panday

Roles Conceptualization, Formal analysis, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliation: Department of Agronomy and Horticulture, University of Nebraska-Lincoln, Lincoln, Nebraska, United States of America

ORCID logo http://orcid.org/0000-0001-8452-3797

Bijesh Maharjan

Roles Methodology, Validation, Writing – review & editing

Affiliation: Department of Agronomy and Horticulture, University of Nebraska-Lincoln, Lincoln, Nebraska, United States of America

Devraj Chalise

Roles Writing – review & editing

Affiliation: Nepal Agricultural Research Council, Lalitpur, Nepal

Ram Kumar Shrestha

Roles Data curation, Resources, Writing – review & editing

Affiliation: Institute of Agriculture and Animal Science, Lamjung Campus, Lamjung, Nepal

Bikesh Twanabasu

Roles Data curation, Methodology, Resources, Software, Visualization, Writing – review & editing

Affiliations Hexa International Pvt. Ltd., Lalitpur, Nepal, Institute for Geoinformatics, Westfalische Wilhelms-Universitat Munster, Munster, Germany

Introduction

Applications of pedometric mapping, also called predictive mapping, i.e., the spatial prediction of soil variables at unobserved locations using statistical inference, have become increasingly important since their initial development in the early 1800s. The utility of such maps was due to the introduction of geostatistics, allowing researchers to accurately interpolate spatial patterns of soil properties [1]. One of the current versions of pedometric mapping, digital soil mapping (DSM), involves the creation and population of spatial soil information systems using field and laboratory observational methods coupled with spatial and non-spatial soil inference systems [2–5]. Field sampling is used to determine the spatial distribution of soil properties, and these surface grids point data are then interpolated to estimate soil properties in areas not sampled [6]. In contrast, existing conventional soil survey methods are relatively slow and expensive; therein, soil databases are neither exhaustive enough nor precise enough to promote an extensive and credible use of soil information within the spatial data [7].

Several statistical models can be used in DSM to develop a relationship between soil properties and environmental variables (often related to soil forming factors such as terrain attributes—altitude, aspect, and slope) rather than from soil observations alone [8], and McBratney et al. [9] thorough reviewed the various models used. Soil nutrients are one of the most important properties governing soil quality, and hence have a significant impact on the variability of soil productivity and crop production. The spatial variability of soil properties can be mapped using an interpolation technique [10]. Many spatial interpolation methods have been developed and several terms have been used to distinguish them, including “deterministic” and “stochastic” [11]. Deterministic interpolation methods such as thiessen, density estimation, inverse-distance-weighted, and splines, provide no assessment of errors, whereas stochastic interpolation and kriging methods do provide prediction error assessments.

Kriging is a geostatic interpolation technique that has proven sufficiently robust for estimating values at non-sampled locations based on sampled data. It provides the best linear unbiased estimates and information on the distribution of the estimation error and shows strong statistical advantages [12]. The use of the geostatistical interpolation technique also reduces the costs of field sampling and laboratory analysis, provided that a given set of soil samples sufficiently represents the study area [13]. However, the reliability of spatial variability maps depends on the adequate sampling data and the accuracy of the spatial interpolation [14].

There is a significantly increased trend in the use of DSM mainly due to recent advances in technology related to quantitative methodologies and geographic information systems. For example, spatial variability of organic matter (OM), pH, and potassium (K) were mapped using kriging by Lopez-Granados et al. [15] in a 40 ha field located in southern Spain. Santos Francés et al. used the kriging interpolation method for the production of spatial distribution of the heavy metal contents in the soils of northern Spain [16] and northern Peru [17], respectively. Balkovič et al. [18] reported that the DSM model represents a complete alternative to classical soil mapping at very fine scales on erosion affecting 37 ha of arable land in Slovakia, even when soil profile descriptions were collected merely by field estimation methods. In northwestern Australia, a DSM soil carbon map at the farm scale was developed from a total of 127 soil sampling locations in an area of 2300 ha [19]. Similarly, Zhang et al. [20] produced spatial variability maps of nitrogen (N), phosphorus (P) and K in winter wheat and summer maize in northeast China. Recently Zhu et al. [21] used an alternative DSM method, individual predictive soil mapping (iPSM), to map OM content in the topsoil layer of an area of 6000 ha in China and observed that iPSM is an effective alternative when existing soil samples are limited in their ability to fully represent the entire study area.

Despite the successful application of DSM in regions around the world, no single study has examined the use of DSM to represent soil nutrient variability in any part of Nepal. The commonly practiced soil fertility assessment is based on a random soil sampling protocol to obtain an average fertility value for a farmer’s field [22]. It ignores spatial variability, or those soil testing results that do not provide randomness of variations from one place to another. Consequently, some parts of the field may receive surplus fertilizer while others may lack nutrients and experience the undesired levels of productivity. The objective of this research was to determine and map the spatial distribution of variability in soil chemical properties of agricultural floodplain lands in the southern part of the Bara district of Nepal. The country’s economy relies heavily on agriculture and any breakthrough in soil mapping would immensely benefit farmers. Information on spatial variability of soil nutrients is also essential for sustainable management of soil fertility.

Materials and methods

Study area

Nepal is located in the south of Asia bordering neighbors of India and China and covers an area of approximately 147,181 km2. It has been divided broadly into three geographic regions: Himalayan, Hilly, and Terai. The study was conducted in the southern part of the Bara district, which falls under Terai region, and included the 23 Village Development Committees (VDCs) and covers 12,516 ha of land as shown in Fig 1. The topographic variation of the study area ranges from 80 to 95 m.

[Figure omitted. See PDF.]

Fig 1. Study area in southern part of Bara district, Nepal which includes 23 Village Development Committees (VDCs).

A large part of the study area lies in Terai region where the topographic variation ranges from 80 to 95 m, and climate is subtropical monsoon. A total of 109 soil samples were taken from the depth of 0–15 cm (topsoil layer) for determination of pH, OM, N, P2O5, K2O, Zn and B status on it.

https://doi.org/10.1371/journal.pone.0206350.g001

There are four seasons in Nepal: pre-monsoon (March to May), monsoon (June to September), post-monsoon (October to November) and winter (December to February). Monsoons are the Nepal’s main source of precipitation, accounting for 85% of the country’s total annual rainfall of 1800 mm, with the remaining 15% occurring in winter [23–24]. During a monsoon, all of the rivers are in spate, with bank-full discharges that cause flooding and inundation in several parts of the Terai region [25]. The study area becomes hottest (37 to 42°C) during the monsoon months compared to the country’s average temperature (28°C) due to seasonal changes and low altitude.

Rice (Oryza sativa L.) is the principal staple food of Nepal, accounting for about 67% of total cereal consumption. Most of the food crops for the entire country are grown in the Terai region, the granary of Nepal. The terai is generally made up of flat terrain with a hot, humid climate. About 80% of the land in this area is occupied by farmlands. Rice -wheat (Triticum aestivum L.)-fallow is the dominant cropping system in the study area followed by rice-wheat/lentils (Lens culinaris)-fallow, rice-wheat-maize (Zea mays L.), and sugarcane (Sacharum officinarum).

The soil association of the study area is developed by the changing river morphology. The soils have predominantly evolved from alluvial deposits and are dominated by sandy loam and silty clay, although clay loam and loamy sand are also present at considerable levels. It was observed that the majority of the area is occupied by land system unit 2b (deep alluvium: <0.5 degree slope, flat, imperfect drainage, sandy loam to silty clay, Aeric, Haplaquepts, Typic, and Fluventic) followed by 2a (deep alluvium: < 0.5 degree slope, depression, poor drainage, loam to silty clay, Aeric, Haplaquepts, and Typic), 3a (deep alluvium: < 1 degree slope, gently undulating, moderate drainage, sandy loam to silty clay, Haplaquepts, Typic, Ustocrepts, and Dystrochrepts), and 2c (stratified alluvium: < 1 degree slope, micro-relief, variable drainage, low areas subject to flooding, sandy loam to silty clay, Typic, and Fluventic).

Soil sampling and analysis

Surface soil samples (0 to 15 cm depth) were collected during May 2013 using a soil auger in the study area. Soil sampling locations were selected to best represent the land use condition in each VDC while considering terrain attributes and drainage conditions. A few VDCs such as Parsurampur and Golganj had only one soil sampling location. By following one soil sample per location, a total of 109 soil samples were collected from the study area, and the details of soil sampling locations are given in Fig 1. A global positioning system receiver with 1 m precision was used to record the longitude and latitude of soil sampling locations. No specific permissions were required for soil sampling in these locations and the field studies did not involve endangered or protected species.

The collected soil samples were air-dried and sieved through a 2 mm sieve for chemical analysis conducted at the Regional Soil Testing Laboratory, Kaski district of Nepal. The soil chemical parameters tested and methods used are given in Table 1. Sodium bicarbonate (NaHCO3) and ammonium acetate (C2H7NO2) were used as the extractants for laboratory analysis of available phosphorus and potassium, respectively.

[Figure omitted. See PDF.]

Table 1. Methods used for testing of soil chemical parameters at Regional Soil Testing Laboratory, Kaski district of Nepal.

https://doi.org/10.1371/journal.pone.0206350.t001

Statistical and geostatistical analysis

Descriptive statistics of soil properties, including mean, standard deviation, coefficient of variation, minimum, maximum, skewness (skew), and kurtosis, were calculated. For all measured soil characteristics, the visual method (histogram, boxplot and normal plot) and values of skew and kurtosis were used and no figures were included for the univariate test of normality in SAS software [33] prior to ordinary kriging.

This study focused on ordinary kriging (the general name for kriging), a linear geostatistical interpolation technique. Kriging estimates were calculated as weighted sums of the adjacent sampled concentrations. It is an improvement over inverse distance weighting (another geostatistical tool) interpolation because prediction estimates in kriging tend to be less biased and are accompanied by prediction standard errors [34]. Details of the kriging formula and calculation are given in Yao et al. [14]. The main application of geostatistics in soil science has been the estimation and mapping of soil attributes out of sampled areas [35].

Regardless of data distribution, kriging can provide the best-unbiased predictor of values at unsampled points, though data that have closer to a normal distribution can provide the best estimates of probability maps [36]. Therefore, it was necessary to normalize the dataset prior to geostatistical analysis because of high skew (Table 2) and the presence of outliers. Since the coefficient of skew was greater than one (except for pH and P2O5), the logarithmic transformation was applied for a kriging analysis (lognormal kriging, hereafter referred to as kriging) to stabilize the variance [35]. The logarithmic transformation resulted in smaller skew and kurtosis for OM, N, K2O, B, and Zn, and the transformed data passed the normality test.

[Figure omitted. See PDF.]

Table 2. Summary statistical overview for selected soil chemical properties of study area (N = 109), including original and log transformed data (for skew and kurtosis).

https://doi.org/10.1371/journal.pone.0206350.t002

A dbf file consisting of data for X and Y coordinates with respect to sampling site location was created in ArcGIS (version 10.2). Several digital map layers were then developed, using kriging in ArcMap, based on each soil chemical property at 1:25000 scale. The ranges for soil pH are classified as strongly acidic (<5.5), moderately acidic (5.5 to 6.2), neutral (6.2 to 7), moderately alkaline (7 to 7.8), and strongly alkaline (>7.8). Similarly, the rating charts for other soil parameters are given in Table 3, which is based on recommendations given by the Soil Management Directorate of the Department of Agriculture for the Terai region of Nepal [22].

[Figure omitted. See PDF.]

Table 3. Range for different soil parameters given by the Soil Management Directorate, Department of Agriculture for Terai region of Nepal.

https://doi.org/10.1371/journal.pone.0206350.t003

The kriging method uses semivariance to estimate the spatial distribution structure of the soil properties [37–38]. Semivariogram modeling and estimation are essential for structural analysis and spatial interpolation, which is akin to fitting a least-squares line in regression analysis [39]. It produces geostatistical parameters, including nugget, structural, sill, and range [38]. The spatial dependency (Sp. D) of soil parameters (the ratio of nugget to sill variances) is expressed as a percentage [40]. To ensure Sp. D, as a rule of thumb the sampling interval (lag) should be less than half of the range of the spatial variation [15]. If the ratio is less than 0.25, the variance has strong Sp. D and if the ratio ranges between 0.25 and 0.75, the variance has moderate Sp. D [41].

Moran’s I Index was used to measure spatial autocorrelation between sample points on the semivariagram cloud, which was evaluated using z-scores. Values greater than 1.96 or smaller than −1.96 are significant at p < 0.05 [42]. Similarly, the mean error (ME) and root mean square error (RMSE) was used for a cross-validation approach (or any given variogram model) to evaluate the accuracy or best fit of the kriging tool [43]. A ME value close to zero indicates that the interpolation method is unbiased. The lowest RMSE value indicates the best fit to the variogram model.

Data analysis

Descriptive statistics and geostatistics were used to analyze the dataset, and descriptive statistics along with a normality test were run in SAS software. All maps were produced using GIS software ArcMap (version 10.2) and its spatial analyst and geostatistical analysis extensions. The structure of spatial variability was analyzed through semivariogram. Spatial distribution was analyzed through kriging interpolation.

Results and discussion

Descriptive statistics for soil chemical properties

The commonly used descriptive statistical summary of the pH, OM, N, P2O5, K2O, Zn, and B is presented in Table 2. The variability was interpreted using the coefficient of variation (CV) and classified into most (CV: >35%), moderate (CV: 15 to 35%) and least (CV: <15%) variable ranges [44]. The CV ranged from 16.0% (in pH) to 97.12% (in K2O). The range of CV for the soil sampling locations suggested different degrees of heterogeneity among the properties studied.

The pH values were ranging from 4.2 to 8 with a mean of 6.4, which was also similar to the median value of 6.4. The concentration of OM was low (ranging from 1 to 2.5%), with a mean of 2.13%. Total N was relatively low (ranging from 0.05 to 0.01%) with a median of 0.09%, though the mean was 0.11%. Available P2O5 (40.08 kg ha-1) and K2O (110.61 kg ha-1) were within their respective medium ranges. Between the two micronutrients measured, Zn was very low (range: <0.25 mg kg-1) with a median of 0.07 mg kg-1 and a mean of 0.08 mg kg-1, while B was at medium (ranging from 0.5 to 1.2 mg kg-1) with a median of 0.99 mg kg-1, though the mean was 1.03 mg kg-1.

Among the soil chemical properties, OM, N, K2O, B, and Zn were found to be not normally distributed due to higher values of skew and kurtosis. Those datasets were then subjected to logarithmic transformation to narrow down the skew and kurtosis values (Table 2) and the transformed datasets were subsequently used in the spatial analysis.

Digital soil maps using kriging

Digital maps of soil chemical properties were produced by using kriging on the log transformed dataset, and the results (shown in Figs 2 through 8) were grouped into various classes based on the range representing their magnitude in the soil. The estimated area of each class is given in Table 4.

[Figure omitted. See PDF.]

Fig 2. Soil pH spatial variability map in southern part of Bara district, Nepal.

Most of the study area was with moderately alkaline (30.69%) followed by moderately acidic (22.91%) and neutral (22.80%) pH. Strongly alkaline was present in about 2.5% of total area but could not see in the variable map.

https://doi.org/10.1371/journal.pone.0206350.g002

[Figure omitted. See PDF.]

Fig 3. Soil OM spatial variability map in southern part of Bara district, Nepal.

Most of study area was with low (48.80%) and very low (26.88%) OM content.

https://doi.org/10.1371/journal.pone.0206350.g003

[Figure omitted. See PDF.]

Fig 4. Soil N spatial variability map in southern part of Bara district, Nepal.

Most of study area was with low (50.86%) for total N content.

https://doi.org/10.1371/journal.pone.0206350.g004

[Figure omitted. See PDF.]

Fig 5. Soil P2O5 spatial variability map in southern part of Bara district, Nepal.

Most of study area was with medium (42.95%) and low (29.88%) for available P2O5.

https://doi.org/10.1371/journal.pone.0206350.g005

[Figure omitted. See PDF.]

Fig 6. Soil K2O spatial variability map in southern part of Bara district, Nepal.

Most of study area was with low (41.19%), followed by very low (30.77%) and medium (21.79%) for available K2O.

https://doi.org/10.1371/journal.pone.0206350.g006

[Figure omitted. See PDF.]

Fig 7. Soil Zn spatial variability map in southern part of Bara district, Nepal.

Almost of the study area was with very low (98.33%) for Zn content.

https://doi.org/10.1371/journal.pone.0206350.g007

[Figure omitted. See PDF.]

Fig 8. Soil B spatial variability map in southern part of Bara district, Nepal.

Most of study area was with low (55.52%) and medium (43.21%) for B content.

https://doi.org/10.1371/journal.pone.0206350.g008

[Figure omitted. See PDF.]

Table 4. Areas under different soil categories based on soil fertility parameters.

https://doi.org/10.1371/journal.pone.0206350.t004

Soil pH.

Soil pH varied from strongly acidic (< 5.5) in 21.12% to strongly alkaline (> 7.8) in 2.48% of the total area (Table 4 and Fig 2). These results are in agreement with those reported in a recent study of soils of the Terai region [22, 45]. The variation in soil pH could be attributed to the nature of the alluvial parent material, micro topography, and the type and history of fertilizer used [46]. The losses of basic cation and other nutrients through erosion and leaching leaves the hydrogen and aluminum ions that can cause soil acidity [47]. Management practices such as crop nutrient uptake and harvest without replenishment [48] and poor crop residue management lowers the pH and leads to low levels of soil OM [28, 49].

Urea (46% N) and di-ammonium phosphate (18% N and 46% P) are the most commonly used fertilizers by Nepalese farmers. The national average for the use of chemical fertilizer has increased dramatically from 16.7 kg ha−1 in 2002 to 67.4 kg ha−1 in 2014 [50]. Of the major fertilizer nutrients, types of N fertilizer containing ammonium-N are the main factors affecting soil pH. As the ammonium-N in fertilizers undergoes nitrification, hydrogen ions are released, which can increase acidity [51], and Tripathi and Shrestha [52] reported an increase of acidity up to 4.1 (in 2000) from 5.6 (in 1997) after the application of fertilizers at four locations in western Nepal.

Plant growth and most soil processes are favored by a specific pH range. The low pH leads to Al and Mn toxicity, along with deficiency and/or unavailability of plant nutrients such as P, Ca, K, Mg, and Mo as observed by Dembele et al. [53] and Tisdale et al. [54]. To produce a sustained crop growth and yield, efforts should be made to increase the pH, which can be addressed through liming and OM management or adoption of the acid tolerant crops.

Soil organic matter.

Soil OM was relatively low (1 to 2.5%) in the majority (48.8%) of the study area, followed by very low (<1%) in 26.88% (Table 4 and Fig 3). The low organic content in the soils can generally be accounted for through the general sparse vegetation and competing use of crop residue as animal feed which then constrains their return to the soil [55–56]. A study conducted in the Dhading district of Nepal in 2003 showed that 37% of crop residue was used to feed livestock, 35% was used as fuel, 15% was burnt on lands, and the remaining 13% was incorporated into the soil through methods other than burning [57].

Another possible reason for low OM is a high soil OM decomposition rate resulting from soil and higher air temperature that decreases soil organic carbon (SOC). The SOC is affected by the addition of farm yard manure (FYM), tillage, and cropping pattern [58–59]. Around 14.56% of the study area revealed medium (2.5 to 5%) levels of soil OM, which could be due to waterlogged conditions, leading to shallow rooting and the confinement of biological activity to the upper soil layer. Similar results were reported by Shrestha et al. [60] for the soils of lowland irrigated rice fields in Nepal.

Among Nepalese farmers, there is an increase in the use of chemical fertilizer in agriculture, though this increase is not being matched by an increase in the use of organic manure (manures, organic fertilizers, compost, or other soil improvers) [61]. The present rate of organic manure application is 2.5 to 3 t ha-1 for soil fertility management [62], with an estimated composition of 0.5% N, 0.2% P, and 1.25% K on a dry weight basis, far below the global average and a rate that may not meet crop demand on a long-term basis [63–64]. As OM decreases, it also decreases available N, P, K, and some micronutrients [65]. Zhao et al. [66] reported that this low level of OM is indication of soil degradation and a high risk of soil erosion. Farmers should be encouraged to add much crop residues to the soil along with manure and compost.

Total nitrogen.

Usually, N has a greater effect on crop growth, crop quality, and yield. However, N was deficient in most of the areas with values <0.1 (low and very low) recorded in 70.5% of the total area (Table 4 and Fig 4). The variation in N content in different parts of the study area may be related to soil management, application of FYM and applied fertilizer to previous crops, etc. [67]. The acute deficiency of N is due to low OM content, increased rate of mineralization, and insufficient application of N fertilizer to nutrient exhausting crops like wheat and maize [46]. The rate of soil OM decomposition and N mineralization holds complex interactions with the microbial population and other environmental factors, mainly soil moisture and temperature. A field with 40 kg N ha-1 of soil nitrate build-up led to the loss of N from the entire field when the soil, which contained moisture levels > 46%, filled pore space at the onset of the monsoon rains in a lowland field in the west central part of Chitwan, Nepal [68].

Available phosphorus.

The available P2O5 was medium in 42.95% and low in 29.88% of the study area (Table 4 and Fig 5). The low level of OM may account for the low level of available P2O5 in the soils. However, the relatively higher availability of P2O5 observed in some areas may be due to the dissolution of Ca-P under neutral soil reaction under cultivated conditions [69–70]. Phosphorous is more directly affected by soil pH than other major plant nutrients such as N, K, and S; for example, at alkaline values, greater than pH 7.5, the HPO42- phosphate ions tend to react quickly with calcium (Ca) and magnesium (Mg) to form less soluble compounds. At acidic pH values, the H2PO4- phosphate ions react with aluminum (Al) and iron (Fe) to again form less soluble compounds [71]. Soils with inherent pH values between 6 and 7.5 are ideal for P availability. Besides pH, the amount of OM and the placement of P fertilizers also control the availability of P2O5, whereas erosion and runoff are associated with its loss from soil. Studies from many developed countries have shown that the use of flue gas desulfurization gypsum, a source of Ca and S, can be used as a soil amendment, especially to reduce soil and soluble P loss from agricultural fields and improve acidic soils [72]. Hence, whether or not farmers attempt to adjust pH levels, it is important to understand methods to increase the availability and use of added nutrients [73].

Available potassium.

The available K2O was at low levels in the majority of the study area (Table 4 and Fig 6). Soil pH also affects the availability of K2O. When soil pH is greater than 7, the greater Ca concentration increases the K availability through the displacement of exchangeable K by Ca. Conversely, when soil pH is less than 5.5, the reduction in Ca concentration reduce the K availability. In addition, low levels of OM due to low clay content, high hydraulic conductivity, and possible nutrient losses through leaching and erosion without replenishment also reduces the K level [74]. Water for irrigation to many of these study areas comes from Nepal’s rivers, which are flooded during monsoon season and carry heavy sediments (for example, mica) a source of exchangeable K [75]. However, due to year-round cropping practices, there is very little time for K to release from sediments and remain in the exchangeable site [76]. This could be another reason why the majority of the study area included low amounts of K2O.

Zinc and boron.

The micronutrient Zn was low and B was at medium level throughout the study area (Table 4 and Figs 7 and 8), possibly due to unfavorable soil pH (moderately alkaline in 30.69% of the total study area), intensive cropping, the use of high yielding varieties, and different fertilizer application strategies practiced by smallholder farmers. The Khaira disease (leaf bronzing) in rice due to Zn deficiency [77–78] and sterility in wheat induced by an inadequate B supply [79–82] are major concerns in the study area. In a study of micronutrient deficiencies in grain legumes, Srivastava et al. [83] found that B severely restricted the growth of lentils (Lens culinaris M.), chickpeas (Cicer arietinum L.), and pigeonpeas (Cajanus cajan L.) in the Terai region. Since rice is the major staple crop in Nepal, farmers use zinc sulphate (ZnSO4) before transplanting or sowing at the time of land preparation, along with a combination of ZnSO4 and lime during the growing stage, if the crop is infected [84].

Farmers’ practice and use of digital soil maps.

In many developing countries including Nepal, soil fertility management recommendations are solely based on soil types and agro-ecological zones. Details about soil pH range and the recommended agricultural lime rate, as well as the recommended doses of chemical fertilizers for specific crops in Nepal are given in Pandey et al. [22]. Despite the advisory recommendation made from research, farmers do not apply balanced doses of fertilizer, and instead use mostly acid forming nitrogenous fertilizers.

Most of the farmers apply FYM to their lands at the same rate as it is produced. The practice for FYM preparation and application is not an improved one because farmers dump FYM in open spaces and expose it to the sun, wind, and rain for several days before ploughing [85]. Farmers also follow nutrient-exhaustive high-yielding crop varieties under intensive cropping all year-round, leading a heavy loss of nutrients after every harvest. Therefore, a balanced rate of chemical fertilizers and organic manures must be applied every year.

Maps that characterize the spatial distribution of each soil property can be produced using kriging to group individual fields into potentially low- and high-productivity areas. Hence, management strategies to enhance soil nutrients could be implemented in the study area by using these DSMs as a guide [86], such as famers following fertilizer recommendations based on buildup and maintenance levels. Normally, nutrient values that are at low levels require relatively higher amount of fertilizer application; therefore, these DSMs may lead to proper understanding of existing farm soils by allowing easier management and maintaining the sustainability of productivity. This research sets a precedent for future DSM in other parts of the country.

Geostatistics for soil chemical properties

Semivariogram analysis.

The semivariogram model and some of the geostatistical parameters of soil chemical properties are shown in Table 5. Based on the lowest root mean square error (RMSE), different theoretical semivariogram models were selected for the significant fit of soil chemical properties [87]. An exponential model provided the best fit to the semivariogram of pH, OM, N, and Zn. The spherical model was the best fit to the semivariogram of K2O and B, whereas Gaussian was the best fit for P2O5. Many findings suggest that the exponential model is the most suitable for assessing spatial variability in soil chemical properties [88–92] because it explains the maximum variability in the spatial dataset [93–94].

[Figure omitted. See PDF.]

Table 5. Semivariance analysis of spatial structure in soil chemical properties.

https://doi.org/10.1371/journal.pone.0206350.t005

Table 6 shows that Sp. D of soil parameters ranged from 0.3 (in N) to 0.92 (in Zn). There was a moderate (in N, OM, P2O5, and pH) and weak (in K2O, Zn, and B) Sp. D. of the kriging model, which could be attributed to external factors such as variable rates of fertilizer application and incorporation of amendments by farmers within a cropped region. The ranges of spatial dependencies were large and vary between 4951 m for OM to 5945 m for P2O5 indicating that the optimum sampling interval varies greatly among different soil properties. Determination of the range values provides an idea of the correlation between different sampling locations, along with the maximum spatial dependence distance between them [95]. Fluctuation in the range with different lag sizes indicates that spatial structure may merely be regarded with a single model for semivariogram [96]. This difference may not be important for semivariance calculation, but it may be important if the purpose is to understand the underlying spatial structure of the data [97].

[Figure omitted. See PDF.]

Table 6. Test of significance of pattern analysis for selected soil chemical properties.

https://doi.org/10.1371/journal.pone.0206350.t006

Spatial autocorrelation.

The analysis of spatial autocorrelation based on Global Moran’s I Index was used to identify the spatial pattern soil chemical properties that may be dispersed, random, or clustered based on feature locations and attribute values simultaneously [42], as presented in Table 6. The hypothesis for the pattern analysis was that the soil chemical properties (including pH, OM and some nutrients) across the study area were randomly distributed.

According to ESRI [98], for the theory of random patterns, when the p-value is very small (in this study, p < 0.05) and the z-score is either very high or very low (1.96 < z and z < −1.96), the spatial pattern is not likely to reflect a random form of distribution. A positive Moran’s I index value indicates the neighboring values are similar, suggesting spatial dependency. A negative Moran’s I index value indicates the neighboring values are dissimilar, suggesting inverse spatial dependence. A Moran’s I index value of zero implies a lack of spatial pattern [99–101]. With the exception of K2O, all other soil variables had a positive Moran’s I index for their spatial pattern (Table 6).

Test of significance for values returned by the analysis of the major soil chemical properties indicated that pH, OM, N and Zn showed clustered distributions in the study area, with low levels clustered at one location and high levels at the other (no figures included). On the other hand, the pattern of distribution of P2O5, K2O and B did not appear significantly different from a random distribution at p < 0.05.

Conclusions

The application of the geostatistical approach, including descriptive statistics and semivariogram analysis, improved the description of spatial variability for soil chemical properties at 0 to 15 cm depth on a field scale. The descriptive statistics showed that most of the measured soil chemical variables were skewed and non-normally distributed and the available K2O data were highly variable (5 to 696 kg ha-1). Geostatistical interpolation identified that exponential, spherical, or Gaussian models provided the best fit to the semivariograms, depending on the soil chemical variable and, in general, showed weak or moderate spatial dependency for all of the variables. The kriging maps of soil chemical properties were found effective in explaining the distribution of soil properties in non-sampled locations based on sampled data. These maps aid farmers in to making efficient management decisions based on their proper understanding of the conditions of existing farm soils. These results show geostatistical analysis using kriging is an effective prediction tool for exploring the spatial variability of soil nutrients, and we recommend this tool for future soil sampling campaigns in Nepal.

Supporting information

[Figure omitted. See PDF.]

S1 File. An excel file include coordinates of 109 sampling locations and report of different soil chemical properties.

https://doi.org/10.1371/journal.pone.0206350.s001

(CSV)

Acknowledgments

The authors would like to thank Hexa International Pvt. Ltd., Lalitpur and En. Geo. Global Pvt. Ltd., Bhaktapur, Nepal for providing access and assistant with the datasets. Also thanks to the National Land Use Project, Ministry of Land Reform and Management, Government of Nepal for their supports. The residents of the study area who supported and provided valuable information are deeply acknowledged.

Our sincere thanks to Bharat Sharma Acharya, Nikita Bhusal, and Ian Rogers for their technical support during the preparation of the manuscript. Finally, we would like to thank anonymous reviewers and an academic editor for their valuable comments and suggestions which helped us on improving this paper.

Citation: Panday D, Maharjan B, Chalise D, Shrestha RK, Twanabasu B (2018) Digital soil mapping in the Bara district of Nepal using kriging tool in ArcGIS. PLoS ONE 13(10): e0206350. https://doi.org/10.1371/journal.pone.0206350

References

1. Webster R. The Development of Pedometrics. Geoderma. 1994;62(1–3):1–15.

2. Lagacherie P, McBratney AB, Voltz M. Digital soil mapping: an introductory perspective. Developments in soil science. Elsevier publication, The Netherlands. 2007.

3. McBratney AB, Minansy B, Malone BP, Sulaeman Y. Digital mapping of soil carbon. 19th World Congress of Soil Science, Soil Solutions for a Changing World, Brisbane, Australia. 2010.

4. Sarmadian F, Keshavarzi A, Odagiu A. Sampling design optimization based on soil-land inference model (SoLIM). ProEnvironment. 2013;6:556–561.

5. Martínez-Graña AM, Goy JL, Zazo C, Silva PG. Soil map and 3D virtual tour using a database of soil-forming factors. Environmental Earth Sciences. 2016;75(21):1402.

6. Sanchez PA, Ahamed A, Carre F, Zhang G. Digital soil map of the world. Science. 2009;325(5941):680–1. pmid:19661405

7. Malone BP, Minasny B, McBratney AB. Using R for digital soil mapping. Springer International Publishing, Switzerland. 2017.

8. Stoorvogel JJ, Bakkenes M, Temme A, Batjes NH, Brink B. S-World: A global soil map for environmental modelling. Land Degradation and Development. 2016;28(1):22–33.

9. McBratney AB, Mendonça Santos ML, Minasny B. On digital soil mapping. Geoderma. 2003;117:3–52.

10. Cambardella CA, Karlen DL. Spatial analysis of soil fertility parameters. Precision Agriculture. 1999;1:5–14.

11. Myers DE. Spatial interpolation: an overview. Geoderma. 1994;62(1–3):17–28.

12. Wang JF, Li LF, Christakos G. Sampling and kriging spatial means: efficiency and conditions. Sensors. 2009;9:5224–5240. pmid:22346694

13. Johnson CJ, Hurley M, Rapaport E, Pullinger M. Using expert knowledge effectively: lessons from species distribution models for wildlife conservation and management. In: Expert Knowledge and its Application in Landscape Ecology (eds Perera AH, Drew CA, Johnson CJ). Springer, New York. 2012;153–171.

14. Yao X, Fu B, Lu Y, Sun F, Wang S, Liu M. Comparison of four spatial interpolation methods for estimating soil moisture in a complex terrain catchment. PLoS One. 2013;8(1):e54660. pmid:23372749

15. Lopez-Granados F, Jurado-Exposita M, Pena-Barragan JM, Garcia-Torres L. Using geostatistical and remote sensing approaches for mapping soil properties. European Journal of Agronomy. 2005;23(3):279–289.

16. Santos-Francés F, Martínez-Graña A, Zarza CÁ, Sánchez AG, Rojo PA. Spatial Distribution of Heavy Metals and the Environmental Quality of Soil in the Northern Plateau of Spain by Geostatistical Methods. International journal of environmental research and public health. 2017;14(6):568.

17. Santos-Francés F, Martínez-Graña A, Alonso Rojo P, García Sánchez A. Geochemical Background and Baseline Values Determination and Spatial Distribution of Heavy Metal Pollution in Soils of the Andes Mountain Range (Cajamarca-Huancavelica, Peru). International journal of environmental research and public health. 2017;14(8):859.

18. Balkovič J, Rampašeková Z, Hutár V, Sobocká J, Skalský R. Digital soil mapping from conventional field soil observations. Soil and Water Research. 2013;8:13–25.

19. Malone BP, Styc Q, Minasny B, and McBratney AB. Digital soil mapping of soil carbon at the farm scale: A spatial downscaling approach in consideration of measured and uncertain data. Geoderma. 2017;290:91–99.

20. Zhang Q, Yang Z, Li Y, Chen D, Zhang J, Chen M. Spatial variability of soil nutrients and GIS-based nutrient management in Yongji County, China, International Journal of Geographical Information Science. 2010;24(7):965–981.

21. Zhu AX, Liu J, Du F, Zhang SJ, Qin CZ, Burt J, et al. Predictive soil mapping with limited sample data. European Journal of Soil Science. 2015;66:535–547.

22. Pandey S, Bhatta NP, Paudel P, Pariyar R, Maskey KH, Khadka J, et al. Improving fertilizer recommendations for Nepalese farmers with the help of soil-testing mobile van. Journal of Crop Improvement. 2018;32(1):19–32.

23. Panday D. Adapting Climate Change in Agriculture: The sustainable way in Nepalese context. Hydro Nepal Special Issue: Conference Proceedings. 2012;91–94.

24. Paudel B, Acharya BS, Ghimire R, Dahal KR, Bista P. Adapting agriculture to climate change and variability in Chitwan: long-term trends and farmers’ perceptions. Agricultural Research. 2014;3:165–174.

25. Adhikari BR. Flooding and inundation in Nepal Terai: issues and concerns. Hydro journal. 2013;12:59–65.

26. Jackson ML. Soil chemical analysis, New Delhi: Prentice Hall of India Pvt. Ltd. 1973.

27. Walkley A, Black IA. An examination of the Degtjareff method for determining organic carbon in soils: Effect of variations in digestion conditions and of inorganic soil constituents. Soil Science. 1934;63:251–263.

28. Bremner DC, Mulvaney JM. Total Nitrogen. In: Methods of Soil Analysis. (eds Page AL, Miller RH. Keaney DR). American Society of Agronomy. 1982;9(2).

29. Olsen SR, Cole CV, Watanabe FS, Dean LA. Estimation of available phosphorus in soils by extraction with sodium bicarbonate. U. S. Department of Agriculture Circular No. 939. Banderis, AD., DH. Barter, and K. Anderson. Agricultural and Advisor. 1954.

30. Toth SJ, Prince AL. Estimation of CEC and exchangeable Ca, K, and Na content of soil by Flame photometer technique. Soil Science. 1949;67:439–445.

31. Berger KC, Truog E. Boron determination in soils and plants. Industrial and Engineering Chenistry, Analytical Edition. 1939;11:540–545.

32. Lindsay WL, Norvell WA. Development of a DTPA soil test for zinc, iron, manganese, and copper. Soil Science Society of America Journal. 1978; 42:421–428.

33. SAS. SAS 9.4 in-database products: User’s guide, fifth edition. 2015.

34. Yasrebi J, Saffari M, Fathi H, Karimian N, Moazallahi M, Gazni R. Evaluation and comparison of ordinary kriging and inverse distance weighting methods for prediction of spatial variability of some soil chemical parameters. Research Journal of Biological Sciences. 2009;4(1):93–102.

35. Goovaerts P. Geostatistics in soil science: state-of-the-art and perspectives. Geoderma, 1999;89(1–2):1–45.

36. Tziachris P, Metaxa E, Papadopoulos F, Papadopoulou M. Spatial Modelling and Prediction Assessment of Soil Iron Using Kriging Interpolation with pH as Auxiliary Information. ISPRS International Journal of Geo-Information. 2017;6(9):283.

37. Zandi S, Ghobakhlou A, Sallis P. Evaluation of spatial interpolation techniques for mapping soil pH. In International Congress on Modeling and Simulation. 2011;1153–1159.

38. Wang YQ, Shao MA. Spatial variability of soil physical properties in a region of the Loess Plateau of PR China subject to wind and water erosion. Land Degradation & Development. 2013;24(3):296–304.

39. Chen S, Guo J. Spatial interpolation techniques: their applications in regionalizing climate-change series and associated accuracy evaluation in Northeast China. Geomatics, Natural Hazards and Risk. 2017; 8(2):689–705.

40. Al-Omran AM, Al-Wabel MI, El-Maghraby SE, Nadeem ME, Al-Sharani S. Spatial variability for some properties of the wastewater irrigated soils. Journal of the Saudi Society of Agricultural Sciences. 2013;12(2):167–75.

41. Orman EE. Improving the prediction accuracy of soil mapping through geostatistics. International Journal of Geosciences. 2012;3:574–590.

42. Moran PA. Notes on continuous stochastic phenomena. Biometrika. 1950;37:17–23. pmid:15420245

43. Johnston K, Ver H, Jay M, Krivoruchko K, Lucas N. Using ArcGIS Geostatistical Analyst. Environmental Systems Research Institute, Redlands, CA. 2001.

44. Wilding LP. Spatial Variability: its documentation, accommodation, and implication to soil surveys. In Nielsen, DR. Bouma, J. (Eds.). Soil Spatial Variability, Pudoc, Wageningen, Netherlands. 1985.

45. Khadka D, Lamichhane S, Shrestha SR, Pant BB. Evaluation of soil fertility status of Regional Agricultural Research Station, Tarahara, Sunsari, Nepal. Eurasian Journal of Soil Science. 2017;6(4):295.

46. Vasu D, Singh SK, Sahu N, Tiwary P, Chandran P, Duraisami VP, et al. Assessment of spatial variability of soil properties using geospatial techniques for farm level nutrient management. Soil and Tillage Research. 2017;169:25–34.

47. Brady A.C., Weil R.R. (2007): The Nature and Properties of Soils (13th ed.) Pearson Prentice-Hall Education Inc. 976.

48. Westarp SV, Sandra Brown HS, Shah PB. Agricultural intensification and the impacts on soil fertility in the Middle Mountains of Nepal. Canadian journal of soil science. 2004 Aug 1;84(3):323–32.

49. Regmi BD, and Zoebisch MA. 2004. Soil fertility status of Bari and Khet land in a small watershed of middle hill region of Nepal. Nepal Agricultural Research Journal 2004; 5:38–44.

50. WBG. The World Bank Group. Fertilizer consumption (kilograms per hectare of arable land). 2017; [Online] Accessed January 22, 2018. http://data.worldbank.org/indicator/AG.CON.FERT.ZS

51. Moore GA. Soilguide (Soil guide): A handbook for understanding and managing agricultural soils. 2001.

52. Tripathi BP, Shrestha SP. Nitrogen content in farm yard manure and its effects on the productivity and soil properties of rice-wheat, upland blackgram and maize-fingermillet systems. 2000. Lumle Working Paper No. 2000/14, Lumle Agricultural Research Center, Kaski, Nepal. 8 pages.

53. Dembele D, Traore K, Quansh C, Jnr E, BA B, Ballo M. Optimizing soil fertility management decision in Mali by remote sensing and GIS. Donnis Journal of Agricultural Research. 2016;3(4):22–34.

54. Tisdale SL, Nelson WL, Beaton JD. Soil and fertilizers. 4th ed. MacMillan Publications Inc, New York. 1985.

55. Paudel GS, Thapa GB. Changing farmers’ land management practices in the hills of Nepal. Environmental Management. 2001;28(6):789–803. pmid:11915967

56. Awasthi KD, Singh BR, Sitaula BK. Profile carbon and nutrient levels and management effect on soil quality indicators in the Mardi watershed of Nepal. Acta Agriculturae Scandinavica Section B-Soil and Plant. 2005 Sep 1;55(3):192–204.

57. Pandit K, Balla MK. An Assessment of soil Fertility Management Issues in Pokhare Khola Watershed, Dhading. Nepal Journal of Science and Technology. 2006;7:89–96.

58. Shrestha RK, Ladha JK, Gami SK. Total and organic soil carbon in cropping systems of Nepal. Nutrient cycling in agroecosystems. 2006 Jul 1;75(1–3):257–69.

59. Ghimire R, Lamichhane S, Acharya BS, Bista P, Sainju UM. Tillage, crop residue, and nutrient management effects on soil organic carbon in rice-based cropping systems: A review. Journal of integrative agriculture. 2017;16(1):1–5.

60. Shrestha BM, Sitaula BK, Singh BR, Bajracharya RM. Soil organic carbon stocks in soil aggregates under different land use systems in Nepal. Nutrient cycling in agroecosystems. 2004;70(2):201–13.

61. Khadka YG, Rai SK, Raut S. Long term effects of organic and inorganic fertilizers on rice under rice-wheat cropping sequence. Nepal Journal of Science and Technology. 2008;9:7–13.

62. NARC. The objectives of soil test and methods to take sampling for test (in Nepali). Lalitpur, Nepal: Nepal Agricultural Research Council Extension Fact Sheet. Soil Science Division. 2013.

63. Bishwakarma BK, Dahal NR, Allen R, Rajbhandari NP, Dhital BK, Gurung DB, et al. Effects of improved management and quality of farmyard manure on soil organic carbon contents in small-holder farming systems of the Middle Hills of Nepal. Climate and development. 2015;7(5):426–36.

64. Regmi AP. Effects of long term application of mineral fertilizers and manure on rice-wheat system. 2000. [Online] Accessed July 14, 2017. http://libcatalog.cimmyt.org/download/cim/67242.pdf

65. Wang J, Raman H, Zhang G, Mendham N, Zou M. 2006. Aluminium tolerance in barely (Horidium vulgarie L.): Physiological mechanisms, genetics and screening methods. Journal of Zhejiang University Science. 2006;7:769–787.

66. Y Zhao, He X, X Huang, Y Zhang, X Shi. Increasing soil organic matter enhances inherent soil productivity while offsetting fertilization effect under a rice cropping system. Sustainability. 2016;8:879.

67. Sherchan DB, and Gurung BD. An Integrated Nutrient Management System for Sustaining Soil Fertility: Opportunities and Strategy for Soil Fertility Research in the Hills. Challenges in Mountain Resource Management in Nepal. Processes, Trends and Dynamics in Middle Mountain Watersheds. Proceedings of a Workshop held in Kathmandu, Nepal 10–12 April, 1995 (Eds. H. Schreier, P. B. Shah, and S. Brown). ICIMOD/IDRC/UBC, Kathmandu. 1995;50–62.

68. Pande KR, Becker M. Seasonal soil nitrogen dynamics in rice‐wheat cropping systems of Nepal. Journal of Plant Nutrition and Soil Science. 2003;166(4):499–506.

69. Barrow CJ. Land degradation: Development and break-down of terrestrial environments. Cambridge University Press, Cambridge. 1991.

70. Pal DK, Wani SP, Sahrawat KL. Vertisols of tropical Indian environments: Pedology and edaphology. Geoderma. 2012;189–190:28–49.

71. IPNI. Soil pH and the availability of plant nutrients. International Plant Nutrition Institute. 2010, 2. [Online] Accessed January 22, 2018. http://www.ipni.net/ipniweb/pnt.nsf/0/97c1b6659f3405a28525777b0046bcb9/$FILE/Plant%20Nutrition%20Today%20Fall%202010%202.pdf

72. Panday D, Ferguson RB, Maharjan B. Flue Gas Desulfurization (FGD) Gypsum as Soil Amendment. In: Soil Amendments for Sustainability: Challenges and Perspectives (eds Rakshit A, Sarkar B, Abhilashis PC. CRC Press, FL. 2018;199–208.

73. McCauley A, Jones C, Jacobsen J. Soil pH and organic matter. Nutrient management module. 2009;8:1–2.

74. Chopart JJ, Azevedo MCB, Mezo LL, Mariaon D. Sugarcane root system depth in three different countries. Proceedings- International Society of Sugarcane Technololgy. 2012;27: 1–6.

75. Bajwa MI. Soil K status, K fertilizer usage and recommendations in Pakistan. Potash Rev. 1994;3:1–8.

76. Karki KB, MentlerA,Blum WEH. Food security and crop productivity in Kathmandu valley, Nepal. In: Proceedings of Internutionul Mbrkshop in Food Security of Urbun und Peri-urban Systems in Developing Countries, 15–18 November 2000, Vienna. Austria. 2000.

77. Mainuri ZG, Owino JO. Linking landforms and land use to land degradation in the Middle River Njoro Watershed. International Soil and Water Conservation Research. 2014;2(2):1–10.

78. Chaudhary RP, Subedi BP, Vetaas OR, Aase TH. Vegetation and society: their interaction in the Himalayas. Proceedings of the workshop on Bergen-Tribhuvan Human Ecology Programme, Kathmandu, Nepal. 2001.

79. Wissuwa M, Ismail AM, Yanagihara S. Effects of zinc deficiency on rice growth and genetic factors contributing to tolerance. Plant Physiology. 2006;142(2):731–741. pmid:16905666

80. Rerkasem B, Jamjod S. Boron deficiency induced male sterility in wheat (Triticum aestivum L.) and implications for plant breeding. Euphytica. 1997;96(2):257–262.

81. Huang L, Pant J, Dell B, Bell RW. Effects of boron deficiency on anther development and floret fertility in wheat (Triticum aestivum L.). Annals of Botany. 2000;85:493–500.

82. Abdel-Motagally FMF, El-Zohri R. Improvement of wheat yield grown under drought stress by boron foliar application at different growth stages. Journal of the Saudi Society of Agricultural Sciences. 2016.

83. Srivastava SP, Johansen C, Neupane, RK, Joshi M. Severe boron deficiency limiting grain legumes in the inner Terai of Nepal. In Proceedings of an International Workshop: Micronutrients in South and South East Asia, 8–11 September, 2004, Kathmandu, Nepal. ICIMOD. Nepal. 2005;67–76.

84. Andersen P. Micronutrient strategies for marginal areas. Revised version of a paper presented at the IGU international conference on geographical marginality: opportunities and constraints, 3–9 February, 2003 Geographi Bergen Arbelder fra Institutt for geografi—Bergen No 256–2003, Nepal. 2003.

85. Thapa GB, Paudel GS. Farmland degradation in the mountains of Nepal: a study of watersheds ‘with’and ‘without’external intervention. Land degradation & development. 2002;13(6):479–93.

86. Antwi M, Duker AA, Fosu M, Abaidoo RC. Geospatial approach to study the spatial distribution of major soil nutrients in the Northern region of Ghana. Cogent Geoscience. 2016;2:1201906.

87. Robertson GP. GS+: Geostatistics for the environmental sciences. Gamma Design Software, Plainwell, Michigan USA. 2008:165.

88. Reza SK, Sarkar D, Baruah U, Das TH. Evaluation and comparison of ordinary kriging and inverse distance weighting methods for prediction of spatial variability of some chemical parameters of Dhalai district, Tripura. Agropedology. 2010;20(1):38–48.

89. Bhunia GS, Shit PK, Maiti R. Comparison of GIS-based interpolation methods for spatial distribution of soil organic carbon. Journal of the Saudi Society of Agricultural Sciences. 2016.

90. Liu D, Wang Z, Zhang B, Song K, Li X, Li J, et al. Spatial distribution of soil organic carbon and analysis of related factors in croplands of the black soil region, Northeast China. Agriculture, Ecosystems and Environment. 2006;113(1–4):73–81.

91. Varouchakis EA, Hristopulos DT. Comparison of stochastic and deterministic methods for mapping groundwater level spatial variability in sparsely monitored basins. Environmental Monitoring and Assessment. 2013;185(1):1–19. pmid:22311559

92. Venteris ER, Basta NT, Bigham JM, Rea R. Modeling spatial patterns in soil arsenic to estimate natural baseline concentrations. Journal of Environmental Quality. 2013;43(3):936–946.

93. Tripathi R, Nayak AK, Shahid M, Raja R, Panda BB, Mohanty S, et al. Characterizing spatial variability of soil properties in salt affected coastal India using geostatistics and kriging. Arabian Journal of Geosciences. 2015;8(12):10693–10703.

94. Lark RM. Estimating variograms of soil properties by the method-of-moments and maximum likelihood. European Journal of Soil Science. 2000;51:717–728.

95. Akpa S, Odeh I, Hartemink A. Digital mapping of soil particle-size fractions for Nigeria. Soil Science Society of America Journal. 2014;78:1953–1966.

96. Silva RA, Siqueira GM, Costa MK, Guedes Filho O, e Silva ÊF. Spatial Variability of Soil Fauna under Different Land Use and Managements. Rev Bras Cienc Solo. 2018;42:e0170121.

97. Chung SO, Sudduth KA, Drummond ST, Kitchen NR. Spatial variability of soil properties using nested variograms at multiple scales. Journal of Biosystems Engineering. 2014;39(4):377–88.

98. ESRI. The principles of geostatistical analysis (3) 2010 [Internet]. 2017. http://maps.unomaha.edu/Peterson/gisII/ESRImanuals/Ch3_Principles.pdf

99. Lloyd CD. Nonstationary models for exploring and mapping monthly precipitation in the United Kingdom. International Journal of Climatology. 2010;30: 390–405.

100. Saxena R, Nagpal BN, Das MK, Srivastava A, Gupta SK, Kumar A, et al. 2012. A spatial statistical approach to analyze malaria situation at micro level for priority control in Ranchi district, Jharkhand. Indian Journal Medical Research. 2012;136(5):776–782.

101. Al-Ahmadi K, Al-Zahrani A. Spatial autocorrelation of cancer incidence in Saudi Arabia. International Journal of Environmental Research and Public Health. 2013;10:7207–7228. pmid:24351742