Lahar events in the last 2000 years from Vesuvius eruptions – Part 3: Hazard assessment over the Campanian Plain

In this study we present a novel general methodology for probabilistic volcanic hazard assessment (PVHA) for lahars. We apply the methodology to perform a probabilistic assessment in the Campanian Plain (southern Italy), focusing on syn-eruptive lahars from a reference size eruption from Somma–Vesuvius. We take advantage of new field data relative to volcaniclastic flow deposits in the target region (Di Vito et al., 2024b) and recent improvements in modelling lahars (de' Michieli Vitturi et al., 2024). The former allowed defining proper probability density functions for the parameters related to the flow initial conditions, and the latter allowed computationally faster model runs. In this way, we are able to explore the effects of uncertainty in the initial flow conditions on the invasion of lahars in the target area by sampling coherent sets of values for the input model parameters and performing a large number of simulations. We also account for the uncertainty in the position of lahar generation by running the analysis on 11 different catchments threatening the Campanian Plain. The post-processing of the simulation outputs led to the production of hazard curves for the maximum flow thickness reached on a grid of points covering the Campanian Plain. By cutting the hazard curves at selected threshold values, we produce a portfolio of hazard maps and probability maps for the maximum flow thickness. We also produce hazard surface and probability maps for the simultaneous exceeding of pairs of thresholds in flow thickness and dynamic pressure. The latter hazard products represent, on one hand, a novel product in PVHA for lahars and, on the other hand, a useful means of impact assessment by assigning a probability to the occurrence of lahars that simultaneously have a relevant flow thickness and large dynamic pressure.

1 Introduction

Lahars are flows of water and entrained sediments that originate from the remobilization of volcaniclastic deposits by water, either from rain, melting ice or snow, or a sudden release by a crater lake (see the companion paper by de' Michieli Vitturi et al., 2024, and references therein). They represent one of the processes causing the highest death toll among volcanic phenomena. According to the analysis by Auker et al. (2013), syn-eruptive lahars are responsible for about 14 % of the fatalities in their database, whereas post-eruptive or inter-eruptive lahars cause another 3 %. Among the most tragic episodes of lahar impact, we recall the lahar generated from Nevado del Ruiz, which buried the city of Armero, killing about 23 000 people, making it the second-worst volcanic disaster of the 20th century (Pierson et al., 1990; Voight et al., 2013; Parra and Cepeda, 1990). Other examples include the series of lahars that hit the surroundings of the Pinatubo volcano in the years after the 1991 eruption (Pierson et al., 1996; Umbal and Rodolfo, 1996; Rodolfo et al., 1996), the Tangiwai disaster at Ruapehu (Manville, 2004) and the lahars from the 2008–2009 Chaiten eruption (Pierson et al., 2013).

A simplified method commonly used so far to describe lahar impact is the LAHARZ model (Schilling, 1998). It is based on the statistical correlation between the inundated area and the mass flow volume inferred from past events. However, in recent years, a few examples of probabilistic hazard assessment for lahars based on more robust statistical treatments, like statistical surrogates or emulation approaches, have been proposed for different volcanoes worldwide, such as Mead and Magill (2017) on Ruapehu (New Zealand), Tierz et al. (2017) on Vesuvius (Italy), and Gattuso et al. (2021) on Vulcano (Italy).

Hazard assessment for lahars needs to consider (i) the identification of potential source regions for volcanic material and for water, including snowcaps and glaciers; (ii) the potential magnitude and characteristics of the flow; (iii) the topography between the source region and the potential targets at risk; (iv) the potential for modification of the flow properties along the path; (v) the frequency of such events in the past; and (vi) the meteorological data at the source region and along the potential path of such flows, especially for extreme events.

As explained in the companion paper by de' Michieli Vitturi et al. (2024), lahars can change character downstream through processes of flow bulking and debulking, generating high variability in both time (i.e. unsteadiness) and space (i.e. non-uniformity) for variables pivotal for hazard assessment, such as particle concentration, granulometry and componentry, bulk rheology, and velocity. The full complexities associated with these processes prevent us from effectively modelling these flows for quantitative hazard assessment purposes. At present, even if we could describe all the underlying physics, a full 3D simulation of all these phenomena would require prohibitive computational costs, and current numerical models describe only some of the observed phenomena or use simplified approaches (see the companion paper by de' Michieli Vitturi et al., 2024). In terms of numerical modelling, a good compromise between the completeness of the physics behind these phenomena and the computational feasibility is represented by the shallow water approach (de' Michieli Vitturi et al., 2024), where model complexity is reduced with a depth averaging of flow properties. This approach approximates the original 3D problem with a 2D model, and it is the one we apply in this work.

In the Campania region, which is largely covered by fallout and pyroclastic density current (PDC) deposits from eruptions of Somma–Vesuvius, Campi Flegrei, and Ischia volcanoes, the signature of several syn- and post-eruptive lahars has been found in the geological record (Di Vito et al., 2024b; Sulpizio et al., 2006; Zanchetta et al., 2004a, b). Furthermore, detailed lists of documented lahars in the 20th century are available in the literature (Fiorillo and Wilson, 2004). Despite such evidence, up to now most of the probabilistic volcanic hazard assessments (PVHAs) for this region have mainly focused on PDCs (e.g. Neri et al., 2008, 2015; Gurioli et al., 2010; Sandri et al., 2018; Tierz et al., 2018) and tephra fallout (e.g. Costa et al., 2009; Selva et al., 2010, 2018; Sandri et al., 2016; Massaro et al., 2023), while systematic quantitative hazard assessments from lahars (see, for example, Jenkins et al., 2022) have been lacking. An exception is provided by Tierz et al. (2017), who applied a Bayesian belief network to assess the effect of different factors (linked to rainfall intensity and volcanoclastic volume) on the probability of different initial volumes of lahars. However, that study did not explore the variability in the hazard assessment related to the initial flow conditions (mostly linked to the flow volume, detachment area, and volumetric solid fraction).

PVHA for lahars requires performing a high number of simulations in order to enable a quantification of the uncertainty linked to model parameters.

Recent technical improvements (e.g. code parallelization) and generalizations (e.g. description of erosion and deposition during the flow) of lahar models, such as that implemented in the most recent version of the IMEX_SfloW2D model described in the companion paper by de' Michieli Vitturi et al. (2024), permit hundreds of simulations from different catchments, with different initial and boundary conditions in reasonable times, necessary for characterizing the intrinsic variability and for the production of hazard maps.

Following several surveying campaigns carried out to characterize lahar deposits in natural exposures, archaeological excavations, and ad hoc trenches in the plain surrounding the Vesuvius edifice and along the Apennine valleys, Di Vito et al. (2024b) present the results of a multidisciplinary study which shows the presence of volcaniclastic deposits (mostly debris and mud flows but also from hyperconcentrated flood flows) even in areas very far from both the Apennine hills and the valleys of Somma–Vesuvius, demonstrating the high mobility of these flows. In particular, Di Vito et al. (2024b) focused on the analysis of the syn- and post-eruptive lahar deposits generated by the two sub-Plinian eruptions of Vesuvius in 472 CE and 1631. Thicknesses, sedimentological features (lithic content, pumice provenance, grain size, boulder entrainment), and vertical and/or lateral continuity of the deposits were reported during the campaigns in order to establish characteristic facies (massive to structured, poorly sorted to better-sorted with respect to the primary pyroclastic deposits and topography) and general flow dynamics (velocity, dynamic pressure, thickness) of those volcaniclastic systems. Results show that the inclusion of fine ash in the whole deposit distribution, the depositional mechanism of the primary pyroclastic deposits (fallout vs. current), and the large-scale topographic effects (plain vs. valley) are the main geological features affecting the size and style of the remobilization that occurred for the two eruptions (Pollena and 1631).

In this work, we take advantage of these new field data analyses and recent improvements in modelling lahar flows (de' Michieli Vitturi et al., 2024) to explore the effect of uncertainty in the flow initial conditions on the invasion of lahars in the Campanian Plain (Fig. 1a) by sampling coherent sets of values for the input model parameters and subsequently performing a considerable number of lahar simulations (1100 in total) needed for the production of hazard maps.

Figure 1 (a) Campanian Plain highlighting the target area of the present study (black rectangle) covered by a 50 m × 50 m grid; the coloured dots respectively show the location of two illustrative points, which are San Marzano sul Sarno (S) and Casilli (C) train station, that will be used as examples in the paper. (b) Shaded topography zoom on the target area highlighting the 11 catchments considered (1–6 and 7–11, respectively, for Vesuvius and Apennine catchments). (c) The dashed red line encompasses the part of the domain where the TIN Italy DEM has been used.

We present a novel general methodology for PVHA for lahar flows but we focus on syn-eruptive lahars, conditional on the occurrence of a reference eruption for Somma–Vesuvius (the medium-magnitude scenario by Cioni et al., 2008; Macedonio et al., 2008; Sandri et al., 2016).

In particular, we account for the following:

• 11 different catchments (Fig. 1b) where lahars could originate and impact the target area of the Campanian Plain from both the Somma–Vesuvius edifice and from the Apennines sectors to the east and south;

• deposits from PDCs (mostly on the Somma–Vesuvius catchments) and from tephra fallout (on the Apennine catchments) from the reference eruption of Somma–Vesuvius;

• the maximum expected rainfall in a few days, taken to be of the order of 500 mm, as extracted from the rainfall record in the last 70 years in the Campania Region (Fiorillo and Wilson, 2004), as well as coherence among the initial values of the flow (initial thickness, detachment area, and volumetric solid fraction), the deposit porosity (water-saturated), and the available water from rainfall and volcaniclastic sediments from the reference eruption. In order to do so, we build up a strategy to sample the input model parameters, which will be illustrated in Sect. 3.1.

The first and last points, in particular, allow us to explore the uncertainty in the position of lahar generation and in the initial flow conditions (in terms of area of detachment, initial volume, volumetric solid fraction).

The goal of the study is to quantify the conditional probability of invasion by at least one lahar originating from the remobilization of tephra deposits due to the reference eruption at Somma–Vesuvius (Sandri et al., 2016). In order to provide useful results for future quantification of the impact of syn-eruptive lahars, we express our results in terms of exceedance probability for selected thresholds of two variables pivotal for hazard assessment, such as maximum flow thickness and dynamic pressure. The domain of interest, covering the Campanian Plain (see Fig. 1a), was discretized for computational reasons on a 50 m × 50 m grid. This resolution represents a good compromise between solution accuracy and computational time required for a simulation (see companion paper, de' Michieli Vitturi et al., 2024). As we mentioned above, we also compute the exceedance probability of selected pairs of critical thresholds in dynamic pressure and flow thickness that jointly are key parameters for lahar impact assessment (e.g. Wilson et al., 2014). In other words, we take into account the flow “history” in every target grid point, and we compute the probability of the flow to simultaneously exceed (in at least one time step of the simulated flow) two values of flow thickness and dynamic pressure, repeating this for several pairs of values.

The paper is organized as follows: first we very briefly summarize the geological information and the features of the model; second, we present the method used for the PVHA, and then we show the results as maps. Finally, we present a discussion and conclusion in order to highlight the main achievements and the current limitations, which will be addressed in future works.

2 Field data and transport model

The code developed by generalizing the IMEX_SfloW2D model for describing lahar flows, described in the companion paper by de' Michieli Vitturi et al. (2024), was calibrated using the field data presented by Di Vito et al. (2024b). The field data were used also to define the available initial deposit from a reference eruption from Somma–Vesuvius in the different catchments from PDC and tephra fallout deposits.

2.1 Remobilizable PDCs and tephra deposits

In this study, for hazard assessment purposes we considered a reference eruption belonging to the medium-magnitude scenario (MMS; Sandri et al., 2016). The initial volumes that can be remobilized as lahars come from PDC and tephra fallout deposits as well as from the available water (rain in our case). PDC deposits are more relevant for Somma–Vesuvius catchments, which are proximal to the source; tephra fall deposits, dispersed by the local wind fields, are mostly relevant for the Apennine catchments, where PDCs from an MMS eruption do not leave any appreciable deposit.

As regards the PDC deposits, we used the field data from the most recent sub-Plinian eruptions, which are the 472 CE (Pollena; Sulpizio et al., 2005) and the 1631 CE (Rosi et al., 1993) eruptions. Cautiously, in each grid cell of a given catchment we considered the maximum thickness between these two PDC events to be available PDC deposits, as mapped by Gurioli et al. (2010).

As regards tephra fall deposits, we rely on the results of the simulations presented in Sandri et al. (2016), where 1000 simulations were performed for the MMS considering variability of the meteorological conditions and ESPs (total erupted mass, mass eruption rate, total grain size distribution). Specifically, we randomly sampled (without repetition) 100 fallout deposits from the 1000 available. Hence we used those deposit distributions for the simulations performed with the generalized IMEX_SfloW2D model.

2.2 Lahar runout

A considerable set of lahar runout estimations inferred from the deposits in the Campanian Plain (see Fig. 2), associated with the 472 and 1631 eruptions (Di Vito et al., 2024b), were used to calibrate the empirical parameters needed for friction, erosion, and deposition terms (de' Michieli Vitturi et al., 2024).

Figure 2 The dots highlight the locations of field data analysed in Di Vito et al. (2024b), where measurements and facies analyses of lahar deposit features were taken or retrieved by an inverse engineering approach.

2.3 Grain size distribution

The grain size distribution (GSD) is another critical input for modelling lahar transport (de' Michieli Vitturi et al., 2024). For this study we used those reconstructed by Di Vito et al. (2024b) based on the field data from the catchments around Somma–Vesuvius (catchments 1 to 6 in Fig. 1b) integrated with those from Pozzelle quarry (Sulpizio et al., 2006) and from Vallo di Lauro (at the base of catchment 8 in Fig. 1b). In particular, both types of GSDs have been reconstructed by fitting a mixture of two Weibull distributions and then averaging the values of the local grain sizes sampled in the field at the base of these catchments separately. The reconstructed GSDs are reported in Fig. 3: for the Apennine catchments (7–11 in Fig. 1b) we used the GSD reconstructed from the field data taken in Vallo di Lauro (Fig. 3b), while for the Vesuvian catchments (1–6 in Fig. 1b) we obviously use the GSD reported in Fig. 3a and reconstructed from field data taken there.

Figure 3 Grain size distribution used for Vesuvian (a) and Apennine (b) catchments. The former has been derived from local samples from the Somma–Vesuvius catchments, while the latter are from samples from Vallo di Lauro (located at the base of catchment 8 in Fig. 1b).

Finally, in order to assess the effect of the uncertainty in the reconstructed GSDs we performed a sensitivity test (see Sect. S1 in the Supplement) of the weight of the coarse and fine populations on the simulated deposits, showing that it is not very critical in terms of simulated deposit, maximum thickness, and dynamic pressure 24 h after the flow mobilization.

2.4 Digital elevation model

For a correct modelling of the areas invaded by lahars it is necessary to use a digital elevation model (DEM) as accurate as possible. To this end we used a digital terrain model (DTM) derived from an airborne lidar survey of 2012 combined with the TIN Italy topography (Tarquini et al., 2007) in a portion of the sub-Apennine areas where the lidar data were not available. The lidar data were provided by the Italian Ministero dell'Ambiente e della Tutela del Territorio e del Mare (MATTM) through a series of ASCII files storing the elevation data in the latitude–longitude WGS84 reference system. The tiles, each covering 1 km², have been processed to create a single elevation matrix georeferenced in the WGS84-UTM-Zone 33N geodetic cartographic system and memorizing the elevation data at 32 bits. The obtained matrix, with spatial resolution of 1 m and vertical accuracy < 30 cm (Pizzimenti et al., 2016), was cleaned from residual anthropic or artificial features and subsequently resampled at 10 m cell size in order to be combined with the TIN Italy model. The resulting matrix (4129 × 5088 cells), covering 1600 km², was used as the topography model for the area (Fig. 1c) for the simulations, discretized at a computational grid of 50 m × 50 m, which was tested by de' Michieli Vitturi et al. (2024) as a good resolution able to reproduce the main features of the flow in reasonable computational times.

3 Methods

3.1 Sampling strategy

In order to explore the natural variability in the processes governing lahar initiation, we first identify the key independent parameters. Given a catchment, the key parameters are the initial flow volume and initial solid volumetric fraction (α_s); the first depends on the initial area of lahar detachment, on the thickness of the available remobilizable deposits, and on the available water from both rain or other external water sources and the water content filling soil pores.

The initial value of α_s is very variable and hard to assess; thus, we uniformly sample it in the range [0.10–0.60], which represents general limits for a lahar, encompassing a wide range of flows from debris flows (solid concentration α_s>0.5) to hyperconcentrated flows (α_s= 0.5–0.1) to muddy streamflows (α_s<0.1) (Vallance and Iverson, 2015; Neglia et al., 2022).

Volcanoclastic deposits from Vesuvius are characterized by an effective porosity (α_d) that has been estimated as 0.37 ± 0.10 (Di Vito et al., 2024b). Thus, α_d is sampled from a Gaussian distribution with a mean of 0.37 and standard deviation equal to 0.03 so that 99 % of the sampled values are within the above estimate.

To identify the possible initial area of lahar formation in a given catchment, we first select the grid points (Fig. 1a, b) falling within the catchment. Then, we rely on three empirical assumptions based on field evidence (Bisson et al., 2014): (i) only grid cells that are “steep enough”, i.e. having a slope larger than 20–30°, can be the site of remobilization. (ii) The steeper a cell, the more likely its deposits will be remobilized; however, (iii) on very steep cells (slope > 40°), deposits tend not to accumulate, and thus we assume no deposit is available for remobilization on such grid cells.

In order to define the potential initial area of lahar generation accounting for points (i) to (iii) above, we assume that remobilization can occur in a grid cell if it is steeper than a given δ_min and less steep than 40°. Then we distinguish between Somma–Vesuvius catchments, where most of the deposits are fine-grained PDC deposits and are thus easier to be remobilized, and Apennine catchments, where deposits are coarser-grained fallout layers more difficult to be remobilized. As pointed out by Pierson et al. (2013) the thickness of fine ash is also important because it can prolong low infiltration capacity and a high runoff rate. To reflect this feature, for each catchment and simulation, the slope threshold δ_min is sampled randomly from a triangular distribution, independently for Somma–Vesuvius and Apennine catchments. For Somma–Vesuvius we considered a lower-bound distribution δ_min of [20–30]° and for the Apennines [20–35]°. In this way, steeper cells are much more likely to be the site of remobilization (to reflect point ii) that can never occur at slopes lower than 20° (to reflect point i) and that definitely occurs above 30 and 35°, respectively (to reflect the difference in the types of deposits). These values are in agreement with Pierson et al. (2013).

Having defined the initial solid volumetric fraction and the area of lahar initiation, we can univocally define the initial flow volume by taking into consideration some physical constraints given by the availability of solid deposits and of water from pores and from rain in the domain of remobilization. The simulated volumes in each catchment (minimum, maximum, mean, and some percentiles) are provided in Table 1. As typical conditions during lahar flows we assume that the deposits are already water-saturated by an extra amount of water from previous rain when the lahar is triggered (e.g. Fiorillo and Wilson, 2004; Di Vito et al., 2024b). Following de' Michieli Vitturi et al. (2024), the thickness of available compacted deposit (i.e. devoid of the water filling its pore), ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{s}}$, can be bounded as

$$\begin{array}{}\text{(1)}& {\stackrel{\mathrm{‾}}{h}}_{\mathrm{s}}\le \min \left\{h_{\mathrm{d}}\left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{d}}\right),{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{s}}\left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{d}}\right)}{\left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{s}}-{\mathit{\alpha }}_{\mathrm{d}}\right)}}{h}_{\mathrm{r}}\right\},\end{array}$$

