Eruption forecasting is one of the key goals for volcanology (Sparks, 2003). Eruption forecasting aims to support decision-makers faced with challenges ranging in timescale from evacuation calls in the short-term to land-use planning in the long-term (Sparks, 2003; Marzocchi and Bebbington, 2012). During periods of volcanic unrest, there is a high degree of uncertainty due to the complexity of volcanic systems and the variability in eruption characteristics, such as eruption vent location and eruptive style and size (Marzocchi et al., 2012). Volcanoes can exhibit a range of geophysical, geochemical and geodetic changes during unrest periods, which can be captured using a range of monitoring techniques as magma ascends towards the surface (Marzocchi et al., 2008; Newhall et al., 2017; Gottsmann et al., 2019). However, the interpretation of these changes can place significant pressure on monitoring volcanologists, especially when critical emergency management decisions depend on them (Marzocchi and Woo, 2007; Marzocchi et al., 2012; Papale, 2017). The high stakes have encouraged the development of quantitative models to support increasingly complex eruption forecasting (Connor et al., 2003; Lindsay et al., 2010; Sandri et al., 2012; Selva et al., 2012; Aspinall and Woo, 2014; Sobradelo et al., 2014; Newhall and Pallister, 2015; Cassisi et al., 2016; Sheldrake et al., 2017). Recent reviews can be found in Poland and Anderson (2020) and Whitehead and Bebbington (2021).

One approach gaining prominence for conducting short-term quantitative eruption forecasting assessments is the use of event trees, first proposed by Newhall and Hoblitt (2002). Event trees provide a graphical representation of a volcanic event with branches connecting nodes representing mutually exclusive and exhaustive behaviours with conditional probabilities based on the previously observed behaviour (Marzocchi et al., 2004, Marzocchi et al., 2008). The idea of an event tree is to provide a generic framework that can be applied to any volcano, where the required input data can be derived from any or all of: knowledge of past eruptions, insights from expert opinions, data from analogous volcanoes, and monitoring data. Event trees can consider the aleatoric and epistemic uncertainty for the volcanic behaviour by applying Bayesian inference (Marzocchi et al., 2004, Marzocchi et al., 2008; Neri et al., 2008; Marzocchi and Bebbington, 2012; Sobradelo et al., 2014; Bartolini et al., 2016). There are two prominent short-term event tree frameworks that incorporate Bayesian inference, Bayesian Event Tree for Eruption Forecasting (BET_EF; Marzocchi et al., 2008) and ST_HASSET (Bartolini et al., 2016). These tools are designed to be developed during quiescence, and thus contrast with other event tree structures for assessing short-term volcanic eruption likelihood, such as the USGS Volcano Disaster Assistance Program (VDAP) event tree approach (Newhall and Pallister, 2015), which is designed to document the volcanologist’s understanding of what is occurring in the midst of a crisis.

The output of an eruption forecasting assessment can inform an evacuation call when combined with cost-benefit analysis (CBA; Marzocchi and Woo, 2007; Woo, 2008; Marzocchi and Woo, 2009; Bebbington and Zitikis, 2016). CBA has been applied in conjunction with BET outputs to support evacuation decision-making based on risk-to-life in the AVF from base-surge phenomena (Sandri et al., 2012) and for other volcanoes (Marzocchi and Woo, 2009; Wild et al., 2019). CBA defines a threshold for an action by weighing the cost of action (C) versus the loss from no action (L) and compares this to the probability of impact (p). When the probability of impact reaches or exceeds the CBA threshold, i.e., p => C/L, the action, e.g., evacuation, is cost-beneficial.

This paper explores the challenges of integrating short-term eruption forecasting with cost-benefit analysis for crisis decision making in areas of distributed volcanism by developing and testing a BET_EF for the AVF. First, we present the model set-up, including the parameter selection methodology. Following this, we conduct a sensitivity analysis of key model parameters. In the absence of observed activity in the AVF, the performance of the developed BET_EF is examined using a published synthetic unrest dataset from multiple eruptive scenarios. Finally, we review the integration of BET_EF with CBA to assess its utility for short-term eruption forecasting and evacuation decision-support for distributed volcanic fields such as the AVF.

