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

Abstract

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.,. 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.,. 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.,. Load from tephra deposits may generate structural damages to buildings e.g.,, disrupt infrastructures roads, rails, power lines, etc.; e.g.,, and impinge on agriculture and livestock e.g.,. The injection of large amounts of ash particles into the atmosphere may also alter climatic conditions e.g., and strongly affect aviation safety e.g.,, leading to heavy economic loss e.g.,.

Volcanoes are intrinsically multi-hazard and multi-risk phenomena, with potential significant interaction among the different phenomena and/or damaging mechanisms . 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.,. 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.,.

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. 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.,, extending the output to Bayesian hazard curves .

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.,.

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.,. Sometimes, the explicit numerical modeling of individual scenarios is required, like in the case of tsunamis or volcanic eruptions e.g.,. 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.,. 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 . To this end, some studies applied methodologies to rank multiple volcanoes according to their hazards e.g., or their population exposure e.g.,. 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.,.

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.,. This approach was very useful to support civil protection emergency plans e.g., Emergency Plan of Vesuvius;, and it is computationally less expensive and more feasible for near-real-time applications e.g.,. However, it had limitations in representing the full volcanic hazard because uncertainty only accounted for wind conditions e.g., 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.,. For this purpose, the analysis of the intra-scenario variability is certainly more complete, allowing a reduction of the epistemic uncertainty e.g.,.

In recent PVHAs at the Neapolitan volcanoes, the importance of quantifying uncertainty has largely been addressed. For instance, 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 , 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, 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, 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 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 ©.

[Figure omitted. See PDF]

In this work, by merging the model of into the new PVHA workflow developed by , 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 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. ), the methodology (Sect. ), and the discussion of the obtained results (Sect. ), followed by the conclusions (Sect. ).

2 Eruptive history of Neapolitan volcanoes

2.1 Somma–Vesuvius

The Somma–Vesuvius volcanic complex (Fig. ) consists of an older volcano, Mt. Somma, dissected by a summit caldera, in which the Vesuvius cone grew after the 472 CE eruption . Four Plinian eruptions associated with caldera collapses repeatedly truncated the Somma volcanic edifice, forming the present-day summit caldera e.g.,. These are the Pomici di Base 22.03 $\pm$ 0.18 cal kyr BP;, Mercato 8.89 $\pm$ 0.09 cal kyr BP;, Avellino 3.90 $\pm$ 0.04 cal kyr BP;, and Pompeii 79 CE; eruptions. Sub-Plinian eruptions occurred at 19.1 cal kyr BP Greenish;, between Avellino and Pompeii eruptions APs;, in 472 CE and in 1631 CE . 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 . During this period, Vesuvius produced a few violent Strombolian eruptions, such as in 1822 CE, 1906 CE, and 1944 CE e.g.,.

2.2 Campi Flegrei caldera

The Campi Flegrei caldera (Fig. ) results from at least two main nested collapses related to the Campanian Ignimbrite CI; 40 ka; e.g., and Neapolitan Yellow Tuff NYT; 15 ka; e.g., eruptions. A magnitude 6.6 eruption (corresponding to a volcanic explosivity index (VEI) $=$ 6) occurred at about 29 ka . 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 . After a long quiescence, the last eruption generated the Monte Nuovo tuff cone in 1538 CE e.g.,. Most of the explosive eruptions were low- or medium-magnitude events with products dispersed over areas ranging from a few to 500 km${}^{\mathrm{2}}$, while the highest-magnitude Plinian events such as Pomici Principali 12.4 cal kyr BP; and Agnano-Monte Spina AMS; 4.5 cal kyr BP; dispersed their products over areas $>$ 1000 km${}^{\mathrm{2}}$ .

2.3 Ischia

Ischia is a volcanic island located in the northwestern part of the Gulf of Naples Fig. ;, characterized by alternating effusive and explosive eruptions e.g.,. 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.,. A poorly defined period of pyroclastic activity predated the large Mt. Epomeo Green Tuff caldera-forming eruption 55 ka;, which was followed by resurgence . 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 . 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.

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., 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 $\mathrm{\Delta }T$. Following and , we quantify by

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}$$ where

• ${\mathit{\nu }}_i$ represents the annual rate of occurrence of eruptions from the volcano $i$ for the exposure time $\mathrm{\Delta }T$. The evaluation of this parameter is discussed in Sect. .

