Probabilistic Volcanic Ash Hazard Analysis (PVAHA) II: assessment of the Asia-Pacific region using VAPAH

Abstract

Volcanic ash is an increasingly common, long-range hazard, impacting on our globalised society. The Asia-Pacific region is rapidly developing as a major contributor to the global population and economy and is home to one-quarter of the world’s active volcanoes. Here we present a regional-scale volcanic ash hazard assessment for the Asia-Pacific using a newly developed framework for Probabilistic Volcanic Ash Hazard Analysis (PVAHA). This PVAHA was undertaken using the Volcanic Ash Probabilistic Assessment of Hazard (VAPAH) algorithm. The VAPAH algorithm considered a magnitude-frequency distribution of eruptions and associated volcanic ash load attenuation relationships for the Asia-Pacific, and integrated across all possible events to arrive at an annual exceedance probability for sites of interest. The Asia-Pacific region was divided into six sub-regions (e.g. Indonesia, Philippines and Southeast Asia, Melanesia/Australia, Japan/Taiwan, New Zealand/Samoa/Tonga/Fiji and Russia/China/Mongolia/Korea) characterised by 276 source volcanoes each with individual magnitude-frequency relationships. Sites for analysis within the Asia-Pacific region were limited to land-based locations at 1-km grid spacing, within 500 km of a volcanic source. The Indonesian sub-region exhibited the greatest volcanic ash hazard in the region at the 100-year timeframe, with additional sources (in Japan, the Philippines, Papua New Guinea, Kamchatka - Russia and New Zealand) along plate boundaries manifesting a high degree of hazard at the 10,000-year timeframe. Disaggregation of the volcanic ash hazard for individual sites of interest provided insight into the primary causal factors for volcanic ash hazard at capital cities in Papua New Guinea, the Philippines and Japan. This PVAHA indicated that volcanic ash hazard for Port Moresby was relatively low at all timeframes. In contrast to this, Jakarta, Manila and Tokyo are characterised by high degrees hazard at all timeframes. The greatest hazard was associated with Tokyo and the PVAHA was able to quantify that the large number of sources impacting on this location was the causal factor contributing to the hazard. This evidence-based approach provides important insights for decision makers responsible for strategic planning and can assist with prioritising regions of interest for more detailed volcanic ash hazard modelling and local scale planning.

Background

Explosive eruptions pose a serious hazard to both society and the environment (Blong 1984; Spence et al. 2005; Horwell and Baxter 2006; Durant et al. 2010; Wilson et al. 2012). Volcanic ash can be transported up to thousands of kilometres from a source, and understanding the spatial distribution of fallout from all potential events is essential for quantifying the risk to a region. The occurrence of volcanic ash fallout may not present an immediate danger to large populations; however, the long lasting nature of the hazard can cause serious disruption to aviation routes, critical infrastructure and agriculture production (Wilson et al. 2012). The potentially widespread impact of volcanic ash hazard, coupled with increased globalisation of communication, transportation and commercial networks, highlights the importance of characterising the long-term implications for strategic planning. Furthermore, regional- to global-scale analyses will enable rapid determination of at-risk zones that would benefit from focused, local hazard and risk analysis studies.

The potential impacts of volcanic ash fallout are widespread, varying and highly dependent on the scale of the eruption and the distance from source (Blong 1984; Heiken et al. 1992; Spence et al. 2005; Horwell and Baxter 2006; Costa et al. 2009; Durant et al. 2010). These impacts include but are not limited to: (1) damage to human settlements and buildings in the form of roof collapse from ash loading; (2) disruption of transportation systems due to a decrease in or loss of visibility, covering of roads/railways by ash or direct damage to vehicles; (3) partial or total destruction of agricultural crops, damage to forestry, decrease in soil permeability and increased surface run-off promoting flooding; (4) destruction of pastures and health risks for livestock (e.g. fluorosis); (5) disruption of communication systems (e.g. equipment and power lines); (6) temporary shutdown of airports due to degraded engine performance, failure of navigation equipment or loss of visibility; (7) volcanic ash leaching, which can lead to chemical and physical changes in the quality of open water supplies; (8) adverse health effects (e.g. irritation of eyes and skin and potential respiratory symptoms) associated with ash inhalation. These potential impacts stress the socio-economic implications of volcanic ash fallout and highlight the relevance of adequate hazard assessment and risk mitigation (Folch et al. 2008a; Folch et al. 2008b; Folch and Sulpizio 2010).

Previous probabilistic ash hazard assessments in the Asia-Pacific region

