Aso volcano, Japan: assessing the 100-year probability of a new caldera-forming eruption based on expert judgements with Bayes Net and Importance Sampling uncertainty analysis

Abstract

The Aso-4 explosive eruption on Kyushu, Japan, 89,500 years ago was one of the biggest eruptions in the last one hundred millennia, with a magnitude of approximately M8. Modern society requires the likelihood of natural events with potentially disastrous consequences to be evaluated, even if probabilities of occurrence are diminishingly small. For some situations, it is not satisfactory to assert an event scenario probability is “negligible” or can be “ignored”. Judicial hearings or litigation may require risk levels to be quantified, in which case, statements of scientific confidence could be decisive. Internationally, e.g., for nuclear site safety evaluations, event likelihoods on order of 10^–7/year are often considered for quantitative assessment. At such hazard levels, this might include evaluating the proposition that a particular volcano can deliver a future super-eruption, a supposition that could be attached to Aso volcano. But, simplistically taking the average recurrence interval between past caldera-forming eruptions at a given volcano is an unreliable guide to the likelihood of a future repeat: each past event represented a unique set of tectonic and magmatic conditions within a continually evolving volcanic system. Such processes are not temporally stationary nor statistically uniform. To evaluate the probability of a new M8 event at Aso, within the next 100 years, we performed a comprehensive stochastic probability uncertainty analysis using a model implemented with advanced computational Bayes Net (BN) software. Our eruption process model is informed by multiple strands of evidence from volcanology, petrology, geochemistry and geophysics, together with estimates of epistemic (knowledge) uncertainty, adduced from reviews of published data, modelling and from expert judgement elicitation. Several lines of evidence characterise the likely structure, magmatic composition and eruptive state of the present-day Aso volcano, which has had numerous smaller eruptions since Aso-4. To calculate the probability of another M8 eruption of Aso, we implemented probabilistic ‘Importance Sampling’ in our model. With this approach, we find the chance of an Aso-4 scale eruption (characterised by mean volume 500 km³ DRE and approximate 90% credible interval [210 ‥ 1200] km³ DRE) is less than 1–in–1 billion in the next 100 years (i.e., < 10^–9 probability). Based on current volcanological understanding and evidence, we believe this probability estimate is robust to within an order of magnitude.

Introduction

The Aso-4 explosive eruption on Kyushu, Japan, 89,500 years ago is one of the biggest global eruptions in the last one hundred millennia, with a magnitude of approximately M8 (Takarada and Hoshizumi 2020). Its effects were widespread across the whole of Japan; a similar scale event today would have huge societal impact because more than 5 million people live in Fukuoka Prefecture and more than 14 million in Kyushu (including Okinawa Prefecture).

The location and setting of Aso volcano are shown in Fig. 1, with the present caldera rim outlined and the extent of mapped Aso-4 eruption ignimbrite deposits depicted.

Fig. 1

Location of Aso Volcano, Kyushu, Japan. The outline of the present caldera rim and distribution of Aso-4 ignimbrite deposits are based on GIS data in Hoshizumi et al. (2023)

Increasingly, society demands, and, in many industries, regulatory requirements insist that probability estimates are made concerning very low probability natural hazards, even those as extremal as Aso-4, which might have extensive or severe consequences. These assessments are needed to inform decisions about protecting people and safety–critical facilities. Thus, scientists can be faced with the obligation to evaluate the likelihood of extraordinarily rare hazard scenarios, perhaps verging on the incredible, when the event probability is vanishingly small on normal human timescales, problematic to evaluate and cannot be tested definitively using conventional scientific approaches.

However, given destruction in proximal areas in such a massive eruption will be total, possible mitigation measures, and regulatory safety concerns, can only rationally relate to hazards and risks in distal areas, and therefore entail spatial evaluation of hazard magnitudes, intensities and uncertainties at distances typically exceeding a few tens of km. Thus, the first step, before considering such possible hazards, is to determine the likelihood of occurrence of a future causative event, i.e., a caldera-scale eruption. At the time this analysis was carried out, Japanese law and regulations stipulated that safety–critical power plants could be operated for a maximum of 60 years. Elsewhere, nuclear decommissioning may be considered part of plant total lifetime. Thus, 100 years is adopted here as a generic evaluation period for the Aso caldera eruption probability assessment.

