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) and fountaining during moderately explosive eruptions. 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 (1-50 vol.%), manifesting strong topographic control. The simulation of their dynamics and mapping of their hazards pose several different problems to researchers and practitioners, 5 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 (where 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 10 are formulated to allow a robust description of wet-dry conditions (thus allowing to accurately track the front propagation) and to implicitly solve the non-linear 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 of 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 February 11th pyroclastic avalanche at Etna volcano (Italy) is finally presented. In the present formulation, 15 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 million of cubic metres) 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. It is worth remarking that this is also the volume threshold identified by Ogburn and Calder (2017) above which modelling of pyroclastic currents becomes more problematic, showing 20 transitional features between the two phenomenologies.

Modelling 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 25 volcanological research and applications (Pitman et al., 2003;Kelfoun and Druitt, 2005;Shimizu et al., 2017), we also adopt here a physical formulation based on depth-averaging, which is appropriate for shallow granular avalanches and it 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 one hand, 30 terrain-following coordinates are often used to express the depth-averaged transport equations. However, on 3D 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ẑ is nonnegligible, 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 revealed to be problematic for strongly stratified and non-homogeneous flows (Bartelt et al., 2016;Kelfoun, 2017;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). 5 In addition to the latters, some additional difficulties arise from the numerical solution of the conservation equations: despite numerical method based on conserv ative, approximate Riemann solvers are robust and well tested 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, where friction 10 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 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 15 (absolute) coordinate system, so that it is possible to include non-hydrostatic terms arising from the steep topographic slopes or from more accurate approximation of the vertical momentum equations. 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. 20 Numerically, the first and most relevant advancement 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 non-linear rheologies. The model is indeed discretized in time with an explicit-implicit Runge-Kutta method where 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- 25 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 where the elements of the Jacobian on the non-linear system are evaluated numerically with a complex step derivative technique. This automatic procedure allows to use different formulations of the friction term without the need of major modifications of the code. In particular, the Voellmy-Salm empirical model is implemented in the present version.

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

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 the Newton's second equation of motion.

Depth-averaged equations
The model we use for the fluid evolution is described by the Saint Venant's equations (Pudasaini and Hutter, 2007;Toro, 2013), 5 coupled with source terms modeling frictional forces. Such partial differential equations assume that vertical flow motion can be neglected and model the changes in space and time of flow depth h(x, y, t) and horizontal velocities u(x, y, t) and v(x, y, t) (averaged across the vertical column) over a topography B(x, y) (we assume here that the topography does not change with time t). Here we write the equations in global Cartesian coordinates, thus the two velocities are defined as the components along the x and y axes, orthogonal to the z axis, which is parallel to the gravitational acceleration g = (0, 0, g). In addition, 10 in this work we assume an 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) With these assumptions, and without considering frictional forces, the 2D inviscid depth-averaged equations in differential form can be written in the following way: The first equation represents conservation of mass, while the other two equations describe conservation of momentum in x and y directions. The last terms in Eqs. (3-4) account for the gravity-induced force, which results from the hydrostatic pressure 20 and depend also on the bathymetry B(x, y).
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. Please note that if we substitute w to h in the temporal derivative of Eq.
(2), this still represents mass conservation, because we assume B not changing with time.
If we introduce now the vector of conservative variables Q = (w, hu, hv) T (where the superscript notation Q (i) is used to 25 denote the i−th component), the governing equations can be written in the compact form: 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: S 1 (Q) = (0, ghB x , ghB y ) T .

5
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 to easily adopt finite-volumes discretization schemes for its numerical solution.

