Machine learning-based prediction of soil aggregate stability for sustainable land management in the eastern Himalayan region of India

1 Introduction

Soil aggregate stability (SAS) is a key indicator of soil health, reflecting the ability of soil particles—composed of sand, silt, and clay—to remain bound together and resist disintegration under external forces such as water and wind erosion (1). As an important indicator of soil structure, SAS strongly influences water infiltration, nutrient movement, aeration, and root penetration, thereby regulating several essential soil functions (2, 3). Stable aggregates create interconnected pore spaces that enhance hydraulic conductivity, water storage, gaseous exchange, and microbial activity, ultimately supporting plant growth and ecosystem functioning (4). Since sustainable soil management is fundamental to environmental security and the well-being of present and future generations (5–7), maintaining SAS is crucial for preserving soil productivity, ecosystem resilience, and long-term agricultural sustainability (2, 8–11).

The formation and persistence of stable soil aggregates are largely governed by interactions between clay particles and soil organic matter (12). Stable soil aggregates provide a well-structured pore network that facilitates the storage and movement of water, air and nutrients, making aggregate stability a widely recognized indicator of soil quality (13, 14). Aggregate formation is a complex process driven by interactions among organic and mineral soil constituents, resulting in the hierarchical arrangement of micro-aggregates within macro-aggregates (13, 15). Soil Organic Carbon (SOC) plays a vital role in soil stability by improving soil structure, water infiltration and root development (16–18). SOC is a key regulator of soil physical, chemical and biological functions, and underpins numerous ecosystem services (19, 20). Through physical protection of soil organic matter from oxidation and microbial decomposition, stable aggregates contribute significantly to carbon sequestration and nutrient retention (3, 21–23). They also improve resistance to compaction and reduce soil degradation and erodibility (2, 24). Aggregate breakdown caused by raindrop impact, runoff, wind or mechanical disturbances leads to the detachment and transport of finer particles, which clog soil pores, decrease infiltration and intensify erosion processes (25, 26). Consequently, aggregate stability directly affects soil mechanical strength, microbial activity, root development, seed germination and overall soil fertility (2, 24, 27).

Besides physico-chemical factors, biological agents also contribute to soil aggregation. Biological soil crusts (biocrusts), composed of cyanobacteria, algae, fungi, lichens and other microorganisms, can enhance aggregate stabilization, particularly in arid and semi-arid environments where vegetation cover is sparse (28). However, in humid mountainous zones such as the eastern Himalayan region, their contribution is likely to be less pronounced because dense vegetation, extensive root systems, microbial activity, and higher organic matter inputs constitute the dominant mechanisms governing soil aggregate formation and SAS (28).

The eastern Himalayan soils in Sikkim state of India are highly susceptible to structural instability, weathering and erosion because of their coarse texture, acidic nature, steep terrain and prolonged exposure to intense rainfall events (29). These conditions frequently accelerate nutrient depletion and land degradation, posing serious challenges to sustainable land use in the region. In such a fragile ecosystem, understanding the factors controlling SAS is essential for assessing erosion susceptibility and developing effective soil conservation strategies. Therefore, evaluating SAS under eastern Himalayan conditions can provide important insights into soil degradation processes and ecological sustainability across mountainous landscapes.

Moreover, SAS is commonly evaluated using indices such as water-stable aggregates (WSA), geometric mean diameter (GMD), mean weight diameter (MWD) and water-dispersible clay (WDC) content (30). Among these, the MWD is considered one of the most reliable indicators for assessing SAS (31). Several soil physical properties, including Soil Organic Carbon (SOC), bulk density (BD) and soil texture—particularly the proportions of sand, silt and clay—strongly influence MWD and soil aggregate formation (25, 32). Salt and clay particles function as important binding agents that help maintain soil structure and stability (33). Therefore, understanding the relationships among these properties is essential for evaluating soil health and implementing sustainable land management practices.

Although substantial progress has been made in the direct measurement of SAS, many conventional techniques remain labor-intensive, expensive, technically demanding and difficult to standardize across studies (25). To overcome these limitations, pedotransfer functions (PTFs) have been developed as practical and cost-effective tools for predicting soil properties using easily measurable variables (34). These PTFs utilize routinely available soil survey data, including soil morphology, texture, structure and pH, to estimate complex soil attributes (35). Consequently, PTF-based approaches offer considerable potential for improving the assessment of SAS, evaluating erosion risks and supporting sustainable soil management in ecologically sensitive regions such as the eastern Himalayas.

Over the past few decades, regression methods have been predominantly used to develop PTFs. Among these, Multiple Linear Regression (MLR) has been widely employed as an empirical model. However, MLR assumes a linear relationship between the ‘input’ and ‘target’ variables, limiting its ability to capture non-linear interactions. In recent years, machine learning (ML) techniques have gained prominence for constructing PTFs due to their ability to model complex, non-linear relationships without requiring a pre-defined relationship between ‘predictors’ and ‘target’ variables (36, 37). Techniques such as Artificial Neural Network (ANN), Random Forest (RF), Support Vector Machine (SVM), etc., have been extensively utilized in the past for predicting SAS indicators like the MWD. As ML frameworks become increasingly integrated into scientific investigations, understanding model behavior and the factors influencing predictions has become as important as achieving high predictive performance (38). Recent studies from the Indian Himalayan region have further demonstrated the potential of ML approaches for predicting difficult-to-measure soil hydraulic properties. For example, Singh et al. (39) developed ML models for predicting quasi-saturated hydraulic conductivity across different land-use systems in the Northwest Himalayas and demonstrated that prediction accuracy improved when models were calibrated using region-specific datasets. Their findings highlight that pedotransfer functions are highly dependent on local soil characteristics and environmental conditions, emphasizing the need for developing localized models rather than relying on generalized equations.