with being h_d the thickness of the water-saturated deposit, h_r the column of available rainwater, α_d the deposit effective porosity, and α_s the initial solid volumetric fraction. Following the considerations reported in de' Michieli Vitturi et al. (2024), we set h_r equal to 500 mm as a conservative value corresponding to the maximum 2 d accumulated rainfall over the Somma–Vesuvius area in the last 70 years in the region, although in principle it can be sampled from a suitable distribution (in this case an assumed Dirac delta distribution).

Table 1 Values (in m³) of the volume of simulated flows in each catchment. Numbering of catchments is the same as in Fig. 1b. In the different columns we provide the minimum (non-zero) simulated volume, the maximum, three percentiles (5th, 50th, and 95th), and the mean.

3.2 PVHA workflow and combination of the model simulation output

After the sampling of the relevant parameters, for our PVHA we follow the step of the workflow illustrated in Fig. 4. For each catchment, we run N_s= 100 simulations with the generalized IMEX_SfloW2D model (see purple box in Fig. 4), each with a different set of initial values for the model parameters. The simulations are run over a sub-domain with a resolution of 50 m × 50 m, cut in order to save computational time to avoid simulating negligible flow thicknesses over very distal grid points from a given catchment.

Figure 4 Workflow describing the sampling of ESP (and their combination through Eq. 1) and IMEX_SfloW2D simulation (purple box); the post-processing of the simulation outputs (brown, pink, and yellow boxes) to apply Eqs. (2) and (3) to achieve the hazard curves, maps and surfaces, and probability maps for a given catchment (green box); and the further processing through Eqs. (4) and (5) to achieve those for any catchment (cyan box). The number of catchments N_c considered in this application is N_c= 11, and the number of simulations for each catchment is N_sim= 100. The darker purple sub-box indicates that, in the present study, we set h_r fixed at 0.5 m, but it can ideally be sampled from climatological distributions.

