Geostatistical modelling of physico-chemical properties of soils impacted by artisanal gold mining and farming: implications for soil environmental quality assessment

ABSTRACT

In this study, geostatistical techniques were used to investigate spatial variability of physico-chemical properties of soil samples from farming- and artisanal gold mining-affected areas. Based on coefficients of variation, chemical properties, except pH, had high variability compared to physical properties. This could be attributed to the local lithology, vegetation cover and land-use. Spatial distribution patterns of silt, clay, total nitrogen, pH and C/N could be controlled by extrinsic factors such as biological activity and land-use patterns, whereas those of electrical conductivity, organic matter (SOM) and sand were governed by both intrinsic and extrinsic factors. Because of the dominance of clay type, the distribution of cation exchange capacity (CEC) was spatially correlated with clay. About 80% of the area had acidic to circumneutral pH, which may increase mobility and toxicity of most cationic heavy metals. Distribution maps showed that microbial activities, expressed as C/N ratios, tend to be weak around farmlands, reflecting the impacts of soil management practices. In contrast, abundance of clay around the river bank is an indication of soil erosion and its redistribution from upland to depositional areas. Mean error and root mean square error of kriging were higher than those of cokriging for SOM, suggesting that cokriging was more reliable for spatial interpolation of this variable. The work reported in this study indicated that that 29%, 77% and 75% of the area had high probabilities to have pH, CEC and SOM contents below which most cationic heavy metals would be mobile and bioavailable in the soils.

KEYWORDS:

• Soil environmental quality
• ordinary kriging
• indicator kriging
• spatial dependency

1. Introduction

In recent years, misuse of fertilizers and pesticides in farming, high population growth and artisanal mining have been releasing huge quantities of heavy metals into soils worldwide (Bamba et al., 2013; Marrugo-Negrete et al., 2017; Sako et al., 2018). Once in the soil environment, these metals can be immobilized or be available to plants and animals and end up in the food web (Ashraf et al., 2017; Kumar et al., 2019). The bioavailability and geochemical behavior of heavy metals in soils are directly and indirectly influenced by soil physico-chemical properties such as soil organic matter (SOM) content, particle size distribution, pH, electrical conductivity (EC) and cation exchange capacity (CEC). Contrary to temperate soils, semi-arid tropical soils are byproducts of intensive chemical weathering over a pedological time scale through which primary minerals (e.g. felsic and ferromagnesian minerals) are converted into clay and residual minerals with near complete removal of base cations (da Silva et al., 2022; Hampl et al., 2023). Not only these cations are essential macronutrients for plant growth and development, their presence also increases soil CEC, which determines its buffering capacity to maintain metallic pollutants against leaching and subsequent migration to the groundwater systems (Aharonov-Nadborny et al., 2017; B. Huang et al., 2015; Puga et al., 2016; Suzuki et al., 2020). Accordingly, several studies have indicated that CEC is one of the key soil quality indicators that are sensitive to anthropogenic activities (e.g. Graber et al., 2017; Khaledian et al., 2017).

Given their large specific surface area, negative charges, ability to hold and exchange positively charged ions, clay minerals are major carriers of heavy metals in soils (Uddin, 2017). Unlike the soils with dominant coarse particles, clayey soils tend to retain heavy metals, and thus reduce their mobility and bioavailability (Zhang et al., 2020). Similarly, soils with high clay and SOM contents tend to have high CEC (Moral & Rebollo, 2017; Rahal & Alhumairi, 2019). Under high SOM contents, dissolved organic matter can be complexed with heavy metals, which, in turn, will improve soil structure and its adsorptive capacity (Jiang et al., 2021). Loss of SOM, due to farming, may destroy soil structure and imped water flow and atmospheric O₂ and soil CO₂ exchanges (Helliwell et al., 2014; Siedt et al., 2020). The subsequent buildup of CO₂ decreases soil pH, increases heavy metal concentrations and exacerbates their migration and toxicity (Al-Wabel et al., 2019; Ouyang et al., 2020). Moreover, the low pH and O₂ contents are likely to inhibit the growth and activity of soil microbial community, leading to a low humus production and formation of unstable soil aggregates (Obalum et al., 2017; Oldfield et al., 2018). Ultimately, this will affect soil structure and its ability to immobilize anthropogenically derived heavy metals. Another factor that may cause soil acidity and heavy metal content is the application of acidifying fertilizers with trace concentrations of heavy metals (Wei et al., 2020). Such agricultural practices are likely to increase heavy metal background levels and reduce soil pH (Ning et al., 2017). In the resulting acidified soils, nutrient loss and heavy metal concentrations will be increased (Anza et al., 2021; Kicińska et al., 2021). Therefore, among soil physico-chemical properties, pH is considered as one of the most important environmental indicators that can be used to estimate heavy metal mobility and toxicity in the soil environment (Maurya et al., 2020; Sungur et al., 2016).

Artisanal gold mining and its associated activities also have profound effects on the environmental quality of the surrounding soils (Compaore et al., 2019; Eludoyin et al., 2017; Tankari Dan-Badjo et al., 2019). That is, the compounded impacts of climate change on agricultural yields and high gold price, have led to an upsurge of artisanal gold mining sites across Burkina Faso in recent years. Notwithstanding artisanal gold mining contributes to poverty alleviation among rural communities through income generation, it causes a number of environmental concerns. These unmonitored activities profoundly affect land-use patterns, aggravate soil degradation and soil and water pollution (Essandoh et al., 2021; Torrance et al., 2021). Due to pedogenic processes, topography, bioclimatic conditions, underlying geology and land-use changes, soil environmental quality indicators usually show a strong spatial dependence (Babur et al., 2021). Therefore, a thorough evaluation of the spatial variability of the indicators at macro- and micro-scales and identification of their potential sources by using adequate modeling tools is crucial for setting up site-specific pollution monitoring strategies. However, conventional statistical techniques such univariate (e.g. coefficient of variation) and multivariate analyses (e.g. correlation analysis and factor analysis), which have been used to model spatial variability, are inadequate because they assume that the measured variables are independent of their spatial distribution (Hengl et al., 2015).

