## Abstract

Prediction of the emplacement of volcanic mass flows (lava flows, pyroclastic density currents, debris avalanches and debris flows) is required for hazard and risk assessment, and for the planning of risk-mitigation measures. Numerical computer-based models now exist that are capable of approximating the motion of a given volume of volcanic material from its source to the deposition area. With these advances in technology, it is useful to compare the various codes in order to evaluate their respective suitability for real-time forecasting, risk preparedness and post-eruptive response. A ‘benchmark’ compares codes or methods, all aimed at simulating the same physical process using common initial and boundary conditions and outputs, but using different physical formulations, mathematical approaches and numerical techniques. We set up the basis for a future general benchmarking exercise on volcanic mass-flow models and, more specifically, establish a benchmark series for computational lava-flow modelling. We describe a set of benchmarks in this paper, and present a few sample results to demonstrate output analysis and code evaluation methodologies. The associated web-based communal facility for sharing test scenarios and results is also described.

## Motivation and background

The on-going demographic congestion around volcanic edifices increases the potential risks and costs that volcanic flows represent, and leads to a growing demand for faster and more accurate predictions of flow extent, both spatially and temporally. As a result, various computer codes have been developed to predict lava-flow footprints and emplacement dynamics. Both physical and numerical implementations differ between models, depending on the aim, whether for the generation of hazard maps (Favalli *et al.* 2012; Tarquini & Favalli 2013), the forecasting of lava-flow emplacement (Hidaka *et al.* 2005; Vicari *et al.* 2007; Proietti *et al.* 2009; Ganci *et al.* 2012) or the evaluation of the potential of human counter-measures to deviate lava flows (e.g. channels generation from explosives, spraying water, building levees: Fujita *et al.* 2009).

Lava flows are a constant threat for surrounding infrastructures and populations (Komorowski *et al.* 2002; Behncke & Neri 2003). There is a clear need to relate eruption conditions and the resulting flow extent. Walker (1973) first established the importance of the effusion rate (*Q*) in the determination of the lava final length (*l*_{f}), a relationship still used for Etna's lava flows (Calvari & Pinkerton 1998) but which does not hold for thermally-insulated systems such as lava tubes, as is evident from the work of Malin (1980) and Pinkerton & Wilson (1994). In these conditions, the lower cooling rate allows the flow to reach a further extent, as has been quantified through analogue experiments and scaling arguments (e.g. Griffiths *et al.* 2003). Finally, the viscous and plastic characteristics of lavas have a major impact on the spatial extent and final shape of lava flows (Walker 1967; Hulme 1974). For example, plastic yield stress determines a minimum thickness of the lava before flow can begin and partly determines the height of the emplaced lava flow: hence, yield stress has been widely used in numerical codes as a flow stopping condition.

At a first-order level, lava flows can be described as gravity currents, spreading under their own weight, slowed down by their internal viscous resistance or by the topography. For a given isoviscous Newtonian lava flow, the mean velocity is (Jeffreys 1925):
(1)
which connects the main driving parameters – the slope (α), flow height (*h*), density (ρ) and gravity acceleration (*g*) – with the resisting forces due to the shear viscosity (η) and the parameter *S*, a shape factor reflecting the effect of the flow cross-section geometry. For example, *S*=3 for an equant rectangular cross-section and decreases when the width becomes much greater than the depth (Sakimoto & Gregg 2001; Lev & James 2014). One of the first uses of Jeffreys equation on lavas occurs in Nichols (1939), and the first thorough description of lava-flow movement and spreading features was made during the 1943–52 Paricutin eruption (Krauskopf 1948). Walker (1973) then first established the importance of the effusion rate (*Q*). Additional equations take into account the lava yield stress and consider the basal profile section below this critical stress as being non-flowing (Moore 1987; Ishihara *et al.* 1990). Lava cooling leads to an increase in magma viscosity, and is controlled by lava thermal conductivity (κ), specific heat (*C*_{p}) and density (ρ). Cooling may eventually be slowed down through latent heat supplied by crystallization. Mechanical stress applied to the system can also lead to a decrease in viscosity through viscous heating and shear thinning (Costa & Macedonio 2003; Cordonnier *et al.* 2009, 2012).

A numerical code for lava emplacement must consider:

topography or slope;

eruptive input conditions – volume effusion rate, vent geometry (point or linear) and effusive temperature;

thermal boundary conditions at the top and bottom of the lava – insulation, convection, radiation and conduction;

physical properties of the lava – density, thermal conductivity, rheological laws (i.e. how viscosity varies as a function of stress, temperature and bubble/crystal content).

None of these properties impacts the lava emplacement in the same way; the magnitude of their effect varies with distance from the vent but also with time (Hulme 1974; Chester *et al.* 1985). Lava flows have several different characteristic flow morphologies, such as Pahoehoe, a'a’, channels, tilted blocks or lava-tube flows (see Fig. 1). Each of these morphologies is associated with a specific cooling rate function of the efficiency of heat loss through radiation, conduction or convection. This cooling rate determines a vertical thermal profile directly affecting the lava viscosity and velocity field, and potentially leading to a new flow morphology. Some magma properties, such as rheology, are mostly known from experimental measurement rather than *in situ* measurements. Numerical codes and theoretical models cannot take the entire complexity of lava properties into account, and simplifications are required to speed up codes. Hence, the need to evaluate the accuracy and robustness of the model results.

