Skip to main content

Probabilistic Volcanic Ash Hazard Analysis (PVAHA) I: development of the VAPAH tool for emulating multi-scale volcanic ash fall analysis


Significant advances have been made in recent years in probabilistic analysis of geological hazards. Analyses of this kind are concerned with producing estimates of the probability of occurrence of a hazard at a site given the location, magnitude, and frequency of hazardous events around that site; in particular Probabilistic Seismic Hazard Analysis (PSHA). PSHA is a method for assessing and expressing the probability of earthquake hazard for a site of interest, at multiple spatial scales, in terms of probability of exceeding certain ground motion intensities. Probabilistic methods for multi-scale volcanic ash hazard assessment are less developed. The modelling framework presented here, Probabilistic Volcanic Ash Hazard Analysis (PVAHA), adapts the seismologically based PSHA technique for volcanic ash. PVAHA considers a magnitude-frequency distribution of eruptions and associated volcanic ash load attenuation relationships and integrates across all possible events to arrive at an annual exceedance probability for each site across a region of interest. The development and implementation of the Volcanic Ash Probabilistic Assessment tool for Hazard (VAPAH), as a mechanism for facilitating multi-scale PVAHA, is also introduced. VAPAH outputs are aggregated to generate maps that visualise the expected volcanic ash hazard for sites across a region at timeframes of interest and disaggregated to determine the causal factors which dominate volcanic ash hazard at individual sites. VAPAH can be used to identify priority areas for more detailed PVAHA or local scale ash dispersal modelling that can be used to inform disaster risk reduction efforts.


Numerous approaches have been adopted in the past to assess volcanic ash hazard at the local-scale (10s km) including observational, statistical, deterministic and probabilistic techniques (Bonadonna et al. 2002a; Bonadonna et al. 2002b; Blong 2003; Bonadonna and Houghton 2005; Costa et al. 2006; Magill et al. 2006; Jenkins et al. 2008; Costa et al. 2009; Folch et al. 2009; Folch and Sulpizio 2010; Simpson et al. 2011; Bear-Crozier et al. 2012; Jenkins et al. 2012a; Jenkins et al. 2012b). However, incomplete historical data on the magnitude and frequency of eruptions worldwide and the difficulties associated with up-scaling computationally intensive volcanic ash dispersal models have limited regional or global scale assessments (100 s km).

Simple assessments of volcanic ash hazard are based on compiling observations of the distribution of volcanic ash from historical eruptions, an approach that is still adopted worldwide (McKee et al. 1985; Barberi et al. 1990; Bonadonna et al. 1998; Costa et al. 2009). These maps discriminate land areas buried by volcanic ash fallout in the past from those that have not. Deterministic methods extend the usefulness of observational methods by utilising the benefits of numerical and computational models and typically consider the causes driving the hazard. The advantage of this approach is that it is computationally straightforward and provides a conservative result, which can be used to maximise safety. The disadvantage is that subjective and implicit assumptions made on the probability of the chosen scenario commonly result in an overestimation of conservative hazard values, whereby the largest possible eruption may be possible but highly unlikely.

Probabilistic methods estimate the probability of occurrence of the hazard at a site according to the location, magnitude and frequency of occurrence of hazardous events around that site. They are flexible and can take into account as much data as you have available. Probabilistic methods produce hazard curves, which provide information on the level of expected hazard for any given timeframe. Incorporating occurrence rate information into hazard analysis is more complex than deterministic, statistical or observational approaches. However, the resulting hazard curve is more useful for prioritising regions where more detailed analysis is needed. Probabilistic approaches to volcanic ash hazard assessment are the focus of this study.

Previous work

