Spatial variability of soil fertility parameters and their correlation with crop productivity using geospatial techniques

Introduction

Uniform fertilizer application across heterogeneous agricultural fields wastes resources, reduces profitability, and contributes to environmental pollution through nutrient losses to groundwater and atmosphere [1] . Soil properties vary spatially due to the combined influences of parent material, topography, climate, biological activity, and management history, creating complex patterns that traditional sampling approaches fail to capture adequately [2] . Understanding this spatial variability represents the foundation for precision agriculture approaches that match nutrient inputs to site-specific requirements. Geographic Information Systems (GIS) combined with geostatistical analysis provide powerful tools for characterizing and mapping soil spatial variability [3] . Kriging interpolation, based on regionalized variable theory, generates optimal predictions at unsampled locations while quantifying prediction uncertainty, advantages unavailable through conventional interpolation methods [4] . Research across diverse agricultural systems has demonstrated that geostatistical approaches can identify spatial patterns in soil fertility that correlate meaningfully with crop yield variation, enabling zone-based management strategies [5] . The Indo-Gangetic plains of northern India, among the world's most productive agricultural regions and the breadbasket of the nation, face mounting challenges from soil degradation, declining organic matter content, and nutrient imbalances following decades of intensive rice-wheat cultivation [6] . But despite the recognized importance of soil fertility management in this region, systematic characterization of spatial variability using modern geospatial techniques remains limited for many districts. This research was therefore conducted with objectives to quantify spatial variability of key soil fertility parameters across an irrigated agricultural area in Karnal district of Haryana, to establish correlations between soil fertility indicators and wheat crop productivity, and to develop fertility zone maps suitable for guiding site-specific nutrient management recommendations.

Research Area Characterization

The research was conducted across 2,400 hectares of irrigated agricultural land in the Nissing block of Karnal district, Haryana, India (29°38'N to 29°44'N latitude, 76°52'E to 77°00'E longitude). This area lies within the alluvial plain of the Yamuna River, characterized by nearly level topography with elevation ranging from 245 to 252 m above sea level. The climate is classified as semi-arid subtropical with hot summers and cold winters, mean annual temperature of 24.2°C, mean annual precipitation of 680 mm concentrated during the monsoon months of July through September, and potential evapotranspiration exceeding 1,800 mm annually [7] . Soils are predominantly Typic Ustochrepts developed from Indo-Gangetic alluvium, with textures ranging from sandy loam to loam. The area has been under continuous irrigated rice-wheat rotation for over four decades, a system that has provided food security but also raised concerns about soil health decline. Canal and tubewell irrigation provide water during both kharif and rabi seasons, though groundwater quality concerns have increased in some pockets due to excessive extraction. Land holdings are fragmented into small parcels averaging 2-4 hectares, with management intensity and fertilizer application rates varying substantially among farmers based on economic capacity and awareness levels [8] .

Material and Methods

Material: Soil sampling was conducted during October 2019 prior to wheat sowing, following harvest of the preceding kharif rice crop. A systematic grid sampling design with 400 m spacing generated 150 sampling points across the research area. At each georeferenced location, composite samples were collected from 0-30 cm depth by combining five subsamples taken within a 2 m radius using a stainless steel auger. GPS coordinates were recorded using Garmin GPSMAP 64s receivers with horizontal accuracy of ±3 m. Samples were air-dried, ground to pass 2 mm sieve, and analyzed at the Soil Testing Laboratory following standard protocols recommended by the Indian Council of Agricultural Research [9] .

Methods