Numerical simulations of lava flows have undergone major developments over the past years (e.g. coupling flow and thermal equations, going from 2D to 3D). These steps were achieved through various methods and solvers (see the following section on ‘Modelling lava-flow emplacement’), each with their own limits and approximations. Codes differ in their physical implementations, numerical accuracy and computational efficiency. Therefore, a common testing ground is essential, so that: (i) developers can test their code efficiency and accuracy under conditions of varying complexity; and (ii) the target audience interested in lava-flow model output, such as civil defence authorities, can correctly evaluate the reliability of model predictions and the applicability of each code (early eruptive phase, thermal cooling of emplaced systems, counter-measures and hazard maps). A meaningful comparison of multiple numerical codes needs to be based on well-defined benchmarks.

A benchmark compares the outputs of models that simulate the same physical process (the same initial and boundary conditions, and the same physical flow properties). Such models may have different physical implementations, mathematical approaches or numerical techniques but the benchmark yields an identification of efficient and accurate strategies. To achieve such an ambitious goal, the complexity of the benchmarks must increase gradually (e.g. from an established initial constant volume to time-dependent flow rates, from isothermal conditions to cooling flows, from Newtonian to non-Newtonian elasto-visco-plastic flows). Our goal is to trigger a community-based effort to set-up such benchmarks, encouraging model developers and those who use the results to share their need and expertise:

to define common calibration and validation test cases for numerical models of lava-flow advance;

to better identify the dominant control parameters on lava-flow advance and velocity;

to evaluate the reliability of model predictions during a volcanic crisis.

## Modelling lava-flow emplacement

Historically, most models of lava-flow emplacement used empirical laws based on natural observations (Yokoyama 1957; Hedervari 1963; Friedman & Williams 1968; Scandone 1979), and attempted to link mass and heat produced. Later similar approaches linked heat flux to effusion rate or mass flux from eruptive events (Pieri & Baloga 1986; Crisp & Baloga 1990; Francis *et al.* 1993; Harris & Stevenson 1997; Wooster *et al.* 1997; Harris *et al.* 1998) or were constrained experimentally (Garel *et al.* 2012). Some work also improved the pioneering studies of Walker (1973) by linking the area covered by lava to the discharge rate (Dragoni & Tallarico 2009; Harris & Baloga 2009; Harris *et al.* 2010), again with experimental support (Griffiths & Fink 1992, 1993).

Existing flow-emplacement models can be divided into different categories, depending on the approach (deterministic v. stochastic/probabilistic models), on the numerical method and on the complexity of the physical modelling. Physical processes being modelled by flow simulation include the motion of the lava (conserving momentum and mass), the changes in rheology due to various factors (e.g. temperature, strain-rate, bubble and/or crystal content) and the thermal evolution of the flow. Conditions defining the simulation set-up include topography data and effusion rate (mass input). We shortly present hereafter the main characteristics of the existing model types (for further details, see the lava-flow modelling section of this volume).

### Existing numerical models

In this subsection, we provide a short overview of the available flow emplacement codes, divided into subdivisions based on the numerical approach. We first discuss deterministic methods, in which a solution is achieved given a unique set of boundaries and initial conditions. We then look at the stochastic methods, the solutions of which result from up to thousands of sets of input parameters (see the subsection on ‘Stochastic models’ later). It is noteworthy that most codes can be run independently in a deterministic or stochastic mode depending on the specific hazard targeted. Table 1 presents a summary of the current existing models simulating lava-flow emplacement, with associated pros and cons. For all approaches, the DEM of the volcano topography and the vent location is required to simulate lava-flow emplacement. Deterministic models also require the temporal evolution of the mass or volume effusion rate at the eruptive vent, as well as any thermal or rheological model required by the model set-up.

#### Channelled models

For a channelled, confined liquid, flow advance occurs downslope in only one direction, as opposed to free flows that can also spread across the slope. Models of channelled lava dynamics are thus simplified regarding the *emplacement* physical processes, and can implement more complex thermal modelling (e.g. crystallization rate).

Velocity profiles for lava flow down a channel depend on the channel dimensions and on the lava rheology model (Jeffreys 1925; Dragoni 1989; Tallarico & Dragoni 2000). Flow can also feed back into the rheology for a strain-rate dependent viscosity (Robertson & Kerr 2012*a*, *b*).

The main code for one-dimensional (1D) channel lava emplacement is FLOWGO (e.g. Harris & Rowland 2001; Harris *et al.* 2011), which models the motion of finite-mass lava parcels in a channel confined by levees (see also Harris *et al.* 2015). FLOWGO models the heat lost by radiation and convection at the surface of the flow, and by conduction through the base. The resulting cooling rate sets a crystallization rate, and rheology is calculated from its crystal- and temperature-dependency. The flow is a 1D line, and stops once the flow has cooled down to rheological properties that equate to a solid.

