Cost-benefit analysis for evacuation decision-support: challenges and possible solutions for applications in areas of distributed volcanism

Abstract

During a volcanic crisis, evacuation is the most effective mitigation measure to preserve life. However, the decision to call an evacuation is typically complex and challenging, in part due to uncertainties related to the behaviour of the volcano. Cost-benefit analysis (CBA) can support decision-makers: this approach compares the cost of evacuating versus the expected loss from not evacuating, expressed as a ‘break-even’ probability of fatality. Here we combine CBA with a Bayesian Event Tree for Short-term Volcanic Hazard (BET_VHst) to create an evacuation decision-support tool to identify locations that are cost-beneficial to evacuate in the event of volcanic unrest within a distributed volcanic field. We test this approach with the monogenetic Auckland Volcanic Field (AVF), situated beneath the city of Auckland, New Zealand. We develop a BET_VHst for the AVF, extending a recently revised Bayesian Event Tree for Eruption Forecasting (BET_EF) to consider the eruptive style, phenomena produced, and the impact exceedance probability as a function of distance. The output of the BET_VHst is a probability of volcanic hazard impact at a given location. Furthermore, we propose amending the weight of the monitoring component within the BET_VHst framework to a transitional parameter, addressing limitations identified in a previous study. We examine how three possible transitional monitoring component weights affect the spatial vent likelihood and subsequent BET_VHst outputs, compared to the current default weight. For the CBA, we investigate four thresholds, based on two evacuation durations and two different estimates for the value of life that determine the cost of not evacuating. The combinations of CBA and BET_VHst are tested using a synthetic unrest dataset to define an evacuation area for each day. While suitable evacuation areas were identified, there are further considerations required before such an approach can be applied operationally to support crisis management.

Introduction

Evacuation is a common strategy to mitigate the loss of life of those acutely at risk from a volcanic eruption (Tobin and Whiteford 2002; Marzocchi and Woo 2007; Moriarty et al. 2007; Wilson et al. 2012). Nevertheless, the decision to evacuate is not always straightforward, and evacuation is logistically complex and costly (Woo 2008). During a volcanic crisis, decision-makers often rely on experience and qualitative interpretation of monitoring observations (Paton et al. 2008; Marzocchi et al. 2012; Newhall and Pallister 2015). However, this is fraught with uncertainty due to the complexity of volcanic systems and the variability in unrest and eruption characteristics. This uncertainty, coupled with often limited data, places substantial pressure on monitoring volcanologists and emergency management officials when making high-stakes decisions during an event (Marzocchi and Woo 2007). In terms of the preservation of life, examples of successful and justifiable evacuations include those made during the 1973 Heimaey (Williams and Moore 1983), 1991 Pinatubo (Newhall 1996; Newhall et al. 1998), and 2010 Merapi (Mei et al. 2013) eruptions. There have also been examples of unrest episodes that did not result in an eruption, such as in 1976 at Guadeloupe and during the 1980s at Campi Flegrei (Woo 2008). Nevertheless, lacking the advantages of hindsight, evacuation may have been justified even if there is no subsequent eruption.

Cost-benefit analysis (CBA) is a standard decision-making method applied across a variety of sectors, including industry and government (for example in engineering projects), and has been proposed for evacuation decision-support during a volcanic crisis (Marzocchi and Woo 2007, 2009; Woo 2008; Sandri et al. 2012; Bebbington and Zitikis 2016). CBA defines a threshold for action by weighing the cost of action (C) versus the loss from no action (L). This ratio is compared to the probability of an adverse state (p). When this probability exceeds the CBA threshold, i.e., p = > C/L, the action, in this case, evacuation, is cost beneficial, providing a binary decision outcome (Katz and Murphy 1997; Woo 2008).

While CBA provides a simplistic framework, populating the different parameters for volcanic crisis management is convoluted, even with regards to defining the ‘cost’ and ‘loss’ of human life (Sobradelo et al. 2015; Woo 2015). More particularly, establishing the probability of an adverse state is challenging across an unrest episode, especially in complex volcanic settings with varying vent locations, eruptive styles and expected phenomena. This has led to 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). One quantitative approach for eruption forecasting gaining prominence is the use of event trees, which provide a graphical representation of volcanic events with nodes representing changes in the volcanic behavioural state connected by branches (Newhall and Hoblitt 2002; Marzocchi et al. 2004, 2008; Newhall and Pallister 2015). Each node consists of conditional probabilities, given the previous state is observed, with the inputs derived from a combination of knowledge of past eruptions, insights from expert opinions, data from analogous volcanoes, and monitoring observations.

Quantitative eruption forecasting models and CBA have many advantages. They can be developed and reviewed in advance of a volcanic event, thereby removing pressure from the volcanic monitoring team during a crisis (Lindsay et al. 2010; Papale 2017). They can draw from the input of a range of experts (Aspinall 2006; Neri et al. 2008; Newhall and Pallister 2015; Constantinescu et al. 2016) and can be documented and auditable, so any subsequent decision is justifiable with evidence (Woo 2008, 2015). The use of CBA has been theoretically assessed in conjunction with event trees to support evacuation decision-making in volcanic crises in a number of past studies (Marzocchi and Woo 2009; Marrero et al. 2012; Sandri et al. 2012; Sobradelo et al. 2015; Wild et al. 2019b).

This paper explores the use of integrating CBA with an eruption forecasting event tree to provide evacuation decision-support in distributed volcanic environments. This builds on the previous study by Wild et al. (2022), who used Bayesian Event Tree for Eruption Forecasting (BET_EF; Marzocchi et al. 2008) and CBA for the Auckland Volcanic Field (AVF), New Zealand to define evacuation zones in an series of evolving hypothetical unrest episodes. That work revealed two fundamental issues in the setup as it applies to distributed volcanism (Wild et al. 2022):

• The BET_EF spatial vent likelihood methodology may not be appropriate for informing quantitative approaches for defining evacuation areas.

• There is a need to consider the volcanic hazard intensity and extent and the associated risk-to-life safety, which BET_EF does not assess.

To address these issues, we first review the input parameters and the outputs from the eruption forecasting framework that are required to apply CBA for evacuation decision-support. Then we review the BET_EF spatial vent likelihood approach for distributed volcanic environments based on the findings of Wild et al. (2022). Following this, we extend the previously developed BET_EF to consider the spatial probability of vent location, style and produced phenomena to assess the associated risk-to-life using the Bayesian Event Tree for Short-Term Volcanic Hazard (BET_VHst; Selva et al. 2014). Finally, we review the integration of BET_VHst with CBA to assess its utility for short-term eruption forecasting and evacuation decision-support for the AVF (and by extension, other areas of distributed volcanism). We conclude by providing a discussion of challenges of applying this approach in areas of distribution and make recommendations for possible ways to overcome these.

Cost-benefit analysis background

The decision to initiate mitigative measures, such as whether, when and where to evacuate, is often difficult for decision-makers due to the systematic complexity of the impacts of such a decision, especially in the midst of a volcanic crisis (Woo 2008). This has led to the proposed use of CBA to quantify the cost of mitigation measures against the potential loss from no action, and to provide decision-makers with an exceedance threshold to determine a yes or no outcome (Vrijling et al. 1995; Katz and Murphy 1997; Woo 2008). CBA when applied to volcanic crisis decision-support aims to reduce the pressure on scientists and emergency management officials by having a predefined threshold of exceedance for an evacuation call.

Marzocchi and Woo (2007) expressed the described CBA approach as a ratio representing the cost to mitigate (C) against the loss given no action (L) in the event of a hazardous event, compared to the probability of impact (${P}_{Evac}^{\left(X\right)}$) exceeding a predefined threshold, which in this application is the probability of volcanic hazard impacting a given location (X). CBA indicates the action becomes cost-beneficial when:

$${P}_{Evac}^{\left(X\right)}>\frac{C}{L}$$

(1)

Woo (2008) showed that, in the case solely considering risk to life, the C/L ratio specifically for evacuation decision-support is:

$$\frac{C}{L}=\frac{N.R}{N.V.E}=\frac{R}{V.E}$$

(2)

where N is the number of people at risk, R is the socio-economic loss per capita for the duration of the evacuation, V is the economic value of life and E is the proportion of those that owe their life to the evacuation call given the hazardous event occurs. The value of life (V) must be quantified, which requires economic and ethical considerations. The latter is particularly important as the value selected will affect subsequent evacuation decisions. Previously documented approaches for evaluating V include the financial sums to relocate people out of hazardous zones (Marzocchi and Woo 2009) and the loss of GDP per capita over the average working life (Jones-Lee 1994; Woo 2015).

