MinVoellmy v1: a lightweight model for simulating rapid mass movements based on a modified Voellmy rheology

The Voellmy rheology has been widely used for simulating snow avalanches and also for rock avalanches. Recently, a modified version of this rheology was proposed. While the conventional version of Voellmy's rheology uses the sum of Coulomb friction and a velocity-dependent friction term, the modified version assigns the two terms to different regimes of velocity. The software MinVoellmy presented here provides the first numerical implementation of the modified rheology in a two-dimensional, depth-averaged model. It consists of MATLAB and Python classes, where simplicity and parsimony were the design goals. In contrast to the majority of the models in this field, MinVoellmy uses a Cartesian coordinate system with the thickness of the fluid measured vertically and the velocity averaged vertically instead of perpendicularly to the bed. Furthermore, MinVoellmy implements a simple upstream scheme, which turns out to be sufficient for rheologies of the Voellmy type. Numerical tests reveal that the modified Voellmy rheology reproduces the empirical relation between runout length, height drop, and volume of large rock avalanches fairly well. Furthermore, there seems to be a large potential for further research on hummocky deposit morphologies and longitudinal striations.

1 Introduction

Modeling of rapid mass movements was pushed strongly by the ideas of Savage and Hutter (1989), who extended the shallow-water equations towards granular media. The shallow-water equations provide a two-dimensional, depth-averaged description of flow processes with a free surface. In their original form, the shallow-water equations assume that the bed and the free fluid surface are almost horizontal. However, this is not the case for typical scenarios of granular flow. In order to overcome this limitation, the Savage–Hutter model provides an extension towards thin layers on gently curved surfaces. Beyond the general formalism, the Savage–Hutter model also includes an approximation for the stresses arising from internal deformation of a medium with a given angle of internal friction.

The idea behind the Savage–Hutter model is adopted by almost all two-dimensional continuum models of granular flow over a given topography. Existing models differ mainly concerning rheology, coordinate system, and approach to reduce numerical diffusion. An overview of some of the available models is given by McDougall (2017).

In the context of snow and rock avalanches, the rheology proposed by Voellmy (1955) is widely used. It assumes a shear stress of

$$\begin{array}{}\text{(1)}& \mathit{\tau }=\mathit{\mu }\mathit{\sigma }+{\displaystyle \frac{\mathit{\rho }g}{\mathit{\xi }}}{v}^{\mathrm{2}}\end{array}$$

at the bed. The first term describes Coulomb friction with a coefficient μ, where σ is the normal stress. The second term was adopted from the respective relation for turbulent flow of water in open channels with a rough bed, where ρ, g, and v are density, gravity, and vertically averaged velocity, respectively. The parameter ξ refers to the roughness of the bed. As detailed by Salm (1993), Eq. (1) does not imply turbulent flow in the sense of a complete mixing of the granular medium, which would be incompatible with the preservation of stratigraphy often found in deposits of rock avalanches (Dufresne et al., 2016). The second term in Eq. (1) can be interpreted as the result of converting a part of the kinetic energy of translation parallel to the bed into random particle motion (see also Hergarten, 2024).

As discussed by Hergarten (2024), Voellmy's rheology cannot predict the long runout typically observed for large rock avalanches without further assumptions. In a nutshell, long runout requires an effective coefficient of friction, $\frac{\mathit{\tau }}{\mathit{\sigma }}$, that is much lower than typical Coulomb friction coefficients μ. However, $\frac{\mathit{\tau }}{\mathit{\sigma }}\ge \mathit{\mu }$ for Voellmy's rheology (Eq. 1), which makes it incompatible with long runout unless artificially low friction coefficients μ are assumed for large rock avalanches.

Hergarten (2024) proposed a modification of Voellmy's rheology to overcome this limitation. Instead of adding the two contributions in Eq. (1), a transition between two distinct regimes of movement in the form

