IMEX_SfloW2D 1.0: a depth-averaged numerical flow model for pyroclastic avalanches

Pyroclastic avalanches are a type of granular flow generated at active volcanoes by different mechanisms, including the collapse of steep pyroclastic deposits (e.g., scoria and ash cones), fountaining during moderately explosive eruptions, and crumbling and gravitational collapse of lava domes. They represent end-members of gravity-driven pyroclastic flows characterized by relatively small volumes (less than about 1 Mm³) and relatively thin (1–10 m) layers at high particle concentration (10–50 vol %), manifesting strong topographic control. The simulation of their dynamics and mapping of their hazards pose several different problems to researchers and practitioners, mostly due to the complex and still poorly understood rheology of the polydisperse granular mixture and to the interaction with the complex natural three-dimensional topography, which often causes rapid rheological changes. In this paper, we present IMEX_SfloW2D, a depth-averaged flow model describing the granular mixture as a single-phase granular fluid. The model is formulated in absolute Cartesian coordinates (whereby the fluid flow equations are integrated along the direction of gravity) and can be solved over a topography described by a digital elevation model. The numerical discretization and solution algorithms are formulated to allow for a robust description of wet–dry conditions (thus allowing us to accurately track the front propagation) and an implicit solution to the nonlinear friction terms. Owing to these features, the model is able to reproduce steady solutions, such as the triggering and stopping phases of the flow, without the need for empirical conditions. Benchmark cases are discussed to verify the numerical code implementation and to demonstrate the main features of the new model. A preliminary application to the simulation of the 11 February pyroclastic avalanche at the Etna volcano (Italy) is finally presented. In the present formulation, a simple semi-empirical friction model (Voellmy–Salm rheology) is implemented. However, the modular structure of the code facilitates the implementation of more specific and calibrated rheological models for pyroclastic avalanches.

1 Introduction

Pyroclastic avalanches are rapid flows of pyroclastic material (volcanic ash, lapilli, pumices and scoriae) that propagate down volcanic slopes under the effect of gravity. They share many phenomenological features with other natural granular avalanches (such as landslides and rock and snow avalanches) and they pose similar modeling, monitoring and risk mitigation challenges (Pudasaini and Hutter, 2007).

Despite the fact that the term pyroclastic avalanche is not widely used among volcanologists, who often adopt the term pyroclastic density current (PDC), there are reasons to prefer the former in some specific cases. The term avalanche derives from the old french avaler, which meant “move down the valley”. Implicit in it is the idea that a pyroclastic avalanche must move along and be confined within the volcanic slopes and be dominantly driven by the longitudinal (i.e., parallel to the ground) component of gravity. Pyroclastic avalanches stop when either the slope is reduced or their momentum is dissipated by friction. This is opposite to a pyroclastic density current, a term deriving from the general concept of density current (Von Karman, 1940; Benjamin, 1968; Simpson, 1999), which is a flow moving under the dominant action of the hydrostatic pressure associated with its density contrast with respect to the atmosphere. Thick density currents are able to propagate inertially even on flat topographies, and the effect of friction is usually negligible. Low-aspect-ratio ignimbrites (Fisher et al., 1993; Dade and Huppert, 1996; Dade, 2003) or flows produced by the collapse of Plinian columns (e.g., Shea et al., 2011) can generally be (and have been) modeled as inertial PDCs for most of their run-out. PDCs indeed generally stop because of the progressive increase in buoyancy due to the combined effect of air entrainment and particle sedimentation (Bursik and Woods, 1996; Esposti Ongaro et al., 2016), which may give rise to the formation of co-ignimbrite eruption columns (Woods and Wohletz, 1991; Engwell et al., 2016), or because they cannot overcome topographic obstacles (Woods et al., 1998). The two phenomenologies often overlap. The basal, concentrated part of a PDC can be strongly topographically controlled and behave more like an avalanche (classically termed a pyroclastic flow). Although modeling of the latter has been done with an approach similar to that presented in this paper (e.g., Pitman et al., 2003, and derived works), it is becoming apparent that finding the appropriate rheological model is problematic (Kelfoun, 2011) and that the basal concentrated and the upper dilute turbulent layers should be coupled dynamically (Doyle et al., 2010; Kelfoun, 2017a; Shimizu et al., 2017). Pyroclastic flows can indeed show unexpected flow transformations, marking a transition from a dominantly frictional to an inertial behavior, in which mobility is drastically enhanced (Fisher, 1995; Komorowski et al., 2013).

We will refer to pyroclastic avalanches in this work for pyroclastic flows that (1) remain confined within the volcano slopes, (2) show evidence of a dense basal granular flow and (3) are controlled by topography (i.e., they mostly move in the direction of the maximum slope). Such conditions generally reduce the applicability of this category to relatively small flows (less than about 1 million m³) generated by mildly explosive activity (e.g., Strombolian) or by the gravitational collapse of basaltic scoria cones or of relatively small viscous and degassed lava domes. This is also the volume threshold identified by Ogburn and Calder (2017) above which modeling of pyroclastic currents becomes more problematic, showing transitional features between the two phenomenologies. The numerical model presented in this paper can be extended to be applicable to more general PDCs or pyroclastic flows sensu lato, but some aspects of the model formulation should be refined in that case. In particular, for what concerns the rheological model, the depositional–erosional features, the entrainment of atmospheric air and the thermodynamic properties. This will be the objective of future studies.

1.1 Modeling and numerical simulation of shallow pyroclastic avalanches

As in classical fluid dynamics, even more so for granular fluids, the choice is between continuum and discrete field representation (Guo and Curtis, 2015). In this work, we prefer the continuum approach, which is more suited to large geophysical systems (that would otherwise require a prohibitive number of discrete particles). As in similar models already considered in volcanological research and applications (Pitman et al., 2003; Kelfoun and Druitt, 2005; Shimizu et al., 2017), we also adopt a physical formulation based on depth averaging, which is appropriate for shallow granular avalanches and is computationally less expensive. Finally, the model is formulated for a single granular fluid. Future developments and implementations will consider multiphase flows as a more accurate representation of pyroclastic avalanches (Dufek, 2015).

Despite these simplifying hypotheses, several difficulties arise in granular avalanche depth-averaged models. On the one hand, terrain-following coordinates are often used to express depth-averaged transport equations. However, on 3-D rough surfaces, they need to be corrected with curvature terms, which introduce problems with irregular topographies, cliffs, obstacles and high curvatures (Denlinger and Iverson, 2004; Fischer et al., 2012). Moreover, on steep slopes, where acceleration along $\widehat{z}$ is non-negligible, the hydrostatic approximation is flawed (Denlinger and Iverson, 2004; Castro-Orgaz et al., 2015; Yuan et al., 2017).

