the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# A control volume finite-element model for predicting the morphology of cohesive-frictional debris flow deposits

### Ying-Chen Wu

### Chi-Yao Hung

### Hervé Capart

### Vaughan R. Voller

To predict the morphology of debris flow deposits, a control volume finite-element model (CVFEM) is proposed, balancing material fluxes over irregular control volumes. Locally, the magnitude of these fluxes is taken proportional to the difference between the surface slope and a critical slope, dependent on the thickness of the flow layer. For the critical slope, a Mohr–Coulomb (cohesive-frictional) constitutive relation is assumed, combining a yield stress with a friction angle. To verify the proposed framework, the CVFEM numerical algorithm is first applied to idealized geometries, for which analytical solutions are available. The Mohr–Coulomb constitutive relation is then checked against debris flow deposit profiles measured in the field. Finally, CVFEM simulations are compared with laboratory experiments for various complex geometries, including canyon–plain and canyon–valley transitions. The results demonstrate the capability of the proposed model and clarify the influence of friction angle and yield stress on deposit morphology. Features shared by the field, laboratory, and simulation results include the formation of steep snouts along lobe margins.

When they transition from steep gullies to milder topography, debris flows typically spread out and slow down to form fresh deposits. By burying houses, bridges, or other assets, these may cause considerable damage to communities and infrastructure (Liu and Huang, 2006; Scheidl et al., 2008; Tai et al., 2019). This is illustrated in Fig. 1 for a case in Taiwan (courtesy of the Chi Po-lin Foundation, 2009), where debris flow deposition near a gully mouth buried the lower stories of multiple buildings. To mitigate debris flow hazards, it is therefore important to anticipate the possible extent and thickness of their deposits.

To simulate the flow and deposition of debris flows, many highly resolved models have been proposed. These typically apply mass and momentum balance equations to flows over non-erodible (O'Brien et al., 1993; O'Brien, 2006; Liu and Huang, 2006; Murillo and García-Navarro, 2012; Pudasaini, 2012; Kowalski and McElwaine, 2013; Gregoretti et al., 2016; Meng and Wang, 2016; Tai et al., 2019; Pudasaini and Fischer, 2020) or erodible substrates (Armanini et al., 2009; Bartelt et al., 2017). Such simulations, however, require detailed hydrological input data and various rheological parameters which may be difficult to obtain and may also differ dramatically from one case to another. In this context, it is worth exploring whether reduced complexity models could predict key features of debris flow deposits with less computational effort and more limited data requirements. Unlike existing models that also attempt to predict the runout at high velocities, we limit our scope and focus on predicting the final deposit morphologies of debris flows, modeled as slow, quasi-static processes.

A class of reduced complexity models developed for fluvial problems rests on defining a constitutive model for the mass flux, which in turn can be used with a mass balance equation (e.g., the Exner equation) to evolve the bed surface elevation. For applications to alluvial fans and river deltas, for instance, some models have been proposed that simply set the mass flux proportional to the current slope at that point (Voller and Paola, 2010; Lorenzo-Trueba and Voller, 2010; Lorenzo-Trueba et al., 2013). More sophisticated approaches employ the device of a critical threshold (Mitchell, 2006; Lai and Capart, 2007), whereby sediment transport occurs only when the bed inclination exceeds a critical slope (Lai and Capart, 2007; Hsu and Capart, 2008; Lai and Capart, 2009). In these models, the critical slope for the fluvial sediment flux can be derived by considering the friction stress at the sediment–water interface (the Shields stress). In some sense, this idea of a critical slope is analogous to the angle of repose governing the shapes of dry sand piles (Kuster and Gremaud, 2006; Giudice et al., 2019) or the morphology of idealized deltas and fans (Ke and Capart, 2015; Zhao et al., 2019; Chen and Capart, 2022).

Mass flux models have also been used to model mud flows. In particular, we refer to the work of Yuhi and Mei (2004), where a flux law was obtained by combining lubrication theory with a cohesive yield stress criteria. Predictions from this model were verified by comparing with analytical solutions which constrain the slope of the deposit, in axisymmetric domains, based on a cohesive yield stress criteria (Coussot et al., 1996; Yuhi and Mei, 2004). Unlike what might be seen in a sand pile or fluvial system close to the threshold, here the slope at a point varies with the thickness of the deposit.

