The stability forecasting of embankment dams during construction represents a critical challenge in geotechnical engineering, where structures are at their most vulnerable to complex, multi-hazard failure mechanisms. This was catastrophically illustrated by the large-scale sliding event at the Megech Dam in Ethiopia, a case study that emphasizes the dire need for models that can not only predict failure but also quantify the confidence of their predictions. In our preceding study [1], we successfully developed a hybrid Finite Element Method (FEM) with an Artificial Neural Network-Long Short-Term Memory (ANN-LSTM) model, achieving high predictive accuracy (R² = 0.94) for forecasting rainfall-induced seepage failures. This framework effectively bridged physical principles with data-driven pattern recognition [1]. However, a significant limitation remained: the model produced deterministic point forecasts, lacking any formal quantification of predictive uncertainty. For operational early-warning systems, a single-valued prediction without a confidence interval offers limited practical utility for risk-based decision-making under the inherent uncertainties of geotechnical systems [2] .

The field of deep learning has increasingly addressed this critical gap through sophisticated Uncertainty Quantification (UQ) methods, now considered essential for deploying reliable models in safety-critical applications. Techniques such as Monte Carlo (MC) Dropout provide a practical Bayesian approximation to capture epistemic uncertainty (model uncertainty due to limited data), while Mixture Density Networks (MDNs) directly model aleatory uncertainty (inherent noise in the data) by predicting the parameters of an output distribution [3]. The transformative value of moving from deterministic to probabilistic forecasting is evident across disciplines, from structural engineering, where it quantifies uncertainty in damage identification [4], to broader computational fields where multi-fidelity approaches are enhancing the tractability of strong UQ [5].

However, the transfer and application of these advanced probabilistic deep learning frameworks to the specific domain of embankment dam stability forecasting remains limited (see for example [6] and [7]). While our previous work and others have established the value of hybrid modeling, most dam safety applications continue to produce deterministic outputs [8]. This severely limits the development of confidence-based early-warning protocols, which require not just a prediction of failure, but a quantified measure of confidence in that prediction to guide proactive interventions, such as the slope trimming and drainage measures implemented post-slide at Megech [9] .

Building directly upon the validated hybrid FEM-ANN-LSTM architecture established in our prior work [1], this research bridges the critical gap between accurate forecasting and actionable risk assessment. We evolve the deterministic framework into a strong, uncertainty-quantified predictive system. The primary objectives are to: (1) enhance the hybrid model by rigorously disentangling and quantifying both epistemic and aleatory uncertainties using MC Dropout and an MDN output layer, and (2) demonstrate its superior performance over deterministic benchmarks by achieving a high Prediction Interval Coverage Probability (PICP) on the historical Megech Dam behavior data. By embedding uncertainty directly into the forecasting workflow, this work enhances interpretability and provides a more decision-ready basis for early-warning systems in high-risk civil infrastructure.

The geotechnical variability of the clay core was statistically characterized, with grain size distributions found to span 0.4–26.8% gravel, 11.0–82.4% sand, 13.4–86.7% silt, and 48.7–89.5% clay. Atterberg limits were determined to range from LL = 31.0–89.5%, PL = 15.9–47.1%, and PI = 16.8–56.3%, while permeability was measured between 1.45×10⁻⁸ and 1.17×10⁻⁵ cm/sec. Free swell tests yielded values of 50–168%, the majority of which exceeded the 90% high-risk threshold.

Key climatic drivers were identified through Principal Component Analysis (PCA) applied to the 13-year rainfall record, reducing the dataset to dominant climate modes. This transformation normalized variability and isolated the multi-year and seasonal patterns most relevant to dam stability, consistent with established hydro-climatic methodologies [12], [13] and [14]. .

Accordingly, the PCA variance was primarily explained by two components: PC1 accounted for 68% of the variance, representing the primary precipitation pattern, while PC2 accounted for 24%, capturing seasonal temperature influences. This analysis was conducted on a dataset spanning from 2008 to 2021, comprising 156 monthly observations.