Previous probabilistic analyses of volcanic ash hazard in the Asia-Pacific region have largely focused on qualitative assessments of the frequency and potential consequences of future eruptions (e.g. Simpson et al. 2011). Quantitative approaches have typically considered a single source or site of interest at a local scale with any determined probabilities conditional on an eruption from that volcano occurring (Bonadonna et al. 2005; Jenkins et al. 2008; Costa et al. 2009; Folch and Sulpizio 2010; Bear-Crozier et al. 2012). Jenkins et al. (2012a; 2012b) developed a probabilistic framework for assessing ash fall hazard on a regional-scale (100’s of km) by up-scaling the ASHFALL dispersion model for regional-scale applications (Hurst 1994; Hurst and Turner 1999). Limitations associated with this approach centre on the computationally intensive nature of up-scaling ash dispersal modelling requiring large numbers of CPUs, significantly increased simulation times and subsequent limitations on the number of events that can be considered (in the order of hundreds to thousands). Dispersion models are not typically optimised for parallel processing on high performance computing infrastructure at a scale needed to resolve the hundreds of thousands or even millions of events considered in a fully probabilistic analysis. Additionally, this method does not produce hazard curves of annual exceedance probability versus volcanic ash for individual sites of interest, which precludes disaggregation of contributing factors for the hazard.

This paper presents an application example for a newly developed methodology at Geoscience Australia, which seeks to overcome these limitations through the adaptation of a well-established approach to Probabilistic Seismic Hazard Assessment (PSHA) for volcanic ash at a regional-scale. The Probabilistic Volcanic Ash Hazard Analysis (PVAHA) framework considers the magnitude-frequency distribution of eruptions and associated volcanic ash load attenuation relationships to produce and statistically emulate volcanic ash hazard for all events across a region of interest. An algorithm developed to facilitate PVAHA, the Volcanic Ash Probabilistic Assessment tool for Hazard (VAPAH), is used here to present an application example for the Asia-Pacific region.

The PVAHA framework

A framework for PVAHA was developed to investigate volcanic ash hazard in terms of the annual exceedance probability for each site across a region of interest (Bear-Crozier et al. 2016). For a full account of the conceptual frame for PVAHA and the development and validation of the VAPAH algorithm, the reader is referred to the companion paper (Bear-Crozier et al. 2016). A brief summary of the PVAHA framework as pertinent to the Asia-Pacific regional case study is provided below.

PVAHA consists of four steps:

• 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 magnitude-frequency relationships of past events.

• 3. For a catalogue of synthetic events, 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.

The VAPAH algorithm

The VAPAH algorithm developed as a mechanism to facilitate PVAHA utilises high performance computing technology in order to undertake assessments at multiple spatial scales. The VAPAH algorithm combines magnitude-frequency relationships, a catalogue of Ash Load Prediction Equations (ALPEs) 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. The reader is referred to the companion paper Bear-Crozier et al. (2016) for a detailed account of the development of the VAPAH algorithm. A brief summary of the VAPAH procedure is provided below and presented in Fig. 1.

Fig. 1

Schematic representation of the operational procedure for the VAPAH tool (Bear-Crozier et al. 2016)

Inputs for the VAPAH algorithm include:

• 1. Identification of volcanic sources for analysis.

• 2. Characterisation of magnitude-frequency relationships for each volcanic source.

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

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

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

VAPAH generates a database file of hazard calculations collectively referred to here as the PVAHA. PVAHA results can be post-processed using VAPAH to generate hazard curves for annual probability of exceedance at sites of interest (or all sites), maps for maximum expected ash hazard at timeframes of interest and histograms which disaggregate the hazard (e.g. location, magnitude, source etc.) for determining the primary causal factors at sites of interest.

Methodology

Preparation of input datasets for VAPAH

Identification of volcanic sources for the Asia-Pacific region

Potential volcanic sources for the Asia-Pacific region were identified using the Smithsonian Institution’s Global Volcanism Program (GVP) database of current and past eruptive activity (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. Sources were grouped by sub-region for a more robust statistical analysis. Six sub-regions were delineated in accordance with the Smithsonian catalogue structure (Fig. 2):

Fig. 2

Map of the Asia-Pacific region study area featuring sub-regions of interest. Potential volcanic sources are indicated by red triangles. The 500 km buffered region used for the extent of the site grid is shown in blue

• A. PNG, Vanuatu, the Solomon Islands and Australia

• B. Philippines and Southeast Asia

• C. Japan and Taiwan

• D. Indonesia

• E. New Zealand, Samoa, Tonga and Fiji

• F. Kamchatka and the Kuriles (Russia), China, Mongolia and Korea

Following the procedure outlined in Bear-Crozier et al. (2016) for source selection, database fields including volcano ID, region, sub-region, volcano type, volcano name, latitude, longitude, eruption year and Volcano Explosivity Index (VEI; an indication of magnitude; Newhall and Self 1982) were captured for each eruption at each source for the Asia-Pacific. These selections were further interrogated and sources classified as submarine, hydrothermal, fumerolic or of unknown type were discarded. Sources with no assigned VEI but designated caldera ‘C’ or Plinian ‘P’ were allocated to the larger magnitude class as arbitrary VEI 4 events. The authors acknowledge that Plinian style, caldera-forming eruptions are commonly associated with eruptions greater than VEI 4; however, VEI 4; is selected as a minimum magnitude representing a conservative estimate for the magnitude for the smaller number of these events given the absence of further information (only 13 events identified where this is the case for the Asia-Pacific region). A total of 276 sources were identified across the six sub-regions for 3997 historical events in the catalogue (Fig. 2; Table 1).

Table 1 Asia-Pacific region sources identified for the PVAHA

The annual rate of occurrence of events for each magnitude was used to characterise the magnitude-frequency relationships. The annual rate of occurrence is defined here as the total number of eruptions from a given volcano divided by the length of time for which the record is assumed to be complete (i.e. Record of Completeness - ROC). The completeness of the record for each source under consideration is highly variable owing mostly to the data sparse nature of this region. Following the method outlined in Bear-Crozier et al. (2016), the ROC for each sub-region established for small magnitude events (VEI 3 or below) and large magnitude events (VEI 4 or above) using the break in slope method (Fig. 3; Table 2; Additional file 1). In this instance, each slope break was manually identified (refer to red line on Fig. 3; Additional file 1); however, recent work has identified ways to automate this process (Mead and Magill 2014).

Fig. 3

Cumulative number of small (VEI 2–3; top) and large (VEI 4–7; bottom) events over time for three of the six sub-regions a PNG, Vanuatu, Solomon Islands, Australia, b Philippines, Southeast Asia, c Japan, Taiwan. Additional plots for sub-regions d Indonesia, e New Zealand, Samoa, Tonga, Fiji and f Kamchatka (Russia), China, Mongolia, Korea) can be found in Additional file 1

