The effect of covariates on Soil Organic Matter and pH variability: a digital soil mapping approach using random forest model

ABSTRACT

This research focuses on understanding the spatial variation of Soil Organic Matter (SOM) and pH levels in the North of Morocco. The study employs a comprehensive approach to enhance predictive modelling, incorporating the Boruta algorithm for effective environmental covariates selection and optimizing model parameters through hyperparameter optimization. Utilizing a Random Forest (RF) model with remote sensing indices and topographic features, the research predicts SOM and pH to identify key contributors to their spatial variability. SOM prediction saw significant success, with a notable correlation to remote sensing indices such as the RVI, NDVI, and TNDVI. These indices, indicative of vegetation health and productivity, emerged as primary influencers of SOM. In comparison, the influence of topographic features like elevation, slope, and aspect was found to be less significant. Conversely, predicting pH was challenging due to the minimal spatial variability within the dataset. Addressing this limitation could involve dataset expansion or alternative models for low-correlated data handling. Despite the RF model’s limited efficacy in pH prediction, an observable correlation between SOM and pH was identified, consistent with prior research. Areas with higher SOM exhibited lower pH values, indicating relative soil acidification from organic matter decomposition. The study’s RF model demonstrated potential in SOM prediction using remote sensing indices, but enhancing pH prediction is essential. Future research may explore dataset expansion, diverse sampling, or testing alternative predictive models for better performance with low-correlated datasets. The study offers valuable insights for advanced predictive model development and enriches understanding of soil management practices.

KEYWORDS:

• Digital soil mapping
• feature selection
• soil organic matter
• hyperparameter optimization
• machine learning

1. Introduction

The importance of soil organic matter (SOM) content and pH as critical indicators of soil quality and health has been extensively acknowledged by various authors (Aboutayeb, El Yousfi, and El Gharras 2020; He et al. 2021; Wulanningtyas et al. 2021). These properties significantly influence mineral availability, soil fertility, and crop yield (Laurent et al. 2020). The increase in SOM is linked to carbon sequestration, leading to enhanced physical, chemical, and biological properties, and improved nutrient bioavailability (Ashworth et al. 2017; Falahatkar et al. 2014; Garosi et al. 2022; Karchegani et al. 2012; Moussadek et al. 2014; Mrabet et al. 2012).

Soil pH and SOM content are crucial indicators of soil quality, with low levels of pH and SOM rendering the soil susceptible to degradation and various challenges like erosion, compaction, pollution, salinization, loss of biodiversity, and essential nutrients (Bünemann et al. 2018; Devkota et al. 2022; Doran et al. 2018; Zornoza et al. 2015). To safeguard soil quality, extensive soil measurements, laboratory analysis (Stone et al. 2016), and digital soil mapping techniques are utilized to estimate spatial variability in these properties (Asgari et al. 2020; Tajik, Ayoubi, and Zeraatpisheh 2020; Zeraatpisheh et al. 2017, 2019).

The crucial role of soil pH and SOM in supporting biological, physical, and chemical activities within the soil necessitates comprehensive monitoring on a large scale. However, spatial variations in these properties across diverse soil areas pose a challenge, as highlighted by various studies (Robertson, Crum, and Ellis 1993; Yu et al. 2017). Addressing these variations intensifies the monitoring process, requiring substantial human, organizational, and financial resources for accurate representation of real-world conditions, as emphasized by Bouslihim et al. (2021b) and Gholizadeh et al. (2015).

Efficiently understanding and preventing soil degradation necessitates a thorough study and mapping of these properties. Digital soil mapping (DSM) emerges as a pivotal approach in this regard, enabling precise predictions of soil SOM and pH beyond the capabilities of conventional soil maps. DSM not only offers a cost-effective solution but also provides insights into map accuracy, empowering informed decision-making (Robertson, Crum, and Ellis 1993; Yu et al. 2017).

Traditional techniques like kriging and inverse distance weighting (IDW) have been widely employed for developing soil parameter maps at varying scales (Bouasria et al. 2021; Tibhirine et al. 2023; Zhao et al. 2019). It’s essential to clarify that while techniques like ordinary kriging rely on spatial relationships among known points, advanced geostatistical approaches like regression kriging incorporate multiple environmental covariates impacting soil formation and distribution (L. Chen et al. 2019).

Machine learning (ML) applications in DSM have significantly expanded opportunities (Bouasria et al. 2023; S. Chen et al. 2022). Extensive research underscores its growing importance at continental and national levels, supported by the availability of free and regular datasets, including satellite images at various resolutions, climate data, and topography (Bouasria et al. 2022; Bouslihim et al. 2021; Hengl et al. 2017, 2021; John, Bouslihim, Bouasria, et al. 2022; John, Bouslihim, Ofem, et al. 2022; Shi et al. 2018).