• Considering for each volcano a set of possible eruptive size classes, $P\left({\mathit{\sigma }}_{ik}\right|E_i)$ represents the probability that, in the case of eruption $E$ from the volcano $i$, it will be of size class $k$ (${\mathit{\sigma }}_{ik}$). For Somma–Vesuvius, we considered three eruptive (explosive) size classes proposed by : small in which the 1906 CE and 1944 CE eruptions are included;, medium in which the 472 CE and 1631 CE eruptions are included;, and large in which the Pompei and Avellino eruptions are included;. For Campi Flegrei, the three explosive eruption size classes adopted in and are considered again: small e.g., Averno 2;, medium e.g., Astroni 6;, and large e.g., Agnano Monte Spina;. Ischia has one eruption size class defined as large, representing the Cretaio Tephra ca. 60 CE;, which is the largest explosive eruption that occurred in the last 3 kyr beside 33 other eruptions of lower intensity . We only considered this class because the Cretaio Tephra is the best-characterized eruption in terms of eruptive source parameters ESPs; and represents the only eruptive scenario which can have a significative impact beyond the island . The parameters related to all eruption size classes are reported in Table . The probability of the different size classes corresponds to node 5 of the BET model , and it is here evaluated jointly to ${\mathit{\nu }}_i$:

• Considering for each volcano a set of possible eruptive vent locations, $P\left({\mathit{\gamma }}_{ij}\right|E_i)$ represents the probability of the vent activation ${\mathit{\gamma }}_{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 . 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 , 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. ). For Campi Flegrei, given the dimensions of the caldera with respect to its proximal locations in the target domain (Fig. ), 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 . This probability corresponds to node 4 in the BET methodology, used for example in and ;

• $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 ${\mathit{\sigma }}_{ik}$ from the vent ${\mathit{\gamma }}_{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 using as weather conditions those retrieved from the ECMWF ERA5 database 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${}^{-\mathrm{1}}$ 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 , , ). 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 , , and .

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 $\mathrm{\Delta }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. () exist, they jointly represent an estimation of the epistemic uncertainty seeand 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 . 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. (), 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.

[Figure omitted. See PDF]

3.2 Tephra dispersion modeling

In this study, we used FALL3D-8.0 , 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 . 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 , , and . In Table we report the summary of the key model parameters which are not dependent on the eruption size class (Appendix ).

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${}^{\circ }$ N in latitude and from 12.31${}^{\circ }$ W to 19${}^{\circ }$ E in longitude, with a 0.03 ${}^{\circ }$ resolution gridded domain and a vertical $\mathit{\sigma }$ coordinate system with a linear decay (Fig. ).

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 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ e.g.,. As discussed in Sect. , 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 , the probability of the different sizes at each volcano $P\left({\mathit{\sigma }}_{ik}\right|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

2$${\mathit{\nu }}_{ik}={\mathit{\nu }}_{i}P\left({\mathit{\sigma }}_{ik}\right|E_i),$$ considering the potential switch between regimes and the different size distributions in each region . To account for this, the contribution of volcano i in Eq. () can easily be rewritten as 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}).$$ The quantification of ${\mathit{\nu }}_{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 . For short $\mathrm{\Delta }T$ (e.g., shorter than 100 years), ${\mathit{\nu }}_{ik}$ varies depending on the initial state of the volcano and the length of $\mathrm{\Delta }T$, due to the persistence of the volcano in its initial state. Indeed, during short $\mathrm{\Delta }T$ the state of the volcanic activity is not likely to change; on the contrary, the longer the $\mathrm{\Delta }T$, the larger is the probability to change the regime. In Table , we report the values of ${\mathit{\nu }}_{ik}$ for Somma–Vesuvius, Campi Flegrei, and Ischia in 50 years. The rates ${\mathit{\nu }}_{ik}$ intrinsically account for the aleatory uncertainty in eruption occurrence, changes in the regime, and consequent eruption size distribution . Alternative models do not exist, so the epistemic uncertainty in these values is not quantified here.

The probabilities $P\left({\mathit{\gamma }}_{ij}\right|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 . 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. ) 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 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$, which roughly correspond to thicknesses from 1 mm to 1.6 m (considering a typical deposit density of 1600 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$).

Generally, values of $\sim \mathrm{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.,. In the Neapolitan area, the typical thresholds of 200, 300, and $\mathrm{400}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ are used for roof collapses e.g.,. According to these references, here we considered two tephra load thresholds: $\mathrm{10}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ ($\sim \mathrm{1}$ cm) and $\mathrm{300}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ ($\sim \mathrm{30}\mathrm{cm}$). In Fig.  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 , 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 (${\mathit{\lambda }}_i$, see Eq. ) to the mean annual frequency obtained by the sum of the impacts of the three volcanoes in 50 years, which is

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)}}.$$ 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 ${\mathrm{10}}^{\mathrm{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 ${\mathit{\sigma }}_{ik}$. As in Eq. (), we express this disaggregation in terms of relative dominance of eruption size class $D_{ik}$, evaluated as

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)}}.$$ 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.  the long-term mean hazard maps considering each volcano individually (Fig. a–c) and the total mean hazard maps obtained by aggregating the effect of each volcano (Fig. d) (see Eq. ) 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. a) we show that, for an ARP of 100 years, tephra loads range from $\sim \mathrm{100}$ to $\mathrm{300}\text{kg}{\text{m}}^{-\mathrm{2}}$ in proximity to the volcano, decreasing to 1–$\mathrm{10}\text{kg}{\text{m}}^{-\mathrm{2}}$ within $\sim \mathrm{250}\text{km}$. At an ARP of 500 years, higher loads are obviously expected (300–$\mathrm{1000}\text{kg}{\text{m}}^{-\mathrm{2}}$ within $\sim \mathrm{100}\text{km}$); at larger distances the maximum tephra load expected with $\text{ARP}=\mathrm{500}$ years is $\mathrm{100}\text{kg}{\text{m}}^{-\mathrm{2}}$. For the remaining part of the domain, the expected tephra load with this ARP is $\sim$ 1–10$\text{kg}{\text{m}}^{-\mathrm{2}}$. 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. b), the expected tephra load with an ARP of 100 years ranges from 1 to 100 kg m${}^{-\mathrm{2}}$ 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 . Higher values of tephra loads (from 1 to 500 kg m${}^{-\mathrm{2}}$) are expected when ARP is 500 years, decreasing with distances larger than $\sim$ 80 km. For an ARP of 1000 years, the expected tephra load is higher in the proximity of Campi Flegrei, reaching $\sim$ 1–10 kg m${}^{-\mathrm{2}}$ within $\sim$ 50 km.

