The only two nodes in Multihaz that have no parents are RI and PV (Figure 2). Although they are defined in terms of qualitative states, they can also be described in a quantitative way by defining thresholds that characterize their states (see section “Multihaz applied to Somma-Vesuvius” and Supplementary Material). These thresholds can be general but the probabilities for each state within the node need to be set according to the volcano studied and the potential variability of RI and PV. As mentioned above, the rainfall regime can change significantly depending on regional climate (e.g., Pierson et al., 1992). Similarly, the spatial distribution and total amount of pyroclastic material available at the different catchments surrounding the volcano will depend on factors, such as (a) aleatory variability in eruption behavior: eruption column height, grain size distribution, PDC volume and mobility; (b) aleatory variability in wind vector at different altitudes; (c) topography of the volcano and surrounding terrain. We quantify all these aspects in our application for Somma-Vesuvius in section “Multihaz Applied to Somma-Vesuvius” and Supplementary Material.

Lahar initiation is described in Multihaz through the parameterization of the CPTs for two nodes: RE and IV (see Figure 2 and Table 1). The former has RI, WR and RD as its parent nodes while both the WR and RD nodes have CPTs derived from their common parent node, RI. Finally, the CPT of the IV node is described as a function of the states of its parent nodes, RE and PV (Figure 2).

The CPTs are parameterized following a general-purpose approach, based on diverse strands of evidence about the physical processes behind the complex interrelationships that exist between the variables modeled by Multihaz. This evidence is extracted from the literature and it is summarized in Table 2. The majority of data is related to the RE node which has the largest CPT within Multihaz (Table 1).

In Figure 4, we show some probabilistic trends that are observed within the CPT of the RE node and confirm that the evidence listed in Table 2 is correctly incorporated into Multihaz. For example, the probability of lahar occurrence, i.e., P(RE > 0), is very similar for long-lasting but low-intensity rainfall events (RI = low, RD = high) and for short-lived but high-intensity rainfall events (RI = high, RD = low), especially when the predominant water-routing mechanism is infiltration or a partition between infiltration and runoff (Figure 4D). The CPT of the RE node also shows the expected behavior when runoff is the dominant water-routing mechanism, with increasing probability of lahars occurring as the rainfall intensity increases (Figure 4D), and the lahars that are initiated are expected to have medium-high volumes of material (Figure 4B; Thouret and Lavigne, 2000; Pierson and Major, 2014). Additionally, the importance of rainfall duration in controlling lahar occurrence and volume when shallow landsliding is the main triggering mechanism (Iverson, 2000; Mead et al., 2016) is observed as peaks in probability in all the graphs (and above all for medium-large remobilization efficiencies) when WR = both and RD = high.

The probabilistic patterns included in the CPT of the IV node (derived from the information presented in Table 2) show that the probability of having lahars, especially of medium-large sizes, increases more or less linearly with increasing RE and PV. However, we ensure that the probability of any size of lahars increases considerably when the PV available is relatively high (e.g., Thouret and Lavigne, 2000). Thus, for example, the probability of having medium or large lahars given that PV = high but RE = low is slightly greater than the probability of having such lahars given that PV = low but RE = high. Likewise, the probability of medium-large lahars given that PV = high and RE = medium is higher than the same probability given that PV = medium and RE = high.

Somma-Vesuvius explosive eruptions capable of producing enough pyroclastic material, dispersed along its flanks and over medial-distal surrounding areas, so as to generate considerable lahar activity, range from violent Strombolian up to sub-Plinian and Plinian eruptions (Cioni et al., 1999, 2008; Sulpizio et al., 2006). The natural variability of these eruptions, in terms of tephra-fallout dispersal and PDC inundation, is quite large (Cioni et al., 2008; Gurioli et al., 2010; Sandri et al., 2016; Tierz et al., 2016a). In this paper, we quantify this aleatory uncertainty in pyroclastic volume over different catchments around Somma-Vesuvius (Figure 5), by computing probabilistic hazard curves (e.g., Rougier et al., 2013) of this variable for tephra fallout and dense PDCs. Our hazard assessment can be expressed as either conditional on the occurrence of a small, medium or large eruption (the sizes defined as in Sandri et al., 2016; Tierz et al., 2016a), or accounting for the aleatory uncertainty in the eruption size (e.g., Sandri et al., 2016).

