Development of machine learning surrogate models for slope stability prediction using AI techniques: a case study of the Meizhou landslide

The results of sieve analysis after conducting the tests at the soil mechanics laboratory indicates that the material at the site, consisting of gravel, sand, silt, coarse clay, and medium clay particles derived from completely weathered granite soil, is characterized by three-dimensional angular particles. These particles exhibit minimal variation in the diameters of their major and minor axes. Additionally, the large frictional resistance between the particles, attributed to their irregular geometry and lack of rounded edges, significantly influences their mechanical behavior. Laboratory testing was conducted to determine the plastic limits of three soil samples with particle sizes smaller than 0.425 mm, following the relevant standards and procedures. The plastic limit moisture contents for samples 1, 2, and 3 were found to be 30.74%, 35.76%, and 24.36%, respectively, while the corresponding liquid limit moisture contents were 84.59%, 103.74%, and 41.32%. These results indicate that all three soil samples, with particle sizes below 0.425 mm, exhibit clay-like properties based on their consistency limits. Based on the Engineering Classification Standard for Soil (GB/T 50145 − 2007), soil samples 1, 2, and 3 are classified as sandy soils, specifically clayey sand and fine-grained sand. Comprehensive measurements of completely weathered granite soils from Hong Kong⁶⁷ reveal that the original material exhibits a total bulk density of 16–21 kN/m³, a dry bulk density of 14–19 kN/m³, an effective internal friction angle ranging from 35° to 44°, an effective cohesion of 5–15 kPa, and a permeability coefficient between 10^{− 5} to 10^{− 7} m/s. For completely weathered granite fill compacted on-site, the total bulk density increases to 19–21 kN/m³, with a dry bulk density of 15–19 kN/m³. The effective internal friction angle ranges from 38° to 42°, effective cohesion decreases to 0–5 kPa, and permeability ranges from 10^{− 6} to 10^{− 7} m/s. These parameters illustrate distinct mechanical and hydraulic behaviors depending on the soil’s degree of weathering and compaction.

A slope stability model of the Meizhou Expressway landslide is developed using GeoStudio, a leading geotechnical engineering software suite. In parallel, AI-assisted coding frameworks, such as ChatGPT-MATLAB integration, have recently been employed to streamline geotechnical model development and enhance automation in numerical simulations⁶⁸. The model incorporates the geological, hydrological, and material properties of the slope to simulate failure mechanisms under various conditions, including heavy rainfall. GeoStudio is chosen for its robust capabilities in slope stability and seepage analysis, offering tools like SLOPE/W for evaluating stability and SEEP/W for modeling water infiltration and pore pressure changes. Its advanced features enable accurate simulations of complex geotechnical scenarios, making it ideal for understanding the intricate processes leading to landslides like the one at Meizhou. Figure 4 shows the slope model developed on GeoStudio.

Slope stability analysis with water level variation

The trapped rain water to the side of the slope as mentioned in Fig. 2 are penetrating in the soil below the road and the water level is changing time to time. Keeping this point in consideration, the slope is analyzed to find out FoS values with varying water level. The overall height from point A to F is 30 m. During the analysis, the top height considered is 25 m and lower height as 13 m. The total number of analyses performed are seventeen with each variation of 1 m interval to find out the FoS at each and every interval. Table 1 shows the mechanical properties of the slope and FoS vales with the variation of the water level. The analysis types are Morgenstern-Price. The staged pseudo-static analysis option is selected as effective strength. The pore water pressure conditions are selected as piezometric surface and the unit weight of water is 9.81 kN/m³ as constant. Figure 5 shows the graph between the FoS with the water level variation.

The graph shows in Fig. 5 illustrates the linear relationship between water level (m) and the FoS. It implies 98.23% of the variation in FoS is explained solely by changes in water level. Such a high R² suggests a very reliable predictive relationship under the modeled conditions.

Machine learning analysis for FoS

