Field-Scale Spatial Modeling of Soil Properties in the Harran Plain Using IDW and Auxiliary Variable-Based Co-Kriging

This study aimed to evaluate the field-scale spatial variability of selected soil properties in a 30-decare agricultural area located in the southeastern Harran Plain, Şanlıurfa, Türkiye. A total of 54 surface soil samples (0–30 cm) were collected using a 25 × 25 m grid design. Soil physical (clay, silt, sand) and chemical properties (pH, electrical conductivity (EC), lime, organic matter (OM), and available phosphorus (P)) were analyzed.Spatial distribution was modeled using a deterministic interpolation method (Inverse Distance Weighting, IDW) and a multivariate geostatistical approach (Co-Kriging, COK). The results indicated pronounced spatial heterogeneity among soil parameters. Compared with IDW, Co-Kriging reduced RMSE values particularly for organic matter (0.34% → 0.33%) and available phosphorus (0.67 → 0.60 kg da⁻¹). The highest prediction improvement was observed for phosphorus, with a 10.44% RMSE reduction relative to IDW. Auxiliary variables (EC for OM and lime for P) were incorporated based on their significant correlation with the target variables.Spatial maps revealed localized distribution patterns for clay, lime, organic matter, and phosphorus, supporting site-specific nutrient management and fertilization zoning. The findings demonstrate that integrating auxiliary variables within a multivariate geostatistical framework enhances prediction performance at the field scale and improves the reliability of spatial soil assessments.