Environmental covariates form fundamental inputs for DSM based on the soil – environment relationship. These include derivatives of digital elevation models (DEM) and remotely sensed spectral data (RS). DEM-derived data represent relief or topographical units, while RS imagery offers valuable information about soil properties through empirical methods (Forkuor et al. 2017). Assessing the influence of environmental covariates on SOM and pH is pivotal for predicting soil properties and understanding their relationships with environmental factors (Lamichhane, Kumar, and Wilson 2019; Samuel-Rosa et al. 2015). Machine learning presents an effective approach for mapping soil parameters based on these covariates.

This study distinguishes itself from previous research by employing a specific set of environmental covariates and remote sensing imagery, enhancing the interpretability of the modelling approach by focusing on relevant predictors. The goal is to provide a framework for leveraging machine learning techniques and high-resolution remote sensing data in predicting and mapping soil properties, thereby improving soil management and conservation practices. The primary objectives include utilizing spectral indices derived from multispectral satellite imagery and terrain attributes from high-resolution digital terrain models to model the spatial variability of SOM and pH using automated hyperparameter tuning of the random forest. Additionally, the study aims to analyse resulting maps to infer soil quality, understand the relationships between environmental covariates and soil properties, and assess their potential impact on land management practices in the Province of Taounate, North Morocco.

The research ultimate aim is to enhance soil management and conservation practices through the effective utilization of machine learning techniques and high-resolution remote sensing data for soil property prediction and mapping.

2. Methods and materials

2.1. Study area

The study area is located in the North of Morocco (Province of Taounate), in the hills north of the Atlas Mountains, 20 km south of Taounate, and 15 km north of Fez; it covers an area of 805 km² between the following coordinates 34.134°/34.392°N and 4.547°/5.103°S (Figure 1). The climate is the Mediterranean. Short periods of low temperatures can occur during the winter, and sometimes it can even snow. In summer, heat waves from the desert can occur. According to the Köppen-Geiger classification, the study area is in the CSA class (Hot-summer Mediterranean climate) with a mean annual temperature of 17.8°C and mean precipitation of 549 mm. The altitudes vary from 135 and 833 m. The area is characterized by a dominance of Forest soils and a combination of Mountain soil & Red-Brown soil. Geologically, the study area is characterized by predominantly sedimentary facies (conglomerate, marl, sandstone, clays, limestones, and dolomites) of Triassic, manifested by alternating marls and limestones and dolomites, which outcrop in the coverage of middle quaternary age.

Figure 1. The geographical location of the study area and soil samples.

2.2. Soil data and descriptive statistics

Grid sampling of 1 km was applied (Figure 1). At each point, 1 kg of disturbed soil from the uppermost 30 cm was collected. In total, 191 samples were collected over two months (November and December 2013). The Walkley-Black method (FAO 2019) was used for organic carbon measurement, and the conversion factor of 1.724 was applied to estimate the percentage of organic matter. Soil pH was determined by measuring the hydrogen ion activity in an aqueous solution with a pH metre. In this study, various general statistical parameters were calculated, including minimum, maximum, mean, median, 1st quartile, and 3rd quartile. Additionally, the Sapiro-Wilk test was chosen to determine if the data follows a normal distribution. This test is a good choice for small data sets and is often used instead of the Shapiro-Francia test, which can be less reliable with small samples. The Sapiro-Wilk test has several benefits: it is simple to use and requires minimal data, it is generally more reliable than other tests for determining normality, and it is widely used and accepted in the research community and is often available in popular statistical software like R.

2.3. Environmental covariates

Several environmental variables influence the distribution of soil characteristics, such as remote sensing indices, topography, climate, and parent materials. However, this type of data is unavailable in many areas, complicating the application of DSM methods. For this reason, we used just two data types freely available worldwide in the present study. Figure 2 shows all the covariates used in this study. First, a Landsat-8 OLI_TIRS_L2SP image acquired on 3 November 2013 with 0% cloud cover for the area of interest and a spatial resolution of 30 m was used to extract twenty-one spectral indices as follows: Red, Green, Blue, Near-Infrared (NIR), Albedo, Brightness Index (IB), Chlorophyll Vegetation Index (CVI), Crop Moisture Index (CMI), Difference Vegetation Index (DVI), Ferrous Mineral Index (FMI), Green Vegetation Index (GVI), Iron Oxide Index (IOI), Log-Ratio Vegetation Index (LRVI), Normalized Difference Vegetation Index (NDVI), Normalized Difference Water Index (NDWI1 & NDWI2), Ratio Vegetation Index (RVI), Soil Adjusted Vegetation Index (SAVI), Soil Color Index (ICS), Soil Redness Index (SRI), Transformed Normalized Difference Vegetation Index (TNDVI). Also, eight terrain attributes such as elevation, slope, profile curvature (PrC), plan curvature (PlC), general curvature (GCurv), Multi-resolution Valley Bottom Flatness (MrVBF), Multi-resolution Ridge Top Flatness (MrRTF) and Topographic Wetness Index (TWI) were extracted from ALOS-PALSAR with a resolution of 12.5 m. Next, all 21 remote sensing indices were resampled to resolution of the topography attributes. Furthermore, the R/RStudio (R Core Team 2021) were using to raster data manipulation using the following packages (raster, terra, rgdal, sp, sf and PROJ).

