Machine Learning Enhances Soil Aggregate Stability Mapping for Effective Land Management in a Semi-Arid Region

Abstract

Soil aggregate stability (SAS) is needed to evaluate the soil’s resistance to degradation and erosion, especially in semi-arid regions. Traditional laboratory methods for assessing SAS are labor-intensive and costly, limiting timely and cost-effective monitoring. Thus, we developed cost-efficient wall-to-wall spatial prediction maps for two fundamental SAS proxies [mean weight diameter (MWD) and geometric mean diameter (GMD)], across a 5000-hectare area in Southwest Iran. Machine learning algorithms coupled with environmental and soil covariates were used. Our results showed that topographic covariates were the most influential covariates in predicting these SAS proxies. Overall, our SAS maps are valuable tools for sustainable soil and natural resource management, enabling decision-making for addressing potential soil degradation and promoting sustainable land use in semi-arid regions.

1. Introduction

Soil aggregate stability (SAS) is one of the main soil physical properties indicating the level of soil resistance to water and wind erosion [1,2,3]. It is also a proxy of soil structure, demonstrating the soil’s ability to control water infiltration and nutrient movement. SAS significantly impacts soil aeration, microbial activity, and tillage performance by improving air circulation, supporting microbial life, and reducing compaction for easier tillage [4].

The spatial distribution of SAS depends on many factors including particle size distribution, plant characteristics (e.g., plant type, growth stages, plant root system, and residues), and land management (e.g., tillage operations), just to name a few [4]. Since SAS is not easily quantified as a single value, certain indices have been proposed as proxies to represent this property: the mean weight diameter (MWD) and the geometric mean diameter (GMD) of soil aggregates [5,6]. These proxies provide insights into soil structure and resilience to erosion [7]. MWD measures the average size of soil aggregates, weighted by their mass. It is often used to evaluate the soil’s physical resistance to external forces such as water erosion and compaction. MWD values greater than 2 indicate more stable and robust soil aggregates [8], contributing to soil health and resistance to degradation [7]. GMD reflects the distribution of aggregate sizes, providing a finer assessment of soil structure by emphasizing the arrangement and connectivity of soil particles [9]. This measure is useful when evaluating soil porosity, aeration, and water infiltration capacity, which are crucial for plant growth and microbial activity [10]. Studies have shown that both MWD and GMD serve as essential proxies for predicting soil health, erosion susceptibility, and nutrient cycling efficiency, thus highlighting their relevance in soil conservation practices and land management efforts [8]. Specific to semi-arid regions, other proxies have been adopted as well [11,12,13,14,15,16]. For example, in some Australian regions with low organic matter, SAS is largely influenced by clay content. Here, wetting-and-drying cycles, especially with the addition of organic amendments like straw, can significantly enhance SAS [17]. Additionally, biological soil crusts (biocrusts)—composed of cyanobacteria, algae, fungi, and other microorganisms—play a crucial role in stabilizing soil aggregates, particularly in arid climates where they act as a biological agent to enhance SAS [18]. However, for humid conditions, their stabilizing effect may be overshadowed by the presence of vascular plants and higher organic matter content [18].

Given the importance of SAS in facilitating land management efforts, it is important to quantify and map its distribution over vast landscapes, and digital soil mapping (DSM) has made this possible and more rapid than solely depending on conventional approaches [19]. DSM has previously been applied to understand the relationship between SAS indices and environmental covariates as well as present their spatial and temporal variations at the landscape scale [20,21]. In DSM, freely accessible environmental covariates such as those from remote sensing data (RS) have frequently been used for such modeling tasks. DSM involves the use of different machine learning algorithms (MLAs) to extrapolate or quantify over space different soil properties using observational data [22,23,24]. However, one reliable and widely applied MLA is the Random Forest (RF) algorithm. RF has demonstrated superior performance in soil property prediction and mapping over the years as it can handle complex relationships between input and target variables, making it suitable for capturing nonlinear patterns in soil datasets [25]. RF is an ensemble approach that combines multiple decision trees eventually resulting in improved accuracy and robustness of model predictions [16,26,27,28,29]. The application of DSM for mapping soils is popular, with varying successes and failures. Despite this, there is still not sufficient information on SAS modeling and mapping for soil management purposes in dryland regions [24,25].

