TephraProb: a Matlab package for probabilistic hazard assessments of tephra fallout

Abstract

TephraProb is a toolbox of Matlab functions designed to produce scenario–based probabilistic hazard assessments for ground tephra accumulation based on the Tephra2 model. The toolbox includes a series of graphical user interfaces that collect, analyze and pre–process input data, create distributions of eruption source parameters based on a wide range of probabilistic eruption scenarios, run Tephra2 using the generated input scenarios and provide results as exceedence probability maps, probabilistic isomass maps and hazard curves. We illustrate the functionality of TephraProb using the 2011 eruption of Cordón Caulle volcano (Chile) and selected eruptions of La Fossa volcano (Vulcano Island, Italy). The range of eruption styles captured by these two events highlights the potential of TephraProb as an operative tool when rapid hazard assessments are required during volcanic crises.

Introduction

Tephra is the most widespread of all volcanic hazards. Variations in deposit thickness and and grain–size, from proximal (i.e. a few kilometres from the vent) to distal locations (i.e. hundreds of kilometres) result in a wide range of potential impacts, from total loss (e.g. human and livestock casualties, destruction of vegetation and crops, collapse of buildings) to less dramatic but highly problematic consequences (e.g. damage to crops, damage to non-structural building elements, disruption of water and electricity supplies; Blong 1984; Jenkins et al. 2014; Wilson et al. 2011). Tephra is also responsible for widespread disruption of civil aviation (Biass et al. 2014; Bonadonna et al. 2012; Guffanti et al. 2009; Scaini et al. 2014). Both the dispersal and sedimentation of tephra can generate complex and multi-faceted impacts at different spatial and temporal scales, affecting physical, social, economic and systemic aspects of societies. These complex impacts require proactive strategies to effectively increase the level of preparedness and mitigate the risk to exposed communities (Birkmann 2006). In order to develop effective proactive risk mitigation strategies, it is necessary to both develop comprehensive probabilistic hazard assessments for tephra dispersal and sedimentation and communicate the outcomes to decision-makers and stakeholders.

Our method for quantifying the hazard from tephra fallout begins with the identification of eruption scenarios based on detailed stratigraphic studies. Eruption scenarios reflect the typical eruption styles at a given volcano and are characterized by ranges of eruption source parameters (ESP) such as plume height, erupted mass, mass eruption rate (MER) and total grain–size distribution (TGSD). Eruption scenarios are constrained by the completeness of the eruptive record, as well as an understanding of the range of past activity observed at analogous volcanic systems (Marzocchi et al. 2004; Ogburn et al. 2016). Highly–studied systems result in eruption scenarios describing precise eruptive episodes, whereas poorly–known systems can often only be characterized by generic eruption scenarios based on analogue systems. Some hazard studies have relied on a deterministic approach, defining ESPs as a single set of best guess values, and hazard as expected tephra accumulation (e.g. isomass maps, kg m⁻²; Barberi et al. 1992). Other strategies rely on a large number of simulations to explore the range of variability in ESPs and atmospheric conditions. In the latter case, the hazard is quantified as a probability to exceed a given tephra accumulation and can be expressed either as conditional to the occurrence of the eruption scenario (e.g. scenario–based hazard assessment; Biass et al. 2016; Bonadonna et al. 2005; Scaini et al. 2014; Volentik and Houghton 2015) or as absolute when the long–term probability of the eruption scenario is also quantified (Bear-Crozier et al. 2016; Jenkins et al. 2012; Marzocchi and Bebbington 2012; Sandri et al. 2014, 2016; Thompson et al. 2015).

The TephraProb package provides a framework for i) accessing, pre–processing, and analyzing model inputs (e.g. retrieval and probabilistic analysis of Reanalysis wind data), ii) building probabilistic eruption scenarios, iii) running a suitable tephra dispersal model and iv) post-processing and exporting model results. The package is written in Matlab and released under an open–source GPLv3 license on VHub and GitHub. This paper provides an in–depth review of the possibilities offered by the TephraProb package. The case studies of Cordón Caulle volcano (Chile) and La Fossa volcano (Vulcano Island, Italy) are used to illustrate the flexibility of TephraProb for working with a wide range of probabilistic eruption scenarios and eruptive styles. The user manual provides a technical description of each functionality and each individual Matlab function. The latest version of the package can be found at https://github.com/e5k/TephraProb.

Case studies

Cordón Caulle volcano

The 2011 eruption of Cordón Caulle volcano (Chile), part of the Puyehue–Cordón Caulle system, is considered as an example of a long–lasting sustained sub–Plinian eruption. Cordón Caulle started erupting on the 4th of June 2011 after a month–long seismic swarm (Bonadonna et al. 2015; Collini et al. 2013; Pistolesi et al. 2015). The initial and most vigorous phase (Unit 1; Bonadonna et al. 2015; Pistolesi et al. 2015) was characterized by plume heights varying between 10–14 km asl during a 24–30 h period with a mean MER of 10⁷ kg s⁻¹. The eruption continued until the 15th of June with MER ≥10⁶ kg s⁻¹ and was followed by long–lasting, low–intensity activity (Bonadonna et al. 2015; Pistolesi et al. 2015). Modelling details of the 2011 eruption of Cordón Caulle can be found in Elissondo et al. (2016).