On the other hand, from a physical point of view, the description of the depth-averaged rheology of the granular fluid was revealed to be problematic for strongly stratified and nonhomogeneous flows (Bartelt et al., 2016; Kelfoun, 2017a; Shimizu et al., 2017). Even in the more recent literature, a unifying model for the rheology of fast granular flows is still lacking (Bartelt et al., 1999; Iverson and Denlinger, 2001; Mangeney et al., 2007; Forterre and Pouliquen, 2008; Kelfoun, 2011; Iverson and George, 2014; Lucas et al., 2014; Delannay et al., 2017).

In addition to the latter, some additional difficulties arise from the numerical solution to the conservation equations: despite the fact that numerical methods based on conservative, approximate Riemann solvers are robust and well tested (Denlinger and Iverson, 2001; Mangeney-Castelnau et al., 2003; Patra et al., 2005; Christen et al., 2010; Toro, 2013), non-hydrostatic terms arising from the vertical momentum equation (Denlinger and Iverson, 2004) can be computationally expensive and/or need particular treatment. In many cases, further difficulties arise in the treatment of source terms (especially for thin flow in which friction dominates), often requiring empirical yielding–stopping criteria that might (and usually do) deteriorate the numerical solution (Charbonnier and Gertisser, 2009; Ogburn and Calder, 2017). Last, but not least, there are just a few open-source codes, and those available are not easy to modify due to the lack of documentation.

In this paper, we present and show verification tests of the new IMEX_SfloW2D (which stands for Implicit–Explicit Shallow Water flow) numerical model for shallow granular avalanches that we designed to address most of the above difficulties. In particular, the model is formulated in a geographical (absolute) coordinate system so that it is possible to include non-hydrostatic terms arising from steep topographic slopes, rapid topographic changes (e.g., local curvature) or from more accurate approximation of the vertical momentum equations. In this version of the model, however, non-hydrostatic terms are neglected. The model can deal with different initial and boundary conditions, but its first aim is to treat gravitational flows over topographies described as digital elevation models (DEMs) in the ESRI ASCII format. The same format is used for the output of the model so that it can be handled very easily with GIS software.

Several different geophysical flow depth-averaged models already exist, which are able to simulate pyroclastic avalanches on DEMs. Among them, two of the most well-established models in the volcanological community are VolcFlow and Titan2D. VolcFlow is a MATLAB finite-difference Eulerian code based on an explicit upwind or double-upwind discretization. Recently, A two-layer depth-averaged version for both the dilute and the concentrated parts of pyroclastic currents has been implemented (Kelfoun, 2017b). Titan2D is a parallel finite-volume Eulerian code based on a high-order slope-limiting Godunov scheme, and it allows for the use of adaptive mesh refinement, reducing the computational cost while maintaining the simulation's accuracy (Patra et al., 2005).

Numerically, the first and most relevant advancement in IMEX_SfloW2D is represented by the implicit treatment of the source terms in the transport equations, which avoids most problems related to the stopping of the flow, especially when dealing with strongly nonlinear rheologies. Any rheological model can in principle be implemented, including formulations not dependent on velocity, as purely frictional or plastic rheologies without the need for introducing an artificial numerical viscosity. The model is indeed discretized in time with an explicit–implicit Runge–Kutta method whereby the hyperbolic part and the source term associated with topography slope are solved explicitly, while other terms (friction) are treated implicitly. The finite-volume solver for the hyperbolic part of the system is based on the Kurganov and Petrova (2007) semi-discrete central-upwind scheme and it is not tied to the knowledge of the eigenstructure of the system of equations. The implicit part is solved with a Newton–Raphson method whereby the elements of the Jacobian on the nonlinear system are evaluated numerically with a complex-step derivative technique. This automatic procedure allows for the use of different formulations of the friction term without the need for major modifications of the code. In particular, the Voellmy–Salm empirical model is implemented in the present version.

The FORTRAN90 code can be freely downloaded and it is designed in a way that users can simply use it without any intervention or they can easily modify it by adding new transport and/or constitutive equations.

2 Physical and mathematical model

In this section we present the governing equations based on the shallow water approximation. We omit their derivation that comes from the manipulation of the mass conservation law and Newton's second equation of motion.

2.1 Depth-averaged equations

The model we use for the flow evolution is described by the Saint-Venant equations (Pudasaini and Hutter, 2007; Toro, 2013) coupled with source terms modeling frictional forces. Saint-Venant equations are partial differential equations suitable when the flow horizontal length scale is much greater than the vertical one, allowing us to disregard vertical flow motion. Here we write the equations in global Cartesian coordinates, and thus the topography can be expressed as B(x,y) (we assume it does not change with time t) and the two velocities are defined as the components along the x and y axes orthogonal to the z axis parallel to the gravitational acceleration $\mathit{g}=(\mathrm{0},\mathrm{0},g)$. We denote the z-averaged horizontal velocities with $u(x,y,t)$ and $v(x,y,t)$ and the flow depth with $h(x,y,t)$. In addition, in this work we assume a hydrostatic pressure distribution, resulting in the following relationship between pressure p, bulk density of the flow ρ (assumed to be constant) and flow depth h:

$$\begin{array}{}\text{(1)}& {\displaystyle }p=\mathit{\rho }gh.\end{array}$$

With these assumptions, and without considering frictional forces, the 2-D inviscid depth-averaged equations in differential form can be written in the following way:

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{\partial \left(hu\right)}{\partial x}}+{\displaystyle \frac{\partial \left(hv\right)}{\partial y}}=\mathrm{0},\text{(3)}& {\displaystyle }& {\displaystyle \frac{\partial \left(hu\right)}{\partial t}}+{\displaystyle \frac{\partial \left(h{u}^{\mathrm{2}}\right)}{\partial x}}+{\displaystyle \frac{\partial \left(huv\right)}{\partial y}}+gh{\displaystyle \frac{\partial \left(h+B\right)}{\partial x}}=\mathrm{0},\text{(4)}& {\displaystyle }& {\displaystyle \frac{\partial \left(hv\right)}{\partial t}}+{\displaystyle \frac{\partial \left(huv\right)}{\partial x}}+{\displaystyle \frac{\partial \left(h{v}^{\mathrm{2}}\right)}{\partial y}}+gh{\displaystyle \frac{\partial \left(h+B\right)}{\partial y}}=\mathrm{0}.\end{array}$$

The first equation represents the conservation of mass (or volume because of the constant density), while the other two equations describe the conservation of momentum in the x and y directions. The last terms in Eqs. (3)–(4) account for the gravity-induced force, which results from the hydrostatic pressure and depends on the slope of the free surface, and thus also on the topography elevation B(x,y). At present, no terms accounting for centrifugal acceleration effects caused by terrain curvature are present in model.

Table 1 List of model variables with notation and units.

As stated above, the variables of the model are $(h,u,v)$, but to deal more easily with the topography we introduce the additional variable $w=h+B$, describing the height of the free surface of the flow. We observe that, if we substitute h with w in the temporal derivative of Eq. (2), this still represents mass conservation because we assume that B is not changing with time.