$$\begin{array}{}\text{(2)}& \mathit{\tau }=\left\{\begin{array}{ccc}\mathit{\mu }\mathit{\sigma }& & v<v_{\mathrm{c}}\\ \frac{\mathit{\rho }g}{\mathit{\xi }}{v}^{\mathrm{2}}& \text{for}& v\ge v_{\mathrm{c}}\end{array}\right.\end{array}$$

was assumed. While a given constant crossover velocity v_c is the simplest idea, Hergarten (2024) also developed a model for the dependence of v_c on the thickness h of the layer. This model was obtained by reinterpreting the random kinetic energy (RKE) model (Buser and Bartelt, 2009; Bartelt and Buser, 2010), which describes the supply of kinetic energy of random particle motion and its consumption. Introducing some simplifying assumptions, the relation

$$\begin{array}{}\text{(3)}& v_{\mathrm{c}}\propto \sqrt[\mathrm{3}]{\mathit{\xi }h}\end{array}$$

was obtained. This approach turned out to predict the scaling relation between volume and runout length of rock avalanches better than the version with constant v_c and is therefore used in the following.

Concerning the implementation in numerical models, numerical diffusion is typically the most serious problem. Numerical diffusion causes a progressive smoothing of sharp fronts and an artificial damping of waves. In the context of granular media, smoothing of fronts is the major problem.

Lagrangian methods are the straightforward approach to avoid numerical diffusion. In contrast to Eulerian methods, they use a coordinate system moving with the particles. In general, however, Lagrangian methods are complicated. There seem to be only two Lagrangian models in this field. While the model DAN3D (McDougall, 2006) implements the concept of smoothed-particle hydrodynamics, which is much simpler than a classical Lagrangian approach, the model AvaFrame com1DFA (Tonnel et al., 2023) additionally uses a grid in order to avoid problems at low particle densities. In turn, the vast majority of the available models uses the Eulerian approach with a fixed coordinate system.

The total variation diminishing non-oscillatory central differencing (TVD-NOC) scheme introduced by Nessyahu and Tadmor (1990) turned out to be powerful in reducing numerical diffusion without introducing strong artificial oscillations. It is implemented in several models, e.g., the comprehensive model r.avaflow (Mergili et al., 2017). In turn, however, it will be shown in Sect. 5.2 that numerical diffusion is not a huge problem in combination with rheologies of the Voellmy type. Practically, even the simple upstream scheme works reasonably well here, which allows for simple and lightweight implementations.

The simplest form of the Savage–Hutter model refers to a channel and uses a coordinate system aligned to the bed with the x coordinate in the principal flow direction. Using a curvilinear coordinate system in this spirit on an arbitrary topography is, however, not feasible. Therefore, simpler approaches are typically preferred. The model RAMMS (Christen et al., 2010) widely used in practical applications uses a coordinate system with the x and y coordinates aligned to the Cartesian axes but locally inclined to become parallel to the topography. However, this coordinate system does not only involve a profile curvature along the axes but is also non-orthogonal. The limitations arising from these properties can be overcome by introducing more or less complicated correction terms in the equations (Fischer et al., 2012).

As an alternative, some models use a Cartesian coordinate system (Bouchut and Westdickenberg, 2004; Denlinger and Iverson, 2004; Hergarten and Robl, 2015; Rauter and Tuković, 2018). Here the challenge is that the velocity at the bed must be parallel to the bed. The approach of Hergarten and Robl (2015) was even simplified in such a way that a solver for the original shallow-water equations could be used. In turn, however, only the balance of the horizontal components of the momentum was considered, which introduces a serious limitation for scenarios with a strong profile curvature.

In the following, some kind of minimum implementation of the modified Voellmy rheology (Eqs. 2 and 3) will be developed. The model uses a Cartesian coordinate system in combination with an approximation to the driving acceleration by gravity. It is designed for simple applications in research but may also be useful in teaching since the code can be fully understood with limited knowledge about numerics and is short enough to be transferred to different programming languages easily.

In turn, it is not intended to compete with comprehensive models such as r.avaflow (Mergili et al., 2017), which even includes direct coupling to a geographic information system and options for multi-phase flow (Pudasaini and Mergili, 2019). Even more important, it should not be used for operational hazard assessment. The model RAMMS widely used in this context not only has a much longer history of continuous development but also includes estimates of its model parameters based on a large number of studies, which are essential for real-world applications.

2 Coordinates and governing equations

The model MinVoellmy presented in this paper uses Cartesian coordinates with the topography of the bed b(x,y), as illustrated in Fig. 1a. The time-dependent model variables are the thickness of the mobile layer $h(x,y,t)$ and the velocity vector $\mathit{v}(x,y,t)$. As in all models based on the theory developed by Savage and Hutter (1989), v is the component of the depth-averaged velocity parallel to the bed.

Figure 1 Definition of the model variables in two dimensions (only in the x–z plane) for (a) the Cartesian approach used here and (b) a coordinate system aligned to the bed as proposed by Savage and Hutter (1989).

Using Cartesian coordinates circumvents several problems arising from non-orthogonal or curvilinear coordinate systems aligned to the topography. In turn, the treatment of the velocity is more complicated. The condition that $\mathit{v}(x,y,t)$ must be parallel to the bed requires

$$\begin{array}{}\text{(4)}& \mathit{v}\cdot \mathit{n}=\mathrm{0}\end{array}$$

with the normal vector

$$\begin{array}{}\text{(5)}& \mathit{n}=\cos \mathit{\beta }\left(\begin{array}{c}-\mathrm{\nabla }b\\ \mathrm{1}\end{array}\right),\end{array}$$

where ∇b is the two-dimensional gradient of the bed and β the slope angle ($\tan \mathit{\beta }=\left|\mathrm{\nabla }b\right|$). While this relation would allow for reducing the velocity vector to two components in the equations, it is kept as a three-component vector here, as already proposed by Rauter and Tuković (2018). As a major difference, however, the thickness is not measured normal to the bed but vertically (Fig. 1).

As illustrated in Fig. 2, considering the vertical thickness also avoids geometrical problems with the thickness normal to the bed. Since the orange lines are not parallel if the bed is curved, an exact balance of mass and momentum must take the curvature into account. Bouchut and Westdickenberg (2004) developed such an approach for general topographies and general coordinate systems (including Cartesian coordinates) without introducing additional approximations to the Savage–Hutter model. The red lines, however, illustrate that each approach with the thickness normal to the bed is geometrically limited. Here, the thickness normal to the bed is greater than the radius of curvature of the bed, which causes an intersection of the lines. In this case, it is not possible to define the thickness normal to the bed consistently.

Figure 2 Different ways to define the thickness. Blue lines describe the definition of h as the vertical thickness proposed in this study. Orange and red lines describe the thickness normal to the bed. The normal thickness becomes inconsistent if the lines intersect (red lines).

Since the Savage–Hutter model is an approximation for thin layers and the problem does not occur if the fluid layer is sufficiently thin, it may be considered irrelevant. However, as pointed out by Hutter et al. (2005), the Savage–Hutter model is often applied to situations that are formally outside the range of validity of the assumptions but still yields reasonable results. In this sense, circumventing geometrical problems for thick layers by considering the vertical thickness may be useful.

In sum, however, considering the vertical thickness is a tradeoff. It simplifies the balance of mass and momentum (Sect. 2.1) and avoids geometrical problems. In turn, we will see in Sect. 2.2 that it requires an approximation for the acceleration by gravity, so it is not fully consistent with the original Savage–Hutter approximation unless the fluid surface is either parallel to the bed or horizontal.

2.1 The balance of mass and momentum

As in all depth-averaged models, the mass balance is taken into account by balancing the fluxes into and out of volumes aligned in the same direction as h is measured (gray areas in Fig. 1). Assuming an incompressible fluid, the mass balance can be replaced by a volumetric balance in which the (constant) density ρ does not appear. The balance equation can be written as an advection equation for the thickness h:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{\partial \left(v_{x}h\right)}{\partial x}}+{\displaystyle \frac{\partial \left(v_{y}h\right)}{\partial y}}=\mathrm{0}.\end{array}$$

Since the volumes are vertical in the approach proposed here, the second and third terms in Eq. (6) refer to horizontal advection with the velocities v_x and v_y. The balance equation is the same as for the original shallow-water equations for an incompressible fluid (e.g., Vreugdenhil, 1994).

The momentum balance can also be written as an advection equation in the form

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial \left(h\mathit{v}\right)}{\partial t}}+{\displaystyle \frac{\partial \left(v_x\right(h\mathit{v}\left)\right)}{\partial x}}+{\displaystyle \frac{\partial \left(v_y\right(h\mathit{v}\left)\right)}{\partial y}}=h\left(\mathit{a}-f{\displaystyle \frac{\mathit{v}}{\left|\mathit{v}\right|}}+c\mathit{n}\right).\end{array}$$

The advection term at the left-hand side is basically the same as Eq. (6) but for the depth-integrated momentum per unit mass hv instead of the thickness h. It is identical to the respective terms in the original shallow-water equations (e.g., Vreugdenhil, 1994). However, v is not approximately horizontal (as required for the original shallow-water equations), and the coordinate system is not aligned to the bed (as in the original Savage–Hutter model). Therefore, hv must be treated as a three-component vector, and Eq. (7) consists of three scalar equations instead of two equations in the original shallow-water equations and in the Savage–Hutter model.

In contrast to Eq. (6), Eq. (7) has a nonzero source term at the right-hand side, which describes the total depth-integrated acceleration. While the left-hand side of Eq. (7) is practically the same in all depth-averaged models, the right-hand side is model-specific. In the approach proposed here, the total acceleration is separated into three parts. The first term, a, is the gravitational acceleration parallel to the bed. Its computation is basically the same as in the original Savage–Hutter model and will be explained in Sect. 2.2. As illustrated in Fig. 3, a always points downslope with respect to the fluid surface.

Figure 3 Directions of velocities and accelerations for different situations. Only situations with the fluid surface parallel to the bed are considered for a two-dimensional geometry (x–z plane). The symbols refer to the terms v, a, $-f\frac{\mathit{v}}{\left|\mathit{v}\right|}$, and cn at the right-hand side of Eq. (7).

