Natural Hazards and Earth System Sciences Modelling shallow debris flows of the Coulomb-mixture type over temporally varying topography

We propose a saturated binary mixture model for debris flows of the Coulomb-mixture type over temporally varying topography, where the effects of erosion and deposition are considered. Due to the deposition or erosion processes, the interface between the moving material and the stagnant base is a non-material singular surface. The motion of this singular surface is determined by the mass exchange between the flowing layer and the ground. The ratio of the relative velocity between the two constituents to the velocity of the solid phase is assumed to be small, so that the governing equations can be reduced to a system of the quasisingle-phase type. A shock-capturing numerical scheme is implemented to solve the derived equation system. The deposition shapes of a finite mass sliding down an inclined planary chute are investigated for a range of mixture ratios. The geometric evolution of the deposition is presented, which allows the possibility of mimicking the development of levee deposition.


Introduction
Debris flows are collections of geological materials (e.g.sand, cobbles, pebbles or even rocks) mixed with water or slurry, which are driven by gravity to slide over complex terrain.Such flows often cause heavy casualties and the destruction of houses and other facilities.There is considerable interest in predicting what the endangered areas are and in forecasting their evolution and behaviour in order to quantify the risks.A comprehensive description of both theoretical and applied aspects of debris flow is available in Takahashi (2007).Efforts for the theoretical modelling of the movement of debris flows have mainly concentrated either on the rheology of the flowing material or on the influences due to the geometry of the basal surface.The present study is devoted to the description of the movement of shallow debris flows over a non-trivial topography, with the processes of entrainment and deposition being integrated into the model.
One of the pioneering trials describing geophysical flows over real complex terrain with a depth-integrated equation system can be dated back to Savage and Hutter (1989), of which a review together with its extensions are given in Pudasaini and Hutter (2007).Pudasaini and Hutter (2003) and Pudasaini et al. (2005b) designed an orthogonal coordinate system, which is suitable for the bed surface of chutes or of a channel-type form.The Particle-Image-Velocimetry (PIV) technique has been adpoted for velocity measurement and the results compared with theoretical predictions is available in Pudasaini et al. (2005a).In their coordinate system, the curvature and the twisting of the channel are included.Bouchut and Westdickenberg (2004) (BW), on the other hand, presented an arbitrary coordinate system for modelling gravity-driven flows of the shallow-water type over general topographic surfaces.Although in BW's coordinate system explicit expressions of basis vectors and curved coordinates along the locally varying basis vectors are only available for special cases, mostly for simple surfaces or locally defined surfaces, it is desirable to have some flexibility in the choice of (possibly curvilinear) coordinates, and is very helpful to describe the local topography for the purpose of numerical simulations.In a similar approach, De Toni and Scotton (2005) proposed a formulation for granular flows over a nontrivial terrain topography.Based on BW's work, Luca et al. (2009b) derived non-Cartesian, topography-based avalanche equations for the gravity-driven flows of ideal and viscous fluids.A thorough discussion of the hierarchy of avalanche models with an arbitrary topography can be found in Luca et al. (2009a).
Published by Copernicus Publications on behalf of the European Geosciences Union.
The aforementioned models are all related to single-phase granular flows.However, since debris flows are often triggered by excessive rainfall or initiated due to excessive groundwater, a more realistic description should include solid-fluid mixtures, as described by Iverson (1997), Iverson and Denlinger (2001), Pudasaini et al. (2005c), Pitman and Le (2005), Pelanti et al. (2008), Luca et al. (2009b) and more recently by Pudasaini (2011) and Luca et al. (2012).Iverson (1997) and Iverson and Denlinger (2001) derived the mass and momentum balance equations by treating the debris flows as a saturated mixture (Coulomb-Mixture theory).Based on the well-documented, grain-fluid measurements, they assumed the ratio of the relative velocity between the two constituents to the solid velocity to be rather small.As a consequence, the force of the interaction between the constituents is annihilated in the equation system.Pudasaini et al. (2005c) proposed a model based on the same concept, but in an orthogonal curved and twisted coordinate system.Notably in Pudasaini et al. (2005c), levee formation was first observed in their numerical simulation.Despite the fact that interactions between the phases have been largely omitted in the previous work, they are generally associated with the constituent pressures and play an important role.Taking into account the effects due to the non-vanishing relative velocity, Pitman and Le (2005) derived a model for a binary mixture, in which each constituent moves according to its own mass and momentum equations.A slightly modified version of Pitman and Le (2005) can be found in Pelanti et al. (2008).Along the same line, Luca et al. (2009b) formulated model equations in the arbitrary coordinate system proposed by Bouchut and Westdickenberg (2004).In Pudasaini (2011), the Mohr-Coulomb plasticity and the solid volume fraction gradient enhanced non-Newtonian viscous stresses were respectively applied to model the solid stress and fluid extra stress, in which the interaction between the two constituents is taken into account.The generalized interfacial momentum transfer is modeled by including the viscous drag, buoyancy force, and the relative acceleration between the solid and the fluid particles (the virtual mass).The generalized drag force covers both the solid-like and fluidlike contributions.The model equations are written in a multi-curvature, curvilinear coordinate system that fits the basal topography.Beyond the fully developed mixture models, Luca et al. (2012) considered a shallow two-layer debris flow with clean water in the upper layer.They derived a depth-integrated equation system for a three-dimensional topography.
Taking into account the mass exchange between the moving layer and static base, Bouchut et al. (2008) proposed a one-dimensional model for a general topography, in which the basal surface is time dependent.The formulism of the unified coordinate (UC) system (see e.g.Hui, 2004Hui, , 2007) ) allows one to construct a deformable mesh system that moves coincidentally with pseudo-particles on the basal surface.Taking advantage of the introduction of the UC system and the concept of general topography in Bouchut and Westdickenberg (2004), Tai and Kuo (2008) proposed a singlephase model with a deformable coordinate system for twodimensional mountain slopes, where the coordinate ζ = 0 always coincides with the varying topographic surface.A single-phase model for a three-dimensional topography was proposed by Tai et al. (2012), and it is capable of successfully resolving the real catastrophic Hsiaolin landslide (Taiwan), see Kuo et al. (2011).The present study is a Coulombmixture extension of Tai et al. (2012)'s work with temporally varying coordinates for the entrainment-deposition processes.Through the numerical investigation of a finite mass flow, the key features are illustrated in which the geometric evolution of the deposition as well as the levee formation are presented.
The rest of the paper is constructed as follows: Sect. 2 summarizes the derivation of the governing equations where the terrain-following coordinate system, the simplified balance equations, the boundary conditions, the manipulation of the stresses, and the erosion-deposition rate are briefly discussed.The numerical investigation of finite mass flows leading to the formation of the deposition heaps is described in Sect.3. Concluding remarks noting the achievements and inferences for further work are recapitulated in Sect. 4.