If we introduce the vector of conservative variables $\mathit{Q}=(w,hu,hv{)}^T$ (where the superscript notation Q⁽ⁱ⁾ is used to denote the ith component), the governing equations can be written in the compact form:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\mathit{Q}}_t+\mathit{F}(\mathit{Q}{)}_x+\mathit{g}(\mathit{Q}{)}_y={\mathit{S}}_{\mathrm{1}}\left(\mathit{Q}\right),\end{array}$$

where the letter-type subscript denotes the partial derivative with respect to the correspondent variable, and the terms appearing in the equation are defined in the following way:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathit{F}\left(\mathit{Q}\right)=(hu,h{u}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}g{h}^{\mathrm{2}},huv{)}^T,\\ \text{(6)}& {\displaystyle }& {\displaystyle }\mathit{g}\left(\mathit{Q}\right)=(hv,huv,h{v}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}g{h}^{\mathrm{2}}{)}^T,{\displaystyle }& {\displaystyle }{\mathit{S}}_{\mathrm{1}}\left(\mathit{Q}\right)=(\mathrm{0},gh{B}_x,gh{B}_y{)}^T.\end{array}$$

We observe that the homogeneous part of Eq. (5) represents the shallow water equations over a flat topography, for which the eigenstructure and hyperbolicity are well studied (Toro, 2013). Furthermore, the homogeneous part is written in a conservative form, allowing us to easily adopt finite-volume discretization schemes for its numerical solution.

2.2 Voellmy–Salm rheology

To properly model shallow pyroclastic avalanches, we have to modify the classic Saint-Venant equations by introducing an additional source term S₂ accounting for friction forces (Pudasaini and Hutter, 2007; Toro, 2013). To this purpose, we implemented in IMEX_SfloW2D, as a prototype nonlinear model, the Voellmy–Salm rheology, which simulates mixture motion as a homogeneous mass flow. This model is commonly used for avalanches and debris flows, but it is also relevant for volcanic gravitational flows.

The terms we consider appear only in the momentum equations ($S_{\mathrm{2}}^{\left(\mathrm{1}\right)}=\mathrm{0}$) and are given by

$$\begin{array}{l}{\displaystyle }{S}_{\mathrm{2}}^{\left(\mathrm{2}\right)}\left(\mathit{Q}\right)={\displaystyle \frac{u}{\sqrt{u^{\mathrm{2}}+v^{\mathrm{2}}}}}\left[\mathit{\mu }h\mathit{g}\cdot \mathit{n}+{\displaystyle \frac{g}{\mathit{\xi }}}(u^{\mathrm{2}}+v^{\mathrm{2}})\right],\\ \text{(7)}& {\displaystyle }{\displaystyle }{S}_{\mathrm{2}}^{\left(\mathrm{3}\right)}\left(\mathit{Q}\right)={\displaystyle \frac{v}{\sqrt{u^{\mathrm{2}}+v^{\mathrm{2}}}}}\left[\mathit{\mu }h\mathit{g}\cdot \mathit{n}+{\displaystyle \frac{g}{\mathit{\xi }}}(u^{\mathrm{2}}+v^{\mathrm{2}})\right].\end{array}$$

The total basal friction in the Voellmy–Salm model (represented by the common term in square brackets in the two right-hand sides of Eq. 7) is split into two components: (1) a velocity-independent dry Coulomb friction proportional to the coefficient μ, the flow thickness and the component of the gravitational acceleration normal to the topography; and (2) a velocity-dependent turbulent friction inversely proportional to the coefficient ξ and commonly considered to represent the effect of granular collisions. For the sake of simplicity, μ is named the friction coefficient and ξ the turbulent coefficient. If the topography is a function of the global coordinates $z=B(x,y)$, then the component of the gravitational acceleration normal to the topography is given by

$$\begin{array}{}\text{(8)}& {\displaystyle }\mathit{g}\cdot \mathit{n}={\displaystyle \frac{g}{\sqrt{\mathrm{1}+B_x^{\mathrm{2}}+B_y^{\mathrm{2}}}}}.\end{array}$$

With the notation introduced above, the final form of the equations modeling pyroclastic avalanches is

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathit{Q}}_t+\mathit{F}(\mathit{Q}{)}_x+\mathit{g}(\mathit{Q}{)}_y={\mathit{S}}_{\mathrm{1}}\left(\mathit{Q}\right)+{\mathit{S}}_{\mathrm{2}}\left(\mathit{Q}\right),\end{array}$$

where the hyperbolic terms and the source terms are defined by Eqs. (6)–(7). In this formulation we have kept the two source terms separated because, as shown in the next section, they require different numerical treatments.

3 Numerical model

IMEX_SfloW2D is based on a finite-volume central-upwind scheme in space and on an implicit–explicit Runge–Kutta scheme for the discretization in time. The main purpose of the code is to run simulations on colocated grids derived from DEMs, and for this reason the standard input files defining the topography are raster files in the ESRI ASCII format, defining a uniform grid of equally sized square pixels whose values (in our case representing the terrain elevation above sea level) are arranged in rows and columns. The procedure to define the elevation values at the face centers and cell centers of the computational grid is represented in Fig. 1. First, the pixel values of the ESRI file (represented with colored squares in Fig. 1) are allocated to the coordinates of the center of the pixels (filled circles in Fig. 1), and then these values are linearly interpolated at the four corners of the computational grid (no-fill circles in Fig. 1). Finally, the elevation values B_j,k at the centers of each cell (filled squares in Fig. 1) are defined as the average value of the four cell corners, while the values at the centers of each face, denoted with $B_{j,k+\frac{\mathrm{1}}{\mathrm{2}}}$ and $B_{j+\frac{\mathrm{1}}{\mathrm{2}},k}$ and represented by the no-fill squares in Fig. 1, are defined as the average value of the two face corners. With this definition, B_j,k will also be the average of the values at the centers of the four faces of the cell (i,j) and this fact plays an important role for a correct numerical discretization of the last two terms in Eqs. (3)–(4), resulting in a scheme capable of preserving steady states of the form $h+B=$ const. This interpolation procedure does not degrade the original DEM, provided that the computational grid size has the same resolution. In fact, if they have the same resolution and the computational grid corners coincide with the center of the DEM's pixels, no interpolation is done.

Figure 1 Computational grids. The colored pixels represent the elevation values of the original DEM. The lines define the edges of the IMEX-SfloW2D computational cells. The elevation values at the centers (filled squares), faces (no-fill squares) and corners (no-fill circles) of the computational cells are obtained by interpolating the pixel values associated with their centers (filled circles).

3.1 Central-upwind scheme