Although ML techniques have recently been applied in the Indian Himalayas for predicting soil hydraulic properties and identifying erosion-prone areas (40, 41), no study, to the best of our knowledge, has developed ML-based pedotransfer functions specifically for predicting SAS or MWD in the eastern Himalayan region. Therefore, the present study addresses this important regional knowledge gap. The present study aimed to: (i) characterize the variability and multicollinearity of soil physico-chemical properties, (ii) develop MWD-based SAS prediction models using MLR, ANN, SVM and RF techniques, and (iii) assess model performance through statistical and graphical evaluation. While MLR offers simplicity and less mathematical complexity, the ML techniques are more computationally intensive but excel at capturing non-linear relationships among soil properties. The goal of this research is to identify the most effective approach for modeling MWD, balancing predictability with computational demand. Ultimately, the findings could guide the integration of these techniques as sub-routines in physical or empirical erosion prediction models for sustainable soil management. The developed pedotransfer functions can provide rapid estimates of SAS, which may subsequently serve as input information for regional soil erosion assessment and conservation planning, complementing ML-based erosion hotspot identification approaches recently demonstrated in Himalayan environments (41).

Furthermore, different state and central organizations measure the basic soil properties for different uses. One such activity of the Government of India is the soil health card, which contains vast soil information. A suitable PTF to predict MWD using this soil information would be an added advantage, which can be used for soil erosion estimation in data-scarce regions of eastern and northeastern India, and for further development of a digital soil erosion potential map of the region. The PTFs should be developed for different regions using localized data, so that they can capture the correct inter-relation between the different parameters. As a suitable PTF for the estimation of MWD is not available for the eastern Himalayan region of India, an attempt has been made in this study to develop a localized PTF which would be beneficial for different stakeholders.

Artificial Neural Network (ANN) has been identified as a robust ML technique for predicting various soil hydrological characteristics, as demonstrated in numerous past studies (26, 42). The algorithm excels at learning complex, non-linear relationships between ‘predictor’ and ‘target’ variables, making it highly adaptable for soil property modeling (43). On the other hand, Random Forest (RF) operates as an ensemble of Classification and Regression Tree (CART) algorithms. This method builds multiple regression trees using randomly selected subsets of ‘input’ variables, and the final prediction is obtained by averaging the outputs of all decision trees (44). In contrast, Support Vector Machine (SVM) represents a robust pattern recognition-based approach designed for both classification and regression tasks. SVM minimizes prediction errors while managing data complexity, making it particularly suitable for precise soil property modeling (45). To enhance model accuracy, cross-validation techniques were employed in constructing PTFs. Cross-validation helps to identify the most effective PTFs by comparing performance metrics across models, ensuring optimal predictions (46, 47).

2 Materials and methods

2.1 Study area

The present study was conducted in the Ranikhola sub-watershed, located in the western part of the East District within the eastern Himalayan range of Sikkim state, India (Figure 1). The study area spans 113.33 km² over Gangtok, Pakyong, Rangpo and Namchi blocks (administrative units), lying between 27.224°N–27.330°N latitudes and 88.491°E–88.706°E longitudes, with elevations ranging from 1,200 m to 1,500 m. The region experiences a humid subtropical climate, characterized by an average maximum temperature of 28 °C in summer and a minimum of 4 °C in winter, and receives an annual rainfall of approximately 2,578 mm, placing it in the humid Himalayan zone. Geologically, the area is dominated by gneiss and schist formations, leading to shallow and nutrient-poor brown clay soils. Dominant land use type involves a mix of tropical forests (97.35 km²; 85.9% of total area) and agricultural activities (10.31 km²; 9.1% of total area), both of which significantly impact local hydrological processes. Intense rainfall, coupled with steep slopes, accelerates soil erosion, rock weathering and soil destabilization. According to the National Bureau of Soil Survey and Land Use Planning (NBSS&LUP), water erosion affects 38.24% of Sikkim’s total area (2,71,405 hectares), with degradation levels ranging from ‘very low’ (<5 tonnes/ha/year) to ‘extremely severe’ (>80 tonnes/ha/year) (48).

2.2 Data collection and laboratory analysis

A total of 126 topsoil samples (0–15 cm depth) were randomly collected from different locations across the study area during winter (January–February 2026), primarily covering agricultural and forest land uses. Out of these, 61 samples were collected from agricultural land and the remaining 65 from forest land. Soil samples were collected in polythene bags following standard procedures, air-dried, and carefully hand-sieved through a 4-mm sieve to remove roots and stones. Eleven physico-chemical properties of soil samples, i.e., %Sand, %Silt, %Clay, Organic Carbon content (%OC), Bulk Density (BD) [g/cm³], pH, Calcium (Ca²⁺) [meq/L], Magnesium (Mg²⁺) [meq/L], Carbonate (CO₃²⁻) [meq/L], Bicarbonate (HCO₃⁻) [meq/L], and mean weight diameter (MWD) [mm], were analyzed using standard methods. The hydrometer method was employed to determine particle size distribution, including the proportions of sand, silt and clay (49). Organic Carbon content (%OC) in soil was measured using the wet oxidation method of Walkley and Black (50), while bulk density (BD) was estimated using the core cutter method (49). Additionally, soil pH, Ca²⁺, Mg²⁺, CO₃²⁻ and HCO₃⁻ ions were analyzed following standard protocols provided by the Indian Council of Agricultural Research-Central Soil Salinity Research Institute (ICAR-CSSRI), Karnal, Haryana state, India (51). A flowchart containing the complete step-by-step methodology is presented in Figure 2.

Soil aggregate stability (SAS) was estimated using the wet sieving method as outlined by Kemper and Rosenau (52). Air-dried soil samples were gently hand-crushed and sieved through a 4-mm mesh, and the portion retained on the sieve was used for analysis. From this retained sample, 50 g of soil was placed on a stack of sieves arranged in decreasing opening sizes (2 mm, 1 mm, 0.5 mm, 0.25 mm and 0.053 mm). The sample was evenly distributed on the 2-mm sieve and submerged in water to allow thorough soaking. Subsequently, the sieves were oscillated vertically at a rate of approximately 30 times per minute in a water tub, following the method described by Yoder (53). After sieving, soil aggregates retained on each sieve were oven-dried at 105 °C for 24 hours and weighed. This process was repeated for all sieve sizes, and the mean weight diameter (MWD) was calculated using the following formula (Equation 1):

(1) MWD=∑i=1NMi×Yi

where Y_i = mean diameter of each fraction size on the i^th sieve (mm), M_i = proportion of the total sample weight on the i^th sieve, and N = number of fractions. For a more meaningful evaluation of SAS, the MWD values were interpreted using the classification thresholds proposed by Bissonnais (54), as presented in Table 1.

