Model description paper 06 Apr 2021
Model description paper  06 Apr 2021
Towards a model for structured mass movements: the OpenLISEM hazard model 2.0a
 ^{1}Faculty of GeoInformation Science and Earth Observation, University of Twente, Enschede, the Netherlands
 ^{2}State Key Laboratory of Geohazard Prevention and GeoEnvironment Protection, Chengdu University of Technology, Chengdu, China
 ^{3}Institute of Mountain Hazards and Environment, Chinese Academy of Sciences, Chengdu, China
 ^{1}Faculty of GeoInformation Science and Earth Observation, University of Twente, Enschede, the Netherlands
 ^{2}State Key Laboratory of Geohazard Prevention and GeoEnvironment Protection, Chengdu University of Technology, Chengdu, China
 ^{3}Institute of Mountain Hazards and Environment, Chinese Academy of Sciences, Chengdu, China
Correspondence: Bastian van den Bout (b.vandenbout@utwente.nl)
Hide author detailsCorrespondence: Bastian van den Bout (b.vandenbout@utwente.nl)
Mass movements such as debris flows and landslides differ in behaviour due to their material properties and internal forces. Models employ generalized multiphase flow equations to adaptively describe these complex flow types. Such models commonly assume unstructured and fragmented flow, where internal cohesive strength is insignificant. In this work, existing work on twophase mass movement equations are extended to include a full stress–strain relationship that allows for runout of (semi)structured fluid–solid masses. The work provides both the threedimensional equations and depthaveraged simplifications. The equations are implemented in a hybrid material point method (MPM), which allows for efficient simulation of stress–strain relationships on discrete smooth particles. Using this framework, the developed model is compared to several flume experiments of clay blocks impacting fixed obstacles. Here, both final deposit patterns and fractures compare well to simulations. Additionally, numerical tests are performed to showcase the range of dynamical behaviour produced by the model. Important processes such as fracturing, fragmentation and fluid release are captured by the model. While this provides an important step towards complete mass movement models, several new opportunities arise, such as application to fragmenting mass movements and block slides.
The Earth's rock cycle involves sudden release and gravitydriven transport of sloping materials. These mass movements have a significant global impact in financial damage and casualties (Nadim et al., 2006; Kjekstad and Highland, 2009). Understanding the physical principles at work at their initiation and runout phase allows for better mitigation and adaptation to the hazard they induce (Corominas et al., 2014). Many varieties of gravitationally driven mass movements have been categorized according to their material physical parameters and type of movement. Examples are slides, flows and falls, consisting of soil, rocks or debris (Varnes, 1987). Major factors in determining the dynamics of mass movement runout are the composition of the moving material and the internal and external forces during initiation and runout.
Within the cluster of existing mass movement processes, a distinction can be made based on the cohesive of the mass during movement. Postrelease, a sloping mass might be unstructured, such as mud flows, where grain–grain cohesive strength is absent. Alternatively, the mass can be fragmentative, such as strongly deforming landslides or fragmenting of rock avalanches upon particle impact. Finally, there are coherent or structured mass movements, such as can be the case for block slides where internal cohesive strength can resist deformation for some period (Varnes, 1987). The general importance of the initially structured nature of mass movement material is observed for a variety of reasons. First, block slides are an important subset of mass movement types (Hayir, 2003; Beutner and Gerbi, 2005; Reiche, 1937; Tang et al., 2009). This type of mass movement features some cohesive structure to the dynamic material in the movement phase. Secondly, during movement, the spatial gradients in local acceleration induce strain and stress that results in fracturing. This process, often called fragmentation in relation to structured mass movements, can be of crucial importance for mass movement dynamics (Davies and McSaveney, 2009; Delaney and Evans, 2014; Dufresne et al., 2018; Corominas et al. 2019). The lubricating effect from basal fragmentation can enhance velocities and runout distance significantly (Davies et al., 2006; Tang et al., 2009). Otherwise, fragmentation generally influences the rheology of the movement by altering grain–grain interactions (Zhou et al., 2005). The importance of structured material dynamics is further indicated by engineering studies on rock behaviour and fracture models (Kaklauskas and Ghaboussi, 2001; Ngekpe et al., 2016; Dhanmeher, 2017).
Dynamics of geophysical flows are complex and depend on a variety of forces due to their multiphase interactions (Hutter et al., 1994). Physically based models attempt to describe the internal and external forces of all of these mass movements in a generalized form (David and Richard, 2011; Pudasaini, 2012; Iverson and George, 2014). This allows these models to be applied to a wide variety of cases, while improving predictive range. A variety of both one, two and threedimensional sets of equations exist to describe the advection and forces that determine the dynamics of geophysical flows.
For unstructured (fully fragmented) mass movements, a variety of models exist relating to Mohr–Coulomb mixture theory. Such mass movements are described as nonNewtonian granular flows with dominant particle–particle interactions, assuming perfect mixing and continuous movement. Examples are debris flows and mudslides, while block slides and rockslides do not fit these criteria. Within these models, the Mohr–Coulomb failure surface is described with zero cohesive strength and only an internal friction angle (Pitman and Le, 2005). Examples that simulated a single mixed material exist throughout the literature (e.g. Rickenmann et al., 2006; Julien and O'Brien, 1997; Luna et al., 2012; van Asch et al., 2014). Twophase models describe solids, fluids, and their interactions; provide additional detail; and generalize in important ways (Sheridan et al., 2005; Pitman and Le, 2005; Pudasaini, 2012; Iverson and George, 2014; Mergili et al., 2017). Recently, a threephase model has been developed that includes the interactions between small and larger solid phases (Pudasaini and Mergili, 2019). Typically, implemented forces include gravitational forces and, depending on the rheology of the equations, drag forces, viscous internal forces and a plasticity criterion. The assumption of zero cohesion in the Mohr–Coulomb material is invalid for any structured mass movement. Some models do implement a nonNewtonian viscous yield stress based on depthaveraged strain estimations (Boetticher et al., 2017; Fornes et al., 2017; Pudasaini and Mergili, 2019). However, this approach lacks the process of fragmentation and internal failure.
For structured mass movements, limited approaches are available. These movements feature some discrete interparticle connectivity that allows the moving material to maintain a elastoplastic structure. Examples here are block slides, rockslides and some landslides (Aaron and Hungr, 2016). These materials can be described by a Mohr–Coulomb material with cohesive strength (Spencer, 2012). Aaron and Hungr developed a model for simulation of initially coherent rock avalanches (Aaron and Hungr, 2016) as part of DAN3D Flex. Within their approach, a rigidblock momentum analysis is used to simulate initial movement of the block. After a specified time, the block is assumed to fragment, and a granular flow model using a Voellmytype rheology is used for further runout. Their approach thus lacks a physical basis for the fragmenting behaviour. Additionally, by dissecting the runout process in two stages (discrete block and granular flow), benefits of holistic twophase generalized runout models are lost. Finally, Greco et al. (2019) presented a runout model for cohesive granular matrix. Their approach similarly lacks a description of the fragmentation process. Thus, within current mass movement models, there might be improvements available from assuming nonfragmented movement. This would allow for description of structured mass movement dynamics.
In this paper, a generalized mass movement model is developed to describe runout of an arbitrarily structured twophase Mohr–Coulomb material. The model extents on recent innovations in generalized models for Mohr–Coulomb mixture flow (Pudasaini, 2012; Pudasaini and Mergili, 2019). Section 2 provides the derivation of the extensive set of equations that describe structured mass movements in a generalized manner. Section 3 validates the developed model by comparison with results from controlled flume runout experiments. Additionally, Sect. 3 shows numerical simulation examples that highlight fragmentation behaviour and its influence on runout dynamics. Finally, in Sect. 4 a discussion on the potential usage of the presented model is provided, together with reflection on important opportunities of improvement.
2.1 Structured mass movements
Gravitational mass flows are triggered when the local driving forces within an often steep section of a slope exceed a critical threshold. The instability of such materials is generally understood to take place along a failure plane (Zhang et al., 2011; Stead and Wolter, 2015). Along this plane, forces exerted due to gravity and possible seismic accelerations can act as a driving force towards the downslope direction, while a normal force on the terrain induces a resisting force (Xie et al., 2006). When internal stress exceeds specified criteria, commonly described using Mohr–Coulomb theory, fracturing occurs, and the material becomes dynamic. Observations indicate material can initially fracture predominantly at the failure plane (Tang et al., 2009; Davies et al., 2006). Full finiteelement modelling of stability confirms no fragmentation occurs at initiation, and runout can start as a structured mass (Matsui and San, 1992; Griffiths and Lane, 1999).
Once movement is initiated, the material is accelerated. Due to spatially nonhomogeneous acceleration, either caused by a nonhomogeneous terrain slope or impact with obstacles, internal stress can build within the moving mass. The stress state can reach a point outside the yield surface, after which some form of deformation occurs (e.g. plastic, brittle, ductile) (Loehnert et al., 2008). In the case of rock or soil material, elastic or plastic deformation is limited, and fracturing occurs at relatively low strain values (Kaklauskas and Ghaboussi, 2001; Dhanmeher, 2017). Rocks and soil additionally show predominantly brittle fracturing, where strain increments at maximum stress are small (Bieniawaski, 1967; Price, 2016; Hušek et al., 2016). For soil matrices, cohesive bonds between grains originate from causes such as cementing, frictional contacts and root networks (Cohen et al., 2009). Thus, the material breaks along either the grain–grain bonds or on the molecular level. In practice, this process of fragmentation has frequently been both observed and studied. Cracking models for solids use stress–strain descriptions of continuum mechanics (Menin et al., 2009; Ngekpe et al., 2016). Fracture models frequently use smooth particle hydrodynamics (SPH) since a Lagrangian, meshfree solution benefits possible fracturing behaviour (Maurel and Combescure, 2008; Xu et al., 2010; Osorno and Steeb, 2017). Within the model developed below, knowledge from fracturesimulating continuum mechanical models is combined with finiteelement fluid dynamic models.
The Mohr–Coulomb mixture models on which the developed model is based can be found in Pitman and Le (2005), Pudasaini (2012), Iverson and George, (2014), and Pudasaini and Mergili (2019). While these are commonly named debris flow models, their validity extends beyond this typical category of mass movement. This is both apparent from model applications (Mergili et al., 2018) and theoretical considerations (Pudasaini, 2012). A major cause for the usage of debris flow as a term here is the assumption of unstructured flow, which we are aiming to solve in this work.
2.2 Model description
We define two phases within the flow, solids and fluids, indicated by s and f, respectively. A specified fraction of solids within this mixture is at any point part of a structured matrix. This structured solid phase, indicated by sc, envelops and confines a fraction of the fluids in the mixture, indicated by fc. The solids and fluids are defined in terms of the physical properties such as densities (ρ_{f},ρ_{s}) and volume fractions (${\mathit{\alpha}}_{\mathrm{f}}=\frac{f}{f+s},{\mathit{\alpha}}_{\mathrm{s}}=\frac{s}{f+s})$. The confined fractions of their respective phases are indicated as f_{sc} and f_{fc} for the volume fraction of confined solids and fluids, respectively (Eqs. 1, 2 and 3).
For the solids, internal friction angle (ϕ_{s}) and effective (volumeaveraged) material size (d_{s}) are additionally defined. We also define ${\mathit{\alpha}}_{\mathrm{c}}={\mathit{\alpha}}_{\mathrm{s}}+{f}_{\mathrm{fc}}{\mathit{\alpha}}_{\mathrm{f}}$ and ${\mathit{\alpha}}_{\mathrm{u}}=(\mathrm{1}{f}_{\mathrm{fc}}){\mathit{\alpha}}_{\mathrm{f}}$ to indicate the solids with confinedfluid and freefluid phases, respectively. These phases have a volumeaveraged density ρ_{sc},ρ_{f}. We let the velocities of the unconfined fluid phase (${\mathit{\alpha}}_{\mathrm{u}}=(\mathrm{1}{f}_{\mathrm{fc}}){\mathit{\alpha}}_{\mathrm{f}}$) be defined as ${u}_{\mathrm{u}}=\left({u}_{\mathrm{u}},{v}_{\mathrm{u}}\right)$. We assume velocities of the confined phases (${\mathit{\alpha}}_{\mathrm{c}}={\mathit{\alpha}}_{\mathrm{s}}+{f}_{\mathrm{fc}}{\mathit{\alpha}}_{\mathrm{f}}$) can validly be assumed to be identical to the velocities of the solid phase, ${u}_{\mathrm{c}}=\left({u}_{\mathrm{c}},{v}_{\mathrm{c}}\right)={u}_{\mathrm{s}}=\left({u}_{\mathrm{s}},{v}_{\mathrm{s}}\right)$. A schematic depiction of the represented phases is shown in Fig. 1.
A major assumption is made here concerning the velocities of both the confined and free solids (sc and s), that have a shared averaged velocity (u_{s}). We deliberately limit the flow description to two phases, opposed to the innovative work of Pudasaini and Mergili (2019) that develop a multimechanical threephase model. This choice is motivated by considerations of applicability (reducing the number of required parameters), the infancy of threephase flow descriptions and finally the general observations of the validity of this assumption (Ishii, 1975; Ishii and Zuber, 1979; Drew, 1983; Jakob et al., 2005; George and Iverson, 2016).
The movement of the flow is described initially by means of mass and momentum conservation (Eqs. 4 and 5).
Here we add the individual forces based on the work of Pudasaini and Hutter (2003), Pitman and Le (2005), Pudasaini (2012), Pudasaini and Fischer (2016), and Pudasaini and Mergili (2019) (Eqs. 6 and 7).
Here f is the body force (a part of which is gravity), M_{DG} is the drag force, M_{vm} is the virtual mass force, and T_{c}andT_{u} are the stress tensors for solids with confined fluids and unconfined phases, respectively. The virtual mass force described the additional work required by differential acceleration of the phases. The drag force describes the drag along the interfacial boundary of fluids and solids. The body force describes external forces such as gravitational acceleration and boundary forces. Finally, the stress tensors describe the internal forces arising from strain and viscous processes. Both the confined and unconfined phases in the mixture are subject to stress tensors (T_{c}, and T_{u}), for which the gradient acts as a momentum source. Additionally, we follow Pudasaini (2012) and add a buoyancy force (p_{c}∇α_{c} and p_{f}∇α_{u}).
2.2.1 Stress tensors describing internal structure
Based on known twophase mixture theory, the internal and external forces acting on the moving material are now set up. This results in several unknowns, such as the stress tensors (T_{c} and T_{u}, described by the constitutive equation), the body force (f), the drag force (M_{DG}) and the virtual mass force (M_{vm}). This section will first describe the derivation of the stress tensors. These describe the internal stress and viscous effects. To describe structured movements, these require a full stress–strain relationship, which is not present in earlier generalized mass movements models. Afterwards, existing derivation of the body, drag and virtual mass force are altered to conform the new constitutive equation.
Our first step in defining the momentum source terms in Eqs. (6) and (7) is the definition of the fluid and solid stress tensors. Current models typically follow the assumptions made by Pitman and Le (2005), who indicate that “the proportionality (and alignment) of the tangential and normal forces that is imposed as a basal boundary condition is assumed to hold throughout the thin flowing layer of material … following Rankine (1857), an earth–pressure relation is assumed for diagonal stress components.” Here, the earth–pressure relationship is a vertically averaged analytical solution for lateral forces exerted by an earth wall. Thus, unstructured columns of moving mixtures are assumed. Here, we aim to use the full Mohr–Coulomb relations. Describing the internal stress of soil and rock matrices is commonly achieved by elastoplastic simulations of the material's stress–strain relationship. Since we aim to model a full stress description, the stress tensor is equal to the elastoplastic stress tensor (Eq. 8).
Here σ is the elastoplastic stress tensor for solids. The stress can be divided into the deviatoric and nondeviatoric contributions (Eq. 9). The nondeviatoric part acts normally on any plane element (in the manner in which a hydrostatic pressure acts equally in all directions). Note that we switch to tensor notation when describing the stress–strain relationship. Thus, superscripts (α and β) represent the indices of basis vectors (x, y, or z axis in Euclidian space), and obtain tensor elements. Additionally, the Einstein convention is followed (automatic summation of nondefined repeated indices in a single term).
Here s is the deviatoric stress tensor and ${\mathit{\delta}}^{\mathit{\alpha}\mathit{\beta}}=[\mathit{\alpha}=\mathit{\beta}]$ is the Kronecker delta.
Here, we define the elastoplastic stress (σ) based on a generalized Hooketype law in tensor notation (Eqs. 10 and 11), where plastic strain occurs when the stress state reaches the yield criterion (Spencer, 2012; Necas, and Hlavácek, 2017; Bui et al., 2008).
Here ${\dot{\mathit{\u03f5}}}_{\mathrm{elastic}}$ is the elastic strain tensor, ${\dot{\mathit{\u03f5}}}_{\mathrm{plastic}}$ is the plastic strain tensor, ${\dot{\mathit{\sigma}}}^{\mathrm{m}}$ is the mean stress rate tensor, ν is Poisson's ratio, E is the elastic Young's modulus, G is the shear modulus, $\dot{s}$ is the deviatoric shear stress rate tensor, $\dot{\mathit{\lambda}}$ is the plastic multiplier rate and g is the plastic potential function. Additionally, the strain rate is defined from velocity gradients as Eq. (12).
By solving Eqs. (9), (10) and (11) for $\dot{\mathit{\sigma}}$, a stress–strain relationship can be obtained (Eq. 13) (Bui et al., 2008).
Here $\dot{e}$ is the deviatoric strain rate (${\dot{e}}^{\mathit{\alpha}\mathit{\beta}}={\dot{\mathit{\u03f5}}}^{\mathit{\gamma}\mathit{\gamma}}\frac{\mathrm{1}}{\mathrm{3}}{\dot{\mathit{\u03f5}}}^{\mathit{\alpha}\mathit{\beta}}{\mathit{\delta}}^{\mathit{\alpha}\mathit{\beta}}$), ψ is the dilatancy angle, and K is the elastic bulk modulus and the material parameters defined from E and ν (Eq. 14).
Fracturing or failure occurs when the stress state reaches the yield surface, after which plastic deformation occurs. The rate of change of the plastic multiplier specifies the magnitude of plastic loading and must ensure a new stress state conforms to the conditions of the yield criterion. By means of substituting Eq. (13) in the consistency condition ($\frac{\partial f}{\partial {\mathit{\sigma}}^{\mathit{\alpha}\mathit{\beta}}}d{\mathit{\sigma}}^{\mathit{\alpha}\mathit{\beta}}=\mathrm{0}$), the plastic multiplier rate can be defined (Eq. 15) (Bui et al., 2008).
The yield criteria specifies a surface in the stress state space that the stress state can not pass and at which plastic deformation occurs. A variety of yield criteria exist, such as Mohr–Coulomb, Von Mises, Drucker–Prager and Tresca (Spencer, 2012). Here, we employ the Drucker–Prager model fitted to Mohr–Coulomb material parameters for its accuracy in simulating rock and soil behaviour and numerical stability (Spencer, 2012; Bui et al., 2008) (Eqs. 16 and 17).
Here I_{1} and J_{2} are tensor invariants (Eqs. 18 and 19).
Here the Mohr–Coulomb material parameters are used to estimate the Drucker–Prager parameters (Eq. 20).
Using the definitions of the yield surface and stress–strain relationship, combining Eqs. (13), (15), (16) and (17), the relationship for the stress rate can be obtained (Eqs. 21 and 22).
In order to allow for the description of large deformation, the Jaumann stress rate can be used, which is a stressrate that is independent from a frame of reference (Eq. 23).
Here $\dot{\mathit{\omega}}$ is the spin rate tensor, as defined by Eq. (24).
Due to the strain within the confined material, the density of the confined solid phase (ρ_{c}) evolves dynamically according to Eq. (25).
Here ϵ_{v} is the total volume strain, ${\dot{\mathit{\u03f5}}}_{\mathrm{v}}\approx {\mathit{\u03f5}}_{\mathrm{1}}+{\mathit{\u03f5}}_{\mathrm{2}}+{\mathit{\u03f5}}_{\mathrm{3}}$, ϵ_{i} is one of the principal components of the strain tensor. Since we aim to simulate brittle materials, where volume strain remains relatively low, we assume that changes in density are small compared to the original density of the material ($\frac{\partial {\mathit{\rho}}_{\mathrm{c}}}{\partial t}\ll {\mathit{\rho}}_{\mathrm{c}}$).
2.2.2 Fragmentation
Brittle fracturing is a processes commonly understood to take place once a material internal stress has reached the yield surface, and plastic deformation has been sufficient to pass the ultimate strength point (Maurel and Cumescure, 2008; Hušek et al., 2016). A variety of approaches to fracturing exist within the literature (Ma et al., 2014; Osomo and Steeb, 2017). Finite element method (FEM) models use strainbased approaches (Loehnert et al., 2008). For SPH implementations, as will be presented in this work, distancebased approaches have provided good results (Maurel and Cumbescure, 2008). Other works have used strainbased fracture criteria (Xu et al., 2010) . Additionally, dynamic degradation of strength parameters have been implemented (Grady and Kipp, 1980; De Vuyst and Vignjevic, 2013; Williams, 2019). Comparisons with observed fracture behaviour has indicated the predictive value of these schemes (Xu et al., 2010; Hušek et al., 2016). We combine the various approaches to best fit the dynamical multiphase mass movement model that is developed. Following Grady and Kipp (1980), we simulate a degradation of strength parameters. Our material consists of a soil and rock matrix. We assume fracturing occurs along the intergranular or interrock contacts and bonds (see also Cohen et al., 2009). Thus, cohesive strength is lost for any fractured contacts. We simulate degradation of cohesive strength according to a volume strain criteria. When the stress state lies on the yield surface (the set of critical stress states within the sixdimensional stressspace), during plastic deformation, strain is assumed to contribute to the fracturing of the granular or rock material. A critical volume strain is taken as a material property, and the breaking of cohesive bonds occurs based on the relative volume strain. Following Grady and Kipp (1980) and De Vuyst and Vignjevic (2013), we assume that the degradation behaviour of the strength parameter is distributed according to a probability density distribution. Commonly, a Weibull distribution is used (Williams, 2019). Here, for simplicity we use a uniform distribution of cohesive strength between 0 and 2c_{0}, although any other distribution can be substituted. Thus, the expression governing cohesive strength becomes Eq. (26).
Here c_{0} is the initial cohesive strength of the material, ϵ_{v0} is the initial volume, $\left(\frac{{\mathit{\u03f5}}_{\mathrm{v}}}{{\mathit{\u03f5}}_{\mathrm{v}\mathrm{0}}}\right)$ is the fractional volumetric strain rate and ϵ_{c} is the critical fractional volume strain for fracturing.
2.2.3 Water partitioning
During the movement of the mixed mass, the solids can thus be present as a structured matrix. Within such a matrix, a fluid volume can be contained (e.g. as originating from a groundwater content in the original landslide material). These fluids are typically described as groundwater flow following Darcy's law, which poses a linear relationship between pressure gradients and flow velocity through a soil matrix. In our case, we assumed the relative velocity of water flow within the granular solid matrix as very small compared to both solid velocities and the velocities of the free fluids. As an initial condition of the material, some fraction of the water is contained within the soil matrix (f_{fc}). Additionally, for loss of cohesive structure within the solid phase, we transfer the related fraction of fluids contained within that solid structure to the free fluids.
Beyond changes in f_{fc} through fracturing of structured solid materials, no dynamics are simulated for influx or outflux of fluids from the solid matrix. The initial volume fraction of fluids in the solid matrix defined by ff_{fc} and sf_{sc} remains constant throughout the simulation. The validity of this assumption can be based on the slow typical fluid velocities in a solid matrix relative to fragmented mixed fluid–solid flow velocities (Kern, 1995; Saxton and Rawls, 2006). While the addition of evolving saturation would extend the validity of the model, it would require implementation of pretransfer functions for evolving material properties, which is beyond the scope of this work. An important note on the points made above is the manner in which fluids are repartitioned after fragmentation. All fluids in fragmented solids are released, but this does not equate to free movement of the fluids or a disconnection from the solids that confined them. Instead, the equations continue to connect the solids and fluids through drag, viscous and virtual mass forces. Finally, the density of the fragmented solids is assumed to be the initially set solid density. Any straininduced density changes are assumed small relative to the initial solid density ($\frac{{\mathit{\rho}}_{\mathrm{c}}}{{\mathit{\rho}}_{\mathrm{s}}}\ll \mathrm{1}$).
2.2.4 Fluid stresses
The fluid stress tensor is determined by the pressure and the viscous terms (Eqs. 29 and 30). Confined solids are assumed to be saturated and constant during the flow.
Here I is the identity tensor, τ_{f} is the viscous stress tensor for fluids , p_{f} is the fluid pressure, η_{f} is the dynamic viscosity of the fluids and 𝒜 is the mobility of the fluids at the interface with the solids that acts as a phenomenological parameter (Pudasaini, 2012).
The fluid pressure acts only on the free fluids here, as the confined fluids are moved together with the solids. In Eq. (30), the second term is related to the nonNewtonian viscous force induced by gradients in solid concentration. The effect as described by Pudasaini (2012) is induced by a solid concentration gradient. In the case of unconfined fluids and unstructured solids (${f}_{\mathrm{sc}}=\mathrm{0},{f}_{\mathrm{fc}}=\mathrm{0}$); this force is identical to the description in the unstructured equations. Within our flow description, we see no direct reason to eliminate or alter this force with a variation in the fraction of confined fluids or structured solids. We only consider the interface between solids and free fluids as an agent that induces this effect, and therefore the gradient of the gradient of the solids and confined fluids ($\mathrm{\nabla}({\mathit{\alpha}}_{\mathrm{s}}+{f}_{\mathrm{fc}}{\mathit{\alpha}}_{\mathrm{f}})=\mathrm{\nabla}{\mathit{\alpha}}_{\mathrm{c}}$) is used instead of the total solid phase (∇α_{s}).
2.2.5 Drag force and virtual mass
Our description of the drag force follows the work of Pudasaini (2012), Pudasaini et al. (2018), where a generalized twophase drag model is introduced and enhanced. We split their work into a contribution from the fraction of structured solids (f_{sc}) and unconfined fluids (1−f_{fc}) (Eq. 31).
Here U_{T,c} is the terminal or settling velocity of the structured solids, U_{T,uc} is the terminal velocity of the unconfined solids, 𝒫 is a factor that combines solid and fluidlike contributions to the drag force, 𝒢 is the solidlike drag contribution, ℱ is the fluidlike drag contribution, and S_{p} is the smoothing function (Eqs. 32 and 34). The exponent j indicates the type of drag: linear (j=0) or quadratic (j=1).
Within the drag, the following functions are defined:
where M is a parameter that varies between 2.4 and 4.65 based on the Reynolds number (Pitman and Le, 2005). The factor 𝒫 that combines solid and fluidlike contributions to the drag is dependent on the volumetric solid content in the unconfined and unstructured materials $\left(\mathcal{P}={\left(\frac{{\mathit{\alpha}}_{\mathrm{s}}(\mathrm{1}{f}_{\mathrm{sc}})}{{\mathit{\alpha}}_{\mathrm{f}}\left(\mathrm{1}{f}_{\mathrm{fc}}\right)}\right)}^{\mathrm{m}}\right)$ with m≈1. Additionally, we assume the factor 𝒫 is zero for drag originating from the structured solids. As stated by Pudasaini and Mergili (2019), “as limiting cases: 𝒫 suitably models solid particles moving through a fluid”. In our model, the drag force acts on the unconfined fluid momentum (u_{uc}α_{f}(1−f_{fc})). For interactions between unconfined fluids and structured solids, larger blocks of solid structures are moving through fluids that contains solids of smaller size.
Virtual mass is similarly implemented based on the work of Pudasaini (2012) and Pudasaini and Mergili (2019) (Eq. 35). The adapted implementation considers the solids together with confined fluids to move through a freefluid phase.
Here ${C}_{\mathrm{DG}}=\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\mathrm{1}+\mathrm{2}{\mathit{\alpha}}_{\mathrm{c}}}{{\mathit{\alpha}}_{\mathrm{u}}}\right)$ is the drag coefficient.
2.2.6 Boundary conditions
Finally, following the work of Iverson and Denlinger (2001), Pitman and Le (2005), and Pudasaini (2012), a boundary condition is applied to the surface elements that contact the flow (Eq. 36).
Here N is the normal pressure on the surface element and S is the shear stress.
2.3 Depthaveraging
The majority of the depthaveraging in this work is analogous to the work of Pitman and Le (2005), Pudasaini (2012), and Pudasini and Mergili (2019). Depthaveraging through integration over the vertical extent of the flow can be done based on several useful and oftenused assumptions, e.g. $\frac{\mathrm{1}}{h}{\int}_{\mathrm{0}}^{h}x\mathrm{d}h=\stackrel{\mathrm{\u203e}}{x}$, for the velocities (u_{u} and u_{c}); solid, fluid, and confined fractions (α_{f}, α_{s}, f_{fc} and f_{sc}); and material properties (ρ_{u}, ϕ and c). Besides these similarities and an identical derivation of depthaveraged continuity equations, three major differences arise.
 i.
Fluid pressure. Previous implementations of generalized twophase debris flow equations have commonly assumed hydrostatic pressure ($\frac{\partial p}{\partial z}={g}^{z}$) (Pitman and Le, 2005; Pudasaini, 2012; Abe and Konagai, 2016). Here we follow this assumption for the fluid pressure at the base and solid pressure for unstructured material (Eqs. 37 and 38).
$$\begin{array}{}\text{(37)}& {\displaystyle}{P}_{{\mathrm{b}}_{s,u}}=\left(\mathrm{1}\mathit{\gamma}\right){\mathit{\alpha}}_{\mathrm{s}}{g}^{z}h\text{(38)}& {\displaystyle}{P}_{{\mathrm{b}}_{\mathrm{u}}}={g}^{z}h\end{array}$$Here $\mathit{\gamma}=\frac{{\mathit{\rho}}_{\mathrm{f}}}{{\mathit{\rho}}_{\mathrm{s}}}$ is the density ratio (not to be confused with a tensor index when used in superscript) (–).
However, larger blocks of structure material can have contact with the basal topography. Due to density differences, larger blocks of solid structures are likely to move along the base (Pailhia and Pouliquen, 2009; Iverson and George, 2014). If these blocks are saturated, water pressure propagates through the solid matrix and hydrostatic pressure is retained. However, in cases of an unsaturated solid matrix that connects to the base, hydrostatic pressure is not present there. We introduce a basal fluid pressure propagation factor $\mathcal{B}({\mathit{\theta}}_{\mathrm{eff}},\stackrel{\mathrm{\u203e}}{{d}_{\mathrm{sc}}},\mathrm{\dots})$ that describes the fraction of fluid pressure propagated through a solid matrix (with θ_{eff} as the effective saturation, and as $\stackrel{\mathrm{\u203e}}{{d}_{\mathrm{sc}}}$ the average size of structured solid matrix blocks). This results in a basal pressure equal to Eq. (39).
$$\begin{array}{}\text{(39)}& \begin{array}{rl}{P}_{{\mathrm{b}}_{\mathrm{c}}}& =\left(\mathrm{1}{f}_{\mathrm{sc}}\right)\left(\mathrm{1}\mathit{\gamma}\right){\displaystyle \frac{\left(\mathrm{1}{f}_{\mathrm{sc}}\right){\mathit{\alpha}}_{\mathrm{s}}}{\left(\mathrm{1}{f}_{\mathrm{fc}}\right){\mathit{\alpha}}_{\mathrm{f}}}}{g}^{z}h\\ & {f}_{\mathrm{sc}}\left(\mathrm{1}\mathit{\gamma}\right)\mathcal{B}{\displaystyle \frac{\left({f}_{\mathrm{sc}}\right){\mathit{\alpha}}_{\mathrm{s}}}{\left({f}_{\mathrm{fc}}\right){\mathit{\alpha}}_{\mathrm{f}}}}{g}^{z}h\end{array}\end{array}$$The basal pressure propagation factor (ℬ) should theoretically depend mostly on saturation level, similarly to the pedotransfer function, as a full saturation means perfect propagation of pressure through the mixture, and low saturation equates to minimal pressure propagation (Saxton and Rawls, 2006). Additionally, it should depend on pedotransfer functions and the size distribution of structured solid matrices within the mixture. For low saturation levels, it can be assumed that no fluid pressure is retained. Combined with an assumed soil matrix height being identical to the total mixture height, this results in ℬ=0. Assuming saturation of a structure's solids results in a full propagation of pressures and ℬ=1.
 ii.
Stress–strain relationship. Depthaveraging the stress–strain relationship in Eqs. (22) and (23) requires a vertical solution for the internal stress. First, we assume any nonnormal vertical terms are zero (Eq. 40). Commonly, Rankine's earth pressure coefficients are used to express the lateral earth pressure by assuming vertical stress to be induced by the basal solid pressure (Eqs. 41 and 42) (Pitman and Le, 2005; Pudasaini, 2012; Abe and Konagai, 2016).
$$\begin{array}{}\text{(40)}& {\displaystyle}{\mathit{\sigma}}^{zx}={\mathit{\sigma}}^{zy}=\phantom{\rule{0.125em}{0ex}}{\mathit{\sigma}}^{yz}=\phantom{\rule{0.125em}{0ex}}{\mathit{\sigma}}^{xz}=\mathrm{0}\text{(41)}& {\displaystyle}\stackrel{\mathrm{\u203e}}{{\mathit{\sigma}}^{zz}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{P}_{{\mathrm{b}}_{\mathrm{s}}},{\mathit{\sigma}}^{zz}{}_{\mathrm{b}}={P}_{{\mathrm{b}}_{\mathrm{s}}}\text{(42)}& {\displaystyle}{K}_{\mathrm{a}}={\displaystyle \frac{\mathrm{1}\mathrm{sin}\left(\mathit{\varphi}\right)}{\mathrm{1}+\mathrm{sin}\left(\mathit{\varphi}\right)}},{K}_{\mathrm{p}}={\displaystyle \frac{\mathrm{1}\mathrm{sin}\left(\mathit{\varphi}\right)}{\mathrm{1}+\mathrm{sin}\left(\mathit{\varphi}\right)}}\end{array}$$Here we enhance this with Bell's extension for cohesive soils (Eq. 45) (Xu et al., 2019). This lateral normaldirected stress term is added to the full stress–strain solution.
$$\begin{array}{}\text{(43)}& \stackrel{\mathrm{\u203e}}{{\mathit{\sigma}}_{xx}}=K{\mathit{\sigma}}_{zz}{}_{\mathrm{b}}\mathrm{2}c\sqrt{K}+{\displaystyle \frac{\mathrm{1}}{h}}{\int}_{\mathrm{0}}^{h}{\mathit{\sigma}}_{xx}\mathrm{d}h\end{array}$$Finally, the gradient in pressure of the lateral interfaces between the mixture is added as a depthaveraged acceleration term (Eq. 44).
$$\begin{array}{}\text{(44)}& {S}_{{x}_{\mathrm{c}}}={\mathit{\alpha}}_{\mathrm{c}}\left({\displaystyle \frac{\mathrm{1}}{h}}\left({\displaystyle \frac{\partial \left(h{\mathit{\sigma}}^{xx}\right)}{\partial x}}+{\displaystyle \frac{\partial \left(h{\mathit{\sigma}}^{yx}\right)}{\partial y}}\right)\right)+\mathrm{\dots}\end{array}$$  iii.
Depthaveraging other terms. While the majority of terms allow for depthaveraging as it was proposed by Pudasaini (2012), an exception arises. Depthaveraging of the vertical viscosity terms is required. The nonNewtonian viscous terms for the fluid phase were derived assuming a vertical profile in the volumetric solid phase content. Here, we alter the derivation to use this assumption only for the nonstructured solids, as opposed to the structured solids where $\frac{\partial {\mathit{\alpha}}_{\mathrm{s}}}{\partial z}=\mathrm{0}$.
$$\begin{array}{}\text{(45)}& \begin{array}{rl}{\int}_{\mathrm{b}}^{\mathrm{s}}{\displaystyle \frac{\partial}{\partial z}}\left({\displaystyle \frac{\partial {\mathit{\alpha}}_{\mathrm{s}}}{\partial z}}\left({u}_{\mathrm{u}}{u}_{\mathrm{c}}\right)\right)\mathrm{d}z& ={\left[{\displaystyle \frac{\partial {\mathit{\alpha}}_{\mathrm{s}}}{\partial z}}\left({u}_{\mathrm{u}}{u}_{\mathrm{c}}\right)\right]}_{\mathrm{b}}^{\mathrm{s}}\\ & =\left(\stackrel{\mathrm{\u203e}}{{u}_{\mathrm{u}}}\stackrel{\mathrm{\u203e}}{{u}_{\mathrm{c}}}\right){\left[{\displaystyle \frac{\partial {\mathit{\alpha}}_{\mathrm{s}}}{\partial z}}\right]}_{\mathrm{b}}^{\mathrm{s}}\\ & =\left(\stackrel{\mathrm{\u203e}}{{u}_{\mathrm{u}}}\stackrel{\mathrm{\u203e}}{{u}_{\mathrm{c}}}\right){\left[{\displaystyle \frac{\partial {\mathit{\alpha}}_{\mathrm{s}}}{\partial z}}\right]}_{\mathrm{b}}^{\mathrm{s}}\\ & ={\displaystyle \frac{\left(\stackrel{\mathrm{\u203e}}{{u}_{\mathrm{u}}}\stackrel{\mathrm{\u203e}}{{u}_{\mathrm{c}}}\right)\left(\mathrm{1}{f}_{\mathrm{sc}}\right)\mathit{\zeta}\stackrel{\mathrm{\u203e}}{{\mathit{\alpha}}_{\mathrm{s}}}}{h}}\end{array}\end{array}$$Here ζ is the shape factor for the vertical distribution of solids (Pudasaini, 2012). Additionally, the momentum balance of Pudasaini (2012) ignores any deviatoric stress (τ_{xy}=0), following Savage and Hutter (1989) and Pudasaini and Hutter (2007). Previously this term has been included by Iverson and Denlinger (2001), Pitman and Le (2005), and Abe and Kanogai (2016). Here we include these terms since a full stress–strain relationship is included.
2.3.1 Basal frictions
Additionally we add the Darcy–Weisbach friction, which is a Chézytype friction law for the fluid phase that provides drag (Delestre et al., 2014). This ensures that without solid phase a clear fluid does lose momentum due to friction from basal shear. This was successfully done in Bout et al. (2018) and was similarly assumed in Pudasaini and Fischer (2016) for fluid basal shear stress.
Here n is Manning's surface roughness coefficient.
2.3.2 Depthaveraged equations
The following set of equations is thus finally achieved for depthaveraged flow over sloping terrain (Eqs. 47–71).
Here X is the shape factor for vertical shearing of the fluid (X≈3 in Iverson and Denlinger, 2001), R is the precipitation rate and I is the infiltration rate.
2.3.3 Closing the equations
Viscosity is estimated using the empirical expression from Julien and O'Brien (1997), which relates dynamic viscosity to the solid concentration of the fluid (Eq. 72).
Here α is the first viscosity parameter and β the second viscosity parameter.
Finally, the settling velocity of small (d<100 µm) grains is estimated by Stokes equations for a homogeneous sphere in water. For larger grains (>1 mm), the equation by Zanke (1977) is used (Eq. 73).
In which U_{T} is the settling (or terminal) velocity of a solid grain, η is the dynamic viscosity of the fluid, ρ_{f} is the density of the fluid, ρ_{s} is the density of the solids and d is the grain diameter (m).
2.4 Implementation in the material point method numerical scheme
Implementing the presented set of equations into a numerical scheme requires considerations of that scheme's limitations and strengths (Stomakhin et al., 2013). Fluid dynamics are almost exclusively solved using an Eulerian finiteelement solution (Delestre et al., 2014; Bout et al., 2018). The diffusive advection part of such scheme typically does not degrade the quality of modelling results. Solid material, however, is commonly simulated with higher accuracy using an Lagrangian finiteelement method or discreteelement method (Maurel and Cumbescure, 2008; Stomakhin et al., 2013). Such schemes more easily allow for the material to maintain its physical properties during movement. Additionally, advection in these schemes does not artificially diffuse the material since the material itself is discretized, instead of the space (grid) on which the equations are solved. In our case, the material point method (MPM) provides an appropriate tool to implement the set of presented equations (Bui et al., 2008; Maurel and Cumbescure, 2008; Stomakhin et al., 2013). Numerous existing modelling studies have implemented in this method (Pastor et al., 2014; Abe and Kanogai, 2016). Here, we use the MPM method to create a twophase scheme. This allows the usage of finiteelement aspects for the fluid dynamics, which are so successfully described by the that method (particularly for water in larger areas; see Bout et al., 2018).
2.4.1 Mathematical framework
The mathematic framework of smooth particles solves differential equations using discretized volumes of mass represented by kernel functions (Libersky and Petschek, 1991; Bui et al., 2008; Stomakhin et al., 2013). Here, we use the cubic spline kernel as used by Monaghan (2000) (Eq. 74).
Here r is the distance, h is the kernel size and q is the normalized distance ($q=\frac{r}{h}$).
Using this function mathematical operators can be defined. The average is calculated using a weighted sum of particle values (Eq. 75), while the derivative depends on the function values and the derivative of the kernel by means of the chain rule (Eq. 76) (Libersky and Petschek, 1991; Bui et al., 2008).
Here ${W}_{ij}=W({x}_{i}{x}_{j},h)$ is the weight of particle j to particle I, and $r=\left\phantom{\rule{0.125em}{0ex}}{x}_{i}{x}_{j}\right$ is the distance between two particles. The derivative of the weight function is defined by Eq. (77).
Using these tools, the momentum equations for the particles can be defined (Eqs. 78–84). Here, we follow Monaghan (2000) and Bui et al. (2008) for the definition of artificial numerical forces related to stability. Additionally, stressbased forces are calculated on the particle level, while other momentum source terms are solved on a Eulerian grid with spacing h (identical to the kernel size).
Here i and j are indices indicating the particle, Π_{ij} is an artificial viscous force as defined by Eqs. (83) and (84), and ${F}_{ij}^{n}{R}_{ij}^{\mathit{\alpha}\mathit{\beta}}$ is an artificial stress term as defined by Eqs. (85) and (86).
Here ϵ_{0} is a small parameter ranging from 0 to 1, α_{Π} and β_{Π} are constants in the artificial viscous force (often chosen close to 1), and u_{sound} is the speed of sound in the material.
The conversion from particles to gridded values and vice versa depends on a grid basis function that weighs the influence of particle values for a grid centre. Here, a function derived from dyadic products of onedimensional cubic Bsplines is used as has been done previously by Steffen et al. (2008) and Stomakhin et al. (2013) (Eq. 84).
2.4.2 Particle placement
Particle placement is typically done in a constant pattern, as initial conditions have some constant density. The simplest approach is a regular square or triangular network, with particles on the corners of the network. Here, we use an approach that is more adaptable to spatially varying initial flow height. The R_{2} sequence approaches, with a regular quasirandom sequence, a set of evenly distributed points within a square (Roberts, 2020) (Eq. 85).
Here x_{n} is the relative location of the nth particle within a grid cell, and ${c}_{\mathrm{p}}={\left(\frac{\mathrm{9}+\sqrt{\mathrm{69}}}{\mathrm{18}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}+{\left(\frac{\mathrm{9}\sqrt{\mathrm{69}}}{\mathrm{18}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}\approx \mathrm{1.32471795572}$ is the plastic constant.
The number of particles placed for a particular flow height depends on the particle volume V_{I}, which is taken as a global constant during the simulation.
3.1 Flume setup
In order to validate the presented model, several controlled experiments were performed and reproduced using the developed equations. The flume setup consists of a steep incline, followed by a nearflat runout plane (Fig. 3). A massive obstacle is placed on the separation point of the two planes. This blocks the path of twofifths of the width of the moving material. For the exact dimensions of both the flume parts and the obstacle, see Fig. 3.
Two tests were performed whereby a cohesive granular matrix was released at the upper part of the flume setup. Both of these volumes had dimensions of $\mathrm{0.2}\times \mathrm{0.3}\times \mathrm{0.25}$ m (height, length, width). For both of these materials, a mixture of highorganiccontent siltyclay soils where used. The material's strength parameters were obtained using triaxial testing (cohesion, internal friction angle Young's modulus and Poisson ratio). The first set of materials properties were c=26.7 kPa and ϕ=28^{∘}. The second set of materials properties were c=18.3 kPa and ϕ=27^{∘}. For both of the events, pre and postrelease elevation models were made using photogrammetry. The model was set up to replicate the situations using the measured input parameters. Numerical settings were chosen as follows: $\mathit{\{}{\mathit{\alpha}}_{\mathrm{s}}=\mathrm{0.5},{\mathit{\alpha}}_{\mathrm{f}}=\mathrm{0.5},{f}_{\mathrm{sc}}=\mathrm{1.0},{f}_{\mathrm{fc}}=\mathrm{1.0},{\mathit{\rho}}_{\mathrm{f}}=\mathrm{1000}$, ${\mathit{\rho}}_{\mathrm{s}}=\mathrm{2400},E=\mathrm{12}\cdot {\mathrm{10}}^{\mathrm{6}}\phantom{\rule{0.125em}{0ex}}\mathrm{Pa},K=\mathrm{23}\cdot {\mathrm{10}}^{\mathrm{6}}\phantom{\rule{0.125em}{0ex}}\mathrm{Pa},\mathit{\psi}=\mathrm{0},{\mathit{\alpha}}_{\mathrm{\Pi}}=\mathrm{1},{\mathit{\beta}}_{\mathrm{\Pi}}=\mathrm{1},X,\mathit{\zeta},j=\mathrm{2},{u}_{\mathrm{sound}}=\mathrm{600},\phantom{\rule{0.125em}{0ex}}\mathrm{d}x=\mathrm{10},{V}_{\mathrm{I}}=\mathrm{0}\phantom{\rule{0.125em}{0ex}}\mathrm{m}/\mathrm{s},h=\mathrm{10},n=\mathrm{0.1},\mathit{\alpha}=\mathrm{1},\mathit{\beta}=\mathrm{10},M=\mathrm{2.4},\mathcal{B}=\mathrm{0},{N}_{R}=\mathrm{15}\phantom{\rule{0.125em}{0ex}}\mathrm{000},{N}_{RA}=\mathrm{30}\mathit{\}}$. Calibration was performed by means of input variation. The solid fraction and elastic and bulk modulus were varied between 20 % and 200 % of their original values with increments of 10 %. Accuracy was assessed based on the percentage accuracy of the deposition (comparison of the modelled vs. the observed presence of material).
3.2 Results
Both the mapped extent of the material after flume experiments, as the simulation results are shown in Fig. 5. Calibrated values for the simulations are {α_{s}=0.45, $E=\mathrm{21.6}\cdot {\mathrm{10}}^{\mathrm{6}}\phantom{\rule{0.125em}{0ex}}\mathrm{Pa},K=\mathrm{13.8}\cdot {\mathrm{10}}^{\mathrm{6}}\phantom{\rule{0.125em}{0ex}}\mathrm{Pa}$ }.
As soon as the block of material impacts the obstacle, stress increases as the moving object is deformed. This stress quickly propagates through the object. Within the scenario with lower cohesive strength, as soon as the stress reached beyond the yield strength, degradation of strength parameters took place. In the results, a fracture line developed along the corner of the obstacle into the length direction of the moving mass. Eventually, this fracture developed to half the length of the moving body and severe deformation resulted. As was observed from the tests, the first material experienced a critical fracture, while the second test resulted in moderate deformation near the impact location. Generally, the results compare well with the observed patterns, although the exact shape of the fracture is not replicated. Several reasons might be the cause of the moderately accurate fracture patterns. Other studies used a more controlled setup where uncertainties in applied stress and material properties were reduced. Furthermore, the homogeneity of the material used in the tests can not completely assumed. Realistically, minor alterations in compression used to create the clay blocks have left spatial variation in density, cohesion and other strength parameters.
4.1 Numerical setup
In order to further investigate some of the behaviours of the model and highlight the novel types of mass movement dynamics that the model implements, several numerical tests have been performed. The setup of these tests is shown in Fig. 6.
Numerical settings were chosen for three different blocks with equal volumes but distinct properties. Cohesive strength and the bulk modulus were varied (see Fig. 6). Remaining parameters were chosen as follows: $\mathit{\{}{\mathit{\alpha}}_{\mathrm{s}}=\mathrm{0.5},{\mathit{\alpha}}_{\mathrm{f}}=\mathrm{0.5},{f}_{\mathrm{sc}}=\mathrm{1.0},{f}_{\mathrm{fc}}=\mathrm{1.0},{\mathit{\rho}}_{\mathrm{f}}=\mathrm{1000}\phantom{\rule{0.125em}{0ex}}\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{3}}$, ${\mathit{\rho}}_{\mathrm{s}}=\mathrm{2400}\phantom{\rule{0.125em}{0ex}}\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{3}},E=\mathrm{1}e\mathrm{12}\phantom{\rule{0.125em}{0ex}}\mathrm{Pa},\mathit{\psi}=\mathrm{0},{\mathit{\alpha}}_{\mathrm{\Pi}}=\mathrm{1},{\mathit{\beta}}_{\mathrm{\Pi}}=\mathrm{1},X,\mathit{\zeta},j=\mathrm{2},{u}_{\mathrm{sound}}=\mathrm{600}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}},\mathrm{d}x=\mathrm{10}m,{V}_{\mathrm{I}}=\mathrm{0}\phantom{\rule{0.125em}{0ex}}\mathrm{m}/\mathrm{s},h=\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{m},n=\mathrm{0.1},\mathit{\alpha}=\mathrm{1},\mathit{\beta}=\mathrm{10},M=\mathrm{2.4},\mathcal{B}=\mathrm{0},{N}_{R}=\mathrm{15}\phantom{\rule{0.125em}{0ex}}\mathrm{000},{N}_{RA}=\mathrm{30}\mathit{\}}$.
4.2 Results
Several time slices for the described numerical scenarios are shown in Figs. 7 and 8.
Fractures develop in the mass movements based on acceleration differences and cohesive strength. For scenario 2A, the stress state does not reach beyond the yield surface, and all material is moved as a single block. Scenario 2B, which features lowered cohesive strength, fractures and the masses separate based on the acceleration caused by slopes.
Fracturing behaviour can occur in MPM schemes due to numerical limitations inherent in the usage of a limited integration domain. Here, validation of real physically based fracturing is present in the remaining cohesive fraction. This value only reduces in the case of plastic yield, where increasing strain degrades strength parameters according to our proposed criteria. Numerical fractures would thus have a cohesive fraction of 1. In all simulated scenarios, such numerical issues were not observed.
Fragmentation occurs due to spatial variation in acceleration in the case of scenario 3A and 3B. For scenario 3A, the yield surface is not reached and the original structure of the mass is maintained during movement. For 3C, fragmentation is induced by lateral pressure and buoyancy forces alone. Scenario 3B experiences slight fragmentation at the edges of the mass but predominantly fragments when reaching the valley, after which part of the material is accelerated to count to the velocity of the mass. For all the shown simulations, fragmentation does not lead to significant phase separation since virtual mass and drag forces converge the separate phase velocities to their mixtureaveraged velocity. The strength of these forces partly depends on the parameters; effects of more immediate phase separation could be studied if other parameters are used as input.
A variety of existing landslide models simulate the behaviour of lateral connected material through a nonlinear, nonNewtonian viscous relationship (von Boetticher et al., 2017; Fornes et al., 2017; Pudasaini and Mergili, 2019; Greco et al., 2019). These relationships include a yield stress and are usually regularized to prevent singularities from occurring. While this approach is incredibly powerful, it is fundamentally different from the work proposed here. These viscous approaches do not distinguish between elastic or plastic deformation and typically ignore deformations if stress is insufficient. Additionally, fracturing is not implemented in these models. The approach taken in this work attempts to simulate a full stress–strain relationship with Mohr–Coulomb type yield surface. This provides new types of behaviour and can be combined with nonNewtonian viscous approaches as mentioned above. A major downside to the presented work is the steep increase in computational time required to maintain an accurate and stable simulation. Commonly, a near 100 times increase in computational time has been observed during the development of the presented model.
The presented model shows a good likeness to flume experiments, and numerical tests highlight behaviour that is commonly observed for landslide movements. There are, however, inherent scaling issues and the material used in the flume experiments is unlikely to form larger landslide masses. The measured physical strength parameters of the material used in the flume experiments would not allow for sustained structured movement at larger scales. There is thus the need for more realscale validation cases. The application of the presented type of model is most directly noticeable for blocktype landslide movements that have fragmented either upon impact with some obstacle or during the transition phase. Of importance here is that the moment of fragmentation is often not reported in studies on fastmoving landslides, potentially due to the complexities involved in knowing the details of this behaviour from postevent evidence. Validation would therefore have to occur in cases where deposits are not fully fragmented, indicating that this process was ongoing during the whole movement duration. The spatial extent of initiation and deposition would then allow validation of the model. Another major opportunity for validation of the novel aspects of the model is the full threedimensional application to landslides that were reported to have lubrication effects due to fragmentation of the lower fraction of flow due to shear.
An important point of consideration in the development of complex multiprocess generalized models is the applicability. As a detailed investigative research tool, these models provide a basic scenario of usage. However, both for research and beyond for applicability in disaster risk reduction decision support, the benefit drawn from these models depends on the practical requirement for parameterization and the computational demands for simulation. With an increasing complexity in the description of multiprocess mechanics comes the requirement for more measured or estimated physical parameters. Inspection of the presented method shows that in principle, a minor amount of new parameters are introduced. The cohesive strength, a major focus of the model, becomes highly important depending on the type of movement being investigated. Additionally, the bulk and elastic modulus are required. These three parameters are common simulation parameters in geotechnical research and can be obtained from common tests on sampled material (Alsalman et al., 2015). Finally, the basal pressure propagation parameter (ℬ) is introduced. However, within this work the value of this parameter is chosen to have a constant value of 1. As a result, the model does require additional parameters, although these are relatively easy to obtain with accuracy.
There are a variety of aspects of the model that could be significantly improved. Here, we list several major opportunities for future research.

Groundwater mechanics. The presented model allows for the a solid or granular matrix to be present within the flow. We have assumed the flows in and out of these matrices are small enough that they can be ignored. In reality, there is a fluid flux in and out of structured solids. This could occur both due to pressure differences and due to stress and strain of the structured solids. Implementing this kind of mechanic requires a dynamic, solidpropertiesdependent, soil water retention curve (Van Looy et al., 2017). An example of MPM soil mechanics with dynamic groundwater implementation can be found in Bandera et al. (2016).

Implementing entrainment and deposition. Current equations for entrainment (erosion with major graingrain interactions) are limited to unstructured mixture flows (Iverson, 2012; Iverson and Ouyang, 2015; Cuomo et al., 2016; Pudasaini and Fischer, 2016). Extending these models to include a contribution from structured solids would be required to implement entrainment in the presented work.

Separation of phases. A major assumption in the presented work is that the velocities of structured solids, free solids and confined fluids are all equal. In reality, there might be separation of structured and freesolid phases. Additionally, we already discussed the possibility of influx and outflux of confined fluids from the solid matrix. Recent innovations in threephase mixture flows might be used to extend the presented work to a three, four or fivephase model by separating free solids and confined fluids or adding a Bingham viscous solid–fluid phase (Pudasaini and Mergili, 2019). However, while this would implement an additional process, it would significantly increase the complexity of the equations (in an exponential manner with relation to the number of phases) and the numerical solutions, which could hinder practical applicability.

Application to large slowmoving landslides. When confined fluids would act as a distinct phase guided by the mechanics of water flow in granular matrix, groundwater pressures and movement through the structured solids could be described. This might enable the model to perform detailed deformation and groundwater simulation of large slowmoving landslides.

Numerical improvements. Numerical techniques for particlebased discretized methods (SPH, MPM) have been proposed in the literature. A common issue is numerical fracturing of materials when particle strain increases beyond the length of the kernel function. Due to this, the connection between particles is lost, and fracturing occurs as an artefact of the numerical method. This issue is partly solved by the artificial stress term that is also used by Bui et al. (2008). Additionally, a geometric subdivide, as used by Xu et al. (2012) and Li et al. (2015), could counter these artificial fractures. Implementing this technique does require additional work to maintain mass and momentum conservation.

Threedimensional solutions. In a variety of scenarios, the assumptions made in the depthaveraged application of flow models are invalid. A common example is the impact of mass movements into lakes or other large water bodies. In such cases, the vertical velocity and concentration variables are not well described by their depthaveraged counterparts. Additionally, the lubrication effect of basal fragmentation of landslides due to shear can not be described without velocity profiles and a vertical stress solution. A full threedimensional application would therefore have the potential to increase understanding of these important processes.
We have presented a novel generalized mass movement model that can describe both unstructured mixture flows and structured movements of Mohr–Coulombtype material. The presented equations are part of the continuous development of the OpenLISEM hazard model, an opensource tool for physically based multihazard simulations. The model builds on the works of Pudasaini (2012) and Bui et al. (2008) to develop a single holistic set of equations. The model was implemented in a GPUbased material point method (MPM) Code. The equations were validated on flume experiments and numerical tests that highlight the new movement dynamics possible with the presented model. The integration of cohesive structure and a full stress–strain relationship for the structured solids allows for movement of blocktype slides as a single whole. Interactions with terrain, other flow masses or obstacles lead to elastoplastic deformation and eventually fragmentation. This type of selfalteration of flow properties is novel with mass movement models. Although the presented equations can provide additional detail for specific mass movement types, applicability of the model for real events need to be investigated as computational costs are significantly increased.
The presented simulation both validates the basic behaviour of the model and highlights the types of flow dynamics made possible by the presented equations. The model's dependency on cohesive strength and internal friction angles matches the flume experiments. The numerical examples show commonly described behaviour for landslide movements. Although the simulations compare well to the flume experiments, validation is required for realscale application to various types of mass movements. Additionally, the presented equations still lack descriptions of processes that might become important. Separating the fluid and solid phases (as in Pudasaini and Mergili, 2019) could improve flow dynamics and phase separation. With added groundwater mechanics, such as those in Bandera et al. (2016), slowmoving landslide simulations might be described.
h  flow height 
s  solid phase 
f  fluid phase 
sc  structured solid phase 
fc  confined fluid phase 
ρ_{f}  density of fluids 
ρ_{s}  density of solids 
α_{f}  volumetric fluidphase fraction 
α_{s}  volumetric solidphase fraction 
f_{sc}  fraction of solids that is structured (confining) 
f_{fc}  fraction of fluids that is confined 
α_{c}  volumetric fraction of solids, structured solids and confined fluids 
α_{u}  volumetric fraction of free fluids (unconfined phase). 
ρ_{sc}  volumeaveraged density of the solids and confined fluids 
u_{u}  velocity of the unconfined phase (free fluids) 
u_{c}  velocity of the solids, confining solids and confined fluids 
u_{s}  velocity of the solids 
f  body force 
M_{DG}  drag force 
M_{vm}  virtual mass force 
T_{c}  stress tensor for the solids, confining solids and confined fluids 
T_{u}  stress tensor for the freefluid phase 
σ  stress tensor 
$\dot{s}$  deviatoric shear stress rate tensor 
δ  Kronecker delta 
${\dot{\mathit{\u03f5}}}_{\mathrm{plastic}}$  plastic strain rate 
${\dot{\mathit{\u03f5}}}_{\mathrm{elastic}}$  elastic strain rate 
λ  plastic multiplier rate 
g  plastic potential function 
${\dot{\mathit{\u03f5}}}_{\mathrm{total}}$  total strain rate 
$\dot{e}$  deviatoric strain rate 
ν  Poisson's ratio 
E  elastic Young's modulus 
G  shear modulus 
K  bulk elastic modulus 
f(I_{1},J_{2})  yield surface or yield criterion 
g(I_{1},J_{2})  plastic potential function 
ψ  dilatancy angle 
I_{1}  first stress invariant 
J_{2}  second stress invariant 
α_{ϕ}  first Drucker–Prager material constant 
k_{c}  second Drucker–Prager material constant 
$\dot{\mathit{\omega}}$  spin rate tensor 
ϵ_{v0}  initial volumetric strain 
ϵ_{v}  volumetric strain 
c_{0}  initial cohesion 
τ_{f}  fluid Cauchy stress tensor 
p_{f}  fluid pressure 
η_{f}  fluids dynamic viscosity 
𝒜  mobility of the fluid at the interface 
𝒞_{DG}  drag coefficient 
U_{T,c}  settling velocity of the solids, structured solids and confined fluids 
U_{T,uc}  settling velocity of the unstructured solids 
ℱ  drag contribution from solidlike drag 
𝒢  drag contribution from fluidlike drag 
S_{p}  smoothing function 
𝒦  absolute total mass flux 
M(Re_{p})  empirical function weakly dependent on the Reynolds number 
𝒫  partitioning parameter for the fluid and solidlike contributions to drag 
m  an exponent for P 
𝒞_{VMG}  virtual mass coefficient 
S  norm of the shear force 
N  normal force on a plane element 
g  gravitational acceleration 
${P}_{{\mathrm{b}}_{s,u}}$  basal pressure from the unconfined phase 
${P}_{{\mathrm{b}}_{u}}$  basal pressure from the free fluids 
${P}_{{\mathrm{b}}_{\mathrm{c}}}$  basal pressure from the solids, structured solids and confined fluids 
ℬ  pressure propagation factor for structured solids 
K_{a}  active lateral earth–pressure coefficient 
K_{p}  passive lateral earth–pressure coefficient 
ζ  shape factor for the vertical gradient in solid concentration 
n  Manning's surface roughness coefficient 
X  shape factor for the vertical fluid velocity profile 
Re_{p}  particle Reynolds number 
N_{R}  Reynolds number 
N_{RA}  interfacial Reynolds number 
H  typical height of the flow 
L  typical length of the flow 
α  first viscosity parameter 
β  second viscosity parameter 
d  grain diameter 
W  kernel weight function 
r  distance 
h  kernel width (not to be confused with the flow height) 
q  normalized particle distance 
Π_{ij}  an artificial viscosity term 
${F}_{ij}^{n}{R}_{ij}^{\mathit{\alpha}\mathit{\beta}}$  an artificial stress term 
ϵ_{0}  a constant parameter for the artificial stress term 
α_{Π} and β_{Π}  constants in the artificial viscous force 
u_{sound}  speed of sound in the material 
N(x)  grid kernel function 
c_{p}  plastic coefficient 
If the state of the stress tensor lies beyond the yield surface, either due to degradation of strength parameters or building numerical errors, a correction must be applied. We implement the correction scheme used by Bui et al. (2008). This scheme considers two primary ways in which the stress can have an undesired state: tension cracking and imperfectly plastic stress.
B1 Tension cracking
In the case of tension cracking, the stress state has moved beyond the apex of the yield surface, as described by Chen and Mizuno (1990). The employed solution in this case is to remap the stress tensor along the I_{1} axis to be at this apex. The apex is provided by the following yield function:
To solve for this condition, the nondeviatoric stress state is increased (since ${I}_{\mathrm{1}}\frac{{k}_{\mathrm{c}}}{{\mathit{\alpha}}_{\mathit{\varphi}}}$ is negative) to lie perpendicular to the apex point on the I_{1} axis as follows:
B2 Imperfect plastic stress
Imperfect plastic stress described the state where the stress tensor lies above the apex but beyond the yield criterion and thus have more stress than is supported by the failure criteria that is set. This criterion is simply the yield surface itself (Eq. B3).
For this state, remapping is done by scaling of the J_{2} value (Eqs. B4, B5 and B6).
The model presented in this article is part of the continued development of the OpenLISEM modelling tools. The most recent set of equations were implemented in the opensource alpha version of OpenLISEM hazard model 2.0a. Here we describe the details of the implementation of the model into software.
C1 Hybrid MPM
We utilize the MPM framework to be able to discretize part of the equations on a Eulerian regular grid and part of the equations on the Lagrangian particles. Our distinct take on this method is the representation of the fluid phase entirely as a finiteelement solution, while solids are simulated as discrete particle volumes. This allows the model to use the major benefits that are present when depthaveraged fluid flow is simulated in a grid. Both numerical efficiency and highaccuracy coupling with hydrology are lacking in particle methods. For the solid phase, nondissipative advection, fracturing and stiffness is a major benefit of the MPM approach. Since our model assumed that confined fluids share their velocity with the solids, we advect the confined fluids as part of the particles. Total fluid volume is then calculated from the free fluids in the finiteelement data and the gridded particle data. A flow chart of the software setup is provided in Fig. 6.
C2 Finiteelement solution
We use a regular Cartesian grid to describe the modelling domain. Terrain and cellboundarybased variables are reproduced using the Monotonic Upstreamcentered Scheme for Conservation Laws (MUSCL) piecewise linear reconstruction (Delestre et al., 2014). For each cell boundary, a left and right estimation of acceleration terms, velocity updates, and new discharges is made. The left estimates use leftreconstructed variables, while the right estimates use rightreconstructed variables. The final average flux through the boundary determines the actual mass and momentum transfer. Local acceleration is averaged from the right estimate of the left boundary and left estimate of the right boundary. An additional benefit of the used scheme is the automatic estimation of continuous and discontinuous terrain. The piecewise linear reconstructions do not guarantee smooth terrain; for sharp locally variable terrain, pressure terms from vertical walls arise that block momentum. These terms allow for better estimation of momentum loss by barriers but can be turned off if required for the simulated scenario.
C3 GPU acceleration using OpenCL/OpenGL
In order to create a more efficient setup, both the finiteelement and particle interactions are performed on the GPU. We utilize the OpenCL API to compile kernels written in Cstyle language. These kernels are compiled at the start of the simulation and thereby allow for easy customization by users. While the usage of OpenCL 1.1 forces the usage of single precision floating point numbers, it allows for a wider range of GPU types to be supported. Finiteelement solutions on the GPU are straightforward, as maps are a basic data storage type for graphical processing units. Particles are stored as single precision floating point arrays. Within the framework of MPM, iteration of particles within a kernel is required for each time step and particle. This effectively means O(n^{2}) operations are required. Significant efficiency improvements are obtained by precalculation sorting. Particles are sorted based on their location within the finiteelement grid. Based on the ID of the grid cell, a bitonic mergesort is performed. This sorting algorithm works seamlessly on parallel architecture and operates as O(nlog ^{2}(n)) (Batcher, 1968). Following this, a raster is allocated to store the first indexed occurrence within the sorted list of particles of that grid cell. Since the kernel used for the presented work extends at most to a full width of two grid cells, we must iterate over all particles present in nine neighbouring grid cells.
A final benefit to the usage of OpenCL is direct access to simulation variables for visualization in OpenGL using the OpenGL/OpenCL interoperability functionality. The builtin viewing window of OpenLISEM hazard 2.0a directly uses the data to draw particles, shapefiles and grid data using customizable shaders written in the OpenGL shader language.
All code and data used within this work are made opensource as part of the continuous development of the OpenLISEM hazard model under the GNU General Public Licence v3.0. The code and the data are hosted on GitHub (https://github.com/bastianvandenbout/OpenLISEMHazard2.0PreRelease, https://doi.org/10.17026/dansxz42tut) (van den Bout, 2020). Both binaries and a copy of the source code are also available on SourceForge, where the manual and compilation guide can similarly be found (https://sourceforge.net/projects/lisem/, van den Bout, 2021). Finally, more information can be found at the blog of Bastian van den Bout and Victor G. Jetten (https://blog.utwente.nl/lisem/, van den Bout and Jetten, 2020)
The software and its user interface are written for Windows, but platformindependent libraries are used and compilation can be performed on other platforms.
Hardware requirements for the usage of the model are a 64bit operating system that can compile all required external libraries (see the manual for a full list and description), a graphical processing unit conforming to at least the OpenCL 1.2 standard, and support for both OpenGL 4.2 and OpenGL/OpenCL interoperability. Additionally, an approximate 500 MB of hard drive space and 750 MB of memory must also be available.
During this work, model derivation and conceptual work was performed by BVDB, VGJ, CJVW, TVA and OM. The presented flume experiments were carried out by CXT, WH and BVDB. Finally, manuscript production, data analysis, and numerical implementation of the methods and equations were carried out by BVDB.
The authors declare that they have no conflict of interest.
This paper was edited by Andrew Wickert and reviewed by two anonymous referees.
Aaron, J. and Hungr, O.: Dynamic simulation of the motion of partiallycoherent landslides, Eng. Geol., 205, 1–11, 2016.
Abe, K. and Konagai, K.: Numerical simulation for runout process of debris flow using depthaveraged material point method, Soils Found., 56, 869–888, 2016.
Alsalman, M., Myers, M., and SharfAldin, M.: Comparison of multistage to single stage triaxial tests, in: 49th US Rock Mechanics/Geomechanics Symposium, Francisco, California, 28 June 2015, ARMA2015767, 2015.
Bandara, S., Ferrari, A., and Laloui, L.: Modelling landslides in unsaturated slopes subjected to rainfall infiltration using material point method, Int. J. Numer. Anal. Met., 40, 1358–1380, 2016.
Batcher, K. E.: Sorting networks and their applications, in: Proceedings of the spring joint computer conference, Atlantic City, New Jersey, 30 April–2 May 1968, 307–314, 1968.
Beutner, E. C. and Gerbi, G. P.: Catastrophic emplacement of the Heart Mountain block slide, Wyoming and Montana, USA, Geol. Soc. Am. Bull., 117, 724–735, 2005.
Bieniawski, Z. T.: Mechanism of brittle fracture of rock: part I – theory of the fracture process, Int. J. Rock. Mech. Min., 4, 395–406, https://doi.org/10.1016/01489062(67)900307, 1967.
Bout, B., Lombardo, L., van Westen, C. J., and Jetten, V. G.: Integration of twophase solid fluid equations in a catchment model for flashfloods, debris flows and shallow slope failures, Environ. Modell. Softw., 105, 1–16, 2018.
Bui, H. H., Fukagawa, R., Sako, K., and Ohno, S.: Lagrangian meshfree particles method (SPH) for large deformation and failure flows of geomaterial using elastic–plastic soil constitutive model, Int. J. Numer. Anal. Met., 32, 1537–1570, 2008.
Chen, W.F. and Mizuno, E.: Nonlinear analysis in soil mechanics, Elsevier, Amsterdam, 1990.
Cohen, D., Lehmann, P., and Or, D.: Fiber bundle model for multiscale modeling of hydromechanical triggering of shallow landslides, Water resources research, 45, W10436, https://doi.org/10.1029/2009WR007889, 2009.
Corominas, J., van Westen, C., Frattini, P., Cascini, L., Malet, J.P., Fotopoulou, S., Catani, F., Van Den Eeckhaut, M., Mavrouli, O., and Agliardi, F.: Recommendations for the quantitative analysis of landslide risk, Bull. Eng. Geol. Environ., 73, 209–263, 2014.
Corominas, J., Matas, G., and RuizCarulla, R.: Quantitative analysis of risk from fragmental rockfalls, Landslides, 16, 5–21, 2019.
Cuomo, S., Pastor, M., Capobianco, V., and Cascini, L.: Modelling the space–time evolution of bed entrainment for flowlike landslides, Eng. Geol., 212, 10–20, 2016.
David, L. G. and Richard, M.: A twophase debrisflow model that includes coupled evolution of volume fractions, granular dilatancy, and porefluid pressure, Ital. J. Eng. Geol. Environ., 43, 415–424, 2011.
Davies, T. and McSaveney, M.: The role of rock fragmentation in the motion of large landslides, Eng. Geol., 109, 67–79, 2009.
Davies, T., McSaveney, M., and Beetham, R.: Rapid block glides: slidesurface fragmentation in New Zealand's Waikaremoana landslide, Q. J. Eng. Geol. Hydroge., 39, 115–129, 2006.
Delaney, K. B. and Evans, S. G.: The 1997 Mount Munday landslide (British Columbia) and the behaviour of rock avalanches on glacier surfaces, Landslides, 11, 1019–1036, 2014.
Delestre, O., Cordier, S., Darboux, F., Du, M., James, F., Laguerre, C., Lucas, C., and Planchon, O.: FullSWOF: A software for overland flow simulation, in: Advances in hydroinformatics, Springer, Singapore, 221–231, 2014.
De Vuyst, T. and Vignjevic, R.: Total Lagrangian SPH modelling of necking and fracture in electromagnetically driven rings, Int. J. Fracture, 180, 53–70, 2013.
Dhanmeher, S.: Crack pattern observations to finite element simulation: An exploratory study for detailed assessment of reinforced concrete structures, Delft University, Delft, the Netherlands, 2017.
Drew, D. A.: Mathematical modeling of twophase flow, Annu. Rev. Fluid Mech., 15, 261–291, 1983.
Dufresne, A., Geertsema, M., Shugar, D., Koppes, M., Higman, B., Haeussler, P., Stark, C., Venditti, J., Bonno, D., and Larsen, C.: Sedimentology and geomorphology of a large tsunamigenic landslide, Taan Fiord, Alaska, Sediment. Geol., 364, 302–318, 2018.
Fornes, P., Bihs, H., Thakur, V. K. S., and Nordal, S.: Implementation of nonNewtonian rheology for Debris Flow simulation with REEF3D, in: Eproceedings of the 37th IAHR World Congress, Kuala Lumpur, Malaysia, 13–18 August 2017, 2017.
Grady, D. E. and Kipp, M. E.: Continuum modelling of explosive fracture in oil shale, Int. J. Rock Mechanics Min., 17, 147–157, 1980.
Greco, M., Di Cristo, C., Iervolino, M., and Vacca, A.: Numerical simulation of mudflows impacting structures, J. Mt. Sci., 16, 364–382, 2019.
Griffiths, D. and Lane, P.: Slope stability analysis by finite elements, Geotechnique, 49, 387–403, 1999.
Hayir, A.: The effects of variable speeds of a submarine block slide on nearfield tsunami amplitudes, Ocean Eng., 30, 2329–2342, 2003.
Hušek, M., Kala, J., Hokeš, F., and Král, P.: Influence of SPH regularity and parameters in dynamic fracture phenomena, Proc. Eng., 161, 489–496, 2016.
Hutter, K., Svendsen, B., and Rickenmann, D.: Debris flow modeling: a review, Continuum Mech. Therm., 8, 1–35, 1994.
Ishii, M.: Thermofluid dynamic theory of twophase flow, NASA Sti/recon Technical Report A, 75, 29657, https://doi.org/10.1007/9781441979858, 1975.
Ishii, M. and Zuber, N.: Drag coefficient and relative velocity in bubbly, droplet or particulate flows, AICHE J., 25, 843–855, 1979.
Iverson, R. M.: Elementary theory of bedsediment entrainment by debris flows and avalanches, J. Geophys. Res.Earth, 117, F03006, https://doi.org/10.1029/2011JF002189, 2012.
Iverson, R. M. and Denlinger, R. P.: Flow of variably fluidized granular masses across threedimensional terrain: 1. Coulomb mixture theory, J. Geophys. Res.Solid Ea., 106, 537–552, 2001.
Iverson, R. M. and George, D. L.: A depthaveraged debrisflow model that includes the effects of evolving dilatancy. I. Physical basis, P. Roy. Soc. A Math. Phy., 470, 20130819, https://doi.org/10.1098/rspa.2013.0819, 2014.
Iverson, R. M. and George, D. L.: Modelling landslide liquefaction, mobility bifurcation and the dynamics of the 2014 Oso disaster, Géotechnique, 66, 175–187, 2016.
Iverson, R. M. and Ouyang, C.: Entrainment of bed material by Earthsurface mass flows: Review and reformulation of depthintegrated theory, Rev. Geophys., 53, 27–58, 2015.
Jakob, M., Hungr, O., and Jakob, D. M.: Debrisflow hazards and related phenomena, Springer, Berlin, Heidelberg, 2005.
Kaklauskas, G. and Ghaboussi, J.: Stressstrain relations for cracked tensile concrete from RC beam tests, J. Struc. Eng., 127, 64–73, 2001.
Kern, J. S.: Evaluation of soil water retention models based on basic soil physical properties, Soil Sc. Soc. Am. J., 59, 1134–1141, 1995.
Kjekstad, O. and Highland, L.: Economic and social impacts of landslides, in: Landslides–disaster risk reduction, SpringerVerlag Berlin Heidelberg, 573–587, 2009.
Li, C., Wang, C., and Qin, H.: Novel adaptive SPH with geometric subdivision for brittle fracture animation of anisotropic materials, Visual Comput., 31, 937–946, 2015.
Libersky, L. D. and Petschek, A. G.: Smooth particle hydrodynamics with strength of materials, in: Advances in the freeLagrange method including contributions on adaptive gridding and the smooth particle hydrodynamics method, Springer, Berlin, Heidelberg, 248–257, 1991.
Loehnert, S. and MuellerHoeppe, D. S.: Multiscale methods for fracturing solids, in: IUTAM symposium on theoretical, computational and modelling aspects of inelastic media, Springer, Dordrecht, the Netherlands, pp. 79–87, 2008.
Luna, B. Q., Remaître, A., Van Asch, T. W., Malet, J.P., and Van Westen, C.: Analysis of debris flow behavior with a one dimensional runout model incorporating entrainment, Eng. Geol., 128, 63–75, 2012.
Ma, G., Wang, Q., Yi, X., and Wang, X.: A modified SPH method for dynamic failure simulation of heterogeneous material, Math. Probl. Eng., 2014, 808359, https://doi.org/10.1155/2014/808359, 2014.
Matsui, T. and San, K.C.: Finite element slope stability analysis by shear strength reduction technique, Soils Found., 32, 59–70, 1992.
Maurel, B. and Combescure, A.: An SPH shell formulation for plasticity and fracture analysis in explicit dynamics, Int. J. Num. Meth. Eng., 76, 949–971, 2008.
Menin, R., Trautwein, L. M., and Bittencourt, T. N.: Smeared crack models for reinforced concrete beams by finite element method, Revista IBRACON de Estruturas e Materiais, 2, 166–200, 2009.
Mergili, M., Fischer, J.T., Krenn, J., and Pudasaini, S. P.: r.avaflow v1, an advanced opensource computational framework for the propagation and interaction of twophase mass flows, Geosci. Model Dev., 10, 553–569, https://doi.org/10.5194/gmd105532017, 2017.
Mergili, M., Emmer, A., Juřicová, A., Cochachin, A., Fischer, J. T., Huggel, C., and Pudasaini, S. P.: How well can we simulate complex hydro‐geomorphic process chains? The 2012 multi‐lake outburst flood in the Santa Cruz Valley (Cordillera Blanca, Perú), Earth Surf. Proc. Land., 43, 1373–1389, 2018.
Monaghan, J. J.: SPH without a tensile instability, J. Comput. Phys., 159, 290–311, 2000.
Nadim, F., Kjekstad, O., Peduzzi, P., Herold, C., and Jaedicke, C.: Global landslide and avalanche hotspots, Landslides, 3, 159–173, 2006.
Necas, J. and Hlavácek, I.: Mathematical theory of elastic and elastoplastic bodies: an introduction, Elsevier, https://doi.org/10.1016/B9780444997548.500084, 2017.
Ngekpe, B., Ode, T., and Eluozo, S.: Application of totalstrain crack model in finite element analysis for punching shear at edge connection, Int. J. Res. Eng. Social Sci., 6, 1–9, 2016.
Julien, P. Y. and O'Brien, J. S.: Selected notes on debris flow dynamics, in: Recent Developments on Debris Flows. Lecture Notes in Earth Sciences, , edited by: Armanini, A. and Michiue, M., vol 64, Springer, Berlin, Heidelberg, https://doi.org/10.1007/BFb0117766, 1997.
Osorno, M. and Steeb, H.: Coupled SPH and Phase Field method for hydraulic fracturing, Proc. Appl. Math. Mech., 17, 533–534, https://doi.org/10.1002/pamm.201710236, 2017.
Pailha, M. and Pouliquen, O.: A twophase flow description of the initiation of underwater granular avalanches, J. Fluid Mech., 633, 115, https://doi.org/10.1017/S0022112009007460, 2009.
Pastor, M., Blanc, T., Haddad, B., Petrone, S., Sanchez Morles, M., Drempetic, V., Issler, D., Crosta, G. B., Cascini, L., Sorbino, G., and Cuomo, S.: Application of a SPH depthintegrated model to landslide runout analysis, Landslides, 11, 793–812, https://doi.org/10.1007/s103460140484y, 2014.
Pitman, E. B. and Le, L.: A twofluid model for avalanche and debris flows, Philos. Trans. A. Math. Phys. Eng. Sci., 363, 1573–1601, https://doi.org/10.1098/rsta.2005.1596, 2005.
Price, N. J.: Fault and joint development: in brittle and semibrittle rock, Elsevier, https://doi.org/https://doi.org/10.1016/C20130054102, 2016.
Pudasaini, S. P.: A general two‐phase debris flow model, J. Geophys. Res.Earth, 117, F03010, https://doi.org/10.1029/2011JF002186, 2012.
Pudasaini, S. P. and Fischer, J.T.: A mechanical erosion model for twophase mass flows, arXiv preprint, arXiv:1610.01806, 2016.
Pudasaini, S. P. and Hutter, K.: Rapid shear flows of dry granular masses down curved and twisted channels, J. Fluid Mech., 495, 193–208, 2003.
Pudasaini, S. P. and Hutter, K.: Avalanche dynamics: dynamics of rapid flows of dense granular avalanches, Springer Science & Business Media, SpringerVerlag Berlin Heidelberg, https://doi.org/10.1007/9783540326878, 2007.
Pudasaini, S. P. and Mergili, M.: A multiphase mass flow model, J. Geophys. Res.Earth, 124, 2920–2942, 2019.
Pudasaini, S. P., Hajra, S. G., Kandel, S., and Khattri, K. B.: Analytical solutions to a nonlinear diffusion–advection equation, ZAMMZ Angew. Math. Phys., 69, 150, https://doi.org/10.1007/s0003301810426, 2018.
Rankine, W. J. M.: II. On the stability of loose earth, Philos. T. R. Soc. Lond., 147, 9–27, 1857.
Reiche, P.: The TorevaBlock: A distinctive landslide type, J. Geol., 45, 538–548, 1937.
Rickenmann, D., Laigle, D., McArdell, B., and Hübl, J.: Comparison of 2D debrisflow simulation models with field events, Comput. Geosci., 10, 241–264, 2006.
Roberts, M.: Evenly Distributing Points in a Triangle, available at: http://extremelearning.com.au/evenlydistributingpointsinatriangle/ (last access: 3 March 2021), 2020.
Savage, S. B. and Hutter, K.: The motion of a finite mass of granular material down a rough incline, J. Fluid Mech., 199, 177–215, 1989.
Saxton, K. E. and Rawls, W. J.: Soil water characteristic estimates by texture and organic matter for hydrologic solutions, Soil Sci. Soc. Am. J., 70, 1569–1578, 2006.
Sheridan, M. F., Stinton, A. J., Patra, A., Pitman, E., Bauer, A., and Nichita, C.: Evaluating Titan2D massflow model using the 1963 Little Tahoma peak avalanches, Mount Rainier, Washington, J. Volcanol. Geotherm. Res., 139, 89–102, 2005.
Spencer, A. J. M.: Continuum mechanics, Courier Corporation, Dover Publications, ISBN13 9780486139470, 2012.
Stead, D. and Wolter, A.: A critical review of rock slope failure mechanisms: the importance of structural geology, J. Struct. Geol., 74, 1–23, 2015.
Steffen, M., Kirby, R. M., and Berzins, M.: Analysis and reduction of quadrature errors in the material point method (MPM), Int. J. Num. Meth. Eng., 76, 922–948, 2008.
Stomakhin, A., Schroeder, C., Chai, L., Teran, J., and Selle, A.: A material point method for snow simulation, ACM Transactions on Graphics (TOG), 32, 1–10, 2013.
Tang, C.L., Hu, J.C., Lin, M.L., Angelier, J., Lu, C.Y., Chan, Y.C., and Chu, H.T.: The Tsaoling landslide triggered by the ChiChi earthquake, Taiwan: insights from a discrete element simulation, Eng. Geol., 106, 1–19, 2009.
Van Asch, T. W., Tang, C., Alkema, D., Zhu, J., and Zhou, W.: An integrated model to assess critical rainfall thresholds for runout distances of debris flows, Nat. Haz., 70, 299–311, 2014.
van den Bout, B.: OpenLISEM Hazard 2.0 alpha, DANS, Faculty of GeoInformation Science and Earth Observation (ITC), University of Twente https://doi.org/10.17026/dansxz42tut, 2020.
van den Bout, B.: openLISEM, available at: https://sourceforge.net/projects/lisem/ (last access: 22 March 2021), SourceForge, 2021.
van den Bout, B and Jetten, V. G.: openLISEM a spatial model for runoff, floods and erosion, available at: https://blog.utwente.nl/lisem/ (last access: 22 March 2021), 2020.
Van Looy, K., Bouma, J., Herbst, M., Koestel, J., Minasny, B., Mishra, U., Montzka, C., Nemes, A., Pachepsky, Y. A., and Padarian, J.: Pedotransfer functions in Earth system science: challenges and perspectives, Rev. Geophys., 55, 1199–1256, 2017.
Varnes, D. J.: Slope movement types and processes, Special Rep., 176, 11–33, 1978.
von Boetticher, A., Turowski, J. M., McArdell, B. W., Rickenmann, D., Hürlimann, M., Scheidl, C., and Kirchner, J. W.: DebrisInterMixing2.3: a finite volume solver for threedimensional debrisflow simulations with two calibration parameters – Part 2: Model validation with experiments, Geosci. Model Dev., 10, 3963–3978, https://doi.org/10.5194/gmd1039632017, 2017.
Williams, J. R.: Application of SPH to coupled fluidsolid problems in the petroleum industry, Videos of Plenary Lectures presented at the IV International Conference on ParticleBased Methods, Fundamentals and Applications (PARTICLES 2015), Technical University of Catalonia (UPC), Barcelona, Spain, 28–30 September 2015, 2019.
Xie, M., Esaki, T., and Cai, M.: GISbased implementation of threedimensional limit equilibrium approach of slope stability, J. Geotech. Geoenviron. Eng., 132, 656–660, 2006.
Xu, F., Zhao, Y., Li, Y., and Kikuchi, M.: Study of numerical and physical fracture with SPH method, Acta Mech. Solida Sinica, 23, 49–56, 2010.
Xu, R., Kang, L., and Tian, H.: A goctree based fast collision detection for largescale particle systems, in: 2012 International Conference on Computer Science and Electronics Engineering, IEEE, 3, 269–273), 2012.
Xu, S. Y., Lawal, A. I., Shamsabadi, A., and Taciroglu, E.: Estimation of static earth pressures for a sloping cohesive backfill using extended Rankine theory with a composite logspiral failure surface, Acta Geotech., 14, 579–594, 2019.
Zhang, L., Zhang, J., Zhang, L., and Tang, W.: Stability analysis of rainfallinduced slope failure: a review, P. I. Civil Eng.Geotec., 164, 299–316, 2011.
Zhou, F., Molinari, J.F., and Ramesh, K.: A cohesive model based fragmentation analysis: effects of strain rate and initial defects distribution, Int. J. Solid. Struct., 42, 5181–5207, 2005.
 Abstract
 Introduction
 A set of mass movement equations incorporating internal structure
 Flume experiments
 Numerical tests
 Discussion
 Conclusions
 Appendix A: List of symbols
 Appendix B: Stress remapping
 Appendix C: Software implementation
 Code and data availability
 Author contributions
 Competing interests
 Review statement
 References
 Abstract
 Introduction
 A set of mass movement equations incorporating internal structure
 Flume experiments
 Numerical tests
 Discussion
 Conclusions
 Appendix A: List of symbols
 Appendix B: Stress remapping
 Appendix C: Software implementation
 Code and data availability
 Author contributions
 Competing interests
 Review statement
 References