The finite-volume method adopted for IMEX_SfloW2D is based on the semi-discrete central-upwind scheme introduced in Kurganov and Petrova (2007), in which the term central refers to the fact that the numerical fluxes at each cell interface are based on an average of the fluxes at the two sides of the interface; the term upwind is employed because, in the flux averaging, the weights depend on the local speeds of propagation at the interface. Following Kurganov and Petrova (2007), the semi-discretization in space leads to the following ordinary differential equation system in each cell:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{d}{\mathrm{d}t}}{\stackrel{\mathrm{‾}}{\mathit{Q}}}_{j,k}\left(t\right)=-{\displaystyle \frac{{\mathit{H}}_{j+\frac{\mathrm{1}}{\mathrm{2}},k}^x\left(t\right)-{\mathit{H}}_{j-\frac{\mathrm{1}}{\mathrm{2}},k}^x\left(t\right)}{\mathrm{\Delta }x}}\\ \text{(10)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{{\mathit{H}}_{j,k+\frac{\mathrm{1}}{\mathrm{2}}}^y\left(t\right)-{\mathit{H}}_{j,k-\frac{\mathrm{1}}{\mathrm{2}}}^y\left(t\right)}{\mathrm{\Delta }y}}+{\stackrel{\mathrm{‾}}{\mathit{S}}}_{j,k}\left(t\right),\end{array}$$

where ${\stackrel{\mathrm{‾}}{\mathit{Q}}}_{j,k}$ denotes the average of the conservative variables Q(x,y) over the control volume (j,k), and H^x and H^y are the numerical fluxes calculated from the value of the variables reconstructed at the cell interfaces.

The choice of the variables to reconstruct at the interface is fundamental for the stability of the numerical scheme. The homogeneous system associated with Eq. (9) admits smooth steady-state solutions, as well as non-smooth steady-state solutions. A good numerical method for the solution to the homogeneous system should accurately capture both the steady-state solutions and their small perturbations (quasi-steady flows). From a practical point of view, one of the most important steady-state solutions is a stationary one:

$$\begin{array}{}\text{(11)}& {\displaystyle }w=h+B=\mathrm{const},u=v=\mathrm{0}.\end{array}$$

This suggests using, as the vector of variables for the linear reconstruction at the interfaces, the vector $\mathit{U}=(w,u,v{)}^T$, denoted as the vector of physical variables of the system, for which the boundary conditions are also prescribed. IMEX_SfloW2D also allows the user to choose a set of variables for the linear reconstruction of the vector of conservative variables $\mathit{Q}=(w,hu,hv{)}^T$. In this case, a correction procedure is required to limit the values of the velocity components at the interfaces when the flow thickness goes to zero, as done in Kurganov and Petrova (2007).

For the reconstruction procedure based on the physical variables, we introduce the notation $\mathbf{\Gamma }:{\mathbb{R}}^{\mathrm{3}}\to {\mathbb{R}}^{\mathrm{3}}$ for the mapping from conservative variables Q to physical variables U and Γ⁻¹ for the inverse mapping from physical to conservative variables. From the average values of the physical variables we can operate a linear reconstruction inside each cell in order to obtain the values at the interfaces sides. In particular, given the local partial derivatives at the cell center (U_x)_j,k and (U_y)_j,k, the one-side values at the east, west, north and south interfaces of the cell (j,k) are given by

$$\begin{array}{l}{\displaystyle }{\mathit{U}}_{j,k}^E={\stackrel{\mathrm{‾}}{\mathit{U}}}_{j,k}+{\displaystyle \frac{\mathrm{\Delta }x}{\mathrm{2}}}({\mathit{U}}_x{)}_{j,k},{\mathit{U}}_{j,k}^W={\stackrel{\mathrm{‾}}{\mathit{U}}}_{j,k}-{\displaystyle \frac{\mathrm{\Delta }x}{\mathrm{2}}}({\mathit{U}}_x{)}_{j,k},\\ {\displaystyle }{\mathit{U}}_{j,k}^N={\stackrel{\mathrm{‾}}{\mathit{U}}}_{j,k}+{\displaystyle \frac{\mathrm{\Delta }x}{\mathrm{2}}}({\mathit{U}}_y{)}_{j,k},{\mathit{U}}_{j,k}^S={\stackrel{\mathrm{‾}}{\mathit{U}}}_{j,k}-{\displaystyle \frac{\mathrm{\Delta }x}{\mathrm{2}}}({\mathit{U}}_y{)}_{j,k}.\end{array}$$

These partial derivatives are calculated using an opportune geometric limiter. In IMEX_SfloW2D it is possible to choose between MinMod, Superbee and Van Leer limiters. We observe that at each cell interface and for each variable, there are two reconstructed values, one from each cell at the two sides of the interface.

During the reconstruction step, particular care should be taken to avoid unrealistic values of the physical variables, such as negative flow thickness or velocities that are too large. For this reason, in the case that one of the reconstructed interface values of w is smaller than the topography B at the same location (thus resulting in a negative thickness h), the relative derivative is further limited to have a zero thickness at such an interface. We remark that the correction is applied to the derivative, and thus the reconstructed value of w at the opposite interface will also be affected. For example, if $w_{j,k}^S<B_{j,k-\frac{\mathrm{1}}{\mathrm{2}}}$, then we take the following derivative in the y direction:

$$\begin{array}{}\text{(12)}& {\displaystyle }{\left(w_y\right)}_{j,k}={\displaystyle \frac{{\overline{w}}_{j,k}-B_{j,k-\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{\Delta }y/\mathrm{2}}},\end{array}$$

which gives the two reconstructions at the S and N interfaces of the (i,j) control volume:

$$\begin{array}{}\text{(13)}& {\displaystyle }{w}_{j,k}^S=B_{j,k-\frac{\mathrm{1}}{\mathrm{2}}},w_{j,k}^N=\mathrm{2}{\overline{w}}_{j,k}-B_{j,k-\frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

As stated above, an additional problem in the reconstruction step is that h can be very small, or even zero, leading to large values of the velocities. Thus, when the physical variables u and v at the interfaces are computed from the conservative variables, a desingularization is applied to avoid division by very small numbers and the corrected values are given by the following formulas:

$$\begin{array}{}\text{(14)}& {\displaystyle }u={\displaystyle \frac{\sqrt{\mathrm{2}}h\left(hu\right)}{\sqrt{h^{\mathrm{4}}+max(h^{\mathrm{4}},\mathit{ϵ})}}},v={\displaystyle \frac{\sqrt{\mathrm{2}}h\left(hv\right)}{\sqrt{h^{\mathrm{4}}+max(h^{\mathrm{4}},\mathit{ϵ})}}},\end{array}$$

where ϵ is a prescribed tolerance.

Finally, once the physical variables are reconstructed at the interfaces, the numerical fluxes in the x direction are given by

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{H}}_{j+\frac{\mathrm{1}}{\mathrm{2}},k}^x=\\ {\displaystyle }& {\displaystyle }{\displaystyle \frac{a_{j+\frac{\mathrm{1}}{\mathrm{2}},k}^{+}\mathit{F}\left({\mathbf{\Gamma }}^{-\mathrm{1}}\right({\mathit{U}}_{j,k}^E),B_{j+\frac{\mathrm{1}}{\mathrm{2}},k})-a_{j+\frac{\mathrm{1}}{\mathrm{2}},k}^{-}\mathit{F}\left({\mathbf{\Gamma }}^{-\mathrm{1}}\right({\mathit{U}}_{j+\mathrm{1},k}^W),B_{j+\frac{\mathrm{1}}{\mathrm{2}},k})}{a_{j+\frac{\mathrm{1}}{\mathrm{2}},k}^{+}-a_{j+\frac{\mathrm{1}}{\mathrm{2}},k}^{-}}}\\ \text{(15)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{a_{j+\frac{\mathrm{1}}{\mathrm{2}},k}^{+}{a}_{j+\frac{\mathrm{1}}{\mathrm{2}},k}^{-}}{a_{j+\frac{\mathrm{1}}{\mathrm{2}},k}^{+}-a_{j+\frac{\mathrm{1}}{\mathrm{2}},k}^{-}}}\left({\mathbf{\Gamma }}^{-\mathrm{1}}\left({\mathit{U}}_{j+\mathrm{1},k}^W\right)-{\mathbf{\Gamma }}^{-\mathrm{1}}\left({\mathit{U}}_{j,k}^E\right)\right),\end{array}$$

where the right- and left-sided local speeds $a_{j+\frac{\mathrm{1}}{\mathrm{2}},k}^{+}$ and $a_{j+\frac{\mathrm{1}}{\mathrm{2}},k}^{-}$ are estimated by

$$\begin{array}{l}{\displaystyle }{a}_{j+\frac{\mathrm{1}}{\mathrm{2}},k}^{+}=max(u_{j,k}^E+\sqrt{g{h}_{j,k}^E},u_{j+\mathrm{1},k}^W+\sqrt{g{h}_{j+\mathrm{1},k}^W},\mathrm{0}),\\ {\displaystyle }{a}_{j+\frac{\mathrm{1}}{\mathrm{2}},k}^{-}=min(u_{j,k}^E-\sqrt{g{h}_{j,k}^E},u_{j+\mathrm{1},k}^W-\sqrt{g{h}_{j+\mathrm{1},k}^W},\mathrm{0}).\end{array}$$

In a similar way, the numerical fluxes in the y direction are given by

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{H}}_{j,k+\frac{\mathrm{1}}{\mathrm{2}}}^y=\\ {\displaystyle }& {\displaystyle }{\displaystyle \frac{b_{j,k+\frac{\mathrm{1}}{\mathrm{2}}}^{+}\mathit{G}\left({\mathbf{\Gamma }}^{-\mathrm{1}}\right({\mathit{U}}_{j,k}^N),B_{j,k+\frac{\mathrm{1}}{\mathrm{2}}})-b_{j,k+\frac{\mathrm{1}}{\mathrm{2}}}^{-}\mathit{G}\left({\mathbf{\Gamma }}^{-\mathrm{1}}\right({\mathit{U}}_{j,k+\mathrm{1}}^S),B_{j,k+\frac{\mathrm{1}}{\mathrm{2}}})}{b_{j,k+\frac{\mathrm{1}}{\mathrm{2}}}^{+}-b_{j,k+\frac{\mathrm{1}}{\mathrm{2}}}^{-}}}\\ \text{(16)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{b_{j,k+\frac{\mathrm{1}}{\mathrm{2}}}^{+}{b}_{j,k+\frac{\mathrm{1}}{\mathrm{2}}}^{-}}{b_{j,k+\frac{\mathrm{1}}{\mathrm{2}}}^{+}-b_{j,k+\frac{\mathrm{1}}{\mathrm{2}}}^{-}}}\left({\mathbf{\Gamma }}^{-\mathrm{1}}\left({\mathit{U}}_{j,k+\mathrm{1}}^S\right)-{\mathbf{\Gamma }}^{-\mathrm{1}}\left({\mathit{U}}_{j,k}^N\right)\right),\end{array}$$

where local speeds in the y direction $b_{j,k+\frac{\mathrm{1}}{\mathrm{2}}}^{+}$ and $b_{j,k+\frac{\mathrm{1}}{\mathrm{2}}}^{-}$ are given by

$$\begin{array}{l}{\displaystyle }{b}_{j,k+\frac{\mathrm{1}}{\mathrm{2}}}^{+}=max(v_{j,k}^N+\sqrt{g{h}_{j,k}^N},v_{j,k+\mathrm{1}}^S+\sqrt{g{h}_{j,k+\mathrm{1}}^S},\mathrm{0}),\\ {\displaystyle }{b}_{j,k+\frac{\mathrm{1}}{\mathrm{2}}}^{-}=min(v_{j,k}^N-\sqrt{g{h}_{j,k}^N},v_{j,k+\mathrm{1}}^S-\sqrt{g{h}_{j,k+\mathrm{1}}^S},\mathrm{0}).\end{array}$$

Following Kurganov and Petrova (2007), the source term S₁ is calculated trough a quadrature formula:

$$\begin{array}{}\text{(17)}& {\displaystyle }{S}_{\mathrm{1}}^{\left(\mathrm{2}\right)}=-g(w_{j,k}-B_{j,k}){\displaystyle \frac{B_{j+\frac{\mathrm{1}}{\mathrm{2}},k}-B_{j-\frac{\mathrm{1}}{\mathrm{2}},k}}{\mathrm{\Delta }x}},\text{(18)}& {\displaystyle }{S}_{\mathrm{1}}^{\left(\mathrm{3}\right)}=-g(w_{j,k}-B_{j,k}){\displaystyle \frac{B_{j,k+\frac{\mathrm{1}}{\mathrm{2}}}-B_{j,k-\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{\Delta }y}}.\end{array}$$

This discretization, coupled with the fact that the elevation B_j,k at the center of the control volumes is defined as the average value of the elevation at the center of the faces, guarantees that the resulting second-order numerical scheme is well balanced (i.e., preserves steady-state solutions) and the solutions are nonnegative (Kurganov and Petrova, 2007).

3.2 Runge–Kutta method

The semi-discrete system of Eq. (10) is solved using an implicit–explicit (IMEX), diagonally implicit Runge–Kutta scheme (DIRK) because such a scheme is well suited for solving stiff systems of partial differential equations, and the governing equations are expected to be stiff given the strong nonlinearities present in the friction terms. In addition, an implicit treatment of these terms allows for a better coupling of the equations and a proper recovery of the stoppage condition without the need to impose additional criteria or arbitrary thresholds.

The family of IMEX methods (Ascher et al., 1997) have been developed to solve stiff systems of partial differential equations written in the form

$$\begin{array}{}\text{(19)}& {\displaystyle }{\mathit{Q}}_t+\mathit{P}\left(\mathit{Q}\right)=\mathit{R}\left(\mathit{Q}\right),\end{array}$$

where in P are lumped all the non-stiff terms (in our case the semi-discretized conservative fluxes F and g and the term S₁), while R denotes the stiff terms of the system (here represented by the friction term S₂). The system of Eq. (19) must be solved for each control volume (j,k), but here, to keep the notation simpler, we omit the subscripts.

An IMEX Runge–Kutta with ν steps consists of applying an implicit discretization to the stiff terms and an explicit one to the non-stiff terms, obtaining

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{Q}}^{\left(l\right)}={\mathit{Q}}^n-\mathrm{\Delta }t\sum _{m=\mathrm{1}}^{l-\mathrm{1}}{\stackrel{\mathrm{̃}}{a}}_{lm}\mathit{P}\left({\mathit{Q}}^{\left(m\right)}\right)+\mathrm{\Delta }t\sum _{m=\mathrm{1}}^{\mathit{\nu }}{a}_{lm}\mathit{R}\left({\mathit{Q}}^{\left(m\right)}\right),\\ \text{(20)}& {\displaystyle }& {\displaystyle }l=\mathrm{1},\mathrm{\dots },\mathit{\nu }\text{(21)}& {\displaystyle }& {\displaystyle }{\mathit{Q}}^{n+\mathrm{1}}={\mathit{Q}}^n-\mathrm{\Delta }t\sum _{m=\mathrm{1}}^{\mathit{\nu }}{\stackrel{\mathrm{̃}}{b}}_m\mathit{P}\left({\mathit{Q}}^{\left(m\right)}\right)+\mathrm{\Delta }t\sum _{m=\mathrm{1}}^{\mathit{\nu }}{b}_m\mathit{R}\left({\mathit{Q}}^{\left(m\right)}\right).\end{array}$$