The temporal progression of settlement was captured through instrumentation monitoring across four measurement sets ( Table 3). The resulting time-series data enabled the identification of distinct behavioral phases through temporal correlation analysis, a fundamental technique for interpreting geotechnical monitoring data [15].

Furthermore, to strongly identify the primary environmental drivers from the complex meteorological record, we derived climate principal components (Table 2) from a 13-year raw rainfall dataset using Principal Component Analysis (PCA). This dimensionality reduction technique is well-established for isolating dominant climate patterns from noisy, long-term data, with advanced implementations like Rotated Spectral PCA (rsPCA) proving effective for identifying dynamical modes in climate systems [16].

The application of such PCA-based methods is foundational in geotechnical and hydrological studies for linking climate forcing to engineering responses [17], and the core principles are further detailed in foundational literature on the method [18].

This integrated multi-source dataset was engineered to create a comprehensive feature space for machine learning by fusing three primary data streams: field and laboratory geotechnical data, climate principal components derived from a 13-year meteorological record and instrumentation monitoring that characterized geotechnical behavior through settlement progression across four measurement sets. The application of temporal correlation analysis to this fused dataset subsequently enabled the identification of distinct behavioral phases, forming a rich, multi-dimensional foundation for predictive modeling. The fusion of these data streams (Table 4) created a unified feature set, which was mapped to defined risk categories (Table 5) to facilitate supervised learning for probabilistic forecasting, following paradigms for risk-informed geotechnical modeling as in the works of [19].

Analysis of the instrumentation data revealed a clear triphasic failure mechanism, a pattern consistent with the progressive failure of geotechnical-structures under evolving hydraulic and mechanical stresses (see for example, [20]). We explicitly formulated this physical sequence not merely as a post-failure description, but as a foundational constraint for the AI model, ensuring its predictions are grounded in recognizable geotechnical behavior and enhancing the interpretability of its outputs.. The registered progression in failure is described in Fig. 2.

This initial phase manifested as a near-linear settlement trend, attributable to secondary consolidation (creep) and the development of desiccation cracks during prolonged construction delays. The behavior is modeled by the linear function in Eq. (1):

Where S₁(t) is the cumulative settlement at time t, β₀ represents the immediate deformation, β₁ is the constant rate of long-term creep, and ε is the measurement error.

This critical phase captured the system's transition to impending failure, characterized by nonlinear acceleration triggered by the saturation of the foundation and the destabilizing toe excavation. The response is modeled by the exponential function in (2).

The exponential growth parameter γ was identified as a key instability indicator, directly quantifying the rate of strength loss, and was subsequently calibrated as the core metric for the early-warning system.

Stabilization (2022–2023). This final phase represents the period of kinematic stabilization following the major translational slide, where the slope reached a new, though failed, equilibrium geometry under residual shear strengths.

Building upon a validated deterministic hybrid FEM-ANN-LSTM baseline [1], this research introduces a probabilistic framework to quantify predictive uncertainty, which is critical for risk-based decision-making.

To transition from deterministic to probabilistic forecasting, we implemented a dual-strategy UQ approach:

To quantify uncertainty in the model itself (e.g., from limited data), we enabled dropout layers during inference. For each input, stochastic forward passes were performed. The epistemic variance was calculated as (3):

Where an individual prediction and is the mean prediction across all samples.

The probabilistic forecasting capability of the hybrid model was evaluated by plotting the temporal settlement predictions against the actual instrumentation data, as shown in Fig. 4. The model's mean prediction was visualized alongside the 68% and 95% confidence bands, which were derived from the total predictive uncertainty (). These intervals were calculated by combining the epistemic uncertainty from the MC Dropout with the aleatory uncertainty from the Mixture Density Network output. The resulting visualization demonstrates that the model not only accurately captured the trend of the actual settlement but also successfully quantified its confidence, with the observed data points predominantly falling within the predicted uncertainty bounds, particularly during the critical acceleration phase.