To fill this fundamental knowledge gap, we aimed to estimate two critical indices of SAS (MWD and GMD) based on two MLAs, RF and k-nearest neighbor (k-NN), coupled with selected environmental covariates, in the Lapuee Plains of Southwestern Iran. Specifically, our study objectives were to (1) develop spatial prediction maps for SAS properties and their associated uncertainties; (2) evaluate the performance of different MLAs in predicting SAS indices; and (3) highlight the benefits of SAS maps for operational purposes.

2. Materials and Methods

Our study follows a sequential workflow described in (Figure 1). The main steps were to (i) compile a database including soil observations of SAS indices for the top soil layer (0 to 30 cm) coupled with the environmental covariate, (ii) select the effective environmental covariates through the recursive feature elimination (RFE) method, (iii) model SAS indices using RF and k-NN MLAs based on two scenarios S1 [all covariates (RS indices, topographic attributes, and soil covariates)] and S2 (considering only topographic and RS covariates), (iv) evaluate the performance of each MLA based on conventional model validation criteria, and (v) map both the SAS indices of interest, MWD and GMD, alongside their associated uncertainties. Our reason for adopting a scenario-based approach was to examine how the same model would perform given a varying number of covariates. Particularly, we would be able to establish whether fewer variables, relative to more, would yield a parsimonious outcome or not.

2.1. Study Area Description

Our study area is located between 29°78′ to 29°88′E longitude and 52°68′ to 52°71′N latitude in the Lapuee Plains of Fars province, Southwestern Iran, and covers 5000 ha (Figure 2). The mean annual temperature and precipitation in this area are 17 °C and 446 mm, respectively, a typical BSh (Semi-Arid Climate) class based on the köppen classification system [30]. The climatic data from the Zarghan Meteorological Station was provided by the Iranian Meteorological Organization covering the period from 2001 to 2021 (Zarghan Meteorological Station, 2001–2021). The study area consists mainly of two soil orders, Entisols and Inceptisols. Soils here have relatively low organic matter ranging between less than 1 and 3% [31]. Also, these soils tend to have low nutrient levels limiting vegetation cover and sustenance. Nonetheless, for areas where farming (cultivation) takes place, farmers integrate organic manure to improve soil aggregation and nutrient levels as well as organic matter input. Soils here also have relatively low moisture due to the climatic conditions (xeric and thermic) [32]. The region has three main landscape types ranging from South to North: Plain, Piedmont, and Hill lands. The Lapuee Plain is largely flat, with elevation ranging between 1591 and 1900 above sea level (Figure 2), and predominantly features flat terrain. Over 85% of the area consists of Plain landscapes, Less than 10% of the area is Piedmont landscape also with Hill land areas covering less than 5% of the total region. (Figure 2).

2.2. Field Study and Laboratory Analysis

The study area comprises primary land-use types: agricultural irrigated, dry farming, and pasture lands. A total of 129 georeferenced points were sampled. At each point, surface soils were collected (0–30 cm depth) and transferred to the laboratory where they were immediately air-dried following routine laboratory procedures.

Since our goal was to estimate SAS, we first used a standard sieving method to separate stable aggregates. For this purpose, the soil samples were each passed through a 6 mm sieve to remove debris. After that, 100 g of each soil sample was placed on top of a set of sieves. The sieves were in the ranges >2, 2, 1, 0.5, 0.25, and <0.053 mm open diameters. The set of sieves was then placed in a water-filled container of a wet sieving apparatus, and vertically shaken for 5 min. After shaking, each fraction of aggregate size was washed with water, oven-dried at 65 °C, and weighed [7,8]. Having obtained the weights of different soil aggregate fractions, we then estimated GMD and MWD. These were calculated as follows:

$$\mathrm{M}\mathrm{W}\mathrm{D}={\sum }_{i=1}^n{w}_i{\overline{d}}_i$$

(1)

$$\mathrm{G}\mathrm{M}\mathrm{D}=\mathrm{e}\mathrm{x}\mathrm{p}\left[\sum _{i=1}^n{w}_i\mathrm{Log}{\overline{d}}_i\right]$$

(2)

${\overline{d}}_i$ is the arithmetic mean of the aggregate size on the ith sieve, w_i is the proportion of stable aggregates on the ith sieve, and n is the number of aggregate size groups.

For better interpretation, especially regarding the MWD, values were classified based on the criteria in [8] (

Table 1. The classification of the mean weight diameter (MWD) of soil aggregates according to their stability and crustability classifications [8].

Class	MWD (mm)	Stability	Crustability
1	<0.4	Very unstable	Systematic crust formation
2	0.4–0.8	Unstable	Crusting frequent
3	0.8–1.3	Medium	Crusting moderate
4	1.3–2.0	Stable	Crusting rare
5	>2.0	Very stable	No crusting

).