Table 1. The classification of the mean weight diameter (MWD) of soil aggregates according to their stability and crustability categories (54).

2.3 Pre-processing of datasets for PTF development

Ten soil physico-chemical properties such as %Sand, %Silt, %Clay, %OC, BD, pH, Ca²⁺, Mg²⁺, CO₃²⁻ and HCO₃⁻ were used as ‘input’ (predictor) variables to predict the ‘target’ (predictand) variable, which is mean weighted diameter (MWD) (2, 26). Pedotransfer functions (PTFs) were developed separately for individual land-use classes. Multicollinearity and collinearity among ‘input’ and ‘target’ variables were assessed using IBM SPSS Statistics v22. The Variance Inflation Factor (VIF) analysis was performed to detect multicollinearity, and variables with VIF values exceeding 10 were excluded from modeling (55). The Pearson correlation test identified highly correlated predictor variables, which were prioritized for model development, while less correlated variables were included through trial-and-error. The selected ‘input’ variables were used in the PTF model development and evaluated using various performance metrics (56). In this study, three supervised machine learning (ML) techniques such as Artificial Neural Network (ANN), Support Vector Machine (SVM) and Random Forest (RF), along with empirical Multiple Linear Regression (MLR), were employed. For model development and validation, the dataset was randomly partitioned into an 80:20 ratio, with 80% used for model training and 20% reserved for model testing. The 10-fold cross-validation technique was performed exclusively on the training subset during model development, while the testing subset was used only for final model evaluation. To account for variability associated with random data partitioning, the procedure was repeated 50 times using different random partitions, and the final performance metrics were calculated as the average across the 50 repetitions.

2.4 PTF development using regression and ML models

2.4.1 Multiple linear regression

Regression analysis is a statistical technique that uses the relationship between two or more predictor (independent) variables to predict a predictand (dependent) variable. The general form of the Multiple Linear Regression (MLR) model is as follows (Equation 2):

(2) Yi=β0+(β1·X1i+β2·X2i+……………….+βp·Xpi)+ϵi

where, Y_i is the dependent variable; X_i is the independent variable; β₀ = intercept; β₁, β₂ …… β_p = unknown regression coefficients (slopes); and ϵ_i = random error term with a mean of zero and constant variance (8).

2.4.2 Artificial neural network

Artificial neural networks (ANNs) excel at forecasting and predicting outcomes, especially in scenarios with non-linear relationships and non-stationary time series, where clear patterns are challenging to discern. Functioning similarly to the human brain, they process non-linear information through interconnected nodes called ‘neurons’. Signals are transferred via links with associated weights, and each neuron’s output is determined by applying an activation function to the net input (57). This computational technique learns from experience, similar to the human brain, and leverages acquired knowledge to enhance output accuracy. The output (Y_k) of neuron ‘k’ in the hidden layer is as follows (Equation 3):

(3) Yk=f[ωk+∑i=n(ωk·xi)]

where, ω_k = bias of the output neuron ‘k’, ω_ki = connection weight between hidden neuron ‘i’ and output neuron ‘k’, x_i = output from neuron ‘i’ in the hidden layer, n = number of neurons in the hidden layer, and ‘f’ is the activation function, also known as the transfer function.

An ANN architecture typically comprises of three layers: ‘input’, ‘hidden’ and ‘output’, with information flowing sequentially through these layers to contribute to the final output. In the ‘hidden’ layer, an activation function is employed, which derives the output from a set of inputs. The sigmoid activation function is most commonly used due to its ability to account for non-linearity between the ‘predictor’ and ‘predicted’ variables. In this study, a Feed-Forward Neural Network (FFNN), specifically a Multi-Layer Perceptron (MLP), was employed, operating with a single-direction flow from the ‘input’ to the ‘hidden’ and finally to the ‘output’ layers. The ANN models were developed using R-programming software v4.6.0, with the ‘neuralnet’ package (8). The ANN architecture was determined through a trial-and-error approach, in which different combinations of ‘hidden’ layers and numbers of ‘hidden’ neurons were evaluated, and the architecture providing the best predictive performance was selected as the final configuration. The number of hidden layers and neurons in each layer was specified using the ‘hidden’ argument of the ‘neuralnet’ function. The ‘input’ layer consisted of the selected soil ‘predictor’ variables, while one ‘output’ neuron was used to predict mean weight diameter (MWD). The final ANN architecture was selected based on the model performance criteria used in this study.

2.4.3 Support Vector Machine

The Support Vector Machine (SVM) is a supervised machine learning algorithm designed for both classification and regression analysis. SVM for regression analysis was used in this study to predict the distinct values. This algorithm partitions the dataset into different classes by creating hyperplanes that are not split linearly (58, 59). The core idea behind SVM is to employ a linear model to construct non-linear class borders by translating input vectors into a high-dimensional feature space in a non-linear manner. In the original space, non-linear decision boundaries can be represented by a linear model in the new space (60).

Consider a dataset: D={xi,yi}i=1N(where, x_i = vector of the ‘predictor’ variables, y_i = vector of the ‘predictand’ variable, n = number of data points). The regression function of the SVM is expressed as follows (Equation 4):

(4) Y=f(x)=wi·ϕi(x)+b

where, ϕ_i(x) = feature of the predictors, w_i = weight vector, and b = bias term. The coefficients (w_i and b) are estimated by minimizing the following regularized risk function (Equations 5–8):

(5) R(ϕ)=12|w|2+C·∑i=1N(ξi+ξi*)

subject to:

(6) {yi−f(xi)}≤(ϵ+ξi)

(7) {f(xi)−yi}≤(ϵ+ξi*)

(8) ξi, ξi*≥0 ∀i

The selection of kernel functions and ideal sets of tuning parameters are the main factors that affect the performance of non-linear support vector regression (NLSVR) models. To minimize the empirical error, a positive regularization parameter (C) is utilized, which is similar to the weight parameter in the optimization problem. It is known as the cost/penalty parameter. Furthermore, ‘ϵ’ is the width of the tube and is measured perpendicular to the hyperplane. This tube is known as an ϵ-insensitive tube. The loss function for the selected error tolerance ‘ϵ’ penalizes training errors using the slack variables. ξi and ξi* (61). Both ‘C’ and ‘ϵ’ are user-determined hyperparameters. The distance between the real values and the tube-corresponding values is represented by two positive slack variables, ξi and ξi* . Points below this tube are called. ξi* , and points above this tube are called ξi (62). An ϵ-insensitive tube is the area limited by (y_i ± ϵ). In the SVM model, the Radial Basis Function (Gaussian) kernel was used for model building. The SVM model was developed in R-programming software v4.6.0, using the ‘e1071’ package.