La Fossa volcano

During the last 1000 years, La Fossa volcano (Vulcano Island, Italy) experienced two sub-Plinian eruptions with similar intensities and magnitudes (i.e. plume heights of 7–8 km asl and erupted masses of 2.1–2.4 ×10⁹ kg) and at least 8 long–lasting Vulcanian cycles (Di Traglia et al. 2013). Each of these Vulcanian cycles is characterized by total fallout masses of 0.1–1 ×10⁷ kg over durations of 3 weeks to 3 years. The most recent eruption of La Fossa in 1888–1890, described in detail by Mercalli and Silvestri (1891) and De Fiore (1922), was characterized by explosions occurring every 4 hours to 3 days and associated with plume heights of 1–10 km asl (Di Traglia 2011; Di Traglia et al. 2013; Mercalli and Silvestri 1891). A detailed description of the probabilistic modelling of these eruption scenarios can be found in Biass et al. (2016).

The Tephra2 model

The Tephra2 model (Bonadonna et al. 2005, 2012; Connor and Connor 2006; Volentik et al. 2009) uses an analytical solution of the advection–diffusion equation to compute the tephra mass accumulation depending on eruptive and wind conditions. Main characteristics of the model include: i) grain–size-dependent diffusion and particle density, ii) a vertically stratified wind, iii) particle diffusion within the rising plume and iv) settling velocities that include variations in the Reynolds number depending on whether particles follow a laminar or turbulent flow regime. Tephra2 requires 3 input files: i) a configuration file that specifies the dynamics of the eruption and the surrounding atmosphere, including location of the vent, plume height, erupted mass, TGSD, range of particle sizes and particle densities, and atmospheric and dispersion thresholds, ii) a table of locations where tephra accumulation is calculated, and iii) a table of wind conditions influencing the horizontal dispersion of tephra, including atmospheric level and the wind speed and wind direction at this level. Note that since Tephra2 solves the advection–diffusion equation analytically, the wind is vertically stratified but horizontally and temporally homogeneous, thus limiting the validity of our approach to a few hundred kilometers around the vent. The user is referred to Bonadonna et al. (2005); Connor and Connor (2006) and Volentik et al. (2009) for a complete description of Tephra2 and to Scollo et al. (2008a,b) and Bonadonna et al. (2012) for a comparison with other volcanic ash transport and dispersion models (VATDM).

Tephra2 simplifies atmospheric turbulence by using two different diffusion laws, one for coarse and another for fine particles (Bonadonna 2006; Bonadonna et al. 2005; Connor and Connor 2006; Volentik et al. 2009). These two diffusion regimes of fallout are based on three empirical parameters, namely the fall–time threshold (FTT) acting as a threshold for the modelling of the diffusion of small and large particles (i.e. power law vs. linear diffusion), a diffusion coefficient used for the linear diffusion law and an apparent eddy diffusivity fixed at 0.04 m ² s⁻¹ for the power law diffusion (Bonadonna et al. 2005; Suzuki 1983). Secondary atmospheric effects (e.g. topography effects) are neglected. The reader is referred to Bonadonna et al. (2005); Connor and Connor (2006); Volentik et al. (2009) and Volentik et al. (2010) for further information regarding the modeling parameters used by Tephra2. The version of Tephra2 used here was modified to accept an additional input file defining the TGSD, which accounts for aggregation processes (e.g. Biass et al. 2014; Bonadonna et al. 2002a). Aggregation is further discussed in Section “Total grain–size distribution”.

The probability P(x) of tephra released from an eruption plume at some vertical height above the vent is controlled by a beta probability density function (PDF):

$$ P(x) = \frac{(1-x)^{\beta-1}x^{\alpha-1}}{B(\alpha,\beta)}, $$

(1)

where x is a dimensionless height normalized to the plume height and B(α,β) represents the beta function defined by two parameters, α and β, both greater than 0. Changing the values of α and β changes the shape of the PDF. When α=β=1, the beta function is identical to a uniform random function and the probability of tephra release is equal along the height of the plume. When α>β, the release pattern is shifted towards the top of the plume; when α=β≠1 tephra release is greatest from the middle of the plume; when α<β the release pattern is shifted toward the base of the plume.

The TephraProb package

The TephraProb package contains a set of Matlab functions integrated as a graphical user interface (GUI) assisting the user with every step required to produce comprehensive hazard assessments for tephra fallout. The GUI provides an interface facilitating programming and computation aspects, but more advanced users can customize the Matlab functions finding assistance from the detailed comments within the code and the user manual. Figure 1 summarizes the main functionality of the code, which consists of four main modules including i) Input Parameters: to retrieve and analyse the various inputs required by probabilistic hazard assessments, ii) Eruption Scenarios: to stochastically sample ESP and run eruption scenarios, iii) Post–processing: to process individual Tephra2 runs into probabilistic outputs and iv) Output: to visualize and export results. The following sections use the two case studies presented in Section “Case studies” to illustrate every step of the tephra fallout hazard assessments.

Fig. 1

Summary of the workflow of the TephraProb package

Input parameters

Calculation points

