## Abstract

The intrusion of debris flow into a river is a common phenomenon that usually occurs in mountainous areas. Debris flow, unlike a landslide, has a complex component that makes its simulation more difficult. In this paper, a two-layer, 2D model is presented for use in the preliminary study of debris flow intruding into a river. In this model, the upper layer is assumed to be a single-phase flow and the lower layer to be a solid–liquid two-phase flow. The drag force between the fluid and solid phases of the debris flow is considered in this model. The basal friction term of the debris flow is considered as the phase-averaged value of the solid shear stresses and fluid shear stresses. A numerical method based on the finite-volume method is proposed to solve the complex model equations involved. Several numerical tests were performed to confirm the feasibility of the proposed numerical method and model. The numerical results yielded a reasonable prediction of debris flow motion and tsunami wave generation and illustrate the complex interplay between the debris flow and tsunami waves. The characteristics of the river have extra effects on the distribution of subaqueous debris flow (e.g. flow velocity) and should not be neglected.

As one of the major natural hazards, debris flow occurs mainly in steep terrain that is covered with loose material. Debris flow is induced by rainfall and indirectly induced by earthquakes. It is precisely because of these environmental influences that debris flow often has huge volumes and travels large distances. Given these characteristics of debris flow, it is common for a river channel in a mountainous area to be blocked by debris flow. This can result in surges and can threaten life and property in low-lying areas. A series of such debris ﬂow events occurred in the town of Yingxiu, transporting huge amounts of sediment to the Minjiang River and forming debris dams, which resulted in the ﬂooding of the newly reconstructed town (Fig. 1a) (Tang *et al*. 2011). Similar examples include the debris flow in the HongChun gully, which was studied by Liao *et al*. (2013); the rainfall-induced debris ﬂows that completely cut off traffic on the Du-wen Highway; and a weir of the Ming River becoming a damming body that diverted the river plug body and flooded dozens of new houses in Yingxiu (Fig. 1b) (Li *et al*. 2014). Given the damage that debris flows can cause, there is an urgent need for a sound theoretical model for and physical understanding of debris flow blockage of a river in a natural environment for prediction of their geomorphological processes and assessment of the risks they pose.

The dynamics of debris flow blocking a river depend on many factors, including the properties of the debris flow and the river, the topography, and the initial and boundary conditions. There is a significant need for reliable methods for predicting such event dynamics. In recent years, two-layer models that reflect the interaction between different layers have been adopted by many researchers to simulate the dynamics of debris flow (Castro *et al*. 2001; Capart & Young 2002; Salmon 2002; Li *et al*. 2013; Kurganov & Miller 2014). In the classical two-layer shallow-water flow model, the two layers can be endowed with distinct densities and velocities but the same rheological properties (Abbott 1979). Assuming that the lower layer is a Bingham fluid, Chen & Peng (2006) proposed a two-layer, 2D, shallow-water computational model and successfully simulated the confluence phenomenon qualitatively. In another work by Chen *et al*. (2007), a two-layer shallow-water model in which the lower layer was modelled as a Herschel–Bulkley fluid was used to simulate mud flow intrusions into quiescent water. Fernandez-Nieto *et al*. (2008) used a two-layer model of the Savage–Hutter type to study submarine avalanches. Their model considers the effects of buoyancy and the centripetal acceleration of the grain movement owing to the curvature of the bottom in the deﬁnition of the Coulomb term. However, these models are more suitable for use in describing single-phase ﬂow. For debris flow, the coexistence and interaction of ﬂuid and solid granular phases within a ﬂowing mass always play key roles in natural gravity-related instabilities.