The second term is the frictional deceleration. Its direction is always opposite to the velocity. Since $\frac{\mathit{v}}{\left|\mathit{v}\right|}$ is a unit vector, f is the absolute value of the deceleration, which depends on the assumed rheology.

In contrast to the first two terms, which are parallel to the bed, the third term (cn) is normal to the bed. It describes the centripetal acceleration required to keep the velocity parallel to the bed. This term is not immediately needed in models with a coordinate system aligned to the bed, in which the two-component velocity vector is parallel to the bed by definition. In such models, the centripetal acceleration is only used for computing the dynamic contribution to Coulomb friction via the normal stress at the bed. In some models, such as the original version of RAMMS, this contribution is even neglected. Other approaches (e.g., Fischer et al., 2012) compute c from the local curvature. However, the concept proposed in this paper uses a simpler approach, which will be presented in Sect. 3.2.

2.2 The gravitational acceleration

Since the expression for a is typically derived in a bed-parallel coordinate system, its computation in Cartesian coordinates is briefly recapitulated in the following. The main idea stems from the Navier–Stokes equations for an inviscid fluid, where

$$\begin{array}{}\text{(8)}& \mathit{\rho }\mathit{a}=-\mathrm{\nabla }p+\mathit{\rho }\left(\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ -g\end{array}\right)\end{array}$$

with the fluid pressure p.

The central approximation of the original shallow-water equations is that a is horizontal and constant along vertical columns, which results in hydrostatic vertical pressure profiles. Similarly, the Savage–Hutter model assumes that a is parallel to the bed and constant along columns normal to the bed. The same is assumed for the Cartesian version proposed here except that the columns are vertical instead of normal to the bed (Fig. 1). So the first condition to be met is $\mathit{a}\cdot \mathit{n}=\mathrm{0}$ with the normal vector n defined in Eq. (5), which implies

$$\begin{array}{}\text{(9)}& \mathrm{\nabla }p\cdot \mathit{n}=-\mathit{\rho }g\cos \mathit{\beta }.\end{array}$$

As a second condition, p=0 at the free surface $s=b+h$. Then, ∇p must be normal to the surface, and thus

$$\begin{array}{}\text{(10)}& \mathrm{\nabla }p=\mathit{\lambda }\left(\begin{array}{c}-\mathrm{\nabla }s\\ \mathrm{1}\end{array}\right)\end{array}$$

with an unknown factor λ. Inserting this relation into Eq. (9) yields

$$\begin{array}{}\text{(11)}& \mathit{\lambda }=-{\displaystyle \frac{\mathit{\rho }g}{\mathrm{1}+\mathrm{\nabla }s\cdot \mathrm{\nabla }b}}.\end{array}$$

Since a and thus also ∇p are vertically constant, p increases linearly downward from the surface (p=0). Then the pressure at the bed is

$$\begin{array}{}\text{(12)}& p_{\mathrm{b}}=\mathrm{\nabla }p\cdot \left(\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ -h\end{array}\right)={\displaystyle \frac{\mathit{\rho }gh}{\mathrm{1}+\mathrm{\nabla }s\cdot \mathrm{\nabla }b}}\end{array}$$

and the bed-parallel acceleration

$$\begin{array}{}\text{(13)}& \begin{array}{rl}\mathit{a}& ={\displaystyle \frac{g}{\mathrm{1}+\mathrm{\nabla }s\cdot \mathrm{\nabla }b}}\left(\begin{array}{c}-\mathrm{\nabla }s\\ \mathrm{1}\end{array}\right)+\left(\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ -g\end{array}\right)\\ & =-{\displaystyle \frac{g}{\mathrm{1}+\mathrm{\nabla }s\cdot \mathrm{\nabla }b}}\left(\begin{array}{c}\mathrm{\nabla }s\\ \mathrm{\nabla }s\cdot \mathrm{\nabla }b\end{array}\right)=-{\displaystyle \frac{p_{\mathrm{b}}}{\mathit{\rho }h}}\left(\begin{array}{c}\mathrm{\nabla }s\\ \mathrm{\nabla }s\cdot \mathrm{\nabla }b\end{array}\right).\end{array}\end{array}$$

As a central property, the absolute value of the acceleration is

$$\begin{array}{}\text{(14)}& \left|\mathit{a}\right|=g{\displaystyle \frac{\left|\mathrm{\nabla }b\right|}{\sqrt{\mathrm{1}+|\mathrm{\nabla }b{|}^{\mathrm{2}}}}}=g\sin \mathit{\beta }\end{array}$$

if the free surface is parallel to the bed (∇s=∇b) and zero if the free surface is horizontal (∇s=0).

Figure 4 shows the one-dimensional version of Eq. (13) for slope angles $\mathit{\beta }=\mathrm{15}{}^{\circ }$ and $\mathit{\beta }=\mathrm{45}{}^{\circ }$, where the angle ϕ describes the slope of the free surface (Fig. 1). As a striking property, a singularity occurs at $\mathit{\varphi }=-\mathrm{75}{}^{\circ }$ for $\mathit{\beta }=\mathrm{15}{}^{\circ }$ and at $\mathit{\varphi }=-\mathrm{45}{}^{\circ }$ for $\mathit{\beta }=\mathrm{45}{}^{\circ }$. It occurs if the denominator in Eq. (12) approaches zero, which is the case if the free surface is normal to the bed ($\mathit{\varphi }=\mathit{\beta }-\mathrm{90}{}^{\circ }$). The pressure p_b grows towards infinity then and even becomes negative after passing the singularity, which is unrealistic. The respective range corresponds to the dashed lines in Fig. 4, which show a positive (downslope) acceleration for an uphill-facing front ($\mathit{\varphi }<\mathrm{0}{}^{\circ }$).

Figure 4 Acceleration by gravity with the original pressure at the bed (Eq. 13) and with the modified pressure (Eq. 17). The angle ϕ describes the slope of the free surface s (Fig. 1). Negative values of ϕ correspond to an inclination opposite to the bed. The bed-parallel acceleration a is $\pm \left|\mathit{a}\right|$ with a positive sign downslope and a negative sign upslope. The dashed lines illustrate the unrealistic range for the pressure derived from Eq. (13) and the dash-dotted lines the respective parts of the curves in a coordinate system aligned to the bed.

One might argue that this situation is far outside the scope of the theory proposed by Savage and Hutter (1989) and that the occurrence of the singularity is irrelevant. Concerning numerical simulations, however, it is a big advantage if the solution is still well-defined beyond the range of the approximations made. At this point, the widely used formulation in a local coordinate system is better than the Cartesian version. If the thickness h is measured perpendicular to the bed, there is no singularity but just $\left|\mathit{a}\right|\to \mathrm{\infty }$ for $\left|\mathrm{\nabla }h\right|\to \mathrm{\infty }$. The allowed range for ϕ in Fig. 4 would be $\mathit{\varphi }\in (\mathit{\beta }-\mathrm{90}{}^{\circ },\mathit{\beta }+\mathrm{90}{}^{\circ })$ instead of $(-\mathrm{90},\mathrm{90}{}^{\circ })$, which means that the dashed line segments to the left of the singularity exist no longer. In turn, the allowed range of ϕ includes the dash-dotted line segments.