Along with the measurements of GMD and MWD, we also assessed other physical properties of the soil known to influence SAS. These properties included bulk density (BD) and Atterberg limits—specifically, the Plastic Limit (PL), Liquid Limit (LL), Shrinkage Limit (SL), and plasticity index (PI)—(BD) was measured via the cylinder method, and all Atterberg limits measured following international ASTM standards and protocols (D 422 and D 4318). These properties provide valuable information about the behavior of soil under different moisture conditions, which is directly linked to the formation and stability of soil aggregates.

2.3. Environmental Covariates

We selected environmental covariates including topographic attributes, soil data, and RS data. As mentioned earlier, soil data were BD and Atterberg limit values. RS data were based on Landsat 8 image satellite (OLI/TIRS), including individual bands and vegetation indices commonly used for predicting SAS, the Normalized Ratio Vegetation Index (NRVI), and the Greenness (Green). Vegetation indices serve as proxies for vegetation cover, while individual spectral bands offer insights into the physical properties of land surface features (i.e., soil texture, moisture content, and mineral composition), pertinent to soils [26,33]. Topographic covariates like slope gradient, aspect, and microtopography influence SAS, influencing water redistribution, microclimatic conditions, vegetation cover, and exposure to environmental factors [16,34]. Steeper slopes are more prone to erosion, while microtopographic features can create localized variations in soil moisture and organic matter content [6,35].

Topographic covariates were derived from the ASTER Global Digital Elevation Model (DEM) with a spatial resolution of 30 m. From the DEM, 33 topographic covariates were obtained using SAGA GIS software, version 7.4.0 [36]. For the spectral data, 20 covariates were obtained including the ratios of visible, near-infrared, short-wave infrared, and thermal bands (2, 3, 4, 5, 6, 7, 10, and 11). To do so, the Landsat 8 images (path/row: 163/39) with spatial 30 m resolution were used which corresponded to the soil sampling period. Pre-processing of remote sensing data was performed, and radiometric correction was performed in ERDAS IMAGINE (version 14) software. Additionally, the five corresponding soil covariates previously mentioned were measured in the laboratory. Finally, all environmental covariates were resampled to a spatial resolution of 30 m using the bilinear resampling method in R statistical software, version 4.2.1. We should highlight here that for the soil covariates, each was first mapped using ordinary kriging to obtain wall-to-wall spatial maps for each property.

2.4. Environmental Covariate Selection for Modeling Purposes

In selecting covariates, we approached it based on two scenarios. In the first scenario (S1), we modeled SAS indices using all the environmental covariates (RS, topographic, and soil covariates) and then adopted a recursive feature elimination (RFE) method based on the “caret” package in R to select the most highly predictive set of environmental covariates for modeling SAS indices (MWD and GMD). More detailed explanations about RFE can be found in [37]. For the second scenario (S2), only the RFE-selected environmental covariates (RS and topographic covariates) were used. This comparison allowed us to assess the model performance while using either many or fewer variables.

2.5. Modeling and Mapping Aspects

We implemented two MLAs, the k-nearest neighbors (k-NN) and random forest (RF), to predict MWD and GMD indices under the two scenarios. The main reason for using these two algorithms is that the RF has proven effective for mapping several soil properties with success in the past [38,39,40,41]. Meanwhile, k-NN, which is considered a simple algorithm, was only used to compare our results with those of the RF.

The RF algorithm is a well-known tree-based approach. It involves many decision trees which are combined together to generate a prediction [42]. Individual trees are trained using bootstrap samples of the data, with each set of original data (in the bag) being used for training. Residuals for internal (out-of-bag) cross-validation of this algorithm address non-linear and complex relationships between soil properties and environmental covariates. The RF algorithm is advantageous because it requires few parameters for training while being easy to implement [43]. It also reduces overall error by using subsets of randomly chosen input variables for each tree [44]. Its predictive power increases with the number of independent trees, which decreases collinearity between them [45]. Two optimizing parameters, including the number of trees (ntree) ranging from 100 to 1000 based on an interval of 50 trees [46] and the number of mtrys, were applied to minimize error (RMSE) using the “randomForest” package in R software, version 4.2.1.