To diagnostically assess the calibration of the model's uncertainty estimates, a scatter plot of the predicted uncertainty against the absolute prediction error was generated, as presented in Fig. 5. This analysis was performed to verify that the model's self-reported uncertainty was a reliable indicator of its actual accuracy. A positive correlation between these two quantities confirms that the model correctly assigned higher uncertainty to predictions that had a larger error, a hallmark of a well-calibrated probabilistic model. This validation was a critical step in establishing that the uncertainty quantification could be trusted for operational decision-making, as it demonstrates the model's ability to signal its own potential inaccuracies in advance.

The epistemic uncertainty inherent in the hybrid model was quantified using Monte Carlo (MC) Dropout, a technique that provides a practical approximation of Bayesian inference. During the inference stage, dropout layers were retained with probabilities of p = 0.3 in LSTM layers and p = 0.5 in dense layers, and for each input sequence, T = 100 stochastic forward passes were performed to sample from the approximate posterior distribution of the model weights. The epistemic uncertainty () was then calculated as the standard deviation across these stochastic predictions, following Eq. (3). This procedure was implemented to capture the model's uncertainty due to limited data and architectural choices. The resulting analysis is presented in Fig. 6, which visually demonstrates the methodology's effectiveness: panel (a) displays the distribution of the multiple forward passes, illustrating the range of possible outcomes; panel (b) charts the temporal evolution of the quantified epistemic uncertainty, clearly showing its significant increase during the identified high-risk precursor periods (Weeks 25–30 and 75–80); panel (c) contrasts this probabilistic output with a single deterministic forecast, highlighting the added information; and panel (d) decomposes the total uncertainty at critical junctures, confirming that epistemic uncertainty became the dominant component leading up to the failure event. This approach ensured that the model's confidence was explicitly quantified, moving beyond a single point estimate to a more informative and reliable probabilistic forecast.

The distributions of MC-Dropout predictions at four representative monitoring weeks—10, 30, 75, and 80—provide clear insight into how uncertainty evolves across the deformation process. At Week 10, the predictive ensemble is narrow and tightly clustered, indicating high model confidence during the stable phase. By Week 30, the distribution widens slightly, reflecting the onset of mild precursor activity. A pronounced increase in spread is observed by Week 75, signaling accelerated deformation and a reduction in the model’s confidence as the behaviour departs from previously observed patterns. By Week 80, just before failure, the distribution becomes substantially broader with significantly higher variance, demonstrating the dominance of epistemic uncertainty during the pre-failure stage. Together, these snapshot distributions illustrate how the stochastic forward passes progressively diverge as instability develops, offering a clear and interpretable measure of emerging risk.

As shown in Fig. 7, the temporal evolution of epistemic uncertainty exhibits distinct shifts that align with the system’s stability state. During the early and midstages, the uncertainty remains relatively low and fluctuates within a narrow band, reflecting a well-constrained model response. However, pronounced increases emerge near the precursor windows—highlighted around Weeks 25–30 and 75–80—indicating reduced model confidence as the embankment transitions toward a less stable condition. This behavior demonstrates that the model effectively captures early signs of destabilization by quantifying the growth of epistemic uncertainty over time.

The probabilistic and deterministic forecasts overlap closely during the stable and moderate deformation phases. After approximately Week 70, when deformation accelerates toward failure, the probabilistic model provides a widening 95% prediction interval that reflects increasing uncertainty, whereas the deterministic model offers only a single trajectory. This divergence highlights the advantage of the probabilistic approach, which not only captures the trend but also quantifies confidence, providing more reliable early-warning information for decision-making.

The Mixture Density Network (MDN) architecture was implemented to explicitly capture the aleatory uncertainty inherent in the geotechnical monitoring data. As demonstrated in Fig. 8, the MDN outputs the full parameters of a 3-component Gaussian Mixture Model—specifically the mixture weights (), means (), and variances (σₖ²)—enabling the model to represent complex, multi-modal probability distributions for the Factor of Safety (Table 6). This approach moves beyond single-point estimates to provide a complete probabilistic characterization of the forecast. The total predictive uncertainty was synthesized by combining the epistemic uncertainty from MC Dropout with the aleatory uncertainty derived from the GMM variance, following. As quantified in Table 6, this decomposition revealed that epistemic uncertainty became the dominant component (72%) during the critical pre-failure period (Week 78), while aleatory uncertainty remained relatively constant. The resulting 95% prediction intervals, calculated as ± 1.96·, provided rigorously calibrated confidence bounds that expanded appropriately during high-risk periods.