2.4.4 Random Forest

Random Forest (RF) is a supervised machine learning algorithm, which was invented by Breiman (63). The RF is a versatile and intuitive ML algorithm that was designed to reduce the limitations of single regression trees (63). RFs are ensembles of Classification and Regression Trees (CART). The RF models can deal with both continuous and categorical variables (64). Regression trees are constructed by applying a bootstrap to all datasets, with splits at each node from the randomly selected predictor variables of the entire predictor set. This helps reduce the problem of overfitting (65). All trees give different outputs, and the final result is calculated by averaging the outputs of all trees (44). The RF model incorporates the variable importance of all input variables. This helps generate PTFs in this study. A less important variable shows little importance for model prediction; thus, it can be removed from model building.

Two tuning parameters were used in the RF model prediction; these are ‘ntree’, which corresponds to the number of trees, and ‘mtry’, which corresponds to the number of variables randomly sampled at each split. The value of ‘mtry’ will not exceed the number of input variables, and if the number of trees (ntree) is greater, the model is over-fitted (64). The RF model was developed in R-programming software v4.6.0, using the ‘randomForest’ package.

2.5 Performance evaluation of developed models

The performances of the ANN, SVM, RF and MLR models were assessed using a 10-fold cross-validation procedure. It provides accurate and unbiased predictions on small sample datasets (66). Cross-validation provides a modeling framework for partitioning dataset into several calibration and validation sets, ensuring that each sample is assigned to validation at least once. Therefore, a 10-fold cross-validation procedure was used to evaluate the prediction performance of all PTF models. Each cross-validation process was repeated 50 times using different random partitions (80:20) which produced a significant randomness effect on the creation of training and testing datasets. It was carried out to ensure robust statistical estimates and account for variability inherent in random data splitting procedures (67, 68). The performances of PTF-based ML models were assessed using suitable graphical and statistical indicators. The statistical evaluation criteria include: (i) Coefficient of Determination (R²), (ii) Root Mean Square Error (RMSE), (iii) Mean Absolute Error (MAE), (iv) Nash-Sutcliffe Efficiency (NSE), and (v) Akaike Information Criteria (AIC) (Equations 9–14). The AIC metric was used as a complementary model-selection criterion alongside R², RMSE, MAE, and NSE to account for the trade-off between model fit and complexity. Lower AIC values indicate a more favorable balance between goodness-of-fit and model complexity (69). AIC has previously been used as an additional criterion for evaluating ML-based models in soil-related modeling applications (70, 71).

These five statistical indicators (i.e., R², NSE, RMSE, MAE and AIC) were calculated for each of the 50 repetitions, and the results were summarized as “mean ± standard deviation” across the repeated train-test sets. To analyze model performance based on R², model performance categorization criteria given by Li et al. (72) was used, i.e., R² < 0.5 (unsatisfactory prediction capacity), 0.5 < R² < 0.75 (satisfactory prediction capacity), and R² > 0.75 (excellent prediction capacity).

(9) R2=[∑i=1n(Oi−O¯)·(Pi−P¯)∑i=1n(Oi−O¯)2·∑i=1n(Pi−P¯)2]2

(10) RMSE=1n∑i=1n(Pi−Oi)2

(11) MAE=1n∑i=1n|Pi−Oi|

(12) NSE=1−[∑i=1n(Oi−Pi)2∑i=1n(Oi−O¯)2]

(13) AIC=n·ln(SSE)−n·ln(n)+2(K+1)

(14) SSE=∑i=1n(Pi−Oi)2

where, O_i = measured values, P_i = predicted values, O¯= average of the measured values, n = number of data points, and K = number of predictor variables.

3 Results

3.1 General characteristics of soil physico-chemical parameters

Statistical analysis was performed on various soil physico-chemical parameter datasets (%Sand, %Silt, %Clay, %OC, BD, pH, Ca²⁺, Mg²⁺, CO₃²⁻, HCO₃⁻ and MWD) for two major land use classes, i.e., agriculture (n = 61), forest (n = 65), and also for pooled dataset (n = 126). General statistical estimates, including minimum, mean, maximum, range, standard deviation, skewness and kurtosis are summarized in Table 2.

Table 2. Descriptive statistics of all soil physico-chemical parameters.

Across the study area, significant variability in the soil physico-chemical parameters was observed. The mean bulk density (BD) of agricultural soil was 1.25 g/cm³, whereas that of forest soil was 1.17 g/cm³. The mean BD value of soil from the forest land is lower than that of agricultural land, which is attributed to the impact of frequent cultivation on soil structure. Loose, well-aggregated, and porous soils generally have lower BD and higher Organic Carbon content (%OC) (73). The pH of both land uses indicates the acidic nature of the soil in the watershed, which is attributed to high rainfall and ion leaching. The %OC in forest soil is higher (4.23%) than that in agricultural soil (3.37%). Plant litter, including dead leaves, stems, and logs, is a major source of forest organic matter (74). Soil organic matter participates in various processes in the formation and stability of soil aggregates (24). Based on the average mean weight diameter (MWD) values presented in Table 2 and the classification proposed by Bissonnais (54) in Table 1, the agricultural (1.09 mm) and forest (1.11 mm) soils fall within the ‘Medium Stable’ soil aggregate stability (SAS) class. This indicates that, on average, the soils of the study area possess moderate resistance to aggregate breakdown and erosion. The MWD of forest soils was higher than that of agricultural soils, indicating greater aggregate stability and resistance to aggregate breakdown and erosion in forest soils (75).

