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

. 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 5 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.


Introduction
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 incur :::: 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.
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.
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 axi-symmetric 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.
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 3-fold, 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.Over the long term, for ::: 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), andSavi et al. (2020).Our focus here, however, is on the formation of fresh deposits, possibly over a pre-existing 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.In Sect.8, finally, ::::::: 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.

Governing equations
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 by z(x, y, t) the time-varying surface elevation during flow, and by z b (x, y) the underlying bed topography.The corresponding profiles are shown on Fig. 2b on a schematic section.
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.

Numerical method
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 andPatankar (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. 3 b).This will result in i = 1, 2, . . .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 z(x, y, t), to be determined.To evaluate the values of z b and 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 , 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 Ax , n Ay 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. 3 c,d.Thus the control volume around node i has j = 1, 2, . . .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, A CV,i ::::::: 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 q j = (q xj , q yj ) 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 n = (∆y/ℓ, −∆x/ℓ).This provides us enough information to fully approximate the discharge in Eq. ( 9) in terms of the current nodal values of zi 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: where Q in,i accounts for source node points where material is added to the domain.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. (

Critical slope
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., |∂ z/∂x| ≤ S c ), the mixture remains in static equilibrium, with τ ≤ σ tan ϕ + τ 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); (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 2-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 Eq : s. : ( 14) and ( 15) into the force balance relations, Eq. ( 13), we obtain the following expression for the shear stress 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 = z − z 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. with ::::::::::: observations.

Flux limiter
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 1 3 (z A − z bA )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.

Analytical solutions
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 everywhere aligned with this transect, then the surface profile will satisfy the simpler equation with coordinate x taken along the transect direction.In this expression, the plus operators denote downhill deposition (z and z b decreasing in the same direction), and the minus operators denote uphill deposition (z and z b decreasing in opposite directions).
Substituting z = z b + H, the equation becomes an ODE for the deposit thickness For the special case in which the bed slope ∂z b /∂x = tan β 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: clarify model properties, verify the numerical method, and calibrate material parameters when comparing model results with field and laboratory data. 1 tanφ As an example, analytical solutions for the centerline profiles of cohesive-frictional deposits over an inclined plane are illustrated in Fig. 5. Deposits of different heights are shown, for material supplied at a single point corresponding to the apex of each deposit ::: For :::: each :::: case ::: and :::::: deposit :::::: height :: H, ::: we :::::: supply ::::::: material :: at : a ::::: single ::::: point :::::::::::: corresponding :: to ::: the :::: apex :: of :::: each ::::::: deposit.B ln (B))/A 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 ϕ centred 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.

Numerical model evaluation
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.

Bed profile
Analytical solution Computational result

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); 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.
x 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 ratio H/h ::::::::: 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 z = 0.1h, normalized by the analytical value r 10 = 1.628 m.

Comparison with field profiles
Coussot et al. (1996) observed six natural debris flow deposits in the French Alps.By categorizing these deposits by their fines fraction (ratio of particles whose diameter is less than 40 µm to total solid volume), they found that debris flow deposits with a low fines fraction (< 1%), at Bourgeat, Le Bez and Ste-Elisabeth, exhibit nearly straight profiles, whereas debris flow deposits with a high fines fraction (10%-15%), at Les Sables, St-Julien, and Mont Guillaume, exhibit significant snout-like toes.Coussot et al., therefore, focused on the latter case to test their model involving only the effect of yield stress.For each deposit with a high fines fraction, 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 τ Y /(ρg).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 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 fines fraction.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 τ Y /(ρg)).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 τ Y /(ρg).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 τ Y /(ρg).
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.

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 bevelled wood panels.Different from alluvial fan experiments (Le Hooke and Rohrer, 1979;Whipple et al., 1998;Delorme et al., 2018;Savi et al., 2020), involving 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.[-] [-] Field data Fitted analytical profile 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.
Two series of tests were conducted: canyon-plain experiments (T01-T04), using the geometry shown in Fig. 9a, and canyonvalley 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 degrees 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% kaolin, 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 respectively to 3 and 6 degrees.For these runs, the topography was scanned with the laser oriented perpendicular to the canyons, and the resulting DEM data have resolution 2 mm x 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 degrees, connect at right angles to a wide trapezoidal channel of longitudinal inclination equal to 3 degrees.Two of the canyons (TC1 and TC2) connect on the right side, and one on the left (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 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 3 tributaries.For run T15, deposits were formed in three separate stages.In the first stage, deposits were formed as in run T13 by supplying * The parameter ranges in documented field cases are collected or calculated from the data of Iverson (1997) and Zhou and Ng (2010).
** Viscosity for the experimental interstitial fluid is estimated from a set of viscometer measurements.

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 centrelines 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 τ Y /(ρg) remain to be determined.For this set of experiments, some variability in material properties was caused by uncontrolled variations in moisture in the kaolin.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 τ Y /(ρg) = 0.115 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 z(x, y, 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.For the simulations, numerical stability is insured by choosing for the time step the constant value ∆t = ∆ℓ 2 /4ν * = 0.004 s.Provided that the diffusivity coefficient ν * is chosen sufficiently large, we have checked that the final shape of the deposit does not vary with the value adopted for ν * , or with the rate at which material is supplied at 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 reasonably well the measured topography of the deposits.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 well the outer boundaries of the deposits, both along the steep canyon and valley sides, and over the mildly inclined floor.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, back to steep snouts at the downstream toes.In both the simulations and experiments, furthermore, 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.

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 τ Y /(ρg) = 0.344 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, and a constant time step ∆t = ∆ℓ 2 /4ν * = 0.00625 s. : .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 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 well the deposit thickness distributions obtained in the different experiments.
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.

Conclusions
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 fines 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 applicable also to deposits with a low fines 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 pre-existing deposits.Using material parameters calibrated from one or more transects, the model is found to reproduce well the measured topography in all cases.
Deposits 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, the model is nevertheless subject to various limitations.

Figure 3 .
Figure 3. Global and local mesh geometry: (a) The discretized domain and elements; (b) a triangular element divided by the segments connecting the centroid and the midpoint of each side; (c) the control volume and the region of support of an internal node; (d) those of a node on the boundary.

Figure 4 .
Figure 4. Force balance of a small piece of material on a fixed bed whose local gradient has a value equal to tan β and direction pointing towards ξ.

Table 1 .
Chen et al. (2022)size on model accuracy and computational time.Avg.elementsize[m]# of elements H/h :::::::: (h − H)/h : (R10max − R10min)/r10 Computational time[s]To further test the model, in the 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/output data discussed in this section are available inChen et al. (2022).