To reflect the nature of debris flow more realistically, numerical simulations of debris flow based on two-phase models have been conducted (George & Iverson 2011; Pudasaini 2012; Rosatti & Begnudelli 2013; von Boetticher *et al*. 2014; Ouyang *et al*. 2015; Bouchut *et al*. 2016). Iverson & Denlinger (2001) considered the flow to be a mixture of fluid and solid masses. They assumed that the relative velocity between the two constituents was small, and thus that the drag force between the solid and fluid phases could be neglected. However, the solid- and fluid-phase velocities may differ substantially in natural debris flow, and thus the drag force has an important effect on the motion of debris flow. Pitman & Le (2005), on the other hand, accounted for the relative velocity between the solid and fluid phases, assuming the fluid phase to be an ideal fluid. Pudasaini (2012) presented a generalized two-phase debris flow model that considers the effects of many physical factors, such as buoyancy, drag force and virtual mass. In this model, the Mohr–Coulomb plasticity is used to approximate the solid stress, and the fluid stress is modelled as a non-Newtonian viscous stress. Bouchut *et al*. (2015) proposed a thin-layer, depth-averaged, two-phase model that employs a dissipative energy balance to describe avalanches of solid–fluid mixtures. However, study of a model that incorporates the two-phase model of debris flow into the two-layer model to describe the dynamic process of debris flow intruding into a river is still lacking.

The main purpose of this paper is to present a depth-averaged two-layer model for two-phase debris flow intruding into a river that combines a two-phase model with a two-layer model. In this model, the interaction between the solid and fluid phases and the interaction between the upper and lower layers are considered in describing the dynamics of debris flow into a river, from rest to ﬁnal deposition. The effect of the river on the movement of subaqueous debris flow is also studied. Moreover, a well-balanced numerical scheme based on a Harten–Lax–van Leer contact (HLLC) scheme is proposed for use in simulating debris ﬂows into rivers over complex domains involving wetting and drying.

This paper is organized as follows. First, a description of the depth-averaged, two-phase, two-layer model is presented. The details of a computationally efficient numerical approach are then described, and the results of numerical simulations are discussed. Finally, concluding remarks are presented.

## The mathematical model

The intrusion of debris flow into a river, as a complicated dynamic process, is considered to be a typical two-layer problem. As illustrated in Figure 2, debris flow (layer 2) on the bottom of the river (layer 1) is flowing through a fixed river bed, and the two flowing layers are assumed to have distinct densities *ρ*_{1} and *ρ*_{2} that are uniform over the corresponding depths *h*_{1} and *h*_{2}. The free surface *δ*_{1}, layer interface *δ*_{2} and basal surface *δ*_{3} separate the atmosphere from the river, the river from the debris flow and the debris flow from the fixed bed, respectively. To each of these three surfaces corresponds a unit normal vector field, which points upward as depicted in Figure 2: *n*_{1} corresponds to *δ*_{1}, *n*_{2} corresponds to *δ*_{2}, and *n*_{3} corresponds to *δ*_{3}. Here, a locally inclined Cartesian-type coordinate system, 0–*xyz*, is introduced by a 2D reference surface that follows the mean downslope chute topography (Gray *et al*. 1999; Ma *et al*. 2015; Meng & Wang 2015; Bouchut *et al*. 2016). At any location, the *x*-axis is oriented in the downslope direction, the *y*-axis lies in the cross-slope direction and the *z*-axis is an outward-directed vector normal to them. As noted by Gray *et al*. (1999), the downslope inclination angle *θ* of the reference surface is not uniquely defined and can be considered as a function of the downslope coordinate *x*. For this complicated two-layer problem, the continuum mechanics theory and modelling equations are described below.

### The upper layer

Here, we consider the upper layer with the following properties: (1) the upper layer is an incompressible and viscous fluid with constant density *ρ*_{1}; (2) at the free surface *δ*_{1} the mass exchange (e.g. owing to evaporation) is not considered; (3) at the layer interface *δ*_{2} the mass exchange is not considered.