Table 2 Sub-region record of completeness (total number of eruptions included in each sub-region)

Magnitude-frequency relationships for Asia-Pacific sources

Events within the ROC for each sub-region were included in the magnitude-frequency relationship calculations (Additional files 2, 3, 4, 5, 6 and 7). Each source was classified into a volcano type category according to morphology and previous eruption style (Caldera, Large Cone, Shield, Lava Dome, Small Cone), following Jenkins et al. (2012a, b). The probability of an event occurring irrespective of magnitude and the conditional probability of an event of particular magnitude were calculated following the methodology outlined in Bear-Crozier et al. (2016) to arrive at the annual eruption probability for each source. In this study, 22 volcanoes were assigned proxy values for annual rate of occurrence (highlighted in bold; Table 1). In this instance, an average value for all volcanoes in a given category (i.e. large cone) was used as a proxy ROC (Bear-Crozier et al. 2016). In this way, sources that have more detailed historical records are used as analogues for those that are underrepresented in the database. The breakdown of eruptions for each sub-region ROC calculation is outlined in Table 2. A subset of annual eruption probabilities for sources across the region are reported in Table 3, and the reader is referred to Additional files 8, 9, 10, 11, 12 and 13 for all remaining annual eruption probabilities for the Asia-Pacific.

Table 3 Subset of sources across the Asia-Pacific and associated annual eruption probabilities

Characterisation of the volcanic ash load attenuation relationship

A catalogue of 1056 synthetic events was developed, with each event assigned an equal weighting for probability of occurrence following the methodology of Bear-Crozier et al. (2016). Events are not source specific but rather represent a suite of synthetic eruptions, which when coupled with magnitude-frequency statistics and prevailing meteorological conditions for a region of interest, can be used to represent a range of potential events at any volcanic source. The relationship for rate of volcanic ash load decay with distance from the source, as a function of magnitude, column height, duration, wind turbulence, direction and speed, is needed for each event as an input to the VAPAH algorithm. The ALPE catalogue reported in Bear-Crozier et al. (2016) is used here to characterise this attenuation relationship. 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. The ALPEs employed here were developed using the 3D time-dependent Eulerian model FALL3D (Folch et al. 2012). The reader is referred to the companion paper Bear-Crozier et al. (2016) for background on generation of the ALPE catalogue used for this study.

A spatial grid of pre-determined resolution clipped to the domain extent

The spatial extent of the domain was customised for the Asia-Pacific with consideration given to available high performance computing resources. The finer the resolution of the grid, the greater the number of sites and therefore the higher the computing demands for random access memory (RAM) and central processing units (CPUs). Definition of the full extent of the Asia-Pacific region resulted in a grid measuring 111° by 107°. The gridded area was then clipped to land-based sites within 500 km of a source (a conservative distal limit imposed by the authors for ground loading of volcanic ash using dispersal models; Fig. 2). The resolution of the grid was 0.01° × 0.01° (approximately 1 km²) and the resulting analysis considered 4,473,375 gridded nodes (sites) across an area of approximately 4,600,000 km² (Fig. 2).

Characterisation of meteorological conditions

Sixty-four years of meteorological data (January 1950 - December 2014) was sourced from the NCEP/NCAR re-analysis for Asia-Pacific region, available at grid intervals of 2.5° globally (Kistler et al. 2001). Wind direction and wind speed vector components are available through this global reanalysis for 17 pressure levels to a height of 40 km. Following the methodology of Bear-Crozier et al. (2016), monthly mean vector components for meridional wind (u-component) and zonal wind (v-component) were extracted at 60 locations across the Asia-Pacific region, from NCEP grid points closest to each volcanic source at the 250 mb pressure level (Tropopause) for a 64-year period (Fig. 4). 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. Prevailing wind direction and wind speed were assigned to each source from the closest NCEP point. A subset of monthly mean wind direction and wind speed aggregated for a 64-year period for the Asia-Pacific region are reported in Table 4, and the reader is referred to Additional file 14 for all remaining meteorological parameters used for this study.