Figure 2. Presentation of environmental covariates 21 remotes sensing indices & 8 terrain attributes.

CMI: Crop Moisture Index, CVI: Chlorophyll Vegetation Index, DVI: Difference Vegetation Index, Elevation (metre); FMI: Ferrous Mineral Index, GCuv: General curvature (Degrees), GVI: Green Vegetation Index, IB: Brightness Index, ICS: Soil Color Index, IOI: Iron Oxide Index, SRI: Soil Redness Index, LRVI: Log-Ratio Vegetation Index, MrRTF: Multi-resolution Ridge Top Flatness, MrVBF: Multi-resolution Valley Bottom Flatness, NDVI: Normalized Difference Vegetation Index, NDWI 1 & 2: Normalized Difference Water Index, NIR: Near-Infrared, PlC: plan curvature (Degrees), PrC: profile curvature (Degrees), RVI: Ratio Vegetation Index, SAVI: Soil Adjusted Vegetation Index, Slope (Degrees); TNDVI: Transformed Normalized Difference Vegetation Index and TWI: Topographic Wetness Index.

2.4. Environmental covariates selection

Before developing the RF model, environmental covariates selection process is highly recommended to minimize the number of variables and select the most relevant ones. This technique reduces the information redundancy generated by certain variables, which would introduce unnecessary noise while improving and optimizing the modelling results. In the present study, the Boruta algorithm was used for environmental covariates selection. This algorithm was performed in R statistical software using the Boruta package (Kursa and Rudnicki 2010). Boruta evaluates the significance of environmental covariates in relation to the target variable, helping to filter out irrelevant or noise-inducing environmental covariates. This step ensures that the selected environmental covariates are directly related to the outcome of interest, aligning with the goals of the study. The Boruta package’s fundamental approach relies on the Random Forest technique. This technique is pivotal in assessing the importance of each environmental covariate by comparing it with the importance of environmental covariates generated through random value shuffling. These ‘shadow’ covariates are created by randomly shuffling the values of the original environmental covariates. This process allows the Boruta package to effectively determine which covariates are truly important for the model. If environmental covariate’s importance is similar to or lower than that of the shuffled environmental covariates, it suggests that the environmental covariates may not provide meaningful information for prediction. Conversely, if the environmental covariate’s importance is significantly higher than that of the shuffled environmental covariates, it indicates its relevance in making accurate predictions. By iteratively applying this process and marking environmental covariates as ‘confirmed’, ‘rejected’, or ‘tentative’. Boruta ultimately helps in selecting a subset of environmental covariates that contribute the most to the model’s performance. This can lead to improved model accuracy, reduced overfitting, and a better understanding of which variables are truly influential in the dataset. Rejected and confirmed attributes and sometimes tentative attributes are determined based on three scores (Z-score) minimum, average, and a maximum of a shadow attribute. Attributes above the maximum Z-score are considered confirmed attributes and can be used to develop the RF model. Furthermore, a shadow attribute refers to a synthetic environmental covariate that is created by permuting or shuffling the values of a real environmental covariate within the dataset. These shadow attributes are used for comparison purposes during the environmental covariate’s selection process. In the Boruta algorithm, the ‘TentativeRoughFix’ function is instrumental in determining the importance of environmental covariates and iteratively refining the selection of important environmental covariates based on their significance compared to the shadow attributes. It helps in identifying the subset of environmental covariates that provide meaningful predictive power for the model while excluding irrelevant or redundant environmental covariates.

2.5. Random forest model and environmental covariates importance

Random Forest (RF) is a non-parametric technique that consists of many individual tree models trained from bootstrap data samples (Breiman, 2001). The output of each respective tree is aggregated to make a single prediction. This technique can also rank the predictor variable’s relative importance based on the regression prediction error of out-of-bag (OOB) predictions. RF perturbs each variable and calculates the importance as the resulting change in the OOB error (Pouladi et al. 2019; Shi et al. 2018; Taghizadeh-Mehrjardi, Nabiollahi, and Kerry 2016). Two important parameters in RF development are the number of trees (ntree) and the number of variables available for selection in each split (mtry) (Houborg and McCabe 2018).