The objective in this part is to model the relationship between FoS and various geotechnical parameters using ML, specifically multiple linear and polynomial regression, ensuring R² > 0.90. The dataset consists of 249 data points (Annex-1) derived through model analysis using GeoStudio. The independent variables are Effective Cohesion (kPa), Effective Friction Angle (φ), Unit Weight (kN/m³), Surcharge Load (kN/m²), Water Table Level (m) and the dependent variable is FoS. FoS values with the variation of Cohesion and keeping all other parameters as constant are mentioned in Table 2, while for all other parameters, the FoS values are mentioned in Annex-1.

The variation of the FoS in response to different geotechnical parameters, as illustrated in Fig. 6, demonstrates distinct trends that reflect the underlying mechanics of soil stability. An increase in cohesion results in a progressive enhancement of FoS, indicating a near-linear positive correlation due to the direct contribution of cohesive forces to shear strength. The response of FoS to changes in the friction angle is markedly non-linear, exhibiting an exponential growth pattern especially beyond 60° highlighting the dominant influence of interparticle resistance on slope stability. In contrast, the relationship between unit weight and FoS displays a diminishing return behavior; FoS rises rapidly at lower densities but plateaus as unit weight continues to increase, suggesting a balance between beneficial normal stress and adverse self-weight effects. Interestingly, FoS appears unaffected by variations in surcharge load under the studied conditions, which may indicate that the influence of external loading was either minimal or counteracted by other factors in the analysis. Collectively, these patterns underscore the complex, often non-linear interactions between strength parameters and stability, emphasizing the need for robust modeling frameworks when predicting FoS under variable geotechnical conditions.

Figure 7 presents Kernel Density Estimation plots comparing the distributions of six key geotechnical parameters before and after data cleaning. 41 extra data points are removed to ensure the dataset is free from influential outliers, statistically sound for analysis, uniformly filtered across features and reflective of realistic geotechnical conditions.

For each parameter, the original dataset (blue line) and the cleaned dataset (red line) are overlaid to visualize the impact of the cleaning process on data distribution. In all cases, the cleaned data exhibit sharper and more peaked density curves, indicating a reduction in variability and outliers. This sharpening effect suggests the removal of anomalous values and a more centralized, representative dataset, which improves the reliability of further probabilistic or stability analyses. Notably, the distribution of FoS becomes more concentrated around a mean value, which is critical for consistent slope stability assessments.

Random forest-based stability prediction model

Random Forest is a powerful ensemble ML technique that constructs multiple decision trees and aggregates their predictions to enhance accuracy and reduce overfitting. In geotechnical engineering, RF is particularly effective for modeling complex, non-linear relationships between soil properties and performance indicators such as the FoS. Its ability to capture variable interactions and assess feature importance makes it highly suitable for parametric stability analysis. RF models are typically developed using data science platforms such as Python (scikit-learn), R, or MATLAB, with Python being the most widely used due to its flexibility and integration with statistical libraries. In this study, RF was employed to develop a robust predictive correlation between FoS and key soil parameters including cohesion, friction angle, unit weight, surcharge load, and water level. Using the full Annex-1 data set and a 300-tree Random-Forest regressor, the ensemble’s mean prediction can be written as the following surrogate equation, obtained by averaging individual tree splits and fitting a second-order surface to the resulting partial-dependence curves in Eq. (1):

Table 3 shows the mean importance of each parameter in this correlation.

The Random-Forest analysis confirms that FoS is governed primarily by material strength parameters cohesion and friction angle jointly explain over 60% of the variance, with a positive first-order contribution and a modest interaction term, indicating that steeper friction gains are realized in low-cohesion soils. Water level exerts the strongest negative influence: the quadratic term shows FoS degradation accelerates once the piezometric surface approaches the failure plane, reflecting reduced effective stress. Unit weight has a secondary negative effect because greater self-weight adds to driving forces, but its interaction with water level slightly offsets that penalty when buoyancy is significant. Surcharge is the least influential variable under the tested ranges, producing a near-linear, weak reduction in stability. The residual quadratic in φ captures the plateau observed beyond 60°, where additional friction yields diminishing returns because the mobilization factor approaches unity. Collectively, the surrogate equation reproduces the ensemble’s non-linear response.

Coefficient of Determination, R² = 0.75, means 75% of the variation in the actual FoS values is explained by the model.

