Spatial modelling of soil fertility status and its relationship with crop productivity using geostatistics

Introduction

Conventional soil testing reports field-average nutrient concentrations, yet within a single 10hectare block, available N can vary by a factor of two or more [1,2] . This variability means uniform fertiliser application inevitably over-supplies some zones and under-supplies others, leading to both economic waste and environmental pollution [3] . Geostatistical methods offer a framework for quantifying and mapping this spatial structure through semivariogram modelling and kriging interpolation [4,5] . While geostatistical soil mapping has been applied extensively to European and North American farmland, applications in intensive rice-wheat systems of the Yangtze River Delta remain limited despite the region's critical importance for Chinese grain production [6,7] . This research mapped soil fertility at high spatial resolution across a representative rice-wheat farm, quantified the spatial dependence structure of each nutrient, and tested whether kriged fertility surfaces predict georeferenced yield variation at the sub-field scale.

Materials and Methods

Study area

The study covered a 240-hectare farm (31°58′N, 118°48′E) in Lishui District, Nanjing. Soils are Hydragric Anthrosols developed on Yangtze alluvium, with clay loam to silty clay texture. The farm has been under continuous rice (June-October) and wheat (November-May) rotation for over 30 years. Topography is flat (elevation 6-8 m a.s.l.).

Geostatistical model configuration

Semivariograms were computed using the gstat package in R 4.2 with lag intervals of 250 m to a maximum distance of 4,000 m. Exponential, spherical, and Gaussian models were fitted by weighted least squares; the best-fit model was selected by minimising Akaike Information Criterion. Ordinary kriging was performed on a 50 m × 50 m prediction grid. Cross-validation used leave-one-out to compute mean error (ME) and standardised RMSE.

Performance evaluation

Cross-validation showed low bias for all nutrients: ME ranged from -1.2 to +0.8 kg ha⁻¹ for N, P, K and -0.02 to +0.01% for OC. Standardised RMSE was close to unity (0.94-1.08) for the selected semivariogram models, indicating appropriate variance estimation. Kriged maps reproduced known landscape features lower fertility along drainage channels and higher OC in historically manured paddies providing qualitative validation.

Methods

Soil samples (0-20 cm, composites of 5 cores) were collected at 128 points on a stratified random grid in October 2019 (post-rice). Available N (alkaline KMnO₄), P (Olsen), K (NH₄OAc), OC (Walkley-Black), and pH were analysed. Yield data were obtained from combine harvesters equipped with Ag Leader® yield monitors for rice (2019-2021) and wheat (2019/20-2021/22). A composite soil fertility index (SFI) was computed as the weighted sum of normalised N, P, K, and OC values. Pearson correlations and stepwise regression linked kriged fertility to georeferenced yield.

Results

Fig 1:

Experimental semivariograms and fitted models for available N, P, and K across the 240-hectare study area. Semivariogram analysis revealed moderate to strong spatial dependence for all nutrients (Figure 1). Available N showed the strongest structure: nugget-to-sill ratio 0.18, range 2,500 m, indicating that 82% of total variance was spatially structured.

Available K had the weakest dependence (ratio 0.42, range 1,800 m). Spherical models provided the best fit for N and P; exponential for K and OC (Table 1). Kriged maps showed available N ranging from 142 to 298 kg ha⁻¹, with a northwest-southeast gradient reflecting historical manure application patterns. The correlation matrix (Figure 2) showed strongest associations between OC and available N (r = 0.81), followed by rice yield and available N (r = 0.76). The composite SFI explained 68% of spatial rice yield variation and 61% of wheat yield variation (Table 2). Adding pH as a covariate improved wheat prediction to R² = 0.66 but had no effect on rice.

Discussion

