A Multi‐Phase Mass Flow Model

Geomorphic mass flows are often complex in terms of material composition and its evolution in space and time. The simulation of those hazardous phenomena would strongly benefit from a multi‐phase model, considering the motion and—importantly—interaction of phases characterized by different physical aspects including densities, frictions, viscosities, fractions, and their mechanical responses. However, such a genuine multi‐phase model is still lacking. Here, we present a first‐ever, multi‐mechanical, multi‐phase mass flow model composed of three different phases: the coarse solid fraction, fine‐solid fraction, and viscous fluid. The coarse solid component, called solid, represents boulders, cobbles, gravels, or blocks of ice. Fine‐solid represents fine particles and sand, whereas water and very fine particles, including colloids, silt, and clay, constitute the viscous fluid component in the mixture. The involved materials display distinct mechanical responses and dynamic behaviors. Therefore, the solid, fine‐solid, and fluid phases are described by Coulomb‐plastic, shear‐ and pressure‐dependent plasticity‐dominated viscoplastic, and viscosity‐dominated viscoplastic rheologies. They are supposed to best represent those materials. The new model is flexible and addresses some long‐standing issues of multi‐phase mass flows on how to reliably describe the flow dynamics, runout, and deposition morphology of such type of phenomena. With reference to some benchmark simulations, the essence of the model and its applicability are discussed.

These materials behave differently due to different material parameters, phase interactions, mechanical responses, and basal boundary conditions. However, three-phase mass flow composition has always been neglected in simulations due to the high level of complexity. Here, we propose a novel multi-phase, multi-mechanical mass flow model, by extending the two-phase viscous fluid and Coulomb solid model (Pudasaini, 2012) to additionally combine it with the Coulomb-viscoplastic material. An extra phase adds several additional mechanical, interfacial, and dynamical problems. However, it addresses many long-standing challenges in multi-phase mass flows. The proposed coupled model considers three phases: viscous fluid, fine-solid, and solid. Determining the interfacial terms is challenging. However, Pudasaini (2012) provides the closures for interphase momentum transfers between fluid and solid. In contrast, the closures between the fluid and fine-solid, and the fine-solid and solid phases will be developed. Similarly, drags and virtual masses between fluid and viscoplastic, and Coulomb and viscoplastic materials will be constructed by exploring the basic idea from Pudasaini (2012).
Some basic simulations shed light on the essential mechanical, dynamical, and depositional characteristics of the new model. This paves the way for the simulation of real events, and the prediction of the flow dynamics of complex multi-phase mass flows.

