Surrogate Modeling Reveals Hidden Failure Risks in Complex Soil-Rock Composite Slopes

Geotechnical engineering researchers have demonstrated that parameter uncertainty in soil-rock composite slopes can produce a 2.767 percent probability of slope failure even when traditional deterministic evaluations calculate an apparently stable factor of safety exceeding 1.0. To overcome the severe computational bottlenecks and dimensional complexity associated with modeling heterogeneous open-pit mine slopes comprising overlying loose silty clay and underlying weathered mylonite bedrock, an engineering research team established a staged reliability evaluation framework. This methodology couples two-dimensional limit equilibrium screening with an explicit second-order response surface surrogate derived from quasi-three-dimensional finite-difference numerical simulations in FLAC3D. By evaluating 100,000 Monte Carlo sampling realizations against the calibrated algebraic surrogate model, the researchers demonstrated an equivalent reliability index of 1.916 under baseline parameter independence. The findings reveal that conventional mean-value stability evaluations mask sudden sliding collapse hazards, providing open-pit mining operations with an efficient protocol for probabilistic safety assessments.

The investigated open-pit mining slope possesses an excavated stepped morphology characterized by a toe elevation of 60 meters and a maximum vertical relief difference of approximately 142 meters. Geotechnical cross-sections expose a dual composite stratigraphy consisting of an uppermost silty clay layer with uneven thickness resting above fractured mylonite bedrock with undulating contact surfaces. Weathering processes have subdivided the underlying bedrock into an upper highly weathered zone and a deeper moderately weathered zone with contrasting mechanical integrity. Site engineering geological surveys established that the highly weathered mylonite has suffered extensive structural destruction, yielding fragmented drill cores occurring as blocky fragments with rock quality designation values ranging from 0 percent to 39 percent, an average of 4.06 percent, and an engineering classification of Grade V. Conversely, the deeper moderately weathered mylonite retains recognizable rock structure dissected by weathering fractures, exhibiting rock quality designations between 27 percent and 57 percent, an average of 40.91 percent, and a Grade IV classification based on comprehensive rock mass ratings.

Deterministic baseline parameters acquired from site investigation reports documented representative unit weights of 19.00 kilonewtons per cubic meter for the silty clay, 25.90 kilonewtons per cubic meter for the highly weathered mylonite, and 26.30 kilonewtons per cubic meter for the moderately weathered mylonite. Strata stiffness was characterized by elastic moduli of 0.010 gigapascals for clay, 0.500 gigapascals for highly weathered rock, and 0.612 gigapascals for moderately weathered rock, with Poisson ratios of 0.35, 0.28, and 0.23. Recognizing spatial heterogeneity in natural geological strata, the researchers applied coefficients of variation of 5 percent to unit weight, 10 percent to internal friction angle, and 25 percent to cohesion across all three geotechnical units. Because incorporating nine uncertain parameters simultaneously into numerical simulations creates severe computational bottlenecks, the team deployed an antecedent sensitivity screening procedure to isolate dominant variables.

To execute computationally inexpensive dimensionality reduction, the team performed sensitivity screening in the two-dimensional Slide software using the simplified Bishop limit equilibrium method with circular slip surface searches across 25 slices. Utilizing a three-level orthogonal array comprising 27 parameter combinations under natural static conditions, multivariate analysis of variance in SPSS Statistics evaluated the main effects of all nine parameters on slope safety factors. The statistical screening revealed that the cohesion, internal friction angle, and unit weight of the moderately weathered mylonite exerted dominant control over slope stability, exhibiting F-values of 182.511, 163.347, and 5.639 with respective p-values below 0.001, below 0.001, and equal to 0.030. These three variables contributed 49.70 percent, 44.48 percent, and 1.54 percent to total safety factor variation, collectively explaining 95.72 percent of the corrected sum of squares with omega-squared effect sizes of 0.494, 0.441, and 0.013. The remaining six parameters of the clay and highly weathered layer showed statistically insignificant influence and were fixed at mean baseline values.

For rigorous mechanical analysis of global slope deformation and failure, the investigators established a 1-meter-thick quasi-three-dimensional continuum model in FLAC3D under plane-strain lateral displacement constraints. Mesh discretization sensitivity was evaluated across coarse, medium, and fine numerical grids containing 4,180, 8,620, and 13,460 zones, which yielded strength reduction safety factor values of 1.164, 1.148, and 1.140, respectively. Advancing from the medium mesh with 10,197 gridpoints to the fine mesh increased zone count by 56 percent while shifting the safety factor by only 0.008, demonstrating diminishing mesh sensitivity.

Discontinuity surveys across the slope sector documented joint trace lengths of 0.4 to 3 meters with spacing densities between 4 and 10 joints per meter, confirming that localized joint sets were not kinematically oriented to trigger planar or wedge failures. This geological condition justified adopting an equivalent continuum Mohr-Coulomb constitutive model across all three strata to capture global rotational deformation. Numerical plastic shear strain increment contours generated in FLAC3D revealed that shear localization initiated at the slope toe within the moderately weathered mylonite before propagating progressively upward through the composite strata. The resulting continuous, curved shear band closely aligned in spatial geometry with the critical circular slip surface identified by the automated Bishop search, corroborating a deep-seated rotational failure mode.