Root Mean Square Error, RMSE ≈ 0.25.

RMSE measures the average magnitude of prediction error. Here, the average error between the model-predicted FoS and the actual FoS is about 0.25 units. Since typical FoS values range from 1.0 to 2.0, this is a very small error, suggesting the model is highly precise.

Equation (2) presents the formula for the Pearson Correlation Coefficient (rₓ_γ). It quantifies the linear relationship between two variables, calculated as the covariance of x and y divided by the product of their standard deviations.

The Pearson correlation matrix is employed to examine how two continuous variables are linearly related. It utilizes a statistical measure known as Pearson’s correlation coefficient (r), which takes values from − 1 to + 1. An r value near + 1 suggests a strong direct relationship, while a value near − 1 reflects a strong inverse relationship; values close to zero imply a weak or nonexistent linear link. In the field of geotechnical engineering, especially in data-driven studies like the present one, this matrix plays a crucial role in exploring the extent of correlation between various input parameters and the FoS. Such analysis helps in determining key contributing variables, refining the feature selection process, and providing deeper insights into the behavior of models such as RF. Table 4 presents the Pearson correlation matrix for the analysis based on RF method.

Table 4 illustrates the Pearson correlation matrix constructed for five geotechnical variables, cohesion (c), friction angle (φ), unit weight (γ), surcharge load (q), water level and their linear association with the FoS. Among all variables, the friction angle (φ) exhibits the strongest positive correlation with FoS (r = 0.838), indicating a dominant linear contribution to slope stability. Cohesion (c) also shows a moderate positive correlation with FoS (r = 0.259), while unit weight (γ), surcharge load (q), and water level display negative correlations of −0.105, −0.161, and − 0.196 respectively, suggesting that increases in these parameters could unfavorably impact FoS. Notably, the water level has a substantial negative correlation with cohesion (r = −0.623), hinting at the degradation of shear strength under higher moisture conditions.

Friction Angle (φ) – Positive effect: A higher friction angle increases the shear resistance along potential failure surfaces, enhancing slope stability and therefore increasing the FoS.

Cohesion (c) – Positive effect: Greater cohesion strengthens the bonding between soil particles, contributing to overall shear strength and improving slope resistance to failure.

Unit Weight (γ) – Negative effect: An increase in unit weight results in a higher self-weight of the soil mass, which amplifies driving forces and thereby tends to reduce the FoS.

Surcharge Load (q) – Negative effect: Additional surcharge increases the external loading on the slope, intensifying the driving forces acting on potential failure planes, which lowers the FoS.

Water Level (WL) – Negative effect: Elevated water levels raise pore water pressures, reducing effective stress and shear strength in the soil, thus negatively impacting slope stability and decreasing the FoS.

These insights align with the mean feature importance rankings from the RF model, emphasizing that while correlation highlights linear dependencies, it complements model-based importance metrics by offering transparency into variable interactions and helping validate the physical relevance of the predictive features.

Ordinary least squares (OLS) regression analysis method

OLS regression is a widely applied statistical approach used to establish a linear relationship between a response variable and one or more predictors. It works by identifying the line that minimizes the total of the squared differences between the observed values and the values estimated by the model. This approach provides the most accurate linear fit across the dataset. OLS relies on certain assumptions, including linearity between variables, constant error variance (homoscedasticity), and the independence of residuals. Due to its straightforward implementation, clear interpretation, and effectiveness in capturing trends, OLS remains a popular tool in disciplines such as geotechnical analysis, environmental studies, and economic forecasting. An empirical Eq. (3) was developed using the OLS regression method to quantify the relationship between the FoS and the selected geotechnical parameters.

XGBoost (extreme gradient boosting) method