Tephra2 requires an ASCII file of locations arranged in a three-column format containing Easting, Northing, and elevation. Since Tephra2 uses an analytical solution to the advection-diffusion-equation that solves the differential equation using a constant elevation boundary condition, calculations must be performed on flat surfaces and topography is neglected (Lim et al. 2008). Two options are available in TephraProb. The first option is to perform calculations on regularly spaced grids allowing for the compilation of probability maps. Although this is often the preferred option, calculations using large grids is time consuming and might not be achievable in some emergency contexts if computing power is limited. The second option is to perform calculations for user–selected points of interest, which results in a series of hazard curves that require a considerably shorter computation time.

One typical problem arising from the use of UTM coordinates is the shift in Easting coordinates when crossing UTM zones and the shift in Northing coordinates when crossing the equator. In both cases, TephraProb includes a correction to produce grids with contiguous UTM coordinates, where the resulting distortion is assumed negligible compared to the maximum distance deemed valid for the application of Tephra2.

Wind data

Wind velocity is a key input parameter in any tephra dispersal model because it controls the advection of the plume and the sedimentation of tephra. Probabilistic modelling requires access to wind datasets spanning at least one decade in order to capture the variability of wind conditions over a region. Most hazard assessments for tephra dispersal and sedimentation rely on Reanalysis datasets, which provide access to decades of atmospheric observations continuously interpolated in space and time. The two most frequently used datasets are the National Oceanic and Atmospheric Administration (NOAA) National Centers for Environmental Prediction NCEP/NCAR Reanalysis (i.e. joint project between the National Centers for Environmental Prediction and the National Center for Atmospheric Research; Kalnay et al. 1996) and the European Centre for Medium-Range Weather Forecasts (ECMWF) Era-Interim datasets (Dee et al. 2011). Both the NOAA NCEP/NCAR and the ECMWF Era-Interim datasets can be accessed via GUIs using TephraProb. The TephraProb package downloads the files in the NetCDF format and converts them into three–column ASCII files including altitude, wind direction (i.e. direction the wind blows towards) and wind velocity. Once the files are downloaded and converted, they can be used as an input to any tephra model and/or to explore the wind conditions over a region of interest. Note that Tephra2 uses a single wind profile per simulation extracted for one set of coordinates and one time. Therefore, each simulation is performed in an atmosphere that is vertically stratified but horizontally and temporally homogeneous.

Figure 2 illustrates the interface and the outputs resulting from wind analysis above Cordón Caulle volcano using the ERA–Interim dataset for the period 2001–2010. Figure 2 a shows median profiles of wind velocity and direction with height for the entire population. Figure 2 b and 2 c show the median wind direction averaged per year and per month, respectively. Such plots are useful to identify potential variability associated with El Niño/La Niña oscillations (Fig. 2 b) or seasonality (Fig. 2 c). Finally, wind roses (Fig. 2 d) show, for a fixed altitude, the probability that wind blows towards a given direction and at a given intensity.

Fig. 2

Wind statistics produced by the TephraProb package as inferred from the ERA–Interim dataset for the 2001–2010 period above Cordón Caulle volcano. a Screenshot of the GUI shows the median wind velocity and median wind direction with height (solid black lines) including the 25^th and 75^th percentiles (error bars); b median wind direction per year; c median wind direction per month; d probability of the wind to blow in given directions and intensities at an altitude of 10 km asl

Eruptive history

Hazard assessments are based upon the assumption that future activity will be similar to past activity. In volcanology, this implies a necessity to develop scenarios of future eruptions based on the most complete possible study of the geological record in order to best constrain the eruptive history of a volcano. In practice, assessing the full spectrum of eruptions of a given system is impossible, due firstly to the fact that geological records are not continuously complete through time and secondly because field studies rarely provide a complete picture of the eruptive history (Kiyosugi et al. 2015). In most cases, eruption scenarios are constructed based on a few studied eruptions representing common eruptive styles and/or often relating underlying geologic processes with analogue eruptive events (Sheldrake 2014). Eruption databases such as the Global Volcanism Program of the Smithsonian Institute (GVP; Siebert et al. 2010; Simkin and Siebert 1994) or the Large Magnitude Explosive Volcanic Eruptions (LaMEVE; Crosweller et al. 2012, Ogburn et al. 2016) can prove useful in acquiring a global picture of eruptive history. TephraProb includes a module to access and explore the Holocene eruptive history as recorded in the GVP database. Note that the GVP module does not constitute a source of input in itself but serves as a support tool to identify eruption scenarios and key ESPs.

Figure 3 a illustrates the GUI of the GVP module, which retrieves eruptive history based on the volcano number (e.g. 357150 for the Puyehue-Cordón Caulle system and 211050 for La Fossa). Key features of the GVP database are preserved (e.g. VEI, confirmation and evidence; Siebert et al. 2010; Simkin and Siebert 1994) and the data can be plotted either as a histogram or as a cumulative distribution. Considering the entire Holocene catalogue, Fig. 3 a shows the frequency of eruptions per VEI class and Fig. 3 b shows the cumulative number of eruptions through time. The various segments characterized by distinctive slopes in Fig. 3 b illustrate the complications arising when unifying geological and historical records to compile a comprehensive eruptive history, discussed by Siebert et al. (2010). On one hand, the geological record is biased toward large eruptions and extends back in time. On the other hand, the historical record is detailed and preserves small and moderate eruptions but does not sufficiently extend back in time to capture eruptions of larger magnitudes. Figure 3 c illustrates the historical segment of the eruptive history of the Puyehue-Cordón Caulle system as inferred from the GVP database (i.e. break in slope at 1880 A.D. on Fig. 3 b). Note that breaks–in–slope associated with completeness might vary as a function of the VEI/magnitude class and should be assessed separately (Biass and Bonadonna 2013; Dzierma and Wehrmann 2010; Jenkins et al. 2012; Mendoza-Rosas and De la Cruz-Reyna 2008).

