 Research
 Open Access
Probabilistic Volcanic Ash Hazard Analysis (PVAHA) I: development of the VAPAH tool for emulating multiscale volcanic ash fall analysis
 A. N. BearCrozier^{1}Email author,
 V. Miller^{1},
 V. Newey^{1},
 N. Horspool^{1, 2} and
 R. Weber^{1}
https://doi.org/10.1186/s1361701600434
© BearCrozier et al. 2016
 Received: 14 July 2014
 Accepted: 17 January 2016
 Published: 27 January 2016
Abstract
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 multiscale 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 magnitudefrequency 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 multiscale 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.
Keywords
 Probabilistic
 Volcanic ash
 Statistical emulator
 Hazard
 Computational modelling
Background
Numerous approaches have been adopted in the past to assess volcanic ash hazard at the localscale (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; BearCrozier 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 upscaling 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 AsiaPacific 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 (BearCrozier et al. 2012); Okataina, New Zealand (Jenkins et al. 2008); SommaVesuvio, 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; BearCrozier et al. 2012)
Regionalscale 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 upscales implementation of the ash dispersal model ASHFALL for regionalscale assessments (Hurst 1994; Hurst and Turner 1999). This approach presented a method for assessing regionalscale 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 regionalscale applications that would require significant highperformance 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 regionalscale. 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 regionalscale named Probabilistic Seismic Hazard Analysis (PSHA; Cornell 1968; McGuire 1995, 2008). PSHA consists of a fourstep 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 highperformance computing have led to the development of eventbased PSHA. Eventbased Probabilistic Seismic Hazard Analysis allows calculation of groundmotion fields from stochastic event sets.
 1.
Seismicity data for the region of interest must be spatially disaggregated into discrete seismic sources.
 2.
For each seismic source the seismicity is characterised with respect to time (i.e. the annual rate of occurrence of different magnitudes).
 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.
This database of groundmotion fields, representative of the possible shaking scenarios that the investigated area can experience over a userspecified 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 fourstep procedure of PSHA for volcanic ash hazard analysis at a regionalscale. The framework named here, Probabilistic Volcanic Ash Hazard Analysis (PVAHA) considers the magnitudefrequency 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 regionalscale probabilistic volcanic ash hazard assessment.
 1.
A description of the proposed framework for PVAHA.
 2.
The procedure used for identification of source volcanoes, development of eruption statistics and calculation of magnitude frequency relationships for each source.
 3.
Derivation and validation of ash load prediction equations derived from volcanic ash dispersal modelling used to inform the PVAHA.
 4.
Development of the VAPAH algorithm.
Methodology
A framework for Probabilistic Volcanic Ash Hazard Analysis (PVAHA)
 1.
Volcanic sources with respect to any given site of interest must be identified.
 2.
For each volcanic source the annual eruption probability must be calculated based on magnitudefrequency relationships of past events.
 3.
For a set of stochastic events (synthetic catalogue) volcanic ash load attenuation relationships must be calculated (derived from conventional ash dispersal modelling).
 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 datasparse regions like the AsiaPacific include incomplete or nonexistent 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 AsiaPacific examples to provide context where needed (Miller et al., 2016). Database fields including volcano ID, region, subregion, 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 subregions 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.
Magnitudefrequency relationships
Having established the record of completeness for each source in the database a procedure analogous to developing earthquake magnitudefrequency 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
Annual eruption probability ‘λ (VEI 2–7)’ for each volcanic source in the Indonesian subregion using a conversion factor of 1.75
Volcano  Category  N (VEI 2–3)^{a}  λ (VEI 2–3)^{a}  N (VEI 4–7)^{b}  λ (VEI 4–7)^{b}  N (VEI 2–3) normalised^{c}  N (VEI 2–7)^{c}  λ (VEI 2–7)^{c} 