The choice of the number of the Runge–Kutta steps ν, of the ν×ν matrices $\stackrel{\mathrm{̃}}{A}=\left({\stackrel{\mathrm{̃}}{a}}_{lm}\right)$ and A=(a_lm), and of the vectors $\stackrel{\mathrm{̃}}{b}=({\stackrel{\mathrm{̃}}{b}}_{\mathrm{1}},\mathrm{\dots },{\stackrel{\mathrm{̃}}{b}}_{\mathit{\nu }})$ and $b=(b_{\mathrm{1}},\mathrm{\dots },b_{\mathit{\nu }})$ differentiates the various IMEX Runge–Kutta schemes. We remark that the explicit discretization of the non-stiff terms requires ${\stackrel{\mathrm{̃}}{a}}_{lm}=\mathrm{0}$ for l≥m, while the implicit treatment of the stiff terms requires a_lm≠0 for some l≥m.

Following Pareschi and Russo (2005), the IMEX scheme used in this work satisfies an additional condition, i.e., a_lm=0 for l>m. This family of IMEX schemes are called direct implicit Runge–Kutta (DIRK) schemes and their use leads to the following implicit problem to solve at each step in the Runge–Kutta procedure:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathit{N}\left({\mathit{Q}}^{\left(l\right)}\right)\equiv {\mathit{Q}}^{\left(l\right)}-\mathrm{\Delta }t\cdot a_{ll}\mathit{R}\left({\mathit{Q}}^{\left(l\right)}\right)-{\mathit{Q}}^n\\ \text{(22)}& {\displaystyle }& {\displaystyle }+\mathrm{\Delta }t\sum _{m=\mathrm{1}}^{l-\mathrm{1}}\left[{\stackrel{\mathrm{̃}}{a}}_{lm}\mathit{g}\left({\mathit{Q}}^{\left(m\right)}\right)+a_{lm}\mathit{R}\left({\mathit{Q}}^{\left(m\right)}\right)\right]=\mathrm{0}.\end{array}$$

To enforce the stopping condition that can result from the application of the total basal friction, the velocity-independent friction term and the velocity-dependent term are computed in two steps. First, the dry Coulomb friction is computed and its value is limited to account for the fact that this force at maximum can stop the flow. Then, the system of nonlinear Eq. (22) in the unknowns Q^(l) is solved using a Newton–Raphson method with an optimum step size control. The method requires the computation of the Jacobian matrix J of the left-hand side of Eq. (22) with the highest possible accuracy, since R(Q^(j)), accounting for the dependence of the friction force on flow variables, can be strongly nonlinear. Following Martins et al. (2003) and La Spina and de' Michieli Vitturi (2012), this can be obtained with the use of complex variables to estimate derivatives. With the complex-step derivative approximation we can approximate the Jacobian J needed for the Newton–Raphson method with an error of the same order as the machine working precision. We simply extend the function N to the complex plane, introducing the new function $\stackrel{\mathrm{̃}}{\mathit{N}}:{\mathbb{C}}^{\mathrm{3}}\to {\mathbb{C}}^{\mathrm{3}}$, and compute the real-valued columns of the Jacobian at Q as

$$\begin{array}{}\text{(23)}& {\displaystyle }\mathbf{J}\cdot {\mathit{e}}_j={\displaystyle \frac{\mathrm{\Im }\left(\stackrel{\mathrm{̃}}{\mathit{N}}\right(\mathit{Q}+i\mathit{ϵ}{\mathit{e}}_j\left)\right)}{\mathit{ϵ}}},\end{array}$$

where $({\mathit{e}}_j{)}_{j=\mathrm{1},\mathrm{\dots },\mathrm{3}}$ represents the standard basis vectors of ℝ³, ℑ(⋅) denotes the imaginary part of complex numbers and ϵ is a real number of the order of the machine working precision. Once the Jacobian is computed, the descent direction of the Newton–Raphson is updated and the descent step is obtained by applying a globally convergent method as described in Press et al. (1996).

Figure 2 1-D Riemann problem with no friction. Numerical solution at three different times: 0 s (a, b), 0.1 s (c, d), 0.2 s (e, f). Flow thickness and bottom topography (blue line) are plotted in the left panels, while the velocity is plotted in the right panels. The different colors represent numerical solutions obtained with different numbers of cells: 100 (dashed red line) and 400 (solid green line).

4 Numerical tests and code verification

In this section we present numerical tests aimed at demonstrating the mathematical accuracy of the numerical model results (the verification step, following Oberkampf and Trucano, 2002). Numerical tests are aimed at proving the following:

1. the capability to manage the propagation of discontinuities;

2. the potential to deal with complex and steep topographies and dry–wet interfaces; and

3. the ability of the granular avalanche to stop, thereby achieving the expected steady state.