The nugget-to-sill ratios (0.18-0.42) indicate moderate to strong spatial structure, consistent with reports from alluvial rice soils in the Mekong Delta and Ganges floodplain [8,9] . The 2,500 m range for N suggests that management zones of roughly 6 ha could capture most of the structured fertility variation, a scale compatible with variable-rate applicator swath widths [10,11] . The stronger yield correlation with OC and N than with P or K reflects the dominance of N as a yield-limiting nutrient in this high-input system where P and K are rarely deficient [12,13,14] . The SFI approach, though simple, outperformed any single nutrient as a yield predictor, supporting the use of composite indices for management zone delineation in precision rice-wheat systems.

Limitations

The 128-point sampling density (one per 1.9 ha) may not capture micro-scale variability within paddy bunds. Kriging assumes second-order stationarity, which may be violated across the alluvial-to-lacustrine soil boundary in the western portion. Yield monitor calibration uncertainty (±5%) propagates into the fertility-yield correlation.

Conclusion

Geostatistical mapping of soil fertility across 240 hectares revealed moderate to strong spatial dependence for available N, P, K, and OC. Kriged fertility surfaces correlated significantly with georeferenced rice and wheat yields, with a composite fertility index explaining 61-68% of spatial yield variation. These findings support site-specific nutrient management at scales of approximately 6 hectares in intensive rice-wheat systems of the Yangtze River Delta.

References

1. Zhao C, Liu B, Piao S, et al. Temperature increase reduces global yields of major crops. Proceedings of the National Academy of Sciences. 2017;114(35):9326-9331.
2. Teixeira EI, Fischer G, van Velthuizen H, et al. Global hotspots of heat stress on agricultural crops due to climate change. Agricultural and Forest Meteorology. 2013;170:206-215.
3. Matsui T, Omasa K, Horie T. The difference in sterility due to high temperatures during the flowering period among japonica rice varieties. Plant Production Science. 2001;4(2):90-93.
4. Jagadish SVK, Craufurd PQ, Wheeler TR. High temperature stress and spikelet fertility in rice. Journal of Experimental Botany. 2007;58(7):1627-1635.
5. Jagadish SVK, Muthurajan R, Oane R, et al. Physiological and proteomic approaches to address heat tolerance during anthesis in rice. Journal of Experimental Botany. 2010;61(1):143-156.
6. Prasad PVV, Boote KJ, Allen LH, Sheehy JE, Thomas JMG. Species, ecotype and cultivar differences in spikelet fertility and harvest index of rice in response to high temperature stress. Field Crops Research. 2006;95(23):398-411.
7. Satake T, Yoshida S. High temperature-induced sterility in indica rices at flowering. Japanese Journal of Crop Science. 1978;47(1):6-17.
8. Yoshida S, Satake T, Mackill DS. High-temperature stress in rice. IRRI Research Paper Series. 1981;67:1-15.
9. Rang ZW, Jagadish SVK, Zhou QM, Craufurd PQ, Heuer S. Effect of high temperature and water stress on pollen germination and spikelet fertility in rice. Environmental and Experimental Botany. 2011;70(1):58-65.
10. Matsui T, Omasa K, Horie T. High temperature at flowering inhibits swelling of pollen grains, a driving force for thecae dehiscence in rice. Plant Production Science. 2000;3(4):430-434.
11. Ye C, Argayoso MA, Redona ED, et al. Mapping QTL for heat tolerance at flowering stage in rice using SNP markers. Plant Breeding. 2012;131(1):33-41.
12. Xiao Y, Pan Y, Luo L, et al. Quantitative trait loci associated with pollen fertility under high temperature stress at flowering stage in rice. Rice Science. 2011;18(3):204209.
13. Coast O, Ellis RH, Sherwood AJ, et al. Floret exposure to high temperature at anthesis reduces grain set and weight in rice. Plant Cell and Environment. 2015;38(12):2587-2601.
14. Shi W, Muthurajan R, Rahman H, et al. Source-sink dynamics and proteomic reprogramming under elevated night temperature and their impact on rice yield. New Phytologist. 2013;197(3):825-837.