In the past, probabilistic analyses of volcanic ash hazard have focused on quantitative assessments of the frequency and potential consequences of eruptions. Simpson et al. (2011) undertook a quantitative assessment of volcanic ash hazard across the Asia-Pacific region using the Smithsonian Institution’s Global Volcanism Program (GVP) database that enabled the straightforward production of magnitude–frequency plots for each country and, to some extent, provinces within countries, in the region. Quantitative approaches have focused on a single source or site of interest at the local scale (10s of km) using tephra dispersal models (e.g. Campi Flegrei, Italy (Costa et al. 2009); Gunung Gede, Indonesia (Bear-Crozier et al. 2012); Okataina, New Zealand (Jenkins et al. 2008); Somma-Vesuvio, Italy (Folch and Sulpizio 2010) and Tarawera, New Zealand (Bonadonna and Houghton 2005). Numerical simulations of volcanic ash fallout generally involve running a deterministic eruptive scenario that represents the most likely event (based on historical investigation and/or modern analogues) over a period of time sufficiently large as to capture all possible meteorological conditions (Magill et al. 2006; Folch et al. 2008a; Folch et al. 2008b; Folch and Sulpizio 2010; Bear-Crozier et al. 2012)

Regional-scale probabilistic volcanic ash hazard assessments are less common (Yokoyama et al. 1984; Hoblitt et al. 1987; Hurst 1994; Hurst and Turner 1999; Magill et al. 2006; Ewert 2007). (Jenkins et al. 2012a; Jenkins et al. 2012b) employed a stochastic simulation technique that up-scales implementation of the ash dispersal model ASHFALL for regional-scale assessments (Hurst 1994; Hurst and Turner 1999). This approach presented a method for assessing regional-scale ash fall hazard, which had not been attempted previously and represents an important step forward in the development of techniques of this kind. However, limitations associated with up scaling conventional ash dispersal modelling methods include the computationally intensive nature of regional-scale applications that would require significant high-performance computing resources, long simulation times and could potentially constrain the spatial resolution, geographic extent and number of sources considered. Hazard curves of annual exceedance probability versus volcanic ash hazard for individual sites of interest are typically not generated and therefore disaggregating the dominant contribution to the hazard at particular sites of interest by magnitude, source or distance is not captured.

Motivation for the current work

Workers in other geohazards fields (earthquake, wind, flood etc.) have faced similar limitations associated with quantifying hazard on the regional-scale. Major developments were made by seismologists working in this space in the 1960’s, with a view to assessing ground motion hazard at multiple sites associated with potential earthquake activity (Cornell 1968). A methodology was developed for quantifying earthquake hazard at the regional-scale named Probabilistic Seismic Hazard Analysis (PSHA; Cornell 1968; McGuire 1995, 2008). PSHA consists of a four-step framework for which uncertainty in size, location and likelihood of plausible earthquakes can be incorporated to model the potential impact of future events (Robinson et al. 2006). This methodology has since been adapted for tsunami (Lin and Tung 1982; Rikitake and Aida 1988; Geist and Parsons 2006; Thio et al. 2007; Thomas and Burbidge 2009; Sørensen et al. 2012; Power et al. 2013) and applied to regional–scale tsunami hazard assessments (e.g. Indonesia; Horspool et al. 2014). Early attempts at partially adapting PSHA to volcanic ash (on a local scale) were reported by Stirling and Wilson (2002) for two volcanic complexes on the North Island of New Zealand (Okataina and Taupo). This study seeks to further advance the adaptation of PSHA for volcanic ash hazard at regional spatial scales.

Probabilistic Seismic Hazard Assessment (PSHA)

For decades seismologists have been faced with the inherent difficulty of assessing and expressing the probability of earthquake hazard at a site of interest in terms of maximum credible intensity (Cornell 1968; McGuire 1995). PSHA was developed to consider a multitude of earthquake occurrences and ground motions and produce an integrated description of seismic hazard representing all events (Cornell 1968; McGuire 1995). It is derived from the early formulation of seismic hazard analysis by Cornell (1968) and Esteva (1968). Initially, PSHA was developed to assess seismic risk at individual sites however over time the methodology was applied systematically to a grid of points yielding a regional seismic probability map with contours of maximum ground motion of equal timeframe (Cornell 1968; McGuire 1995). Traditional PSHA considers the contribution of magnitude and distance to the hazard and selects the most likely combination of these to accurately replicate the uniform hazard spectrum (McGuire 1995). Advances in seismic hazard analysis and the proliferation of high-performance computing have led to the development of event-based PSHA. Event-based Probabilistic Seismic Hazard Analysis allows calculation of ground-motion fields from stochastic event sets.

The four-step procedure for event-based PSHA reported by Musson (2000) is summarised below and presented in Fig. 1:

Fig. 1
figure 1

Schematic representation of the four stage procedure for PSHA: 1. Source; 2. Recurrence; 3. Ground Motion and; 4. Probability of exceedance (modified after TERA (1980) and Musson (2000))

  1. 1.

    Seismicity data for the region of interest must be spatially disaggregated into discrete seismic sources.

  2. 2.

    For each seismic source the seismicity is characterised with respect to time (i.e. the annual rate of occurrence of different magnitudes).

  3. 3.

    A stochastic event set is developed which represents the potential realisation of seismicity over time and a realisation of the geographic distribution of ground motion is computed for each event (taking into account the aleatory uncertainties in the ground motion model).

  4. 4.

    This database of ground-motion fields, representative of the possible shaking scenarios that the investigated area can experience over a user-specified time span, are used to compute the corresponding hazard curve for each site. Hazard curves are computed for each event individually and aggregated to form probabilistic estimates.

This paper presents a methodology developed at Geoscience Australia, which modifies the four-step procedure of PSHA for volcanic ash hazard analysis at a regional-scale. The framework named here, Probabilistic Volcanic Ash Hazard Analysis (PVAHA) considers the magnitude-frequency distribution of eruptions and associated volcanic ash load attenuation relationships and produces an integrated description of volcanic ash hazard for all events across a region of interest. An algorithm was developed to facilitate a PVAHA named here, the Volcanic Ash Probabilistic Assessment tool for Hazard (VAPAH). This approach builds on the previous work of Stirling and Wilson (2002), Simpson et al. (2011) and Jenkins et al. (2012a; 2012b) towards the development of tools and techniques for conducting regional-scale probabilistic volcanic ash hazard assessment.

An assessment for the Asia-Pacific region was undertaken during the development of the VAPAH algorithm. The reader is referred to the companion paper (Miller et al. 2016) for a detailed workflow and discussion of the Asia-Pacific region case study. This paper focuses on the adaptation of PSHA for volcanic ash, the development of the PVAHA framework and the VAPAH algorithm itself. Results from sub-regions of the Asia-Pacific study are only included here as needed to illustrate concepts and to describe the advantages and disadvantages of the overall approach. This manuscript is divided into four sections:

  1. 1.

    A description of the proposed framework for PVAHA.

  2. 2.

    The procedure used for identification of source volcanoes, development of eruption statistics and calculation of magnitude frequency relationships for each source.

  3. 3.

    Derivation and validation of ash load prediction equations derived from volcanic ash dispersal modelling used to inform the PVAHA.

  4. 4.

    Development of the VAPAH algorithm.


A framework for Probabilistic Volcanic Ash Hazard Analysis (PVAHA)

A probabilistic framework for assessing volcanic ash hazard at multiple spatial scales (PVAHA) adapted from PSHA is presented here (Fig. 2). The modified four-step procedure is outlined below:

Fig. 2
figure 2

Schematic representation of the modified PSHA procedure for Probabilistic Volcanic Ash Hazard Assessment (PVAHA); 1. Sources; 2. Magnitude-frequency relationships; 3. Ash load attenuation with distance and; 4. Probability of exceedance

  1. 1.

    Volcanic sources with respect to any given site of interest must be identified.

  2. 2.

    For each volcanic source the annual eruption probability must be calculated based on magnitude-frequency relationships of past events.

  3. 3.

    For a set of stochastic events (synthetic catalogue) volcanic ash load attenuation relationships must be calculated (derived from conventional ash dispersal modelling).

  4. 4.

    Calculation of the annual exceedance probability versus volcanic ash hazard for each stochastic event at each site across a region of interest.

Database development and data completeness

Following the procedure analogous to Jenkins et al. (2012a) a database of volcanic sources and events for the region of interest was prepared using entries from the GVP catalogue. The Smithsonian Institution’s Global Volcanism Program (GVP) catalogue of Holocene events was used to identify volcanic sources for analysis (Siebert et al. 2010). The GVP reports on current eruptions from active volcanoes around the world and maintains a database repository on historical eruptions over the past 10,000 years. We acknowledge this database is not a complete record and does contain gaps in the eruption record. Factors that contribute to these gaps, particularly in data-sparse regions like the Asia-Pacific include incomplete or non-existent historical records, poor preservation of deposits or lack of accessibility to geographically remote sources. However, the GVP database is widely recognised as the most complete global resource currently available and represents the authoritative source for information of this kind.

Other sources of data can be used to augment an analysis of this kind including volcano observatory archive data and other databases including but not limited to the Large Magnitude Explosive Volcanic Eruptions database (LaMEVE; Crosweller et al. 2012). The procedure for creating a database of volcanic sources and events for a PVAHA is described below using the Asia-Pacific examples to provide context where needed (Miller et al., 2016). Database fields including volcano ID, region, sub-region, volcano type, volcano name, latitude, longitude, eruption year and Volcano Explosivity Index (VEI; Newhall and Self 1982) were captured for each eruption at each volcano in a region of interest. These entries were further examined and volcanic sources classified as submarine, hydrothermal, and fumerolic or of unknown type were discarded.

Before calculating magnitude frequency relationships each source in the database, the record must be assessed for completeness (Simpson et al. 2011). The eruption record consists of all known events for each volcano in the database. Different magnitude eruptions have different time periods for which the record is considered complete, and these periods may vary significantly across a region (Simpson et al. 2011). Additionally, larger eruptions are better preserved in the record than smaller eruptions and this has important implications for data completeness (Jenkins et al. 2012a). With this in mind, events in the database are grouped into sub-regions defined by geographic boundaries already adopted by the GVP catalogue for consistency and further subdivided into magnitude classes, VEI 2–3 for smaller magnitude events and VEI 4–7 for larger magnitude events following the methodology of Jenkins et al. (2012a). Jenkins et al. (2012a) approach for calculating individual record of completeness for each magnitude class was based on reporting by Simkin and Siebert (1994) who declared smaller magnitude eruptions globally complete from the 1960’s and larger magnitude eruptions globally complete over the last century.

Sources with no assigned VEI but designated caldera ‘C’ or Plinian ‘P’ are allocated to the larger magnitude class as arbitrary VEI 4 events. This does not include all the remaining caldera and Plinian eruptions in the database, which were assigned a specific VEI (typically in range VEI 4–7). It’s important to note here that only those events classified ‘C’ or ‘P’ with no assigned VEI were arbitrarily allocated to VEI 4. We acknowledge that Plinian style, caldera-forming eruptions are commonly associated with eruptions greater than VEI 4 however VEI 4 is selected as a minimum magnitude representing a conservative estimate for the magnitude for the smaller number of these events given the absence of further information (only 13 events identified where this is the case for the Asia-Pacific region). The record of completeness (ROC) was then assessed for each sub-region using the ‘break in slope’ method by plotting the cumulative number of eruptions against time for each magnitude class for each sub region (Simpson et al. 2011). Similar to Jenkins et al. (2012a) completeness was identified by a linear increase in the cumulative number of eruptions per unit time. The reader is referred to Miller et al. (2016) for the individual ROC values for the Asia-Pacific. An example is provided here for the Indonesia sub-region (Fig. 3).

Fig. 3
figure 3

Example Record of Completeness plots for events in the sub-region of Indonesia; (a) Record of completeness (red line) for small magnitude eruptions; (b) Record of completeness (red line) for large magnitude eruptions

Magnitude-frequency relationships

Having established the record of completeness for each source in the database a procedure analogous to developing earthquake magnitude-frequency distributions for PSHA is adopted here for assessing the annual rate of occurrence for eruptions of different magnitudes at each source (Musson 2000). Where traditional probabilistic techniques focus on a single volcano for which the hazard is estimated independent of the probability of the eruption occurring, the framework reported here is based on the premise that the ash fall hazard associated with a given site may represent a maximum expected hazard from multiple sources. By extension of this, each source is likely to have varying eruption probabilities, styles and magnitudes and therefore traditional approaches must be modified to accommodate for this heterogeneity (Connor et al. 2001; Bonadonna and Houghton 2005; Jenkins et al. 2008; Jenkins et al. 2012a). In order to calculate the annual eruption probability for each volcanic source at each magnitude the probability of an event of any magnitude occurring and the conditional probability of an event of a particular magnitude occurring must be calculated first.

Probability of an event of any magnitude

Firstly, the annual eruption probability for each volcanic source (λ) must be determined by dividing the total number of events (N) by the time period for which the catalogue is thought to be complete (T):

$$ \uplambda = \mathrm{N}/\mathrm{T} $$

Equation 1 must be solved for each volcanic source in the small and large magnitude classes separately using the associated record of completeness values calculated in the previous section (i.e. λ (VEI 2–3) and λ (VEI 4–7); Table 1). In order to arrive at the likelihood of an event of ‘any magnitude’ occurring (i.e. λ (VEI 2–7) at each volcanic source (analogous to PSHA) the total number of events (N) and ROC (T) for each magnitude class in a sub-region must be aggregated into a single value for each. This is achieved by normalising the occurrence interval for the small magnitude class (240 years) to the occurrence interval of the large magnitude class (420 years), with an assumed constant eruption rate. The conversion factor for the number of events is calculated by taking the ratio of the large to the small magnitude completeness periods (e.g. Indonesia; 420:240 = 1.75). The normalised number of small magnitude events ‘N (VEI 2–3) normalised’ is then calculated for each source by multiplying the conversion factor (e.g. Indonesia = 1.75) by the number of small magnitude events ‘N (VEI 2–3; Table 1). The final value for ‘N (VEI 2–7) for each source is calculated by adding the number of small magnitude events ‘N (VEI 2–3) normalised’ to the number of large magnitude events ‘N (VEI 4–7; Table 1)’. Sources with no events during the time period for which their record is deemed complete (e.g. Besar) are assigned records from analogous volcanoes (those of the same type category) following the method of Jenkins et al. (2012a) in order to provide some insight on likely eruption behaviour in the absence of empirical data. Equation 1 can then be solved for the annual eruption probability of each source at any magnitude ‘λ(VEI 2–7)’ using the calculated values for ‘N (VEI 2–7)’ and ‘T’ (e.g. Indonesia = 420 years).

Table 1 Annual eruption probability ‘λ (VEI 2–7)’ for each volcanic source in the Indonesian sub-region using a conversion factor of 1.75

Conditional probability of an event of a particular magnitude

In the previous step, the probability of an event of any magnitude occurring λ(VEI 2–7) at each source was ascertained. The next step is to ascertain the conditional probability of an event being a ‘particular’ magnitude (e.g. VEI 2, 3, 4, 5, 6 or 7). This calculation is performed using the database of events and a classification scheme for volcano morphology (shape) used by Jenkins et al. (2012a). Firstly, all source volcanoes in the database (e.g. Asia-Pacific) are assigned a type category related to their morphology and previous eruption style. Five type categories are considered; lava dome, small cone, large cone, shield and caldera; Table 1). The conditional probability of an event at each magnitude is equal to the number of events for a type category at a particular magnitude divided by the total number of events in the magnitude class (i.e. small magnitude class (VEI 2–3) - Table 2; large magnitude class (VEI 4–7) - Table 3) in the database.

Table 2 Conditional probability of an event VEI 3 or less occurring for each of the five volcano types in the database (e.g. Asia-Pacific)
Table 3 Conditional probability of an event VEI 4 or greater occurring for each of the five volcano types in the database (e.g. Asia-Pacific)

Annual probability of an event

The annual probability of an event of a given magnitude for each source, needed for the PVAHA, can now be calculated by multiplying the annual probability of an event of any magnitude for a source by the probability that the event will be a particular magnitude (e.g. Indonesia; Table 4). Metadata is developed to preserve the distinction between annual probability values based on historical data versus analogues (e.g. Besar) so that uncertainty associated with these assumptions is carried through the remaining PVAHA. These magnitude-frequency relationship calculations are repeated for all sub-regions in the database.

Table 4 Annual probability of an event of a particular magnitude occurring for each volcanic source in the Indonesian sub-region

Emulating volcanic ash load attenuations relationships

Earthquake hazard is measured in terms of the level of ground motion that has a certain probability of being exceeded over a given time period (McGuire 1995, 2008). Ground motion prediction equations (GMPEs) or attenuation relationships are used to provide a means of predicting the level of ground shaking and its associated uncertainty at any given site or location (Fig. 4). GMPEs are based on an earthquake magnitude, source-to-site distance, local soil conditions and fault mechanism and are an integral part of PSHA. A process was developed here for adapting the GMPE approach used for seismic hazard to volcanic ash hazard. The process involves derivation of a mathematical expression for volcanic ash load attenuation with distance from source, for an event, using a gridded hazard footprint generated by an ash dispersal model. The resulting equation, named here an Ash Load Prediction Equation (ALPE), statistically emulates the volcanic ash attenuation relationship (Fig. 4). Where up scaling of conventional volcanic ash dispersal modelling techniques would be computationally intensive and time consuming, ALPEs can be used to emulate generalised volcanic ash hazard (derived from dispersal models) for any given event(s) at any location(s), from any volcanic source of interest as a function of distance of the site from each source.

Fig. 4
figure 4

a Example plot of predicted peak ground acceleration with changing distance using two European GMPE’s for different magnitude ranges modified after Akkar and Bommer (2007) and Bommer et al. (2007); (b) Example plot of volcanic ash load attenuation relationships using three ALPEs generated through dispersal modelling for this study (ALPE 100 - VEI 3; ALPE 150 - VEI 4; ALPE 201 - VEI 5)

The procedure for calculating ALPEs for a PVAHA is described below and includes the following:

  1. 1.

    Development of a synthetic catalogue of events

  2. 2.

    Volcanic ash dispersal modelling

  3. 3.

    Derivation of an ALPE

  4. 4.

    Validation of an ALPE

Development of synthetic catalogue of events

Similar to approach taken for PSHA, a synthetic catalogue of events is developed as a basis for the generation of the ALPEs (one ALPE per event) needed for the PVAHA. A relationship for the rate of volcanic ash load decay with distance from the source, as a function of magnitude, column height, duration, wind turbulence, direction and speed, must be established for each event of interest. The dispersal of volcanic ash through the atmosphere produces deposits at ground level that diminish gradually in load (kg/m2) with distance from the source but in directions controlled by the wind. Consequently, ash load attenuation is a complex function of distance and azimuth from source (Stirling and Wilson 2002). Synthetic events developed here were based on the development of a logic-tree data structure (Fig. 5). The purpose of this structure was to capture all possible variations in volcanological conditions and to quantify the uncertainty associated with the inputs for each event (Bommer and Scherbaum 2008). The influences of site-specific meteorological conditions are considered separately at a later stage in the procedure.

Fig. 5
figure 5

Schematic representation of logic tree structure adopted for creation of synthetic catalogue of events; Node 1: Eruption column height (1000 - 40000 m); Node 2: Eruption Duration (1–12 h); Node 3: Source term (Plinian, Sub-Plinian, Vulcanian, Strombolian or Violent Strombolian)

Fig. 6
figure 6

Schematic representation of a FALL3D volcanic ash load hazard footprint; (b) Example plot of volcanic ash load attenuation relationship derived from the ALPE extracted from the generic hazard footprint in (a)

A simplified, schematic representation of the logic tree data structure used is presented in Fig. 5.

Input parameters included:

  • eruption column height (in meters; between 1 000 and 40 000),

  • eruption duration (in hours between 1 and 12) and;

  • eruption style (Strombolian, Vulcanian, Sub-Plinian and Plinian).

A total of 1056 events were developed and assigned an equal weighting for probability of occurrence. Events are not volcano specific but rather represent a suite of synthetic eruptions, which when coupled with magnitude-frequency statistics and prevailing meteorological conditions for a region of interest, can be used to assess a range of potential events at any volcanic source. It is important to note that assignment of weightings can be modified where prevalent eruption behaviour for a region of interest is well known.

Volcanic ash dispersal modelling

The volcanic ash dispersal model FALL3D was used here to computationally model volcanic ash fall hazard footprints needed for the calculation of ALPEs. FALL3D is a time-dependant Eulerian model that solves the advection–diffusion-sedimentation (ADS) equation on a structured terrain-following mesh. FALL3D outputs time-dependant deposit load at ground level as a hazard footprint of changing ash load with distance from source (Folch et al. 2012). It is acknowledged that FALL3D is one of a number of suitable ash dispersal models that could be utilised for an analysis of this kind.

An assessment of the mass eruption rate, height and shape of the eruption column is made for each stochastic event in the catalogue. These parameters together describe the eruptive source term needed to simulate the dispersal of volcanic ash using FALL3D. The source term can be defined as either a 1-D buoyant plume model (Bursik 2001) or as an empirical relationship (Suzuki 1983). The empirical relationship (Suzuki, 1983) used here estimates the mass eruption rate (MER) given an eruption column height (H) using known best-fit relationships of MER versus H (Sparks et al. 1997). A generalised total grainsize distribution (TGSD) is used to account for a range of eruption potential eruption styles which includes minimum and maximum grainsize (phi), average grainsize (phi), sorting, density range (kg/cm3) and sphericity of clasts. FALL3D was used to simulate 1056 events in the synthetic catalogue.

Derivation of ash load prediction equations

A script was developed for extracting the volcanic ash attenuation relationship (changing ash load with distance) for each hazard footprint generated by the dispersal model in the synthetic catalogue (Fig. 6). Each ALPE represents a single event (1056 total). When coupled with magnitude-frequency statistics and prevailing meteorological conditions for a region of interest, each ALPE can be used to statistically emulate the expected volcanic ash hazard from an event of this kind at any location of interest from any volcanic source as a function of distance of the site from the source. Not unlike GMPEs used to conduct PSHA, the generation and application of ALPEs will have a considerable influence on the outcome of the PVAHA. The ALPEs developed here use the dispersal model FALL3D however other dispersal models could be used (e.g. ASHFALL (Hurst 1994; Hurst and Turner 1999), HAZMAP (Costa et al. 2009), or TEPHRA (Bonadonna and Houghton 2005)) and the authors would encourage the development of ALPEs using a range of dispersal models currently available to build on and compare/contrast with the current work.


In order to determine the uncertainty of the ALPEs or the degree to which they accurately reproduce simulated ash fallout generated by FALL3D and observed deposit data gathered from field studies of historical eruptions, a validation is presented for the F2 Plinian fall deposit generated by the 1815 eruption of Tambora, on the island of Sumbawa, Indonesia. FALL3D has already been widely validated against several tephra deposits and airborne ash cloud observations from different eruptions (Costa et al. 2006; Macedonio et al. 2008; Scollo et al. 2008a; Scollo et al. 2008b; Scollo et al. 2009; Corradini et al. 2011; Costa et al. 2012; Kandlbauer et al. 2013; Costa et al. 2014; Kandlbauer and Sparks 2014). An inversion simulation for the Tambora F2 deposit data is presented here using FALL3D by comparing observed measurements for the F2 deposit at 27 sample locations with simulated ash measurement generated by FALL3D at all corresponding sample sites (Sigurdsson and Carey, 1989). We then generate an ALPE from the FALL3D data following the methodology outlined above and use this equation to calculate the expected ash fallout for the F2 deposit. We present difference plots to compare observed versus simulated (FALL3D), simulated versus calculated (ALPE) and calculated versus observed estimates (Sigurdsson and Carey 1989) for volcanic ash load and comment on performance of the ALPE for emulating the Tambora F2 Plinian fall deposit.

Validation of FALL3D for the Tambora F2 Plinian fall deposit

To inverse model the Tambora F2 Plinian fall deposit, the validation method of Costa et al. (2012) is adopted. Costa et al. (2012) constrained the eruption dynamics and ash dispersal characteristics associated with the Campanian Ignimbrite (39 ka) eruption in Italy and later the Youngest Toba Tuff (75 ka) super-eruption in Indonesia. This approach combines time-dependant meteorological fields for the region, a spectrum of volcanological parameters (erupted mass, mass eruption rate (MER), column height and total grainsize distribution) and over 100 simulations of the ash dispersal model, FALL3D. Optimal values of the input parameters are obtained by best-fitting measured Tambora F2 tephra thicknesses over the entire dispersal area (27 locations) and minimising the deviation of regression. The explored range of input parameters for the Tambora-F2 eruption is reported in Table 5.

Table 5 The explored range and best-fit input parameters modelled dispersion of the Tambora-F2 fall deposit using FALL3D

Ten years of wind data (January 2000 to December 2010) were obtained from the National Centres for Environmental Protection (NCEP) and Atmospheric Research (NCA) global reanalysis project (Kistler et al. 2001). The NCEP/NCAR reanalysis data archive contains six hourly data at 17 pressure levels ranging from 1000 to 10 mb with a 2.5° horizontal resolution. The methodology of Costa et al. (2012) used here assumes that this collection of modern meteorological fields can statistically represent a proxy for those at the time of the Tambora F2 eruption (~200 years ago). Vertical meteorological profiles were extracted from the wind data at a gridded location closest to the source and interpolated to the FALL3D computational grid. In order to reduce the computational requirements, vertical, but not horizontal, changes in wind conditions with distance from the source were accounted for here at six hourly intervals using linear temporal and spatial interpolations. The computational domain was discretised by a horizontal grid step of Δx = Δy = 1.35 km and a vertical step of Δz = 1 km. The computational domain extended from 9° S to 7° N and from 117° E to 118° W.

The distribution of mass within the column was calculated using an empirical parameterisation based on that of Suzuki (1983) and Pfeiffer et al. (2005). In order to account for aggregation processes (clustering of fine ash particles) an aggregation model similar to that of Cornell et al. (1983) and used by Costa et al. (2012; 2014) was adopted here. This aggregation model assumes that 50 % of the 63–44 μm (4–4.5Ø) ash, 75 % of the 44–31 μm (4.5-5Ø) ash and 95 % of the sub-31 μm (<5Ø) ash fell as aggregated particles. The diameter and density of the aggregates were determined by best fit in the simulations. An approach developed for best fitting the spatial variation between recorded and simulated tsunami heights, the Aida indices Aida (1978) was adopted by Costa et al. (2014) for tephra and is similarly used here to measure the reliability of modelled results. The Aida index K represents the geometric average of the distribution and the second index k is the associated standard deviation of the distribution:

$$ \log\ \mathrm{K}=1\mathrm{n}\mathrm{i}=1\mathrm{n} \log\ {\mathrm{K}}_{\mathrm{i}} $$
$$ \log\ \mathrm{k} = {\left[1\mathrm{n}\mathrm{n}=1\mathrm{n}{\left( \log\ \mathrm{K}\mathrm{i}\ \right)}^2\hbox{-} {\left( \log\ \mathrm{K}\ \right)}^2\right]}^{\frac{1}{2}} $$

where n is the total number of measurements and Ki = Mi/Hi is the ratio of measured thickness (load) at Mi,,the i-the location and Hi is the simulated thickness (load) at the same site. In keeping with the approach for tsunami and Costa et al. (2014) we consider the simulated tephra thickness results satisfactory when:

$$ 0.95<\mathrm{K}>1.05\ \mathrm{and}\ \mathrm{k}<1.45 $$

The best-fit results from FALL3D are reported in Table 5 and indicate that the MER was ~1.2 x 108 kg/s, the eruption column height was ~33 km and the total mass deposits as fallout was ~ 1.2 x 1014 kg. These results are in good agreement with estimates made by Sigurdsson and Carey (1989) based on field observations. The corresponding simulated, best-fit, deposit thicknesses are depicted in a difference plot and reported in Fig. 7. The correlation coefficient between log (measured thickness) and log (simulated thickness) is 0.87 for the best meteorological fit. All simulated thicknesses are between 1/5 and 5 times the observed thicknesses and the reliability of the best-fit results are further emphasised by the Aida index values; reflecting a geometric average, K = 0.98 and a geometric standard deviation, k = 1.29.

Fig. 7
figure 7

a Comparison between the thicknesses (cm) from best fit FALL3D simulations and field data for the F2 Plinian fall event, Tambora at each of the 27 sample points (modified after Sigurdsson and Carey, 1989); (b) Difference plot for best fit FALL3D simulation against observed thicknesses converted to load (kg/m2; Sigurdsson and Carey, 1989). The solid line represents a perfect agreement and the dotted and dashed black lines mark the region that is different from the observed by a factor of 10 (1/10) and 5 (1/5) respectively; (c) Difference plot for calculated ash load (ALPE) against simulated load (FALL3D) (d) Difference plot for calculated load (ALPE) against observed load (kg/m2)

Validation of an ALPE for the Tambora F2 Plinian fall deposit

Following the procedure outlined above the volcanic ash attenuation relationship for the best-fit FALL3D simulation was considered and an ALPE was derived for this event. For the purposes of this validation all volcanic ash thicknesses were converted to load (kg/m2) the preferred unit of measurement for PVAHA. Validation of the ALPE was a two-step process involving first, verification that the ALPE could statistically emulate the simulated deposit load generated by FALL3D in previous step and secondly that when compared with the observed data of Sigurdsson and Carey (1989) for this event, ash load values were in good agreement with field measurements at the 27 sample localities. Following the methodology of Costa et al. (2012) the modelled results are considered to be in good agreement with the measured observations when they are between 1/5 and 5 times the observed thickness (Fig. 7). A difference plot indicates good agreement (within 1/5 to 5 times), between the simulated best-fit (FALL3D) load (kg/m2) and the observed load (converted from thickness; (Sigurdsson and Carey 1989) plotted here at 27 sample localities (Fig. 7b). To verify that the calculated load from the ALPE closely approximates the simulated load from FALL3D (from which it was derived), a difference plot was generated. As expected, the calculated load and the simulated best-fit load are in good agreement (Fig. 7c). Finally, to complete the validation process the calculated load must be compared with the observed field data. A difference plot is generated and the calculated load (ALPES) and the observed load are found to be in good agreement (Fig. 7d).


Volcanic Ash Probabilistic Assessment tool for Hazard (VAPAH)

An algorithm was developed to facilitate the fourth step of the PVAHA framework named here, the Volcanic Ash Probabilistic Assessment tool for Hazard (VAPAH). VAPAH utilises a scripted interface and high performance computing technology in order to undertake assessments at multiple spatial scales. The VAPAH algorithm reads in magnitude-frequency relationships, an ALPE catalogue and global scale meteorological conditions for a region of interest and integrates across all possible events to arrive at a preliminary annual exceedance probability for each site across the region of interest. Other algorithms of this kind have been developed for probabilistic earthquake hazard assessment (e.g. the Earthquake Risk Model (EQRM; (Robinson et al., 2005) however this algorithm is the first of its kind specifically designed for volcanic ash hazard. Inputs for the VAPAH algorithm include:

  1. 1.

    Identification of volcanic sources for analysis.

  2. 2.

    Characterisation of magnitude-frequency relationships for each volcanic source.

  3. 3.

    Characterisation of the volcanic ash load attenuation relationship (ALPE catalogue).

  4. 4.

    A spatial grid of pre-determined resolution clipped to the domain extent (default = auto-generated).

  5. 5.

    Characterisation of meteorological conditions - prevailing wind direction (degrees) and wind speed (m/s)

Identification of volcanic sources, characterisation of magnitude-frequency relationships for each source and development of an ALPE catalogue have all been described previously. Characterisation of the meteorological conditions for the region of interest are discussed below using the Indonesian sub-region as an example.

Characterisation of meteorological conditions

The VAPAH algorithm requires an estimation of the prevailing wind direction (degrees) and speed (m/s) at a single pressure level for each source. Meteorological data can be sourced from direct observations (e.g. weather balloons, anemometers and wind vanes) or modelled data available at multiple spatial scales depending on the purpose of the PVAHA (e.g. National Centres for Environmental Protection and Atmospheric Research re-analysis (NCEP/NCAR), Weather Research and Forecasting Model (WRF) and the Australian Community Climate and Earth System Simulator ACCESS). These estimates for wind direction and wind speed are considered to be highly prevalent at each location by the algorithm and are assigned a probability weighting according to a Gaussian distribution with a standard deviation assigned by the user. An example is provided below for the deriving these variables for the Indonesia sub-region using one potential source of meteorological data.

Sixty-four years of meteorological data (January 1950 - December 2014) was sourced from the NCEP/NCAR re-analysis for the Indonesia sub-region, available at grid intervals of 2.5° globally. Among other variables, wind direction and wind speed vector components are available for 17 pressure levels to a height of 40 km. Monthly mean vector components for meridional wind (u-component) and zonal wind (z-component) were extracted at 16 locations across Indonesia, from NCEP grid points closest to each volcanic source at the 250mb pressure level (Tropopause) for a 64 year period (Fig. 8). Monthly mean wind direction (degrees) and wind speed (m/s) were derived from the u and v wind components and aggregated first for each year and then for the 64 year period using the freeware meteorological analysis and plotting tool, WRPLOT (Table 6). Prevailing wind direction and wind speed were assigned to each source from the closest NCEP point.

Fig. 8
figure 8

Monthly mean wind direction (deg) and wind speed (m/s) for the 250mb pressure level (Tropopause) aggregated for a 64 year period (1950–2014) for four of the 16 NCEP grid points used for the Indonesian sub-region. Wind rose diagrams depict the frequency of winds blowing ‘from’ a particular direction over the 64-year period. The length of each spoke is related to the frequency that the wind blows from a particular direction (e.g. N, S E or W) over the 64-year period and each concentric circle represents a different frequency (e.g. 10 %, 20 %, 30 %) emanating from zero at the centre

Table 6 Monthly mean wind direction and wind speed aggregated for a 64-year period (1950–2014) for 16 NCEP grid points across the Indonesian sub-region

VAPAH algorithm procedure

The operational procedure for the VAPAH algorithm is presented in Fig. 9. A configuration file is used to customise the extent of the assessment (Attachment 1). The configuration file reads in a series of CSV (comma separated value) files including:

Fig. 9
figure 9

Schematic representation of the operational procedure for the VAPAH tool. * VAPAH will auto-generate a spatial grid using user-specified parameters in the configuration file if one is not provided as a pre-processed input

  1. 1)

    the ALPE catalogue,

  2. 2)

    the volcano sources (included prevailing wind speed and direction)

  3. 3)

    the sites of interest (pre-processed spatial grid)

The user can then customise the run further through the configuration file. If no sites file is provided by the user, geographic coordinates for the region of interest and the resolution can be input for auto-generation of the spatial grid by VAPAH. The user must define the timeframes of interest (e.g. 100, 500, 1000 etc.) in the configuration file. Mean wind direction and mean wind speed for each source are pre-configured (see previous section) however wind direction and wind speed probabilities vary seasonally, particularly in equatorial regions. The configuration file captures this uncertainty through the wind direction distribution parameter (e.g. normal or bimodal), the wind direction distribution standard deviation parameter (degrees), the wind speed distribution parameter (e.g. normal or bimodal distribution), the wind speed distribution standard deviation parameter (m/s) and the number of wind speeds. Finally, hazard thresholds such as maximum distance from source, maximum ash load value at source and minimum sum of ash load needed to generate hazard curves or histograms can be set here by the user.

The VAPAH algorithm can be run in serial or parallel computing environments (i.e. on one or many processors simultaneously) but is optimised for high performance computing platforms utilising thousands of CPUs. Simulation time will vary according to the number of events, the number of sources, the resolution of the hazard grid and the distribution and standard deviation of meteorological conditions. A script is used to execute the procedure as follows:

  1. 1.

    The first site is located on the hazard grid and the distance is calculated in kilometres between this site and the first source in the catalogue.

  2. 2.

    The distance value is used to evaluate the first ALPE for the first synthetic event in the catalogue in order to derive the expected ash load (kg/m2) at that site for the first event.

  3. 3.

    The calculated ash load for the first event and its associated probability derived from the magnitude-frequency relationships for the first source in the catalogue are written to a results file.

  4. 4.

    The algorithm repeats (2) and (3) for each ALPE for the first source.

  5. 5.

    The algorithm then moves to the next source in the catalogue and repeats steps (2), (3) and (4) until all ALPES have been assessed for the first source for the first site.

  6. 6.

    The algorithm will then calculates the cumulative probability of each event (e.g. each instance where a volcanic ash load was recorded from one or more sources) for the first site and generates a hazard curve of volcanic ash load (kg/m2) versus annual probability of exceedance and the maximum expected ash hazard for timeframes of interest (specified by the user in the configuration file). This process will capture all instances where the first site might experience multiple ash hazards from more than one source.

  7. 7.

    The algorithm then loops to the next site and repeats steps (1) – (6) until all sites have been evaluated.

VAPAH results

VAPAH generates a database file of hazard calculations collectively referred to here as the PVAHA. Like, PSHA the PVAHA results can be post-processed using VAPAH to generate hazard curves for annual probability of exceedance at sites of interest (or all sites), maps for maximum expected ash hazard at timeframes of interest and histograms which disaggregate the hazard (e.g. location, magnitude, source etc.) for determining the primary causal factors at sites of interest. Examples of each for the Indonesia sub-region are reported in Fig. 10). Hazard curves report the annual exceedance probability versus volcanic ash load for a site of interest. These curves capture all instances where the site experienced ash loading from events originating from one or more sources (potentially thousands of events) that ash load will exceed a particular value (Fig. 10). The expected maximum ash load for timeframes of interest is also calculated. By aggregating the hazard calculations for each site, hazard maps can be generated which display the maximum expected ash load (kg/m2) at each site across the region for a timeframe of interest (e.g. 1-in-100 year event; Fig. 10). It’s important to clarify that a 1-in-100 year event does not suggest that the maximum expected ash load for a site will occur regularly every 100 years, or only once in 100 years but rather, given any 100 year period, the maximum expected ash load for a particular site may occur once, twice, more, or not at all.