All the numerical tests presented here are available on the Wiki page of the code (https://github.com/demichie/IMEX_SfloW2D/wiki, last access: 30 January 2019), together with animations of the output and scripts to reproduce the results.

4.1 One-dimensional test with discontinuous initial solution and topography

The first example is a 1-D test for a Riemann problem with a discontinuous topography, as presented in Andrianov (2004) and (Kurganov and Petrova, 2007). No frictional forces are considered in this test, aimed only at showing the capability of the numerical scheme used for the spatial discretization to properly model the propagation of strong discontinuities in both flow thickness and velocity.

The domain is the interval [0;1] and bottom topography B is the following step function:

$${\displaystyle }B\left(x\right)=\left\{\begin{array}{ll}\mathrm{2},& x\le \mathrm{0.5},\\ \mathrm{0.1},& x>\mathrm{0.5}.\end{array}\right.$$

The gravitational constant is g=2 and the initial data are

$${\displaystyle }\left(w\right(x,o),u(x,\mathrm{0}\left)\right)=\left\{\begin{array}{ll}(\mathrm{2.222},-\mathrm{1}),& x\le \mathrm{0.5},\\ (\mathrm{0.8246},-\mathrm{1.6359}),& x>\mathrm{0.5}.\end{array}\right.$$

The initial solution (Fig. 2, top panels) presents a discontinuity at x=0.5 and the exact solution at t>0 consists of a left-going rarefaction wave and a right-going shock wave. We present here the numerical solution obtained with a three-step IMEX Runge–Kutta scheme, whereby the reconstruction from the cell centers to the faces is applied to the physical variables. The solutions obtained by discretizing the domain with 100 cells (Fig. 2, dashed red lines) and 400 cells (Fig. 2, solid green lines) are compared. The numerical results at two different times t>0 (Fig. 2, middle and bottom panels) show the capacity of the numerical scheme to properly model the propagation of both the rarefaction and shock waves, and the good description of the shock with a small number of cells highlights the low numerical diffusivity of the central-upwind finite-volume numerical scheme implemented for the spatial discretization of the governing equations.

Figure 3 Numerical solution of a 1-D Riemann problem with friction at four different times. The solid blue line represents the topography, while the solid green line represents the free surface of the wet region.

4.2 One-dimensional problem with dry–wet interface and friction

This example is a 1-D test for a system with friction, as presented in Kurganov and Petrova (2007). As in the previous test, an initial discontinuity is present, but this time representing the interface from a “wet” (presence of flow) and a “dry” (no flow or zero thickness) region, with the terminology borrowed from the common use of Saint-Venant equations in hydrology. Thus, there is an additional numerical difficulty involving the capability of the numerical solver to propagate these discontinuities without creating regions with negative flow thickness. For this problem the bottom topography presents both smooth regions and a step, and it is defined as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }B\left(x\right)=\\ {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}\mathrm{1},& x\le \mathrm{0},\\ {\cos }^{\mathrm{2}}\left(\mathit{\pi }x\right),& \mathrm{0}\le x<\mathrm{0.4},\\ {\cos }^{\mathrm{2}}\left(\mathit{\pi }x\right)+\mathrm{0.25}\left(\cos \right(\mathrm{10}\mathit{\pi }(x-\mathrm{0.5}))+\mathrm{1}),& \mathrm{0.4}\le x<\mathrm{0.5},\\ \mathrm{0.5}{\cos }^{\mathrm{4}}\left(\mathit{\pi }x\right)+\mathrm{0.25}\left(\cos \right(\mathrm{10}\mathit{\pi }(x-\mathrm{0.5}))+\mathrm{1}),& \mathrm{0.5}\le x<\mathrm{0.6},\\ \mathrm{0.5}{\cos }^{\mathrm{4}}\left(\mathit{\pi }x\right),& \mathrm{0.6}\le x<\mathrm{1},\\ \mathrm{0.25}\sin \left(\mathrm{2}\mathit{\pi }\right(x-\mathrm{1}\left)\right),& \mathrm{1}\le x<\mathrm{1.5},\\ \mathrm{0},& x\ge \mathrm{1.5}.\end{array}\right.\end{array}$$

The initial conditions are defined by

$${\displaystyle }w(x,\mathrm{0})=\left\{\begin{array}{ll}\mathrm{1.4},& x\le \mathrm{0},\\ B\left(x\right),& x>\mathrm{0},\end{array}\right.u(x,\mathrm{0})=\mathrm{0}.$$

For this simulation, the gravitational constant has a value g=1 and a simpler friction term is employed (−κ(h)u, with $\mathit{\kappa }\left(h\right)=\mathrm{0.001}(\mathrm{1}+\mathrm{10}h{)}^{-\mathrm{1}}$). Also, for this test the numerical solution is obtained with a three-step IMEX Runge–Kutta scheme, whereby the reconstruction from the cell centers to the faces has been applied to the physical variables and the domain $[-\mathrm{0.25};\mathrm{1.75}]$ has been discretized with 400 cells. The boundary conditions are prescribed to model a closed domain and also to check that the total mass contained in the domain is kept constant by the numerical discretization schemes. The numerical solution at four different times (t=0, t=0.2, t=2 and t=40 s) is presented in Fig. 3, showing the capability of the model to deal with the propagation of dry–wet interfaces and to reach a steady solution for which the horizontal gradient of $w=h+B$ is null.

Figure 4 Numerical solution of a 1-D problem with Voellmy–Salm friction at three different times. In (a) the topography (blue line) and the profile of the material (green line) are plotted. In (b) the corresponding velocities are shown.

4.3 One-dimensional pyroclastic avalanche with Voellmy–Salm friction

This example is a 1-D test for the Voellmy–Salm rheology, with a pile of material initially at rest released on a constant slope topography. This test is aimed at checking if the model is able to preserve an initial steady condition, when the tangent of the pile free surface slope is smaller than the Coulomb friction coefficient, μ, and to properly simulate the stopping of the flow, i.e., when inertial and gravitational forces are smaller than total basal friction.

For this test, the domain is the interval [0;500] and the center of the pile coincides with the center of the domain. The gravitational constant is g=9.81 ms⁻² and the parameters of the Voellmy–Salm rheology are μ=0.3 and ξ=300 ms⁻².

We present the results for a numerical simulation with a topography with a constant slope of 13^∘ and an initial pile of material with a relative slope (i.e., with respect to the topography) of 20^∘. Thus, for this test, the initial condition is not steady and the pile of material starts to move along the slope until it reaches the stoppage condition.

The numerical solution at three different times (t=0, t=5 and t=40 s) is presented in Fig. 4. For this simulation, the reconstruction technique with limiters has been applied to the physical variables and an IMEX four-step Runge–Kutta has been adopted, while the domain has been discretized with 400 cells. The plot of the numerical solution at the intermediate time clearly shows that the front of the flow has a larger propagation velocity than the rest of the flow. On the other hand, comparing the left and right middle panels, we observe that part of the tail is not moving at all. After a few tens of seconds from the release, the flow reaches a steady condition, as shown by the plot of the velocity at 40 s (bottom right panel in Fig. 4). This test highlights the capability of the model to reach a steady condition not only when the gradient of $w=h+B$ is null, as shown in the previous example, but also with a flow with a positive slope below a critical condition.

Figure 5 Numerical simulation of two-dimensional pyroclastic avalanche with Voellmy–Salm friction. The contour plots on the bottom plane of each panel represent constant values of the thickness of the flow, whose free surface is represented in blue. The outermost contour on the bottom plane corresponds to a thickness of 0.06 m, and thus the thinner portion of the flow is not represented by the contour plot. A visual comparison between the two bottom plots highlights the fact that a steady condition has been reached.

4.4 Two-dimensional pyroclastic avalanche with Voellmy–Salm friction

This test extends the simulation with a Voellmy–Salm rheological model presented in the previous section from one to two dimensions, with an example of an avalanche of finite granular mass sliding down an inclined plane merging continuously into a horizontal one. The initial conditions and the topography of this tutorial are the same as in Example 4.1 from Wang et al. (2004). The computational domain is the rectangle $[\mathrm{0};\mathrm{30}]\times [-\mathrm{7};\mathrm{7}]$, in which a hemispherical shell holding the material together with a maximum thickness of 1.85 m is suddenly released so that the bulk material starts to slide on an inclined flat plane at 35^∘ (for $\mathrm{0}\le x\le \mathrm{17.5}$) into a horizontal run-out plane (for x≥21.5) connected by a smooth transition. Here we do not use the same rheological model as in the original paper of the example, but a Voellmy–Salm rheology is applied with μ=0.3 and ξ=300 ms⁻².

The numerical solution at four different times (t=0, t=7.5, t=20 and t=25 s) is presented in Fig. 5, in which both the three-dimensional flow shape over the topography and thickness contours are presented. For this simulation, the reconstruction technique with limiters has been applied to the physical variables and an IMEX two-step Runge–Kutta has been adopted, while the domain has been discretized with 150×100 cells.

Figure 6 Map of the avalanche deposit from the February 2014 pyroclastic avalanche and the lava flow.

As shown by the plots presented in Fig. 5, the model is able to simulate the propagation of the flow with no numerical oscillations or instabilities, without the need for artificial numerical diffusion. As the front reaches the maximum run-out and horizontal spreading (top left panel), the tail of the flow is still accelerating and the avalanche body starts to contract. Comparing the two panels at the bottom, it is evident that after about 20 s the flow has stopped propagating, with the deposit located at the transition region between the inclined and horizontal zones. This simulation took 238 s on an Intel Core i5-3210M CPU at 2.50 GHz.

Figure 7 Results of IMEX_SfloW2D numerical simulations of a pyroclastic avalanche overlapped on the hill-shaded relief of the Etna summit and the boundaries of the 11 February 2014 event. Numerical parameters as follows: (a) μ=0.1, ξ=500; (b) μ=0.2, ξ=500; (c) μ=0.2, ξ=100; (d) μ=0.3, ξ=500; (e) μ=0.4, ξ=500; (f) μ=0.4, ξ=5000. Avalanche volume is equal to 0.5×10⁶ m³.

5 Simulation of a pyroclastic avalanche at Etna volcano

On 11 February 2014, a hot pyroclastic avalanche was generated at the New Southeast Crater (NSEC) of Etna, triggered by the instability and collapse of its eastern flank where several vents had been actively effusing lava flows towards Valle del Bove since 22 January. The avalanche propagation was recorded by the INGV (Istituto Nazionale di Geofisica e Vulcanologia) monitoring system and, in particular, it was filmed by the thermal IR camera from Monte Cagliato, located on the east slope of the Valle del Bove at about 7 km from the NSEC, and the Catania CUAD visible camera (ECV) about 26 km south of the summit. A 500 m wide avalanche front propagated about 2.3 km along the steep slopes of the Valle del Bove before stopping at the break in slope at the valley bottom. At the same time, a voluminous buoyant ash cloud was generated by elutriation of the finest ash from the avalanche and rapidly dispersed in the north-northeast direction by an intense wind. The event is accurately reported by Andronico et al. (2018). The flank collapse left a detachment niche that allowed us to estimate the avalanche volume to range between 0.5 and 1.0×10⁶ m³. Due to difficult weather and environmental conditions, the presence of active lava flows in the region and the persistence of the Strombolian activity, it was unfortunately not possible to obtain a detailed map of the deposit thickness. Comparison of numerical results will therefore be limited to the flow boundary and run-out.

IMEX_SfloW2D numerical simulations have been performed over the 2014 digital elevation model of Mount Etna (De Beni et al., 2015). We have developed a MATLAB tool to modify the original topography by excavating a detaching volume with an ellipsoidal shape oriented towards the local maximum slope and with a prescribed total volume. The initial avalanche volume is defined by the difference between the original and the modified topography. It is also worth remarking that the DEM did not account for the presence of the thick lava flow that had been emplaced in the days before the avalanche event (Fig. 6), which likely controlled the avalanche path by confining it along its southern edge (Andronico et al., 2018).

The avalanche rheological parameters have been varied in ranges consistent with previous studies of geophysical granular avalanches (e.g., Bartelt et al., 1999). For quasi-static granular flows, μ is physically related to the basal friction of the granular material. However, its value for rapid avalanches cannot be easily defined a priori. We thus simply consider μ and ξ as empirical model parameters whose values need to be calibrated. In particular, $\mathrm{0.2}<\mathit{\mu }<\mathrm{0.5}$ and $\mathrm{300}<\mathit{\xi }<\mathrm{5000}$. Results of computations are reported graphically in Fig. 7.

The misfit between the simulated and observed flow path is due to the presence of the lava flow depicted in Fig. 6 that was not included in the DEM. Overall, the fit is reasonable in terms of maximum run-out and area covered by the modeled deposit for runs (d) $(\mathit{\mu }=\mathrm{0.3},\mathit{\xi }=\mathrm{500})$ and (f) $(\mathit{\mu }=\mathrm{0.4},\mathit{\xi }=\mathrm{5000})$. It is interesting to note that the value $\mathit{\mu }=\mathrm{0.3}-\mathrm{0.4}$ corresponds to fairly low friction coefficients of dry granular materials (corresponding to repose angles of 16–20^∘). This observation is commonly acknowledged as frictional weakening for rapid granular flows (e.g., Lucas et al., 2014).

A thorough discussion about the optimal choice of the rheological model and parameters for pyroclastic avalanches would require an extensive comparison with similar phenomena that occurred at Etna (few of which have been documented so far) and at other analogous volcanoes, for which more accurate measurements will be needed in the future to achieve a better calibration of the model. This, is however, clearly beyond the scope of the present work.

6 Conclusions

We have presented the physical formulation, numerical solution strategy and verification tests of the new IMEX_SfloW2D numerical model for shallow granular avalanches. The numerical code is available open-source and freely downloadable from a GIT repository, where the users can also find the documentation and example tests described in this paper. The main features of the new model make it suited for research and application to geophysical granular avalanches, in particular the following.

• The flexible discretization and numerical solution algorithm (not tied to knowledge of the eigenstructure of the system of equations) allows for the easy implementation of new transport equations.

• The formulation in Cartesian geographical coordinates is suited for running on digital surface models (read in standard ESRI ASCII grid format) and for the integration of non-hydrostatic terms, even on steep slopes.

• The conservative and positivity-preserving numerical scheme allows for a robust and accurate tracking of 1-D and 2-D discontinuities, including wet–dry interfaces and flow fronts.

• The implicit coupling of nonlinear rheology terms allows for the simulation of steady-state equilibrium solutions and, in particular, favors flow stopping without the need for any ad hoc empirical criteria.

• The numerical procedure to evaluate the Jacobian of the nonlinear system (based on a complex-step derivative technique) allows for an easy implementation and testing of new rheological models for complex geophysical granular avalanches.