Agung  Large cone  3  0.0125  1  0.00238  5.25  6.25  0.0148 
Awu  Large cone  13  0.0541  2  0.00476  22.75  24.75  0.0589 
Batur  Caldera  26  0.1083      45.5  45.5  0.1083 
Besar  Large cone              0.0420^{a} 
Cereme  Large cone  4  0.0166      7  7  0.0166 
Dieng V.C  Large cone  22  0.0916      38.5  38.5  0.0916 
Dukono  Large cone  3  0.0125      5.25  5.25  0.0125 
Galunggung  Large cone  3  0.0125  2  0.00476  5.25  7.25  0.0172 
Gamalama  Large cone  57  0.2375      99.75  99.75  0.2375 
Gamkonora  Large cone  11  0.0458  1  0.00238  19.25  20.25  0.0482 
Gede  Large cone  23  0.0958      40.25  40.25  0.0958 
Guntur  Large cone  22  0.0916      38.5  38.5  0.0916 
Ijen  Large cone  10  0.0416      17.5  17.5  0.0416 
Kelut  Large cone  13  0.0541  6  0.01428  22.75  28.75  0.0684 
Krakatau  Caldera  41  0.1708  1  0.00238  71.75  72.75  0.1732 
LokonEmpung  Large cone  25  0.1041      43.75  43.75  0.1041 
Mahawu  Large cone  6  0.025      10.5  10.5  0.025 
Makian  Large cone  6  0.025  3  0.00714  10.5  13.5  0.0321 
Merapi  Large cone  54  0.225  2  0.00476  94.5  96.5  0.2297 
Raung  Large cone  59  0.2458  4  0.00952  103.25  107.25  0.2553 
Rinjani  Large cone  15  0.0625      26.25  26.25  0.0625 
Salak  Large cone  5  0.0208      8.75  8.75  0.0208 
Sangeang Api  Large cone  16  0.0666      28  28  0.0666 
Semeru  Large cone  59  0.2458      103.25  103.25  0.245833 
Serua  Large cone  7  0.0291  1  0.00238  12.25  13.25  0.0315 
Sinabung  Large cone  1  0.0041      1.75  1.75  0.0041 
Suoh  Caldera      1  0.00238  0  1  0.0023 
Tambora  Large cone  6  0.025  1  0.00238  10.5  11.5  0.0273 
Teon  Large cone  1  0.0041  1  0.00238  1.75  2.75  0.0065 
Tongkoko  Large cone  4  0.0166  1  0.00238  7  8  0.0190 
Wurlali  Large cone  1  0.0041      1.75  1.75  0.0041 
Conditional probability of an event of a particular magnitude
Conditional probability of an event VEI 3 or less occurring for each of the five volcano types in the database (e.g. AsiaPacific)
Category  Events  VEI 0  VEI 1  VEI 2^{a}  VEI 3 

Caldera  387  0.007  0.037  0.075  0.008 
Large cone  2422  0.023  0.171  0.522  0.079 
Lava dome  35  0.000  0.010  0.001  0.001 
Shield  168  0.001  0.008  0.042  0.005 
Small cone  33  0.000  0.001  0.007  0.004 
Conditional probability of an event VEI 4 or greater occurring for each of the five volcano types in the database (e.g. AsiaPacific)
Category  Events  VEI 4  VEI 5  VEI 6  VEI 7 

Caldera  49  0.118  0.024  0.024  0.003 
Large cone  209  0.500  0.177  0.042  0.007 
Lava dome  8  0.017  0.010  0.000  0.000 
Shield  18  0.049  0.007  0.007  0.000 
Small cone  4  0.014  0.000  0.000  0.000 
Annual probability of an event
Annual probability of an event of a particular magnitude occurring for each volcanic source in the Indonesian subregion
Volcano  P(VEI 0)  P(VEI 1)  P(VEI 2)  P(VEI 3)  P(VEI 4)  P(VEI 5)  P(VEI 6)  P(VEI 7) 