Confined flow models are very fast to run as they do not involve calculating the fluid motion. However, important limitations arise from the assumption of one-dimensionality forcing the local channel width to directly correlate with ground slope. The choice of the rheological model is also crucial to link thermal evolution and flow velocity. Moreover, in order to match the effusion rate input, the FLOWGO model needs to assume at-vent channel dimensions, which greatly influence the model flow results (Wantim *et al.* 2013).

#### Cellular automata

The cellular automata (CA) approach is one of the most popular used to model lava-flow emplacement, with codes such as MAGFLOW (Vicari *et al.* 2007, 2009; Del Negro *et al.* 2008; Herault *et al.* 2009; Bilotta *et al.* 2012; Ganci *et al.* 2012), SCIARA (Crisci *et al.* 1986, 1998, 2008; Barca *et al.* 1994, 1999, 2003, 2004; Barca *et al.* 2004; Avolio *et al.* 2006) and FLOWFRONT (Young & Wadge 1990; Wadge *et al.* 1994). In CA, the computational domain is represented by a 2D grid of cells. Each cell is characterized by properties such as altitude, lava height and temperature, and lava-flow advance and cooling are described through the evolution of cell properties. This evolution depends on the state of the neighbouring cells, with in/out fluxes computed for each cell at each time step. For example, the flow of lava from one cell to another depends on the ‘slope’ (altitude difference between two adjacent cells) and on the flow density, and the resulting new lava temperature will reflect the different in/out fluxes. At first, CA codes used plastic rheology to stop the flow (Ishihara *et al.* 1990). Today, solidification takes place when the cell temperature gets below a predefined temperature (Gregorio & Serra 1999). Cellular automata allow for the modelling of compound lava flows and more generalized lava-flow fields as cell altitude is updated by adding the lava thickness if solidification occurred.

Cellular automata models are fast to run, and each of them has developed specific strategies to transfer mass, energy and momentum between neighbouring cells: see Rongo *et al.* this volume, in review; Cappello *et al.* this volume, in press for the examples of the SCIARA and MAGFLOW codes. However, 2D models will always lack the detailed vertical structure of the lava flow, which is important in coupling surface and basal heat losses to the bulk rheology evolution.

#### Depth-averaged models

Depth-averaged methods use the shallow-water equations (SWE) hypothesis introduced by Barre de Saint-Venant (1871), which assumes that the horizontal length scale is much greater than the vertical one and neglects the vertical component (i.e. homogeneous flow properties are assumed throughout the vertical section). Newtonian, but also some more complex rheologies such as Bingham (Liu & Mei 1990) and Herschel–Buckley, can be implemented (Balmforth *et al.* 2002, 2006, 2007).

Depth-averaged codes for volcanological applications include the model developed by G. Macedonio (Costa 2005; Macedonio *et al.* 2005), the VOLCFLOW code (Kelfoun & Druitt 2005) detailed in Kelfoun & Vallejo Vargas (2015) and the implementation made in RHEOLEF (Bernabeu *et al.* 2013; Saramito *et al.* 2013). The depth-average formulation leads to a dimension reduction as variations with depth are neglected, and are thus relatively fast to compute. Drawbacks include the necessity to assume a constant viscosity value along the vertical profile, which limits a detailed feedback of rheology into flow dynamics.

#### Nuclear-based models

An important subset of numerical flows formulating spreading under thermal cooling conditions originates from the field of nuclear reactor safety. Since the Three Mile Island and Chernobyl accidents, nuclear safety groups have established the necessity to control any potential accident in the reactor containment in order to minimize off-site consequences. The most urgent task during a severe nuclear reactor meltdown is to rapidly spread the fusion melt product (corium) over the widest possible surface in order to quench it as rapidly as possible. Thus, this shallow fluid process shares many similarities with lava flows, and simulation codes can be used by both communities.

Many experimental platforms have been employed over the years to study corium spreading. Several computer codes have been jointly developed to confirm the physical processes engaged in corium spreading: CORFLOW (Wittmaack 1997), CROCO (Michel *et al.* 2000), LAVA (Allelein *et al.* 2000), MELTSPREAD (Farmer *et al.* 1990), RASPLAV (Bolshov *et al.* 1995), SAMPSON (Hidaka *et al.* 2002) and THEMA (Spindler & Veteau 2006). A scaling-based model has been developed by Dinh *et al.* (2000). These codes are an evolution of the CA method by discretizing the height of each cell into several vertical cells, yet none models the true free surface.

Addition of a slope or topography is the only additional feature to implement to switch these nuclear-engineering codes to lava-flow emplacement models. The mass, momentum and energy conservation are resolved for this 3D system, which allows the development of complex emplacement features related to solidification, such as lava lobes or tubes. The code LavaSIM (Hidaka *et al.* 2002, 2005; Proietti *et al.* 2009), for example, has been successfully adapted to simulate lava-flow advance (see Fujita & Nagai, this volume, in press).

#### Generic 3D computational fluid dynamics codes