Contrasting with fluvial and mud flows, for debris flows it is believed that both friction angle and yield stress can affect the morphology of deposits (O'Brien et al., 1993; Mangeney et al., 2010; Murillo and García-Navarro, 2012; Pudasaini, 2012; Gregoretti et al., 2016; Tai et al., 2019; Pudasaini and Fischer, 2020). The study of Coussot et al. (1996) emphasizes this point. Using only a yield stress criterion, these authors derived solutions for deposit profiles which they compared with surveyed debris flow transects. This model was found to work well for cohesive debris flow deposits with high clay content. For lower clay content, however, deposit inclinations are more consistent with control by the saturated angle of friction (Takahashi, 1991). For debris deposits mixing coarse and fine material, therefore, it appears necessary to consider both a yield stress and a saturated friction angle, as in the well-known Mohr–Coulomb model for cohesive-frictional materials.

The objective of the current work is threefold, first, we will develop a mass flux expression that considers both friction angle and yield stress in setting the critical slope under a quasi-static assumption. Secondly, we will use this mass flux in an unstructured control volume finite-element method (CVFEM) solution of the Exner mass balance equation to arrive at, for a given input mass, predictions of the final deposit location and shape. Finally, we will assess the predictive performance of this model by comparing predictions with available closed-form expressions, experimental measurements, and field observations.

In line with our objectives, we note that, in general, alluvial and debris fans build up over time in more complex ways than those immediately addressed by our proposed model and experiments. For example, channel formation, migration, and avulsion are expected to significantly affect fan evolution, especially for large-scale debris flow fans. For alluvial fan experiments devoted to these processes, the reader is referred to Le Hooke and Rohrer (1979), Whipple et al. (1998), Delorme et al. (2018), and Savi et al. (2020). Our focus here, however, is on the formation of fresh deposits, possibly over a preexisting fan surface, by unchannelized debris flows. For such conditions, illustrated by Fig. 1, we hope to formulate and verify a simplified model that could later be extended to more general conditions.

The paper is structured as follows. Section 2 presents the governing equations that form the core of our model. The CVFEM algorithm developed to obtain numerical solutions is then described in Sect. 3. Section 4 describes how we incorporate a Mohr–Coulomb constitutive relation into this framework. In Sect. 5, we explain how to supplement our CVFEM with a flux limiter to model flow over non-erodible surfaces. In Sect. 6, we check simulations against available analytical solutions. In Sect. 7, we verify our model by comparing results with field data and laboratory experiments. Finally, in Sect. 8, we discuss the contribution and limitations of our work, emphasizing how our model can help understand the influence of material properties on the morphology of debris flow deposits.

To write governing equations, we consider a debris mixture depositing over a fixed substrate of arbitrary topography. An example is shown in
Fig. 2a: supplied upstream of a steep triangular channel, the mixture flows into a trapezoidal channel of mild inclination, where it
spreads out and slows to a complete stop. We denote the time-varying surface elevation during flow using $\stackrel{\mathrm{\u0303}}{z}(x,y,t)$ and the
underlying bed topography using *z*_{b}(*x*,*y*). The corresponding profiles are shown in Fig. 2b on a schematic section.

To capture the deposition process and predict the final deposit morphology, we express mass conservation using the Exner equation (Exner, 1920, 1925)

where $\mathit{q}=({q}_{x},{q}_{y})$ is the volumetric flux (volume transferred per unit width and time), ∇⋅** q** with $\mathrm{\nabla}=(\partial /\partial x,\partial /\partial y)$ is the divergence of this flux,

*δ*is the Dirac delta function, ${\mathit{x}}_{\mathrm{s}}=({x}_{\mathrm{s}},{y}_{\mathrm{s}})$ is the location of the source, and

*Q*

_{in}is the inflow source volumetric flux. For simplicity, we assume that the flow is sufficiently slow to be regarded as quasi-static, allowing inertia effects to be neglected. At each location (

*x*,

*y*), the flux

**is assumed to be aligned with the direction of steepest descent according to**

*q*The diffusivity *ν*, however, is not assumed constant but instead depends on the local surface slope $\left|\right|\mathrm{\nabla}\stackrel{\mathrm{\u0303}}{z}\left|\right|$ according to the formula

where *ν*^{∗} is a real and positive constant, and *S*_{c}(*x*,*y*) is a critical slope dependent on material properties and on the local
instantaneous thickness of the flow layer. This dependence of *S*_{c} on the flow layer thickness is derived in Sect. 4. Combining
Eqs. (1)–(3), we see that we obtain a non-linear diffusion process with a diffusivity *ν* that depends on
the difference between the magnitude of the local gradient and the critical slope *S*_{c}. With this model, the flux is only non-zero when the
local slope $\left|\right|\mathrm{\nabla}\stackrel{\mathrm{\u0303}}{z}\left|\right|$ exceeds the critical slope *S*_{c}. By contrast, models that consider momentum effects can produce local
deposit slopes that are smaller than critical slopes (e.g., Tregaskis et al., 2022). In our model, on the deposit surface where the flow slows down to a complete stop, the flux ** q** vanishes as the local slope $\left|\right|\mathrm{\nabla}\stackrel{\mathrm{\u0303}}{z}\left|\right|$ decreases from a value that exceeds the critical
slope

*S*

_{c}to exactly the critical slope

*S*

_{c}, imposing the mathematical condition that

everywhere on the final deposit surface. Make particular note that the critical slope developed in our model (see Sect. 4) will involve the sum of two components, a constant friction slope and a yield stress term that will be an inverse function of the deposits thickness; thus, the final slope over the predicted deposited debris flow may not take a constant value.

To incorporate the above flux definition (Eqs. 2 and 3) in an Exner balance, our model includes three main components. First, we need a numerical method to solve the governing mass balance equation with the proposed flux model. Second, we need to derive an appropriate expression for the critical slope. In doing this we will consider both a friction angle and a yield stress. Third, we need to provide a limiter in our evolution algorithm to avoid fluxing out from a control volume more than the amount of material available.

To solve the Exner equation as formulated above, we adopt the control volume finite-element method (CVFEM), a method first proposed by Winslow (1966) and later extended by Baliga and Patankar (1980, 1983), Voller (2009), and Tombarevic et al. (2013). The CVFEM is a useful tool for this application because it couples the finite-element flexibility of fitting the domain geometry with the explicit mass balance of the control volume.

The application of the CVFEM to model debris flow deposits over an existing topography starts by identifying a 2-D planar problem domain (*x*,*y*) and
then covering this domain with a mesh of connected, non-overlapping, plane geometric elements. In our case, we use a rectangular domain and cover it
with an unstructured mesh of linear triangle elements (Fig. 3a). Each triangular element is associated with three vertex
node points (locally labeled *A*,*B*, and *C*) (Fig. 3b). This will result in $i=\mathrm{1},\mathrm{2},\mathrm{\dots}N$ node points in the domain,
each storing values for the fixed bed substrate elevations *z*_{b}(*x*,*y*), assumed given, and for the time-dependent flow surface elevations
$\stackrel{\mathrm{\u0303}}{z}(x,y,t)$ to be determined. To evaluate the values of *z*_{b} and $\stackrel{\mathrm{\u0303}}{z}$ at internal points in an element we use the classic
finite-element interpolation based on linear shape functions. In this way, at a point (*x*,*y*) in a given element we approximate the bed substrate
elevation as

and the flow surface elevation as

where the shape functions, *n*_{A}, *n*_{B}, and *n*_{C}, which are linear functions in *x* and *y*, take a unit value at nodes *A*, *B*, and *C*, respectively, and
vanish along the element sides opposite the labeled node, i.e, sides *B*−*C*, *C*−*A*, and *A*−*B*, respectively. Thus, the CVFEM discretization provides
piece-wise linear approximations of the bed substrate and flow surfaces. In particular, we note that in any element *j* in our domain we can readily
approximate the surface gradient by

where, ${n}_{{A}_{x}},{n}_{{A}_{y}}$, etc., are the derivatives of the shape functions. Due to the linear nature of the shape functions, we note this approximation renders a constant value for the slope in each element.

To move on, we construct an additional geometric element on our grid of triangular elements. We join the midpoint of each element side to the centroid
of each element, generating a set of connected non-overlapping control volumes around each node *i* in the domain; see
Fig. 3c and d. Thus, the control volume around node *i* has $j=\mathrm{1},\mathrm{2},\mathrm{\dots}m$ elements connected to it (the region of
support), and each of these elements contains two faces of the control volume. To discretize our governing equation, Eq. (1), we
integrate the equation over the control volume, use the divergence theorem, and make an explicit finite-difference approximation in time to arrive at
a discrete equation for the surface elevation at each node point and time step,

where *A*_{CV,i} is the area of the control volume, *Q*_{in,i} is the source flux at node *i*, and

is the net discharge out of the control volume across the two faces in element *j*, e.g., sides *S*_{AB} and *S*_{AC} in
Fig. 3b.

With an appropriate constitutive equation for determining the critical slope (see discussion below), we can use our approximations for the deposit
slope in the element, Eq. (7) to, through Eq. (2), arrive at an approximation for the flux
${\mathit{q}}_{j}=({q}_{{x}_{j}},{q}_{{y}_{j}})$ in element *j*; we should expect this value to be constant over the element. Further, if we use Δ*x* and Δ*y* to
express the change in the *x* and *y* values along a face as we move counter-clockwise around node *i* (see Fig. 3b), we can express the constant outward normal on a face with length ℓ as $\mathit{n}=(\mathrm{\Delta}y/\mathrm{\ell},-\mathrm{\Delta}x/\mathrm{\ell}$). This provides us enough
information to fully approximate the discharge in Eq. (9) in terms of the current nodal values of ${\stackrel{\mathrm{\u0303}}{z}}_{i}$ in the element (for
full details refer to Voller, 2009). On making this approximation for each element in the support of node *i* and rearranging
Eq. (8), we arrive at the following update for the surface elevation:

We note that when a node *i* is on the domain boundary (see Fig. 3d), we set the discharge across the control volume faces that coincide with the boundary to zero. Hence, Eq. (10) provides us with an explicit means of updating the nodal values of the surface elevation at time *t*+Δ*t* from the known values at time *t*. To speed up computations while still ensuring numerical stability, we use a dynamic time step
$\mathrm{\Delta}t=\mathrm{0.2}\mathrm{\Delta}{\mathrm{\ell}}^{\mathrm{2}}/max\left({\mathit{\nu}}_{\text{ele}}\right)$, where Δℓ is the average element size and *ν*_{ele} is the element diffusivity
given by ${\mathit{\nu}}_{\text{ele}}={\mathit{\nu}}^{\ast}\left(\right|\left|\mathrm{\nabla}{\stackrel{\mathrm{\u0303}}{z}}_{\text{ele}}\right||-{S}_{\mathrm{c},\mathrm{ele}})/\left|\right|\mathrm{\nabla}{\stackrel{\mathrm{\u0303}}{z}}_{\text{ele}}\left|\right|$. To let material diffuse rapidly to surrounding elements when the element slope exceeds the critical slope, we set ${\mathit{\nu}}^{\ast}=\mathrm{100}max\left({Q}_{\text{in}}\right)$.

In the previous sections, we assumed that flow occurs when the surface slope exceeds a critical slope, or, upon assuming that the direction of
steepest descent coincides with the *x* axis,

To set this critical slope, we adopt a Mohr–Coulomb failure criterion. For flow to occur, the shear stress *τ* at the base must then satisfy

where *σ* is the normal stress, *ϕ* is the saturated friction angle dependent on the solid fraction, the void fraction, and the fine content
in the fluid (Takahashi, 1991), and *τ*_{Y} is the yield stress. When the deposit surface slope is less than or equal to the critical
slope (i.e., $\left|\partial \stackrel{\mathrm{\u0303}}{z}/\partial x\right|\le {S}_{\mathrm{c}}$), the mixture remains in static equilibrium, with $\mathit{\tau}\le \mathit{\sigma}\mathrm{tan}\mathit{\varphi}+{\mathit{\tau}}_{\mathrm{Y}}$. In the limiting state, we can therefore use a force balance to derive an expression for the critical slope.

In the CVFEM model, we express this force balance element by element under the following two simplifying assumptions: (i) the surface slope in an
element is uniform (a direct consequence of our choice of linear elements), and (ii) the flow thickness in an element is also uniform. This latter
restriction is needed to keep expressions simple but will still allow us to apply the model to flows of variable thickness. Under these assumptions,
we can simply consider a two-dimensional force balance in the (*ξ*,*η*) coordinate system aligned with the surface inclination, as illustrated in
Fig. 4. Force balance in the normal and tangential directions can then be expressed as

where *ρ* is the density of the mixture, *g* the gravitational acceleration, *h* the oblique layer thickness in the *η* direction, and
*β* the bed inclination angle. To move forward, we note, by our assumptions, that

and that the vertical and oblique thicknesses are related by

Thus, on substituting Eqs. (14) and (15) into the force balance relations, Eq. (13), we obtain the following expression for the shear stress

which is an expression that matches the derivation made by Yuhi and Mei (2004). Finally, on substituting this shear stress into the Mohr–Coulomb criterion, we arrive at a model for the critical slope

The critical slope in each element can therefore be determined by setting values for the saturated friction angle and yield stress, taking into
account the local vertical layer depth $H=\stackrel{\mathrm{\u0303}}{z}-{z}_{\mathrm{b}}$. What distinguishes our expression from previous suggestions for the critical
slope Liu and Mei (1989), Coussot et al. (1996) and Yuhi and Mei (2004) is the appearance of the friction angle in addition to the yield
stress. We emphasize that this combination of the friction angle and bed thickness in the definition of the critical slope is an essential ingredient
in our model. This affects the slope of the final deposit as follows. Where the thickness of the deposit is large, for instance close to the source,
the final slope will approach the constant value tan *ϕ*. Towards the margins, by contrast, where the deposit thickness decreases, steeper slopes
and snout-like features will produced, consistent with observations.

In our CVFEM model, we assume a non-eroding bed substrate. This will require the use of a “flux limiter” to ensure mass conservation in an element over each time step of the calculation. Over a time step, we cannot flux out more material than what is available at the beginning of the time step.

With reference to the selected element in Fig. 3b, we note that one-third of the element area *A*_{ABC} contributes to the
control volume around node *A* and thus, at the start of a time step, the material available for fluxing from this sub-section of the control volume
will be $\frac{\mathrm{1}}{\mathrm{3}}({\stackrel{\mathrm{\u0303}}{z}}_{A}-{z}_{\mathrm{b}A}){A}_{ABC}$. In this way, over a time step Δ*t*, the maximum discharge that can be fluxed out from
this section, contributing to the inflows to nodes *B* and *C*, is given by

From this, following the time step calculation of the flux *Q*_{A} across faces *S*_{AB} and *S*_{AC}, we can provide a limiter by setting

where the limiting factor ≤ 1 is calculated as

Similar limiters must likewise be applied to the outflows from nodes *B* and *C*. In practice, to ensure that fluxes balance out, we apply a single
value of the limiting factor

to each element in the solution domain.

As the flow spreads and slows, it will eventually come to a complete stop and freeze in place. At each point of the resulting deposit, the limit equilibrium condition, Eq. (4), will then be satisfied. If, say because of symmetry, the surface gradient along a certain transect is aligned everywhere with this transect, then the surface profile will satisfy the following simpler equation:

with coordinate *x* taken along the transect direction. In this expression, the plus operators denote downhill deposition ($\stackrel{\mathrm{\u0303}}{z}$
and *z*_{b} decreasing in the same direction), while the minus operators denote uphill deposition ($\stackrel{\mathrm{\u0303}}{z}$ and *z*_{b} decreasing in
opposite directions). Substituting $\stackrel{\mathrm{\u0303}}{z}={z}_{\mathrm{b}}+H$, the equation becomes an ODE for the deposit thickness:

For the special case in which the bed slope $\partial {z}_{\mathrm{b}}/\partial x=\mathrm{tan}\mathit{\beta}$ is constant, Eq. (23) becomes a first-order
autonomous ODE that can be integrated analytically. In implicit form, the resulting depth profile *H*(*x*) is given by

where

In the above expressions, *H*(*x*_{0}) is the boundary condition at *x*_{0}, which can be any point within the depositing region. Note that *A* will be zero
for frictionless material deposits on a horizontal plane or frictional materials depositing downhill when the friction slope equals the bed slope,
and *B* will be zero when there is no yield stress. In what follows, these analytical solutions will be used for three purposes: to clarify model
properties, to verify the numerical method, and to calibrate material parameters when comparing model results with field and laboratory data.

As an example, analytical solutions for the centerline profiles of cohesive-frictional deposits over an inclined plane are illustrated in
Fig. 5. For each case and deposit height *H*, we supply material at a single point corresponding to the apex of each deposit. To
facilitate comparison, the source locations are adjusted to let the deposits have the same toe location. These locations *x*_{s} are determined
using the formula ${x}_{\mathrm{s}}=(AH-B\mathrm{ln}(|AH+B|\left)\right)/{A}^{\mathrm{2}}-C$, where $A=-\mathrm{tan}\mathit{\beta}+\mathrm{tan}\mathit{\varphi}$, $B={\mathit{\tau}}_{\mathrm{Y}}/\left(\mathit{\rho}g\right)$, $C=(-B\mathrm{ln}(B\left)\right)/{A}^{\mathrm{2}}$. In all cases, the material properties are the same, and the origin is taken at the downstream end of each deposit. This representation is
chosen to highlight two important features of the solutions. First, the shape of the deposit toe does not change with the size of the deposit and
depends only on the bed slope and material properties. Secondly, the different material properties affect separate features of the profiles. The yield
stress *τ*_{Y} controls the scale of the steep snouts, where the deposit thickness reaches zero, whereas the friction slope tan *ϕ* sets
the deposit inclination far away from the snouts, where the deposit thickness becomes large.

It follows from these properties that a single profile of sufficient length through the toe of a deposit is sufficient to calibrate the material
properties of the model. This is very useful as it greatly facilitates model application to field and experimental cases. A second implication is
that, for deposits of large size compared to the scale of the snouts, deposit shapes may be well approximated by surfaces of constant slope. For the
deposits of Fig. 5, setting the yield stress to zero would produce upright cones of slope tan *ϕ* centered at the apex of each
deposit. In general, however, the morphology of deposits will be affected by both the yield stress and the friction angle.

In this section, we evaluate the CVFEM numerical model by comparing results with analytical solutions. This provides an opportunity to show how model results depend on material parameters, for some additional simple cases. We also examine how mesh geometry and size affect the accuracy and performance of the model.

## 7.1 Comparison with analytical solutions

To verify our CVFEM algorithm we consider deposits formed by supplying material from a point source onto three idealized geometries: (i) a horizontal
plane, (ii) an axisymmetric conical basin of slope tan *β*=0.05, and (iii) an inclined plane of constant slope (tan *β*=0.02). The CVFEM model
for each of these cases operates in Cartesian coordinates and will produce 3D deposit shapes. Thus, to compare with analytical solution profiles we
need to select appropriate transects. For the horizontal plane and conical basin cases, we examine radial profiles (see
Fig. 6a, b, e, f, i, and j). For the inclined plane, we select two profiles through the source point: a longitudinal
profile in the direction of the base slope and a transverse profile orthogonal to this direction (see Fig. 6c, d, g,
h, k, and l). For the longitudinal profile (Fig. 6c, g, and k), we can use the analytical solution in
Eq. (24) as the exact solution. For the transverse profile (Fig. 6d, h, and l), the transect is
not a true symmetry axis. Nevertheless, the analytical solution obtained by setting tan *β*=0 can be used as an approximate solution. For each
case, we impose a fixed thickness of the deposit at the origin for both analytical solutions and numerical solutions.

To show how parameters affect results and check the numerical model under different assumptions, we compare numerical and analytical solutions for
three groups of material properties: (1) tan *ϕ*>0, *τ*_{Y}=0; (2) tan *ϕ*=0, *τ*_{Y}>0; and (3) tan *ϕ*>0,
*τ*_{Y}>0. We find excellent agreement between the computational results and the analytical solutions in each case regardless of the choice
of parameters (Fig. 6) and therefore verify the proposed CVFEM algorithm.

In the cases with constant friction stress and no yield (tan *ϕ*>0, *τ*_{Y}=0, Fig. 6a–d), the simulated
final deposits have constant surface slopes equal to the friction slope, which is consistent with physical and computational models for sand piles
(Kuster and Gremaud, 2006; Giudice et al., 2019).

In the cases with only yield stress (tan *ϕ*=0, *τ*_{Y}>0, Fig. 6e–h), we obtain piles with mild
slopes in the central regions and steep slopes along the margins of the deposit, resulting in toes that have a snout-like profile. This matches the
analytical solutions and models proposed by Coussot et al. (1996) and Yuhi and Mei (2004) for slow mud flows (fluids with a Bingham plastic
rheology). By considering the yield stress, it is therefore possible to reproduce the snout-like toes observed along the margins of many debris flow,
mud flow, and snow avalanche deposits (Johnson, 1970; Pudasaini and Hutter, 2007).

Finally, in the cases with both friction and yield stress (tan *ϕ*>0, *τ*_{Y}>0, Fig. 6i–l), we note
that snout-like profiles are again obtained at the toes. Away from the toes, however, the deposit slope now tends toward a finite inclination,
controlled by the friction angle. Overall the results in Fig. 6 clearly demonstrate how the friction angle and yield
stress affect deposit shapes.

## 7.2 Influence of mesh geometry and size

By using triangular elements as building blocks, the CVFEM model can be applied to either structured or unstructured meshes. In
Fig. 7, we show how model results are affected by mesh geometry and size. For these calculations, we again consider a
simple test case in which material supplied at the origin deposits over a horizontal substrate, under the combined influence of friction angle and
yield stress (tan *ϕ*>0, *τ*_{Y}>0). For these tests a prescribed volume of material is supplied, by controlling the accumulated
discharge supplied at the source.

Three different meshes are considered: a structured mesh, built from triangular elements laid out in a row-column pattern (Fig. 7a); an unstructured mesh, constructed by the mesh generation algorithm of Engwirda (2014) (Fig. 7b); and a fine unstructured mesh, constructed by the same algorithm (Fig. 7c). The corresponding model results are shown in Fig. 7d–f, representing the calculated topography by elevation contours.

In Fig. 7d, clear directional errors can be seen when results are computed on the structured mesh. In this case, the
deposits contours visibly protrude along the *x* and *y* directions. Such errors can be reduced by using an unstructured mesh
(Fig. 7e) and by calculating on a finer grid (Fig. 7f). By doing so, the calculated
contours become closer to the expected circular pattern.

By performing tests on progressively finer meshes, we can also check the convergence of our CVFEM algorithm. For this purpose, we consider two
predictive measures to assess grid convergence. The first is the normalized modeling error $(h-H)/h$ between the calculated deposit height *H* and the
analytical value *h*=0.241 m. Noting that even unstructured meshes can introduce some bias (in particular when the mesh is coarse), our second
measure is the difference between the maximum and minimum radii associated with the contour $\stackrel{\mathrm{\u0303}}{z}=\mathrm{0.1}\phantom{\rule{0.125em}{0ex}}h$, normalized by the analytical
value *r*_{10}=1.628 m.

In Table 1, we list these height and radius measures for different mesh sizes, as characterized by the average length of element edges and by the number of elements of the mesh. As the mesh is refined, we see that both measures converge to 0. In particular, the normalized modeling error $(h-H)/h$ clearly converges to the first order with respect to element size. In Table 1, we also report the computational time in seconds needed to run these simulations on an i5-9500 Intel processor. We emphasize that the use of the dynamic time step in our solution contributes significantly to its efficiency. Preliminary versions of the code used a constant time step selected by $\mathrm{\Delta}t=\mathrm{0.25}\mathrm{\Delta}{\mathrm{\ell}}^{\mathrm{2}}/{\mathit{\nu}}^{\ast}$. This approach produces the same predictions as those reported here but requires over an order of magnitude more CPU time.

To further test the model, in this section we present comparisons with field and laboratory data. Measured profiles for the toes of debris flow deposits are first exploited to verify the applicability of the critical slope and Mohr–Coulomb model to field cases. Comparisons with new laboratory experiments are then made to check the ability of the CVFEM model to predict the overall morphology of cohesive-frictional deposits. The calibration and CVFEM numerical model code and input and output data discussed in this section are available in Chen et al. (2022).

## 8.1 Comparison with field profiles

Coussot et al. (1996) observed six natural debris flow deposits in the French Alps. By categorizing these deposits by their fine fractions
(ratio of particles whose diameter is less than 40 µm to total solid volume), they found that debris flow deposits with a low fine
fractions (< 1 %), at Bourgeat, Le Bez, and Ste-Elisabeth, exhibit nearly straight profiles, whereas debris flow deposits with a high fine
fractions (10 %–15 %), at Les Sables, St-Julien, and Mont Guillaume, exhibit significant snout-like toes. Therefore, Coussot et al. (1996) focused on the latter case to test their model involving only the effect of yield stress. For each deposit with a high fine fractions, they documented two profiles, frontal and lateral, which they sought to fit by calibrating two parameters: the bed slope tan *β* and the yield stress over specific weight ${\mathit{\tau}}_{\mathrm{Y}}/\left(\mathit{\rho}g\right)$. For each site, they calibrated these parameters separately for the frontal and lateral profiles. The lengths of the profiled deposits were in the range 2 to 15 m, and the corresponding thicknesses were in the range 1.5 to 3 m.

Straight profiles, characterized by a constant slope, can be reproduced in our model by setting the yield stress to zero and the saturated friction
slope tan *ϕ* equal to the deposit surface slope. We therefore need to check whether our model can reproduce also the snout-like profiles observed
for the case of high fine fractions. By taking both friction angle and yield stress into account, we can test whether analytical profiles can
reproduce the field profiles using only one set of parameters per site. For this purpose, we assume that the frontal and lateral profiles at the same
site share the same material properties (tan *ϕ* and ${\mathit{\tau}}_{\mathrm{Y}}/\left(\mathit{\rho}g\right)$). For the frontal profile, we treat the substrate bed
slope (tan *β*) as unknown, while for the lateral profile we assume that the bed slope is zero (tan *β*=0).

To estimate the three parameters, we fit the analytical solution given by Eq. (24) to the two measured profiles. Assigning the
measured toe position as the boundary condition (*x*_{0}=*x*_{toe} and *H*(*x*_{0})=0), we obtain a predicted profile for given values of the frontal
substrate bed slope tan *β*, the saturated friction angle tan *ϕ*, and the ratio of yield stress over specific weight ${\mathit{\tau}}_{\mathrm{Y}}/\left(\mathit{\rho}g\right)$. Then, on minimizing the root-mean-square error (RMSE) between predicted and measured fan profiles we arrived at best-fit estimates for
tan *β*, tan *ϕ* and ${\mathit{\tau}}_{\mathrm{Y}}/\left(\mathit{\rho}g\right)$.

In Fig. 8, we compare the resulting profiles with the field data, normalizing both axes by the length scale ${\mathit{\tau}}_{\mathrm{Y}}/\left(\mathit{\rho}g\right)$. From the figure, we see that our critical slope model based on the Mohr–Coulomb constitutive law provides close fits to the field
observations in the cases of Les Sables (Fig. 8a and b) and Mont Guillaume (Fig. 8e and f) and an
acceptable fit in the case of St-Julien (Fig. 8c and d). The ability to use the same parameters (friction angle and yield
stress) to fit both frontal and lateral profiles indicates that, while it is significant during the flowing stage, inertia may only play a limited
role in determining the final deposit morphologies. For the debris flow deposits in Les Sables and St-Julien, the additional parameter (tan *ϕ*)
plays an important role in determining the deposit morphology and provides the degree of freedom needed to describe each pair of profiles for the
same site using the same set of parameters. For Mont Guillaume, calibration produces a low value for the saturated friction angle, indicating that the
yield stress and bed slope are sufficient to represent the deposit morphology. This may be due to the high clay content at this site.

Depending on scale and material composition, either the friction angle or the yield stress alone may be sufficient to characterize certain debris deposits in the field. Both influences, however, must be considered for intermediate cases and to encompass the range of possible behaviors in a single description. However, other effects could also cause or contribute to the formation of steep snouts in debris flows, this could include, for example, pore pressure loss at the front and margins (e.g., Iverson, 1997; Iverson and Vallance, 2001; Savage and Iverson, 2003; Iverson et al., 2010; Gray, 2018) and the frictional hysteresis of the angular sand particles (e.g., Félix and Thomas, 2004; Mangeney et al., 2007; Edwards et al., 2017; Rocha et al., 2019; Edwards et al., 2019). These effects are not currently included in our model.

To go beyond transect comparisons, in the next section we will use laboratory experiments to test the ability of our CVFEM model to simulate the complete morphology of cohesive-frictional deposits.

## 8.2 Experimental design and conditions

To investigate the morphology of cohesive-frictional deposits in well-controlled conditions but more complex geometries, we conducted new laboratory experiments at the Hydrotech Research Institute of National Taiwan University. As illustrated in Fig. 9, these experiments were conducted in faceted flumes assembled from beveled wood panels. Different from alluvial fan experiments (Le Hooke and Rohrer, 1979; Whipple et al., 1998; Delorme et al., 2018; Savi et al., 2020), which involve water and cohesionless sediment, here the deposits are built from mixtures of sand, kaolinite, and water, mixed together thoroughly to behave as a cohesive-frictional material. To produce varied deposits, controlled volumes of these mixtures were supplied upstream of steep V-shaped canyons and conveyed by these canyons to zones of milder topography where they could spread, slow, and freeze in place. Water-soluble dyes were added to distinguish the materials supplied to different canyons. Finally laser scanning (Ni and Capart, 2006; Lobkovsky et al., 2007) was used to acquire high-resolution maps of the substrate and deposit topography.

As illustrated by the photographs of Fig. 9d–f, the experiments generate rather idealized deposits, which nevertheless reproduce various features exhibited by debris flow deposits in the field. These include steep snouts along lobe margins and cusped weld lines where separate lobes come into contact. Surface folds, indicative of viscoplastic behavior, can be observed at various locations (see, for instance, the lobe in the foreground of Fig. 9e), similar to the folds visible in some areas of the field deposit shown in Fig. 1.

The following two series of tests were conducted: canyon–plain experiments (T01–T04), using the geometry shown in Fig. 9a, and
canyon–valley experiments (T11–T15), using the geometry shown in Fig. 9b. For the canyon–plain experiments (runs T01–T04), two V-shaped canyons connect to a wide U-shaped plain that has a planar floor and vertical walls. The canyon thalwegs have an inclination of
18.8^{∘} relative to the planar floors. The experiments were designed so that the whole flume could be tilted away from horizontal, in the
longitudinal direction of the tributary channels. In each run, a mixture of 61.4 wt % silica sand (*d*_{50}=0.6 mm),
8.8 wt % kaolinite, and 29.8 wt % water was used to deposit a fan into an initially empty and clean flume.

For run T01 the flume floor was horizontal, and two equal volumes of mixture were poured simultaneously upstream of the two canyons. For run T02 the
inclination was the same, but the volumes supplied to the two tributary channels TC1 and TC2 were in a ratio of 1 to 2. The continuous mass input was
arranged to start and end at the same time. Runs T03 and T04 were identical to run T02 apart from different flume tilt angles, set to 3
and 6^{∘}, respectively. For these runs, the topography was scanned with the laser oriented perpendicular to the canyons, and the resulting DEM data have a
resolution of 2 mm × 2 mm.

For the canyon–valley experiments (T11–T15), the flume had a more complex configuration, illustrated in Fig. 9b. Three
V-shaped canyons, having thalweg inclinations equal to 14^{∘}, connect at right angles to a wide trapezoidal channel of longitudinal inclination
equal to 3^{∘}. Two of the canyons (TC1 and TC2) connect on the right side, and one connects on the left side (TC3) slightly downstream. In all runs the
initial state of the canyon was clean wood, but the main channel was covered by a 2 cm thick layer of unconsolidated silica sand
(*d*_{50}=0.6 mm). For run T11 a controlled volume of mixture was supplied to tributary TC1 only. For run T12 different volumes were supplied
simultaneously to tributaries TC1 and TC2 and arranged to start and end at the same times. For run T13 different volumes were supplied to
tributaries TC1 and TC3, and for run T14 different volumes were supplied simultaneously to all three tributaries.

For run T15, deposits were formed in three separate stages. In the first stage, deposits were formed as in run T13 by supplying different volumes to tributaries TC1 and TC3. In the second stage, a constant water discharge was supplied to the main channel for 20 min, eroding the first-stage deposits. The resulting topography was scanned to provide initial conditions for the third stage, in which new volumes of material were supplied to tributaries TC1 and TC3. This provides an opportunity to examine the formation of fresh deposits onto a preexisting deposit surface. For all canyon–valley experiments, the topography was scanned with the laser oriented orthogonal to the main channel and parallel to the canyons, and the resulting DEM data have a resolution of 5 mm × 5 mm.

^{∗} The parameter ranges in documented field cases are collected or calculated from the data of Iverson (1997) and Zhou and Ng (2010). ${}^{\ast \ast}$ Viscosity for the experimental interstitial fluid is estimated from a set of viscometer measurements.

In Table 2, we present the range of parameter values covered by the laboratory experiments (runs T11–T15) and compare them to typical values for natural debris flows (Iverson, 1997; Zhou and Ng, 2010). From the table, we can see that the experiments exhibit smaller Froude and Reynolds numbers; hence, inertia effects are smaller in the experiments than in field cases. Nevertheless, the Bagnold number (ratio between collisional and viscous forces), the Savage number (ratio between collisional and frictional forces), and friction number (ratio between frictional and viscous forces) in the experiments are within the range of values encountered for natural debris flows.

In the next sections, the data from these different experiments will be used to calibrate model parameters and compare CVFEM simulation results with the topography measurements acquired in each case.

## 8.3 Comparison with canyon–plain experiments

To apply the CVFEM method to the canyon–plain experiments (runs T01–T04), we first determine model parameters from longitudinal deposit profiles,
measured along the centerlines of the deposits from each canyon (see example profile locations in Fig. 10a). The calibration method
used is the same as the one applied to the field profiles, except that the substrate slope tan *β* is known from the flume geometry; hence, only
the material parameters tan *ϕ* and ${\mathit{\tau}}_{\mathrm{Y}}/\left(\mathit{\rho}g\right)$ remain to be determined. For this set of experiments, some variability in material
properties was caused by uncontrolled variations in moisture in the kaolinite. For this reason, we use all eight of the available measured profiles together to estimate a pair of parameters that best fit the whole series of experiments. The resulting estimates are tan *ϕ*=0.063 and
${\mathit{\tau}}_{\mathrm{Y}}/\left(\mathit{\rho}g\right)=\mathrm{0.115}\phantom{\rule{0.125em}{0ex}}\mathrm{cm}$.

Initial and boundary conditions are set up as follows. An unstructured mesh of average element size Δℓ=4 mm is generated over the
problem domain. The flume topography measured before each experiment is then used to set the substrate and initial surface elevations
*z*_{b}(*x*,*y*) and $\stackrel{\mathrm{\u0303}}{z}(x,y,\mathrm{0})$. To input the deposits, constant discharge sources are placed at the vertices closest to the upstream ends
of the two channel thalwegs (*x*,*y*) = (0,10) cm and (*x*,*y*) = (0,34.8) cm, respectively. The rates of these discharges are set
to ensure that at the end of the chosen simulation time the volumes supplied match the measured experimental volumes for each source.

In Fig. 10, we compare simulation results with the experimental measurements for the four runs T01–T04. Qualitatively and quantitatively, the simulations are found to predict the measured topography of the deposits reasonably well. As indicated by the contours, both the simulations and experiments produce steep snouts along lobe margins, well-defined cusps along weld lines, where two lobes come into contact, and saddle points along these same weld lines.

In planform (Fig. 10a–h), the model is able to reproduce the outer boundaries of the deposits well, both along the steep canyon and valley sides and over the mildly inclined floor. This agreement holds for both the symmetric (equal volumes supplied to the two canyons) and asymmetric cases (unequal volumes). The model also reproduces the gradual elongation of the deposit lobes as the flume inclination is increased.

In profile (Fig. 10i–l), model results also compare well with the measurements. The model is able to capture the observed deposit slope variations, from steep upstream of the canyons, to mild over the thick lobes, and back to steep snouts at the downstream toes. In both the simulations and experiments, the deposits become gradually shallower as the flume slope is increased.

Nevertheless, there are some discrepancies between the CVFEM model and the experiments. Within the canyons and at canyon outlets, the model produces narrower and shallower deposits than the experimental results. This could be due to the geometrical simplifications used to derive the critical slope model, in which the basal substrate was assumed approximately parallel to the surface. There are also some mismatches in planform length and width, possibly due to the previously mentioned moisture variations between runs. This is especially notable for the distal parts of run T04.

## 8.4 Comparison with canyon–valley experiments

For the canyon–valley experiments (T11–T15), the moisture was better controlled; hence, the material composition was more nearly identical for all
runs. We can therefore use the longitudinal profile for the single deposit produced in run T11 (red line in Fig. 11a) to calibrate
the parameters for all cases. The resulting estimates for the material parameters, tan *ϕ*=0.118 and ${\mathit{\tau}}_{\mathrm{Y}}/\left(\mathit{\rho}g\right)=\mathrm{0.344}\phantom{\rule{0.125em}{0ex}}\mathrm{cm}$,
are used for all CVFEM simulations of this series.

To simulate these runs, we use an unstructured mesh of average element size Δℓ=5 mm. Like before, for each case we obtain the initial condition by sampling the measured pre-event topography at the mesh nodes. For runs T11–T14, we prescribe point sources of constant discharge at the vertices where canyon thalwegs intersect the domain boundaries (red points in Fig. 11). For run T15, the deposits partly buried the canyons; hence, line sources are used instead at cross sections along the domain boundary (red lines in Fig. 11j). The discharge for these various sources are again set to match the volumes of the individual deposits.

To compare measured and simulated results, topographic contours and deposit thickness maps for the different cases are presented in Fig. 11. Overall, good agreement is observed between the CVFEM simulations and the experiments. Because the main channel dips to the left, the deposit lobes acquire an asymmetric, distorted shape, which is well reproduced by the simulations. In both the experiments and the simulations, steep snouts are produced along the outer and side margins of the deposits, where they connect with the valley bed and sides. For runs T12 to T14, the weld lines obtained where different lobes come into contact are also accurately predicted. Using a single set of material parameters, the simulations also reproduce the deposit thickness distributions obtained in the different experiments well.

Similar to the canyon–plain experiments, some discrepancies are nevertheless observed between the simulations and experiments. The simulated fans are
slightly wider (*x* direction) and shorter (*y* direction) than their experimental counterparts. This could be due to momentum, neglected in our CVFEM
model, allowing the experimental mixture to flow out further in the canyon direction.

Finally, the T15 experiment (Fig. 11i and j) allows us to test our model for the case of fresh deposits onto a preexisting deposit of complex shape. This case is similar to the 2009 Xinfa debris flow shown in Fig. 1, where the inundation of houses suggests a flow of around 2–3 m deep occurring on a much larger (around 40 m tall) preexisting debris flow fan. We see that the experimental deposit exhibits similar features to the Xinfa debris flow deposits, in particular the well-defined snouts of the secondary deposits on top of the preexisting deposit. In the experiments, the earlier deposit may deform slightly due to the new deposition, but we neglect this complication and take it as a new rigid boundary in the simulations. For this challenging case, the CVFEM model again provides an excellent overall prediction of the thickness, extent, and morphology of the secondary deposits. In both the experiments and simulations, the fresh deposits do not completely cover the preexisting lobes. The fresh material stops over these lobes at some locations, flowing further at other locations to form new secondary lobes. The corresponding margins again feature well-defined snouts.

In this paper, we proposed a novel computational model to simulate the morphology of debris flow deposits. The numerical algorithm uses the control volume finite-element method (CVFEM) to discretely approximate fluxes over a finite-element mesh and explicitly enforce mass balance over prescribed control volumes. Unlike fluvial and mud flow deposits, debris flow deposits are affected by both cohesion and friction. To set the critical slope at which flow starts or stops, we therefore adopted a Mohr–Coulomb criterion that includes both a yield stress and a friction angle.

We verified the CVFEM algorithm by comparing computational results to analytical solutions in idealized cases, obtaining excellent agreement. Comparisons with field profiles were then performed to check that our critical slope model based on the Mohr–Coulomb relation can reproduce the key features of debris flow deposits. For deposits characterized by a high fine fraction, the inclusion of a yield stress allows our model to reproduce the blunted snouts observed at deposit toes. Accounting for a friction angle, on the other hand, allows our model to match the trailing slope observed away from the toes and makes the model also applicable to deposits with a low fine fraction, which feature more even slopes.

Finally, comparisons with new laboratory experiments were conducted to test the ability of our CVFEM model to predict the extent, thickness, and morphology of cohesive-frictional deposits in more complex geometries. The conditions considered include supply by single and multiple sources, and deposition over faceted substrates and preexisting deposits. Using material parameters calibrated from one or more transects, the model is found to reproduce the measured topography well in all cases. Deposit features captured accurately by the model include steep snouts along the margins of primary and secondary lobes and cusped weld lines where different lobes come into contact.

Although good agreement was obtained for the different comparisons, we do recognize some possible limitations. First, the model cannot simulate the dynamic evolution of debris flows, since it is only designed for computing the quasi-static morphology of debris flow deposits and relies on a hypothetical diffusivity and pseudo time steps. Second, the model neglects flow momentum and basal erosion; hence, it does not apply to rapid or erosive debris flows (Armanini et al., 2005). Besides, as noted above, the model does not include other effects that may lead to the formations of snouts and channel levees, such as pore pressure loss and frictional hysteresis. Likewise, it does not account for the thixotropic behavior whereby deposits gradually solidify to form a new substrate for fresh deposits (Murata, 1984; Roussel, 2006). Finally, our model and experiments do not include processes like channel formation, migration, and avulsion that also affect the evolution over time of debris and alluvial fans (Le Hooke and Rohrer, 1979; Whipple et al., 1998; Delorme et al., 2018; Savi et al., 2020).

Despite these current limitations, we have shown that a critical slope model accounting for yield stress and friction angle can simulate deposit morphology with excellent efficiency using dynamic time steps. Aside from computation time, another key consideration is the work involved in calibrating model parameters. In this regard, an important advantage of our proposed simple model is that its parameters can be calibrated directly from topography profile data. As done in the paper for the experimental cases, all model parameters can be acquired from a single long profile through observed deposits. It is therefore not necessary to run the three-dimensional model multiple times to adjust model parameters by trial and error. More complex models, by contrast, typically require multiple iterations or must rely on other sampling and material analysis to acquire parameter values. The model can easily calibrate parameters for a broader range of conditions than has been considered previously. To simulate such deposits in complex geometries, moreover, the control volume finite-element method (CVFEM) was found to provide a promising numerical approach and could possibly be extended in the future to more general processes or other geomorphic systems.

The MATLAB codes of the CVFEM model and parameter calibrations, input data (experimental pre-event and post-event topography data), and model output data are available at Chen et al. (2022), https://doi.org/10.5281/zenodo.7324739.

The algorithm for constructing unstructured mesh used in this paper is an open-access MATLAB package built by Engwirda (Engwirda, 2014) available at https://github.com/dengwirda/mesh2d (Engworda, 2019).

The algorithms for linearly interpolating triangulation and plotting contours for triangular mesh used in this paper are open-access MATLAB codes built by Hanselman (Hanselman, 2021a, b) available at https://www.mathworks.com/matlabcentral/fileexchange/38925-linearly-interpolate-triangulation and https://www.mathworks.com/matlabcentral/fileexchange/38858-contour-plot-for-scattered-data.

TYKC performed numerical simulations, contributed to designing the research methodology, analyzed data, created figures, and wrote the first draft. YCW conducted laboratory experiments, contributed to designing the research methodology, and analyzed data. CYH discussed the results and contributed to the writing, visualization, and funding acquisition. HC contributed to designing the research methodology and experiment design and setup; discussed the results; advised the project; acquired funding; and contributed to the visualization, writing, and editing of the manuscript. VRV contributed to designing the research methodology; discussed the results; advised the project; provided resources; and contributed to the visualization, writing, and editing of the manuscript.

The contact author has declared that none of the authors has any competing interests.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Useful feedback regarding the implementation and applications of our proposed scheme was provided by Philippe Delandmeter (Fugro Engineers Belgium).

This research has been supported by the Ministry of Science and Technology, Taiwan (grant nos. MOST110-2636-M-005-001 and MOST108-2911-I-002-564).

This paper was edited by Tom Coulthard and reviewed by Stefan Hergarten and Chris Johnson.

Armanini, A., Capart, H., Fraccarollo, L., and Larcher, M.: Rheological stratification in experimental free-surface flows of granular–liquid mixtures, J. Fluid Mech., 532, 269–319, https://doi.org/10.1017/S0022112005004283, 2005. a

Armanini, A., Fraccarollo, L., and Rosatti, G.: Two-dimensional simulation of debris flows in erodible channels, Comput. Geosci., 35, 993–1006, https://doi.org/10.1016/j.cageo.2007.11.008, 2009. a

Baliga, B. R. and Patankar, S. V.: A new finite-element formulation for convection-diffusion problems, Numer. Heat Transfer, 3, 393–409, https://doi.org/10.1080/01495728008961767, 1980. a

Baliga, B. R. and Patankar, S. V.: A control volume finite-element method for two-dimensional fluid flow and heat transfer, Numer. Heat Transfer, 6, 245–261, https://doi.org/10.1080/01495728308963086, 1983. a

Bartelt, P., Bieler, C., Bühler, Y., Christen, M., Deubelbeiss, Y., Graf, C., McArdell, B., Salz, M., and Schneider, M.: RAMMS – rapid mass movement simulation, A numerical model for debris flows in research and practice, User Manual v1.7.0, Debris Flow, WSL Institute for Snow and Avalanche Research SLF, http://ramms.slf.ch/ramms/downloads/RAMMS_DBF_Manual.pdf (last access: 19 October 2022), 2017. a

Chen, T.-Y. K. and Capart, H.: Computational morphology of debris and alluvial fans on irregular terrain using the visibility polygon, Comput. Geosci., 169, 105228, https://doi.org/10.1016/j.cageo.2022.105228, 2022. a

Chen, T.-Y. K., Wu, Y.-C., Hung, C.-Y., Capart, H., and Voller, V. R.: Model codes and data for “A control volume finite element model for predicting the morphology of cohesive-frictional debris flow deposits”, Zenodo [data set], https://doi.org/10.5281/zenodo.7324739, 2022. a, b

Coussot, P., Proust, S., and Ancey, C.: Rheological interpretation of deposits of yield stress fluids, J. Non-Newton. Fluid, 66, 55–70, https://doi.org/10.1016/0377-0257(96)01474-7, 1996. a, b, c, d, e, f, g

Delorme, P., Devauchelle, O., Barrier, L., and Métivier, F.: Growth and shape of a laboratory alluvial fan, Phys. Rev. E, 98, 012907, doi10.1103/PhysRevE.98.012907, 2018. a, b, c

Edwards, A. N., Viroulet, S., Kokelaar, B. P., and Gray, J. M. N. T.: Formation of levees, troughs and elevated channels by avalanches on erodible slopes, J. Fluid Mech., 823, 278–315, https://doi.org/10.1017/jfm.2017.309, 2017. a

Edwards, A. N., Russell, A. S., Johnson, C. G., and Gray, J. M. N. T.: Frictional hysteresis and particle deposition in granular free-surface flows, J. Fluid Mech., 875, 1058–1095, https://doi.org/10.1017/jfm.2019.517, 2019. a

Engwirda, D.: Locally optimal Delaunay-refinement and optimisation-based mesh generation, PhD thesis, University of Sydney, http://hdl.handle.net/2123/13148 (last access: 20 April 2015), 2014. a, b

Engworda, D.: MESH2D: Delaunay-based mesh generation in MATLAB, GitHub [code], https://github.com/dengwirda/mesh2d (last access: 27 March 2021), 2019. a

Exner, F. M.: Zur Physik der Dünen, Akad. Wiss. Wien, Math. Naturwiss. Klasse, 129, 929–952, 1920. a

Exner, F. M.: Über die Wechselwirkung zwischen Wasser und Geschiebe in Flüssen, Akad. Wiss. Wien, Math. Naturwiss. Klasse, 134, 165–204, 1925. a

Félix, G. and Thomas, N.: Relation between dry granular flow regimes and morphology of deposits: formation of levées in pyroclastic deposits, Earth Planet. Sc. Lett., 221, 197–213, https://doi.org/10.1016/S0012-821X(04)00111-6, 2004. a

Giudice, L. A., Giammanco, G., Fransos, D., and Preziosi, L.: Modeling sand slides by a mechanics-based degenerate parabolic equation, Math. Mech. Solids, 24, 2558–2575, https://doi.org/10.1177/1081286518755230, 2019. a, b

Gray, J. M. N. T.: Particle segregation in dense granular flows, Annu. Rev. Fluid Mech., 50, 407–433, https://doi.org/10.1146/annurev-fluid-122316-045201, 2018. a

Gregoretti, C., Degetto, M., and Boreggio, M.: GIS-based cell model for simulating debris flow runout on a fan, J. Hydrol., 534, 326–340, https://doi.org/10.1016/j.jhydrol.2015.12.054, 2016. a, b

Hanselman, D.: Linearly Interpolate Triangulation, Mathworks [code], https://www.mathworks.com/matlabcentral/fileexchange/38925-linearly-interpolate-triangulation (last access: 30 November 2021), 2021a. a

Hanselman, D.: Contour Plot for Scattered Data, Mathworks [code], https://www.mathworks.com/matlabcentral/fileexchange/38858-contour-plot-for-scattered-data (last access: 30 November 2021), 2021b. a

Hsu, J. P. C. and Capart, H.: Onset and growth of tributary-dammed lakes, Water Resour. Res., 44, W11201, https://doi.org/10.1029/2008WR007020, 2008. a

Iverson, R. M.: The physics of debris flows, Rev. Geophys., 35, 245–296, 1997. a, b, c

Iverson, R. M. and Vallance, J. W.: New views of granular mass flows, Geology, 29, 115–118, https://doi.org/10.1130/0091-7613(2001)029<0115:NVOGMF>2.0.CO;2, 2001. a

Iverson, R. M., Logan, M., LaHusen, R. G., and Berti, M.: The perfect debris flow? Aggregated results from 28 large-scale experiments, J. Geophys. Res.-Earth, 115, F03005, https://doi.org/10.1029/2009JF001514, 2010. a

Johnson, A. M.: Physical processes in geology: A method for interpretation of natural phenomena; intrusions in igneous rocks, fractures, and folds, flow of debris and ice, Freeman, Cooper, ISBN 10:0-87735-319-0, ISBN 13:978-0-87735-319-5, 1970. a

Ke, W.-T. and Capart, H.: Theory for the curvature dependence of delta front progradation, Geophys. Res. Lett., 42, 10–680, https://doi.org/10.1002/2015GL066455, 2015. a

Kowalski, J. and McElwaine, J. N.: Shallow two-component gravity-driven flows with vertical variation, J. Fluid Mech., 714, 434–462, https://doi.org/10.1017/jfm.2012.489, 2013. a

Kuster, C. M. and Gremaud, P. A.: Accurately Computing the Shape of Sandpiles, in: Multiscale Optimization Methods and Applications, Springer, 305–312, https://doi.org/10.1007/0-387-29550-X_15, 2006. a, b

Lai, S. Y. J. and Capart, H.: Two-diffusion description of hyperpycnal deltas, J. Geophys. Res.-Earth, 112, F03005, https://doi.org/10.1029/2006JF000617, 2007. a, b

Lai, S. Y. J. and Capart, H.: Reservoir infill by hyperpycnal deltas over bedrock, Geophys. Res. Lett., 36, L08402, https://doi.org/10.1029/2008GL037139, 2009. a

Le Hooke, R. B. and Rohrer, W. L.: Geometry of alluvial fans: Effect of discharge and sediment size, Earth Surf. Proc. Land., 4, 147–166, https://doi.org/10.1002/esp.3290040205, 1979. a, b, c

Liu, K. F. and Huang, M. C.: Numerical simulation of debris flow with application on hazard area mapping, Computat. Geosci., 10, 221–240, https://doi.org/10.1007/s10596-005-9020-4, 2006. a, b

Liu, K. F. and Mei, C. C.: Slow spreading of a sheet of Bingham fluid on an inclined plane, J. Fluid Mech., 207, 505–529, https://doi.org/10.1017/S0022112089002685, 1989. a

Lobkovsky, A. E., Smith, B. E., Kudrolli, A., Mohrig, D. C., and Rothman, D. H.: Erosive dynamics of channels incised by subsurface water flow, J. Geophys. Res.-Earth, 112, F03S12, https://doi.org/10.1029/2006JF000517, 2007. a

Lorenzo-Trueba, J. and Voller, V. R.: Analytical and numerical solution of a generalized Stefan problem exhibiting two moving boundaries with application to ocean delta formation, J. Math. Anal. Appl., 366, 538–549, https://doi.org/10.1016/j.jmaa.2010.01.008, 2010. a

Lorenzo-Trueba, J., Voller, V. R., and Paola, C.: A geometric model for the dynamics of a fluvially dominated deltaic system under base-level change, Comput. Geosci., 53, 39–47, 2013. a

Mangeney, A., Bouchut, F., Thomas, N., Vilotte, J.-P., and Bristeau, M.-O.: Numerical modeling of self-channeling granular flows and of their levee-channel deposits, J. Geophys. Res.-Earth, 112, F02017, https://doi.org/10.1029/2006JF000469, 2007. a

Mangeney, A., Roche, O., Hungr, O., Mangold, N., Faccanoni, G., and Lucas, A.: Erosion and mobility in granular collapse over sloping beds, J. Geophys. Res., 115, F03040, https://doi.org/10.1029/2009JF001462, 2010. a

Meng, X. and Wang, Y.: Modelling and numerical simulation of two-phase debris flows, Acta Geotech., 11, 1027–1045, https://doi.org/10.1007/s11440-015-0418-4, 2016. a

Mitchell, N. C.: Morphologies of knickpoints in submarine canyons, Geol. Soc. Am. Bull., 118, 589–605, https://doi.org/10.1130/B25772.1, 2006. a

Murata, J.: Flow and deformation of fresh concrete, Matériaux et Constructions, 17, 117–129, https://doi.org/10.1007/BF02473663, 1984. a

Murillo, J. and García-Navarro, P.: Wave Riemann description of friction terms in unsteady shallow flows: Application to water and mud/debris floods, J. Comput. Phys., 231, 1963–2001, https://doi.org/10.1016/j.jcp.2011.11.014, 2012. a, b

Ni, W. J. and Capart, H.: Groundwater drainage and recharge by networks of irregular channels, J. Geophys. Res.-Earth, 111, F02014, https://doi.org/10.1029/2005JF000410, 2006. a

O'Brien, J. S.: 2-Dimensional Flood Routine Model Manual, Version 2006.01, FLO-2D Inc., Nutrioso, AZ, USA, 2006. a

O'Brien, J. S., Julien, P. Y., and Fullerton, W.: Two-dimensional water flood and mudflow simulation, J. Hydraul. Eng., 119, 244–261, 1993. a, b

Pudasaini, S. P.: A general two-phase debris flow model, J. Geophys. Res.-Earth, 117, 002186, https://doi.org/10.1029/2011JF002186, 2012. a, b

Pudasaini, S. P. and Fischer, J.-T.: A mechanical model for phase separation in debris flow, Int. J. Multiphas. Flow, 129, 103292, https://doi.org/10.1016/j.ijmultiphaseflow.2020.103292, 2020. a, b

Pudasaini, S. P. and Hutter, K.: Avalanche dynamics: dynamics of rapid flows of dense granular avalanches, Springer Science & Business Media, ISBN 13:978-3-540-32686-1, ISBN 10:3-540-32686-3, 2007. a

Rocha, F., Johnson, C., and Gray, J.: Self-channelisation and levee formation in monodisperse granular flows, J. Fluid Mech., 876, 591–641, https://doi.org/10.1017/jfm.2019.518, 2019. a

Roussel, N.: A thixotropy model for fresh fluid concretes: Theory, validation and applications, Cement Concrete Res., 36, 1797–1806, https://doi.org/10.1016/j.cemconres.2006.05.025, 2006. a

Savage, S. and Iverson, R.: Surge dynamics coupled to pore-pressure evolution in debris flows, in: Proc. 3rd Int. Conf. on Debris-Flow Hazards Mitigation: Mechanics, Prediction, and Assessment, edited by: Rickenmann, D. and Chen, C., Citeseer, Millpress, Rotterdam, Netherlands, 503–514, https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=27c51b90e843fe21f2ccd1709de44fba611535de (last access: 30 August 2022), 2003. a

Savi, S., Tofelde, S., Wickert, A. D., Bufe, A., Schildgen, T. F., and Strecker, M. R.: Interactions between main channels and tributary alluvial fans: channel adjustments and sediment-signal propagation, Earth Surf. Dynam., 8, 303–322, https://doi.org/10.5194/esurf-8-303-2020, 2020. a, b, c

Scheidl, C., Rickenmann, D., and Chiari, M.: The use of airborne LiDAR data for the analysis of debris flow events in Switzerland, Nat. Hazards Earth Syst. Sci., 8, 1113–1127, https://doi.org/10.5194/nhess-8-1113-2008, 2008. a

Tai, Y., Heß, J., and Wang, Y.: Modeling Two-Phase Debris Flows With Grain-Fluid Separation Over Rugged Topography: Application to the 2009 Hsiaolin Event, Taiwan, J. Geophys. Res.-Earth, 124, 305–333, https://doi.org/10.1029/2018JF004671, 2019. a, b, c

Takahashi, T.: Debris Flow, IAHR Monograph, Balkema, Rotterdam, ISBN 10:9054101040, ISBN 13:978-9054101048, 1991. a, b

Tombarevic, E., Voller, V., and Vušanovic, I.: Detailed CVFEM algorithm for three dimensional advection-diffusion problems, CMES-Comp. Model. Eng., 96, 1–29, 2013. a

Tregaskis, C., Johnson, C. G., Cui, X., and Gray, J. M. N. T.: Subcritical and supercritical granular flow around an obstacle on a rough inclined plane, J. Fluid Mech., 933, A25, https://doi.org/10.1017/jfm.2021.1074, 2022. a

Voller, V. R.: Basic control volume finite element methods for fluids and solids, vol. 1, World Scientific, ISBN 13:978-981-283-498-0, ISBN 10:981-283-498-2, 2009. a, b

Voller, V. R. and Paola, C.: Can anomalous diffusion describe depositional fluvial profiles?, J. Geophys. Res.-Earth, 115, F00A13, https://doi.org/10.1029/2009JF001278, 2010. a

Whipple, K., Parker, G., Paola, C., and Mohrig, D.: Channel dynamics, sediment transport, and the slope of alluvial fans: Experimental study, J. Geol., 106, 677–693, https://doi.org/10.1086/516053, 1998. a, b, c

Winslow, A. M.: Numerical solution of the quasilinear Poisson equation in a nonuniform triangle mesh, J. Comput. Phys., 1, 149–172, https://doi.org/10.1016/0021-9991(66)90001-5, 1966. a

Yuhi, M. and Mei, C. C.: Slow spreading of fluid mud over a conical surface, J. Fluid Mech., 519, 337–358, https://doi.org/10.1017/S0022112004001478, 2004. a, b, c, d, e

Zhao, M., Salter, G., Voller, V. R., and Li, S.: Can the growth of deltaic shorelines be unstable?, Earth Surf. Dynam., 7, 505–513, https://doi.org/10.5194/esurf-7-505-2019, 2019. a

Zhou, G. G. and Ng, C. W.: Dimensional analysis of natural debris flows, Can. Geotech. J., 47, 719–729, https://doi.org/10.1139/T09-134, 2010. a, b

- Abstract
- Introduction
- Governing equations
- Numerical method
- Critical slope
- Flux limiter
- Analytical solutions
- Numerical model evaluation
- Comparisons with field and laboratory data
- Conclusions
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References

- Abstract
- Introduction
- Governing equations
- Numerical method
- Critical slope
- Flux limiter
- Analytical solutions
- Numerical model evaluation
- Comparisons with field and laboratory data
- Conclusions
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References