Fig. 4

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

Table 4 Monthly mean wind direction and wind speed aggregated for a 64-year period (1950–2014) for a subset of NCEP grid points across the Asia-Pacific region

Application of the VAPAH algorithm for the Asia-Pacific region

Configuring and running VAPAH for the Asia-Pacific region

The operational procedure for the VAPAH algorithm is summarised in Fig. 1, and the reader is referred to the companion paper Bear-Crozier et al. (2016) for further detail. A configuration file was used to customise the assessment for the Asia-Pacific region. The configuration file calls a series of CSV (comma separated value) files including:

• 1. The ALPE catalogue developed in Bear-Crozier et al. (2016).

• 2. The volcano sources and assigned magnitude-frequency relationships, including prevailing wind speed and direction derived through NCEP reanalysis.

• 3. The sites of interest consisting of a pre-processed spatial grid at 0.01° resolution clipped to a 500 km radius around volcanoes of interest.

The configuration file was further customised for the Asia-Pacific region, details are presented in Table 5.

Table 5 VAPAH configuration file parameters for Asia-Pacific study

For this study, 276 sources were simulated across a 0.01° site grid with one degree resolution wind direction spacing and 16 individual wind speed directions. The VAPAH algorithm was initialised in parallel using the Raijin Fujitsu Primenergy HPC located at the National Computing Infrastructure (NCI) at the Australian National University (ANU). Memory requirements fluctuated between 32 and 248 GB depending upon the number of sites (this varied spatially). The assessment considered over 1 billion potential events and was completed on 32 cores (CPUs) over a period of approximately 8 weeks.

Results

PVAHA of the Asia-Pacific

VAPAH generated a database file of hazard calculations collectively referred to as a PVAHA for the Asia-Pacific region. Hazard calculations were aggregated for each site and maps were generated which display the maximum expected ash hazard (kg/m²) at each site across the Asia-Pacific at the 100-, 1000- and 10,000-year timeframes (e.g. 1-in-100 year event; Fig. 5). Parts of the Indonesian, Japan and Russian Kamchatka Peninsula (Russian) sub-regions exhibit the greatest volcanic ash hazard occurrence across the Asia-Pacific for all observable timeframes. The remaining three sub-regions, in particular the Philippines, PNG and the north island of New Zealand, exhibit higher degrees of volcanic ash hazard at the longer timeframes. The Indonesian sub-region is subjected to a prevailing easterly wind (blowing towards the west) and has 74 sources contributing to the hazard at this regional scale, more than any other sub-region. A number of these sources have high eruption recurrence rates across all VEIs in the range and 10–20 % of the Indonesian sub-region land mass is subjected to >1 kg/m² of volcanic ash loading, even at the shorter 100-year timeframe (Fig. 5). This is observed predominantly on the island of Java. For higher ash loads, capable of causing building damage and roof collapse (100–500 kg/m²), the hazard is typically more localised proximal to individual sources with the greatest hazard in east Java where a dense concentration of sources (Tengger Caldera, Semeru, Lamongan centres) combine to create a higher load value zone 80–120 km in diameter. Japan (Northeast Honshu Arc), the Philippines (Mayon, Kanlaon, Bulusan volcanoes), PNG (Bismarck Volcanic Arc), New Zealand (central Taupo Volcanic Zone) and the Kamchatka peninsula (dominated by Kluchevskoi volcano) also feature localised, proximal volcanic ash loads of >100 kg/m² at the lower timeframe (Fig. 5).

Fig. 5

Estimated max ash load for the Asia-Pacific region at the a 100-year timeframe, b 1000-year timeframe and c 10,000-year timeframe. The colour bar indicates the maximum ash load in kg/m²

As the timeframe increases, the level of volcanic ash hazard in the Asia-Pacific region is also found to increase. At the 1000-year timeframe, widespread regions of Indonesia, Japan, the Philippines, PNG, the Kamchatka peninsula and New Zealand show elevated estimates of maximum ash load above the threshold of 1 kg/m² (Fig. 5). At the 10,000-year timeframe, all countries along the island arcs and plate boundaries of the region show widespread and high volcanic ash hazard.

It should also be noted that there are areas of high ash hazard observed in some of the smaller Pacific island nations, such as Vanuatu. The hazard here is more difficult to observe at the regional-scale, given the small area of these island land masses and the larger regional groupings in this analysis. Whilst the total area likely to be impacted by these levels of ash may be small, it will comprise a larger percentage of the habitable land mass in these nations.

Disaggregation for sites of interest