The k-NN algorithm is a simple and easy-to-implement method used for both classification and regression tasks [47]. It determines the closest probability of occurrence and prediction values based on the optimal selection of nearest neighbors, aiming to minimize model errors [48]. The k-NN was tuned using the “caret” package in R software, with the optimal number of ‘k’ parameters selected to calibrate the model. This approach was chosen because of its simplicity and effectiveness. Based on these models, we further generated wall-to-wall maps of MWD and GMD together with their associated uncertainty estimates. For the uncertainty estimates, the bootstrap method using n = 100 bootstraps was applied. Specifically, we generated 100 maps for each of the SAS indices and averaged them.

Having randomly divided the dataset into calibration (80%) and validation (20%) sets, three common model evaluation metrics including the coefficient of determination (R²), root mean square error (RMSE), and normalized root mean squared error (nRMSE) were used to assess model accuracy (equations below).

$${\mathrm{R}}^2=1-{\displaystyle \frac{{\sum _{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{a}}_{\mathrm{i}}-{\mathrm{b}}_{\mathrm{i}}\right)}^2}{{\sum _{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{a}}_{\mathrm{i}}-\overline{\mathrm{b}}\right)}^2}}$$

(3)

$$RMSE=\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\left(a_i-b_i\right)}^2}$$

(4)

$$\mathrm{n}\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}={\displaystyle \frac{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}{\overline{\mathrm{b}}}}$$

(5)

where a_i: observed value, b_i: predicted value, $\overline{\mathrm{b}}$: mean of the observed values, n: total number of observations.

3. Results and Discussion

3.1. Summary Statistics of Soil Variables and SAS Indices

A statistical summary of SAS indices is presented in (

Table 2. Summary statistics of soil aggregate stability (SAS) indices along with the other soil covariates for all datasets (n = 129 points).

Aggregate Stability Indices	Unit	Min	Max	Average	SD	CV (%)
MWD	Mm	0.24	2.25	1.64	1.01	61.5
GMD	Mm	0.009	0.81	0.48	0.15	31.2
BD	g·cm−3	0.36	1.95	1.14	0.46	40.3
SL	%	10.1	45	23.2	7.59	32.6
PL		18.6	49.2	39.4	6.04	15.3
LL		30	73	53.1	9.28	17.4
PI		2.4	31.9	14.2	5.44	38.1

Note: MWD: Mean weight diameter; GMD: Geometric mean diameter; Min: Minimum; Max: Maximum; SD: Standard deviation; CV%: coefficient of variation percentage, SL: Shrinkage Limit, PL: Plastic Limit, LL: Liquid Limit, PI: Plasticity Index.

). The minimum and maximum values of the MWD and GMD vary from 0.24 to 2.25 and 0.009 to 0.81 mm, respectively. The mean values of MWD and GMD were 1.64 and 0.48 mm, respectively. The minimum and maximum values of the BD vary from 0.36 to 1.95 g·cm⁻³ and its mean value was 1.14 g·cm⁻³. According to Hazelton and Murphy [49], based on the mean value of BD in this study, soils belong to a low-class BD rating (1 < BD < 1.3 g·cm⁻³). Meanwhile, the minimum and maximum values of the Shrinkage Limit (SL), Plastic Limit (PL), Liquid Limit (LL), and Plasticity Index (PI) vary from 10.1 to 45.0, 18.6 to 49.2, 30.0 to 73.0, and 2.4 to 31.9%, respectively. Moreover, the coefficient of variation (CV%) values for MWD, BD, and PI showed high variability, while GMD, SL, PL, and LL varied moderately among the classes proposed by [45].

3.2. The Importance of Environmental Covariates and Modeling Performance

