Assessing minimum pyroclastic density current mass to impact critical infrastructures: example from Aso caldera (Japan)

We describe a method for calculating the probability that a distal geographic location is impacted by a pyroclastic density current (PDC) of a given size, considering the key related uncertainties. Specifically, we evaluate the minimum volume and mass of a PDC generated at the Aso caldera (Japan) that might affect each of five distal infrastructure (marker) sites, with model input parameter uncertainties derived from expert judgment. The 5 marker sites are all located 115–145 km from the caldera; as these lie in well-separated directions, we can test the effects of the different topographic shielding effects in each case. To inform our probabilistic analysis, we apply alternative kinetic energy assessment approaches, i.e., rock avalanche and density current dynamics. In the latter formulation, the minimum mass needed to reach the markers ranges between median values of $\sim \mathrm{153}\times {\mathrm{10}}^{\mathrm{12}}$ and $\sim \mathrm{465}\times {\mathrm{10}}^{\mathrm{12}}$ kg (M 7.2–7.7), depending on the site. Rock avalanche dynamics modeling indicates that a ∼3-times greater mass would be required to reach the marker sites with 50 % probability, while the hypothetical scenario of a relatively dilute distal ash cloud would require ∼3-times less mass. We compare our results with the largest recorded Aso eruption, showing that a catastrophic eruption, similar to Aso-4, ≈ M8, would present a significant conditional probability of PDCs reaching the marker sites, in the density current formulation and contingent on uncertainty in the erupted mass and on marker site direction.

1 Introduction

Catastrophic caldera-forming eruptions can generate gigantic clouds of ash and aerosols, extensive fallout deposits, and extremely mobile pyroclastic density currents (PDCs; Baines and Sparks, 2005; Costa et al., 2014; Black et al., 2015; Oppenheimer and Donovan, 2015; Self, 2006, 2015). For instance, the 18.8 Ma Peach Spring Tuff (AZ, USA) pyroclastic flows traveled >170 km (Roche et al., 2016) and pyroclastic flows from other large eruptions have propagated >100 km from their source vents (Wilson, 1991; Wilson et al., 1995; Streck and Grunder, 1995). However, such catastrophic caldera-forming eruptions have never been observed directly. Understanding the dynamics of these large-scale PDCs is thus limited, being based on interpretations of their deposits, analogue experiments, numerical modeling, and on extrapolations derived from the observation of smaller scale events (e.g., Fisher et al., 1993; Cas et al., 2011; Shimizu et al., 2019). Consequently, there are large uncertainties that require the use of techniques to consider them in probabilistic terms (e.g., Neri et al., 2015; Rutarindwa et al., 2019; Tadini et al., 2022a; Bevilacqua et al., 2021, 2022).

Figure 1 (a) Map of Kyushu Island (South Japan), where the Aso caldera is located (inside red rectangle). The positions of the five marker sites are indicated by dotted gray circles. A sketch showing the location of Kyushu Island within Japan is also displayed in the insert. (b) Zoom on the Aso caldera. Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) 1 arcsec 30 m (NASA, 2013).

Several volcanic complexes in the world have produced catastrophic caldera-forming eruptions in the past and may do so again in the future (e.g., Lowenstern et al., 2006; Larsen et al., 2007; Aoki, 2008; Williams, 2012; Newhall et al., 2018; Suñe-Puchol et al., 2019). The expected recurrence rate of M7-8 eruptions at local and global scales has been evaluated from different datasets and through various statistical approaches (Decker, 1990; Simkin, 1993; Mason et al., 2004; Deligne et al., 2010; Brown et al., 2014; Kiyosugi et al., 2015; Rougier et al., 2016, 2018; Papale, 2018; Papale et al., 2021). Inevitably, all such analyses present substantial uncertainties.

In Japan, the largest eruption in the last 100 ka was Aso-4, which created the current Aso caldera (Fig. 1). It was responsible for the emission of about several hundred km³ dense rock equivalent (DRE) (Takarada and Hoshizumi, 2020; Aspinall et al., 2021; Rougier et al., 2022), and thus could be classed an M 8 “super-eruption”, or very close to that scale. The Aso caldera is located in the densely populated Kyushu Island (∼12.7 million people; October 2021), the volcano is currently active and poses significant challenges for risk mitigation (Tajima et al., 2017; Cigolini et al., 2018). One concern is the high runout distance of PDCs from this caldera which, in the past, have extended up to ≈160 km in some directions, e.g., for the Aso-4 eruption (Ono and Watanabe, 1985).

In this work, through first-order integral PDC models based on different assumptions, we evaluate the minimum volume and mass of a PDC generated at the Aso caldera that has the potential to affect any of five selected target “marker” sites (MS). These are populated cities or critical infrastructure positions within 160 km from the caldera, located in well-separated directions, i.e., MS1, MS2, MS3, MS4 and MS5 (see Fig. 1; Kato et al., 2016). For this, we compared four alternative kinetic energy models with variable input parameters based on structured expert judgment. Our constraints of PDC volume and mass were then compared with the main characteristics of the eruption of Aso-4. Note that our study is not dependent on accurately knowing the volumes of the Aso-4 or its individual flow units. Indeed the volumes have large uncertainties so we can use such estimates to constrain the range of volumes of interest in the modeling and then back check that the volumes with a high likelihood of reaching the observed distances are consistent with the known facts.

This work is organized in six sections. In Sect. 2 we describe the main characteristics of the volcanism of the Aso caldera. In Sect. 3 we briefly outline the first-order integral PDC models used in this analysis and the simplified strategy adopted to account for the shielding effects of topography. In Sect. 4, we present the main results derived from this study, including data elicited by expert judgment (Sect. 4.1), the minimum PDC volume and mass that would be able to reach the marker sites (Sect. 4.2) and the probability likelihoods that this threshold could be exceeded by PDCs with the characteristics of Aso-4 PDCs (Sect. 4.3). Then we present a brief sensitivity analysis of our numerical results (Sect. 4.4). Finally, in Sect. 5 we discuss the differences and limitations of the simplified models adopted to evaluate the resulting PDC runout distances and the probability of occurrence of caldera-forming eruptions on the scale of Aso-4. Supplement S1 contains technical notes on the derivation of the kinetic models and Supplement S2 is the email-log of the discussions that took place during the elicitation sessions.

2 Geological framework

2.1 The Aso caldera