Auckland, located in New Zealand’s North Island (Figure 1), is the country’s largest city (population ∼1.6 million; Statistics New Zealand 2018) and produces 37.9% (NZ$₂₀₁₈ 107.8 billion; Statistics New Zealand 2019) of the nation’s GDP. Auckland is also situated upon the basaltic monogenetic Auckland Volcanic Field (AVF). The monogenetic nature of the AVF poses a significant challenge to forecasting the next eruption location (Allen and Smith, 1994; Lindsay et al., 2010). The AVF has an estimated 53 eruptive centres formed in the last 193 kyr, with the last eruption occurring 550–600 years ago (Allen and Smith, 1994; Lindsay et al., 2011; Leonard et al., 2017; Hopkins et al., 2020). Most AVF eruptions have been <0.1 km² in volume; however, the most recent two eruptions, namely Mt Wellington and Rangitoto, have been significantly larger than average (Kereszturi et al., 2013). The magma source supplying the AVF is thought to originate from 70 to 90 km depth based on seismic tomography revealing a zone of anomalously low P-wave velocities (Horspool et al., 2006) as well as geochemical evidence (McGee et al., 2013). The mantle-crust boundary is thought to be at 25–30 km depth (Horspool et al., 2006). This is an important parameter in the context of eruption forecasting, as it is considered that seismicity might not be observed until the magma ascends into the crust (Sherburn et al., 2007).

Given the monogenetic nature of the AVF, there are significant uncertainties around the spatial location of the next eruption and its likely volcanic hazards (Magill et al., 2005; Sherburn et al., 2007; Lindsay et al., 2010; Ashenden et al., 2011; Bebbington and Cronin, 2011; Bebbington, 2013, Bebbington, 2015; Kereszturi et al., 2017). Therefore, a critical issue with evaluating the hazard during the next AVF eruption is identifying the vent location. The AVF is thought to be bounded N-S by the mantle source geometry at depth, and E-W by faults associated with the Dun Mountain Ophiolite Belt (Spörli and Eastwood, 1997; Le Corvec et al., 2013). Runge et al. (2015) defined the spatial extent of the AVF using an ellipsoid around past vents (Figure 1), and recommend that a buffer is applied to this boundary given the most recent eruption was anomalous both in terms of polygenetic behaviour and comparatively large erupted volume. Previous studies have looked at the spatial distribution of the next vent location in the AVF using probabilistic approaches for long-term assessments (Magill et al., 2005; Bebbington and Cronin, 2011; Bebbington, 2013, Bebbington, 2015).

The geologic record indicates two primary eruptive styles in the AVF: magmatic and phreatomagmatic (Allen and Smith, 1994; Lindsay et al., 2010; Kereszturi et al., 2014). Approximately 83% of past eruptions have initiated with an explosive phreatomagmatic phase due to magma interacting with groundwater and/or surface water (Morrissey et al., 2000; Kereszturi et al., 2014; Ang et al., 2020); this eruptive style forms maar craters and tuff rings and produces base-surges (which can extend up to 6 km from the vent in a large eruption; Brand et al., 2014), ballistics and tephra fall. Past AVF magmatic eruptions were primarily Hawaiian or Strombolian in style, forming scoria cones, ballistics and lava flows. At least 60% of past AVF eruptions have transitioned from a phreatomagmatic phase into a magmatic phase. The likelihood of future transitions varies spatially across the field, depending on near surface geology and hydrology (Kereszturi et al., 2014, 2017). Due to the eruption hazards, primarily base-surge phenomena, the Auckland Volcano Field Contingency Plan (Auckland Council, 2015) outlines two evacuation zones: the primary evacuation zone, which extends radially from the vent uncertainty zone; and the secondary evacuation zone, which extends 2 km from the primary. Although the two-zones allow for prioritization, in a future crisis both zones are required to evacuate (Auckland Council, 2015).

The likely expected magma ascent time is poorly constrained in the AVF. Blake et al. (2006) estimate magma ascent rates of 0.03–6 ms⁻¹ for the AVF, based on modeling dyke propagation, xenolith settling speeds and xenocryst host reaction speeds from analogous volcanoes. In the AVF, seismicity associated with magma ascent is unlikely to be detected until it reaches the crust-mantle boundary (Sherburn et al., 2007). If detection occurs at this boundary, i.e., 25–30 km depth, these ascent rates would yield ascent times from this depth of <2 h to 12 days. Brenna et al. (2018) evaluated diffusion gradients in AVF xenocrysts and derived magma ascent rates of 0.01–0.03 ms⁻¹, producing ascent times of 9–35 days. Based on the current understanding of the AVF, in which there is no evidence of magma storage in the crust, magma ascent rates are thought to be too quick to support stalling for significant periods (Lindsay et al., 2010; Mazot et al., 2013; Hopkins et al., 2016, Hopkins et al., 2020). Brenna et al. (2018) also suggested that once magma had reached 1–2 km depth, an eruption is likely in less than 12 h. Estimating ascent times is complex, as a single magma batch may ascend at variable rates (Hopkins et al., 2020).