Alternatively, researchers have applied geostatistical techniques to interpolate the spatial variability of soil physico-chemical properties and to predict the values at unsampled locations (e.g. Bhunia et al., 2018; Ghong et al., 2023). The techniques can also be used to quantify spatial characteristics of agricultural lands so that the different landscape services (e.g. food supply, clean water and beautiful sceneries) can be mapped and sustainably managed (Aredehey et al., 2018; Hilal et al., 2024). The description of how spatial relationships change with the distance of separation between sampling points (lag distances) are accomplished through variogram analysis, which decomposes the spatial components of the variables into structural components that are used for interpolation. Ordinary kriging (OK) and cokriging (CK) are the interpolation techniques widely used to estimate spatial variability of soil physico-chemical properties (e.g. Delbari et al., 2019; Mousavi et al., 2017; Pirestani et al., 2024). Ordinary kriging is a geostatistical interpolation technique that uses statistical properties of sampled points (Isaaks & Srivastava, 1989), whereas CK uses one or more secondary variables to account for spatial cross-correlation between the primary and secondary variables (Goovaerts, 1998). As a result, CK can produce much more accurate and reliable estimation results than OK (Belkhiri et al., 2020; Khosravi & Balyani, 2019).

Soil physico-chemical properties like most environmental data often contain highly skewed values. In which case, applying indicator kriging (IK), a nonparametric type of conditional kriging and OK, may be a more appropriate geostatistical technique (Antunes & Albuquerque, 2013; Bradaï et al., 2016). As a result, IK model shows the quantified probability (0–1) of a variable that may exceed or not exceed a specified threshold value. That is, the expected value at an unsampled location derived from indicator data is equivalent to the cumulative distribution of the variable. Consequently, IK has been used to produce the probability maps of estimate values that exceed given limits such as regulatory threshold values (Belkhiri et al., 2017; Liang et al., 2018).

The main objectives of the work reported here were as follows: 1) to evaluate spatial variability of selected physico-chemical properties of surface soil samples in an area largely affected by cultivation and artisanal gold mining, 2) to determine the plausible sources of the spatial variability and 3) to identify areas with high environmental risks. The study is based on the hypothesis that land-use patterns such as farming and artisanal gold mining are the major driving factors of spatial variability of soil environmental quality. We have thought that a combination of geostatistical techniques (i.e. OK, CK and IK) will be appropriate for assessing spatial variability of the soil physico-chemical properties. The findings of the present study can help environmental geochemists in predicting behavior and fate of inorganic pollutants in highly anthropized semi-arid tropical soils at a relatively low cost.

2. Materials and methods

2.1. Study site and fieldwork

The study area (~42 km²) is located in southwestern Burkina Faso within the Bougouriba River drainage basin (Figure 1b). With a slop of 0.2 m/km and total length of 406 km, the Bougouriba River is the main watercourse in the area (Moniod et al., 1977). The area belongs to the South Sudanian climatic zone with predominant savannah vegetation. The regional annual rainfall varies between 900 and 2300 mm, whereas the average annual temperature is between 20.8 and 34.7°C (Forkuor et al., 2017). As far as the geology is concerned, the area is underlain by volcano-sedimentary rocks composed of andesite, tonalite and basalt with granitoid intrusion (Figure 1c). Under warm and humid climatic conditions, these geological formations are highly weathered, producing acidic soils with kaolinite as the dominant clay type. Hence, about 62% of the soil in the study area is classified as leached tropical ferruginous, 28% as latosols and 10% as hydromorphic soils with low SOM content (Hounkpatin et al., 2018). Artisanal gold mining is the main economic activity in the area. This unregulated activity includes explosive excavation and ore crushing, panning and mercury amalgamation, which have dramatic impacts on environmental soil quality. Subsistence (maize, sorghum, and yam) and cash crop (cotton) agricultures are the second largest land-use practice (Figure 1b).

Figure 1. Maps showing the study area (a), land occupation (b) and lithology with 132 topsoil sampling points from southwestern Burkina Faso (c).

In the present study, a total of 132 surface soil (0–30 cm) samples including 17 termite mounds were randomly collected using a stainless-steel scoop and stored in clean unused polyethylene bags. Back in the laboratory, the samples were dried at room temperature and sieved to 2 mm.

2.2. Soil physico-chemical property analyses

Both pH_H2O and EC of the samples were measured in 1 g soil to 2.5 mL distilled water using a calibrated electronic multimeter, whereas CEC was determined by the neutral ammonium acetate method (Bain & Smith, 1987). Particle size distribution of the samples was characterized through sieving for silt and sand, and through hydrometric method for clay content. Soil organic carbon content was determined using the Walkley and Black dichromate oxidation method (Nelson & Sommers, 1982). The method consisted of adding concentrated H₂SO₄ to a mixture of soil and aqueous K₂Cr₂O₇. After heating and oxidizing, the residual dichromate is back-titrated using ferrous sulfate (FeSO). The difference in added FeSO₄ is compared with a blank titration and determines the amount of easily oxidizable OC. The percentage of OC was calculated as follows (Eq. 1):