The Random Forest (RF) algorithm inherently provides an estimation of environmental covariates importance by assessing the mean decrease in node impurity, quantified through metrics such as Gini impurity or mean squared error (MSE), resulting from environmental covariates-based tree splits within the forest. In essence, environmental covariates that significantly reduce variance are considered more pivotal. The environmental covariates importance mechanism within RF assigns a numerical value to signify the extent of each environmental covariate’s impact on enhancing the model’s predictive accuracy. By leveraging this built-in capability, RF automatically ranks environmental covariates by their importance, offering insights into the key determinants behind modelled outcomes. In accordance with this approach, we analysed the environmental covariates’ importance outcomes generated by our RF model to understand the principal environmental variables shaping the spatial patterns of SOM and pH.

2.6. Hyperparameter optimization, model performances and uncertainty

Hyperparameter optimization involves searching for the best parameter values that govern the behaviour of a machine learning algorithm. In the present study, a combined grid search was applied to test several combinations between mtry 1 to 20 and 1 to 14 for SOM and pH, respectively. For ntree, five values were tested (500, 1000, 1500, 2000), using three repeats of 5-fold cross-validation (CV). The best combination was chosen based on the repeated CV’s root mean square error (RMSE) results. The optimal values of ntree and mtry were applied to develop our SOM and pH prediction models. The performance of both models was assessed using a 5-fold cross-validation approach. This approach divides the initial database (191 samples) into several subsets (5 in this case); each time, four subsets are used to train the model, and the fifth subset is used to validate the model. The mean result of the five tests was used to evaluate the final performance of the model. In this sense, three statistical indices, including coefficient of determination (R²), mean absolute error (MAE), and root mean square error (RMSE), were used to evaluate the accuracy of the RF model for the prediction of SOM and pH, and its equations are as follows:

(1) $R^2=1-\frac{{\sum }_i{\left(O_i-P_i\right)}^2}{{\sum }_i{\left(O_i-\bar{O}\right)}^2}$(1)

(2) $\mathit{RMSE}=\surd \left[\sum \frac{{\left(P_I-O_i\right)}^2}{\mathit{n}}\right]$(2)

(3) $\mathit{MAE}=\frac{1}{n}\ast \sum \left|O_i-P_i\right|$(3)

The prediction interval coverage probability (PICP) 90% is a metric that quantifies the accuracy of the uncertainty assessment in the model’s predictions. This metric is calculated as the percentage of cross-validation observations that are contained within the 90% prediction interval. A PICP of close to 90% indicates that the uncertainty in the model’s predictions has been properly assessed. However, a PICP greater than 90% suggests that the uncertainty was underestimated, while a PICP less than 90% indicates that the uncertainty was overestimated (Poggio et al., 2021; Shrestha and Solomatine, 2006). To visualize the uncertainty as a map, the 90th prediction interval (PI90) is also calculated as the difference between the 95th and 5th quantiles of the predictions (Equation 4). A lower PI90 value indicates a more certain prediction, while a higher PI90 value signifies higher uncertainty in the model’s predictions.

(4) $\mathrm{PI}90=\mathrm{q}0.95-\mathrm{q}0.05$(4)

Moreover, the correlation map was performed between pH and SOM raster grids using the ‘corLocal’ function within the ‘raster’ R package. The local correlation coefficient for the two rasters was computed using a focal neighbourhood window of 5 × 5pixels and the Pearson method.

3. Results

3.1. Descriptive statistics

Table 1 shows the summary statistics for the 191 samples. Table 1 shows the values for SOM range from 0.505% to 5.094%, with a first quartile of 1.787, a median of 2.361%, a mean of 2.649%, a third quartile of 3.631, and a maximum of 5.094%. Similarly, the values for pH range from 6.840 to 7.840, with a first quartile of 7.290, a median of 7.400, a mean of 7.393, a third quartile of 7.490, and a maximum of 7.840. The average SOM value found in the study area is higher than the 1.34% value reported by Bouasria et al. (2021) at Sidi Bennour, located in the Doukkala plain (western Morocco). On the other hand, Bouslihim et al. (2021a) reported a higher mean value of organic matter (3.76%) in the Settat-Ben Ahmed plain (Central Morocco).

Table 1. Summary statistics for the 191 samples.

	Min.	1^st Qu.	Median	Mean	3^rd Qu.	Max.
