On the interaction of dilatancy and friction in the behavior of fluid-saturated sheared granular materials: a coupled Computational Fluid Dynamics--Discrete Element Method study

TL;DR

The paper tackles how dilatancy and friction interact with pore-fluid pressures to control failure in fluid-saturated granular media under subaerial and subaqueous conditions. It uses a 3D CFD–DEM approach (MFIX–DEM) to simulate dense and loose packings, coarse-grains microscale data to RVEs, and defines apparent friction and normalized pore pressure to map rheology across inertial and viscous regimes via the combined $K$-framework with $K=I_v+\alpha I_n^2$. An analytical poromechanics solution is used to benchmark steady-state excess pore pressure during breaching, revealing good agreement at steady state and transient deviations during onset due to frictional transients. The results advance physics-based hazard modeling by linking poromechanical feedbacks to failure modes, informing predictions of landslides and submarine slope stability through a scalable, mechanistic framework that integrates microscale grain dynamics with mesoscale rheology.

Abstract

Frictional instabilities in fluid saturated granular materials underlie natural hazards, including submarine landslides and earthquake initiation. Experiments show distinct failure behaviors under subaerial and subaqueous conditions due to coupled deformation, interparticle friction, and particle fluid interactions. We use three-dimensional coupled computational fluid dynamics, discrete element method (CFD - DEM) to investigate collapse and runout of dense and loose granular assemblies in both environments. Parametric analyses show that pore pressure evolution controls failure mode in saturated settings (fast vs slow sliding), consistent with prior laboratory experiments and lattice Boltzmann discrete element simulations: dense assemblies stabilize via dilation, whereas loose assemblies compact rapidly and transiently fluidize. At mesoscale, we coarse grain particle contact statistics and Eulerian fluid fields to define apparent friction and normalized pore pressure, and organize inertial and viscous responses using log10(In/Iv). Spatiotemporal analyses of these coarse grained fields reveal strain rate dependent behavior governed by evolving porosity and effective stress. In both environments, friction in failure shear zone is rate-strengthening with respect to inertial number (In, for dry) and viscous number (Iv, for fluid-saturated). We further utilize mesoscale stress framework to compare evolution of pore pressure in CFD - DEM of subaqueous slope collapse with an analytical solution for development of failure front, using inputs derived from numerical triaxial DEM tests on same assemblies. The analytical model reproduces steady-state excess pore pressures and captures fluid-particle coupling, but mismatch near failure onset suggests transient frictional effects. These results support physics-based hazard models and improve mechanistic understanding of saturated granular failure.

Figures (15)

• Figure 1: Model setup for subaerial and subaqueous conditions: (a) grain configuration, (b) grain size distribution, (c) first step of the model setup, where grains are poured into the simulation box and allowed to settle (initial height, $H_i = 50d_p$ or 0.05 m, and initial length, $L_i = 62.5d_p$ or 0.06 m), with aspect ratio, A = $H_i$/$L_i$ = 0.8, (d) second step of the model setup, where failure is induced by removing the supporting wall, (e) initial state of the grains in subaqueous conditions before failure (length, $L_x = 261.25d_p$ or 0.23 m; height, $L_y = 71d_p$ or 0.071 m; and width, $L_z = 25d_p$ or 0.025 m of the computational domain). Non-slip boundary conditions with frictional walls having the same coefficient of friction (0.5) as the grains, and (f) final state of the grains in subaqueous conditions after failure (final height, $H_f$, and final length, $L_f$, of the deposit).
• Figure 2: Comparison of failure and runout dynamics using velocity for subaqueous and subaerial models with dense ($\phi$ = 0.579) and loose ($\phi$ = 0.444) configurations at two time steps: 0.04 s and 0.4 s. Subaqueous model exhibit slower failure and reduced runout distances compared to subaerial model, due to particle–fluid interactions. Panels (a) and (b) show the dense subaqueous model with a minimal runout of 0.054 m (61% shorter than the maximum), while panels (c) and (d) represent the loose subaqueous model, which displays greater displacement and a more extended runout of 0.062 m (59% shorter). Panels (e) and (f) display the dense the subaerial model, characterized by abrupt failure with a runout of 0.15 m (maximum), and panels (g) and (h) illustrate the loose subaerial model, which also exhibits a runout of 0.15 m (maximum).
• Figure 3: (a) Front-position evolution for a dense subaqueous column collapse. Orange line: CFD–DEM ($\phi=0.579$, $A=0.80$, 3D, laminar). Squares: Rondon et al. (2011) flume data ($\phi=0.60$, $A\approx0.67$, quasi-2D). Both show rapid advance followed by a plateau; offsets reflect geometry and sidewall effects. (b) Normalized runout, $L/L_i$ for dense subaerial and subaqueous MFIX–DEM runs compared with Rondon data. Our subaqueous dense result ($L_f/L_i=2.92$ at $\phi=0.579$ and $3.21$ at $\phi=0.444$) falls within the experimental envelope ($L_f/L_i=2.05$–3.08 for $\phi=0.55$–0.60). Subaerial values ($L_f/L_i=3.69$ at $\phi=0.579$ and $4.01$ at $\phi=0.444$) are much larger than quasi-2D flumes ($L_f/L_i=2.6$ at $\phi=0.60$), reflecting reduced sidewall dissipation in 3D setups Jop2005. The loose case ($\phi=0.444$) lies below experimental packing ranges and densifies rapidly, preventing pore-pressure-driven mobility. Importantly, the relative ordering (subaerial $>$ subaqueous) remains consistent with experiments Rondon2011Polania2024.
• Figure 4: Depiction of shear zone development and pore pressure evolution at two time steps (0.1 s and 0.8 s) for dense ($\phi$=0.579) and loose ($\phi$=0.444) subaqueous models. Panels (a) and (b) show shear zones, defined by zero volumetric strain within the critical state framework Wood1991, for the dense configuration at the two respective time steps, while panels (e) and (f) show the corresponding results for the loose configuration. Panels (c) and (d) illustrate pore pressure changes in the dense configuration, which begins with lower initial values that increase over time; in contrast, panels (g) and (h) show the loose configuration, which starts with higher initial pore pressures that decrease with continued deformation.
• Figure 5: (a) and (b) show the dense ($\phi$ = 0.579) and loose ($\phi$ = 0.444) subaqueous models, respectively, with points A and B located within the shear zone of the dense configuration, and points C and D within the shear zone of the loose configuration. Panels (c) and (d) illustrate the dense and loose subaerial models, respectively, with points E and F located within the shear zone of the dense configuration, and points G and H within the shear zone of the loose configuration.