VAPAH was used to disaggregate the ash hazard as a function of volcanic source for four capital cities of interest; Jakarta (Indonesia), Port Moresby (PNG), Manila (Philippines) and Tokyo (Japan) in order to ascertain the primary causal factors contributing to the hazard at these locations. The estimated maximum ash load hazard is further broken down by VEI for three timeframes (100, 1000, and 10,000 years). Figure 6 features the site locations for each city on enlarged hazard maps for Indonesia, PNG, the Philippines, and Japan at the 1000-year timeframe.

Fig. 6

Hazard maps of volcanic ash load at 1 km resolution for the 1000 year timeframe for a Indonesia, b Japan, c the Philippines and d Papua New Guinea

Jakarta, Indonesia

A map of maximum expected volcanic ash load on a 1-km grid for Indonesia is presented here for a 1-in-1000 year event in Fig. 6. A hazard curve for the annual probability of exceedance versus volcanic ash load (kg/m²) generated for a site in Jakarta, Indonesia and maximum expected ash load (kg/m²) at timeframes of interest derived from the data used to generate the hazard curves are presented in Fig. 7. Histograms disaggregating the percentage contribution to the hazard by source and magnitude for this site for a 1-in-100 year, 1-in-1000 year and 1-in-10,000 year event are also presented in Fig. 7. Seventeen sources in the Indonesian sub-region contribute to the hazard for Jakarta (Dempo, Besar, Suoh, Krakatau, Gagak, Salak, Gede Tangkubanparahu, Guntur, Galunggung, Ciremai, Slamet, Dieng V.C, Sundoro, Merbabu and Merapi; Fig. 8). The maximum calculated ash load for the site in Jakarta at the 100-year timeframe was found to be ~0.013 kg/m². The density of sources within the vicinity of this site resulted in a somewhat even distribution of percentage contribution from each source. Salak was the primary contributor at ~10.5 % (predominantly for events VEI 3 and VEI 4) with other primary contributions from a cluster of sources at Gagak, Gede and Tangkubanparahu each contributing 8–10 % of the ash load hazard for the site at this timeframe. At the 1000-year timeframe, the estimated ash load increased to 0.13 kg/m². Other primary contributors included Gagak (11 %) and Gede (10.5 %) with five additional sources (Krakatau, Tangkubanparahu, Guntur, Galunggung and Ciremai) contributed to >5 % of the ash load at the site for this timeframe. The expected maximum ash load at the 10,000-year timeframe increased to ~11 kg/m². Salak, Gagak and Gede are again major contributors to the hazard for this timeframe (25, 23 and 21 % respectively) and the more distal sources of Krakatau, Tangkubanparahu, Papandayan, Guntur, Galunggung and Ciremai have a greater percentage contribution to the hazard for this site at this longer timeframes (collectively ~35 % contribution; Fig. 8).

Fig. 7

a Annual exceedance probability curve for a site in Jakarta, Indonesia. Histograms disaggregating the % contribution to ash load hazard of magnitude (VEI) for each source at the site in Jakarta at the b 100, c 1000 and d 10,000 year timeframe

Fig. 8

a Annual exceedance probability curve for a site in Port Moresby, Papua New Guinea. Histograms disaggregating the % contribution to ash load hazard of magnitude (VEI) for each source at the site in Port Moresby at the b 100, c 1000 and d 10,000 year timeframe

Port Moresby, PNG

A map of maximum expected volcanic ash load on a 1-km grid for Papua New Guinea is presented for a 1-in-1000 year event in Fig. 6. A hazard curve for the annual probability of exceedance versus volcanic ash load (kg/m²) generated for a site in Port Moresby, Papua New Guinea and maximum expected ash load (kg/m²) at timeframes of interest derived from the data used to generate the hazard curves are presented in Fig. 8. Histograms disaggregating the percentage contribution to the hazard by source and magnitude for this site for a 1-in-100 year, 1-in-1000 year and 1-in-10,000 year event are also presented in Fig. 8. Six sources in the PNG sub-region contribute to the hazard for Port Moresby and include Long Island, Ritter Island, Langila, Lamington, Victory, and Waiowa (Fig. 8). The maximum calculated ash load for Port Moresby at the 100-year timeframe was found to be approaching zero (~1.1 × 10⁻⁹⁴ kg/m²). Lamington was the major relative contributor to the hazard at this site at the 100-year timeframe with events at all magnitudes represented totalling 33 % of the ash load anticipated for this site. At the 1000-year timeframe, the estimated ash load increased to ~1.0 × 10⁻⁶ kg/m². Lamington was again the primary contributor (~52 % of the hazard) dominated by VEI 3 and VEI 4 eruptions contributing to >30 % of the expected ash load. Victory and Waiowa are also major contributors for this timeframe (23 and 24 % respectively). The expected maximum ash load at the 10,000-year timeframe increased to ~0.2 kg/m². Again, the nearby sources of Lamington, Victory and Waiowa are major contributors to the hazard for this timeframe (30, 23 and 17 % respectively) however the more distal sources of Long Island, Ritter Island and Langlia also represent a hazard to Port Moresby for this longer timeframe (collectively ~30 % contribution; Fig. 8).