Many fields, in addition to volcanology and nuclear engineering, including civil and mechanical engineering, metallurgy and biophysics share the requirement to predict flow advance. The field of computational fluid dynamics (CFD) has evolved dramatically over the past few decades, and nowadays many advanced CFD tools are available for the community. Tools range from commercial software packages to community-driven or government- supported open-source libraries. Adopting such libraries and packages to use for lava-flow models can require additional capabilities that are not always built-in, such as crystallization effects, complex topography, variable effusion rates, complex rheology including both viscous and brittle components, and the strong temperature dependence of physical properties. Nonetheless, we discuss here two advanced CFD packages, one open-source and one commercial, applicable to lava-flow modelling.

OpenFOAM (Open Field Operation And Manipulation) toolbox (http://www.openfoam.com, Jasak

*et al.*2007) is a finite-volume-method-based open-source software package produced by OpenCFD Ltd. OpenFOAM solvers can handle complex fluids, chemical reactions, turbulence, heat transfer, solid mechanics and electromagnetics. The code is fully parallelized using OpenMPI, and has straightforward interfaces with external meshing, and pre- and post-processing tools. OpenFOAM can be used as a standard simulation package with its existing capabilities, and users can easily add new equations, solvers and applications (e.g. solidification of lava can be implemented). Adaptive mesh refinement is a standard utility in OpenFOAM, and is convenient to model moving flows and interfaces, which can arise in lava-flow simulations.Flow3D is a commercial software aimed largely at engineers, with particular emphasis on computational fluid dynamics. Flow3D is produced and distributed by Flow Science Inc. At its core, Flow3D is based on Volume-of-Fluid algorithms, combined with interface tracking tools (level set). It is therefore well suited for the free surface flow that takes place during lava emplacement. All manner of heat transfer can be simulated by Flow3D, as well as porous, two-phase and viscous flows. Flow3D is a popular choice for a range of industrial applications, and the user group is large and diverse. The two main downsides for Flow3D are its slowness and high price. In exchange, the users get a fully developed and tested modelling environment, with a graphic user interface and product support.

#### Meshless and ‘bottom-up’ methods

All the above codes are based on a mesh-dependent description of the flow, which can be costly in terms of computational resources due mainly to re-meshing procedures. In addition, they are commonly so-called ‘top-down’ methods that describe macroscopic systems discretized with partial differential equations where certain quantities may not be conserved. Hence, meshless methods hold a great potential for fluid dynamics modelling. Brittle material can also be incorporated into such models more easily. They commonly form ‘bottom-up’ schemes based on a discrete microscopic model that conserves the desired quantities through its construction. However, ‘bottom-up’ methods have remained poorly used in Geosciences. The two major approaches are the Lattice Boltzmann methods (LBM: Parmigiani *et al.* 2011, 2013) and smoothed particles hydrodynamics (SPH: Herault *et al.* 2011): only the latter has been so far been applied to lava-flow dynamics. For each of these methods, additional terms such as the Boltzmann discretization for the LBM or kernel size limitation for the SPH inherently generate numerical diffusion, and require careful benchmarking. Smoothed particle hydrodynamics, meshless methods with a bottom-up scheme, are based on individual flowing particles, the characteristics of which (e.g. velocity and temperature) are defined by their interaction with the surrounding particles. This interaction is weighted, based on the particle distance, and is considered null beyond a given threshold in order to increase code efficiency. SPH are naturally 3D and can be used to simulate volcanic mass flows. One example of SPH simulations is presented in Herault *et al.* (2011) and in Bilotta *et al.* this volume, in review. We also present in the section ‘First benchmark results’ additional results from the code NB3D. This code, currently under development by B. Cordonnier, is an SPH algorithm with similar physical grounds to the ones described in Herault *et al.* (2011) and implemented within a MATLAB® structure (results are presented in the ‘First benchmark results’ section). For this paper, the code has been set to solve 1 326 240 particles under a Quintic Wendland type kernel. The time-step algorithm is symplectic and is associated every 10 time steps with a zero-order density filter to reduce oscillations of the particles pressure field (Shepard type).

### Stochastic models

In a stochastic process, the system's state is not described by unique values but, rather, by probability distributions. While, in principle, most numerical codes could be used as stochastic models, in the sense that they can be run several times with slightly changed parameters, only codes with strongly simplified physics are commonly used this way due to the computing time required. The most simplified stochastic models for lava-flow emplacement are reduced to gravity currents that flow towards the steepest slope (Dobran 1992). These models are only determined by the digital elevation model (DEM), and must account for potential errors between the available DEM and the real topography, due, for example, to acquisition uncertainties, vegetation changes, erosion or lava-flow emplacement since the last DEM acquisition. By slightly changing the elevation model (randomly or based on previous runs), the stochastic models establish the probability distribution for the lava-flow path, and which route is the most likely. For the most simplified version, stochastic models are not bound by criteria able to halt the flow advance. They are unable to present a time-evolution or downslope limit of the flow but may be used to generate hazard flow maps and to evaluate the impact of counter-measures to deviate lava flows. In more evolved versions, such as CA-based stochastic models, time evolution of the lava front is estimated.

