A Highly Efficient Depth-Averaged Material Point Method for Geophysical Mass Flows

• Kammholz is credited with the original draft, software, formal analysis, methodology, investigation, and visualization.Source: CRediT authorship contribution statement · p. 1
• Zhao, Blatny, and Pauze are credited with software and methodology; Zhao and Blatny also contributed to validation and supervision, respectively.Source: CRediT authorship contribution statement · p. 1
• Gaume is credited with review and editing, supervision, conceptualization, and funding acquisition.Source: CRediT authorship contribution statement · p. 1

1. Introduction

• Reliable probabilistic hazard mapping under climate uncertainty requires computationally efficient models for ensembles spanning wide parameter ranges.Source: §1 · p. 1
• Depth-averaged mass-flow modeling progressed from bed-aligned Savage–Hutter equations to 1D operational tools and diverse 2D Eulerian and Lagrangian models.The introduction identifies examples including RAMMS and AvaFrame com1DFA.Source: §1 · p. 2
• Eulerian models can diffuse flow margins and require shock-handling or stopping techniques, while SPH models need stabilization and neighbor searches; depth-resolved models are costly for catchment-scale ensembles.Several established models also remain commercial or closed-source.Source: §1 · p. 2
• MPM represents the material with moving particles and computes forces on a stationary grid, combining margin tracking without numerical diffusion with fixed stencils that avoid neighbor searches.The introduction notes that blended or higher-order transfers address dissipation or noise, and B-splines address spurious accelerations and cell-crossing errors.Source: §1 · p. 2
• One to two orders of magnitude faster than established operational 2D depth-averaged models, the open-source framework combines global Cartesian coordinates, a separate friction step, and Coulomb, Voellmy, and μ(I) laws.The paper describes the resulting approach as a stable, robust, and flexible solution for geophysical mass-flow simulation and probabilistic hazard assessment.Source: §1 · pp. 2–3

2.1.1. Depth-Averaged Conservation Equations

• The model assumes incompressible, shallow flow with no entrainment, bed-parallel uniform velocity, negligible bed-normal acceleration, and isotropic internal stress.It solves for flow height and horizontal velocity in global Cartesian coordinates, with vertical velocity recovered from the bed-parallel-flow condition.Source: §2.1.1 · p. 3

2.1.2. Stress Closure

• Isotropic stress closure sets the earth-pressure coefficient to K = 1 and assumes hydrostatic normal stress before depth-averaging.The resulting stress tensor is isotropic; the paper describes this as treating the granular mass like a fluid without internal friction.Source: §2.1.2 · p. 4

2.1.3. Basal Friction

• Coulomb friction uses a basal friction coefficient, while Voellmy adds a velocity-dependent turbulent term governed by ξ.The μ(I)-rheology makes friction rate-dependent; μ1 and μ2 are its quasi-static and dynamic limits, and I0 controls their transition.Source: §2.1.3 · pp. 4–5
• A separate friction step reduces trial-velocity magnitude without changing its direction, preventing nonphysical reversal and allowing flow to stop.Source: §2.1.3 · p. 5, §2.2 · p. 7

2.2. Depth-Averaged Material Point Method

• Depth-averaged MPM represents each particle as a vertical granular column; a regular Cartesian grid follows the DEM resolution, and particle positions, velocities, heights, areas, and densities carry the flow state.Release areas can come from polygon shapefiles or GIS rasters, and outputs include GIS rasters and PLY point clouds.Source: §2.2 · p. 5
• Particle sampling may be regular or Poisson-disk; the authors recommend Poisson-disk sampling to avoid the directional bias of regular spacing.They typically sample four particles per cell; quadratic B-splines are used for particle-to-grid transfer.Source: §2.2 · pp. 5–6
• Each time step is CFL-constrained, transfers particle mass and momentum to the grid, updates grid velocities, applies friction, and transfers updated velocities back to particles using blended PIC/FLIP.The authors typically choose α = 0.95 for the PIC/FLIP blend to suppress high-frequency noise while preserving energy.Source: §2.2 · pp. 5–7

3.1.1. Constant Slope

• Sliding-block simulations agree with the reference solution, and the slope-angle sweep reproduces the transition from no flow to flow and the analytical steady-state behavior.Coulomb friction produces indefinite acceleration, whereas the velocity-dependent laws approach their analytical steady-state velocities.Source: §3.1.1 · p. 8, Figure 2 · p. 8, Figure 3 · p. 9

3.1.2. Dam-Break Problem

• For the flat dam break, flow height and velocity closely match the analytical solution in the subcritical regime, with slight height underprediction in the supercritical regime; the inclined case shows similar accuracy.Grid refinement reduces relative L2 error until deviations localized at the wave front dominate and remain nearly unchanged under further refinement.Source: §3.1.2 · pp. 9–10, Figure 4 · p. 10, Figure 5 · p. 10