For Ischia, Fig. c shows that non-negligible tephra load (between 10 and 300 kg m${}^{-\mathrm{2}}$) 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${}^{-\mathrm{2}}$ for distances between $\sim$ 100–170 km.

The total long-term mean hazard maps (Fig. d) are obtained by summing the three individual hazard curves (Eq. ) 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${}^{-\mathrm{2}}$; loads of $\sim$ 200–100 kg m${}^{-\mathrm{2}}$ are expected within $\sim$ 150 km, decreasing up to 1–10 kg m${}^{-\mathrm{2}}$ at $\sim$ 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. ).

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 (; Table ). 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 ${\mathit{\nu }}_{ik}$ corresponding to each eruption size class for a time exposure of 50 years (Table ).

Figure 3

Long-term hazard maps reporting the mean hazard intensity (tephra load, $\text{kg}{\text{m}}^{-\mathrm{2}}$) 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 omitted. See PDF]

Figure 4

Disaggregation per volcanic source. (a) Map showing, in each point, which volcano contributes the most to the hazard of exceeding $\mathrm{300}\text{kg}{\text{m}}^{-\mathrm{2}}$ 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 }_{\mathrm{10}}$ scale) producing a tephra load greater than $\mathrm{300}\text{kg}{\text{m}}^{-\mathrm{2}}$ is exceeded.

[Figure omitted. See PDF]

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 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ 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 }_{\mathrm{10}}$ number of simulations that produce a tephra ground load exceeding 300 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$.

[Figure omitted. See PDF]

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${}^{-\mathrm{2}}$ 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 }_{\mathrm{10}}$ number of simulations that produce a tephra ground load exceeding 300 kg m${}^{-\mathrm{2}}$.

[Figure omitted. See PDF]

4.2 Hazard disaggregation analysis

4.2.1 Volcanic source disaggregation

The disaggregated contribution of each volcanic source in overcoming the threshold of 300 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ in 50 years is shown in Fig. a. Somma–Vesuvius (green area) dominates the largest part of the sub-domain considered, up to distances of $\sim \mathrm{90}\mathrm{km}$ (green area; Fig. a), while Campi Flegrei and Ischia show their dominance within a radius of a few kilometers from their respective vents (red and blue areas; Fig. a). 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 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ (Fig. b, 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 }_{\mathrm{10}}$ scale) produce a tephra load greater than $\mathrm{300}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$: 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 $\mathrm{300}$ and $\mathrm{10}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$

