A new depth-averaged model for ﬂ ow-like landslides over complex terrains with curvatures and steep slopes

Flow-like landslides are one of the most catastrophic types of natural hazards due to their high velocity and long travel distance. They travel like ﬂ uid after initiation and mainly fall into the ‘ ﬂ ow ’ movement type in the updated Varnes classi ﬁ cation (Hungr et al., 2014). In recent years, depth-averaged models have been widely reported to predict the velocity and run-out distance of ﬂ ow-like landslides. However, most of the existing depth-averaged models present di ﬀ erent shortcomings for application to real-world simulations. This paper presents a novel depth-averaged model featured with a set of new governing equations derived in a mathematically rigorous way based on the shallow ﬂ ow assumption and Mohr-Coulomb rheology. Particularly, the new mathematical for- mulation takes into account the e ﬀ ects of vertical acceleration and curvature e ﬀ ects caused by complex terrain topographies. The model is derived on a global Cartesian coordinate system so that it is easy to apply in real- world applications. A Godunov-type ﬁ nite volume method is implemented to numerically solve these new governing equations, together with a novel scheme proposed to discretise the friction source terms. The hy- drostatic reconstruction approach is implemented and improved in the context of the new governing equations, providing well-balanced and non-negative numerical solutions for mass ﬂ ows over complex domain topo- graphies. The model is validated against several test cases, including a ﬁ eld-scale ﬂ ow-like landslide. Satisfactory results are obtained, demonstrating the model's improved simulation capability and potential for wider appli- cations.