The model was trained on 110 weeks of pre-failure data using a combined loss function (4).

Where: is the Negative Log-Likelihood loss for the MDN and is a FEM-based constraint that penalizes predictions violating fundamental geotechnical principles. The hyperparameter λ = 0.1 was determined through cross-validation to balance data quality and reduce physically unlikely predictions by 61% compared to the pure data-driven approach while maintaining minimal validation loss (0.136). Training employed the Adam optimizer with an initial learning rate of 0.001, automatically reduced by upon validation loss plateau at epochs 45 and 82. Early stopping with patience of 25 epochs prevented overfitting, with the final model converging efficiently in 150 epochs (Table 2). Cross-validation across 5 folds demonstrated consistent performance ( = 0.136 ± 0.003), confirming the strength of the selected hyperparameters and training strategy.

In this regard the entire optimization process is visualized in, which systematically breaks down the training efficacy: panel (a) quantifies the impact of the physics weight (λ), showing the 61% reduction in physically unlikely predictions at the optimal value of λ = 0.1, a balance between data fidelity and physical plausibility consistent with principles of physics-informed neural networks [21]. Panel (b) details the training dynamics, tracing the stable convergence of the total loss and its components, and explicitly marks the adaptive learning rate reductions at epochs 45 and 82 that refined the final convergence. The model's reliability is further substantiated in panel (c) through a comprehensive 5-fold cross-validation table, confirming the minimal variance in performance. Finally, panel (d) decomposes the training into three distinct phases, highlighting the progressive refinement where the physics constraint loss () was systematically reduced from 0.38 to 0.24, demonstrating a training progression analogous to curriculum learning strategies in deep learning [22]. Together, these panels provide a transparent and multi-faceted validation of the training protocol, ensuring the model's predictions are both accurate and physically credible.

The quantitative results of this optimized training protocol are consolidated in Table 7 and Table 8. Table 7 provides definitive evidence for the hyperparameter selection, confirming that λ = 0.1 delivered the optimal balance by achieving a validation loss of 0.136 and a 61% reduction in physical violations. The performance metrics in Table 8 affirm the model's quality, demonstrating excellent convergence ( = 0.136), the absence of overfitting ( gap of 0.015), and strong physical consistency ( = 0.24). These numerical outcomes validate the training strategy and confirm the model's strength for reliable application.

This work presents a probabilistic deep-sequence forecasting framework that integrates hybrid ANN–LSTM modeling with dual-source uncertainty quantification to improve early detection of embankment dam instability. Applied to the Megech Dam failure sequence, the model reproduced the complete deformation trajectory—from stable behavior to accelerated movement and collapse—while providing calibrated uncertainty bounds that support risk-aware interpretation. Decomposing epistemic and aleatory components revealed that epistemic uncertainty increased sharply during the acceleration phase, indicating a departure from the system’s learned behavior and aligning closely with physical failure processes.

A central outcome of this study is the introduction of the instability metric γ, which yields a quantifiable and operational early-warning signal. The γ > 0.08 threshold corresponded precisely with observed crack initiation in October 2021 and provided an actionable lead time of roughly 10–12 weeks before failure. This predictive capability contrasts strongly with deterministic models, which capture deformation trends but cannot express confidence or identify impending instability.

The probabilistic framework further demonstrated strong performance through near-ideal PICP, narrow MPIW during stable periods, and low CRPS values.

In general, the findings show that uncertainty-aware deep learning offers a strong and interpretable approach for real-time geotechnical monitoring. By linking uncertainty evolution with physical degradation processes, the framework enables more reliable early-warning decisions and improves understanding of failure progression. Future developments will extend the method to multi-hazard conditions, real-time sensor integration, and broader classes of geotechnical systems, supporting scalable and data-driven infrastructure risk management.