Based on the output generated by the RFE method, 15 out of the total 58 environmental covariates were selected as the most important predictors of SAS. Carefully selecting a subset of covariates has been shown to improve model outcomes compared with using all available covariates [50]. We selected several covariates for analysis. The first group consisted of five indices obtained from RS data: The Normalized Ratio Vegetation Index (NRVI), Greenness (Green), Difference Vegetation Index (DVI), Soil Adjusted Vegetation Index (SAVI), and Transformed Vegetation Index (TVI). The second group included five parameters derived from topographic attributes: Channel Network Base Level (CNBL), Watershed Basins (WB), Valley Depth (VD), Diffuse Insolation (Diff), and Standardized Height (SH). Finally, we incorporated five soil covariates: BD, PL, SL, LL, and PI (see Table 3). The relative importance of the selected environmental covariates while predicting GMD and MWD based on the two scenarios is presented in (Figure 3 and Figure 4). Based on the ranking of environmental covariates, topographic attributes showed the greatest importance in predicting SAS indices (Figure 3 and Figure 4), which validated the dominant influence of topography on the spatial variability of SAS indices in the study area. To further elaborate on the influence of topographic factors on soil aggregate stability (SAS), factors like Channel Network Base Level, Watershed Basins, Valley Depth, and Diffuse Insolation play crucial roles by influencing soil moisture, erosion patterns, nutrient cycling, and microclimate conditions. Channel Network structure regulates soil moisture and erosion, essential for SAS [51], while Base Level influences water table depth and aggregation processes [51]. Additionally, Watershed Basins shape water flow and nutrient cycles, affecting soil structure [52], and Valley Depth creates microclimatic conditions, thus altering soil cohesion [53]. Meanwhile, Diffuse Insolation modulates soil temperature and moisture, affecting microbial activity critical for soil aggregation [54]. Overall, among all the soil and environmental covariates, the two topographic covariates, CNBL and Diff, had the highest influence on the modeling process for SAS indices (Figure 3 and Figure 4). The CNBL has the highest influence on both GMD and MWD in the study area (Figure 3 and Figure 4). CNBL, representing topography, plays a significant role in SAS indices, aligning with SOC patterns observed in the area due to soil erosion and deposition [55,56].

Topographic attributes, which describe the land physiography, are widely applied in DSM [57,58,59,60]. We observe that these attributes drive SAS variability in the area, influencing soil moisture regimes, organic matter, microbial diversity and activity, and plant cover, as well as surface features known to influence runoff and erosion. All these directly or indirectly influence the formation and stability of soil structure. The variability and distribution of SAS indices have also been shown to be dependent on the abovementioned topographic covariates [6,61]. Meanwhile, soil erosion results in the deposition of surface soil, changes in micro-topographic characteristics, and the composition of soil properties, all of which affect SAS stability [29,62,63,64]. Moreover, research by Parker et al. [53] highlighted that elevation stands out as the most crucial terrain factor influencing the spatial distribution of vegetation. Vegetation, in turn, acts as a source of SOM, directly contributing to the formation and stabilization of soil aggregates [6,54]. Also, topographic covariates are important factors in controlling the effects of climate, water availability, lithological dissimilarity, and biota in soil formation [65,66]. Topographic and climatic variables control the state of soil water, and the dynamics of plant litter mineralization, erosion, and deposition [46,67]. Several studies have confirmed the importance of topographic covariates such as elevation, slope, and aspect [2,57] in the prediction of soil properties. Millard and Richardson [68] also showed the importance of topographic attributes in predicting and mapping soil properties. Bouslihim et al. [21] specifically stated that topography and geology are the most crucial parameters in the spatial prediction of SAS. The results based on the two scenarios indicate that different environmental covariate combinations affect the prediction of SAS indices. The relative importance analysis (Table 3) showed that the presence of three sets of environmental covariates, shown in S1, increased the RF and k-NN performances. The results validate the significance and benefit of soil covariates in influencing the variability of SAS indices in the study area. Soil covariates were ranked lower in importance in the RF algorithm compared with the k-NN algorithm because RF evaluates variable significance based on their contribution to improving node purity across multiple decision trees, which could dilute the impact of individual variables. In contrast, k-NN directly uses soil covariates to calculate distances between data points, making these covariates more prominent in determining predictions. This difference arises because RF averages the importance across many trees, while k-NN relies on immediate distance metrics, highlighting the soil covariates’ significance more clearly (Figure 3a,c and Figure 4a,c). Therefore, soil covariates showed a high ability to predict the changes in some soil secondary attributes, such as SAS indices, and demonstrated acceptable results in improving modelling performance [24]. Kamamia et al. [6] found that soil covariates have an important role in MWD prediction. Our results, in agreement with Zeraatpisheh et al. [24], show that soil covariates have the most important effect in predicting SAS compared with RS indices. DVI and SH had relatively greater effects when predicting MWD, whereas DVI, SH, and WB covariates had lesser effects when predicting GMD (Figure 3 and Figure 4).

