Assessing long-term tephra fallout hazard in southern Italy from Neapolitan volcanoes

Nowadays, modeling of tephra fallout hazard is coupled with probabilistic analysis that takes into account the natural variability of the volcanic phenomena in terms of eruption probability, eruption sizes, vent position, and meteorological conditions. In this framework, we present a prototypal methodology to carry out the long-term tephra fallout hazard assessment in southern Italy from the active Neapolitan volcanoes: Somma–Vesuvius, Campi Flegrei, and Ischia.

The FALL3D model (v.8.0) has been used to run thousands of numerical simulations (1500 per eruption size class), considering the ECMWF ERA5 meteorological dataset over the last 30 years. The output in terms of tephra ground load has been processed within a new workflow for large-scale, high-resolution volcanic hazard assessment, relying on a Bayesian procedure, in order to provide the mean annual frequency with which the tephra load at the ground exceeds given critical thresholds at a target site within a 50-year exposure time. Our results are expressed in terms of absolute mean hazard maps considering different levels of aggregation, from the impact of each volcanic source and eruption size class to the quantification of the total hazard. This work provides, for the first time, a multi-volcano probabilistic hazard assessment posed by tephra fallout, comparable with those used for seismic phenomena and other natural disasters. This methodology can be applied to any other volcanic areas or over different exposure times, allowing researchers to account for the eruptive history of the target volcanoes that, when available, could include the occurrence of less frequent large eruptions, representing critical elements for risk evaluations.

1 Introduction

Volcanic eruptions produce a large variety of hazards, which widely span spatial distributions and impact the environment and society. The most frequent and widespread one is tephra fallout, occurring in more than 90 % of all eruptions (e.g., Newhall and Hobblit, 2002; Jenkins et al., 2015). Although large explosive eruptions impacting populated areas are relatively infrequent, the tephra fallout can result in considerable direct and indirect impacts on society, which are a function of the hazards; exposure of assets and living beings; and vulnerability of critical infrastructure (e.g., Spence et al., 2005; Menoni et al., 2012; Wilson et al., 2012, 2017; Jenkins et al., 2015). The effects of tephra dispersal and deposition vary according to the size of the eruption, the distance from the vent, the occurrence of aggregation phenomena, and the wind field at the time of the eruption. Lapilli mostly fall from the buoyant plume within tens of kilometers, while ash travels much further, sometimes draping the landscape with ephemeral deposits that can be remobilized long after the eruption due to resuspension by winds and possible generation of volcaniclastic flows (e.g., Blong, 1984; Sulpizio et al., 2006; Alderton and Elias, 2020). Load from tephra deposits may generate structural damages to buildings (e.g., Spence et al., 2005; Costa et al., 2009; Mingari et al., 2022), disrupt infrastructures (roads, rails, power lines, etc.; e.g., Swords-Daniels, 2011; Wilson et al., 2012; Jenkins et al., 2015), and impinge on agriculture and livestock (e.g., Blong, 1984; Wilson et al., 2012; Sulpizio et al., 2014; Loughlin et al., 2015). The injection of large amounts of ash particles into the atmosphere may also alter climatic conditions (e.g., Rampino and Self, 2000) and strongly affect aviation safety (e.g., Casadevall, 1994; Miller and J., 2000; Sulpizio et al., 2012; Biass et al., 2014), leading to heavy economic loss (e.g., Folch and Sulpizio, 2010; Tesche et al., 2012; Wilson et al., 2014, 2017; Bonadonna et al., 2021).

Volcanoes are intrinsically multi-hazard and multi-risk phenomena, with potential significant interaction among the different phenomena and/or damaging mechanisms (Zuccaro et al., 2013; Selva, 2022). Forecasting the type and size of the next eruption from a given volcano can be achieved only statistically through a probabilistic quantification based on past eruptive activity, analogue volcanoes, and other geological information and should account for the natural variability of physicochemical processes and the potential scarcity of data, as well as the limited knowledge of volcanic systems and their processes (e.g., Marzocchi and Bebbington, 2012; Selva et al., 2018). Such an approach sustains the so-called probabilistic volcanic hazard assessment (PVHA) method based on Bayesian procedures for transforming all the volcanological information into probability distributions, allowing aleatory and epistemic uncertainties to be properly accounted for (e.g., Marzocchi et al., 2004, 2007, 2010).

In the case of tephra fallout, the main uncertainties are due to the high variability in eruptive source parameters (e.g., total erupted mass, mass eruption rate, time duration, vent position, total grain-size distribution) and meteorological conditions. Marzocchi et al. (2010) presented a Bayesian event tree workflow (BET_VH) to calculate the probability of any kind of long-term volcanic hazard from a general prior event to subsequent events, weighting each one with its own probability of occurrence. Input incorporates results from numerical simulations aimed at exploring the potential impact of the volcanic phenomena over a selected area, while output is in the form of hazard maps showing the probabilities associated with exceeding critical hazard thresholds at given locations. Such Bayesian event tree (BET) models have largely been used to assess both long- and short-term hazard assessments (e.g., Selva et al., 2010, 2014; Sandri et al., 2012, 2014; Thompson and Lindsay, 2015; Constantinescu et al., 2016), extending the output to Bayesian hazard curves (Tonini et al., 2015).

PVHA is often computed for single volcanoes. On the contrary, in many other fields, the hazard is computed by integrating over all the potential sources, as the interest is not in the source (e.g., probability of an earthquake), but in the occurrence of the dangerous phenomenon (e.g., ground shaking) from whatever source (e.g., Cornell, 1968).

Other geohazards, such as for probabilistic tsunami hazard assessment (PTHA) or probabilistic seismic hazard assessment (PSHA) methodologies, provide a framework for assessing the exceedance probability of a given measure of the intensity of the phenomena (e.g., tsunami wave height, peak ground acceleration) at a particular location within a given time window.

As for volcanic eruptions, historical catalogs are usually incomplete, and therefore a computational scheme is usually adopted, based on the combination of probabilistic source models and empirical or numerical models of propagation of the hazardous phenomena (e.g., Grezio et al., 2017; Gerstenberger et al., 2020). Sometimes, the explicit numerical modeling of individual scenarios is required, like in the case of tsunamis or volcanic eruptions (e.g., Selva et al., 2016; Grezio et al., 2017). In all cases, a common description based on the quantification of the mean annual rates of exceedance is possible, making possible an explicit comparison among the different hazards and consequent risks (e.g., Grunthal and Wahlstrom, 2006; Marzocchi et al., 2012). Besides volcanic activity, other geohazards (such as landslides, meteorological events) could be treated within the same computational framework.

In volcanology, this is complicated by the fact that volcanic systems may be very different to each other and may experience different phases of activity, complicating the potential homogeneous integration of different volcanoes in the same hazard quantification (Selva et al., 2022). To this end, some studies applied methodologies to rank multiple volcanoes according to their hazards (e.g., Aspinall et al., 2011; Auker et al., 2015) or their population exposure (e.g., Small and Naumann, 2001; Freire et al., 2019). Jenkins et al. (2022) quantified and ranked exposure to multiple volcanic hazards for 40 volcanoes in southeast Asia, showing which of these volcanoes may have been overlooked because of not being frequently (or recently) active but with the potential to affect large numbers of people and assets. Other studies accounted for multi-hazards (i.e., lava flows, pyroclastic density currents, fallout, tsunami waves, earthquakes) for a single volcano, typically for volcanic islands (e.g., Becerril et al., 2014; Selva et al., 2019).