While CBA has its benefits, there are some limitations; namely, its rigid parameters do not allow input of intangible losses and costs, such as psychosocial well-being and stress resulting from an evacuation. Furthermore, the cost of an evacuation likely extends beyond the immediate loss of GDP per capita. There could be wider economic costs from an evacuation, such as loss of industry due to the displacement of staff, or an evacuation zone isolating facilities, which could extend the economic losses to a nation’s GDP, and even result in global supply chain issues, similar to those observed due to loss of production during the COVID-19 lockdowns. As a result, such costs and losses need to be factored into a CBA. Additional costs are also incurred during an evacuation, for example, for accommodation, transport, and financial assistance for the displaced population (Woo 2015). As such, these values need to be factored into R.

CBA has been applied in conjunction with eruption forecasting frameworks to support evacuation decision-making in a number of studies (Marzocchi and Woo 2009; Sandri et al. 2012; Sobradelo et al. 2015; Wild et al. 2019b, 2022). All these examples used event-tree forecasting approaches, either retrospectively or for future planning, to demonstrate the methodology for evacuation decision-making.

Bayesian event tree for forecasting eruptions and volcanic hazards

Overview and past applications

The Bayesian Event Tree for Eruption Forecasting (BET_EF; Marzocchi et al. 2008) was the first in the suite of BET’s (Fig. 1). At each node, BET_EF uses the combination of a non-monitoring component, formed using a Priori belief and past data, and a monitoring component, based on exceedance(s) of monitoring parameter thresholds.

Fig. 1

Schematic diagram of BET_VHst framework. The dashed outline presents the nodes that comprise the BET_EF framework. Any branch that terminates with “clone” is identical to that of the top branch for that node

At each node (k), the conditional probability distribution (${\theta }_{k})$ is a linear combination of the two-probability distribution functions for the monitoring $({\theta }_{k}^{\{\text{{\rm M}}\}})$ and non-monitoring $({\theta }_{k}^{\{\stackrel{-}{\text{{\rm M}}}\}})$ components (Marzocchi et al. 2008):

$${\theta }_{k}={\gamma }_{k}\left[{\theta }_{k}^{\left\{\text{{\rm M}}\right\}}\right]+\left(1-{\gamma }_{k}\right)\left[{\theta }_{k}^{\left\{\stackrel{-}{\text{{\rm M}}}\right\}}\right]$$

(3)

where ${\gamma }_{k}$ is a value in the interval [0, 1] that represents the weighting of the monitoring component. At all nodes, ${\gamma }_{k}$ is determined based on the degree of unrest $\left(\eta \right)$. At Node 1, observing a single anomalous monitoring parameter sets $\eta =1$, thereby the model deems the volcano to be in unrest and the probability of unrest to 1. At Nodes 2 and 3, ${\gamma }_{k}=\eta$, which means when the volcano is in unrest (i.e., $\eta =1$), the output probabilities at each of these nodes are fully informed by the monitoring component.