Although the numerical difference is small, both values fall within the ‘Medium Stable’ SAS category, indicating comparable levels of SAS under the two land-use systems. Welch’s t-test indicated statistically significant differences in pH, (Ca²⁺ + Mg²⁺), %OC and BD between agricultural and forest soils at the 0.1% significance level (Figure 3). In contrast, %Sand, %Silt, %Clay, MWD, CO₃²⁻, and HCO₃⁻ did not exhibit any statistically significant differences between the two land-use types. Nevertheless, the lower MWD observed in agricultural soils is consistent with their comparatively lower soil %OC (3.37%) than forest soils (4.23%). Soil organic matter plays a vital role in the formation and stabilization of soil aggregates by promoting particle binding and enhancing aggregate cohesion. Therefore, the comparatively lower %OC in agricultural soils may have contributed to the observed reduction in aggregate stability, highlighting the importance of maintaining appropriate soil %OC for improving soil structural quality. Furthermore, sand fraction (%Sand) is found more in agricultural soil (61.34%), while clay fraction (%Clay) is found more in forest land (30.11%). Therefore, a high %Sand and low %Clay enhance soil vulnerability to erosion because clay particles hold soil particles together and make them resistant to detachment (76).

The skewness coefficients of agricultural soil parameters [%Sand, %Silt, %OC, BD, (Ca²⁺ + Mg²⁺), CO₃²⁻, and HCO₃⁻] range from –0.47 to 0.44, indicating moderate symmetry. However, pH, %Clay and MWD values deviate from a Normal Distribution. For forest land use, the skewness coefficients of %Silt, %OC, pH, (Ca²⁺ + Mg²⁺), MWD and HCO₃⁻ parameters show moderate symmetry. The mean value was considered as the center of distribution for these parameters (77). Furthermore, the watershed soils were classified into four major USDA (United States Department of Agriculture) textural classes using the USDA triangle plot, namely, (i) sandy clay, (ii) sandy clay loam, (iii) sandy loam, and (iv) loamy sand (Figure 4). The spatial variation maps of all the soil physico-chemical parameters are illustrated in Figures 5a–j.

3.2 Correlation between soil physico-chemical parameters

The effects of different soil physico-chemical properties on MWD are illustrated using correlation matrices for both agricultural and forest land uses, as well as for the pooled dataset. The scatter plot matrix involving agricultural soil properties is shown in Figure 6. The diagonal of the matrix displays density plots for various soil parameters. In the lower section, scatter plots depict relationships between different parameters, whereas the upper section showcases correlation values among the studied soil parameters. Correlation analysis highlights %Clay, %OC and BD as the most relevant soil parameters influencing SAS or MWD, with respective coefficients of 0.221, 0.258 and –0.506, respectively. For forest land use, as depicted in Figure 7, %Sand, %Clay, %OC, BD, pH and (Ca²⁺ + Mg²⁺) emerged as major parameters influencing MWD, exhibiting significant correlations of –0.434, 0.453, 0.451, –0.529, –0.417 and –0.386, respectively. The correlation matrix for the pooled dataset underscores %Sand, %Clay, %OC and BD as major influencers of MWD, with significant correlation values of –0.205, 0.307, 0.302 and –0.482, respectively (Figure 8). Notably, these findings align with previous studies emphasizing the significant impact of %Clay, %OC and BD on soil aggregate stability (SAS) or MWD (15, 32, 78).

The analysis of scatter plots from different land use datasets revealed a consistent pattern in which %Clay, %OC and BD emerged as three major soil parameters affecting SAS or MWD. This observation is in line with the existing literature underscoring the importance of %Clay and %OC as key binding agents, which are influenced by soil pH, polyvalent cations, and their mobility (1, 15, 79, 80). The correlation analysis also revealed a negative correlation between pH and MWD. This indicates that SAS decreases with an increase in soil pH. Past studies have also justified this phenomenon, as at lower pH (<7), microbial activity increases, promoting macro-aggregation and stability. Cations have higher solubility and mobility, facilitating their association with soil %Clay and %OC (81, 82).

3.3 MWD prediction using MLR, ANN, RF and SVM

3.3.1 Multiple linear regression-based PTFs

The Multiple Linear Regression (MLR)-based pedotransfer function (PTF) models were developed through the outcomes of correlation analysis and model performance assessment. Collinearity (between two variables) or multicollinearity (involving more than two variables) can affect model accuracy (83), and errors in the fitting of regression models. Identifying and addressing these issues involves selecting suitable predictors and eliminating those affecting model output (56). Parameters exhibiting collinearity were excluded and examined for Variance Inflation Factor (VIF) analysis. In the present study, the scatter plot matrix for agricultural land use identified multicollinearity between BD, pH and (Ca²⁺ + Mg²⁺), and collinearity among %Sand and %Clay. Consequently, %Sand and (Ca²⁺ + Mg²⁺) were thus excluded from the ‘input’ parameter dataset. Forest soil parameters showed no multicollinearity or collinearity among them. Moreover, the pooled dataset’s scatter plot matrix revealed multicollinearity between %Silt, %OC and CO₃²⁻, and among pH, (Ca²⁺ + Mg²⁺) and HCO₃⁻, leading to the removal of CO₃²⁻ and HCO₃⁻ from ‘input’ parameter dataset. According to the VIF analysis, none of the ‘input’ parameters showed a VIF larger than 10 for agricultural and forest lands, as well as for the pooled dataset (Table 3).

Table 3. Multicollinearity between ‘input’ and ‘target’ variables for different land uses.

The MLR models were fitted to predict MWD across the studied land uses. The PTFs using %Clay, %OC, BD, and pH exhibited better performance for agricultural soil dataset, yielding R² values of 0.43 and 0.14 for training and testing phases, respectively. Likewise, MLR models were developed for forest areas and pooled dataset, based on optimal predictor combinations. The resulting MLR models for agricultural and forest land uses, and pooled dataset are expressed in (Equations 15–17), respectively.

1. Agricultural Land Use:

(15) MWD=1.863−0.0239×(%OC)−2.2945×BD+0.00169×(%Clay)+0.366×pH

ii. Forest Land Use:

(16) MWD=3.2656−0.0625×pH−0.092×(Ca2++Mg2+)+0.051×(%OC)+0.896×BD −0.0137×(%Silt)−0.147×CO32−−0.0163×(%Sand)−0.011×(%Clay)

iii. Pooled Dataset:

(17) MWD=2.109+0.044×pH−0.012×(%OC)−1.086×BD−0.00003×(%Silt)+0.004×(%Clay)