Figure a shows the largest relative dominance among the three eruption size classes considered for Somma–Vesuvius and their contributions separately, while Fig. b–d show, respectively, the dominance in overcoming the tephra threshold of 300 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ 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 $\sim \mathrm{100}\mathrm{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 $\sim \mathrm{20}\mathrm{km}$. Distal areas are dominated by the large size class up to distances of $\sim \mathrm{100}\mathrm{km}$ (Fig. a). In Fig.  we show the same information for Campi Flegrei. The large size class homogeneously dominates up to distances of $\sim \mathrm{100}\mathrm{km}$ (Fig. a). We note that the relative dominance $D_{i\mathrm{1}}$ of the small class is not visible because it is ubiquitously lower than $\mathrm{0.1}$ (Fig. b). Although the medium class provides a contribution from 10 % to 40 % in exceeding the selected threshold of 300 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ in the proximity of the caldera (Fig. c), this is nowhere higher than that shown by the large class ($>\mathrm{50}$ %; Fig. d). 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.  we report an example of hazard disaggregation per eruption size class for Campi Flegrei in 50 years, considering a tephra threshold of 10 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$. In this case, results show that the caldera is dominated by the small and medium classes (up to distances of $\sim \mathrm{50}\mathrm{km}$), while the large class nearly completely covers the rest of the southeastern domain, reaching distances of $\sim \mathrm{250}\mathrm{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 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ due to the occurrence of eruption in 50 years. Dashed black contours represent the isolines of the ${\log }_{\mathrm{10}}$ number of simulations that produce a tephra ground load exceeding 300 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$.

[Figure omitted. See PDF]

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 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ due to the occurrence of eruption in 50 years. Dashed black contours represent the isolines of the ${\log }_{\mathrm{10}}$ number of simulations that produce a tephra ground load exceeding the threshold of 10 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$.

[Figure omitted. See PDF]

Figure 9

Map showing, in each point, which eruption size class of Somma–Vesuvius contributes the most to the hazard of overcoming 10 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ 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 }_{\mathrm{10}}$ number of simulations that produce a tephra ground load exceeding 10 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$.

[Figure omitted. See PDF]

Figure 10

Hazard curves showing the mean annual frequency of exceeding a set of intensity thresholds (tephra load, in $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$) 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).

[Figure omitted. See PDF]

4.2.3 Visualization of the overall disaggregation

In order to have an overall picture of hazard disaggregation, in Fig.  we show the results of the relative dominance per volcanic source and eruption size class for the threshold of 300 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$. As already observed in Fig. a, 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 $\sim \mathrm{30}\mathrm{km}$, as indicated in Fig. a. 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 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ conditional to the occurrence of an eruption of medium size on land . 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 $\sim \mathrm{10}\mathrm{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.  we show the disaggregation by volcanic source and eruptive size class corresponding to the tephra threshold of 10 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$. 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 $\sim \mathrm{50}\mathrm{km}$ from the vent and the large class over $\sim \mathrm{50}\mathrm{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 ). On the contrary, Campi Flegrei shows a reduced-area largest dominance because the mean annual frequency to overcome the 10 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ 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.  we show the relative dominance $D_{ik}$ of each eruption size class of Somma–Vesuvius in overcoming the tephra threshold of 10 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ in 50 years, noting that the small class contributes to this hazard intensity only within distances of $\sim \mathrm{60}\mathrm{km}$ from the vent, with a relative dominance of $\sim \mathrm{60}\mathit{\%}$ only in the proximal areas (Fig. a), while the medium class contributes within distances of $\sim \mathrm{100}\mathrm{km}$ (including the Gulf of Naples, the northern part of the peninsula Sorrentina, and Ischia; Fig. b). 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. c). In Fig. b–d, we report the same type of results in overcoming the threshold of 300 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ 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 $\mathit{\lambda }$ 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 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$. 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.  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. a) 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 $\sim$ 10–20 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$, but for higher thresholds (from $\sim \mathrm{20}$ to 1000 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$) Campi Flegrei shows a higher hazard curve (Fig. b). This is in accordance with , 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 , and the city center is closer to the eruptive vents of Campi Flegrei and more directly downwind e.g.,. Beyond the 1000 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ 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. c), being located downwind. The results shown in Figs. – 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 ($\sim \mathrm{34}\mathit{\%}$ in 50 years vs. $\sim \mathrm{3}\mathit{\%}$ and $\sim \mathrm{12}\mathit{\%}$ in 50 years, respectively; . Wind directions statistically prevail towards the east in the Neapolitan area e.g., with moderate seasonal variations e.g.,: this aspect affects the model results, making tephra load dispersion oriented (e.g., Figs. –). 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 . 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., ), conditional to the occurrence of eruptions at one volcano e.g.,. In this sense, our PVHA is, for the first time, fully compatible with other hazard studies like the ones for seismic events . 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 . 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., and the Campanian Ignimbrite e.g.,, are not taken into account because of their negligible probability of occurrence over an exposure time of 50 years e.g.,. 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.,. 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 .

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 , 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 . 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 ($\mathrm{kg}{\mathrm{s}}^{-\mathrm{1}}$) are reported.