For each source catchment i, we compute the probability (given N_s simulations) of exceeding a given threshold h_j in maximum flow thickness in a target grid point x as

$$\begin{array}{}\text{(2)}& p_{i,h_j}\left(x\right)={\displaystyle \frac{\sum _{k=\mathrm{1}}^{N_{\mathrm{s}}}{\mathit{\theta }}_{ik}}{N_{\mathrm{s}}}},\end{array}$$

where θ_ik is 1 if, in the kth simulation from catchment i, the maximum simulated flow thickness in x was larger than h_j or 0 otherwise. The set of thresholds in flow thickness used is composed of 21 values: from 0.1 to 1 m at a step of 0.1 m, from 1.5 to 4 m at a step of 0.5 m, from 6 to 10 m at a step of 2 m, and the last values at 15 and 20 m.

Similarly, for each source catchment i, we compute the probability of simultaneously exceeding a pair (h_j,P_l) of threshold values in flow thickness and dynamic pressure in a target grid point x as

$$\begin{array}{}\text{(3)}& p_{i,h_j,P_l}\left(x\right)={\displaystyle \frac{\sum _{k=\mathrm{1}}^{N_{\mathrm{s}}}{\mathit{\theta }}_{ik}}{N_{\mathrm{s}}}}.\end{array}$$