For tephra fallout, we derive the hazard curves of pyroclastic volume at each catchment from data in Sandri et al. (2016). For dense PDCs, we extend the procedure presented by Tierz et al. (2014) as follows. First, we apply a ratio of 0.5 between the maximum flow depth simulated by Titan2D and the potential thickness of dense PDCs, similarly to findings by Charbonnier and Gertisser (2012) and Ogburn (2015). Secondly, we account for enhanced deposition on low-angle slopes in comparison to high-angle slopes (Doyle et al., 2010; Gurioli et al., 2010) by applying a correction similar to Doyle et al. (2010). However, instead of setting up a slope threshold above which there is no PDC deposition, we use the minimum non-zero value of slope over the study area such that the final thickness is a local-slope-weighted fraction of the initial thickness (the greater the local slope of a point, the smaller the fraction), with (2)Hf=H0·min(δ)​/δi where H_f is the final PDC-deposit thickness estimate, H₀ is the thickness estimate after applying the first correction above, δ_i is the (non-zero) slope at a given point, and min(δ) is the minimum (non-zero) slope in the study area. We note that our PDC modeling is restricted to dense PDCs so, given that dilute PDCs are not simulated, a certain thickness of PDC deposit is expected to be missing from our probabilistic quantification of the PDC volume. Nevertheless, dilute PDCs at Somma-Vesuvius may predominantly invade medial and distal sectors away from the central crater (e.g., Gurioli et al., 2010). These sectors are typically located over the Campanian Plain where conditions that favor lahar generation (e.g., steep slopes) are not likely to be met.

Finally, for those catchments in the proximal areas of Somma-Vesuvius where both tephra-fallout and dense-PDC deposition are expected to occur, we randomly sample the hazard curves of tephra fallout and dense PDCs and sum up the samples to obtain the final quantification of the aleatory uncertainty for the variable PV (i.e., hazard curves of PV). This step assumes that tephra-fallout and PDC propagation are independent, although volumes are correlated through the eruption sizes. Using the hazard curves for PV, we can parameterize the prior table of the node, at each catchment. We set two (general-purpose) thresholds: PV ≤ 10⁴ m³ for PV = low; PV > 10⁶ m³ for PV = high; and 10⁴ < PV ≤ 10⁶ m³ for PV = medium.

Figure 6 shows some examples of the hazard curves of tephra-fallout and dense-PDC volumes at catchments on the flanks of Somma-Vesuvius. We observe quite different curves depending on (a) the location of the catchment with respect to the central crater (we remark that the predominant wind in the area, and therefore the tephra-fallout transport, is generally toward the east) and the topographic barrier of Mount Somma, and (b) the size of the eruption considered. Generally, catchments located beyond Mount Somma (Figures 6A–C) do not experience significant dense-PDC deposition during small eruptions (i.e., pyroclastic volume is < 10² m³). On the central and oriental sectors beyond Mount Somma (Figures 6A,B), greater accumulations of pyroclastic material are expected to occur due primarily to tephra fallout rather than dense PDCs. At the westernmost end of Mount Somma however (Figure 6C), we see that greater pyroclastic volumes can result from dense-PDC deposition rather than from tephra fallout. On the southern flank (Figure 6D), small-size dense PDCs can deposit a substantial amount of pyroclastic volume (even though these volumes are lower than those produced by tephra fallout) and medium and large dense PDCs have the potential to deposit considerably greater amounts of pyroclastic volume in comparison with tephra fallout. We also notice the influence of the eruptions of different sizes on the hazard curves of eruptions of any size. For instance, the hazard curve for dense PDCs of any size in catchment 331 (black dashed line, Figure 6D) shows a step-wise shape. This indicates that (1) small-size eruptions have a high weight (i.e., higher probability of occurrence, Marzocchi et al., 2004; Sandri et al., 2016) in defining the any-size hazard curve since the latter follows the small-size curve over its domain, and (2) the gap between the end of the small-size hazard curve and the values of medium-sized dense PDCs with a exceedance probability far from unity results in the plateau observed in the hazard curve for an eruption of any size (Figure 6D).