SOM (%)	0.505	1.787	2.361	2.649	3.631	5.094
pH	6.840	7.290	7.400	7.393	7.490	7.840

In this research, it is possible to use a normality test to check if the data follows a normal distribution and then use the Random Forest for prediction. In this situation, the normality test would be used first to evaluate the data distribution. However, it would not have a direct connection to the use of the Random Forest algorithm. Figure 3 depicts the distribution of SOM and pH data, with the histogram indicating that SOM is not normally distributed. A statistical test will be used to provide a more meaningful analysis of the sample.

Figure 3. Correlation matrix plot with a distribution of SOM (a) and pH data (b).

The Shapiro-Wilk test found that there is not enough evidence to reject the hypothesis that the data are normally distributed, since the p-value of 0.4612 is larger than the alpha level of 0.05. However, the SOM data has a p-value of 2.263 × 10⁻⁶, indicating that it is not normally distributed.

The distribution shows that 25% of pH values are strictly below 7.2. On the other hand, 75% of the organic matter values are higher than 1.9. The pH value ranged from near neutrality to slightly alkaline conditions. This value may be attributed to the predominant soil type and the region’s farming practices. The low level and high variability of SOM are mainly due to land use patterns, such as crops, grassland, and farming practices.

The normality of the values is shown in Figure 3, where pH values show a normal distribution of data. However, the distribution is skewed to the left for the organic matter values. Furthermore, Figure 3 shows the matrix correlation between dependent and independent variables. A perfect correlation was identified between SOM, TNDVI, NDVI and SAVI. Other positive correlations were detected between Red & Albedo, NIR & Blue and IB & Albedo. On the other hand, MRVBF, MRRTF, PIC, PrC, and Slope have a low correlation (or even non-linearity) with almost all the other inputs.

3.2. Environmental covariates selection

Figure 4 shows the results of variable selection by Boruta’s algorithm. The blue box plots represent the Z score of the minimum (shadowMin), average (shadowMean), and maximum (shadowMax) of a shadow attribute. The red and green box plots represent the Z-scores of the rejected and confirmed attributes, respectively. For the SOM (Figure 4a), between 29 variables (21 spectral indices and 8 terrain attributes), 20 variables were considered relevant, and all between them are spectral indices except CMI. Besides, other terrain attributes have been rejected and considered irrelevant. For pH (Figure 4b), in the first step, 11 variables were considered relevant (Slope, MRVBF, MRRTF, elevation, LRVI, SAVI, TNDVI, NDVI, RVI, TWI and DVI), 14 variables were rejected. The remaining 4 variables (IOI, NDWI2, IB and GVI) are designated as tentative. The ‘TentativeRoughFix’ function was used to resolve this issue and the 4 variables were consequently rejected. Compared to the SOM, the terrain attributes (Slope, MRVBF, MRRTF and elevation) were ranked as the most relevant, while the spectral indices (RVI, NDVI and TNDVI) for the SOM were found to be the most relevant.

Figure 4. Result of environmental covariates selection by Boruta’s algorithm for a) SOM and b) pH prediction.

3.3. Environmental covariates importance

Here we showed the importance plots as a form of a ‘post hoc analysis’ as it offers insights into the contribution of each Environmental covariates to the RF performance. This understanding is valuable for interpretation of the machine learning models and presents to stakeholders about the most critical variables in the model for predicting SOM or pH. The RF’s importance measures of the predictor variables revealed that remote sensing indices have the most significant contribution in explaining the spatial distribution of SOM (Figure 5a), more precisely, the RVI, NDVI, and TNDVI, which directly relate to vegetation cover.

Figure 5. Environmental covariates importance for a) SOM prediction and b) pH prediction.

Healthy vegetation generally indicates a higher potential for organic matter input to the soil, thus areas with high vegetation indices, such as RVI, NDVI, and TNDVI, might be expected to have higher SOM content. While the indices LRVI, NDWI1, SAVI, NDWI2, and FMI were of less importance, they still significantly contributed to SOM prediction. NDWI1 and NDWI2, which are related to water content, could indicate waterlogged soils where decomposition rates might be slower, leading to higher organic matter content.

On the other hand, as seen in Figure 5b, the variable importance in the pH prediction was significantly different. Contrary to SOM, most remote sensing indices’ contribution was low. Instead, topographic environmental covariates played a significant role, particularly the slope, followed by MRVBF and MRRTF.

Topographic features represent different aspects of the landscape’s physical form, which can significantly influence soil formation and properties. Steeper slopes might experience more erosion or faster water movement, reducing organic matter accumulation, and thus influencing SOM and pH. MRVBF (Multi-resolution Valley Bottom Flatness) and MRRTF (Multi-resolution Ridge Top Flatness) can influence water flow and accumulation patterns, further impacting SOM and pH distribution.