The scatter plots of ‘predicted’ and ‘measured’ MWD values for agricultural and forest land uses, and pooled dataset are depicted in Figures 9a–f. The model statistical performance indicators for the training and testing phases are summarized in Table 4. The MLR model demonstrated unsatisfactory prediction accuracy (R² < 0.50), which is consistent with the criteria proposed by Li et al. (72).

Table 4. MLR model performance indicators for different land uses during ‘training’ and ‘testing’ phases.

The numerical values of model performance statistics are shown as: mean ± standard deviation, estimated from 50 repetitions of train-test split and 10-fold cross-validations. Bold values indicate best PTF model for a particular land-use pattern, identified separately for training and testing phases.

3.3.2 Random Forest-based PTFs

The Random Forest (RF)-based PTF models were developed using a 10-fold cross-validation procedure, and their predictive performance was evaluated using RMSE, MAE, R², NSE and AIC statistical indicators. The optimal models were selected based on combinations of predictor variables that yielded the highest predictive accuracy. For the agricultural soil dataset, the RF model achieved an R² of 0.89 and an RMSE of 0.08 during training, and an R² of 0.66 with an RMSE of 0.11 during testing. For the forest soil dataset, the model attained an R² of 0.87 and an RMSE of 0.05 in the training phase, while the testing phase yielded an R² of 0.69 and an RMSE of 0.08. The pooled (combined) dataset produced an R² of 0.87 and an RMSE of 0.07 during training, and an R² of 0.58 with an RMSE of 0.15 during testing. The performance statistics for all RF models are presented in Table 5, indicating satisfactory predictive capability across different land-use categories (0.50 < R² < 0.75). To further improve model performance, a trial-and-error hyperparameter optimization approach was employed. The optimal parameter combinations were: ntree = 300 and mtry = 5 for agricultural soils, ntree = 60 and mtry = 5 for forest soils, and ntree = 70 and mtry = 5 for the pooled dataset (Table 6). Scatter plots comparing ‘predicted’ and ‘measured’ MWD values for the training and testing datasets are presented in Figures 10a–f.

Table 5. RF model performance indicators for different land uses during ‘training’ and ‘testing’. phases.

The numerical values of model performance statistics are shown as: mean ± standard deviation, estimated from 50 repetitions of train-test split and 10-fold cross-validations. Bold values indicate best PTF model for a particular land-use pattern identified separately for training and testing phases.

Table 6. Calibrated (tuned) hyperparameters of different machine learning models.

The RF model output generated variable importance plots, which were used to identify the most influential predictor variables for developing the PTFs. The variable importance plots for the agricultural, forest, and pooled datasets are presented in Figures 11a–c. These plots rank the predictor variables according to their contribution to model performance, thereby distinguishing the most and least influential soil physico-chemical properties. Based on the variable importance analysis, PTFs incorporating %Silt, %Clay, %OC, BD, and pH consistently exhibited superior predictive performance across the agricultural, forest, and pooled datasets.

3.3.3 Support Vector Machine-based PTFs

The Support Vector Machine (SVM)-based PTF models were developed by optimizing several hyperparameters, including Kernel Function, Cost/Penalty Parameter (C), Epsilon (ϵ) and Gamma (γ) (Table 6). Model performance was evaluated using a 10-fold cross-validation procedure and assessed based on RMSE, MAE, R², NSE, and AIC. Among the evaluated kernel functions [Linear, Sigmoid, Radial Basis Function (RBF), and Polynomial], the RBF kernel produced the highest predictive accuracy, consistent with the findings of Bhattacharya et al. (8). The remaining hyperparameters were optimized through a trial-and-error approach to further improve model performance. For the agricultural soil dataset, the SVM model achieved an R² of 0.71 and an RMSE of 0.13 during training, and an R² of 0.61 with an RMSE of 0.11 during testing. For the forest soil dataset, the corresponding R² and RMSE values were 0.95 and 0.03 for the training phase, and 0.77 and 0.07 for the testing phase. The pooled (combined) dataset yielded an R² of 0.98 and an RMSE of 0.03 during training, and an R² of 0.72 with an RMSE of 0.12 during testing. Here, the gap between the training and testing performance shows a moderate level of overfitting. This may be due to heterogeneous pooled datasets wherein soils come from diverse land-use systems, thus having substantially varying predictor-response relationships. Hyperparameter optimization was performed using a bounded trial-and-error approach combined with 10-fold cross-validation, where the Kernel Function, C, ϵ and γ were systematically evaluated within pre-defined ranges before selecting the optimal hyperparameter combination (Table 6). The complete performance statistics of SVM model are presented in Table 7. The SVM model demonstrated acceptable (satisfactory) predictive capability (0.50 < R² < 0.75) for all datasets and exhibited improved performance for PTFs based on %OC, BD, %Clay, %Silt and pH. Scatter plots comparing ‘predicted’ and ‘measured’ MWD values for the training and testing datasets are presented in Figures 12a–f.

Table 7. SVM model performance indicators for different land uses during ‘training’ and ‘testing’ phases.

The numerical values of model performance statistics are shown as: mean ± standard deviation, estimated from 50 repetitions of train-test split and 10-fold cross-validations. Bold values indicate best PTF model for a particular land-use pattern identified separately for training and testing phases.

3.3.4 Artificial Neural Network-based PTFs

The optimal Artificial Neural Network (ANN) architectures for the agricultural, forest, and pooled datasets are presented in Figures 13a–c. Various combinations of input variables were systematically evaluated for each dataset, and the optimal network architecture, including the number of hidden layers and neurons, was determined through a trial-and-error approach. For the agricultural soil dataset, the most suitable PTF model was based on %Sand, %Clay, %OC, BD and pH. The optimal ANN architecture (5–2–1), comprising five input neurons, one hidden layer with two neurons, and one output neuron, converged after 54,173 training iterations and achieved a classifier error of 0.52% (Figure 13a). For the forest soil dataset (Figure 13b), the optimal PTF incorporated seven input variables (%Sand, %Silt, %Clay, %OC, BD, pH and CO₃²⁻). The best-performing ANN architecture (7–4–1–3–1), consisting of three hidden layers with 4, 1, and 3 neurons, respectively, converged after 57,530 training iterations with a classifier error of 0.05%. For the pooled dataset (Figure 13c), the optimal PTF was developed using four input variables (%Clay, %OC, BD and pH). The optimal ANN architecture (4–4–4–2–1), comprising three hidden layers with 4, 4, and 2 neurons, converged after 39,600 training iterations and yielded a classifier error of 0.46%.