Introduction
Flow-like landslides are one of the most catastrophic types of natural hazards due to their high velocity and long travel distance. They travel like fluid after initiation and mainly fall into the 'flow' movement type in the updated Varnes classification (Hungr et al., 2014). Numerical models have been widely used to predict the dynamics of flow-like landslides and hence quantify the run-out distance and flow velocity to facilitate risk assessment and management. Due to their much simplified formulation and less computational demand compared with the fully 3D models, depth-averaged models have been widely reported and successfully applied to simulate granular flows, including flow-like landslides. Savage and Hutter (1989) made the first attempt to develop a depth-averaged model for granular flows based on the Mohr-Coulomb internal rheology law and constant Coulomb bed friction. Their approach has since been adopted and extended by numerous researchers to develop granular flow models (e.g. Hungr, 1995;Iverson, 1997;Gray et al., 1999;Iverson and Denlinger, 2001;Denlinger and Iverson, 2001;Gray et al., 2003;Mcdougall and Hungr, 2004;Denlinger and Iverson, 2004;Pudasaini and Hutter, 2003;Pudasaini et al., 2005;Mangeney et al., 2007;Luca et al., 2009Luca et al., , 2012Gray and Edwards, 2014;Edwards and Gray, 2014;Iverson and George, 2014;George and Iverson, 2014).
For geophysical granular flows such as avalanches, landslides and debris flows, a challenging task is to simulate the real-world events taking into account the effect of complex 3D topographies. A major difference between the flow-like landslides and water flows, such as river flows or overland flood waves, is that flow-like landslides usually take place on steep slopes rather than nearly horizontal and flat ground surface. This poses a major challenge in developing depth-averaged models. As shown in Fig. 1, over a steep slope, the vertical acceleration (a v ) of a particle or mass element is non-zero. Consequently, the pressure distribution along the vertical direction can no longer be trivially calculated in the same way as the conventional shallow water equations defined on the Cartesian coordinates. In addition, there exists a centrifugal force (a c ) along the direction normal to the bed when it is curved. Incorporating the vertical acceleration and centrifugal force into the depth-averaged models is essential for accurate simulation of granular flows over complex terrains. This has been a challenge since the very first depth-averaged model was developed.
When deriving a depth-averaged model, the flow direction is normally assumed to be parallel to the domain surface, so that the normal pressure term is trivial to calculate. As a result, surface-fitted curvilinear coordinate systems have been widely adopted for developing depth-averaged granular flow models. The first of such attempts was made by Gray et al. (1999) through the introduction of an orthogonal curvilinear coordinate system to develop their model; but their model assumes that predominant topographic variation only occurs in one direction and so faces a major difficulty in adapting to real-world terrains with more complex topographic features. Later on, Pudasaini and Hutter (2003) introduced a model on a non-orthogonal curvilinear coordinate system for avalanches in arbitrary curved and twisted channels. Such a non-orthogonal curvilinear coordinate system has also been adopted in other granular flow models (e.g. Pudasaini et al., 2005Pudasaini et al., , 2007 and achieved certain level of success. A curvilinear coordinate system relies on the thalweg along the bed to define its main axis. For a real-world complex topography, however, it is usually difficult to define the downslope direction or the thalweg due to large variations of the topography in different directions. Furthermore, to facilitate the simulation of geophysical flows (e.g. landslides, debris flows and avalanches) in the real world, Digital Elevation Models (DEMs) are commonly used to describe the terrain topography, which is commonly defined on a Cartesian coordinate system. Transformation of topographic data from the Cartesian coordinate system to the curvilinear coordinate system must be applied, inevitably increasing computational effort and leading to a loss of accuracy, particularly in the cases where the topographies are featured with abrupt changes. To avoid this, Bouchut and Westdickenberg (2004) introduced a shallow water flow model on an arbitrary coordinate system for simulations over topographies with small curvatures. Their overall governing equations were derived on a fixed Cartesian coordinate system while the variables were defined on a local reference coordinate system aligning with the local topography, i.e. the flow depth is normal to and the velocities are parallel to the bed. This model was later extended and applied to granular flows by Mangeney et al. (2007). In their model, the terrain topographies may be directly described by a DEM, but coordinate transformation is still needed to provide the initial depth along the vertical axis. For real-world applications with complex topographies, defining flow depth normal to bed can be inconvenient because performing such coordinate transformation is not only time-consuming but also sometimes difficult, if not impossible, especially when discontinuous topographies arise.
In practice, it is desired to develop a model based on a global Cartesian coordinate system in which the vertical axis is aligned with the gravity direction so that DEMs can be directly used to support model setup without the need of coordinate transformation. But a global Cartesian coordinate system based model may also have its limitations, as the determination of pressure/stress terms becomes a much more difficult task within such a configuration. Granular flows commonly happen on steep inclined slopes. The flow acceleration along the vertical direction has a magnitude comparable to the gravity and is not negligible (Iverson, 2014). The vertical normal pressure thus becomes more difficult to calculate. The centrifugal force caused by bed surface curvature is also no longer trivial to quantify because the velocity variables are not defined parallel to the terrain surface.
In order to utilise a global Cartesian coordinate system and meanwhile maintain solution accuracy, Denlinger and Iverson (2004) subtracted the vertical acceleration from the gravity acceleration to account for the non-hydrostatic pressure effect; their model provided better results than a hydrostatic model when applying to a granular dam break test. More recently, Castro-Orgaz et al. (2014) suggested that the Boussinesq-type models which retain the non-hydrostatic pressure to a certain level and have been successfully applied in modelling shallow water waves may be also implemented for simulating gravitydriven granular flows. They also pointed out that Denlinger and Iverson's model can actually be categorized as a Boussinesq-type model for granular flow. The model by Castro-Orgaz et al. (2014) has recently been further simplified by Yuan et al. (2017) and implemented in an existing code. These Boussinesq-type non-hydrostatic models generally predict better results than their hydrostatic counterparts. However, more sophisticated numerical schemes must be used to solve the Boussinesq-type governing equations due to the presence of additional higher-order derivative terms. Therefore, a Boussinesq-type model is computationally much more demanding and usually less stable than a model solving the shallow water equations or similar depth-averaged formulations.
Other simpler global Cartesian coordinate system based models have also been reported (e.g. Juez et al., 2013;Hergarten and Robl, 2015). These models simply modified the original shallow water type equations by including a projection factor to the pressure and source terms, determined by the consideration of bed or surface topography gradients according to heuristic geometric arguments. These models produce very similar results to those models based on local curvilinear coordinate systems, which has been confirmed by the authors' previous study (Xia et al., 2015). Compared with the Bousinessq-type models, the numerical implementation of these shallow water type models can be much easier to achieve and many well-documented numerical schemes developed for shallow flow hydrodynamics can be directly used. However, although these models have proven to be successful for certain applications, they have not been fully justified in a mathematically rigorous way and all of them do not consider the effect of the centrifugal force induced by bed curvatures which may become significant for applications involving complex topographies.
Numerous robust numerical schemes have been reported in the literature for solving the shallow water equations in the context of hydrodynamic simulation. In the last two decades, particular attention has been paid to develop shock-capturing numerical schemes to support accurate and stable simulation of shallow flow hydrodynamics over dry terrains with complex topographies (e.g. Gray et al., 2003;Audusse et al., 2004;Liang and Marche, 2009;Hou et al., 2014). Such numerical schemes are generally required to maintain the C-property (i.e. preserving the lake at rest solution at the discrete level) and include proper numerical techniques to handle wet/dry interface and discretise the friction source terms. However, some of these issues (e.g. C-property and friction term discretisation) have not been thoroughly considered and resolved in the context of flow-like landslide modelling, which calls for more research efforts (e.g. Mangeney et al., 2007;Juez et al., 2013;Zhai et al., 2015).
In order to correctly take into account the effects of large slope gradients and curvatures, but meanwhile allow the users to directly take advantages of DEM data, this work presents a new depth-averaged model based on a global Cartesian coordinate system with the following highlights: 1. New depth-averaged equations are derived through depth-integration and asymptotic analysis, taking into account the effects of vertical acceleration and centrifugal force; the resulting equations are hyperbolic and rotationally invariant, and mathematically preserve the lake at rest solution. 2. The hydrostatic reconstruction method is improved and implemented to properly handle wetting and drying and maintain the C-property in the context of a second-order Godunov-type finite volume scheme. 3. A splitting method is proposed to discretise the friction source terms for accurate and stable simulations.
The rest of the paper is organised as follows: Section 2 presents the depth-averaged governing equations for which the detailed derivation is provided in Appendix A; Section 3 introduces the proposed numerical scheme; the model is then validated by carefully selected test cases in Section 4, followed by brief conclusions in Section 5. A glossary of notations is provided at the end of the manuscript.

Governing equations
This section introduces the new depth-averaged governing equations derived on a global Cartesian coordinate system and presents the relevant mathematical properties.

The depth-averaged governing equations
The depth-averaged equations can be derived from the three-dimensional governing equations by assuming the Mohr-Coulomb rheology. The detailed derivation is presented in Appendix A, and the final depth-averaged equations in a matrix form are written as where the vector terms are given by (2) where (x, y) define the two-dimensional Cartesian coordinates, g is the acceleration due to gravity, h is the depth, b is bed elevation, and u and v are the x − and y −direction depth-averaged velocities. The above equations appear to be similar to the shallow water equations, which may be beneficial in terms of directly adopting many existing numerical methods originally developed for the shallow water equations. Compared with the conventional shallow water equations, in addition to the friction terms there are three major differences. Firstly, the gravity terms have an additional factor of 1/ϕ 2 that reduces the gravity effect. This factor is only related to the bed topography and is independent of either the coordinate system or the velocity direction. It is an essential condition to ensure rotational invariance of the above depth-averaged equations. The inclusion of this factor is theoretically important for the governing equations to properly describe the effects of complex topography in a Cartesian coordinate system. As a consequence, a flow modelled by the current equations may move slower than that predicted by the conventional shallow water equations. Secondly, the term v T Hv is included to account for the effect of centrifugal force. The effect of centrifugal force increases the normal pressure and hence the friction force, so that the movement modelled by the new equations may become slower. But the centrifugal force can also lead to faster flow movement than models without considering curvature in certain situations. Considering flow moving from an inclined slope into a horizontal plane, it is apparent that the velocity is always aligned to the slope direction. As the flow moves onto the horizontal plane, the horizontal velocity component is larger than it is on the inclined slope because the velocity has become fully aligned with the x-direction. However, if the centrifugal force is not considered, the vertical velocity component will disappear when it moves into the horizontal plane, leading to reduced velocity. Explicitly including the centrifugal force will retain the full magnitude of the velocity because it is the centrifugal forces that change the direction of the flow. Thirdly, the second terms  .40), are included to mathematically preserve the lake at rest condition as will be proved in the next subsection. Numerical experiments will later demonstrate that neglecting these terms will lead to inaccurate results when reproducing experimental granular flows.
Relevant formulations have been reported in the literature (e.g. Juez et al., 2013;Hergarten and Robl, 2015), but they do not include the new centrifugal force term and the extra terms in S b . Compared with the more complicated Boussinesq-like models (e.g. Iverson, 2004 andCastro-Orgaz et al., 2014), the models solving the current depth-averaged governing equations will be computationally much less demanding due to the use of a much simplified formulation and subsequently the use of less sophisticated numerical schemes.
Although they are derived based on the Mohr-Coulomb rheology, the new governing equations can easily incorporate with friction laws of varying coefficients, such as the velocity dependent friction law proposed by Pouliquen and Forterre (2002). In principle, they can be also extended to include the more complex two-phase rheologies.

Properties of the new depth-averaged equations
In this section, the new depth-averaged governing equations will be proved to possess three important properties, i.e. mathematical preservation of the lake at rest solution, hyperbolicity and rotational invariance, which are essential to guarantee accurate numerical solutions in the context of implementing a Godunov-type numerical scheme.

Mathematical preservation of the lake at rest solution
The lake at rest conditions may be defined as 0 , a n d 0 .
A model satisfying the above conditions is often referred to as wellbalanced (Greenberg and Leroux, 1996), or preserving the C-property (Bermúdez et al., 1998). Preservation of the lake at rest solution is a necessary condition to ensure numerically stable and physically correct simulation results. It is the simplest yet an important case of any steady stationary configurations. It primarily concerns the balance between the depth slope and the bed slope for shallow water flow simulations. But it is also relevant to landslide modelling, for which the friction slope is of more interest because the friction slope may also be generalised as a bed slope.
To prove that the model preserves the C-property, we may substitute Eq.(6) into Eqs.(1)-(3), where the continuity equation is automatically satisfied. Herein we focus on the x-direction momentum equation and the y-direction equation can be proved in a similar way. After eliminating all of the terms containing zero velocities, the x-direction momentum equation reduces to After applying the chain rule and performing simple manipulations, the equation can be further simplified to become Since we have h + b = const, the above equation is apparently satisfied, confirming that the new governing equations mathematically preserve the lake at rest solution.

Hyperbolicity
The Jacobian matrix corresponding to the flux terms is given by where the flux vector is f(q) = f(q)n x + g(q)n y , and with n x and n y being the two Cartesian components of the unit vector, and = c gh ϕ is the wave speed. The eigenvalues of J are x y x y x y 1 2 3 (11) Obviously, these eigenvalues are distinct real values when h≠0. Therefore the corresponding homogeneous equations of Eqs.
(1)-(3) are strictly hyperbolic when h≠0. This implies that the new depth-averaged mass flow governing equations can be numerically solved using a range of Godunov-type numerical schemes that have been widely developed for the shallow water equations.

Rotational invariance
Rotationally invariant equations are independent of the choice of coordinate directions, therefore facilitating the implementation of a more robust numerical scheme and ensuring more reliable simulation results. Adopting the following rotational matrix in which θ is an arbitrary angle. It is straightforward to verify that the fluxes of the depth-averaged equations satisfy and therefore rotationally invariant. Similarly, the source terms can be also proved to be rotationally invariant and so the overall formulation strictly satisfies the rotational invariance requirement.

Numerical method
In this work, the hyperbolic conservation laws formed by the newly derived depth-averaged governing equations are solved using a Godunov-type finite volume scheme to allow automatic shock-capturing numerical solutions. In order to effectively maintain the Cproperty and handle wetting and drying, particular attention is paid to implement a hydrostatic reconstruction approach (Audusse et al., 2004) with effective improvements to discretise the source terms in a divergent form. A fractional splitting method is used to evaluate the friction source terms to reinforce physically-sound numerical solutions.

Finite volume discretisation
A finite volume numerical scheme solves the integrated form of the governing equations and the resulting semi-discretised equation for an arbitrary cell 'i' is given by where 'k' is the index of the cell edges (N = 4 for the Cartesian grids adopted in this work) of cell 'i', l k is the length of cell edge 'k', Ω i is the cell area, and F k contains the interface fluxes at cell edge 'k'.
Herein, the overall solution procedure adopts a fractional splitting method to update the flow variables to a new time level, i.e. + q i n 1 , where n denotes the time level. In the first step, the flow variables are updated against only the interface fluxes and slope source terms using a second-order Runge-Kutta time marching scheme where the Runge-Kutta coefficient is defined as In the second step, the flow variables are fully updated by also taking into account the friction source terms using the following timemarching formula Details of the calculation for each term will be given in the following sub-sections.

Interface flux calculation
The interface fluxes in Eq. (14) are calculated by solving the Riemann problems defined locally at each cell interface, i.e.
where q L and q R are the Riemann states, simply the face values of the flow variables at either side of the interface obtained through projection onto the local coordinate direction defined at the cell interface in which n and n ⊥ are the basis vectors (see Fig. 2), with n being (1, 0), (0, 1), ( −1, 0), (0, −1) and n ⊥ being (0, 1), ( −1, 0), (0, −1), (0, 1) respectively for the east, north, west and south cell interfaces. For a second-order numerical scheme, the linear reconstruction as originally proposed by van Leer (1979) is used to estimate the face values of the flow variables from their cell-centred values: X. Xia, Q. Liang Engineering Geology 234 (2018) 174-191 where f denotes a flow variable, 'j' represents a neighbour of cell 'i', the subscripts 'L' and 'R' represent 'left' and 'right', ∇f is the unlimited gradient of f, which is calculated by an upwind difference method, r is the distance vector from cell centre to face centre, Ψ(r) is a limiter function to ensure monotonicity during the interpolation to suppress spurious oscillations in the numerical solutions. In this work, the simple minmod slope limiter is adopted, which is defined on a rectangular Cartesian grid for cell 'i' as in which the subscripts + and − denote the downstream and upstream cells, respectively. With the face values, the final left and right Riemann states are obtained by implementing the hydrostatic reconstruction method as proposed by Audusse et al. (2004) to avoid negative flow depth. The first step is to define a single bed elevation at the cell interface under consideration: where b L and b R are calculated by finding the difference between the linearly reconstructed face values of flow surface elevation (i.e. s = h + b) and depth. Then the left and right Riemann states of the depth are redefined as which are then used to reconstruct other Riemann states The Riemann states are then used to define the local Riemann problems, which are numerically solved to evaluate the interface fluxes.
Due to the similarity between the new mass flow governing Eqs. (1)-(3) and the classic shallow water equations, a standard approach developed for the shallow water equations can be used after certain modifications. Assuming that 1/ϕ 2 is constant across the cell interface (within the area enclosed by dashed line in Fig. 2), the new depthaveraged equations are locally reduced to the shallow water equations after defining a modified gravity acceleration; the modified gravity acceleration across the cell interface under consideration may be simply determined by finding the average between the two neighbouring cells: where the gradient terms of b can be calculated using the slope gradients constrained by a slope limiter to avoid excessively large gradients where the topography is discontinuous.
Therefore solving Eqs.(1)-(3) essentially becomes seeking numerical solution to the shallow water equations with a modified but still constant gravity acceleration. A Godunov-type finite volume scheme established for the shallow water equations is directly applied here. Due to its ease to incorporate wetting and drying treatment, an HLLC approximate Riemann solver (Toro, 2001) is implemented to evaluate interface fluxes, taking the x-direction flux f as an example, in which f L = f(q L ) and f R = f(q R ) are calculated from the left and right Riemann states, S L , S R and S M are the characteristic wave speeds, and f *L and f *R are the fluxes in the left and right middle regions of the HLLC solution structure, and calculated as follows with the HLL fluxes f * provided by the following formula The formulae for the left and right characteristic wave speeds S L and S R are The middle characteristic wave speed S M is calculated as Finally, the flux F k on a cell face 'k' can be obtained by projecting f back to the global coordinates, it is given as where f 1 , f 2 and f 3 are the three components of f.

Discretisation of slope source terms and C-property
Apart from the flux calculation, the slope source terms in S b must also be discretised properly to maintain the C-property. The formulation as proposed in the original hydrostatic reconstruction method by Audusse et al. (2004) cannot be directly implemented in the new mass flow simulation framework due to the difference in governing equations. A modified scheme must be used.
To discretise the slope source terms, we first integrate S b over the entire domain Ω of an arbitrary cell, i.e. X. Xia, Q. Liang Engineering Geology 234 (2018) In a spatially second-order scheme, the gradient term ∇s is constant within a cell, leading to the following relationship: where b M is defined at the centre of the domain, essentially the cell centred value. After applying the Green-Gauss theorem, the discretised equation for the slope source term vector S b can be obtained as where g k , calculated by Eq. (25), is the modified gravity acceleration at cell edge 'k'; = g g ϕ / i i 2 is defined at the cell centre; s k is the interpolated surface elevation at cell edge 'k'; h Lk is the hydrostatic reconstructed flow depth at cell edge 'k' obtained according to Eq.(23); a i is calculated according to Eq.(4) with the second-order derivatives approximated using central differences. Unlike the original hydrostatic reconstruction method, the above slope source term discretisation scheme does not involve any additional cell centred source terms to maintain the consistency of the overall numerical scheme, removing extra effort in numerical implementation.
The C-property of the overall numerical scheme can be proved as follows. For the lake at rest condition: s = const and v = 0, the centrifugal force term v T Hv vanishes, thus a i = g i , and s k is equal to s i at all cell edges; hence 1 2 2 will also disappear. As a result, the discretised S bi becomes It is trivial to verify that Eq.(40) exactly balances the summation of flux at all cell edges. Therefore the C-property is maintained.

Discretisation of friction source terms
If the soil is static and the forces induced by the pressure gradient and slope gradient are smaller than the friction force, a static resistant force smaller than the friction will exist to balance exactly the other forces to maintain the static soil condition. Friction forces alone can only stop the flow but can never reverse the flow. In order to represent this physical reality in the numerical model, the friction source term must be restricted as follows: where q = (0,hu,hv) T contains the flow variables updated using the second-order Runge-Kutta scheme as introduced previously; S* f is calculated from Eq.(3) with the updated water depth h and velocity vector v at the cell centres.
Due to the use of the above restriction measures, the current friction discretization scheme will never revert the flow. It is also straightforward to verify that the current scheme maintains the static soil condition. If friction is larger than other forces, the velocity will be slowed down to nil in the discretised friction source term formulation in Eq. (18). Therefore, the physical nature of the friction forces is correctly interpreted in the numerical scheme.

Results and discussion
In this section, the new mass flow model as presented in the previous sections is validated against several test cases, including an idealised uniform flow on an inclined frictional slope, three experimental granular flow tests and a real-world flow-like landslide event.

A uniform but unsteady flow on an inclined frictional slope
A uniform but unsteady flow on an inclined frictional slope is considered here to test the model's capability in capturing the effect of a large slope. Being a simple test, the analytical solution may be derived. This test has been previously considered by other researchers (e.g. Hergarten and Robl, 2015). Herein we only provide an overview with essential details. The flow occurs on a slope with an angle of θ and friction coefficient of μ. The flow depth is h and the depth-averaged velocity ũ is parallel to the slope surface. The flow acceleration can be then obtained from the Newton's second law   However, if the conventional shallow water equations without the modified gravity acceleration are used, the flow acceleration will be predicted as where a factor of cos 2 θ has been missed out compared with the correct solution. For real-world landslides, the slope can be as steep as 30°-40°, and the conventional shallow water equations will therefore overestimate the acceleration by up to 30%, leading to unacceptable errors in the velocity and run-out distance predictions.

USGS granular flow experiment
An experimental granular flow test reported by Denlinger and Iverson (2001) is used herein to further validate the current depthaveraged landslide model. The experiment was undertaken in USGS. In Fig. 6. Granular flow over a smooth two-dimensional bump: simulated and measured surface profiles at different moments.
X. Xia, Q. Liang Engineering Geology 234 (2018) 174-191 the experiment, 290 cm 3 of dry quartz was placed behind a gate on an inclined flume. The flume consists of three parts: an inclined slope of 31.4°, a curved section with a 10 cm radius and a horizontal run-out surface at the end. The width of the flume is 20 cm. A sketch of the experiment setup is shown in Fig. 3. The friction angle between the grains and the bed is 29°as suggested by Denlinger and Iverson (2001). During the experiment, the gate was suddenly opened to release the grains, which flowed rapidly downhill through the flume and finally deposited at the horizontal surface. The simulated flow profile at different output times is shown in Fig. 4. After the gate is opened, the mass front rapidly moves down the slope. When the mass front reaches the horizontal plane, it begins to deposit to form a pile, which is consistent with the observation during the experiment. The measured and predicted flow extents, defined by the front and rear positions of flow at different times, are compared in Fig. 5. Satisfactory agreement has been achieved, demonstrating the capability of the current model in reproducing this laboratory test. In order to investigate the influence of vertical acceleration and curvature effects, we respectively turn off the gravity correction factor, i.e. set 1/ ϕ 2 = 1, and omit centrifugal forces, and plotted the corresponding front and rear positions also in Fig. 5. Clearly, the vertical acceleration and curvature have significant influences on the results. Negligence of vertical acceleration leads to faster mass movement and the deposition of flow front at a further position. On the contrary, negligence of curvature effect causes slower mass motion and nearer deposition position of the flow material because the velocity component normal to the bed surface disappears when the direction of the bed slope changes. Viroulet et al. (2017) conducted an experiment to study the dynamics of dense granular flows over a two-dimensional symmetric bump. In their experiment, grains were released from the upstream end of an inclined 5 cm wide chute with a curved bump. They observed that even a small amount of grains initially deposited at the upstream side of the bump could slow down the flow and gradually form a steady shock. They also demonstrated that a terrain fitted avalanche theory is able to predict accurately the location of the steady shock with different slope angles; however a standard avalanche theory on Cartesian coordinates is not able to make correct predictions, suggesting that the local slope and curvature of the bump play an important role in the granular flow dynamics. Therefore, this experiment is chosen to further test confirm the capability the new depth-averaged model in correctly representing the effects of bed curvature.

Experimental granular flow over a smooth two-dimensional bump
In this simulation, the friction between the grains and side walls must also be considered. The formula suggested in Viroulet et al. (2017) is adopted and modified to reflect the different definition of depth where μ b and μ W are the basal and side wall friction coefficients respectively, and W is the width of the chute. We use μ b = tan23°, μ W = tan7.5°and W = 5 cm, identical to those used in Viroulet et al. (2017). As shown in Fig. 6, an erosional shock wave forms after the granular flow reaches the initial deposit at about t = 0.4 s, which then gradually moves upstream and develops into a steady shock. The current numerical predictions are in good agreement with the experiment observation and comparable with the results produced by a terrain fitted depth-averaged model (see Viroulet et al., 2017 for details). The rear section of the simulated flow depth profile is sharp, as opposed to the smooth transitions observed in the experiment. This is because the shock waves are idealised as discontinuous flow depth and velocities without explicitly considering any diffusing mechanism such as viscosity ( Fig. 7 also presents the simulation results with these terms being omitted. The simulations clearly fail to predict the locations of the shock wave, confirming that these new terms are indeed essential for accurate predictions. Successful reproduction of this experiment indicates that the effects of large slope steepness and curvature have been correctly represented in the current depth-averaged model.

Experiment granular flow with an obstacle
The granular flow experiment reported by Gray et al. (2003) is further simulated using the new model. The topography on which the grains travel is featured with a pyramid shaped obstacle on an inclined flume as shown in Fig. 8. At the beginning of the experiment, grains were released from upstream with a constant velocity and depth for 10 s until the controlling gate was shut off.
During the simulation, the friction angle between the flume and the granule is chosen as 32°, the same as the one used in Gray et al. (2003). The simulated flow patterns are shown in Fig. 9 for different output times, compared with the snapshots taken at the experiment. As the grains travel downstream, the flow front touches and is then blocked by the obstacle. The Coulomb basal friction balances the internal pressure  X. Xia, Q. Liang Engineering Geology 234 (2018) 174-191 gradients caused by the obstacle to create a stationary zone. At the meantime, a detached shock wave is formed upstream and separates the flow into two streams because the inflow velocity is greater than the wave speed gh in this case. Such a phenomenon is physically similar to a 'hydraulic jump' in a transcritical water flow. The simulated results are observed to be consistent with the experimental records. Particularly, the location of the shock wave front and the extent of the stationary zone both match well with those observed during the experiment, demonstrating that the new model is able to capture complex flow characteristics, such as shock waves formed around complex topographies.

2015 Heifangtai landslide, China
Heifangtai is a loess terrace located in Yongjing County, 42 km to the west of Lanzhou city, Gansu province, China. Various reasons including fluvial erosion at the base of slopes and rainfall infiltration may cause slope instabilities. According to previous field studies ( Xu et al., 2012, Fig. 9. Experimental granular flow with an obstacle: simulated (left column) and measured (right column; from Gray et al., 2003) flow patterns at different output times. To better illustrate the boundary of flow, depth smaller than 2 mm is plotted as void.
X. Xia, Q. Liang Engineering Geology 234 (2018) 174-191 2014), this region is prone to landslides particularly as a result of intensive agricultural irrigation and the subsequent rise of ground water level. The loess terrace can be as high as 100 m with steep slopes up to 30°, creating a large gravitational potential that may drive a landslide to travel rapidly for a long distance. Peng et al. (2017) conducted a comprehensive field survey of landslides in this region. Among all those landslides in this region, the one happened near to Jiaojia village in 2015, which is indexed as JJ4 in Peng et al. (2017), is particularly interesting because of its rich morphological features of both bed and deposit. As shown in the pre-event satellite image (Fig. 10), the topography on which the landslide mass travels is complex and featured with several trenches (TCs in Fig. 10) caused by previous slope failures and soil erosion and three stages (STs in Fig. 10) formed by construction works. The landslide traveled a long distance, reaching the bottom of the third stage (ST3), and the lateral extent are well defined by the trenches. Peng et al. (2017) carried out a quantitative survey of the deposit depth (see Fig. 7 of their paper for details), based on which we have sketched the areas where the deposit depth is significantly thicker than elsewhere on the post-event satellite image (Fig. 11). The sketches show that the deposits mainly concentrate at the bottom of the landslide scarp and along the trenches. This landslide's long travel distance (∼ 400 m) and small deposit depth (between 1 m and 10 m) indicate that the depth-averaged models including the one presented in this paper are applicable. Therefore, this landslide event is chosen to test our model's simulation capability, especially the numerical stability in the presence of highly irregular topographic features. For this simulation, we have chosen an empirical velocity-dependent friction law proposed by Pouliquen and Forterre (2002) and summarised in Gray and Edwards (2014) to match both the runout extent and the deposit morphology. The friction coefficient is calculated as where μ 1 and μ 2 are the static and dynamic friction coefficients re- is the Froude number, β is a dimensionless coefficient and ℒ is a critical length comparable to the particle diameter. The parameters being used in this simulation are summarised in Table 1. Simulation results at t = 15 s, 30 s, 45 s and 60 s are plotted in Fig. 12 to present the dynamic processes of the landslide. At the    Xia, Q. Liang Engineering Geology 234 (2018) 174-191 beginning, the soil moves down rapidly through the channels; after 30 s, the landslide reaches its maximum runout distance; then it begins to deposit and completely settles down after 60 s. The simulation has well captured both the extent and morphology of the final deposit. Similar to the satellite image and post-event survey, the deposit predicted by the simulation reaches the bottom of ST3, is laterally confined by the trenches, and mainly concentrates at the bottom of the landslide scarp and along the trenches. From the post-event satellite image (Fig. 11), the soil at the right hand side of the flow direction seemed to choose its path to travel through the low-lying area in-between TC3 and TC4 and finally reached the bottom of ST3 to form a tongue shape deposit after meeting with the flow at the left hand side; a fraction of the volume is deposited on its original course without enough bulk momentum to reach the bottom of ST1. The simulation reasonably reproduces the trend of the flow travelling leftwards. However, the amount of soil being diverted is not large enough, so the soil at the right hand side still reaches the bottom of ST1 and the tongue shape of the deposit at the bottom of ST3 has not been well predicted. A possible explanation is that, in reality, the flow may have entrained and widened TC4 and subsequently enabled more soil to travel through TC4. The effect of bed entrainment is currently not taken into account in the model. Generally, this test confirms the simulation capability of the present depth-average landslide model and demonstrates its applicability for simulating flow-like landslides over terrains with complex topographies.

Conclusions
This paper presents a new depth-averaged model for simulating flow-like landslides. Based on the shallow flow assumption and the Mohr-Coulomb rheology, new depth-averaged governing equations are derived in a mathematically rigorous way through asymptotic analysis in a global Cartesian coordinate system. A correcting coefficient taking into account the vertical acceleration due to large slope gradient, the centrifugal acceleration due to bed curvatures, and two additional terms to preserve the lake at rest solution are naturally included in the formulation without any ad-hoc assumptions. The final governing equations are rotationally invariant, hyperbolic and mathematically well-balanced to preserve the lake at rest solution. A second-order Godunov-type finite volume numerical scheme is then implemented to solve the governing equations. To ensure preservation of the C-property and non-negative flow depth during a simulation, the hydrostatic reconstruction approach is implemented with necessary modifications in the context of the new governing equations. A simple fractional scheme with a novel discretisation method for the friction terms is proposed to properly simulate the static resistance and stop condition. The new model is therefore physically sound and numerically robust to support landslide simulations over complex real-world terrains.
The new model has been validated against several test cases and is able to produce the exact solution to the uniform flow on a frictional slope and satisfactory results for three laboratory-scale test cases, including the more complex cases of granular flow interacting with a curved bump or a pyramid-shape obstacle. Finally, the model is applied to simulate a field-scale landslide with complex topographic features, successfully capturing the major characteristics of the landslide dynamics and predicting both final deposit extent and morphology informed by post-event survey. This confirms the simulation capability of the model and its potential for a wider range of applications.

Glossary of notations
a v , a c vertical and centrifugal accelerations x, y and z three dimensions in space t time q conservative variable of the depth-averaged equations f(q) x-direction flux g(q) y-direction flux S b source terms related to topography S f frictional source terms u and v depth-averaged velocities (three-dimensional velocities in Appendix A) in x and y dimensions respectively w velocity along the z direction u v , and w depth-averaged velocities in x, y and z dimensions respectively h flow depth b bed elevation s flow surface elevation g gravitational acceleration ϕ denominator in the unit vector normal to bed topography a modified gravitational acceleration v the vector of depth-averaged velocities H Hessian matrix U the vector of three-dimensional velocities X. Xia, Q. Liang Engineering Geology 234 (2018) 174-191 ∇ Gradient operator T Cauchy stress tensor g the vector of the gravity field p isentropic pressure τ deviatoric stress tensor I identity tensor L characteristic length H characteristic deptĥ b non-dimensional bed elevation s non-dimensional flow surface elevation h non-dimensional flow deptĥ u ,v and  w non-dimensional velocities ρ density τ xx , τ yy , τ zz , τ xy , τ yz , τ xz components of the deviatoric stress tensor τ xx ,τ yy ,τ zz ,τ xy ,τ yz ,τ xz components of the non-dimensional deviatoric stress tensor μ friction coefficient ε ration between characteristic depth H and characteristic length basal deviatoric stress tensor K lateral stress coefficient F(q) flux vector n x , n y Cartesian components of the unit vector J Jacobian matrix λ 1 , λ 2 , λ 3 eigenvalues of the Jacobian matrix R rotational matrix Ω cell area l length of cell edge K Runge-Kutta coefficients Δt time step length n, ⊥ n unit vectors normal and perpendicular to cell interfaces f arbitrary flow variables S L , S R , S M left, right and middle characteristic wave speeds ũ depth-averaged velocity parallel to the slope surface W width of the chute μ b , μ W friction coefficients for the bed and wall μ 1 , μ 2 Static and dynamic friction coefficients in the velocity-dependent friction law β dimensionless number in the velocity-dependent friction law F r Froude number ℒ critical length in the velocity-dependent friction law

Acknowledgment
This work is partly funded by the NERC SINATRA and TENDERLY projects (grant no. NE/K008781/1). The authors thank Professor Fuchu Dai for his help and advice on the 2015 Heifangtai landslide. Thanks also go to Dr Yanbo Cao and Dr Haisong Liu for their helpful discussion. Three anonymous reviewers have given constructive comments, which have helped substantially improve the paper.

Appendix A. Derivation of the depth-averaged governing equations
On a rectangular Cartesian coordinate system (x, y, z) with the vertical axis z parallel to the direction of gravity, the three-dimensional governing equations for the single-phase granular flows are given by where ∇ is the gradient operator, U = (u,v,w) T defines the velocity field, ρ is the bulk density, t denotes the time, ⊗ is the dyadic product, T is the Cauchy stress tensor, g = (0,0, −g) T gives the gravity field. According to Savage and Hutter (1989), the moving mass may be assumed to satisfy an internal Mohr-Coulomb frictional rheology. The Cauchy stress tensor may be decomposed into an isentropic pressure p and a deviatoric stress τ: where I is the identity tensor and τ = (τ xx ,τ xy ,τ xz ,τ yy ,τ yz ,τ zz ) contains the stress components in different directions. X. Xia, Q. Liang Engineering Geology 234 (2018) 174-191 = u v w gH u v w ( , , ) ( , , ), (A.5) in which s and b respectively denote the free surface and bed elevations, as illustrated in Fig. A.13, h is the depth of the mass flow (s = h + b), L and H are characteristic length and thickness; the corresponding aspect ratio must satisfy ε = H/L ≪ 1 to reinforce the shallow flow assumption, and the quantities with hats represent the scaled variables. Herein, the downhill flow speed is assumed to be of the same order as the gravity wave celerity gH which is different from gL as used in Savage and Hutter (1989). The net difference between the gravity force and the resistance is of a lower order compared with the gravity force itself; thus the flow velocity can only be accelerated to the order of gH but not gL , the same scaling strategy has also been used in other models (e.g. Gray and Edwards, 2014).
The scaled stress tensor components may be subsequently given as ̂̂̂̂̂= p τ τ τ τ τ τ ρgH p μτ μτ μτ μτ μτ μτ ( , , , , , , ) ( According to the Mohr-Coulomb assumption, the deviatoric tensor is friction-related. Because the coordinate is fixed but the displacement direction is arbitrary, all of the deviatoric tensor components could be friction-related at certain sitiuations, and therefore they are all scaled to O(μρgH), where μ is the friction coefficient. Replacing the variables in Eqs. (A.1)-(A.2) with the non-dimensionalised variables, we can obtain leading to the following depth-integrated equations: and σ body = (0,0, −h) T is the gravity-related body force. Fig. A.14 illustrates an example of the stress vectors acting on a soil column. The projection of Eqs. (A.19)-(A.21) onto the basal normal direction can be therefore written as 14. An example of the stress vectors acting on a soil column.
For shallow flows, it is generally assumed that the velocity is (nearly) parallel to the bed surface, the relationship between the three components of the velocity u can be written in the following formal form: In Eq. (A.28), the first term at the right hand side (RHS) indeed contains the centrifugal forces, which have been derived naturally from the above X. Xia, Q. Liang Engineering Geology 234 (2018) 174-191 depth-integrating procedure without making any ad-hoc assumptions. Now the basal normal traction (as illustrated in Fig. A.14) can be easily in which the term ⋅ σ n b in Eq. (A.25) has become part of the residual term O(ε) as evidenced from Eq. (A.23). According to the Mohr-Coulomb rheology assumption, the basal tangential stress vector acting on the bed (also illustrated in Fig. A.14) is given as where U b is the basal velocity vector. The vertical velocity component w b in U b is related to the horizontal velocity components U b and v b through A popular method to determine the lateral stress coefficient was provided by Savage and Hutter (1989). However the approach cannot be directly used here because it was originally derived on a local coordinate system with the vertical axis normal to the bed. A rotated stress tensor should be applied to obtain a new lateral stress coefficient, but the resulting formula becomes complicated. Indeed the validity of the lateral pressure coefficient introduced by Savage and Hutter (1989) is still in debate and the more recent molecular dynamic simulations (Silbert et al., 2003;GDR Midi, 2004) suggest that the lateral pressure coefficient is actually closer to unity. In this work, the lateral stress coefficient is taken as unity (i.e. K = 1) for simplicity but more complicated formula may be derived and adopted in the future. The unity lateral stress coefficient has also been adopted by many other researchers to develop their models (e.g. Bouchut and Westdickenberg, 2004;Pudasaini and Hutter, 2003;Gray et al., 1999). Substituting