Voellmy-Salm rheology
To shallow pyroclastic avalanches, we have to modify the classic Saint Venant's equations introducing an additional source 10 term S 2 accounting for friction forces (Pudasaini and Hutter, 2007;Toro, 2013). To this purpose, we implemented in IMEX-SfloW2D, as a prototype non-linear model, the Voellmy-Salm rheology, which simulates the mixture motion as a homogeneous mass flow. This model is commonly used for avalanches and debris flows, but it fits also for volcanic gravitational flows.
The terms we consider appear only in the momentum equations (S (1) 2 = 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 Eqs. (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 20 turbulent friction, inversely proportional to the coefficient ξ. 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: With the notation introduced above, the final form of the equations modelling pyroclastic avalanches is: 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 treatments.
IMEX-SfloW2D is based on a finite-volume central-upwind scheme in space and on a implicit-explicit Runge-Kutta scheme for the discretization in time. The main purpose of the code is to run simulations on co-located 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 5 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 associated 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 10 B j,k+ 1 2 and B j+ 1 2 ,k and represented by the no-fill squares in Figure 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 to preserve steady states. 15 The finite-volume method here adopted is based on the semi-discrete central-upwind scheme introduced in Kurganov and Petrova (2007), where 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, while the term upwind is employed because, in the flux averaging, the weights depend on the local speeds of propagation at the interface.

Central-Upwind scheme
Following Kurganov and Petrova (2007), the semi-discretization in space leads to the following ordinary differential equa- 20 tions system in each cell: where 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 25 homogeneous system associated to Eq. (9) admits smooth steady-state solutions, as well as non-smooth steady-state solutions.
A good numerical method for the solution of 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:

30
This suggest to use, 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 as set of variables for the linear reconstruction 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 Γ : R 3 → R 3 for the mapping from conservative variables Q to physical variables U, and Γ −1 for the inverse mapping from physical to conservative vari-5 ables. 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: 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, for each variable, there are two reconstructed values, one from each cell at the two sides of the interface.
During the reconstruction step, a particular care should be taken in order to avoid unrealistic values of the physical variables, such as negative flow thickness or velocities too large. For this reason, in the case one of the reconstructed interface values of w 5 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 interface. We remark that the correction is applied to the derivative and thus also the reconstructed value of w at the opposite interface will be affected. For example, if w S j,k < B j,k− 1 2 , then we take the following derivative in the y-direction: 10 which gives the two reconstructions at the S and N interfaces of the (i, j) control volume: As stated above, an additional problem in the reconstruction step is that h can be very small, or even zero. Thus, when the physical variables u and v at the interfaces are computed from the conservative variables, a desingularization is applied in order to avoid the division by very small numbers and the corrected values are given by the following formulas: 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: 20 where the right-and left-sided local speeds a + j+ 1 2 ,k and a − j+ 1 2 ,k are estimated by: In a similar way, the numerical fluxes in the y-direction are given by: Following Kurganov and Petrova (2007), the source term S 1 is calculated trough a quadrature formula: 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 non-negative.

Runge-Kutta method 10
The semi-discrete system of equations (10) is solved using an Implicit-Explicit (IMEX), Diagonally-Implicit Runge-Kutta scheme (DIRK), because such schemes are well-suited for solving stiff systems of partial differential equations, and the governing equations are expected to be stiff given the strong non-linearities present in the friction terms. In addition, an implicit treatment of these terms allows for a better coupling of the equations and to properly recover the stoppage condition without the need to impose additional criteria or arbitrary thresholds. 15 The family of IMEX methods (Ascher et al., 1997) have been developed to solve stiff systems of partial differential equations written in the form: where in P are lumped all the non-stiff terms (in our case the semi-dicretized conservative fluxes F and G and the term S 1 ), while R denotes the stiff terms of the system (here represented by the friction term S 2 ). The system of equations (19) must be 20 solved for each control volume (j, k), but here, to keep the notation simpler, we omit the subscripts.
An IMEX Runge-Kutta with ν steps scheme consists of applying an implicit discretization to the stiff terms and an explicit one to the non-stiff terms, obtaining: for some l ≥ m.
Following Pareschi and Russo (2005), the IMEX scheme used in this work satisfy an additional condition, i.e. a lm = 0 5 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 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 10 limited to account for the fact that this force at maximum can stop the flow. Then, the system of nonlinear equations (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 equation (22) with the highest possible accuracy, since R(Q (j) ), accounting for the dependence of the friction force on flow variables, can be strongly non-linear. 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. 15 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:Ñ : C 3 → C 3 and compute the real-valued columns of the Jacobian at Q as where (e j ) j=1,...,3 are the standard basis vectors of R 3 , (·) denotes the imaginary part of complex numbers, and is a real 20 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).

Numerical tests and code verification
In this section we present numerical tests aimed at demonstrating the mathematical accuracy of the numerical model results

25
(the verification step, following Oberkampf and Trucano, 2002). Numerical tests are aimed at proving: 1. the capability to manage the propagation of discontinuities; 2. the potential to deal with complex and steep topographies and dry/wet interfaces; 3. the ability of the granular avalanche to stop, achieving the expected steady state. All the numerical tests presented here are also available on the wiki page of the code (https://github.com/demichie/IMEX_ SfloW2D/wiki), together with animations of the output and scripts to reproduce the results.

One dimensional test with discontinuous initial solution and topography
The first example is a 1D 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 5 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 (Figure 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 3 step IMEX Runge-Kutta scheme, where the reconstruction from the cell centers to the faces is applied to the physical variables.

One dimensional problem with dry/wet interface and friction
This example is a 1D 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" and a "dry" region. Thus, there is an additional numerical difficulty involving the capability of the numerical solver to propagate this discontinuities without creating regions with negative flow thickness. For this problem the bottom topography presents both smooth regions and a step, and it Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd-2018-224 Manuscript under review for journal Geosci. Model Dev. is defined as follows: cos 2 (πx) + 0.25(cos(10π(x − 0.5)) + 1), 0.4 ≤ x < 0.5, 0.5 cos 4 (πx) + 0.25(cos(10π(x − 0.5)) + 1), 0.5 ≤ x < 0.6, 0.5 cos 4 (πx), 0.6 ≤ x < 1, 0.25 sin(2π(x − 1)), 1 ≤ x < 1.5, 0, x ≥ 1.5. The initial conditions are defined by: For this simulation, the gravitational constant has value g = 1 and the friction coefficient is κ = 0.001(1 + 10h) −1 . Also for this test the numerical solution is obtained with a 3 step IMEX Runge-Kutta scheme, where 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, 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=0s, t=0.2s, t=2s and 5 t=40s) 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 where the horizontal gradient of w = h + B is null.

One dimensional pyroclastic avalanche with Voellmy-Salm friction
This example is a 1D 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 both preserve an initial steady condition, when the tangent of 10 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 the center of the pile coincides with the center of the domain. The gravitational constant is g = 9.81 and the parameter of the Voellmy-Salm rheology are µ = 0.3 and ξ = 300.
We present the results for a numerical simulation with a topography with a constant slope of 13°and an initial pile of material 15 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 start to move along the slope, until it reach the stoppage condition.
The numerical solution at three different times (t=0s, t=5s and t=40s) is presented in Fig. 4. For this simulation, the reconstruction technique with limiters has been applied to the physical variables and an IMEX 4-steps Runge-Kutta has been adopted, while the domain has been discretized with 400 cells. The plot of the numerical solution at the intermediate time 20 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 reach a steady condition, as shown by the plot of the velocity at 40s (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. The numerical solution at four different times (t=0s, t=7.5s, t=20s and t=25s) is presented in Fig. 5, where both the three-5 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 2-steps Runge-Kutta has been adopted, while the domain has been discretized with 150 × 100 cells.
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 of artificial numerical diffusion. As the front reach the maximum runout and    . 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 ( Figure 6), which likely controlled the avalanche path by confining it along its Southern edge (Andronico et al., 2018).

10
The avalanche rheological parameters have been varied in ranges consistent with previous studies of geophysical granular avalanches (e.g., Bartelt et al., 1999). In particular, 0.2 < µ < 0.5 and 300 < ξ < 5000 . Results of computations are reported graphically in Figure 7. The misfit between the simulated and observed flow path is due to the presence of the lava flow depicted in Figure 6 that was not included in the DEM. Overall, the fit is reasonable in terms of maximum runout and extension of the deposit for runs d) (µ = 0.3, ξ = 500) and f) (µ = 0.4, ξ = 5000).

15
A thorough discussion about the optimal choice of the rheological model and parameters for pyroclastic avalanches would  other analogous volcanoes, where 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.

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 5 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 the research and application to geophysical granular avalanches. In particular: -The flexible discretization and numerical solution algorithm (not tied to the knowledge of the eigenstructure of the system of equations) allows to easily implement new transport equations.
-The formulation in Cartesian geographical coordinates is suited for running on Digital Surface Models (read in standard 10 ESRI ASCII grid format) and for integration of non-hydrostatic terms, even on steep slopes.
-The conservative and positivity preserving numerical scheme allows a robust and accurate tracking of 1D and 2D discontinuities, including wet-dry interfaces and flow fronts.
-The implicit coupling of non-linear rheology terms allows the simulation of steady state equilibrium solutions and, in particular, favours the flow stopping without the need of any ad-hoc empirical criterion. -The numerical procedure to evaluate the Jacobian of the non-linear system (based on a complex step derivative technique) allows an easy implementation and testing of new rheological models for complex geophysical granular avalanches.
Pre-processing 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 (github.com/demichie/IMEX_SfloW2D/wiki)

5
where detailed informations on how to run the simulations are given.
Author contributions. MdMV has developed the numerical algorithm and implemented the 1D 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 2D. AA has tested the model and helped with comparison with similar depth-averaged models in volcanology.
Competing interests. No competing interests are present.