The fundamental balance laws for mass and momentum may be written as follows (Ma *et al*. 2015; Yavari-Ramshe & Ataie-Ashtiani 2015):
(1)where **u**_{1} = (*u*_{1}, *v*_{1}, *w*_{1}) is the velocity field of the upper layer and **P**_{1} is the corresponding stress tensor of the upper layer. These variables are functions of space (*x*, *y*, *z*) and time *t*. **g** is the gravitational acceleration. However, the 3D theoretical model is unmanageable and computationally complex. According to previous studies (Molls & Chaudhry 1995; Wu 2004; Costa & Macedonio 2005; Liang *et al*. 2007; Iverson & George 2014), a depth-averaged approximation is commonly used to simplify the complex governing equations owing to the disparity of the length scales *H*/*L* of river flow in *x* and *z* directions, where *H* is the typical thickness and *L* is the typical topography-parallel length. Moreover, these conservation equations are subject to kinematic and boundary conditions, which express relationships between rates of boundary elevation change and the velocity components of the material adjacent to the boundary at the free surface and the base (Gray 2001; Iverson & Ouyang 2015). The following kinematic conditions are considered:
(2)
(3)where **u**(*δ*_{1}) and **u**(*δ*_{2}) are the velocities of upper layer at the top and bottom boundary, respectively, which range from zero to the average stream velocity, which is dependent on the flow characteristics (Le & Pitman 2009; Johnson *et al*. 2012). *z*(*δ*_{1}) and *z*(*δ*_{2}) denote the position of the top and bottom boundaries of upper layer. A stress-free boundary condition is imposed at the free surface *δ*_{1}:
(4)For the layer surface *δ*_{2}, a friction boundary condition is imposed (Fernandez-Nieto *et al*. 2008):
(5)where **f*** _{in}* is a friction term between the two layers. By depth-averaging the system of conservation equations, and kinematic and boundary conditions (Iverson & Ouyang 2015; Liu & He 2015; Yavari-Ramshe & Ataie-Ashtiani 2015), the model equations of the upper layer can be expressed as
(6)
(7)
(8)where

*g*,

_{x}*g*and

_{y}*g*represent the components of the gravitational acceleration in the reference surface induced downslope, cross-slope and in the normal direction, respectively;

_{z}*z*=

_{m}*h*

_{2}+

*z*denotes the basal surface for the upper flow;

_{b}*z*represents the river basal surface.

_{b}### The lower layer

As a classical two-phase flow, both the solid and fluid phases of a debris flow play an important role in determining its dynamic characteristics. Thus, debris flows are considered in this study to be mixtures of solid and fluid phases with separate evolutions of their volume fractions and velocities. Here, we consider a debris flow with the following assumptions: (1) the solid and ﬂuid components of debris flows are individually incompressible and with constant densities; (2) at the layer surface *δ*_{2} mass exchange is not considered; (3) at the basal interface *δ*_{3} mass exchange is not considered; (4) the fluid phase is viscous; the normal solid stress in the *z*-direction at any height is equal to the (buoyancy) reduced weight of the solid material overburden.

As given by Anderson & Jackson (1967), Pitman & Le (2005) and Pudasaini (2012), the mass and momentum conservation equations for the solid and fluid constituents are
(9)
(10)
(11)
(12)where *ρ _{s}* and

*ρ*are the solid and fluid densities, respectively; the velocities are

_{f}**u**

*= (*

_{s}*u*,

_{s}*v*,

_{s}*w*) for the solid phase and

_{s}**u**

*= (*

_{f}*u*,

_{f}*v*,

_{f}*w*) for the fluid phase;

_{f}*α*and

_{s}*α*= (1–

_{f}*α*) denote the volumetric fractions for the solid and fluid constituents, respectively; the mixture density of debris flow

_{s}*ρ*=

_{2}*α*+

_{s}ρ_{s}*α*;

_{f}ρ_{f}**ζ**

*and*

_{s}**ζ**

*are the partial stresses of the solid phase and fluid phase, respectively;*

_{f}**f**

*represents the sum of the buoyant force*

_{ce}**f**

*and all the nonbuoyant components of the resultant forces*

_{b}**f**

*exerted by the fluid on the solid phase,*

_{m}**f**

*=*

_{ce}**f**

*+*

_{b}**f**

*. As detailed by Jackson (2000) and Pudasaini (2012), the term*

_{m}**f**

*combines the drag, lift and virtual mass forces. Here, the force*

_{m}**f**

*is expressed simply by the drag force*

_{m}**f**

*. According to Anderson & Jackson (1967), the buoyancy force*

_{d}**f**

*is expressed as –*

_{b}*α*∇

_{s}*p*, where

*p*denotes the fluid pressure. For the solid phase, the value of

**ζ**

*relates to the solid volume fraction, and can be expressed as*

_{s}**ζ**

*= –*

_{s}*α*

_{s}**T**

*, where*

_{s}**T**

*is the stress tensor of the dry granular material (Meng & Wang 2015). For the fluid phase, the common relation*

_{s}**ζ**

*= –*

_{f}*p*

**I**+

*α*

_{f}**τ**

*is used, where*

_{f}**τ**

*is the additional stress of the pure fluid and*

_{f}**I**is the unit tensor. According to Pudasaini (2012),

**τ**

*is expressed as viscous-fluid stresses. Finally, by substituting these expressions into (10) and (12), the momentum conservation equations are rewritten in the following forms: (13) (14)The kinematic conditions for the lower layer are given as follows: (15) (16)where*

_{f}**u**

_{(s, f)}(

*δ*

_{2}) and

**u**

_{(s, f)}(

*δ*

_{3}) are the velocities of the lower layer at the top and bottom boundary, respectively.

*z*(

*δ*

_{3}) denotes the position of the bottom boundary of the lower layer. The friction boundary condition is also imposed at the layer surface

*δ*

_{2}for the solid and fluid phase (Bouchut

*et al*. 2016): (17)where

**T**

_{2}=

**ζ**

*+*

_{s}**ζ**

*is the stress tensor of the lower layer. However, the friction between the upper layer and the fluid phase of the lower layer is neglected. Thus, (17) can be expressed as (18)For the solid phase, it also satisfies the Coulomb friction boundary condition at the basal surface*

_{f}*δ*

_{3}(Meng & Wang 2015): (19)where

*φ*is the Coulomb friction angle of the basal surface. For the fluid phase, considering the interaction between the fluid phase and bed surface, a friction law is imposed on the bottom (Benkhaldoun

_{bed}*et al*. 2011; Li & Duffy 2011): (20)where

*n*is the Manning friction coefficient. It should be noted that any grid cell with

_{b}*h*

_{2}< 10

^{−5}m is assumed to be filled with only the upper layer and the basal friction of the lower layer is set to zero. Moreover, because of the complexity of the components of debris flow, the basal friction term of debris flow is considered to be the phase-averaged value of the solid shear stresses and fluid shear stresses (as illustrated in Fig. 3). The phase-averaged basal friction term

**T**

*can be written as follows: (21)where*

_{zb}**T**

*and*

_{szb}**T**

*are the basal shear stresses of the solid and fluid phases, and are described by (19) and (20), respectively. Similarly to the deduced process of the upper layer equations, the depth-averaged procedure is also applied to simplify the lower layer equations. The analogous derivation of these equations has been proposed in earlier papers (Pitman & Le 2005; Hutter & Luca 2012; Pudasaini 2012; Meng & Wang 2015; Bouchut*

_{fzb}*et al*. 2016). By expressing quantities containing the variables

*α*,

_{s}*α*and

_{f}*h*in terms of the conserved quantities

_{2}*h*=

_{s}*α*and

_{s}h_{2}*h*=

_{f}*α*, the mass-balance equations for the solid and fluid phases of the lower layer are as follows: (22) (23)where

_{f}h_{2}*h*and

_{s}*h*indicate the depth of solid and fluid phases in the lower layer, respectively. Neglecting the shear stress of the solid phase because it is much lower than the normal stress, the depth-averaged momentum conservation equations for the solid and the fluid phases can be expressed as follows: (24) (25) (26) (27)where

_{f}*r*=

*ρ*/