Terrain-following coordinate system
The modern geographic information system (GIS) delivers digital representations of ground surface topography.These data can be interpolated as explicit descriptions of non-trivial terrains in horizontal-vertical, Cartesian coordinates.The unit normal vector of a topographic surface can thus be defined.Letting the x-and y-axes lie on the horizontal plane and the z-axis point upward in the vertical direction, the topographic surface S can be given by z − z b (x,y,t) = 0. Following Bouchut and Westdickenberg (2004), the unit normal vector is determined by where ∂ x z b and ∂ y z b are the partial derivatives with respect to the horizontal coordinates and c = 1 represents a point on the topographic surface S, one can define an arbitrary coordinate system O ξ ηζ by allowing to be the natural basis of the tangent space to S, and can be decomposed as (see Fig. 1) (4) Thus, the tangential vectors of the local natural bases read It is straightforward to prove that Luca et al. (2009b).As mentioned in Bouchut and Westdickenberg (2004), this system is an arbitrary coordinate system, so that the coordinates are generally nonorthogonal.The readers are referred to Pudasaini et al. (2003) for a systematic application of non-orthogonal coordinates in avalanche flow modeling.
The unified coordinate method (e.g., Hui, 2004Hui, , 2007) introduces a transformation from the Cartesian coordinates (x,y,z) to some arbitrary coordinates (ξ,η,ζ), or in the component form is the velocity of the local coordinate (mesh) and Ω is the Jacobian matrix of the coordinate transformation.With definition (4), the Jacobian matrix is given by where With the small curvature and shallow flow depth assumptions, i.e. ζ∂ x s = O(ǫ 1+α ), the Jacobian matrix and the tangential vectors of local natural basis reduce to and respectively, where ǫ = H/L ≪ 1 is the aspect ratio of shallow flows and α ∈ (0,1).Approximations ( 9) and ( 10) greatly reduce the complexity of the coordinate transformation for deriving the model equations; see Sect.2.4.

Simplified equations of mass and momentum conservation
We consider here saturated binary mixtures.The true densities of both of the constituents are assumed to be constant, i.e. ρs =const.and ρf =const., and their partial densities in the mixture are determined by their corresponding volume fractions, i.e., the density of the mixture ρ = ρ s + ρ f = ν s ρs + ν f ρf , and the volume fractions are subjected to ν s + ν f = 1.
With the assumption of no mass exchange between the constituents, the balance equations of mass and momentum read for the solid constituent and for the fluid constituent.In this notation g is the gravitational acceleration, and v s/f and Ts/f represent the velocity and stress tensor of the solid/fluid constituent, respectively.In (11) 2 and (12) 2 , ι denotes the interaction force between the two constituents.With the help of the condition ν s + ν f = 1 and the assumption of constant true densities, (11) 1 and (12) 1 merge to become Following Pudasaini et al. (2005c), with the introduction of the fluid specific discharge q = ν f (v f − v s ) the momentum equation of the fluid phase (12) 2 can be recast to to be the unit normal vector.With the definitions in Eqs. ( 2) and (3), one can easily prove that τ ζ = n.For a point at a distance ζ above the topographic surface, its position vector, r, can be decomposed as (see Fig. 1) Thus, the tangential vectors of the local natural bases read It is straightforward to prove that g ζ ≡ g ξ × g η /||g ξ × g η || = n, see e.g.Luca et al. (2009b).As mentioned in Bouchut and Westdickenberg (2004), this system is an arbitrary coordinate system, so that the coordinates are generally non-orthogonal.The readers are referred to Pudasaini et al. (2003) for a systematic application of non-orthogonal coordinates in avalanche flow modeling.
The unified coordinate method (e.g.Hui, 2004Hui, , 2007) ) introduces a transformation from the Cartesian coordinates (x,y,z) to some arbitrary coordinates (ξ,η,ζ ), or in the component form where Q = Q x êx + Q y êy + Q z êz is the velocity of the local coordinate (mesh) and ˜ is the Jacobian matrix of the coordinate transformation.With definition (4), the Jacobian matrix is given by where I d is the 2 × 2 unit matrix, and respectively, where = H /L 1 is the aspect ratio of shallow flows and α ∈ (0,1).Approximations Eqs. ( 9) and ( 10) greatly reduce the complexity of the coordinate transformation for deriving the model equations; see Sect.2.4.