While evaluating the performance of our RF and k-NN algorithms (Table 4, also see Supplementary Figure S1), the results based on the two scenarios (S1 and S2) showed that S1 was the best scenario for predicting both GMD and MWD. We link this improvement in model performances to the integration of soil covariates which are crucial and should be considered in SAS prediction. Altogether, both algorithms succeeded in modelling SAS. Overall, the RF algorithm resulted in higher accuracy for the SAS indices (R² = 0.70 for GMD and R² = 0.72 for MWD) than the k-NN algorithm (R² = 0.62 for GMD and R² = 0.63 for MWD). Kamamia et al. [6] predicted GMD and MWD indices via the RF algorithm with R² values of 0.50 and 0.39, respectively, less than when compared with the current study. The main advantages of RF are its high capability to generate robust estimates, its non-linear capabilities to relate the response and independent variables, and its partial interpretability associated with the ability to derive relative importance estimates of each covariate applied in the model [69]. Applying the RF algorithm, Zeraatpisheh et al. [24] also reported comparable R² values of 0.50 to 0.75 while predicting SAS indices in another semi-arid region of Iran. In contrast, compared with our results, higher R² values (0.8 > R² > 0.65) were also reported by Bouslihim et al. [21] in Morocco while also applying the RF algorithm.

3.3. Spatial Prediction and Uncertainty Maps

Predicted maps of SAS indices (MWD and GMD) based on both RF and k-NN algorithms are presented in (Figure 5).

While adopting the best scenario, here S1, we see that both MLAs had comparable overall spatial patterns of MWD and GMD where high, relative to low, lands show low and high SAS predictions, respectively. Nonetheless, observing the spatial patterns in detail, we see that the k-NN contains several artifacts that are mainly introduced by the covariates used here (see Supplementary Figures S2 and S3). The regions with the highest values of MWD (>2 mm) and GMD (0.70 mm) involve the irrigation of croplands. In contrast, the lowest values of MWD (1.2 mm) and GMD (0.2 mm) were found at the northern and southern boundaries, where pastures with sparse vegetation were significantly impacted by high erosion rates. The central parts of the study region have mostly flat landscapes (lowland) and all the dispositional material (alluvial) is assumed to have eroded and been deposited from the higher area (upland) to settle here. Therefore, the accumulation of the fine soil texture components such as clay could have led to more SOC stock and eventually increased SAS levels in the soils. Some physio-chemical properties such as clay, SOC, and calcium carbonate equivalent (CCE) are binding agents that improve soil structure by forming large and stable aggregates [70]. The study area is mainly on the Piedmont plain and Hill land landscapes. By comparing these two landscapes, we observe that for the highlands, our SAS indicators show that soil stability is lower. However, for lowlands where land-use activities (e.g., crop rotation) take place, we observe higher soil stability. Variability in soil structure and quality is linked to variations in land-use types, which are directly correlated with soil aggregate stability [71]. Furthermore, the type of agricultural land, how much fertilizer is applied, and the type of crops can influence SAS [72], which is also relevant to our study area. Agricultural activities are known to encourage higher rooting density and reduce erosion rates. Under these situations, the level of SOC and the percentage of stable aggregates are expected to be high. Our results follow those of Kamamia et al. [6] and Ostovari et al. [70]. According to the land management in the study area, it seems that increasing SAS in agricultural lands may relate to agricultural practices that involve adding crop residues to the soil [73]. This process promotes the formation of macroaggregates, which then contribute to the development of macroaggregates through the influence of transient organic binding agents such as polysaccharides, roots, and fungal hyphae [74]. Additionally, bioturbation by earthworms and soil micro and macrofauna plays a crucial role in this aggregation process. Earthworms and soil fauna enhance soil structure by physically mixing the soil and organic matter, which accelerates the formation of macroaggregates, and their subsequent consolidation into macroaggregates. Their activity also facilitates the distribution of organic binding agents throughout the soil, further supporting the development and stability of soil aggregates [75]. In sparse vegetation pastures, the sparse vegetation may result in weaker soil structures and low structural stability, as vegetation plays a crucial role in protecting the soil from erosion and maintaining aggregate stability [76].

Additionally, a histogram analysis of MWD and GMD across the entire study area revealed that frequency distribution followed a normal distribution with a mean MWD of 1.64 mm and a mean GMD of 0.48 mm (Figure 6) (equal to the mean MWD value of the entire data set) (Table 2). Le Bissonnais [8] classified the stability of aggregates into the five classes of “very unstable”, “unstable”, “medium”, “stable”, and “very stable”, when the values of MWD are <0.4, 0.40–0.80, 0.4–0.80, 0.8–1.3, 1.3–2, and >2 mm, respectively. Based on this classification, most of the soils in our study area are classified as “stable” and “very stable”, which implies they are not prone to a high risk of soil degradation.