The Neapolitan area represents one of the highest volcanic risk areas in the world due to the presence of three active and potentially explosive volcanoes and the extremely high exposure (over 3 million people living in the metropolitan area; https://www.cittametropolitana.na.it/, last access: 16 June 2023). Previous studies assessed tephra fall hazard from Neapolitan volcanoes, combining field data and numerical simulations and calculating conditional ash load probability maps for one or a few reference volcanic scenarios (e.g., Barberi et al., 1995; Cioni et al., 2003; Macedonio et al., 2008; Costa et al., 2009; Folch and Sulpizio, 2010; Sulpizio et al., 2012). This approach was very useful to support civil protection emergency plans (e.g., Emergency Plan of Vesuvius; Dipartimento di Protezione Civile, 2019), and it is computationally less expensive and more feasible for near-real-time applications (e.g., Selva et al., 2014). However, it had limitations in representing the full volcanic hazard because uncertainty only accounted for wind conditions (e.g., Macedonio et al., 2008, 2016) but not for the full range of eruptive style, magnitude, and frequency. Furthermore, a single eruptive scenario is not appropriate to achieve an unbiased PVHA (e.g., Selva et al., 2010; Sandri et al., 2016). For this purpose, the analysis of the intra-scenario variability is certainly more complete, allowing a reduction of the epistemic uncertainty (e.g., Selva et al., 2014, 2018; Sandri et al., 2016; Titos et al., 2022; Martinéz Montesinos et al., 2022).

In recent PVHAs at the Neapolitan volcanoes, the importance of quantifying uncertainty has largely been addressed. For instance, Sandri et al. (2016) proposed a stratified sampling procedure to fully explore the intra-size-class aleatory variability, comparing the results with the classical approach based on reference volcanic scenarios, in the case of proximal and medial areas and large tephra loads from Somma–Vesuvius and Campi Flegrei.

It is worth noting that with the term “eruption size class” we referred to the broad range of possible eruptive sizes identified by the total erupted mass, which is used to define the eruption magnitude. Following Sandri et al. (2016), we consider splitting the eruptive size range into a few classes that can be linked to representative members like the classical approach used in past studies. These classes ideally span the general range continuously, whereas representative members, by definition, discretize it. In this regard, the intra-size-class variability represents the aleatory variability due to combinations of parameters characterizing eruptions which belong to the same eruptive class. Moreover, Selva et al. (2018) provided unbiased tephra fallout hazard estimations at Campi Flegrei by using ensemble modeling of alternative choices related to the treatment of submarine eruptive vents and tephra total grain-size distributions considering different mass fractions of ash and percentages of ash aggregates. Since high-performance computational resources are more and more available, Martinéz Montesinos et al. (2022) proposed a new workflow that was applied to Campi Flegrei, aiming at performing robust and unbiased short- and long-term PVHA on a large-scale (thousands of km) and high-resolution (about 2 km) domain, to be used by civil protection agencies, aviation companies, and other stakeholders. However, all these studies provide output PVHA for single volcanoes and do never merge into a single quantification the contribution of the three active volcanoes in the Neapolitan area: Somma–Vesuvius, Campi Flegrei, and Ischia.

A breakthrough study in providing the estimation of the present state of the Neapolitan volcanoes has been proposed by Selva et al. (2022) through the development of a simple physics-based statistical model that satisfactorily fits the eruptive history of all the Neapolitan volcanoes, accounting for potential changes in their eruptive regime. The model is compatible with existing data (including isolated events and long repose periods) and accounts for two activity regimes (high–low) able to describe the temporal modulations in eruptive activity. Thus, the model can provide a homogeneous quantification of the probability of eruption, which takes into account the state of the volcano and the possible transitions. Such estimations do not account for monitoring data or potential state of unrest, as these types of information have been recorded only in the last few decades, while they are obviously missing along the rest of the eruptive history, but they provide a homogeneous input for long-term volcanic hazard.

Figure 1 Map showing the computational domain (southern Italy) used in the numerical simulations with a magnification of the Neapolitan area. The Somma–Vesuvius, Campi Flegrei, and Ischia volcanoes are indicated as yellow stars. Original map data are available from https://www.openstreetmap.org#map=9/-2.0475/28.5535 (last access: 16 June 2023) OpenStreetMap contributors 2022. Distributed under the Open Data Commons Open Database License (ODbL) v1.0 ©(OpenStreetMap contributors, 2022).

In this work, by merging the model of Selva et al. (2022) into the new PVHA workflow developed by Martinéz Montesinos et al. (2022), we provide a prototypical methodology for long-term tephra fallout hazard assessment on a large-scale domain (southern Italy) at 3 km resolution, combining the impacts of the three active Neapolitan volcanoes. Since the statistical evaluation of the fallout hazard requires us to take into account a wide spectrum of different volcanic scenarios, we explored the intrinsic variability of the explosive eruptions by performing a large number of numerical simulations of tephra dispersion from each volcano, considering the wind patterns over the last 30 years. All simulations were combined to calculate the averaged return period (ARP) of overcoming different tephra load thresholds within a 50-year exposure time. Then, we applied hazard disaggregation (Bazzurro and Cornell, 1999) to quantitatively evaluate the importance of the different volcanoes and different eruptive sizes in the different target areas. This technique, which is widely used in seismic and tsunami hazard analyses in volcanology, can provide important clues about many choices made whenever specific eruptive scenarios were chosen or when high/low priority was given to one or to another volcanic source or size. In this regard, the outcomes of this study add new insights into volcanic risk assessment in southern Italy. In the following, we present the eruptive history of the Neapolitan volcanoes (Sect. 2), the methodology (Sect. 3), and the discussion of the obtained results (Sect. 4), followed by the conclusions (Sect. 5).

2 Eruptive history of Neapolitan volcanoes

2.1 Somma–Vesuvius

The Somma–Vesuvius volcanic complex (Fig. 1) consists of an older volcano, Mt. Somma, dissected by a summit caldera, in which the Vesuvius cone grew after the 472 CE eruption (Santacroce et al., 2008). Four Plinian eruptions associated with caldera collapses repeatedly truncated the Somma volcanic edifice, forming the present-day summit caldera (e.g., Cioni et al., 1999; Santacroce et al., 2008). These are the Pomici di Base (22.03 ± 0.18 cal kyr BP; Santacroce et al., 2008), Mercato (8.89 ± 0.09 cal kyr BP; Mele et al., 2011), Avellino (3.90 ± 0.04 cal kyr BP; Sevink et al., 2011), and Pompeii (79 CE; Sigurdsson and Carey, 1985) eruptions. Sub-Plinian eruptions occurred at 19.1 cal kyr BP (Greenish; Cioni et al., 2003), between Avellino and Pompeii eruptions (APs; Andronico and Cioni, 2002), in 472 CE (Sulpizio et al., 2005) and in 1631 CE (Rosi et al., 1993). The most recent cycle of Vesuvius activity between 1631 CE and 1944 was characterized by recurrent summit and lateral lava effusions associated with semi-persistent and mild explosive activities, interrupted by quiescence periods lasting from months to a maximum of 7 years (Santacroce, 1987; Cioni et al., 2008). During this period, Vesuvius produced a few violent Strombolian eruptions, such as in 1822 CE, 1906 CE, and 1944 CE (e.g., Arrighi et al., 2001; Cole and Scarpati, 2010).

2.2 Campi Flegrei caldera

The Campi Flegrei caldera (Fig. 1) results from at least two main nested collapses related to the Campanian Ignimbrite (CI; 40 ka; e.g., Orsi et al., 1996; De Vivo et al., 2001; Giaccio et al., 2017) and Neapolitan Yellow Tuff (NYT; 15 ka; e.g.,  Deino et al., 2004) eruptions. A magnitude 6.6 eruption (corresponding to a volcanic explosivity index (VEI) = 6) occurred at about 29 ka (Albert et al., 2019). In the last 15 kyr, volcanism occurred within the NYT caldera in three main epochs of activity (Epoch I, II, and III) dated between 15–9.5, 8.6–8.2, and 4.8–3.8 ka, respectively (Orsi et al., 2004). After a long quiescence, the last eruption generated the Monte Nuovo tuff cone in 1538 CE (e.g., Di Vito et al., 1999; Costa et al., 2022). Most of the explosive eruptions were low- or medium-magnitude events with products dispersed over areas ranging from a few to 500 km², while the highest-magnitude Plinian events such as Pomici Principali (12.4 cal kyr BP; Sulpizio et al., 2010) and Agnano-Monte Spina (AMS; 4.5 cal kyr BP; de Vita et al., 1999; Sulpizio et al., 2010) dispersed their products over areas > 1000 km² (Orsi et al., 2009).

2.3 Ischia

Ischia is a volcanic island located in the northwestern part of the Gulf of Naples (Fig. 1; Orsi et al., 1996), characterized by alternating effusive and explosive eruptions (e.g., Selva et al., 2019). Volcanism at Ischia dates back to more than 150 ka and continued with centuries to millennia of quiescence, until the most recent eruption occurred in 1302 CE (e.g., Vezzoli and Barberi, 1988; Sbrana and Toccaceli, 2011; Sbrana et al., 2018). A poorly defined period of pyroclastic activity predated the large Mt. Epomeo Green Tuff caldera-forming eruption (55 ka; Civetta and Gallo, 1991), which was followed by resurgence (Orsi et al., 1991). Volcanism continued up to 33 ka with a series of explosive eruptions of trachytic magmas from the present southwestern and northwestern periphery of the island. After 5 kyr of quiescence volcanism resumed at about 28 ka with the eruption of trachybasaltic magma along the southeastern coast and then continued sporadically until 18 ka. The most recent period of activity began at about 10 ka, but it is mainly concentrated in the past 2.9 kyr, with almost all the vents in the eastern part of the island (Selva et al., 2019). At least 35 effusive and explosive eruptions took place, emplacing lava domes, high-aspect ratio lava flows, tuff cones, tuff rings, and variably dispersed fallout and pyroclastic density current deposits. The most significant recent explosive eruption is the Cretaio eruption that occurred in the 1st century CE(Orsi et al., 1996; Selva et al., 2019).

3 Methodology

3.1 PVHA strategy

Similarly to seismic or tsunami hazard (PSHA or PTHA, respectively), PVHA should, in principle, aggregate the contribution of the different sources, quantifying the exceedance probability or the mean annual frequency of a specific seismic threshold at a specific location in a given exposure time window (e.g., Grezio et al., 2017; Gerstenberger et al., 2020,  and references therein). Translated into the volcanic context, we implemented the PVHA related to a set S of volcanoes by calculating the mean annual frequency with which a certain hazardous phenomenon quantified by the intensity measure z (tephra load at the ground, in our case) exceeds the threshold Z at a geographical location x during a given exposure time window ΔT. Following Martinéz Montesinos et al. (2022) and Selva et al. (2022), we quantify by

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathit{\lambda }\left(z>Z;x,\mathrm{\Delta }T\right)& =\sum _{i\in S}{\mathit{\nu }}_i\sum _{j\in {\text{vent}}_i}\sum _{k\in {\text{size}}_i}P\left({\mathit{\sigma }}_{ik}|E_i\right)\\ & P\left({\mathit{\gamma }}_{ij}|E_i\right)P\left(z>Z|E_i,{\mathit{\sigma }}_{ik},{\mathit{\gamma }}_{ij}\right),\end{array}\end{array}$$

where

• ν_i represents the annual rate of occurrence of eruptions from the volcano i for the exposure time ΔT. The evaluation of this parameter is discussed in Sect. 3.3.
  

• Considering for each volcano a set of possible eruptive size classes, P(σ_ik|E_i) represents the probability that, in the case of eruption E from the volcano i, it will be of size class k (σ_ik). For Somma–Vesuvius, we considered three eruptive (explosive) size classes proposed by Sandri et al. (2016): small (in which the 1906 CE and 1944 CE eruptions are included; Arrighi et al., 2001; Cole and Scarpati, 2010), medium (in which the 472 CE and 1631 CE eruptions are included; Sulpizio et al., 2005; Rosi et al., 1993), and large (in which the Pompei and Avellino eruptions are included; Macedonio et al., 2008; Sulpizio et al., 2010). For Campi Flegrei, the three explosive eruption size classes adopted in Sandri et al. (2016) and Selva et al. (2018) are considered again: small (e.g., Averno 2; Di Vito et al., 2011), medium (e.g., Astroni 6; Mele et al., 2020), and large (e.g., Agnano Monte Spina; de Vita et al., 1999). Ischia has one eruption size class defined as large, representing the Cretaio Tephra (ca. 60 CE; Orsi et al., 1992, 1996), which is the largest explosive eruption that occurred in the last 3 kyr beside 33 other eruptions of lower intensity (de Vita et al., 2010; Selva et al., 2019, 2022; Primerano et al., 2021). We only considered this class because the Cretaio Tephra is the best-characterized eruption in terms of eruptive source parameters (ESPs; Primerano et al., 2021) and represents the only eruptive scenario which can have a significative impact beyond the island (Selva et al., 2019). The parameters related to all eruption size classes are reported in Table A1. The probability of the different size classes corresponds to node 5 of the BET model (Marzocchi et al., 2007, 2010), and it is here evaluated jointly to ν_i:
  

• Considering for each volcano a set of possible eruptive vent locations, P(γ_ij|E_i) represents the probability of the vent activation γ_ij, conditional to the occurrence of the eruption E from the volcano i. For Somma–Vesuvius, being a stratovolcano, a single vent location is fixed on the present-day crater, since the recurrence of the recent eruptions started from there, and the likelihood of different vent positions is an order of magnitude smaller (Sandri et al., 2018). At Ischia, the majority of eruptive vents are preferentially distributed along the NS direction; therefore we fixed a single vent corresponding to the source of the Cretaio Tephra (Primerano et al., 2021), assuming that the spatial probability on vent locations would not affect the tephra ground load at distal areas on land in the target domain (Fig. 1). For Campi Flegrei, given the dimensions of the caldera with respect to its proximal locations in the target domain (Fig. 1), tephra dispersion is simulated from an equally spaced grid of potential vents having a different probability of opening (conditional to the occurrence of an eruption) based on the geological and geomorphological features of the past eruptive vents opened in the last 5000 years (Selva et al., 2012). This probability corresponds to node 4 in the BET methodology, used for example in Selva et al. (2012) and Martinéz Montesinos et al. (2022);
  

• $P(z>Z|E_i,{\mathit{\sigma }}_{ik},\mathit{\gamma }ij)$ represents the probability that a tephra load z exceeds the threshold Z at the target location x due to the eruption E of size class σ_ik from the vent γ_ij from the volcano i, conditional to the occurrence of such an eruption. To calculate this probability, simulations of tephra dispersion and deposition were performed with FALL3D-8.0 (Folch et al., 2020) using as weather conditions those retrieved from the ECMWF ERA5 database (Hersbach et al., 2018) with a spatial resolution of $\mathrm{0.03}{}^{\circ }\times \mathrm{0.03}{}^{\circ }$ and a temporal resolution of 3 h. These data were randomly sampled considering 50 d yr⁻¹ from 1 January 1991 to 31 December 2020 so that the simulation days are homogeneously distributed within the 30-year time span. Adopting the approach proposed in (), each simulation represents a volcanic scenario characterized by a set of eruptive source parameters that are randomly sampled from their probability density functions (Tables A1, A2, A3). For each volcano, 1500 volcanic scenarios were assigned to each eruption size class. In this way, the total number of simulations is 10 500 simulations (4500, Somma–Vesuvius; 4500, Campi Flegrei; 1500, Ischia). This probability corresponds to nodes 7 and 8 in the BET methodology used, for example, in Selva et al. (2012), Tonini et al. (2015), and Martinéz Montesinos et al. (2022).

In this way, the mean annual frequency $\mathit{\lambda }(z;x,\mathrm{\Delta }T)$ to have a tephra load z overcoming a threshold Z at location x during the 50-year exposure time ΔT from the Neapolitan volcanoes is computed. The curve $\mathit{\lambda }(z>Z;x,\mathrm{\Delta }T)$ as a function of the threshold Z represents the final result of PVHA in each target point x, and it is usually referred to as the hazard curve. This is replicated all over the computational domain, which here covers southern Italy. Whenever alternative quantifications of one or more of the parameters in Eq. (1) exist, they jointly represent an estimation of the epistemic uncertainty (see Marzocchi et al., 2021, and references therein).

The specific contribution of each volcano, each eruptive size at one volcano, or whatever other definition of a group of potential eruptions may be quantitatively evaluated using the hazard disaggregation strategy discussed in Bazzurro and Cornell (1999). This method consists of evaluating the probability that the occurrence of a given exceedance (z>Z) is caused by a given group of sources (G), $P\left(G\right|z>Z)$, based on the Bayes rule. Since this evaluation reduces to the quantification of specific ratios of the addends of Eq. (1), it can always be quantified, and thus it nowadays represents a standard post-processing tool for probabilistic hazard, especially seismic hazard.

Figure 2 Workflow of the PVHA strategy used to provide the calculation of the annual frequency of exceeding a specific intensity threshold value at the target site in 50 years associated with tephra fallout from the Neapolitan volcanoes and the relative hazard curves and maps.

3.2 Tephra dispersion modeling

In this study, we used FALL3D-8.0 (Folch et al., 2020), the Eulerian tephra dispersion model based on the so-called advection–diffusion–sedimentation (ADS) equations, for simulating dispersion of volcanic tephra, gas, and radionuclides, with a wide range of possible model parameterization options (e.g., eruptive parameters, source model, ash aggregation, domain discretization), including the possibility to describe the gravitational spreading of the umbrella region (Costa et al., 2013). In the used version (v8.0) the parallelization strategy, input/output (I/O), model pre-process workflows, and memory management have been improved, leading to a better code scalability, efficiency, and an overall capability to handle much larger problems. The outputs are time-dependent loads at the ground and atmospheric ash concentrations.

In this work, simulations were performed on the HPC cluster ADA of the Istituto Nazionale di Geofisica e Vulcanologia, Bologna (Italy). We run one simulation per sampled day in order to cover all possible daily and seasonal meteorological conditions over the last 30 years.

The eruptive source parameters represent the main volcanological inputs in the model, which are described in detail in Tables A1, A2, and A3. In Table A4 we report the summary of the key model parameters which are not dependent on the eruption size class (Appendix A).

The total time spent for each simulation varies according to the duration of the eruption, which depends on the eruption size class, although the model runs for 24 additional hours after the source term is switched off, given the large size of the target domain: this is necessary in order to ensure that most of the remaining airborne material has sufficient time to settle or to leave the computational domain through the lateral boundaries. We considered a spherical projection computational domain from 36.31 to 42.61^∘ N in latitude and from 12.31^∘ W to 19^∘ E in longitude, with a 0.03 ^∘ resolution gridded domain and a vertical σ coordinate system with a linear decay (Fig. 1).

3.3 Construction of hazard curves and maps

In this study, hazard intensity represents tephra accumulated on the ground per unit area, typically expressed in kg m⁻² (e.g., Sandri et al., 2016). As discussed in Sect. 3.1, we quantify the hazard as the mean annual frequency of exceeding a specific intensity threshold value at the target site x in 50 years, due to the activity of the three Neapolitan volcanoes. As discussed in Selva et al. (2022), the probability of the different sizes at each volcano P(σ_ik|E_i) may significantly vary through time, as the regime of the volcanoes changes. For example, at Somma–Vesuvius during open-conduit periods, small sizes are favored, while large explosive eruptions are mainly expected during close-conduit periods. On the contrary, it is possible to quantify the average annual eruption rate of the eruption size class k of the volcano i as

$$\begin{array}{}\text{(2)}& {\mathit{\nu }}_{ik}={\mathit{\nu }}_{i}P\left({\mathit{\sigma }}_{ik}\right|E_i),\end{array}$$

considering the potential switch between regimes and the different size distributions in each region (Selva et al., 2022). To account for this, the contribution of volcano i in Eq. (1) can easily be rewritten as

$$\begin{array}{}\text{(3)}& {\mathit{\lambda }}_i(Z;x)=\sum _{j\in \text{vent}}\sum _{k\in \text{size}}{\mathit{\nu }}_{ik}P\left({\mathit{\gamma }}_{ij}\right|E_i\left)P\right(Z>z|E_i,{\mathit{\sigma }}_{ik},{\mathit{\gamma }}_{ij}).\end{array}$$

The quantification of ν_ik was based on the recorded activities of the Neapolitan volcanoes in the last 2 kyr for Somma–Vesuvius, 15 kyr for Campi Flegrei, and 3 kyr for Ischia, assuming a non-homogeneous Poisson model because volcanoes randomly oscillate through time between two discrete regimes with different size distributions (Selva et al., 2022). For short ΔT (e.g., shorter than 100 years), ν_ik varies depending on the initial state of the volcano and the length of ΔT, due to the persistence of the volcano in its initial state. Indeed, during short ΔT the state of the volcanic activity is not likely to change; on the contrary, the longer the ΔT, the larger is the probability to change the regime. In Table C1, we report the values of ν_ik for Somma–Vesuvius, Campi Flegrei, and Ischia in 50 years. The rates ν_ik intrinsically account for the aleatory uncertainty in eruption occurrence, changes in the regime, and consequent eruption size distribution (Selva et al., 2022). Alternative models do not exist, so the epistemic uncertainty in these values is not quantified here.

The probabilities P(γ_ij|E_i) and $P(z>Z|E_i,{\mathit{\sigma }}_{ik},{\mathit{\gamma }}_{ij})$ were calculated using the BET workflow at node 4 (vent opening) and nodes 7 and 8 (tephra reaching point z – tephra load overcoming threshold Z), respectively. The BET model formally accounts for potential epistemic uncertainty in the probability at each node. However, this may be severely underestimated whenever alternative approaches do produce significantly different results (Marzocchi et al., 2021). For this reason, the epistemic uncertainty is neglected, and only the average probabilities are evaluated.

Simulations were post-processed for each eruption size class in order to quantify the mean annual frequencies of exceeding a given accumulation of tephra on the computational domain (see Fig. 1) in 50 years. The relative hazard maps were built by assessing a mean hazard curve in each point, considering 16 intensity thresholds from 1 to 1600 kg m⁻², which roughly correspond to thicknesses from 1 mm to 1.6 m (considering a typical deposit density of 1600 kg m⁻³).

Generally, values of ∼1 and $\sim \mathrm{10}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ are referred to as non-conservative and conservative bounds for airport disruption (e.g., Folch and Sulpizio, 2010; Sulpizio et al., 2012). In the Neapolitan area, the typical thresholds of 200, 300, and 400 kg m⁻² are used for roof collapses (e.g., Pareschi et al., 2000; Orsi et al., 2004; Spence et al., 2005). According to these references, here we considered two tephra load thresholds: 10 kg m⁻² (∼1 cm) and 300 kg m⁻² (∼30 cm). In Fig. 2 we summarize the workflow of the PVHA strategy used in this study.

3.4 Hazard disaggregation

From a multi-volcano hazard assessment perspective, establishing which volcano contributes the most to the hazard at a given location, assuming that a threshold is overcome, is crucial; equally relevant is assessing which eruption size class provides the largest contribution exceeding a given threshold. These issues can be addressed throughout hazard disaggregation analysis as detailed in the following.

The hazard disaggregation scheme permits us to post-process hazard results to display the relative contributions of the different source to the hazard. The contribution of each individual source depends on many factors, like the size and the position of the event, the annual rates of each specific size and position, and the propagation. These are all ingredients of the hazard, but their combination and their balance are not trivial at all. For example, is the hazard impact higher in one unlikely big event or a more likely smaller event with favorable wind? Is the hazard impact higher in one unlikely event in a location usually upwind or a more likely event in a location usually downwind but with a rare variable wind direction? To answer these questions, the only solution is to post-process the hazard combination and balance the different contribution.

3.4.1 Hazard disaggregation by volcanic source

The hazard disaggregation was originally developed to analyze the importance of the different seismic source regions on the seismic probabilistic hazard assessment (SPHA). In this study, it is used to calculate the contribution to the total hazard in each target point of each volcanic source. This intrinsically accounts for both the probability of eruption and the probability that each specific eruption propagates from the source to the target point.

At first, the disaggregation is evaluated for individual volcanoes. In this case, we evaluate $P\left(E_i\right|z>Z)$, which is the probability that, given the observation of the exceedance of a given threshold Z, this has been caused by an eruption E from the volcano i. To simplify its interpretation, we name this probability dominance D_i in a given spatial point x.

Following Bazzurro and Cornell (1999), this probability can be calculated by simply evaluating the ratio of the mean annual frequency of overcoming a given threshold Z in x from that volcano (λ_i, see Eq. 3) to the mean annual frequency obtained by the sum of the impacts of the three volcanoes in 50 years, which is

$$\begin{array}{}\text{(4)}& D_i=P\left(E_i\right|z>Z)={\displaystyle \frac{{\mathit{\lambda }}_i(z>Z,x)}{{\sum }_i{\mathit{\lambda }}_i(z>Z,x)}}.\end{array}$$

To be meaningful, this calculation must be carried out upon a statistically significant number of simulations overcoming such a tephra load threshold; therefore we focus our attention only on a sub-domain, neglecting the distal areas that are reached by a number of simulations less than 10^1.5 (ca. 30 simulations).

3.4.2 Hazard disaggregation by eruption size class

This kind of disaggregation is performed by calculating, for Somma–Vesuvius and Campi Flegrei only (since Ischia has only one explosive size), the hazard contribution of the different eruption size classes considered. In this case, we first quantify the probability $P\left({\mathit{\sigma }}_{ik}\right|z>Z;E_i)$ of overcoming the tephra threshold Z due to the occurrence of an eruption E from the volcano i having the eruption size class σ_ik. As in Eq. (4), we express this disaggregation in terms of relative dominance of eruption size class D_ik, evaluated as

$$\begin{array}{}\text{(5)}& D_{ik}=P\left({\mathit{\sigma }}_{ik}\right|z>Z;E_i)={\displaystyle \frac{{\mathit{\lambda }}_{ik}(z>Z;x)}{{\sum }_k{\mathit{\lambda }}_{ik}(z>Z,x)}}.\end{array}$$

Notably, the summation in the dividend is not extended to i, meaning that each volcano is treated separately.

4 Results and discussion

4.1 Long-term hazard maps

In this section, we report the results of the PVHA and different levels of disaggregation, from the total hazard assessment of the volcanic sources to the evaluation of the impact due to each volcano obtained by aggregating their eruption size classes.

Generally, hazard curves may be expressed in terms of ARP, which is the inverse of the mean annual frequency. For the sake of simplicity, the most likely values of the ground tephra load occur frequently and may be associated with relatively short ARPs (e.g., 100 years), while a larger Z is exceeded more rarely, corresponding to longer ARPs.

In Fig. 3 the long-term mean hazard maps considering each volcano individually (Fig. 3a–c) and the total mean hazard maps obtained by aggregating the effect of each volcano (Fig. 3d) (see Eq. 4) are shown. Maps are obtained by cutting hazard curves at three different ARPs of 100, 500, and 1000 years, respectively, and reporting in each point of the domain the corresponding tephra load.

At Somma–Vesuvius (Fig. 3a) we show that, for an ARP of 100 years, tephra loads range from ∼100 to 300 kg m⁻² in proximity to the volcano, decreasing to 1–10 kg m⁻² within ∼250 km. At an ARP of 500 years, higher loads are obviously expected (300–1000 kg m⁻² within ∼100 km); at larger distances the maximum tephra load expected with ARP=500 years is 100 kg m⁻². For the remaining part of the domain, the expected tephra load with this ARP is ∼ 1–10 kg m⁻². For an ARP of 1000 years, the boundaries of tephra load isolines are larger than the previous ones, as expected, covering almost the entire computation domain.

At Campi Flegrei (Fig. 3b), the expected tephra load with an ARP of 100 years ranges from 1 to 100 kg m⁻² only in the proximity of the caldera. This is partly due to the uncertainty in the vent position that “blurs” the resulting hazard maps, although offering a more realistic degree of knowledge on future eruptions (Sandri et al., 2016; Massaro et al., 2022). Higher values of tephra loads (from 1 to 500 kg m⁻²) are expected when ARP is 500 years, decreasing with distances larger than ∼ 80 km. For an ARP of 1000 years, the expected tephra load is higher in the proximity of Campi Flegrei, reaching ∼ 1–10 kg m⁻² within ∼ 50 km.

For Ischia, Fig. 3c shows that non-negligible tephra load (between 10 and 300 kg m⁻²) is only expected on the island when ARP is 100 years; for an ARP of 500 and 1000 years, the highest tephra loads on the Italian mainland are expected around Ischia and the Campi Flegrei area, decreasing up to a range of 1–10 kg m⁻² for distances between ∼ 100–170 km.

The total long-term mean hazard maps (Fig. 3d) are obtained by summing the three individual hazard curves (Eq. 1) and cutting the resulting total hazard curve. This corresponds to assuming the independence between different volcanoes. For an ARP of 100 years, areas surrounding the volcanic sources are exposed to tephra loads from 300 to 1000 kg m⁻²; loads of ∼ 200–100 kg m⁻² are expected within ∼ 150 km, decreasing up to 1–10 kg m⁻² at ∼ 250 km. For ARPs of 500 and 1000 years, the boundaries of the tephra load isolines become larger with respect to the previous ones, as expected.

We stress that the fringed edges of the distal tephra load are due to a numerical artifact depending on the low number of simulations reaching the targets on the computational domain (Fig. 3).

Somma–Vesuvius shows the largest expected tephra load for all ARPs because its mean annual rate of eruptions is higher by 1 order of magnitude than those of Campi Flegrei and Ischia (Selva et al., 2022; Table C1). Indeed, Campi Flegrei and Ischia are responsible for an expected tephra load $Z=\mathrm{10}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ only in the proximity of their vents for longer ARPs (500–1000 years). We also remark that the long-term hazard maps show the aleatory uncertainty since, for each volcano, we used the mean annual rate ν_ik corresponding to each eruption size class for a time exposure of 50 years (Table C1).

Figure 3 Long-term hazard maps reporting the mean hazard intensity (tephra load, kg m⁻²) in 50 years as a function of the ARP (100, 500, and 1000 years) due to the impact of (a) Somma–Vesuvius (aggregated eruption sizes: small, medium, large), (b) Campi Flegrei (aggregated eruption sizes: small, medium, large), and (c) Ischia (single eruption size: large). (d) Long-term hazard maps in which the combined impact of the three volcanic sources is taken into account. For each intensity threshold (color bar label), the maps show the areas where that value is expected to be exceeded after 24 h from eruption onset with that ARP.

Figure 4 Disaggregation per volcanic source. (a) Map showing, in each point, which volcano contributes the most to the hazard of exceeding 300 kg m⁻² due to the occurrence of an eruption in 50 years. In panels (b), (c), and (d), the relative dominance of exceeding the threshold due to the occurrence of an eruption from Somma–Vesuvius, Campi Flegrei, and Ischia, respectively, is shown. The dotted black isolines encompass the areas where the labeled number of simulations (in log ₁₀ scale) producing a tephra load greater than 300 kg m⁻² is exceeded.

Figure 5 Disaggregation per eruption size class. (a) Map showing, in each point, which eruption size class of Somma–Vesuvius contributes the most to the hazard of overcoming 300 kg m⁻² due to the occurrence of eruption in 50 years. The relative dominance of overcoming the threshold is accounted for in (b) small, (c) medium, and (d) large size classes. Dashed black contours represent the isolines of the log ₁₀ number of simulations that produce a tephra ground load exceeding 300 kg m⁻².

Figure 6 Disaggregation per eruption size class. (a) Map showing, in each point, which eruption size class of Campi Flegrei contributes the most to the hazard of exceeding 300 kg m⁻² due to the occurrence of an eruption in 50 years. The relative dominance of exceeding the threshold is accounted for in (b) small, (c) medium, and (d) large size classes. Dashed black contours represent the isolines of the log ₁₀ number of simulations that produce a tephra ground load exceeding 300 kg m⁻².

4.2 Hazard disaggregation analysis

4.2.1 Volcanic source disaggregation

The disaggregated contribution of each volcanic source in overcoming the threshold of 300 kg m⁻² in 50 years is shown in Fig. 4a. Somma–Vesuvius (green area) dominates the largest part of the sub-domain considered, up to distances of ∼90 km (green area; Fig. 4a), while Campi Flegrei and Ischia show their dominance within a radius of a few kilometers from their respective vents (red and blue areas; Fig. 4a). It is worth noting that this considers both the different probabilities of eruptions and the different probabilities of dispersions for the various eruptive sizes. Therefore, these results show the much higher probability of eruption at Somma–Vesuvius led to a general dominance of this volcano, which is overcome by the other volcanoes only upwind (prevalent winds here flow toward southeast), where dispersion from Somma–Vesuvius is less likely. This tendency is also highlighted by mapping, for each volcano separately, the actual value of D_i for overcoming the threshold of 300 kg m⁻² (Fig. 4b, c, d). We also remark that the dotted black isolines represented in this figure and the following figures encompass the areas in which different numbers of simulations (in log ₁₀ scale) produce a tephra load greater than 300 kg m⁻²: the larger this number, the more constrained this probability is, since it is evaluated on a larger number of simulations.

4.2.2 Size class disaggregation considering the tephra thresholds of 300 and 10 kg m⁻²

Figure 5a shows the largest relative dominance among the three eruption size classes considered for Somma–Vesuvius and their contributions separately, while Fig. 5b–d show, respectively, the dominance in overcoming the tephra threshold of 300 kg m⁻² in each target point of the domain in 50 years. As for the previous analysis, for such a threshold, the significant number of simulations is restricted to distances up to ∼100 km from the volcanic sources. We note the small size class is dominant in the proximity of the vent, while the medium dominates the adjacent area up to distances of ∼20 km. Distal areas are dominated by the large size class up to distances of ∼100 km (Fig. 5a). In Fig. 6 we show the same information for Campi Flegrei. The large size class homogeneously dominates up to distances of ∼100 km (Fig. 6a). We note that the relative dominance D_i1 of the small class is not visible because it is ubiquitously lower than 0.1 (Fig. 6b). Although the medium class provides a contribution from 10 % to 40 % in exceeding the selected threshold of 300 kg m⁻² in the proximity of the caldera (Fig. 6c), this is nowhere higher than that shown by the large class (>50 %; Fig. 6d). It is important to stress that a variability in the results of this disaggregation is expected if different thresholds are selected. In this regard, in Fig. 6 we report an example of hazard disaggregation per eruption size class for Campi Flegrei in 50 years, considering a tephra threshold of 10 kg m⁻². In this case, results show that the caldera is dominated by the small and medium classes (up to distances of ∼50 km), while the large class nearly completely covers the rest of the southeastern domain, reaching distances of ∼250 km.

Figure 7 Disaggregation per volcanic source and eruption size class. Map showing, in each point, which eruption size class of Somma–Vesuvius, Campi Flegrei, and Ischia contributes the most to the hazard of overcoming 300 kg m⁻² due to the occurrence of eruption in 50 years. Dashed black contours represent the isolines of the log ₁₀ number of simulations that produce a tephra ground load exceeding 300 kg m⁻².

Figure 8 Disaggregation per volcanic source and eruption size class. Map showing, in each point, which eruption size class of Somma–Vesuvius and Campi Flegrei contributes the most to the hazard of overcoming 10 kg m⁻² due to the occurrence of eruption in 50 years. Dashed black contours represent the isolines of the log ₁₀ number of simulations that produce a tephra ground load exceeding the threshold of 10 kg m⁻².

Figure 9 Map showing, in each point, which eruption size class of Somma–Vesuvius contributes the most to the hazard of overcoming 10 kg m⁻² due to the occurrence of eruption in 50 years. The relative dominance of overcoming the threshold is accounted for in (a) small, (b) medium, and (c) large size classes. Dashed black contours represent the isolines of the log ₁₀ number of simulations that produce a tephra ground load exceeding 10 kg m⁻².

Figure 10 Hazard curves showing the mean annual frequency of exceeding a set of intensity thresholds (tephra load, in kg m⁻²) at three target sites: (a) Pozzuoli, (b) the municipality of Naples, and (c) Torre del Greco due to the impact of Somma–Vesuvius (blue curve) and Campi Flegrei (orange curve) in 50 years. The position of the target sites is expressed in UTM coordinates (UTM zone: 33T).

4.2.3 Visualization of the overall disaggregation

In order to have an overall picture of hazard disaggregation, in Fig. 7 we show the results of the relative dominance per volcanic source and eruption size class for the threshold of 300 kg m⁻². As already observed in Fig. 5a, the eruptive size classes of Somma–Vesuvius (green shaded colors) dominate in the largest part of the domain, from the small to the large class at increasing distances from the vent. Notably, the dark red area corresponding to the large class of Campi Flegrei covers the whole caldera up to distances of ∼30 km, as indicated in Fig. 6a. Interestingly, this area almost coincides with the Yellow Zone for tephra fallout hazard taken as reference by the national Italian Civil Protection Department (https://www.protezionecivile.gov.it/it/, last access: 20 June 2023), which is defined as the area having a probability larger than 5 % to overcome 300 kg m⁻² conditional to the occurrence of an eruption of medium size on land (Dipartimento di Protezione Civile, 2019). Therefore, our results confirm that Campi Flegrei mainly dominates in this area, even if that choice was not based on a quantitative evaluation of the dominance. More than 800 000 inhabitants live in this area, located in Naples and the surrounding municipalities. Intuitively, Ischia Island is dominated by the Ischian events. This occurs up to distances of ∼10 km, showing that only in this area, with Ischia being upwind with respect to other volcanoes, Ischian fallouts dominate even if the volcano produces significantly smaller eruptions. The results drastically change when varying the tephra load thresholds. In Fig. 8 we show the disaggregation by volcanic source and eruptive size class corresponding to the tephra threshold of 10 kg m⁻². Here, the entire computational domain is reached by a significant number of simulations with the exception of the northwestern sector and a large part of Sicily. At this load threshold, Somma–Vesuvius provides the largest dominance (green shaded colors), having the small and medium classes able to dominate up to distances of ∼50 km from the vent and the large class over ∼50 km along the preferential towards the east wind direction. This result is expected for Somma–Vesuvius, which shows the highest mean annual rates in 50 years (Table C1). On the contrary, Campi Flegrei shows a reduced-area largest dominance because the mean annual frequency to overcome the 10 kg m⁻² threshold due to an eruption of the large size class (dark red area) is less than that shown by the medium and large size classes of Somma–Vesuvius. Ischia, modeled with a single eruption size class, dominates only locally. In Fig. 9 we show the relative dominance D_ik of each eruption size class of Somma–Vesuvius in overcoming the tephra threshold of 10 kg m⁻² in 50 years, noting that the small class contributes to this hazard intensity only within distances of ∼60 km from the vent, with a relative dominance of ∼60 % only in the proximal areas (Fig. 9a), while the medium class contributes within distances of ∼100 km (including the Gulf of Naples, the northern part of the peninsula Sorrentina, and Ischia; Fig. 9b). At further distances we observe the large class with radial lobes of relative dominance, likely due to the randomness of the sampling of wind conditions (Fig. 9c). In Fig. 6b–d, we report the same type of results in overcoming the threshold of 300 kg m⁻² for Campi Flegrei.

4.3 Implications for the PVHA in the Neapolitan urban area and for multi-hazard evaluation

The proposed PVHA is firstly presented as a combination of the relative contributions of each volcano in terms of mean annual frequency λ of exceeding different tephra load thresholds at a specific location, considering a reference exposure time of 50 years. Then, disaggregation analyses show which volcano, and further which eruption size class, impacts the most on a hazard, i.e., has the largest probability of causing the exceedance of tephra load thresholds of 10 and 300 kg m⁻². Even if the disaggregation results are relative to the selected tephra load threshold, this result is general, and it is valid for all thresholds. The complete information for all thresholds is reported only in the entire hazard curves. In Fig. 10 we provide the hazard curves, showing the mean annual frequency of exceeding the tephra load thresholds for three sites along the dominant wind direction: Pozzuoli, Naples, and Torre del Greco. As expected, within the Campi Flegrei caldera (e.g., in Pozzuoli, Fig. 10a) the major hazard is due to the Campi Flegrei volcano for all thresholds, also considering that Pozzuoli is located upwind with respect to Somma–Vesuvius. The center of the city of Naples represents a transition case: the hazard curves seem to be overlapped up to tephra loads of ∼ 10–20 kg m⁻², but for higher thresholds (from ∼20 to 1000 kg m⁻²) Campi Flegrei shows a higher hazard curve (Fig. 10b). This is in accordance with Selva et al. (2012), indicating that hazard exposure of Naples due to the Campi Flegrei caldera, even in quiet periods, is higher than for Somma–Vesuvius, given that expected eruption size classes are comparable (Marzocchi et al., 2004; Orsi et al., 2009), and the city center is closer to the eruptive vents of Campi Flegrei and more directly downwind (e.g., Selva et al., 2010). Beyond the 1000 kg m⁻² threshold, the inversion occurs (Somma–Vesuvius dominates), but this trend is not reliable. In the case of Torre del Greco, the highest hazard curve is due to Somma–Vesuvius for all the thresholds (Fig. 10c), being located downwind. The results shown in Figs. 4–8 are due to the superimposition of two effects: the prevalent wind direction and probability of eruption. The probability of eruption is much higher at Somma–Vesuvius than at the other volcanoes (∼34 % in 50 years vs. ∼3 % and ∼12 % in 50 years, respectively; Selva et al. (2022). Wind directions statistically prevail towards the east in the Neapolitan area (e.g., Macedonio et al., 2016) with moderate seasonal variations (e.g., Costa et al., 2009): this aspect affects the model results, making tephra load dispersion oriented (e.g., Figs. 3–8). This means, for example, that the western side of the domain has a lower probability to be reached by a tephra fallout due to Somma–Vesuvius. It is worth noting that some modeling limitations need to be taken into account, since they significantly affect the model results, especially at greater distances from the volcanic sources. Although the number of the simulated volcanic scenarios per eruption size class (1500) represents a good compromise between the quality of the results and computational costs, it restricts the area of reliable results in terms of the number of simulations exceeding the tephra load thresholds at an increasing distance from the vents. However, this equally affects all three volcanoes, even if this effect is more evident for Somma–Vesuvius, which has the highest annual rates. It is also worth noting that the model neglects the epistemic uncertainty that, in the future, could be treated by substituting single estimates with ensembles of alternative estimates (Marzocchi et al., 2021; Selva et al., 2022). Considering this as an intrinsic limitation of the methodology, along with the incomplete geological records, this work provides a prototypal methodology for long-term multi-volcano hazard assessment focused on the full quantification of the natural variability of the modeled phenomena, the aleatory uncertainty. The results obtained in this study produce, for the very first time, a complete PVHA, considering all the volcanoes in order to produce the total unconditional hazard. In this sense, our results significantly improve those from the previous PVHA for tephra fallout conditional to the occurrence of specific representative eruptive scenarios from Neapolitan volcanoes (i.e., Macedonio et al., 2008; Costa et al., 2009), conditional to the occurrence of eruptions at one volcano (e.g., Orsi et al., 2009; Selva et al., 2010, 2018; Sandri et al., 2016). In this sense, our PVHA is, for the first time, fully compatible with other hazard studies like the ones for seismic events (Gerstenberger et al., 2020). The reason is that the natural variability of such phenomena is now fully explored by merging the results of a large number of numerical simulations for three volcanic sources, taking into account the impact of low-probability but high-consequence events. In particular, we estimate the mean annual frequency of exceeding specific tephra load thresholds at selected sites in 50 years (but extendible to any time window), considering the mean annual eruption rate from the model of Selva et al. (2022). Notably, this allows us to also consider that, for each volcano, the sequences of eruptive clusters with different eruptive size distributions may occur, also accounting for isolated eruptions and long-term repose times. It is worth noting that the very large magnitude eruptive events, i.e., for Campi Flegrei, events like the Neapolitan Yellow Tuff (e.g., Lirer and Munno, 1975; Scandone et al., 1991; Orsi et al., 2004) and the Campanian Ignimbrite (e.g., Barberi et al., 1978; Costa et al., 2013), are not taken into account because of their negligible probability of occurrence over an exposure time of 50 years (e.g., Bevilacqua et al., 2016; Orsi et al., 2009). However, they would need to be considered if we selected a larger exposure time window. The proposed methodology also represents a clear step forward into the multi-hazard and risk perspective with respect to the previous PVHAs in response to the need for homogeneous model definitions for an effective comparison among volcanoes and for producing coherent multi-volcano long-term hazards and multiple risk quantifications and risk rankings (i.e., Marzocchi et al., 2015). In doing this, a further challenge is to coordinate the efforts of hazard scientists and decision-makers in order to maximize the PVHA potential benefit for society (Selva, 2022).

5 Conclusions

In this study, we provide a new long-term PVHA of tephra load associated with explosive eruptions from the Somma–Vesuvius, Campi Flegrei, and Ischia volcanoes. By defining a set of different eruption size classes for each volcano (small, medium, and large), we created a synthetic dataset of tephra ground loads composed by a total of 10 500 numerical simulations (1500 for each eruption size class) that consider a meteorological variability over the last 30 years. The hazard evaluation was performed through the HPC workflow recently developed by Martinéz Montesinos et al. (2022), accounting for the uncertainty in the eruptive source parameters, vent opening (only for Campi Flegrei), and the mean annual rate of eruption for each eruption and size class (Selva et al., 2022). In this way, we obtained a hazard model from which we derived a set of hazard maps for southern Italy, showing the threshold tephra load that would be exceeded with selected ARPs (i.e., 100, 500, 1000 years) within a 50-year exposure time. As expected, the south-southeast (SSE) area of the computational domain is mainly threatened by the tephra fallout hazard, especially in the proximity of the Neapolitan area. By performing hazard disaggregation, we quantified the relative dominance of the three volcanic sources and their eruptive size classes, showing that Somma–Vesuvius gives a major contribution to the total tephra load hazard for most of southern Italy, compared to Campi Flegrei and Ischia. This is mainly due to its greater mean eruption rate, which is an order of magnitude larger than the ones of the other two volcanoes. Campi Flegrei dominates in the area of the city center of Naples, in which the lower mean annual rate of Campi Flegrei is compensated by the low probability of eruption dispersing tephra toward the west in the Neapolitan area. Ischia, instead, dominates only locally. From a wider methodological perspective, this study aims to improve the PVHA in areas in proximity of multiple active volcanoes, like Naples, homogeneously cumulating the effect of all existing volcanoes. This produces a robust approach that allows easier comparisons among the different possible eruptive scenarios (e.g., volcano, size, etc.), as well as with those used for seismic phenomena and other natural disasters. Moreover, this approach can be applied to compute tephra fallout hazard in different areas or for different time spans, accounting for less frequent events that can still be significant for hazard assessment.

Appendix A: definition and sampling of the eruptive source parameters

Here we report some information about the definition and sampling of the eruptive source parameters (ESPs) for each volcano.

Table A1 ESPs for Somma–Vesuvius, Campi Flegrei, and Ischia. For each eruption size class, TEM and the mass of the fallout phase (kg), the duration of the fallout phase (h), and MER (kg s⁻¹) are reported.

Table A2 Parameters of the beta distribution best fitting the field-based total grain size distributions (TGSDs) for Campi Flegrei (Mele et al., 2020). In this case, the reported parameters are referred to as the juvenile, lithics, and crystals; p and (1−p) are, respectively, the coarse and fine sub-population weights; μ₁, μ₂ and σ₁, σ₂ are the mean and standard deviations of the two Gaussian distributions in Φ units.

Table A3 Parameters of the beta distribution best fitting the field-based TGSDs for Campi Flegrei (Mele et al., 2020). In this case, the reported parameters are referred to as the juvenile, lithics, and crystals; p and (1−p) are, respectively, the coarse and fine sub-population weights; μ₁, μ₂ and σ₁, σ₂ are the mean and standard deviations of the two Gaussian distributions in Φ units.

Table A4 Some parameters and models used to run numerical simulations in FALL3D v.8.0.

The total erupted mass (TEM) and duration of the fallout phase were randomly sampled from uniform distributions with different ranges for each eruption size class (Sandri et al., 2016; Primerano et al., 2021; Table A1). In this study, we refer to the mass of the fallout phase representing 80 % and 25 % of the total erupted mass (TEM) in the case of Somma–Vesuvius and Campi Flegrei, respectively (Sandri et al., 2016). In the case of Ischia, we assume that the TEM associated with the Cretaio Tephra eruption is almost entirely due to a fallout phase (Primerano et al., 2021).

For Somma–Vesuvius and Campi Flegrei, the mass eruption rate (MER) is assumed to be constant during the eruption, and it was obtained by calculating the ratio between the mass of the fallout phase and its duration. For Ischia, we use the empirical relationship of Mastin et al. (2009), sampling the eruptive column height from a beta distribution (α=2 and β=3, scaled by a factor of 10 km and translated by 5 km such that $\stackrel{\mathrm{‾}}{H}=\mathrm{9}\mathrm{km}$) within the interval 5–15 km (as in Primerano et al., 2021). It is worth noting that in the BET workflow, each simulation is weighted according to a power law which depends on the mass of the fallout phase of each eruption size class (Sandri et al., 2016).

Since the dispersion of tephra is strongly influenced by the geometry of the dispersed particles, an accurate characterization of the whole size range of erupted particles is necessary to assign the associated mass and describe the tephra distribution in the proximal and medial areas around the volcano. In this study, TGSD is randomly sampled by a range of values for each size class by assuming a bi-Gaussian distribution (Eq. A1) for the variable Φ (grain size, considering that $d={\mathrm{2}}^{-\mathrm{\Phi }}$ is the particle diameter in mm):

$$\begin{array}{}\text{(A1)}& f\left(\mathrm{\Phi }\right)=p{\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathrm{1}}\sqrt{\mathrm{2}\mathit{\pi }}}}{e}^{-\frac{(\mathit{\varphi }-{\mathit{\mu }}_{\mathrm{1}}{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}+(\mathrm{1}-p){\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathrm{2}}\sqrt{\mathrm{2}\mathit{\pi }}}}{e}^{-\frac{(\mathit{\varphi }-{\mathit{\mu }}_{\mathrm{2}}{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}},\end{array}$$

where p and 1−p are the fine and coarse sub-population weights, and μ₁, σ₁ and μ₂, σ₂ are the mean and standard deviation of the two Gaussians, respectively. The five parameters of the distribution (p, μ₁, σ₁, μ₂, σ₂) are defined for each eruption size class according to what is reported in the literature (Poret et al., 2020; Mele et al., 2020; Primerano et al., 2021) and randomly sampled from beta distributions (Tables A2, A3). In the case of Ischia, TGSD is estimated following Costa et al. (2016).

For Somma–Vesuvius and Ischia, we consider a bulk class of particles; for Campi Flegrei, the availability of literature data (Mele et al., 2020) allowed Martinéz Montesinos et al. (2022) to consider the representative fractions for juvenile, lithic, and crystal components. Juvenile particles are modeled using a bi-Gaussian distribution, while lithics and crystals through a Gaussian distribution (Mele et al., 2020).

Tephra particle density is set as follows: from 900 to 2500 kg m⁻³ for Somma–Vesuvius (see Sandri et al., 2016, and references therein) and from 740 to 2460 kg m⁻³ for Ischia (Primerano et al., 2021). For Campi Flegrei, we distinguished the density of lithics (2500 kg m⁻³) and crystals (2800 kg m⁻³) as in Martinéz Montesinos et al. (2022), considering the estimates provided by Mele et al. (2020).

In Table A4 we report the models and other physical parameters used to run the simulations. For each volcano, the emission source term is modeled considering the Suzuki option (Suzuki, 1983; Pfeiffer et al., 2005), which assumes a mushroom-like vertical distribution of emission points depending on two dimensionless parameters λ and A. The parameter λ, introduced by Pfeiffer et al. (2005), controls the distribution of the emitted mass around the maximum, while A controls the vertical location of the maximum of the emission profile. For our applications, λ is set to 1 and A randomly sampled in the range [3,5].

The aggregation of tephra particles affects the sedimentation dynamics and deposition (e.g., Durant et al., 2009; Folch and Sulpizio, 2010). FALL3D, besides the option based on the model of Costa et al. (2013), includes some simple a priori aggregation parameterizations consisting of empirically based predefined fractions of aggregating classes being transferred to one or more classes of aggregates within the eruptive column (i.e., aggregation is performed before transport). We use the Cornell aggregation model (Cornell et al., 1983), selecting a single tephra bin for Somma–Vesuvius (diameter of 200 µm, density sampled within a uniform distribution, 100–600 kg m⁻³) and Ischia (diameter of 200 µm, density sampled within a uniform distribution, 50–500 kg m⁻³). For Campi Flegrei, we use the percentage model (Sulpizio et al., 2012), selecting two classes of aggregates for each eruption size class. For accretionary lapilli the diameter is set to 2000 µm, and density is sampled in the range 1000–2000 kg m⁻³. For other aggregates, the diameter is set to 200 µm, and density varies in the range 100–600 kg m⁻³ (as in Martinéz Montesinos et al., 2022).

Appendix B: disaggregation per eruption size class for Campi Flegrei

Figure 6 reports the size class disaggregation analysis (see Eq. 5) for Campi Flegrei with a ground load threshold of 300 kg m⁻². Here we present the same information for the threshold of 10 kg m⁻².

Figure B1 Disaggregation per eruption size class. (a) Map showing, in each point, which eruption size class of Campi Flegrei contributes the most to the hazard of overcoming 10 kg m⁻² due to the occurrence of eruption in 50 years. D_ik of overcoming the threshold is accounted for by (b) small, (b) medium, and (d) large size classes. Dashed black contours represent the isolines of the log ₁₀ number of simulations that produce a tephra ground load exceeding the threshold of 10 kg m⁻².

Appendix C: mean annual eruption rates

The mean annual rates presented in this appendix have been used in Eq. (3) and have been obtained with the method presented in Selva et al. (2022).

Table C1 Mean annual rates (ν_ik) calculated for each eruption size class k of the Somma–Vesuvius, Campi Flegrei, and Ischia volcanoes in 50 years.

Table C2 Mean annual rates (ν_ik) calculated for each eruption size class k of the Somma–Vesuvius, Campi Flegrei, and Ischia volcanoes in 5 years.

Table C3 Mean annual rates (ν_ik) calculated for each eruption size class k of the Somma–Vesuvius, Campi Flegrei, and Ischia volcanoes in 500 years.