Fig. 10
figure 10

a Hazard map of volcanic ash load (kg/m2) for Indonesia at 1 km resolution for the 1000 year timeframe; (b) Annual exceedance probability curve versus ash load for a site in Jakarta and (c) Histogram disaggregating the % contribution to ash load hazard of magnitude (VEI) for each source at the site in Jakarta for the 1000 year timeframe

Disaggregation of the hazard calculations can ascertain which events dominate the hazard at a particular site. The ability to disaggregate the primary causal factors contributing to the hazard at a given site (i.e. magnitude, source, distance, ash load etc.) is an inherent strength of this approach. The user can specify disaggregation parameters in the configuration file prior to undertaking an assessment and the VAPAH algorithm will generate histograms as part of the PVAHA results set (Fig. 10). Alternatively, histograms can be generated using a pre-generated database of events from assessments already computed. This functionality is useful for demonstrating what the percentage contribution to hazard of different volcanic sources on a site located in for example ‘Jakarta’ might be, and of those volcanic sources what proportion are VEI 2, 3, 4, 5, 6 and 7 (Fig. 10).


The PVAHA methodology presented here integrates across all possible events and expected volcanic ash loads to arrive at a combined probability of exceedance for a site of interest. The analysis incorporates the relative frequencies of occurrence of different events and a wide range of volcanic ash dispersal characteristics. Like PSHA, this quantitative method of estimating volcanic ash hazard has the advantage of providing consistent estimates of hazard and can be prepared for one site or many, all in the same region but in significantly different geographic orientations with respect to potential hazard sources (McGuire 1995, 2008).

