Voellmy-type mixture rheologies for dilatant, two-layer debris flow models

We formulate and test different Voellmy-type mixture rheologies that can be introduced into two-layer debris flow models. The formulations are based on experimental data from the Swiss Illgraben test site as well as on mathematical constraints in steady flow conditions. In agreement with the ideas of Iverson, we show that the uniform, fixed rheological models cannot accurately represent the changing frictional resistance when debris flows undergo spatial and temporal changes in solid–fluid composition. Indeed, the experimental results of Illgraben indicate that flow friction decreases with increasing volumetric fluid concentration; however, the degree of reduction depends on both the pore pressure and the solid particle agitation. The interplay between these processes makes friction in debris flows highly nonlinear and difficult to quantify. Changing the friction according to the flow composition must be carefully executed, because it can lead to numerical instabilities, which is a recurrent problem in two-layer debris flow models. We test the different rheological formulations using a real event documented with differential topographic data collected using unmanned aerial vehicles (UAVs). The model is able to reproduce the correct erosion pattern and exhibit the right density profile. The event includes de-watering at the front and deposition of sediment, which causes a change from debris flow to debris flood or hyperconcentrated flow, which indicates that two completely different flow states can be modeled with a single Voellmy-type mixture rheology.


Introduction
The increasing application of two-phase debris flow models in hazard engineering (Kowalski 2008;Pudasaini 2012;Kowalski and McElwaine 2013;Iverson and George 2014;George and Iverson 2014;Takahashi 1981) has renewed interest and subsequently research into the debris flow rheology problem, see Fig. 1.Presently debris flow calculations are performed using fixed rheologies using frictional parameters that vary strongly depending on the debris flow type (Graf et al. 2019).This limits the predictive power of model calculations and their application in hazard engineering.The evolving behavior of debris flows is simply too complex to be modeled by a fixed rheology (with regards to space and time) that ignores spatial and temporal variations of fluid content, pore pressure, and particle agitation.Iverson (2003) or more recently Hungr (2000) termed the impossibility of finding a fixed rheology, the rheology myth, and proposed Coulomb mixture theories to model debris flows of variable solid-fluid composition.
In this paper, we investigate two Coulomb mixture theory formulations of the Voellmy type.Our use of Voellmy-type models is motivated by the fact that they are widely applied in debris flow practice (at least in Switzerland, Graf et al. 2019), and therefore, empirical parameter ranges are known.Because Voellmy-type models contain only two parameters for each layer, the total number of parameters remains small.Our primary motivation, however, is to show that improved simulation results in terms of predictive power are possible.In this work, we aim at demonstrating, using shear and normal force measurements at the Illgraben test site (McArdell and Sartori 2020;McArdell et al. 2007), that there is also a good experimental foundation for this rheological approach.We investigate the influence of rheology on debris flow composition, de-watering, entrainment, and prediction of channel outbreaks (Cuomo et al. 2014(Cuomo et al. , 2016)).
The general mathematical form of the Voellmy-type mixture formulations is: where ⃗ S i denotes the frictional resistance of the i-th layer (i = 1 stands for the mixed debris flow layer and i = 2 for the free water layer, for more details, see "Two-layer debris flow model with dilatancy" section).As in the case for single-layer/layer Voellmy-type models, the frictional resistance is split into two parts: the Coulomb friction part ⃗ S ,i of the i-th layer and the hydraulic, velocity-dependent friction ⃗ S ,i .Here, we adopt the Swiss convention of denoting the Coulomb friction with the greek symbol and the hydraulic friction with the symbol .Coulomb friction depends on the normal stress N i and is given by: and S ,i is given by: Where i is the turbulent friction coefficient for the layer i and ⃗ v i is the velocity of the i-th layer.The vector ⃗ e i is a unit vector along the flowing direction of the i-th layer, ⃗ e i = ⃗ v i ‖v i ‖ .To test the rheological formulations, it is necessary to implement them into a two-layer debris flow model.For this purpose, we use the approach recently proposed by Meyrat et al. (2021).We choose to work with a depth-averaged model since it has been shown to accurately reproduce the spatial and temporal flow composition measured in actual debris flows at the Swiss Illgraben test station (McArdell et al. 2007;Schlunegger et al. 2009;Badoux et al. 2009). (1) Landslides 20 • (2023)   Review Paper This includes the simulation of the measured density distribution from leading edge of the debris flow to the flow tail (Meyrat et al. 2021).It contains both a free-fluid layer and mixture layer (solid-bonded fluid) and therefore allows the determination of excess pore pressures as well as an agitation-driven dilatancy terms arising from the shearing of the granular solid.As such, it fulfills the basic debris flow model requirements formulated by several authors including Takahashi (1981), Iverson and George (2014) and Iverson (1997).The rheological formulations presented here are formulated in general terms such that they can be applied in any two-layer model.Moreover, we wish to separate as much as possible the proposed rheological formulations from the specific two-layer debris flow model that we apply.The natural state of a debris flow is unsteady, especially when considering flow in steep, twisting torrents in complex mountain terrain.Bed erosion also serves to disrupt any possibility of steady flow.Nonetheless, we decide to test the Voellmy-type mixture rheologies both in steady and unsteady conditions.Indeed, a steadystate analysis of the mixture formulations is necessary, especially to test the overall model performance (and the correctness of the numerical implementation since the model results must converge to analytically derived steady solutions).Since steady solutions exist in most, if not all, two-layer debris flow models, the steady-state analysis facilitates a comparison between different two-layer model approaches.It also helps clarifying the underlying physics of the different mixture formulations.
In the following section, we present the governing equations of our two-layer approach.We discuss the physical ideas underpinning the model.This presentation serves to identify the friction terms in the momentum balance equations.In "Voellmy-type mixture rheologies" section, we present the two basic rheological formulations based on pore pressure and dilatancy concepts.In "Steady state solutions" section, we investigate the numerical stability of the proposed formulations.Documented real event simulations are found in "Case study: Ritigraben, Switzerland" section.The paper is rounded-off with some concluding remarks and an outlook to future work in "Discussion and conclusions" section.