_{s}*ρ*is the density ratio,

_{f}*ε*=

*H*/

*L*is the aspect ratio of debris flow and

*k*

_{(x,y)}represents the lateral earth pressure coefficients. A formula proposed by Gray

*et al*. (1999) was adopted to illustrate an active (or passive) state of stress if an element of material is elongated (or compressed): (28) (29)

*k*

_{(x,y)}can be defined as a function of the internal and basal angles of friction: (30)where

*φ*is the internal friction angle of the flow material; ‘−’ corresponds to an active state, and ‘+’ corresponds to a passive state.

_{int}### The interaction force

The interaction force between two layers is considered to depend primarily on the velocity difference between two layers (Fernandez-Nieto *et al*. 2008) and can be expressed as follows:
(31)where *χ* is a positive constant. For the lower layer, various empirical relations have been proposed in the literature for the drag force between the solid and fluid phases (Pitman & Le 2005; Pelanti *et al*. 2008; Pudasaini 2012; Bouchut *et al*. 2015). In this study, the following simple drag interaction formula was adopted (Pitman & Le 2005):
(32)where *v _{T}* is the terminal velocity of a typical solid particle falling in the ﬂuid under gravity and

*m*is related to the Reynolds number of the flow and varies from

*c.*4.25 to 2.4 as the Reynolds number increases from zero to inﬁnity (Bouchut

*et al*. 2015).