Manila, the Philippines

A map of maximum expected volcanic ash load on a 1-km grid for the Philippines is presented for a 1-in-1000 year event in Fig. 6. A hazard curve for the annual probability of exceedance versus volcanic ash load (kg/m²) generated for a site in Manila in the Philippines and maximum expected ash load (kg/m²) at timeframes of interest derived from the data used to generate the hazard curves are presented in Fig. 9. Histograms disaggregating the percentage contribution to the hazard by source and magnitude for this site for a 1-in-100 year, 1-in-1000 year and 1-in-10,000 year event are also presented in Fig. 9. Six sources contribute to the expected ash hazard at the site in Manila (Bulusan, Mayon, Taal, Pinatubo, Ambalatungan Group, Cagua; Fig. 9). The maximum calculated ash load for this site in Manila at the 100-year timeframe is ~0.06 kg/m². Taal was responsible for the majority of the hazard at the 100-year timeframe with events at all magnitudes represented (mostly VEI 3 and VEI 4) totalling ~27 % of the ash load anticipated for this site. Pinatubo was the second largest contributor (~23 %) again with all magnitudes represented (predominantly VEI 3 and VEI 4). Bulusan, Mayon, the Ambalatungan Group and Cagua collectively represented the remaining 50 % of the ash load contribution. At the 1000-year timeframe the expected maximum ash load increased to ~5.13 kg/m². Again, Taal was the primary ash load contributor (~57 %) for this longer timeframe and Pinatubo was again the second greatest contributor (~36 %). The contributions from both were predominantly VEI 4 and above events. Minimal contributions from other sources were observed for this timeframe (Fig. 9). The expected ash load increased to ~688 kg/m² at the 10,000-year timeframe. Taal and Pinatubo were the only contributors to the ash load hazard for the site in Manila at this timeframe again accounting for approximately 70 and 30 % of the total anticipated ash load respectively, with events at VEI 4 and above only (Fig. 9).

Fig. 9

a Annual exceedance probability curve for a site in Manila, the Philippines. Histograms disaggregating the % contribution to ash load hazard of magnitude (VEI) for each source at the site in Manila at the b 100, c 1000 and d 10,000 year timeframe

Tokyo, Japan

A map of maximum expected volcanic ash load on a 1-km grid for Japan is presented here for a 1-in-1000 year event in Fig. 6. A hazard curve for the annual probability of exceedance versus volcanic ash load (kg/m²) generated for a site in Tokyo, Japan and maximum expected ash load (kg/m²) at timeframes of interest derived from the data used to generate the hazard curves are presented in Fig. 10. Histograms disaggregating the percentage contribution to the hazard by source and magnitude for this site for a 1-in-100 year, 1-in-1000 year and 1-in-10,000 year event are also presented in Fig. 10. Over 30 sources contribute to the ash hazard for Tokyo at the timeframes of interest (Isu-Tobu, Hakone, Fuji, On-take, Haku-san, Norikura, Yake-dake, Niigata-Yake-yama, Tate-yama, Myoko, Asama, Kusatsu-Shirane, Haruna, Akagi, Hiuchi, Nikko-Shirane, Takahara, Nasu, Numazawa, Bandai, Adatara, Azuma, Zao, Kurikoma, Chokai, Akita-Komaga-take, Iwate, Akita-Yake-yama, Oshima, Nii-jima, Kozu-shima, Miyake-jima, Hachijo-jima, Aoga-shima; Fig. 10). The maximum calculated ash load for a site in Tokyo at the 100-year timeframe is ~0.19 kg/m². The density of sources within the vicinity of this site resulted in a somewhat even distribution of percentage contribution from each source with the possible exceptions being Akita-Komaga-take, Iwate and Akita-Yake-yama, which collectively represent <4 % of the hazard for the site. Hakone represented the primary contributor (~5.5 %). Five sources contributed >4 % to the total ash accumulation for this site in Tokyo across all magnitudes for this timeframe (Isu-Tobu, Hakone, Fuji, Akagi and Oshima; Fig. 10). At the 1000-year timeframe, the estimated maximum ash load for the site increased to ~735 kg/m². Hakone and Fuji were the primary contributors (25 and 16 % respectively) by events VEI 4 or greater in magnitude (Fig. 10). Secondary contributors included Akagi (10 %), Oshima (9 %) and Isu-Tobu (8 %). All other sources contributed less than 5 % to the maximum expected ash load for this site for this timeframe. The expected ash load for this site in Tokyo increased to ~6.0 × 10³ kg/m² at the 10,000-year timeframe and Hakone again was the dominant contributor at 42 % of the total anticipated ash load predominantly from events VEI 4 and above (Fig. 10). Fuji was the second largest contributor (18 %). Additional sources contributed 5–10 % of the total expected ash load for this timeframe (Fig. 10).

Fig. 10