Laboratory analyses included organic carbon by Walkley-Black wet oxidation, available nitrogen by alkaline permanganate method, available phosphorus by Olsen bicarbonate extraction, available potassium by ammonium acetate extraction, pH in 1:2.5 soil-water suspension, and electrical conductivity in saturation extract. Quality control included duplicate analyses for 10% of samples and certified reference materials. Wheat grain yield data were obtained from farmer records and field measurements at harvest in April 2020 and April 2021, with yields georeferenced to corresponding soil sampling locations through village land records and GPS verification [10] . Geostatistical analysis was performed using ArcGIS Pro 3.0 Geostatistical Analyst extension. Exploratory data analysis examined distributions, identified outliers, and assessed normality. Semivariograms were calculated using classical Matheron estimator and fitted with spherical, exponential, and Gaussian models, with optimal model selection based on cross-validation statistics. Ordinary kriging generated continuous prediction surfaces with 50 m cell resolution. Spatial dependence was classified based on nugget-to-sill ratio: strong (<0.25), moderate (0.25-0.75), or weak (>0.75). Cross-validation assessed prediction accuracy through root mean square error and standardized mean error. Correlation analysis examined relationships between soil parameters and crop yield, while cluster analysis delineated fertility management zones [11] . 1). Organic carbon ranged from 0.28% to 1.18% with mean of 0.64%, reflecting the generally depleted status of intensively cultivated Indo-Gangetic soils following decades of ricewheat rotation. The coefficient of variation (CV) of 32.8% for organic carbon indicates moderate to high spatial variability. Available phosphorus exhibited the highest variability (CV = 47.6%), likely reflecting differential fertilizer application histories among land parcels. In contrast, pH showed low variability (CV = 4.1%) typical of buffered calcareous soils of northwestern India. Electrical conductivity varied considerably (CV = 66.7%) with some locations approaching salinity thresholds, indicating localized problems in lower-lying areas with poor drainage.

Results

Comprehensive Interpretation

Semivariogram analysis revealed distinct spatial structures for all measured parameters (Table 2). Organic carbon exhibited the strongest spatial dependence with nugget ratio of 0.24, indicating that 76% of total variance is spatially structured and predictable through kriging. The range parameter of 1,840 m suggests Correlation analysis demonstrated that wheat yield was most strongly associated with organic carbon (r = 0.847, p<0.01), followed by available nitrogen (r = 0.782) and available phosphorus (r = 0.694) (Table 3, Figure 1). Electrical conductivity showed significant negative correlation with yield (r = -0.524), confirming that even moderate salinity levels affect productivity in some portions of the research area. Cluster analysis based on the integrated fertility index delineated five management zones: Zone A (high fertility, 23.4% of area), Zone B (medium-high, 18.6%), Zone C (medium, 12.8%), Zone D (medium-low, 24.2%), and Zone E (low, 21.0%). The spatial distribution of organic carbon across these zones is illustrated in Figure 3, showing clear differentiation in central tendency and variability among fertility categories.

Discussion

The substantial spatial variability documented in this research confirms that uniform fertilizer recommendations are inappropriate for the heterogeneous soils of the Indo-Gangetic alluvial plain. The coefficient of variation values exceeding 30% for organic carbon and 45% for available phosphorus fall within ranges reported for intensively cultivated irrigated systems across Punjab and Haryana, where management history creates variability superimposed on natural soil heterogeneity [3, 5]. The strong spatial dependence of organic carbon, with 76% of variance spatially structured, suggests that this parameter can be mapped effectively using the sampling density employed, enabling zone-based management recommendations.

The strong correlation between organic carbon and crop yield (r = 0.847) reflects the multiple functions of soil organic matter in supporting productivity, including nutrient supply, water retention, and soil structure maintenance [12] . In the semi-arid subtropical environment of Haryana, where high temperatures during wheat grain filling accelerate organic matter decomposition, maintaining adequate soil carbon presents particular challenges that continuous ricewheat rotation has failed to address. The correlation coefficients obtained in this research exceed those reported in some temperate region studies, possibly because soil organic matter limitations are more acute and therefore more directly yield-limiting under the harsh climatic conditions of northwestern India [13] .