Here θ_ik is 1 if, in the kth simulation for catchment i, the pair (h_j,P_l) has been overcome in x at least once or 0 otherwise. In this case the set of threshold pairs in flow thickness and dynamic pressure simultaneously overcome is composed of 36 values, consisting of all the possible combinations of the values 0, 0.1, 0.5, 1, 2, and 5 m of thickness and 0, 0.5, 1, 2, 5, and 30 kPa for dynamic pressure. The value equal to 0 in one of the two target variables that allows computing the probability maps of exceeding the six thresholds in the other variable. In other words, by selecting a threshold value equal to 0 in one variable and a given threshold t_i in the other variable, we can visualize the probability of a maximum value larger than or equal to t_i in the latter variable.

We then combine all the N_c= 11 catchments together by computing

$$\begin{array}{}\text{(4)}& p_{h_j}\left(x\right)=\mathrm{1}-\prod _{i=\mathrm{1}}^{N_{\mathrm{c}}}\left(\mathrm{1}-p_{i,h_j}\left(x\right)\right)\end{array}$$

and

$$\begin{array}{}\text{(5)}& p_{h_j,P_l}\left(x\right)=\mathrm{1}-\prod _{i=\mathrm{1}}^{N_{\mathrm{c}}}\left(\mathrm{1}-p_{i,h_j,P_l}\left(x\right)\right),\end{array}$$