The predictive performance of the ANN models was evaluated using RMSE, MAE, R², NSE, and AIC for both the training and testing datasets (Table 8). For the agricultural soil dataset, the ANN model achieved R² values of 0.67 and 0.57, with corresponding RMSE values of 0.14 and 0.12 during the training and testing phases, respectively. For the forest soil dataset, the model attained R² values of 0.92 and 0.91, and RMSE values of 0.04 and 0.05 for the training and testing phases, respectively. Similarly, the pooled dataset yielded R² values of 0.75 and 0.70, with RMSE values of 0.09 and 0.13 during the training and testing phases, respectively. Overall, the ANN models demonstrated satisfactory predictive performance across all datasets (0.50 < R² < 0.75), with the forest soil model exhibiting the highest accuracy. Scatter plots comparing the predicted and measured MWD values for the training and testing datasets are presented in Figures 14a–f.

Table 8. ANN model performance indicators for different land uses during ‘training’ and ‘testing’ phases.

The numerical values of model performance statistics are shown as: mean ± standard deviation, estimated from 50 repetitions of train-test split and 10-fold cross-validations. Bold values indicate best PTF model for a particular land-use pattern identified separately for training and testing phases.

4 Discussion

Comparative evaluation of the machine learning (ML) models for soil mean weight diameter (MWD) prediction across different land-use datasets revealed distinct model-specific strengths. The Random Forest (RF) model exhibited the highest predictive accuracy for the agricultural land dataset, achieving the lowest RMSE, MAE, and AIC values in both the training and testing phases. In contrast, the Artificial Neural Network (ANN) model demonstrated superior performance for the forest land dataset, while the Support Vector Machine (SVM) model provided the most accurate predictions for the pooled dataset. Overall, all machine learning models, except the Multiple Linear Regression (MLR), exhibited robust predictive capability, highlighting the effectiveness of non-linear approaches for estimating soil aggregate stability (SAS) or MWD across diverse land-use conditions.