Two-layer debris flow model with dilatancy
We model debris flows with a two-layer, depth-averaged formulation (Meyrat et al. 2021), see Fig. 2. The debris flow contains two material components: a solid component (subscript s) consisting of coarse granular sediments (e.g., boulders, cobbles, and gravel), associated with a density s , and a fluid component (subscript f) consisting of fine sediment likely to behave as suspended sediment (e.g., sand, silt, clay), hereafter referred to as the muddy fluid content, the density of which is denoted by f .
The solid and fluid components are divided into two layers.The first layer (subscript 1) contains the granular solid material and a part of the fluid.The fluid is contained in the interstitial space between particles and is assumed to be bonded to the solid particles.The second layer (subscript 2) is formed by the fluid which can flow independently from the first layer.The total mass per unit of area (all the quantities in the following are defined per unit area A) of the debris flow M is the sum of the mass of layer 1 (solid particles and inter-granular fluid) and mass of layer 2 (free fluid), where h f ,1 and h f ,2 represent the amount of fluid in the first (interstitial fluid) and second (free fluid) layer, respectively, and h f defines the total amount of fluid in the flow.We write h s instead of h s,1 because the solid always belongs to the first layer only.M f ,i = f h f ,i is the fluid mass contained in the i-th layer.These assumptions lead to a system of three depth-averaged mass balance equations, (4) where ⃗ ∇ is the divergence operator in Cartesian coordinates.The vectors ⃗ v 1 and ⃗ v 2 represent the depth-averaged velocity of the first and second layers, respectively.For notation purposes, we do not denote the depth-averaged quantities differently in the following of this paper.The right-hand side of the inter-granular fluid and free fluid equations contains the term Q (Meyrat et al. 2021), which is the mass exchange rate between the inter-granular and free fluid driven by dilatant actions in the solid matrix (see following).The equations relevant for the first phase (Eqs.5 and 7) also contain the erosion rate, denoted by E. The density e characterizes the density of the entrained mass.In order to compute the erosion rate E, we began with the model introduced in Frank et al. (20152017) which is based on careful field measurements (Schürch et al. 2011) and well-documented events.The model was subsequently adapted for the two-layer model.In the modified model, the erosion rate E is no longer uniform (as in Frank et al. 2015Frank et al. , 2017)), but a function of the flow composition.If E s represents the erosion rate for a purely solid flow and E f for a completely fluid layer, the erosion rate E is expressed as follows: where is the volumetric fluid fraction of the first layer, defined as: The term E is absent in Eq. ( 6), because the second layer, i.e., the free fluid, cannot entrain solid material by definition.Solid material must be entrained by the first layer.Unlike in Cuomo et al. (20142016) in which the erosion rate is a function of the flow depth ( 5) solid and inter-granular fluid (6) and velocity, the solid and fluid erosion rates E s and E f are constant and uniform.However, the maximum potential erosion depth is computed from the shear stress (Frank et al. 2015(Frank et al. , 2017)), which depends on the flow height, velocity, and composition.We emphasize that the pseudo-variable h 1 does not represent the physical height of the first phase (the one that a sensor would measure).This variable is introduced to simplify the mass and momentum conservation equations.The pseudo-height h 1 is defined as the sum of the solid and bonded fluid mass (equal to first layer mass), normalized by the solid density s : This relation can be found by noting that the total mass of the first layer M 1 can be expressed as the sum of the solid material M s and the interstitial fluid mass M b .As the heights h s and h f ,1 represent the height of the solid and inter-granular fluid components in the first layer, respectively (these quantities can be considered as the solid and fluid volumetric parts of the debris flow mass), we find that, Where h sb = h s + h f ,1 represents the physical height of the first layer and M 2 = M f the mass of the second layer, see Fig. 2. All the fluid which flows above the first layer is considered free, and subsequently, it can escape the matrix of solid particles, allowing the debris flow to de-water.As the solid particles settle and the interstitial space between particles collapses, fluid can escape from the solid part of the flow.
At first instance, the division of masses using a pseudo-height h 1 appears cumbersome.h 1 is however associated with a constant and uniform density s (Eq.10), in contrast with h sb which is associated with the "real" first layer density 1 which evolves with the dilatant action of the solid matrix.Therefore, using h 1 simplifies the momentum balance equations by the density.They are given for the two layers (Savary and Zech 2007; Mandli 2011) by: (10) Fig. 2 A sketch of the debris flow model.The flow consists of solid material in the form of boulders (granules) as well as two types of fluid.The bonded/interstitial fluid is located between the boulders and is "fixed" to the solid.The free fluid is located above the granular solid/fluid mixture and moves independently.In the numerical simulations, the debris flow is divided into a series of volumes V (black rectangles) Landslides 20 • (2023)