These conceptual relationships between variables highlight the importance of comprehensive environmental data in predicting soil properties and illustrate the value of machine learning in handling complex, non-linear relationships between variables for effective digital soil mapping.

3.4. Model performances and uncertainty

Figure 6 presents the tuning hyperparameter results for the RF model regarding the prediction of SOM and pH. From the same Figure, the best combination of ntree and mtry is provided, leading to the best performance of the model based on the RMSE results. For SOM, the results show that ntree = 500 and mtry = 3 yield the lowest RMSE value (0.366). For pH, the lowest RMSE value (0.15) was obtained by a value of ntree = 2000 and mtry = 2.

Figure 6. RF model tuning for a) SOM prediction and b) pH prediction.

Generally, the default value of mtry for the RF regression model is p/3 (p represents the number of input variables) (Genuer and Poggi 2020). Based on this rule, the mtry for the present case must be 20/3 = 7 for SOM and 11/3 = 4 for pH. The optimum values of mtry obtained by the automated hyperparameter tuning applied in the present study are different from the default value of the RF model for both cases (3 for SOM and 2 for pH). Furthermore, a set of implications can be derived from Figure 6. First, mtry equals 1 always gives the highest value of RMSE for SOM and pH. Also, increasing the number of variables (p) contributes directly to a decrease in RMSE until the optimum value is reached for either SOM or pH, and exceeding these values leads to a significant increase in the error. Generally, the default value proposed by the RF model (mtry equals 7 and 4 for SOM and pH, respectively) gives higher RMSE values than those obtained by the automated hyperparameter tuning (mtry equals 3 and 2 for SOM and pH, respectively). Also, varying the ntree values does not change the RMSE results significantly, and the results for each p remain close.

The optimal mtry and ntree for SOM and pH were applied with a five-fold cross-validation to develop the two models. Table 2 and scatter plots (Figure 7) represent the performance of the RF model for predicting SOM and pH. The performance results show that the RF performed well in predicting SOM (R² = 0.90, RMSE = 0.37%, MAE = 0.28%). The prediction output obtained here was higher than the reports of Žížala et al. (2022), Wiesmeier et al. (2011), and Zhang et al. (2021), respectively, for R². On the other hand, the automated hyperparameter in the RF model performed poorly in pH prediction, producing (R² = 0.26, RMSE = 0.15, MAE = 0.11).

Figure 7. Scatter plots for observed vs. RF-model predicted values for a) SOM and b) pH.

Table 2. RF model performances (5-fold cross-validation) for SOM and pH prediction.

Soil properties	RMSE	MAE	R^2
SOM (%)	0.37	0.28	0.90
pH	0.15	0.11	0.26

In general, automated hyperparameter tuning can help improve the performance of RF models by identifying the optimal combination of hyperparameters that result in the best model performance, as assessed by evaluation metrics such as R², RMSE, and MAE. However, it is essential to note that the performance of RF models can also be influenced by additional factors, including the complexity of the underlying relationships between predictor and response variables, among others. Given this, removing redundant variables before the model is calibrated through automated hyperparameter tuning may lead to a higher performance level for the RF model.

The 90th prediction interval (PI90) map for SOM (Figure 8a) confirms the RF model’s prediction of SOM uncertainty. The map shows a small prediction uncertainty across the study area, with only a few pixels with high uncertainty. This suggests that the model’s predictions are accurate and consistent, and that the model has captured the underlying patterns in the data.

Figure 8. 90th prediction interval map for a) SOM and b) pH.

Conversely, the PI90 map of pH (Figure 8b) is noisier, although a substantial portion of pixels still demonstrates low prediction uncertainty. There are a few areas with high uncertainty, especially in the area drained by the river. This suggests that there are several factors contributing to the variability in soil pH in this region, and the model’s inability to capture all requisite information for accurate predictions in these areas highlights the intricate nature of the pH dynamics.

Despite the poor performance of the RF model in predicting pH, the PI90 map shows that the uncertainty across large parts of the study area is small. This means that the model is confident in its predictions, even though they may not be accurate. This phenomenon is unlikely attributed to overfitting, considering the implementation of a robust 5-fold cross-validation approach that ensures comprehensive exposure to the dataset during training. Instead, it is more likely that the model is not capturing the underlying relationship between the input environmental covariates and the output variable (pH) or the insignificance of spatial pH variations.

3.5. Spatial prediction of SOM and pH

