Articles | Volume 12, issue 1
Model description paper
01 Feb 2019
Model description paper |  | 01 Feb 2019

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

Mattia de' Michieli Vitturi, Tomaso Esposti Ongaro, Giacomo Lari, and Alvaro Aravena

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 Mm3) 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 Hutter2007).

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 Karman1940; Benjamin1968; Simpson1999), 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 Huppert1996; Dade2003) 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 Woods1996; Esposti Ongaro et al.2016), which may give rise to the formation of co-ignimbrite eruption columns (Woods and Wohletz1991; 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 (Kelfoun2011) and that the basal concentrated and the upper dilute turbulent layers should be coupled dynamically (Doyle et al.2010; Kelfoun2017a; 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 (Fisher1995; 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 m3) 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 Curtis2015). 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 Druitt2005; 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 (Dufek2015).

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 Iverson2004; Fischer et al.2012). Moreover, on steep slopes, where acceleration along z^ is non-negligible, the hydrostatic approximation is flawed (Denlinger and Iverson2004; 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; Kelfoun2017a; 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 Denlinger2001; Mangeney et al.2007; Forterre and Pouliquen2008; Kelfoun2011; Iverson and George2014; 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 Iverson2001; Mangeney-Castelnau et al.2003; Patra et al.2005; Christen et al.2010; Toro2013), non-hydrostatic terms arising from the vertical momentum equation (Denlinger and Iverson2004) 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 Gertisser2009; Ogburn and Calder2017). 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 (Kelfoun2017b). 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 Hutter2007; Toro2013) 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 g=(0,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:

(1) p = ρ g h .

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:


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 1List of model variables with notation and units.

Download Print Version | Download XLSX

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 Q=(w,hu,hv)T (where the superscript notation Q(i) is used to denote the ith component), the governing equations can be written in the compact form:

(5) Q t + F ( Q ) x + g ( Q ) y = S 1 ( Q ) ,

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:


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 (Toro2013). 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 S2 accounting for friction forces (Pudasaini and Hutter2007; Toro2013). 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 (S2(1)=0) and are given by


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

(8) g n = g 1 + B x 2 + B y 2 .

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

(9) Q t + F ( Q ) x + g ( Q ) y = S 1 ( Q ) + S 2 ( Q ) ,

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 Bj,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 Bj,k+12 and Bj+12,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, Bj,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 1Computational 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:


where Qj,k denotes the average of the conservative variables Q(x,y) over the control volume (j,k), and Hx and Hy 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:

(11) w = h + B = const , u = v = 0 .

This suggests using, as the vector of variables for the linear reconstruction at the interfaces, the vector 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 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 Γ:R3R3 for the mapping from conservative variables Q to physical variables U and Γ−1 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 (Ux)j,k and (Uy)j,k, the one-side values at the east, west, north and south interfaces of the cell (j,k) are given by


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 wj,kS<Bj,k-12, then we take the following derivative in the y direction:

(12) w y j , k = w ¯ j , k - B j , k - 1 2 Δ y / 2 ,

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

(13) w j , k S = B j , k - 1 2 , w j , k N = 2 w ¯ j , k - B j , k - 1 2 .

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:

(14) u = 2 h ( h u ) h 4 + max ( h 4 , ϵ ) , v = 2 h ( h v ) h 4 + max ( h 4 , ϵ ) ,

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


where the right- and left-sided local speeds aj+12,k+ and aj+12,k- are estimated by


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


where local speeds in the y direction bj,k+12+ and bj,k+12- are given by


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


This discretization, coupled with the fact that the elevation Bj,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 Petrova2007).

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

(19) Q t + P ( Q ) = R ( Q ) ,

where in P are lumped all the non-stiff terms (in our case the semi-discretized conservative fluxes F and g and the term S1), while R denotes the stiff terms of the system (here represented by the friction term S2). 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


The choice of the number of the Runge–Kutta steps ν, of the ν×ν matrices Ã=(ãlm) and A=(alm), and of the vectors b̃=(b̃1,,b̃ν) and b=(b1,,bν) differentiates the various IMEX Runge–Kutta schemes. We remark that the explicit discretization of the non-stiff terms requires ãlm=0 for lm, while the implicit treatment of the stiff terms requires alm≠0 for some lm.

Following Pareschi and Russo (2005), the IMEX scheme used in this work satisfies an additional condition, i.e., alm=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:


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 Ñ:C3C3, and compute the real-valued columns of the Jacobian at Q as

(23) J e j = ( N ̃ ( Q + i ϵ e j ) ) ϵ ,

where (ej)j=1,,3 represents the standard basis vectors of 3, ℑ(⋅) 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 21-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 Trucano2002). 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 (, 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 Petrova2007). 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:


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


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 3Numerical 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:


