Abstract

Soil fertility is a determining factor for agricultural productivity and food security. This study aims to map the spatial variability of soil fertility parameters (nitrogen, pH, and carbon) in the Sudan-Guinea and Guinea zones of Benin using geostatistical approaches. To achieve this objective, fifty-two (52) soil samples were collected using a stratified random sampling design and analyzed at the Soil, Water and Environmental Science Laboratory (LSSEE) of the Agricultural Research Center based in Agonkanmey, Abomey-Calavi. Total nitrogen, organic carbon content, and water pH were determined using conventional methods. Four kriging methods (simple, ordinary, universal, and indicator) were compared after fitting three variogram models (spherical, Gaussian, and exponential). The results reveal that the spherical model provided the best fit to the data for nitrogen, pH, and carbon. Indicator kriging proved optimal for nitrogen (RMSE = 0.47), simple kriging for pH (RMSE = 0.40), and ordinary kriging for carbon (RMSE = 0.41). Spatial analysis indicates low nitrogen content (0.08% ± 0.03%), slightly acidic pH (6.63 ± 0.22), and carbon content (0.79% ± 0.17%). The carbon stock was estimated at 1.62 t∙C∙ha⁻¹, revealing a moderate sequestration potential. These results provide a scientific basis for differentiated soil fertility management in these agro-ecological zones.

Keywords

Geostatistics, Fertility, Nitrogen, Organic Carbon, Kriging, Benin

Share and Cite:

1. Introduction

Population growth and the intensification of food needs exacerbate the risks of food insecurity in developing countries, particularly in Benin. This situation results notably from the progressive degradation of agricultural soils due to inappropriate farming practices. Among these practices is shifting slash-and-burn cultivation, which destroys flora, organic matter, as well as soil fauna and microfauna. The excessive use of synthetic chemical fertilizers on nutrient-poor soils leads to a decline in the quality of agricultural products. These practices affect not only production but also the sustainability of soils and agricultural systems [1]. Soil is a medium rich in nutrients (micro and macro elements) whose properties vary in space and time [2]. It also constitutes a significant carbon reservoir. Soil represents one of the most important carbon sinks globally (alongside oceans and forests). By helping to reduce the amount of atmospheric CO, carbon sinks influence the planetary climate and thus all environmental components dependent on climate. The amount of organic carbon stored in the topsoil layer is estimated at 700 Gt (gigatonnes) worldwide [3]. According to [4], for decades, the most studied soil nutrients in precision agriculture have been nitrogen (N), phosphorus (P), and potassium (K), due to the importance of the roles they play. According to [5], nitrogen (N) promotes growth, endosperm and chlorophyll formation, while phosphorus (P) is useful in the transformation of energy substances, the formation of roots, flowers, and fruits, and potassium (K) for plant health and stress resistance. The spatial autocorrelation of soils means that two soil samples located in geographically closer locations tend to have more similar properties than samples considered sufficiently distant [6]. Knowledge of soil nutrient content is important in crop production as its variability affects crop yield as well as the nutritional value of plant species [7] [8]. Conventional soil fertility mapping approaches use categorization of soils at sampling points, assuming that variability within identified units is homogeneous [9]. Unlike conventional mapping approaches that assume spatial homogeneity, geostatistics makes it possible to quantify and model the spatial autocorrelation of soil properties. This approach offers increased accuracy in estimating values at unsampled points while providing a measure of the uncertainty associated with predictions [10] [11]. Geostatistics is a true statistical approach for processing spatial or spatio-temporal data, as it relies on both statistical and mathematical methods to minimize prediction uncertainties [12], thus providing more reliable estimated values [13]. In geostatistics, several methods have been developed, including kriging. The different existing kriging techniques include simple, ordinary, universal, and indicator kriging. For decades, these methods have been used variously in environmental studies to assess the spatial variability of water quality [14], rainfall [15], and atmospheric pollutants [16]. This study aims to map the spatial variability of soil fertility parameters (nitrogen, pH, and carbon) in the Sudan-Guinea and Guinea zones of Benin using geostatistical approaches for mitigating the effects of climate change.