Generally, SOM content ranges between 2.2 to 2.6% in Chaouia Region (Aboutayeb, El Yousfi, and El Gharras 2020) and can reach a value of 6.6% in some specific areas (Bouslihim, Rochdi, and Paaza 2021). For the current study (Province of Taounate), soils with 3% SOM content could be explained by the presence of trees, which contribute to providing soil in terms of organic matter. This finding draws our attention to restitute SOM to the soil using organic amendments and stabilized biowaste like composts. This practice improves bio-physicochemical soil properties and thus enhances its fertility (El-Mrini et al. 2021), leading to increased crop yields.

The maps in Figure 9 illustrate the spatial distribution of SOM and pH. These maps were performed using the best model for each parameter and all input variables in a raster format with a resolution of 12.5 m.

Figure 9. Spatial prediction of a) SOM and b) pH.

The distribution of SOM values varies between 0.66 and 4.77% and does not follow a specific pattern; different ranges are mixed in each area (Figure 9a). According to Hounkpatin et al. (2022) classification, these soils range between Class IV characterized by low fertility and severe limitations (SOM ranges between 0.5 to 1.0%). Class I considered highly fertile soils where SOM exceeds 2%.

Nevertheless, we can generally divide the study area into two parts. The eastern part includes most soils with a low percentage of organic matter, and the western part includes most soils with a high percentage. A close analysis of this distribution, especially by linking it to vegetation distribution, shows that most bare land in the eastern region has low percentages of organic matter. As for the lands containing trees (without specifying their type), the percentages of organic matter are high and can reach 4.77%, which is aligned with the results of Tesfaye et al. (2016) and highlights the effect of land use on the distribution of several soil parameters and specifically organic matter. In terms of acidity, soils are neutral to slightly alkaline. The pH values range from 7.11 to 7.69, as shown in Figure 9b. These findings agree with those of Aboutayeb et al. (2020). They mentioned a pH varying from 7.23 to 7.28 for conventional and No-till systems, respectively, in the Chaouia Region in Northwest Morocco.

Visually, the correlation between SOM and pH is clearly observed in the mapping results (Figure 10). The areas with higher SOM content are also characterized by lower pH values, indicating relative soil acidification. This correlation is consistent with the findings of previous studies, which have shown that the decomposition of organic matter produces organic acids that can lower the pH of the soil (Hong, Gan, and Chen 2019). This highlights the importance of considering the interplay between different soil quality parameters when managing and protecting our soils. For example, efforts to increase SOM content through organic matter management practices may also lead to a decrease in soil pH and potential acidification. Therefore, a holistic approach is needed to balance the benefits and potential drawbacks of different management strategies.

Figure 10. Local correlation between OM and pH predictions.

4. Discussion

Soil pH and organic matter are two parameters that strongly influence different soil functions and the availability of nutrients that can be useful for plants (Bai et al. 2018). Alkaline pH could negatively affect the bioavailability of mineral nutrients. For example, it can reduce phosphorus and ferrous absorption, which cause plant deficiencies and compromises crop yields. Furthermore, several studies demonstrated the effect of pH on the availability of nutrients for plants (McCauley, Jones, and Jacobsen 2009). For instance, potassium is highly available under different pH conditions; phosphorus availability is limited in a range of 5.5 and 7.5. Also, other elements, such as calcium and magnesium, are available in a narrow pH range (between 6.5 and 8.5). Thus, it indicates that this parameter is very sensitive and that small changes in pH can have significant effects such as pesticide efficacy, organic matter decomposition, and nutrient availability (Hock 2012; Leifeld, Zimmermann, and Fuhrer 2008; Neina 2019). Overall, and as mentioned earlier, the range of soil acidity in the study area is appropriate and consistent with most of the soil fertility indicators that may be needed to improve soil fertility.

On the other hand, the importance of SOM cannot be ignored; it plays a crucial role in several soil properties either i) biological: it provides a source of energy and a reservoir of nutrients (N, P, S) and contributes to soil/plant system resilience; ii) chemical: it contributes to cation exchange capacity, buffers pH changes, nutrient availability, reduces toxic cation concentrations or iii) physical: it influences water retention properties, structural stability, and plasticity (Murphy 2015; Obalum et al. 2017; Yang, Chen, and Yang 2019). The SOM directly or indirectly affects various soil functions and services, such as food production and climate change mitigation (Bot and Benites 2005; Lal 2004; Minasny et al. 2017).