To further increase the R² value and minimize the RMSE in predicting the FoS, a more sophisticated and in-depth analysis of the Annex 1 dataset is required. While traditional methods like OLS and RF offer baseline and moderately complex modeling capabilities, they may fall short in capturing the full range of nonlinear interactions and feature dependencies present in geotechnical data. Therefore, to enhance predictive accuracy and develop a more reliable correlation, we will employ the XGBoost method. XGBoost is known for its superior performance due to its ability to handle nonlinearities, incorporate regularization to prevent overfitting, and manage feature interactions automatically through ensemble learning. Unlike simpler regression models, XGBoost constructs decision tree ensembles in a stage-wise manner, allowing it to iteratively correct residual errors and achieve high model fidelity. Its robustness, scalability, and proven track record in regression tasks make it an ideal choice for developing a high-performance correlation from the Annex-1 data. Equation (4) presents the correlation developed using XGBoost method.

Where, R² = 0.88 and RMSE = 0.12.

This newly developed second-order polynomial correlation, fitted to the XGBoost model output, demonstrates exceptional predictive capability with an R² value of 0.88 and an RMSE of only 0.12. Compared to the previously derived RF and OLS models, which had R² values of 0.75 and 0.83 respectively, this surrogate provides a significantly more accurate approximation of the FoS. The enhanced performance is attributed to the XGBoost model’s ability to capture complex nonlinear interactions, which are then preserved and symbolically represented through the polynomial regression. This makes the new correlation both interpretable and highly precise, offering a robust tool for geotechnical design and analysis.

It is evident that none of the model predictions perfectly align with this ideal line. The RF and OLS predictions generally overestimate the FoS values across the range, while the XGBoost model exhibits significant non-linear deviations, especially at higher values. These discrepancies indicate that the derived correlation equations require further adjustment and calibration to improve their accuracy and bring the model predictions into closer agreement with the actual FoS values ultimately achieving better alignment with the 1:1 reference line. This is essential for improving model reliability in geotechnical design and safety analysis. The adjusted linear correlations 5, 6 and 7 shown below were developed to improve the agreement between predicted and actual FoS values by calibrating the original model outputs using simple linear regression post-processing. Specifically, the predicted FoS values from the RF, OLS, and XGBoost models were individually regressed against the corresponding actual FoS values from the dataset to establish direct linear mapping equations of the form; FoS_actual = a⋅FoS_pred + b. Following are the adjusted correlations for all the three methods:

1. Random Forest (RF):

2. OLS Regression:

3. XGBoost (Polynomial):

This calibration step effectively adjusts the model outputs to better reflect the real trend and scale observed in the actual measurements. The resulting equations significantly improve the alignment with the ideal 1:1 line, as evidenced by the high coefficients of determination (R² values exceeding 0.95) and low RMSE values as shown in Fig. 8.

Table 5 shows the details of R² and RMSE for all the three methods, and the Fig. 9 shows the percent improvement in R² and reduction in RMSE.

The bar graph presents a comparative analysis of the performance of three regression methods, RF, OLS, and XGBoost based on their predicted and adjusted R² and RMSE values. The upper panel illustrates that R² values increase progressively from RF to XGBoost, indicating improved model fit, with adjusted R² values showing further refinement, particularly for the XGBoost method which achieves the highest value of 0.96. The lower panel displays a corresponding decrease in RMSE values, demonstrating enhanced prediction accuracy after adjustment. Notably, the adjusted RMSE for RF reduces drastically to 0.010, reflecting a significant calibration improvement. This graphical representation highlights the superior predictive capability and robustness of the XGBoost model, followed by OLS and RF, especially after post-processing or model tuning. The RF method shows the greatest improvement in R², increasing from 0.75 to 0.90. This significant enhancement reflects the strong effect of post-processing or model adjustment in boosting the model’s explanatory power. Although XGBoost achieves the highest final R² value (0.96), its improvement margin is smaller (from 0.88 to 0.96), indicating that it was already well-calibrated in its raw prediction. Similarly, RF exhibits the largest decrease in RMSE, dropping from 0.25 to 0.10, an improvement of 0.15. This suggests a notable reduction in prediction error after adjustment. XGBoost, while having the lowest final RMSE (0.040), shows a relatively smaller improvement (0.08), again due to its already optimized performance prior to adjustment. RF benefits the most from adjustment in terms of both R² and RMSE. Similarly, XGBoost demonstrates the best absolute performance with highest R² and lowest RMSE, but shows less relative improvement.