Volcano	Size class	TEM (kg)	Fallout mass (kg)	Fallout duration (h)	MER ($\mathrm{kg}{\mathrm{s}}^{-\mathrm{1}}$)
Somma–Vesuvius	Small	${\mathrm{10}}^{\mathrm{10}}$–${\mathrm{10}}^{\mathrm{11}}$	$\mathrm{8}\times {\mathrm{10}}^{\mathrm{9}}$–$\mathrm{8}\times {\mathrm{10}}^{\mathrm{10}}$	$\mathrm{11.11}$–$\mathrm{87.60}$	$\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{4}}$–$\mathrm{2}\times {\mathrm{10}}^{\mathrm{6}}$
	Medium	${\mathrm{10}}^{\mathrm{11}}$–${\mathrm{10}}^{\mathrm{12}}$	$\mathrm{8}\times {\mathrm{10}}^{\mathrm{10}}$–$\mathrm{8}\times {\mathrm{10}}^{\mathrm{11}}$	$\mathrm{6.24}$–$\mathrm{11.11}$	$\mathrm{2}\times {\mathrm{10}}^{\mathrm{6}}$–$\mathrm{2}\times {\mathrm{10}}^{\mathrm{7}}$
	Large	${\mathrm{10}}^{\mathrm{12}}$–${\mathrm{10}}^{\mathrm{13}}$	$\mathrm{8}\times {\mathrm{10}}^{\mathrm{11}}$–$\mathrm{8}\times {\mathrm{10}}^{\mathrm{12}}$	$\mathrm{6.14}$–$\mathrm{6.24}$	$\mathrm{2}\times {\mathrm{10}}^{\mathrm{7}}$–$\mathrm{3.6}\times {\mathrm{10}}^{\mathrm{8}}$
Campi Flegrei	Small	${\mathrm{10}}^{\mathrm{10}}$–${\mathrm{10}}^{\mathrm{11}}$	$\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{9}}$–$\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{10}}$	$\mathrm{3.48}$–$\mathrm{27.36}$	$\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{4}}$–$\mathrm{1.9}\times {\mathrm{10}}^{\mathrm{5}}$
	Medium	${\mathrm{10}}^{\mathrm{11}}$–${\mathrm{10}}^{\mathrm{12}}$	$\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{10}}$–$\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{11}}$	$\mathrm{1.95}$–$\mathrm{3.48}$	$\mathrm{1.9}\times {\mathrm{10}}^{\mathrm{6}}$–$\mathrm{3.6}\times {\mathrm{10}}^{\mathrm{7}}$
	Large	${\mathrm{10}}^{\mathrm{12}}$–${\mathrm{10}}^{\mathrm{13}}$	$\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{11}}$–$\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{12}}$	$\mathrm{1.92}$–$\mathrm{1.95}$	$\mathrm{3.6}\times {\mathrm{10}}^{\mathrm{7}}$–$\mathrm{3.6}\times {\mathrm{10}}^{\mathrm{8}}$
Ischia	Large	${\mathrm{10}}^{\mathrm{11}}$–$\mathrm{3}\times {\mathrm{10}}^{\mathrm{11}}$	${\mathrm{10}}^{\mathrm{11}}$–$\mathrm{3}\times {\mathrm{10}}^{\mathrm{11}}$	2–620	$\mathrm{1.3}\times {\mathrm{10}}^{\mathrm{5}}$–$\mathrm{1.3}\times {\mathrm{10}}^{\mathrm{7}}$

Table A2

Parameters of the beta distribution best fitting the field-based total grain size distributions (TGSDs) for Campi Flegrei . In this case, the reported parameters are referred to as the juvenile, lithics, and crystals; $p$ and $(\mathrm{1}-p)$ are, respectively, the coarse and fine sub-population weights; ${\mathit{\mu }}_{\mathrm{1}}$, ${\mathit{\mu }}_{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{1}}$, ${\mathit{\sigma }}_{\mathrm{2}}$ are the mean and standard deviations of the two Gaussian distributions in $\mathrm{\Phi }$ units.

Volcano	Size class	$p$	${\mathit{\mu }}_{\mathrm{1}}$	${\mathit{\sigma }}_{\mathrm{1}}$	${\mathit{\mu }}_{\mathrm{2}}$	${\mathit{\sigma }}_{\mathrm{2}}$
Somma–Vesuvius	Small	$[\mathrm{0.5},\mathrm{0.6}]$	$[-\mathrm{2.8},-\mathrm{1.4}]$	$[\mathrm{1.3},\mathrm{1.8}]$	$[-\mathrm{0.16},\mathrm{3.2}]$	$[\mathrm{1},\mathrm{2.8}]$
	Medium	$[\mathrm{0.6},\mathrm{0.8}]$	$[-\mathrm{1.3},\mathrm{0.5}]$	$[\mathrm{2},\mathrm{3}]$	$[\mathrm{4},\mathrm{5.2}]$	$[\mathrm{0.7},\mathrm{1.5}]$
	Large	$[\mathrm{0.7},\mathrm{0.8}]$	$[-\mathrm{0.5},\mathrm{0.5}]$	$[\mathrm{1},\mathrm{3}]$	$[\mathrm{4},\mathrm{6}]$	$[\mathrm{0.5},\mathrm{1.5}]$