which respectively yield the probability of exceeding the maximum flow thickness h_j and the probability of exceeding simultaneously the pair (h_j,P_l) in flow thickness and dynamic pressure in the target grid point x from at least one catchment, given the remobilization of the volcanoclastic deposits from a medium-sized eruption at Somma–Vesuvius.

4 Results

The ranges in the initial volume of simulated lahars are shown in Table 1. We can see that some catchments are more prone to generating large initial volumes, specifically catchments 1 to 4 (on Somma–Vesuvius) and 7, 8, and 9 (on the Apennines) which have large maximum and/or mean simulated volumes. The latter catchments are characterized by a large extension and the former by the availability of a large thickness of deposits from PDCs. The catchments on Somma–Vesuvius and those in the northeastern Apennine section (numbers 7 and 8 in our case) were also identified in Tierz et al. (2017) as those able to generate larger-initial-volume lahars in the case of a medium-sized eruption. We also notice that in some Apennine catchments (numbers 7 to 11) some simulations do not have significant deposits from tephra fallout to be remobilized (null initial volumes, probably because we simulated deposits from eruptions under wind fields not directed towards those catchments).

Following the approach described in the previous section, for each target grid point we can build a hazard curve (e.g. Tonini et al., 2015) for the maximum flow thickness from each different catchment (Eqs. 2 and 3) and one from at least one catchment (Eqs. 4 and 5). As an example, we show these results on the grid point located in San Marzano sul Sarno (S in Fig. 1a), an inhabited town in the Campanian Plain, in Fig. 5: in Fig. 5b–e we show the hazard curves only relative to the catchments that are able to generate a lahar reaching this target point and do not show the hazard curves for the other catchments, whereas in Fig. 5f we show the hazard curve from any catchment. These hazard curves can be cut at different thresholds in terms of exceedance probability or thickness thresholds and mapped on the whole target domain to obtain hazard or probability maps, respectively, such as those in Fig. 6a–k (at 5 % exceedance probability) and Fig. 7a–k (for a flow maximum thickness of 1 m).