## Computational schemes

The proposed numerical program employs a Godunov-type finite-volume scheme to discretize the model equations. After the discretization, an HLLC Riemann approximate solver is applied to manage the discontinuous problem at the grid cell interface. To improve the calculation accuracy of the HLLC scheme, the MUSCL (monotonic upstream-centred scheme for conservation laws) approach is used to reconstruct the conserved variables of the model before calculation. Moreover, a well-balanced numerical scheme is also used to handle the reconstructed variables before calculation to adapt the complex domains involving wetting and drying. Because the model equations are divided into two 1D sub-equations by using an operator-splitting technique, each sub-equation is calculated twice in a consecutive way to improve time accuracy. With this numerical program, second-order accuracy both in time and space can be achieved. Further details are described in the following sections.

In matrix form, the conservation law of the 2D, two-phase, two-layer model equations can be written as follows:
(33)where **w**, **f**, **m**, **s** and **q** are vectors representing the conserved variables, the fluxes in the *x* and *y* directions, and the source terms in the *x* and *y* directions, respectively. The explicit time-marching conservative finite-volume form of (33) is as follows:
(34)where the superscript *k* is the time level, the subscript *i* is the cell index and Δ*t* is the time step. As shown in Figure 4, **f*** _{et}* and

**f**

*are the fluxes through the western and eastern cell interfaces, respectively;*

_{wt}**m**

*and*

_{nt}**m**

*are the fluxes through the northern and southern cell interfaces, respectively; Δ*

_{st}*x*and Δ

*y*are the side lengths of a grid cell in the

*x*and

*y*directions, respectively. To simplify calculation, two separate 1D problems in the

*x*and

*y*directions can be obtained by dividing (34) based on the operator-splitting technique (Liang

*et al*. 2006; Toro 2009) and are expressed as follows: (35) (36)Then, the solution at the next time step can be obtained by an efficient step to obtain the second-order time accuracy (Ouyang

*et al*. 2013) as follows: (37)Here, we assume that

*L*and

_{x}*L*represent the computational procedures and are used for solving (35) and (36) in the

_{y}*x*and