(1) $\mathit{OC}(\mathrm{\%})=\mathit{M}\times \frac{\left(V_1-V_2\right)}{m_s}\times 0.3\times \mathit{CF}\times 100$(1)

Where M is the molarity of FeSO₄ solution of blank, V₁ is the volume (mL) of FeSO₄ required in blank titration, V₂ is the volume (mL) of FeSO₄ required in sample titration, m_s (g) is the mass of the soil samples and CF is the correction factor used to account for carbon that does not oxidize in the procedure. The soil organic matter content was then obtained by multiplying the percentage of OC by the Bemmelen factor of 1.724 according to the following equation (Eq 2; Jackson, 1958):

(2) $\mathit{SOM}\mathit{\hspace{0.17em}}\left(\mathrm{\%}\right)\mathit{\hspace{0.17em}}=\mathit{\hspace{0.17em}}\mathit{OC}\mathit{\hspace{0.17em}}\left(\mathrm{\%}\right)\mathit{\hspace{0.17em}}\times \mathit{\hspace{0.17em}}1.724$(2)

The total nitrogen content (TN) was determined using the Kjeldahl method (Bremner & Sparks, 1996). A sub-sample of each sample was digested with concentrated H₂SO₄ by using CuSO₄ as a catalyst through which any nitrogen in the soil sample was converted into ammonia borate. Total nitrogen is then estimated by titrating ammonia borate with hydrochloric acid. Internal soil samples with known OC and TN content were analyzed along with the samples, and the recovery was better than 90%. Furthermore, three replicates were carried out for each analysis, and the relative standard deviations were less than 10%.

2.3. Exploratory data analysis

Soil physico-chemical property data were explored through a series of exploratory data analysis techniques. Thus, arithmetic mean and median were used to express the central tendency, whereas the spread of the data was evaluated through standard deviation (SD), maximum and minimum values. Traditionally, coefficient of variation (CV) has been used to evaluate the overall variability of soil physico-chemical properties (Castrignanò et al., 2000). The CV of a soil property is considered low if it is less than 15%, moderate if it is between 15 and 30% and high if it is greater than 35% (Wilding & Drees, 1978). Likewise, skewness and kurtosis were used to evaluate the data normality and asymmetry. Outliers and high skewness can affect negatively the spatial continuity of the variogram function. In order to provide stable variograms and kriging results, a transformation is considered as essential. The Pearson’s correlation analysis was used to investigate the relationship among the soil physico-chemical properties. The quantile-quantile (Q-Q) plots were used to test the normality of soil physico-chemical properties. In this study, pH, EC, TN and C/N did not show normal distribution, and they were log-transformed. All statistical analyses were performed using SPSS 20.0, PAST3 and Origin 2000 Pro.

2.4. Geostatistical analysis

2.4.1. Ordinary kriging and cokriging

Geostatistical techniques have been widely used to predict the spatial distribution of soil physico-chemical variables across unsampled areas based on the assumption that concentrations of a variable in closer samples are more spatially correlated than in those further away (Matheron, 1963). Kriging, in honor of Danie Krige, who first formulated this form of interpolation in 1951, has been extensively used as an important interpolation method in environmental studies. In the present study, OK was used to produce predictive maps of soil physico-chemical properties for unsampled locations. Some authorized theoretical models such as linear, exponential, spherical and Gaussian are usually used to fit the experimental semivariograms (Webster & Oliver, 2007). The semivariogram is expressed by using the Matheron’s classical estimator (Eq. 3).

(3) $\mathit{\gamma }(\mathit{h})=\frac{1}{2\mathit{N}(\mathit{h})}\sum _{\mathit{i}=1}^{\mathit{N}(\mathit{h})}[\mathit{z}(\mathit{x})-\mathit{z}(x_i+\mathit{h}){]}^2$(3)

Where Z(x) is the value of Z at location x_i and Z(x_i + h) is the value of Z at a location apart from xi by distance h and N(h) is the number of pairs of separated by h. Based on the exploratory data analysis, non-skewed data (less than 0.5), do not require transformation, whereas those with skewness value greater than 0.5 must be transformed before variogram analysis. There was no obviously anisotropy in semivariograms of soil physico-chemical properties in four different direction classes of 45°each. As a consequence, isotropic models were adjusted to fit to the data set of each individual parameter and to auto-fit one of the models (i.e. spherical, exponential, Gaussian or linear), and then select the best matching one.

The spatial structure of each model can be assessed using three variograms parameters: 1) nugget effect (C₀), which is the value of the semivariogram when the distance is 0. The nugget describes the spatial variation caused by either natural or anthropogenic factors. 2) Sill (C₀ + C), which is the value approaching the data variance, and it reflects the intensity of the total variation of regionalized variables, and higher values indicate higher degrees of variation. 3) The range (a) reflects the spatial autocorrelation of a regionalized variable and larger values indicate weaker spatial variability of the regionalized variable and stronger spatial correlation (Kashem & Singh, 2001; Matheron, 1963). The nugget to sill ratio also known as spatial dependence ratio (SDR) that is expressed as a percentage ((C₀/C₀ + C) ×100) was calculated for individual soil physico-chemical properties (Rosemary et al., 2017). A SDR value less than 25% indicates a strong spatial dependence, whereas values between 25% and 75% and higher than 75% indicate moderate and weak spatial dependences, respectively (Cambardella et al., 1994). Ordinary kriging was extended to cokriging when a strong correlation existed between a primary variable and a secondary variable.

