 Research
 Open Access
 Published:
Probabilistic Volcanic Ash Hazard Analysis (PVAHA) I: development of the VAPAH tool for emulating multiscale volcanic ash fall analysis
Journal of Applied Volcanology volume 5, Article number: 3 (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.
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.
The fourstep procedure for eventbased PSHA reported by Musson (2000) is summarised below and presented in Fig. 1:

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.
An assessment for the AsiaPacific 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 AsiaPacific 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 subregions of the AsiaPacific 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.
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)
A probabilistic framework for assessing volcanic ash hazard at multiple spatial scales (PVAHA) adapted from PSHA is presented here (Fig. 2). The modified fourstep procedure is outlined below:

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.
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, calderaforming 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 AsiaPacific region). The record of completeness (ROC) was then assessed for each subregion 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 AsiaPacific. An example is provided here for the Indonesia subregion (Fig. 3).
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
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):
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 subregion 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).
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. AsiaPacific) 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.
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 magnitudefrequency relationship calculations are repeated for all subregions in the database.
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, sourcetosite 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.
The procedure for calculating ALPEs for a PVAHA is described below and includes the following:

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
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/m^{2}) 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 logictree 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 sitespecific meteorological conditions are considered separately at a later stage in the procedure.
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, 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
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) supereruption in Indonesia. This approach combines timedependant 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 bestfitting 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 TamboraF2 eruption is reported in Table 5.
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.55Ø) ash and 95 % of the sub31 μ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:
where n is the total number of measurements and Ki = Mi/Hi is the ratio of measured thickness (load) at Mi,,the ithe 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:
The bestfit results from FALL3D are reported in Table 5 and indicate that the MER was ~1.2 x 10^{8} kg/s, the eruption column height was ~33 km and the total mass deposits as fallout was ~ 1.2 x 10^{14} kg. These results are in good agreement with estimates made by Sigurdsson and Carey (1989) based on field observations. The corresponding simulated, bestfit, 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 bestfit results are further emphasised by the Aida index values; reflecting a geometric average, K = 0.98 and a geometric standard deviation, k = 1.29.
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)
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 magnitudefrequency 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.
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.
Sixtyfour years of meteorological data (January 1950  December 2014) was sourced from the NCEP/NCAR reanalysis for the Indonesia subregion, 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 (ucomponent) and zonal wind (zcomponent) 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.
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:

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.
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.
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
VAPAH generates a database file of hazard calculations collectively referred to here as the PVAHA. Like, PSHA the PVAHA results can be postprocessed 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 subregion 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/m^{2}) at each site across the region for a timeframe of interest (e.g. 1in100 year event; Fig. 10). It’s important to clarify that a 1in100 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.
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
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.
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.
References
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):511530.
Barberi F, Macedonio G, Pareshci MT, Santacroce R. Mapping tephra fallout risk: an example from Vesuvius, Italy. Nature. 1990;344:142–4.
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.
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):21522170.
Bommer JJ, Scherbaum F. The Use and Misuse of Logic Trees in Probabilistic Seismic Hazard Analysis. Earthquake Spectra. 2008;24(4):997–1009.
Bonadonna C, Houghton BF. Total grainsize distribution and volume of tephrafall deposits. Bull Volcanol. 2005;67:441–56.
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.
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.
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.
Cornell CA. Engineering seismic risk analysis. Bull Seismol Soc Am. 1968;58(5):1583–606.
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.
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.
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.
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.
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.
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.
Esteva L. Bases para la formulacion de decisiones de disen ̃o sismico. Mexico: Universidad Autonoma Nacional de Me ́xico; 1968.
Ewert JW. System for ranking relative threats of US volcanoes. Nat Hazard Rev. 2007;8(4):112–24.
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.
Folch A, Cavazzoni C, Costa A, Macedonio G. An automatic procedure to forecast tephra fallout. J Volcanol Geotherm Res. 2008a;177:767–77.
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.
Folch A, Costa A, Macedonio G. FALL3D: A computational model for transport and deposition of volcanic ash. Comput Geosci. 2009;35:1334–42.
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.
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 nuclearpower plants in the Pacific Northwest. US Geological Survey. US Numbered Series. Series Number 87297.
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.
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.
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)
Jenkins S, Magill C, McAneney J, Blong R. Regional ash fall hazard I: a probabilistic assessment methodology. Bull Volcanol. 2012a;74(7):1699–712.
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.
Kandlbauer J, Sparks R. New estimates of the 1815 Tambora eruption volume. J Volcanol Geotherm Res. 2014;286:93–100.
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.
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.
Lin IC, Tung CC. A preliminary investigation of tsunami hazard. Bull Seismol Soc Am. 1982;72(6A):2323–37.
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.
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.
McGuire RK. Probabilistic seismic hazard analysis and design earthquakes: closing the loop. Bull Seismol Soc Am. 1995;85(5):1275–84.
McGuire RK. Probabilistic seismic hazard analysis: Early history. Earthquake Eng Struct Dyn. 2008;37(3):329–38.
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.
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.
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.
Pfeiffer T, Costa A, Macedonio G. A model for the numerical simulation of tephra fall deposits. J Volcanol Geotherm Res. 2005;140:273–94.
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.
Rikitake T, Aida I. Tsunami hazard probability in Japan. Bull Seismol Soc Am. 1988;78(3):1268–78.
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):453459.
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:15731585.
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 coignimbrite tephra fall from the. Bull Volcanol. 1989;51(4):243–70.
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.
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.
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).
Sparks RSJ, Bursik M, Carey S, Gilbert JS, Graze LS, Sigurdsson H, et al. Volcanic Plumes. Chichester: Wiley and Sons; 1997.
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.
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.
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/CR1582.
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.
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/s0044501408938.
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.
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.
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.
Author information
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 AsiaPacific 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 AsiaPacific case study. VN developed the scientific coding behind the VAPAH tool and assisted in the execution of the AsiaPacific case study. NH contributed to the conceptual framework for the PVAHA and RW contributed to the statistical analysis for the AsiaPacific 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 (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.
About this article
Received
Accepted
Published
DOI
Keywords
 Probabilistic
 Volcanic ash
 Statistical emulator
 Hazard
 Computational modelling