BET_EF was extended to include the assessment of volcanic hazards and their spatial extent in the Bayesian Event Tree for Volcanic Hazards (BET_VH; Marzocchi et al. 2010). BET_VH adds three nodes for assessing the occurrence of hazardous phenomena (at Node 6), their spatial extent (at Node 7), and the exceedance of the hazard threshold (at Node 8). However, BET_VH merges Nodes 1 to 3 into a single node representing the probability of eruption (P(Erup) and does not allow for the input of monitoring parameters. As a result, BET_VH is a long-term hazard assessment model with limited use within crisis management. The BET_VH framework has been implemented in the application PyBETVH (Tonini et al. 2015) in which Nodes 7 and 8 are merged into one node to output the probability of hazard threshold exceedance at a given location. The BET_VH framework has been applied in a number of past studies, including for tephra fall hazard at Campi Flegrei, Italy (Selva et al. 2010), at Okataina Volcanic Centre (Thompson et al. 2015) and at Mt. Taranaki (Wild et al. 2019b), both in New Zealand, for multi-hazards for El Misti, Peru (Sandri et al. 2014) and for base-surge phenomena in the AVF (Sandri et al. 2012).

BET_EF and BET_VH have been amalgamated to form the Bayesian Event Tree for Short-Term Volcanic Hazard (BET_VHst; Fig. 1) to allow for short-term probabilistic volcanic hazard assessment (Selva et al. 2014). BET_VHst takes Nodes 1 to 5 from the BET_EF with Nodes 6 and 7/8 from BET_VH, which allows monitoring data to be input, and can thus support crisis management. BET_VHst has been applied retrospectively to the MESIMEX simulation to assess tephra dispersal hazard at Vesuvius (Selva et al. 2014). In addition, BET_VHst has been proposed as an automated near-real time tool for assessing tephra dispersion pre-eruption at Vesuvius, Campi Flegrei and Etna (Selva et al. 2015a, b).

Consideration for spatial vent likelihood at node 4

Node 4 is used to assess the spatial vent likelihood for an eruption. In BET, two styles of volcanic geometries are considered: a ‘cone’ with a central vent and four flanks; and a ‘field’ defined by a grid used for assessing volcanic fields or calderas. As with Nodes 1 to 3, Node 4 is calculated using non-monitoring components as prior models and past vent locations, and monitoring components. The input for monitoring components is a file containing conditional probabilities of an eruption at each location. How the monitoring component is informed is volcanic context-specific and various approaches can be applied, such as the localisation of any of the monitored anomalous activities (Constantinescu et al. 2016; Tonini et al. 2016), or, where monitoring capability permits, a more specific approach such as weighting the likelihood by the inverse of the hypocentral earthquake depth (Lindsay et al. 2010; Wild et al. 2019b). Unlike Nodes 2 and 3, at Node 4, even when $\eta =1$, both the monitoring and non-monitoring components are used to assess the outputs (Marzocchi et al. 2008):

$${\gamma }_{k}=Min\left(0.5, \eta \right)$$

(4)

The justification for this is to not rely solely on the interpretation of the unrest signal (derived from the monitoring component) and to mitigate false localisation due to monitoring anomalies, such as the migration of seismicity (Marzocchi et al. 2008). Wild et al. (2022) recently identified some issues with the weighting (4) at Node 4 during the application of BET_EF in the AVF (a field geometry), when combined with CBA to inform evacuation zones. While the outputs of absolute eruption probabilities at a location were highest where there was increased unrest, when aggregating these probabilities within a 5 km radius of a location to determine the evacuation zone, informed by the Auckland Volcanic Field Contingency Plan (Auckland Council 2015), areas far from the focus of unrest were shown as cost-beneficial to evacuate. Furthermore, the evacuation area grew as the P(Erup) increased, eventually encompassing the entire field. This contradicts the expected behaviour, in which we would envision the area of the evacuation zone to decrease when unrest parameters converge, and thus confidence in the vent location increases.

The reason for this contradictory behaviour is the maximum weighting of 0.5 at Node 4 for the monitoring component (Eq. 4). This is illustrated in Wild et al. (2022), who showed for distributed volcanism the effect that the high value (0.5) of the non-monitoring component has on output spatial vent probability. They found that with a 500 × 500 m grid, when using a uniform prior distribution, the non-monitoring component of spatial probability of occurrence in a cell is ~ 1.5 × 10^− 4. Within the 5 km radius from a possible vent location, the sum of all vents’ non-monitoring components is ~ 4.7 × 10^− 2. Application of the code thus only requires a P(Erup) (i.e., the product of Nodes 1–3) of ~ 0.31 to include every cell in the field as being cost-beneficial to evacuate, regardless of the location of monitoring observations (using the 6-month CBA threshold of 0.0143 from Sandri et al. 2012). Wild et al. (2022) recommended that the weighting of monitoring and prior components at Node 4 should be revised for distributed volcanism. Here we add that any weighting change also needs to prevent even further inflation of spatial likelihoods when there are minimal unrest observations (e.g., seismicity) to inform the monitoring component for vent localisation. For example, when there are only a few seismic events at depth, the vent likelihood based on the monitoring component is shown as more confident and more spatially confined than when there is lots of seismicity, even at a shallower depth, which is also counterintuitive. A transitional weighting of the monitoring component as unrest increases should prevent this issue; in other words, as unrest increases, the output from Node 4 would transition from being informed predominantly by the prior to being predominantly driven by monitoring observations. We test this approach here.

Case study: evacuation decision from an AVF eruption

Here, we combine CBA with an eruption forecasting event tree to evaluate when and where to evacuate in the lead up to an eruption in a distributed volcanic environment. We use the case-study of the Auckland Volcanic Field (AVF), a monogenetic volcanic field (Fig. 2) situated beneath Auckland, a city in Aotearoa New Zealand of 1.6 million people (Statistics New Zealand 2018) that is responsible for 37.9% (NZ$₂₀₁₈107.8 billion; Statistics New Zealand 2019) of the nation’s GDP. Sandri et al. (2012) developed an initial BET_VH for base-surge phenomena for the AVF, building on the 2008 Exercise Rūaumoko scenario (Auckland Region CDEM Group, 2008; Horrocks 2008; Brunsdon and Park, 2009) and extending the previous BET_EF from Lindsay et al. (2010). They also combined the outputs with CBA to present a preliminary evacuation decision-support tool for the AVF.

Fig. 2

Auckland Volcanic Field vent locations with “tight” and 5 km buffer boundaries from Runge et al. (2015) to illustrate the current understanding of the AVF extent. The 5 km buffer is a conservative estimate of the spatial extent of the AVF

In addition to the problem discussed above regarding the ratio between monitoring and non-monitoring components at Node 4 (i.e., when assessing the spatial vent likelihood), Wild et al. (2022) also identified another limitation when integrating CBA with BET_EF to define evacuation zones in the AVF. Namely, applying the 5 km radial area from the AVF contingency plan (Auckland Council 2015) assumes an equal risk-to-life within that zone, i.e., the same risk for someone at the vent as someone at a 5 km distance, which does not reflect reality for the spatial extent of the hazard nor the likelihood of casualty. Indeed, an understanding of the spatial variability in eruptive style and consequent hazardous phenomena in the AVF is required to fully assess risk-to-life-safety to inform ${P}_{Evac}^{\left(X\right)}$ within a CBA. As such, it is desirable to extend the BET_EF to a BET_VHst for AVF.

In the following subsections, we outline the methodology applied to combine CBA with a BET_VHst to inform an evacuation decision-support approach for the AVF.

Cost-benefit analysis thresholds

To apply the CBA approach presented in Eq. 2, there is a need to define R, V, and E. As part of defining R, we consider, pre-eruption, two possible evacuation durations to review the influence these have on the CBA threshold and the subsequent formed evacuation zone. In reality, however, the actual evacuation duration will vary depending on the volcanic and societal context, likely ranging from weeks to over a year, as observed in past global events (Woo 2008). First, we consider a duration of 3 months. This follows Woo (2008), who considered a 3-month-long unnecessary evacuation an appropriate average disruption duration to assess the cost, C, for an unrequired evacuation (assessed on a probability-weighted basis considering the variability in duration and rates of compliance). As an alternative, we use an evacuation duration of 1 month, which might be indicative of an evacuation duration in the event of a shallow intrusion with no eruption in the AVF. It is presently thought that the lead up to an eruption in the AVF is likely to be ~ 9 - ~35 days, based on detection of seismicity occurring once magma has crossed the mantle-crust boundary at 25–30 km depth (Horspool et al. 2006), with ascent rates of 0.01–0.03 m/s (Brenna et al. 2018). The conceptual model for volcanism in the AVF includes ascent rates from the crust-mantle boundary to the surface that are considered too fast for magma stalling to occur in the upper crust (Mazot et al. 2013; Hopkins et al. 2016, 2020). This suggests that if an eruption were to eventuate following detection of unrest, it would most likely occur within a month of any associated evacuation occurring. Sandri et al. (2012) presented a duration of 6 months for the evacuation period, which was suggested by Auckland emergency management officials in the past, based on the length of the 1973 Heimaey eruption and associated evacuation. However, the Heimaey evacuation length was based on post-onset eruption duration. In Auckland, once the vent is known, an exclusion zone can be established, allowing those who initially evacuated from areas to return if they do not end up being within the area impacted by the new volcano. However, it is simply not known how long an AVF eruption would go on for, and who of those affected by the eruption will be able to return. Other issues around the dynamic complexity of the evacuation timing and duration are examined by Bebbington and Zitikis (2016); here we consider only the initial decision of whether to evacuate, and what the evacuation zone should be.

The CBA is deciding between an ultimately unwarranted (no eruption) evacuation and a failure to evacuate preceding an eruption occurring. The 1 and 3 month evacuation durations are combined with the GDP per capita in Auckland of NZ$72,000 (Statistics New Zealand 2021) to assess the per capita loss due to an evacuation. The GDP per capita represents an average across the total population. To account for additional costs of an evacuation, we consider the expense of housing and/or financial support for those that do not have alternate accommodation, such as a secondary home or the ability to stay with friends or family. Based on a recent survey of Aucklanders, 35% of respondents would evacuate to a friend’s or relative’s house (Thakur et al. 2022), and it is assumed that the remaining proportion would require accommodation support. Here, to demonstrate the approach, we select a value of NZ$100 per day per person for the duration of the evacuation for the evacuees modelled as requiring financial assistance and/or accommodation. This accommodation support cost is included here for the proportion of evacuees (65%) requiring accommodation support, in order to demonstrate how such additional costs can be considered as part of the overall cost of an evacuation (C) within CBA; however, the actual value could differ (e.g., Whitehead 2003). This leads to R values of NZ$23,931 (72,000 * 3/12) + (100 * 0.65 * 365 * 3/12) and NZ$7,977 per person (72,000 * 1/12) + (100 * 0.65 * 365 * 1/12), for an evacuation duration of 3 months and 1 month, respectively.

For this study, we consider two values of V to demonstrate the influence this value has on evacuation decisions. The first is a value of NZ$4.46 million, based on the value-of-statistical life (VOSL) from the Ministry of Transport (2021), as applied in the New Zealand Treasury CBA tool (New Zealand Treasury 2021); this is the metric applied by Sandri et al. (2012). Additionally, we consider the lower value of 25-times the GDP per capita proposed by Woo (2008), which is NZ$1.8 million. Of note, V is not dependent on the evacuation duration. E is estimated as 0.65, based on a survey finding that 65% of people would wait for an official evacuation call to be made before evacuating, whereas 35% would self-evacuate in advance (Thakur et al. 2022). This results in four CBA thresholds, presented in Table 1.

Table 1 AVF CBA thresholds (estimated probability of volcanic impact) calculated in this study based on two evacuation durations and two V values

Risk-to-life-safety: development of an Auckland volcanic field bayesian event Tree for volcanic hazard short term

After evaluating various available approaches, BET_EF was selected as an appropriate model for generating short-term volcanic hazard data for a crisis decision-support tool in the AVF context (Wild et al. 2020). It can assess both the eruption probability and spatial vent likelihood in one model and is available as a software application, which was considered necessary for potential operational uptake. However, given the spatial variability in eruptive style and hazardous phenomena when assessing the risk-to-life-safety at a location, and the need for this to inform a CBA for evacuation decision-support, the previous BET_EF (Wild et al. 2022) is here extended to the BET_VHst framework for the AVF (hereinafter, BETVHst_AVF). This allows for complete alignment with Eq. 16 from Marzocchi and Woo (2009):

$${P}_{Evac}^{\left(X\right)}={p}_{1}{p}_{2}{p}_{3}\sum _{k}{p}_{4}^{\left(k\right)}{p}_{5}^{\left(S\right)}{p}_{6}^{\left(S\right)}{p}_{7}^{\left(k,S, X\right)}$$

(5)

where p_1to7 are outputs from BETVHst_AVF, k is a vent within the field, S is for eruption style size, and X is a location within the defined study extent. This returns the probability ${(P}_{Evac}^{\left(X\right)}$) of needing to evacuate location X, subsequently used with Eq. 1.

Node 1 to 4: eruption and spatial vent probability

The monitoring and prior model parameters from the recently updated BETEF_AVF (Wild et al. 2022) are applied here for Nodes 1 to 4. The uniform prior is used at Node 4 to generate the outputs. Seismicity is used as input for the monitoring component for estimating the vent location, whereby each event is weighted inversely proportional to its depth in an isotropic Gaussian kernel density. The field extent remains as the bounding box around the AVF boundary + 5 km buffer (Runge et al. 2015), and the field is divided into 500 m spaced cells, forming a 53 × 79 grid, i.e., 26 × 39 km.

Node 5: eruption style/size

While previous studies have illustrated the spatial variability in eruptive style across the AVF (Sandri et al. 2012; Hayes et al. 2019; Ang et al. 2020), this was not considered in the original implementation of BETEF_AVF by Wild et al. (2022), as the Auckland Volcanic Contingency Plan states an evacuation is required for everyone within 5 km of the inferred eruption centre/area (Auckland Council 2015), irrespective of eruption style. However, as the risk to life depends on the style of eruption, here we consider three eruptive styles, typical phreatomagmatic (TP), large phreatomagmatic (LP) and magmatic/effusive (M), to align with Sandri et al. (2012). While it is acknowledged that styles can transition during an AVF eruption (Kereszturi et al. 2014, 2017), only the initial eruptive phase is considered here as our study focuses on the risk-to-life safety and associated evacuation decision-making in the lead-up to an eruption.

To assess the spatial variability in phreatomagmatic eruptions (p_phr), Eq. 6 from Ang et al. (2020) is applied at each grid location of the BETEF_AVF extent. It is taken that if the eruption is not phreatomagmatic, it is magmatic; therefore, p_mag = 1 – p_phr. The probability, given it is phreatomagmatic, for TP and LP is 0.7 and 0.3, respectively, taken from the prior in Sandri et al. (2012). An ‘equivalent data’ value of 10 is applied, representing a significant degree of confidence.

Node 6: phenomena

The phenomenon that will result in the most significant risk to the population from a phreatomagmatic eruption is a base surge (Magill et al. 2005). However, for magmatic/effusive eruptions, there are several vent proximal volcanic products that can result in causalities such as edifice and/or fissure formation, scoria cone development, lava spatter and ballistics. To simplify the complexity and improve computational efficiencies, a single ‘footprint’ considering a “need to evacuate” represents the risk of casualty to a person at a proximal location, regardless of the volcanic eruptive style and products produced. A conditional probability of 1 for this node is applied for each eruption style, as without action, they would all result in casualty for an exposed person.

We assess the need to evacuate locations across Auckland using the same 500 m grid spacing as the AVF grid used to assess spatial vent likelihood (Node 4), but have extended it by 5 km to account for hazard run out distances beyond the field boundary, forming a 73 × 99 grid, i.e., 36 × 49 km.

Node 7–8: overcoming thresholds at target areas

Base surge is the primary hazard of concern for impact on life safety during a phreatomagmatic eruption. Here we have applied the same base-surge run-out probability functions for the two eruptive styles as were presented in Sandri et al. (2012). This yields the probability of base surge run out exceeding a given distance, and is used to assess the exposure to the phenomenon at a location, given the designated eruptive style. The ‘equivalent data values’ are 5 and 7 for the TP and LP, respectively, from Sandri et al. (2012). While there are examples from past eruptions that show that such phenomena could be survivable at their distal edges, albeit likely with severe body burns and respiratory issues requiring extensive medical treatment (e.g., 2019 Whakaari New Zealand eruption; Baker et al. 2021 and 2010 Merapi Baxter et al. 2017), there are not sufficient data to produce a reliable/robust survivability function. Furthermore, there are differences in susceptibility based on whether people are sheltered (e.g., inside buildings) or exposed in the open (Jenkins et al. 2013). It is assumed that all those exposed outside would not survive, due to either intense temperature exceeding 200°C or surge concentration exceeding 0.1 kg/m³, resulting in asphyxiation and respiratory issues (Baxter et al. 1998). Interpretation of the Brand et al. (2014) numerical modelling of surge parameters for both LP and TP scenarios shows that at least one of these thresholds is still exceeded at the end of the run-out distance. Given the limited data, the variability and associated uncertainty in whether people are inside or outside, we have applied the conservative assumption that all exposed people are outside, resulting in a binary impact as a function of hazard exceedance of base-surge run-out (Fig. 3).

Fig. 3

Probability of exposure to base surge as a function of distance from the vent. Distribution for “typical phreatomagmatic” and ‘large phreatomagmatic” after Sandri et al. (2012) Fig. 4. Magmatic distribution is based on past eruptions – refer to the text for more detail. In our analysis, exposure to base surge results is assumed to result in casualty (binary impact)

Analogous past eruptions inform the magmatic eruptive style function in the absence of empirical data for the AVF. The 1973 Heimaey eruption started with a fissure, formed during the first hours of the eruption, with a length of 1.5 km onshore but potentially extending to 3.5 km offshore, with fire fountaining and lava splatter proximal to the opening (Thorarinsson et al. 1973; Williams and Moore 1983). During the initial phase of the 1943 Parícutin eruption, a cinder cone reaching 50 m high was formed due to lapilli and bombs falling 300-400m from the vent, while lava travelled 700 m from the vent in the first 24 hours (Foshag and González-Reyna 1956; Rees 1979). The phenomena associated with these analogous eruptions align with expert-derived feasible magmatic scenarios for the AVF (Hayes et al. 2018, 2019). Impact from magmatic style eruptive phenomena is assumed to result in fatality. As such, only the hazard exceedance of magmatic style reaching a distance from the vent is required (Fig. 3). An ‘equivalent data value’ of 1 is used to represent the maximum uncertainty for the input in the BET_VHst, given that the function is based on limited past analogue eruptions.

Weight of the monitoring component at BET node 4

As discussed above, when assessing the spatial vent likelihood in the case of distributed volcanism it could be advantageous to modify the weight of the monitoring component at Node 4 of BET_EF and BET_VHst. To demonstrate the proposed approach of a transitional weighting for the monitoring component at the Node 4 parameter, three interval ranges for the weight will be considered:

• 0.2 to 0.8,

• 0.1 to 0.9,

• 0 to 1.

The scaled weight should be positively correlated with the increase in eruption likelihood, because the purpose of this parameter is to reflect escalating unrest and to modify how Node 4 outputs are assessed, which is not to be confused with the Node 4 inputs.

To demonstrate this approach, the parameter selected is the Node 3 parameter 2 from BETVHst_AVF (refer to Wild et al. 2022, for more information). This parameter assesses the number of volcano-tectonic earthquakes with M_L>2 not attributed to a mainshock/aftershock sequence within the AVF GeoNet monitoring extent, with a lower and upper threshold of 10 (normal activity) and 100 (highly suggestive of an impending eruption), respectively. This parameter is selected as the presence of VT events is considered to suggest magma movement in the crust, with an increase being an indicator of escalating volcanic unrest (Lindsay et al. 2010). Additionally, seismicity is constantly monitored through the Auckland seismic monitoring network, whereas there is currently no permanent geochemistry or geodesy monitoring equipment within the AVF.

To examine the effect of the change in monitoring weights on the spatial vent likelihood outputs, we will use the BETEF_AVF developed by Wild et al. (2022) and apply the AVF synthetic unrest sequence for the Birkenhead Scenario from Hayes et al. (2018). This will be conducted for the three monitoring component weight ranges and compared to the default BET_EF Node 4 weight of 0.5, when the volcano is determined to be in unrest (i.e., $\eta =1$).

AVF evacuation decision-support testing

The synthetic unrest sequence from the Birkenhead scenario (Hayes et al. 2018, 2019) is used to review the performance of the integration of CBA with the BET_VHst to inform a decision-support approach for the AVF. The Birkenhead scenario pre-eruptive sequence occurs across a 15-day period, with observable phenomena limited to just seismic activity and ground cracking, with 901 earthquakes starting at 29 km before shallowing. Refer to Hayes et al. (2018) for more information on the scenario’s pre-eruptive sequences.

At the time of writing, there is no downloadable BET_VHst application. To implement the BET_VHst framework, the output of Nodes 1 to 4 are assessed using BETEF_AVF, implemented using the PyBETUnrest tool (Tonini et al. 2016), with the BET_VH described above. In the BET_EF framework, the a and b prior distribution parameters are used to assess the degree of unrest at Nodes 1 to 3. As in BETEF_AVF, to present the outputs here, a = U(0, 2) and b = U(0.75, 1) are used to calculate the outputs for up to Nodes 1 to 4 from PyBETUnrest. The PyBETUnrest code is modified to add a transitional monitoring component weight at Node 4. This code base is run for each of the three Node 4 transitional monitoring component weight intervals for the Birkenhead unrest sequence, along with the default 0.5 weight used in Wild et al. (2022). PyBETVH (Tonini et al. 2015) is used to produce outputs for Nodes 5 to 7/8. The Node 1–3 probability of eruption and Node 4 conditional vent location likelihood inputs in BET_VH are the respective outputs from BETEF_AVF. This produces a series of outputs representing the probability of fatality at a given location for each of the Node 4 monitoring weights.

The four presented CBA thresholds (Table 1) are used with the outputs of the BETVHst_AVF for the Birkenhead scenario. This will identify which locations are deemed cost-beneficial to evacuate based on each threshold, thus identifying an evacuation zone. For comparison, we generate the evacuation zone following the AVF Contingency Plan of 5 km around the vent uncertainty area (Auckland Council 2015). In this instance, we will define the required vent uncertainty area as either the 50th or 95th percentile output area produced from Node 4 from BETEF_AVF.

Using the developed evacuation zones across the unrest sequence in the scenario, the population and subsequent evacuation clearance time can be evaluated by applying the methodology presented in Wild et al. (2021) in which the exposed population is defined as those within the given evacuation zone. The population data used are from the Statistics New Zealand Statistical Area Two (SA2) 2018 census dataset (Statistics New Zealand 2018), which represents a “semi-suburb” resolution, considered an appropriate size for operational evacuation zone boundaries. The Wild et al. (2021) evacuation clearance time methodology is divided into a pre-travel phase of 36 h, which is spatially independent, and the travel time phase, which is spatially dependent on the evacuation zone. The travel time is modelled using a high-level geospatial approach, combining private transport ownership data with road network data and vehicle carrying capacity.

Despite an evacuation call, there is some residual risk for the non-evacuated population beyond the evacuation zone. ‘Residual statistical fatalities’ can thus be assessed as the product of the population at a location (using the SA2 population dataset and excluding those within the evacuation zone) and the BETVHst_AVF estimated probability of impact at a location. The output of this calculation shows the potential fatalities if the eventual vent location generates hazardous phenomena that extend into areas beyond where the BETVHst_AVF output does not exceed the CBA threshold. When run for every day of the sequence, the influence of the evacuation zone and monitoring weights on the residual risk to the surrounding population if an eruption hypothetically occurred on each day of the sequence can thus be determined.

Results

We calculated the absolute probabilities of an eruption at a location, i.e., Node 4, from BETEF_AVF for the Birkenhead scenario using the four considered monitoring component weights described in Sect. 4.3 (see supplementary material). As in Wild et al. (2022), the vent outbreak eventuated in the centre of the identified likely vent hotspot (warm-coloured cells) for all transitional monitoring component weights. However, compared to the 0.5 monitoring component weight, the estimated vent location is more tightly constrained in the transitionally weighted applications, given the greater emphasis on the monitoring component in later stages of the eruption.

Figure 4 presents evacuation zones identified for selected days of the Birkenhead Scenario comparing the use of BETVHst_AVF, which produces the probability of impact at a location, with BETEF_AVF, which sums the probability of an eruption from all possible vents within 5 km of a location, ignoring eruptive style and phenomena (as in Wild et al. 2022). To present the results, the 0.1 to 0.9 interval for the transitional monitoring component weight is used for assessing the spatial vent likelihood and the 3-month evacuation duration with V = NZTA VOSL CBA threshold (Table 1). To further facilitate comparison, the evacuation zone following the AVF Contingency Plan, i.e., 5 km buffer around the vent uncertainty area, is also presented for comparison. The vent uncertainty area is represented by the 50th and 95th percentile output area from BETEF_AVF to represent the vent uncertainty. As shown, the evacuation zones are smaller using BETVHst_AVF than BETEF_AVF, due to the decay in probabilities of phenomena exceedance as the distance from the vent increases. The AVF contingency plan defined evacuation zone formed using the 95th percentile output vent uncertainty is significantly larger than the others. However, the evacuation zones defined by the AVF Contingency Plan defined evacuation zone formed using the 50th percentile output vent uncertainty are similar to the output of the BETVHst_AVF from 7 days prior to the eruption.

Fig. 4

Evacuation zone extents formed using BETVHst_AVF and BETEF_AVF with 0.1 to 0.9 transitional weighting at Node 4. Additionally, evacuation zones are presented by applying the AVF Contingency Plan using the 50th and 95th percentile vent uncertainty area for the Birkenhead scenario + 5 km evacuation zone buffer. All examples use the CBA threshold of p = 0.00826. The base probability grid is the spatial distribution of the absolute short-term average probability of volcanic hazard impact at a location from BETVHst_AVF.

Figures 5 and 6 present for the default 0.5 and 0.1 to 0.9 Node 4 monitoring component weightings, respectively, the output spatial probability of needing to evacuate a location due to risk-to-life-safety based on the BETVHst_AVF, and the evacuation zones defined with each of the four presented CBA thresholds (Table 1). The outputs for the transitional weights of 0.2 to 0.8 and 0 to 1 are presented in the supplementary material. As expected, for both Node 4 monitoring components, the 1-month evacuation duration with V = NZ$4.46 million with the smallest CBA threshold yields the largest evacuation zone. However, what is counterintuitive is that, for the three smaller CBA thresholds with the 0.5 Node 4 weighting, the evacuation zones extend to include areas beyond the focus of unrest as the eruption approaches. The increase beyond the monitoring area is influenced by the likely eruption style, as those locations first included have a higher likelihood of being a phreatomagmatic style.

Figure 7 presents the derived evacuation zone area for each of the Node 4 monitoring component weightings for each of the four CBA thresholds across the Birkenhead scenario unrest sequence. For the 0.5 monitoring component weight, it is shown that in all cases except the smaller CBA threshold, the area increases as the unrest sequence develops. For the three transitional Node 4 monitoring weight intervals, except for the 0.2 to 0.8 set with the 1-month and V = NZTA VOSL CBA threshold, areas stabilise around 10 days before eruption. For the case of a 1-month duration evacuation with V = NZTA VOSL CBA threshold, the evacuation area spiked on day 11 before the eruption when the monitoring component weight was ~ 0.5 for all the interval sets, before decreasing with the increased importance of the monitoring component within the BETEF_AVF to assess the vent location.

Fig. 5

Evacuation zones defined for each of the CBA approaches using the spatial distribution of the absolute short-term average probability of volcanic hazard impact at a location from BETVHst_AVF where the Node 4 weight of the monitoring component is 0.5 for selected days leading up to the eruption in the Birkenhead scenario

Fig. 6

Evacuation zones defined for each of the CBA approaches using the spatial distribution of the absolute short-term average probability of volcanic hazard impact at a location from BETVHst_AVF where the Node 4 weight of the monitoring component transitions from 0.1 to 0.9 for selected days leading up to the eruption in the Birkenhead scenario

Fig. 7

Evacuation zone areas defined for each of the CBA approaches using the spatial distribution of absolute short-term average probability with vent location monitoring component weights of 0.5, 0.2 to 0.8, 0.1 to 0.9, and 0 to 1 for the eruption in the Birkenhead scenario

The population exposure (i.e., the number needing to evacuate) and corresponding evacuation clearance time are assessed based on the defined evacuation zone for each day across the sequence (Fig. 8). For the two 3-month CBA thresholds, the population exposure and evacuation clearance times follow similar trends for all weighting options for the Node 4 monitoring component. The population exposure is the least for the 0.5 monitoring component weight as the evacuation zone is the smallest compared to the other weights. The difference in population exposure between the two 3-month CBA thresholds is driven by the evacuation zone area (Fig. 7).

In contrast, the two 1-month CBA thresholds result in very different exposure outputs for the four monitoring component weightings. When considering the CBA threshold with V = NZ$1.8 million, the three transitional monitoring component weightings all follow the same trend, with minimal differences in population and evacuation clearance times, steadying from 10 days prior to eruption. However, with the 0.5 monitoring component weighting, the population evacuation constantly trends up because of the growing evacuation zones (Fig. 7). For the smallest CBA threshold, i.e., 1-month duration with V = NZ$4.46 million, there is a diverging split in population exposure (Fig. 8) as the evacuation area grows for the 0.5 and 0.2 to 0.8 transitional vent weights as more locations are identified as cost beneficial to evacuate as P(Erup) increases, despite vent likelihood converging to a location.

Figure 9 presents the ‘residual statistical fatalities’ across the Birkenhead scenario unrest sequence for the different evacuation zones and V values evaluated in this study. On the day before the eruption, these residual statistical fatalities range from ~ 300 to ~ 1,800 (using the transitional monitoring weight of 0.1 to 0.9 at Node 4), and from ~ 25 to ~ 5,900 (using the monitoring weight of 0.5 at Node 4), based on the evacuation zones presented in Fig. 4.

Fig. 8

Population exposure and median evacuation clearance times for each of the CBA approaches for all four considered Node 4 monitoring component weights based on the evacuation zone on that day of the scenario unrest sequence

Fig. 9

Residual statistical fatalities for each of the presented evacuation zone boundaries across the Birkenhead scenario unrest sequence for a 3-month evacuation duration. Results are shown for Node 4 monitoring component weights of 0.1 to 0.9 (left) and 0.5 (right)

Discussion

Spatial-temporal eruption forecasting using BET

A review of the weight of the monitoring component at Node 4 was recommended by Wild et al. (2022) to help assess the spatial vent likelihood for distributed volcanism. As expected, irrespective of the weighting approach used, vent likelihood converged on the same region in our Birkenhead synthetic unrest scenario at Node 4 (Eq. 3). However, the weight of the monitoring and prior components did indeed influence the specific spatial vent likelihood and any outputs at subsequent Nodes of a BET, given subsequent outputs are conditional on the vent location probability. Therefore, with a larger weight for the non-monitoring component, the likelihood of possible vent locations beyond the focus of unrest area is increased. This is because the component representing the prior has a greater influence, as the estimated eruption probability (Node 3) increases the likelihood of eruption everywhere.

Marzocchi et al. (2008) applied a maximum weight for the monitoring component of 0.5 (Eq. 4). This was a deliberate decision to not have outputs fully informed by the monitoring inputs (in contrast to the previous Nodes 1 to 3) in order to account for potential uncertainty and to mitigate false vent localisation that could arise due to possible lateral ascent of magma or migration of seismicity. However, BET_EF was initially developed primarily for stratovolcanoes, where a central vent is commonly the most likely future vent. In contrast, distributed volcanoes typically have many possible vent locations, resulting in the vent location probability being significantly more spatially distributed in the prior model. This is even more pertinent for volcanoes such as the AVF that have diffuse spatial patterns and/or erratic spatio-temporal trends to inform the prior model (Bebbington and Cronin 2011; Bebbington 2013). Given the likely variability across the field of both the prior and the monitoring component, outputs may not converge as significantly onto a likely vent or general area if only 0.5 weighting is given to the monitoring component. This was illustrated in the case of the AVF when applying a uniform prior (Wild et al. 2022). To address this, we undertook a sensitivity analysis of the influence that the Node 4 monitoring component weight has on BET outputs and evacuation zones thus identified. While setting the monitoring component weight to 1 provides the greatest vent likelihood localisation from a BET_EF, we agree with Marzocchi et al. (2008) that this should not be done. Some prior contribution should be retained to account for the influence of the subsurface structure and its influence on vent likelihood (e.g., to account for the possibility of lateral magma movement as indicated above) and to remove the risk of false localisation by relying only on the monitoring observations.

However, the monitoring component weight at Node 4 needs to be increased to provide evacuation zones that stabilise prior to eruption, and so we reviewed the use of a transitional weight between two values. This means that in the initial phases, when vent uncertainty is large and a limited number of monitoring observations, the prior models dominate the spatial posterior. The transition from a prior-dominated to a monitoring-observation-dominated weight as unrest increases reduces the influence of a single seismic event at depth on the output vent likelihood. In the Birkenhead scenario example above, the parameter selected to inform the monitoring component transition is based on the number of VT events. The weight of the monitoring component does not increase from the lower threshold until 12 days before the eruption. Prior to that time, the spatial probability output from BETVHst_AVF remains relatively uniform for the three interval ranges (Sect. 4.3), in contrast to the outputs to the fixed 0.5 monitoring component weight (Eq. 4). The upper threshold is exceeded 9 days from eruption when there are more than 100 VT events, considered a significant enough number that using these would allow for confidence in modelling vent localisation. A similar approach can be reached by applying judgement to modify the spatial monitoring input in Node 4. However, this introduces manual intervention and subjectivity around when that level of confidence is reached, which in the midst of an event loses one of the main benefits of a quantitative forecasting approach such as a BET, namely to remove the pressure from monitoring volcanologists.

While here we have used the number of VTs to inform the transitional parameter for the monitoring component weight, other strategies could be adopted, e.g., cumulative energy release or number of LP events. To demonstrate the approach for this study, we selected a parameter already within the BET_EF at Node 3. However, a composite parameter including other observations could be advantageous and warrants further investigation. Possibilities include observation of deformation or localised gas recording, or perhaps the estimated eruption probability, which in a BET is a result of synthesis of all the monitoring data. However, non-seismic monitoring observations are challenging to process in real-time, for the purpose of tracking moving in the subsurface. Additionally, P(Erup) is already used to produce the absolute probability of eruption at a vent location and ideally the spatial vent likelihood should be to some degree decoupled.

From the results of this study, we recommend a movement to a transitional parameter at Node 4 for the monitoring component weight for distributed volcanism. The weighting split between the two thresholds should be determined based on the confidence in spatial monitoring capability and understanding of the volcano. Here we have presented three ranges of thresholds to test the sensitivity this has on Node 4 outputs based on a synthetic unrest sequence to review the influence on outputs at the later Nodes within a BET_VHst. As with the Node 1 to 3 monitoring parameters, informed selection is required to select the parameter(s) and upper and lower thresholds during the setup of a BET_EF or BET_VHst.

As demonstrated across both this study and Wild et al. (2022), the sensitivity analysis for key parameters requires further review from monitoring scientists before a BET could be operationally applied for the AVF. Furthermore, across this study, a single synthetic scenario has been used for testing, and other scenarios, including failed eruptions should be reviewed. Even if operationally applied, its use would likely be helpful to support volcanologists and decision-makers during a time of crisis, rather than be prescriptive.

Applying cost-benefit analysis and BET_VHst as an evacuation decision-support for the AVF

Previous studies integrated a BET with CBA to support evacuation decision-making (Marzocchi and Woo 2009; Sandri et al. 2012; Wild et al. 2019b). Here we constructed BETVHst_AVF and used the outputs representing the probability of fatality for a person at a location combined with CBA to inform evacuation zones during the lead-up to an AVF eruption. The extension of the BETEF_AVF to BETVHst_AVF provides the output probability of adverse volcanic hazard impact (representing the probability of fatality in this case study) at a location, which is a more appropriate value for p within the CBA to define the need to evacuate, compared to the cumulative P(Erup) from a vent within a 5 km radius that was presented in Wild et al. (2022). In all cases, the model performed well for the reviewed Node 4 weightings, and the four reviewed CBA thresholds, as the evacuation zones derived were all smaller than those formed using the BETEF_AVF, thereby reducing the number of people who would need to evacuate to a more realistic number.

The Node 4 monitoring weight impacts the size of the evacuation zone, as there remains a residual vent likelihood in the areas beyond focus of unrest due to the effect of the prior. This results in a greater BETVHst_AVF maximum probability of volcanic hazard outputs with the increase in monitoring weight, which subsequently affects the evacuation zone based on the CBA threshold applied. In the case of the 3-month duration with V = 1.8 million, the smallest CBA threshold, across the assessed Node 4 monitoring weights, all evacuation zones are consistent in shape. The 0.5 monitoring component weight produces the smaller evacuation zone, as there remains an inflated vent likelihood beyond the focus of unrest. This is due to the greater influence of the prior, but the probabilities at locations beyond the focus of unrest are not significant enough to exceed the CBA thresholds. In contrast, the evacuation zone derived using the smaller CBA thresholds is shown as more sensitive to the BETVHst_AVF Node 4 monitoring component weight. When the monitoring component weight is set to 0.5, the evacuation zones increase in size, similar to the results observed when using BETEF_AVF to identify where to evacuate (Wild et al. 2022). This is because the prior inflates the vent likelihood exceeding the CBA threshold in areas beyond the focus of unrest, informed by the seismic monitoring observations. This occurred from 10 to 6 days from eruption for V as NZ$1.8 million and NZ$4.46 million, respectively. However, instead of the evacuation zone expanding radially, as in Wild et al. (2022), the direction of expansion is informed by the eruptive style at Node 5. The zone grows where the initial eruptive style is more likely to be phreatomagmatic (Ang et al. 2020). This is because phreatomagmatic phenomena are modelled to have a greater run out (Fig. 3), and as such, locations where the dominant eruptive style is phreatomagmatic have a greater probability of impact. As the P(Erup) increases, more locations are identified as cost-beneficial to evacuate beyond the focus of unrest, thereby increasing the evacuation zone (Fig. 10). The Birkenhead scenario is limited to seismic unrest phenomena, thereby restricting the Z-value at Nodes 1 to 3, preventing a greater P(Erup) value (Wild et al. 2022), but a different unrest sequence (for example, the Rūaumoko scenario, Fig. 4 in Wild et al. (2022) can have pre-eruptive values of P(Erup) nearer 1.

With the transitional monitoring component weights, the counterintuitive expanding evacuation zone with decreasing time to eruption is less often observed. With both the 0.1 to 0.9 and 0 to 1 intervals, the evacuation zones remain focused around the focus of unrest. For the 0.2 to 0.8 range, at the smallest CBA threshold, i.e., 1-month duration with V = NZ$1.8 million, the evacuation zone starts to include areas beyond the focus of unrest 8 days from the eruption. This suggests that the smaller the CBA threshold probability for evacuation, the higher the maximum value of h at Node 4 needs to be in order to produce a stable evacuation zone.

Fig. 10

Evacuation zones defined for each of the CBA approaches using the spatial distribution of the absolute short-term average probability for Node 4 monitoring component weight of 0.5 for (A) P(Erup) = 0.7 and (B) P(Erup) = 1 and Node 4 monitoring component weight of 0.1 to 0.9 for (C) P(Erup) = 0.7 and (D) P(Erup) = 1 for requiring to evacuate a location from BETVHst_AVF on the day before the eruption in the Birkenhead scenario

Comparing the evacuation zones applying BETVHst_AVF with CBA against those applying the AVF Contingency Plan using the 50th percentile vent uncertainty found they converge as the eruption nears (Fig. 4). In other words, as the scenario’s unrest sequence escalates, the vent likelihood starts to converge on an estimated eruptive area. However, at 13 days before the eruption, both tested AVF Contingency Plan evacuation zones are large in spatial extent, as there is little confidence in the likely vent location. In contrast, BETVHst_AVF does not define an evacuation zone until 10 days before the eruption, as the P(Erup) has not reached a high enough value for the CBA to be exceeded. While here we have selected for review purposes the 50th and 95th percentile vent uncertainty area using BETEF_AVF output, there is no current documented approach within the AVF Contingency Plan on defining a vent uncertainty area. As such, this requires expert judgement to determine when the eruption likelihood is significant, and when there is high enough confidence to define a vent uncertainty area. Using a quantitative approach, as shown here, removes the pressure on the advisory group and decision-makers to produce finely balanced judgements during a crisis. The fact that the AVF Contingency Plan evacuation zone using the 50th percentile vent uncertainty area and BETVHst_AVF with the CBA threshold selected yield similar evacuation zones suggests some merit in using BETVHst_AVF with CBA as a check to help inform decision-makers in a future AVF event.

While the presented decision-support approach combining BETVHst_AVF with CBA performs well in dynamically identifying evacuation zones, there are limitations. The applied CBA approach requires as one of its inputs a nominal (false positive) evacuation duration when in reality this is unknown until the end of an unrest period, and of course it will be difficult to determine when a ‘false-alarm’ “unrest” period ends if no eruption ensues. As shown, different evacuation durations affect the extent of the evacuation zone area.

BETVHst_AVF, in its present form, considers a uniform likelihood of hazard exceedance in all directions as a function of distance, irrespective of topography and built environment, which can affect the hazard extent and intensity of phenomena. Base-surge modelling for each grid location would enhance the input resolution for Node 7/8 of the BETVHst_AVF, allowing for output hazard information such as dynamic pressure and temperature. This would improve the assessment of building fragility and likelihood of fatality when combined with data such as construction type. Here we have taken a conservative approach, as New Zealand has a low risk tolerance regarding life safety, as demonstrated in policy and legislation (e.g., Auckland Council 2014; CDEM Act 2002; MCDEM 2015) and in the standard requirements for construction and land-use planning (e.g., Saunders et al. 2013; Ministry for Business, Innovation and Employment 2021). In areas with different risk tolerances, one could explore the consequences of shelter for those at distal margins of hazard footprints, although we emphasise there is limited reliable data on survival rates for at the extent of base-surges, both in and out of buildings.

Considerations for cost-benefit analysis as an evacuation decision-support approach

CBA provides a systematic framework that provides a binary yes-no decision. Nevertheless, as demonstrated in our application to AVF, there are a number of limitations, both in terms of the approach and parameter selection. Firstly, the CBA approach applied in this study is a static criterion. This means it is only cost-beneficial to evacuate when p is exceeded. In the presented hypothetical Birkenhead scenario, a stable evacuation zone is formed early enough so that the required population has sufficient time to clear. However, this of course is only apparent in hindsight in this presented scenario, and we stress that an evacuation based on the CBA approach should be called at the time it becomes cost-beneficial, rather than waiting for any reduction or stabilisation of evacuation zones. In some circumstances, such as a much shorter sequence or where there is significant uncertainty, this could potentially leave insufficient time to evacuate those at risk, as the required clearance time is not considered. Bebbington and Zitikis (2016) presented a dynamic CBA approach which considers time-dependent factors such as the uncertainty around the eruption onset time and the distribution of time to evacuate. However, the Bebbington and Zitikis (2016) approach would need to be modified to reflect the complexity to incorporate the spatial component required for the AVF.

CBA lacks the ability to include intangible costs that come with an evacuation. Such factors that should be included as part of decision-making include the mental and cultural costs of an evacuation, and the removal of the sense of community cohesion that comes with being displaced (Gottsmann et al. 2019). Furthermore, displacing a population can have secondary health implications, as observed following the 1991 Pinatubo evacuation, with a measles outbreak occurring in one of the evacuation centres (Bautista 1996; Floret et al. 2006). All of these factors require consideration as either a ‘cost’ or ‘loss’ from a decision, but are often not directly attributable to an economic value, limiting their inclusion with a CBA. This has led to the consideration of broader holistic decision-making approaches, such as Multi-criteria analysis (MCA), which consider more than what is directly attributable to an economic ‘cost’ or ‘loss’ (Ackerman 2008). These holistic approaches do not align with a fully quantifiable approach, such as probabilistic eruption forecasting, and have not been explored for use in volcanic crisis evacuation management. However, as with CBA, they can be developed and reviewed before an event.

In the presented approach for assessing CBA, the proportion of people that owe their life to the evacuation (E) is a fixed value. Here we set E to 0.65 based on Thakur et al. (2022), as we estimate 35% would self-evacuate in advance of the official evacuation call, following the issuance of advisory messages. However, this issue is complex, with several societal and environmental factors likely to influence behaviour in variable ways across any eruption sequence. For example, it is expected that both reported and felt earthquakes will influence Auckland residents’ behaviour (Coomer et al. 2015). Risk perception and understanding of the volcanic hazards and likely impacts will affect individual and household behaviour, such as when they will initiate evacuating (Coomer et al. 2015), as will social cues, such as family or neighbours evacuating (Baker 1991). As the unrest sequence develops, but before the mandatory evacuation order is issued, the proportion that has self-evacuated is likely to increase, thus decreasing the value of E. Furthermore, shadow evacuations, described as those self-evacuations that occur from outside the evacuation zone, have also been observed in past crises (Dow and Cutter 2002; Dash and Gladwin 2007) and are expected in a future AVF event (Coomer et al. 2015; Thakur et al. 2022). This will increase pressure on evacuation resources, congestion on transport routes and subsequently the clearance times, but is not captured within the CBA approach applied here. On the other hand, non-compliance with a mandatory evacuation order can also occur, and has been observed in past volcanic crises (Tayag et al. 1996; Mei and Lavigne 2012; Elissondo et al. 2016; Bird and Gísladóttir 2018; Jumadi et al. 2018; Lavigne et al. 2018). This can often be attributed to a reluctance to leave property, possessions and pets, inexperience with the type of hazard and its impacts, and religious grounds (Blong 1984; Lindell and Perry 1992; Cola 1996; Tobin and Whiteford 2002). Non-compliance to an evacuation order is not considered within CBA, as it cannot be accounted for in either the cost of evacuating (C) nor lives saved from the evacuation call (L). Notwithstanding all of the above, past work suggests that there would be a very high compliance rate in Auckland for those within an evacuation zone (Horrocks 2008; Coomer et al. 2015).

There are ethical considerations for the development of a CBA for crisis management. An example of this is the value-of-human life parameter (V). Reviewing the examples in which CBA is applied to volcanic crisis management, we note that this value changes based on the approach used to inform its selection. There are ethical concerns around the ability to quantify the value of human life in financial terms, with some weighting the importance of different demographics such as age or location differently (Baker et al. 2008; Fischhoff 2015). Furthermore, the approach selected to select V affects the derived CBA threshold. For example, the GDP per capita is NZ$57,000 per annum for Waikato (Statistics New Zealand 2021), another volcanic region in New Zealand. With V set to NZ$4.46 million from the NZTA VOSL, it results in a smaller CBA threshold, thereby implicitly stating the risk tolerance for Waikato is lower than for Auckland, as an evacuation would be called for a location with a lower likelihood of fatality. This shows the benefit of alternative approaches to calculate V, such as the average national contribution of GDP over a lifetime (e.g., 25 years). When applying the GDP over a lifetime parameter for V, for a 3-month evacuation, R/V = 0.04 (not considering the accommodation costs) and the CBA threshold becomes 0.04E, thereby removing the GDP value considerations, which forms an equitable and consistent threshold. The willingness-to-pay approach to relocate is an alternate metric for V (Marzocchi and Woo 2009; Marzocchi et al. 2012), but this can be dependent on the financial means available and political priorities at the time. However, attributing a financial value to a human life is incompatible with Kantian ethics as it undermines the intrinsic worth, dignity, autonomy, and universal moral principles that prioritize treating individuals as ends in themselves rather than means to an end. This raises a challenge in applying CBA for decision-making involving humans or other elements which do not have a financial value (Ackerman 2008). Further research is also warranted to explore the applicability of the CBA approach applied in potentially related contexts, such as adventure tourism. In the latter, the costs and losses are often borne by different entities.

Another way of illustrating the challenges associated with the V value is by looking at the different approaches for defining evacuation zones presented in Fig. 4, and back-calculating what it would be given their different extents. By taking the median value that intercepts the evacuation zone boundary from the BETVHst_AVF, taken to represent the threshold required to call an evacuation, a value for p can be input into the CBA equation to output corresponding V values. This demonstrates the variability in V for each of the evacuation zone approaches (Table 2). While not directly applied to inform the evacuation zone, when defining the vent uncertainty area by the 95th percentile area, V is 25-times greater than the NZTA VOSL and 63-times greater than the 25-year average GDP per capita. In contrast, with a 50th percentile vent uncertainty with the AVF Contingency Plan, V is only 75% of the NZTA VOSL. This illustrates the variability and subjective nature of the value of life parameter under the different selected evacuation zone approaches and for different contextual settings.

Table 2 V for day 1 in Fig. 7 using a 3-month evacuation duration

The cost parameter for the proposed evacuation, i.e. value for R in Eq. 2, is limited to only considering the socio-economic loss per capita for the evacuation duration, noting a detailed economic assessment of the costs and losses is beyond the scope of this study. However, other economic costs need to be more broadly considered as a result of an evacuation call, such as the longer-term broader regional economic impacts from a future eruption (e.g., McDonald et al. 2017). Woo (2008) noted the scale independence fails when a site of significant economic importance is impacted by the evacuation. For example, in the Birkenhead scenario, the proposed evacuation zone included the Auckland Port, a vital import and export location for New Zealand, which adds NZ$1.4 billion per year to the GDP, with benefits that extend beyond Auckland (Maralani and Wilson 2019). Furthermore, operations at the port could be impacted, as staff or freight movement could be impeded by an evacuation zone including or restricting access to the port. The CBA approach does not allow the inclusion of these broader economic costs associated with an evacuation call. Nevertheless, given their relative importance, bringing these costs into any CBA decision-making framework would improve the tool’s effectiveness. Adding these additional costs to the numerator of Eq. 2 also requires including the number of evacuees (N). Given the spatial population distribution, this results in an inconsistent CBA threshold across Auckland. A new AVF vent that is located offshore would likely result in the need to evacuate a much smaller proportion of the population (Wild et al. 2021). If the new vent uncertainty area was near the port, this would increase the influence of the economic loss and thus require a greater threshold to call the evacuation for those people compared to evacuating a more populous area. While this could be beneficial in practice, as it could be assumed that evacuating fewer people would require less time and resources, it is still questionable in terms of equity to have a higher threshold for one region than for another. Other considerations for assessing R could include individual or household evacuation destinations, and whether evacuees can remain working in their new location. If an individual evacuates intra-regionally and is someone who can work remotely, it is considered there is no loss of GDP for both New Zealand and Auckland. In contrast, someone who works on-site and is unable to do their job may leave Auckland. Following the 2011 Christchurch Earthquake, the New Zealand Government provided support payments to those individuals and employers affected (New Zealand Government 2011; Fischer-Smith 2013), demonstrating potential funding that may be given to support AVF evacuees; this would inflate the costs of the evacuation. In this study, to include the cost of financial support to evacuees for accommodation within C, a fixed proportion of people requiring support (65%) was applied, irrespective of socio-demographic factors of evacuees, which vary across Auckland. Alternatively, the proportion of people requiring financial accommodation support could be assessed at a finer resolution based on a number of elements, such as age, and/or income (from census data) or by using the Social Deprivation Index (Wood et al. 2017; Atkinson et al. 2020). This would potentially provide a more accurate sheltering cost for an evacuated population. However, while these elements could be incorporated into assessing R, this could create spatially variable CBA thresholds across Auckland, which again would be ethically questionable. This demonstrates the importance of achieving balance, within a CBA, between detailing nuance of the economic parameters and maintaining equity.

Another potential ethical challenge in using CBA in this context is that it assumes that some residual risk is accepted, and as such, a proportion of casualties is tolerated. While the locations in which the CBA threshold is not exceeded are by definition in areas with a reduced probability of fatality, i.e., farther from the vent, there are still residual ‘statistical fatalities’ (Fig. 9). This demonstrates the considerations required for complex volcanic settings such as volcanic fields, and those with long repose periods, where an eruption and events leading up to it have not been previously observed, putting potentially a significant number of people considered “safe” at risk. The use of CBA and the acceptance of fatalities may contradict and not align with the risk appetite in many contexts. In the case of New Zealand legislation and low societal and political risk tolerance, fatalities are not considered acceptable (Deligne et al. 2018). This demonstrates the ethical and political balance between the risk tolerance setting and the conservatism from preventing needless evacuations, albeit only known in hindsight.

The use of CBA approach applied here, while having its limitations, provides a balance between regional nuance while remaining equitable for forming an Auckland per-capita risk threshold for determining when and where to evacuate. This is achieved by applying regional or national values, representing the average across the population, irrespective of the number of people to evacuate, and the socio-economic or demographic differences at a more local scale. However, some groups may be more vulnerable during an evacuation. Groups that cannot self-evacuate (such as those in hospitals, care homes or prisons) require greater logistics, thereby increasing evacuation clearance times. Evacuations for these populations would likely need to commence earlier than for the general population in Auckland (Wild et al. 2019a), which could be achieved by using a lower CBA threshold to initiate evacuations earlier and potentially larger region. In this paper, we present an example parameter selection for the Woo (2008) CBA approach to demonstrate the application for evacuation decision-support, when combined with BET_VHst for the AVF. As discussed, this CBA approach has limitations, particularly in the selection of suitable parameters. Engagement with a range of experts, including emergency managers, policy makers and lawyers, may lead to the identification of alternative values or other factors to consider, to inform parameter selection within the CBA. Alternatively, a completely different threshold could be used with outputs from BETVHst_AVF to inform the decision of when and where to evacuate (e.g., Rogers 1975). As an example, in New Zealand and Australia, a threshold of 1 × 10^− 4 is recommended for annualised individual risk to life, where exceedance represents an intolerable risk for land-use planning purposes (Australian Geomechanics Society 2007; Taig et al. 2011; Deligne et al. 2018). Given this is smaller than the lowest CBA threshold in this assessment, if it were used as an evacuation decision threshold, it would produce larger evacuation zones than those presented here.

Given the discussed limitations with the use of CBA, we suggest it requires further review from decision-makers before being operationally applied. However, while there are known limitations, the use of CBA to support crisis decision-making removes the longstanding geoethical challenges and decision pressures on monitoring scientists and volcanologists, who traditionally focus on the spatio-temporal eruptive likelihood and advise decision-makers, rather than make the decisions themselves (Papale 2017; Peppoloni et al. 2023). As such, while this study was exploratory, looking at the integration of BET_VHst with CBA to support short-term evacuation decision-making in the AVF, given its potential, we recommend ongoing broader discussions with stakeholders to continue to explore operational applications during a future crisis.

Summary and conclusions

This paper presents the amalgamation of CBA with a BET_VHst for the AVF to be used as a decision-support tool for where and when to evacuate for a future crisis. For the CBA, we reviewed four thresholds, based on two evacuation durations and two values of human life (V). The BET_VHst was extended from the BET_EF developed in Wild et al. (2022) to produce probabilities of volcanic hazard impact across the region. The combination of these models was tested with a synthetic unrest dataset to derive evacuation zones in the lead-up to an eruption.

Based on our development and application of CBA with BETVHst_AVF, we offer the following key general conclusions regarding the integration of these models as a tool for supporting decision-makers calling an evacuation in a distributed volcanic environment:

• There is a potential benefit in changing the Node 4 monitoring component weight within BET_EF and BET_VHst when assessing the spatial vent likelihood for volcanic fields. The change in monitoring weight reduces the likelihood of the vent in regions beyond the focus of unrest. The proposed change resolves the issue identified in Wild et al. (2022) in which the evacuation zones grew as the probability of eruption increases. We propose that, in cases of distributed volcanism, a transitional weight be used as unrest escalates, to prevent false localisation with minimal monitoring observations.

• Extending the previously developed short-term eruption forecasting framework BET_EF to consider, within a BET_VHst framework, the volcanic style and phenomena produced, their extent and the likelihood of causing casualty produces outputs that are more practical for crisis management. In the context of this study, BET_VHst produces a probability of volcanic hazard impact at a given location, which aligns with the input for p for CBA. The evacuation zones were reduced in size due to the decay in the likelihood of the hazard exceedance (and associated fatality) as a function of distance from the vent, compared to the binary approach using BETEF_AVF.

• While the approach for combining BET_VHst with CBA can support decision-makers, there is complexity around the selection of the CBA parameters, which can alter the resulting thresholds and thus the evacuation zone. A balance needs to be struck between the calculated economic losses and maintaining equity across a region. The different thresholds can significantly affect the number of people and locations to evacuate.

A tool combining CBA with BET_VHst can be developed and reviewed in advance of a crisis and run in near-real time to provide an initial discussion starting point to support crisis management decision-makers when challenged with when and where to call an evacuation. Expert input from scientific, legal, and other domains is required for selecting the weight of the monitoring component and political and ethical considerations to inform the CBA component. These models need to be volcano specific and should be reviewed regularly, given advances in the understanding of a given volcano and forecasting capability, as well as the selection of the CBA parameters. As such, further work remains before such an approach can become operational.

This section presents the spatial probability of volcanic hazard impact at a location with the 0.2 to 0.8 and 0 to 1 interval sets for Node 4 monitoring component, with evacuation zones defined with each of the four presented CBA thresholds in Figs. 5 and 6 respectively. This follows the same approach as explained in the main paper.