Examples of probabilistic models for lava-flow emplacement include DOWNFLOW (e.g. Favalli *et al.* 2005; Tarquini & Favalli 2011), ELFM (Damiani *et al.* 2006), LAZSLO (Bonne *et al.* 2008) and recent models for random pahoehoe lava-lobe emplacement (Glaze & Baloga 2013; Hamilton *et al.* 2013). The DOWNFLOW code is presented in Tarquini & Favalli, this volume, in press, and represents the most efficient probabilistic lava flow model to date. It requires a very short computation time and can predict final lava-flow emplacement by taking into account probabilistic lava-flow length (Tarquini & Favalli 2013).

### Trade-off in code selection

When selecting a simulation code, accuracy is not the only consideration; speed and efficiency are also important. There is often a trade-off between speed and physical completeness. The more physical processes that a model takes into account, the longer it will take to produce a solution. Hence, if models need to solve mass, momentum and energy conservation equations, the computing time may greatly differ depending on the physical modelling, and on the vertical and lateral resolution.

Naturally, a code solving most of the active physical processes (e.g. all styles of heat loss, coupled with rheological changes) will be complex and will take a significant time to produce a result. Hence, it will not aim for a stochastic approach. In contrast, a stochastic model, similar to the Monte Carlo methods, calculates a considerable number of solutions with slightly modified input parameters and simplified flows. Stochastic codes are optimized for fast calculation and commonly only consider a few major parameters. They are capable of quickly adapting to incorporate the latest observations during an eruptive crisis such as DEM evolution, eruptive flow rate or adjustment of physical properties of lavas. However, they may offer a limited insight into the dynamics of the flow and the physical processes taking place within it. Deterministic codes, given their relative slowness, are more commonly used during the preparation and assessment phase or the post-eruptive period of volcanic hazards.

Hence, there is a balance to be found between a reasonable computation time that can be used in near-real time during an eruptive crisis (and updated depending on the observed lava-flow advance) and the accuracy of the flow dynamics predicted by the model. Note also that stochastic models can be run using a low-resolution DEM (<30 m) and are thus a reasonable choice for poorly monitored volcanic areas. This is why benchmarks are needed in order to estimate the minimum complexity of physical and thermal modelling required for models of lava-flow emplacement.

## Benchmarking fundamentals

A benchmark couples a numerical code with a specific solution, which may be an analytical solution, an experimental measurement or a natural observation. Analytical solutions are, by definition, exact and offer a unique result that numerical codes must aim to approximate most closely. As such, analytical solutions test the correctness of solvers and computation algorithms. They do not, however, reflect the complexity that real situations may exhibit. Experimental measurements offer well-constrained initial and boundary conditions, known material properties and accurate measurement of flow evolution. Experimental results are arguably the best choice for testing code accuracy, as most discrepancies are not related to uncertainties in input parameters. For example, extensive experiments were performed in the field of nuclear energy research, modelling reactor failure and corium flow. Some of these experimental results are available as benchmark data sets (Greene *et al.* 1988; Suzuki *et al.* 1993; Fieg *et al.* 1996; Cognet *et al.* 2001; Piluso *et al.* 2001; Journeau *et al.* 2003) and are used to confirm that the physics implemented in the numerical simulations is sufficient. Finally, natural cases are the ultimate target of the numerical simulations as they contain all the factors at play. However, they present many unknowns and measurement errors, which makes it difficult to determine whether the source of errors between computed solutions and natural observations is numerical or linked to data acquisition.

As discussed earlier, lava-flow behaviour reflects the interplay of various physical processes, the relative importance of which depends on the emplacement conditions and flow regimes. All the simulation codes model only a subset of these parameters and usually apply them with specific applications in mind (e.g. hazard maps, real-time forecasting). Thus, no single benchmark could represent all of the numerical codes in a consistent fashion. An exhaustive examination of all possible combinations of parameters and conditions produces a very large number of test scenarios, too large to systematically test codes on. A list of these benchmarks with the code capabilities to resolve them and possible solutions are given in Appendix A, but this paper limits the discussion to a representative subset, including four benchmarks with increasing complexity (BM) and one code comparison test case (CCTC) based on a natural lava flow:

BM1: Dam-break flow test, an initial reservoir of isoviscous Newtonian ‘lava’ is released and spreads on a flat surface. Lava is not interacting thermally with its surroundings.

BM2: Newtonian ‘lava’ is extruded from a point source at a constant known rate, and flows down an inclined plane with a known slope. Lava is isoviscous and does not interact thermally with its surroundings.

BM3: Fluid with a temperature-dependent rheology spreading outwards from a point source on a flat plane. The fluid interacts thermally with the environment. Ultimately cooling will lead to solidification and/or crystallization.

BM4: Lava is extruded onto an inclined plane, where it interacts with an obstacle and splits. Rheology is that of natural basaltic lavas. Lava interacts thermally with its environment.

CCTC1: Based on a natural flow, using conditions such as effusion rate and topography, measured in the field. Attempt to best fit the natural footprint.

### BM1: Dam-break flow