The Aso caldera (dimensions 25 km N–S, 18 km E–W) is one of the largest calderas in Japan (Ono and Watanabe, 1985), with an area of 379 km², an average wall height of 500 m and a maximum wall height of 900 m (Matumoto, 1943). It is located in central Kyushu, southwest Japan (Fig. 1), where the Philippine Sea Plate subducts beneath the Amurian Plate (e.g., Petit and Fournier, 2005; DeMets et al., 2010).

The Aso caldera was formed by four successive large magnitude eruptions: Aso-1 (0.3 Ma, M 6.7), Aso-2 (0.14 Ma, M 6.7), Aso-3 (0.12 Ma, M 7.4) and Aso-4 (86.8–87.3 ka; M 8.1–8.4; Aoki, 2008; Takarada and Hoshizumi, 2020), where M=log ₁₀ [erupted mass (kg)] − 7 (Hayakawa, 1993; Pyle, 2000). Each of Aso-1–4 caldera cycles contains abundant dacitic to rhyodacitic pumice and volcanic ash, with mafic scoria included in the last stage of each eruption cycle (Hasenaka, 2016). The four caldera-forming eruptions of Aso produced >600 km³ DRE of pyroclastic deposits in total (>1200 km³ bulk), including hundreds of cubic kilometers of pyroclastic flow deposits (e.g., Matsumoto et al., 1991; Ono et al., 1977; Machida et al., 1985; Machida, 1999; Machida and Arai, 2003; Kaneko et al., 2007). In particular, Aso-1 erupted ≥40 km³ DRE (≥100 km³ bulk) of pyroclastic material, Aso-2 produced ∼22 km³ DRE (50 km³ bulk), and Aso-3 emitted about 100 km³ DRE (>150 km³ bulk; Matsumoto et al., 1991; Crosweller et al., 2012; Kaneko et al., 2015). Aso-4, which represents the largest eruption of Aso caldera, produced several hundred km³ DRE (please see next subsection for more details about Aso-4). The PDCs generated from these eruptions reached runout distances from the caldera of ∼30 km for Aso-1, ∼30 km for Aso-2, ∼70 km for Aso-3, and ∼166 km for Aso-4 (Ono and Watanabe, 1983; Takarada and Hoshizumi, 2020).

In addition to the caldera-forming events, Aso volcano has been the source of numerous smaller magnitude eruptions, both preceding and following Aso-4:

• Inter-caldera volcanism produced silicic pumice and lava flows, which may be related to the initial phases or “precursory” eruptions to caldera-forming events (Hasenaka, 2016). For instance, the deeply dissected stratovolcano Nekodake, which lies in the central eastern sector of Aso caldera, may have been formed about 150 ka BP, although recent work would suggest a more recent age (Shinmura et al., 2021).

• Post-caldera volcanism initiated soon after the Aso-4 eruption and has been dominated by smaller scale events with a wide compositional range (from basalt to rhyolite). At least 17 central cones were constructed during this period (Ono and Watanabe, 1985; Miyoshi et al., 2005; Miyabuchi, 2009), and many more edifices have been detected beneath the present central cones (Uto et al., 1994; Hoshizumi et al., 1997). Nakadake is the youngest and most active post-caldera volcano (see Fig. 1). The total thickness of scoria and ash-fall deposits is in the order of 100 m, including 36 silicic pumice-fall deposits. Most of the silicic eruptions produced deposits of <0.1 km³, while the Nojiri pumice (84 ka; ∼1 km³) and Kusasenrigahama pumice (30 ka; 2.2 km³), were one order of magnitude larger (Miyabuchi, 2011).

The total volume of post-caldera tephra has been calculated by Miyabuchi et al. (2003, 2004) and Miyabuchi (2009, 2011) as 18.1 km³ DRE. The total volume of the edifices of the post-caldera cones (including buried edifices) is estimated at ∼112 km³ DRE (Komazawa, 1995) of which approximately 24 km³ DRE are exposed at the surface (Miyoshi et al., 2012). Thus, the average rate of magma discharge during the post-caldera stage is about 1.5 km³ kyr⁻¹, similar to the rate at other Japanese quaternary volcanoes (Miyabuchi, 2009). However, an analysis of these published volumes shows that the post-Aso-4 eruption rate has experienced a significant waning trend. The bulk (∼115 km³) of the post-Aso-4 magma erupted before approximately 70 ka, suggesting elevated eruption rates of 7 km³ kyr⁻¹ between 70–87 ka. During the period 22–50 ka approximately 7 km³ of magma erupted at a rate of ∼0.3 km³ kyr⁻¹, with silicic magma slightly more abundant than mafic magma. Since 22 ka, predominantly (98 %) mafic magma has erupted at a rate of approximately 0.1 km³ kyr⁻¹.

2.2 The Aso-4 eruption and its products

The Aso-4 eruption (86.8–87.3 ka) is the largest eruption during the last 100 kyr in Japan. The PDC deposits of this eruption reached ∼166 km to the NNE from the source, being widely distributed in the northern and central zones of Kyushu Island and also in the SW part of Honshu, invading the Yamaguchi prefecture (>160 km from the source; Matumoto, 1943; Watanabe, 1978; Ono and Watanabe, 1985; Takarada and Hoshizumi, 2020).

The pyroclastic flow deposits of Aso-4 have been divided into two cycles, each one characterized by a progression from silicic to more mafic magma, with evidence of magma mixing (Lipman, 1967; Watanabe, 1979; Ono and Watanabe, 1983; Kaneko et al., 2007; Ishibashi et al., 2018). The ash of Aso-4 consists of rhyodacitic glass shards with brown hornblende and orthopyroxene phenocrysts (Machida and Arai, 1983). Plagioclase and magnetite phenocrysts occur in all ejecta from the Aso-4 eruption as well (Kaneko et al., 2007). The pyroclastic flows of Aso-4 also contain abundant lithic fragments (Ono et al., 1977) with regionally differing lithologies. This lateral heterogeneity of lithic fragments may suggest that the magmas of the Aso-4 eruption were erupted from either a ring fissure or multiple vents (Kaneko et al., 2007), as has been inferred for other explosive super-eruptions (e.g., Costa and Martí, 2016); this may also be explained by surface erosion across different lithic populations.