In Figure 7, we provide an example of the kind of forecast that Multihaz can produce. In particular, we select three medial-distal catchments to highlight the changes in these predictions as a function of (1) the location of the catchment relative to the main dispersion axis for tephra fallout, and (2) the area of the catchment. The selected catchments are: #195 (Figure 7A; area ~2 km², NNE of the vent), #225 (Figure 7B; area ~ 10 km², NNE of the vent), and #277 (Figure 7C; area ~2 km², ENE of the vent). In this example, we show Multihaz models built in the event of a large-size eruption occurring. However, the eruption size has a strong influence on the Multihaz behavior, as we show later in this section for all the catchments and as demonstrated in the previous section for the PV variable over proximal catchments (Figure 6).

The location of the catchment has a strong role on the potential lahar hazard, as the prevailing wind direction leads to greater tephra-fallout accumulation to the E and SE rather than NNE of Somma-Vesuvius (Figure 7). For example, the probability of large volumes of pyroclastic material being stored in catchment #277 (Figure 7C) is more than 30% higher than this probability at catchment #195 (Figure 7A). The area of the catchment can also be important; if we compare a small catchment situated on the tephra dispersion axis (catchment #277, Figure 7C) with a larger catchment situated further from the dispersion axis, but still downwind from the volcanic vent (catchment #225, Figure 7B), we observe that the probability of PV = low is similar between the two catchments. However, catchment #277, on the tephra dispersion axis, has a probability of PV = high almost 25% higher (Figure 7C), illustrating the importance of location for tephra accumulation in medial-distal catchments.

If we instantiate (i.e., we set one state of the node to be true) the RI node to “high” and the WR node to “both”, we find that the probabilities at the RD and RE nodes mimic each other for all three catchments (Figure 7), since they all share the same parameterization of the CPTs for RD (as determined from regional meteorological behavior and inversely related to RI) and RE (since we have assumed that the tephra properties are consistent across the catchments). However, in terms of initiation volume, we notice that there are important differences in the probability of small-, medium-, and large-sized lahars occurring. In the case of catchments #225 and #277, the differences due to catchment size are offset by the location of the catchments relative to the primary tephra dispersion axis. Nevertheless, we find that the probability of having large lahars is still around 3% greater for the smaller catchment #277 that is located on the on dispersion-axis (Figures 7B,C).

In this study, we parameterize 1092 Multihaz models (273 catchments × 4 eruption sizes that may precede the lahar triggering) that give a wide spectrum of probabilistic hazard assessments for lahars of specific volumes occurring at Somma-Vesuvius. However, exploring all these models (and their coupling with LaharFlow) would be extremely costly in computational terms. We focus our preliminary application of the proposed probabilistic multi-hazard framework on three specific scenarios designed to allow assess the roles of (i) the size of an eruption occurring at Somma-Vesuvius, (ii) the specific value of RE triggering the rainfall lahars, and (iii) the maximum value of IV that a specific catchment can produce, given (i) and (ii). The multi-hazard scenarios are:

Scenario 1: a large eruption occurs at Somma-Vesuvius and lahars are triggered by means of high remobilization efficiencies (i.e., RE = high).

We recall that Multihaz for each catchment shares a common parameterization of the WR, RD, RE, and IV nodes but it has its own parameterization for the RI and PV nodes (see section “Multihaz Applied to Somma-Vesuvius” and Supplementary Material). Figure 8 shows the probabilities for the different states of the IV node, for each of the explored scenarios. To aid in visualizing the results, we divide the catchments into four zones (see Figure 8B): (1) medial-distal catchments (~20–30 km from the vent) toward the SE and E; (2) medial-distal catchments toward the NE; (3) medial catchments (~15–20 km from the vent) toward the E; and (4) proximal catchments (within a radius of about 3–4 km from the vent).

Scenario 1 (Figure 8A) is characterized by an almost constant probability for the occurrence of medium-size lahars in each catchment while the probabilities for small and large lahars are both smaller and tend to mirror each other (the states of each node are mutually exclusive and exhaustive, so the sum of their probabilities must be 1, and P(IV = zero | RE ≠ zero) = 0 in our parameterization of Multihaz). In Figure 8D, we display the spatial distribution of P(IV = high) under scenario 1. The northwesternmost catchments of zone 2 have the lowest probabilities of initiating large-volume lahars. In zone 4, however, the probability of small-size lahars is reduced and it is almost zero for the catchments located on the southern flank of Somma-Vesuvius (Figure 8A). In these catchments, IV = high is the most likely state (with probabilities above 50%) under scenario 1 (Figures 8A,D).

Scenario 2 (Figure 8B) is characterized by small and medium lahars having relatively similar probabilities of occurrence in zones 1 and 3. The probability of large lahars increases slightly over zone 3 and more significantly over zone 4, similarly to that observed in scenario 1 (Figures 8A,B). Nonetheless, large lahars are never the most likely outcome in scenario 2.

Scenario 3 is strongly dominated by all catchments having IV = low as their most likely state (Figure 8C). The probability of having medium lahars increases over zone 4 coinciding with a very slight increase in the probability of large lahars, up to 4–5% (Figure 8C). The jump in probability of occurrence for small and medium lahars inside zone 4 marks the separation between the catchments located on the north flank of the volcano, where small lahars are much more likely in scenario 3, and those on the south flank, where small lahars are only slightly more probable than medium lahars (Figure 8C).

The results shown above can be combined with simulations of lahar dynamics to examine the possible impact of such lahars in probabilistic terms. For each of the scenarios, we perform LaharFlow simulations for each of 24 simulation groups each containing multiple lahar sources in a subset of the 273 catchments (see Supplementary Material). We define the hazard domain for each simulation group as a grid of 100 × 100 m of spatial resolution that extends from the source areas of lahars through the areas over which they propagate. The simulations allow us to compute, at each grid point of the hazard domain, three values of lahar flow depth and speed and three values of probability associated with such an event, according to the Multihaz assessments. To do this, we must compute the probability of each source within each simulation group. We assume the sources in each catchment are independent and derive these probability values by multiplying the probability of each catchment to generate the volume of lahars simulated with LaharFlow, given the scenario: (3)pi,j=∏kNjPk(IVmax|scei) where p_i,j is the probability assigned to the hazard footprint computed with LaharFlow for scenario i and simulation group j; N_j is the number of catchments that form the simulation group j; P_k(IV_max|sce_i) is the probability of the catchment k to trigger lahars with the maximum possible volume for the catchment (IV_max, the volume simulated with LaharFlow), given the scenario i (sce_i). These P_k values come from the Multihaz model implemented in each catchment and for an eruption of the size indicated by a specific scenario.

The probabilistic independence between the P_k values holds if we assume that: (1) the triggering rainfall is homogeneous over the catchments of the same simulation group, (2) the catchments of the same simulation group are located close enough to each other that changes in available tephra volume (i.e., PV) are only due to the catchment area and not to the spatial distribution of PV, and (3) the spatial distribution of the catchment area does not show any clustering pattern, that is: small and large catchments do not tend to cluster together, spatially. We recognize that some of these assumptions may not hold in some cases. For example, assumption (2) might be reasonable for tephra fallout but it may be questionable for dense PDCs in the presence of topographic barriers, such as Mount Somma. Relaxing these assumptions would require the development of another probabilistic model to assess the dependencies among the P_k values for each simulation group and this is beyond the scope of this paper.

The lahar flows simulated at Somma-Vesuvius and surroundings tend to be channelized following the main (narrow) valleys that descend from their initiation points but there are also cases in which steep short valleys lead to lahars that converge and accumulate material over larger flat areas, like in the area near Sarno. We summarize the observed patterns in lahar propagation in Figure 9. We select four simulation groups, which are characterized by (a) including catchments with large areas, (b) representing one of the zones defined in section “Multihaz Probabilistic Assessments,” and (c) having densely-populated areas nearby. Thus, the simulation groups are the following:

The total number of catchments belonging to these four simulation groups is 62, almost one quarter of the total number of catchments (273). In Figure 9, we show the hazard footprints of the maximum flow-depth of lahars computed for scenario 1, for which the most voluminous lahars are produced.

In simulation group 3 (Figure 9A), the flows propagate following the principal valleys of the catchments from SE to NW, primarily, and they spread out after reaching the flat area where Castellammare di Stabia is located. The longest runout of lahars in the simulation group is about 7–8 km (Figure 9A). The maximum values of flow depth tend to be around 10 m while the maximum values of speed are >25 m/s, which are attained close to the initiation of some lahars, but are around 5–10 m/s over the flat lowlands.

In simulation group 19–20 (Figure 9B), the flows propagate following the principal valleys in some cases, especially on the eastern and southeastern sectors of the hazard domain (Figure 9B). The largest flow, formed over the southeastern sector, can develop cross-section widths of few hundreds of meters while keeping flow depths of almost 10 m (Figure 9B). The flow crosses the area where Mugnano di Cardinale, Avella, and other municipalities are located and stops after ~15 km of propagation. The maximum speeds recorded in the simulation can be 30 m/s at the initiation areas but they decrease to about 5 m/s over the flat areas.

In simulation group 22 (Figure 9C), the flows propagate short distances along quite straight valleys before reaching the flat area on the Campanian Plain where Episcopio and Sarno are located. The lahars in the eastern part of the hazard domain converge over the Campanian Plain to produce an extended inundated area where flow depths can be close to 10 m (Figure 9C). The maximum run-outs of the whole simulation group are not greater than 3 km. In terms of lahar speed, maximum values can be above 30 m/s near to the initiation points, are sustained around 10–15 m/s over the initial part (~1 km) of the Campanian Plain and then gradually decrease to below 5 m/s.

In simulation group 24 (Figure 9D), the flows are triggered from the upper part of the Somma-Vesuvius edifice, including the north flank beyond Mt Somma, and propagate radially developing braided patterns when they have reached the mid slopes as the topography is not as confining as in other simulation groups (Figure 9D). The majority of the lahars simulated on the south flank reach the sea after having traveled for ~5–6 km. The lahars on the other flanks of the volcano have similar maximum runouts (Figure 9D). The maximum flow depths are around 10 m in some proximal sectors where flows might get channelized. Maximum flow speeds can reach 25 m/s and be maintained around 10-12 m/s over proximal-medial sectors (up to ~4 km). On distal sectors (~4–6 km), the lahars slow down to speeds of 5 m/s or lower (Figures 9D,10A).

Combining the probabilistic assessments performed by Multihaz with the hazard footprints of lahars computed with LaharFlow, we are able to extract discrete probabilistic output distributions for two hazard variables of particular interest: lahar flow depth and speed. In Figure 10, we show such a probabilistic quantification of the lahar hazard for the flow speed and simulation group 24 (Somma-Vesuvius flanks). We compute values of the hazard variable (in this example the lahar speed) for each scenario and our BBN provides the likelihood of each initial volume of lahars which, with additional assumptions about probabilistic independence, allows us to compute the probability value relative to each lahar speed (p_i,j in Equation 3; we apply the formula for the catchments generating lahars that actually interact before impacting three selected locations; red filled circles in Figure 10). The final result are discrete probabilistic output distributions of the lahar speed at each location (Figure 10D). Note that we can select any point within our hazard domain and visualize its output probability distribution for lahar speed, and also for lahar flow-depth or any other variable that could be extracted from the simulations.

We find that the fastest lahars are expected to occur for flows simulated under scenario 1, with the maximum speed recorded at San Sebastiano al Vesuvio (~6–7 m/s; p_{1, 24} ~ 12%). Scenario 3 is associated with lahars of maximum speeds between 0 and 2 m/s for the three selected points and probabilities always below 6% (Figure 10D). Scenario 2 is considered the most likely of the three scenarios (p_{2, 24} ~[3–18]%). The lahars simulated here show speeds around 2 m/s at Torre Annunziata (p_{2, 24} ~ 3%), 5 m/s at San Sebastiano al Vesuvio (p_{2, 24} ~ 18%), and 0.3 m/s at Somma Vesuviana (p_{2, 24} ~ 4%). Interestingly, we see that the range of lahar speed at Somma Vesuviana is narrower compared to the ranges at the other two locations. Thus, the maximum lahar speed at Torre Annunziata for scenario 1 is slightly higher than that at Somma Vesuviana but, for scenario 3, the former is 0 m/s and the latter is around 0.3 m/s. In summary, we observe that relatively slow lahars may be expected to impact Somma Vesuviana (e.g., scenarios 2 and 3) while lahars may reach Torre Annunziata only during medium-large events (e.g., scenarios 1 and 2) but these lahars could have higher speeds than those impacting areas northwards of Mount Somma (Figure 10).

In order to model the rainfall-triggering of lahars, we develop a simplified but very versatile probabilistic model that can be easily adapted to many volcanoes, provided that there are probabilistic descriptions available at least for the regional rainfall intensity around the volcano and the spatial distribution of tephra deposits around the volcano following an eruption. Another requirement for building Multihaz at other volcanoes may be the calibration of the parameters of the WR node in order to model tephra deposits with variable properties, e.g., infiltration capacity related to the median grain size distribution (e.g., Pierson and Major, 2014; Jones et al., 2015).

The triggering process is an important ingredient of a lahar hazard model, as it controls the location and timing of lahar events, and influences the magnitude of the flow. However, the physical processes leading to initiation remain poorly understood. While some studies consider erosive runoff to be the principal initiation mechanism (e.g., Collins and Dunne, 1986; Lavigne et al., 2000a; Pierson and Major, 2014; Jones et al., 2015) others deem that shallow-landsliding is the primary generation mechanism for lahars (e.g., Iverson, 2000; Volentik et al., 2009; Mead et al., 2016). Discrepancies in interpretation of volcanological phenomena can lead to contrasting model assessments (e.g., Hincks et al., 2014) so it is important that the selected model is able to incorporate the range of possible behaviors. It is likely that both initiation processes (and perhaps others) can generate lahars and further research could elucidate the regimes in which each is active. Until a deeper understanding of the initiation processes is available, models should capture this source of epistemic uncertainty. Multihaz attempts to model the two distinct initiation mechanisms for lahars through the parameterization of the RE node as a function of its parents: RI, RD and, especially, WR (see Tables 1, 2; and Figures 2–4). Nonetheless, Multihaz could be adapted to model only erosive runoff or shallow-landsliding, for instance, by permanently instantiating the WR node to “runoff” or “both,” respectively (although WR = ifl also implies the possibility of shallow landslides triggering lahars; see Figures 2, 3). Alternatively, the structure of Multihaz could be changed to reflect new fundamental physical insights into lahar initiation, for example removing the WR node and adding other variables that are thought to directly control the triggering for each initiation mechanism, e.g., hydraulic diffusivity (e.g., Mead et al., 2016).

Multihaz is a first attempt to probabilistically model (and forecast) lahar initiation following explosive volcanic eruptions using an integrated multi-hazard framework, and we consciously construct a relatively simple BBN. Inevitably, this leads to some limitations of our model. Thus, the parameterization of Multihaz is preliminary and the performance of the model is not tested against real data. However, at certain volcanoes (e.g., Mount Merapi, Indonesia; Lavigne et al., 2000a; De Bélizal et al., 2013), large datasets of lahar observations may be used to either calibrate the Multihaz parameters or test the performance of the network (Figure 11). That is, using combinations of node states that have been observed (e.g., RI = high, RD = medium, IV = medium), it is possible (i) to refine/update, through Bayesian inference, the probability values (i.e., parameters) in Multihaz; or (ii) to calculate the likelihood of observing the aforementioned combinations, given our paramerization of Multihaz (e.g., Murphy, 2001; Koller and Friedman, 2009). Furthermore, a variety of observations could be collected to calibrate or test Multihaz at any volcano. For instance (1) direct observations of the water-routing mechanism, given rainfall intensity, through slope plots (e.g., Chinen, 1986; Leavesley et al., 1989); (2) discrimination between event/no event, given rainfall intensity and duration, through acoustic flow monitoring (e.g., Lavigne et al., 2000b; Jones et al., 2015); (3) estimation of initiation volume and/or remobilization efficiency through field analyses and remote-sensing techniques (e.g., Bremer and Sass, 2012; Amici et al., 2013; Pierson et al., 2013).

Even though the Multihaz probabilistic assessments for the initiation volume are an advance on previous studies in which the outputs are binary (e.g., event or no event, as in Van Westen and Daag, 2005; Berti et al., 2012; Jones et al., 2015), our BBN model is unable to give a detailed temporal analysis of the triggering (e.g., Iverson, 2000; Mead et al., 2016) or even of the spatial changes in lahar susceptibility with time (e.g., Fan et al., 2017). Further development of our multi-hazard framework could include the implementation of dynamic Bayesian networks (e.g., Koller and Friedman, 2009) to model the temporal aspects of the lahar initiation. Additionally, LaharFlow is capable of modeling sources that are active at different times so such a dynamic BBN could be coupled to the flow simulator to provide spatio-temporal hazard footprints of lahars.

Secondary hazardous volcanic processes, such as lahars, represent a further challenge for hazard assessment in comparison with primary processes, such as tephra fallout and PDCs. This is largely due to the difficulty in quantifying the cascade effect among the processes. That is, the aleatory (and epistemic) uncertainty in the primary processes conditions the spatial and temporal variability in the secondary process (e.g., Pierson and Major, 2014). Our proposed multi-hazard framework provides a viable solution to deal with such situation.

Firstly, we believe that constraining the spatial distribution of tephra deposits in probabilistic terms and incorporating this information into Multihaz is one of the major strengths of our multi-hazard framework. The coupling between the areas that are expected to be affected by tephra deposition, given an eruption, and the areas that may subsequently act as lahar sources is indeed rarely seen in the literature (e.g., Volentik et al., 2009). However, there are strong links observed between the spatial distribution of tephra accumulation and the magnitude and frequency of lahars at several volcanoes after recent eruptions (e.g., Pierson et al., 1992; Pierson and Major, 2014). In addition, the importance of quantifying the aleatory uncertainty in tephra accumulation is highlighted by our results, which account for the possible occurrence of unlikely (but still possible) events (e.g., large-size lahars at catchment #195 after a large eruption; see Figure 7A). Indeed, low-probability, high-consequence situations are often very relevant not only in volcanic hazards but in any kind of natural hazard (e.g., Rougier et al., 2013; Mignan et al., 2014).

Secondly, our multi-hazard framework fundamentally relies on the strategy of separating the calculation of the probability of a specific event (e.g., medium-size lahars at catchment #195, given the occurrence of a large eruption at Somma-Vesuvius; see Figure 7) and the hazard footprint that ensues from this specific event (e.g., Spiller et al., 2014). This is a great advantage because it implies that changes in the probability distributions of Multihaz, for instance due to a modification in the prior table of RI, do not inevitably require many more simulations to be performed.

The final probabilistic products that we obtain, after having run the lahar simulator, are discrete output distributions for lahar speed (Figure 10) or flow depth. These are important outputs since they link the hazard-intensity measure and its probability of occurrence at a given point of the hazard domain, in a similar way to hazard curves. However, our probabilistic output distributions still consider too few “realizations” of the hazard (they are based on 3 values, one per scenario) and, therefore, we do not have enough information about the exceedance probabilities over the whole possible domain of the hazard-intensity measure. This limitation could be overcome by building a BBN model that used continuous PDFs at each of its nodes (e.g., Hanea et al., 2006) as well as by propagating uncertainty through performing a larger number of lahar simulations and/or using sophisticated uncertainty quantification techniques (e.g., Spiller et al., 2014).

This study presents an example of lahar hazard assessment at Somma-Vesuvius that includes novel aspects: (1) an explicit and quantitative coupling between lahar triggering and propagation, and (2) a probabilistic assessment of the lahar initiation volumes. Previous studies have mostly dealt with the analysis of the stratigraphical sequences and historical data of syn- and inter-eruptive lahars (e.g., Rosi et al., 1993; Lirer et al., 2001; Sulpizio et al., 2006). We note that if a precise reconstruction of the spatial distribution of some syn-eruptive lahars was available, this kind of data could represent a first step to test the reliability of our integrated framework, for instance by comparing these lahar deposits to our LaharFlow simulations (e.g., Tarquini and Favalli, 2011; Tierz et al., 2016b).

Other hazard studies of water-sediment flows at Somma-Vesuvius have focused on the detailed description of the susceptibility of different catchments (both around the volcano and in very distal sectors) to act as sources for these volcaniclastic flows according to the spatial distribution of past pyroclastic deposits and the hydrogeomorphological characteristics of these catchments (e.g., Bisson et al., 2010, 2014). Even though our maps are not directly comparable (e.g., our hazard metric—lahar initiation volume—is different from those of Bisson et al., 2010, 2014—instability or disruption proneness), we note that the spatial distribution of the most hazardous catchments is roughly similar. However, with respect to the analyses presented by Bisson et al. (2010, 2014), our study crucially incorporates the aleatory uncertainty in RE and PV. Thus, some catchments ENE of Somma-Vesuvius, which according to morphological features only may not be classified as highly hazardous (e.g., Figure 8 in Bisson et al., 2014), do seem to have the potential to trigger large-volume lahars in our study (Figure 8D) because they are located downwind the main tephra-dispersion axis (Figure 7).

In terms of lahar propagation at Somma-Vesuvius, Favalli et al. (2006) presented a study focusing on the “maximum expected (eruptive) event” at the volcano which, according to the emergency plan implemented by the Italian civil protection at that time, was a sub-Plinian I eruption equivalent to the 1631 AD eruption (e.g., Rosi et al., 1993; Cioni et al., 2008). In contrast, our model includes the aleatory uncertainty in tephra-fallout dispersion and dense-PDC deposition. Nevertheless, our modeling of tephra dispersal does not explicitly consider the occurrence of long-lasting violent Strombolian/ash-emission eruptions (Cioni et al., 2008). This type of eruption could produce significant amounts of tephra near the Somma-Vesuvius cone for long periods of time and, therefore, lahar production might be promoted in the mid to long term (e.g., Jones et al., 2015).

In contrast to Favalli et al. (2006), our results do not identify the Acerra-Nola Plain (Figure 5) as the most-hazardous area in terms of rain-triggered lahars. Due to the expected spatial distribution of the tephra-fallout deposits (Figure 7), the closest catchments to the central part of the plain (approximately catchments #200–225 in Figure 8, zone 3) show the lowest probability, among all catchments, of producing large lahars (Figures 8A,D). Therefore, even when considering an extreme scenario (large eruption combined with high remobilization efficiencies: scenario 1), the LaharFlow runs, for instance, of the simulation group 19-20 (Figure 9B), do inundate the eastern part of the Acerrra-Nola Plain but do not propagate as far as the flows simulated by Favalli et al. (2006). We argue that this can be also related to the quite high value of rainfall intensity used by Favalli et al. (2006): peak intensity at 108 mm/h during 0.5 h of rainfall. According to the hydrological report presented by AdBCC (2015), the maximum rainfall intensity for an infinitesimally-short rainfall (i.e., RD → 0) at the area of the Acerra-Nola Plain, is just above 85 mm/h (we remind that our parameterization of RI through yearly maxima does not influence the hazard analysis of scenario 1 because IV is independent of RI, given RE; see Figure 2).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.