The dam-break flow benchmark describes the transient evolution of the profile of a fixed volume of fluid initially confined in a box geometry, one of sides of which is suddenly removed, similar to the rupture of a dam. BM1 aims at evaluating the adequacy of the code solution of the momentum equation under isothermal conditions. This test is extremely popular in the fluid mechanics discipline, and is considered as a standard measurement when characterizing consistency and yield stress of fluid rheology: hence, numerous experimental data sets are available. This test of rupture is also often used in numerical simulations of transient flow processes with geological applications, such as avalanches or mud flows, and can be used as a benchmark for lava-flow numerical codes. Analytical solutions exist if the assumption of hydrostatic pressure is made so that the problem reduces to a 1D problem (shallow-water equations). This may be generalized to a 2D horizontal plane problem to provide the basis for practical numerical simulations.

Dam-break tests are characterized by two asymptotic flow regimes at short and long times. The description of this test is fully dimensionless and based on existing rheological standards. The reservoir box is defined by its height (*H*) and length (*L*). The width of the box equals the length (*W*=*L*). The aspect ratio (*L*/*H*) of 6.6, of the common standard, is taken (see Fig. 2); a code aiming to reproduce lavas must use the long-time asymptotic regime solution of this dam-break test. Here we test a Newtonian rheology, and the asymptotic solution of the front velocity with time is described using the exponent power law of 0.5 for short times and 0.2 for long times (see equation 2). Note, however, that non-Newtonian solutions for plastic fluids also exist and could be used for a future benchmark. For numerical codes that use dimensional input we recommend the following parameters: *L*=6.6 [m]; *H*=1 [m]; total box length, *L*_{T}=11*L*=72.6 [m]. The fluid should have a density of ρ=2700 [kg m^{−3}] and a viscosity η=10^{4} [Pa s]. Simulations should run for *t*= 1.65×10^{6} [s], by which the expected run-out distance is *L*_{R}=10*L*=66 [m].

Results are compared to the solution proposed in Balmforth *et al.* (2007) and Saramito *et al.* (2013). The final front evolution is described by:
(2)
with the characteristic time *T*=(*L*/*H*)^{2} (η/ρ*gH*). Similarly, the height evolution of the fluid at the dam position may be described for longer times in the reservoir section as:
(3)

### BM2: Inclined viscous isothermal spreading

This benchmark is based on the analytical solution derived by Lister (1992). A Newtonian liquid of viscosity η and density ρ is injected at a constant supply rate, *Q*. The isothermal flow spreads onto an inclined plane of slope α (Fig. 3a). This benchmark, combined with BM1, provides useful results on how well fast-computing models handle a simple viscous flow geometry.

In order to test whether a code can reproduce the scaling laws derived by Lister (1992), we propose the following set of parameters that follow the set-up of the analogue experiments (using silicone oil) by Lister (1992): point source, inclined plane at a slope of 2.5° from horizontal, with a fluid supply rate of *Q*=1.4810^{−6} m^{3} s^{−1} and a kinematic viscosity of η/ρ=11.3× 10^{−4} [m^{2} s^{−1}]. Figure 3b presents the experimental observed solution for flow contour. The theoretical solution derived by Lister (1992) yields scaling solutions for the flow front advance *L*_{d}(*t*), and cross-slope extent *y*_{P} at ‘long times’ (once the flow is asymmetric):
(4)
(5)

The up-slope extent (*L*_{u}) is predicted to rapidly reach a steady value. Code outputs would be tested against their ability to reproduce the scaling of equations (4) and (5) as a function of time (*t*) and of supply rate (*Q*).

### BM3: Axisymmetric cooling and spreading

This benchmark is the first to deal with a non-isothermal flow. It is an intermediate step between the ‘spreading’ benchmarks (BM1 and BM2) and the ‘complex’ lava benchmarks (BM4 and CCTC1). Temperature is only considered as a passive tracer in this benchmark, with no feedback between thermal structure and rheology. BM3 is based on analogue experiments in which hot silicone oil is injected at a constant supply rate and cools as it spreads in the ambient air (Garel *et al.* 2012). Figure 4 shows the evolution of flow contours and surface temperatures with time. Simulation inputs are shown in Table 2.

Here we compare the evolution of dimensionless radial surface temperature profile during experiment C14 of Garel *et al.* (2012) (see Fig. 4). Dimensionless temperature is (*T*−*T*_{a})/(*T*_{0}−*T*_{a}), with *T*_{0} and *T*_{a} the source and ambient temperatures, respectively. At long times, temperature profiles are the same, indicating the establishment of a thermal steady state. Radial flow advance, *R*(*t*), follows the theory of Huppert (1982) with .

### BM4: Split flow experiment

This test case is based on experimental results obtained at the Syracuse University Lava Project laboratory (lavaproject.syr.edu). The material used was natural basalt, heated above the liquidus and fully degassed. The flow interacted with a triangular steel obstacle located at the flow centre, 45 cm from the start of the sandy bed. The obstacle side walls are each 14 cm long, and the opening angle between them is 90°. Figure 5 shows snapshots from the experiment, as well as the position of the flow front over time. Simulation inputs are shown in Table 3.

Measurements from the experiment that can be compared with numerical simulation outputs are the change over time of: (1) the downslope distance of the flow front(s); (2) the maximum width; and (3) the average flow surface temperature.