Having verified the failure mechanism and screened the parameter space, the researchers implemented a three-factor, three-level Box-Behnken Design in Design-Expert software to construct a second-order response surface surrogate. The experimental design matrix comprised 12 edge-midpoint runs and five replicated center-point evaluations across the three retained parameters of the moderately weathered mylonite. Multivariate regression analysis generated a quadratic polynomial surrogate that achieved an F-value of 805.59 with a p-value below 0.0001, confirming exceptional statistical significance. Model diagnostics demonstrated a coefficient of determination of 0.9990, an adjusted R-squared of 0.9978, and a predicted R-squared of 0.9846, with the narrow 0.0132 discrepancy reflecting outstanding predictive stability. The surrogate exhibited a low standard deviation of 0.0047, a coefficient of variation of 0.41 percent, and an adequate precision ratio of 110.263, indicating high signal-to-noise performance. Zero pure-error variance across five replicated FLAC3D center points reflected the deterministic nature of continuum calculations, while Box-Cox diagnostics confirmed that untransformed safety factors accurately captured response surface behavior.

Topological examination of three-dimensional response surfaces and two-dimensional contour plots clarified complex nonlinear coupling among the controlling geotechnical parameters. Cohesion and internal friction angle displayed approximately parallel linear contour distributions without significant interaction, matching an interaction p-value of 0.6079 in analysis of variance testing. Conversely, the interaction between cohesion and unit weight proved statistically significant with an F-value of 18.45 and a p-value of 0.0036, reflecting how gravitational self-weight alters sliding forces and normal effective stresses along shear surfaces. Calculating the local marginal effect of unit weight at the center point yielded a partial derivative of negative 0.00737, confirming that increased rock weight diminishes stability despite an uncentered positive first-order algebraic coefficient. Independent generalization testing evaluated 20 validation combinations generated using Latin Hypercube space-filling principles, including 12 near-limit-state points and 8 global parameter configurations. Comparing surrogate predictions against direct FLAC3D runs across these 20 cases yielded a root mean square error of 0.0142, a mean absolute error of 0.0114, and maximum relative errors within 2.55 percent.

The validated explicit response surface equation was coupled with a vectorized Monte Carlo solver in Python to compute probabilistic reliability metrics across 100,000 sampling iterations. To prevent physically unrealistic realizations or extrapolation beyond calibrated design ranges, the three parameters were modeled as truncated normal distributions spanning 106.73 to 177.88 kilopascals for cohesion, 28.91 to 35.33 degrees for friction angle, and 24.99 to 27.62 kilonewtons per cubic meter for unit weight. Under baseline parameter independence, 2,767 realizations produced safety factors below 1.0, defining a failure probability of 2.767 percent and an equivalent reliability index of 1.916. The simulated distribution exhibited a mean safety factor of 1.144, a standard deviation of 0.075, a skewness of negative 0.052, and an excess kurtosis of negative 0.513. Tracking running failure estimates confirmed sample convergence with a binomial standard error of 0.0519 percentage points and an exact 95 percent Clopper-Pearson confidence interval of 2.666 percent to 2.871 percent. Furthermore, independent limit equilibrium verification comparing nine-variable and three-variable Monte Carlo models yielded failure probabilities of 6.510 percent and 6.630 percent, demonstrating that removing six minor variables produced only a 0.120 percentage point difference.

To explore dependence between cohesion and internal friction angle without sufficient paired field measurements, the researchers executed a Gaussian Copula sensitivity analysis across correlation scenarios of 0, negative 0.2, negative 0.4, and negative 0.6. The simulated failure probability systematically dropped from 2.767 percent to 1.566 percent, 0.703 percent, and 0.137 percent, showing that negative correlation reduces the likelihood of simultaneously low strength values. Methodologically, the staged framework eliminated repetitive continuum simulations, requiring only 17 design runs and 20 validation calculations in FLAC3D before conducting 100,000 Monte Carlo evaluations. The investigators emphasized that the current reliability assessment is limited to natural static conditions, as groundwater variations, rainfall infiltration, and earthquake ground motions were excluded from the analysis. These external loading processes can alter effective stresses and shear strength, prompting the research team to propose extending the dimensionality reduction response surface framework to coupled hydro-mechanical and seismic stability evaluations in future open-pit mine engineering studies.

Subject of Research: Probabilistic reliability analysis and surrogate modeling of soil-rock composite slopes

Article Title: Reliability analysis of soil-rock composite slopes based on dimensionality reduction response surface method

Article References: Feng, L., Chen, S., Sun, J., Jia, B., Zhao, B., Gu, Y., Li, Z., Zhang, K., Xue, Y., & Zhou, Y. (2026). Reliability analysis of soil-rock composite slopes based on dimensionality reduction response surface method. Frontiers in Materials, 13, Article 2007637. https://doi.org/10.3389/fmats.2026.2007637

Image Credits: AI Generated

DOI: 10.3389/fmats.2026.2007637

Keywords: slope stability, soil-rock composite slopes, response surface method, Monte Carlo simulation, limit equilibrium method, FLAC3D, geotechnical engineering, reliability analysis, Reliability, analysis, soil-rock, composite