a Annual exceedance probability curve for a site in Tokyo, Japan. Histograms disaggregating the % contribution to ash load hazard of magnitude (VEI) for each source at the site in Japan at the b 100, c 1000 and d 10,000 year timeframe

Discussion

Benefits of the proposed framework for PVAHA

The PVAHA framework is advantageous for quantifying the uncertainty in modelling the volcanic ash hazard. The synthetic catalogue of events represents the full spectrum of input parameters and together with the associated probabilities this enables the uncertainty assumptions to be incorporated into the final results. The event-based approach, which calculates hazard curves of annual probability of exceedance versus volcanic ash load at individual sites also allows for disaggregation of the hazard results (e.g. by magnitude, distance or source). This ability to easily interrogate and understand the contributing factors would be a useful tool for determining impact at a site of interest (e.g. a densely populated city or piece of critical infrastructure). It would also provide input to prioritising funding or developing mitigation strategies for sites with multiple sources contributing to the anticipated hazard where emergency planning would be inherently more complex. In addition the results can be used for identifying priority areas for more detailed dispersion modelling at the local-scale.

The PVAHA developed for the Asia-Pacific region using VAPAH highlighted the degree of volcanic ash hazard across the region as a function of timeframe. The Indonesian sub-region was identified as the greatest contributor to ash hazard for the 100-year timeframe. The remaining five sub-regions located along plate boundaries presented considerable hazard at the longer timeframes. This PVAHA, which employed all available eruption data for the Asia-Pacific region, allowed us to quantitatively assess the relative ash hazard across the region and identify key contributing sources to the hazard at sites of interest providing an evidence base for prioritising further, more detailed, local-scale analyses. Indonesia, Papua New Guinea, the Philippines and Japan were all countries which exhibited a high density of active volcanoes. The volcanic ash hazard maps developed in this study highlighted the potential exposure of communities in these countries to volcanic ash hazard. In order to understand the impact from volcanic ash hazard for a given population, their location and their corresponding infrastructure must be taken into account.

The PVAHA for the four capital city sites of interest (Jakarta, Port Moresby, Manila and Tokyo) indicate that Port Moresby has the lowest ash hazard of the four. At the 100-year and 1000-year timescales, it was demonstrated that Tokyo had the highest anticipated hazard of the four capital cities. At the 10,000-year timescale, the hazard for Tokyo was of the same order as Manila despite the lesser number of sources contributing to the ash hazard for Manila. The number of sources and distance from source greatly influenced the resulting ash hazard for each capital city and has important implications for prioritising planning for more local scale hazard analysis (e.g. conventional ash dispersal modelling). For example if workers intended to customise monitoring strategies for a maximum likelihood event (e.g. using the 10,000-year timeframe), an understanding of the primary causal factors (e.g. which sources contribute more often and at what magnitude) is critical for developing evacuation maps and monitoring strategies. This study suggests that for Manila, a hazard map could be developed in the first instance with more detailed modelling from Taal volcano, as it is likely that 70 % of the ash hazard contributed at this timescale would be from this source. In addition, if monitoring strategies needed to be prioritised, more resources could be employed for Taal volcano. In contrast, for Tokyo there is 25 % likelihood that maximum ash hazard would come from Hakone; however, there is still a significant likelihood that Fuji would be a major contributor. This information could be used to provide evidence that resourcing for monitoring and planning should be distributed amongst a number of contributing sources for this site and that prioritisation of detailed local-scale dispersal modelling might target Hakone and Fuji when considering the hazard for Tokyo for this timeframe.

Limitations of the proposed framework for PVAHA