GNS Science conducts volcanic monitoring in New Zealand as part of the GeoNet monitoring program. Auckland’s current seismic monitoring network consists of 10 short-period seismographs, seven of which are buried to reduce anthropogenic noise, and one broadband station (van Wijk et al., 2021). In 2007 it was thought that recorded seismicity might have to be as shallow as 5 km before being useful to estimate potential vent location (Sherburn et al., 2007). There is no permanent geochemistry or geodesy monitoring equipment within the AVF. At the time of writing, there had been no observed volcanic unrest episodes in the AVF.

A previous BET_EF for the AVF was developed by Lindsay et al. (2010) and run as part of the 2008 Exercise Rūaumoko, a national Civil Defence exercise conducted in Auckland, in which the lead-up phase to an AVF eruption was simulated (Auckland Region CDEM Group, 2008; Horrocks, 2008; Brunsdon and Park, 2009). While the Lindsay et al. (2010) BET_EF performed well as part of the exercise, there have been significant subsequent improvements in the scientific understanding of the AVF, which include new constraints on its geographic extent (Runge et al., 2015), spatial vent likelihood (Bebbington and Cronin, 2011; Bebbington, 2013, Bebbington, 2015), eruptive style and the age of past eruptions (Kereszturi et al., 2013, Kereszturi et al., 2014; Leonard et al., 2017; Ang et al., 2020). These improvements have supported the development of seven additional AVF scenarios (Hayes et al., 2018, Hayes et al., 2019). In addition, the Lindsay et al. (2010) consultation with monitoring scientists to inform parameter selection was minimal and informal. In contrast, this study has taken a more formal approach to elicit expert opinion.

Sandri et al. (2012) applied CBA in conjunction with BET outputs in the AVF to support evacuation decision-making based on risk-to-life safety from base-surge phenomena. That study applied the evacuation-specific CBA approach by Woo (2008), which defines the cost of evacuation as the average socio-economic loss per-capita for the evacuation duration, and the loss from no action as the product of the statistical value of life and the proportion of the population that evacuates that owes its life to the decision. Sandri et al. (2012) produced three thresholds based on the evacuation lasting 2 weeks, 2 months and 6 months, and used a proportion of evacuees that owe their life to the evacuation call of 0.5. This made use of the spatial vent likelihood outputs from Lindsay et al. (2010) to produce evacuation zones for the Exercise Rūaumoko sequence.

Bayesian Event Tree for Eruption Forecasting (BET_EF), originally developed by Marzocchi et al. (2008), is an event tree structure developed to assess short-term probabilities of eruption, thereby allowing it to support crisis management (Figure 2). BET_EF was selected as the most appropriate tool for estimating short-term volcanic hazard for crisis decision-support in the AVF context after evaluating a range of available approaches (Wild et al., 2020).

The BET_EF defined structure comprises five nodes (Figure 2):

Node 1: there is either unrest or no unrest, in the time interval (t ₀ , t ₀ + τ), where t ₀ is the present time, and τ is the time window considered;

At each node, there are three types of inputs:

• Priori beliefs

Priori beliefs are derived from theoretical models or expert judgment, such as the recurrence rate of unrest or eruption at the volcano. At each node, the priori mean belief input is accompanied by an “equivalent number of data” parameter to reflect the weighting, where a value of 1 indicates low confidence and results in a wider distribution to represent the uncertainty. Past data relates to previous observations of activity at the volcano, for example, the number of observed unrest periods.

Each branch is assigned a probability employing a Bayesian approach, represented by a probability density function, formed by the priori information and weighting values, and past data to produce a distribution representing the uncertainty (both aleatoric and epistemic). The more available data, the more uncertainty is reduced, and subsequently, the variability in each node’s calculated output can potentially be reduced. The output from an event tree is an estimated probability for each outcome node.

Key inputs into a BET_EF for short-term eruption forecasting are monitoring data parameters and threshold values for Nodes 1, 2, and 3 (Marzocchi et al., 2008). In previous implementations of the BET_EF, parameters and thresholds have been determined through expert judgement by evaluating available and theoretical streams of monitoring data (e.g., seismic tremor, gas concentrations) and considering what behavior would be considered anomalous at each node. When setting up a BET_EF, the user can either provide one or two thresholds. A single threshold is reached when the parameter is observed, regardless of the magnitude, e.g., “detection of SO₂.” For some parameters, two thresholds are established: one for “background behavior” and another for “anomalous behaviour.” In this case, when observed input data is between the two thresholds, the gradual transition between the two states is modeled using a fuzzy procedure. An example could be “number of long-period earthquakes,” where the lower bound of 3 would indicate potentially anomalous activity for unrest, whereas the upper bound of 8 would be a strong indicator of unrest. Two thresholds can be beneficial as it can be challenging to determine whether the recorded behaviour indicates anomalous diversion from background levels (Gottsmann et al., 2019).