Fig. 3

Eruptive history of Puyehue-Cordón Caulle volcanic system as recorded in the GVP database. a Screenshot of the GUI shows the number of eruptions per VEI class for all eruptions recorded in the GVP, where U stands for undefined and O for other. Alternative plots include, cumulative number of eruptions considering (b) all eruptions and (c) the most recent complete segment of the record, (d) and probability of eruption through time considering a Poisson model for all eruptions (black line, using the record shown in (c)), eruptions of VEI 2 (red line, completeness at 1880 A.D.), VEI 3 (green line, completeness at 1920 A.D.) and VEI 5 (orange line, completeness at 5000 B.C.)

If the eruptions are assumed to occur i) with a constant rate and ii) independently of the time since the last event, it is possible to obtain a first–order estimate of probability of occurrence of an eruption in a given time window by considering a Poisson process (e.g. Biass and Bonadonna 2013; Borradaile 2003). We aim at quantifying the probability that a repose interval T is smaller or equal to an arbitrary time window t:

$$ F(t) = P(T \geq t) $$

(2)

which, in the simplest case of a Poisson process, results in an exponential distribution:

$$ F_{exp}(t) = 1-e^{-\lambda t}, $$

(3)

where t is a forecasting window (years) and λ is the eruption rate (number of eruptions per year) defined over a complete section of the eruptive catalog. As an example, the black line in Fig. 3 d shows the probability of eruption through time for a case where Fig. 3 c represents a complete historical record. In contrast, the colored lines plot probabilities for eruptions of VEI 2 (red line; completeness at around 1880 A.D.), VEI 3 (green line; completeness at around 1920 A.D.) and VEI 5 (orange line; completeness at around 5000 B.C.).

Note that the aim of TephraProb is scenario–based hazard assessments, i.e. based upon the conditional probability that the associated eruption scenario occurs; two caveats should be considered when using the GVP module. First, the GVP database is not a direct input to the probabilistic modelling in TephraProb and should be used as a support tool to develop eruption scenarios and identify critical ESPs with the full knowledge of the limitations and assumptions behind such databases (Biass and Bonadonna 2013; Dzierma and Wehrmann 2010; Jenkins et al. 2012; Mendoza-Rosas and De la Cruz-Reyna 2008; Siebert et al. 2010; Simkin and Siebert 1994). Second assessing the long–term probability of a future eruption is not trivial and should be achieved using a rigorous probabilistic framework in order to quantify and propagate various sources of uncertainties on final estimates (Bear-Crozier et al. 2016; Connor et al. 2003; Jenkins et al. 2012; Marzocchi and Bebbington 2012; Sandri et al. 2016; Sheldrake 2014; Thompson et al. 2015). Probabilities resulting from a Poisson process should not be viewed as more than a first–order estimate, and only when the hypotheses of stable rate and independent events are satisfied (Borradaile 2003).

Eruption scenarios

Based upon considerations described in Section “Eruptive history”, probabilistic strategies were introduced to account for biases associated with both small and large eruptions. On one hand, they help account for parts of the geological record that are not accessible or that have been removed. On the other hand, they allow for the consideration of larger eruptive events associated with lower probabilities of occurrence that have not taken place during the recent history but cannot be excluded. Probabilistic eruption scenarios rely on a large number of simulated events to explore the uncertainty associated with variable parameters. For tephra hazard, the two main variable parameters are magnitude/style and atmospheric conditions of a future eruption. The variability of a future eruption is commonly described by sets of critical ESPs defined as ranges and/or probability distributions, where the shape of a distribution reflects the degree of certainty for a particular scenario. The maximum degree of uncertainty for a given ESP will be typically described by a uniform distribution, whereas a logarithmic distribution favors the occurrence of smaller values compared to larger ones, and a Gaussian distribution preferentially samples around a central value. Similarly, some regions of the world experience seasonal wind trends, making a seasonal analysis of tephra hazard necessary. Identifying the degree of variability allowed during stochastic sampling requires subjective choices that are influenced by both the degree of knowledge of a volcanic system and the purpose of the hazard assessment.

Total grain–size distribution

Plume height, erupted mass and TGSD are the three most important ESPs for modeling tephra accumulation with Tephra2. In TephraProb, the sampling of plume heights and erupted masses are inherent to the probabilistic eruption scenarios described below, but the stochastic scheme for sampling TGSDs is similar for all scenarios and is therefore described first. TGSDs in TephraProb are assumed Gaussian (in ϕ units) and are defined by ranges of median diameters and sorting coefficients. It is then possible to modify the TGSD to account for aggregation using the empirical approach proposed by Bonadonna et al. (2002a) and Biass et al. (2014) based on observations from Cornell et al. (1983), Bonadonna et al. (2011), and Bonadonna et al. (2002b) by defining a range of values for the empirical aggregation parameter and a maximum diameter affected by aggregation processes.

The empirical aggregation parameter introduced in TephraProb follows Bonadonna et al. (2002a) and Cornell et al. (1983) and represents a weight fraction of particles that form aggregates. In our approach, a fraction of mass equal to the empirical aggregation parameter is removed from all bins equal to or finer than the maximum diameter (Biass et al. 2014; Bonadonna et al. 2002a). The total mass removed is then equally redistributed into bins ≤−1ϕ and >the maximum diameter (where > means coarser than). Although this aggregation scheme is a simple answer to a complex concept (e.g. Brown et al. 2012; Gilbert and Lane 1994; James et al. 2002, 2003; Rose and Durant 2011; Van Eaton et al. 2015) and is based uniquely on observations (Bonadonna et al. 2002b, 2011; Cornell et al. 1983), it is computationally efficient and suitable for hazard assessment purposes (Biass et al. 2014, 2016; Bonadonna et al. 2002a). This scheme has already been validated in previous studies with good agreement with field observations (e.g. Bonadonna et al. 2002a; Bonadonna and Phillips 2003; Cornell et al. 1983). Since no physics is involved in this method, it is not possible to draw a distinction between dry or wet aggregates. Instead, following the nomenclature proposed by Brown et al. (2012), using a maximum diameter of 5 ϕ suggests that particles <63 microns are influenced by aggregation processes (e.g. ash clusters, coated particles and poorly-structured pellets; Bonadonna et al. 2011), whereas a maximum diameter of 4 ϕ extends aggregation up to <125 microns (e.g. pellets with concentric structures (e.g. accretionary lapilli) and liquid pellets (i.e. mud rain); Gilbert and Lane 1994; Van Eaton et al. 2012, 2015).

At each probabilistic run, values of median diameters and sorting coefficients are sampled and used to create a Gaussian distribution. If wanted, the aggregation scheme is applied by randomly sampling an empirical aggregation parameter. The TGSD is then written as a text file and passed as an input to Tephra2.

General eruption scenarios

Probabilistic eruption scenarios implemented in TephraProb allow for the stochastic sampling of eruption source parameters and/or wind conditions, and various probability distributions are proposed for the sampling of ESPs (Table 1; Figs. 4 and 5). The two most general eruption scenarios implemented in TephraProb are (Bonadonna 2006):

• Eruption Range Scenarios (ERS), designed to assess the probability of accumulation of a critical tephra accumulation based on the statistical distribution of possible eruption source parameters and wind conditions (i.e. stochastic sampling of one set of ESPs and one wind profile at each run);

  Fig. 4

  

  Workflows defined for the sub–Plinian and Plinian eruption scenarios implemented in TephraProb. Get refers to ESPs deterministically defined, Set defines a variable and Sample indicates a stochastic sampling. The index i refers to the run number, where the total number of runs is given by nr. For long–lasting eruptions, index j refers to the simulation number of run i, where the total number of simulations for run i is given by nw

  Fig. 5

  

  Workflows defined for the Vulcanian eruption scenarios implemented in TephraProb. The nomenclature follows Fig. 4

  Table 1 Summary of probabilistic eruption scenarios implemented in TephraProb. All Plinian–type scenarios can be modeled as Long–lasting if the eruption duration is longer than the wind sampling interval

• One Eruption Scenarios (OES), designed to assess the probability of accumulation of a critical tephra accumulation based on the statistical distribution of wind conditions and a deterministically–defined eruption (i.e. stochastic sampling of one wind profile at each run while ESPs remain constant).

Derivations of these two scenarios are implemented in TephraProb in order to allow additional subjective choices to the stochastic sampling (Table 1):

• Wind Range Scenarios (WRS) constrain the sampling of wind conditions within pre–defined radial sector around the volcano. Such scenarios are useful to assess the hazard at specific sites considering a specific wind scenario (e.g. Volentik et al. 2009). The probability of occurrence of the wind scenario can itself be assessed using the wind module of TephraProb;

• Fixed Date Scenarios (FDS) fix the eruption starting date, thus assessing the probability of tephra accumulation based on the variability of eruption source parameters only. Rather than expressing the hazard of a future eruption, such scenarios are useful to assess the probability of a studied fallout event considering the statistical distribution of wind conditions (e.g. Elissondo et al. 2016). Outcomes of Fixed Date Scenarios can be compared to Eruption Range Scenarios and field–based isomass maps.

Eruptive styles

Different assumptions are made in TephraProb to calculate the erupted mass of sustained sub–Plinian and Plinian eruptions and nonsustained Vulcanian explosions, giving rise to Plinian–type and Vulcanian–type eruptions (Table 1; Figs. 4 and 5). Note that i) in our modelling scheme, sub–Plinian eruptive styles are hereafter considered as a part of Plinian–type eruptions and ii) both Plinian and Vulcanian–type styles can be modeled using the eruption scenarios described in Section“General eruption scenarios ”.

Plinian–type eruptions

Plinian–type eruption scenarios are defined by ranges and distribution shapes of plume height, eruption duration and eruption start date (Fig. 4). Additionally, the user also defines a range of erupted tephra mass. The sampling scheme, developed by Biass et al. (2014) and summarized in Fig. 4, is then applied. At each run, an eruption date, a plume height and an eruption duration are first sampled. The eruption date is used to retrieve the corresponding wind profile and combined with the plume height to calculate the MER following Degruyter and Bonadonna (2012):

$$ MER = \pi\frac{\rho_{a0}}{g'}\left(\frac{\alpha^{2}\bar{N}}{10.9}H^{4} + \frac{\beta^{2}\bar{N}^{2}\bar{v}}{6}H^{3}\right) $$

(4)

where ρ _a0 is the reference density of the surrounding atmosphere (kg m⁻³), g ^′ the reduced gravity at the source (m s⁻²), α is the radial entrainment coefficient, $\bar {N}$ is the average buoyancy frequency (s ⁻¹), H is the plume height (m above the vent), β is the wind entrainment coefficient and $\bar {v}$ the average wind velocity across the plume height (m s⁻¹). Note that Reanalysis datasets are only used to calculate $\bar {v}$, and the calculation of other parameters (e.g. $\bar {N}$, ρ _a0) follows from Degruyter and Bonadonna (2012). The MER is used with the eruption duration to calculate the erupted tephra mass. If the resulting value falls within the user–defined mass range, then the run is sent to the model to be executed, or else the sampling process is re–started. Note that TephraProb also includes an option to sample plume heights and erupted tephra mass independently.

Vulcanian–type eruptions

In contrast, Vulcanian–type eruption scenarios are modeled as thermals (i.e. instantaneous release of mass) (Fig. 5). At each run, a plume height is sampled and used to calculate the mass of a thermal with the relationship (Bonadonna et al. 2002a; Druitt et al. 2002; Woods and Kienle 1994):

$$ H = 1.89Q^{0.25} $$

(5)

where Q=f×M×C×Δ T is the excess thermal mass of the thermal injection, f is the solid mass fraction capable of loosing heat to the plume, M (kg) is the plume mass, C (J kg ⁻¹ K ⁻¹) is the pyroclasts specific heat and Δ T (K) is the initial temperature contrast between the erupted mixture and the surrounding air. For Soufrière Hills volcano (Montserrat), f is taken as 0.8, C is 1100 J kg ⁻¹ K ⁻¹, Δ T is 800 K (Bonadonna et al. 2002a; Druitt et al. 2002), and the relationship between the Vulcanian plume height and the mass of the plume can be expressed as:

$$ H = 55M^{0.25} + H_{V} $$

(6)

where H is the plume height (m asl), M is the plume mass (kg) and H _V the vent height (m asl).

Long–lasting eruption scenarios

The hazard associated with long–lasting eruptions is complicated by the fact that wind conditions might vary during the course of an eruption. Long–lasting eruptions can be associated either with sustained plumes (e.g. 2010 eruption of Eyjafjallajökull, Iceland), or with repetitive emissions of limited to moderate amounts of ash (e.g. ongoing eruption of Soufrière Hills volcano, Montserrat, West Indies; ongoing eruption of Sakurajima volcano, Japan). In TephraProb these two eruptive styles are modelled as long-lasting Plinian and long-lasting Vulcanian eruptions, respectively (Figs. 4 and 5).

Long–lasting Plinian eruptions

Long–lasting Plinian scenarios are constrained by the temporal resolution Δ w of available wind data, which is 6 h for most Reanalysis datasets but can be varied within TephraProb. In cases where the eruption duration defined by the user is longer than the temporal resolution between wind profiles, the eruption is divided into Δ w–duration periods. The corresponding wind profile for each period is retrieved and the Eruption Range Scenario sampling scheme is performed, assuming that ESPs remain steady throughout the duration of the period (Fig. 4). If the sums of the masses of independent periods fall within the initial mass range, each period becomes a separate Tephra2 simulation, and the final accumulation of a long–lasting eruption is the sum of all periods.

Long–lasting Vulcanian eruptions

Long–lasting Vulcanian scenarios describe Vulcanian cycles and require the identification of the total duration of the cycle and the repose interval between explosions in order to apply the workflow shown in Fig. 5. At each run (i.e. each Vulcanian cycle), repose intervals and associated plume heights are constantly sampled until the sum of the repose intervals exceeds the eruption duration. The mass of each explosion is calculated with Eq. 6 and used in a separate Tephra2 simulation. The final accumulation of a Vulcanian cycle is the sum of all explosions. Note that TephraProb does not account for any relationship describing the influence of the repose interval on the corresponding plume height. More details on modeling Vulcanian cycles with TephraProb can be found in Biass et al. (2016).

Seasonality

Because some regions at low latitudes experience significant changes in wind patterns between dry and rainy seasons, a seasonality option is implemented in TephraProb to assess the variability in hazard as a function of the time of year. If the seasonality option is enabled, the user can define start and end dates of the rainy season and three scenarios will be run using winds for i) all months of the year, ii) months of the rainy season and iii) months of the dry season.

Post-processing

Post-processing functions in TephraProb process the individual outputs of a given eruption scenario into probabilities of exceeding given thresholds of tephra accumulations. Following Bonadonna (2006), we quantify the probability of hazardous thresholds of mass accumulations using:

$$P[M(x,y) \geq M_{T} \mid eruption] $$

where M(x,y) is the tephra mass accumulation (kg m⁻²) accumulated at given locations and M _T a mass accumulation threshold. For a given eruption scenario, the probability P _M at coordinates x,y is calculated based on the number of times a given threshold of accumulation is reached and divided by the total number of runs N _R :

$$P_{M}(x,y) = \frac{\sum_{i=1}^{N_{R}}n_{i}}{N_{R}} $$

$$n_{i} = \left\{ \begin{array}{l l} 1 & \quad \text{if \(M_{i}(x,y) \geq\) threshold \(\mid\) eruption}\\ 0 & \quad \text{otherwise.} \end{array} \right.$$

Outputs

Two scenarios for each volcano in Section “Case studies” were used as examples for TephraProb. First, the climatic phase (i.e. 4th of June) of the 2011 eruption of Cordón Caulle served as a case-study to model long–lasting sub-Plinian/Plinian eruptions. Eruption scenarios include a Long–Lasting Eruption Range Scenario (i.e. variable ESPs and wind conditions) and a Long–Lasting Fixed Date Scenario (i.e. variable ESPs and wind conditions set for the 4th of June 2011). ESPs for each scenario were defined from Bonadonna et al. (2015) and Pistolesi et al. (2015) and are summarized in Table 2. Second, La Fossa volcano was used to illustrate short–lasting sub–Plinian/Plinian eruptions and long–lasting Vulcanian cycles. The sub–Plinian eruptions were simulated using a One Eruption Scenario, whereas Vulcanian cylces were simulated using a Vulcanian Long–Lasting Eruption Scenario. Table 2 summarizes ESPs for both scenarios based on Mercalli and Silvestri (1891); Bianchi (2007) and Di Traglia et al. (2013).

Table 2 ESPs for the Cordón Caulle and the La Fossa case studies

Four variables need to be considered when displaying the probability of exceeding a given tephra accumulation, including geographic coordinates (Easting, Northing or longitude, latitude), a threshold of tephra accumulation and its associated exceedance probability. Since typical maps are limited by three dimensions, it is necessary to fix at least one degree of freedom. TephraProb can produce three main types of outputs.

Probability maps

Probability maps fix a threshold of tephra accumulation (kg m⁻²) to contour the associated spatial probability of exceedance, based upon the conditional probability that the given eruption scenario occurs (e.g. Fig 6). In Fig. 6, the minimum displayed probability is set to 0.1 and the red line contours the extent of the computational grid.

Fig. 6

Hazard maps show the conditional probability of exceeding a threshold of tephra accumulation given the occurrence of the associated eruption scenario: a Long–lasting Eruption Range Scenario for Cordón Caulle volcano for a tephra accumulation of 10 kg m⁻²; b Long–lasting Fixed Date Scenario for Cordón Caulle for a tephra accumulation of 10 kg m⁻²; c One Eruption Scenario for a sub–Plinian eruption of La Fossa for a tephra accumulation of 10 kg m⁻²; d Vulcanian Long–Lasting Eruption Scenario for a Vulcanian eruption of La Fossa for a tephra accumulation of 100 kg m⁻². ESPs for all scenarios are summarized in Table 2. The minimum displayed probability is 0.1 and the red line contours the extent of the calculation grid. The vent is indicated by a red triangle

Figure 6 a–b show probability maps compiled for an accumulation of 10 kg m⁻² (i.e. critical threshold for crops) for the Long–Lasting Eruption Range Scenario and the Long–Lasting Fixed Date Scenarios of Cordón Caulle based on the climatic phase of the 2011 eruption (i.e. Unit 1 in Bonadonna et al. 2015; Pistolesi et al. 2015, Table 2). These maps show the importance of wind patterns in the hazard assessment of tephra fallout and, put in perspective of Fig. 2, show that the opening phase of the 2011 eruption (i.e. 4–5th of June) occurred in wind conditions of low probability but with high consequences for the town of San Carlos de Bariloche, which experienced accumulations of about 5 kg m⁻² (Bonadonna et al. 2015; Pistolesi et al. 2015).

Figure 6 c–d show the probabilities of exceeding tephra accumulations of 10 and 100 kg m⁻² following sub–Plinian OES (Fig. 6 c) and Vulcanian V–LLERS (Fig. 6 d) eruption scenarios at La Fossa volcano, respectively. Note that in the case of long–lasting Vulcanian scenarios, no syn–eruptive erosion between explosions is considered and Fig. 6 d should therefore be regarded as a probability map of the maximum accumulation.

Hazard curves

Hazard curves fix the geographical location to display the probability of exceeding any tephra accumulation at a given point (e.g. Bonadonna 2006). Figure 7 shows hazard curves for San Carlos de Bariloche using both Cordón Caulle scenarios described in Table 2. The dashed grey line on Fig. 7 corresponds to the accumulation of 5 kg m⁻² reached in San Carlos de Bariloche during the 2011 eruption. Probabilistic analyses show that such an accumulation had an ∼2 % probability of occurrence considering the statistical distribution of wind directions of the past 10 years.

Fig. 7

Hazard curves showing the conditional probability of exceeding any tephra accumulation given the occurrence of the associated eruption scenario at San Carlos de Bariloche. The dashed grey lines show an accumulation of 5 kg m ⁻² corresponding to the accumulation reached in San Carlos de Bariloche during the 2011 eruption of Cordón Caulle volcano

Probabilistic isomass maps

Probabilistic isomass maps fix a probability threshold to represent the typical tephra accumulation given a probability of occurrence of the hazardous phenomenon conditional on the occurrence of the associated eruption scenario (e.g. Biass and Bonadonna 2013; Biass et al. 2016). Figure 8 shows isomass contours (kg m⁻²) for a LLERS eruption of Cordón Caulle for probabilities of occurrence of 25 % (Fig. 8 a) and 75 % (Fig. 8 c). The choice of the probability threshold (which may be regarded as an acceptable level of hazard), is critical, and must be determined in collaboration with decision makers. The resultant probabilistic isomass maps should be presented in a style and format that meets the maps’ communication goals and decision makers’ needs. However, their use facilitates the incorporation of probabilistic approaches into pre–event impact assessments and are useful when communicating typical eruption scenarios to stakeholders.

Fig. 8

Probabilistic isomass maps showing the typical tephra accumulations (kg m⁻²) for a conditional probability of occurrence of the hazard of 25 % (a) and 75 % (b) for a Long–lasting Eruption Range Scenario scenario at Cordón Caulle volcano

Discussion and conclusions

The TephraProb package aligns with recent efforts to bridge gaps between academic and operational contexts (Bartolini et al. 2013; Bear-Crozier et al. 2012; Felpeto et al. 2007). It provides a set of flexible tools for the assessment of tephra hazard and designed to compile comprehensive assessments in a wide range of conditions (e.g. poor constraints on scenarios, rapid assessments), that can be used in the context of probabilistic hazard assessments or separately (e.g. accessing and analyzing wind data or the GVP database). TephraProb integrates various levels of computational requirements and allows hazard assessments to be performed on single CPU computers (i.e. hazard curves only), multi–core personal computers (i.e. grids of moderate resolutions) and computer clusters (i.e. grids of fine resolutions). Resulting hazard assessments are conditional on the occurrence of the associated scenarios, and can serve as direct inputs to probabilistic frameworks such as Bayesian event trees to assess the long–term probability of tephra accumulation (Thompson et al. 2015) or to fit in multi–hazards assessments (Sandri et al. 2014). In order to facilitate the integration of further analyses, each output of TephraProb is saved in a variety of formats (e.g. ASCII columns format, ASCII ArcMap rasters).

The TephraProb package is currently tailored to generate Tephra2–friendly configuration files. However, probabilistic strategies implemented in the package are independent of the adopted model and could be modified to work with any VATDM. For instance, early versions of the macros implemented in TephraProb were applied by Biass et al. (2014) and Scaini et al. (2014) to assess the hazard and the related impact resulting from the eruption of critical Icelandic volcanoes on the European air traffic using the Eulerian model FALL3D (Costa et al. 2006; Folch et al. 2009). Similarly, these probabilistic approaches have already been successfully applied by Bonadonna et al. (2002a); Scollo et al. (2007); Jenkins et al. (2012) and Jenkins et al. (2014) based on both the HAZMAP and ASHFALL models (Armienti et al. 1988; Hurst and Turner 1999; Hurst 1994; Macedonio et al. 2005,1988).

The user manual of TephraProb, submitted as a supplementary file, provides in–depth technical descriptions of all the functions of the package. Each function is thoroughly commented in order to allow customization for more advanced users. In particular, functions to calculate the MER (Eq. 4; Degruyter and Bonadonna 2012) or to calculate the mass of a thermal (Eq. 6; Bonadonna et al. 2002a) contain guidance for modifying all empirical parameters. In addition, all users are free to modify the code for their own needs following the terms and conditions of the GNU GPL3 license.

Intentionally, no default scenario is provided in order to not promote TephraProb as a black box. When volcanoes with little documented history are considered, we encourage the user to rely on a global understanding of eruption processes combined with the use of global databases (Crosweller et al. 2012; Siebert et al. 2010; Simkin and Siebert 1994) to identify analogue eruptions. When eruption scenarios are based upon detailed field studies, the TError code (Biass et al. 2014) can serve as a systematic tool to identify the probability distributions of ESPs. In either case, the user manual of TephraProb provides a list of empirical parameters for Tephra2 for eruptions characterized by the inversion technique of Connor and Connor (2006).

This version of TephraProb is published on GitHub as a basis for further development based on inputs from the scientific community. Identified directios of development include i) the implementation of a variety of models available in the literature to quantify ESPs such as Sparks (1986); Wilson and Walker (1987); Mastin et al. (2009); Woodhouse et al. (2013) and Mastin (2014), ii) probabilistic inversion schemes to systematically assess the likelihood of past events (Elissondo et al. 2016) and iii) better integration of Reanalysis datasets to calculate the required atmospheric parameters (Degruyter and Bonadonna 2012).