A caldera eruption is one of a class of rare natural threats that can have global impacts; another is the risk of a near-Earth asteroid (NEA) impact. Because available observations do not yet allow the orbit trajectories of NEA objects to be determined well-enough to preclude collision with the Earth, and because there are many such bodies, the possibility exists of a major destructive asteroid impact in the future, albeit this a very low probability event. To quantify the threat, the Caltech Jet Propulsion Laboratory (JPL) maintains its Sentry Impact Risk Page (https://cneos.jpl.nasa.gov/sentry/) which lists known NEA objects with calculations of their impact probabilities in next 100 years (i.e., the same timescale adopted for the present study). In the current JPL Sentry list, impact likelihoods for individual NEA objects are enumerated from 0.1 down to 1.3 × 10^–10 probability (in next 100 years).

Because deterministic or qualitative assessments are impossible to integrate into a comparative risk framework for prioritising protection and mitigation measures, in the civil nuclear industry many safety–critical decisions rely heavily on similar probabilistic hazard and risk assessments (Hill et al. 2009; Hill 2018). In the case of Japan, Nuclear Regulation Authority regulations (NRA, 2019) require the evaluation of volcanic hazards to ensure "…the possibility that a beyond-design-basis volcanic event could affect a nuclear power plant in operation is assessed as not sufficiently small, the site is considered unsuitable for a nuclear power plant.". This said, there is no clear guidance as to what constitutes "sufficiently small" but, internationally, event likelihoods on order of 10^–7/year are often considered for assessment in nuclear safety evaluations (e.g., IAEA (International Atomic Energy Agency), 2012, 2016).

Licensing of a geological repository for nuclear waste in the USA, which had a planned performance lifetime of 1,000,000 years (NRC (US Nuclear Regulatory Commission), 2010), states that only features, events and processes (FEPs) that "are estimated to have less than one chance in 100,000,000 per year of occurring" do not require consideration, i.e., equivalent to no more than 10^–8 per year exceedance probability. Also, licensing of new nuclear power plants in the USA requires consideration of event sequences between 10^–4/year to 5 × 10^–7/year as "beyond design basis events", which are analysed to develop confidence that the proposed facility has additional capacity to safely withstand rare events (NRC (US Nuclear Regulatory Commission), 2020; Nuclear Energy Institute, 2019). In addition, event sequences with probability less than 5 × 10^–7/year probability are evaluated for potentially significant failures of critical components, and for defence–in–depth considerations.

Here, guided by Japan NRA regulations, we use Aso-4, the largest Quaternary eruption from Aso volcano, as representing a maximal magnitude event, and seek to determine the likelihood that an event of such magnitude could occur or be exceeded in the next 100 years. This must be done before the risks associated with the postulated event can be deemed acceptably low or unacceptably high and requiring additional analysis. Motivated, in part, by specific elements of recent legal challenges to regulatory decision-making about extreme events in Japan, including a defined future timescale, we exemplify our hazard assessment approach by addressing the question: what is the likelihood of another caldera forming M8 eruption from Aso volcano in the next 100 years?

Evaluating the likelihood of a future caldera-scale eruption, however, represents a significant challenge in volcanology. Caldera eruptions are rare (i.e., only four past events have occurred at Aso volcano in the past 300,000 years), and these eruptions each represent a unique set of tectonic and magmatic conditions within the volcanic system. Consequently, the recurrence rate of caldera-forming eruptions at any particular volcano cannot be derived accurately from the timing of past events, because each of these events was controlled by a set of evolving geological factors. In essence, this means that such events cannot replicate themselves exactly, which is the implicit assumption underlying a simplistic average recurrence rate.

Expert judgment is used to evaluate scientific hypotheses that are untestable when data are very sparse and frequentist statistical methods, such as magnitude-recurrence relations, cannot reliably inform decision-making. The science of volcanology may, one day, gain enough understanding to characterise the state of the caldera system such that its future activity can be anticipated with some confidence. Nevertheless, decisions regarding the potential for a future Aso4-scale eruption necessitate scientific questions are addressed now, using such information and knowledge as is available. To meet this need, the only rational approach is to conduct a transparent and traceable analysis of current knowledge. Faced with significant epistemic uncertainties on how calderas develop eruptions and with volcanology’s inability to directly measure the current state of the magma system for accurate forecasting, the only viable recourse is to rely on probabilistic judgments provided by volcano experts.

Because such a probability estimate is made to inform decisions about safety issues, any dependable approach should be based on a rationale that enhances the credibility of numerical probability results. In this contribution, we describe an analysis to estimate the probability that a future explosive eruption of Aso volcano − occurring within the next 100 years − might produce a total eruption volume that equals or exceeds the volume erupted in the Aso-4 event, 89.5 kyr BP.

Implicitly, our basis for quantifying the probability an eruption of Aso volcano in the next 100 years could exceed the Aso-4 total erupted volume is to presuppose the maximal probability of an eruption of Aso of any size within 100 years, can be taken as 1 (or very close to 1).The many small-volume mafic eruptions that have occurred historically at Aso volcano (e.g., Kawaguchi et al. 2021) favour this supposition.. Consequently, estimating the probability of an eruption on the scale of Aso-4 is conditional on the likelihood there will be a sufficient volume of eruptible magma present in the volcano to sustain such a large eruption (and, crucially, the uncertainty on the Aso-4 volume). In essence, we are compounding the probability of any eruption in the next 100 years (approximately unity) with the conditional probability its scale will match the Aso-4 case: thus, the latter probability dictates the enumeration of this joint probability.

For calculating the latter, conditional probability, we use a Bayes Net (BN) approach (e.g., Jenson 2001; Cooke et al. 2007; Fenton and Neil 2013), configured as a model designed to evaluate this likelihood. Almost all the volcanological, geological and geophysical aspects of the problem involve substantial scientific uncertainties; these are represented in our BN model to the extent they are informed by observational data, theory or expert judgment (Scourse et al. 2014).

Siting and designing technological facilities that need to be located in regions susceptible to major tectonic events require evaluation of the full range of knowledge and appraisal of plausible alternative models and interpretations, all within a probabilistic framework. This challenge has been clearly demonstrated by the extreme effects of the March 2011 Tohoku earthquake and is nowhere more problematic than in siting facilities with hazard potentials that last for thousands of years, such as geological repositories for radioactive waste. The use of formalised expert elicitation to help derive credible impact scenarios and their likelihoods of occurrence was trialled for the first time in the Japanese geological disposal programme (Scourse et al. 2014). That study looked at the methodology for eliciting expert judgement under uncertainty and explores the broader possibilities of this approach for tectonic hazard forecasting.

As a generality, the likelihood of an extreme or unparalleled scenario, such as a caldera-forming eruption, can be challenging to consider and make judgements about, even probabilistically. Some scientists prefer to restrict their thinking to known precedents and what those might be able to tell us about future hazard probabilities in a quasi-frequentist sense. However, while a particular scenario of concern may have happened in the past, usually there are too few examples to provide a reliable basis for conventional statistical frequentist analysis. Other scenarios, which may be considered volcanologically possible (even if highly improbable), may have occurred but did not leave evidence in the geological record.

Not every Earth scientist is prepared to tackle this challenge of considering the probability of an apparently unprecedented scenario. With respect to the present discussion of Aso volcano, our discussion and expert judgement elicitations have been conditioned explicitly and solely on the likelihood of a caldera-forming eruption at Aso in the next 100 years. Thus, this is one distinct scenario at one specific volcano; other future timescales for Aso (or for other volcanoes), would require further work, elicitations and analyses. The present findings cannot be extrapolated to longer intervals (although interpolation to shorter intervals might be defensible).

Here, we present a re-appraisal of the Aso-4 eruption magnitude, based on recently published or new analyses of field data, with the magnitude estimate expressed with metrics of its uncertainty. Then, we take advantage of fundamental advances in stochastic uncertainty modelling to refine the basis for calculating the probability of an eruption on the scale of Aso-4 in the next 100 years. We argue that this probability can be quantified in support of societal risk assessment and decision making.

Approach

The Aso-4 eruption volume/magnitude is the critical scaling parameter against which we test for the probability that a near-future explosive eruption of Aso could equal (or exceed) the Aso-4 volume. An extensive review of the knowledge base concerning Aso volcano and its history and the publication of an important paper discussing Aso-4 eruption deposit volumes (Takarada and Hoshizumi 2020) are key to a re-determination of Aso-4 magnitude and, critically, the uncertainty associated with that magnitude estimate. Our approach involves the application of a new statistical modelling technique (Rougier et al 2022) for integrating contoured tephra isopach (thickness) data to estimate the total tephra deposit volume and associated uncertainties.

With respect to petrology, recently published investigations on the composition and origins of Aso-4 magma (Ushioda et al. 2020) and post-Aso-4 magma (Kawaguchi et al. 2021) help constrain our understanding of the evolution of Aso volcano and provide a framework for interpreting trends of post-Aso-4 silicic and mafic magmatism. Most importantly, the post Aso-4 magma system appears to represent a new and distinct phase of magmatic activity at Aso volcano. The system appears to form small volumes of predominantly mafic magma, with small amounts of silicic magma, and does not show petrological evidence of generating large volumes of silicic magma arising from extensive mush zones in the deeper crust (e.g., Kaneko et al. 2007; Miyoshi et al. 2011, 2012).

This information, and other published petrological data, was used to develop a conceptual model for the current Aso volcano magma system in an International Group preliminary workshop, held in 2019; the main components of the model are summarised in Fig. 2.

Fig. 2

Simplified diagrammatic representation of Aso volcano magma system for the purposes of this assessment, with three putative reservoirs of uncertain size and magma composition, and alternative schematic pathways to the surface. The same model elements are implemented in various combinations in the Bayes Net probability calculation model, described in text and shown in Figs. 3 - 8

Consideration of available geophysical data, together with the fresh petrological insights, mentioned above, were subsequently synthesized through expert elicitation to inform parameters for the Bayes Net uncertainty modelling of reservoir sizes, compositions, and eruptible volumes.

Benefit is also taken from new advances in performing Importance Sampling analysis (Rubin 1987), here applied within the BN numerical framework. This enables stochastic modelling calculations to be undertaken relating to events and situations with hitherto unquantifiable extremely low probabilities of occurrence; coupled with these stochastic modelling developments, advantage is taken of the increased capacity to process massive spreadsheet datasets afforded by the Data Model and Power Query capabilities in the current version of Excel software (Microsoft 365 Apps for Enterprise version 2104).

In this contribution, the basis for calculating the substantive probability question is informed by model parameter distributions derived from judgements provided by an International group of volcanologists and by a cohort of Japanese specialists (Domestic experts). The International Group undertook the first, initial scoping assessment and recommended that Domestic Group colleagues should be invited to collaborate with any subsequent updating of the assessment, with their knowledge of the most recent information and new work published in the Japanese literature. Although they could disagree with the International Group, individual Domestic Group experts endorsed the basic 2 or 3 reservoirs concept, previously proposed.

These group judgements were obtained using a structured elicitation procedure to pool the judgments of each group (Cooke 1991; Dias Morton Quigley 2018; Hanea Nane Bedford French 2021). Additional file 1 records the names of participating scientists and their affiliations; note, however, that the judgements are not ascribed here to specific individuals (although that information is retained confidentially for purposes of scientific traceability, if needed).

The judgements of the two groups were also combined into a Supergroup to provide a counterpart, collective semi-epistemic basis for contextualising the Aso-4 scale eruption probability findings. The Supergroup findings are deemed semi-epistemic here because, due to COVID-19, it was not possible to hold workshops to ensure all experts had access to the same information and the opportunities to discuss jointly conceptual models, hypotheses and data.

Aso eruption probability Bayes Net design

The fundamental concept of our Aso BN model is a magma-volume accounting framework with two main elements: the total volume of magma that was erupted in Aso-4, 89.5 kyr BP, and the estimated total volume of eruptible magma that may reside currently within the volcano (or could become present within the next 100 years). The lines of evidence informing these aspects are many, varied, and all entail uncertain variables and parameters which characterise the current Aso volcanic system; these are treated as random variables (RV) in our BN framework.

For the purposes of the present assessment, we assume these variables can be regarded as static –albeit very uncertain– over the specified 100-year future timeframe: such a magmatic system evolves significantly only over much longer timescales. This assumption affords an outlook period for which we judge we can have the most confidence when it comes to extrapolating the state of Aso into the near future. As discussed below, we do not believe it is at all likely that the Aso system would change from its current state (mafic volcanism) to major caldera unrest (silicic system) in less than 100 years, but we do consider the possibility of such an extraordinary switch within the elicited judgments.

Clearly, magma composition is a critical element: the Aso-3 and Aso-4 large magnitude explosive eruptions were almost wholly silicic in composition, as were the smaller-magnitude Aso-1 and Aso-2 eruptions. There is no evidence from the geological record that mafic volcanoes are capable of large explosive eruptions of the kind that create large calderas, such as that formed by Aso-4. Nevertheless, eruptions younger than 89.5 kyr BP at Aso volcano are predominantly mafic in composition.

The issue here is not one of a simple binary representation of alternatives: mafic or silicic. Our uncertainty distributions for the two compositions seek to express likelihoods that the magma is either predominantly silicic or predominantly mafic. In the BN calculations, we repeat re-sample the uncertainty distribution of one magma type at random, taking the complement of this sample as the corresponding probability for predominance of the other type; thus, while individual pairs of probability samples sum to 1, the BN resampling process captures, statistically, the relative spreads in the predominance of each magma type, reflecting joint intrinsic uncertainties.

The question might be asked: although the magma system of Aso is predominantly mafic at present, could it evolve to silicic within 100 years? Aside from the problem of whether such a fundamental alteration is physically feasible, we are not aware of any volcano where such a rapid and major change has been documented. This said, silicic and mafic magmas are quite commonly erupted together or in sequence at certain volcanic centres but, for all cases we know about, magmas with such profoundly distinct compositions are stored in different parts of the crust or upper mantle.

While our Aso BN model includes assessments of the key, complementary likelihoods that the presently stored magma is either predominantly mafic or predominantly silicic, for completeness we also incorporate the additional (very remote) possibility that stored mafic magma could evolve to silicic inside the next 100 years. Thus, in our Bayes Net model, we allow for a finite, but minute, probability of some compositional non-stationarity in that one (or more) of the magma reservoirs might, just conceivably, evolve rapidly from mafic to silicic, within 100 years.

The key point is our BN assessment is based on an appraisal of existing knowledge and understanding to create a current hazard “snapshot” for Aso volcano. Unless substantive new knowledge emerges or until something volcanologically significant happens at the volcano, the findings of this assessment would likely hold good for the next 100 years. But, if some important change were to occur, in knowledge or activity, then updating the assessment might be warranted.

The Bayes Net numerical analysis uses stochastic sampling of the multiple variable uncertainty distributions (and other variables and factors which feed into them) to ascertain the probability that the available eruptible magma volume could equal or exceed the estimated Aso-4 erupted volume. This is implemented using the UNINET software package, developed originally by TU Delft in The Netherlands (Cooke et al. 2007; now maintained by LightTwist Software https://lighttwist-software.com/uninet/—accessed 9 June 2022). UNINET is an advanced analytical graphical user interface (GUI) program for high-dimension stochastic uncertainty modelling, multivariate data mining and machine learning, using Bayes Nets, probability vines and dependence trees (Hanea et al. 2006, 2010; Cooke et al. 2007; Kurowicka and Cooke 2011).

The framework of our complete Aso BN comprises three parts:

Sub-net [a]: a subordinate BN with nodes to model lines of evidence contributing to the estimation of the composite total volume of the Aso-4 deposits (shown schematically in Figs. 3 & 4);

Fig. 3

Aso BN Sub-net [a] model: total Aso-4 erupted volume, synthesised from different deposits volumes; calculations implemented in UNINET (Cooke et al. 2007; software maintained by LightTwist Software https://lighttwist-software.com/uninet/; accessed 9 June 2022)—see text for key to node labels. The UNINET package was originally designed for numerical research into stochastic uncertainty modelling and supports only basic functional graphics

Fig. 4

Aso-4 deposits volume synthesis: BN node inputs and results PDFs in summary graphical form. Numbers at the bottom of each panel indicate the node mean value and standard deviation, with volumes expressed as cubic km; detailed PDF statistics can be extracted from UNINET for further calculations or for formal presentation purposes. Note that deposit volumes are expressed as Dense Rock Equivalent (DRE), to allow for direct comparison to magma volumes. The hints of bimodality in some of the top-level nodes arise from the different contributing sub-type PDFs and are not considered meaningful in any physical sense

Sub-net [b]: a subordinate BN for the nodes that represent lines of evidence informing the estimation of the total volume(s) of eruptible magma available in reservoirs in the present-day volcano, which includes the next 100 years (shown schematically in Figs. 6 & 7);

Sub-net [c]: a top-level sub-net of nodes for performing Importance Sampling from two partial statistical cumulative distribution functions (CDFs), derived from the results of computations with Sub-nets [a] and [b]: i.e., for current eruptible reservoir magma volume available in present-day Aso volcano and for total Aso-4 erupted volume, respectively (Fig. 8).

The latter two partial CDFs are defined by sample volumes in the uncertainty range where the upper tail of the total eruptible volume PDF from Sub-net [b] overlaps the lower tail of the PDF for the Aso-4 total volume distribution from Sub-net [a] (for illustration, see Fig. 9). The two CDFs are subjected to joint importance sampling to estimate the probability of an Aso-4-scale eruption in the next 100 years, as discussed below.

With respect to these three sub-nets, the importance sampling Sub-net [c] is, fundamentally, simply a basic calculational node that operates on just two uncertainty distributions, i.e., the syntheses resulting from evaluating Sub-nets [a] and [b] independently. Quantification of uncertainty nodes for Sub-net [a] is discussed later.

Sub-net [b], related to the eruptible volume of magma in the system, presents the greatest challenge when it comes to framing and quantifying uncertainty distributions for the various lines of evidence; to encompass at least some elements of epistemic uncertainty, this has been addressed in the present study by conducting separate structured elicitations with two different groups of experts. This aspect is discussed in the next section.

Aso BN Sub-net [a] – estimation of composite total volume of Aso-4 deposits

In Figs. 3 and 6, nodes representing input data, variables or parameter distributions are shown as plain ellipses; calculational (functional) nodes are shown as ellipses with end blocks. The latter contain equations or dependence conditions for operating on variable samples from data nodes to which they are linked by arrows (‘arcs’ in BN terminology). Top level nodes, representing the key outputs from the BN networks are coloured for identification purposes.

Green nodes, labelled < Unifn > , are simple Uniform probability density functions (PDF) on the interval [0‥1], which are sampled randomly and independently to differentially activate alternative relative weights where these are needed in linked calculational nodes. These Uniform samples are randomised at each sampling / calculation iteration.

BN node sample conditioning using Uniform PDFs – i.e., for nodes in the model to which < Unifnn > nodes are linked directly – can be explained by an example. The < Total_composite_vol > distribution for tephra (Figs. 3 and 4) is calculated based on a sampling regime for its three child nodes, which represent alternative assessments of Aso-4 tephra volumes by different authors or methods. Rather than simply mixing and averaging samples from these three different distributions, at each BN iteration one of the three alternative samples is chosen according to the value returned for that iteration from < Unif2 > , in the range [0‥1]. Within < Tephra_composite_vol > node, the following formula is applied:

$$\mathrm{if}(\mathrm{Unif}2 >= 0.6,\mathrm{ TH}\_\mathrm{tephra}\_\mathrm{vol},\mathrm{ if}(\mathrm{Unif}2 >= 0.5,\mathrm{ Bicub}\_\mathrm{tephra}\_\mathrm{vol},\mathrm{ BF}\_\mathrm{tephra}\_\mathrm{vol}))$$

Depending on the value drawn from < Unif2 > at each BN iteration, this formula selects alternative samples from: < TH_tephra_vol > distribution (Takarada & Hoshizumi 2020); < Bicub_tephra_vol > distribution (bicubic spline fitting), and < BF_tephra_vol > (basis function estimator), with relative proportions: 0.4; 0.1 and 0.5, respectively. These weights were ascribed by the International Group to reflect, in their judgment, comparative evidential strengths of the three alternatives. In this way, a composite uncertainty distribution is created for the variable Aso-4 tephra volume, ensuring uncertainties associated with the different Aso-4 tephra volume estimates are included in the analysis. However, this step is not necessary for nodes comprising just a single uncertain variable.

While the relative weights in this example were assigned initially by the International group, different weights could be chosen and applied; UNINET is ideal for such sensitivity testing. As a general principle for BN stochastic modelling, deploying separate Uniform distributions for conditioning samples from different variable nodes ensures that the latter samples are independent of one another and not prey to systematic or latent dependencies.

In Figs. 3 and 4, node labels comprise combinations of terms and abbreviations, as follows:

• App4_ relates to alternative geophysical data or values, determined as start points for initial scoping tests of the BN and retained for comparative purposes

• BF_ denotes basis function determined deposit volumes (Rougier et al. 2022)

• Bicub_ indicates tephra volume obtained from bicubic spline fitting

• Nonweld is the abbreviation for non-welded pyroclastic density current (PDC) / ignimbrite data

• PDC stands for “pyroclastic density current”

• TH_ denotes data or values taken from Takarada and Hoshizumi (2020)

• Unif is “Uniform”, i.e., a Uniform statistical distribution on interval [0‥1]

• _vol or Vol_ are shorthand for “volume”

Figures 4 and 7 show the same sub-nets, corresponding to Figs. 3 and 6, with selected node uncertainty distributions depicted as histograms for illustration purposes. Summary indicators of node distribution mean value and standard deviation are displayed at the bottom of each histogram panel. Note, however, node uncertainty distributions are not necessarily Gaussian, or even approximately so: some have irregular uni/bimodal shapes, some are heavy-tailed. This can give rise to high standard deviations on the mean when scientific uncertainty about a variable is substantial or skewed, as is often the case in volcanology.

For each Bayes Net node variable, samples are generated from within an interval that is prescribed to be volcanologically plausible for that variable, usually informed by deriving an extended spread that covers, within modest added margins, the full range of judgments expressed in expert elicitation; in this way, no latent or deliberate censoring is needed to avoid non-physical or scientifically egregious sample values.

For input nodes derived for the present study, volume estimates are expressed as Normal distributions with means and standard deviations to match analysis results, where appropriate. For parameters which can take only positive real values (e.g., physical volumes), lognormal PDFs may be invoked simply to avoid negative volume samples being drawn from lower tails of such distributions.

For parameter spreads based on Takarada & Hoshizumi (2020; i.e., nodes labelled: TH_xxxx), these are represented by triangular distributions to accord with the way those authors reported their findings. In the BN model, the triangular distributions are extended below and above their minimum and maximum values allowing the reported endpoint values to be actively sampled as real values. In effect, this is tantamount to ascribing arbitrary – but realistic–uncertainties to their minimum and maximum values for stochastic sampling purposes. To implement this intrinsic range adjustment to their variables, the distributions are here extended beyond Takarada and Hoshizumi’s bounding values by 10% below the relevant minima to 10% greater than their maxima. Some simple tests suggest these extension values are not critical and appear reasonable for the present application and circumstances.

The results of computing the various contributory Aso-4 erupted volumes can be read off the BN chart Fig. 4; the key composite (summed/combination) results are summarised on Table 1. Note that deposit volumes are expressed as cubic km Dense Rock Equivalent (DRE).

Table 1 Aso-4 deposits volumes determined from stochastic re-sampling uncertainty distributions in Bayes Net sub-net model (Fig. 4)

Quantile details of the top-level node Total_composite_vol PDF distribution (from 20,000,000 BN samples) are shown below in Fig. 5. This corresponds to the top node histogram with the same label in Fig. 4.

Fig. 5

BN PDF solution for Aso-4 total composite deposits volume (in cubic km DRE). Mean estimated volume and its standard deviation, and three distribution percentile values are reported; the smallest and largest sample values are indicated on the x-axis. The hint of bimodality arises from the contributions of different types of deposits and their PDFs, this said, the double peak is not considered meaningful in any physical sense

Note, the smallest sample volume is about 212 km³ DRE, substantially smaller than the Takarada and Hoshizumi (2020) minimum estimate of 465 km³ DRE. In addition, the greatest BN sample volume, shown on Fig. 5 (about 1202 km³ DRE), exceeds Takarada and Hoshizumi (2020) maximum bound of 962 km³ DRE. In other words, the BN model generates some bounding lower and upper tail volumes in our sample space that more than span their extrema.

When converted to equivalent eruption magnitudes, Fig. 5 DRE volumes equate to:

$$\text{Mean magnitude M8.1}\pm\text{0.08 magnitude units}$$

$$\mathrm M7.7\,(\mathrm{minimum});\mathrm M8.0\,(5^{\mathrm{th}}\%\mathrm{ile});\mathrm M8.1\,(\mathrm{median});\mathrm M8.2\,(95^{\mathrm{th}}\%\mathrm{ile});\mathrm M8.5\,(\mathrm{maximum}).$$

The LaMeve database (Crosweller et al. 2012) reports the Aso-4 eruption magnitude at M7.7; Takarada and Hoshizumi (2020) assert M8.1 – 8.4. The magnitude uncertainty distribution shown as Aso4_magnitude in Fig. 4, enumerated with the BN model of Figs. 3 and 4, spans these published magnitude values; the BN mean Aso-4 magnitude is the same as the low end of Takarada and Hoshizumi’s range.

This engenders confidence that our Aso-4 volume assessment, as determined here, is an appropriate reference criterion for the purpose of estimating the probability of an eruption of this scale in the next 100 years.

Aso BN Sub-net [b] – eruptible volume estimation

In Fig. 6 (and Fig. 7), many node labels comprise a mix of abbreviations and terms, as follows:

Avail is short for “available”

Comb means “combination”

EruptFrac is shorthand for “eruptible fraction”

Interm is short for “Intermediate”, as in intermediate depth reservoir

Nr_ represents “number of”

Reserv is “reservoir” abbreviated

Unif… is “Uniform”, i.e., separate Uniform statistical distributions on interval [0‥1]

_vol or Vol_ are shorthand for “volume”

Wts signifies “weights”

Fig. 6

Reservoirs / Eruptible volumes synthesis: BN model branches—see text for keys to node titles. As with Figs. 2 and 3, above, each node holds a statistical probability density function (PDF) expressing the variable uncertainty and arrows denote conditional uncertainty dependencies between nodes (i.e., ‘parent–child’ links). During development of a BN model, legacy PDFs from earlier runs can be retained in the UNINET data file and are available for sensitivity testing and for exploring the impacts of weighted alternates

Fig. 7

Aso BN Sub-net [b] model of Fig. 5, with certain node PDFs displayed in summary graphical form. Values at the bottom of each panel indicate the node mean value and its standard deviation; more detailed node PDF statistics can be exported from UNINET for detailed post-processing

Figure 7 shows node PDF distributions for the reservoir volumes and eruptible magma volumes in that sub-network as histograms. For most variable nodes in Fig. 7, the quantile outcomes from expert elicitation are shown in sample histogram form, corresponding to piece-wise CDFs (not shown); UNINET provides summary plots such as these purely for quick visualisation and checking purposes, and auto-scaling can obscure details in some instances.

The eruptible volume sub-net part of the BN model (i.e., Figs. 6 and 7) includes both two- and three-reservoir possibilities for feeding an eruption. Those nodes (viz., < 2ReservCombWts > and < 3ReservCombWts >) are depicted as discrete value distributions.

Also included are uncertainty distributions for: the reservoirs’ present and near-future volumes; the proportions of these volumes that are eruptible, and the likelihoods that individual reservoirs are silicic or mafic. Further model nodes accommodate the likelihoods of scenarios where one or more of the currently mafic reservoirs might transform to silicic magma within 100 years.

Aso BN primary model: remarks

Thus, with this model formulation, the two main trunks of the BN model –one for estimating the Aso-4 total erupted volume and the other quantifying the volume of magma eruptible in the next 100 years– jointly encompass a wide range of plausible scenarios, whilst staying within appropriate boundary volcanological conditions.

In the Aso BN model, all reservoir-related parameters and their uncertainty distributions are presumed independent of each other. Therefore, their parametrisations and, more critically, their uncertainties are treated as uncorrelated. If different variable uncertainties are treated as independent, such random uncertainties tend to cancel out when evaluating known empirical relations or solving equations with stochastic sampling. However, it is possible certain variables in the Aso BN model, and their uncertainties, may be correlated (or anti-correlated) with other variables or variable uncertainties. In principle, and with plentiful data, it would be feasible to derive inter-variable correlations and to introduce these into UNINET to account properly for correlated sampling issues. That issue remains a topic for further research.

But, in the present case, relevant statistical data are very sparse and uncertain, and an alternative option is to obtain and define inter-variable correlations from expert judgement elicitation. However, while this has been done in one or two published studies under very favourable circumstances (e.g., Bamber et al 2019), the process is complex, requires the application of advanced, high dimensionality copula calculus and can be taxing for analysts to implement numerically. At its present stage of development, such an approach, relying as it does on elicited inter-variable correlations which are challenging for experts to enumerate, must be regarded as exploratory, and currently falls outside the scope of the present study. Moreover, if features like this were added to the Aso BN, which are novel and scientifically innovative, potential exists for them to be queried and for the credibility of the modelling results to be challenged.

Aso BN model: practicalities enumerating Aso-4-scale eruption probability

In this analysis, the following fundamental presumption is made: the Aso BN model calculates the conditional probability that the volume of volcanic material produced by Aso volcano in any eruption, within the next 100 years, might equal or exceed that produced during the Aso-4 eruption, 89.5 kyr BP.

Thus, the probability of an eruption of this scale in the next 100 years can be defined as:

$$Pr[Aso\!\!\,-\!\!\,4-scale\,eruption\,in\,100\,yrs] =Pr[an\,eruption\,in\,100\,yrs] *Pr[volume >=Aso\!\!\,-\!\!\,4\,volume|\, expert\,judgements\,\& \,|an\,eruption\,occurs]$$

where the notation | within the probability brackets Pr [..|..] indicates a GIVEN condition or conditions associated with the particular conditional probability.

It is very reasonable to assume that Pr[an eruption in 100 yrs], i.e., the probability of ‘at least one eruption (of any size) in next 100 years’ is asymptotic to 1 (indeed, Aso volcano erupted on 20 October 2021 while this manuscript was in preparation). Therefore, the upper bound probability for an Aso-4-scale eruption in the next 100 years can be enumerated as:

$$Pr[Aso-4-scale\,eruption\,in\,100\,yrs] = 1 *Pr[volume>=Aso-4volume\,|\,expert\,judgements\,\&\,|\,an\,eruption\,occurs]$$

i.e., numerically identical to the righthand term, viz. the Aso-4-scale eruption conditional probability.

Thus, the (conditional) eruption volume exceedance probability, enumerated by the BN model computation, is taken to characterise, with numerical equivalence, the fundamental issue: that is, what is the (conditional) probability of an Aso-4-scale eruption in the next 100 years?

In the discussion that follows, this conditionality is implicit with respect to the probability of a future eruption of Aso of the requisite size.

Importance sampling

To estimate the probability of an Aso-4-scale eruption in the next 100 years, we introduce the technique of Importance Sampling (Rubin 1987) into our BN analysis as an innovative extension of the numerical framework. This technique is not the same as the more familiar Acceptance-Rejection sampling, first proposed by von Neumann (1951). The latter is a procedure for generating samples virtually for a target probability distribution which, usually, has an arbitrary density function envelope. Rejection sampling works by drawing uniformly-distributed samples from a simpler proposal distribution –often a standard form PDF– and then applying a rule to accept just those samples which fall within the target PDF envelope.

Attempting to resolve a very low event probability (e.g., well below 10^–6 probability) by ‘brute force’ Monte Carlo simulation can require upward of a billion samples to quantify the result and its associated uncertainty; put bluntly, manipulating such massive numbers of samples is generally impractical. Instead, Importance Sampling inference (e.g., Rubin 1987) can be used to estimate a posterior density or expectation in a state- or parameter estimation probabilistic model that is too complex to treat analytically as a whole problem, or to solve directly with conventional numerical software. This is where the flexibility of the UNINET software can be leveraged in the present analysis.

In the present study, a separate sub-net is deployed which, selectively, uses just those fractional parts of the generic distributions in the main BN model that are informative and relevant to the probability quantification problem, and can be derived to characterise the relevant statistical uncertainty space. That is, we analyse only those samples from the extreme upper tail of the potential eruption volumes distribution which are large enough to exceed –or are close to exceeding– corresponding low-end stochastic samples from the estimated Aso-4 erupted volume distribution. The importance sampling approach, therefore, is a numerically feasible replacement for the unattainable full-scale problem sampling method.

Figure 8 shows the Aso BN Sub-net [c] which implements the requisite Importance Sampling.

Fig. 8

Aso BN Sub-net [c] model: Importance Sampling model for estimating the probability of an Aso-4-scale eruption in the next 100 years, implemented in UNINET (see text for an explanation of importance sampling). XS_eruptible_vol stands for excess eruptible volume, i.e., eruptible volume sample value is equal to or greater than the corresponding Aso-4 volume PDF sample; node < positive_vol_exceedances > counts the number of times this test is met. Node < Number_Samples > is a constant indicating the overall number of samples taken in the main BN and is the denominator rescaling for the count obtained with < Eruptible_ge_Aso4 > test

Here, importance sampling inference is applied as follows. With the dual branched Aso BN model, outlined in Sect. 3 above –i.e., comprising the Aso-4 volume estimator in Fig. 4, and the available eruptible volume estimator net (Fig. 7)– the upper tail of the PDF for available eruptible volume just overlaps the lower tail of the PDF for node < Total_composite_vol > (i.e., the Aso-4 eruption volume estimate, per Fig. 4); from that pair of BN models, this overlap ‘uncertainty space’ spans the volume range: 213 – 350 km³ DRE.

When the two BN sub-nets are run jointly with 20 × 10⁶ sampling iterations each, just 73 samples are found in the overlap range for < Total_eruptible_vol > (blue distribution: upper tail in Fig. 9), while 23,138 samples are obtained for < Total_composite_vol > (red distribution: lower tail Fig. 9). The latter ≈23 k samples, in the uncertainty space of interest for the Aso-4 volume distribution still represent only about 0.1% of the original 20 million samples for the whole BN!

Fig. 9

Importance sampling: close-up of the uncertainty space where two PDF distributions overlap: in this case, the upper tail of available eruptible magma volume PDF (blue) versus the lower tail of Aso-4 eruption volume PDF (pink); full distributions are not shown to magnify the tiny overlap zone

UNINET’s conditionalized sampling option is used to export those samples which fall inside the overlap range, drawn from the two BN sub-net top-level nodes < Total_eruptible_vol > and < Total_composite_vol > , respectively. These samples are processed in a spreadsheet to define cumulative density functions (CDF) with which to characterise the corresponding quantile statistics and tail properties of the two variable distributions. The resulting CDFs are input back into the importance sampling sub-net of the BN, for the target exceedance probability calculations. This probability is determined by multiplying the importance samples’ exceedance probability test distribution < Eruptible_ge_Aso-4 > jointly with the two sample size ratios; thus, the equation for enumerating the probability of an Aso-4-scale eruption in 100 years is:

$$\begin{array}{l}\Pr[Eruptible\_ge\_Aso\!\!\,-\!\!\,4\,|\,importance\,sampling] * (73/Number\_Samples) * (23138/Number\_Samples)\\\mathrm{where }\,Number\_Samples\,\text{drawn from the main BN is}\,20\times10^6.\end{array}$$

One million iterations of the importance sampling BN sub-net are sufficient to allow the required exceedance probability to be calculated to an appropriate precision.

Without importance sampling, the main Aso BN model – with maximum 20 million samples allowed by the UNINET software – fails to find a single instance where the potential future eruption volume exceeds the smallest volume quantified in the Aso-4 erupted volume distribution. This demonstrates that the probability of an Aso-4-scale eruption in the next 100 years is definitively lower than 5 × 10^–8 but cannot enumerate exactly how much smaller.

Thus, to resolve this probability numerically, recourse to the advanced and innovative importance sampling technique, just outlined, is indispensable. When the BN model top-level volume nodes overlap, samples are jointly pooled in importance sampling mode, as just described, and the mean (conditional) probability of an Aso-4-scale eruption in the next 100 years is enumerated about 5.6 × 10^–10. The associated standard deviation on this mean –derived from all the parameter uncertainties included in the Aso BN– is: ± 1.4 × 10^–9. Numerical uncertainty analysis indicates the corresponding 99% confidence level Aso-4-scale eruption probability in the next 100 years is not greater than 4.2 × 10^–9.

These results are determined numerically from computations with our BN framework, with the latest geological, geophysical and volcanological evidence about Aso contributing to model parameter distributions. For most purposes, and to avoid spurious precision, the probability values reported above should be read rounded to the nearest single significant digit, i.e., 6 × 10^–10 (mean probability) and 4 × 10^–9 (99% confidence).

Thus, to resolve the unique and exceptional challenge of estimating the probability of an Aso-4-scale eruption in the next 100 years, we have adopted an advanced importance sampling technique and implemented it, as an innovation, directly within the Aso BN model, as just outlined.

When a full stochastic sampling run with importance sampling is successfully completed, UNINET has an option to provide an automated report recording the elements of the BN model and the statistical parameters of input and output node distributions; a copy of the final run report for our Aso BN model is reproduced in Additional file 2.

Discussion of findings

In the Aso BN model, and subject to the node parameterisations we have adopted, the minimum sample value for the estimated Aso-4 total erupted volume is 212 km³ DRE. For a new eruption (i.e., one within next 100 years) to match this smallest conceivable volume for the Aso-4 eruption –without involving a deep magma reservoir– would require at least 86% evacuation of the combined maximum eruptible volumes from the shallow magma reservoir (i.e., 95 km³ DRE) and from the intermediate magma reservoir volume (i.e., 150 km³ DRE). Alternatively, 100% evacuation of the shallow reservoir maximum eruptible volume and 78% of the intermediate reservoir maximum eruptible volume would match the modelled minimum Aso-4 total erupted volume. Both these scenarios, however, depend jointly on the shallow and the intermediate magma reservoirs being fully charged with magma that is wholly silicic in composition.

Clearly, if the ‘target’ volume for a future eruption is substantially greater than the sum of these two reservoir maxima, i.e., 95 + 150 = 245 km³ DRE (~ magnitude M7.8), then silicic magma would need to be erupted from the deep reservoir for the prospective eruption to match or exceed the Aso-4 erupted volume. In respect of this scenario, the critical prerequisite is that the shallow and intermediate reservoirs are both currently fully charged with silicic magma to their (barely credible) maximal volumes, as determined by the reservoir volume extrema derived by expert elicitation. If either, or both, of these reservoirs has substantially smaller capacity – or only much lower fractions of reservoir volumes are, for some reason, available for fuelling an eruption (as would be expected from volcanological considerations) – then a very substantial additional volume of extra stored silicic magma would need to be present at depth, currently residing in the lower crustal or upper mantle region. Moreover, to plausibly match the scale of an event approaching the mean volume of the Aso-4 eruption, a large proportion of this present-day deep magma must also be eruptible.

From petrological and geophysical considerations, the existence of the requisite massive volume of silicic magma –stored at lower depths below the volcano– is not deemed credible at this time. In addition, the compositions of the shallow and intermediate magma systems are judged to be predominantly mafic, with only a small likelihood of predominantly silicic magmas existing in these reservoirs.

As reported above, it is here asserted that the (conditional) probability of an Aso-4 scale eruption in the next 100 years is less than 10^–9. Of course, this finding is contingent on the assumptions that underpin our Aso BN model. However, our collective view is the parameter distributions used in the BN model are based on sound, current volcanological reasoning. Therefore, it seems very improbable that an exceedance probability an order of magnitude greater than our result – or even more – could be obtained without asserting strongly different parameter distributions which, at the same time, must remain consistent with available volcanological information and our current understanding of caldera-forming eruption processes.

One important attribute of the present analysis is that the calculations are fully auditable in numerical terms and traceable back to the volcanological properties embedded in our BN model. Thus, with this defined and transparent implementation, we believe the required probability estimate is well founded.

However, this is not to say the treatment of all possible volcanological issues has been totally comprehensive or that every possible nuance has been considered. Two fundamental, and possibly decisive factors, in the context of estimating the probability of a caldera-forming eruption of Aso in the next 100 years are the questions of how many effective reservoirs of magma should be included in our conceptual model and how much eruptible magma is available to power such an eruption.

Thus, to explore these issues in more detail, the present study has been organised to take advantage of the involvement of Japanese volcanologists who have extensive knowledge of Aso volcano.

Updating Aso BN eruptible volumes with expert judgements

The work reported here is based on elicitations of selected key node parameters from two separate expert groups. One group comprised a team of seven domestic experts in Japan (Domestic experts), and the other a group of five international volcanologists (International group) who had been involved in preceding exploratory studies and analyses (names and affiliations are listed in Additional file 1). Both groups had access to similar data sets, although the International group was able to review only incompletely much of the relevant Japanese-language literature. In addition, because of the Covid-19 pandemic, the Domestic experts group elicitation was conducted later than that of the International group elicitation and after the Aso BN model had been composed. As a consequence, however, the Domestic experts were able to consider summary reports prepared by the International group for their elicitation. Due to the exigencies of the pandemic, there were only very limited scientific interactions between the two groups for these two parallel-track elicitations.

In the specific context of estimating the probability of a future Aso-4-scale eruption, this was the first time Japanese colleagues’ judgements were acquired concerning selected eruption BN model parameter distributions, using a structured expert elicitation of the Domestic experts group. The elicited views of the International group represented updates to a previous assessment exercise, modified, where individual scientists deemed it appropriate, by newer scientific knowledge or information about Aso volcano.

It should be noted that, in effect, the two expert groups followed separate approaches to evaluating parameter distributions by expert elicitation. Thus, when the judgements of both groups are melded into joint distributions, these may not have fully captured epistemic uncertainties concerning the state of knowledge about Aso volcano; the combined Supergroup combination distributions are likely to incorporate some different elements of uncertainty, arising from having to conduct separate elicitations, of necessity. Additionally, each group is likely to have considered some subtle variations in available information sources. Thus, this component of overall epistemic uncertainty may not be fully encompassed, and this should be recognised when appraising the Supergroup joint distributions and associated findings.

In many areas of science, continuous data from uncertain variables are often encountered that are not easily- or well-characterized by familiar continuous probability distributions, such as the Gaussian, lognormal, Weibull, etc. Such data may come from measurements, models, assessments or elicitations. However, the true underlying process, or processes, that generates this data are either unknown or do not lend themselves to convenient equations for computing probability density functions (PDF), cumulative distribution functions (CDF) or quantiles (i.e., inverse CDF).

Introduced by Keelin (2016), the Metalogistic distributions are a family of continuous probability distributions that directly address the need to handle such data in hazard and risk assessment work. Compared to traditional distributions, a Metalogistic distribution offers the following advantages: shape flexibility for matching data from almost any source (e.g., empirical data, expert-assessed data, or simulated data); choice of boundedness (unbounded, semi-bounded, or bounded); relative simplicity of the equations; ease of fitting to data with ordinary least squares; closed form quantile function that facilitates simulation.

For the present study, uncertainty distributions –derived from the elicitations of the two Aso volcano expert groups– are presented as fitted Metalogistic curves (see Additional file 4). This allows the judgements of each group to implemented separately for sensitivity runs of the Aso BN and also conflated jointly from all twelve participating experts (i.e., seven from Japan and five internationally), under the communal title ‘Supergroup’.

As discussed below (and in detail, in Additional file 3), these sensitivity tests include two alternative configurations of magma reservoirs under the current Aso volcano, one assuming two effective magma reservoirs, the other envisaging three. From the elicitations, expert opinion on the merits of each configuration, while based on similar scientific evidence and interpretations, was divided. This dichotomy impinges on the weights that can be accorded to these two alternative volcano configurations in the Aso BN probability calculations.

For the elicited Aso BN model variables, some were parameterised by the experts simply as single point values only, while others were characteristics in term of their medians and uncertainties with two additional quantiles. The corresponding distribution fits are reproduced in Additional file 4 and presented there graphically as cumulative distribution functions (CDF) or as probability density functions (PDF), as appropriate.

Evaluating Aso eruption probability with Bayes Net sampling test criterion

Based on our conceptual model of Aso volcano, in its current state and expected state for the next 100 years, we have tested whether random samples –drawn from the distribution characterising the uncertainty associated with the potential volume of explosively eruptible silicic magma– ever equal or exceed random samples drawn from an uncertainty distribution for the actual volume of silicic magma ejected during the 89.5 kyr Aso-4 eruption. The elements of the Aso BN sub-net determining the Aso-4 volume distribution were outlined earlier in this paper; the key BN model metric here is the minimum erupted magma volume estimate for the Aso-4 eruption, i.e., about 213 km³ DRE.

Thus, any sample drawn from the uncertainty distribution for eruptible silicic magma volume must equal 213 km³, at least, to be considered for a positive exceedance outcome; given the Aso-4 volume distribution scales with much greater volumes than the eruptible volume distribution, this means that only samples from the upper tail of the latter have any chance of meeting the volume exceedance comparison test.

This testing is done by extracting 20 million random samples from all the uncertainty nodes in the BN magma reservoir model, keeping count of the number of times the total eruptible volume meets the criterion; for forensic purposes, all corresponding BN node sample values are also recorded: i.e., those that are associated with a positive test outcome.

When attempting to estimate posterior densities or expectations in probabilistic state- or parameter estimation problems that are too complex to treat analytically, as here and as noted earlier, Importance Sampling can be used as a numerically feasible surrogate for full distribution sampling. Full distribution sampling is often precluded by the need for massive computing power to implement –and track– many billions of samples to achieve numerically stable results, e.g., for extremely low probability event occurrence scenarios.

Following the approach adopted for the present analysis, a separate sub-net BN (Fig. 8) performs Importance sampling probability calculations on those elements of the Aso BN model parameter distributions which correspond to the region of uncertainty space that needs to be quantified, i.e., where samples of potential eruptible silicic magma eruption volume may match or exceed counterpart samples from the low end of the Aso-4 eruption volume uncertainty distribution.

The probability results from these tests, determined by Importance sampling with the Aso BN model are summarised now, per expert group. First, the results of the International group analysis:

$$BN\,maximal\,sample\,value-\,eruptible\,silicic\,magma\,volume: 337\,km^3 (\approx M\,7.9)$$

$$Probability\,of\,exceeding\,\mathit{Aso}\!\!\,-\!\!\,4\,volume\,in\,next\,100\,years:2.8\times10^{-9}$$

The maximal eruptible volume and its corresponding exceedance probability reflects parameter updating by the group, e.g., in respect of the weights for the two-reservoir and three-reservoir options and for weights for the different combinations of the two-reservoir scenario (i.e., shallow reservoir alone activated; intermediate alone; or both) and for the seven different combinations of the three model reservoirs (shallow; intermediate; deep – see Fig. 6). The “fatness” of the upper tail in the summed volume distribution is mainly influenced by the weights ascribed by the group to the triple combination of the three reservoirs option.

The corresponding results for the Domestic experts are as follows.

$$BN\,maximal\,sample\,value-eruptible\,silicic\,volume: 136\,km^3 (\approx M\,7.5)$$

$$Probability\,of\,exceeding\,Aso-4\,volume\,in\,next\,100\,years:\,indeterminate$$

The exceedance probability is indeterminate because the two distributions –eruptible volume and actual Aso-4 volume– do not overlap; based on this group’s BN parameters and from the least possible exceedance probability evaluated by importance sampling, it can be asserted that their implied exceedance probability is substantially lower than 10^–12 in the next 100 years.

The results from combining both groups of experts into one Supergroup (joint domestic and international groups collective judgements) are:

$$BN\,maximal\,sample\,value-eruptible\,silicic\,volume: 206\,km^3 (\approx \,M\,7.7)$$

$$Probability\,of\,exceeding\,Aso-4\,volume\,in\,next100\,years:\,indeterminate$$

Again, the exceedance probability is indeterminate for the model parameters provided by the Supergroup combination.

Influential Bayes Net parameters

Magma reservoir existence probabilities and reservoir sizes, and the corresponding maximal available eruptible silicic magma volumes, are clearly very influential when determining the probability result generated by the Aso BN model. This said, these factors are not well-constrained and certain features of their uncertainty distributions are somewhat convoluted.

For instance, in the responses of both the Domestic experts and the International group, there are suggestions of a fundamental bimodality in respect of weightings for the two-reservoir versus three-reservoir alternatives (see Figs. 10A and 10B for two-reservoir model weights).

Some experts in each group ascribe higher relative weights to the two-reservoir scenario, whereas others express higher relative weights for the three-reservoir alternative; thus, the collective responses are evidently bimodal and distinct from one another by group and, as such, do not concede much overlap across these differing judgements.

In ascending order of weights, the five individual International group experts ascribed central weights (coded as WR2 in the elicitation questionnaire) to the two-reservoir option, as follows: 20%; 25%; 30%; 70%; 80% (Fig. 9). In other words, three experts ascribed relatively low weights to the two-reservoir option, while the other two experts accorded it much higher weights.

A similar low/high split is present among the responses of the seven Domestic experts (Fig. 10b), but with a small majority opting for higher weights for the two-reservoir alternative.

Fig. 10

Metalogistic cumulative distribution function (CDF blue line) for experts’ WR2 weights (yellow markers) for the two-reservoir model option: a International Group experts. b. Domestic Group expert. c Supergroup’ combination of Domestic Group and International Group experts (see Figs. 10a & b)

If the two groups are combined into the Supergroup, with twelve experts in total, the balance of low/high dichotomy weights turns out to be six-six (Fig. 10c); more details on these elicitations are presented in Additional files 3 and 4.

It appears that the experts’ attributions of weights to this particular model variable may reveal an underlying issue of judgement coherency. One plausible explanation for this apparent dichotomy is that it reflects, correctly, fundamental limitations to our volcanological understanding. In which case, it would be justified to adopt a bimodal fit to the collective judgements for this underlying epistemic uncertainty. This said, more extensive discussion among the experts could ensure both groups had a common understanding of the Bayes Net conceptual model.

On the other hand, these contrasting distributions might arise from significantly different interpretations by each expert of Aso magmatism and processes, thus giving rise to apparent disjointedness among the views of the participating experts. If so, some of this may be due to the different, and separate, circumstances under which the two groups were elicited, rather than exposing fundamental epistemic uncertainty per se; such incongruity issues can be complex and difficult to resolve. In the present case, they concern mainly parameters for which single-point values were elicited (see, for instance, Issue 1 values and plots in Additional file 4).

Moreover, with sparse data or data of variable quality or with apparently incoherent elicitation judgements, trying to fit some form of familiar distribution (e.g., Normal; lognormal; Beta …) can entail tricky numerical exigencies. For the present hazard analysis, the appeal of the Metalogistic distribution is that it can accommodate multi-point quantile-based uncertainty estimates, with the option of defining different distribution bounds (see also Additional file 3). When processing tens of elicitation items, as here, the Metalogistic approach offers helpful expediency and an element of data treatment consistency.

In some circumstances, fitting a specific type of distribution to data or elicitation items may be rational but, as a general rule, such alternative distributional characterisations often have little effect on bottom-line hazard estimates.

Sensitivity testing: alternative weightings for two-reservoir versus three-reservoir models

This section briefly discusses the sensitivity of Aso eruption probability findings in relation to the apparent bimodality in the weights for the two-reservoir versus three-reservoir model. This analysis explores the implications of the apparent bimodality when estimating the probability of an Aso-4-scale eruption in the next 100 years.

Whereas the substantive reference calculations evaluate a range of epistemic and aleatory uncertainties in developing the eruption likelihood estimate, for sensitivity testing certain epistemic or aleatory uncertainties are systematically removed from the calculations in order to ascertain the significance of those uncertainties for the results. Consequently, probability estimates calculated under sensitivity testing should not, themselves, be interpreted as plausible likelihood values for a future Aso-4-scale eruption -they act only as indicators of the sensitivity of the probability calculation to uncertainties in selected key parameter values.

Nevertheless, because the elicited judgements are such that the Aso eruption probability is numerically indeterminate for the separate cumulative distributions of weights from the Domestic experts or the combination Supergroup, tractable sensitivity testing of the two-reservoir versus three-reservoir issue can be performed only with the International Group’s responses (i.e., Fig. 10b).

The International group Metalogistic probability density function (PDF) for the two-reservoir model weights, as shown in Fig. 11, is the basis for the BN model sensitivity test. The shape of this unimodal curve is obliged to reflect the spread and wide separation of individual weighting judgements, shown in Fig. 11, and is constrained to fall in the bounded range [0 ‥ 100%].

Fig. 11

Metalogistic probability density function (PDF) for basis case International Group WR2 weights (blue line), with sub-regions (within yellow lines) selected for sensitivity testing (see text)

In Fig. 11, two equal width sub-regions are marked by vertical yellow lines for sensitivity testing using UNINET sample-based conditioning: i.e., stochastic BN sampling is restricted to weights in the ranges 20%—30%, and 75%—85%, respectively. These are taken as representing a proxy “low weight” mode (i.e., 25%) and a proxy “high weight” mode (80%), each with the same uncertainty spread. The central lower modal value, 25%, signals that the experts offering WR2 weights in this lower band think the two-reservoir model is approximately only one-third as likely to be appropriate as the three-reservoir model, whereas proponents for the higher WR2 weight mode, centred near 80%, view the two-reservoir model to be about four times more likely to apply than the three-reservoir alternative.

With the BN model, multiple sample values are drawn separately and randomly from either of the two constrained weight regions (and corresponding WR3 complement weights) and the expected (mean) probability of an Aso-4-scale eruption in the next 100 years is re-computed, together with indicative uncertainty distribution statistics.

These are tabulated on Table 2 and compared with the overall result from the International Group’s base case distribution.

Table 2 Results of sensitivity tests on two-reservoir and three-reservoir model options and weightings

For the test when the two-reservoir weight WR2 is assumed to have ‘high’ weighting, within the range 75% to 85%, the expected eruption probability is approximately three times lower than for the (International group) reference basis case. The uncertainty on this estimate (Stdev) is marginally smaller than for the reference case. The 95^th percentile probability is also lower, but to a more limited extent than the reduction in mean probability.

For the two-reservoir ‘low’ weight scenario, the expected eruption probability is roughly 1.2 × greater than the reference case result, and the standard deviation on the mean is marginally greater. The 95^th percentile probability is also greater than that of the reference case.

To provide some additional context, two further sensitivity tests are undertaken (rows 4 & 5 on Table 2). The first of these assumes there is, at most, only a 5% chance the two-reservoir configuration is preferred. The mean probability and the 95^th percentile probability then are greater than those of the basis case or of the two preceding sensitivity tests.

The last example (row 5 on Table 2) considers the situation where the two-reservoir model is assumed to have zero weight and the three-reservoir model is the only option; for this scenario, the mean eruption probability and the corresponding 95^th percentile probability are about 22 × and 30 × greater than the reference case, respectively. Arguably, these might be regarded as approximate upper bound values for the Aso-4-scale eruption probability in the next 100 years (albeit the probability remains very small).

This pair of sensitivity tests provide a quantified perspective on the potential impact on eruption probability calculations of selecting one or the other of the ‘two-schools-of-thought’ dichotomy regarding weightings for the two-reservoir and three-reservoir alternatives, evident in the judgements of both groups and the combined ‘supergroup’, as shown in Figs. 10a-c.

As mentioned earlier, the reason for this dichotomy is unresolved, and may warrant further investigation and discussion among the experts to advance our understanding and assessment of Aso’s volcanology. For this particular hazard assessment, however, the issue is moot: the eruption probability determined with the two-reservoir alternative is so much smaller than that of the three-reservoir model.

Thus, for present purposes, the apparent dichotomy of views about the two-reservoir or three-reservoir alternatives – evident in both groups – can be accommodated in the full BN analysis. This said, the implications of adopting either as the modal preference are investigated by sample-based conditioning in UNINET and the results of this analysis are reported in Additional file 3 Sect. 4.

Further discussion in Additional file 3 Sect. 5 contextualises the issue in terms of maximal reservoir volumes and their BN sampling probabilities and expands on the limitations and uncertainty implications of eliciting single point values for some variables.

Synoptic discussion

As noted above in the section on Evaluating Aso eruption probability, the reference basis overall probability for an Aso-4-scale eruption in the next 100 years is 2.8 × 10^–9, derived from the International Group’s judgements. The theoretical asymptotic minimum joint probability that can be enumerated for such an event—by Importance Sampling with the present BN model—is determined by the BN parameters themselves and conditional on finding at least one single sample that meets the exceedance criterion with importance sampling.

Empirically, the smallest computable probability with this BN model is 1.25 × 10^–22; that minimum probability emerges if only one exceedance sample is found when 20 million samples are taken from each trial distribution and combined with 20 million importance samples drawn from within the uncertainty distribution space where the eruptible volume may be found to equal or exceed the Aso-4 eruption volume.

With each of the different groups’ parameter distributions, the joint probability that all three reservoirs exist, and that their individual maximal silicic magma volumes could be erupted simultaneously, is inferred to be close to, but less than, the numerical minimum probability that is resolvable with our BN model. Hence, for the groups’ different reservoir volumes, the related probabilities of occurrence of erupting the summed maxima from all three reservoirs jointly is in each case < 10^–21, as noted in column 5 on Table 3.3, Sect. 5 Additional file 3.

On that Table, the highest overall joint probability is the product of 0.123 and 1.25 × 10^–21 from the International Group judgements (columns 5 and 6). Thus, given our present state of knowledge, and as far as the extreme scenario of erupting maximal silicic magma volumes from all reservoirs simultaneously is concerned, the corresponding probability of occurrence in the next 100 years can be regarded as zero for all rational intents and purposes.

This extreme scenario provides a perspective on the biggest imaginable – but volcanologically incredible – volume of silicic magma that could be erupted from Aso volcano on a timescale of about 100 years from now.

For the situation where the potential eruptible silicic magma volume might, just possibly, match or exceed the absolute lowest volume estimate for the Aso-4 deposits, only the International Group BN parameters allow joint conditions for a quantifiable probability to be calculated. Counterpart probability calculations are precluded, out of hand, for the BN node judgements of the Domestic experts or for those obtained by combining the judgements of both groups in the Supergroup.

Extreme event probabilities, such as the likelihood of an Aso-4-scale eruption in the next 100 years are based mainly on elicitations of experts’ judgements about the likelihood of process scenarios that could give rise to such an event. These can include circumstances that may have happened in the past (but are, usually, too few in number to provide a reliable basis for frequentist analysis), scenarios which have not been recognized, or that have left no discernible evidence in the geological record but may be considered volcanologically possible (even if highly improbable).

With respect to the Aso problem, the BN discussed here is conditioned on the one specific scenario of an Aso-4-scale eruption and its likelihood of occurrence in the next 100 years. If one wished to consider other eruption scenarios or different future timescales then these would entail modifications to the conceptual BN model, re-framed analysis approaches and additional elicitation questions. For instance, the present findings cannot be simply extrapolated to longer intervals as elements of the present BN model are exclusively couched in terms of the next 100 years (although interpolation to shorter intervals might be defensible).

Over the course of this study, the component node PDFs in our Aso BN were populated with information provided by the various contributing experts and expert groups. Thus, some nodes carry older, prior information while others contain newer interpretations, worked up or published during the project. Ideally, the compositions of the input nodes would have been documented with information sources, contributors, dates and reviews. However, opportunities for in-depth collaborative efforts in this regard were limited by COVID-19 pandemic travel restrictions. Thus, the present applied research explored mainly extreme event probability estimation, breaking new ground with novel methodological advances; it does not represent a fully-fledged engineering analysis and the quality assurance that normally entails, for example. This said, there are, unquestionably, some scientific aspects that could be further developed in greater detail, more effectively capturing the depth of current volcanological knowledge and the breadth of associated uncertainties.

With the current Aso BN model, run in UNINET with importance sampling and relying on the appraised lines of volcanological evidence, the mean probability of an Aso eruption in the next 100 years is enumerated at about 6 × 10^–10. The standard deviation on this mean – derived from all the parameter uncertainties now included in our model – is: ± 1.4 × 10^–9. This numerical uncertainty analysis indicates the corresponding 99% confidence level eruption probability is not greater than 4 × 10^–9 in the next 100 years.

This upper confidence probability is radically smaller than the value 3 × 10^–4 probability in 100 years which would be an inaccurate interpretation of the fact that there has been just one Aso-4-scale eruption within a continuous time period of 300 kyr (i.e., the approximate duration of caldera-forming eruptions at Aso volcano). The latter, frequentist basis for inferring a mean century rate for Aso-4-scale eruptions of Aso, fails to account for the significant evolutionary changes in the petrology of the magma system that have accompanied each of the Aso caldera-forming eruptions, most notably the Aso-4 eruption. Given the magmatic system changed significantly with Aso-4, a frequentist approach to eruption recurrence rate estimation –spanning prior different compositional epochs– is void.

Conclusions

The formal approach, followed here for estimating the conditional probability of an Aso-4-scale eruption, is founded on a rationale which supports the credibility of the quantitative results obtained.

Naturally, such hazard and risk estimates are entirely conditional on the available data, the design of the BN, the robustness and comprehensiveness of component models and any critical assumptions that, inevitably, underpin any analysis. Given our knowledge is, unavoidably, imperfect concerning the stochastic nonlinear processes that drive caldera-forming eruptions, such probabilistic findings cannot claim to be absolute ‘scientific truth’. In fact, this is a truism: no single volcano ever exactly repeats its eruptions: volcanoes are evolutionary complex systems by their very nature and therefore precise prediction per se is not possible.

This said, our understanding and judgement of Aso volcano’s capability to produce a magnitude M8 caldera-forming eruption can be expected only to change minimally over a timespan of years or decades, if at all. In that sense, the present analysis –for the next 100 years and based on our current understanding of Aso’s magmatic system– can be regarded as a stable, robust estimate of that likelihood.

This is not to say the chance of a future Aso-4 eruption in some other future 100-year time window will remain at or below 3 × 10^–9 probability; simply, this is the most dependable probability estimate that can be suggested currently for Aso volcano for the next 100 years, with key uncertainties quantified commensurate with present knowledge and understanding. And if we were to apply the same hazard assessment approach to a different caldera system, it is conceivable we could reach a different conclusion about the 100 years’ probability of an M8 eruption from that other volcano.

More speculatively, had an analysis like this been done, using the same methodology, a century before the Aso-4 eruption occurred, it seems certain an entirely different probability for an ensuing M8 eruption would have been obtained for the subsequent 100 years. It must be supposed that, in the years before Aso-4 89.5 kyr BP, the volcano was manifesting major precursory signs, given its prior history and given the state it must have evolved into just before that caldera-forming eruption. Present-day Aso is an entirely different proposition: it is not presenting with major ominous signs of unrest, collateral dynamic or deformational ‘distress’; and its activity is restricted to relatively minor mafic eruptions from the much smaller central cone complex.

Any elicitation for parameterising a volcanological model entails ‘taking a snapshot’ of current knowledge and understanding. Subsequently, new information likely will continue to emerge. For example, in a paper published after the present study was completed and the present contribution had been submitted (we are grateful to a reviewer for drawing this to our attention), a new petrological interpretation of the evolution of Aso volcano magmas has been proposed by Miyagi et al. (2023). As with any new information, its potential significance for the existing assessment can be difficult to assess on an immediate basis. Some fresh information might provide reasons to update existing parameter ranges, such as a substantially revised volume estimate for the Aso-4 eruption. Sensitivity analyses can be conducted to determine if some new information is potentially significant, given the specific context, and, if so, would support the need for an assessment update. In contrast, other new information might present alternative conceptual approaches which are more difficult to assess for potential significance.

The current elicitation developed a ‘multiple reservoirs’ model that integrated existing petrological, volcanological, and geophysical information, and used expert judgment to assess, inter alia, the likelihoods of different model states and different magma pathway combinations. Determining, in isolation, the potential significance of further information (e.g., the detailed petrological analyses reported by Miyagi et al. 2023) on this conceptual model is challenging: without reconvening the expert groups, it cannot be known how they would now weigh this ‘new’ evidence within – or against – the totality of other petrological, volcanological and geophysical information. This would be a non-trivial exercise and still might have only marginal impact on the substantive issue: the probability of an M8 + eruption of Aso volcano in the next 100 years.

The main conclusion to be drawn from this Aso BN model analysis is that, notwithstanding inevitable epistemic and aleatory uncertainties, the current state of the volcano’s magmatic system is believed to be such that, fundamentally, the likelihood of an Aso-4-scale eruption in the next few decades is infinitesimal; on this basis, consideration of detailed societal risk-exposure evaluations and mitigation measures does not appear warranted.