2.4.2. Indicator kriging

Indicator kriging estimates the probability of exceeding specific threshold values, z_k, at a given location (Goovaerts, 1997). In the geostatistical techniques, the probability of attribute value is not greater than a specific critical limit (Zk) was estimated using IK (Goovaerts, 1997). The threshold value or critical limit is selected arbitrary based on the case study at hand. In indicator kriging, the data, z(x), are transformed into an indicator variable as follows (Eq. 4):

(4) $\mathit{i}(\mathit{x},z_k)=\mathit{f}(\mathit{x})=\left\{\begin{array}{c}1,\mathit{ifz}\left(\mathit{x}\right)\le z_k\\ 0,\mathit{otherwise}\end{array}\right.$(4)

At an unsampled location, the probability that z(x) ≤z_k can be estimated using a linear combination of neighboring indicator variables. In IK, data values are transformed into indicator values. The transformation of the original data on values of probabilities was carried out as follows: the values lower than the median value were assigned 0 values and values equal or greater than the median value was transformed into 1. That is, the original values which exceed the selected threshold value are coded 1, and those below the threshold value are coded 0.

Then, the semivariogram γ_i (h) is used to quantify the spatial correlation of the indicator codes, i (x_i, z_k), and it is written as follows (Eq. 5):

(5) ${\gamma }_i(\mathit{h})=\frac{1}{2\mathit{N}}\sum _{\mathit{i}=1}^{\mathit{N}(\mathit{h})}{\left[i(x_i;z_k)-\mathit{i}(x_i+\mathit{h};z_k\right]}^2$(5)

Where h is the distance between both locations, x_i and x_i + h, and N (h) is the number of pairs for x_i and x_i+h. The indicator kriging estimator, i (x_i; z_k) at the location x₀ can be calculated by (Eq. 6):

(6) $\mathit{i}\ast (x_0;z_k)=\sum _{\mathit{i}=1}^N{\lambda }_i\mathit{i}(x_i;z_k)$(6)

And the indicator kriging system is given the sum of λ_i as follows (Eq. 7)

(7) $\sum _{\mathit{i}=1}^N{\lambda }_i(x_j-x_i)={\gamma }_i(x_0-x_i)-\mathit{\mu }$(7)

Where λ_i is the weight coefficient, γ_i is the semi-variance of the indicator codes at the respective lag distance, and µ is the Lagrange multiplier.

The recommended international references for SOM and CEC contents (ITA, 1977; LCA, 2008) above which soils are likely to retain cationic nutrients and heavy metals (2% and 25 cmol/kg for SOM and CEC, respectively) were used to indicator variograms and probability maps of sensitive areas. That threshold values of pH (pH < 6.5) at which most cationic heavy metals (e.g. Pb, Cr and Cd) and oxyanion-forming elements such as Se, As and Sb (pH > 7.4) are mobile in soils are also used to delineate area at risk of heavy metal pollution.

2.4.3. Cross-validation

Accuracy of the kriging was evaluated through the leave-one-out cross-validation technique (Kerry et al., 2012). This technique is primarily used in scenarios where the objective is to predict accurately how a model will perform in practice. The cross-validation was conducted by visual inspection and optimization semivariogram parameters such as lag size (400) and lag number (12). The mean error (ME), mean standardized error (MSE), root mean square error (RMSE), root mean square standardized error (RMSSE) and average standard error (ASE) were the five validation parameters used to determine the best fitting models. ME values close to zero indicates unbiased predictions of soil physico-chemical properties, whereas any deviation of prediction values from measured values is evaluated by RMSE. An RMSE value close to that of ASE indicates that the assessment of prediction errors has been properly performed (Erdogan Erten et al., 2022; Nikroo et al., 2010). However, if RMSE is less than ASE, the variability of prediction errors is overestimated and the reverse case indicates an underestimation of error variability. The value of RMSSE should be close to 1. In the case where RMSE is greater than 1, the variability of prediction errors is underestimated, whereas a RMSSE value less than 1, then the variability of errors is overestimated (Ghazi et al., 2014; Nas & Berktay, 2010). When more than one model could be chosen for a variable, the choice was based the Akaike Information Criterion (Akaike, 1973). Geostatistical Analyst Extension within ArcMap 10.8 software was used to interpolate the spatial variability of soil physico-chemical properties and to produce the maps. For a given variable, all kriged maps were kept on the same scale in order to allow easier comparison.

3. Results and discussion

3.1. General soil physico-chemical properties

The spatial variability of soil physico-chemical properties is the result of intrinsic and extrinsic factors such as geological and pedogenic processes, fertilizer application, livestock grazing and mining (Bitencourt et al., 2016; Eljebri et al., 2018). Among the 10 soil physico-chemical properties used in this study, OC, SOM, TN, CEC, clay and silt data had skewness values less than 0.5 (Table 1; Figure 2), whereas EC, pH, C/N and sand were positively skewed (0.60–3). On the normal Q-Q plots, skewed variables deviated considerably from the normal distribution pattern (Figure 2a,b,f). After log-transformation, EC, pH, C/N and sand were normally distributed. In overall, soil physico-chemical properties exhibited high variabilities as highlighted by CV values (11–123%). However, there was a remarkable difference between physical properties and chemical properties. That is, the soil chemical properties, except pH (CV = 11%), exhibited relatively higher variabilities compared to physical properties (i.e. sand, silt and clay content). These high variabilities could be attributed to the local geology, residual fertilizers and vegetation cover (Ayele et al., 2019; Bogunovic et al., 2017). Nevertheless, the CV values of sand, silt and clay were greater than 30%, suggesting their high variabilities. This could be an indication of degradation of soil structure through tillage, deforestation and widespread artisanal gold mining in the study area (Bouasria et al., 2021; Eludoyin et al., 2017). With the application of acidifying fertilizers, intensive agriculture is also likely to impact other soil properties such as soil pH, EC and SOM content.

Figure 2. Normal Q-Q plots of soil physico-chemical properties collected in southwestern Burkina Faso: a) pH, b) EC, c) OC, d) SOM, e) N, f) C/N, g) CEC, h) clay, i) silt and j) sand.

Table 1. Univariate statistical summary of physico-chemical properties of topsoil samples from southwestern Burkina Faso

Variable	N	Unit	Min.	Max.	Mean	SD	Median	Skewness	Kurtosis	CV (%)
pH	132		5.1	9.6	6.7	0.7	6.6	1.2	3.5	11
EC	132	µS/cm	9.90	625.00	85.02	104.49	48.60	3.14	11.0	123
OC	132	%	0.19	1.93	0.99	0.36	1.03	−0.04	−0.43	36
SOM	132	%	0.33	3.33	1.71	0.62	1.77	−0.04	−0.43	36
TN	132	%	0.024	0.212	0.089	0.037	0.091	0.40	0.10	41
C/N	132		5.93	17.98	11.47	2.08	11.21	1.11	1.76	18
CEC	132	cmol/kg	4.92	38.58	20.30	7.70	20.92	−0.09	−0.71	38
Clay	132	%	6.59	59.53	31.68	11.31	33.06	−0.23	−0.71	36
Silt	132	%	9.80	51.96	27.66	8.98	28.16	0.23	−0.35	32
Sand	132	%	9.10	79.69	40.66	14.94	38.51	0.61	−0.06	37

SD: standard deviation and CV: coefficient of variation = (SD/Mean) x100.

The soil pH in the present study varied from 5.1 to 9.7 (6.7 ± 0.7) with the majority of the samples (up to 75%) having pH less than 7, indicating that the soils were mostly acidic to circumneutral (Penn & Camberato, 2019). The low soil pH values are generally inherent to tropical soils that are the byproducts of intensive chemical weathering, depletion of alkali and alkaline earth elements (Latifah et al., 2017; Russell et al., 2017). The absence of these elements in the soil is corroborated by the relatively low soil EC (Rangel-Peraza et al., 2017). Thus, the recorded soil EC ranged from 9 to 625 µS/cm with an average value of 85.0 ± 102.5 µS/cm. According to Richards’s salinity classification (Richards, 1954), the measured EC values were low with a large difference between the median and average values (85.02 µS/cm versus 48.60 µS/cm). However, a moderate positive correlation (r = 0.53; p < 0.05) between pH and EC (Figure 3) suggests that soil pH controls ion abundance in soil solution (Penn & Camberato, 2019). The low pH values may also affect the rates of biological activities, and thereby soil C/N ratios.

Figure 3. Pearson correlation heatmap of soil physico-chemical properties at p = 0.05.

Overall, the soil C/N ratios ranged from 5.9 to 17.98 with 75% of the samples having C/N ratios less than 12. According to the LCA guidelines (LCA, 2008), C/N ratios comprised between 8 to 12 indicate good organic matter decomposition, with ratios greater than 12 imply that activities of soil microbial community is reduced and the organic matter decomposition rate is slow. The areas with high C/N ratios (~25%) greater than 12 are characterized by low biological activities, which could be attributed to low dissolved oxygen content and pH (Z. Xu et al., 2020; Zhao et al., 2018). The SOM content of the present study ranged from 0.33 to 3.3% (1.7 ± 0.6%) with 50% of soil samples having moderate SOM content (~2%; ITA, 1977). This is an indication of good biological activities and organic matter decomposition rates. Soil organic matter plays a prominent role in mobility and biogeochemical cycling of heavy metals in the soil environment (Kwiatkowska-Malina, 2018; Zhu et al., 2018). Thus, in the presence of high SOM content in the upper horizon, certain heavy metals, particularly Pb, tend to be accumulated in the soil (Y. Li et al., 2017). Consequently, SOM could be a limiting factor in immobilizing heavy metals in the present soils. As expected, the soil OC contents were relatively lower than that of SOM, and it ranged from 0.19 to 1.93% (0.99 ± 0.3%). About 75% of the soil samples had low TN content (<0.11%; Batjes, 2014). The highly positive correlations (Figure 3), observed in the present work, between SOM and OC contents and between OC and TN content suggests that OC and TN contents are largely controlled by the supply of readily decomposable organic materials (Bongiorno et al., 2019). Furthermore, the significant correlation between OC and TN contents is an indication of soil C/N stabilization in the plow layer (Kopittke et al., 2020; Y. Wang et al., 2018). Similar close relationships have been observed in terrestrial ecosystems with different land-use or management practices (Ferreira et al., 2016; T. Wang et al., 2016).

On the textural triangle, the majority of the soil samples was projected in the clayey and silty-clayey textural classes (Figure 4). Under intense and long-term chemical weathering of crystalline bedrocks, one could expect a predominance of sand and residual Fe and Al oxide minerals (Souza et al., 2018). However, the observed textural classes could be a result of tillage, water and erosional redistribution of fine-grained particles from nearby upland areas to deposition areas (Figure 6a) in the study area (Berhe & Torn, 2016; Nie et al., 2016). The soil CEC, which varied broadly across the study area (4.9–38.6 cmol/kg), showed a strong positive correlation with clay (r = 0.94, p < 0.05). This indicates that CEC is solely related to the local clay type (Pincus et al., 2017; Ramos et al., 2018). The recorded CEC values fall within the range of kaolinite (3–57 cmol/kg; Kabata-Pendias, 2011).

Figure 4. Soil textural triangle of the samples based on the International Union of soil Sciences (IUSS Working Group, 2015).

Figure 6. Spatial distribution maps of soil physico-chemical properties collected in southwestern 884 Burkina Faso derived ordinary kriging. a) elevation map, b) pH, c) EC, d) OC, e) SOM, f) TN,885 g) C/N, h) CEC, i) clay, j) silt and k) sand. pH, EC, C/N and sand were log-transformed.

3.2. Structural characteristics of soil physico-chemical properties

Isotropic semivariogram models were fitted to the raw and log-transformed soil physico-chemical property data (Figure 5). That is, the Gaussian model was best-fitted for pH, OC, C/N, CEC and sand, whereas exponential model was more adequate for EC, SOM and clay, and the spherical model for TN and silt. For all samples, RMSE values were close to 1 (1.00–0.96) and MSE values were low (−0.009742–0.0174), which indicated that the fitted variograms and subsequent kriging standard errors were satisfactory (Table 2). Thus, OK could successfully interpolate pH, OC, SOM, C/N, CEC and silt with minimum RMSE (0.04–9.3) and ASE (0.04–9.3) values, whereas EC, clay and sand showed larger uncertainty with maximum RMSE (98, 11 and 15) and ASE (139, 11 and 17) values. Furthermore, the low ME values, except for EC, indicated a good match between observed and predicted soil physico-chemical properties.

Figure 5. Experimental semivariograms for a) log-transformed pH, b) log-transformed EC, c) OC, d) SOM, e) TN, f) log-transformed C/N, g) CEC, h) clay, i) silt and j) log-transformed sand based on cross validation analysis.

Table 2. Best-fitted variogram parameters and cross-validation prediction errors for ordinary kriging (OK) of soil physico-chemical properties

Variable	Method	Model	Nugget effect	Sill	Range	SDR (%)	ME	MSE	RMSSE	RMSE	ASE	Spatial dependance
			C0	C + C0	(m)
logpH	OK	Gaussian	0.007	0.011	3950	63	0.0040	0.0038	1.00	0.6	0.58	Moderate
logEC	OK	Exponential	0.458	0.827	6009	55	−4.7293	−0.1198	1.58	98.5	80.69	Moderate
OC	OK	Gaussian	0.105	0.130	3019	81	0.0055	0.0157	1.00	0.34	0.34	Weak
SOM	OK	Exponential	0.293	0.392	4370	75	0.0061	0.0103	1.00	0.59	0.59	Moderate
TN	OK	Spherical	0.001	0.001	2804	88	0.0006	0.0174	0.98	0.04	0.04	Weak
logC/N	OK	Gaussian	0.026	0.030	2502	88	0.0070	0.0001	1.00	1.94	1.95	Weak
CEC	OK	Gaussian	14.318	60.272	530	24	0.0005	−0.0042	0.98	7.3	7.4	Strong
Clay	OK	Exponential	96.763	119.27	1440	81	0.0287	0.0009	0.99	10.97	11.05	Weak
Silt	OK	Spherical	75.674	79.044	1032	96	0.0003	−0.0009	0.96	8.89	9.33	Weak
logSand	OK	Gaussian	0.099	0.151	796	65	−0.0443	−0.0409	0.96	15.29	16.99	Moderate

SDR = dependence ratio ((C₀/C+C₀) ×100); ME= mean error; MSE= mean standardized error; RMSSE= root mean square standardized error.

RMSE= root mean square error; ASE= average standard error.

Five variables (pH, TN, C/N, OC and sand) showed very low nugget effects (0.007–0.105). In contrast, clay, silt and CEC (96, 75 and 14.3) showed the highest nugget effects, indicating a random spatial distribution behavior of the variables at short range (Song et al., 2018). The large nugget effect suggests that the lag distance did not characterize the spatial variations, and thus an additional sampling of the variable at smaller lag distances and in larger number is required. Another factor that may explain the large nugget effect is the subscale variability of the variable (Guedes et al., 2020). In other words, the random distribution may be attributed to the presence of extremely high or extremely low values that are outliers or to sampling design.

The range was used to reflect the distance of spatial correlation for the soil physico-chemical properties (Bitencourt et al., 2016). Electrical conductivity, TN and clay were autocorrelated in longer ranges than pH, OC, SOM and silt. This is consistent with the high CV values of EC, TN and clay. The SDR of the individual soil properties were in the following decreasing order (Table 2): The SDR of the individual soil properties were in the following decreasing order (Table 2): silt (96%)> TN=C/N(88%)>OC=clay (81%)>SOM (75%)>sand (65%)> pH (63%)>EC(55%)> CEC (24). All the ratios, except that of CEC, were greater than 25%, implying that the spatial dependence was weak to moderate. With SDR values greater than 75%, silt, TN, C/N, clay and OC had weak spatial dependence, implying that the spatial distribution of these soil properties was only controlled by extrinsic factors such as biological activity and land-use pattern (Eljebri et al., 2018; Phefadu & Kutu, 2016). Conversely, SDR values of EC, SOM, sand and pH comprised between 25 and 75%, and were therefore considered to be moderately dependent. The spatial distribution of these chemical properties can be attributed to both intrinsic and extrinsic factors. Thus, EC, pH and sand had SDR slightly over 50%, which suggests that the extrinsic factors, such as tillage, vegetation cover and biological activities were the most controlling factors in the spatial distribution of these variables. Soil CEC was the only soil property with strong spatial dependency in the study area. This implies that the soil CEC might have been controlled by intrinsic factors such as the dominant clay type (i.e. kaolinite).