Simplified equations of mass and momentum conservation
We consider here saturated binary mixtures.The true densities of both of the constituents are assumed to be constant, i.e. ρs =const.and ρf =const., and their partial densities in the mixture are determined by their corresponding volume fractions, i.e. the density of the mixture ρ = ρ s + ρ f = ν s ρs + ν f ρf , and the volume fractions are subjected to ν s + ν f = 1.With the assumption of no mass exchange between the constituents, the balance equations of mass and momentum read for the solid constituent and for the fluid constituent.In this notation, g is the gravitational acceleration, and v s/f and Ts/f represent the velocity and stress tensor of the solid/fluid constituent, respectively.In Eqs.(11) 2 and (12) 2 , ι denotes the interaction force between the two constituents.With the help of the condition ν s + ν f = 1 and the assumption of constant true densities, (11) 1 and (12) 1 merge to become Following Pudasaini et al. (2005c), with the introduction of the fluid specific discharge q = ν f (v f − v s ) the momentum equation of the fluid phase (12) 2 can be recast to In Iverson (1997), it is justified that As shown in Pudasaini et al. (2005c), with this condition and after some algebraic manipulation, equations ( 13) and the addition of (11) 2 and ( 14) form an equation system of the quasi-single-phase, where ρ = ρ s + ρ f is the density of the mixture and Ttotal = Ts + Tf indicates the total stress in the mixture.In another words, in this simplified theory, the assumption Henceforth, for brevity, we will use v instead of v s unless otherwise stated.

Boundary conditions and mixture behaviour laws
The flowing layer is bounded by the upper free surface and the lower basal surface.The free surface is assumed to be a material surface for both constituents, so that there is no mass exchange across it.If the free surface is described by F h = 0, the kinematic boundary condition reads At the basal surface, erosion/deposition is considered.If F b = 0 represents the basal surface with a velocity v int , the kinematic boundary condition reads Choosing the unit normal vector n pointing into the flowing layer, the erosion-deposition function (volume flux) can be defined by Here, it is worth noting that E G = E L only occurs when v • n = 0, i.e. there is pure sliding at the basal surface.This also indicates that the densities above and below the basal interface must be different if v • n = 0 (see e.g.Tai and Kuo, 2008).To obtain a terrain-following coordinate system, we choose the local coordinate velocity Q = (v int • n)n to move coincidentally with erosion or deposition, i.e.
We consider the Coulomb-mixture approach proposed by Iverson (1997), Iverson and Denlinger (2001) or Pudasaini et al. (2005c).The total mixture stress is assumed to be Ttotal = Ts + Tf = Te + pI + ν f Tf vis , where Te is the inter-granular effective stress, p is the pore fluid pressure, and Tf vis is the deviator stress due to the viscous interstitial fluid, and I is the identity tensor.In the cited works, the fluid stress is defined as Tf = −pI + ν f Tf vis and the solid stress is postulated to be Ts = Te .Iverson and Denlinger (2001) have argued that a linear variation of fluid pressure is appropriate and suggested a parameter, f ∈ (0,1), to express the fluid pressure as a fraction of the total basal normal stress, i.e.
The linear variation of fluid pressure implies that the parameter f remains constant in the depth-wise direction.Pudasaini et al. (2005c) mentioned that f is a continuous parameter that depends on factors such as the flow thickness and the time and diffusion of the basal pore pressure in the stream-wise direction of the flow body.The linear distribution assumption of the pore fluid pressure and the constant value of f in the flow depthwise direction suggest that the solid and fluid stresses are in the form of throughout the flow depth.The interstitial fluid is assumed to obey the conventional law of incompressible Newtonian fluids, and the fluid pressure is expressed by a fraction f of the total normal stress.As a consequence, the viscous stress of the fluid constituent is thus proposed to be where µ f is the viscosity of the pore fluid and D = 1 2 (∇v + ∇v T ) is the stretching tensor.
The solid constituent is assumed to follow Mohr-Coulomb behaviour: The granular material is an incompressible continuum and the slip plane appears inside the bulk of the material when the internal state of stress equals the yield stress, i.e. the yield occurs when ||t t || = tanφ||t n ||, where t t and t n are the shear and normal stresses, respectively, acting on a planar element inside the granular material, and φ is the internal friction angle of the medium.Some other variants of the Mohr-Coulomb flow law can be found in (e.g.Luca et al., 2009a;De Toni and Scotton, 2005;Greve et al., 1994 et al., 1999;Iverson and Denlinger, 2001;Pudasaini and Hutter, 2003;Gray et al., 2003).In the present study, we would like to exempt ourselves from the complexity of the Mohr-Coulomb constitutives and to illustrate the key features of the model.Thus, we focus on a simpler isotropic (hydraulic) model for the solid constituent.
For the dynamic boundary conditions, the flow surface is assumed to be traction free for both of the constituents, i.e.
Ts n h = 0 and Tf n h = 0, where n h is the unit normal vector of the flow surface.At the basal surface, the Coulomb friction law is applied to the solid constituent, viz.
with δ being the angle of basal friction and n the unit normal vector pointing into the flow body.For the fluid constituent, the no-slip condition is applied at the basal surface.