Figure 5 Hazard curves for the maximum flow thickness from catchments 3, 4, 8, and 9 (shown in panel a) from San Marzano sul Sarno (point S) (the hazard curves are respectively shown in panels b, c, d, e) and from any catchment (f).

Figure 6 Hazard map at the 5 % probability level, showing the maximum flow thickness with the probability to be exceeded given a lahar from specific catchments (a–k) and from any catchment (l).

Figure 7 Probability map for the maximum flow thickness of 1 m, showing the exceedance probability given a lahar from specific catchments (a–k) and from any catchment (l).

In Figs. 6l and 7l we show the same hazard and probability maps for at least one lahar from any catchment, as derived from Eq. (5).

Similarly to hazard curves, we can think of a “hazard surface” in terms of pairs of threshold values in flow thickness and dynamic pressure simultaneously overcome at every grid point. These hazard surfaces aim at an easier visualization for impact assessment: in reality, in order to evaluate the impact of a lahar on the built environment, it is the joint effect of these two parameters that gives a more relevant measure of the impact.

An example for a Vesuvian location (the Casilli train station, point C) is shown in Fig. 8 from specific catchments (Fig. 8b and c) and from any catchment (Fig. 8d). This location has been selected as it is impacted by more than one catchment, which is not common on the studied domain. By cutting these surfaces at specific pairs of values for the flow thickness and flow dynamic pressure, we achieve maps showing the probability of simultaneously exceeding the pair of threshold values selected.

Figure 8 Hazard surface for the simultaneous exceeding of different pairs of values in maximum flow thickness and flow dynamic pressure at the Casilli (C) train station. Panels (b) and (c) respectively show the probability to simultaneously exceed the threshold pairs at least once during the flow, given a lahar from catchments 3 and 4 (location of catchments shown in panel a). Panel (d) shows the hazard surface from all the catchments.

In Fig. 9 we show the example for 0.5 m of flow thickness and 1 kPa in flow dynamic pressure from each different catchment (Fig. 9a–k) and from any catchment, as derived from a catchment (Fig. 9l).

Figure 9 Probability map for the simultaneous exceeding of 0.5 m in maximum flow thickness and 1 kPa in flow dynamic pressure, showing the probability to simultaneously exceed the two thresholds at least once during the flow, given a lahar from specific catchments (a–k) and from any catchment (l).

Typically, the output of probabilistic hazard studies is a portfolio of scientific products (different hazard and/or probability maps at different thresholds in terms of exceedance probability and/or intensity, respectively, corresponding to different average return times). In this respect, the present study is no different, and in Sect. S2 we provide the whole set (apart from those already given in Figs. 6, 7, and 9) of the following:

• probability maps for the maximum flow thickness (and at the other 20 thresholds in thickness considered except the one in Fig. 7 – Figs. S7 to S26);

• hazard maps for the maximum flow thickness at the exceedance probability thresholds of 1 %, 2 %, 10 %, 50 %, and 90 % (Figs. S27 to S31); and

• probability maps for the simultaneous exceeding of flow thickness and dynamic pressure threshold pairs (at the other 35 threshold pairs considered except the one in Fig. 9, Figs. S32 to S65).

The visual inspection of the portfolio of probabilistic maps for lahar maximum thickness shows that, in the case of a reference size eruption at Vesuvius and heavy rain, flows of maximum thickness of half a metre (e.g. Fig. S11) are to be expected in the southwestern, northwestern, northern, northeastern, and eastern sectors around Vesuvius, having conditional probabilities very close to 1 even in valleys down to 10–15 km from the volcano summit. In the southwestern sector they would likely reach the shoreline, in agreement with Tierz et al. (2017), with decimetric maximum thickness. In such sectors, the hazard maps (e.g. Fig. S30) tell us that a flow of metric maximum thickness has a 50 % chance to be exceeded in the case of a primary lahar from the reference eruption. In the areas threatened by Apennine catchments, the highest conditional probability for such a scale in maximum flow thickness is found in valleys, especially in the Vallo di Lauro, Sarno, and Castellammare di Stabia areas, with values up to 50 %. This is in agreement with the findings by Tierz et al. (2017), who found metric flow depths when considering lahars triggered from Vallo di Lauro and Avella catchments.