Review Paper
Where the vectors ⃗ v 1 and ⃗ v 2 represent the depth-averaged velocity of the first and second layers, respectively, and b ∶= b(x, y) denotes the bottom topography.The symbol ⊗ denotes the tensor product, and I is the two-dimensional unity matrix.The left side of Eqs. ( 13) and ( 14) is the total variation of the momentum with respect to time, including the effect of gravitation and the influence of each layer on the other (Savary and Zech 2007;Mandli 2011).The right side represents the change in momentum due to external forces (excluding gravitation).⃗ S i corresponds to the shearing forces acting on the i-th layer (we divide it by the density s for the first layer and f for the second one, because Eqs. ( 13) and ( 14) are defined per unit of density).Different Voellmy-type mixture models for the shearing stresses ⃗ S i are provided in the next chapter.Contrarily to Pudasaini and Fischer (2016) and Pudasaini and Krautblatter (2021), we do not account for any momentum production or softening of the basal topography due to erosion.Consequently, erosion is not directly present in the momentum balance equations, Eqs. ( 13) and ( 14).However, it influences the mass balance and the mixture layer composition, which governs the shear resistance of the flow.Therefore, although we do not explicitly consider momentum production, entrainment changes the entire flow dynamics.
The vector ⃗ P is the rate of momentum exchange associated with the mass exchange Q, The mass exchange Q is a result of the dilatant actions of the solid particles in layer 1.Under interactions with the rough bed of the channel, the solid matrix can expand its volume during flowing, (13) leading to different flowing configurations, associated with different densities (see left sketch of Fig. 3).Note that, even if the solid mass is conserved, according to Eq. ( 11), the first layer density can vary if the inter-granular fluid concentration evolves: When the flow is at rest, the first layer is in the so called co-volume configuration, the height of which is denoted h 0 , which corresponds to the configuration where the solid matrix is completely collapsed, see right sketch on Fig. 3.However, potential energy is required to rise up the center of mass of the solid matrix.We call it the configurational energy, D, and it can be expressed, considering buoyancy, as: where z denotes the vertical displacement of the center of mass from the collapsed configuration to a dilated configuration, and gravity g z is the slope-perpendicular gravity component.Changes in the center of mass of the solid matrix imply changes in the void space.When the void space increases, free fluid fills the additional space between the particles.Conversely, if the solid matrix collapses, interstitial bonded fluid is transformed into free fluid and will be squeezed out of the first layer.Therefore, changes in the solid matrix center of mass will always be accompanied by fluid mass exchanges between the interstitial bonded fluid and the free fluid.These exchanges are denoted Q and are responsible for the evolution of the flow composition.
We postulate that the configurational energy D is governed by a simple production (parameter ) and decay (parameter ) term (Meyrat et al. 2021 2021).We define three variables associated with the solid mass, h s , h 0 , and h sb .The height h s is the volume of the solid particles in the flow; h 0 represents the reference height of the non-dilated mass; we call it the co-volume (it is different from h s , because we consider that, even in the nondilated configuration, the void space is not zero), whereas h sb represents the dilated height of the solid mass, "Two-layer debris flow model with dilatancy" section.We refer to h 0 as the co-volume, in an analogy to Van der Waals work on non-ideal gasses with large molecules and cohesion (Visser 1989) The quantity W f represents the shear work.That is, the change in the energy of configuration D (dilatation) is directly related to the shear work, in accordance to Reynolds (1885).The parameter defines the fraction of the shear work that produces a dilatation, whereas the parameter defines how quickly the dilatation collapses in the absence of shear, due to energy dissipation caused by shearing between particles.The balance between the production of D and its decay defines the degree of saturation in the debris flow, as this defines the amount of void space (inter-granular water) in the moving solid.For more details, the reader is referred to Meyrat et al. (2021).

Voellmy-type mixture rheologies
The two parameters Voellmy formulation is popular in natural hazard mitigation because the Coulomb friction determines the debris flow runout, as it defines the critical slope angle at which the flow begins to decelerate tan = (Graf et al. 2019).The hydraulic friction defines the steady flow speed of the movement.Thus, with only two parameters, we can define the approximate runout distance and steady flow velocity.Default values for single-layer/layer debris flows are known (approximately = 0.20 and = 200 m/s 2 ).
However, these values can change significantly in practical casestudies, because they depend on the fluid content of the debris flow and torrent conditions, such as grain size distribution.With the aim of reducing the range of possible values of the two Voellmy parameters, we now introduce two models, the Voellmy mixture model (denoted as VM), "Voellmy mixture model (VM model)" section and the Voellmy mixture model including pore pressure effects (denoted PP), "Voellmy mixture model including pore pressure effects (PP model)" section.In these models, the rheology is not uniform and constant anymore but evolves with the flow composition.Before going to the details of the different models, here, we make some general comments: The rheological formulation of any two-layer/phase model with an entirely fluid layer/phase (the free fluid one here) has to ensure the following consistency conditions: where is the volumetric fluid fraction, as defined in Eq. ( 9).These limit conditions ensure the following: if the first layer density becomes close to the density of the second layer (meaning that it contains no solid material), the frictional shear forces acting on the first layer must be similar to the second layer.The density of the first layer can vary, from a rocky density ( ≈ 2000-2200 kg/m 3 , debris flow front) to the density of the muddy fluid ( ≈ 1000-1300 kg/m 3 , debris flow tail).Therefore, if the first layer acquires the density of the muddy fluid (for example, at the tail of the flow), the consistency condition enforces a rheological uniformity between the layers.Because bed entrainment can suddenly change the solid and fluid concentrations of the debris flow anywhere, the consistency condition physically constrains the effects of bed erosion on the evolving flow properties of the debris flow.
The assumption that the debris flow rheology is not constant and uniform but evolving with the flow composition is supported by the experimental data from Illgraben.The Illgraben debris flow test site (Fig. 4) is located near Leuk, Canton Wallis, Switzerland (Schlunegger et al. 2009;Badoux et al. 2009;McArdell and (19) Sartori 2020).Since 2005, the Illgraben torrent has been instrumented with a rectangular force plate (area A = 8 m 2 ) that measures shear and normal stresses at the base of a passing debris flow.
A laser sensor located above the plate measures the total debris flow height h as the flow passes over the plate.The force plate is located at the end of a 5 km long torrent (in orange on Fig. 4) that has an average slope of 5 o at the downstream end of the channel at the location of the force plate.The torrent is fed by a large (9.5 km 2 ), steep catchment zone (in blue on Fig. 4), which supplies the channel.The Illgraben sub-catchment supplies most of the water and sediment, whereas the Illbach catchment supplies mainly water (McArdell et al. 2007).
From the normal stress and height measurements, it is possible to compute the flowing density or equivalently the fluid volume fraction (Meyrat et al. 2021).It is therefore possible to investigate the link between the flow composition and the intensity of the shear forces.The ratio between the shear and the normal stress (S/N) is a good representation of the shear force intensity.In Fig. 5, the ratio S/N is plotted as a function of the volumetric fluid concentration (which is simply equal to 100 × f ).In this figure, we can see clearly that the ration S/N depends strongly on the flowing fluid saturation.
In order to take into account the effects of the evolving debris flow composition in the first layer on its rheology, we define four Review Paper rheology coefficients: s and s which are relevant for the densest configuration, i.e, the co-volume configuration while f and f if the flow is composed of fluid only.As the second layer is always composed only by fluid, its coefficients are uniform and constant, i.e., 2 = f and 2 = f .

