Data analysis

Using SAS 9.2 software package [51], the descriptive statistics of studied soil parameters were obtained. The relations among the studied soil parameters were evaluated by Pearson’s correlation coefficient analysis. The datasets pertaining to AS, AZn, AFe, ACu, AMn, and AB were checked for normal distribution [52], prior to geostatistical analysis, and data transformations were carried out. The semivariogram (Eq 1) analysis was carried out using ArcGIS 10.6 software (Esri, Redlands, California, United States) for characterizing spatial structure of AS, AZn, AFe, ACu, AMn, and AB [39]. The semivariogram draws a plot between the variance (of spatially separated points of each dataset) and the lag distance. The semivariogram at h distance interval is γ(h). The sample pair value at h distance interval is N(h). The symbols z(xi) and z(xi+h) denote the sample points separated by a distance h.

The best fitted semivariogram models were chosen by cross-validation and depending upon the value of root mean square error (RMSE) (Eq 2).

Where, z(xi, yi) is observed value, z _* (xi, yi) is predicted value and n is number of observations.

The widely used technique of ordinary kriging (OK) was used in order to develop distribution maps. OK estimates the values at unsampled locations based on optimal, unbiased and linear estimation [53]. It is built on the assumptions of second-order stationary and spatial autocorrelation. The accuracy of interpolated maps was verified by goodness-of-prediction-criterium (G value expressed in %) [35]. The positive G value indicates more accuracy of interpolated map from the samples than the area average.