Flows of maximum thickness of the order of 1–2 m (e.g. Figs. 7 and S17) have a conditional probability close to 1 only in the bottom of valleys from the northwest to southeast (clockwise) sectors around Somma–Vesuvius, up to approximately 5 km. In the other Vesuvian valleys around the volcano they have 10 %–50 % conditional probabilities, whereas in the Apennine valleys the maximum conditional probability is about 30 % for this order of flow thickness, again in the Vallo di Lauro, Sarno, and Castellammare di Stabia areas. Flows of maximum thickness of the order of 5 m or more are conditionally unlikely (probability less than 10 %) and confined to valley bottoms, at least at the resolution used in this study.

Figure 10 Probability maps for the simultaneous exceeding of (a) 0.5 m in maximum flow thickness and 2 kPa in flow dynamic pressure as well as (b) 1 m and 5 kPa. The maps show only points where the probability to simultaneously exceed the two thresholds at least once during the flow is larger than 5 %, given a lahar from any catchment.

When accounting for the simultaneous exceeding of conjoint thresholds in flow thickness and dynamic pressure, we can see that for the Vesuvian area the maximum conditional probability for flows of 0.5–1 m and 0.5–5 kPa is found downhill in valley bottoms (e.g. Figs. 9, 10, S45–S47, S50–S53), with the highest values being located well below 500 m altitude: these areas are densely populated and built. In these areas, these ranges of conjoint thresholds have a conditional probability close to 1 in the bottom of the valleys. High dynamic pressure values (e.g. 5 kPa or larger – Figs. S41, S42, S46, S47, S52, S53) are overcome only in steep slopes, where presumably the flow speeds up. Their conditional probability is quite low in Apennine steep flanks, where such dynamic pressure thresholds are overcome with very thin flows (of the order of 0.1 m, Fig. S42). In steep flanks of Vesuvian catchments 1 to 4, high dynamic pressure flows (5 kPa or larger) are more probable than in the Apennines, still associated with flows of 0.1–0.5 m (Figs. S42 and S47). In Di Vito et al. (2024b), a reverse engineering approach is used to invert the occurrence of external clasts (bricks, walls, limestone fragments, etc.) into the volcaniclastic deposits to estimate the local flow dynamic pressure, velocity, and thickness. For flows with an estimated thickness of 0.5–1 m, the most characteristic range of dynamic pressure is 4–8 kPa, corresponding to a representative range of velocity of 2–4 m s⁻¹. This means that the higher the velocity (and dynamic pressure) the higher the capability of the flow to entrain accidental clasts (or damage infrastructures). Comparing such estimates with the probabilistic maps from the present study, we see that the probability to overcome these combinations of thresholds in thickness and dynamic pressure (0.5–1 m; 2–30 kPa, Figs. S46–S48, S52–S54) is statistically significant, being larger than 5 % in several grid points. In particular, the pairs of lower thresholds (e.g. 0.5 m and 2 kPa, Fig. S46) have a probability larger than 5 % in steep valleys throughout the domain and in some locations at the mouth of the narrowest and steepest valleys (e.g. the Vesuvian ones or in Vallo di Lauro). The pairs of higher thresholds (e.g. 1 m and 5 kPa, Fig. S53) have a probability larger than 5 % only in the steepest valleys on Vesuvius slopes. As a concluding remark, we can state that the estimated combined values of flow thickness and dynamic pressure from field data in Di Vito et al. (2024b) are well-captured by our probabilistic maps and do not represent outliers or unlikely values with respect to them.

As mentioned above, the probability maps shown in Figs. S32 to S37 (where the threshold in maximum flow thickness is 0 m) in practice show the probability of flows with maximum dynamic pressure exceeding the six thresholds considered (see Sect. 3.2). It is important to highlight that where high values of dynamic pressure are overcome, it may be due to simulated flows of negligible thickness: this is the reason why we decided to consider the simultaneous exceeding of non-zero thresholds in thickness and dynamic pressure (e.g. in Figs. 9 and S38–S65). We think that this type of information could be of great importance when incorporating probabilistic hazard into impact and quantitative risk assessment (e.g. Zuccaro et al., 2008; Zuccaro and De Gregorio, 2013).

Comparing the different catchments, the catchments threatening the largest areas (e.g. the area invaded by at least 10 cm thick flows with exceedance probability of 10 %, Fig. S29) are the Apennine sectors 7 (Vallo di Avella), 8 (Vallo di Lauro), and 9 (upward of the towns of Nocera Inferiore, Gragnano, and Castellammare di Stabia, a very densely inhabited area). Overall, the Vesuvian sectors 2 (north), 3 (northeast), and 4 (east) show the largest maximum expected thickness in hazard maps, given an exceedance probability (e.g. Figs. 6, S27–S31), and the largest probability of given maximum thickness in probability maps (e.g. Figs. 7, S7–S26). In a few words, the smaller conditional hazard ubiquitously found in the Apennine areas, compared to the Vesuvian ones, is due to the smaller quantity of available sediments to be remobilized (only from medio-distal fallout, whereas on Somma–Vesuvius both proximal fallout and PDC deposits are available) and their coarser grain size. In our simulation strategy scheme, the latter feature has been assumed to have higher resistance to remobilization.