This study aimed to develop Random Forest (RF) models for predicting Soil Organic Matter (SOM) and pH using remote sensing indices and topographic Environmental covariatess. The analysis of variable importance was conducted to determine the most influential parameters in explaining the spatial variation of SOM and pH. The results showed that remote sensing indices, such as RVI, NDVI, and TNDVI, had the most significant contribution in explaining the spatial distribution of SOM. Our study results are similar to those reported by Zhang et al. (2021), who mapped regional SOM using an RF model with sentinel-2A and MODIS data, contrary to Wiesmeier et al. (2011). They reported land use as the most important predictor. Several environmental factors influence SOM spatial distribution, including climatic conditions, geological units, plant cover data, and land use data (Zhang et al. 2021). It is likely that the remote sensing data, specifically the RVI and NDVI, are the most contributing factors in the RF model’s performance because they are good indicators of vegetation health and productivity. Therefore, vegetation health is related to soil’s amount of organic matter (Van Geel et al. 2019), so it makes sense that these indices would be good predictors of soil organic matter (SOM). The RVI is a ratio between the red and near-infrared bands. It can be used to monitor vegetation health, as healthy vegetation absorbs more near-infrared radiation and reflects more red radiation. NDVI is a measure of the difference between near-infrared and red bands, and it is also a good indicator of vegetation health and productivity. Both indices are sensitive to chlorophyll content, which is a good indicator of the amount of photosynthesis in the plant. A higher chlorophyll content means that the plant is healthy and actively photosynthesizing, which is directly linked to the amount of carbon uptake (Croft et al. 2015), which is an important component of SOM (Bhattacharyya et al. 2022; Van Geel et al. 2019). Topographic Environmental covariatess, such as elevation, slope, and aspect, may also influence SOM, but to a lesser degree than remote sensing data. For example, elevation can affect the amount of water available to plants and the amount of solar radiation the soil receives (Zapata‐Rios et al. 2016), affecting vegetation growth and organic matter accumulation. However, the effect of topographic Environmental covariatess on SOM may be more subtle and may be overshadowed by the more direct relationship between vegetation health and SOM.

Model performance for pH was weak compared to studies from the Czech Republic, New Zealand, USA, Quebec, and India, respectively (Ramcharan et al. 2017; Reddy et al. 2021; Roudier et al. 2020; Sylvain, Anctil, and Thiffault 2021; Žížala et al. 2022). It is possible that the insignificant spatial variability of pH could be a contributing factor to the poor performance of the RF model. The pH values are all relatively similar (a minimum and maximum of 6.84 and 7.84, respectively). There may not be enough variation in the data for the model to learn from and make accurate predictions. The model will not be able to identify patterns and relationships between the input Environmental covariatess and the output variable.

Another thing to consider is that we are working on a small dataset. This problem can be resolved by increasing the size of the dataset or gathering more samples, which may help improve the model’s performance. Also, the application of other models that could be more robust to the low-correlated dataset can be tested. The same problem was reported by Zeraatpisheh et al. (2022). They reported that missing variables’ variability could cause poor prediction performance. Dharumarajanet al. (2017) also showed poor performance in pH prediction with the RF model (R² = 0.30). The pH distribution does not depend on the vegetation distribution as organic matter. However, despite the poor performance, it is more related to the topography, and the effect of slope and other topographic Environmental covariatess on soil acidity distribution appears. The findings in our study corroborate with the report of Zhang et al. (2019) and Seibert et al. (2007) that the spatial variability of soil pH is influenced by annual temperature range and terrain models. Thus, it is possible to observe the accumulation of low-value land in river channels.

However, low SOM soil types (e.g. Luvisols) and pastoral areas are distributed mainly in the east part of the study area. Therefore, the predicted SOM in that region was relatively low. Moreover, cambisols are typically found in low-relief places, where soil migration driven by external factors such as runoff and running water resulted in soil deposition in the area’s eastward direction, resulting in a more significant SOM. As a result, SOM content in the Taounate region tended to be higher in the west and lower in the east direction. Zhang et al. (2021) came up with contrary conclusions for Songnen Plain in China. The low uncertainty suggests that the model is overconfident in its predictions, even though they are not accurate.

5. Conclusion

In conclusion, this study demonstrates the potential of using high-resolution remote sensing data and machine learning models for digital soil mapping. The Random Forest algorithm, when combined with a collection of environmental covariates, was able to accurately predict soil organic matter content with an R² of 0.90 but performed poorly in pH prediction due to insignificant spatial variability in the study area. Also, vegetation indices play a significant role in the SOM prediction, while the topographic Environmental covariatess significantly influence soil acidity (pH). Despite the limitations in pH predictions, this study provides valuable insights into the spatial variability of soil quality parameters. It highlights the need for further research to improve the region’s predictions of pH and other soil quality parameters.

Acknowledgements

The authors would like to express their gratitude to the anonymous reviewers for their constructive feedback.

Disclosure statement

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

Data availability statement

The datasets used for the current study are available from the corresponding author [YB] on reasonable request.

Additional information

Funding

The author(s) reported there is no funding associated with the work featured in this article.