Parameters and Variables
Let the solid, fine-solid, and fluid phases be denoted by the suffix s, fs, f, respectively. The three phases are designated by their distinct material properties: The fluid phase is governed by its true density f , viscosity f , and isotropic stress distribution; the fine-solid and solid phases are characterized by their true densities fs , s ; internal friction angles fs , s ; basal friction angles fs , s ; anisotropic stress distribution, K s (lateral earth pressure coefficient); and the viscosity of the fine-solid fs . These material properties and mechanical responses, and induced relative motions between the phases lead to three distinct mass and momentum balance equations for the solid, fine-solid, and the fluid phases, respectively. Let u s = ( , and s , fs , f denote the velocities with their components along the flow directions (x, y, z), and the volume fractions for the solid, fine-solid, and fluid constituents. Similarly, f is gravity; −T s , −T fs , −T f are the negative Cauchy stress tensors; fs , f are the extra stresses; and p fs and p f are the pressures. This implies T fs =−p fs I+ fs and T f =−p f I+ f ; and C DG and C vm constitute the interfacial force densities-the drags and the virtual mass forces. In what follows, the superscript-pair represents the considered phases, for example, C s, DG means the drag force exerted by fluid on solid. As we will see later, the final model equations are also characterized substantially by the distinct basal boundary conditions applied to the different phases in the mixture.

Effective Phase Fractions, Densities, Ambient Pressures, and Buoyancies
The total mixture density m composed of water (w); colloids, including clay, silt, and very fine particles (c); fine-solid (fs); and solid (s) can be written in different ways. Let f denotes the mixture composed of water where ∶= c + w , =̃∕ ,̃= c c + w w ; and the fraction of the fine-solid and fluid mixture is fsf ∶= fs + c + w . So the density of the mixture of fine-solid and fluid is given by s =̃s ∕ s ,̃s = s s + . Moreover, the hold-up identity, s + fs + c + w = 1, or s + fs + f = 1, is satisfied.
Correct understanding of the ambient pressures of different phases in the mixture is important. Since we consider the total debris mixture as three-phase material, fluid (water plus colloids), fine-solid, and solid, the ambient pressures on the fine-solid and the solid particles must be understood and described properly. The fluid molecule is surrounded by the fluid itself. Similarly, for the fine-solid particles, the ambient pressure is due to the fluid pressure. However, the situation for the solid particles is different and depends on competence, which is the suspension capacity of the solid-type material in the mixture of the fluid-type material.
Depending on the competence of the fine particles in the fluid, a certain fraction of the fine particles can be suspended and loaded in the fluid that increases the ambient pressure for the remaining fine particles larger than the diameter of the competent particle. This fraction of the fine-solid is the competence (parameter), and is denoted by fs . By definition fs ∈ (0, 1). Given the yield strength of the fluid matrix surrounding the fine-solid, fs can be estimated (Hampton, 1975(Hampton, , 1979. However, in many mass flows, fs appears to be closer to 1 than to 0. Figure 1 shows how the finer (or the weaker) material supports the coarser (or the stronger) material in the mixture. This is how the buoyancy force works in multi-phase mixtures.
With this physical understanding, the relevant density for the scaling of the background fluid for the fine-solid is f . Nevertheless, the relevant density of the background fluid for the solid particles must be considered as follows. First, we need to redefine fsf with respect to the competence parameter fs : fsf = fs fs + f . With this construction, the relevant density of the ambient mixture surrounding a solid particle is s =̃s ∕ s , but now,̃s = s s s + . This definition is consistent with the fact that as fs → 0, fsf fsf → f f , and as fs → 1, fsf fsf → fs fs + f f , which is the situation of the fully competent fine-solid. The latter results in highly buoyant or suspended (coarse) solids in debris flows (McArdell et al., 2007).
With respect to the competence fs , the new and legitimate definitions for the effective densities, and volume fractions of the solid, fine-solid, and the fluid phases are provided by the assignments: Similar assignments can also be constructed for the relevant phase viscosities (Mooney, 1951;Rutgers, 1962). In general, the left-hand quantities and variables for the phase densities and the volume fractions in (2) should be used while developing the model equations (the mass and momentum balances) for the mixture flows. Even when fs is a constant, the simulations obtained with this approach may produce different 10.1029/2019JF005204 results than the quantities on the right. Otherwise, the results might be substantially different. However, this adds lots of complexities. So, without loss of generality, we simplify the situation and consistently use the right-hand quantities and variables in the assignments (2). Nevertheless, it is important to keep track of fs while defining the initial phase volume fractions and densities of the three constituents in the mixture; and presenting and interpreting the final results in terms of the original phase densities ( s , s , ) and the phase volume fractions ( s , s ,

)
. Later, we will see that the above assignments suggest different legitimate scalings for the ambient pressures of the fluid, fine-solid, and solid phases in the mixture.
It is important to note that the definition (2) plays a crucial role in obtaining the proper interfacial force densities, including drags and virtual mass forces; hydraulic and basal pressure gradients; buoyancies and basal frictional forces, all (or most) of which are related to the buoyancy-reduced normal loads of the fine-solid and solid particles in the mixture. A further important point regarding the assignments in (2) is that a fraction fs fs of the fine-solid mass fs is transferred to the fluid, so this portion of the volume fraction of the fine-solid is reduced. This also consistently upgrades the effective density of the fluid. However, the density of the remaining portion of the fine-solid is unchanged. Strictly speaking, this is a mass transfer from fine-solid to fluid, and this can change in time and space, depending on the state of the fine-solid, that is, the changing competence of the fine-solid, and on how the yield strength of the fluid mixture evolves. So, a formally better way to deal with such a complex situation of interfacial mass transfer, and thus the associated momentum transfer, is required. One way to achieve this is by following and expanding the two-phase mechanical model for rock-ice avalanches (Pudasaini & Krautblatter, 2014) that deals with internal mass and momentum exchanges between solid-type material (rock and ice), and the fluid phases in the mixture. However, here we follow the simple redefinition (2) of the densities and volume fractions so that, at this stage, we do not need to explicitly deal with very complex mass and momentum exchanges between the phases.

A Three-Dimensional Three-Phase Mass Flow Model
Following Drew (1983) and Pudasaini (2012), we consider the phase-averaged balance equations for mass and momentum conservations and, for simplicity, assume that surface tension is negligible; the solid, fine-solid, and fluid components are incompressible; and no phase change occurs. For phase changes in (debris) mass flows see, for example, Pudasaini and Krautblatter (2014).
Extending the two-phase model (Pudasaini, 2012), the three-dimensional, three-phase mass flow model can be written as with s + fs + f = 1. Here, (3a)-(3c) are the mass balance equations for the solid, fine-solid, and fluid phases, respectively. And (4a)-(4c) are the corresponding momentum balance equations. It is important to note that in (4), the solid, fine-solid, and fluid stresses are accompanied by the respective volume fractions and that T s and fs , f are not coupled. Furthermore, s , fs , and f appear inside the differential operators. The drag terms, C DG , sum-up to zero and so do for the virtual mass terms, C vm . These distinct features of the field equations, consideration of distinct rheological models for the phases, and the way we develop and implement the generalized drag, virtual mass, enhanced non-Newtonian viscous stresses and the distinct phase-basal boundary conditions, ultimately lead to a set of genuinely new multi-phase, multi-mechanical model equations in a well-structured conservative form.

Constitutive Equations
Constitutive equations are required for the generalized interfacial force densities included through the drags C DG on the particulate phases (solids and fine-solids) and the virtual mass forces C vm due to the relative acceleration of the solid with respect to the fine-solid and fluid; and fine-solids with respects to the fluid; the solid stress tensor T s ; the fine-solid and fluid viscous stress tensors fs and f ; and the pressures p s , p fs and p f . Modeling these terms is challenging (Biesheuvel & Spoelstra, 1989;Drew, 1983;Ishii & Mishima, 1984;Jeffrey, 1973;Mokeyev, 1977;Pitman & Le, 2005;Pudasaini, 2012;Richardson & Zaki, 1954;van Wijngaarden, 1976;Zuber, 1964). Furthermore, the basal boundary conditions must be applied in accordance with the physics observed for the solid, fine-solid, and fluid phases (Domnik et al., 2013;von Boetticher et al., 2016).

Interfacial Force Density I. Generalized drag forces:
The two-phase generalized drag model proposed by Pudasaini (2012) and enhanced by Pudasaini (2018Pudasaini ( , 2020 is (pairwise) expanded to three phases as follows

PUDASAINI AND MERGILI 2925
The drags between the solid, fine-solid, and fluid phases are shown in Figure 2. In (5a),  s, ∈ (0, 1) combines the fluid-like,  s, = s ( s ) 3 Re s, p ∕180, and solid-like,  s, = M s, −1 , drag contributions between solid and fluid components in three-phase mass flows;  s, T is the terminal velocity of a particle falling through the fluid. Moreover, = 1 or 2 should be selected according to whether linear or quadratic drag coefficients are used, and M s, = M s, (Re s, p ) depends on the particle Reynolds number Re s, , 2005;Richardson & Zaki, 1954). Furthermore, d s is particle diameter, s = ∕ s is the fluid to solid density ratio, and s = ∕ s is the fluid to solid fraction ratio. The drag coefficients explicitly depend on the essential physical parameters, for example, the volume fractions of the solid and fluid, the solid and fluid densities, terminal velocity of solid particles, particle diameter, and fluid viscosity. The term  s, describes fluid flow through a solid skeleton (e.g., granular-rich debris flows), whereas  s, models solid particles moving through fluid. As  s, and  s, may have different degrees of importance, the generalized drag coefficient is expressed by a linear combination, with  s, , of these contributions. As limiting cases:  s, = 0 suitably models solid particles moving through a fluid, and  s, = 1 better suits for flows of fluids through dense packings of grains (Kytömma, 1991;Ouriemi et al., 2009;Pitman & Le, 2005;Pailha & Pouliquen, 2009). (5a) is called the smoothing function, where  s, = | s u s + u | is determined by the corresponding mixture mass flux per unit mixture density, typically  s, = 10 m/s (Pudasaini, 2018(Pudasaini, , 2020.  s, is a function of the solid volume fraction  s, = m s , where m is a positive number, close to 1. The emergence of   in (5a) is crucial for the broad structure of the enhanced drag that removes the singularity from the existing drag coefficients and offers a great opportunity to simulate a wide spectrum A similar structure can be developed for the virtual mass forces with the relative accelerations between the phases, that is, of flows. Without this smoothing term, the drag (in the previous models, for example, Ishii & Hibiki, 2011;Pailha & Pouliquen, 2009;Pitman & Le, 2005;Pudasaini, 2012) increases exponentially with solid volume fraction and thus produces singularity. Smoothly varying drags, for any fractions of solids and fine-solids, are necessary to properly and fully simulate multi-phase mass flows. Otherwise, either we need to restrict the flow simulations to the dilute regimes or the simulation results cannot be applied in solving the real flow problems due to the singular behavior of the drag terms in the dynamical equations. With this, (5a) is called the enhanced generalized drag in mixture mass flows. This fully describes the drag for any values of the solid volume fraction s . Similar discussions hold for the other drags C s, s DG and C s, DG with corresponding functions for .

II. Virtual mass forces:
In general, the solid particles may accelerate relative to the fine-solids and fluid, and the fine-solids may accelerate relative to the fluid. This results in the acceleration of parts of the ambient fluid and the fine-solids. These relative accelerations induce three additional forces, called the virtual mass forces (Biesheuvel & Spoelstra, 1989;Cook & Harlow, 1984;Dong, 1978;Lamb, 1952;Mokeyev, 1977;Zuber, 1964). Pairwise expanding the two-phase expression of Pudasaini (2012), we model the virtual masses (C vm ) for three phases by The physically based, fully analytical, and well-bounded two-phase virtual mass force coefficient (Pudasaini, 2019): can be extended for three-phase flows, where  vm is the virtual mass number, and and n are some numerical parameters. This model covers any distribution of the dispersive phase (dilute to dense distribution of the solid particles). As justified in Pudasaini (2019), the physically most relevant values for the parameters can be  0 vm = 10, = 0.12 and n = 1. The novel virtual mass force (7) is general and evolves automatically as a function of solid volume fraction. The other virtual mass force coefficients  s, s and  s, can be constructed from (7). Note that the model (7) is developed by considering the solid and fluid mass fluxes from the general two-phase mass flow model (Pudasaini, 2012). Then, by solving the virtual mass force induced solid and fluid mass fluxes, the fully analytical model (7) for the virtual mass force is obtained.

Stresses I. Solid stresses:
Solid stresses are assumed to satisfy the Mohr-Coulomb plasticity criterion (Pudasaini, 2012;Savage & Hutter, 1989). The lateral Cauchy stress components can be written in terms of the normal pressure (stress) and earth pressure coefficients.

II. Pressures:
The fundamental understanding of the ambient fluid pressures, explained in section 3.2, is important. The solid and fine-solid pressures p s , and p fs can be related to the reference fluid pressure p f . This will be made clear at section 4.1.

III. Fluid extra stresses:
By utilizing the two-phase models (Drew, 1983;Ishii, 1975;Pudasaini, 2012), the phase-averaged viscous fluid stresses for fine-solids and fluid phases in three-phase mixtures can be modeled by extending non-Newtonian fluid rheologies as follows: (8b) is the Newtonian-type viscous stress, where D f is the strain-rate tensor for fluid, and D is the corresponding deviatoric stress tensor, is the non-Newtonian viscous stress that depends on the solid-and fine-solid-volume fraction gradients which enhances the apparent viscous stress of fluid, f . Moreover, the superscript e in e and e s indicates the effective viscosities and will be explained at IV below. A similar description holds for the non-Newtonian stress contribution for the fine-solid in (8a). Unlike for fs , which contains only one non-Newtonian term associated with  s,s , contains two non-Newtonian terms associated with  ,s and  , s . This is so, because the kinetic energy of the fluid is enhanced by the motion of the solid and the fine-solid particles. The coefficient  ,s is called the mobility of the fluid at the interface with solid. Similar analyses hold for  s,s , and  , s . The importance of enhanced extra stress has been exclusively explained for the two-phase debris mixture in Pudasaini (2012). However, for the new three-phase model, it is important to realize that in (8b), the fluid shear stress is enhanced (or reduced) by the gradients of the solid and fine-solid concentrations.

IV. The fluid and fine-solid viscosities:
The fluid phase is modeled by the shear-rate-dependent Bingham-type fluid. So the effective fluid (kinematic) viscosity is given by where is the yield strength (usually a parameter) normalized with the relevant density and ||D f || is given by the second invariant (I I D ) of the strain-rate tensor for fluid: ||D || = √ I I D with, . We model the fine-solid phase as a pressure and shear-rate-dependent Coulomb-viscoplastic material (Domnik et al., 2013). So the effective kinematic viscosity for fine-solid is given by

10.1029/2019JF005204
where s is the pressure-and friction-dependent yield strength of the fine-solid (normalized with the relevant density), which can be modeled as s = sin s p s (Domnik et al., 2013;Khattri & Pudasaini, 2019b;von Boetticher et al., 2016von Boetticher et al., , 2017. p fs is the pressure of the fine-solid. Furthermore, ||D fs || is defined in terms of the second invariant of the strain-rate tensor for the fine-solid. In (9)-(10), e s and e are the kinematic viscosities of the Coulomb-viscoplastic and Bingham plastic fine-solid and fluid-type materials, respectively. These are fundamentally different rheologies. If ||D f || → 0, ||D fs || → 0, (9) and (10) , and , where the parameters r y are introduced for a smooth transition between yielded and unyielded regions (Domnik et al., 2013;Papanastasiou, 1987;von Boetticher et al., 2016).

Plastic Yielding
It is important to note that both the fine-solid and fluid phases yield plastically if the measures of the deviatoric stress tensors overcome the strengths of the materials (Domnik et al., 2013;Khattri & Pudasaini, 2019b;Prager & Drucker, 1952):

Free-Surface and Basal Conditions
Kinematic free-surface and bottom boundary conditions are satisfied for all the phases (Pitman & Le, 2005;Pudasaini, 2012). The free surfaces are traction-less. However, fundamentally different basal boundary conditions are required to physically reliably describe the interactions of the mixture material with the basal surface. Following the well-established laws, the Coulomb basal slip for the solid phase and the no-slip basal condition for the fluid phase can be applied. However, the fine-solids interactions with the base should better be modeled with the pressure-, rate-, and friction-dependent Coulomb-viscoplastic basal slip conditions (Domnik & Pudasaini, 2012;Domnik et al., 2013). The applied boundary conditions will be elaborated further at section 4.2.

System of Equations
In (3)- (4), there are 12 equations, but 14 unknowns. The system can be closed if the pressures of any of the two-phases can be related to the pressure of the third phase. This will be discussed at section 4.1, and at section 4.4C. Then, in principle, the system (3)-(4) can be solved numerically. There are strong couplings between the solids, fine-solids, and the fluid components both through the interfacial momentum transfers due to the viscous drags and virtual mass forces, the enhanced non-Newtonian viscous fluid stresses, and the novel descriptions of buoyancies and ambient pressures (section 3.2).

A Reduced Two-Dimensional Three-Phase Mass Flow Model
In general, the three-dimensional model equations (3)-(4) require huge computational efforts when applied to natural-scale geophysical mass flows. This is particularly so for multi-phase flows (Bout et al., 2018;Mergili et al., 2017;Mergili, Emmer, et al., 2018;Mergili, Frank, et al., 2018;von Boetticher et al., 2016von Boetticher et al., , 2017. However, to a large extent, flows under consideration can be long and wide compared to their depth. So depth-averaging can be utilized. For this purpose, the basal-and the free-surface of the flow are denoted by b = b(t, x, y) and s = s(t, x, y), and h = s − b is the flow depth (in the direction normal to the substrate surface). Although depth-averaging processes are great challenges in developing shallow flow models, particularly for multi-phase flows, such problems can be resolved, for example, by following the methods described in Pudasaini (2012). Depth-averaging involves a sequence of steps that rigorously transforms the three-dimensional equations into a relatively simple set of two-dimensional equations. This drastically reduces the computational cost and complexities. Here, we mention the main aspects needed to develop the reduced depth-averaged model equations. After that, we present the final set of equations for the flowing mass for the solids, fine-solids, and fluid components.

Scaling Analysis
The scaling analysis suggests to discard the physically insignificant terms (Pudasaini, 2012;Pudasaini & Hutter, 2007). With this, the reduced system is much simpler to solve numerically in an efficient way (Bout et al., 2018;Mergili et al., 2017;Mergili, Emmer, et al., 2018;Mergili, Frank, et al., 2018) than the full three-dimensional model. Another important aspect is the legitimate pressure scaling as suggested by the structure and composition of the ambient fluid pressure for each considered phase (section 3.  (2012), and the assignments (2), the buoyancy-reduced normal loads for the fine-solids and solids, respectively, scale with: s , and , where s = ∕ s , and s = ∕ s . Without the buoyancy effects, the factors s , and ( 1 − s ) would be 1. For notational brevity, the hats will now be dropped from the respective terms. Using the scaling equations as in Pudasaini (2012), the mass and momentum equations (3)-(4) can be written in non-dimensional form.
In what follows, g x , g y , g z are the components of gravitational acceleration. Here, the equations will be written in a locally inclined Cartesian-type topography (Pitman & Le, 2005;Pudasaini, 2012). However, the detailed basal topographic information can be included, for example, through the gradients of the basal topographies, and also in the Coulomb friction terms (Fischer et al., 2012;Mergili et al., 2017). Furthermore, with some algebra, the virtual mass terms can be combined with the inertial terms. Once the model equations are derived in non-dimensional form that suggests neglecting the terms of higher orders in small parameters (Pudasaini, 2012) (e.g., typical flow depth/extend) such terms are neglected. Afterwards, for ease of application, all the model equations are written back in their original physical dimensions.

Boundary Conditions
Following the single-phase and two-phase modeling approaches (Fernandez-Nieto et al., 2008;Iverson & Denlinger, 2001;Pitman & Le, 2005;Pudasaini, 2012;Pudasaini et al., 2005;Savage & Hutter, 1989), we assume that the solid-, fine-solid, and fluid phase constituents separately satisfy the kinematic free-surface and the bottom boundary conditions. However, different dynamic bottom boundary conditions should be applied for different constituents. Following the standard considerations (Pudasaini, 2012;Pudasaini & Hutter, 2007), Coulomb sliding for the solid and no-slip condition for the fluid are satisfied at the flow base. However, the fine-solid will be modeled in different ways.

I. No-slip:
This is a standard condition at the base for a usual fluid. During the depth-averaging process, two viscous terms in the momentum equations remain for the possible shearing of the u and v velocity components in the xz-and yz-planes perpendicular to the slip surface (xy-plane). Usually, these shearings are modeled by assuming parabolic velocity profiles and no-slip conditions at the bottom: where, [ * ] b indicates that " * " is evaluated at the bed, h is the total flow depth of the mixture (composed of fluid, fine-solid, and solid), and the parameters u s , v s are zero if there is no substantial shearing (plug-like flow), and about three for parabolic velocity profiles (Iverson & Denlinger, 2001;Pudasaini, 2012).
II. The pressure-and rate-dependent dynamic basal slip: This is a natural and more suitable basal slip condition for granular-type Coulomb-viscoplastic material (Domnik et al., 2013). The pressure-and rate-dependent dynamic basal slip for fine-solid is given by Here, C F is the Coulomb friction coefficient. Whether I or II is more appropriate still has to be evaluated.

Hydrostatic/Lithostatic Pressures; Uniform Velocity
The depth-averaged models essentially assume that the momentum transfer in the direction normal to the slope can be assumed negligible as compared to the momentum transfers in the slope-parallel directions. However, with our novel approach, we have essentially different pressures for the fluid, fine-solid, and solid phases (see section 4.4C). Furthermore, the phase velocities are largely assumed to be uniform through the flow depth, with some possible momentum corrections (Iverson & Denlinger, 2001;Pitman & Le, 2005;Pudasaini & Hutter, 2007).
With these procedures, the three-dimensional, three-phase mass flow equations (3)-(4) can be formally transferred into a simple, depth-averaged, two-dimensional, three-phase mass flow model.

Journal of Geophysical Research: Earth Surface
10.1029/2019JF005204

The Model Equations
The mass balance equations for the solid, fine-solid, and fluid phases are, respectively The x-directional depth-averaged momentum conservation equations for the solid, fine-solid, and fluid phases are, Due to symmetry, the y-directional momentum equations for the solid, fine-solid, and fluid phases can be written similarly here and in all the following considerations. This is achieved by formally utilizing the replacements: x ↔ y and u ↔ v, whenever necessary, both for variables and associated parameters. Below, we systematically present models for all the source terms and forces in the momentum equations.

A. The x-directional source terms are
These expressions are more general than those in Pudasaini (2012) for two-phase flow, for example, the fine-solid and fluid viscosities e s and e are now generally evolving functions and appear inside the gradients.
B. The drag coefficients are given by:

Journal of Geophysical Research: Earth Surface
with functions and parameters given at section 3.4.1I. These coefficients extend the generalized drag from two-phase (Pudasaini, 2012) to three-phase mass flows.
C. The effective basal pressures for the fluid, fine-solid, and the solid are This follows from the discussions at section 4.1 and formally expands the two-phase derivation in Pudasaini (2012) to three phases. Equation (17) shows distinct scalings for the fluid, fine-solid, and solid phases in the three-phase mixture flow. This also indicates that the fine-solid and solid pressures are reduced due to the respective buoyancies by the factors s and . Such pressure scalings are novel, legitimate, and influence the flow dynamics substantially, providing better solutions to mixture flow dynamics.

D. The x-directional hydraulic pressure coefficients for solid, fine-solid, and fluid are
where K is the earth pressure coefficient (Pudasaini & Hutter, 2007). The buoyancy-reduced normal load of the solid particles, , is due to the fluid composed of water and very fine particles and the fine-solids, and thus s is the corresponding mixture density normalized by the solid density. The structures in (18) basically emerge from (17). Similar statements hold for fine-solids.
E. The virtual mass induced mass and momentum enhancements for the solid phase due to fluid and the fine-solid are denoted by u vm s and uu vm s , uv vm s , and can be written as which are much more complex here for three phases as compared to two-phases (Drew, 1983;Drew & Lahey, 1987;Ishii & Mishima, 1984;Ishii & Zuber, 1979;Pudasaini, 2012;Zuber, 1964). As explained earlier, the virtual mass coefficients  can be constructed from (7). Similarly, the virtual mass force induced mass and momentum enhancements for the fine-solid and fluid phases can be obtained as follows: and respectively, where, s s = s ∕ s , s = s ∕ and s = s ∕ are the fraction ratios. By consistently replacing u by v in (19)-(21), we obtain the virtual mass induced mass and momentum enhancements in the y-direction.
F. The viscous stresses associated with e s and e in (15b)-(15c) are related to the Newtonian-type viscous stresses which are much more complex here for three phases as compared to the classical Newtonian stresses. This is due to their pressure-, rate-, yield strength-, and friction-dependency, and become more clearer below at H. G. The x-directional fluid-type basal shear stresses in the xz-plane are given (see section 4.2), either by the no-slip condition (for both the fluid and fine-solid): or by the no-slip condition for fluid and the Coulomb-slip condition for fine-solid: with the Coulomb friction coefficient C F u s = −u s ∕|u s | tan s , where fs is the basal friction angle for the fine-solid. The parameters u and u s in (22) and (23)

H. The effective fluid and fine-solid kinematic viscosities are given by
where and s are the corresponding yield stresses, r y are the parameters for regularization, and s = sin s p s , see section 3.4.2IV. In the viscosities (24), the depth-averaged norm of D f is obtained as A similar structure holds for the norm of D fs .

I. Flow and no-flow regions:
The yield criteria depend on the measures of the deviatoric stress tensors and the material strengths for both the fine-solid and fluid phases (Domnik et al., 2013;Prager & Drucker, 1952) and are given by . The yield criteria help to precisely distinguish the flow and no-flow regions and thus are very useful for exactly defining the final deposits. This new feature is a great advantage of the presented multi-phase model that was not available in previous mixture and two-phase mass flow models (Iverson & Denlinger, 2001;Pitman & Le, 2005;Pudasaini, 2012).
J. The x-directional enhanced non-Newtonian viscous stress contribution (denoted by nN) for fine-solid due to the non-uniform distribution of the solid particles in the fine-solid is given by

10.1029/2019JF005204
which results from (8a). Similarly, the enhanced non-Newtonian viscous stress contribution for fluid due to the non-uniform distribution of the fine-solid and solid particles in the fluid is given by which is derived from (8b). Due to the evolving fine-solid and fluid viscosities, and the additional fine-solid component in the mixture, the non-Newtonian viscous contributions (26)-(27) are more general than the non-Newtonian viscous contribution in existing models (Drew, 1983;Ishii, 1975;Ishii & Hibiki, 2006;Pudasaini, 2012).

A closed system of equations:
The model (13)-(14) constitutes a set of nine equations for mass and momentum balances (including the y-components) for three-phase mixture mass flows in nine unknowns, namely, the solid, fine-solid, and fluid phase velocities in the downslope ( u s , u s , u ) , and cross-slope

and the respective phase depths
. Note that h s + h fs + h f = h. The model is written in a well-structured form of partial differential equations and can be solved numerically once appropriate initial and boundary conditions are prescribed. Several important features of the model equations (13)-(14) are discussed in section 3.2 and supporting information.

Reduction to Existing Models
By setting the fine-solid fraction to zero ( fs → 0), the new mass flow model (13)-(14) reduces exactly to the general two-phase mass flow model of Pudasaini (2012). This offers a great advantage to directly expand the existing computational tools, such as r.avaflow (Mergili et al., 2017), which has been successfully applied to simulate natural events (Mergili, Emmer, et al., 2018;Mergili, Frank, et al., 2018) assuming two-phase mixtures of solid particles and viscous fluid (Pudasaini, 2012). Moreover, we note that the Pudasaini (2012) model extends and unifies the previously existing single-phase and two-phase models (Iverson & Denlinger, 2001;Pitman & Le, 2005;Pudasaini et al., 2005;Savage & Hutter, 1989). Hence, the multi-phase model constructed here presents the most extensive, adaptable, new generation of mass flow model.

Benchmark Numerical Simulations
In order to apprehend the rapidly changing behavior of the flow variables, the model equations such as (13)-(14) are solved in conservative variables with high-resolution numerical schemes (Nessyahu & Tadmor, 1990;Pudasaini & Hutter, 2007;Tai et al., 2002). This allows to extend the numerical strategy from two-phases (Pudasaini, 2012) to three phases.
The new multi-phase model is implemented in the geographic information system-based software tool r.avaflow 2.0, representing a further development of r.avaflow (Mergili et al., 2017). We demonstrate the model performance by using two benchmark tests imposed on regular, computer-generated topographies: a bell-shaped mound with simultaneous release of multiple flows characterized by fundamentally different initial material compositions (Test 1) and a valley head with landslides of different materials successively impacting a water body (Test 2). Whereas Test 1 illustrates the influence of material composition on the flow dynamics in an idealized geometry, Test 2 is representative of a large, extremely rapid rock avalanche and a large, but slower, viscous earth flow impacting a reservoir.

Bell-Shaped Mound with Multiple Flows: Test 1
First, we consider a circular, bell-shaped hill defined by a cosine function, with a height of 400 m and a maximum slope angle of approximately 14 • . This choice of the geometry allows to simulate most of the dynamics and depositions within a reasonably small area. Eight mass flow release areas of different material compositions are defined by distorted hemi-ellipsoidal shapes with radii of 100 m in horizontal, and 30 m in vertical direction. Mass flow A represents a completely dry flow of coarse material, B a fine-solid, highly viscous flow (e.g., earth flow), and C a flood of pure water. The mass flows D-H are composed of more than one type of material: D could represent an earth flow mixed with water from heavy rain or surface runoff, E a debris flow of coarse material and water, F a mixture flow of fine earth and coarse blocks (e.g., till material), G a flow of coarse material with low water content, and H a mixture flow with equal fractions of all three materials, which could represent a multi-phase debris flow. Each release area is characterized by its specific amount of solid, fine-solid, and fluid material (Figure 3a). The advantage of such a simulation set-up is that several fundamentally different initial materials (by compositions) can be considered, and flow dynamics, runout, and depositions can be analyzed in a simple and compact frame. This allows a direct comparison among the different simulations.
All flows are released at t = 0 s, whereas the simulation stops at t = 120 s. The distinct behavior of the individual flows is illustrated in Figure 3 in bird's eye view. We note that the slight asymmetry of those flows including viscous material is related to the numerical discretization. Whereas Figures 3b-3e show the flow heights at selected time intervals, Figure 3f illustrates the maximum flow height throughout the entire simulation. The influence of the material composition on the flow dynamics and mobility is clearly seen. The purely fluid flow is by far the most mobile. Results in panel (f) indicate that even 10% of fluid mixed with otherwise solid material results in a substantial increase in mobility, compared to pure solid. This effect becomes much stronger when considering a 50% solid-50% fluid mixture flow. In this latter case, partial separation of the phases occurs (Figure 3d), with the fluid phase rather moving to the front. The fine-solid, due to its high viscosity, is much less mobile than the fluid material. The flow of pure fine-solid is more mobile than the flow of pure solid material but hardly expands laterally. Also a mixture of 50% of fine-solid with 50% of fluid material shows a limited degree of mobility. This initial material composition first leads to an escape of the fluid at the front, and then on the lateral sides; still with limited mobility of the fine-solid that remains mainly along the central line. A mixture of one third of solid, fine-solid, and fluid each results in separation of the phases, with solid, and-particularly-fluid, at the front, and the fine-solid fraction of the flow remaining behind. It appears that viscosity may play a crucial role in mixing or preventing from mixing between different phases. These are crucially new understandings in multi-phase mass flow simulation.
The maximum flow heights through the entire simulations represented in Figure 3f reveal several mechanically important features of the flow dynamics. The landslide masses with different initial material compositions evolve differently both in runout, extent, and the actual material compositions. With respect to the evolving geometries, these can broadly be categorized into three classes (with some overlaps).

I. Limited motion and deformation (A, B, F):
The flows of pure solid, fine-solid, and a mixture of both show the least degree of mobility and deformation, because they are composed of frictional and/or highly viscous materials. These materials consume lots of energy against deformation and motion. The highly viscous fine-solid material deforms and propagates in downslope direction, but its extent remains largely constrained in cross-slope direction throughout the entire simulation. This type of behavior is commonly observed for highly viscous material movements. The mobility of these materials would increase when increasing the slope angle or decreasing friction and/or viscosity.

II. Expanding fans (C, D, E, G, H):
Importantly, fluid is involved in all of these initial material compositions and plays a considerable to dominant role. Due to the fluid only material in C, the fan is the widest, and the runout distance is the longest. This is quite understandable, because there is no friction, and viscosity is minimal. In D, E, and H, fluid and solid and/or fine-solid contribute equally. The flow moves more slowly than in C, and the lateral expansion is reduced. This is due to the frictional and viscous dissipation associated with the solid and fine-solid components. Importantly, as the fluid is mechanically weak as compared to the solid and fine-solid, the fluid escapes quickly both in the longitudinal and lateral directions. So the fluid material generates an envelope at the sides and the front of the solid and/or fine-solid dominated central surge of the debris mixture. Such behaviors are also observed in solid-fluid mixture mass flows (Kattel et al., 2018). In the simulation G with 90% solid and 10% fluid, a clear fan and frontal surge and lobe is developed. However, both the lateral extent, and typically the runout distance, are largely reduced as compared to C and E. As the initial landslide is largely solid dominated, it forms a triangular fan, but as the fluid contribution is minimal, the fan expansion and runout are reduced substantially, and so is the flow mobility. Such behaviors are observed in solid or solid dominated mass flows (Pudasaini et al., 2005). However, we note that the fine-solid in the fine-solid and fluid mixture moves slower than the fine-solid only. This can be explained because, first, the amount of fine-solid in the fine-solid and fluid mixture is reduced by 50%, compared to the fine-solid only flow. This lowers the hydraulic pressure of the fine-solid. Second, as the fluid is weaker than the fine-solid, it deforms and moves quickly to the front of the flow. So the fluid surge forms a frontal support for the fine-solid, resulting in the enhanced lateral expansion of the fine-solid, compared to the fine-solid only simulation. This leads to the slower frontal motion of the fine-solid in this simulation than in the fine-solid only flow.

III. Fingering and lobe formation (D, H):
The situation becomes more interesting when substantial fine-solid and fluid fractions are involved in the mixture-(viscous) fingerings and lobes are generated. Emergence of fingerings and lobes is new in multi-phase flow simulation. Clear and narrow fluid-dominated fingerings are observed in the mixture with fine-solid and fluid. Principally centrally down-hill motion of the fine-solid appears to force the mechanically weaker fluid to quickly escape laterally, which allows the fluid to rush down the slope forming longitudinally elongated multiple lobes on either sides of the flank. However, in the mixture of solid, fine-solid, and fluid, the lobes and fingerings are less pronounced along the down-hill direction, but they expand substantially across the slope. In this material composition, both the fine-solid and solid contribute to the formation and dynamics of fingering and lobes.
This analysis clearly reveals that the flow dynamics is strongly controlled by the material composition. These results are in line with the experimental evidence (de Haas et al., 2015) and also with the expected behavior of the flow materials with the given parameters. However, more detailed research is needed to better understand the dynamics of phase-separation, fingerings, and lobes.

Landslide-Reservoir Interactions: Test 2
Second, we consider a mountain landscape characterized by several elements: a U-shaped valley with an amphitheater-like valley head; a trapezoid-shaped dam with a height of 50 m in direction transversal to the valley; a reservoir impounded by the dam, with a maximum depth of 40 m and a volume of 16.8 million m 3 (fluid material); two landslide release areas of frictional solid and viscous fine-solid materials, in the shape of a hemisphere with a radius of 100 m and a volume of 2.1 million m 3 each (Landslide A and Landslide B). The drop height from the centers of the landslide release areas into the reservoir is approximately 220 m (Landslide A) and 170 m (Landslide B). The inclination of the valley slopes ranges between 30 • and 40 • near the landslide release areas, whereas it reduces to approximately 16 • near the reservoir level.
In contrast to the previous example described in section 5.1, initially, all three phases are separate. The arrangement of the elements is illustrated in Figure 4a, whereas the main difference in parameters now is that the viscosity of the fine-solid is higher than in Test 1. Landslide A is released immediately at the start of the simulation, whereas Landslide B is released after a delay of t = 240 s when the impact wave generated by Landslide A has largely settled. This setting allows, first, to analyze the dynamics of solid, fine-solid, and fluid separately, and second, the complex interactions among the phases, as solid and fine-solid impact the fluid reservoir. The simulation stops at t = 600 s. Figures 4b-4i illustrate the evolution of the mixture flow at different points in time, whereas Figure 4j displays the maximum flow height throughout the entire simulation. Figures S1a-S1f (in supplementary information) show vertical profiles and the corresponding fields of kinetic energy for the initial phase of the simulation (generation of the impact wave).
Landslide A reaches a maximum velocity of 88.1 m/s. It enters the reservoir before t = 10 s, causing a displacement wave that overtops the dam at about t = 50 s. Partial mixing between the solid particles in landslide and the water in the reservoir occurs. The flood wave at output hydrograph O1 ( Figure S1g) shows a strong peak of 83,000 m 3 /s at t = 80 s, and then quickly alleviates to a much lower level. The solid and fine-solid components of the discharge at O1 remain negligible, meaning that realistically almost the entire  Figure S1. landslide material remains in the reservoir. The very viscous, fine-solid Landslide B is released later. It represents a clay-rich, highly cohesive earth flow which slowly dives into the reservoir, without causing a major displacement wave. No clear signal of Landslide B is seen in the output hydrograph O1 ( Figure S1). Mixing of solid, fine-solid, and fluid material in the reservoir occurs to some extent, but the main impact area of Landslide B remains clearly discernible until the end of the simulation, showing that water is pushed away rather than incorporated in the fine-solid material. Even though the solid and fluid materials could be mixed, the fine remains relatively un-mixed with solid and fluid in this simulation. Such phenomena are realistic. This is mainly due to the highly viscous material considered here, which is only slowly penetrated by other types of material. In summary, also this benchmark test yields the expected results. Rigorous test runs against field data or laboratory experiments will be necessary for a quantitative evaluation, which, however, is out of scope of the present work. A further important aspect is the influence of the viscosity of the fine-solid component. In Figure 3, the deformation of fine-solid was relatively large. Whereas in Figure 4, one order of magnitude increase in the viscosity of the fine-solid results in substantially reduced and relatively slow deformation. This clearly indicates that proper choice of the parameters in real multi-phase flow simulations is crucial.

Discussion and Future Perspectives
We have presented a novel, multi-mechanical, and multi-phase model that better captures the dynamics of complex mass flows consisting of coarse particles, fine particles, and viscous fluid. This work should be seen as a first step in the direction of simulating real multi-phase mass flows. Several important aspects can further be considered when simulating multi-phase flows: • Interaction processes are the key drivers to cause transitioning from one dynamic mass flow regime to another, eventually leading to a cascading mass flow event (Mergili, Emmer, et al., 2018;Mergili, Frank, et al., 2018). Interactions with the basal sliding surface and with water bodies are most relevant for rapid gravity-driven mass flows. Even though advanced approaches to dynamically simulate entrainment and also deposition do exist in theory (Pudasaini & Fischer, 2016a), they still have to be implemented in r.avaflow (Mergili et al., 2017) in a way to also be suitable for real-world applications. Erosion of the bed materials enters the mass and momentum balance equations as additional mass and momentum production terms on the right-hand sides of (13a)-(13c) and (14a)-(14c), respectively. Thereby, nonlinear coupling and feedback mechanisms accounting for erosion, deposition, phase-separation, and particle sorting should be considered (Conway et al., 2010;de Haas et al., 2015;Johnson et al., 2012;Major & Iverson, 1999;McArdell et al., 2007;Pudasaini & Fischer, 2016b;Schneider et al., 2011;Steinkogler et al., 2015). • Phase transformations are important aspects. Theoretical concepts do exist also for transformations between phases (see, e.g., Pudasaini & Krautblatter, 2014). More research is necessary to use these concepts in the multi-phase model. The same is true for the frequency-dispersion and turbulence (Kafle et al., 2019). • The present multi-phase simulation tool (r.avaflow 2.0) still builds on the shallow flow assumption which is a reasonable approximation in many real-world cases, but not always (Domnik & Pudasaini, 2012;von Boetticher et al., 2016). The impact of a large landslide in a reservoir, for example, calls for a full 3D consideration as presented at section 3. The design of interfaces to use 2D wherever admissible and 3D wherever necessary remains a challenge for the future. The 1D-2D coupling approach in Domnik et al. (2013) can be extended to address such a problem. • The dynamic change of flow parameters depending on material composition needs to be investigated. Following Pudasaini and Krautblatter (2014), we hypothesize that key model parameters such as the internal and basal friction angles of the solid and fine-solid material are characterized by spatiotemporal variations depending on the material composition. Well-founded physically based rules to consider these variations may greatly improve our ability to more realistically simulate real-world mass flows. • A concept for block sliding would be useful to enable the simulation of process components dominated by basal sliding and no or little internal deformation, such as rock slides. In this way, the interaction of rock slides with reservoirs could be investigated, which can be different from the interaction of mass flows (such as rock avalanches) with reservoirs, shown in the benchmark example in section 5.2. The most prominent example of a rock slide impacting a reservoir is the 1963 Vajont event (Genevois & Tecca, 2013). • Further, the applicability of the proposed model to slow-moving, highly viscous flows such as earth flows, lava flows, or even the flow of glaciers has to be investigated. Including appropriate flow transformation mechanisms, this could, for example, allow for the simulation of complex volcanic process chains such as lava dome growth followed by dome collapse and the development of pyroclastic surges. Whereas, in principle, the model would largely serve for this type of process chains, the combination of different time scales might represent a challenge in terms of practical applications.
• Finally, field data and experimental data have to be utilized to systematically evaluate and enhance the computational model. Even though natural debris flows are monitored in a few debris flow torrents (Hürlimann et al., 2003;Marchi & Tecca, 2013;McCoy et al., 2010), these initiatives face many difficulties because of the short-lasting, infrequent, and destructive flows. There is a great need to address these problems in order to enable a more reliable evaluation of the proposed multi-phase flow model.

Conclusions
The new multi-phase, multi-mechanical model presented in this work represents a major step forward with regard to our ability to simulate cascading landslides, particularly in high-mountain areas where various materials including rock, debris, ice, mud, and water may interact dynamically. The model opens up opportunities to more realistically back-calculate documented events in order to build a basis for predictive simulations of possible future events, so as to establish a better foundation for the design of risk reduction strategies. Thereby, the possibility to consider and combine more than two types of material (phases) with their specific characteristics in an arbitrary way greatly increases our flexibility. Building on the Pudasaini (2012) two-phase flow model, the proposed multi-phase model considers the conservation of mass and momentum of each phase. Taylor-made models of the material properties of each phase, interfacial dynamic interactions, and the forces acting between the flow and the basal surface allow for a realistic representation of complex flow-dominated landslides. Whereas, in principle, the model would allow to consider an arbitrary number of phases, we only deal with three phases in the present work as the complexity of interactions would increase in a nonlinear way with each additional phase added.
In order to demonstrate the scope of the proposed multi-phase model, we have presented two benchmark simulations. The results of these simulations clearly confirm the physical plausibility of the model and disclose patterns that were never before reproduced with simulation models (e.g., fingering of mixed flows). The spatiotemporal patterns produced by the complex interactions between the phases are in line with the expected patterns, indicating a reasonable degree of physical robustness of the proposed model. However, as a first step, these tests were performed on synthetic case studies only. Future work has to confirm the practical applicability with tests against experimental data and real-world case studies, out of scope here.
In the future, the model could also profit from extensions such as entrainment, phase-separation, flow transformation, and full three-dimensional implementation.