At Node 1, observing a single anomalous parameter sets the probability to 1. Monitoring parameters at Nodes 2 and 3 include a weight. A parameter with an associated weight of 2 means that it is considered to have the same weight as two parameters with a weight value of 1. For Nodes 2 and 3, (Eq. 1) is used to calculate the total degree of anomaly (Z) at a node: Z = ∑ z i ⋅ w i (1) where z is the anomaly and w is the weight for parameter i.

The degree of anomaly is used to calculate the conditional probability distribution for a node. The mean, P, of a Beta distribution is estimated using P = 1 − a . e − b Z (2) where a and b represent the minimum bound (1−a) and the sensitivity to monitoring (b) used to form the Beta distribution. Prior distributions are input for these two parameters, which are updated within the BET_EF framework using data from available observed unrest periods at the volcano in question to form posterior distributions.

BET_EF has been retrospectively applied to the volcanic activity between 2001–2005 at Mt. Etna (Brancato et al., 2011), and between 1981–2009 at Campi Flegrei, Italy (Selva et al., 2012) and to the 1,631 eruption of Mt. Vesuvius (Sandri et al., 2009). The framework has also been used in the emergency management exercise MESIMEX for Mt. Vesuvius, Italy (Marzocchi et al., 2008), Exercise Rūaumoko for the Auckland Volcanic Field (AVF), New Zealand (Lindsay et al., 2010) and the simulation of Cotopaxi, Ecuador (Constantinescu et al., 2015). In all these exercises, estimated probabilities from the BET_EF were consistent with the observed (or simulated) timing and location of eruption and aligned with monitoring scientists’ judgement during the exercise, demonstrating the merits of the framework. An application of BET during the VUELCO Simulation Exercise in Dominica, Lesser Antilles (Constantinescu et al., 2016) demonstrated the importance of users of BET outputs being familiar with probabilistic frameworks in advance of any operational use to ensure the outputs are trusted and applied appropriately.

The BET_EF framework is executed within the PyBETUnrest application (Tonini et al., 2016). This implements both the BET_EF and the BET_UNREST framework (Rouwet et al., 2014; Tonini et al., 2016), the latter of which also considers non-magmatic eruptive phenomena, e.g., hydrothermal eruptions (refer to Tonini et al., 2016 for more detail). PyBETUnrest is an update to the old application form of the BET_EF framework. PyBETUnrest is implemented using the Python programming language, allowing code to be easily modified. Previous testing of the PyBETUnrest application found some discrepancies in how the BET_EF framework is applied between PyBETUnrest and the original application (Wild et al., 2020). In particular, Eq. 2 was not calculated using probability distributions (prior or posterior) but instead used single values of a and b, meaning the procedure is no longer Bayesian, as in the original BET_EF application. In addition, the same values of a and b were used for all three nodes. In the original BET_EF application, prior distributions for a and b are U (0.5, 1) and U (0, 2). For this study, the PyBETUnrest code was modified to reflect the original Bayesian functionality and thus allow uncertainty to be assigned and estimated in the sensitivity analyses. The inputs and outputs from Lindsay et al. (2010) were used to validate the modified PyBETUnrest against the original BET_EF.

To set up the BET_EF for the AVF (hereinafter, BETEF_AVF), both monitoring and non-monitoring components need to be considered for parameter selection. The non-monitoring components were informed by previously published studies for the AVF, following Lindsay et al. (2010).

To inform the monitoring parameters, a scientific consensus approach was applied; this is an approach that has been used to define inputs in past studies for developing volcanic event trees (e.g., Lindsay et al., 2010; Selva et al., 2012; Constantinescu et al., 2016; Tierz et al., 2020) and applied operationally by the USGS’s VDAP team to inform volcanic event-trees in past crises (Newhall and Pallister, 2015). There are alternate methods for collecting expert opinion data, such as Cooke’s classical model (Cooke, 1991), that have been applied in volcanology (e.g., Aspinall, 2006; Neri et al., 2008; Hincks et al., 2014; Bebbington et al., 2018). However, in our case, given the limited number of participants with expertise in each given topic area (e.g., seismic unrest in the AVF), a quantitative elicitation method may have led to very distributed and non-useable data. Furthermore, volcanologists in New Zealand and internationally have operational experience with reaching a consensus for estimating volcanic activity (Potter et al., 2014; Fearnley and Beaven, 2018), while the closely connected nature of the volcanological community suggested a significant risk of a narrow perspective or cognitive bias if a formal elicitation was conducted (Morgan, 2014). As such, the consensus approach within subject matter areas, combined with open discussion, moderated through facilitation, was deemed appropriate for our study.