UC formulation and leading-order model equations
With the unified coordinates, there are two coordinate systems, and every vectorial quantity can be expressed using either systems.For index operation in these coordinate systems, we use i,j ∈ {x,y,z} and m, n ∈ {ξ,η,ζ } for the vectorial components unless otherwise specified.For instance, we have T = T ij êi ⊗ êj = T nm g n ⊗g m = im êi ⊗g m for the same stress tensor expressed with different dyadic combinations of the Cartesian basis êi and the UC basis g n .With the notations they are related by see e.g.Tai et al. (2012).Following Tai and Kuo (2008) or Tai et al. (2012), the simplified quasi-single-phase equation system (15) can be rewritten in the UC formulation, for i ∈ {x,y,z} and m ∈ {ξ,η,ζ }, where 25) are essentially the evolution of J and J v i with respect to the topography-fitted mesh system.The main benefit of this formulation is the uncomplicated computation of the time derivative ∂ t (J v i ), because the velocity v i is expressed in the Cartesian coordinates whose orientations are fixed in space.
The equation system ( 25) is further exploited with scaling analysis to isolate the physically non-significant terms in the governing equations.In the scaling analysis, we first define the aspect ratio as = H /L 1, which is the ratio of a typical flow thickness H to a typical length-scale L on the tangent space related to the basal surface.The variables are non-dimensionalized in a similar way to that done in Gray et al. (1999), Pudasaini et al. (2005b), Tai and Kuo (2008) or Tai et al. (2012).The non-dimensionalization parameters are defined as follows: where ϑ 0 = O( β ) represents the typical magnitude of the friction coefficient and 0 < β < 1.Here, ∂ ξ n and ∂ η n respectively denote the curvatures in the ξ -and η-directions.The curvatures are characterised by the inverse of the curvature radius R. With the small curvature assumption K = L/R = O( α ) and 0 < α < 1, one obtains ζ ∂ x s = O( 1+α ), which leads to see Tai and Kuo (2008) or Tai et al. (2012) for details.The viscous effect of the pore fluid is non-negligible in the present study and it is non-dimensionalized by The velocity component (v x , v y , v z ), the erosion-deposition rate E G/L and the stress tensor ( im ) are scaled by With the help of scalings Eqs. ( 26)-( 28), integrating the equations in (25) from the base ζ = 0 to the free surface ζ = h yields the depth-averaged mass conservation and the depth-integrated momentum equation where the Leibniz rule, the kinematic boundary conditions ( 16)-( 18), the traction free condition at the free surface, and the local coordinate velocity ( 19) have been applied (see also Tai et al., 2012).The over-lines denote the depth-averaged values of the over-lined physical quantities and the subscript "b" indicates the physical quantities at the basal surface.
The velocity component v i can be decomposed into two parts, , where v i represents the velocity component tangential to the topographic surface and the other indicates the velocity component normal to the topographic surface (for details see Tai et al., 2012).With this decomposition, one obtains the depth-averaged value where the coordinate transformation is defined by The velocity components on the right hand side, v m 0 = v m + O( 1+γ ), m ∈ {ξ,η} and γ = min(α,β), are the velocity in the tangential space of the terrain-fitted UC coordinates and the velocity component normal to the topographic surface v ζ 0 vanishes.By virtue of the introduction of the representative mean values, v i 0 and v m 0 , the advection terms in (30) can be approximated by where is the momentum correction factor, i ∈ {x,y,z} and m ∈ {ξ,η} (see Tai et al., 2012).
Multiplying the momentum Eq. ( 15) 2 by ˜ and after some algebraic manipulations, the ζ -component reduces to in the dimensionless form, whereas the second term is the contribution of the local curvature.With traction free boundary condition, Eq. ( 34) can be depth-integrated and yields the total normal stress on the basal surface and the depth-averaged total normal stress The stress term iζ s,b in (30) represents the basal Coulomb friction (23).With the help of ( 35) and the fraction of the total normal stress for the solid constituent at the base, the Coulomb friction of the solid constituent reads Because the normal stresses of the solid constituent are assumed to be isotropic as in hydraulic theories (the simplification used in the present study), i.e.T ξ ξ s , and because f is introduced for the decomposition of the total normal stress, the depth-averaged solid stress term in (30) may be recast as Adapting the Navier-Stokes equations with the no-slip boundary condition at the basal surface, the fluid stresses in (30) can then be approximated by where with i ∈ {x,y,z} and m ∈ {ξ,η}.The dimensionless parameter N R = ρH √ gL(ν f µ f ) −1 is a dynamic scaling factor analogous to the Reynolds number for Newtonian fluids (cf.Iverson and Denlinger, 2001;Pudasaini et al., 2005b).
With the derived stress and the depth-averaged approximations, the equation system ( 29)-( 30) can be rearranged into the forms, for the mass and momentum conservation, where N b = hc − KA c is the normal pressure on the basal surface.As mentioned in Gray et al. (1999), the term KA c tanδ = O( α+β ) is assumed to contribute to the leading-and first-order balances, i.e. α + β ≤ 1.In case of 1 < α + β < 1 + γ or α + β ≥ 1 + γ , this term can be neglected and the momentum balances are accurate to order O( α+β ), or O( 1+γ ), respectively.Note that the momentum Eq. ( 41) is comprised of three equations, but with (32), only two of them are independent.We choose the first two (i ∈ {x,y}) as the representative momentum equations.Equations ( 40)-( 41) are the resultant model equations.These equations essentially describe the evolution of the mass conservation and the Cartesian components with respect to the topography-fitted mesh system.The first term on the right hand side (RHS) of ( 41) accounts for the entrainment-deposition processes and the second term (RHS) is the Coulomb friction of the solid constituent to the flow.Since the x-and y-axes lie on the horizontal plane, the gravity components are zero, g x = g y = 0, and so is the fifth term on the RHS of (41).Thus, it appears that there is no gravityacceleration in the formulation.This, however, is not true, because the gravity acceleration in this formulation is incorporated into the model equations via the reacting basal pressure N b , the third term on the RHS of (41).This is different than the models by Pudasaini and Hutter (2003) and Pudasaini et al. (2005b,c), in which the gravitational force is decomposed into three non-zero components along three mutually orthogonal coordinates that follow the channel.In these models, the basal Coulomb friction is also enhanced by the curvature and the twist of the channel.
The formulation of ( 40)-( 41) has several advantages over the existing models of the Coulomb-mixture type.First, due to the fact that the orientations of the x-and y-axes are fixed in space, the complexity for computing the evolution of the conservative variables is greatly reduced.Second, even with deposition/erosion, the coordinate surface ζ = 0 always coincides with the topographic surface, so that the mesh automatically fits the moving topography with the progression of time.

Erosion-deposition rate
The angle of repose is one of the most important features of granular media.With respect to this angle and the BCRE model proposed by Bouchaud et al. (1994), Tai and Kuo (2008) suggested that the erosion-deposition processes can be divided into three states for dry granular material relative to the neutral angle (the angle of repose of the material, θ n ), and a speed threshold (corresponding to a kinetic energy threshold).These three states are: immobile basal surface, deposition, and erosion, which take place at the base.In the present study, we consider a saturated mixture, in which case the interstitial fluid may affect the value of the neutral angle.Rondon et al. (2011) investigated granular collapses in a fluid and found there to be a roughly linear proportion between the initial solid volume fraction and the final deposition slope.In Coulomb-mixture theory, the solid normal stress at the base is assumed to be a fraction, 1 − f , of the total basal normal stress.Similar to the dry case, the local natural angle of the mixture is thus assumed to be linearly proportional to the fraction of the solid normal stress, i.e. θ nl = (1 − f )θ n .Therefore, the three states of the erosion-deposition mechanism correspond to In ( 42), θ indicates the local inclination angle, ||v || represents the tangential speed along the topography and the speed threshold v th is given by v th = α v (θ − θ nl ) 2 , where α v is an empirical constant.Utilizing the same concept as in Tai et al. (2012), the local inclination angle θ is defined to be the principal inclination angle determined by the ratio of the vertical velocity component, v z 0 , to the tangential speed over the topography, viz. and Following Tai and Kuo (2008) and Tai and Lin (2008), the deposition rate is proposed to be a function of the local thickness of the flowing layer and the difference between the inclination angle and the neutral angle, where h = h + α h √ h with an adjustment coefficient α h that manifests the erosion-deposition rate when the flow is very thin.The three states are implemented by the function e r , where H( the erosion-deposition rate and to distinguish the deposition from stagnancy, the function f reg v 0 , v th , e α is introduced in such a way that with a transition steepness factor e α .Based on (18), we relate , where the parameter α E is essentially the density ratio between the flowing layer and the stationary bed, if one takes into account the jump condition of the mass conservation at the interface.

Numerical investigation
The resultant model Eqs.( 40)-( 41) comprises a nonlinear system, for which a high resolution non-oscillatory central (NOC) scheme (see e.g.Jiang and Tadmor, 1998;Tai et al., 2001;Liu et al., 2007) is applied to solve the steep gradients of the variables, J b h, J b hv x and J b hv y .For the determination of the location of the moving mesh, a trapezoidal integration is applied.At time t k , let the position of the (i,j )−cell in the Cartesian coordinates be denoted by x k i,j = (x k i,j , y k i,j , z k i,j ) T , and let ) T be the cell velocity.The new position of the cell is then updated by where the time step t k = t k+1 − t k is determined by the Courant-Friedrichs-Lewy (CFL) condition.Note that there are alternative updating schemes for determining the coordinate position, and that some errors might be introduced during the simple updating process of (48).Nevertheless, the errors are small and acceptable (see e.g.Hui et al., 1999).For further details of the numerical implementation, the readers are referred to Tai and Lin (2008); Tai et al. (2012).
In the numerical investigation, the equations are solved in dimensionless form.We consider a finite mass of a saturated mixture sliding down an inclined flat chute (inclination angle 40 • onto a horizontal plane.The transition zone lies between ξ = 17.5 and 24.5 (dimensionless).The initial shape of the mass is a parabolic cap with a base radius of 2.4 (dimensionless).Its maximum height is 0.5.The center of the cap is located at (4.6, 0.0) on the ξ -η-surface, and the tangential components of the initial velocity are given by  et al., 2001;Liu et al., 2007) is applied to solve the steep gradients of the variables, J b h, J b hv x and J b hv y .For the determination of the location of the moving mesh, a trapezoidal integration is applied.At time t k , let the position of the (i,j)−cell in the Cartesian coordinates be denoted by x k i,j = (x k i,j , y k i,j , z k i,j ) T , and let ) T be the cell velocity.The new position of the cell is then updated by where the time step ∆t k = t k+1 − t k is determined by the Courant-Friedrichs-Lewy (CFL) condition.Note that there are alternative updating schemes for determining the coordinate position, and that some errors might be introduced during the simple updating process of (48).Nevertheless, the errors are small and acceptable (see e.g.Hui et al., 1999).For further details of the numerical implementation, the readers are referred to Tai and Lin (2008); Tai et al. (2011).
In the numerical investigation, the equations are solved in dimensionless form.We consider a finite mass of a saturated mixture sliding down an inclined flat chute (inclination angle 40 • onto a horizontal plane.The transition zone lies between ξ = 17.5 and 24.5 (dimensionless).The initial shape of the mass is a parabolic cap with a base radius of 2.4 (dimensionless).Its maximum height is 0.5.The center of the cap is located at (4.6, 0.0) on the ξ-η-surface, and the tangential  constituent).The basal friction coefficient of the solid constituent against the chute surface in ( 23) is taken to be µ = tanδ = tan23 • .Once the deposition takes place, the material accumulates on the basal surface and the friction angle is then changed to θ n .Parameter N R is set to be 3 × 10 5 as in Pudasaini et al. (2005c).In the computation, the values of the parameters of the erosion-deposition rate E G in (45) are identical to those in Tai and Lin (2008), which are α v = 1.0, α h = 0.1, α e = 2.0, e α = 20 and α E = 0.9.The initial mesh size is set to be ξ = η = 0.2; the aspect ratio = 1; the superbee slope limiter is used for cell-reconstruction of the physical quantities; and the CFL number is selected to be 0.4.The momentum correction factor is set to be 1.0 in the present simulation for ease of comparison.The readers are referred to Tai et al. (2012) for a more detailed discussion of the factor .In the result figures (Figs. 2 to 7), the flowing mass is represented by the sequential contours of flow thickness in which the levels of the contour lines are 0.001, 0.01 and from 0.05 to 0.49 at increments of 0.04.The thick red lines indicate the outlines of the deposition heaps and there is a transition zone lying between the vertical dash-dotted lines.
Figure 2 shows the thickness contour plots of the dry mass ( f = 0) moving down the chute without the entrainmentdeposition process.In the inclined section of the chute, the flowing mass accelerates in the down-slope ξ -direction and slightly extends in the transverse η-direction.Once the front s flows of the Coulomb-mixture type over temporally varying topography The local neutral angle of the mixture is determined by θ nl = (1 − Λ f )θ n with θ n = 34 • (the angle of repose of the solid constituent).The basal friction coefficient of the solid constituent against the chute surface in ( 23) is taken to be µ = tanδ = tan23 • .Once the deposition takes place, the material accumulates on the basal surface and the friction angle is then changed to θ n .Parameter N R is set to be 3 × 10 5 as in Pudasaini et al. (2005c).In the computation, the values of the parameters of the erosion-deposition rate E G in (45) are identical to those in Tai and Lin (2008), which are α v = 1.0, α h = 0.1, α e = 2.0, e α = 20 and α E = 0.9.The initial mesh size is set to be ∆ξ = ∆η = 0.2; the aspect ratio ǫ = 1; the superbee slope limiter is used for cell-reconstruction of the physical quantities; and the CFL number is selected to be 0.4.The momentum correction factor ̟ is set to be 1.0 in the present simulation for ease of comparison.The readers are referred to Tai et al. (2011) for a more detailed discussion of the factor ̟ .
In the result figures (Figs. 2 to 7), the flowing mass is represented by the sequential contours of flow thickness, in which the levels of the contour lines are 0.001, 0.01 and from 0.05 to 0.49 at increments of 0.04.The thick red lines indicate the outlines of the deposition heaps and there is a transition zone lying between the vertical dash-dotted lines.
Figure 2 shows the thickness contour plots of the dry mass (Λ f = 0) moving down the chute without the entrainment- in the computation; see also the comparison with experimental observations in Tai and Kuo (2008).
Figures 4 to 7 show the results computed with Λ f = 0.1 to Λ f = 0.4 in increments of 0.1.Apparently, there are no significant differences among the outlines of the flow body at t = 3.0 when they are flowing on the inclined section which the levels of the contour lines are 0.001, 0.01 and from 0.05 to 0.49 at increments of 0.04.The thick red lines indicate the outlines of the deposition heaps and there is a transition zone lying between the vertical dash-dotted lines.
Figure 2 shows the thickness contour plots of the dry mass (Λ f = 0) moving down the chute without the entrainmentdeposition process.In the inclined section of the chute, the flowing mass accelerates in the down-slope ξ-direction and slightly extends in the transverse η-direction.Once the front reaches the horizontal plane, the basal friction decelerates the mass and brings the front to rest.At t = 17 the entire flow body is "almost" at the state of rest.There is no clear criterion to determine the duration of the motion as argued in Tai and Kuo (2008).In contrast to Fig. 2, Fig. 3 depicts the results of the same condition but with the entrainmentdeposition procedure.Similarly, the flow front slows down on the horizontal section.The deposition heap begins to form at the front part slightly before t = 12, and develops backwards.The whole flow body comes to the state of rest at t = 16.526, when the total volume of the flowing layer is less than 10 −6 of the initial volume V int = 4.520.No significant difference is found between the shapes of the deposit of the sliding material body between Figs. 2 and 3.However, with the deposition procedure in the latter case, much more details of the flow are revealed, especially in relation to the formation of the deposit.These details illustrate the advantages to be gained by adopting the entrainment-deposition procedure in the computation; see also the comparison with experimental observations in Tai and Kuo (2008).Figures 4 to 7 show the results computed with Λ f = 0.1 to Λ f = 0.4 in increments of 0.1.Apparently, there are no significant differences among the outlines of the flow body at t = 3.0 when they are flowing on the inclined section (ξ ≤ 17.5).This may be due to the short distance from the initial position, so the spreading of the mixture is still in the beginning stages.However, as in Pudasaini et al. (2005c), the influence of the fluid constituent on the sliding mass becomes significantly distinguishable when the mass is deposited on the horizontal plane.For the values of λ f = 0.0 to Λ f = 0.4, the duration of sliding increases from 16.526 to 23.749 dimensionless time units.In addition, the maximum run-out distance increases and the deposition heap extends in the longitudinal direction as the value of Λ f increases.This is because a larger value of Λ f provides higher fluid-like mobility which overcomes the basal friction caused by the solid constituent moving against the basal surface.This matches the physical intuition.In all Λ f = 0 cases, the deposition heap begins to form around t = 12.However, it is worth noting that the deposition heap develops from the front for Λ f = 0.0 (see Fig. 3), whilst it is found to form initially at the rear flanks and develops forwards to the front for Λ f ≥ 0.2 (Figs. 5 to 7).
Our attention is further drawn to the finding that the fluid constituent leads to levee formation at the rear.This levee reaches the horizontal plane, the basal friction decelerates the mass and brings the front to rest.At t = 17 the entire flow body is " almost" at the state of rest.There is no clear criterion to determine the duration of the motion as argued in Tai and Kuo (2008).In contrast to Fig. 2, Fig. 3 depicts the results of the same condition but with the entrainmentdeposition procedure.Similarly, the flow front slows down on the horizontal section.The deposition heap begins to form at the front part slightly before t = 12, and develops backwards.The whole flow body comes to the state of rest at t = 16.526, when the total volume of the flowing layer is less than 10 −6 of the initial volume V int = 4.520.No significant difference is found between the shapes of the deposit of the sliding material body between Figs. 2 and 3.However, with the deposition procedure in the latter case, many more details of the flow are revealed, especially in relation to the formation of the deposit.These details illustrate the advantages to be gained by adopting the entrainment-deposition procedure in the computation; see also the comparison with experimental observations in Tai and Kuo (2008).
Figures 4 to 7 show the results computed with f = 0.1 to f = 0.4 in increments of 0.1.Apparently, there are no significant differences among the outlines of the flow body at t = 3.0 when they are flowing on the inclined section (ξ ≤ 17.5).This may be due to the short distance from the initial position, so the spreading of the mixture is still in the beginning stages.However, as in Pudasaini et al. (2005c)  formation cannot be obtained in the single-phase simulation (see Fig. 2 or 3).With the Coulomb-mixture theory, Pudasaini et al. (2005c) were the first to obtain levee formation in their numerical investigation.They called it the "reverse Barchan dune" type.In their theory, the deposition procedure was not taken into account.With the introduction of the deposition process, our simulation results agree with the results obtained in that former study.Furthermore, the deposition process enables the possibility of analyzing the geometric evolution of the deposition heap in greater detail.In our simulation with Λ f = 0.0, deposition initiates from the margins.When the deposition is initiated at the margins of the rear part, the central part moves at a higher speed in the down-slope direction, which postpones the deposition.The development of the levee deposition along the transverse ηdirection at ξ = 33.0 is illustrated in Fig. 8.At t = 15, deposition heaps begin to form at the side flanks.While the mass-flux decreases and passes through this cross-section, the levee deposition heaps bound the flowing mass, developing fast toward the center (ζ = 0).The thickness of the flowing layer then decreases rapidly depositing a growing heap.The mass flux vanishes shortly after t = 18.0.After this stage, the flowing mass freezes and the final deposition profile remains in a stagnant state.