Ischia	Large	$[\mathrm{0.06},\mathrm{0.34}]$	$[-\mathrm{2.36},\mathrm{1.28}]$	$[\mathrm{0.8},\mathrm{2}]$	$[\mathrm{3.50},\mathrm{6.35}]$	$[\mathrm{1.05},\mathrm{1.87}]$

Table A3

Parameters of the beta distribution best fitting the field-based TGSDs for Campi Flegrei . In this case, the reported parameters are referred to as the juvenile, lithics, and crystals; $p$ and $(\mathrm{1}-p)$ are, respectively, the coarse and fine sub-population weights; ${\mathit{\mu }}_{\mathrm{1}}$, ${\mathit{\mu }}_{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{1}}$, ${\mathit{\sigma }}_{\mathrm{2}}$ are the mean and standard deviations of the two Gaussian distributions in $\mathrm{\Phi }$ units.

CF eruption size class	Particle	$p$	${\mathit{\mu }}_{\mathrm{1}}$	${\mathit{\sigma }}_{\mathrm{1}}$	${\mathit{\mu }}_{\mathrm{2}}$	${\mathit{\sigma }}_{\mathrm{2}}$
Small	Juvenile	$[\mathrm{0.3},\mathrm{0.5}]$	$[\mathrm{2.7},\mathrm{4}]$	$[\mathrm{0.9},\mathrm{1.5}]$	$[\mathrm{4.9},\mathrm{5.4}]$	$[\mathrm{1.2},\mathrm{1.6}]$
	Lithics	$[\mathrm{0.3},\mathrm{0.5}]$	$[-\mathrm{2},-\mathrm{0.5}]$	$[\mathrm{1.4},\mathrm{1.7}]$	$[\mathrm{4.9},\mathrm{5.4}]$	$[\mathrm{1.2},\mathrm{1.6}]$
	Crystals	$[\mathrm{0.3},\mathrm{0.5}]$	$[\mathrm{0.1},\mathrm{0.6}]$	$[\mathrm{0.7},\mathrm{1.2}]$	$[\mathrm{4.9},\mathrm{5.4}]$	$[\mathrm{1.2},\mathrm{1.6}]$
Medium	Juvenile	$[\mathrm{0.2},\mathrm{0.4}]$	$[\mathrm{2.7},\mathrm{4}]$	$[\mathrm{0.9},\mathrm{1.5}]$	$[\mathrm{4.9},\mathrm{5.4}]$	$[\mathrm{1.2},\mathrm{1.6}]$
	Lithics	$[\mathrm{0.2},\mathrm{0.4}]$	$[-\mathrm{2},-\mathrm{0.5}]$	$[\mathrm{1.4},\mathrm{1.7}]$	$[\mathrm{4.9},\mathrm{5.4}]$	$[\mathrm{1.2},\mathrm{1.6}]$
	Crystals	$[\mathrm{0.2},\mathrm{0.4}]$	$[\mathrm{0.1},\mathrm{0.6}]$	$[\mathrm{0.7},\mathrm{1.2}]$	$[\mathrm{4.9},\mathrm{5.4}]$	$[\mathrm{1.2},\mathrm{1.6}]$
Large	Juvenile	$[\mathrm{0},\mathrm{0.3}]$	$[\mathrm{2.7},\mathrm{4}]$	$[\mathrm{0.9},\mathrm{1.5}]$	$[\mathrm{4.9},\mathrm{5.4}]$	$[\mathrm{3.5},\mathrm{5.5}]$
	Lithics	$[\mathrm{0},\mathrm{0.3}]$	$[-\mathrm{2},-\mathrm{0.5}]$	$[\mathrm{1.4},\mathrm{1.7}]$	$[\mathrm{4.9},\mathrm{5.4}]$	$[\mathrm{3.5},\mathrm{5.5}]$
	Crystals	$[\mathrm{0},\mathrm{0.3}]$	$[\mathrm{0.1},\mathrm{0.6}]$	$[\mathrm{0.7},\mathrm{1.2}]$	$[\mathrm{4.9},\mathrm{5.4}]$	$[\mathrm{3.5},\mathrm{5.5}]$

Table A4

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