Agung  3.4E04  2.5E03  7.8E03  1.2E03  7.4E03  2.6E03  6.2E04  1.0E04 
Awu  1.4E03  1.0E02  3.1E02  4.7E03  2.9E02  1.0E02  2.5E03  4.1E04 
Batur  7.8E04  4.1E03  8.1E03  8.2E04  1.3E02  2.6E03  2.6E03  3.8E04 
Besar^{a}  9.7E04  7.2E03  2.2E02  3.3E03  2.1E02  7.5E03  1.8E03  2.9E04 
Cereme  3.8E04  2.8E03  8.7E03  1.3E03  8.3E03  3.0E03  6.9E04  1.2E04 
Dieng V.C  2.1E03  1.6E02  4.8E02  7.3E03  4.6E02  1.6E02  3.8E03  6.4E04 
Galunggung  4.0E04  2.9E03  9.0E03  1.4E03  8.6E03  3.1E03  7.2E04  1.2E04 
Gamalama  5.5E03  4.1E02  1.2E01  1.9E02  1.2E01  4.2E02  9.9E03  1.6E03 
Gamkonora  1.1E03  8.2E03  2.5E02  3.8E03  2.4E02  8.5E03  2.0E03  3.3E04 
Gede  2.2E03  1.6E02  5.0E02  7.6E03  4.8E02  1.7E02  4.0E03  6.7E04 
Guntur  2.1E03  1.6E02  4.8E02  7.3E03  4.6E02  1.6E02  3.8E03  6.4E04 
Ijen  9.6E04  7.1E03  2.2E02  3.3E03  2.1E02  7.4E03  1.7E03  2.9E04 
Kelut  1.6E03  1.2E02  3.6E02  5.4E03  3.4E02  1.2E02  2.9E03  4.8E04 
Krakatau  1.3E03  6.5E03  1.3E02  1.3E03  2.0E02  4.2E03  4.2E03  6.0E04 
LokonEmpung  2.4E03  1.8E02  5.4E02  8.3E03  5.2E02  1.8E02  4.3E03  7.2E04 
Mahawu  5.7E04  4.3E03  1.3E02  2.0E03  1.3E02  4.4E03  1.0E03  1.7E04 
Makian  7.4E04  5.5E03  1.7E02  2.6E03  1.6E02  5.7E03  1.3E03  2.2E04 
Marapi  5.7E03  4.2E02  1.3E01  2.0E02  1.2E01  4.4E02  1.0E02  1.7E03 
Merapi  5.3E03  3.9E02  1.2E01  1.8E02  1.1E01  4.1E02  9.6E03  1.6E03 
Raung  5.9E03  4.4E02  1.3E01  2.0E02  1.3E01  4.5E02  1.1E02  1.8E03 
Rinjani  1.4E03  1.1E02  3.3E02  5.0E03  3.1E02  1.1E02  2.6E03  4.3E04 
Salak  4.8E04  3.6E03  1.1E02  1.7E03  1.0E02  3.7E03  8.7E04  1.4E04 
Sangeang Api  1.5E03  1.1E02  3.5E02  5.3E03  3.3E02  1.2E02  2.8E03  4.6E04 
Semeru  5.7E03  4.2E02  1.3E01  2.0E02  1.2E01  4.4E02  1.0E02  1.7E03 
Serua  7.3E04  5.4E03  1.6E02  2.5E03  1.6E02  5.6E03  1.3E03  2.2E04 
Sinabung  9.6E05  7.1E04  2.2E03  3.3E04  2.1E03  7.4E04  1.7E04  2.9E05 
Suoh  1.7E05  8.9E05  1.8E04  1.8E05  2.8E04  5.8E05  5.8E05  8.3E06 
Tambora  6.3E04  4.7E03  1.4E02  2.2E03  1.4E02  4.8E03  1.1E03  1.9E04 
Teon  1.5E04  1.1E03  3.4E03  5.2E04  3.3E03  1.2E03  2.7E04  4.5E05 
Tongkoko  4.4E04  3.3E03  9.9E03  1.5E03  9.5E03  3.4E03  7.9E04  1.3E04 
Wurlali  9.6E05  7.1E04  2.2E03  3.3E04  2.1E03  7.4E04  1.7E04  2.9E05 
Emulating volcanic ash load attenuations relationships
 1.
Development of a synthetic catalogue of events
 2.
Volcanic ash dispersal modelling
 3.
Derivation of an ALPE
 4.
Validation of an ALPE
Development of synthetic catalogue of events
A simplified, schematic representation of the logic tree data structure used is presented in Fig. 5.

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

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

