Assessing future vent opening locations at the Somma‐Vesuvio volcanic complex: 2. Probability maps of the caldera for a future Plinian/sub‐Plinian event with uncertainty quantification

In this study, we combine reconstructions of volcanological data sets and inputs from a structured expert judgment to produce a first long‐term probability map for vent opening location for the next Plinian or sub‐Plinian eruption of Somma‐Vesuvio. In the past, the volcano has exhibited significant spatial variability in vent location; this can exert a significant control on where hazards materialize (particularly of pyroclastic density currents). The new vent opening probability mapping has been performed through (i) development of spatial probability density maps with Gaussian kernel functions for different data sets and (ii) weighted linear combination of these spatial density maps. The epistemic uncertainties affecting these data sets were quantified explicitly with expert judgments and implemented following a doubly stochastic approach. Various elicitation pooling metrics and subgroupings of experts and target questions were tested to evaluate the robustness of outcomes. Our findings indicate that (a) Somma‐Vesuvio vent opening probabilities are distributed inside the whole caldera, with a peak corresponding to the area of the present crater, but with more than 50% probability that the next vent could open elsewhere within the caldera; (b) there is a mean probability of about 30% that the next vent will open west of the present edifice; (c) there is a mean probability of about 9.5% that the next medium‐large eruption will enlarge the present Somma‐Vesuvio caldera, and (d) there is a nonnegligible probability (mean value of 6–10%) that the next Plinian or sub‐Plinian eruption will have its initial vent opening outside the present Somma‐Vesuvio caldera.


Introduction
Somma-Vesuvio (SV) is one of the most studied and high-risk volcanoes in the world. Its eruptive history has been investigated through many studies since the first eyewitness account of the famous A.D. 79 Pompeii eruption by Plinius the Younger [e.g., Sigurdsson et al., 1985;Cioni et al., 1992Cioni et al., , 2008. Volcanic risk is very high here because surrounding areas are very densely inhabited, with over one million people directly threatened by the potential for devastating ash fallout, pyroclastic density currents, and lahars [DPC, 1995[DPC, , 2014Cioni et al., 2008; see Figure 1].
Over its history, the volcano has exhibited a large variability in eruptive styles as well as a moderate but significant spatial variability in vent locations. In Tadini et al. [2017], a detailed reconstruction of vent locations in past events is presented, together with an assessment of sources of associated spatial uncertainties. For instance, it is now established that, especially over the last 2 ka, many volcanic vents opened outside the boundary of the present SV caldera [Santacroce and Sbrana, 2003;Principe et al., 2013;Tadini et al., 2017], mainly on the western and southern flanks [Acocella et al., 2006;Paolillo et al., 2016]. In these cases, • New vent opening probability maps for a future Plinian/sub-Plinian eruption of Somma-Vesuvio volcano are developed • The maps incorporate uncertainty estimates and uniquely include the special case of Somma-Vesuvio caldera enlargement • New volcanological/structural data sets (from the companion paper) are combined to create spatial probability density functions Supporting Information: • Supporting Information S1 • Table S1 • Table S2 • Table S3 • Table S4 • Data Set S1 • Data Set S2 • Data Set S3 • Data Set S4 • Data Set S5 • Data Set S6 Correspondence to: A. Tadini, alessandro.tadini@ingv.it Citation: Tadini, A., et al. (2017), Assessing future vent opening locations at the Somma-Vesuvio volcanic complex: 2. Probability maps of the caldera for a future Plinian/sub-Plinian event with uncertainty quantification, J. Geophys. Res. Solid Earth, 122, 4357-4376, doi:10.1002 however, volcanic activity involved just effusive lava emissions or small-scale explosive eruptions from parasitic vents (e.g., Strombolian to Violent Strombolian events ).
As far as Plinian and sub-Plinian eruptions are concerned, field data indicate that their volcanic vents significantly varied in location but were always located within the present SV caldera. In particular, according to Cioni et al. [1999], the present outline of the SV caldera (an area about 4 km × 3 km) is the result of multistage collapses after four major Plinian eruptions, with minor contributions from three sub-Plinian I eruptions . During caldera-forming eruptions, the initial vent was always located inside the previously formed caldera, or along its border, but field evidence suggests a subsequent migration of the vent as a result of collapses, with consequent enlargement of the caldera itself. This variability appears to have had a strong influence on the distribution of eruptive products around the volcano, particularly those associated with pyroclastic density currents (PDCs). Evidence for such effects are found in the reconstruction of PDC deposits [e.g., Gurioli et al., 2010] as well as in the outcomes of transient 3D numerical simulations of column-collapse scenarios, which highlight a major effect of vent location and proximal volcano topography on the area invaded by PDCs [Esposti Ongaro et al., 2008a, 2008b.
Based on this knowledge, it is evident that a probability map of vent opening within or around the present SV caldera, together with estimates of associated uncertainties, would represent a crucial, objective basis for the assessment of volcanic hazards from a future explosive eruption at SV volcano. Elsewhere, vent opening probability maps, sometimes including information on associated uncertainty, have been produced to assess potential future activity in monogenetic volcanic fields [Connor et al., 2012], on volcanic islands [Alberico et al., 2008;Marti and Felpeto, 2010;Becerril et al., 2013], at composite stratovolcanoes [Cappello et al., 2012], and at calderas [Alberico et al., 2002;Selva et al., 2012;Bevilacqua et al., 2015]. The only approach which considered vent opening variability at SV is that proposed by Selva et al. [2014], who used a vent opening area solely for the evaluation of tephra fall hazard. In this latter case, the authors defined a vent opening area consisting of a circle of 6 km radius subdivided into five distinct areas, each with constant probability of vent opening. That area included several circum-Vesuvian towns, like Torre del Greco and Ottaviano among others, and neglected any volcanological and structural information about past vents and faults. In contrast, numerical modeling for tephra fallout and PDC hazard scenarios from SV explosive eruptions have, thus far, simply assumed a central vent corresponding to the present Gran Cono edifice [e.g., Neri et al., 2007;Esposti Ongaro et al., 2008a, 2008bMacedonio et al., 2016].
The aim of this manuscript is to present the first long-term (also referred to as the background or base-rate) probability maps of vent opening for the SV caldera that incorporates information and related uncertainty about some key parameters of the volcanic system. In doing so, a single map for a specific scenario is not presented, but rather a set of three maps are developed: a spatial distribution map of the mean probability of vent opening location, together with two subsidiary maps of the associated confidence percentiles, based on epistemic uncertainty quantifications derived from a structured elicitation of specialists, using alternative elicitation procedures to combine judgments. This set of three maps express the probability of vent opening conditional on the occurrence of either a new Plinian or sub-Plinian eruption in the foreseeable future, realized by considering the known eruptive record of SV over the last about 20 ka, as well as the distribution of key structural features of the caldera.
For uncertainty quantification, we relied on just two main contributory sources of uncertainty, which can be classed as epistemic (i.e., knowledge-related) uncertainty and natural physical variability (sometimes called aleatoric uncertainty). Epistemic uncertainty is related to the limitations in knowledge we have about the volcanic system which are present in both the data sets and their interpretation and could be reduced by obtaining more or better data. Physical variability is the intrinsic randomness in any evolving natural process which, in our specific case, will limit the forecast accuracy of the next vent location (and is not reduced or removed by acquiring evermore data from the past). While a strong dichotomy between uncertainty on the knowledge of the system (epistemic) and uncertainty on its future behavior (physical/ aleatoric) is not universally accepted, it is a simplifying, pragmatic paradigm adopted in many natural hazard assessments [e.g., Marzocchi and Bebbington, 2012;Selva et al., 2012]-most notably in probabilistic seismic hazard assessment, where there is an extensive and exhaustive literature. For recent reviews in relation to several natural hazards, see, e.g., Beven et al. [2015Beven et al. [ , 2016. We adopt these two main classes because they are consistent with the "doubly stochastic modeling" concept which comes from statistics literature [e.g., Cox and Isham, 1980;Daley and Vere-Jones, 2008;Harte and Vere-Jones, 2005] and has been recently applied to volcanology [e.g., Jaquet et al., 2008Jaquet et al., , 2012. In such a model, the probability values representing the spatial aleatoric variability (or uncertainty) related to the vent opening process are themselves potentially affected by nontrivial epistemic uncertainty: this dual classification was also adopted within doubly stochastic models applied to volcanic systems as in Bevilacqua et al. [2015], Neri et al. [2015], and Bevilacqua et al. [2016] and is fully described in Bevilacqua [2016]. While a more detailed uncertainty classification is possible [e.g., Tierz et al., 2016aTierz et al., , 2016b, with multiple classes of uncertainties instead of just two, the benefits of adding such additional complexity to the analysis structure are regarded by us as marginal, at best.
By accounting for the possible variability of vent location [e.g., as for the Campi Flegrei caldera; Selva et al., 2010;Neri et al., 2015], the maps derived here represent an important new information basis for producing long-term, or background, probabilistic hazard maps for the main phenomena associated to medium-large scale eruptions (e.g., Plinian and sub-Plinian events) at SV. In case of future unrest, these maps will offer a baseline configuration from which maps of vent position probability can be updated, informed by real-time monitoring data [Selva et al., 2012;Bevilacqua et al., 2015].

Methods
The method adopted here assumes that the probability of new vent opening can be computed as a weighted linear combination of the spatial distributions of key physical variables of the system that reflect, or can influence, this process. In other words, we presume ab initio that a new vent will likely open close to previous ones and that its location will be influenced by the presence of geological structures, such as major faults. A very Journal of Geophysical Research: Solid Earth 10.1002/2016JB013860 similar method was applied by Bevilacqua et al. [2015] for the definition of vent opening maps at Campi Flegrei caldera (Italy), whereas similar approaches, but involving some different techniques, have been applied elsewhere, for instance, by Marti and Felpeto [2010], Cappello et al. [2012], Connor et al. [2012], Selva et al. [2012], and Bartolini et al. [2013] for mapping vent opening at other effusive and explosive volcanoes.
In this study we used data from literature as well as new information presented in Tadini et al. [2017]. The data sets considered in that analysis were (1) the distribution of eruptive vents which produced past Plinian and sub-Plinian events; (2) the distribution of vents which produced explosive events such as Violent Strombolian (VS) and Continuous Ash Emission (AE) eruptions; (3) the spatial distribution of effusive vents and eruptive fissures in the past; and (4) the spatial distribution of deformation structures. Given the volcanological and geological data available and based on our present understanding of the volcano, the distributions of these data, here assumed to reflect the physical variability of vent opening processes, emerged as candidates for being most closely correlated with vent opening potential, with major faults being indicative of sectors of crustal weakness inside the caldera.
This said, it is acknowledged that the probability of new intra-caldera vent opening could be correlated with other system variables or with complex processes that are not considered, due to lack of sufficient knowledge about them. To account for any contribution from these neglected factors and to represent missing information, we include a supplementary conservative spatially uniform distribution inside the SV caldera. For instance, with respect to the possible fault effect, it is important to point out that (i) a rising feeder dike could be captured by a preexisting fault but only if it is favorably oriented [e.g., Gaffney et al., 2007;Le Corvec et al., 2013] and (ii) fault zones involve the development from a main fault plane of a broader volume of distributed brittle deformation, called the "damage zone" [Kim et al., 2004], which might be important for defining the surface expression of an opening dyke [Mazzarini et al., 2013[Mazzarini et al., , 2016. Once the various data sets were defined [see Tadini et al., 2017], Gaussian kernels (Appendix A) with different bandwidths were applied to each variable to produce associated continuous-probability density maps.
An important aspect of the study of Tadini et al. [2017] is the quantification of some of the main sources of uncertainty associated with the available data. In particular, uncertainty in the location of past eruptive vents from analysis of deposits or ash dispersal was considered, as were the number and distribution of past events, documented in the chronicles as having occurred, but for which evidence of vent locations has not survived (so-called "lost vents"). Also mentioned in Tadini et al. [2017] was the uncertainty of weights to be ascribed to the different variables which, when combined linearly, serve to define the spatial probability of vent opening. The effects of this latter source of uncertainty are additionally detailed here, expressing median, 5th percentile, and 95th percentile values of the vent opening probability density maps.
For this part of the study, the structured expert judgment elicitation consisted of two sessions involving 15 participating specialists (included in the author list) with different levels of experience, scientific backgrounds, and from a variety of institutions. The main goal of the elicitation was to achieve transparent, robust, and shared distributional estimates for the weights to be attributed to the different variables considered. This objective was achieved by taking advantage of three well-established elicitation methods: Cooke's Classical Model [CM; Cooke, 1991;Aspinall, 2006], the Expected Relative Frequency method [ERF; Flandoli et al., 2011], and the basic Equal Weights (EW) rule [Bevilacqua, 2016]. A detailed account of the two elicitation sessions together with a description of the different pooling methods is provided in Appendix B. This approach differs from other studies in which mutually exclusive weights were deterministically assigned to variables with unknown values by the authors, rather than by a panel of experts [e.g., Selva et al., 2012;Bartolini et al., 2013], or through a cross-validation technique [e.g., Cappello et al., 2012].

Data Sets Description and Associated Spatial Density Distributions
In this section, the main features of the data sets considered in the study are briefly summarized. We refer to Tadini et al. [2017] for a more detailed description of the data sets, including the basis on which the uncertainty sources were quantified. For each data set, the corresponding continuous density functions representing the conditional probability of vent opening associated with the data sets are also presented. Maps were developed by applying a Gaussian kernel with the specific bandwidth adopted for each data set, as reported in Appendix A.
Journal of Geophysical Research: Solid Earth The six sub-Plinian eruptions include both sub-Plinian I and sub-Plinian II types, as defined by Cioni et al. [2008] and discussed in Tadini et al. [2017].
For this data set, we consider uncertainty areas centered on vent positions derived from the interpretation of isopachs/isopleths and caldera structural information; see Tadini et al. [2017]. In this companion paper, various uncertainty area extents were assumed, depending on the eruption type or the features of the specific event: (i) for Plinian eruptions, the uncertainty areas were taken as equal to the area of the related caldera collapse, as defined in Cioni et al. [1999]; (ii) for the sub-Plinian A.D. 1631 and A.D. 512 eruptions, the uncertainty area was defined as a circle that encloses the area of the Gran Cono; (iii) for the A.D. 472 Pollena eruption, the extent was assumed similar to that of the A.D. 1631 eruption, but centered on the Pollena vent location inferred by Sulpizio et al. [2005] and reshaped according to the SV present caldera outline; (iv) for the Greenish Pumices eruption, uncertainty area extent was assumed the same as the Pomici di Base eruption but centered on the Greenish Pumices vent location and reshaped according to the SV present caldera outline; and (v) for the AP1 and AP2 eruptions, the extent was the same as the Avellino eruption.
Based on the volcanological and structural evidences discussed in Tadini et al. [2017], each uncertainty area is also here assumed to enclose 100% probability of the location of the corresponding past vent, and within each of these areas, probability is assumed uniformly distributed. As described in Appendix A, Gaussian kernel are applied to these latter uniform distributions (Figure 2c and 2d); the bandwidth (h) is assumed equal to the mean minimum distance between the centroids of the circles/ellipses (220 m), since it is related to the spatial spread of the observed past vents. The distribution of the Plinian data set ( Figure 2c) shows that the density is quite smoothly spread over all the caldera, with maximum cell values located about 1 km west of the Gran Cono crater. Conversely, for the sub-Plinian data set (Figure 2d), the density shows significantly higher values with two peaks concentrated about 500 m northwest and about 2 km west of the Gran Cono.

Vent Location of Small-Scale Explosive Eruptions
The VS to AE eruptions have been included jointly into the small scale explosive eruptions data set, which comprises the vent locations of 32 events that span a wide temporal window between the Avellino eruption (4.3 ka BP) and the last eruption of Vesuvio of AD 1944 . The 32 vents of the VS eruptions and AE activity were mainly concentrated in the area of Gran Cono. However, due to a relative paucity of field data (possibly related to the scale of these eruptions), we cannot exclude that some events of a similar magnitude and intensity may have been lost in the stratigraphic record.
For the VS to AE data set, three different uncertainty areas for spatial localization of vents were defined as a function of the temporal frame considered (see Tadini et al. [2017], for a detailed explanation). These are (i) the present crater of Gran Cono enclosing the last 10 VS to AE events that occurred between the AD 1631 and AD 1944 eruptions ( Figure 3a); (ii) the area of the Gran Cono representing the uncertainty area of 22 vents of eruptions that occurred between the Avellino and the A.D. 1631 eruptions (Figure 3b); and (iii) the portion of SV caldera that accounts for all the VS to AE eruptions that may have been lost from the stratigraphic record ( Figure 3c). This portion of the SV caldera corresponds to the extent of the caldera just before the occurrence of the Avellino eruption.
As with the Plinian and sub-Plinian eruption cases, the uncertainty areas of these three maps are assumed to enclose 100% probability for the location of associated past vents and that such probability is distributed uniformly within the areas. However, in contrast with the Plinian and sub-Plinian data set, in this case each area refers to a group of eruptions, whose vent locations cannot be reconstructed, and not to a single event with an imprecise vent location. Figures 3d, 3e, and 3f show the density functions corresponding to the three areas in Figures 3a, 3b, and 3c, respectively, obtained by application of symmetric Gaussian kernels. Kernel bandwidth is assumed fixed at 100 m, corresponding to the cell size of the grid adopted.

Vent Location of Effusive Eruptions
Following the description reported in Tadini et al. [2017], the effusive eruptions data set is subdivided in two separate subdata sets: parasitic vents and eruptive fissures.

Parasitic Vents
Among the 99 parasitic vents cited in Tadini et al. [2017], 46 (from the period 1631-1944) were located within the SV caldera and have been directly used in this work ( Figure 4a). Each of these mapped parasitic vents was associated with a circular uncertainty area of 75 m radius (calculated from the mean radius of parasitic vents from different areas; Tadini et al. [2017]), i.e., an area taken to enclose their location. Also in this case, the uncertainty areas of these three maps are assumed to enclose 100% probability of the location of associated parasitic vents.

10.1002/2016JB013860
The major discrepancy between the number of vents cited in historical accounts and the number of mapped vents prompts consideration of the issue of so-called "lost vents." With the information available, it was possible in Tadini et al. [2017] to identify the approximate positions of some lost vents, but at no finer resolution than being situated somewhere within one of the different sectors of the SV caldera ( Figure 4a). For vent probability mapping, the arithmetic means of the minimum and maximum values calculated in Tadini et al.
[2017] are assumed: lost effusive vents are, therefore, four within sector A, five within sector B, and 21 with an unknown locality somewhere within the whole SV caldera (Figure 4a).
The related density function map for the effusive vent subdata set ( Figure 4b) has been produced by combining the map obtained after the application of the Gaussian kernel to the 46 mapped vents with a map where the lost vent locations are assumed uniformly distributed within their relative sectors. As with the large explosive eruptions data set, the kernel bandwidth is assumed equal to the mean minimum distance between the centroids of the circles (185 m), since bandwidth is related to the spatial spread of observed past vents. Figure 4b displays a map where the peaks of probability highlight the distribution of vents aligned along single eruptive fissures, and in this specific case, two directions of elongation of peak areas can be envisaged (NW-SE and NE-SW).

Eruptive Fissures
Only six of 32 eruptive fissures reported by Tadini et al. [2017] can be mapped in the field, while Acocella et al. [2006] inferred the locations of the rest from historical accounts. Since there is large uncertainty in eruptive fissure locations, the companion paper defined an uncertainty area reflecting a high density of  (Figure 4c), which roughly corresponds to the extent of the Gran Cono itself. The probability of eruptive fissures occurring within this area is considered uniform, as the general pattern of eruptive fissures around the Gran Cono is radial and almost equally distributed in all directions [Tadini et al., 2017].
The resulting related density function, obtained using a bandwidth for the kernel of 100 m, which equals the cell-size resolution (Figure 4d), indicates a uniform probability plateau, smoothly decaying toward the edges of the high-density area.

Deep Faults (Buried)
The "Deep faults" data set comprises faults interpolated mostly within the Mesozoic Quaternary carbonate basement [Tadini et al., 2017]. Among the faults that cross the SV caldera described in Tadini et al. [2017], some have been identified only after extrapolation from seismic sections [Bruno et al., 1998], a procedure that possibly limits the reliability of the data. Therefore, to provide a more self-consistent data set, we consider only faults that either have been interpolated reliably from seismic profiles or that have been identified in at least two different bibliographic sources. Thus, this data set comprises just three faults, two being fault planes extrapolated by Bruno et al. [1998] from seismic reflection surveying, but also confirmed by Ciaranfi et al. [1981] from combined gravimetric and seismic data; the other is a fault plane interpolated from seismic profiling by Bruno and Rapolla [1999].
Following Tadini et al. [2017], in order to take into account the spatial uncertainty related to fault positions, we consider an uncertainty area with a width of 150 m around each fault plane ( Figure 5a). For deep faults, the kernel bandwidth is assumed equal to the average damage zone (DZ) width which, in this study, is taken as proportional (linearly) to fault displacement. The latter is assumed equal to 190 m, an estimate based on the Scholz [2002] fault length (L) relationship for faults located within carbonate rocks (i.e., DZ = 3 × 10 À2 × L, with DZ and L in m). The map of Figure 5b shows probability densities that are focused quite close to mapped fault planes; however, these peak probability density values are only slightly higher than those from other data set maps (see sections 3.1 and 3.3.1).

Results
With the spatial density maps described above, we applied the structured expert elicitation techniques described in section 2 and Appendix B to ascribe a weight to each map and combine them into a joint probability map. An alternative uniform distribution over the whole caldera area was also adopted to represent the possibility that there may be no correlation between the vent opening distribution and the various mapped variables considered here.
As previously stated, three different expert judgment procedures were applied to the elicitation data: (i) the CM, (ii) the ERF method, and (iii) the EW rule. The output probability percentiles of the Decision Maker are represented by triangular distributions in the ERF case and by Beta distributions in the other two cases (CM and EW). The approximation with Beta distributions was obtained by choosing shape parameters that minimize the absolute errors on the three elicited percentiles (minimizing the sum of absolute errors gave consistent results-differences were less than 1% in the weights estimates).

Weights for the Variables
To simplify the quantification of weights for each spatial distribution, the experts were not asked to judge these directly. Instead, as in Bevilacqua et al. [2015], a simple hierarchical logic tree is defined (see Figure 6), where most of the target questions quantify the relative importance, or relevance, of one variable or feature of the system versus others. We follow a Monte Carlo simulation approach, multiplying the individual estimates over the branches of the logic tree to obtain Beta probability distributions for the nine linear weights.
In the logic tree of Figure 6, each branch represents a pair or triplet of target questions, with values elicited from the experts. Tables S1 and S2, in the supporting information, report respectively the initial and the revised questionnaires with target items and the solutions obtained from pooling judgments with the Classical Model (i.e., the so-called Decision Maker).
First, the experts were asked to evaluate the first pair of tree-branching probabilities: that the next Plinian/sub-Plinian eruption at SV will have its initial vent inside (versus outside) the present outline of the caldera. After that, when considering the location of the next future vent, the experts assessed the relative relevance of the information based on past volcanic activity compared to that for the distribution of deep faults and to a homogeneous distribution (Uniform map); this is at the second level of the tree in Figure 6. At the next level down on the tree, the experts compared the relative importance of the volcanic features' distributions associated with the different eruptive styles, i.e., Plinian/sub-Plinian versus VS/AE eruptions and versus effusive eruptions. Last, the perceived relative importance of past eruptive styles data sets was evaluated (i.e., Plinian versus sub-Plinian eruptions, VS to AE of different periods, parasitic vents versus eruptive fissures). Table 1 reports the median values and the 5th and 95th percentiles of the uncertainty distributions for the probability of having an initial vent inside or outside the caldera for the three different scoring rules considered. In this case, the median values are preferred as they are directly elicited from the experts (in the CM/EW cases), and they sum closer to 100% than the corresponding mean values. In the following, the mean values Figure 6. Logic tree of target questions. Structure for assigning weights to different items assessed by expert elicitation. Shaded (orange) colored boxes represent variables/data sets whose probability density maps have been linearly combined to produce the joint probability map.  Table 2 reports the mean values and the 5th and 95th percentiles of the linear weights for the nine contributing probability maps, obtained by Monte Carlo sampling from uncertainty distributions on the logic tree framework of elicited items; these results are the basis for producing the mean and percentile vent opening maps, discussed in the following paragraphs.
The weights assigned to different data sets/variables are found to be consistent between the three scoring rules (mean values differ by less than 4% in all the cases). With respect to the uncertainty ranges, these outcomes show that the ERF and CM models each produce narrower uncertainty intervals than the EW (equal weights) combination for mean estimates. The maps based on the CM results are adopted as reference here because they represent a rational, objective consensus on the key context of uncertainty quantification (tests of the CM and ERF methods have shown the former is generally more reliable for parameter uncertainty quantification, the latter for parameter central tendency accuracy; see Flandoli et al. [2011] for more details on this subject). This procedure has largely been adopted by the scientific community [e.g., Woo, 1999;Bevilacqua et al., 2015].
The probability density functions of the weights for the various data sets, as derived from the CM method (the reference one), are displayed in Figure 7. It is worth pointing out that (a) the higher mean weights have been attributed to the vent location of Plinian and sub-Plinian eruptions, along with the uniform map (equal to about 16.3%, 25.8%, and 19.9% respectively), and the resulting density functions have a quasinormal distribution with slightly wider uncertainty bounds; (b) density distributions of the three VS to AE eruptions variables have lower mean weights (i.e., 6.6%, 7.4% and 8.1%), with a distribution skewed toward the upper bound of uncertainty (95th percentile); a similar situation (even if with even lower mean values) is observed also for the effusive eruption variables; (c) the deep faults variable has an extremely skewed distribution, with a low mean value (i.e., 9.4%) but large 95th percentile value.

Sensitivity of Variable Weights to Group Composition
Sensitivity assessments on the elicitation outcomes were performed also with respect to individuals' experience and main scientific expertise, by partitioning participants into two pairs of subgroups: Group A1 (Seniors -10 experts), Group A2 (Juniors-5 experts), and also Group B1 (Geologists-10 experts) and Group B2 (Modelers-5 experts). Results are displayed in Tables S3 and S4 in the supporting information. The Decision Maker from the CM method applied within the B1 (Geologists) subgroup shows wider uncertainty ranges and a higher probability value for initial vent opening outside SV caldera compared to the reference (i.e., whole group Classical Model) Decision Maker (similar results to the Geologists group derived from the A2 subgroup: Juniors). Other differences between subgroups A1/A2 (Seniors/Juniors) and B1/B2 (Geologists/ assigned less weight to the sub-Plinian eruptions data set with respect to the other subgroups and, principally, gave higher weights both to the Parasitic Vents data set associated with effusive eruptions and to the probability of initial vent opening outside SV caldera; (iv) experts in the Geologists subgroup also assigned a lower 5th percentile weight to the Uniform probability map compared to the other subgroups, suggesting they may have greater confidence than other subgroups when it comes to the volcanological information used in the study. All the differences between groups are, however, small, demonstrating a general agreement within the group of experts on the topics being elicited. Figure 8 shows the vent opening location probability maps, corresponding to the 5th, mean, and 95th percentiles of the probability density values, obtained by weighting the density function maps of the different data sets according to the CM (taken here as reference model). Maps of these three percentiles are available respectively in Data Set S1, Data Set S2, and Data Set S3 in the supporting information.

Vent Opening Probability Maps 4.3.1. Probability Distribution Within the Caldera Boundary
The strong similarities between the maps obtained by the application of the CM, ERF, and EW methods (see Figure S1 from the supporting information) indicate that the outcomes of the approach are particularly informative about the group's consensus.  Alternative mean maps were produced (using the CM and EW methods) with respect to the above-mentioned subgroups of experts. Due to strong similarities with the reference maps, they are not presented here but are available in Figure S2 in the supporting information. All maps are presented on grids with cells 100 m on a side, and each probability density is expressed as percentage probability per cell (in square hectometers-hm 2 ) conditional to the occurrence of a Plinian/sub-Plinian eruption within the SV caldera (spatial integration of the mean maps is close to 100%).
Although the maps of Figure 8 show that vent opening location probability is quite diffuse over the whole SV caldera, it is clear that, for each map, the location of each maximal value roughly corresponds with the present crater summit (peak cell probability values range from 0.15% up to 0.69%). Moreover, the influence of Deep Faults is more pronounced in the 95th percentile maps-an effect due to skewness to higher values in the weighting distribution attributed to this data set by the experts (see Table 2 and Figure 7).

Vent Position During Caldera Enlargement
As previously mentioned, the effect of possible caldera enlargement in the case of a next large-scale explosive eruption has been analyzed based on the following assumptions: (a) caldera enlargement is possible due to the occurrence of a Plinian eruption; (b) caldera enlargement can be neglected in case of a sub-Plinian (or lower magnitude) eruption; (c) caldera enlargement can be modeled by assuming migration of the initial vent position of a Plinian eruption out to the limit of a collapse area related to the eruption.
Given these assumptions, vent position following a caldera collapse has been modeled with a simplified kernel function, called here a "collapse kernel" (Figure 9). This kernel function value is piecewise constant with respect to the distance from the origin, but null outside a circle with a radius of 1.3 km, which is the mean size of collapsed areas in the four previous Plinian, caldera-forming eruptions. The kernel sums to 1, and its density is distributed such that 60% is located up to 0.25 km from the center, 5% from 0.25 km to 1.05 km, and 35% from 1.05 km to 1.3 km. The collapse kernel function is convolved with the values of the initial vent opening map, in case of a Plinian eruption. These limits should suitably model a caldera collapse scenario where most of the erupted material is discharged through the central conduit area and along ring faults.
Relying on the findings of Neri et al. [2008], the probability P of having a Plinian eruption (conditional on having a Plinian or a sub-Plinian one) is here estimated assuming independent Beta distribution [as in Neri et al., 2008]; these results are significantly skewed and have [5th, 50th, 95th] uncertainty percentiles equal to [0.5%, 9.5%, 40%] and mean value 13.5%. During the Monte Carlo simulation for the "initial vent" opening probability map, we also sample the caldera collapse probability value P, and we apply the collapse kernel function only to that case. Specifically, the individual samples for the final vent opening maps ( Figure 10) are the linear combinations of the maps after "collapse kernel" application ( Figure S4 in the supporting information) and the maps without caldera collapse (Figure 8), multiplied by their respective probabilities: P and 1 À P.
The maps reported in Figure 10 thus represent the vent opening mean probability map and its companion credible range probabilities, such that these also include the possibility of caldera enlargement (associated only with Plinian events), conditional on the occurrence of a large eruption (Plinian or sub-Plinian) at SV. Probability percentiles of caldera enlargement (i.e., vent opening spatial probability integral outside the caldera boundary) are [0.8%, 1.9%, 5.8%], with a mean value of 2.4%. Maps from Figure 10 are available respectively in Data Set S4, Data Set S5, and Data Set S6 in the supporting information.

Expert Elicitation Procedure, Scoring Methods, and Vent Opening Probability Maps
During and after the first elicitation sessions, some technical and understanding issues were raised by some participating experts. After a discussion of the initial elicitation results, the elicitation was repeated for those items that had produced the most skewed responses (i.e., outside caldera probability and deep fault weight estimation) or produced debatable bimodal solutions. Other skewness-related issues are apparent in the graphical representation of Decision Maker percentile estimates with continuous density functions (with a strongly skewed distribution, the mean value is inevitably shifted toward the wider uncertainty side, with respect to the corresponding modal or median values).
Notwithstanding such conceptual challenges and considering that each of the pooling methods has its own advantages and drawbacks, there is overall agreement among the different scoring methods, which reflects a general accord among all the experts about the mean values, with strong similarity in the maps presented. Weights assigned to different data sets/variables are consistent between the three scoring rules (mean values differ by less than 2% for almost all cases)-a situation that suggests a fundamental consensus exists among the experts and is captured by the elicitation process and relevant Decision Maker solution. With respect to the uncertainty ranges, outcomes show, as expected, that the ERF model computes smaller uncertainty intervals on central tendency, the CM provides intermediate credible intervals values, and the EW model produces the widest uncertainty distributions, as is usual (for verification, see Colson and Cooke [2017]).
Two main observations can be made with respect to the judgments provided by the experts: (1) a higher weight is attributed to the Plinian/sub-Plinian eruption versus other data sets (i.e., VS to AE and effusive); and (2) the probability of initial vent opening outside the SV caldera is evaluated in the nonnegligible range 6-10%. The first outcome might be rationally explained by the fact that the initial aim of the elicitation was to quantify parameters and uncertainties with which to define a first vent opening probability map specifically for the next Plinian/sub-Plinian eruption at SV, and therefore, the elicited experts relied more heavily on data derived from similar past activity.
The second finding highlights a very challenging potential scenario, and the implications for volcanic hazard assessment could be serious (see section 5.2). The elicitation and modeling results indicate a mean probability value for this occurrence of 6% from the CM method; 10% using the ERF method, or 12% using the CM method with subgroup B1-Geologists. Three possible reasons have been hypothesized by the experts to consider that a vent for a future Plinian/sub-Plinian eruption could plausibly open outside the caldera, despite the fact that the volcanic history of SV tends to discount it having happened in the past [Tadini et al., 2017]. These reasons are as follows: (a) the topic is highly complex, involving multiple elements of process-related aleatoric and epistemic uncertainties, and thus, the experts adopted a precautionary approach that would admit the future possibility of such a scenario; (b) there is the possibility, supported by geophysical data, that a magma reservoir or sill, larger than the present edifice, exists 8-10 km below  [Auger et al., 2001]; and (c) there is an on-going debate, within the volcanological community, about a possible past sector collapse that involved a major part of the SV complex (discussed in Tadini et al. [2017]) that might be repeated in the future. With respect to the last point, such a major structural perturbation in the future could, theoretically, cause the volcano system to be "beheaded," resulting in subsequent vent migration. For the SV case, however, the sector collapse evidence has been strongly questioned [e.g., Sulpizio et al., 2008]. Moreover, Neri et al. [2008] evaluated this scenario as very unlikely (mean probability of occurrence conditional on the occurrence of a magmatic unrest of the volcano is 0.03%). In addition, whereas some shifts of vent locations in past high magnitude/intensity eruptions have been interpreted as related to activity on or around the rim structures of the previous caldera [Cioni et al., 1999;Tadini et al., 2017], some of the expert panel may have considered an alternative possibility that those vents were not controlled by preexisting structure but emerged independently outside (and not along the rim) of a preexisting caldera.
It is also worth reiterating that the "Deep faults" data set significantly influences vent opening probability maps, especially with regard to higher percentiles, because, as discussed above, the relevant PDFs are skewed with longer tails to higher values for this data set. Notwithstanding that the mean values of the ascribed weights do not reach extremely high values (i.e., only about 10%), the contribution of this data set to the density values is also increased because it is the most spatially concentrated of all the data sets, and the corresponding Gaussian kernel function does not spread the probability out to offset that concentration. Given the limited quality of these data (two of the three faults that comprise the data set have been extrapolated from adjacent seismic profiles), this aspect needs to be carefully considered when evaluating the influence of this data set on the final vent opening location probability maps. However, it is undoubtedly the case that large crustal discontinuities are likely to be a factor affecting vent distribution, and this suggests that, if possible, refining this data set is desirable.
To properly evaluate the net effect of the deep faults data set on the vent opening maps in case of a future Plinian/sub-Plinian eruption, we also produce maps using the CM and EW methods, in which the contribution of the deep faults data set is removed and its weight equally distributed over the other data sets. The rationale for this further analysis is that prevailing NW-SE and NE-SW trends (i.e., similar to the structure trends in the "Deep faults" data set) are recognizable without reference to the "Deep faults" data set; as mentioned above, in terms of reliability the latter is the most debatable of the data sets we utilize here. Results shown in Figure S3 in the supporting information highlight that, while significantly less evident with respect to Figure 8, the NE-SW and NW-SE trends can be still identified, particularly in the Piano delle Ginestre area and to the SE side of the Gran Cono [Tadini et al., 2017]. These trends might indirectly reflect the presence of buried main structures, with the same trend controlling the volcanic activity.

Implications for Volcanic Hazard Assessments
The vent opening location probability maps convey important implications in terms of hazard assessment for future Plinian and sub-Plinian eruptions. The mean map of the CM method indicates that less than 50% (90% credible interval: 36.0-55.6% probability) of the cumulative probability of vent opening location in the next Plinian/sub-Plinian eruption is found within the central area of the existing Gran Cono (Sector A of Figure 4a). There is about a 30% probability (90% credible interval: 22.4% and 36.7% probability) that the initial vent will be within the area of the Piano delle Ginestre (Sector D of Figure 4a), i.e., between about 1-3 km west of the Gran Cono crater. To give this context, PDCs from the Avellino Plinian eruption [4.3 ka BP; Cioni et al. [2008]], the initial vent of which was located in the Piano delle Ginestre region, reached the area presently occupied by the city of Naples (Figure 1). Considering the Valle del Gigante area (Sector B of Figure 4a), the mean cumulative probability that the vent of the next Plinian/sub-Plinian eruption will open in this area is around 16% (90% credible interval: 11.9% and 19.4% probability). It is pertinent for assessing potential hazard impacts that a vent opening in this area, even if associated with a sub-Plinian eruption, will engender a high likelihood that resulting PDCs could invade areas and municipalities located on the North flank of the volcano (Figure 1), overtopping the Mount Somma rim [Esposti Ongaro et al., 2008a, 2008b. Finally, the cumulative vent opening location probability for Sector C of Figure 4a (the Valle dell'Inferno area) yields a nonnegligible probability value of about 10% (90% credible interval: 7.5% and 13.7% probability).
These findings need to be considered also against the background of salient aspects of the spatial distribution of past volcanic activity at SV (as discussed in Tadini  Plinian and sub-Plinian eruptions were located within or near the outline of the then-existing SV caldera. Some Plinian/sub-Plinian eruptions exploited a preexisting vent, formed during an earlier eruption (e.g., the AP1/AP2 sub-Plinian eruptions exploited the same vent as the Avellino Plinian eruption), but not all followed this pattern (e.g., the Pompeii eruption did not exploit the same vent as AP1/AP2 nor the A.D. 512 sub-Plinian eruption used the same vent as the Pollena eruption). (b) There is an apparent tendency for SV volcanic activity to be centralized in the area of the Gran Cono edifice, but this is mainly evident only from the VS to AE eruptions data set, which might be incomplete. (c) The parasitic vents eruption data set does not highlight any particular temporal evolution of the spatial distribution of vents, at least over the past 2000 years. (d) The Gran Cono edifice and some groups of Plinian/sub-Plinian eruption sites are located at the intersection of two pairs of extrapolated deep faults.
As discussed in section 5.1, another noteworthy outcome of the present study relates to the relative probability that the next Plinian/sub-Plinian eruption might have its initial vent outside the present SV caldera. This is a challenging finding from this work and should be further investigated.
With respect to the SV conduit condition, that is, whether "open" or "closed," there is little or no other data available to inform considerations of future vent opening location from this perspective. For the case of an open conduit condition [see Acocella et al., 2006], we may conjecture (based on reasonable assumptions) that the areas more prone to vent opening outside the SV caldera, in a new eruption of any scale, would be those facing the sea (i.e., in the SW/SE sector), rather than in the N sectors (where the volcanic features have older ages), and that any new vent might be located within a distance of a few kilometers of the SV crater rim, approximately corresponding to the locations of parasitic vents in the Middle Ages activity (see Figure 5 in Tadini et al. [2017]). In contrast, for the case of a closed conduit condition (as at present), there are virtually no grounds for speculating about where a new vent might form, if outside the caldera. These particular issues pose an interesting and specific challenge for future research. In particular, it will be necessary to carefully evaluate (a) the extent of the area for possible sites of new vent openings and (b) the structures and volcanic features that need to be carefully considered and properly investigated.
Finally, while sites of effusive parasitic vents and eruptive fissures outside the SV caldera can be considered, the potential effect of regional/local structures on the distribution of future vent opening locations might also play a significant role. From this standpoint, structural data sets for areas outside the caldera appear to be more dependable than those considered for constructing our maps; the bulk of structural information outside the SV caldera area, as reported by Bruno et al. [1998], Bruno and Rapolla [1999], and Tadini et al. [2017, comprises interpolated data from seismic profiling and not extrapolated outward like the faults in the "Deep faults" data set.

Concluding Remarks
This study has produced the first long-term (base-rate) vent opening probability maps for the summit caldera area of the SV volcanic complex in the case of the next eruption being Plinian or sub-Plinian. The procedure implemented a structured and quantitative treatment of epistemic and physical uncertainties and their influences on analyses to determine where a new vent will first open in such an eruption. The vent opening location probability density maps, here presented as a mean map and associated uncertainty maps corresponding to the 5th and 95th percentile confidence levels, were obtained by linearly combining, with ascribed weights, nine different volcanological and geological data sets, representing the distribution of past vents for different eruptive categories and known fault information, together with an imprecise uniform distribution map to provide an "ignorance" reference.
Data set weights were defined through a procedure of structured elicitation from a group of experts, who all have experience of the SV volcanic complex, from various disciplines. Results from the elicitation were obtained using three different pooling procedures. Outcomes of this exercise include (a) the definition of continuous probability density functions, obtained through the application of symmetrical Gaussian kernels with appropriate bandwidths, for each individual data set/variable selected from the data presented in Tadini et al. [2017]; (b) the quantification of the weights, and their associated uncertainties, to be assigned to alternative vent location probability maps when linearly combined, based on performance-scored expert judgments; and (c) the application and intercomparison of different expert scoring methods, different Journal of Geophysical Research: Solid Earth 10.1002/2016JB013860 compositions of subgroups of experts, as well as the consideration of different sets of volcanological data (e.g., maps obtained with and without the contribution of information on deep faults).
Inspection of these probability maps, given that the next eruption at SV is Plinian/sub-Plinian, shows the following: 1. The mean probability of vent opening in the area of the present edifice (i.e., the Gran Cono area, assumed circular with a diameter of 1 km) is greater than elsewhere, but less than 50%. Uncertainty bounds around this mean value, expressed as 5th and 95th percentiles, correspond to a credible range between about 36% and 56% probability when the CM weights are used. 2. There is a significant probability, i.e., almost 30% as mean value, that the western portion of the SV caldera (i.e., Piano delle Ginestre area) could host the next vent opening. Uncertainty values correspond to a 90% credible range from about 22% to 37%, again referring to the CM findings; 3. There is a 2.4% mean estimated probability that the caldera will enlarge during the next Plinian/sub-Plinian event. 4. Despite the past evidence of SV history, there is at least a 6% probability that the next high-intensity eruption will have its initial vent outside the present outline of the SV caldera. This outcome, and its potentially significant hazard implications, demands further investigation.
All the findings of the present study appear robust with respect to the adoption of different scoring methods for combining expert judgments, with respect to different subsets of the group of experts and to the selection of basic volcanological variables considered in the analysis. Assessing how uncertainty in the potential location of a new vent opening at SV could influence predictions of tephra fallout or PDC inundation hazard mapping, if the next future eruption is sub-Plinian or Plinian, cannot be projected simply and directly from the findings reported here. Potential impacts of other factors, such as local topography, are many and complex, and thus further detailed investigation and modeling are warranted, as done, for example, for the nearby Campi Flegrei caldera [e.g., Selva et al., 2010;Neri et al., 2015].

Appendix A: Kernel Density Estimation
The probability density maps associated with each volcanic data set were generated from the application of Gaussian spatial density kernel to the uniform distribution assumed within the uncertainty areas enclosing the volcanic features [Connor and Hill, 1995;Connor and Connor, 2009].
Kernel density estimation [Silverman, 1986] builds up a continuous probability density from a number of discrete samplings or from a noncontinuous density function. The two main steps in spatial density estimation are the choices of the kernel function and of its bandwidth, or smoothing parameter. The kernel function can be any positive function K that integrates to 1 [Weller et al., 2006], and in general, given a finite sample X i , i = 1, …, N, a kernel density estimator can be defined as follows: where h is the bandwidth. In our study K is assumed equal to a two-dimensional radially symmetric Gaussian kernel, as with many kernel estimators used in geologic hazard assessments (the Gaussian distribution arises in problems of heat and mass transfer, such as those that might be expected in volcanic systems involving diffusion processes [Weller et al., 2006]).
The bandwidth h is typically selected using specific theoretical and empirical methods developed for optimizing consistency with data [Cappello et al., 2012;Becerril et al., 2013]. For instance, in the case of past vents locations, the bandwidth was associated with the average spatial spreading between a vent and others in the data set. Here h is assumed independent of actual location but does depend on the specific feature being considered (for instance, past vents or faults). An important characteristic of this study is that past locations of feature do not comprise simple points, but areas of uncertainty of different extent; each area can cover several cells in our computational grid, some of them completely, others only partially. Therefore, for each cell, it was taken into account the fraction of each uncertain area that it contains, and then we apply the appropriate kernel convolution to this value. In addition, it was also assumed that the kernel convolution does not spread the probability outside the SV caldera boundary.
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013860 An advantage of this approach is that the spatial density estimate will remain consistent with the spatial distribution of past volcanic features. A disadvantage of a symmetrical kernel function is that it does not explicitly allow for specific directionality of local geological and structural boundaries, or in other volcanological information. Additional details about the approach can be found in Bevilacqua et al. [2015].
Appendix B: Performance-Based Expert Judgment As described in the main text, the procedure of weights assignment to different data sets took advantage of two elicitation sessions. The first plenary session involved multiple presentations and discussion of the topics of concern-i.e., the "target" questions and a calibration questionnaire with appropriate "seed" questions. To facilitate the experts in providing responses to the target items, before the plenary session, a small compendium was provided to the participants that summarized the most important features about single data sets, and also the questionnaire with draft versions of the target questions. These were reviewed when the group convened and, if deemed necessary, modified before individuals then answered them confidentially and independently. During this first session, the calibration seed questions-used for scoring individuals' performances in the Cooke Classical Model Structured Expert Judgment method-were also answered.
These questions were based on carefully researched aspects of SV volcanism, other Italian volcanoes, and about explosive volcanism in general, in relation to which the experts were not expected to know precisely the true values but could be expected to provide meaningful credible ranges to capture them. The overall statistical accuracy and information bandwidth of individual's distribution provided the empirical performance basis for differential weighting of their judgments on the target items.
After this first session, a document was provided to participants that summarized, anonymously, all their responses to the target items, together with the preliminary vent opening probability maps that stemmed from combining their judgments with performance weights, and notes about some ambiguities present in the group's responses. Equal weights solutions to target items were also reported, to give context to the pooled findings. Following further group discussion, and in order to generate a final version of weights to be assigned to variables/data sets, a slightly revised questionnaire was prepared and then sent out to participants to complete in a follow-up, remote elicitation. When this second elicitation was completed, the preliminary probability maps were amended accordingly.
In general, during the performance-based elicitation procedure, statistical accuracy (e.g., calibration) and informativeness scores are derived for each expert from a set of subject matter "seed items" [see, e.g., Cooke, 1991;Aspinall, 2006]. These seed items comprise factual questions the true values of which an expert is not expected to know precisely, but an expert is expected to be able to provide meaningful credible intervals that capture those values reliably and informatively, by informed reasoning. Each expert's accuracy and information scores are combined to produce a performance-based weight for application to their responses -as one member of the group of experts-to questions on variables for which estimates are needed, called "target items." In this study, usual practice was followed, and individual expert's judgments were elicited on 5th percentile, median, and 95th percentile values for seed and target questions, in this way obtaining elemental uncertainty distribution markers for each variable. In this form, target item responses are pooled together with the experts' calibration weights to produce a group-synthesized uncertainty distribution, often called a DM solution.
In this work three alternative expert weighting schemes were applied: the CM [Cooke, 1991], the ERF model [Flandoli et al., 2011], and the basic EW model, and then the different DM results obtained were compared. The general purpose of this approach is to provide robust results by comparing outcomes from selective pooling methods (i.e., the Classical Model DM) with more inclusive methods where more of the experts are combined with modulated weights (as in the ERF method) or all experts are given the same weights (i.e., the EW model). More specifically, the CM and ERF methods involve adopting different performance scoring rules and different pooling algorithms, depending on whether robust statistical quantification of uncertainty is the goal (in which case, the Classical Model is favored), or accuracy in pointwise (mean) value estimation is desired (for a discussion of these aspects, see Flandoli et al. [2011] and Bevilacqua [2016]).
Here, the Classical Model has been used for the computation of the final probabilistic vent opening location maps.
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013860