The MLR model exhibited substantially lower predictive accuracy than the machine learning models for SAS estimation, corroborating the findings of 26Besalatpour et al. (60. Although MLR has been successfully applied for broad-scale prediction of soil properties, its performance generally declines at local scales owing to the greater spatial heterogeneity and complexity of soil attributes (84). The relatively limited dataset from a geographically confined region used in this study likely contributed to the poor performance of MLR. Furthermore, regression-based pedotransfer functions (PTFs) rely on pre-defined linear relationships between ‘predictor’ and ‘response’ variables (85), which restricts their ability to capture the complex non-linear interactions governing SAS.

In contrast, the ANN model does not require an explicit functional relationship between ‘input’ and ‘output’ variables, enabling it to effectively model non-linear patterns. Consequently, ANN consistently outperformed MLR, reducing the RMSE by 0.05, 0.06, and 0.08 for the agricultural, forest, and pooled datasets, respectively. Among all evaluated models, ANN achieved the highest predictive accuracy for the forest land dataset, in agreement with Shalmani et al. (43), who also reported the superiority of ANN over regression-based PTFs for predicting SAS in forest and pasture soils. The enhanced performance of ANN can be attributed to its ability to learn complex non-linear relationships (86), which are characteristic of the highly heterogeneous soils of the Himalayan region. Nevertheless, the predictive capability of ANN remains dependent on both the quality and the selection of input variables (87). Similar improvements in SAS prediction using ANN over conventional regression models have also been reported by Rivera and Bonilla (15) and Besalatpour et al. (26), highlighting the robustness of ANN for modeling complex soil processes (88).

In the present study, the RF model achieved the highest predictive accuracy for the agricultural land dataset, outperforming the results reported by Bouslihim et al. (56) for SAS prediction. Bouslihim et al. (56) evaluated RF and MLR models for predicting MWD using soil properties and remote sensing data in Morocco and demonstrated that both approaches produced satisfactory predictions, with RF showing superior performance. Similarly, Zeraatpisheh et al. (89) compared several ML models for predicting Soil Organic Carbon (SOC) and SAS indices, including water-stable aggregates (WSA), geometric mean diameter (GMD), and MWD, in northern Iran and reported that RF consistently outperformed the other models for SAS prediction, whereas SVM yielded the best results for SOC estimation.

The superior performance of RF observed in this study is consistent with its well-established ability to capture complex, non-linear interactions among environmental variables while maintaining high predictive accuracy and resistance to overfitting (12, 90). Supporting this, Nazeri et al. (12) reported that RF models integrating topographic, geological, land use/land cover, remote sensing, and soil variables achieved the highest validation accuracy for MWD prediction among several ML algorithms. Likewise, Bouslihim et al. (14) successfully mapped multiple MWD indices, including mechanical breakdown (MWDmb), fast wetting (MWDfw), slow wetting (MWDsw), and the overall mean (MWDmean), with RF, obtaining high predictive accuracy for both calibration (R² > 0.92) and validation (R² > 0.65) datasets.

The enhanced performance of RF can be attributed to its ensemble-learning framework, which effectively models complex predictor–response relationships while reducing model variance and overfitting. In addition, RF provides robust estimates of predictor importance, unlike conventional MLR, which typically retains only strongly correlated variables through stepwise selection and may exclude potentially informative predictors under multicollinearity (56, 91). Boateng et al. (92) further demonstrated that RF and Support Vector Regression (SVR) combine high computational efficiency with superior predictive performance relative to k-Nearest Neighbor (kNN) and Neural Network models. Collectively, these findings reinforce the robustness and generalizability of RF as a reliable approach for predicting SAS across diverse environmental conditions.

A few past studies have also reported about the superior performance of SVM over conventional regression-based models for predicting SAS (8, 60). Consistent with these findings, the Radial Basis Function (RBF)-kernel SVM outperformed other kernel functions evaluated in the present study, corroborating previous reports on the prediction of soil properties and SAS (8, 93, 94). Notably, SVM achieved the highest predictive accuracy for the pooled dataset, suggesting its strong generalization capability when trained on larger and more heterogeneous datasets. This performance can be attributed to the structural risk minimization principle of SVM, which enables effective learning of complex non-linear relationships while maintaining good generalization (93, 95). Nevertheless, the predictive performance of SVM remains highly dependent on appropriate hyperparameter optimization, particularly the Cost (C), Epsilon (ϵ), and Gamma (γ) parameters, as well as the representativeness of the training dataset (96).

The present findings are consistent with the growing adoption of Artificial Intelligence (AI) for rapid and scalable assessment of SAS. For example, Fajardo et al. (97) demonstrated the potential of smartphone-based imaging for efficient SAS estimation, whereas the present study employed ML models driven by readily available soil physico-chemical properties. The consistently high performance of RF, ANN, and SVM highlights the capability of ML models to generate accurate and reliable estimates of SAS, making these approaches well suited for rapid and cost-effective regional-scale soil quality assessment, erosion risk evaluation, and precision land management. Furthermore, recent advancements indicate that integrating ML models with synthetic data augmentation can mitigate data scarcity, improve model robustness, and further enhance predictive performance (20), thereby strengthening the applicability of AI-based PTFs for large-scale soil monitoring.

No single ML algorithm consistently performs best across all datasets and environmental conditions, as widely acknowledged in the literature. In the present study, however, the Support Vector Machine (SVM) model exhibited the highest predictive accuracy for the pooled dataset, indicating its strong generalization capability across diverse land-use conditions. The superior performance of the SVM model for the pooled dataset demonstrates its ability to capture the relationship between soil properties and MWD across both agricultural and forest land-use systems. Importantly, the model maintained meaningful predictive performance on the independent testing data despite the heterogeneity associated with these contrasting land-use conditions. These results provide evidence of the potential of SVM as a generalizable approach for MWD prediction under diverse soil and land-use conditions. However, because the present study represents a limited geographical extent and sample size, further validation using larger and geographically distributed datasets is needed to establish the transferability of the model across the broader Eastern Himalayan region.

5 Conclusion

The present study developed Machine Learning (ML)-based pedotransfer functions (PTFs) for estimating soil aggregate stability (SAS) across agricultural and forest land-use systems in a Sikkim Himalayan watershed. The results demonstrated clear land-use-specific differences in model performance. The Random Forest (RF) model achieved the highest predictive accuracy for the agricultural soil dataset using an optimized combination of input variables (%Silt, %Clay, %OC, BD, and pH), together with tuned hyperparameters (ntree = 300, and mtry = 5). In contrast, the Artificial Neural Network (ANN) model performed best for the forest soil dataset, utilizing %Sand, %Silt, %Clay, %OC, BD, pH, and CO₃²⁻ as predictors with an optimized architecture comprising three hidden layers (4, 1, and 3 neurons). These findings highlight the importance of adopting land-use-specific ML models to improve the accuracy and reliability of SAS prediction, as well as for undertaking critical land management decisions.

Notably, the Support Vector Machine (SVM) model with a Radial Basis Function (RBF) kernel achieved the highest predictive accuracy for the pooled dataset, demonstrating strong generalization across diverse land-use conditions. Overall, the ML models consistently outperformed the conventional Multiple Linear Regression (MLR) approach, confirming the advantages of non-linear techniques for SAS prediction. The developed ML-based PTFs, based on routinely measured soil physico-chemical properties, provide a rapid, scalable, cost-effective, and reliable alternative to direct SAS determination. Consequently, these models have considerable potential for regional-scale soil quality assessment, soil erosion risk evaluation, and the development of targeted soil conservation and sustainable land management strategies in the eastern Himalayan region.

The findings of this study have important implications for soil erosion mitigation and sustainable land management in mountainous regions, where steep slopes and intense rainfall make soil stability a critical concern. Therefore, accurate estimation of SAS can support the identification of erosion-prone areas and the development of targeted soil-and-water conservation strategies, thereby improving agricultural sustainability and ecosystem resilience. Beyond the eastern Himalayan region, this methodology is readily transferable to other mountainous landscapes with similar environmental conditions, offering a robust decision-support tool for evidence-based soil monitoring, soil conservation, and long-term sustainable land management.

6 Limitations and scope for future research

Although the developed machine learning-based pedotransfer functions demonstrated satisfactory predictive performance, several limitations should be acknowledged. First, the study was based on 126 topsoil samples (0–15 cm) collected during a single winter season from one sub-watershed of the eastern Himalayas, which may not fully capture seasonal variability, deeper soil horizons, or the spatial heterogeneity of the region. Second, only ten readily measurable physico-chemical soil properties were considered, whereas biological, mineralogical, topographic, climatic, and land management factors that also influence soil aggregate stability were not included. Third, external validation using independent datasets from other Himalayan landscapes was not conducted, limiting the broader applicability of the developed models.

Future studies should incorporate multi-season and multi-depth sampling across diverse physiographic settings to improve model robustness and transferability. Integration of terrain attributes, remote sensing products, climatic variables, soil biological indicators, and management practices may further enhance prediction accuracy. The application of advanced explainable artificial intelligence, ensemble and deep learning techniques, together with uncertainty quantification, should also be explored. Finally, coupling the developed soil aggregate stability models with watershed-scale soil erosion and sediment transport models would facilitate regional erosion risk assessment, sustainable land-use planning, and climate-resilient soil conservation strategies across the Himalayan region.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Layla M. San-Emeterio, Swedish University of Agricultural Sciences, Sweden

Reviewed by: Shahid Ul Islam, Baba Ghulam Shah Badshah University, India

Harendra Singh, University of Petroleum and Energy Studies, India

Footnotes

AI, Artificial Intelligence; AIC, Akaike Information Criteria; ANN, Artificial Neural Network; BD, Bulk Density; MAE, Mean Absolute Error; ML, Machine Learning; MLR, Multiple Linear Regression; MWD, Mean Weight Diameter; NSE, Nash-Sutcliffe Efficiency; %OC, Soil Organic Carbon Content; PTF, PedoTransfer Function; R², Coefficient of Determination; RF, Random Forest; RMSE, Root Mean Square Error; SAS, Soil Aggregate Stability; SOC, Soil Organic Carbon; SVM, Support Vector Machine; USDA, United States Department of Agriculture.