Parameters	Somma–Vesuvius	Campi Flegrei	Ischia
Aggregation model	Cornell	Percentage	Cornell
Source type	Suzuki $A\in [\mathrm{3},\mathrm{5}],\mathit{\lambda }=\mathrm{1}$	Suzuki $A\in [\mathrm{3},\mathrm{5}],\mathit{\lambda }=\mathrm{1}$	Suzuki $A\in [\mathrm{3},\mathrm{5}],\mathit{\lambda }=\mathrm{1}$
Vent latitude (${}^{\circ }$)	$\mathrm{40.82}$	$\mathrm{40.82}$	$\mathrm{40.73}$
Vent longitude (${}^{\circ }$)	$\mathrm{14.42}$	$\mathrm{14.13}$	$\mathrm{13.92}$
Vent altitude (m a.s.l.)	$\mathrm{1281}$	$\mathrm{458}$	$\mathrm{251}$
Terminal velocity model	Ganser	Ganser	Ganser
Horizontal turbulence model	CMAQ	CMAQ	CMAQ
Vertical turbulence model	Similarity	Similarity	Similarity

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 (; Table ). 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 . In the case of Ischia, we assume that the TEM associated with the Cretaio Tephra eruption is almost entirely due to a fallout phase .

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 , sampling the eruptive column height from a beta distribution ($\mathit{\alpha }=\mathrm{2}$ and $\mathit{\beta }=\mathrm{3}$, scaled by a factor of $\mathrm{10}$ km and translated by 5 km such that $\stackrel{\mathrm{‾}}{H}=\mathrm{9}\mathrm{km}$) within the interval 5–15 km as in. 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 .

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. ) for the variable $\mathrm{\Phi }$ (grain size, considering that $d={\mathrm{2}}^{-\mathrm{\Phi }}$ is the particle diameter in mm): 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}}}},$$ where $p$ and $\mathrm{1}-p$ are the fine and coarse sub-population weights, and ${\mathit{\mu }}_{\mathrm{1}}$, ${\mathit{\sigma }}_{\mathrm{1}}$ and ${\mathit{\mu }}_{\mathrm{2}}$, ${\mathit{\sigma }}_{\mathrm{2}}$ are the mean and standard deviation of the two Gaussians, respectively. The five parameters of the distribution ($p$, ${\mathit{\mu }}_{\mathrm{1}}$, ${\mathit{\sigma }}_{\mathrm{1}}$, ${\mathit{\mu }}_{\mathrm{2}}$, ${\mathit{\sigma }}_{\mathrm{2}}$) are defined for each eruption size class according to what is reported in the literature and randomly sampled from beta distributions (Tables , ). In the case of Ischia, TGSD is estimated following .

For Somma–Vesuvius and Ischia, we consider a bulk class of particles; for Campi Flegrei, the availability of literature data allowed 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 .

Tephra particle density is set as follows: from 900 to 2500 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ for Somma–Vesuvius seeand references therein and from 740 to 2460 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ for Ischia . For Campi Flegrei, we distinguished the density of lithics (2500 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$) and crystals (2800 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$) as in , considering the estimates provided by .

In Table  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 , which assumes a mushroom-like vertical distribution of emission points depending on two dimensionless parameters $\mathit{\lambda }$ and $A$. The parameter $\mathit{\lambda }$, introduced by , 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, $\mathit{\lambda }$ is set to 1 and $A$ randomly sampled in the range $[\mathrm{3},\mathrm{5}]$.

The aggregation of tephra particles affects the sedimentation dynamics and deposition e.g.,. FALL3D, besides the option based on the model of , 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 , selecting a single tephra bin for Somma–Vesuvius (diameter of 200 $\mathrm{\mu }$m, density sampled within a uniform distribution, 100–600 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$) and Ischia (diameter of 200 $\mathrm{\mu }$m, density sampled within a uniform distribution, 50–500 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$). For Campi Flegrei, we use the percentage model , selecting two classes of aggregates for each eruption size class. For accretionary lapilli the diameter is set to 2000 $\mathrm{\mu }$m, and density is sampled in the range 1000–2000 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$. For other aggregates, the diameter is set to 200 $\mathrm{\mu }$m, and density varies in the range 100–600 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ as in.

Appendix B disaggregation per eruption size class for Campi Flegrei

Figure  reports the size class disaggregation analysis (see Eq. ) for Campi Flegrei with a ground load threshold of 300 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$. Here we present the same information for the threshold of 10 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$.

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 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ 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 }_{\mathrm{10}}$ number of simulations that produce a tephra ground load exceeding the threshold of 10 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$.

[Figure omitted. See PDF]

Appendix C mean annual eruption rates

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

Table C1

Mean annual rates (${\mathit{\nu }}_{ik}$) calculated for each eruption size class $k$ of the Somma–Vesuvius, Campi Flegrei, and Ischia volcanoes in 50 years.

Somma–Vesuvius	Small	Medium	Large
	$\mathrm{6.0}\times {\mathrm{10}}^{-\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}}$	$\mathrm{2.08}\times {\mathrm{10}}^{-\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}}$	$\mathrm{7.13}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{a}}^{-\mathrm{1}}$