3.3. Spatial distribution of soil physico-chemical properties

Semivariogram parameters were used to produce spatial distribution maps of the individual soil physico-chemical properties through OK (Figure 6). In the study area, about 87% of the total surface area is covered by slightly acidic to circumneutral soils (pH = 6.1–7.3; Figure 6b) in farming dominant northern and eastern parts, whereas neutral to alkaline soils (pH = 7.4–8.7) occupied about 7% of the area n northwest, and acidic soils (~6%) were found mainly in the southeast. In addition to anthropogenic activities, the spatial distribution patterns of pH could be partially attributed to the local lithology, low inputs of fresh organic material, removal of basic nutrients and continuous use of ammonium fertilizers (Khadka et al., 2018; Khan et al., 2019).

Figure 6. (Continued).

Spatial distribution of EC varied widely across the study area with the high values (69–167 µS/cm) confined in northwest, and mainly occupied by farmlands, whereas the rest of the study area was characterized by low EC values (17–69 µS/cm; Figure 6c). The pedogeochemical signature of these soils could be directly related to the prevailing agricultural practices in the area. Excessive application of inorganic fertilizers could considerably increase soil background concentrations in certain ions. Thus, the underlying causes of the spatial distribution of EC are extrinsic (e.g. land-use) factors.

The soil samples with moderate SOM content were collected in the central and southeastern parts of the study area (Figure 6e), whereas the soil samples with low SOM content were mainly collected in farmland. Clearly, the spatial distribution pattern of soil SOM was primarily controlled by farming. It is noteworthy that soils samples with moderate TN content were collected in areas characterized by natural vegetation cover, whereas those with low TN content were located in farmlands (Figure 6f), indicating the influence of tillage on the soil organic matter content. Moreover, soil samples with moderate to high C/N ratios (<12; LCA, 2008) were spread around artisanal gold mining sites and farmland (Figure 6g), whereas those with very high ratios (>12) were mainly collected around human settlements and the textile factory. Therefore, it can be concluded that biological activities and organic matter decomposition rates were directly associated with land-use patterns. Again, this indicates that particular attention should be paid to farmland and artisanal gold mining sites as far as inorganic pollutant mobility is concerned.

Patches of soils with high CEC content (25–37.4 cmol/kg) were scattered throughout the study area (Figure 6h). Moreover, the clay distribution showed a significant spatial correlation with CEC (Figure 6h,i), which is consistent with previous studies that indicated that soil CEC is directly dependent on clay types. Soils with a high clay content (35–43%) were collected around the artisanal gold mining sites and along the river bank (Figure 6i). Based on the significant positive correlations between CEC and clay and between SOM, TN and OC, the spatial distribution of CEC and SOM was mapped using CK. Soil CEC and SOM were used as primary variables, whereas clay, OC and TN as secondary variables. Nonetheless, the ME and RMSE values of OK for SOM were lower than those of CK, suggesting that CK was more accurate and reliable in predicting SOM. In contrast, ME and RMSE values of OK for CEC were lower than those of CK. Hence, in the present study CK was less accurate in estimating CEC and clay (Table 3; Figure 7a,b). Therefore, there is no need to determine TN and OC in future monitoring programs as SOM can be used instead.

Table 3. Best-fitted variogram parameters of ordinary kriging (OK) and cokriging (CK) of CEC and SOM as primary variables and clay. TN and OC as secondary variables

	Method	Model	ME	MSE	RMSSE	RMSE	ASE
SOM	OK	Exponential	0.0061	0.0103	1.00	0.59	0.59
SOM	CK	Gaussien	0.0041	0.0062	0.97	0.58	0.59
CEC	OK	Gaussian	0.0005	−0.0042	0.98	7.3	7.4
CEC	CK	Exponential	0.0339	0.0093	1.09	2.69	2.60

Figure 7. Spatial distribution maps of a) CEC with clay as auxiliary variable and b) SOM with OC and TN as auxiliary variables based cokriging.

3.4. Delineation of environmentally high risk areas

Identification of areas at high risk of heavy metal pollution is urgently required for a thorough assessment of their behavior and fate in highly impacted rural soils. Soil physico-chemical properties (e.g. pH, CEC and SOM) play a critical role in mobility and bioavailability of heavy metals in the soil environment (Antoniadis et al., 2017; Gao et al., 2016). Several studies have reported that pH contributes to metal mobility from adsorption sites to soil solution and their subsequent uptake by organisms (Jović et al., 2017; Pikuła & Stępień, 2021; D. M. Xu et al., 2021). This is likely due to the twofold effects of pH on mobility and speciation of most heavy metals (Jarvela-Reijonen et al., 2016; Korfali & Karaki, 2018). Estimated probability maps of soils physico-chemical properties were produced by IK using published threshold values (LCA, 2008). A Map of land occupation was used to highlight spatial distribution the different soil physico-chemical properties at risk of heavy metal mobility and toxicity (Figure 8a).

