Hydro-mechanical Model for Slope Stability Assessment: A polygonal stabilization-free discretization

Abstract: Rainfall-induced landslides are governed by the interaction between subsurface water flow and soil mechanics, requiring robust numerical methods for the simulation of variably saturated porous media. In this work, we consider a semi-coupled hydro-mechanical model based on Richards' equation and linear elasticity and propose a numerical framework based on a stabilization-free Virtual Element Method for its spatial discretization. The proposed approach naturally accommodates general polygonal meshes while avoiding problem-dependent stabilization terms, whose design may become challenging when heterogeneous and strongly non-linear coefficients are involved. The approach is combined with a mass-lumping technique to improve stability in the treatment of the storage term and with Nitsche's method to weakly impose seepage-face and infiltration boundary conditions, allowing for the automatic switching between Neumann and Dirichlet conditions. Time integration is performed using the backward Euler scheme, while non-linearities are handled through a Picard iteration. Numerical experiments demonstrate the stability and robustness of the proposed methodology and show its effectiveness in simulating rainfall infiltration and evaluating slope stability through the Local Factor of Safety.

• The paper introduces a multiphysics model that integrates Richards’ equation with linear elasticity to compute a local factor of safety without preset failure surfaces.
• It employs a stabilization-free Virtual Element Method on polygonal meshes, achieving optimal convergence and robust handling of boundary transitions via Nitsche’s method.
• Simulations under varied rainfall conditions validate the method’s accuracy in predicting infiltration, saturation, and potential slope instabilities.

Hydro-Mechanical Slope Stability Modeling Using Stabilization-Free VEM

Model Formulation and Theoretical Framework

The paper formulates a semi-coupled multiphysics model for rainfall-induced slope instability, integrating variably saturated flow governed by Richards’ equation and a linearly elastic mechanical response. The definition and evaluation of the Local Factor of Safety (LFS)—grounded in effective stress and the Mohr-Coulomb criterion—are central, enabling pointwise assessment of soil stability without prior failure surface assumptions. The pressure-head boundary transitions between infiltration and seepage are addressed via a generalized boundary enforcement mechanism based on Nitsche’s method, enabling automatic switching between Neumann and Dirichlet conditions driven by local hydraulic state.

Figure 1: The Mohr circle and Mohr-Coulomb envelope underpin the geometrical computation of LFS, which is key to distributed slope failure diagnosis.

The effective stress tensor accounts for partial saturation via Bishop’s parameter and pressure-head-dependent suction. The modeling of soil-water retention and hydraulic conductivity employs both Mualem–van Genuchten and Brooks-Corey parameterizations, supporting simulation of complex, multi-layered media.

Numerical Discretization: Stabilization-Free Virtual Element Method

The principal numerical innovation lies in the use of a stabilization-free Virtual Element Method (SFVEM) for spatial discretization. SFVEM naturally accommodates general polygonal meshes—critical for irregular, layered geomorphic domains—by constructing element-wise polynomial projections without introducing mesh- or problem-dependent stabilization terms required for classical VEM. This circumvents the instability and tuning challenges encountered with strongly nonlinear or heterogeneous coefficients, e.g., in nonlinear Richards’ flow.

The mass-lumping technique enhances stability for the storage term in the Richards equation, suppressing spurious oscillations in evolving infiltration fronts. Time integration leverages the backward Euler scheme, and nonlinearities are handled by a Picard iterative procedure with robust convergence properties.

Nitsche’s method is implemented for weak boundary condition enforcement under potentially time-varying Dirichlet/Neumann regimes, yielding automatic and smooth transitions as pressure and hydraulic flux evolve during infiltration and ponding scenarios.

Numerical Experiments and Performance Evaluation

Extensive benchmarking includes single-physics and semi-coupled tests:

• Test 1 (Richards’ Equation): Convergence tests on both non-uniform quadrilateral and Voronoi meshes exhibit optimal rates, with L2 and H1 errors decaying polynomially in mesh size. Picard iterations demonstrate mesh-independent robustness, requiring on average 8 iterations per timestep.

  Figure 2: Discrete solution for pressure-head on a fine Voronoi mesh, demonstrating spatial accuracy in complex geometric domains.

• Test 2 (Linear Elasticity): The method achieves expected linear convergence in H1 and quadratic in L2, confirming approximation quality for the mechanical subsystem.
• Test 3 (Infiltration): Systematic analysis with rain-to-saturated-permeability ratios below, equal to, and above unity shows the automatic boundary switching via Nitsche’s method. For high rainfall intensity, rapid transition to saturated boundary conditions and the enforcement of pressure control are observed, reflecting physical infiltration dynamics.

  Figure 3: Numerical solution evolution along the vertical profile during sub-saturating rainfall, exhibiting Neumann boundary persistence.

  

  Figure 4: Boundary switching is triggered as surface saturation is achieved, with Picard residual spikes marking the transition.

  

  Figure 5: Immediate saturation under intense rainfall, with Dirichlet conditions dominating, rapidly converging to equilibrium.

• Test 4 (Hydro-Mechanical Stability): Simulation on a two-layered slope domain under prolonged and intense rainfall demonstrates infiltration, saturation propagation, and LFS reduction. The method captures the emergence and spatial expansion of potential failure zones, aligning closely with established multiphysics finite element solvers.

  Figure 6: Domain geometry and investigation zone for slope stability simulation.

  

  Figure 7: Detailed mesh and boundary condition layout for layered slope stability assessment.

  

Figure 8: Water content profile post-slow rainfall, showing gradual infiltration and redistribution.

Figure 9: LFS map post-slow rainfall, with no instability detected.

Figures tracking water content and LFS over time reveal:

• Initial hydrological equilibrium is achieved after slow rainfall; intense rainfall subsequently drives surface saturation.
• LFS near unity marks initial stability; saturation leads to $\rm{LFS}\ll1$ in the upper slope and toe, diagnosing unstable regions.
• The numerical results converge with reference solutions (e.g., COMSOL), validating the SFVEM framework.

Practical and Theoretical Implications

The SFVEM scheme obviates the need for ad hoc stabilization design, directly addressing the practical challenge of mesh-generation and stability in complex geomorphological domains. The approach enables high-fidelity, distributed slope stability analysis in multi-layered, non-convex or irregular geometries. The robust weak imposition of evolving boundary conditions supports realistic modeling of rain-driven infiltration, ponding, and seepage in hydromechanical coupling.

Theoretically, the stability and coercivity of the discretized forms are established independent of mesh size or geometry, provided polynomial moment conditions are met. This generalizes VEM applicability to strongly nonlinear multiphysics scenarios, setting a foundation for future extensions to coupled flow-plasticity, fractured media, and stochastic parameterizations.

Future Directions

Potential avenues for advancement include:

• Extension to 3D polygonal or polyhedral meshes, supporting field-scale applications.
• Incorporation of nonlinear plasticity and unsaturated shear strength models.
• Integration with sensor/model fusion frameworks for real-time landslide risk assessment.
• Parallelization and scalable solvers for large-scale hydro-mechanical simulations.

Conclusion

The study demonstrates an efficient, stable, and mesh-agnostic numerical framework for hydro-mechanical slope stability modeling via stabilization-free VEM, validated across benchmarks and realistic, semi-coupled rainfall-driven landslide scenarios. The approach offers both theoretical and practical advantages in multiphysics geomechanics, establishing a versatile tool for advanced slope stability assessment and future research in coupled unsaturated soil mechanics (2607.04778).