The methodology is highly customisable allowing for the flexible integration of ALPEs generated using numerous ash dispersal models and eruption statistics derived from a variety of sources (Whelley et al. 2015). Capacity to build in flexibility in input assumptions highlights the power of this approach for quantifying uncertainty. The VAPAH algorithm effectively replaces the need for computationally intensive and time-consuming ash dispersal modelling that required to survey the number of events, spatial scale and resolution achievable using the technique. Similar to PSHA, the PVAHA framework proposed here is not intended to replace conventional ash modelling techniques. PVAHA considers all possible events, from all possible sources for a region of interest and produces a broad-brush, first approximation of the hazard (like PSHA) which can be updated and re-run regularly. It provides a quantitative mechanism for constraining the causal factors of ash hazard globally and can be used to underpin prioritisation of sources for further local scale dispersal modelling work.

While the PVAHA procedure, through utilisation of the VAPAH algorithm, is primarily concerned with aggregating the hazard contributions from all sources, disaggregating the volcanic ash hazard has two important implications for the usefulness of the technique over other approaches to multi-scale assessments. Firstly, the causal factors, which dominate the volcanic ash hazard for each site including magnitude, distance and source, are captured in a results file that can be easily interrogated (a limitation of dispersal modelling outputs typically generated as gridded data). Second, disaggregation can be used to identify priority areas (sites or sources) from the multitude of volcanic events, for subsequent, more detailed analysis at the local scale that could be used to inform decision-making (e.g. targeted ash modelling at Merapi for Yogyakarta hazard assessment). The benefit of disaggregating the analysis is a better overall understanding of the contributing factors to volcanic ash hazard for a region and evidence-based targeting of disaster risk reduction efforts.