Concluding remarks
In this study, a depth-integrated equation system with deforming coordinates for a Coulomb-mixture over general topography of small curvatures is introduced.The combination of the unified coordinate (UC) method and BW's arbitrary coordinate system enables the possibility that the coordinates can be set to move coincidentally with the temporally varying influence of the fluid constituent on the sliding mass becomes significantly distinguishable when the mass is deposited on the horizontal plane.For the values of λ f = 0.0 to f = 0.4, the duration of sliding increases from 16.526 to 23.749 dimensionless time units.In addition, the maximum run-out distance increases and the deposition heap extends in the longitudinal direction as the value of f increases.This is because a larger value of f provides higher fluid-like mobility which overcomes the basal friction caused by the solid constituent moving against the basal surface.This matches the physical intuition.In all f = 0 cases, the deposition heap begins to form around t = 12.However, it is worth noting that the deposition heap develops from the front for f = 0.0 (see Fig. 3), whilst it is found to form initially at the rear flanks and develops forwards to the front for f ≥ 0.2 (Figs. 5 to 7).
Our attention is further drawn to the finding that the fluid constituent leads to levee formation at the rear.This levee formation cannot be obtained in the single-phase simulation (see Figs. 2 or 3).With the Coulomb-mixture theory, Pudasaini et al. (2005c) were the first to obtain levee formation in their numerical investigation.They called it the " reverse Barchan dune" type.In their theory, the deposition procedure was not taken into account.With the introduction of the deposition process, our simulation results agree with the results obtained in that former study.Furthermore, the deposition process enables the possibility of analyzing   formation cannot be obtained in the single-phase simulation (see Fig. 2 or 3).With the Coulomb-mixture theory, Pudasaini et al. (2005c) were the first to obtain levee formation in their numerical investigation.They called it the "reverse Barchan dune" type.In their theory, the deposition procedure was not taken into account.With the introduction of the deposition process, our simulation results agree with the results obtained in that former study.Furthermore, the deposition process enables the possibility of analyzing the geometric evolution of the deposition heap in greater detail.In our simulation with Λ f = 0.0, deposition initiates from the margins.When the deposition is initiated at the margins of the geometric evolution of the deposition heap in greater detail.In our simulation with f = 0.0, deposition initiates from the margins.When the deposition is initiated at the margins of the rear part, the central part moves at a higher speed in the down-slope direction, which postpones the deposition.The development of the levee deposition along the transverse η-direction at ξ = 33.0 is illustrated in While the mass-flux decreases and passes through this crosssection, the levee deposition heaps bound the flowing mass, developing fast toward the center (ζ = 0).The thickness of the flowing layer then decreases rapidly depositing a growing heap.The mass flux vanishes shortly after t = 18.0.After this stage, the flowing mass freezes and the final deposition profile remains in a stagnant state.

Concluding remarks
In this study, a depth-integrated equation system with deforming coordinates for a Coulomb-mixture over general topography of small curvatures is introduced.The combination of the unified coordinate (UC) method and BW's arbitrary coordinate system enables the possibility that the coordinates can be set to move coincidentally with the temporally varying basal surface.This temporally varying, basal surface mimics the erosion-deposition processes.The evolution of the basal surface is postulated to be a function of the basal friction coefficient, sliding velocity, local thickness of the flowing layer, and a kinetic energy threshold.
The key features of the proposed model are illustrated by numerical simulations, in which a finite mass of a solid-fluid mixture slides down an inclined flat chute to be deposited on a horizontal plane.The levee formation, first obtained in Pudasaini et al. (2005c), is reproduced.By the introduction of the entrainment-deposition process, it is possible to determine the duration of the flow motion and investigate the evolution of the deposition heap in greater detail.The duration and the maximum run-out distance increase as the value of parameter f increases, although the deposition takes place around the same time point.Unlike with the dry granular flow, deposition of the sliding mixture initially forms at the rear flanks and then develops forwards to the flow front for f ≥ 0.2.In the transverse direction of the flow, the deposition is initiated from the side margins and then develops towards the flow center.It is not only the Coulomb friction which brings the moving mass to the state of rest, but the viscous force at the basal surface which also plays an important role.The effect of the fluid viscous drag becomes significant when the flow thickness tends towards zero.This is conjectured to be the cause of the levee deposition(formation).
The assumption of the small relative velocity between the solid and fluid constituents is nevertheless a strong restriction to many real debris flows.In addition, the fraction parameter f is assumed to be a constant throughout the flow body during the simulation, and the erosion-deposition rate is based on a simple intuitive approach.With these simplifications, the present work sheds light on the roles that the entrainment and deposition processes can play in this process.This topic is still one of the most poorly understood ones in debris-flow science (cf.Matthias and Hungr, 2005), calling for future study and improvement.
curvature and shallow flow depth assumptions, i.e. ζ ∂ x s = O( 1+α ), the Jacobian matrix and the tangential vectors of local natural basis reduce to