Campi Flegrei	Small	Medium	Large
	$\mathrm{3.47}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{a}}^{-\mathrm{1}}$	$\mathrm{1.27}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{a}}^{-\mathrm{1}}$	$\mathrm{3.9}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{a}}^{-\mathrm{1}}$
Ischia	Small	Medium	Large
	–	–	$\mathrm{1.89}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{a}}^{-\mathrm{1}}$

Table C2

Mean annual rates (${\mathit{\nu }}_{ik}$) calculated for each eruption size class $k$ of the Somma–Vesuvius, Campi Flegrei, and Ischia volcanoes in 5 years.

Somma–Vesuvius	Small	Medium	Large
	$\mathrm{5.4}\times {\mathrm{10}}^{-\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}}$	$\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}}$	$\mathrm{8.6}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{a}}^{-\mathrm{1}}$
Campi Flegrei	Small	Medium	Large
	$\mathrm{3}.\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{a}}^{-\mathrm{1}}$	$\mathrm{8.8}\times {\mathrm{10}}^{-\mathrm{5}}{\mathrm{a}}^{-\mathrm{1}}$	$\mathrm{4.7}\times {\mathrm{10}}^{-\mathrm{5}}{\mathrm{a}}^{-\mathrm{1}}$
Ischia	Small	Medium	Large
	–	–	$\mathrm{1.7}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{a}}^{-\mathrm{1}}$

Table C3

Mean annual rates (${\mathit{\nu }}_{ik}$) calculated for each eruption size class $k$ of the Somma–Vesuvius, Campi Flegrei, and Ischia volcanoes in 500 years.

Somma–Vesuvius	Small	Medium	Large
	$\mathrm{7.8}\times {\mathrm{10}}^{-\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}}$	$\mathrm{2.2}\times {\mathrm{10}}^{-\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}}$	$\mathrm{3.2}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{a}}^{-\mathrm{1}}$
Campi Flegrei	Small	Medium	Large
	$\mathrm{9.1}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{a}}^{-\mathrm{1}}$	$\mathrm{3.2}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{a}}^{-\mathrm{1}}$	$\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{a}}^{-\mathrm{1}}$
Ischia	Small	Medium	Large
	–	–	$\mathrm{4.3}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{a}}^{-\mathrm{1}}$

Code and data availability

NetCDF files containing the long-term hazard map data are available in the Supplement. The codes used for the analysis are available on request by contacting the authors.

The supplement related to this article is available online at: https://doi.org/10.5194/nhess-23-2289-2023-supplement.

Author contributions

AC, LS, JS, and RS: conceptualization and design of the study, analysis and/or interpretation of data, and drafting of the manuscript. BMM, SM, MS, LS, JS, MTL, and AC: coding of the scripts, software and analysis, and interpretation of data. SM, MS, BMM, AC, JS, and LS: drafting of the manuscript. BG, MM, EP, MN, RI, GN, and PD: contributions to the interpretation of the results. AC: coordination of the study. All the authors revised and accepted the submitted version of the manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

We thank the Copernicus Climate Change Service ERA5, fifth generation of ECMWF atmospheric reanalyses of the global climate, and the Copernicus Climate Change Service Climate Data Store (CDS) for the access to the meteorological data. MS thanks Stefano Cacciaguerra for the academic and personal support. Giovanni Macedonio and Arnau Folch are warmly acknowledged for very useful and fruitful discussions. We also thank the CSI group of INGV Bologna for the technical support for using the cluster.

Financial support

The research leading to these results was supported by the Italian Civil Protection Department of the Presidency of the Council of Ministers within the Contratto concernente l'affidamento di servizi per il programma per il supporto al rafforzamento della Governance in materia di riduzione del rischio sismico e vulcanico ai fini di protezione civile nell'ambito del PON Governance e Capacità Istituzionale 2014–2020 (grant no. CIG6980737E65). The simulations for Campi Flegrei were supported by the EU project (ChEESE (grant no. 823844)), and the computational resources were provided within the multi-year PRACE Project Access Volcanic Ash Hazard and Forecast (grant no. 2019215114). This part of the work was granted access to the HPC/AI resources of TGCC (grant no. 2020-RA2020235114) made by GENCI. SM was also supported by the PON Research and Innovation 2014–2020 project referring to research contracts on green topics (grant no. CODICE CUP H95F21001440006). This research was also supported by the RETURN Extended Partnership and received funding from the European Union's NextGenerationEU (National Recovery and Resilience Plan – NRRP, Mission 4, Component 2, Investment 1.3 – D.D. 1243 2/8/2022, PE0000005). Manuel Stocchi was also supported by the PNIR – Programma Nazionale Infrastrutture di Ricerca with the CIR01_00013 project.

Review statement

This paper was edited by Maria Ana Baptista and reviewed by Ricardo Ramalho and Karoly Nemeth.