*y*directions, respectively. Each computational procedure is composed of a Godunov-type finite-volume scheme combined with an HLLC approximate Riemann solver and MUSCL approach. However, as the properties of the different layers have shown significant differences, the equations for the upper layer and lower layer are solved separately but with the same computational procedure at the same time step. For

*L*, the internal flux (e.g.

_{x}**f**

*) is computed as follows: (38)where*

_{wt}**f**

*and*

_{l}**f**

*are the interface fluxes on both sides of the east cell interface, which are calculated from the left and right Riemann states*

_{r}**w**

*and*

_{l}**w**

*;*

_{r}*S*,

_{l}*S*and

_{m}*S*represent the speeds of the left, middle and right waves, respectively, for a local Riemann problem; and

_{r}**f**

*and*

_{*l}**f**

*represent the left and right sides of the contact wave, respectively. For the upper layer,*

_{*r}**f**

*and*

_{*l}**f**

*can be expressed as follows: (39)where*

_{*r}*v*and

_{l}*v*are the left and right tangential velocity components of the Riemann states. Coupling the MUSCL approach with the HLLC scheme allows us to compute the interface functions

_{r}**w**

*and*

_{l}**w**

*to ensure accuracy and avoid spurious oscillations (van Leer 1979; Zoppou & Roberts 2000). The*

_{r}**w**

*and*

_{l}**w**

*functions can be expressed as follows: (40)where (41)The function M is a min–mod flux limiter that can be written as follows: (42)Then, the fluxes*

_{r}**f**

_{*}in the middle region are needed to calculate

**f**

*and*

_{*l}**f**

*by using the constructed variables and can be obtained from the Harten–Lax–van Leer (HLL) formula (Harten*

_{*r}*et al*. 1983): (43)Considering the dry-bed condition from the two-rarefaction approximate Riemann solver, the wave speeds are calculated as follows (Fraccarollo & Toro 1995; Soares-Frazão & Zech 2011; Balsara 2012): (44) (45) (46)where

*c*is the speed of gravity waves;

*u*,

_{l}*u*,

_{r}*h*,

_{l}*h*are the components of the left and right Riemann states for a local Riemann problem;

_{r}*u*

_{*}and

*h*

_{*}are the components of the middle Riemann states, which are calculated as follows: (47)It should be noted that the difference of calculation process between the two layers is mainly reflected by the speed of gravity waves, which govern ﬂuxes across cell edges in space and time. For the upper layer, the speed of gravity waves

*c*

_{1}is the standard shallow water expression,

*c*

_{1}

^{2}=

*g*

_{z}h_{1}. For the lower layer, the speeds of gravity waves for solid and fluid phases,

*c*and

_{s}*c*, can be given by (48)The form of

_{f}*c*indicates that the speed of gravity wave for the solid phase is influenced not only by lateral stress transfer but also by the fluid phase. Compared with the standard expression for the speed of gravity wave for shallow water, an extra term in the expression for

_{s}*c*has also been influenced by the fluid phase. Moreover, a well-balanced numerical scheme proposed by Liang & Borthwick (2009) was used for simulating frictional shallow ﬂows over complex domains involving wetting and drying. However, because the corresponding basal surface is different for each layer, this scheme should be applied to each flux computation of two layers. Similarly, this procedure is also applied to

_{s}*L*to solve (36). A completed computational flow diagram is shown in Figure 5. The solution at the next time step can be obtained by operating the procedures each twice with the time step calculated by the stability criterion (Brufau

_{y}*et al*. 2004; Simpson & Castelltort 2006). The time step should satisfy the demand for dynamic computing of two layers simultaneously: (49)where

*cfl*is the Courant number, the value of which should be less than unity, and

*η*is the ratio of the area of the grid to its perimeter.

## Test results and discussion