Takarada and Hoshizumi (2020) presented estimates of the volume and mass of PDC and fallout deposits of Aso-4, by evaluating the magnitude (M 8.1–8.4). They estimated a volume of Aso-4 tephra fall deposits of 240–370 km³ DRE (6.0–9.3×10¹⁴ kg), while the volume of Aso-4 PDC deposits was estimated at about 225–590 km³ DRE (5.6–14.8×10¹⁴ kg). These estimates implied a total volume of the Aso-4 eruption of 465–960 km³ DRE (1.2–2.4×10¹⁵ kg).

Nevertheless, the volume of magma erupted in Aso-4 and the uncertainty in volume estimates remain open questions. Aspinall et al. (2021) undertook a preliminary quantification of this uncertainty by combining several different volume estimates of PDC and fallout deposits of Aso-4, through a Bayesian belief network approach. They defined composite volumes, determined by weighted combinations of different estimates of the various products, including those from Takarada and Hoshizumi (2020), and hypothesized a lower bulk volume for the ash fall deposit and therefore a smaller total magnitude with 90 % confidence (i.e., M 7.9–8.2). In a subsequent related assessment of the Aso-4 tephra fall deposit, Rougier et al. (2022) applied advanced statistical methods and re-evaluated its bulk volume to have been in the range 220–370 km³, which equates to the range 95–160 km³ DRE (2.2–3.7×10¹⁴ kg). As the volumes of co-ignimbrite tephra fall deposits and PDC deposits are typically comparable to one another (Sparks and Walker, 1977), the reduced estimate tephra fall deposits also suggests a somewhat reduced PDC volume.

The exposure area of Aso-4 PDC deposits is about 1000 km², but Takarada and Hoshizumi (2020) evaluated that the original area covered by these deposits was about 34 000 km². However, in some directions the runout distance of Aso-4 pyroclastic flows largely exceeded 130 km, while in others the deposits could not be traced that far. The Aso-4 pyroclastic flow deposits have been likely conditioned by the presence of topographical obstacles such as the Kyushu Mountains.

As discussed above, there is much uncertainty about the total PDC volume and consequently the volumes of individual PDCs generated by the Aso-4 eruption. However, from the perspective of our study these uncertainties are not important as all studies agree that the total PDC volume is likely to significantly exceed 100 km³, while the volume of the largest and most extensive individual PDC emplacement units is also likely to be a significant fraction of the total volume.

3 Methods

3.1 Box model integral formulations

Our analysis relies on the implementation of four different versions of the box model integral formulation for axisymmetric gravity-driven particle currents, based on the pioneering work of Huppert and Simpson (1980) and with theory detailed in Bonnecaze et al. (1995) and Hallworth et al. (1998). We focus on models with analytical solutions, adopting input ranges based on expert judgment (see Sect. 4.1). This enables us to utilize a very fast model inversion approach in the uncertainty quantification process. We note that we are not using “reduced” models (i.e., statistical surrogates, e.g., Rutarindwa et al., 2019; Yang et al., 2020). Further details on the physical equations we adopted as well as the mathematical expression of the analytical solutions, can be found in the Supplement S1. In particular, we focus on two different physical formulations.

3.1.1 Model 1: rock avalanche dynamics with constant stress over the flow basal area

This model for energy dissipation is described in Dade and Huppert (1998) and assumes that the entire amount of solid material falls from a prescribed height and expands laterally while preserving its initial volume. In this model the elicited input parameters are: (a) collapse height (H), (b) flow density (ρ_c), and (c) equivalent stress (τ, one third of the constant stress value, see Supplement S1). The dynamics of a similar assumption on the basal stress has been further explored in Kelfoun et al. (2009) and Kelfoun (2011), using depth-averaged models.

3.1.2 Model 2: density current dynamics with particle deposition

This model is described in Dade and Huppert (1995) for the simulation of oceanic turbidity currents. Dade and Huppert (1996) adopted it to simulate emplacement of the large-scale AD 232 Taupo ignimbrite (New Zealand; Holdaway et al., 2018). The input parameters of model 2 are: (a) initial solid fraction (ϕ₀), (b) velocity of settling of the solid particles (w_s), (c) density of solid particles (ρ), (d) density of ambient air (ρ_a), (e) density of interstitial gas (ρ_i), and (f) the Froude number (Fr). All the input ranges of these parameters were elicited with the exception of the Froude number, which was sampled uniformly in [1.0, 1.2]; this range captures expected conditions for PDC (Esposti Ongaro et al., 2016).

In contrast to Dade and Huppert (1996), in this study we assume the instantaneous release of a fixed volume of pyroclastic material, as in Neri et al. (2015), Bevilacqua (2016), Bevilacqua et al. (2017) and Aravena et al. (2020, 2022). This assumption is considered reasonable because the time scale of the eruption of the PDC volume is shorter than the emplacement time. As the model accounting for the full total grain size distribution (TGSD) does not have analytical solutions (e.g., Biagioli et al., 2019; Tadini et al., 2021a) and, however, a TGSD of the studied PDCs would not be easy to constrain, we assumed monodispersed solid particles with an effective diameter. In particular, we adopted the Sauter diameter, which is the particle diameter where the volume/surface area ratio is the same as that of the particle distribution (Crowe et al., 1998; Neri et al., 2022), to provide a reasonable approximation to the dynamics of the full TGSD.

In this work we also considered three variants of Model 2:

• a.

  Model 2a. This variant includes interstitial gas, which is thermally buoyant with respect to the surrounding cold air. The flow stops propagating when the solid fraction ϕ(t) becomes lower than a critical value ϕ_cr, and the remaining mixture of gas and particles lifts off possibly generating a ignimbrite cloud. Note that the thermal properties and flow volume remain constant for the duration of the flow.

• b.

  Model 2b. The equations of this variant are equivalent to model 2a, but an alternative input of w_s is adopted. Instead of being elicited, this range is based on the results of particle terminal velocity derived from ash fallout studies at the scale of the Sauter diameter for analog flows (Armienti et al., 1988; Bonadonna and Phillips, 2003; Dioguardi et al., 2017). These atmospheric flows are Mt. St. Helens (Costa et al., 2016), the Campanian Ignimbrite (Costa et al., 2012; Martí et al., 2016), and the Youngest Toba Tuff (Costa et al., 2014). Moreover, this variant implements an alternative input range of ϕ₀, based on mass discharge rate (MDR) modeling (Costa et al., 2018). Therefore, this formulation represents the hypothetical scenario of a more dilute ash cloud than the elicited flow.