In turn, passing the singularity in the Cartesian version causes an unrealistic behavior. If an uphill-facing front ($\mathit{\varphi }<\mathrm{0}{}^{\circ }$) becomes steeper, there is an increasing outward acceleration (a<0) at first. At a certain steepness, however, the acceleration changes its direction (a>0), causing the material to pile up rapidly. So passing the singularity needs to be inhibited technically, e.g., by imposing a positive lower limit d_min to the denominator in Eq. (12):

$$\begin{array}{}\text{(15)}& p_{\mathrm{b}}={\displaystyle \frac{\mathit{\rho }gh}{max(\mathrm{1}+\mathrm{\nabla }s\cdot \mathrm{\nabla }b,d_{min})}}.\end{array}$$

However, the maximum acceleration at a downhill-facing front would be limited even for $\mathit{\varphi }\to \mathrm{90}{}^{\circ }$. So there is an asymmetry in the acceleration in the form that uphill-facing fronts will typically cause a higher acceleration than downhill-facing fronts. As a consequence, a small-scale roughness of the free surface will cause an uphill acceleration in total. This issue will be investigated numerically in Sect. 5.1.

In order to overcome this strong limitation, the model MinVoellmy uses a simplified expression for the pressure at the bed, which assumes that the free surface is parallel to the bed. So ∇s is replaced by ∇b in Eq. (12), and thus

$$\begin{array}{}\text{(16)}& p_{\mathrm{b}}={\displaystyle \frac{\mathit{\rho }gh}{\mathrm{1}+|\mathrm{\nabla }b{|}^{\mathrm{2}}}}=\mathit{\rho }gh{\cos }^{\mathrm{2}}\mathit{\beta }.\end{array}$$

This approximation can be interpreted as the hydrostatic pressure caused by the normal component of gravity, gcos β, for a layer of a thickness of hcos β (perpendicular to the bed). The acceleration (Eq. 13) then simplifies to

$$\begin{array}{}\text{(17)}& \mathit{a}=-g{\cos }^{\mathrm{2}}\mathit{\beta }\left(\begin{array}{c}\mathrm{\nabla }s\\ \mathrm{\nabla }s\cdot \mathrm{\nabla }b\end{array}\right).\end{array}$$

This modification transfers the good properties of the original approach in a slope-aligned coordinate system to a Cartesian coordinate system. In particular, a is linear in ∇s at constant ∇b. This linearity ensures that a small-scale roughness of the free surface causes no acceleration in total and that a horizontal surface (∇s=0) is stable.

In turn, however, the curves of the two expressions for the acceleration (Eqs. 13 and 17) are not tangential to each other at ϕ=β but cross each other with different slopes. Thus, Eq. (17) is not a first-order approximation to Eq. (13) concerning the difference ϕ−β for ϕ≈β. The acceleration is predicted correctly if the free surface is parallel to the bed, but the effect of an inclination of the free surface relative to the bed is overestimated. It will be shown that this overestimation has a minor effect on avalanche fronts even on steep slopes. It may, however, have a stronger effect on the tails of avalanches and on the propagation of waves at the free surface. There might be options to achieve a better approximation for ϕ≈β preserving the good properties of the equation. For the sake of simplicity and parsimony, however, the simple approximation is used in the following.

2.3 Friction

The friction term in Eq. (7) refers to the acceleration of a vertical column. Here it has to be taken into account that the shear stress τ at the bed acts on the inclined bed area, which is locally by a factor of $\frac{\mathrm{1}}{\cos \mathit{\beta }}$ greater than the horizontal area. Then the deceleration is

$$\begin{array}{}\text{(18)}& f={\displaystyle \frac{\mathit{\tau }}{\mathit{\rho }h\cos \mathit{\beta }}}\end{array}$$

with the shear stress τ from Eq. (2). For v≥v_c, we obtain

$$\begin{array}{}\text{(19)}& f={\displaystyle \frac{g}{\mathit{\xi }h\cos \mathit{\beta }}}|\mathit{v}{|}^{\mathrm{2}}\end{array}$$

without further assumptions. In the range of Coulomb friction (v<v_c), however, the normal stress at the bed must be specified. For simplicity, effects of internal deformation of the granular medium (typically expressed in terms of the so-called earth pressure coefficients) are neglected in the following. If the friction arising from the centripetal acceleration is also neglected, the normal stress σ is given by the pressure at the bed p_b, which leads to

$$\begin{array}{}\text{(20)}& f=\mathit{\mu }{\displaystyle \frac{p_{\mathrm{b}}}{\mathit{\rho }h\cos \mathit{\beta }}}.\end{array}$$

The increase in friction due to the centripetal acceleration c (Eq. 7) can easily be included in the form

$$\begin{array}{}\text{(21)}& f=\mathit{\mu }\left({\displaystyle \frac{p_{\mathrm{b}}}{\mathit{\rho }h\cos \mathit{\beta }}}+c\right)=\mathit{\mu }\left(g\cos \mathit{\beta }+c\right)\end{array}$$

in combination with the modified pressure at the bed (Eq. 16). On concave profiles, however, f may even become negative. In this situation, the mobile material would detach from the bed, which is not captured by models of this type. Since an acceleration by friction is unrealistic, negative values of f should be replaced by 0, and thus