Our uncertainty maps for SAS indices based on the two models are presented in (Figure 7 and Figure 8). The least amount of uncertainty was observed from the central boundary to the northern part of the study region, where MWD and GMD are (0.003 mm _ 0.005 mm) and (0.002 mm), respectively (Figure 7 and Figure 8). Given the low level of uncertainty in this area, the high density of observation points, and the strong correlation between the environmental covariates and surface SAS indices, the RF algorithm was chosen as the superior method for predicting these indices. Baltensweiler et al. [77] also reported the effect of the observed data density on predictive model uncertainty, showing that increased data density generally leads to reduced uncertainty in predictions. However, in some parts of the northern and southern regions, the amount of uncertainty was slightly more, which may be a consequence of sedimentation conditions in this area. Furthermore, the uncertainty map showed that the results of the S2 scenario have slightly higher uncertainty than that of the S1 scenario based on the applied statistical indices (Table 4 and Figure 7 and Figure 8). Since the S2 scenario excluded the soil covariates and only used the other environmental covariates as predictors, this may have resulted in lower prediction accuracy which led to higher associated uncertainty here when compared with the S1 scenario.

The results of the model uncertainty S1 scenario showed that the uncertainty for MWD and GMD ranged from 0.005 to 0.15 mm and 0.002 to 0.09 mm, respectively. Shahini Shamsabadi et al. [78] stated that to make better decisions about the quality of prediction maps, the results of model uncertainty should also be considered alongside the model validation statistics. The quantification of model uncertainty is an essential step in DSM if the spatial prediction maps of soil properties are to be applied to management and decision-making procedures [79]. Overall, the uncertainty analysis showed that the maps of MWD and GMD performed well under both scenarios. Therefore, the current map products can be considered useful for decision-making related to agricultural management and ecosystem services monitoring, a significant benefit of DSM already highlighted in several studies [80,81].

Table 3. Environmental predictors were used to predict soil aggregate stability (SAS) indices in the present study (detailed descriptions can be seen in Table S1).

Environmental Variables	Source	Abbreviation	Definition	Scorpan Factors	Reference
Topographic	DEM	CNBL	Channel Network Base Level	Topographic (r)	[54,82,83,84,85]
		WB	Watershed Basins	Topographic (r)
		VD	Valley Depth	Topographic (r)
		Diff	Diffuse Insolation	Topographic (r)
		SH	Standardized Height	Topographic (r)
Soil	LAB	BD	Bulk density	Soil (s)	[86]
		PI	Plasticity Index	Soil (s)
		SL	Shrinkage Limit	Soil (s)
		LL	Liquid Limit	Soil (s)
		PL	Plastic Limit	Soil (s)
RS indices	Landsat 8 images	NRVI	Normalized Ratio Vegetation Index	Organism (o, s)	[51,87,88,89]
		Green	Greenness	Organism (o, s)
		DVI	Difference Vegetation Index	Organism (o)
		SAVI	Soil Adjusted Vegetation Index	Organism (o, s)
		TVI	Transformed Vegetation Index	Organism (o)

Note: r, s, o, and p represent relief, soil, organism, and parent materials in the SCORPAN model, respectively.

Table 3. Environmental predictors were used to predict soil aggregate stability (SAS) indices in the present study (detailed descriptions can be seen in Table S1).

Environmental Variables	Source	Abbreviation	Definition	Scorpan Factors	Reference
Topographic	DEM	CNBL	Channel Network Base Level	Topographic (r)	[54,82,83,84,85]
		WB	Watershed Basins	Topographic (r)
		VD	Valley Depth	Topographic (r)
		Diff	Diffuse Insolation	Topographic (r)
		SH	Standardized Height	Topographic (r)
Soil	LAB	BD	Bulk density	Soil (s)	[86]
		PI	Plasticity Index	Soil (s)
		SL	Shrinkage Limit	Soil (s)
		LL	Liquid Limit	Soil (s)
		PL	Plastic Limit	Soil (s)
RS indices	Landsat 8 images	NRVI	Normalized Ratio Vegetation Index	Organism (o, s)	[51,87,88,89]
		Green	Greenness	Organism (o, s)
		DVI	Difference Vegetation Index	Organism (o)
		SAVI	Soil Adjusted Vegetation Index	Organism (o, s)
		TVI	Transformed Vegetation Index	Organism (o)

Note: r, s, o, and p represent relief, soil, organism, and parent materials in the SCORPAN model, respectively.

Table 4. Statistical criteria of the model performance to estimate aggregate stability indices using machine learning (ML) algorithm for two scenarios.

		MWD (mm)			GMD (mm)