Limitations associated with the proposed PVAHA framework for application in the Asia-Pacific include the inherent uncertainties of dispersion models, meteorological data resolution, determination of the record of completeness for data sparse regions like the Asia-Pacific and the subsequent use of proxy magnitude-frequency relationships in the absence of observed behaviour for a source. Mitigating the impact of these limitations on the PVAHA output was addressed in a number of ways. There are two key limitations associated with the underlying dispersion models, which have been addressed by the PVAHA framework. Firstly, the ability of dispersal models to simulate ash fall on the ground at ever increasing distances (100’s km) is not well validated and secondly ash load values tend to be overestimated in the very proximal (2 km) range. To account for the first limitation, a 500 km from source distal limit was imposed for estimating ground loading of volcanic ash using conventional dispersal models (i.e. the accuracy of ground load estimates at sites >500 km from a source cannot be determined confidently). As a consequence, credible volcanic ash accumulation beyond the 500 km limit may not have been captured, and this had important implications for the Asia-Pacific assessment. For example the ash hazard results from Western New Guinea (Irian Jaya) and Kalimantan show no ash hazard at any timeframe purely as a function of location >500 km from any source. However, historically instances of ash accumulation at ground level, more than 1000 km from source are rare with the exception of very large very infrequent events (e.g. VEI 8 event at Toba, Indonesia (74,000 ka; Costa et al. 2014) and the VEI 7 event at Tambora (1815; Sigurdsson and Carey 1989). At longer time scales the footprint for a given threshold of ash loading might be expected to increase beyond the imposed 500 km limit (e.g. a Tambora 1815-style event; Sigurdsson and Carey 1989), and it is acknowledged that the distribution of ash at these longer timeframes (e.g. 10,000 years) in this PVAHA might be underestimated. To account for the second limitation, results in the first 2 km around a source are not utilised, to mitigate the effects of overestimating the proximal hazard. At the 2–10 km distance the likelihood of overestimation of the hazard is significantly reduced and would not affect a regional-scale analysis, therefore these values are included but any impact/risk assessment would take this into consideration.

Meteorological conditions are considered one of the major controls on short to medium range ash fall deposits, particularly for small to moderate magnitude eruptions (Mortazavi 2009) Whilst the effect of wind direction has been recognized for some time, the complex interactions of eruption plumes and meteorological conditions are still an active area of research (e.g. Woodhouse et al. 2013). To adequately represent the prevailing wind conditions at a source for the purpose of fully probabilistic analysis it is necessary to interrogate the full range of meteorological data over a period of years and for temporal ranges that sufficiently characterize the seasonal variations, in order to incorporate the uncertainty in any particular set of conditions at the source. For a regional-scale study it would require both time and computational resources to fully capture ‘all wind’ scenarios, rather than just subsets of data. In this study, prevailing wind speed and direction were derived from global NCEP monthly mean estimates at 60 locations across the Asia-Pacific region at the 250 mb pressure level (Tropopause). Meteorological data for the 64-year period (1950–2014) was aggregated for the time period and each source was assigned a prevailing wind speed and direction together with a standard deviation intended to capture seasonal variations in these variables. Meteorological conditions remain fixed horizontally and vertically at each source. According to an analysis of the variation in wind speed and direction horizontally using NCEP data for Indonesia, Jenkins et al. (2012a) showed that mean wind direction and speed downwind does not vary by more than one standard deviation over the maximum extent considered and therefore accounting for horizontal variations is not essential. However, Jenkins et al. (2012a) did find that there are greater vertical variations in wind conditions at one location. It is a limitation of this technique that only one pressure level is considered for each source (250 mb - Tropopause). This limitation is considered acceptable given the coarse resolution of this regional scale analysis and the intended purpose of the PVAHA framework as a mechanism for providing an evidence base for prioritising further analysis using tested and accepted methods of local scale hazard modelling (e.g. conventional dispersion modelling using fully 3D wind fields).

The magnitude-frequency calculations are dependent upon a best estimate for eruption database record of completeness for the Asia-Pacific region. Volcanoes with long repose periods might be under-represented and uncertainty in estimating the completeness of eruption data will affect the probabilities. The PVAHA undertaken here for the Asia-Pacific was highly sensitive to the accuracy of the statistical characterisation of each source and susceptible to gaps in the eruption database. The preservation of historical fall deposits in the rock record is inherently biased by the size of the eruption where larger eruptions are more likely to be preserved in the record (Simpson et al. 2011). The record completeness for a sub-region was also influenced by country specific factors such as the population and record keeping procedures. These limitations could be managed through the refinement of magnitude-frequency relationships for sources in the Asia-Pacific as new data and observations become available. This might include sustained programs of field studies and laboratory analysis aimed at better characterisation of individual source histories and behaviour. The use of proxy volcanoes when deriving the frequency-magnitude values are another limitation. In the Asia-Pacific regional study, 22 sources were assigned proxy probabilities for some VEI relations indicating a need for further targeted data collection (all other sources were assigned probabilities based on known historical behaviour). Future studies could undertake a more detailed assessment for each source however this is currently limited by the lack of geological and historical records.

Conclusions

Explosive eruptions pose a serious hazard to both society and the environment and the volcanic ash hazard associated with these events can cross international borders. Even small amounts of ash can cause disruption to civil aviation routes, critical infrastructure and other exposed elements. Countries without volcanoes still need to be aware of the hazard that might exist from neighbouring regions or those with interconnected networks and economies. Given the extensive nature of the impacts from ash hazard and the increasing interconnectivity and globalisation of communication, transport and economies, it is important to understand the hazard on the breadth of scales used for long-term strategic planning. Quantifying volcanic ash hazard at the regional scale typically requires simplification of computationally intensive approaches to ash dispersion modelling using coarser input datasets. Here we present an application example of the PVAHA framework using the VAPAH algorithm for the Asia-Pacific region. The VAPAH algorithm combined magnitude-frequency relationships for Asia-Pacific sources, an ALPE catalogue and global scale meteorological conditions for the region and integrated across all possible events to arrive at a preliminary annual exceedance probability for each site. This PVAHA, which employed all available eruption data for the Asia-Pacific region, allowed us to quantitatively assess the relative ash hazard across the region and identify key contributing sources to the hazard at sites of interest providing an evidence base for prioritising further, more detailed, local scale analyses. The results were disaggregated to highlight contributing factors for four capital cities in order to demonstrate how the PVAHA could inform evidence-based decision making and emergency planning.