$$\begin{array}{}\text{(22)}& f=\left\{\begin{array}{ccc}\mathit{\mu }max\left(g\cos \mathit{\beta }+c,\mathrm{0}\right)& & v<v_{\mathrm{c}}\\ \frac{g}{\mathit{\xi }h\cos \mathit{\beta }}|\mathit{v}{|}^{\mathrm{2}}& \text{for}& v\ge v_{\mathrm{c}}\end{array}\right..\end{array}$$

3 Numerical implementation

The numerical implementation uses a regular grid with a constant spacing δx in the x direction and δy in the y direction. The variables are the thickness h and the depth-integrated momentum per unit mass hv as a three-component vector. Both variables are considered at the nodes (x,y).

3.1 Mass balance

The mass balance (Eq. 6) is an advection equation without source terms. An explicit Euler scheme is used in combination with an upstream discretization. The question of why this simple scheme is sufficient for the Voellmy rheology considered here will be addressed in Sect. 5.2. Applying the explicit Euler scheme to Eq. (6) yields

$$\begin{array}{}\text{(23)}& h(t+\mathit{\delta }t)=h\left(t\right)-\mathit{\delta }t\left({\displaystyle \frac{\partial \left(v_{x}h\right(t\left)\right)}{\partial x}}+{\displaystyle \frac{\partial \left(v_{y}h\right(t\left)\right)}{\partial y}}\right).\end{array}$$

First, the values v_x and v_y are computed at the nodes from hv and h. Then the values of v_x are interpolated linearly to the points $(x\pm \frac{\mathit{\delta }x}{\mathrm{2}},y)$. Depending on the sign of this velocity, h at either of the nodes (the upstream point) is adopted to $(x\pm \frac{\mathit{\delta }x}{\mathrm{2}},y)$, and a central difference quotient is used for the x derivative at (x,y). The same procedure is applied to the y derivative.

3.2 Momentum balance

The momentum balance (Eq. 7) is separated into several steps.

Step 1: advection

The first step is basically the same as for h, except that the result is not (hv)(t+δt) but an intermediate value

$$\begin{array}{}\text{(24)}& (h\mathit{v}{)}^{\prime }=(h\mathit{v}\left)\right(t)-\mathit{\delta }t\left({\displaystyle \frac{\partial \left(v_x\right(h\mathit{v}\left)\right(t\left)\right)}{\partial x}}+{\displaystyle \frac{\partial \left(v_y\right(h\mathit{v}\left)\right(t\left)\right)}{\partial y}}\right).\end{array}$$

Step 2: centripetal acceleration

The second step addresses the last term at the right-hand side of Eq. (7) and computes a second intermediate value

$$\begin{array}{}\text{(25)}& (h\mathit{v}{)}^{\prime \prime }=(h\mathit{v}{)}^{\prime }+\mathit{\delta }thc\mathit{n}.\end{array}$$

The centripetal acceleration c is obtained from the condition that v must be parallel to the bed ($(h\mathit{v}{)}^{\prime \prime }\cdot \mathit{n}=\mathrm{0}$), which yields

$$\begin{array}{}\text{(26)}& \mathit{\delta }thc=-(h\mathit{v}{)}^{\prime }\cdot \mathit{n}.\end{array}$$

Since the centripetal acceleration should not change the absolute value of the velocity, (hv)^′′ is then rescaled in such a way that $\left|\right(h\mathit{v}{)}^{\prime \prime }|=|\left(h\mathit{v}{)}^{\prime }\right|$.

Step 3: gravitational acceleration

This step generates the next intermediate values

$$\begin{array}{}\text{(27)}& (h\mathit{v}{)}^{\prime \prime \prime }=(h\mathit{v}{)}^{\prime \prime }+\mathit{\delta }th\mathit{a}.\end{array}$$

The gravitational acceleration a parallel to the bed (Eq. 13 or 17) requires the gradient of the free surface $s=b+h$. In order to avoid checkerboard problems, one-sided difference quotients are used, although central difference quotients would yield a better accuracy theoretically. If central difference quotients were used, a checkerboard pattern with one value of s at the black fields of a checkerboard and another value at the red fields would yield no acceleration and thus be stable. One-sided difference quotients in the direction of steepest increase (or least steep descent) in s turned out to be most robust.

Step 4: friction

The final step has the form

$$\begin{array}{}\text{(28)}& \left(h\mathit{v}\right)(t+\mathit{\delta }t)=(h\mathit{v}{)}^{\prime \prime \prime }-\mathit{\delta }thf{\displaystyle \frac{\mathit{v}}{\left|\mathit{v}\right|}}.\end{array}$$

Since friction is opposite to the velocity, (hv)(t+δt) is obtained by rescaling (hv)^′′′.

For v<v_c, the amount in hv consumed by friction is $\mathit{\delta }th\mathit{\mu }max\left(g\cos \mathit{\beta }+c,\mathrm{0}\right)$ according to Eq. (22). If this amount exceeds the actual amount $\left|\right(h\mathit{v}{)}^{\prime \prime \prime }|$, it is assumed that the movement stops. Accordingly, the length of the vector (hv)(t+δt) is

$$\begin{array}{}\text{(29)}& \left|\right(h\mathit{v}\left)\right(t+\mathit{\delta }t\left)\right|=max\left(\left|\right(h\mathit{v}{)}^{\prime \prime \prime }|-\mathit{\delta }th\mathit{\mu }max\left(g\cos \mathit{\beta }+c,\mathrm{0}\right),\mathrm{0}\right),\end{array}$$

where the outer maximum function takes the stopping criterion into account.

For v≥v_c, we have to take into account that friction depends on velocity. Since f→∞ for h→0 according to Eq. (22), an implicit scheme is used here in order to avoid an additional limitation to δt. This means that v in Eq. (22) is expressed in terms of (hv)(t+δt), which yields

$$\begin{array}{}\text{(30)}& \left|\right(h\mathit{v}\left)\right(t+\mathit{\delta }t\left)\right|=\left|\right(h\mathit{v}{)}^{\prime \prime \prime }|-{\displaystyle \frac{g\mathit{\delta }t}{\mathit{\xi }{h}^{\mathrm{2}}\cos \mathit{\beta }}}|\left(h\mathit{v}\right)(t+\mathit{\delta }t){|}^{\mathrm{2}}.\end{array}$$

This is a quadratic equation in $\left|\right(h\mathit{v}\left)\right(t+\mathit{\delta }t\left)\right|$. It is solved by

$$\begin{array}{}\text{(31)}& \left|\right(h\mathit{v}\left)\right(t+\mathit{\delta }t\left)\right|=\sqrt{{\mathit{\gamma }}^{\mathrm{2}}+\mathrm{2}\mathit{\gamma }\left|\right(h\mathit{v}{)}^{\prime \prime \prime }|}-\mathit{\gamma }\end{array}$$

with

$$\begin{array}{}\text{(32)}& \mathit{\gamma }={\displaystyle \frac{\mathit{\xi }{h}^{\mathrm{2}}\cos \mathit{\beta }}{\mathrm{2}g\mathit{\delta }t}}.\end{array}$$

4 Software description

At present, MATLAB and Python implementations of MinVoellmy are available under the GNU General Public License. None of them requires specific packages, except for NumPy for the Python version. Each version contains separate classes for the one- and two-dimensional versions. The implementation is minimalistic. The classes contain a constructor and a method step for performing a forward time step but neither methods for input/output nor graphics components.

The constructor requires six (one-dimensional) or seven (two-dimensional) mandatory arguments:

• Two arrays are b and h for the bedrock elevation b and the initial thickness h.

• The grid spacing is dx and dy (only for the two-dimensional version).

• The physical parameters are μ (mu), ξ (xi), and v_c (vc). The latter describes the crossover velocity at a thickness h=1 m, and the actual value of v_c is computed from Eq. (3). Values v_c≤0 switch to the conventional Voellmy rheology. These three parameters may be either scalar values or arrays.

Further optional arguments are as follows:

• A minimum thickness is defined as h_min (hmin). It is assumed that material can move only if $h>h_{min}$ (default is 0). In order to increase efficiency, the two-dimensional version restricts the computation to a rectangle around the active region in each time step.

• The minimum value is d_min (dmin) for the denominator in Eq. (15) in case the original expression for the pressure is used. The simplified expression (Eq. 16) is used for $d_{min}=\mathrm{0}$, which is strongly recommended (default is 0).

• A logical value cent is used to define whether the effect of the centripetal acceleration on Coulomb friction is taken into account (default is true).

• The gravitational acceleration is g (g) (default is 9.81).

The method step for the forward time step updates the thickness h and the Cartesian components uh, vh (only in the two-dimensional version) and wh of the momentum vector hv. The time increment δt (dt) is the only mandatory parameter. Optionally, an upper limit cfl for the Courant number,

$$\begin{array}{}\text{(33)}& C={\displaystyle \frac{\left|v_x\right|\mathit{\delta }t}{\mathit{\delta }x}}+{\displaystyle \frac{\left|v_y\right|\mathit{\delta }t}{\mathit{\delta }y}},\end{array}$$

can be defined. In this case, δt is reduced automatically if C> cfl, and the reduced value of δt is returned. This option can be used for adjusting δt dynamically to the velocity. According to the Courant–Friedrichs–Lewy (CFL) criterion, C>1 makes the explicit scheme unstable. Setting cfl to a sufficiently small value, e.g., 0.5, avoids instability of the advection terms. However, the CFL criterion does not capture the acceleration and friction terms, so δt must not be too large even when using the optional argument cfl. In particular, the transition in friction at v_c may require quite small time increments δt.

5 Numerical tests

In this section, several tests addressing the fundamental properties of the modified rheology and the numerical approach are presented. Unless stated explicitly, the parameter values μ=0.75, ξ=250 m s⁻², and v_c=5 m s⁻¹ (at a thickness of h=1 m) are used.

5.1 The modified pressure

The first numerical test refers to the modification applied to the pressure at the bed in Sect. 2.2. The scenario is a pile of 100 m width and 50 m height placed on a slope with $\mathit{\beta }=\mathrm{30}{}^{\circ }$. In order to point out the differences between the approximations for the pressure (Eqs. 15 and 16) more clearly, it is assumed that the material is fluidized from the beginning (μ=0, v_c=0). The simulations are performed on a rather fine grid with δx=1 m and a time increment $\mathit{\delta }t={\mathrm{10}}^{-\mathrm{3}}$ s.

A simulation with a coordinate system aligned to the bed is used as a reference. Technically, g in the acceleration term has to be multiplied by cos β here, and ∇b has to be set to zero. In turn, an additional downslope acceleration gsin β has to be taken into account.

As shown in Fig. 5, the different simulations start similarly with some wave-like structures downslope of the crest, owing to the opposite directions of the acceleration at the crest. For the scenarios with the original pressure and small lower limits d_min of the denominator in Eq. (15), these waves propagate uphill rapidly. The pile even moves uphill in total for $d_{min}\le \mathrm{0.1}$, which is completely unrealistic. This artifact was already discussed in Sect. 2.2.

Figure 5 Collapse of a granular pile on a slope with $\mathit{\beta }=\mathrm{30}{}^{\circ }$. The colored curves refer to the original pressure with different lower limits (Eq. 15) and the black line to the modified pressure (Eq. 16). The gray shaded area shows the result of a simulation with a coordinate system aligned to the slope as a reference. Curves are shifted vertically, and labels at the z axis refer to the first and last plot.

In turn, the modified pressure (Eq. 16) overestimates the acceleration at the steep downslope front. As a consequence, spreading of the downslope front is stronger than in the reference scenario. The overall shape is, however, similar for t≥4 s. Due to the faster spreading, the layer is slightly thinner with the modified pressure, resulting in a lower velocity. So the lead over the reference scenario decreases through time.

The scenario with a strong limitation of the denominator ($d_{min}=\mathrm{1}$) also stays close to the reference scenario for some time. However, the tail persists, while the front is too slow. These results confirm the theoretical arguments given in Sect. 2.2 that the asymmetry in the acceleration term causes artifacts. The version with the original pressure appears to be unsuitable even with the strong limitation of the denominator. In turn, the version with the modified pressure works well, except for the spreading of the downslope front being too fast in the beginning.

5.2 Numerical diffusion

Numerical diffusion is typically the most challenging problem in computational fluid dynamics. In order to investigate the effect of numerical diffusion, the movement of a body of granular material moving down a straight slope with a given slope angle β is considered. Let us focus on steady-state solutions in the sense that the entire body moves at a constant velocity v₀≥v_c without changing its shape. Then the driving acceleration $\left|\mathit{a}\right|$ must be balanced by friction at each point, so

$$\begin{array}{}\text{(34)}& \left|\mathit{a}\right|=f={\displaystyle \frac{g}{\mathit{\xi }h\cos \mathit{\beta }}}{v}_{\mathrm{0}}^{\mathrm{2}}\end{array}$$

according to Eq. (19). The simplest solution is an infinite layer with a constant thickness h₀. Then the surface is parallel to the bed, and the acceleration is $\left|\mathit{a}\right|=g\sin \mathit{\beta }$ (Eq. 14) for both pressure models, so Eq. (34) yields

$$\begin{array}{}\text{(35)}& h_{\mathrm{0}}={\displaystyle \frac{v_{\mathrm{0}}^{\mathrm{2}}}{\mathit{\xi }\cos \mathit{\beta }\sin \mathit{\beta }}}.\end{array}$$

Although this solution does not describe a front moving downslope, it helps to write Eq. (34) in a more convenient form:

$$\begin{array}{}\text{(36)}& \left|\mathit{a}\right|=g\sin \mathit{\beta }{\displaystyle \frac{h_{\mathrm{0}}}{h}}.\end{array}$$

In combination with the modified pressure (Eq. 17), the acceleration is

$$\begin{array}{}\text{(37)}& \begin{array}{rl}\left|\mathit{a}\right|& =g{\cos }^{\mathrm{2}}\mathit{\beta }\left(-{\displaystyle \frac{\partial s}{\partial x}}\right)\sqrt{\mathrm{1}+{\left({\displaystyle \frac{\partial b}{\partial x}}\right)}^{\mathrm{2}}}\\ & =g\cos \mathit{\beta }\left(-{\displaystyle \frac{\partial s}{\partial x}}\right)=g\cos \mathit{\beta }\left(\tan \mathit{\beta }-{\displaystyle \frac{\partial h}{\partial x}}\right),\end{array}\end{array}$$

which leads to the differential equation

$$\begin{array}{}\text{(38)}& {\displaystyle \frac{\partial h}{\partial x}}=\tan \mathit{\beta }\left(\mathrm{1}-{\displaystyle \frac{h_{\mathrm{0}}}{h}}\right).\end{array}$$

This equation can be solved analytically in the form x(h) instead of h(x). It is recognized by computing the derivative that the solution is

$$\begin{array}{}\text{(39)}& x={\displaystyle \frac{h+h_{\mathrm{0}}\ln \left(\mathrm{1}-\frac{h}{h_{\mathrm{0}}}\right)}{\tan \mathit{\beta }}}+\text{const.}\end{array}$$

This analytical solution is the reference for the numerical test. It describes a front that becomes vertical ($\frac{\partial h}{\partial x}\to \mathrm{\infty }$) for h→0, while h→h₀ for $x\to -\mathrm{\infty }$. It is plotted as a blue line without markers in Fig. 6 for h₀=10 m.

Figure 6 Front of a layer with a vertical thickness h₀=10 m after traveling a horizontal distance of 10 km. Curves without markers refer to the analytical solutions. Markers refer to the nodes of the numerical solutions with δx=1 m and δx=10 m. All curves are centered horizontally in such a way that the mean thickness over an interval from −250 to 250 m is $\frac{h_{\mathrm{0}}}{\mathrm{2}}$.

For the numerical representation, it is assumed that h=h₀ and $\left|h\mathit{v}\right|=h_{\mathrm{0}}{v}_{\mathrm{0}}$ at the left-hand boundary. This leads to a front propagating with the velocity v₀ (parallel to the bed). Propagation is simulated over a horizontal distance of 10 km in order to approach the steady-state shape for δx=1 m and δx=10 m. In order to minimize the effect of δt, the small value $\mathit{\delta }t={\mathrm{10}}^{-\mathrm{3}}$ s is still used here.

As expected, the deviation from the analytical solution increases with increasing grid spacing. The artificial widening of the front is, however, limited to a few times δx. More important, this widening is limited in time. An initially sharp front is widened rapidly, but width and shape stabilize soon.

This behavior arises from a specific property of Voellmy-type rheologies. The frictional acceleration (Eq. 19) increases with decreasing thickness. So material running ahead of the front is decelerated and thus overrun by the front. As a consequence, the shape of the front predicted by Eq. (39) is stable. The simple upstream scheme introduces numerical diffusion and widens the front, but the stability of the front counteracts this process. This is the reason why numerical diffusion is not a serious problem in combination with Voellmy-type rheologies, in contrast to many other applications of the shallow-water equations.

For completeness, Fig. 6 also shows the results obtained for the original pressure (Eq. 13) and for the model with a coordinate system parallel to the bed. For the latter, the results are also transferred to the Cartesian coordinate system with h measured vertically. The analytical solution is the same for both versions. Using Eq. (13) instead of Eq. (17) for the acceleration yields the differential equation

$$\begin{array}{}\text{(40)}& {\displaystyle \frac{\partial h}{\partial x}}={\displaystyle \frac{\tan \mathit{\beta }\left(\mathrm{1}-\frac{h_{\mathrm{0}}}{h}\right)}{\mathrm{1}-{\sin }^{\mathrm{2}}\mathit{\beta }\frac{h_{\mathrm{0}}}{h}}},\end{array}$$

which can be solved in the form

$$\begin{array}{}\text{(41)}& x={\displaystyle \frac{h+h_{\mathrm{0}}{\cos }^{\mathrm{2}}\mathit{\beta }\ln \left(\mathrm{1}-\frac{h}{h_{\mathrm{0}}}\right)}{\tan \mathit{\beta }}}+\text{const.}\end{array}$$

While the solution is similar to Eq. (39), the respective front is steeper than with the modified pressure and is even overhanging in the Cartesian coordinate system used here. The numerical solution with the coordinate system aligned to the bed approximates the front better than in Cartesian coordinates, in particular if the slope is steep. This result is not surprising since the front is normal to the bed at the bottom, while the Cartesian version is limited by a vertical front.

Overall, however, the differences between all versions are rather small. In combination with the findings of the previous section, this finding justifies the major approximations introduced in the model MinVoellmy. First, the modified pressure provides a reasonable approximation and makes a treatment in Cartesian coordinates with the thickness measured vertically feasible. Second, the simple upstream scheme for the advection terms works reasonably well in combination with the Voellmy rheology.

5.3 Comparison to the conventional Voellmy rheology

Figure 7 shows the results for a slope with $\mathit{\beta }=\mathrm{45}{}^{\circ }$ combined with a horizontal plane. The source area is a segment of an ellipse with an aspect ratio of 4:1 and a vertical wall at the upper edge at a height of 1900 m above the runout plane. Grid spacing is δx=10 m and time increment δt=0.01 s. In combination with the modified pressure, which is also used in the following examples, the Cartesian approach remains stable, although the applicability of the Savage–Hutter theory to the thick detached body is not given here, regardless of the rheology and the numerical scheme.

Figure 7 Snapshots of a one-dimensional simulation for (a) the modified rheology with μ=0.75 and (b) the conventional Voellmy rheology with Coulomb friction reduced by a factor of  10 (μ=0.075). Orange areas indicate v<v_c for the modified rheology (a) and that the Coulomb friction term is greater than the v² friction term for the conventional rheology (b). Blue areas correspond to the opposite situation, and black areas are at rest. The triangles depict the center of mass and the front. Curves are shifted vertically, and labels at the z axis refer to the first and last plot.

For the scenario with the conventional Voellmy rheology, the Coulomb friction term has to be reduced by a factor of 10 (so μ=0.075) in order to achieve a similar runout length. The two versions differ strongly already during the phase of mobilization. Owing to the artificially reduced coefficient of friction μ, the body is mobilized much faster in the conventional scenario than with the modified rheology. In order to avoid the fast mobilization, Aaron and Hungr (2016) extended the model DAN3D by considering a more or less rigid block during the first phase of movement. This might not be necessary for the modified rheology. However, it should be kept in mind that we are outside the range of applicability of the Savage–Hutter theory here. So the mobilization seems to be more realistic with the modified rheology but is not necessarily physically correct.

The thickness overshoots at the lower edge of the detachment area. This overshooting is inevitable at a sharp kink in topography. Due to the finite acceleration, the absolute value of the velocity changes only gradually across the kink. Thus, the abrupt change in flow direction at the kink causes an abrupt drop in the horizontal component of the velocity. Conservation of mass requires an increase in vertical thickness h, which causes the overshooting. This effect is, however, neither unique to the approach used here nor a serious problem.

Further differences between the two versions occur during the runout in the horizontal plane. For the modified rheology, friction increases instantaneously when the velocity drops below v_c, and movement stops soon. As a consequence, a large part of the mass comes to rest within a narrow time span between t=80 s and t=100 s. This region expands towards both sides and is bounded by two narrow ranges with Coulomb friction. The material upstream of this region piles up, which results in a hill in the deposits.

In turn, there is a rather thin layer propagating with little friction, which finally dominates the runout length. Despite the difference in μ, the maximum runout length is even bigger than for the scenario with the conventional Voellmy rheology, while the distance traveled by the center of mass is similar. This result confirms the findings of Hergarten (2024) that the modified Voellmy rheology allows for a long runout without assuming an artificially low coefficient of friction μ.

5.4 Two-dimensional simulations

This section mainly addresses the effect of the orientation of the grid in two dimensions. The starting point is the situation considered in the previous section but extended into the y direction by an ellipsoidal volume with an aspect ratio of $\mathrm{4}:\mathrm{1}:\mathrm{1}$. For this scenario, the total detached volume is V=0.1 km³. In addition, a smaller volume of V=0.02 km³ and a larger volume of V=0.5 km³ are considered, all with the same aspect ratio and the same height of 1900 m above the horizontal plane. Grid spacing is $\mathit{\delta }x=\mathit{\delta }y=\mathrm{10}$ m. The time increment is δt=0.1 s with an additional upper limit C≤0.25 for the Courant number (Eq. 33). All simulations were run over a total time span of 500 s.

Figure 8 shows the final deposits of the three volumes for a slope aligned to the coordinate axes and for a slope aligned to the diagonal line in the x–y plane. The small-scale topography of the deposits depends on the orientation of the slope. In particular, striations occur if the velocity is parallel to one of the coordinate axes. Since the principal flow direction is parallel to the x axis for the axis-parallel scenario, this effect is stronger here than for the diagonal setup.

Figure 8 Deposits obtained from volumes of V=0.5 km³, V=0.1 km³, and V=0.02 km³ with axis-parallel and diagonal orientations. Only deposits with a thickness h≥0.1 m are shown. The red lines are the outlines of the detachment area at the slope, which is almost circular in the x–y plane. The close-ups illustrate the partly hummocky topography and striations.

As already recognized in Sect. 5.3, the formation of hummocky deposit morphologies arises from the discontinuity in friction at v_c. This discontinuity also allows for the formation of striations since the governing equations do not include transverse diffusion of momentum. Owing to the advective characteristics of Eq. (7), different flow lines are in principle decoupled. So the velocity may differ among parallel flow lines without any effect. However, transverse numerical diffusion results in an exchange of momentum between different flow lines if the velocity is not parallel to any of the coordinate axes. This effect is apparently strong enough to suppress the formation of striations.

The occurrence of longitudinal striations is not unrealistic (e.g., Shreve, 1966; Pietrek et al., 2020), and the modified Voellmy rheology may open a door towards understanding their origin. However, proceeding in this direction requires an extension of the model by transverse diffusion of momentum due to particle collisions in order to find out whether the tendency towards forming striations is strong enough to overcome the diffusion of momentum.

The overall shape of the deposits is, however, not affected strongly by the orientation of the grid. This visual impression from Fig. 8 is confirmed by the profiles plotted in Fig. 9. Overall, these results suggest that there is no need to align the coordinate system to the principal orientation of the slope when simulating real-world scenarios since the local topography will override effects of the orientation.

Figure 9 Profiles across the deposits shown in Fig. 8. The profiles were obtained by integrating the height above the bed along the direction perpendicular to the respective profile, which can also be interpreted as volume per profile length. Solid lines refer to the axis-parallel alignment of the slope and dashed lines to the diagonal alignment.

5.5 Fahrboeschung ratios

Explaining the long runout of large rock avalanches was the main goal of developing the modified Voellmy rheology. The relation between runout and volume is often expressed in terms of the Fahrboeschung ratio $\frac{H}{L}$, also called Heim's ratio, where H is the total height drop and L the maximum horizontal runout length. Both properties are measured from the upper edge of the detached volume. Scheidegger (1973) found a power-law relation

$$\begin{array}{}\text{(42)}& {\displaystyle \frac{H}{L}}\propto V^{-\mathit{\alpha }}\end{array}$$

with α=0.16, which was confirmed later by Legros (2002).

The dashed lines in Fig. 10 show the relation obtained numerically with the volumes from the previous section extended to V=0.001, 0.002, 0.005, 0.01, …, 0.5 km³. While the power-law relation between volume and Heim's ratio is reproduced qualitatively well for V≥0.01 km³, the effect of volume on Heim's ratio is overestimated. Fitting power laws yields α=0.25 for the two-dimensional version and even α=0.31 for the one-dimensional version.

Figure 10 Heim's ratio as a function of the volume. Dashed lines refer to a total height of H=1900 m as considered in Sect. 5.3 and 5.4. Solid lines were obtained taking into account the correlation between H and V observed by Legros (2002). Black lines represent power laws fitted for V≥0.01 km³.

This finding is in line with the results of Hergarten (2024), who considered the same rheology in a lumped-mass model and found a strong influence of the absolute height drop H. While larger volumes tend to have a larger height drop in nature, the resulting increase in L is weaker than the increase in H. So Heim's ratio increases with increasing H. In order to take this effect into account, the solid lines in Fig. 10 describe the same scenarios as before but with H∝V^0.09 as found empirically by Legros (2002). The factor of proportionality was chosen to keep H=1900 m from the previous simulations for V=0.1 km³, which yields a range from H=1255 m for V=0.001 km³ to H=2338 m for V=1 km³. In contrast to the previous simulations, V=1 km³ is now possible without the detachment area reaching the foot of the slope.

Including the correlation of H and V reduces the exponent α considerably. For the two-dimensional version, α=0.20 is obtained, which is still higher than observed in nature. The residual deviation can at least partly be attributed to the spatial scaling assumed here. For simplicity, the detached volume was assumed to scale isotropically (so with V^⅓ in each direction). As explained by Hergarten (2024), however, the thickness should be the primary control on runout rather than V itself. Larsen et al. (2010) found the relation V∝A^1.40 between volume and area, which means that the increase in thickness is weaker than V^⅓ in reality. Since there are several dependencies beyond this scaling relation, the value α=0.20 seems to be good enough for a first test. Overall, these findings confirm the results obtained by Hergarten (2024) for a lumped mass and suggest that the two-dimensional scenario without lateral confinement yields similar scaling properties as the simple one-dimensional lumped-mass model.

6 Conclusions

In this study, a simple and lightweight numerical implementation of the modified Voellmy rheology proposed by Hergarten (2024) was presented. The simplicity of the implementation arises from two features. First, a fully Cartesian description is used, where the thickness of the mobile layer is measured vertically instead of perpendicularly to the bed. This concept harmonizes well with a simplified expression for the pressure at the bed. As a second feature, a simple upstream scheme is used for the advection terms. While upstream schemes typically suffer from numerical diffusion, it was shown that numerical diffusion is not a serious problem in combination with Voellmy-type rheologies.

The results obtained from the numerical tests confirm the findings of Hergarten (2024) for the simple lumped-mass model. Furthermore, the modified Voellmy rheology may open doors towards understanding hummocky deposit morphologies and longitudinal striations. In view of the simplicity of the numerical implementation, these preliminary results are promising.

In turn, however, the purpose of the recent implementation must be kept in mind. The MinVoellmy software is lightweight but minimalistic. It does not offer any user interface or methods for importing and exporting data, so it cannot be operated without programming some parts in MATLAB or Python, which limits its field to research and teaching. It is also not designed to compete with comprehensive models such as r.avaflow, which offer additional options such as multi-phase flow.

In particular, the recent implementation is not designed for operational hazard assessment. This restriction is owing not only to the short time span of development and limited testing but mainly to the parameter values. The parameters μ and ξ of the modified Voellmy rheology are similar to those of the conventional Voellmy rheology in their meaning, but their numerical values are not the same. So calibrations of other models, such as the extensively tested parameter values from the widely used model RAMMS, cannot be transferred to the modified rheology directly. Furthermore, knowledge about the additional parameter v_c is still limited.