Scenarios (S)	ML Algorithm	R2	RMSE	NRMSE	R2	RMSE	NRMSE
S1 (RS + Topographic and Soil covariates)	RF	0.72	0.16	0.09	0.7	0.09	0.18
	k-NN	0.63	0.19	0.11	0.62	0.07	0.14
S2 (RS + Topographic covariates)	RF	0.69	0.13	0.07	0.69	0.12	0.25
	k-NN	0.58	0.22	0.13	0.59	0.1	0.20

Note: RF: Random Forest; k-NN: k-nearest neighbor; MWD: mean weight diameter; GMD: geometric mean diameter; RS: remote sensing; R²: coefficient of determination; RMSE: root mean square error; NRMSE: normalized root means square error.

Table 4. Statistical criteria of the model performance to estimate aggregate stability indices using machine learning (ML) algorithm for two scenarios.

		MWD (mm)			GMD (mm)
Scenarios (S)	ML Algorithm	R2	RMSE	NRMSE	R2	RMSE	NRMSE
S1 (RS + Topographic and Soil covariates)	RF	0.72	0.16	0.09	0.7	0.09	0.18
	k-NN	0.63	0.19	0.11	0.62	0.07	0.14
S2 (RS + Topographic covariates)	RF	0.69	0.13	0.07	0.69	0.12	0.25
	k-NN	0.58	0.22	0.13	0.59	0.1	0.20

Note: RF: Random Forest; k-NN: k-nearest neighbor; MWD: mean weight diameter; GMD: geometric mean diameter; RS: remote sensing; R²: coefficient of determination; RMSE: root mean square error; NRMSE: normalized root means square error.

3.4. Strengths and Limitations of the Research

Innovative use of machine learning: This study employs advanced machine learning algorithms (RF and k-NN) to create spatial prediction maps for MWD and GMD. This approach enhances the efficiency and accuracy of SAS assessments compared with traditional, labor-intensive methods. Practical application for sustainable management: The resulting digital SAS maps serve as valuable tools for sustainable soil and natural resource management. They facilitate timely interventions to combat soil degradation, particularly in semi-arid regions, thus promoting sustainable land-use practices. Limited generalizability: The findings are specific to the Lapuee Plain and may not be directly applicable to other regions with different soil and environmental conditions. Data collection constraints: The time-consuming and expensive nature of experimental and laboratory methods for measuring SAS indices may have limited the scope and scale of the study, potentially impacting the representativeness of the results. Model transferability: While the ML algorithms used in the study were promising, their performance in predicting SAS indices in other regions with distinct soil-forming factors remains uncertain and requires further exploration. Incomplete consideration of factors: This study focused on topographic covariates, remote sensing indices, and soil covariates, potentially overlooking other influential factors that could affect soil aggregate stability, such as land-use practices or climate variability.

4. Conclusions

In conclusion, this study examined the topsoil variations of soil aggregate stability (SAS) indices in the Lapuee plains of Shiraz, Iran. Due to the time-consuming and expensive nature of experimental and laboratory methods for measuring SAS indices (MWD and GMD), this study aimed to predict these over the given area. The results showed that scenario S1, which incorporated topographic attributes (CNBL, WB, VD, Diff, and SH), remote sensing (RS) indices (NRVI, Green, DVI, SAVI, and TVI), and soil covariates (BD, PL, SL, LL, and PI), resulted in better accuracy (Table 4). We see that some of the topographic attributes emerged as the most influential factors in explaining the spatial variation of SAS indices, followed by soil covariates. Our results confirmed that the topographic covariates, by controlling the state of soil water, the dynamics of mineralization of plant litter, erosion, and deposition, can be used to predict soil properties even in low-relief areas. The performance of the two MLAs was evaluated based on their prediction accuracy and associated uncertainty. Both algorithms demonstrated acceptable performance, but the RF algorithm slightly outperformed the k-NN algorithm, showing higher accuracy and lower uncertainty as indicated by the model validation statistics used.

The high-resolution maps of SAS indices generated in the Lapuee plain in Southwestern Iran are important for understanding the soil’s ability to resist erosion, maintain fertility, and support agricultural productivity. By assessing SAS, researchers, and land managers can identify areas at risk of degradation, especially in the northern and southern parts of the study area, implement targeted conservation practices, and make appropriate and scientific-based managerial decisions to sustain soil health and productivity. This knowledge is crucial for developing effective soil management strategies and most importantly for preserving soils.