3.2. Validation

• The simulated granular-column deposit and sloping free surface agree well with the experimental profile, while the circular top view indicates grid-orientation independence.The calibrated Coulomb coefficient was μ = 0.61; the paper attributes its high value to the isotropic stress closure and notes that the initial collapse violates the shallowness assumption.Source: §3.2 · p. 11, Figure 7 · p. 12

3.3.1. Back-calculation

• All three friction laws produce very similar maximum flow heights and affected areas in the Brienz/Brinzauls back-calculation, while Coulomb gives slightly higher velocities along the flow path.Voellmy friction yields near-zero velocities at flow margins because its turbulent term grows as flow height decreases and rapidly halts shallow flow.Source: §3.3.1 · pp. 11–12, Figure 9 · p. 14

3.3.2. Probability Map

• The 1632-simulation Voellmy ensemble at 5 m resolution completed in 76 s, and inundation probability fell downslope from 100% on the main slope to 12% at the church.The intermediate reported probabilities were about 80% at the first road and 40% at the second road.Source: §3.3.2 · pp. 13–14, Figure 10 · p. 15

Verification and validation

• Rigid-block tests reproduce transient and steady motion, dam-break tests reproduce rarefaction-head and wet–dry-front speeds, and the granular-column deposit matches experiment without grid-orientation bias.At steep slopes, explicit rate-dependent friction slightly underestimates terminal speed; the reported relative error is about 0.5%.Source: §4 · pp. 15–16

Friction-law assumptions

• The μ(I) implementation combines basal resistance derived for a Bagnold velocity profile with momentum transport assuming a uniform, plug-flow profile; field-scale friction parameters are effective calibration values.The Bagnold and uniform profiles have momentum shape factors of 5/4 and 1, respectively.Source: §4 · p. 16

Method trade-offs

• The MPM tracks moving margins without Eulerian wet–dry cutoffs, but particle–grid transfers introduce interpolation error and PIC adds numerical dissipation; unlike SPH, it needs no neighbor searches.The authors note that MPM height updates may have inferior local conservation properties compared with Eulerian finite-volume methods.Source: §4 · p. 16

Computational performance

• One to two orders of magnitude faster than RAMMS and AvaFrame com1DFA in the tested configurations, the method also completes about 1,000 typical simulations in one minute.The authors caution that speed depends on configuration, resolution, and hardware, and that the comparisons are not a comprehensive benchmark.Source: §4 · pp. 15–16

Brienz/Brinzauls application

• The Brienz/Brinzauls simulation used about 190,000 particles at 1 m resolution; its deposit was slightly underpredicted at the toe and too wide, and it reached terrain protected by forest in the observed event.The simulation and Aaron et al. (2026) both support a single-stage failure within about 2 min, but differ in release volume and parameter calibration.Source: §4 · pp. 16–17

Probabilistic assessment

• The probability map varies only Voellmy parameters, while the ensemble finishes within 80 s and could be extended to release-volume and release-area uncertainties.The authors also identify event-calibration weights as a possible extension.Source: §4 · p. 17

Limitations and future work

• The current model neglects pore pressure, so its effect on debris-flow mobility can only be represented through effective friction parameters.Further work also identifies curvature, entrainment, alternative stress closures, and additional field-data validation as areas to address.Source: §4 · p. 17

5. Conclusion

• The method reproduced sliding-block and dam-break analytical solutions and matched a laboratory granular-column deposit without grid-orientation bias.Its separate friction step brought flow to rest without artificial stopping criteria.Source: §5 · p. 17
• At 5 m resolution, throughput exceeded 20 simulations per second, one to two orders of magnitude faster than AvaFrame com1DFA and RAMMS.The 2023 Brienz/Brinzauls back-calculation at 1 m resolution also ran faster than real time on a single CPU core.Source: §5 · p. 17
• After calibration to observed runout, Coulomb, Voellmy, and μ(I) friction yielded very similar results, though Coulomb produced a slightly wider affected area.Runout observations alone are insufficient to discriminate between the friction laws.Source: §5 · p. 17
• The model is restricted to dry, single-phase flows; physically sound debris-flow simulations require a two-phase formulation with pore-pressure feedback.Future work also identifies curvature, entrainment, anisotropic stress, and flow-dependent momentum shape factors as simplifications to relax.Source: §5 · p. 17

Funding

• Yidong Zhao and Lars Blatny acknowledge Swiss National Science Foundation support through grants CRSK-2_237763 and P500PT_230265, respectively.Source: Funding · p. 18