Figure 8. Estimated probability maps of soil physico-chemical properties indicating areas at risk of exceeding or not exceeding environmental threshold values with relation to farmlands. a) land occupation, b) pH < 6.5, c) pH > 7.4, d) CEC < 25 cmol/kg and e) SOM < 2 %.

About 29% of the sampled area, particularly, in northern and eastern parts had high probability (50–100%) of not exceeding the pH value (<6.5) at which most heavy metals are not mobile and toxic (Figure 8b; Table 4; Sintorini et al., 2021). By contrast, about 2% of the soils in the northwestern part (Figure 8c) of the study had high probability (50–75%) of exceeding the threshold value (pH = 7.4) above which the mobility of most oxyanion-forming elements such as Cr, Se and As can be enhanced (Kodirov et al., 2018; Verbeeck et al., 2019). This suggests that about 31% of the total areas was at risk of heavy metal mobility (metallic and metalloid elements).

Table 4. Percentage areas at risk of not exceeding (SOM. CEC and pH < 6.5) or exceeding (pH > 7.4) threshold values with respect to the total surface of the study area

Probability (%)	pH < 6.5	pH > 7.4	CEC <25 cmol/kg	SOM <2%
0–25	28.16	83.14	3.67	2.44
25–50	43.01	14.43	19.12	21.95
50–75	19.18	2.31	33.18	43.9
75–100	9.56	-	43.91	31.71

As far as CEC is concerned, 77% of the total area was unlikely (50–100%) to have CEC values above 25 cmol/kg (Figure 8d), suggesting that soils in these areas are characterized by low cation retention capacity. The adsorption capacity of heavy metals onto soils is a function of surface charges, ionic strength and clay mineralogy (Anjolaiya, 2015; Sarkar et al., 2017; Uddin, 2017). In contrast to temperate soils with dominant illite (2:1 clay types) and permanent structural charges, highly weathered tropical soils with dominant kaolinite (1:1 clay types) are characterized by pH-variable surface charges (Chen et al., 2017; Colzato et al., 2018; Singh et al., 2016). That is, the adsorptive capacity of soils in areas with relatively high CEC (>25 cmol/kg) may be reduced at small declines in pH.

Similarly, 75.61% of the study area had high probability (50–100%) of not exceeding the SOM threshold (<2%; Figure 8e). These soils are therefore considered to be relatively poor in SOM and would not be capable of retaining cationic elements including potentially toxic heavy metals (M. Huang et al., 2016; Lwin et al., 2018). In that sense, protonation (high H⁺) and deprotonation (high OH⁻) of soil surfaces are greatly influenced by SOM content as organic matter balances both acid and base addition to the soil solution. Less than 1% of the total surface area had low probability (0–25%) to be below the threshold value of 2% that will allow heavy metal retention. Soil organic matter is deprotonated at high pH and therefore can adsorb cationic heavy metals, whereas at low pH, SOM is protonated and previously adsorbed cations will be released into soil solution. This protonation will pose serious bioenvironmental concerns as previously retained heavy metals may be released through desorption into the soil solution (M. Huang et al., 2019). Based on low pH, low CEC low and SOM content of large portion of the study area, it can be deduced that the studied soils would be vulnerable to anthropogenically derived heavy metals.

Since agricultural soils are characterized by high physical and chemical properties, homogenous management of agricultural fields may lead to overapplication of fertilizers in high-yielding areas and underapplication of fertilizers in low-yielding areas (Bhatt et al., 2018; Q. Li et al., 2021; Nadarajan & Sukumaran, 2021). After establishing spatial dependence through geostatistical modeling, isarithmic maps for spatial distribution of CEC, pH, SOM and clay can be generated. These maps can be used to categorize areas with low, medium and high physico-chemical properties for site-specific soil fertility management (Silva et al., 2023). Not only quantifying spatial variability of physico-chemical properties will increase agricultural yield, but it will also allow better management of natural resources and environmental pollution control. That is, the site-specific field management minimize landscape degradation, eutrophication, soil and water resource contamination through excess inorganic fertilizers.

4. Conclusions

Exploratory data analysis showed that soil chemical properties, except pH, were more variable compared to physical properties. This is due to both intrinsic and extrinsic factors such lithology, chemical weathering, vegetation cover, biological activities and land-use pattern. The pronounced low pH across the study area is likely associated with lithology and pedogenic process and agricultural practices. Due to the crystalline nature of the local lithology, one could expect dominant sandy soil rather than the observed clayey and silty soils. This discrepancy reflects transport and deposition of fine-grained particle from human-impacted agricultural and mining areas to lowlands along the river bank. Cokriging was the best method for modeling spatial variability of SOM, OC and TN. Indicator kriging showed that most of the soil samples had probability of exceeding threshold values for pH, CEC and SOM at which most cationic heavy metals will be mobile and bioavailable. The study also showed that through a combination of geostatistical approaches, the routine soil physico-chemical properties can be used to predict areas at risk of heavy metal pollution so that robust pollution control strategies can be set up.

Availability of data

The data generated and/or analyzed to support the findings of the study is available from the corresponding author, Sako, A., up on reasonable request.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The work was supported by the Higher Education doctoral fellowship (PAES) of Burkina Faso Higher Education doctoral fellowship. K. Coulibaly was supported by the Burkina Faso’s Higher Education doctoral fellowship (PEAS). Valuable comments and suggestions from tow anonymous reviewers and the Editor-in- chief of the journal have greatly improved the early version of the paper.