### CCTC1: Natural case

Most codes simulating lava-flow emplacement have been compared with selected natural cases during their development. However, only a few natural flows have been carefully monitored and characterized during their emplacement to be used as benchmark tests. We selected one of the best documented flows to date, the LSF1 lava flow of the 2001 eruption of Mount Etna, Italy (Coltelli *et al.* 2007). The flow covered a total area of 1.95×10^{6} m^{2}, had a total volume of 21.40×10^{6} m^{3} and an average height of 11 m.

The LSF1 vent was located on the south flank of Etna (UTM coordinates: 500506 East, 4173306 North, zone 33; DMS coordinates: 37°42′ 25.171″N–15°0′20.665″E). The DEM prior to the 2001 eruption, courtesy of C. Del Negro and his team (INGV Catania), is provided under UTM coordinates and is available at the community depository (see the section ‘Future plans: a public facility for code evaluation’ later in this paper).

The effusion rate is extracted from Coltelli *et al.* (2007) (see Table 4) and the simulation parameters established are shown in Table 5.

The apparent viscosity (η_{ap} [Pa s]) is the joint effect of crystals plus solidifying melt, and is given by the equation from Krieger & Dougherty (1959):
(6)
where the pure melt viscosity (η_{m} [Pa s]) is linked to the temperature through the Vogel–Fulcher–Tammann equation:
(7)
where *a*=−4.827 (−6.429, −3.225), *b*=5997 (3553, 8441) and *c*=−330.3 (−448.5, −212.2), coefficients with 95% confidence bounds, *T* is the temperature in °C, Φ is the crystal fraction, Φ_{m} the maximum loose packing crystal fraction (0.65), and *B* the intrinsic viscosity parameter (3.5). It is noteworthy that these parameters may not be optimal and, with regard to the crystal shape at Mt Etna, a combination of Φ_{m}=0.35 and *B*=6 may be more realistic. However, the apparent viscosity for such crystal suspensions and for various dynamics (shear-thinning effect) is not yet clearly constrained. In addition, the lower apparent viscosity obtained here could be justified by the additional fraction of exsolved volatiles in the lava. We also arbitrarily opted for a linear increase in the crystal fraction between the liquidus and solidus (1200 and 1000 °C, respectively).

## First benchmark results

The main result of a benchmark is the adequacy between the numerical and theoretical/natural analysis. In the case of a spreading flow benchmark, the main parameters compared are the transient front position (run-out) and/or the inundation area. Dynamic pressure and temperature could be additional outputs to evaluate transient behaviour and thermal cooling. This section presents simulations results as a guideline for some of the benchmarks described previously (BM1 and CCTC1).

BM1 was tested with the codes VOLCFLOW (K. Kelfoun), FLOW3D (E. Lev), OPENFOAM (E. Lev) and NB3D (B. Cordonnier); CCTC1 with DOWNFLOW (M. Favalli), FLOWGO (A. Harris), MAGFLOW (C. Del Negro, A. Cappello & G. Ganci) and VOLCFLOW (K. Kelfoun). On BM1, we plot the dimensionless front position against the dimensionless time for each of the codes tested (Fig. 6). As mentioned above, a code simulating lava flows only needs to reproduce the long-time behaviour of the dam-break test. At long times, all the codes exhibit the expected asymptotic behaviour characterized by a front advancing with time, with a power law of exponent 1/5. However, some significant differences may be observed with the different codes. VOLCFLOW with its SWE-type algorithm logically produces a good match for both short and long time periods. Divergences exist for OPENFOAM and NB3D for short time periods. Concerning NB3D, because the specific SPH algorithm formalism of the equation of state differs from the traditional incompressible formalism, this result may be improved by adjusting the parameters used.

To examine results for CCTC1, we compared the area covered by the numerical solution to the footprint of the LSF1 event from 18 to 28 of July 2001. The quality of the fit between a model simulation result and reality is estimated from a Boolean combination of the 2001 lava footprint (*F*) and computed emplacement (*C*) (Bilotta *et al.* 2012), based on the following definitions: the intersection (*F*∩*C*) provides the common ground between the calculation and the real emplacement (the green patches in Fig. 7). The union (*F*∪*C*) is the total area covered by both the calculation and the real emplacement (the red patches in Fig. 7). Finally, the difference (*F*−*C*) provides the part that is exclusive to the simulation (the blue patches in Fig. 7). From these numbers we establish three values with which to evaluate the accuracy of the simulation:

: gives the fraction of the intersection over the lava footprint. It provides a quantification of how much of the footprint is filled by the calculation.

: gives the fraction of the footprint towards the union of the estimated and observed area. It provides a quantification of the calculation overestimation of the flow; a ratio of 1 meaning no overflow.

: gives the intersection area over the total surface. It combines the advantages of the two previous numbers and can be taken as the most representative number of the code accuracy.

The reason for having three numbers instead of one single number is linked to the diversity of the codes tested. As some are not meant to be transient or reproduce the lava footprint, but just generate hazard maps or study specific processes, it is necessary to assess the purpose of these codes.