eruption style (Strombolian, Vulcanian, SubPlinian 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 magnitudefrequency 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 timedependant Eulerian model that solves the advection–diffusionsedimentation (ADS) equation on a structured terrainfollowing mesh. FALL3D outputs timedependant 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 1D 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 bestfit 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/cm^{3}) 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 magnitudefrequency 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.
Validation
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
The explored range and bestfit input parameters modelled dispersion of the TamboraF2 fall deposit using FALL3D
Modelled dispersion parameters  Explored range  TamboraF2^{a} 

Tephra mass (10^{12} kg)  0.11.5  1.2 
Tephra volume (km^{3})  0.11.5  1.0 
Tephra volume DRE (km^{3})^{b}  0.11.5  0.4 
Duration (hours)  02.8  2.8^{c} 
MER (10^{8} kg/s)  110  1.3 
TGSDmaxima (in Ø units)  08  8 
Column height (km)  3043  30 
Suzuki coefficients A (−)^{d}  14  4 
Density of aggregates (kg/m3)^{e}  100600  300 
Diameter of aggregates (in Ø units)^{e}  23  2.5 
Average deposit density (kg/m3)^{f}  Assumed  1100 
Aida Indices K/k (−)^{g}  Calculated  0.98/1.29 
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.
Validation of an ALPE for the Tambora F2 Plinian fall deposit
Following the procedure outlined above the volcanic ash attenuation relationship for the bestfit 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/m^{2}) the preferred unit of measurement for PVAHA. Validation of the ALPE was a twostep 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 bestfit (FALL3D) load (kg/m^{2}) 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 bestfit 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).
Results
Volcanic Ash Probabilistic Assessment tool for Hazard (VAPAH)
 1.
Identification of volcanic sources for analysis.
 2.
Characterisation of magnitudefrequency relationships for each volcanic source.
 3.
Characterisation of the volcanic ash load attenuation relationship (ALPE catalogue).
 4.
A spatial grid of predetermined resolution clipped to the domain extent (default = autogenerated).
 5.
Characterisation of meteorological conditions  prevailing wind direction (degrees) and wind speed (m/s)
Identification of volcanic sources, characterisation of magnitudefrequency 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 subregion 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 reanalysis (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 subregion using one potential source of meteorological data.
Monthly mean wind direction and wind speed aggregated for a 64year period (1950–2014) for 16 NCEP grid points across the Indonesian subregion
NCEP point  Longitude  Latitude  Wind direction (deg)  Wind speed (m/s) 

0  97.5  5  254  9.78 
1  97.5  2.5  262  9.90 
2  100  0  269  9.96 
3  102.5  −2.5  266  9.59 
4  105  5  262  9.28 
5  107.5  −7.5  257  7.72 
6  110  −7.5  259  7.57 
7  112.5  −7.5  261  7.52 
8  115  −7.5  262  7.49 
9  117.5  −7.5  263  7.39 
10  120  −7.5  265  7.18 
11  122.5  −7.5  269  6.95 
12  130  −7.5  267  6.29 
13  130  −5  263  8.63 
14  127.5  0  267  11.65 
15  125  2.5  253  11.41 
VAPAH algorithm procedure
 1)
the ALPE catalogue,
 2)
the volcano sources (included prevailing wind speed and direction)
 3)
the sites of interest (preprocessed 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 autogeneration 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 preconfigured (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.
 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.
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/m^{2}) at that site for the first event.
 3.
The calculated ash load for the first event and its associated probability derived from the magnitudefrequency relationships for the first source in the catalogue are written to a results file.
 4.
The algorithm repeats (2) and (3) for each ALPE for the first source.
 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.
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/m^{2}) 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.
The algorithm then loops to the next site and repeats steps (1) – (6) until all sites have been evaluated.
VAPAH results
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 pregenerated 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).
Discussion
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 timeconsuming 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 broadbrush, first approximation of the hazard (like PSHA) which can be updated and rerun 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 multiscale 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 decisionmaking (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 evidencebased 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 magnitudefrequency 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
 1.
Determination of the record of completeness for a subregion 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 underrepresented).
 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.
