the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
MinVoellmy v1: a lightweight model for simulating rapid mass movements based on a modified Voellmy rheology
Stefan Hergarten
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 velocitydependent 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 twodimensional, depthaveraged 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.
 Article
(4413 KB)  Fulltext XML
 BibTeX
 EndNote
Modeling of rapid mass movements was pushed strongly by the ideas of Savage and Hutter (1989), who extended the shallowwater equations towards granular media. The shallowwater equations provide a twodimensional, depthaveraged description of flow processes with a free surface. In their original form, the shallowwater 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 twodimensional 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
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
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
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 smoothedparticle 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 nonoscillatory central differencing (TVDNOC) 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 nonorthogonal. 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 shallowwater 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 multiphase 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 realworld applications.
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 timedependent 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 depthaveraged velocity parallel to the bed.
Using Cartesian coordinates circumvents several problems arising from nonorthogonal 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
with the normal vector
where ∇b is the twodimensional gradient of the bed and β the slope angle ($\mathrm{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 threecomponent 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.
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 depthaveraged 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:
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 shallowwater equations for an incompressible fluid (e.g., Vreugdenhil, 1994).
The momentum balance can also be written as an advection equation in the form
The advection term at the lefthand side is basically the same as Eq. (6) but for the depthintegrated momentum per unit mass hv instead of the thickness h. It is identical to the respective terms in the original shallowwater equations (e.g., Vreugdenhil, 1994). However, v is not approximately horizontal (as required for the original shallowwater 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 threecomponent vector, and Eq. (7) consists of three scalar equations instead of two equations in the original shallowwater equations and in the Savage–Hutter model.
In contrast to Eq. (6), Eq. (7) has a nonzero source term at the righthand side, which describes the total depthintegrated acceleration. While the lefthand side of Eq. (7) is practically the same in all depthaveraged models, the righthand side is modelspecific. 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.
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 twocomponent 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 bedparallel 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
with the fluid pressure p.
The central approximation of the original shallowwater 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
As a second condition, p=0 at the free surface $s=b+h$. Then, ∇p must be normal to the surface, and thus
with an unknown factor λ. Inserting this relation into Eq. (9) yields
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
and the bedparallel acceleration
As a central property, the absolute value of the acceleration is
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 onedimensional 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 uphillfacing front ($\mathit{\varphi}<\mathrm{0}{}^{\circ}$).
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 welldefined 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 dashdotted line segments.
In turn, passing the singularity in the Cartesian version causes an unrealistic behavior. If an uphillfacing 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):
However, the maximum acceleration at a downhillfacing front would be limited even for $\mathit{\varphi}\to \mathrm{90}{}^{\circ}$. So there is an asymmetry in the acceleration in the form that uphillfacing fronts will typically cause a higher acceleration than downhillfacing fronts. As a consequence, a smallscale 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
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
This modification transfers the good properties of the original approach in a slopealigned coordinate system to a Cartesian coordinate system. In particular, a is linear in ∇s at constant ∇b. This linearity ensures that a smallscale 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 firstorder 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}}{\mathrm{cos}\mathit{\beta}}$ greater than the horizontal area. Then the deceleration is
with the shear stress τ from Eq. (2). For v≥v_{c}, we obtain
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 socalled 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
The increase in friction due to the centripetal acceleration c (Eq. 7) can easily be included in the form
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
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 depthintegrated momentum per unit mass hv as a threecomponent 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
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
Step 2: centripetal acceleration
The second step addresses the last term at the righthand side of Eq. (7) and computes a second intermediate value
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
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
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, onesided 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. Onesided 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
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}t\phantom{\rule{0.125em}{0ex}}h\mathit{\mu}max\left(g\mathrm{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
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
This is a quadratic equation in $\left\right(h\mathit{v}\left)\right(t+\mathit{\delta}t\left)\right$. It is solved by
with
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 twodimensional 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 (onedimensional) or seven (twodimensional) mandatory arguments:

Two arrays are
b
andh
for the bedrock elevation b and the initial thickness h. 
The grid spacing is
dx
anddy
(only for the twodimensional 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 twodimensional 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 twodimensional 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,
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.
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^{−2}, and v_{c}=5 m s^{−1} (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 wavelike 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.
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 steadystate solutions in the sense that the entire body moves at a constant velocity v_{0}≥v_{c} without changing its shape. Then the driving acceleration $\left\mathit{a}\right$ must be balanced by friction at each point, so
according to Eq. (19). The simplest solution is an infinite layer with a constant thickness h_{0}. Then the surface is parallel to the bed, and the acceleration is $\left\mathit{a}\right=g\mathrm{sin}\mathit{\beta}$ (Eq. 14) for both pressure models, so Eq. (34) yields
Although this solution does not describe a front moving downslope, it helps to write Eq. (34) in a more convenient form:
In combination with the modified pressure (Eq. 17), the acceleration is
which leads to the differential equation
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
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_{0} for $x\to \mathrm{\infty}$. It is plotted as a blue line without markers in Fig. 6 for h_{0}=10 m.
For the numerical representation, it is assumed that h=h_{0} and $\lefth\mathit{v}\right={h}_{\mathrm{0}}{v}_{\mathrm{0}}$ at the lefthand boundary. This leads to a front propagating with the velocity v_{0} (parallel to the bed). Propagation is simulated over a horizontal distance of 10 km in order to approach the steadystate 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 Voellmytype 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 Voellmytype rheologies, in contrast to many other applications of the shallowwater 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
which can be solved in the form
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.
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 Twodimensional 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^{3}. In addition, a smaller volume of V=0.02 km^{3} and a larger volume of V=0.5 km^{3} 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 smallscale 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 axisparallel scenario, this effect is stronger here than for the diagonal setup.
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 realworld scenarios since the local topography will override effects of the orientation.
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 powerlaw relation
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^{3}. While the powerlaw relation between volume and Heim's ratio is reproduced qualitatively well for V≥0.01 km^{3}, the effect of volume on Heim's ratio is overestimated. Fitting power laws yields α=0.25 for the twodimensional version and even α=0.31 for the onedimensional version.
This finding is in line with the results of Hergarten (2024), who considered the same rheology in a lumpedmass 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^{3}, which yields a range from H=1255 m for V=0.001 km^{3} to H=2338 m for V=1 km^{3}. In contrast to the previous simulations, V=1 km^{3} 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 twodimensional 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 twodimensional scenario without lateral confinement yields similar scaling properties as the simple onedimensional lumpedmass model.
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 Voellmytype rheologies.
The results obtained from the numerical tests confirm the findings of Hergarten (2024) for the simple lumpedmass 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 multiphase 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.
All codes are available in a Zenodo repository at https://doi.org/10.5281/zenodo.10304665 (Hergarten, 2023b) and can be redistributed under the GNU General Public License. This repository also contains data obtained from the numerical simulations. Interested users are advised to download the most recent version of the MinVoellmy software from http://hergarten.at/minvoellmy (Hergarten, 2023a). No data sets were used in this article.
The author has declared that there are no competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
The author would like to thank Fabian Walter for his very constructive comments and Thomas Poulet for the editorial handling.
This openaccess publication was funded by the University of Freiburg.
This paper was edited by Thomas Poulet and reviewed by Fabian Walter and one anonymous referee.
Aaron, J. and Hungr, O.: Dynamic simulation of the motion of partiallycoherent landslides, Engin. Geol., 205, 1–11, https://doi.org/10.1016/j.enggeo.2016.02.006, 2016. a
Bartelt, P. and Buser, O.: Frictional relaxation in avalanches, Ann. Glaciol., 51, 98–104, https://doi.org/10.3189/172756410791386607, 2010. a
Bouchut, F. and Westdickenberg, M.: Gravity driven shallow water models for arbitrary topography, Commun. Math. Sci., 2, 359–389, https://doi.org/10.4310/CMS.2004.v2.n3.a2, 2004. a, b
Buser, O. and Bartelt, P.: Production and decay of random kinetic energy in granular snow avalanches, J. Glaciol., 55, 3–12, https://doi.org/10.3189/002214309788608859, 2009. a
Christen, M., Kowalski, J., and Bartelt, P.: RAMMS: Numerical simulation of dense snow avalanches in threedimensional terrain, Cold Reg. Sci. Technol., 63, 1–14, https://doi.org/10.1016/j.coldregions.2010.04.005, 2010. a
Denlinger, R. P. and Iverson, R. M.: Granular avalanches across irregular threedimensional terrain: 1. Theory and computation, J. Geophys. Res.Earth, 109, F01014, https://doi.org/10.1029/2003JF000085, 2004. a
Dufresne, A., Bösmeier, A., and Prager, C.: Sedimentology of rock avalanche deposits – Case study and review, Earth Sci. Rev., 163, 234–259, https://doi.org/10.1016/j.earscirev.2016.10.002, 2016. a
Fischer, J.T., Kowalski, J., and Pudasaini, S. P.: Topographic curvature effects in applied avalanche modeling, Cold Reg. Sci. Technol., 74–75, 21–30, https://doi.org/10.1016/j.coldregions.2012.01.005, 2012. a, b
Hergarten, S.: MinVoellmy [software], http://hergarten.at/minvoellmy, last access: 8 December 2023a. a
Hergarten, S.: MinVoellmy v1: a lightweight model for simulating rapid mass movements, Zenodo [code and data set], https://doi.org/10.5281/zenodo.10304665, 2023b. a
Hergarten, S.: Scaling between volume and runout of rock avalanches explained by a modified Voellmy rheology, Earth Surf. Dynam., 12, 219–229, https://doi.org/10.5194/esurf122192024, 2024. a, b, c, d, e, f, g, h, i, j
Hergarten, S. and Robl, J.: Modelling rapid mass movements using the shallow water equations in Cartesian coordinates, Nat. Hazards Earth Syst. Sci., 15, 671–685, https://doi.org/10.5194/nhess156712015, 2015. a, b
Hutter, K., Wang, Y., and Pudasaini, S. P.: The Savage–Hutter avalanche model: how far can it be pushed?, Philos. T. Roy. Soc. A, 363, 1507–1528, https://doi.org/10.1098/rsta.2005.1594, 2005. a
Larsen, I. J., Montgomery, D. R., and Korup, O.: Landslide erosion controlled by hillslope material, Nat. Geosci., 3, 247–251, https://doi.org/10.1038/ngeo776, 2010. a
Legros, F.: The mobility of longrunout landslides, Engin. Geol., 63, 301–331, https://doi.org/10.1016/S00137952(01)000904, 2002. a, b, c
McDougall, S.: A new continuum dynamic model for the analysis of extremely rapid landslide motion across complex 3D terrain, PhD thesis, University of British Columbia, https://doi.org/10.14288/1.0052928, 2006. a
McDougall, S.: 2014 Canadian Geotechnical Colloquium: landslide runout analysis – current practice and challenges, Can. Geotech. J., 54, 605–620, https://doi.org/10.1139/cgj20160104, 2017. a
Mergili, M., Fischer, J.T., Krenn, J., and Pudasaini, S. P.: r.avaflow v1, an advanced opensource computational framework for the propagation and interaction of twophase mass flows, Geosci. Model Dev., 10, 553–569, https://doi.org/10.5194/gmd105532017, 2017. a, b
Nessyahu, H. and Tadmor, E.: Nonoscillatory central differencing for hyperbolic conservation laws, J. Comput. Phys., 87, 408–463, https://doi.org/10.1016/00219991(90)902608, 1990. a
Pietrek, A., Hergarten, S., and Kenkmann, T.: Morphometric characterization of longitudinal striae on Martian landslides and impact ejecta blankets and implications for the formation mechanism, J. Geophys. Res.Planet., 125, e2019JE006255, https://doi.org/10.1029/2019JE006255, 2020. a
Pudasaini, S. P. and Mergili, M.: A multiphase mass flow model, J. Geophys. Res.Earth, 124, 2920–2942, https://doi.org/10.1029/2019JF005204, 2019. a
Rauter, M. and Tuković, Ž.: A finite area scheme for shallow granular flows on threedimensional surfaces, Comput. Fluids, 166, 184–199, https://doi.org/10.1016/j.compfluid.2018.02.017, 2018. a, b
Salm, B.: Flow, flow transition and runout distances of flowing avalanches, Ann. Glaciol., 18, 221–226, https://doi.org/10.3189/S0260305500011551, 1993. a
Savage, S. B. and Hutter, K.: The motion of a finite mass of granular material down a rough incline, J. Fluid Mech., 199, 177–215, https://doi.org/10.1017/S0022112089000340, 1989. a, b, c, d
Scheidegger, A. E.: On the prediction of the reach and velocity of catastrophic landslides, Rock Mech., 5, 231–236, https://doi.org/10.1007/BF01301796, 1973. a
Shreve, R. L.: Sherman Landslide, Alaska, Science, 154, 1639–1642, https://doi.org/10.1126/science.154.3757.1639, 1966. a
Tonnel, M., Wirbel, A., Oesterle, F., and Fischer, J.T.: AvaFrame com1DFA (v1.3): a thicknessintegrated computational avalanche module – theory, numerics, and testing, Geosci. Model Dev., 16, 7013–7035, https://doi.org/10.5194/gmd1670132023, 2023. a
Voellmy, A.: Über die Zerstörungskraft von Lawinen, Schweiz. Bauzeitung, 73, 212–217, https://doi.org/10.5169/seals61891, 1955. a
Vreugdenhil, C. B.: Numerical Methods for ShallowWater Flow, Springer, Berlin, Heidelberg, New York, https://doi.org/10.1007/9789401583541, 1994. a, b