Given an initial release volume, the Voellmy-Salm (VS) avalanche model is solved forward in time for height $H$ and velocity $\bf{U}$. Due to the models strong similarity with mathematical models for shallow water flow, existing solution techniques can be adopted. However, acceleration and friction terms require an appropriate numerical treatment of the source terms.
Due to computational efficiency it is advantageous to restrict the original digital elevation model (DEM) to a smaller region of interest. For each avalanche simulation we specify a flow path surrounding polygon, such that no interaction of the flow with the boundary is expected. We refer to this region as the computational domain $\mathcal D$ of the problem.
The DEM directly induces a two-dimensional surface mesh in three-dimensional space, in which each quadrilateral finite volume is determined by four grid points. The discretized computational domain $\mathcal D$ is surrounded by two rows of ghost cells for the implementation of boundary conditions.
The governing differential equations (the first three equations from the Mathematical Model) compose a system of non-linear hyperbolic equations which can be written concisely as
\(\partial_t{V} + \nabla \cdot F (V ) = S ( V ) \)
\(\begin{eqnarray} & V=\begin{pmatrix} H \\ H U_x \\ H U_y \\ H R \end{pmatrix} \qquad S \left( V \right):=\begin{pmatrix} \dot Q \\ S_{gx} - S_{fx} \\ S_{gy} - S_{fy} \\ \dot P - \dot D \end{pmatrix} &\\ &F \left( V \right)=\begin{pmatrix} H U_x & H U_y \\ H U_x^2 + g_z k_{a/p} \frac{H^2}{2} & H U_x U_y \\ H U_x U_y & H U_y^2 + g_z k_{a/p} \frac{H^2}{2} \\ H R U_x& H R U_y \end{pmatrix}& \end{eqnarray}\)
where $\bf V(x,y,t)$ is the phase vector containing the four unknown state variables and $\bf F \left( \bf V \right)$ is the flux function. The vector $\bf S \left( \bf V \right)$, located on the right-hand side, contains source terms (mass entrainment $\dot Q$, gravitational acceleration $S_{gx}$ and $S_{gy}$, production of random energy $\dot P$) and sink terms (frictional forces $S_{fx}$ and $S_{fy}$ and decay of random kinetic energy $\dot D$).
A cell-centered, positivity conserving HLL Finite Volume scheme is formulated over quadrilateral control volumes, where the dependent variables of the system are stored at the center of the cell. The System of non-linear hyperbolic equations is then integrated over the $i^{th}$ control volume $C_i$ and by applying Gauss' theorem, the volume integral of the flux function $\bf F \left( \bf V \right)$ transforms into a surface integral.
$$
\int_{C_i} \partial_t \bf{V} \, d \bf{x} + \oint_{\partial C_i} \bf{F} \left(\bf{V}\right)\cdot n_{i} \, d \sigma = \int_{C_i} \bf{S}\left(\bf{V}\right) d \bf{x}
$$
where $n_{i}$ denotes the outer normal to control volume $C_i$. The surface integral is discretized based on a HLL flux. To guarantee second order convergence in space and time, a linear reconstruction of the averaged, conserved cell quantities is applied. The time-integration is given by a Runge-Kutta Heun scheme. The cell-averaged quantities of the state variables on a time-discrete level are denoted by $\bf{V^{(n)}_{i}}$, and $\bf{V^{(n)}}$ denotes the whole grid function on time level $t^{(n)}$. Then, concisely, the scheme is written as
$$\bf{V_i}^{(*)} = \bf{V_i}^{(n)} + \frac{\triangle t}{A_{C_i}} \, \triangle \bf{F_i^{(HLL)}}\left(\bf{V}^{(n)}\right)$$
$$\bf{V_i}^{(**)} = \bf{V_i}^{(*)} + \frac{\triangle t}{A_{C_i}} \, \triangle \bf{F_i^{(HLL)}}\left(\bf{V}^{(*)}\right)$$
$$\bf{V_i}^{(n+1)} = \frac{1}{2}\,\left( \bf{V_i}^{(n)} + \bf{V_i}^{(**)}\right)$$
where $\triangle t:=t^{(n+1)}-t^{(n)}$, $A_{C_i}$ denotes the area of $C_i$ and $\triangle \bf{F_i^{(HLL)}}$ stands for an approximation of the surface integral. In our case it is given by a midpoint quadrature rule
$$
\triangle \bf{F_i^{(HLL)}}\left(\bf{V}^{(n)}\right):= - \sum_{j=1}^4 \bf{F_{ij}^{(HLL)}}\left(\bf{V}^{(n)}\right) n_{ij} \triangle X %\Gamma_{ij}
$$
$n_{ij}$ being the outer normal of the control volumes edge's within the surface. The HLL flux $\bf{F_{ij}^{(HLL)}}\left(\bf{V}^{(n)}\right)$ approximates the local Riemann problem at edge $j$ of cell $i$. It assumes the local solution as to consist of two waves with one intermediate state. The wave speed estimates are denoted by $s_{\text{R}}$ for the fast wave and $s_{\text{L}}$ for the slow wave, respectively. With $\bf{V_{L}}^{(n)}$ and $\bf{V_{R}}^{(n)}$ being the reconstructions of $\bf{V}^{(n)}$ on the left and right sides of the interface $ij$, the numerical flux is given by
$$
\bf{F_{ij}^{(HLL)}}\left(\bf{V}^{(n)}\right)= \begin{cases} \bf{F}\left(\bf{V_{L}}^{(n)}\right)& 0 \leq s_{\text{L}} \\[0.4cm]
\frac{ s_{\text{R}}\,\bf{F}\left(\bf{V_{L}}^{(n)}\right) -s_{\text{L}}\,\bf{F}\left(\bf{V_{R}}^{(n)}\right) + s_{\text{R}}s_{\text{L}} \, \left( \bf{V_{R}}^{(n)}-\bf{V_{L}}^{(n)}\right)}{s_{\text{R}}-s_{\text{L}}} & s_{\text{L}} \leq 0 \leq s_{\text{R}} \\[0.4cm]
\bf{F}\left(\bf{V_{R}}^{(n)}\right) & s_{\text{R}} \leq 0
\end{cases}
$$
The wave speeds $s_{\text{R}}$ and $s_{\text{L}}$, are defined as proposed by Toro (1992), and hence account for dry states (interfaces with zero height on one side). To ensure stability the time-step is dynamically chosen according to the CFL-condition and slope limitation is applied to the linear reconstruction of the cell values by means of a Minmod limiter. The empirical convergence has been extensively tested in \cite{Kowalski08}. For a periodic test function and a central difference linear reconstruction, we showed a straight second order convergence (EOC $\approx 1.94$), whereas with a Minmod limitation it was somewhat smaller (EOC $\approx 1.67$).