An expert opinion workshop was thus conducted to define monitoring parameters for the BETEF_AVF. The workshop aimed to establish monitoring parameters and thresholds for the BETEF_AVF for Nodes 1 to 3. The workshop, including preparation, was conducted in four steps:

1. Identify workshop participants

Participants were identified following engagement with the Determining Volcanic Risk in Auckland (DEVORA) research programme management group (Step 1). Specific areas of expertise were targeted amongst the monitoring team from GeoNet and volcanologists from New Zealand universities to ensure representation of participants with expertise in both the AVF and volcanic monitoring. A total of 28 individuals accepted the invitation to participate. Before the workshop, participants attended a virtual presentation on the current understanding of the conceptual model of volcanism in the AVF, what a BET_EF is, and how the data collected from the workshop would be incorporated in BETEF_AVF (Step 2).

At the workshop itself, participants were divided into three groups: seismology, geochemistry and geodesy, based on their respective expertise. Each group was provided with the Lindsay et al. (2010) monitoring parameters and thresholds and asked to record any updates or additional parameters, thresholds or weights for each of Nodes 1–3. Following this, an all-workshop discussion was held to allow all participants to discuss and provide any feedback (Step 3).

A virtual follow-up discussion was conducted to present the results, illustrated through the Exercise Rūaumoko scenario and against the Lindsay et al. (2010) parameters, thresholds and weights (Step 4). This was an opportunity for participants to provide any further feedback. The output of this process was subsequently used to inform the monitoring parameter inputs for the BETEF_AVF. These are collated in Table 1.

This paper presents the set-up and application of a revised BET_EF for the AVF to assess the short-term eruption probability and spatial vent likelihood when applied to eight hypothetical unrest scenarios. While not designed to supplant monitoring and advisory groups, output from such models as BETEF_AVF can provide an initial starting point for discussions between volcanic science advisory group members regarding the current unrest and state of the volcano, and can be prepared in advance of a future crisis (Lindsay et al., 2010; Papale, 2017). Additionally, expert elicitation workshops during the development of a BET_EF are opportunities to bring together monitoring and research volcanologists to discuss and share insights to develop a common conceptual understanding of a volcano, and to highlight areas for future research and monitoring enhancements. For example, many of the parameters captured and subsequently applied in BETEF_AVF are not routinely monitored for in the AVF, given it is a distributed volcanic field with no target for geochemical and geodetic monitoring. However, expert opinion revealed that these parameters could be important indicators of unrest, magmatic unrest and eruption, at nodes 1–3, respectively. These insights are valuable in that they can inform the direction of future upgrades of the monitoring network.

Based on our development and application of the BETEF_AVF we offer the following key general conclusions regarding BET_EF as a tool for volcanic eruption forecasting:

• There is value in reviewing previous BET_EF and similar frameworks as the understanding of a volcano or monitoring capability is improved. This can include evaluation and updating the prior datasets, as well as the input monitoring parameters and thresholds. In our case, we conducted our evaluation and update of the latter via an expert opinion workshop using a consensus approach.

Furthermore, our application of BETEF_AVF and comparison with past work in the AVF by Lindsay et al. (2010) and Sandri et al. (2012) also yielded some interesting insights into potential challenges when applying the tool in areas of distributed volcanism:

• The fixed BET_EF Node 4 50:50 weighting of monitoring: prior inputs (Marzocchi et al., 2008; Sandri et al., 2012) may not be suitable for distributed volcanic fields. The spatial vent location probabilities were higher near the eruption outbreak location. However, when aggregating the spatial eruption probability at each vent across the grid to inform the probability of impact, the evacuation CBA threshold was unrealistically exceeded far from monitoring observations

In summary, this study highlighted several key challenges when applying BET_EF for crisis decision-support to distributed volcanic fields, especially where there have been no past monitored eruptions. We demonstrated that BETEF_AVF performed well based on the synthetic unrest dataset from scenario unrest sequences, yielding increasing probabilities of magmatic unrest and eruption as the unrest sequence progressed, and generating spatial probability maps that highlighted higher probability in the areas of eventual vent formation. Although these were synthetic eruption sequences and thus not real, we believe the results still provide important insights. In particular, we highlight the value of the process of continuous refinement of parameters and thresholds by monitoring and research scientists, and how this contributes to a shared understanding of the conceptual model for a particular volcano.