In a few points (such as S and C, respectively shown in Figs. 5 and 8) the hazard is due to more than one catchment, but this appears to be an exception: overall, most of the inspected domain is threatened by one catchment.

5 Discussion and conclusions

The present study aims at providing a probabilistic assessment of lahar hazard conditional on the occurrence of a medium-sized reference eruption from Somma–Vesuvius. For the first time, such a hazard assessment has considered several sources of uncertainties previously overlooked, such as uncertainties in the initial volumes and initial detached areas of lahars, as well as effects of erosion and deposition processes. This has been possible thanks to both the model formulation (de' Michieli Vitturi et al., 2024) and the availability of relatively fast computational resources at INGV. Both factors enabled the simulation, in a reasonable time, of 100 different scenarios from each of the 11 catchments examined here, totalling a larger-than-ever number of simulations for a lahar hazard assessment.

Where possible, constraints on the range of parameters values, or constraints on the probabilistic distribution describing relevant parameters, were obtained by comparison with field data (for example, the GSD or the sediment porosity α_d). For other parameters, only loose constraints were possible, and a maximum-ignorance uniform probability distribution has been used to describe such cases.

Under the hazard–risk separation principle (Jordan et al., 2014; Marzocchi et al., 2021), we remark that it is the role of volcanologists to quantify uncertain scientific information in a way that can be used to mitigate risk.

The present study does not tackle the variability related to the different possible eruption sizes at Somma–Vesuvius, as we have focused on the medium-sized eruption only (Sandri et al., 2016), and we provide the quantification of lahar hazard in the case of such an event. Future research will be devoted to extending the present analysis to other eruptive sizes, especially to larger ones (i.e. Plinian events similar to Pompeii 79 CE or Avellino and Mercato eruptions; see Gurioli et al., 2010), whose lahar hazard potential is significant (Tierz et al., 2017).

However, there are a number of papers in the literature quantifying the probability of a Somma–Vesuvius eruption of medium size similar to the one accounted for in this study (e.g. Marzocchi et al., 2004; Selva et al., 2022). According to Selva et al. (2022), an estimate of the probability of at least one eruption from Somma–Vesuvius, given the current low-activity period, is about 34 % in 50 years (Fig. 6 in Selva et al., 2022), and the conditional probability of a medium-sized event (i.e. a VEI = 4) is about 20 % (Fig. S7 in Selva et al., 2022). Consequently, a gross estimate of the probability of an eruption of medium size at Vesuvius is about 7 % in 50 years. To our knowledge, the probability of syn-eruptive lahar generation, given an eruption of medium size, has never been quantified. It may be broadly estimated on the basis of the frequency of lahar triggering observed at analogue volcanoes (Tierz et al., 2019) with similar climates. Such development is beyond the goals of this study and is foreseen as a future research project.

Finally, we remark that, in the present study, two limitations are related to (i) the number of simulations performed from each catchment (100) and (ii) the assumption of available rainfall, which we have fixed at 50 cm. As for the former limitation, it is mainly dictated by the availability of computational resources. Such a number has been chosen to carry out the simulations in a reasonable time, and it may be overcome in the future as more performing codes and resources are available. At the same time, the exploration of lahar hazard in the case of other possible eruptive size classes from Somma–Vesuvius (Sandri et al., 2016) could be examined in future developments of this work. As for the latter limitation, we have also inferred that possible rain due to condensation of magma-exsolved water vapour in the umbrella region that we have roughly estimated in an extra 5 to 10 cm of water does not significantly increase such an upper limit. We analysed the data from the rain gauge operating at the historical building of the Osservatorio Vesuviano on Mt. Vesuvius since 1940 (Ricciardi et al., 2007), finding that this amount represents approximately the maximum recorded accumulated rain over 2 d since about 1950 in the Campanian region. This value is similar to the maximum rainfall among the episodes of lahars reported by Fiorillo and Wilson (2004). In this perspective, this limit has been taken as conservative. In fact, it implies that the simulated lahars have larger initial volumes, since more water is available, compared to cases with a smaller amount of rainfall. However, due to climate change, the 2 d intensity of rain may be larger over this area in the coming decades. Furthermore, the comparison of the 2 d accumulated rainfall at Vesuvius with the occurrence of lahar cases in Campania shows very little correlation in time (Cantelli, 2021). Following these considerations, we acknowledge that relaxing this assumption, on the one hand, would allow simulating more frequent and smaller-sized lahars that appear to have occurred in the last 50 years even with smaller rainfall intensity (e.g. the Sarno event in 1998) and, on the other hand, would allow accounting for potentially larger rainfall intensities expected in the future by ongoing climate change, which is the subject of future work.