The initial conditions are defined by


For this simulation, the gravitational constant has a value g=1 and a simpler friction term is employed (κ(h)u, with κ(h)=0.001(1+10h)-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 [-0.25;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 4Numerical 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−2 and the parameters of the Voellmy–Salm rheology are μ=0.3 and ξ=300 ms−2.

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 5Numerical 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 [0;30]×[-7;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 0x17.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−2.

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 6Map 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 7Results 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×106 m3.


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×106 m3. 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, 0.2<μ<0.5 and 300<ξ<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) (μ=0.3,ξ=500) and (f) (μ=0.4,ξ=5000). It is interesting to note that the value μ=0.3-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.

Code availability

The numerical code, benchmark tests and documentation are available at (last access: 30 January 2019). Preprocessing scripts (to change the grid resolution and the numerical schemes) and post-processing scripts (to plot the solution variables and to create animations) are also available. Furthermore, each example has a page description on the model Wiki (, last access: 30 January 2019), where detailed information on how to run the simulations is given. The digital object identifier (DOI) for the version of the code documented in this paper is (de' Michieli Vitturi and Lari2019).

Author contributions

MdMV has developed the numerical algorithm and implemented the 1-D model and numerical tests. TEO has contributed to the model formulation and application in the context of volcanological applications. GL has implemented the numerical scheme and algorithm in 2-D. AA has tested the model and helped with the comparison to similar depth-averaged models in volcanology.

Competing interests

The authors declare that they have no conflict of interest.


The work has been supported by the Italian Department of Civil Protection, INGV-DPC agreement B2 2016, task D1. We warmly thank Daniele Andronico, Emanuela De Beni and Boris Behncke for useful discussions about the Etna 2014 event and for providing field data and the digital elevation model. TEO would like to thank Perry Bartelt and Betty Sovilla (SLF, Switzerland) for stimulating discussions on granular avalanches and density currents. The authors are also grateful to Olivier Roche and Karim Kelfoun for their careful review and useful and constructive comments. Mattia de' Michieli Vitturi and Giacomo Lari gratefully acknowledge Andrea Milani, whose contribution to the career and numerical work of the two authors was of great significance.

Edited by: Simone Marras
Reviewed by: Karim Kelfoun and Olivier Roche


Andrianov, N.: Testing numerical schemes for the shallow water equations, Tech. rep., available at: (last access: 30 January 2019), 2004. a

Andronico, D., Di Roberto, A., De Beni, E., Behncke, B., Antonella, B., Del Carlo, P., and Pompilio, M.: Pyroclastic density currents at Etna volcano, Italy: The 11 February 2014 case study, J. Volcanol. Geoth. Res., 357, 92–105,, 2018. a, b

Ascher, U. M., Ruuth, S. J., and Spiteri, R. J.: Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations, Appl. Numer. Math., 25, 151–167, 1997. a

Bartelt, P., Salm, L. B., and Gruberl, U.: Calculating dense-snow avalanche runout using a Voellmyfluid model with active/passive longitudinal straining, J. Glaciol., 45, 242–254,, 1999. a, b

Bartelt, P., Buser, O., Valero, C. V., and Bühler, Y.: Configurational energy and the formation of mixed flowing/powder snow and ice avalanches, Ann. Glaciol., 57, 179–188,, 2016. a

Benjamin, T. B.: Gravity currents and related phenomena, J. Fluid Mech., 31, 209–248, 1968. a

Bursik, M. I. and Woods, A. W.: The dynamics and thermodynamics of large ash flows, B. Volcanol., 58, 175–193,, 1996. a

Castro-Orgaz, O., Hutter, K., Giraldez, J. V., and Hager, W. H.: Nonhydrostatic granular flow over 3-D terrain: New Boussinesq-type gravity waves?, J. Geopys. Res., 120, 1–28,, 2015. a

Charbonnier, S. J. and Gertisser, R.: Numerical simulations of block-and-ash flows using the Titan2D flow model: examples from the 2006 eruption of Merapi Volcano, Java, Indonesia, B. Volcanol., 71, 953–959,, 2009. a

Christen, M., Kowalski, J., and Bartelt, P.: RAMMS: Numerical simulation of dense snow avalanches in three-dimensional terrain, Cold Reg. Sci. Technol., 63, 1–14, 2010. a

Dade, W. B.: The emplacement of low-aspect ratio ignimbrites by turbulent parent flows, J. Geophys. Res., 108, 1–9,, 2003. a

Dade, W. B. and Huppert, H. E.: Emplacement of the Taupo ignimbrite by a dilute turbulent flow, Nature, 381, 509–512, 1996. a

De Beni, E., Behncke, B., Branca, S., Nicolosi, I., Carluccio, R., D'Ajello Caracciolo, F., and Chiappini, M.: The continuing story of Etna's New Southeast Crater (2012–2014): Evolution and volume calculations based on field surveys and aerophotogrammetry, J. Volcanol. Geoth. Res., 303, 175–186, 2015. a

Delannay, R., Valance, A., Mangeney, A., Roche, O., and Richard, P.: Granular and particle-laden flows: from laboratory experiments to field observations, J. Phys. D, 50, 1–40,, 2017. a

de' Michieli Vitturi, M. and Lari, G.: IMEX_SfloW2D v. 1.0.0 (Version 1.0.0), Zenodo,, 2019. a

Denlinger, R. P. and Iverson, R. M.: Flow of variably fluidized granular masses across three-dimensional terrain: 2. Numerical predictions and experimental tests, J. Geophys. Res., 106, 553–566,, 2001. a

Denlinger, R. P. and Iverson, R. M.: Granular avalanches across irregular three-dimensional terrain: 1. Theory and computation, J. Geophys. Res., 109, F01014,, 2004. a, b, c

Doyle, E. E., Hogg, A. J., Mader, H. M., and Sparks, R. S. J.: A two-layer model for the evolution and propagation of dense and dilute regions of pyroclastic currents, J. Volcanol. Geoth. Res., 190, 365–378,, 2010. a

Dufek, J.: The Fluid Mechanics of Pyroclastic Density Currents, Annu. Rev. Fluid Mech., 48, 459–485,, 2015. a

Engwell, S., de' Michieli Vitturi, M., Esposti Ongaro, T., and Neri, A.: Insights into the formation and dynamics of coignimbrite plumes from one-dimensional models, J. Geophys. Res.-Sol. Ea., 121, 4211–4231, 2016. a

Esposti Ongaro, T., Orsucci, S., and Cornolti, F.: A fast, calibrated model for pyroclastic density currents kinematics and hazard, J. Volcanol. Geoth. Res., 327, 257–272,, 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,, 2012. a

Fisher, R. V.: Decoupling of pyroclastic currents: hazards assessments, J. Volcanol. Geoth. Res., 66, 257–263,, 1995. a

Fisher, R. V., Orsi, G., Ort, M., and Heiken, G.: Mobility of a large-volume pyroclastic flow–emplacement of the Campanian ignimbrite, Italy, J. Volcanol. Geoth. Res., 56, 205–220,, 1993. a

Forterre, Y. and Pouliquen, O.: Flows of Dense Granular Media, Annu. Rev. Fluid Mech., 40, 1–24,, 2008. a

Guo, Y. and Curtis, J. S.: Discrete Element Method Simulations for Complex Granular Flows, Annu. Rev. Fluid Mech., 47, 21–46,, 2015. a

Iverson, R. M. and Denlinger, R. P.: Flow of variably fluidized granular masses across three-dimensional terrain: 1. Coulomb mixture theory, J. Geophys. Res.-Sol. Ea., 106, 537–552,, 2001. a

Iverson, R. M. and George, D. L.: A depth-averaged debris-flow model that includes the effects of evolving dilatancy. I. Physical basis, P. R. Soc. A, 470, 20130819,, 2014. a

Kelfoun, K.: Suitability of simple rheological laws for the numerical simulation of dense pyroclastic flows and long-runout volcanic avalanches, J. Geophys. Res., 116, B08209,, 2011. a, b

Kelfoun, K.: A two-layer depth-averaged model for both the dilute and the concentrated parts of pyroclastic currents, J. Geophys. Res.-Sol. Ea., 122, 4293–4311,, 2017a. a, b

Kelfoun, K.: A two-layer depth-averaged model for both the dilute and the concentrated parts of pyroclastic currents, J. Geophys. Res.-Sol. Ea., 122, 4293–4311, 2017b. a

Kelfoun, K. and Druitt, T. H.: Numerical modeling of the emplacement of Socompa rock avalanche, Chile, J. Geophys. Res.-Sol. Ea., 110, B12202,, 2005. a

Komorowski, J.-C., Jenkins, S., Baxter, P. J., Picquout, A., Lavigne, F., Charbonnier, S., Gertisser, R., Preece, K., Cholik, N., Budi-Santoso, A., and Surono: Paroxysmal dome explosion during the Merapi 2010 eruption: Processes and facies relationships of associated high-energy pyroclastic density currents, J. Volcanol. Geoth. Res., 261, 260–294,, 2013. a

Kurganov, A. and Petrova, G.: A second-order well-balanced positivity preserving central-upwind scheme for the Saint-Venant system, Commun. Math. Sci., 5, 133–160, 2007. a, b, c, d, e, f

La Spina, G. and de' Michieli Vitturi, M.: High-resolution finite volume central schemes for a compressible two-phase model, SIAM J. Sci. Comput., 34, B861–B880, 2012. a

Lucas, A., Mangeney, A., and Ampuero, J. P.: Frictional velocity-weakening in landslides on Earth and on other planetary bodies, Nat. Commun., 5, 1–9,, 2014. a, b

Mangeney, A., Bouchut, F., Thomas, N., Vilotte, J. P., and Bristeau, M. O.: Numerical modeling of self-channeling granular flows and of their levee-channel deposits, J. Geophys. Res., 112, F02017,, 2007. a

Mangeney-Castelnau, A., Vilotte, J.-P., Bristeau, M.-O., Perthame, B., Bouchut, F., Simeoni, C., and Yerneni, S.: Numerical modeling of avalanches based on Saint Venant equations using a kinetic scheme, J. Geophys. Res.-Sol. Ea., 108, 2527,, 2003. a

Martins, J. R. R. A., Sturdza, P., and Alonso, J. J.: The complex-step derivative approximation, ACM T. Math. Software, 29, 245–262,, 2003. a

Oberkampf, W. L. and Trucano, T. G.: Verification and validation in computational fluid dynamics, Prog. Aerosp. Sci., 38, 209–272, 2002. a

Ogburn, S. E. and Calder, E. S.: The Relative Effectiveness of Empirical and Physical Models for Simulating the Dense Undercurrent of Pyroclastic Flows under Different Emplacement Conditions, Front. Earth Sci., 5, 63–26,, 2017.  a, b

Pareschi, L. and Russo, G.: Implicit–explicit Runge–Kutta schemes and applications to hyperbolic systems with relaxation, J. Sci. Comput., 25, 129–155, 2005. a

Patra, A. K., Bauer, A., Nichita, C., Pitman, E. B., Sheridan, M., Bursik, M., Rupp, B., Webber, A., Stinton, A., Namikawa, L., and Renschler, C. S.: Parallel adaptive numerical simulation of dry avalanches over natural terrain, J. Volcanol. Geoth. Res., 139, 1–21, 2005. a, b

Pitman, E. B., Nichita, C. C., Patra, A., Bauer, A., Sheridan, M., and Bursik, M.: Computing granular avalanches and landslides, Phys. Fluids, 15, 3638–3646,, 2003. a, b

Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P.: Numerical recipes in Fortran 90, vol. 2, Cambridge university press Cambridge, 1996. a

Pudasaini, S. P. and Hutter, K.: Avalanche dynamics: dynamics of rapid flows of dense granular avalanches, Springer Science & Business Media, 2007. a, b, c

Shea, T., Gurioli, L., Houghton, B. F., Cioni, R., and Cashman, K. V.: Column collapse and generation of pyroclastic density currents during the A.D. 79 eruption of Vesuvius: The role of pyroclast density, Geology, 39, 695–698,, 2011. a

Shimizu, H. A., Koyaguchi, T., and Suzuki, Y. J.: A numerical shallow-water model for gravity currents for a wide range of density differences, Progress in Earth and Planetary Science, 4, 1–13,, 2017. a, b, c

Simpson, J. E.: Gravity currents: In the environment and the laboratory, Cambridge University Press, 1999. a

Toro, E. F.: Riemann solvers and numerical methods for fluid dynamics: a practical introduction, Springer Science & Business Media, 2013. a, b, c, d

Von Karman, T.: The engineer grapples with nonlinear problems, B. Am. Math. Soc., 46, 615–683, 1940. a

Wang, Y., Hutter, K., and Pudasaini, S. P.: The Savage-Hutter theory: A system of partial differential equations for avalanche flows of snow, debris, and mud, ZAMM-Z. Angew. Math. Me., 84, 507–527, 2004. a

Woods, A. W. and Wohletz, K.: Dimensions and dynamics of co-ignimbrite eruption columns, Nature, 350, 225–227,, 1991. a

Woods, A. W., Bursik, M. I., and Kurbatov, A. V.: The interaction of ash flows with ridges, B. Volcanol., 60, 38–51, 1998. a

Yuan, L., Liu, W., Zhai, J., Wu, S. F., Patra, A. K., and Pitman, E. B.: Refinement on non-hydrostatic shallow granular flow model in a global Cartesian coordinate system, Comput. Geosci., 109, F01014,, 2017. a

Short summary
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) and fountaining during moderately explosive eruptions. We present IMEX_SfloW2D, a depth-averaged flow model describing the granular mixture as a single-phase granular fluid. Benchmark cases and preliminary application to the simulation of the 11 February pyroclastic avalanche at Mt. Etna (Italy) are shown.