Voellmy mixture model (VM model)
One method to account for the reduction of friction with increasing fluid concentration (Fig. 5) is to express the rheology coefficients and as a linear function of fluid fraction .Therefore, we assume that the shearing coefficients are given by the weighted average of s and s of the flow composition.Formally, Equation (20) combined with Eqs.(1), (2), and (3) satisfy the consistency condition given by Eq. ( 19).In general, we assume that for the muddy fluid (i = 2), ( 20) Thus, we have for the mixture stresses and Voellmy mixture model including pore pressure effects (PP model) Effective stress and pore pressure concepts have been applied to model debris flow motion (Pudasaini 2012;Iverson and George 2014;Iverson 1997Iverson , 2005)).If we let p be the fluid pressure, then in depth-averaged models, the effective stress N eff is given by the difference between the normal pressure and the fluid pressure, On the assumption that the pressure p is hydrostatic, the Coulomb part of the shear stress is, (26) Note here that we do not include the normal stress induced by free fluid, i.e., N f = f g z h f ,2 , in the computation of the pore pressure.It is not required because N f should be added to both normal stress (N) and pore fluid pressure (p).Therefore, the free fluid normal stress contribution will vanish in the computation of ⃗ S ,1 in Eq. ( 27) as only the subtraction of N by p is relevant.Another interesting feature of the model is that the consistency condition given by Eq. ( 19) is satisfied if and only if 2 = 0 .Therefore, In this formulation, Eq. ( 21) is still valid; therefore, the -terms are the same as in the preceding approach,

Steady state solutions
The purpose of this section is to derive steady-state solutions of the two-layer debris flow equations with dilatancy.Because steadystates represent the balance between the driving and resisting forces acting on the debris flow body, the solutions are dependent on the different rheological formulations.For mixture models, the rheological formulations are, in turn, a function of the fluid fraction .In our model, the fluid fraction is dependent on the available pore space between the churning rocks and boulders in the debris flow body and therefore intimately related to flow dilatancy.The demonstration that steady state solutions even exist for a two-layer mixture model for a specific fluid volume fraction is important.It indicates that states of steady flow are associated with states of constant flow density (zero changes in dilatancy and then in pore space) and therefore unchanging propagation speeds.The reverse statement appears also to be true: with changing speeds, the flow densities cannot remain constant, and therefore the evolution of debris flow speed, from the tail to the front of the debris flow, is always related to a change in flow density (and therefore a change in fluid volume fraction ).Changing flow compositions are always accompanied by mass transfers between the bonded and free fluid components.The rheological mixture models must allow for steady states to exist, if only theoretically.Steady-state results can also be used to test the numerical implementation of the mixture models.
In this section, we study the mathematical solution of a sliding rigid block along an inclined plan with a constant slope angle .The block has a fluid volume fraction .A steady state has many properties, one of them is that the acceleration is zero, i.e., ̇⃗ v = 0 .For our purpose, we study the flow in a single direction x.Then the equation of motion of a sliding block is, with respect to the axis êx and êz , with S the shear force in the x-direction.This implies, ( 27) where ⃗ v is the velocity of the center of mass (i.e., always parallel to the slope direction êx ).Note that as S and N are stresses, indicating that the mass m in Eqs. ( 30) and ( 32) is defined per unit of the basal area A, i.e., m = 1 h 1 .This is not a problem for the numerical simulation and for this mathematical analysis, since the basal area A is in both cases constant and uniform.Note that in Eq. ( 32), and have to be changed for ( ) and ( ) for the VM Model and N has to be changed for N eff for the PP Model.The Froude number of the flow block is a function of the rheological formulation.It can be written for the VM model ( Fr l ) and the PP model ( Fr p ) as, with ( ) and ( ) given by Eq. ( 20).The entire derivation is presented in Meyrat et al. (2021).
To test the model implementation, we compare calculated and analytical Froude numbers.For the simulations, we use a smooth plane with a constant angle; this way, the flow can reach a steady state.A block release is used to initiate the flow, and erosion is not taken into account.In order to compute the convergence of the numerical simulations to the analytical derivation, we first compute the absolute value of the relative difference between the numerical outputs and the mathematical value of the Froude number, averaged on each cell containing mass, we call it "percentage difference."If F n and F a denote the numerical and analytical value of the Froude number and if N is the number of cells containing material: with i running over the cells containing mass.In Fig. 6a (VM model) and Fig. 6b (PP model), we plot the percentage difference as a function of the time.It is clear that the numerical outputs converge to the analytical value.Another aspect of a steady state is that the gravitational force is in an exact balance with the shearing force, independent of the formulation of the shearing processes we consider.This fact is depicted in Fig. 6c (VM model) and Fig. 6d (PP model).The red line represents the gravitational stress, and the colored markers (the blue markers represent the front of the flow while the yellow markers represent the tail, see Fig. 7 plot c and d) depict the shearing stresses.Except for the leading edge of the flow, the two curves match perfectly.This difference in the front is linked with numerical effects due to the shock-wave form of the front of the debris flow.
This equivalence between shearing and gravitational stress (or force) in a steady state can be written in mathematical terms as follow: Landslides 20 • (2023)

