| 16 Oct 2026
faDebrisFOAM: debris-flow modeling with coupled terrain evolution in OpenFOAM
Abstract. Debris-flow numerical models are widely used for hazard assessment and engineering design, yet many models still treat bed evolution in a simplified or uncoupled manner. This study presents faDebrisFOAM, an open-source depth-averaged solver developed in OpenFOAM® that couples debris-flow motion, entrainment, deposition, and terrain evolution within a finite-area framework. The model represents the mixture as an equivalent fluid of variable density and incorporates mass exchanges with an erodible bed layer through entrainment and deposition laws. The Terrain Modification Algorithm (TMA), developed in this work, updates the computational grid at each time step, allowing topographic changes to feed back directly into the governing equations and simulated flow dynamics. The model is evaluated through three applications: validation against large-scale US Geological Survey (USGS) flume experiments, testing in an idealized parabolic geometry, and application to the 2015 debris-flow event in the La Mesilla catchment in northern Chile. The solver reproduces the measured evolution of flow thickness in the USGS flume experiments with good agreement. The idealized application demonstrates that terrain evolution is not a passive output but an active control on run-out, velocity, flow depth, and volume redistribution. The La Mesilla case further shows that coupled terrain modification promotes surge self-generation, self-channelization, and channel-plug formation, which are absent or weakened when the TMA is disabled. These results indicate that coupled terrain evolution can modify flow confinement, stress distribution, and preferential routing throughout the simulation. faDebrisFOAM therefore provides an open-source framework for investigating morphodynamic feedbacks in debris-flow propagation.
This is an interesting work that addresses an outstanding question in debris-flow modelling: how bed evolution feeds back into flow dynamics during an event and, of particular interest to me, how to incorporate this into a numerical model. The implementation of terrain evolution within a finite-area framework on curved, unstructured surfaces is, to my knowledge, novel. The manuscript is generally well written, and the applications illustrate the intended processes well. However, I have concerns about the numerical consistency of the moving-surface formulation, which I detail below.
Major comments
Time derivative on an evolving surface. The time-varying cell area S_B appears in all transient terms (Eqs. 1–3 and 5) and has to be incorporated correctly into the time-discretization scheme. This is not a trivial modification of the OpenFOAM finite-area solver, and the authors should discuss it in detail. In particular:
Is the time derivative evaluated with the old-time areas (S_B^n, S_B^{n−1}) as in a moving-mesh formulation, and how is this combined with the second-order scheme (l. 246) and variable time steps? This might be build into the time derivative scheme but that has to be checked and confirmed with the OpenFOAM source code.
The geometry is updated only after the flow equations have converged (l. 255). This means S_B, n_b, g_s and g_n are lagged by one time step. It is therefore unclear which geometry enters the time derivatives at each step.
Since the mesh deformation moves vertices vertically on an inclined surface, the area swept by the edges has a component tangential to the surface. In a moving-mesh finite-area/finite-volume method, this normally requires a mesh-motion flux in the convective terms to satisfy the space (geometric) conservation law. This is e.g. demonstrated in Tukovic and Jasak (2012) for the FVM (not FAM). Is this accounted for?
After the update of n_b, u is no longer tangential to the new surface, and h (measured normal to the surface) refers to a different direction. How are u and h mapped onto the new geometry, and are mass and momentum preserved in this step?
If approximations were made to avoid a fully consistent implementation, this is probably acceptable given the uncertainties involved in debris-flow modelling. However, it is important that they are stated clearly and maybe mentioned as future work to improve on.
Minor comments
I am a bit surprised not to see Rauter and Köhler (2020, “Constraints on entrainment and deposition models in avalanche simulations from high-resolution radar data”) cited. It seems the closest work to the authors’ approach that I know of (finite area combined with entrainment and deposition).
Some videos would be great.
Code availability: The code is stated to be available “upon reasonable request” (l. 466). I am not sure this complies with the Copernicus/NHESS guidelines. Beyond the manuscript, I think the authors’ contribution is very valuable. It would be a great opportunity to contribute the code to the OpenFOAM avalanche repository (https://develop.openfoam.com/Community/avalanche), where solvers, tools, and friction/entrainment models can be shared and combined.
l. 115, Eq. 4: p_b is a scalar, but the right-hand side is a vector (n_b times a scalar).
l. 121: ∇ is described as “the differential operator”, but the finite-area method uses a specific surface operator, not the ordinary one.
l. 124–125: g_s = (I − n_b n_b)·g and g_n = (n_b·n_b)g use inconsistent operators, and n_b·n_b = 1. Since g_n is used as a scalar in Eq. 4, it should presumably be g_n = n_b·g.
Sect. 2.4, “Solver internal performance”: The title does not quite fit the contents.
Fig. 5: There may be better ways to visualize the geometry, e.g. with elevation contour lines in Fig. 7.
Fig. 7: The scenario description given in Fig. 11 is very helpful and would be helpful in Fig. 7 as well.
l. 321: “Volume balance error remains below 0.02 %” is a good sign but deserves more explanation. What was measured, and where and when? Does it refer to the flow, the bed, or both, and is it per time step or cumulative?