Conclusions: This research has demonstrated that soil fertility parameters exhibit substantial spatial variability across the irrigated agricultural landscape of Karnal district, Haryana, with coefficients of variation ranging from 4% for pH to over 66% for electrical conductivity. The geostatistical analysis successfully characterized spatial structure, with organic carbon showing strong spatial dependence suitable for accurate kriging interpolation. The documented correlations between soil fertility indicators and wheat productivity provide quantitative relationships that can inform variable-rate fertilizer application strategies aligned with precision agriculture principles. The delineation of five fertility management zones based on integrated fertility index calculations offers a practical framework for site-specific nutrient management. The 45% of the research area falling in low to medium-low fertility categories requires enhanced organic matter inputs and balanced fertilization to restore productivity potential, with particular attention to building soil organic carbon through residue retention and organic amendments. Conversely, the 23% classified as high fertility may receive reduced fertilizer inputs without yield penalty, generating both economic savings and environmental benefits through reduced nutrient loading to groundwater. The methodology established in this research can be extended to adjacent districts of the Indo-Gangetic plains, with the semivariogram parameters providing guidance for optimal sampling density in future soil health surveys. The strong correlation between organic carbon and yield suggests that soil carbon mapping may serve as an efficient proxy for delineating productivity zones where comprehensive soil testing resources are limited. Implementation of zone-based management recommendations arising from this research should be accompanied by monitoring to validate yield responses and refine fertility zone boundaries over time. These findings support the Government of India's Soil Health Card scheme by demonstrating how site-specific recommendations can improve both productivity and sustainability of India's most important agricultural region.

Acknowledgement

Acknowledgements Funding SourcesThe authors acknowledge financial support from the Indian Council of Agricultural Research (ICAR) through the National Innovations on Climate Resilient Agriculture project. Additional resources were provided by the Haryana State Department of Agriculture and Farmers Welfare.

References

1. Mulla DJ, McBratney AB. Soil spatial variability. In: Sumner ME, editor. Handbook of Soil Science. Boca Raton: CRC Press; 2000. p. A321-A352.
2. Webster R, Oliver MA. Geostatistics for Environmental Scientists. 2nd ed. Chichester: John Wiley & Sons; 2007. p. 1-315.
3. Kerry R, Oliver MA. Variograms of ancillary data to aid sampling for soil surveys. Precis Agric. 2003;4(3):261-278.
4. Goovaerts P. Geostatistics for Natural Resources Evaluation. New York: Oxford University Press; 1997. p. 1-483.
5. Cambardella CA, Moorman TB, Novak JM, Parkin TB, Karlen DL, Turco RF, et al. Field-scale variability of soil properties in central Iowa soils. Soil Sci Soc Am J. 1994;58(5):1501-1511.
6. Ladha JK, Dawe D, Pathak H, Padre AT, Yadav RL, Singh B, et al. How extensive are yield declines in long-term rice-wheat experiments in Asia? Field Crops Res. 2003;81(2-3):159-180.
7. Pathak H, Ladha JK, Aggarwal PK, Peng S, Das S, Singh Y, et al. Trends of climatic potential and on-farm yields of rice and wheat in the Indo-Gangetic Plains. Field Crops Res. 2003;80(3):223-234.
8. Jat ML, Gathala MK, Saharawat YS, Tetarwal JP, Gupta R, Yadav OP, et al. Double no-till and permanent raised beds in maize-wheat rotation of northwestern Indo-Gangetic plains of India. Field Crops Res. 2013;149:291-299.
9. Sparks DL, Page AL, Helmke PA, Loeppert RH, Soltanpour PN, Tabatabai MA, et al., editors. Methods of Soil Analysis. Part 3: Chemical Methods. Madison (WI): Soil Science Society of America; 1996. p. 11390.
10. Corwin DL, Lesch SM. Apparent soil electrical conductivity measurements in agriculture. Comput Electron Agric. 2005;46(1-3):11-43.
11. Isaaks EH, Srivastava RM. Applied Geostatistics. New York: Oxford University Press; 1989. p. 1-561.
12. Lal R. Soil carbon sequestration to mitigate climate change. Geoderma. 2004;123(1-2):1-22.
13. Oldfield EE, Bradford MA, Wood SA. Global metaanalysis of the relationship between soil organic matter and crop yields. Soil. 2019;5(1):15-32.
14. McBratney AB, Mendonça Santos ML, Minasny B. On digital soil mapping. Geoderma. 2003;117(1-2):3-52.