All events are assumed to follow a memoryless Poisson process meaning the probability of an event occurring today is not contingent on whether or not an event occurred yesterday.
 4.
A 2 km radii is applied to each source area and volcanic ash load estimates within this zone are not utilised due to overestimations 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.
This procedure and the VAPAH tool are intended to be one in a range of tools and techniques used to provide a consistent fitforpurpose 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 highperformance 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.
Conclusions
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 groundshaking intensity for timeframes of interest. The PVAHA methodology presented here modifies the fourstep procedure of PSHA for volcanic ash hazard assessment at multiple spatial scales. This technique considers a magnitudefrequency 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 upscale 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/m^{2}). 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.
Declarations
Acknowledgments
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.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), 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.
Authors’ Affiliations
References
 Aida I. Reliability of a tsunami source model derived from fault parameters. J Phys Earth. 1978;26(1):57–73.Google Scholar
 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):511530.Google Scholar
 Barberi F, Macedonio G, Pareshci MT, Santacroce R. Mapping tephra fallout risk: an example from Vesuvius, Italy. Nature. 1990;344:142–4.View ArticleGoogle Scholar
 BearCrozier A, Kartadinata N, Heriwaseso A, Nielsen O. Development of pythonFALL3D: a modified procedure for modelling volcanic ash dispersal in the AsiaPacific region. Nat Hazards. 2012;64(1):821–38.View ArticleGoogle Scholar
 Blong R (1981) Some effects of tephra falls on buildings. In: Tephra Studies. Springer, pp 405–420Google Scholar
 Blong R. Building damage in Rabaul, Papua New Guinea, 1994. Bull Volcanol. 2003;65(1):43–54.Google Scholar
 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):21522170.Google Scholar
 Bommer JJ, Scherbaum F. The Use and Misuse of Logic Trees in Probabilistic Seismic Hazard Analysis. Earthquake Spectra. 2008;24(4):997–1009.View ArticleGoogle Scholar
 Bonadonna C, Houghton BF. Total grainsize distribution and volume of tephrafall deposits. Bull Volcanol. 2005;67:441–56.View ArticleGoogle 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.View ArticleGoogle 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–537Google Scholar
 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–516Google Scholar
 Bursik M. Effect of wind on the rise height of volcanic plumes. Geophys Res Lett. 2001;18:3621–4.View ArticleGoogle 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–636Google Scholar
 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.View ArticleGoogle Scholar
 Cornell CA. Engineering seismic risk analysis. Bull Seismol Soc Am. 1968;58(5):1583–606.Google Scholar
 Cornell W, Carey S, Sigurdsson H. Computersimulation of transport and deposition of the Campanian Y5 Ash. J Volcanol Geotherm Res. 1983;17(1–4):89–109.View ArticleGoogle 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.View ArticleGoogle Scholar
 Costa A, Macedonio G, Folch A. A threedimensional Eulerian model for transport and deposition of volcanic ashes. Earth Planet Sci Lett. 2006;241:634–47.View ArticleGoogle 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.View ArticleGoogle 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) http://dx.doi.org/10.1029/2012GL051605.
 Costa A, Smith VC, Macedonio G, Matthews NE. The magnitude and impact of the Youngest Toba Tuff supereruption. Front Earth Sci. 2014;2:16.View ArticleGoogle 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.View ArticleGoogle 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.View ArticleGoogle Scholar
 Folch A, Sulpizio R. Evaluating longrange volcanic ash hazard using supercomputing facilities: application to SommaVesuvius (Itay), and consequences for civil aviation over the Central Mediterranean Area. Bull Volcanol. 2010;72:1039–59.View ArticleGoogle Scholar
 Folch A, Cavazzoni C, Costa A, Macedonio G. An automatic procedure to forecast tephra fallout. J Volcanol Geotherm Res. 2008a;177:767–77.View ArticleGoogle 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.View ArticleGoogle Scholar
 Folch A, Costa A, Macedonio G. FALL3D: A computational model for transport and deposition of volcanic ash. Comput Geosci. 2009;35:1334–42.View ArticleGoogle 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–19Google Scholar
 Geist EL, Parsons T. Probabilistic Analysis of Tsunami Hazards*. Nat Hazards. 2006;37(3):277–314.View ArticleGoogle 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.Google Scholar
 Hoblitt RP, Miller CD, Scott WE (1987) Volcanic hazards with regard to siting nuclearpower plants in the Pacific Northwest. US Geological Survey. US Numbered Series. Series Number 87297.Google Scholar
 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.View ArticleGoogle 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.Google Scholar
 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.View ArticleGoogle Scholar
 Jenkins S, Magill C, McAneney J, Hurst AW (2008) Multistage volcanic events: tephra hazard simulations for the Okataina Volcanic Centre, New Zealand. J Geophys Res 113 (F04012)Google Scholar
 Jenkins S, Magill C, McAneney J, Blong R. Regional ash fall hazard I: a probabilistic assessment methodology. Bull Volcanol. 2012a;74(7):1699–712.View ArticleGoogle Scholar
 Jenkins S, McAneney J, Magill C, Blong R. Regional ash fall hazard II: AsiaPacific modelling results and implications. Bull Volcanol. 2012b;74(7):1713–27.View ArticleGoogle Scholar
 Kandlbauer J, Sparks R. New estimates of the 1815 Tambora eruption volume. J Volcanol Geotherm Res. 2014;286:93–100.View ArticleGoogle 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/s0044501307083.View ArticleGoogle Scholar
 Kistler R, Collins W, Saha S, White G, Woollen J, Kalnay E, et al. The NCEPNCAR 50year reanalysis: Monthly means CDROM and documentation. Bull Am Meteorol Soc. 2001;82(2):247–67.View ArticleGoogle Scholar
 Lin IC, 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.View ArticleGoogle 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.View ArticleGoogle 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.View ArticleGoogle 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.View ArticleGoogle Scholar
 Miller V, BearCrozier A, Newey V, Horspool N, Weber R. Probabilistic Volcanic Ash Hazard Analysis (PVAHA) II: Assessment of the AsiaPacific region using VAPAH. J Appl Volcanol. 2016; doi:10.1186/s1361701600443.
 Musson R. The use of Monte Carlo simulations for seismic hazard assessment in the UK. Annali di Geofisica. 2000;43(1):19.Google Scholar
 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.View ArticleGoogle 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.View ArticleGoogle 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.View ArticleGoogle 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.Google Scholar
 Robinson, D, Dhu, T, Schnieder, J. Practical probabilistic seismic risk analysis: a demonstration of capability. Seismol Res Lett. 2006;77(4):453459.Google Scholar
 Scollo S, Folch A, Costa A. A parametric and comparative study of different tephra fallout models. J Volcanol Geotherm Res. 2008a;176:199–211.Google Scholar
 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:15731585.Google Scholar
 Siebert L, Simkin T, Kimberley P (2010) Volcanoes of the World. 3rd edn. Smithsonian Institution, Washington DC. University of California, BerkeleyGoogle Scholar
 Sigurdsson H, Carey S. Plinian and coignimbrite tephra fall from the. Bull Volcanol. 1989;51(4):243–70.View ArticleGoogle 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 AsiaPacific region: probabilistic hazard assessment, population risks, and information gaps. Nat Hazards. 2011;57:151–65.View ArticleGoogle 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.View ArticleGoogle 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.Google Scholar
 TERA (1980). Seismic hazard analysis: a methodology for the Eastern United States. US Nuclear Regulatory Commisson Report No. NUREG/CR1582.Google Scholar
 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.View ArticleGoogle Scholar
 Thomas C, Burbidge D (2009) A Probabilistic Tsunami Hazard Assessment of the Southwest Pacific Nations. Geoscience Australia Professional Opinion 2009/02 Google Scholar
 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/s0044501408938.View ArticleGoogle Scholar
 Wilson TM, Stewart C, SwordDaniels 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.View ArticleGoogle Scholar
 Yokoyama I, Tilling RI, Scarpa R (1984) International Mobile EarlyWarning System (s) for Volcanic Eruptions and Related Seismic Activities. Paris; UNESCO FP/21068201 (2296). pp. 102.Google Scholar