Addressing uncertainty

The quality and value of the resulting assessment is controlled by the quality of the input models and it is critical that the uncertainties in parameter values, as well as those associated with the dispersal model itself are suitably accounted for. Uncertainty is addressed here through the development of a suite of ALPEs that account for the full spectrum of input parameters and the associated uncertainty of the dispersal model in use. This process also allows multiple competing hypotheses on models and parameters to be incorporated into the analysis (i.e. ALPEs based on other ash dispersal models). Integrating synthetic catalogue of events with magnitude-frequency relationships derived for each source and prevailing meteorological conditions for the region and carrying those probabilities through to the analysis outputs using the VAPAH algorithm also mitigates the uncertainty in assumptions made for annual eruption probability. Sensitivity analyses should be periodically carried out for all parameters and models and updated as new data and information become available in order to refine the resulting analysis.

Limitations, assumptions and caveats

The PVAHA methodology presented here incorporates a number of assumptions and is subject to limitations on what is produced and how the information can be interpreted and used. All assumptions are made explicit and are open to review and refinement with new evidence. Key assumptions made, limitations and caveats on the resulting assessment include the following:

  1. 1.

    Determination of the record of completeness for a sub-region is difficult to identify and can have a significant impact on the probability values derived for those sources (e.g. volcanoes with long repose periods might be under-represented).

  2. 2.

    The probability of an event is based on the type of source (e.g. caldera, large cone etc.) and this study assumes a static source type (i.e. a caldera remains a caldera) however morphologies are typically dynamic and evolve over the history of a source (e.g. large cones can become calderas). This has implications for estimating the probability of events.

  3. 3.

    All events are assumed to follow a memory-less Poisson process meaning the probability of an event occurring today is not contingent on whether or not an event occurred yesterday.

  4. 4.

    A 2 km radii is applied to each source area and volcanic ash load estimates within this zone are not utilised due to over-estimations of proximal deposits (a feature of the dispersal model and considered acceptable due to the general absence of population, buildings and infrastructure within 2 km of a source (i.e. on the edifice slopes)

  5. 5.

    This procedure and the VAPAH tool are intended to be one in a range of tools and techniques used to provide a consistent fit-for-purpose approach to hazard assessment across multiple spatial scales.

Future directions

The PVAHA methodology presented here is fully customisable and can be modified to reflect advances in our understanding of the dynamics of volcanic ash dispersal, improvement of statistical analysis techniques for historical eruptive events and ever increasing capabilities in high-performance computing. There is no limit to number of ash dispersal models which could be used to generate ALPEs for consideration and this allows multiple competing hypotheses on models and parameters to be incorporated into the analysis.

The VAPAH algorithm currently addresses hazard, however the modular nature of the tool supports a framework for risk analysis. For example, a python module for damage (i.e. building damage, infrastructure damage or agricultural crop damage), containing vulnerability functions for volcanic ash could be developed and implemented as part of the VAPAH algorithm. Vulnerability functions are defined here as the relationship between the potential damage to exposed elements (e.g. buildings, agricultural crops, critical infrastructure, airports) and the amount of ash load (Blong 1981; Casadevall et al. 1996; Blong 2003; Spence et al. 2005; Guffanti et al. 2010; Wilson et al. 2012). Through integrating the hazard module (presented here) with a damage module, the conditional probability of damage (or loss in dollars) for an exposed element could be calculated for a given threshold of volcanic ash load. The resulting damage curves could be integrated with an exposure data module (e.g. population density, building footprints and crop extents) for the region of interest and the potential impact of events could be quantified in a risk framework.


Significant advances have been made in the field of probabilistic natural hazard analysis in recent decades leading to the development of PSHA, a method for assessing and expressing the probability of earthquake hazard at a site of interest in terms of maximum credible ground-shaking intensity for timeframes of interest. The PVAHA methodology presented here modifies the four-step procedure of PSHA for volcanic ash hazard assessment at multiple spatial scales. This technique considers a magnitude-frequency distribution of eruptions, associated volcanic ash load attenuation relationships and prevailing meteorological conditions and integrates across all possible events to arrive at a combined probability of exceedance for a site of interest. The analysis incorporates the relative frequencies of occurrence of different events and a wide range of volcanic ash dispersal characteristics. This quantitative procedure can provide rigorous and consistent estimates of volcanic ash hazard across multiple spatial scales in less time and using far fewer computational resources than those needed to up-scale conventional ash dispersal modelling.

The VAPAH algorithm, developed here to facilitate this procedure calculates the probability of exceeding a given ash load for a site of interest and generates a hazard curve of annual probability of exceedance versus volcanic ash load (kg/m2). VAPAH also calculates the maximum expected ash load at a given site for timeframes of interest. The database of results obtained for a grid of sites can be aggregated to generate maps of expected volcanic ash load at different timeframes or disaggregated by event in order to determine the percentage contribution to the hazard by magnitude, distance, source etc. This has important implications for understanding the causal factors, which dominate volcanic ash hazard at a given site, and identifying priority areas for more detailed, localised modelling that can be used to inform disaster risk reduction efforts.


  • Aida I. Reliability of a tsunami source model derived from fault parameters. J Phys Earth. 1978;26(1):57–73.

  • Akkar, S, Bommer, JJ. Empirical prediction equations for peak ground velocity driven from strong motion records in Europe and the Middle East. Bulletin of Seismological Society of America. 2007;97(2):511-530.

  • Barberi F, Macedonio G, Pareshci MT, Santacroce R. Mapping tephra fallout risk: an example from Vesuvius, Italy. Nature. 1990;344:142–4.

    Article  Google Scholar 

  • Bear-Crozier A, Kartadinata N, Heriwaseso A, Nielsen O. Development of python-FALL3D: a modified procedure for modelling volcanic ash dispersal in the Asia-Pacific region. Nat Hazards. 2012;64(1):821–38.

    Article  Google Scholar 

  • Blong R (1981) Some effects of tephra falls on buildings. In: Tephra Studies. Springer, pp 405–420

  • Blong R. Building damage in Rabaul, Papua New Guinea, 1994. Bull Volcanol. 2003;65(1):43–54.

  • Bommer, JJ, Stafford, PJ, Alarcon, JE, Akkar, S. The influence of magnitude ranges on ground motion prediction. Bulletin of the Seismological Society of America. 2007;97(6):2152-2170.

  • Bommer JJ, Scherbaum F. The Use and Misuse of Logic Trees in Probabilistic Seismic Hazard Analysis. Earthquake Spectra. 2008;24(4):997–1009.

    Article  Google Scholar 

  • Bonadonna C, Houghton BF. Total grain-size distribution and volume of tephra-fall deposits. Bull Volcanol. 2005;67:441–56.

    Article  Google Scholar 

  • Bonadonna C, Ernst G, Sparks RSJ. Thickness varations and volume estimates of tephra fall deposits: the importance of particle Renolds number. J Volcanol Geotherm Res. 1998;81:173–87.

    Article  Google Scholar 

  • Bonadonna C, Macedonio G, Sparks RSJ (2002a) Numerical modelling of tephra fallout associated with dome collapses and Vulcanian explosions: application to hazard assessment on Montserrat. In: Druitt TH, Kokelaar BP (eds) The eruption of Soufrière Hills Volcano, Montserrat, from 1995 to 1999. Geological Society London Memoir, pp 517–537

  • Bonadonna C, Mayberry GC, Calder ES, Sparks RSJ, Choux C, Jackson P, et al. (2002b) Tephra fallout in the eruption of Soufrière Hills Volcano, Montserrat. In: Druitt TH, Kokelaar BP (eds) The eruption of Soufrière Hills Volcano, Montserrat, from 1995 to 1999. Geological Society London Memoir, pp 483–516

  • Bursik M. Effect of wind on the rise height of volcanic plumes. Geophys Res Lett. 2001;18:3621–4.

    Article  Google Scholar 

  • Casadevall TJ, Delos Reyes P, Schneider DJ (1996) The 1991 Pinatubo eruptions and their effects on aircraft operations. Fire and Mud: eruptions and lahars of Mount Pinatubo. University of Washington Press. Philippines:625–636

  • Connor CB, Hill BE, Winfrey B, Franklin NM, Femina PCL. Estimation of volcanic hazards from tephra fallout. Nat Hazard Rev. 2001;2(1):33–42.

    Article  Google Scholar 

  • Cornell CA. Engineering seismic risk analysis. Bull Seismol Soc Am. 1968;58(5):1583–606.

    Google Scholar 

  • Cornell W, Carey S, Sigurdsson H. Computer-simulation of transport and deposition of the Campanian Y-5 Ash. J Volcanol Geotherm Res. 1983;17(1–4):89–109.

    Article  Google Scholar 

  • Corradini S, Merucci L, Folch A. Volcanic ash cloud properties: Comparison between MODIS satellite retrievals and FALL3D transport model. IEEE Geosci Remote Sens Lett. 2011;8(2):248–52.

    Article  Google Scholar 

  • Costa A, Macedonio G, Folch A. A three-dimensional Eulerian model for transport and deposition of volcanic ashes. Earth Planet Sci Lett. 2006;241:634–47.

    Article  Google Scholar 

  • Costa A, Dell'Erba F, Di Vito MA, Isaia R, Macedonio G, Orsi G, et al. Tephra fallout hazard assessment at Campi Flegrei caldera (Italy). Bull Volcanol. 2009;71:259–73.

    Article  Google Scholar 

  • Costa A, Folch A, Macedonio G, Giaccio B, Isaia R, Smith V (2012) Quantifying volcanic ash dispersal and impact of the Campanian Ignimbrite super‐eruption. Geophys Res Lett 39 (10)

  • Costa A, Smith VC, Macedonio G, Matthews NE. The magnitude and impact of the Youngest Toba Tuff super-eruption. Front Earth Sci. 2014;2:16.

    Article  Google Scholar 

  • Crosweller HS, Arora B, Brown SK, Cottrell E, Deligne NI, Guerrero NO, et al. Global database on large magnitude explosive volcanic eruptions (LaMEVE). J Appl Volcanol. 2012;1(1):1–13.

    Article  Google Scholar 

  • Esteva L. Bases para la formulacion de decisiones de disen ̃o sismico. Mexico: Universidad Autonoma Nacional de Me ́xico; 1968.

    Google Scholar 

  • Ewert JW. System for ranking relative threats of US volcanoes. Nat Hazard Rev. 2007;8(4):112–24.

    Article  Google Scholar 

  • Folch A, Sulpizio R. Evaluating long-range volcanic ash hazard using supercomputing facilities: application to Somma-Vesuvius (Itay), and consequences for civil aviation over the Central Mediterranean Area. Bull Volcanol. 2010;72:1039–59.

    Article  Google Scholar 

  • Folch A, Cavazzoni C, Costa A, Macedonio G. An automatic procedure to forecast tephra fallout. J Volcanol Geotherm Res. 2008a;177:767–77.

    Article  Google Scholar 

  • Folch A, Jorba O, Viramonte J. Volcanic ash forecast - application to the May 2008 Chaiten eruption. Nat Hazards Earth Syst Sci. 2008b;8:927–40.

    Article  Google Scholar 

  • Folch A, Costa A, Macedonio G. FALL3D: A computational model for transport and deposition of volcanic ash. Comput Geosci. 2009;35:1334–42.

    Article  Google Scholar 

  • Folch A, Costa A, Basart S (2012) Validation of the FALL3D ash dispersion model using obsrvations of the 2010 Eyjafjallajokull volcanic ash clouds. Atmos Environ in press:1–19

  • Geist EL, Parsons T. Probabilistic Analysis of Tsunami Hazards*. Nat Hazards. 2006;37(3):277–314.

    Article  Google Scholar 

  • Guffanti M, Casadevall TJ, Budding K (2010) Encounters of aircraft with volcanic ash clouds; A compilation of known incidents, 1953–2009. US Geological Survey. US Numbered Series. Series Number 545.

  • Hoblitt RP, Miller CD, Scott WE (1987) Volcanic hazards with regard to siting nuclear-power plants in the Pacific Northwest. US Geological Survey. US Numbered Series. Series Number 87-297.

  • Horspool N, Pranantyo I, Griffin J, Latief H, Natawidjaja D, Kongko W, et al. A probabilistic tsunami hazard assessment for Indonesia. Nat Hazard Earth Syst Sci. 2014;14(11):3105–22.

    Article  Google Scholar 

  • Hurst AW (1994) ASHFALL - A computer program for estimating volcanic ash fallout (Report and User Guide). Institute of Geological & Nuclear Sciences science report 94 (23). pp 23.

  • Hurst AW, Turner R. Performance of the program ASHFALL for forecasting ashfall during the 1995 and 1996 eruptions of Ruapehu volcano. N Z J Geol Geophys. 1999;42:615–22.

    Article  Google Scholar 

  • Jenkins S, Magill C, McAneney J, Hurst AW (2008) Multi-stage volcanic events: tephra hazard simulations for the Okataina Volcanic Centre, New Zealand. J Geophys Res 113 (F04012)

  • Jenkins S, Magill C, McAneney J, Blong R. Regional ash fall hazard I: a probabilistic assessment methodology. Bull Volcanol. 2012a;74(7):1699–712.

    Article  Google Scholar 

  • Jenkins S, McAneney J, Magill C, Blong R. Regional ash fall hazard II: Asia-Pacific modelling results and implications. Bull Volcanol. 2012b;74(7):1713–27.

    Article  Google Scholar 

  • Kandlbauer J, Sparks R. New estimates of the 1815 Tambora eruption volume. J Volcanol Geotherm Res. 2014;286:93–100.

    Article  Google Scholar 

  • Kandlbauer J, Carey S, Sparks RS. The 1815 Tambora ash fall: implications for transport and deposition of distal ash on land and in the deep sea. Bull Volcanol. 2013;75(4):1–11. doi:10.1007/s00445-013-0708-3.

    Article  Google Scholar 

  • Kistler R, Collins W, Saha S, White G, Woollen J, Kalnay E, et al. The NCEP-NCAR 50-year reanalysis: Monthly means CD-ROM and documentation. Bull Am Meteorol Soc. 2001;82(2):247–67.

    Article  Google Scholar 

  • Lin I-C, Tung CC. A preliminary investigation of tsunami hazard. Bull Seismol Soc Am. 1982;72(6A):2323–37.

    Google Scholar 

  • Macedonio G, Costa A, Folch A. Ash fallout scenarios at Vesuvius: Numerical simulations and implications for hazard assessment. J Volcanol Geotherm Res. 2008;178:366–77.

    Article  Google Scholar 

  • Magill CR, Hurst AW, Hunter LJ, Blong RJ. Probabilistic tephra fall simulation for the Auckland Region, New Zealand. J Volcanol Geotherm Res. 2006;153(3–4):370–86.

    Article  Google Scholar 

  • McGuire RK. Probabilistic seismic hazard analysis and design earthquakes: closing the loop. Bull Seismol Soc Am. 1995;85(5):1275–84.

    Google Scholar 

  • McGuire RK. Probabilistic seismic hazard analysis: Early history. Earthquake Eng Struct Dyn. 2008;37(3):329–38.

    Article  Google Scholar 

  • McKee C, Johnson RW, Lowenstein PL, Riley SJ, Blong RJ, De Saint OP, et al. Rabaul Caldera, Papua New Guinea: volcanic hazards, surveillance, and eruption contingency planning. J Volcanol Geotherm Res. 1985;23:195–237.

    Article  Google Scholar 

  • Miller V, Bear-Crozier A, Newey V, Horspool N, Weber R. Probabilistic Volcanic Ash Hazard Analysis (PVAHA) II: Assessment of the Asia-Pacific region using VAPAH. J Appl Volcanol. 2016; doi:10.1186/s13617-016-0044-3.

  • Musson R. The use of Monte Carlo simulations for seismic hazard assessment in the UK. Annali di Geofisica. 2000;43(1):1-9.

  • Newhall CG, Self S. The volcanic explosivity index/VEI/- An estimate of explosive magnitude for historical volcanism. J Geophys Res. 1982;87(C2):1231–8.

    Article  Google Scholar 

  • Pfeiffer T, Costa A, Macedonio G. A model for the numerical simulation of tephra fall deposits. J Volcanol Geotherm Res. 2005;140:273–94.

    Article  Google Scholar 

  • Power W, Wang X, Lane E, Gillibrand P. A probabilistic tsunami hazard study of the auckland region, part I: propagation modelling and tsunami hazard assessment at the shoreline. Pure Appl Geophys. 2013;170(9–10):1621–34.

    Article  Google Scholar 

  • Rikitake T, Aida I. Tsunami hazard probability in Japan. Bull Seismol Soc Am. 1988;78(3):1268–78.

    Google Scholar 

  • Robinson D, Fulford G, Dhu T (2005) EQRM: Geoscience Australia's Earthquake Risk Model: Technical Manual Version 3.0. Geoscience Australia Record 2006/01. Geoscience Australia. Canberra. pp. 148.

  • Robinson, D, Dhu, T, Schnieder, J. Practical probabilistic seismic risk analysis: a demonstration of capability. Seismol Res Lett. 2006;77(4):453-459.

  • Scollo S, Folch A, Costa A. A parametric and comparative study of different tephra fallout models. J Volcanol Geotherm Res. 2008a;176:199–211.

  • Scollo S, Tarantola S, Bonadonna C, Coltelli M, Saltelli A (2008b) Sensitivity analysis and uncertainty estimatio for tephra dispersal models. J Geophys Res. doi:10.1029/2006JB004864.

  • Scollo S, Prestifilippo M, Spata G, D'Agostino M, Coltelli M. Monitoring and forecasting Etna volcanic plumes. Nat Hazards Earth Syst Sci. 2009;9:1573-1585.

  • Siebert L, Simkin T, Kimberley P (2010) Volcanoes of the World. 3rd edn. Smithsonian Institution, Washington DC. University of California, Berkeley

  • Sigurdsson H, Carey S. Plinian and co-ignimbrite tephra fall from the. Bull Volcanol. 1989;51(4):243–70.

    Article  Google Scholar 

  • Simkin T, Siebert L. Volcanoes of the World: A Regional Directory, Gazetteer, and Chronology of Volcanism During the Last 10,000 Years. Tucson, Ariz: Geoscience; 1994. p. 349.

    Google Scholar 

  • Simpson A, Johnson RW, Cummins P. Volcanic threat in developing countries of the Asia-Pacific region: probabilistic hazard assessment, population risks, and information gaps. Nat Hazards. 2011;57:151–65.

    Article  Google Scholar 

  • Sørensen MB, Spada M, Babeyko A, Wiemer S, Grünthal G. Probabilistic tsunami hazard in the Mediterranean Sea. J Geophys Res Solid Earth. 2012;1978–2012:117 (B1).

    Google Scholar 

  • Sparks RSJ, Bursik M, Carey S, Gilbert JS, Graze LS, Sigurdsson H, et al. Volcanic Plumes. Chichester: Wiley and Sons; 1997.

    Google Scholar 

  • Spence RJS, Kelman I, Baxter PJ, Zuccaro G, Petrazzuoli S. Residential building and occupant vulnerability to tephra fall. Nat Hazards Earth Syst Sci. 2005;5(4):477–94.

    Article  Google Scholar 

  • Stirling M, Wilson C. Development of a volcanic hazard model for New Zealand: first approaches from the methods of probabilistic seismic hazard analysis. Bull NZ Soc Earthq Eng. 2002;35(4):266–77.

    Google Scholar 

  • Suzuki T. A theoretical model for dispersion of tephra. In: Shimozuru D, I Y (eds) Arc volcanism: Physics and tectonics. Tokyo: Terra Scientific Publishing Company; 1983.

  • TERA (1980). Seismic hazard analysis: a methodology for the Eastern United States. US Nuclear Regulatory Commisson Report No. NUREG/CR-1582.

  • Thio HK, Somerville P, Ichinose G. Probabilistic analysis of strong ground motion and tsunami hazards in Southeast Asia. J Earthquake Tsunami. 2007;1(02):119–37.

    Article  Google Scholar 

  • Thomas C, Burbidge D (2009) A Probabilistic Tsunami Hazard Assessment of the Southwest Pacific Nations. Geoscience Australia Professional Opinion 2009/02

  • Whelley P, Newhall C, Bradley K. The frequency of explosive volcanic eruptions in Southeast Asia. Bull Volcanol. 2015;77(1):1–11. doi:10.1007/s00445-014-0893-8.

    Article  Google Scholar 

  • Wilson TM, Stewart C, Sword-Daniels V, Leonard GS, Johnston DM, Cole JW, et al. Volcanic ash impacts on critical infrastructure. Phys Chem Earth Parts A/B/C. 2012;45:5–23.

    Article  Google Scholar 

  • Yokoyama I, Tilling RI, Scarpa R (1984) International Mobile Early-Warning System (s) for Volcanic Eruptions and Related Seismic Activities. Paris; UNESCO FP/2106-82-01 (2296). pp. 102.

Download references


This research was undertaken with the assistance of resources from the Australian National Computational Infrastructure (NCI), which is supported by the Australian Government. The authors gratefully acknowledge our partners at Global Volcano Model (GVM) for numerous scientific discussions throughout the development of this manuscript. The authors also wish to thank D. Robinson, J. Griffin, A. Jones and J. Schneider for their constructive reviews that greatly improved the quality of this manuscript. This paper is published with the permission of the Chief Executive Officer of Geoscience Australia.

Author information

Authors and Affiliations


Corresponding author

Correspondence to A. N. Bear-Crozier.

Additional information

Competing interest

The authors declare that they have no competing interests.

Authors’ contributions

ABC contributed to the conceptual framework for the PVAHA, the technical development of the VAPAH tool, executing the Asia-Pacific case study and was responsible for the preparation of this manuscript. VM contributed to the conceptual framework for the PVAHA, the technical development of the VAPAH tool and executing the Asia-Pacific case study. VN developed the scientific coding behind the VAPAH tool and assisted in the execution of the Asia-Pacific case study. NH contributed to the conceptual framework for the PVAHA and RW contributed to the statistical analysis for the Asia-Pacific case study. All authors read and approved the final manuscript.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Bear-Crozier, A.N., Miller, V., Newey, V. et al. Probabilistic Volcanic Ash Hazard Analysis (PVAHA) I: development of the VAPAH tool for emulating multi-scale volcanic ash fall analysis. J Appl. Volcanol. 5, 3 (2016).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: