Modelling rapid mass movements using the shallow water equations in Cartesian coordinates

We propose a new method to model rapid mass movements on complex topography using the shallow water equations in Cartesian coordinates. These equations are the widely used standard approximation for the flow of water in rivers and shallow lakes, but the main prerequisite for their application – an almost horizontal fluid table – is in general not satisfied for avalanches and debris flows in steep terrain. Therefore, we have developed appropriate correction terms for large topographic gradients. In this study we present the mathematical formulation of these correction terms and their implementation in the open-source flow solver GERRIS. This novel approach is evaluated by simulating avalanches on synthetic and finally natural topographies and the widely used Voellmy flow resistance law. Testing the results against analytical solutions and the proprietary avalanche model RAMMS, we found a very good agreement. As the GERRIS flow solver is freely available and open source, it can be easily extended by additional fluid models or source areas, making this model suitable for simulating several types of rapid mass movements. It therefore provides a valuable tool for assisting regional-scale natural hazard studies.


Introduction
Rapid mass movements such as avalanches, debris flows, and lahars are globally abundant surface processes in steep mountainous areas characterized by high process velocities and large masses of granular material involved (e.g.Kirschbaum et al., 2010).Therefore, rapid mass movements represent first-order threats whenever their process domains intersect with populated areas and infrastructure for transport (streets, railway lines), energy supply (power plants, pipelines, electricity lines), or tourism (e.g.ski resorts).While most of the villages prone to rapid mass movements have already implemented hazard mitigation strategies by a combination of permanent and temporal preventive measures, the latter are progressively developed in remote mountainous regions where historical records on rapid mass movements are sparse, so that the level of threat is ambiguous.In combination with field mapping and remote sensing, numerical models describing the motion of granular material on general topography are the primary tool to evaluate the potential impact of rapid mass movements on infrastructure in these remote places.Run-out distance, flow and depositional depth, velocity, and momentum of a granular flow resulting from physically based numerical models represent key parameters to (a) delineate hazard zones on the regional scale, (b) locate ideal corridors and construction areas for new infrastructure, and (c) develop mitigation strategies for protecting planned and already existing infrastructure against these natural hazards (Hsu et al., 2010;Keiler et al., 2006).To fulfill these tasks, codes have to be equipped with advanced numerical techniques to reach the required computational performance (e.g.adaptive mesh refinement), have to provide an interface to geographic information systems (GIS) and should be controllable by a scripting language to perform Monte Carlo simulations and parameter studies for entire valleys and hundreds of process domains.
Several state-of-the-art codes describe granular flow on general topography but are either not open source (e.g.FLATModel, Medina et al., 2008) or restricted to simple Coulomb-type rheology (e.g.Titan2D, Sheridan et al., 2005).
Published by Copernicus Publications on behalf of the European Geosciences Union.S. Hergarten and J. Robl: Modelling rapid mass movements using the shallow water equations The implementation of a new rheology model even in opensource scientific codes by the user is in general not practicable, because a deep knowledge on the specific code and fluid dynamics is required.Other recent codes such as Flow-R (Horton et al., 2013) replace the equations of continuum mechanics by more empirical, grid-based algorithms and thus require a higher degree of calibration.
The Voellmy rheology (Voellmy, 1955) is frequently used to describe debris flows and dense snow avalanches and is implemented in proprietary software products such as RAMMS (Christen et al., 2010), SAMOS-AT (Sampl and Zwinger, 2004;Sailer et al., 2008;Granig and Jörg, 2012), and ELBA+.The latter has been extensively used by the Austrian avalanche and torrent control, but peer-reviewed publications on technical details are still missing.Astonishingly, there are no open-source codes that fulfill all requirements mentioned above to describe granular flow with a Voellmy rheological model on complex topography.However, there are several open-source packages for a wide range of fluid dynamical problems that provide state-ofthe-art flow solvers and a variety of numerical accessories like automatic meshing routines or adaptive mesh refinement, such as OPENFOAM (Weller and Weller, 2008, http: //www.openfoam.com),CLAWPACK (LeVeque et al., 2011;Berger et al., 2011, http://clawpack.github.io),and GERRIS (Popinet, 2009, http://gfs.sourceforge.net).Besides many other applications, these codes are routinely used to predict the propagation of tsunamis in ocean basins (Popinet, 2012) or to model the extent of inundation areas during flooding (An and Yu, 2012) by solving the nonlinear shallow water equations being the standard approximation for the flow of water in rivers and shallow lakes.These numerical packages are operated by highly flexible parameter files that allow the implementation of new fluid rheology models without writing additional source code, so that it is tempting to describe rapid mass movements with one of these fluid dynamics software packages.
In this spirit, two-dimensional models for rapid mass movements on a given topography are similar to the shallow water equations.In both concepts, vertically averaged velocities are considered, and the rheology of the medium and effects of turbulence are taken into account in form of a friction term depending on flow depth and velocity.However, the widely used shallow water equations are only applicable if the water table is almost horizontal.This condition is in general not satisfied for mass movements in steep terrain, so that more elaborate approaches are required here.
These approaches can be subdivided into two major classes according to the coordinate system used.The different coordinate systems are illustrated in Fig. 1a for the example of a channel not parallel to any of the Cartesian coordinate axes x, y, and z.The first group of models focuses on flow in a given channel and uses a curvilinear coordinate system with a z axis being always normal to the surface, while the x coordinate follows the so-called thalweg in downslope direction (red).A formulation for granular flow in general curved and twisted channels was provided by Pudasaini and Hutter (2003) and later extended and applied to data (Pudasaini et al., 2005(Pudasaini et al., , 2008)).Recently, an implementation of this concept called r.avalanche in the open-source GIS software suite GRASS was presented by Mergili et al. (2012).However, it imposes significant simplifications on the thalweg concept, in particular that it is a straight line in map view, so that it is more suitable for flow on slopes than in predefined channels.The alternative concept also assumes that one coordinate, z , is normal to the surface, while the horizontal projections of the x and y coordinates approximately follow the original Cartesian axes x and y (blue).The software RAMMS (Christen et al., 2010) implements the simplest version of such a local coordinate system by neglecting the surface curvature.Potential limitations arising from this approximation were discussed by Fischer et al. (2012), presenting an extension taking the surface curvature into account for the price of more complicated differential equations.A general formulation for arbitrary coordinate systems was provided by Bouchut and Westdickenberg (2004).
In this paper we introduce a different approach based on the original shallow water equations in Cartesian coordinates.As these equations are only valid for an almost horizontal fluid surface, appropriate corrections must be applied to both acceleration and friction for large slope gradients.We express both corrections in terms of a modified friction term.In the following section, an expression for this friction term is derived, and in Sect. 4 the approach is validated by comparing several scenarios with the established avalanche model RAMMS.

Theory
The shallow water equations (e.g.Vreugdenhil, 1994) provide a two-dimensional approximation for the flow of water (or any liquid or granular medium).They refer to vertically averaged horizontal velocities and assume an almost horizontal water table, so that the vertical component of the velocity can be neglected, and the vertical pressure distribution is hydrostatic.Under these conditions, the horizontal pressure gradient and thus the horizontal acceleration is proportional to the gradient of the water table.In their native form directly referring to the acceleration, the shallow water equations read with and The model variables are the vertically averaged horizontal velocity v h (a two-component vector) and the vertical (not normal to the surface) flow depth h v .Both variables depend on the spatial coordinates x and y and on time.The symbol ∇ denotes the two-dimensional gradient operator, and z b (x, y) is the topography, so that s is the negative gradient of the water table H (x, y).The parameters g and ρ are the gravitational acceleration and the bulk density respectively.The second term on the right-hand side is a friction term in direction opposite to the velocity.Here it is written in terms of a basal shear stress τ , but this does not imply that friction in fact only occurs at the bottom of the fluid layer.For turbulent flow of water, τ is proportional to the square of the velocity, but arbitrary functions involving velocity and flow depth may be used.
In the literature, the shallow water equations are often transformed to the so-called advective form where div is the vector divergence operator and v T h is the transpose of v h (a row vector instead of a column vector).The vector h v v h is often termed momentum in this context.The second term on the left-hand side of Eq. (4) describes the advective transport of momentum with the transport velocity v h , while the right-hand side can be interpreted as a source of momentum.
For theoretical and numerical considerations, a third form of the shallow water equations is frequently used.It is de-noted conservative form and reads where I is the identity matrix.However, the native form directly referring to accelerations (Eq. 1) is more convenient for our purposes.
In each case (Eqs. 1, 4 or 5), the shallow water equations must be combined with the equation of continuity, describing the conservation of volume.
If the gradient of the water table is large, the corresponding acceleration term in Eq. ( 1) overestimates the real acceleration for two reasons: (i) the real acceleration acts in direction parallel to the surface, while Eq. ( 1) involves only its horizontal component; and (ii) the absolute value of the gradient of the water table, |s|, corresponds to the tangent of the slope angle ϕ of the fluid surface, while the downslope acceleration on an inclined plane is in fact proportional to sin ϕ.This effect is also taken into account in all models based on curvilinear coordinate systems discussed above.Compensation of each of these errors requires a multiplication of the gravitational acceleration by a factor cos ϕ, so that the acceleration term in Eq. (1) has to be reduced by a factor cos 2 ϕ in total.
The friction term also requires a correction for large gradients, namely a multiplication by cos ψ, where ψ is the inclination angle of the velocity (Fig. 1b).This angle is in general smaller than ϕ and only equal to it for flow in downslope direction.It is given by Furthermore, the vertical flow depth h v must be replaced by the flow depth normal to the surface that is by a factor cos ϕ smaller.Returning to the vertical flow depth h v requires the division of the friction term by cos ϕ.
With these three modifications to the right-hand side, the shallow water equations turn into As our approach shall be compatible with the original shallow water equations, the acceleration term shall remain linear, so that the reduction of the acceleration must be mimicked by an additional friction term: However, the original shallow water equations only allow a friction term in direction of the velocity.Therefore we only consider the projection of the friction term on the velocity: while the component normal to the flow direction is neglected.With this approximation, Eq. ( 9) turns into In the context of dense snow avalanches, the Voellmy rheology (Voellmy, 1955) is the most frequently used constitutive law for the friction term.It combines a velocity-independent Coulomb friction term with a term proportional to the square of the velocity as it is mostly used for turbulent flow: Here, σ denotes the normal stress at the bottom of the fluid layer.Assuming that the flow bed is parallel to the fluid surface, it is given by Otherwise, ϕ in Eq. 16 and thus one of the factors cos ϕ in Eq. 17 should refer to the flow bed instead of the fluid surface.We will come back to this aspect in Sect.3. The second term on the right-hand side of Eq. ( 15) is independent of the direction of the velocity, so that we obtain Inserting this expression in Eq. ( 14) finally yields with This set of equations differs from the original shallow water equations (Eq. 1) only by the more complicated friction term.
For considerations based on the conservative form (Eq. 5) it is readily transformed to with the same function f (h v , v h , ϕ, ψ).

Implementation
Our approximation can easily be implemented in any continuum fluid dynamics software which is able to solve the shallow water equations for a given bed topography and allows the implementation of arbitrary friction terms.We use the software package GERRIS (http://gfs.sourceforge.net)which is freely available and has been in development for more than 10 years.It provides highly developed numerics, and applications of GERRIS have been presented in numerous publications.
GERRIS uses the conservative form of the shallow water equations where h v and the momentum q = h v v h are the variables.Arbitrary friction terms such as the one in Eq. ( 21) can be implemented by operator splitting.The time step from t to t + δt is split up into two half steps.In the first half step, an interim solution q is computed by solving the shallow water equations without friction.In the second half step, the "real" momentum at the time t + δt is computed from q by applying the friction term only, i.e. by solving the differential equation where the solution at the time t is the interim momentum q.
As this equation does not contain any spatial derivatives of the momentum, it is degenerated to a set of ordinary differential equations.Furthermore it does not alter the direction of q, so that it is in principle even scalar.Applying a mixed Euler scheme with an explicit discretization of the arguments of f and an implicit discretization to the remaining term q, i.e.
yields the solution The angles ϕ and ψ should be computed from the gradient of the surface of the flowing medium (Eq.2) according to Eqs. ( 7) and ( 8) in each time step.Strictly speaking, the angle ϕ occurring in the term µ cos ϕ cos ψ in Eq. ( 20) should refer to the flow bed as discussed in Sect. 2. However, the shallow water equations are in principle only valid as long as the gradient of the flow depth is small, i.e. as long as the flow surface is almost parallel to the bed topography.We therefore suggest an implementation where all angles are computed from the fluid surface H , denoted GERRIS H in the following.
However, computing all angles from the bed topography z b should not make a fundamental difference.The respective implementation is denoted GERRIS Zb in the following.In this context it should be kept in mind that the difference between both variants only concerns the corrections for large gradients.In both implementations as well as in the original shallow water equations and in the general formulation of Bouchut and Westdickenberg (2004), the acceleration due to gravity being the driving force depends on the gradient of the fluid surface, ensuring that the fluid table of a lake is horizontal.So the difference should become significant only for large gradients of the fluid table in combination with large variations in flow depth.As illustrated in Sect.4.1, using GERRIS H may even generate artefacts when this gradient is not computed from the same discretization scheme that is used internally for solving the shallow water equations, which will likely occur when staggered grids or finite volume discretizations are used.Due to this aspect it may be even advisable to use the implementation GERRIS Zb .

Validation
In this section we compare our approximation to the established model RAMMS (version 1.5.01).We use the simplest version based on Voellmy's rheology, neglecting entrainment (Christen et al., 2010) and do not use the recently introduced features of extending Voellmy's rheology by cohesion and taking into account the effect of surface curvature on the frictional force considered by Fischer et al. (2012).However, all these extensions can in principle be adjusted to our formulation based on the shallow water equations in Cartesian coordinates.
Compared to the reference model RAMMS as defined above, our approach introduces two approximations.The most serious one consists of considering only the horizontal component of the velocity.While the accelerations due to the slope gradient and due to friction are corrected accordingly, the horizontal velocity would remain constant in absence of gravity or friction.As a consequence, the total velocity increases artificially on a convex slope and decreases on a concave slope even without gravitational acceleration.The other simplification concerns the projection of the correction terms on the velocity vector (Eq.12).This means that the longitudinal acceleration (i.e. in flow direction) is corrected appropriately for large slope gradients, while the transversal acceleration is directly adopted from the original shallow water equations.
In the following we investigate three scenarios defined with regard to these approximations.In the first example, flow down a planar slope is considered.This scenario should be described well by both RAMMS and by our approach.The second set of tests refers to slopes with a strongly curved part in order to examine whether the first approximation has a serious effect.Finally, we consider a more complex topography as an example closer to real-world applications.

Constant flow depth on a planar slope
The movement of an avalanche with a constant flow depth on a planar slope in one dimension can be described by an analytical solution (e.g.Pudasaini and Hutter, 2007).For this purpose we use the velocity, v, parallel to the slope and Lagrangian coordinates, which means that v is the velocity of a given particle and not at a given position.Then the equation of motion is the same as for a rigid body: where τ is the frictional shear stress.According to the arguments leading from Eq. ( 15) to Eq. ( 18), this shear stress amounts to so that The steady-state solution of this equation (i.e. the asymptotic velocity v ∞ ) is readily obtained by setting the left-hand side to zero: With this, Eq. ( 27) turns into The time-dependent solution of this equation is with the characteristic time describing how slowly the velocity approaches v ∞ .
To test whether our approach reproduces this behaviour correctly, we consider a planar ramp with ϕ = 30 • inclination in x direction with Voellmy parameters µ = 0.2 and ξ = 1000 m s −2 .The release zone is defined by a rectangular area of 350 m × 400 m (in horizontal projection) at the upper edge of the ramp with a release height of h = 1 m measured normal to the topography.Figure 2a shows the topography and the flow depth after 20 s, obtained from the simulation with GERRIS H where the frictional terms (i.e.ψ in Eq. 21) are computed from the surface of the flowing mass.numerical experiments: the realization GERRIS H discussed above, the alternative approach GERRIS Zb where the bed surface is used to compute the friction term, and the reference model RAMMS.Only minor differences among the three models are encountered.The avalanche develops a characteristic tail with a rapidly declining flow depth in upslope direction, while the initial flow depth of h = 1 m is still preserved in the main body.The avalanche has already reached the steady-state velocity of about 18 m s −1 in the main body predicted by Eq. ( 28), while the velocities in the tail are lower as a consequence of the reduced flow depth.The avalanche front of all three experiments is characterized by a slight increase of flow depth and flow velocity relative to the main avalanche body.This artefact is in general small, but most pronounced for GERRIS H , while RAMMS and GERRIS Zb show nearly identical profiles at the avalanche front.The slightly stronger artefact occurring in GERRIS H presumably arises from our simple implementation of the gradient of the fluid surface required for computing the angles ϕ and ψ required in Eq. ( 24).
Here we use the standard gradient of the fluid surface provided by GERRIS that is computed from simple symmetrical difference quotients.The sophisticated numerics implemented in GERRIS itself used for maintaining a sharp front is not incorporated here, so that finally the driving term of the shallow water equations and the friction term use different schemes of discretization, causing artefacts at the avalanche front where the fluid surface is strongly curved.However, we found in all considered examples that these small artefacts are stable and do not grow over time, so that they are not a serious problem at all.As a second test, we consider the velocity of the accelerating fluid layer against the time-dependent analytical solution (Eq.30) for different initial flow depths (h), turbulence (ξ ) and dry friction (µ) parameters of the Voellmy rheology, and hillslope angles (ϕ) (Fig. 3).Similarly to the results on the avalanche profiles, the almost perfect agreement between the velocity predicted by Eq. ( 30) and all sets of numerical experiments verifies the ability of our approach at least for planar slopes.Small deviations occur shortly after the release scale with the time step size of the numerical model and could be reduced by forcing the flow solver towards smaller time increments.However, these initial small deviations disappear rapidly when approaching the terminal velocity, so that a higher temporal resolution at the expense of increasing computational time does not justify this insignificant benefit in practical applications.

The effect of profile curvature
While the tests performed in the previous section only concern the technical correctness of the theory and its implementation, the following numerical experiments address the validity and the limitations of the approximations made.
In order to explore the effects of profile curvature on our approach considering only the horizontal components of the velocity vector (and calculating the total velocity from those), we have performed a series of numerical experiments on curved synthetic topographies and compare the results of our approach with those of RAMMS.The first experiment describes an avalanche on a concave flow path defined by a 30 • dipping ramp and a 5 • inclined run-out zone with a smooth (parabolic) transition between both (Fig. 4).Here and in the following, the curvature of the smooth transition zone is defined in such a way that an avalanche entering from the upper ramp with the terminal velocity according to Eq. ( 28) is exposed to an centrifugal acceleration of about 1 m s −2 (which is neglected in both RAMMS and our approach but considered in detail by Fischer et al., 2012).
The behaviour when travelling along the upper ramp is the same as in the example considered in Fig. 2. At t = 40 s, the avalanche is characterized by a long tail, while the bulk mass of the avalanche remains undeformed and has reached the steady-state velocity of about 18 m s −1 (Fig. 4c, d).At this time the avalanche front approaches the curved transition zone to the gently dipping run-out zone leading to a strong deceleration and thickening.At t = 60 s, the frontal part of the fluid layer is more than 3 times thicker than it initially was.The flow velocity has decreased below 10 m s −1 everywhere (Fig. 4e, f).At t = 80 s, the avalanche thickens further in the run-out zone and grows laterally normal to the flow path and in upslope direction as additional mass from the slower avalanche tail becomes incorporated in the deposit (fluid without significant motion).At this stage, significant flow velocities are confined to the steep flow path section where the avalanche tail is still in motion (Fig. 4g, h).While Fig. 4 shows that the avalanche behaves as expected qualitatively, Fig. 5 provides a quantitative comparison with the reference model RAMMS.The flow depths (Fig. 5a-c) and the velocities (Fig. 5d-f) predicted by both models agree almost perfectly in the domain interesting for hazard assessment, i.e.where the avalanche moves at a significant velocity.The same applies to the run-out distance when the avalanche front finally comes to rest.
Noticeable differences between the results of GERRISbased models and RAMMS only occur in the final phase of the avalanche when the front has almost come to rest.While the main body of the avalanche is characterized by a single maximum in the thickness in the GERRIS-based simulations, RAMMS predicts a bimodal avalanche profile.This difference is also reflected in the shape of the final deposits at least when the RAMMS simulation is stopped automatically using the default settings (i.e. when the momentum has decreased to 5 % of its maximum value).However, the difference arises in a phase where only the long tail of the avalanche moves at a considerable velocity, so that material is pushed from behind on the main avalanche body that it already almost resting.So this difference is probably not related to the different way of treating velocities but rather to the different numerical schemes used in RAMMS and GERRIS; apart from this, it is unimportant for practical purposes.
In Fig. 6, the opposite situation involving a convex topography is considered.Again the release zone is located on a 30 • steep slope, but in contrast to the previous experiment the slope steepens in a smooth transition to 45 • .The transition from ϕ 1 = 30 • to ϕ 2 = 45 • causes the fluid layer to accelerate rapidly from 18 m s −1 to the new terminal velocity of 21.7 m s −1 at a reduced flow depth of 0.83 m.As an effect of our approximation considering only the horizontal components of the velocity vector, the new terminal velocity is reached slightly earlier by GERRIS than by RAMMS.This effect becomes more pronounced by sharp terrain transitions (Fig. 7).In this example, the fluid moves at the terminal velocity of v ∞ = 18 m s −1 in the upper region, corresponding to a horizontal velocity v h = v ∞ cos ϕ 1 = 15.6 m s −1 .As the horizontal velocity persists at the transition in our ap- proach, the velocity immediately after the transition amounts to v = v h cos ϕ 2 = 22 m s −1 , which is even slightly above the new terminal velocity.However, the avalanche of RAMMS also reaches the terminal velocity rapidly after entering the steeper region for both the smooth and the sharp transition.
A quantitative estimate on the range where our approximation affects the flow velocity after the slope has changed can be derived from the theoretical considerations made in Sect.4.1.Instead of the velocity as a function of time, we now consider the velocity as a function of the travelled distance.Using Eq. ( 29) we obtain Equation ( 33) implies that v approaches the terminal velocity v ∞ exponentially with decay length In the example considered above, L amounts to about 42 m, so that the avalanche indeed needs a very short travelling distance to approach the terminal velocity.Considering the more realistic situation of an steady-state avalanche with a constant influx instead of a constant flow depth, i.e. hv = const, leads to basically the same result where only the factor 2 in the denominator turns into a factor 3. Thus, the avalanche should practically approach the terminal velocity even more rapidly than stated above.Returning to Fig. 6, the only significant difference between the results of the two GERRIS-based models and RAMMS occur at a late stage (t > 80 s) where the direct effect of our approximation should have almost vanished.As the avalanche becomes more and more stretched at the backward side, the region with a constant thickness becomes shorter until it finally vanishes.In the simulation with GERRIS Zb , this leads to a rapid decrease in flow depth and consequently in velocity, so that the avalanche decays rapidly.In contrast, RAMMS keeps a sharper distinction between the region of constant flow depth and consequently maintains the original flow depth and velocity for a longer time.In return, rather strong waves occur at the transition to the tail being visible in both flow depth and flow velocity.Such waves are in principle generated by GERRIS Zb , too, but with a significantly smaller amplitude than in RAMMS.In return, GERRIS H generates even stronger oscillations than RAMMS.
Similarly to the differences among the models found for the concave topography, the differences found here are presumably not related to our approximation of considering only the horizontal velocity and computing the total velocity afterwards but rather to the different discretization schemes used in RAMMS and GERRIS.GERRIS itself obviously uses a numerical scheme that is well suited for reducing oscillations at the transition to the avalanche tail, but our simple implementation of the gradient of the fluid surface used in GERRIS H cannot compete with this scheme.In summary, the numerical experiments with GERRIS and the comparison with the leading avalanche model RAMMS performed in this section demonstrate the ability of our approach to model avalanches even on curved topography.The effects of our approximations cause only minor deviations, and in particular their impact on predictions of run-out distance, flow depth, and velocity is practically negligible.The better stability of both the avalanche front and the transition to the tail provides arguments in favour of GERRIS Zb compared to GERRIS H .

Flow over complex topography
The thalweg of rapid mass movements on a real topography is in general curved and twisted.We therefore challenge our approach with the complex topography of a typical alpine avalanche flow path and test the results of our approach against RAMMS.In contrast to the previous examples that are in their basic structure one-dimensional, the second approximation made in our theory also becomes relevant here: beyond considering only the horizontal component of the velocity in the equations, our approach only applies corrections for large slopes to the longitudinal component of the acceleration.
The hypothetical avalanche is located in the Felbertal, a typical glacially shaped alpine valley with large open flanks between ridges and the tree line representing classic snow avalanche release zones.Deeply incised, curved and twisted gullies canalize the avalanche in one or several branches with locally extreme flow depths.These gullies route the granular fluid to the nearly flat valley floor, representing the run-out zone of the avalanche (Fig. 8).
We compare the maximum values (at each point, taken over the entire simulation) of flow depth, momentum, and velocity of the modelled avalanche.We set the release height to h = 1 m and define spatially constant parameters for the Voellmy flow resistance law (µ = 0.2, ξ = 2000 m s −2 ).All simulations are performed on a quadrilateral grid with a spatial resolution of 4 m and terminate when the momentum of the fluid drops below 5 % of the maximum momentum (default of RAMMS).
Generally, the deviations among GERRIS Zb , GERRIS H , and RAMMS are small, and the primary features of the avalanche agree well between the two GERRIS approaches and RAMMS.This includes flow depths, run-out distances, flow velocities, and momentum.
However, a closer examination reveals some deviations between the different numerical approaches.RAMMS shows a more pronounced tendency to overflow counter hillsides and to keep the flow direction even uphill.This is clearly documented at the lower third of the avalanche track where the avalanche is split into two branches.Here the orographic right flow path is characterized by a considerable uphill section.In this domain, the results of RAMMS show higher values in the maximum flow depth compared to the two GERRIS approaches (Fig. 8a-c).The modelled avalanches in RAMMS overflow larger areas, causing a wider flow path than predicted by the GERRIS experiments.This is recognized most clearly in the s-shaped gully section.The broader flow path and tendency to flow uphill, observed in avalanches modelled with RAMMS relative to those modelled with GERRIS Zb and GERRIS H , are caused by larger values in the momentum (Fig. 8d-f) arising from slightly higher flow velocities especially in the gully section of the thalweg (Fig. 8g-i).
In contrast to the small deviations found in the previous examples, this effect is a direct consequence of applying only corrections to the longitudinal acceleration.When an avalanche follows a narrow and strongly curved (in map view) gully, the transversal (centripetal) acceleration preventing the fluid from leaving the gully is overestimated by the shallow water equations, similarly to the longitudinal acceleration.However, in contrast to the longitudinal acceleration, the overestimation of the transversal acceleration is not corrected, so that the tendency of the avalanche to stay within the gully is stronger than in RAMMS.A further deviation is observed in the run-out zone where the shapes of the avalanche deposit of the RAMMS simulations differ considerably from those of the GERRIS simulations, while the run-out distances and base areas of the avalanche deposits of the three models agree well.The discrepancy in the geometry of the deposits is similar to that found when considering the run-out on the simple concave topography in Sect.4.2, and it is presumably related to the different numerical schemes used in RAMMS and GERRIS rather than to our approximations.
Finally, the differences between the realizations GERRIS H , where the corrections in the friction term are based on the fluid surface, and GERRIS Zb , where the corrections are computed from the original topography, are very small in this example.So the stronger (although not serious) artefacts occurring at the avalanche front in GERRIS H (Sect. 4.1) and the higher stability of GERRIS Zb at the transition to the tail remain the only noticeable difference between both GERRIS-based approaches.These differences suggest that using the version GERRIS Zb may in general be preferable to GERRIS H .

Conclusions
The examples considered in Sect. 4 show that granular avalanches can be simulated using the shallow water equations directly in Cartesian coordinates even in steep terrain when an appropriate additional friction term is included.This finding allows the utilization of software that was originally designed for other purposes, namely modelling the flow of water in rivers, lakes, and oceans.
Compared to software packages explicitly developed for modelling avalanches, a wealth of state-of-the-art fluid dynamics software packages potentially being adjustable for this purpose is available.Some of them are even freely available.Therefore, research on avalanches can easily profit from the enormous effort that has already been spent in developing numerical codes in fluid dynamics.The implementations presented in this paper are based on the software GERRIS, but this should only be seen as an example.Apart from being freely available and providing state-of-the-art numerics, GERRIS allows the implementation of our method with a moderate effort.However, this should not imply that GER-RIS is the best software for this purpose.
The examples investigated for validation have only revealed minor deviations from the proprietary model RAMMS used as a reference, in particular with regard to the properties of avalanches relevant for hazard assessment.For the artificial, basically one-dimensional geometries investigated in Sect.4.1 and 4.2, the agreement between our implementations based on GERRIS and RAMMS is excellent.The small differences between the approaches encountered here are probably not related to the approximations introduced in our theory but arise from different numerical schemes used for solving the equations of motion.So the approach presented in this study can fully compete with proprietary software for mass movements flowing on open flanks or if large volumes and high flow depths occur and small gullies do not influence the flow characteristics strongly.For such geometries, the deviations between the results of our approach and those of RAMMS are far off from having any implications on the mitigation strategy based on predicted properties of the modelled avalanche.The example based on a complex topography (Sect.4.3) reveals still rather small, but perhaps not always negligible, differences between our approach and RAMMS.RAMMS solutions show larger flow depths in avalanche tracks with prominent uphill sections and expanded overflowed areas in steep, twisted gullies.In contrast to the differences discussed above, these deviations arise from the approximations discussed in Sect. 3 and represent a small intrinsic model limitation that is inevitable when using the shallow water equation in a Cartesian coordinate system with a friction term acting only in direction opposite to the velocity.
However, when discussing differences between models on such a small level we should keep in mind that all these modelling approaches involve a considerable inherent uncertainty compared to other flow processes such as the flow of water in lakes and oceans.These uncertainties start with the basic assumption of the granular medium as a single layer continuum and the rheology (e.g.Voellmy's friction law).They continue with the determination of the relevant parameters for dry snow avalanches and do not stop at the determination of the release zone in the form of spatial position, extent, and involved volumes (fracture depth).Even the resolution and quality of the applied digital elevation model can highly influence the avalanche path (Bühler et al., 2011), and taking into account further processes such as entrainment introduces an additional uncertainty in the parameters.Assessing these uncertainties quantitatively goes beyond the scope of this paper, but, in summary, they are obviously larger than the small deviations between the models.
The differences between the two proposed implementations based on GERRIS are also small.The version GERRIS Zb where the correction terms are computed from the original topography is less prone to artefacts at the avalanche front and at the transition to the tail than the version GERRIS H using the fluid surface, without revealing significant drawbacks anywhere.We therefore suggest computing the friction terms from the topographic slope instead of the fluid surface.
The approach is suitable to describe snow avalanches flowing on general topography but is also applicable to other rapid mass movements such as debris flows.Debris flows are characterized by lower flow velocities and lower flow path gradients compared to snow avalanches, so that effects of our approximations become even less significant.Besides model set-ups with a pre-defined release volume, a huge number of scenarios involving different release zones characterized by discharge time series can be easily implemented within the GERRIS parameter file.In principle the initiation of surface runoff can be defined at each mesh element, so that flooding and debris flow simulations based on precipitation time series for storm events are possible without preceding precipitation runoff models.Flow resistance laws and their parameterizations are also defined in the parameter file, so that the implementation of other rheological models (e.g.Bingham fluid) is straightforward and requires no specific coding skills.
We propose that our approach based on GERRIS is suitable for regional-scale dense snow avalanche studies on complex terrain and probably also for other types of rapid mass movements.However, dimensioning of permanent protection measures requires numerical models that have been calibrated by the backward analysis of numerous monitored real world avalanches as, for example, performed by the SLF at the Vallée de la Sionne.In principle, the parameters of the Voellmy fluid model calibrated for RAMMS are fully compatible with our approach, and spatially variable parameter values can be easily implemented in the GERRIS parameter file.This also applies to extensions of the flow law such as the cohesion term implemented in the recent version of RAMMS.Thus, a more or less complete compatibility with RAMMS can be achieved.However, as we did not perform backward analysis calculations to calibrate the fluid model for our approach and we tested the compatibility only for a few examples, the modelling results should be taken with caution when mitigation strategies and protection measures are developed.For such applications, the compatibility with established and extensively tested software packages should be ascertained for the given conditions, or at least a careful backward analysis of the specific avalanche should be conducted.
Figure 1.(a)Different coordinate systems used for modelling mass movements in a channel not parallel to any of the Cartesian coordinate axes.Green: Cartesian coordinates x, y, and z (this study).Blue: Coordinates aligned to the surface (x and y parallel to the surface, z normal to the surface).Red: Coordinates aligned to the surface (x, y, z) where the x axis follows the thalweg (dashed red line).(b) Illustration of the slope angle ϕ of the fluid surface (Eq.7) and the inclination angle ψ of the velocity (Eq.8).

Figure 2 .
Figure 2. Comparison of two different numerical solutions of GERRIS with RAMMS for granular flow over a 30 • dipping ramp for Voellmy parameters µ = 0.2 and ξ = 1000 m s −2 and an initial flow depth h = 1 m .(a) Two-dimensional representation of the ramp geometry and the flow depth of the fluid layer measured normal to the surface after 20 s.The white dashed-dotted line indicates the position of the longitudinal profiles shown in (b).(b) Longitudinal profiles of three different numerical models showing flow velocity (solid lines) and flow depth (dashed-dotted lines).The gray dashed line indicates the theoretical terminal velocity according to Eq. (28).GERRIS-based models (blue and green lines) apply gradients of fluid surface (GERRIS H ) and gradients of the topography (GERRIS Zb ) respectively, and black lines are results RAMMS as reference.

Figure 3 .
Figure 3. Test of the numerical results of GERRIS Zb (symbols) against the one-dimensional analytical time-dependent solution for a Voellmy fluid flowing on a flow path with a constant slope (Eq.30, solid lines) for several parameter sets for (a) flow depth h, (b) turbulence parameter ξ , (c) dry friction parameter µ, and (d) slope angle ϕ.

Figure 4 .
Figure 4. Two-dimensional representation of numerical solutions of GERRIS Zb for flow depth and flow velocity partly concave slope at four different time steps.The fluid layer accelerates at an inclined ramp (ϕ 1 = 30 • ) and runs out in a gently dipping surface (ϕ 2 = 5 • ) with a smooth (parabolic) transition starting at x = 0.A topographic profile of the geometry is shown in the inset of (a).The white dashed-dotted line indicates the position of the longitudinal profiles shown in Figs. 5 and 6.

Figure 5 .Figure 6 .
Figure 5.Time series of longitudinal profiles plotted every 5 s for flow depth and flow velocity of the scenario shown in Fig. 4. Profiles are based on numerical solutions of (a, d) GERRIS Zb (b, e) GERRIS H , and (c, f) RAMMS.The insets show the topographic profile.

Figure 7 .
Figure 7. Time series of longitudinal profiles plotted for the scenario considered in Figs. 6 and 7 but with a sharp transition between the planar slope segments.

Figure 8 .
Figure 8. Maximum values of the characteristic properties of a hypothetical avalanche flowing along a curved and twisted thalweg on a real topography (Felbertal, Austria), based on the numerical solutions of GERRIS Zb , GERRIS H , and RAMMS for the same model set-up.Fluid rheology (µ = 0.2, ξ = 2000 m s −2 ) and spatial resolution for the three different numerical models are the same.(a-c) maximum flow depth, (d-f) maximum momentum, (g-i) maximum flow velocity.