Overall, results strongly differ depending on the code objectives. For example, DOWNFLOW, dedicated to identifying areas of potential inundations, covers 92.28% of the real lava footprint (Φ_{1}) but has lower accuracy for the other measured values. Conversely, FLOWGO aims to constrain the extent of the lava flow, and the LSF1 footprint covers 98.48% of the calculated area (Φ_{2}). For codes aimed at simulating the true emplacement of the flow, such as MAGFLOW or VOLCFLOW, Φ_{1} and Φ_{2} are both high, leading to a reasonable Φ_{3}. Note that among them VOLCFLOW uses an averaged viscosity of 10^{4} Pa s for the whole simulation.

Yet unexpected, but remarkable, results are the similarities observed between MAGFLOW and VOLCFLOW concerning the under- and overestimated areas (blue and red patches, respectively). It denotes that, despite fundamentally different numerical schemes, the two codes provide similar solutions. It already highlights the possible importance of input parameter errors, and the need for higher resolution monitoring and finer measurements of the physical properties of lavas.

## Evaluating code efficiency

The optimal numerical code provides accurate results at cheap computing costs. This is especially important for codes used for mitigation of natural hazards during an eruptive crisis. Therefore, we consider code efficiency as part of the benchmark definition. Comparing several codes under a single benchmark requires a common measure of code efficiency independent of the capacities of the infrastructure used for the test. The most direct approach would be to run all of the codes on the same machine, and compare the computing time required to complete the simulation. As this is unrealistic in a reality where codes are developed by different groups and for different computing platforms, alternative measures are the total number of CPU operations, and how well the numerical code uses the available resources (i.e. number of operations per second, or FLOP, over the theoretical maximal number of operations per second of the machine). The measurement of these two last parameters allows for an evaluation of how well the code might perform on a different platform. For example, a code will not perform better on a more powerful multicore system if not parallelized. We propose to compare benchmark results through three parameters:

*T*_{1}, computing time;*T*_{2}, total number of operations;*T*_{3}, efficiency (mean FLOPS over theoretical maximal FLOPS).

*T*_{1} characterizes how fast a numerical scheme provides an answer during an eruptive crisis, *T*_{2} is the real computing cost of a code and *T _{3}* is the capability of the code to be ported over larger platforms. It is noteworthy that numerical codes are commonly not meant to be compared to each other in terms of efficiency, and the parameters proposed are a compromise between common outputs of codes and a more qualitative measurement independent of the computing architecture used. However, the most objective comparison would be to run every calculation on a similar architecture. For example,

*T*

_{1}is a straightforward measurement for every code and reflects its present capabilities for real natural hazard situations but is extremely dependent on the architecture used.

*T*

_{2}is a more objective evaluation of the required computing for every code and is a qualitative comparison, but even this evaluation could be subject to slight variations between different architectures as some hardware might implement some operations for transcendental functions (e.g. logarithm and exponential functions). Finally,

*T*

_{3}evaluates the efficiency of the code to use the available resources but a serial code (non-parallelized) may not perform well under multicore architecture.

Code performances are provided in Figure 8. As expected, one may observe that the general accuracy of the code (Φ_{3}) increases with the total computing time required and/or the total number of operations (*T*_{2}). An additional test to accurately compare the code would require each code to be run at comparable levels of degrees of freedom to compare similar number of operations or close computing times. Within these first results, some codes are penalized when taken out of their applicability domain. As, for example, DOWNFLOW aimed to evaluate the most likely total inundation area, or FLOWGO aimed to evaluate the path and total flow advance, where flow width only represents the lava channel and not the inundation area. As a conclusion, however, every code performs well within the limit of their scope (inundation area, total extent or global footprint recovery).

## Future plans: a public facility for code evaluation

This paper presents a first attempt at defining a community-driven formal evaluation procedure for numerical codes simulating lava flows. We describe five test scenarios, building from a basic isoviscous dam-break model, then adding complexity in the form of temperature, rheology and topography, simulating laboratory experiments, until the ultimate test of simulating a well-documented natural lava flow.

However, a benchmark is only useful if the results of code testing are shared among the developers and the community using the code. For this purpose, we developed a simple platform for sharing benchmark set-ups and results. We utilize the community facility of VHUB, and created a working group for the benchmarking of volcanic mass flows: http://vhub.org/groups/benchmarking_models. To date, the site holds detailed descriptions of the benchmarks described above (both input parameters and solutions). In the future, example solutions and results of these benchmarks by various codes will also be summarized on the site. We envision that this will resemble the public facility used by the computer vision community (http://vision.middlebury.edu/flow/) to keep an up-to-date status of all models of lava-flow emplacement, and their relative accuracy and performance. Code developers will be able to submit their results, and code users will be able to make an informed choice of a code that is appropriate for their needs. Finally, we hope to add additional benchmarks in the future: for example, the flow of materials with non-Newtonian rheology (see Appendix A).

## Appendix A

A list of benchmarks, the code potential to solve them and the potential solutions are given in Table A1 and Figure A1.

- © 2016 The Author(s). Published by The Geological Society of London. All rights reserved