Review Paper
which implies that the derivative of S with respect to time is zero.Therefore, the shear work rate Ẇf , in Eq. ( 18), can be written as: Where x represents the position of the flow block on the x-axis.
Using the expression for Ẇf , Eq. ( 18) can be solved exactly.We find in the limit of t → ∞: combining with Eq. ( 17) and defining Γ = , we obtain, for the steady-state flow height of a dilatant debris flow: Note that this formula is valid for both rheological formulations.The difference between them is in the steady state's value of the velocity, (37) S = mg sin( )

Dilatancy
given by Eqs. ( 34) and ( 35).In our model, we have an additional assumption, which is that a solid matrix cannot dilate if there is no more free fluid above, i.e., h f ,2 = 0 .Therefore, Eq. ( 40) becomes: Here, Γ = is the ratio between the production and the decay of the configurational energy, see Eq. ( 18).We can see that in a steady state, only the ratio between them is important, and governs the time needed to approach a steady state.This formula holds for the two formulations.It is related to the fact that in a steadystates, the shearing forces are always equal to the driving forces, independently of their specific expressions.In Fig. 7, we show the convergence of the numerical height to the analytical solution.In plot Fig. 7a (VM model) and Fig. 7b (PP model), we plot the absolute value of the relative difference between the numerical and the analytical value of the height (Eq.36).We can see that the difference approaches zero after some time, which represents the fact that both values of the height coincide.In plot Fig. 7c (VM model) and Fig. 7d (PP model), we choose a location in the domain sufficient distanced (41)  37), in a steady state, the shear stress (colored circles) does not depend on the rheological formulation and is equal to the gravitational stress (red lines).The color represents the time; the front is in blue, and as we move towards the tail, the color becomes yellow, see Fig. 7.The very first part of the flow (blue markers) is not exactly in a steady state, due to numerical effects related to the shock-wave form of the front.The simulations are done on a smooth plane with constant angle, using a block release Landslides 20 • (2023)   from the release zone to ensure the flow can reach a steady state.We compare the calculated height from the numerical code and the analytical steady state solution.Both the numerical and analytical curves match well, which indicates that the flow has indeed reached a steady state in its height1 .Steady states are indeed associated with no changes in density.Note that we assume that we always have enough fluid to fill the void space between the churning particles.Analytically, we do not consider the influence of one layer on the other.However, this effect is taken into account when we perform the numerical simulations.The fact that the numerical output converges to the analytical value shows that these effects are negligible.Therefore, the analytical derivation can be considered as valid even for a dilatant, two-layer model with solid-fluid interactions.

Case study: Ritigraben, Switzerland
After checking the mathematical consistency of the computer implementation, we now model a recent event at a different catchment.For this purpose, we choose the torrent of Ritigraben, located between Gräschen and St Niklaus, in Wallis, Switzerland2 .The torrent is 3.75 km long with an elevation drop of approximately 1500 m (Fig. 8).At the top of the Ritigraben torrent lies a rock glacier, which provides sediments to the upper torrent catchment.Subsequently, sizeable debris flow events occur almost every year.Many studies have already been performed at this debris flow site (Lugon and Stoffel 2010;Kenner et al. 2017;Lugon and Monbaron 1998;Stoffel et al. 2005;Mutter and Phillips 2012).
The event of the 7th of August 2021 is of particular interest to us because digital elevation models (DEM) were obtained 2 days before and 2 days after the event.This allows to estimate the mass balance for this particular debris flow (Bindereif 2022), accurately, constraining the initial and entrainment volumes.To compute the volume eroded by the debris flow, the longitudinal profile of the torrent has been divided into several bins of 20 m length each (Frank et al. 2015;Quillard et al. 2017), see Fig. 9.For each bin, the eroded and deposited material was computed by taking the difference between the two digital elevation models.A map of the erosion (negative difference, red) and deposited (positive difference, blue)

Review Paper
volumes can be made (Fig. 9).We find that the total eroded volume reached approximately 10'000 m 3 .In some parts of the torrent, no erosion or deposition was observed.It should be noted that even if the entire eroded volume can be computed, it is still not possible to define separately the volume of solid-fluid entrained.
To initiate the flow, we assume a block release of 1000 m 3 (650 m 3 of solid and 350 m 3 of fluid) (Bindereif 2022) and a hydrograph, to simulate the river flow in the channel.As the event occurred just after a heavy rainfall, the river discharge, normally low, cannot be neglected, especially for an event of this size.Video recordings were used to estimate the fluid volume, see video.We use a 7 m 3 /s discharge rate of river fluid.We performed two sets of simulations, with the VM model and the PP model.When we compare the two simulations, they were always performed with the same value of the model parameters.A list of the parameter values can be found in Table 1, Appendix section.Figure 10a depicts the numerical distribution of erosion along the debris flow channel.Figure 10b, represents the erosion computed for each of the bins shown on Fig. 9.The red line represents the measured data (drone flights), while the blue (VM Model) and green (PP model) curves are obtained from numerical simulations.Both models provide similar results; both are in good agreement with the field data.Both models give a total erosion volume close to 10 � 000 m 3 , which is the one obtained from the filed data.Note that we allowed erosion in our model only in the torrent sections where erosion was observed by field observations.
Figures 11a (VM model) and Fig. 11b (PP model) show the simulated spatial evolution of the flow density.For both rheological formulations, the front (red region) has the highest density, ≈ 2000 kg/m 3 , which is consistent with the fact that it is composed of large blocks, while the tail (blue region) has the smallest density ≈ 1300 kg/m 3 , which is the density of a muddy fluid.Even if the computed front and tail densities are similar for both rheological models, the spatial evolution differs slightly.This fact will be discussed in "Discussion and conclusions" section.
The spatial variation of the friction coefficients and is shown in Fig. 12. Figures 12a (VM Model) and 12b (PP Model), resp.12c (VM Model) and 12d (PP Model), show a decrease in as we go from the front to the tail for an erosion density of 2000 kg/ m 3 (a and b) and 1800 kg/m 3 (c and d, respectively).The values of given by both models are significantly different in the front but become closer as we go towards the tail.We find that the value of strongly depends on the erosion density.The spatial evolution of is depicted in Fig. 12e (VM Model) and 12d (PP Model).As expected, the value of increases when going from front to tail.A clear decrease in the shearing stresses is observed from front to tail.This is due to the decrease in density, and then an increase in fluid content of the flow, from front to tail.Note that is given by Eq. ( 21) for both model and by Eq. ( 20) for the VM model and by the ratio of the Coulomb shear stress, Eq. ( 27), to the normal stress for the PP Model.
In Fig. 13, the ratio of the total shear stress (green), the Coulomb shear stress (red) and the turbulent shear stress (blue), to the normal stress, is plotted as a function of the volumetric volume concentration for two different locations in the debris flow torrent (Fig. 8).The first point (plots a and b) is located in the flat part of the channel In Fig. 14, the maximum solid (Fig. 14a) and fluid (Fig. 14b) heights in the run-out zone are shown.We can see that even if both fluid and fluid flow together in the steep channel, it is not the case anymore when the flow reaches the run-out zone.Indeed, the core of the flow stops and the fluid content is washed-out the solid matrix and can continue to flow downwards the valley.This is what is called by some authors phase separation, Pudasaini and Fischer (2020).

Discussion and conclusions
Measurements from the Illgraben debris flow test site indicate that frictional resistance depends on the fluid-solid composition of the flow (Fig. 5).Non-changing flow compositions are found at the front (high solid concentration of boulders and rocky debris moving with quasi-constant velocity) and at the tail (high concentration of fluid, also moving close to a constant velocity).These flow compositions persist over long periods of debris flow motion and mathematically represent steady flow conditions.A challenging problem in debris flow science is to model both flow compositions within the same model framework (and the same debris flow), without ad-hoc changes of resistance parameters.For this purpose, we have investigated the use of Voellmy-type mixture theories implemented within multi-phase/layer debris flow models.
In a first step, it was necessary to demonstrate the steady state properties of a two-layer debris flow model.The rheological mixture models implemented within the framework of two-layer simulation tools must allow for the possibility of two entirely different flow compositions in the same debris flow.Figure 7 is of particular importance because it reveals the fact that the flow body can indeed reach a constant density in steady state, in which no more mass or momentum exchanges between the layers are possible.We therefore demonstrated that constant steady state velocities (given by the friction parameters) are associated with constant flow heights and therefore constant flow densities.A particularly important formula is,

Review Paper
which relates the total flow height of the debris flow to the slope angle, flow velocity, and parameters governing granular dilatation.This formula provides us with a method to experimentally determine the ratio between two important model parameters governing the production and decay of the configurational (dilatational) energy, Γ = .Because in our model formulation the interstitial fluid content is controlled by the granular dilatations of the solid phase, the frictional resistance is intimately related to this ratio.In the calculation example of a real debris flow, we applied the value Γ = 1 s. (42) The primary difference between the two different mixture models is the dependency of the solid Coulomb friction parameter on the flow density (Eqs.24 and 27).As the flow density (solid composition) varies from front to tail of the debris flow, we expected the two formulations to provide strongly varying results.Surprisingly, the numerical differences are small.We note that the difference in density between the debris flow front ( ≈ 2500 kg/m 3 ) and tail ( ≈ 1200 kg/m 3 ) is between 1300 and 1500 kg/m 3 .Interestingly, this range is numerically close to the muddy fluid density f (≈ 1300-1500 kg/m 3 ).Therefore, the difference between the two approaches will only be in the front of the flow, where the density, used for the VM model, is much larger.Assuming the use of the same values of all friction parameters ( s , f , s , and f ), the computed shear stresses will be larger at the front of the flow for the VM For a given erosion density, the VM model exhibits larger values of in the front than the PP model but these values become similar at the tail.Note that the turbulent friction decreases analogously for both models when going from front to tail because the decrease is governed by the same relationships.
For both mixture formulations, the calculated density values at the leading edge of the flow are similar.They appear to be governed almost entirely by the erosion density.The calculated densities at the tail of the flow are also similar, approximately a value corresponding to a muddy fluid.However, the spatial evolution of the calculated densities differs strongly between the limit and PP models.For the VM model, the density decreases faster behind the front.This implies that the solid fraction is more concentrated at the front in comparison to the PP model.Due to its larger friction, the front acts a bit more like a dam in the VM model.Consequently, the flow height in the front is slightly larger.This fact could have important impacts when simulating overflows and channel outbreaks.In Fig. 15, we depict two overflow cases.Even if this overflow is

Review Paper
relatively small, the overflow is nonetheless significantly larger for the VM model than for the PP model.Therefore, even if the general behavior is similar, minor differences in flow height (and streamwise composition) could play a significant role in hazard analysis.
Our model comparison to the Illgraben measurements also helped identifying a possible problem with the experimental observations.To model the measurements, it was necessary to apply very high values of the -friction.These high -values (low friction values) cannot be applied to model debris flows because they lead to unrealistically high flow velocities.In Illgraben, the shear stress is measured with a shear plate (McArdell et al. 2007;Schlunegger et al. 2009).However, the turbulent -stress depends on the roughness of the terrain (McArdell et al. 2007;Salm 1993).The force plate consists of smooth steel without any additional roughness (Salm (1993) and is smoother than the natural sediment bed of the channel.Therefore, it is possible that we do not measure the entire shear stress because one part of the turbulent -stress is missing in the measurements, as if the flow was sliding over the plate.This hypothesis is supported by the fact that small values of the shear stress are measured for non-zero values of the normal stress (Fig. 5).
In Fig. 13, we plot the ratio of the shear to the normal stress, as a function of the volumetric fluid concentration at two different track locations (i.e., the same plot as Fig. 5) for both mixture models, see Fig 8 .The first location is on an almost flat slope, comparable to the location where the data from Illgraben are obtained (less than 10 • steep), while the second location corresponds to steeper terrain (around 30 • steep).These plots show that the decrease of total friction (green dots) is because of the decrease in Coulomb friction (red dots).The turbulent -friction (blue dots) is almost constant.Note that volumetric fluid concentration ≥ 95 , correspond to negligible flow heights and therefore should be interpreted with precaution.This can be understood by referring back to the steady state analysis.When reducing the shear, the velocity increases until the flow reaches a steady state.This leads to the quasi constant behavior of the turbulent term with increasing fluid content.This relationship (Fig. 5) is similar to the behavior of the Coulomb stress (red dots) in Fig. 13.These results strengthen the hypothesis that when measuring the shear stress experimentally, it might be possible that a large part of the turbulent stress normally acting on the flow is absent, due to the smoothness of the force plate.In future, shear plates should be artificially roughened to represent real terrain.
The Ritigraben event of the 7 th August 2021 was of particular interest because drone flights were performed before and after the event.The pre-event flight was carried out two days before the event; the post-event flight was performed 2-days after the debris flow.Bed erosion and deposition could therefore be quantified.Figure 10b depicts the calculated erosion pattern using the two different mixture models.The general erosion pattern is quite similar for both formulations and can catch the general behavior of the field data, which is remarkable, regarding the fact that we never take into account the local structure of the bed channel, which is obviously of major importance in erosion processes.The entire eroded volume ( ≈ 10'000 m 3 ) is computed correctly using standard erosion parameters.The importance of erosion can be understood by comparing Fig. 12a with 12c and 12b with 12d.The difference between them are the erosion density (2000 kg/m 3 for Fig. 12a, c and 1800 kg/m 3 for Fig. 12b  d).It is clear that the shearing processes are strongly dependent of the erosion density.Therefore, the erosion does not only influence the mass balance of a debris flow, but its entire flow dynamics (Cuomo et al. 2014(Cuomo et al. , 2016;;Pudasaini and Fischer 2016;Pudasaini and Krautblatter 2021), through the erosion density, which can modify the solid/fluid ratio of a debris flow.For instance, the simulation performed with an erosion density of 1800 kg/m 3 reaches the runout zone more than one and half minutes earlier than the simulation executed with a density of 2000 kg/m 3 (on a complete flowing time of approximately 10 min).
Another aspect of the model is the ability to simulate dewatering and phase separations.Indeed, when the flow reaches the run-out zone, it decelerates and the energy coming from the shear work is not sufficient anymore to maintain the solid matrix in a dilated configuration.Therefore, the flow collapses and deposits (friction increases), as shown in Fig. 14a.During the collapse of the mixture layer, interstitial water is squeezed out and transformed into free fluid, which can escape from the solid matrix (Fig. 14b).If the free fluid is ejected downstream of the deposit (zone 2 on Fig. 14b), it continues to flow downward.If the fluid is squeezed out upward from the deposit (zone 1 on Fig. 14a), it is stopped by the deposit which acts like a dam.This can be of major importance when simulating debris flow events.Indeed, in Ritigraben and many other torrents, the risk of a lake formation upstream the solid deposit which can potentially release afterward and inundate the valley downstream has been widely studied.Therefore, the possibility to be able to simulate these lake formations with our model could help in the future hazard analysis.
To summarize, there exists a strong interaction between the streamwise evolution of a debris flow and the modeling of the frictional resistance.In debris flow hazard mitigation, the problem is often to choose a consistent set of parameters that govern the frictional resistance.In many case studies, however, a single set of fixed rheology parameters is not able to predict accurately the motion of a debris flow because the streamwise evolution of debris flow structure is not considered.Therefore, large modeling uncertainties remain.Fortunately, we have measurement data from the Illgraben test site which allows us to test the application of Voellmy-type mixture models within the framework of layered debris flow approaches.We have demonstrated that a single set of mixture parameters can represent the entire behavior of a debris flow, provided the model formulation tracks the fluid and solid components separately.Although we can envision the application of more complicated rheological models in the future, our results suggest that Voellmy-type models can be applied as a first step, reducing the uncertainties of practitioners considerably.

Fig. 1
Fig. 1 Debris flows inundated the village of Bondo, Switzerland after the Piz Cengalo rock/ice avalanche collapse in 2017

Fig. 3
Fig.3Sketch of two different debris flow configurations, containing overall of the same amount of solid and fluid mass.The left panel is in a dilated configuration, occurring while flowing.The right configuration is the reference configuration, called the co-volume configuration, typically occurring when the flow is at rest.The different heights and further physical parameters are described inMeyrat et al. (2021).We define three variables associated with the solid mass, h s , h 0 , and h sb .The height h s is the volume of the solid particles in the flow; h 0 represents the reference height of the non-dilated mass; we call it the co-volume (it is different from h s , because we consider that, even in the nondilated configuration, the void space is not zero), whereas h sb represents the dilated height of the solid mass, "Two-layer debris flow model with dilatancy" section.We refer to h 0 as the co-volume, in an analogy to Van der Waals work on non-ideal gasses with large molecules and cohesion(Visser 1989)

Fig. 4
Fig.4Map of the Illgraben test site.The catchment zone is colored in blue, while the channel is drawn in orange.The check dams are also shown by red lines.The blue star represents the starting point of the hydrogaph, used for the numerical simulation, and the green star is the location of the force plate.The original map can be found inMcArdell and Sartori (2020)

Fig. 5
Fig. 5 Experimentally measured shear-normal stress ratio (S/N) as a function of the measured fluid volumetric concentration for four different debris flow events.Each dot represents an experimental measurement.The color provides additional temporal information, and the blue represents the debris flow front and the yellow the tail

Fig. 6 a
Fig. 6 a (VM model) and b (PP model) Convergence of the simulation results to the analytical value of the Froude number(Eq.34, Fig. 4a and  Eq.35, Fig. 4b) in a steady state.An averaging is performed over the cells containing mass, which means that the entire flow is in a steady state (except the really front, see next).c (VM model) and d (PP model) According to Eq. (37), in a steady state, the shear stress (colored circles) does not depend on the rheological formulation and is equal to the gravitational stress (red lines).The color represents the time; the front is in blue, and as we move towards the tail, the color becomes yellow, see Fig.7.The very first part of the flow (blue markers) is not exactly in a steady state, due to numerical effects related to the shock-wave form of the front.The simulations are done on a smooth plane with constant angle, using a block release

Fig. 7 a
Fig. 7 a (VM model) and b (PP model) Convergence of the simulation results to the analytical value of the height of the debris flow (Eq.40) in a steady state.An average is done on each cell containing mass, which means that the entire flow is in a steady state, and it is assumed that we always have enough fluid to filled the void space in between the particles.c (VM model) and d (PP model) Flow height on a particular cell taken in the middle of the flow.The red line represent the numerical outputs, and the colored markers are obtained with Eq. (40).The simulations are done on a smooth plane with constant angle, using a block release

Fig. 8
Fig.8Map of the Ritigraben debris flow site.The original map can be found inBindereif (2022)

Fig. 11
Fig. 10 a Erosion pattern obtained with the VM model, b eroded volume per bins as a function of the position of the bins toward the channel.The red curve is obtained from field data, the blue one with the VM model and the green line from the PP Model.The erosion computed with both models is in good agreement with the measurements

Fig. 12
Fig. 12 Spatial evolution of the shearing coefficients and .Panel a (VM Model) and b (PP Model) show a decrease of from the front to the tail.The erosion density is 2000kg∕m 3 .Panel c (VM Model) and d (PP Model) are the same as a and b but for an erosion density of 1800kg∕m 3 .One can see the value of depends strongly on the erosion density.Therefore, the erosion does not only influence the mass ance of the flow but the entire flow dynamics.Panel e VM Model and f PP Model represent the increase in , which induces a decrease in the turbulent friction, from front to tail of the flow

Fig. 13 Fig. 14
Fig. 13Ratio of the total shear stress (green), the Coulomb shear stress (red) and the turbulent shear stress (blue) to the normal stress, as a function of the volumetric fluid content. is done for two different locations, see Fig.8(top and bottom) and for the two models (left and right).Location 1 (plots a and b) is in the flat part of the channel (  < 10 °) and then comparable to the Illgraben measurement, while location 2 (c and d) is in the steep part of the bottom channel ( ≈ 30°)

Fig. 15
Fig. 15 Overflow obtained with the numerical simulation.Despite the small size of this overflow, one can see that the mixture model predicts larger overflow than the PP model

Table 1
List of the parameter values used for the simulation *For the definitions of these values, refer toFrank et al. (20152017) Δx 1 m