Fig. 2 .
Fig. 2.Contour plots of the flowing mass down the inclined flat plane onto the horizontal area, where the material is assumed to be dry (Λ f = 0), without adopting the deposition procedure.

Fig. 2 .
Fig. 2.Contour plots of the flowing mass down the inclined flat plane onto the horizontal area, where the material is assumed to be dry ( f = 0), without adopting the deposition procedure.

Fig. 3 .
Fig.3.Contour plots of the flowing mass and the evolution of the deposition heap, where Λ f = 0 and the deposition procedure is adopted in the computation.

Fig. 3 .Fig. 4 .
Fig. 3. Contour plots of the flowing mass and the evolution of the deposition heap ,where f = 0 and the deposition procedure is adopted in the computation.Y.C. Tai and C.Y. Kuo: Debris flows of the Coulomb-mixture type over temporally varying topography 9

Fig. 5 .
Fig.5.Contour plots of the flowing mass and the evolution of the deposition heap, where Λ f = 0.2 with the deposition procedure being adopted in the computation.

Fig. 4 .Fig. 4 .
Fig. 4.Contour plots of the flowing mass and the evolution of the deposition heap, where f = 0.1 with the deposition procedure being adopted in the computation.

Fig. 5 .
Fig.5.Contour plots of the flowing mass and the evolution of the deposition heap, where Λ f = 0.2 with the deposition procedure being adopted in the computation.

Fig. 5 .
Fig.5.Contour plots of the flowing mass and the evolution of the deposition heap, where f = 0.2 with the deposition procedure being adopted in the computation.

Fig. 6 .Fig. 7 .
Fig.6.Contour plots of the flowing mass and the evolution of the deposition heap, where Λ f = 0.3 with the deposition procedure being adopted in the computation.

Fig. 8 .
Fig. 8. Cross-sectional view of the flowing mass and evolution of the deposition heap along η at ξ = 33.0.The black dotted lines at ζ = 0 represent the basal surface of the chute, the red lines indicate the interface of the deforming deposition heap and the blue solid lines mark the free surface of the flowing layer.

Fig. 6 .
Fig.6.Contour plots of the flowing mass and the evolution of the deposition heap, where f = 0.3 with the deposition procedure being adopted in the computation.

Fig. 6 .Fig. 7 .
Fig.6.Contour plots of the flowing mass and the evolution of the deposition heap, where Λ f = 0.3 with the deposition procedure being adopted in the computation.

Fig. 7 .Fig. 6 .Fig. 8 .
Fig. 7. Contour plots of the flowing mass and the evolution of the deposition heap, where f = 0.4 with the deposition procedure being adopted in the computation.10 Y.C. Tai and C.Y. Kuo: Debris flows of the Coulomb-mixture type over temporally varying topography

Fig. 8 .
Fig. 8. Cross-sectional view of the flowing mass and evolution of the deposition heap along η at ξ = 33.0.The black dotted lines at ζ = 0 represent the basal surface of the chute, the red lines indicate the interface of the deforming deposition heap, and the blue solid lines mark the free surface of the flowing layer.

Figt
= 15, deposition heaps begin to form at the side flanks.