2. Materials and Methods

2.1. Description of the Sudan-Guinea and Guinea zones

The study was conducted in the Sudan-Guinea and Guinea agro-ecological zones of Benin (Figure 1), covering a total area of approximately 48,347 km². The Sudan-Guinea zone, located between 7˚30'N and 9˚45'N, has an annual rainfall of 900 to 1100 mm. The annual temperature varies between 21.2˚C (average minimum) and 32.5˚C (average maximum). The rainfall pattern in the Sudan-Guinea zone is unimodal (May-October) and most often distributed over an average of 113 days. Relative humidity ranges from 31% to 98% in this zone. Average sunshine amounts to 2305 hours per year. This Sudan-Guinea transition zone, extending the zone with Guinean affinities, is the domain of open forest mosaics, possibly with dry dense forests, interspersed with wooded and shrubby savannas and traversed by gallery forests. This zone consists of weakly evolved, low-fertility mineral soils and ferruginous soils on a crystalline bedrock of variable fertility. In this zone, plants such as Daniellia oliveri, Parkia biglobosa, and Terminalia glaucescens are found on well-drained soils, Anogeissus leiocarpus, Acacia campylacantha, and Terminalia macroptera on hydromorphic soils, and Isoberlinia doka and Detarium microcarpum on soils over iron crusts or shallow rocks [17]-[19].

Figure 1. Map of the geographical location of the study areas.

As for the Guinea zone, it is located between 6˚30'-7˚30'N, 1˚40'-2˚45'E of Benin, and has an annual rainfall of 1200 mm. The annual temperature varies between 25˚C and 29˚C with a bimodal rainfall pattern and relative humidity ranging from 69% to 97%. The Guinea zone is the domain of coastal forests, mangroves, semi-deciduous forests, and semi-deciduous swamp forests. This zone consists of deep ferrallitic soils, rich in vertisols, humus, and minerals. In this zone, plants such as Syzygium guineense, Chrysobalanus icaco, Diospyros tricolor, Rhizophora racemosa, and Avicennia germinans are found on halomorphic and hydromorphic sandy soils. Conversely, on ferrallitic soils are found Strombosia pustulata, Chytranthus macrobotrys, Hannoa klaineana, Landolphia incerta, Pouchetia africana, Mansonia altissima, and Pterygota macrocarpa. Similarly, Xylopia rubescens, Mitragyna ledermannii, Uapaca paludosa, Cynometra vogelii, Grewia malacocarpa, and Grewia barombiensis are developed on hydromorphic soils [20].

2.2. Sampling

A stratified random sampling approach was adopted to ensure that the main sources of spatial variability in soil properties were adequately represented. The strata were defined based on agro-ecological zones and land use types, as these factors strongly influence soil nitrogen, pH, and organic carbon dynamics in the study area. Accordingly, the study area was divided into two main strata corresponding to the Sudan-Guinea and Guinea zones, which differ in rainfall regime, vegetation cover, and land-use intensity.

Within each stratum, sampling points were randomly distributed at a density of one point per 930 km². A total of 52 sampling points were established—30 in the Sudan-Guinea zone and 22 in the Guinea zone—and georeferenced using GPS to ensure spatial tracking and relocation.

At each point, a composite soil sample (≈500 g) was prepared by mixing five subsamples collected within a 10 m radius around the central point. Sampling was performed at a depth of 0 - 20 cm using a Dutch auger. All samples were placed in labeled polyethylene bags and transported to the laboratory for physicochemical analyses.

2.3. Determination of Physico-Chemical Parameters

Soil samples were first air-dried for 72 hours and then sieved through a 2 mm mesh to isolate the fine earth [21]. Total nitrogen was determined by the modified Kjeldahl method [22], organic carbon by the Walkley-Black method [23], and water pH using a pH meter with a soil-water ratio of 1:2.5 [24]. Equation (1) is the formula to calculate organic matter (OM):