• c.

  Model 2c. This variant assumes a “cold” interstitial gas equivalent to ambient air. Thermal buoyancy effects are absent, and the flow stops when ϕ(t=0 and the entire initial solid fraction has been deposited. Under this assumption, a longer runout than model 2a is expected.

Several recent papers have focused on pyroclastic flow emplacement by more sophisticated modeling; e.g., thermodynamics modeling of cooling effects and air entrainment effects (e.g., Bursik and Woods, 1996; Fauria et al., 2016), two layer systems (Burgisser and Bergantz, 2002; Doyle et al., 2007; Kelfoun, 2017; Valentine, 2020), the development of dense and dilute regimes within the same flow (Esposti Ongaro et al., 2011, 2020; Kelfoun and Gueugneau, 2022; Neri et al., 2022), the effects of pore pressure on basal friction (Roche et al., 2021; Aravena et al., 2021), and the build-up of coherent turbulent structures and gravity waves (Lube et al., 2020; Brosch et al., 2021). However, since more advanced models require more detailed data to fit additional input parameters, simplified models are suitable for probabilistic modeling in uncertain frameworks, like the present Aso-4 case study, especially for capturing extreme probabilities.

For instance, some of these models adopt a two layer approach in which a concentrated basal flow is overlaid by a dilute turbulent cloud. The two layer exchange mass during emplacement through particle settling and entrainment processes, making such models computer-intensive and currently unsuitable for probabilistic modeling. Models 1 and 2 can be seen as end members that treat the two layers separately and so are likely to yield a wide range of run-out distances that include those found in a hybrid two layer model. In Sect. 5 we further discuss the interpretation of our results in terms of such models.

3.2 Mass inversion strategy and simplified testing of topographic shielding effects

In a first step, our numerical simulations do not consider the possible shielding effect of topographic obstacles. By using the elicited probability distributions of model inputs and the analytical solutions of the kinetic energy models, for any given collapsing volume, we estimate the conditional probability for the occurrence of a PDC with runout distance on a flat topography equal or larger than the distance between the source and each marker site (i.e., ∼130 km for MS1, ∼135 km for MS2, ∼145 km for MS3, ∼115 km for MS4, and ∼120 km for MS5; see Fig. 1).

As the minimum volumes have a different density in the different models, we computed the equivalent mass to set up a more consistent comparison between them. In particular, model 1 is single-phase and the calculated volume represents the bulk solid material at the elicited flow density. On the other hand, ,models 2a–c are two-phase and include solid particles and interstitial gas, but the volume calculated only represents the solid phase, with a density that is also elicited and generally greater than the flow density of model 1. Nevertheless, because of the possible topographic effects in the propagation of PDCs, the simulation of 115–145 km runout in the absence of topographic effects is not a sufficient condition to assess the likelihood of marker site incursions for those distances from the source. Therefore, in a second step, we tested the inclusion of the topography in numerical simulations, using the program BoxMapProb (Aravena et al., 2020), which can be found in https://github.com/AlvaroAravena-/BoxMapProb (last access: 18 August 2022); this approach adopts the formulation of model 2c. To do this, we followed the “energy conoid” approach, based on the assumption of nonlinear, monotonic decay of flow kinetic energy with distance adopted by Neri et al. (2015) and Bevilacqua et al. (2017). Note that we did not follow the tree-branching algorithm in Aravena et al. (2020) because strong channelizing effects were not envisaged for the Aso case (see Supplement S2).

Figure 2 Probabilistic PDC invasion maps for four sets of simulations (N=1000 each run) with different values of collapsing volume and vent positions sampled uniformly within Aso caldera. The other input parameters are derived from expert judgment (see Table 1). These simulations were performed using the program BoxMapProb (Aravena et al., 2020; excluding tree-branching effects), which adopts the formulation of model 2c. Dotted gray circles mark the positions of the five marker sites. Contour lines show the invasion probability values in the lower left box. SRTM – DEM 1 arcsec 30 m (NASA, 2013).

In Fig. 2 we show some examples of PDC invasion probability maps as a function of the DRE volume of pyroclasts. We compared the kinetic energy of the current front and the potential energy associated with the obstacles encountered. Our approach is most sensitive to the shielding effect of the topography close to the marker sites, and not to the large-scale topography around the source site (Aso caldera), that might also significantly affect the flow dynamics (e.g., Todesco et al., 2006; Esposti Ongaro et al., 2020). Therefore, according to the literature sources about Aso-4, we only simulated the mass of PDC outflowing from the caldera. As the energy calculation is performed axisymmetrically along every radial direction, we randomly varied the center of propagation inside the Aso caldera (see Fig. 1) and averaged the results to include uncertainty related to the geographical coordinates of the source.

Figure 3 shows the results of a set of complementary simulations performed in this way for marker sites MS1–MS3. For example, we estimated that the inputs able to produce a runout distance of 150 km on a flat surface ensures that ∼50 % of the simulated PDCs invade the marker site MS1 when the topography is considered (Fig. 3a). Using these indications to emulate the effect of topography, for MS1 we included in this work the results associated with a runout distance of 195 km, i.e., the value at which 95 % of the simulations reach the marker site. Similarly, we determined that sets of inputs able to produce runout distances of 235, 225, 155 and 170 km on a flat surface will ensure 95 % of the simulations reach MS2, MS3, MS4 and MS5, respectively, when topography is considered (Fig. 3b and c). See Sect. 5.2 for a discussion on the possible limitations of this strategy.

Figure 3 Cumulative curves PDC invasion probability of the marker sites (a) MS1, (b) MS2, (c) MS3, for a set of simulations performed using the program BoxMapProb (Aravena et al., 2020) as a function of the parameter L_max, which represents the runout distance on a flat surface. To perform these simulations, vent positions are uniformly sampled within the Aso caldera, the volume of pyroclasts is uniformly sampled between 10⁹ and 10¹² m³, and the other input parameters are derived from expert judgment (see Table 1 – 5000 samples). We observe that L_max of 195, 235 and 225 km ensure that most of the simulations (>95 %) invade the MS1, MS2, MS3, respectively.

Finally, note that for simplicity in all cases we modeled the sea as a flat topography (Neri et al., 2015; Bevilacqua et al., 2017). However, it is a debatable point whether model 1 flows can go over the ocean. In fact, the ability of a flow or surge to pass over water is a complex matter, and it will depend on many factors, such as the angle of entrance into the water, flow speed, density, and temperature (Cas and Wright, 1991; Carey et al., 1996; 2000; Allen and Cas, 2001; Dufek and Bergantz, 2007; Dufek, 2016). These effects are hardly represented in simplified kinetic models like those adopted in this study but, in our opinion, the results obtained provide a first approximation of the distal runout impacts of such catastrophic eruptions. From a safety perspective, the assumption that PDCs can travel over water is both supported by observations of past (smaller magnitude) eruptions and could be exceedingly conservative, in the sense that the assumption will not underestimate risk.

4 Results

4.1 Input ranges based on expert judgment

We based our input range estimation on structured expert judgment (Cooke, 1991; Aspinall, 2006). The elicitation solutions of models 1 and 2a–c are indicated in Table 1, where we also include the modified input ranges for model 2b based on MDR modeling and Sauter diameter of analogs (see Sect. 3.2). Note that volume is not considered an input parameter but it is varied to calculate the probability of any volume chosen in a wide range to reach the distance of each marker site. For judgment aggregation, we implemented the equal weights combination rule (Bevilacqua et al., 2015; Tadini et al., 2017, 2021b, 2022b). We did not apply performance-based scores because of the relatively small number of experts participating and because the overheads and time demands involved in implementing a formal elicitation protocol were not warranted in this case. Two elicitation sessions were organized remotely by sending out the questionnaires and collecting the responses by email. Supporting information S2 demonstrates that a comprehensive exchange of emails took place remotely, in which the models adopted, the variables elicited and their parameterizations were thoroughly discussed.

Table 1 Elicitation solutions for model 1 and model 2 parameters, and modified input range of model 2b based on mass discharge rate (MDR) modeling and Sauter diameter of analogs (see Sect. 3.2). The reported values express the percentiles of the probability distribution obtained by equal weights pooling of experts' judgments (also called the solution decision maker DM).

^a These results approximately impose the input values in the range (min, median) assessed by the joint DM. ^b Froude number was sampled uniformly in [1.0, 1.2].

More detailed input schemes and calibration approaches could be attempted in further research (Bevilacqua et al., 2019; Patra et al., 2020; Aravena et al., 2022), but these approaches would not be straightforward because of the lack of well-preserved distal deposits. Note that we assume the model inputs consist of an array of independent variables. A discussion of this limitation is provided in Sect. 5.3 and future research might explore the effects of possible correlations between them, especially on the computation of extreme model values. Another possibility is to develop simplified two layer models amenable to probabilistic modeling calibrated by the more computer-intensive numerical models to allow simplified parameterizations of mass exchange between the two layers.

4.2 Minimum PDC volume and equivalent mass to reach the marker sites

The volume inversion results on a flat topography are presented in Tables 2, S3 (in the Supplement), S5, S7, S9 and express, for a given percentage P, the minimum PDC volume that has such a probability P of reaching MS1, MS2, MS3, MS4 or MS5, based on values of runout distance equal to 130, 135, 145, 115 and 120 km, respectively. In all the cases, we converted the minimum volumes into PDC mass by considering the elicited values of density (ρ_c and ρ for model 1 and models 2a–c, respectively; see Table 1). The estimates of the minimum PDC mass (MinMass) able to invade each marker site, obtained with the four described models on a flat topography, are also included in Tables 2, S3, S5, S7, S9.

Table 2 Numerical results of the minimum PDC volume and mass needed to reach MS1, with and without consideration of topographic effects (see Sect. 4.4). The density is different in the two models 1 and 2 (see Table 1). Estimates are based on 100 000 statistical samples.

^* Modified ϕ₀ and w_s based on MDR modeling and Sauter diameter of analogs (see Table 1).

This strategy to calculate MinMass, which is a conservative choice, according to model 2a produces values of $\sim \mathrm{17.4}\times {\mathrm{10}}^{\mathrm{12}}$ kg (DRE equivalent 10.2 km³) if we consider a 5 % probability to reach the distance of MS1, $\sim \mathrm{95.9}\times {\mathrm{10}}^{\mathrm{12}}$ kg (DRE equivalent 55.9 km³) if we consider a 50 % probability and $\sim \mathrm{390}\times {\mathrm{10}}^{\mathrm{12}}$ kg (DRE equivalent 226 km³) if we consider a 95 % probability, always considering a flat topography (Table 2). For MS2 and MS3, the associated values of MinMass according to model 2a are 11 % and 33 % higher than those of MS1, respectively. The estimates for MS4 and MS5 (Tables S7 and S9) are 28 % and 19 % lower, respectively.

Tables 2, S3, S5, S7, S9 also include the results based on larger values of runout distance, i.e., 195, 235, 225, 155 and 170 km for MS1, MS2, MS3, MS4 and MS5, respectively, obtained by using again the analytical solutions of kinetic energy models on a flat surface. This increase in the marker runout distance is introduced so as to emulate the effect of the topography, i.e., runout distance reduction derived from the effect of distal obstacles (see Sect. 3.2). Otherwise, the analytical solutions of kinetic energy models would have been invalidated by numerically considering the topographic shielding effects in each simulation.

These results are less conservative than those presented previously and indicate that the minimum PDC mass able to affect the marker sites for model 2a, is ∼3-times greater than the previous estimates based on shorter values of runout distance. In particular, for MS1 the MinMass is $\sim \mathrm{51.2}\times {\mathrm{10}}^{\mathrm{12}}$ kg (DRE equivalent 30.1 km³) with 5 % probability of incursion, $\sim \mathrm{283}\times {\mathrm{10}}^{\mathrm{12}}$ kg (DRE equivalent 165 km³) if with 50 % probability and $\sim \mathrm{1150}\times {\mathrm{10}}^{\mathrm{12}}$ kg (DRE equivalent 668 km³) if with 95 % probability. The estimates of the minimum PDC mass needed to invade MS2 and MS3 (Tables S3 and S5) are 65 % and 46 % higher than those of MS1, respectively. The estimates for MS4 and MS5 (Tables S7 and S9) are 46 % and 31 % lower, respectively.

Note that our four approaches are based on significantly different physical assumptions and therefore their results are presented separately. In fact, model combination schemes could be attempted, but our expert elicitation did not indicate any model preference (see Supplement S2). Bayesian model averaging (Bevilacqua et al., 2017, 2018) was not practicable in this case because of the lack of enough well-preserved distal outcrops. A more detailed discussion of all the models and the associated MinMass inversion results is provided in Sect. 5.1. In particular, model 1 requires ∼3-times greater MinMass values to reach the same runouts as model 2a with 50 % probability, while model 2b has ∼3-times smaller values and model 2c requires ∼30 % smaller mass values than model 2a.

We finally note that these MinMass estimates represent a minimum threshold required to just reach the marker sites; however, the distal impact of these PDC could differ radically between the models, i.e., from PDCs with high particle concentrations and high dynamic pressures to co-ignimbrite flows more similar to dilute ash clouds (see Sect. 5.1).

4.3 PDC invasion probability at marker sites based on the documented eruption size of Aso-4

In Fig. 4 we show the cumulative distribution of the mass required to reach marker sites MS1–MS3 (MinMass) according to the four models described in the previous sections with (right panel) and without (left panel) consideration of the effect of the topography. Figures 5 and S11 (in the Supplements) how the probability density functions of MinMass with and without topographic effects, respectively. In these figures, we also include very uncertain estimates of the mass of the PDCs deposited during the largest caldera-forming eruption at Aso, i.e., Aso-4. For Aso-4 we used the range of PDC mass estimated by Takarada and Hoshizumi (2020), considering only the mass of PDC over-flow outside the caldera, i.e., 420–1180×10¹² kg, mean value 780×10¹² kg. These values of mass correspond to 168–471 km³ DRE, mean value 314 km³ DRE.

Figure 4 Cumulative distributions of the variable MinMass, calculated using models 1, 2a, 2b, 2c. MinMass represents the mass of pyroclasts in a PDC flow required to invade the different marker sites. (a, c, e) are related to maximum runout distances of 130, 135 and 145 km, respectively (the distance between Aso caldera and the marker sites MS1–MS3). (b, d, f) are related to maximum runout distances of 195, 235 and 225 km, respectively, allowing the flow to overcome possible topography shielding effects near the marker sites. Estimates of the mass associated with the PDCs produced during the largest caldera-forming eruption of Aso are included: Aso-$\mathrm{4}/\mathrm{10}$, plausibly representing the Aso-4T unit, in green; Aso-4, total PDC overflow, in yellow (Takarada and Hoshizumi, 2020). Plots are based on 100 000 statistical samples.

Figure 5 Probability density functions of the variable MinMass, calculated using models 1, 2a, 2b, 2c. MinMass represents the mass of pyroclasts in a PDC flow required to invade the marker sites MS1–MS3, related to maximum runout distances allowing the flow to overcome possible topography shielding effects near the marker sites. Estimates of the mass associated with the PDCs produced during the largest caldera-forming eruption of Aso are included: Aso-$\mathrm{4}/\mathrm{10}$, possibly representing the Aso-4T unit, in green; Aso-4, total PDC overflow, in yellow (Takarada and Hoshizumi, 2020). Plots based on 100 000 statistical samples.

As the mass of the longest runout PDC of Aso-4 may have been a significant fraction of the PDC flows outside of the caldera, we first assume that the total volume of the long runout PDC is equal to the total outflow PDCs of the eruption. In fact, by considering the entire Aso-4 mass we assume that a single flow contains most of the erupted mass of the outflow deposits, and that this mass was erupted within a time scale less than the PDC emplacement time. This assumption is further discussed in Sect. 5.1.

However, Aso-4 contained several emplacement units and the longest runout PDC could be only one of these, e.g., Aso-4A, Aso-4B, and Aso-4T (Lipman, 1967; Watanabe, 1979; Kaneko et al., 2007). These units contain evidence of at least 10–20 PDCs in total (Takarada and Hoshizumi, 2020). In particular, Aso-4T (Tosu unit) is the most widely distributed low-aspect-ratio ignimbrite unit within Aso PDCs, but its mass remains unknown (Suzuki-Kamata and Kamata, 1990). Therefore, we tentatively also tested a mass of $\mathrm{1}/\mathrm{10}$ of the total PDC overflow outside the caldera, i.e., 42–118×10¹² kg, mean value 78×10¹² kg (17–47 km³ DRE, mean value 31 km³ DRE), plausibly representative of Aso-4T mass. Note that this assumption is not aimed at constraining the actual mass of Aso-4T, but it provides an illustrative example of a PDC that does not contain most of the erupted mass of the outflow deposits.

The conditional invasion probabilities obtained according to all models and for each marker site are listed in Tables 3, S4, S6, S8, and S10. In particular, these results indicate that a PDC with a mass of $\mathrm{1}/\mathrm{10}$ of the total PDC overflow of Aso-4, without considering the shielding effects of topography, would likely affect MS1 according to model 2a, with a conditional probability of 23 %–56 %. After including topographic shielding effects (see Sect. 3.2), this conditional probability is reduced to 3.1 %–21 %. Slightly lower conditional invasion probabilities are computed for MS2 and MS3, with values of 0.5 %–9.3 % and 0.9 %–12 %, respectively, slightly higher for MS4 and MS5, i.e., 11 %–40 % and 7 %–32 %. In contrast, the conditional incursion probabilities based on the total outflow PDCs of the eruption are significantly higher, i.e., 46 %–100 % with 90 % confidence, depending on the marker site. Note that Table 3 also shows the results based on the volume estimates of Aspinall et al. (2021), which are lower than those based on Takarada and Hoshizumi (2020) but still significant, i.e., 49 %–87 % instead of 66 %–96 % for TS1.

Table 3 Numerical results of the probability that a PDC derived from a caldera-forming eruption similar to Aso-4 reaches MS1, with and without topographic effects. For each scenario, we present the values of the cumulative curves displayed in Fig. 4 at the central point of the variation range of the PDC mass, while between parentheses we include the results at the extremes of these variation ranges. We test the mass of PDC overflow outside of the caldera of Aso-4 and $\mathrm{1}/\mathrm{10}$ of that estimate, tentatively representing the Aso-4T unit.

^a TE: Topographic effects. ^b Modified ϕ₀ and w_s based on MDR modeling and Sauter diameter of analogs (see Table 1).

4.4 Sensitivity analysis

We use a Bayesian belief network (BBN; e.g., Hincks et al., 2014) related to each model to characterize uncertainties on the MinMass and MinVol estimates (Figs. 6 and S12), including the existing correlation coefficients (marked with linking arc arrows) between the input variables and the calculated MinMass and MinVol. This analysis was performed with the UNINET software (Ababei, 2016). The histograms in each BBN also illustrate the shape of the uncertainty distributions for the intermediate variables (g_c, g_p and ϕ_cr) and for the minimum mass (MinMass) and volume (MinVol) needed to reach the first marker site, based on the elicited distributions of the relevant input variables.

Figure 6 Sensitivity analysis of the evaluation of functional variable MinMass in the different models used in this work. Gray boxes indicate the input parameters, and arrow labels indicate correlation coefficients in relation to model inputs and functional variables. BBN based on 100 000 statistical samples.

Larger correlation values highlight variables to which our probability calculations are most sensitive. Negative correlation means that an increase of that input parameter decreases the MinMass (or MinVol) required to reach the marker site. Positive correlations produce the opposite effect. The correlation coefficients are based on the elicited input distributions and on the 130 km maximum runout distance, i.e., valid for MS1 with no consideration of topographic effects. We tested other runout thresholds but observed negligible differences. In model 1, both the height of collapse and the stress coefficient are significantly relevant, the former negatively, the latter positively, while density does not affect the MinMass estimate.

In all versions of model 2, a significant effect on modeling results comes from the sedimentation velocity and from the initial concentration of solid particles, both positively correlated with MinMass, while the solid density and the Froude number are less relevant. Further research aimed at developing better constraints on the uncertainty distributions of the most influential input parameters would likely better constrain the resolution of probability results.

5 Discussion

5.1 First-order integral PDC models

Most recent models of PDCs have highlighted the stratified nature of such flows (Valentine, 2020; Roche et al., 2021; Neri et al., 2022). Nevertheless, in our study we just focused on two simple end member types, suitable for a first-order quantitative assessment of the hazard (see Supplement S2). Although our simplified models do not include the development of dense and dilute regimes within the same flow (Burgisser and Bergantz, 2002; Doyle et al., 2007; Kelfoun and Gueugneau, 2022), models 1 and 2 may represent the two end members, concentrated and dilute flows, respectively, thus shedding some light on more complex assumptions. For example, our results are consistent with Roche et al. (2021), where their model includes both dense and dilute regions, and interactions. In fact, by analyzing data from well-documented explosive eruptions those authors observed that at a given mass discharge rate dilute PDCs have runouts generally longer than those of their concentrated counterparts.

From a risk point of view, our analysis solely derives probability estimates of PDCs reaching the marker sites. A PDC represented by model 2 may do so as a dilute and low temperature flow. Structural risks from these relatively dilute and “cold” PDCs may be similar to risks from tephra fall and would differ radically from the greater thermal and mechanical loads that might be associated with dense PDCs.

Note that for simplicity, rather than modeling a mass discharge rate for unknown eruption duration, we assumed an instantaneous release of a fixed volume of pyroclastic material, as in Neri et al. (2015), Bevilacqua et al. (2017) and Aravena et al. (2020). Further research could follow a similar probabilistic approach by inverting minimum mass discharge rates from the runout distance rather than the MinMass, but this strategy would introduce challenging questions on how volume, eruption duration, runout and flux are quantitatively related. Further model development could also include eliciting the proportion of magma that is erupted as a single volume within a time scale less than the emplacement time.

5.2 Simplified testing of topographic shielding effects

The strategy that we adopted to emulate the effect of the topography on our results relies on the energy conoid approach. Thus, the energy calculation is calculated axisymmetrically along every radial direction, and in some cases the inferred topographic effect might be too intense, as flows can actually converge around obstacles. Also, this approach cannot include partial blocking by topographic obstacles, and returning waves (e.g. Todesco et al., 2006; Esposti Ongaro et al., 2020). Although we reduced the spurious effects related to the exact geographical coordinates of the PDC source by randomly varying the center of propagation inside the Aso caldera (see Fig. 1), only a 2D or 3D model would be capable to capture the more complex source-dependent features of the fluid dynamics.

Whereas the MinMass estimates obtained on a flat surface are probably underestimated, i.e., conservative, our simplified strategy to account for topographic shielding effects might instead overestimate MinMass in some cases. In fact, our approach invokes runout thresholds such that the inputs that enable the PDC to reach a marker on a flat surface will also project the marker site invasion likelihood with 95 % confidence when the topography is considered by the simplifying expedient adopted. Therefore, lower runout thresholds would produce intermediate MinMass values between those obtained on a flat surface and those that account for topographic effects with 95 % confidence. However, these intermediate values would be associated with a lower chance of marker site incursion when the topography is considered and a relatively low kinetic energy budget at the marker site. These bounding effects may be further studied through statistical methods that analyze the edge of the inundated region once 2D or 3D models are applied (Hyman et al., 2019).

5.3 Expert judgment approach and the input sampling strategy

Our results are valid under the assumption that the BBN samples drawn from the elicited inputs are independent, but this might not be valid in some cases. For example, in model 1 the combination of high collapse height H and low equivalent friction τ could be problematic and, likewise, high density is hard to reconcile with a high collapse height (see Table 1 and Supplement S2). Similar questions may arise about other combinations of input values. While the median values of our estimates are not likely to be significantly affected by such statistical dependence effects, the application of interparameter correlations or probabilistic copula techniques would become really important for making robust probabilistic estimates of extreme events, say with less than 5 % probability (Bamber et al., 2019; Barons and Aspinall, 2020; Werner et al., 2021). Further model inversion analysis could evaluate which extreme events are in fact drawn from combinations of independent parameters that are physically debatable.

An indirect and interesting finding of our study is the identification of input variables that are most difficult to constrain. For example, in model 1, while for dense debris avalanches or dome collapse PDCs the variable H can be well defined, some of the experts found it difficult to constrain height H values for PDCs generated from explosion columns. In fact, the dense underflow is thought to be formed by segregation of particles and so it is essentially unknown how initial potential energy is partitioned between the dense and less dense components of the PDC: the effective height H is not likely to be related to overall column height. There would be a similar issue with the H variable in the energy cone formulation, which was not implemented in this study (Bevilacqua et al., 2022; Aravena et al., 2022).

As the emplacement of a low-aspect-ratio ignimbrite has never been observed, we tested several models and a wide range of input parameters. Our purpose was to capture the uncertainty affecting the plausible dynamics of similar flows, and constrain the minimum mass (MinMass) required to reach marker site distances; we did not endeavor to produce an optimized reconstruction of parameters for a flow identical to Aso-4 PDCs, as such. Nevertheless, one aspect that could be addressed in further research is calibration with large scale PDCs which originated from other volcanoes, e.g., Taupo, AD 232, (Dade and Huppert, 1996, 1997; Wilson, 1985a, b, 1997), or Krakatau, 1883 (Carey et al., 1996, 2000), where volumes as well as runouts are relatively well-established. For example, the initial solid fraction of the PDC of Taupo, AD 232, was estimated as 0.3 % by Dade and Huppert (1996), by model fitting. Nevertheless, while discussing density above the ground, Wilson (1985b, 1997) observed tens percentage solid. Although they both do not consider that column collapse creates highly variable solid fractions, these differences outline the dissimilar conditions ranging from density current dynamics with particle deposition to one possibly similar to a rock avalanche with constant friction.

It might be feasible to repeat our analysis approach to find the runout model parameter PDFs which best fit those observations. However, since these parameters would be related to different source conditions, they might need detailed discussion and perhaps additional expert judgment sessions before being invoked as counterparts to this case study. Similarly, we note that for our Aso-4T analysis, while we did not elicit the volume, we relied on the same runout model parameters that were elicited for the total Aso-4 PDC overflow outside of the caldera; eliciting and setting different parameters in the model might improve the robustness of the Aso-4T probability estimates.

5.4 Probability of occurrence of large-scale events in Japan

The PDC invasion hazard at the marker sites is ultimately related to the recurrence rate of catastrophic caldera-forming eruptions at the Aso volcano. We estimated PDC invasion probabilities that are all conditional on the occurrence of large-scale eruptions. In practice, we are dealing with extremely low probabilities of occurrence, but extremely high impact events. Although this aspect is beyond the objectives of this paper, a full magnitude frequency model for the Aso caldera could be integrated into our PDC invasion probability estimates to obtain hazard assessments at 100 or 1000 years. We briefly review the main historical results available.

For eruptions > M 7, global studies (Decker, 1990; Simkin, 1993; Mason et al., 2004; Deligne et al., 2010; Brown et al., 2014) suggest a minimum global average recurrence rate of ∼1–2 times per 1 kyr (Newhall et al., 2018). Once record incompleteness is corrected for (Kiyosugi et al., 2015; Rougier et al., 2016), the value approaches 2 events per 1 kyr. At individual volcanoes, >M 7 recurrence rates range from 1 per 10 kyr to 1 per 100 kyr (Newhall et al., 2018). For eruptions > M 8, Rougier et al. (2018) estimated a return period of 17 ka (95 % confidence interval 5.5–48 ka) from analysis of the LaMEVE database (Crosweller et al., 2012).

More recently, Aspinall et al. (2021), implementing a comprehensive stochastic uncertainty analysis, estimated the probability in the next 100 years of an eruption from the Aso volcano on the same scale as Aso-4. They set up a Bayesian belief network informed by multiple strands of evidence from volcanology, petrology, geochemistry and geophysics, together with reviews of published data, models, and inputs from expert elicitation. Based on the current and expected state of the Aso volcanic system in the short term, they assessed the probability of a similar eruption ≥M8 at the Aso volcano in the next 100 years to be no greater than 10⁻⁸– 10⁻⁹ (i.e., 10⁻¹⁰–10⁻¹¹ on an annualized basis).

6 Conclusions

We describe a new method for estimating the probabilities that distal geographic locations might be invaded by a large pyroclastic density current (PDC) from a high magnitude explosive eruption by expressing the problem in terms of PDC flow mass and related uncertainties. The ranges of key parameters were obtained by expert elicitation performed in two remotely organized sessions.

Our analysis relied on different versions of the box model integral formulation for axisymmetric gravity-driven particle currents (Huppert and Simpson, 1980), and a simplified testing of topographic shielding effects. We focused on models based on an underlying assumption of instantaneous release of pyroclastic material with fixed volume of pyroclastic material, which permit analytical solutions, enabling us to utilize a very fast model-inversion approach. In particular, we adopted the energy conoid approach to generate regional PDC invasion maps, based on the comparison of the mass-dependent kinetic energy of the flow with the potential energy needed to overcome topographic obstacles along radial directions. We considered four different models based on two physical sets of assumptions:

• a.

  Model 1 assumed that the entire amount of solid material originates from a prescribed height above the volcano and flows as a granular current slowed down by a constant stress, i.e., a rock avalanche dynamics model.

• b.

  Models 2a–c were two-phase density current dynamics with particle deposition and thermally buoyant gas. Three variants of this approach were tested, i.e., the first implemented the expert judgments, the second tested the hypothetical scenario of a more dilute ash cloud, and the third assumed a “cold” gas phase in the generation of PDC flow.

We included analytical solutions of these kinetic models on a flat topography and with a simplified approach for testing topographic shielding effects on PDC flow runout. Finally, we used a Bayesian belief network related to each inversion model to probabilistically evaluate the uncertainties on the mass required to invade the marker sites, estimating correlation coefficients between the input variables and the calculated mass, and showing that the greatest influence of correlation effects is exerted by sedimentation velocity and the initial concentration of solid particles.

Our probabilistic approach was applied to the case of the Aso caldera, Japan, which is located at distances between ∼115 and ∼145 km from 5 marker sites. We found that a rock avalanche dynamics model would require a ∼3-times greater mass to reach the same marker sites with 50 % probability than the elicited density current dynamics model; a hypothetical scenario of a more dilute ash cloud would require a ∼3-times less mass, and a “cold” density current would require ∼30 % lower mass. Our strategy to consider topography with a simplifying procedure produced the finding that ∼3-times greater mass values are needed to overcome topographic controls than is indicated by PDC runout estimates obtained for a flat surface. We noted that in the rock avalanche model the identification and constraining of the input variables by expert judgment was more difficult than in the density current models.

According to the density current formulation and the emulation of topographic effects, we showed that an axisymmetric PDC with a mass similar to the Aso-4 total PDC overflow outside the caldera would likely reach the sites of the three distal critical infrastructures with probabilities in the range 46 %–100 % depending on the marker site and the uncertainty on the erupted mass. By tentatively testing a mass of $\mathrm{1}/\mathrm{10}$ of the total PDC overflow of Aso-4 outside the caldera, plausibly representative of the single largest PDC flow unit, we obtained a probability range of 0.5 %–40 %. These numbers are based on the PDC mass described in Takarada and Hoshizumi (2020) and only present an illustrative example of application of our methodology to very uncertain estimates.

Our methodology provides a rational basis for assessing the probability of PDC invasion at critical geographic locations when there is major uncertainty about the actual or forecast extent of flow runout from a major magnitude ignimbrite eruption scenario.