Some numerical simulations were performed to confirm the feasibility of the numerical scheme and model presented in this paper. In all of the simulations, we set *g* = 9.8. To reduce the computational complexity, the density of the upper layer was set to be the same as the density of the fluid phase of the debris flow. A small value was used for the Courant number, *cfl* = 0.5, to ensure the stability of the computation.

### A tidal wave flow over steps

A tidal wave flow over a bed with two vertical steps was considered to validate the capability of the numerical scheme to deal with discontinuous bed topography. This case was initially proposed by the EU CADAM project, and was later investigated by Zhou *et al*. (2002). The channel is 1500 m long and frictionless, and the bed profile is defined by
(50)Bermudez & Vazquez (1994) provided the asymptotic analytical solution of the flow depth and velocity as follows:
(51)The initial flow depth and velocity were *h*(*x*,0) and *u*(*x*,0). During the simulation, the flow depth and velocity of the left end of the channel changed over time, as expressed by *h*(0,*t*), to drive the flow. A closed-boundary wall is imposed at the right end of the channel as *u*(*L*,*t*) = 0. Comparisons of the numerical predictions and analytical solutions for *t* = 10800 s and *t* = 32400 s are shown in Figure 6. Good agreement is observed and demonstrates that the well-balanced numerical scheme is able to describe a debris ﬂow into a river over complex and even discontinuous bed topography.

### Subaqueous debris flow

In this subsection, a simulation of a subaqueous debris flow in an inclined channel was performed. This numerical test was initially proposed by Fernandez-Nieto *et al*. (2008) and further used by Kurganov & Miller (2014), and is used here to verify the feasibility of the presented model. The initial conﬁguration is shown in Figure 7a. The channel was 10 m long, and the bottom *z _{b}* and the initial data are given by
(52)To compare with the numerical results from Kurganov & Miller (2014), first the flow is set as pure solid, and the interlayer friction term and the lateral earth pressure are also neglected, as

**f**

*= 0 and*

_{in}*k*= 1. The Coulomb friction angle

_{x}*φ*= 25°, and an open boundary condition is imposed on each side of the channel. The values of other relevant parameters are shown in Table 1. The computed results at times

_{bed}*t*= 0.1, 0.5, 1, 1.5 and 5 s with Δ

*x*= 0.05 are shown in Figure 7b–f. Once the mass flow is released from rest on the inclined channel, it spreads out in both directions and leads to the generation of surface waves that deform together with the flow shape. Although there is a slight difference between the two results, which may be caused by the numerical scheme (e.g. the treatment of the wet and dry fronts), the change trends of the flow shape and wave propagation are in good agreement. However, there are some significant differences of flow characteristics between pure solid flow and debris flow, mainly focusing on the presence of a fluid phase and the drag force between the solid and fluid phases. To investigate these differences, the final stationary interface of debris flow that we obtain for three solid volume fractions with

*α*= 0.8,

_{s}*α*= 0.5 and

_{s}*α*= 0.3 are shown in Figure 8. The values of the other parameters are the same as before. It can be found that the larger the value of solid volume fraction, the higher is the mobility of subaqueous debris flow. This change rule is contrary to that for subaerial debris flow, and can be considered as the result of the interaction between the upper water layer and the subaqueous debris flow.

_{s}### The intrusion of two-phase debris flow into a river

The model was next used to simulate a 2D two-phase debris flow intruding into a river to demonstrate the potential of the model to simulate a practical case. A simple sketch of the flume geometry is shown in Figure 9. The flume we chose for this case was a straight, rectangular channel 35 m in length and 8 m in width. Two portions of it were sloping (*x* < 20 m, *θ* = 11.3° and *x* > 30 m, *θ* = 45°) and the rest was flat. Flowing water with a depth of 1 m made contact with the slope at *x* = 15 m and *x* = 31 m. The debris flow had a volume of 2.51 m^{3}, initially contained in a Perspex cap 2 × 2 × 0.4 m in size, and was released from rest on the middle of the slope at *y* = 4 m (centre line). The debris flow consisted of 60% solid and 40% fluid. Here, we assumed that the river has a flow velocity only in the *y* direction, with *v*_{1} = 5 m s^{−1}. The values of the relevant parameters are shown in Table 2. An open boundary condition was imposed on each side of the flume. Figure 10 shows the numerical results for the debris flow motion and tsunami wave generation at times *t* = 1, 3, 4, 6 and 8 s after release of the debris flow. In Figure 10 the upper subpanels represent the debris flow shape and the lower subpanels represent the evolution of the river surface mainly shown at *x* = 15 m and *x* = 31 m. Once the cap was open, the debris flow accelerated and spread out rapidly in the downslope direction at *t* = 1 s. As soon as the debris flow intruded the river, a wave was generated in the same direction as the debris flow at *t* = 3 s. However, differing from still water, except for the water resistance the river provides an extra force, which affects the subaqueous debris flow in the flow direction. Because of the effects of water resistance and this extra force, the front of the debris flow started to accumulate in the upstream part of the river and spread in the downstream part of the river. At *t* = 4 s, the debris flow continued to move into the river in a downslope direction, and part of the subaqueous debris flow had moved along the direction of river flow. In addition, the generated waves also propagated in the downslope direction and downstream. At *t* = 6 s, the front of the debris flow accumulated in the downslope direction, but the spread of subaqueous debris still continued in the flow direction of the river, owing to the persistent effect of the river. As the debris completely flowed into the water at *t* = 8 s, part of it was washed away by the river along the flow direction of river, and a deviation in the distribution of subaqueous debris flow in the flow direction of the river and in the opposite direction was found. It also should be noted that the disturbance of waves was mainly focused on the location of debris flow accumulation. The general patterns of the debris flow motion and associated tsunami waves showed that the river has an extra effect on the distribution of subaqueous debris flow caused by the flow velocity. The area of influence of the debris flow intruding into the river was mainly focused on the downstream part of the river. The numerical results also demonstrate the applicability of the model.

## Conclusions

In this paper, a two-layer model for simulation of the intrusion of a two-phase debris flow into a river is presented. Unlike the common two-layer model, the lower layer of this model is considered to be a solid–liquid two-phase flow based on the complex component of debris flow. The depth-averaged equations of the model are obtained by assuming that the characteristic horizontal length scale of the debris flow is greater than the characteristic vertical length scale. Based on mass and momentum conservation, the interaction force between the two layers and the drag force between the solid and fluid phases are considered to make the two layers behave as a whole. In addition, based on the characteristics of the debris flow, the basal friction term is considered to be the phase-averaged value of the solid shear stresses and fluid shear stresses. To solve the complex and hyperbolic model equations, we also propose a highly accurate computational method based on the finite-volume method and coupled with a well-balanced numerical scheme. Several numerical simulations were conducted to investigate the feasibility of the proposed model. Comparisons of the analytical solutions and numerical results confirm the validity of the proposed model. The numerical results show that the larger the value of solid volume fraction, the higher is the mobility of subaqueous debris flow. It can be considered as the result of the interaction between the upper layer and the subaqueous debris flow. Moreover, the river has an extra effect on the distribution of subaqueous debris flow caused by the flow velocity, and should not be neglected.

Although the results presented in this study demonstrate the applicability of the proposed model to simulation of the intrusion of debris flow into river, further improvements to the model are needed. Some additional factors that may influence the dynamics of debris flow, such as erosion and solid–liquid separation, need to be considered. The model also needs to be validated by comparing its predictions with observations of the coupling of debris flow and tsunami wave generation.

## Acknowledgements

The authors thank the anonymous reviewers for helpful suggestions.

## Funding

Financial support from the NSFC (Grant Nos 41272346 and 41101008) and the Science and Technology support plan of Sichuan province (2016SZ0067) is acknowledged.

- © 2018 The Author(s)