$\text{MO}\left(\%\right)=\text{C}\left(\%\right)\times 1.724$ (1)

Available phosphorus was determined according to the method of [25], exchangeable cations by extraction with ammonium acetate at pH = 7, and potassium by flame photometry.

2.4. Determination of Soil Carbon Stock

To assess the organic carbon stock of the soil according to the area, three parameters were estimated to determine the amount of carbon on the horizon: the bulk density of the horizon, the concentration of organic carbon in a given soil mass for the horizon, and the thickness of the soil layer [26]-[29]. Equation (2) is for assessing the amount of carbon in the soil [30]:

$C=\left(A_h\times D{A}_h\times P_h\times \left(1-F_h\right)\right)\ast S$ (2)

where $C$ = soil organic carbon stock in the Sudan-Guinea and Guinea zones (expressed in ton C∙ha⁻¹); $A_h$ = concentration of organic carbon in a given soil mass for the horizon, obtained by laboratory analyses (g∙C∙kg⁻¹∙soil); $D{A}_h$ = Bulk Density of the horizon (ton∙soil∙m⁻³); $P_h$ = soil sampling depth or thickness of the soil layer (m); $F_h$ = percentage by volume of coarse fragments/100 and $S$ = area occupied by the soils.

Bulk density ($D{A}_h$ ) was determined using the core method, which consists of collecting undisturbed soil cores using a cylindrical metal ring of known volume (100 cm³) at each sampling depth. The samples were oven-dried at 105˚C for 48 hours until constant weight, and the bulk density was calculated as the ratio of the oven-dry mass to the volume of the core [31]. For soils where direct measurements were not possible, bulk density was estimated using pedotransfer functions based on soil texture and organic matter content following the procedure of [32].

2.5. Statistical Analyses

An exploratory data analysis was performed for the first time to characterize the variables related to the different measured parameters, according to the practices described by [33]. Descriptive statistics were calculated, including minimum, maximum, mean, standard deviation, coefficient of variation, skewness, and kurtosis coefficient [34]. Since geostatistical methods are sensitive to deviations from normality, a bias that can affect spatial prediction [35], the normality of the data was verified using the Shapiro-Wilk test (α = 0.05). These analyses focused on the variables nitrogen (N), pH, and soil carbon and were performed using R and ArcGIS 10.1 software.

2.5.1. Selection of the Variogram Model

Several variogram models were compared, following the recommendations of various authors [8] [36]-[39]. The models tested included spherical (Equation (5)), exponential (Equation (6)), and Gaussian (Equation (7)) structures [39] [40]. For each fertility parameter, three models were fitted and compared based on cross-validation criteria, notably the Mean Absolute Error (MAE, Equation (3)) and the Root Mean Square Error (RMSE, Equation (4)). The optimal model was selected according to the RMSE criterion, which is considered more robust for assessing the precision of geostatistical predictions [41]. The formulas used are as follows:

$\text{MAE}=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\left|Z\left(x_i\right)-Z^{*}\left(x_i\right)\right|$ (3)

$\text{RMSE}=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left(Z\left(x_i\right)-Z^{*}\left(x_i\right)\right)}^2}$ (4)

where $Z\left(x_i\right)$ is the observed value, $Z^{*}\left(x_i\right)$ the predicted value, and $n$ the number of observations.

The mathematical expressions of the variogram models used are:

Spherical: $g\left(h\right) =\{\begin{array}{l}{C}_o+C \left(\frac{3h}{2a} - \frac{h^3}{2a^3}\right),\text{ if }h\le a\\ C_0 + C,\text{ if }h>a\end{array}$ (5)

Exponential: $g\left(h\right) = C_0 + C \left(1-{\text{e}}^{-h/a}\right)$ (6)

Gaussian: $g\left(h\right) = C_0 + C \left(1-{\text{e}}^{-\left(\frac{h}{a}\right)}{}^2\right)$ (7)

with h: lag distance, C₀: nugget effect, C: sill, a: range.

2.5.2. Evaluation of Kriging Method Performance

Four kriging methods (simple, ordinary, universal, and indicator) were applied to the soil fertility parameters [34] [42] [43].

Ordinary Kriging (OK): Assumes a constant but unknown mean over the entire study area. This is the most commonly used method for environmental data (Equation (8)).

Model: $Z\left(x_0\right) =\mu +ϵ\left(x_0\right)$ (8)

Assumption: μ is an unknown constant.

Simple Kriging (SK): Corresponds to ordinary kriging with a known mean. It is slightly more powerful than ordinary kriging. However, the main difficulty lies in determining the mean, making the assumption of an exactly known mean unrealistic (Equation (9)).

Model: $Z\left(x_0\right) =\mu +ϵ\left(x_0\right)$ (9)

Assumption: μ is a known constant.

Universal Kriging (UK): Accounts for a trend or drift in the data by incorporating a deterministic model of spatial variation (Equation (10)) [43].

Model: $Z\left(x_0\right)=\mu \left(x_0\right)+ϵ\left(x_0\right)$ (10)

Assumption: $\mu \left(x_0\right)$ is a deterministic function.

Indicator Kriging (IK): This method is used when the data do not follow a normal distribution or when they contain extreme values. It transforms the data into binary variables (0 or 1) based on a specific threshold, allowing non-parametric modeling of the probability of exceeding this threshold (Equation (11)).

Model: $I\left(x_0\right)={\mu }_I+ϵ\left(x_0\right)$ (11)

Assumption: The indicator variable $I\left(x_0\right)$ follows a binomial distribution.

Their performance was compared using two assessment criteria: the mean absolute error (MAE) and the root mean square error (RMSE) [8]. A semi-variogram cloud was generated for each parameter to assess the spatial variability between neighboring observations. The methods were ranked according to the minimum values of MAE and RMSE, with optimal performance corresponding to a mean prediction error (ME) close to zero and a root mean square standardized error (RMSS) close to 1 [44].

3. Results

3.1. Statistical Description of the Studied Soil Characteristics

Table 1 presents the general trend of all parameters and indicates that out of 52 soil samples covering the two agro-ecological zones of Benin, nitrogen averaged 0.08% ± 0.03% and pH averaged 6.63 ± 0.22. Furthermore, the average carbon content in these regions was 0.79% ± 0.17%. The standard deviations around the parameters indicate that each exhibits low dispersion around the mean.

Table 1. General trend of parameters.

N = Nitrogen, pH = Hydrogen Potential, MO = Organic Matter, C = Carbon.

The positive skewness coefficients for N (0.56) and C (1.86) reveal right-skewed distributions, suggesting the presence of localized areas with high contents. The Shapiro-Wilk test (p < 0.05) revealed that none of the parameters followed a normal distribution. Therefore, several transformation methods belonging to the Box-Cox family were applied to the data. After comparing the results, the logarithmic transformation proved to be the most adequate for correcting the normality of all parameters with a lambda value equal to 0 (Table 2).

Table 2. Data distribution before and after logarithmic transformation.

N = Nitrogen, pH = Hydrogen Potential, C = Carbon, Sk = Skewness, Kurt = Kurtosis.

3.2. Spatial Variability of Soil Fertility

3.2.1. Selection of the Variogram Model

The fitting of three theoretical variogram models (spherical, Gaussian, and exponential) to the experimental data reveals differentiated behaviors depending on the studied parameters (Table 3). The comparative analysis of the two criteria (MAE and RMSE) reveals the superiority of the spherical model for all three studied parameters. For nitrogen, although the exponential model shows a slightly lower MAE (0.025 vs 0.030), the spherical model displays a clearly better RMSE (1.010 vs 1.060), indicating less dispersion of prediction errors. For pH, the spherical model clearly dominates with an RMSE of 0.938 compared to 0.952 for the exponential and 1.561 for the Gaussian, despite a slightly higher MAE (0.200 vs 0.152 for the Gaussian). Regarding carbon, the three models show relatively close performances, but the spherical model remains competitive with acceptable MAE (0.114) and RMSE (1.291) values. The predominance of the spherical model suggests a spatial dependence structure with a finite range beyond which observations become independent. This structure indicates that the variability of fertility parameters follows a spatial pattern with well-defined transition zones, characteristic of pedogenetic processes influenced by topography, geology, and localized agricultural practices in the Sudan-Guinea and Guinea zones of Benin.

Table 3. Choice of variogram model.

MAE: Mean Absolute Error, RMSE: Root Mean Square Error, N = Nitrogen, pH = Hydrogen Potential, C = Carbon.

3.2.2. Performance of Kriging Methods

The comparison of the four kriging techniques reveals a specificity of performance depending on the nature of the studied parameter (Table 4). This variability in method efficiency reflects differences in the statistical distribution and spatial structure of each variable. For nitrogen, indicator kriging (IK) stands out with an RMSS of 0.94 (close to the theoretical optimum of 1), although presenting a higher RMSE (0.47). This apparent contradiction is explained by the highly asymmetric distribution of nitrogen (skewness coefficient = 0.56), where indicator kriging excels in predicting extreme values and probabilities of exceeding critical thresholds.

Simple kriging (SK) optimizes the prediction of pH with an RMSE of 0.40 and an RMSS of 1.03, reflecting the quasi-normal distribution of this variable after logarithmic transformation. The performance of SK suggests a relatively stable local mean in space, a fundamental assumption of this method. For organic carbon, ordinary kriging (OK) and simple kriging (SK) show equivalent performances (RMSE = 0.41, RMSS = 0.91), indicating an intermediate spatial structure between the local stationarity of SK and the non-stationarity managed by OK. This equivalence suggests a moderate variability of the local mean of carbon over the study area.

Table 4. Comparison of kriging method performance.

RMSE: Root Mean Square Error, RMSS: Root Mean Square Standardized Error, UK: Universal Kriging, OK: Ordinary Kriging, SK: Simple Kriging, IK: Indicator Kriging, N = Nitrogen, pH = Hydrogen Potential, C = Carbon.

3.2.3. Spatial Patterns of Soil Fertility

The analysis of prediction maps confirms distinct spatial patterns for each fertility parameter (Figure 2). The nitrogen map (Indicator Kriging) shows strong heterogeneity. The majority of the territory is in blue and green, indicating low to very low concentrations (below 0.66). Areas with higher values (yellow and red) are localized and dispersed, suggesting point enrichments. The carbon map (Ordinary Kriging) reveals a more structured distribution than that of nitrogen. Average values (light blue and orange) are predominant, but areas of enrichment (dark blue and red) as well as areas of low content (green) are observed. These enrichment zones seem to correspond to basins or convergence zones, which is consistent with the accumulation of organic matter. The distribution is less “spotty” than that of nitrogen. The pH map (Simple Kriging) shows a well-structured and gradual spatial variability. The transition zones are smooth, ranging from acidic values (red, below 6.57) to neutral to slightly basic values (green and blue, above 6.6). Large acidic zones (red) in the south and more neutral zones in the north and center can be clearly identified.

3.3. Assessment of Soil Carbon Stock

Table 5 presents the spatial distribution of soil organic carbon content classes and the corresponding areas. Based on these classes, the average soil carbon stock was estimated at 1.62 tons of carbon per hectare (t∙C·ha⁻¹).

4. Discussion

The geostatistical approach developed in this study presents several advantages

Figure 2. Spatial variability of soil fertility parameters.

Table 5. Spatial distribution of organic carbon classes.

over classical soil fertility mapping methods. Unlike the study by [1], which required 2264 soil observations, 771 soil profiles, and 2148 analyzed samples to characterize soil fertility in Southern Benin, the geostatistical approach of the present study enabled reliable mapping of spatial variability with only 52 strategically distributed samples in the Sudan-Guinea and Guinea zones. The fitting of the spherical model for nitrogen, pH, and carbon, confirms the observations of [40] on the effectiveness of these theoretical models for geostatistical predictions in soil science. The superiority of indicator kriging for nitrogen (RMSE = 0.47) is explained by the highly asymmetric distribution of this parameter (skewness coefficient = 0.56), characteristic of nutrients subject to leaching and volatilization losses under the tropical conditions of Benin. The results obtained agree with those of [45] and [46], who demonstrated the superiority of indicator kriging and ordinary kriging for the spatial prediction of fertility parameters. This convergence reinforces the robustness of the geostatistical approach for differentiated soil fertility management in West Africa.

The analysis reveals a concerning state of soil fertility in the Sudan-Guinea and Guinea zones of Benin. The average nitrogen contents (0.08% ± 0.03%) are well below the threshold of 0.12% considered acceptable for food crops [25]. This nitrogen deficiency, observed throughout the study area, probably results from the intensification of shifting slash-and-burn cultivation practices that destroy soil organic matter and accelerate nitrogen mineralization. The slightly acidic pH values (6.63 ± 0.22) fall within the optimal range for most tropical crops (6 - 7), but they reveal a trend towards acidification in some areas, particularly in the south of the study region. This progressive acidification can limit the availability of phosphorus and exchangeable bases, thus compromising long-term agricultural productivity. The organic carbon contents (0.79% ± 0.17%) correspond to average values for tropical soils, but they remain insufficient to maintain optimal soil structure and resilience to climatic stresses. The carbon stock estimated at 1.62 t∙C·ha⁻¹ reveals a moderate sequestration potential, lower than the objectives of the “4 per 1000” initiative, which aims to increase global soil organic carbon stocks by 0.4% per year. This international initiative, launched during the 2015 Paris Climate Conference (COP21), promotes the enhancement of soil carbon sequestration as a strategy to improve soil fertility and mitigate climate change by offsetting a significant share of anthropogenic greenhouse gas emissions [47].

The observed spatial variability suggests the need to adopt a differentiated rather than uniform management approach. Areas identified as deficient in nitrogen require priority organic inputs, while sectors with low pH require calcium amendment. The produced mapping constitutes a decision-making tool for optimizing resource allocation and reducing input costs. The revealed spatial patterns indicate a combined influence of topographic, geological, and anthropogenic factors. Carbon accumulation zones generally correspond to lowlands and convergence zones, while deficient sectors are located on slopes subject to erosion and intensive agricultural practices.

Although the spatial interpolation provided useful insights into the spatial variability of soil nitrogen, pH, and carbon stocks, the relatively low sampling density (52 samples over 48,347 km²) may limit the resolution and accuracy of the resulting maps. Consequently, fine-scale spatial patterns or localized variations might not be fully captured, suggesting the need for denser sampling or multi-scale analyses in future studies to improve prediction reliability.

5. Conclusion

This study demonstrates the effectiveness of the geostatistical approach for mapping the spatial variability of soil fertility in the Sudan-Guinea and Guinea zones of Benin. Indicator kriging proved optimal for nitrogen, simple kriging for pH, and ordinary kriging for organic carbon. The results reveal an overall deficient nutritional status with low nitrogen contents, slightly acidic pH, and average organic carbon contents. The carbon stock indicates a moderate sequestration potential requiring an improvement strategy. The observed spatial variability justifies the adoption of differentiated fertility management with targeted interventions according to identified deficits. This approach would optimize the efficiency of amendments and sustainably improve agricultural productivity. The produced maps constitute operational tools for decision-makers and producers in planning soil fertility restoration interventions. Extending this approach would contribute to the development of a national strategy for sustainable soil resource management in Benin.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References