Lahar risk assessment from source identification to potential impact analysis: the case of Vulcano Island, Italy

Abstract

Lahars are rapid flows composed of water and volcaniclastic sediments, which have the potential to impact residential buildings and critical infrastructure as well as to disrupt critical services, especially in the absence of hazard-based land-use planning. Their destructive power is mostly associated with their velocity (related to internal flow properties and topographic interactions) and to their ability to bury buildings and structures (due to deposit thickness). The distance reached by lahars depends on their volume, on sediments/water ratio, as well as on the geometrical properties of the topography where they propagate. Here we present the assessment of risk associated with lahar using Vulcano island (Italy) as a case study. First, we estimated an initial lahar source volume considering the remobilisation by intense rain events of the tephra fallout on the slopes of the La Fossa cone (the active system on the island), where the tephra fallout is associated with the most likely scenario (e.g. long-lasting Vulcanian cycle). Second, we modelled and identified the potential syn-eruptive lahar impact areas on the northern sector of Vulcano, where residential and touristic facilities are located. We tested a range of parameters (e.g., entrainment capability, consolidation of tephra fallout deposit, friction angle) that can influence lahar propagation output both in terms of intensity of the event and extent of the inundation area. Finally, exposure and vulnerability surveys were carried out in order to compile exposure and risk maps for lahar-flow front velocity (semi-quantitative indicator-based risk assessment) and final lahar-deposit thickness (qualitative exposure-based risk assessment). Main outcomes show that the syn-eruptive lahar scenario with medium entrainment capability produces the highest impact associated with building burial by the final lahar deposit. Nonetheless, the syn-eruptive lahar scenario with low entrainment capacity is associated with higher runout and results in the highest impact associated with lahar-flow velocities. Based on our simulations, two critical infrastructures (telecommunication and power plant), as well as the main road crossing the island are exposed to potential lahar impacts (either due to lahar-flow velocity or lahar-deposit thickness or both). These results show that a risk-based spatial planning of the island could represent a valuable strategy to reduce the volcanic risk in the long term.

Introduction

Volcanic eruptions are associated with a variety of primary hazards, such as tephra fallout, Pyroclastic Density Currents (PDCs), toxic gas emissions and lava flows, as well as secondary hazards such as landslides, tsunamis and lahars. All these hazards can strike settled areas nearby active and dormant volcanoes at different spatial and temporal scales. In this framework, lahars represent one of the most impactful volcanic hazards that can potentially strike areas up to hundreds of kilometres from active volcanoes and cause a high number of fatalities as demonstrated by the event that occurred in 1985 at Nevado del Ruiz volcano (Colombia) (Lowe et al. 1986). Lahar is an Indonesian term describing water-saturated flows; both liquid and solid interactions influence their behaviour and distinguish them from debris avalanches (Vallance and Iverson 2015), or from jökulhlaups released from glaciers in Iceland (Gudmundsson 2015). The term lahar was used for the first time by Esher in 1922 (Neall 1976) to indicate a mudflow containing debris and blocks of chiefly volcanic origin (Van Bemmelen 1949) or a flowing mixture of rock debris and water of volcaniclastic materials on the flank of a volcano (Gary et al. 1974). Lahars can be produced either during or short after an eruption (syn-eruptive) and for years to decades after an eruption (post-eruptive) and have the potential to cause significant damage to buildings, public facilities, critical infrastructures as well as losses of human life and disruption to critical services (e.g. Sulpizio et al. 2006; Wilson et al. 2014; Jenkins et al. 2015; Mead et al. 2017).

The destructiveness of lahars is a function of their speed, thickness, runout, sediment fraction and grain-size (e.g., Jenkins et al. 2015), which largely control their kinetic energy and the thickness of deposits. Large lahars commonly reach velocities >20 m/s on the lower flanks of volcanoes and can maintain velocities higher than 10 m/s for more than 50 km from their source when confined to narrow canyons (Pierson et al. 2014; Thouret et al. 2020). In medial to distal areas, lahars are typically confined in river valleys and flood plains, while become rapidly less energetic with increasing distance from vent due to a decrease of flow confinement and slope steepness (Spence et al. 2004; Jenkins et al. 2015; Vallance and Iverson 2015).

Their associated impacts rapidly evolve from physical damage (in flow paths, river valleys and in proximal areas) to disruption due to burial further down floodplains where lahars can deposit material with low energies (Baxter et al. 2005). Nonetheless, empirical, experimental and theoretical data describing vulnerability to lahars is currently rather limited (e.g., Blong 1984; Wilson et al. 2014; Jenkins et al. 2015; Dagá et al. 2018; Thouret et al. 2020). Buildings can be damaged through a number of mechanisms induced mainly by: (1) hydrodynamic pressure, (2) hydrostatic pressure, (3) collisional forces of boulders acting as missiles, (4) burial of buildings and other types of infrastructures under high and widespread thickness of deposits, and (5) indirect damage caused by chemical and biological processes such as seeping induced weakness of mortar (Wilson et al. 2014; Mead et al. 2017; Thouret et al. 2020). Two of the main physical parameters typically considered to assess the potential impact of lahars are, therefore, the flow depth (i.e., hydrostatic pressure) and the dynamic pressure (e.g., Wilson et al. 2014). Where enough information is available on building strength (e.g., tensile strength of the masonry wall, wall thickness and wall width), depth-pressure curves can be derived to assess building loss (Mead et al. 2017; Thouret et al. 2020).

A detailed description of various modelling approaches used in the literature to simulate the runout path and the inundation is available in Thouret et al. (2020). Some of these models identify inundation areas with no information on flow velocity and thickness (e.g. LAHARZ; Schilling 1998). More recently, computational fluid dynamics models have been integrated with flow velocity and thickness, and, therefore, can provide more insights into the potential impact to buildings and infrastructures, such as FLO-2D (O’Brien et al. 1993), TITAN2D (Pitman and Le 2005), LaharFlow (Tierz et al. 2017). Amongst this latter generation of models, the Smoothed Particle Hydrodynamic (SPH) modelling approach was successfully implemented by Pastor et al. (2009) in the GeoFlow_SPH software and applied to case studies of fluidized mass movements (e.g. Cascini et al. 2014; Cuomo et al. 2014; Pastor et al. 2014). Thanks to multiple rheological input laws (viscous, Bingham, viscoplastic, constant stress), SPH has been successfully applied to reproduce various lahar behaviours (e.g. Dumaisnil et al. 2010; Lafarge et al. 2012; Cascini et al. 2014; Procter et al. 2020) and characteristics (e.g. trajectory, inundation areas, flow velocity, and deposit thickness for Huiloac Gorge lahars at Popocatépetl volcano; Haddad et al. 2010, 2016). SPH was also used to develop fragility functions in the form of critical depth–pressure curve and quantify the potential impact on residential buildings in the city of Arequipa (Peru) as a result of future lahars associated with El Misti volcano (Mead et al. 2017).

Regardless of the variety of models developed to describe lahar inundation and the various applications to assess the associated hazard, only few studies have investigated the associated vulnerability (e.g. Jenkins et al. 2015; Mead et al. 2017; Thouret et al. 2020), which is key to quantitative impact analyses and sustainable development plans around volcanic areas. In this paper, we present a new methodology for the assessment of lahar impact using Vulcano island (Italy) as a case study. Because the island is periodically affected by intense rain events, the remobilization of tephra deposits remains an active process (e.g. Di Traglia et al. 2013; Baumann et al. 2019). Galderisi et al. (2013) already showed how the north sector of Vulcano island is the most exposed to lahar generation and, based on the characteristics of buildings and urban fabrics, is characterised by various degrees of vulnerability. However, the potential impact on the buildings could not be evaluated by Galderisi et al. (2013) as the hazard assessment was carried out using the model LAHARZ. Here we describe all the steps to assess risk associated with lahar inundation (i.e. potential impact) based on the assessment of hazard (flow velocity and deposit thickness computed with the SPH model), exposure (i.e. distribution of elements at risk) and physical vulnerability (i.e. susceptibility to physical damage).

Geographical and geological setting of the study area

Vulcano Island, located in the southern part of the Tyrrhenian Sea (Southern Italy), is the southernmost of the seven emerged Aeolian Islands. It is part of the Lipari Municipality and has a surface of about 20 km² with a total of 1,282 permanent residents; almost 70% of the population is located in the Porto-Vulcanello area and 30% in Piano. In the month of August (peak of the high season), up to about 28,000 tourists per month visit the island (Bonadonna et al. 2021). The main economic activities and tourist attractions are located in the northern sector of the Island. Most critical infrastructures and facilities are also located in the north including the Ponente and Levante harbours, the fuel power plant, key stations of the water system (including the recently built waste-water plant), the telecommunication centre, one church, the Red Cross, the only pharmacy of the island and the Medical Centre (Fig. 1). The main roads play a strategic role, since they represent the access to the harbours (Porto Levante and Porto Ponente) and connect the Porto area with Vulcanello and Piano areas (Fig. 1).

Fig. 1

Vulcano Island map showing critical facilities and infrastructures. The largest settlements are indicated. The inset shows the geographic location

Vulcano is a composite volcanic edifice located at the southern extremity of the central segment of the Aeolian Islands, along the Tindari Letojanni (TL) fault. The island is the result of four main eruptive cycles (Primordial Vulcano, Lentia - Mastro Minico, La Fossa and Vulcanello). Two multi-collapsed structures were formed during the geological history of the Island: the “Caldera del Piano” and “Caldera de La Fossa”, where the actual active La Fossa cone is located. La Fossa is a 391 m-high cone with a basal diameter of 1 km, that was built since 5.5 ka (up to historical times) through recurrent hydromagmatic to Vulcanian explosive activities, which produced PDCs and tephra fallout deposits alternating with viscous lava flows (De Astis et al. 2013). The last eruption of La Fossa cone occurred in 1888–1890. The last 1,000 years of eruptive activity were investigated through a stratigraphic reconstruction by Frazzetta et al. (1983), De Astis et al. (2013) and Di Traglia et al. (2013). Di Trapani et al. (2011) also reported that the recent deposits were intensely eroded by running water through the generation and deposition of small lahars. During the last hundreds of years, lahars occurred at La Fossa cone by the remobilisation of pyroclastic material deposited by several Vulcanian and phreatic explosive events occurred after the XIII century. Stratigraphic sections from exposure along the main creeks around La Fossa cone show different sequences of deposits associated to lahar events with bed thickness between 0.1 and 0.6 m in the middle and lower catchments of the cone (Baumann et al. 2019). These loose and incoherent deposits cover the Tufi Varicolori layer, an impermeable substratum that favours the remobilisation of the overlying deposits through lahar generation (Ferrucci et al. 2005). According to meteorological data of the National Institute of Geophysics and Volcanology (INGV-Palermo) obtained from a station located in the Lentia caldera rim (data collected between 2010 and 2014), heavy precipitations (torrential events) occurred three times in 2010 and 2011, two times in 2012, and one time in 2013 and 2014. These rainfall events can last 2–3 hours to 3 days and typically occur in September, October, November, December, and, more rarely, in May. The presence of layered, fine-grained tephra (lapilli and ash), and the steep slopes without a significant vegetation cover favour the remobilization of volcanic deposits as lahars (Ferrucci et al. 2005; Di Traglia 2011; Di Traglia et al. 2013). The Porto Levante area has several flood-control engineering structures designed for rainwater collection (Fig. 2). At the base of the La Fossa Cone, along the main watershed on the west side of the Pietre Cotte lava flow, three drainage grids have been installed to filter water from coarse sediments that come down from the top of the volcanic cone. In October 2017, these structures were in bad condition with some plugging by sediment and organic material. In the Porto Levante area, two drainage tunnels are also present to carry rainfall runoff from the grids at the base of volcanic cone to the sea (Fig. 2).

Fig. 2

Mitigation measures for floods, debris flows and lahars at Vulcano porto area. A-B) Drainage grids at the base of the Pietre Cotte Lava Flow; C) Drainage Tunnels at Levante Harbour

The main objective of this paper is the assessment of risk of syn-eruptive lahar inundation associated with the occurrence of a future, long-lasting Vulcanian eruption of La Fossa cone of similar magnitude of that of 1888-90. Among all explosive eruptions, this scenario represents the most frequent and voluminous events of La Fossa volcano. In fact, at least five of such long-lasting episodes occurred in the last 1000 years, corresponding to annual frequencies of 5x10^-3 events per year (Selva et al. 2020).

Methods

The methodology presented here follows a logical sequence from the generation of the hazardous phenomena to its impact on exposed buildings. First, we estimate the conditions leading to lahar generation. Second, we model lahars associated with various scenarios using GeoFlow_SPH (Smoothed Particle Hydrodynamic) model (Pastor et al. 2009) and quantify their hazards using various Hazard Impact Metrics (HIM). Third, the distribution of elements at risk is characterized based on the resulting lahar inundation areas. Fourth, we estimate the vulnerability of the exposed built environment by calibrating the vulnerability matrix developed by Papathoma-Köhle et al. (2012) to Vulcano using in-situ and remote surveys of the buildings. Finally, two risk maps are compiled. The first risk map considers the flow front velocity and the vulnerability score of each exposed building, whereas the second risk map uses the final deposit thickness and the height of buildings to assess the fraction of building burial. Hazard and exposure data were included in a Geographic Information System (GIS, ESRI ArcMap) for geospatial analyses.

Step 1: Hazard assessment

Lahar source volume estimation

The estimation of lahar source volume was based on both probabilistic modelling of tephra sedimentation and modelling of the volume of deposit that could be potentially remobilised by rainfall. Biass et al. (2016) presented a probabilistic modelling of the potential tephra accumulation at Vulcano associated with a variety of eruptive scenarios using TEPHRA2, including a long-lasting Vulcanian cycle lasting up to 3 years. From this, probabilistic isomass maps were compiled, which quantify the tephra accumulation that occurred at a given percentile over all simulations of the probabilistic scenario. For instance, the distribution of tephra accumulations on a probabilistic isomass map for a probability threshold of 25% implies that only 25% of the total simulation produced larger accumulations. Baumann et al. (2019) have used these associated probabilistic isomass maps to estimate what scenario could produce the largest volume of unstable tephra deposit over the 22 catchments in the north sector of La Fossa cone (Fig. 3). Using the TRIGRS model (Baum et al. 2002) and exploring ranges of initial conditions (e.g. the saturated conductivity (Ks), cohesion (c), friction angle, total unit weight of the soil (ϒ), the saturated diffusivity (D0), saturated and residual water content, and the Gardner parameter), results for the Vulcanian scenario, suggest the largest unstable volume is reached after 18 months of eruption (considering no deposit remobilisation). Here, the following methodology is conditional to the choice of using a single scenario of tephra accumulation representing the worst-case scenario for unstable tephra deposit following a 1888-90–type Vulcanian eruption. Using a probabilistic isomass map compiled for a 25% probability of occurrence, Baumann et al. (2019) used the model TRIGRS to estimate the potential triggering source at 22 upper catchments located in the northern sector of La Fossa cone (Fig. 3).

Fig. 3

Tephra-fallout deposit on the northern sector of the La Fossa cone that could be remobilised by a syn-eruptive lahar as simulated with TRIGRS by Baumann et al. (2019) (i.e., equivalent to 18 months of Vulcanian eruption; see main text for details). Deposit thickness is indicated with different colours

SPH modelling of lahar propagation

Smoothed Particle Hydrodynamic (SPH) is a Lagrangian meshless method firstly introduced by Lucy (1977) and Gingold and Monaghan (1977) for applications in astrophysics and which discretises a propagating mass into a set of “moving particles”. The application of SPH was later expanded by Pastor et al. (2009) and Cuomo et al. (2014, 2016, 2019) for modelling the propagation of flows such as water saturated sediments and other mixture of multiphase materials (see also Braun et al. 2018; Cascini et al. 2020). The governing equations are reported in Pastor et al. (2009). In the SPH discretisation the unknowns and their derivatives are linked to the “moving particles”. Particularly, the unknowns are the velocity of the flow (v) and the pore water pressure at the base of the flow (${p}_w^b$), which can be both decomposed as the sum of two components related to: i) propagation and ii) consolidation along the normal direction to the ground surface. The vertical distribution of pore water pressure from ground surface is approximated using a quarter-cosinus shape function, with a zero value at the surface and zero gradient at the basal surface as given by Eq. 1:

$$\frac{{dp^b}_w}{dt}=\frac{\pi^2}{4{h}^2}{c}_v{p}_w^b$$

(1)

where ${p}_w^b$ is the excess of pore water pressure to the hydrostatic value, h is the propagation thickness and c_v is the consolidation coefficient.

The consolidation process reduces the pore water pressure while increasing the normal stress and the shear resistance at the flow base. The initial pore water pressure is taken into account through two input parameters: ${h}_w^{rel}$, which is the relative water height resulting from the ratio of water table height to deposit thickness, and ${p}_w^{rel}$, which is the relative pore-water pressure resulting from the ratio of pore-water pressure to liquefaction pressure inside the source area. Estimates of both parameters can be obtained from the analysis of the triggering stage as explained below. The importance of pore pressure dissipation during debris-flow propagation has been demonstrated in the literature (Pastor et al. 2009; Xu et al. 2012, Iverson and George 2014).

The frictional rheological model was assumed for the flowing mass, based on the high percentage of solid particles in the flow. The greater the flow thickness, the larger the basal shear resistance with a dependence on both friction angle and pore water pressure at the base of the flow, as in Eq. 2:

$${\tau}_b=-\left(\rho^{\prime}\cdot g\cdot h-{p}_w^b\right){\tan}{\varnothing}_b\cdot {\operatorname{sgn}}\left(\overline{v}\right)$$

(2)

where τ_b is the basal shear stress of the flow, ρ ′ s is the submerged density defined as ρ ′  = (1 − n)(ρ_s − ρ_w) where n is the deposit porosity, ρ_s is the particle density, ρ_w is the water density, g is the gravity acceleration, φ_b is the basal friction angle, sgn is the sign function and $\overline{v}$ is the depth-averaged flow velocity. The assumption to use the depth-averaged variables comes from the shallow thickness of the propagating debris flows compared to their length. Propagation and deposition of material are computed based on two fundamental equations: the balance of mass of the mixture and the balance of the linear momentum (i.e. mass multiplied by velocity) of the mixture. In doing that, the main driving force is the gravity while the resistance force is the shear stress at the base of the propagating mass. Deposition occurs in a given computation point when the velocity reduces to zero. SPH is a meshless method; as a result, integration is done over the moving particles (called nodes) that carry all the information regarding depth, velocity, displacements from the initial position and pore water pressure. The spatial resolution of the model depends on the resolution of the topographic model.

The equation of mass conservation of the flow accounts for the capability of the flow to entrain further material from the ground surface. The entrainment rate (e_r) depends on the thickness (h) and depth-averaged velocity of the flow (v) as well as the local slope angle (θ) of the ground surface where the flow moves on (Blanc et al. 2011) as reported in Eq. 3:

$${e}_r=h\cdot v\cdot K\cdot {\left( tan\ \theta\ \right)}^{2.5}$$

(3)

The erosion coefficient (K) depends on the characteristics of both the flow and the ground surface, and it is usually calibrated with site-specific information (Cuomo 2020). Based on that, the erosion depth h_er (m) is computed as the cumulative variation of the elevation of any point of the topography, as reported in Eq. 4. The “erosion time” (t_er) in a given point, defined by Cuomo et al. (2014) as the time lapse during which the flow depth and velocity are different than zero so that erosion rate is different than zero, can be retrieved as a model output. In each point of the Digital Surface Model the passage of the flowing material causes a different amount of eroded thickness during a different erosion time:

$${h}_{er}={\int}_0^{t_{er}}{e}_r dt={\int}_0^{t_{er}}h\cdot v\cdot K\cdot {\left( tan\theta \right)}^{2.5} dt$$

(4)

This kind of mathematical approach is usually classified as hydro-mechanical coupled because the consolidation process affects the basal pore water pressure, which influences the basal shear stress and, in turn, velocity and thickness of the flow, which are primary factors for entrainment.

Digital surface model and digital terrain model of the area

Both a Digital Surface Model (DSM), which captures both the natural and built/artificial features of the environment, and a Digital Terrain Model (DTM), which represents the Earth surface including vector features of the natural terrain, such as ridges, have been used in this study. Both the DSM and DTM are derived from a 2-m resolution Lidar point cloud acquired in 2017 by Ministero dell’Ambiente for the entire island, with an error on the vertical accuracy lower than ± 15 cm (MATTM; www.minambiente.it). For the numerical analysis and in order to minimize the computational cost of SPH, the DSM and DTM have been cut in the northern sector with a total of 295,693 cells, while 11,832 initial tephra thickness points, 3 m-spacing, corresponding to an initial lahar source volume of 24,857 m³ over an area of 119,048 m² were considered (Baumann et al. 2019). The DSM was used to explore the lahar propagation with SPH, while the combination of DSM and DTM was used the assess the height of buildings (see sections below). The modelling focused on the northern sector of La Fossa, where the main slope consists of a bottom substratum of red, impermeable ash deposits overlain by layers of loose, grey to black ash deposits with a thickness up to 0.26 m. The main morphological features of this slope are represented by the Pietre Cotte lava flow (an XVIII century rhyolitic flow which reached the foot of the slope) and two coalescent craters (Forgia Vecchia) related to phreatic explosions that occurred in the XV-XVI centuries.

Lahar hazard scenarios

The frictional rheological behaviour of the potential lahars was explored through a parametric analysis, exploring ranges of tephra porosity (0.5 to 0.7) and the internal friction angle (30° to 40°) obtained from direct shear laboratory tests of tephra samples associated with the 1888-90 eruption (Baumann et al. 2019). Assuming the density of solid particles as 2.7 g cm^-3 and fluid density as 1 g cm^-3, we obtained an apparent friction angle varying from 11° to 22°, i.e. tanφ from 0.2 to 0.4 (Eq. 2). These values are consistent with the range 0.34-0.40 obtained from the back-analysis of past debris flows that remobilised pyroclastic deposits a few kilometres from Vesuvius volcano (Cuomo et al. 2016). Furthermore, the relative pore water pressure (${p}_w^{rel}$; i.e. ratio of pore water pressure to liquefaction pressure) was set to 0.4 and relative water height (${h}_w^{rel}$; i.e ratio of the height of the water table to the deposit thickness inside the source areas) was set to 1.0, based on the slope stability analysis of the tephra-fallout deposit that would accumulate on La Fossa cone as a result of a Vulcanian eruptive cycle (Baumann et al. 2019). Interestingly, these values of h_w^rel and p_w^rel are also similar to those given by Cuomo et al. (2016). Finally, the consolidation coefficient (c_v) and the erosion coefficient (K) were varied within ranges taken from literature (10^-3 to 10^-2 for c_v and 0.019 to 0.007 for K) (Cuomo et al. 2016). In order to be conservative in our modelling approach, the range of K was extended by one order of magnitude (i.e. down to 0.0007). As a result, 8 scenarios were identified (Table 1). The cases of fresh deposits (loose material with low friction angle) are considered in the scenarios V3 and V6, where small and high entrainment capability of the flow is also considered, respectively. On the other hand, the same high entrainment capability of the flow is considered in the scenarios V1 and V2, but with different shear resistance of the flow at the base, being the two cases representative, respectively, of old (dense materials with high friction, deposited some time ago and consolidated along the slope) or intermediate age deposits (Table 7 in Appendix A). It is important to note that each scenario is associated with one deposit remobilisation event which results in individual lahar from different catchments and channels. The material is released from 22 polygons (catchments), in the north site of the La Fossa cone (Fig. 3), all at once in the three considered scenarios with an initial volume of 24,857 m³. Each scenario was ran for 60 seconds, which represent the duration for the lahars to reach the maximum runout regardless of the scenario. A time step of 0.1 s was and the output results plotted at each 5 s. Thicknesses <1 mm were not computed.

Table 1 Key lahar rheological parameters used in the numerical simulations with SPH for 8 different scenarios (refer to Table 7 in Appendix A for the associated physical meaning of the combination of the various parameters). ϕ_b: basal friction angle of the propagating flow, h_w^rel: relative water height, p_w^rel: ratio of pore water pressure to liquefaction pressure; c_v: consolidation factor, K: empirical parameter for the bed entrainment law of Blanc et al. (2011). * The calibration was performed by using some results of Cuomo et al. (2016) and Baumann et al. (2019)

Step 2: Exposure analysis

The elements at risk considered in this study are residential buildings, critical infrastructures and facilities in the northern sector of the island (Fig. 1). The exposure was estimated in ArcMap as the intersection of buildings and infrastructures with the inundation extent of the various scenarios. In fact, various propagation outcomes are expected for the different scenarios because the lahar runout depends on the associated rheology (e.g., basal friction angle, consolidation coefficient or erosion parameter) and topography.

Step 3: Physical vulnerability of residential buildings and infrastructures

A semi-quantitative physical vulnerability analysis of was performed using the parameters proposed by Papathoma-Köhle et al. (2012), each of which was assigned with a vulnerability score (Table 2; low = 0.25, medium = 0.50, high = 0.75, very high = 1.0). These scores are then used to calculate the total physical vulnerability index for each building obtained by summing the single score of each parameter. The characteristics of exposed buildings were obtained based both on a field survey carried out in 2017 and on a remote survey using Google Earth and 3D Street View images (e.g. Potere 2008; Dell’Acqua et al. 2013; Papathoma-Köhle et al. 2017). The existence of a protection measure such as a concrete wall around the building was assumed to decrease the vulnerability from 1.0 (very high vulnerability) to 0.25 (low vulnerability). At the same time, the existence of a basement with an underground floor increases the vulnerability from 0.5 to 1.0. The level of maintenance of the building ranges from poor (high vulnerability; 1) to good (low vulnerability; 0.25). Very high vulnerability (1.0) was assumed for those structures built in perpendicular direction with respect to the lahar flow direction, with values comprises between 61° and 90° and high vulnerability (0.75) for values varying from 30° to 60°. The height and the number of floors were measured in the field for most of the critical (i.e. at the foot of the volcano) investigated buildings while, for the buildings investigated based on remote strategies, height was retrieved from the difference between the DSM and the DTM. At the base of the Pietre Cotte lava flow, the height of protective walls and of the drainage channel, close to the houses was measured in the field.

Table 2 Physical vulnerability score parameters for buildings

Step 4: Risk assessment

Various strategies of risk assessment exist in literature that include qualitative, semi-quantitative or quantitative analyses (e.g. Lee and Jones 2014; Liu et al. 2015; Van Westen and Greiving 2017; Bonadonna et al. 2021; Poljansek et al. 2021). In particular, quantitative analyses require a quantitative hazard assessment as well as case specific fragility functions that can relate the hazard intensity to the potential impact (Thouret et al. 2020). In addition to fragility functions, physical vulnerability can also be assessed based on vulnerability indicators and damage categories (e.g., Jakob et al. 2012; Wilson et al. 2014; Godfrey et al. 2015; Papathoma-Köhle et al. 2019). Due to the lack of vulnerability functions specific to the Vulcano case study, we opted here for two alternative approaches of risk assessment: an indicator-based risk assessment (semi-quantitative) and an exposure-based risk assessment (qualitative). For both risk approaches, the hazard of lahars is based on a scenario-based hazard assessment of the initial source deposit (i.e., Biass et al. 2016; Baumann et al. 2019).

Here, we consider two HIM that are lahar-flow-front velocity and deposit thickness. Although other HIM such as the dynamic pressure can arguably better capture the dynamics of impact, associated fragility curves require detailed characteristics of exposed infrastructures that are currently unavailable for Vulcano. In addition, flow velocity and deposit thickness are the most common HIMs used to relate lahar hazard to building damage in post-event impact assessments (e.g., Fuchs et al. 2007; Akbas et al. 2009; Totschnig et al. 2011; Lo et al. 2012; Papathoma-Köhle et al. 2012; Totschnig and Fuchs 2013; Kang and Kim 2016; Zhang et al. 2018; Thouret et al. 2020). As a result, our approach produces a first-order estimation of the potential impacts of lahars on the built environment using the fragility curves proposed by Wilson et al. (2014). To despite being inherently associated with large uncertainties, this approach reflects the difficulty of quantifying volcanic risk in data-poor contexts. Although the development of more accurate vulnerability models is required, this is beyond the scope of the present paper.

The risk associated with the lahar-flow-front velocity is assessed following a semi-quantitative (indicator-based) approach based on the following equation:

$$Risk= Hazard\times Exposure\times Vulnerability$$

(5)

Lahar flow velocities (i.e. Hazard) are classified in a score ranging from 1 (low hazard) to 4 (high hazard) (Table 3) based on the review of observed damage on buildings from eruptions in the last 100 years provided by Wilson et al. (2014). Exposure is the Boolean presence of buildings inside a given hazard zone. Considering that the study is done at the building-scale, therefore, the exposure is equal to 1. The vulnerability value is obtained for each exposed building from the total vulnerability index (Table 2). The risk results from the multiplication between the hazard scoring, the vulnerability index and the exposure.

Table 3 Hazard scores based on the velocity categories of Wilson et al. (2014)

Finally, the exposure-based risk associated with the thickness of the final lahar deposit is also compiled based on the burial fraction of the exposed buildings in the lahar inundation area using the following formula:

$$Fraction\ of\ buried\ building\ height=\frac{l_{th}}{h_b}\times 100$$

(6)

with h_b being the building height and l_th the lahar deposit thickness.

Results

Lahar hazard assessment

Lahar inundation areas and deposit characteristics (scenarios V1, V3 and V6)

Out of all the numerical results, those obtained for the scenarios V1, V3 and V6 (Table 7 in Appendix A) are discussed in detail. In particular, V1 has the minimum porosity (0.5) combined with the maximum internal friction angle (40°; i.e. tanφ_b=0.4), which represents old materials and high entrainment capability. V3 has the maximum porosity (0.7) with the minimum internal friction angle (30°; i.e. tanφ_b=0.2), which represents freshly deposited materials and high entrainment capability. V6 represents an intermediate case with the maximum porosity combined with the minimum internal friction angle, or the minimum porosity combined to the maximum internal friction angle; in both cases tanφ_b=0.3, and also associated with low entrainment capability.

The maximum runout distance for V1 scenario is 555 m on the west side of the La Fossa cone, next to the Palizzi Valley (Fig. 4A and Table 4). The remobilized lahar volume, at the end of the 60 seconds of simulation is 260,122 m³, one order of magnitude higher than the initial volume (24,857 m³) with an averaged lahar-deposit thickness of 2.3 m. The V3 scenario is characterized by the highest value of remobilised material (481,570 m³), with an average lahar deposit thickness of about 3.0 m (Table 4). The maximum runout distance of 694 m is shown at the southern sector of the La Fossa cone, next to the Palizzi Valley (Fig. 4B). The V6 scenario is associated with the highest runout distance of about of 717 m, with a final remobilized lahar volume of 41,924 m³ and an average lahar deposit thickness of 0.3 m (Fig. 4C and Table 4).

Fig. 4

Final (at 60 seconds) lahar deposit thickness and runout (m; indicated in white boxes). A): scenario V1; B): scenario V3; C): scenario V6 (Table 5)

Table 4 Main results for the three simulated scenarios (V1, V3 and V6 in Tables 1 and 2). V_in and V_fin indicate initial and final volume

Lahar-flow velocity and deposit thickness distribution (scenarios V3 and V6)

The lahar-flow velocity was investigated for scenarios V3 and V6 only (in four time-steps each 5 seconds), as the lahar runout associated with scenario V1 does not reach the built environment (Fig. 4). The highest percentage of simulated SPH nodes with velocities >5 m/s in V3 scenario was reached at 20 s (93.5%) (Fig. 5a-d and Table 8 in Appendix A). The V6 scenario is associated with higher lahar-flow velocities with respect to the V3 scenario (Fig. 6a-d and Table 8 in Appendix A). The highest percentage of simulated SPH nodes with velocities >5 m/s was reached at 15 seconds (98.6%).

Fig. 5

Lahar hazard maps based on flow velocity classes (according to the strategy of Wilson et al. (2014 modified); Table 3) for the V3 scenario at 20 sec (a), 25 sec (b), 30 sec (c) and 35 sec (d) of the SPH simulation. The dashed lines show the maximum runout for V3 scenario; the main critical infrastructures are also showed (i.e. telecommunication and electrical power plant)

Fig. 6

Lahar hazard maps based on flow velocity classes (adapted from the strategy of Wilson et al. (2014); Table 3) for the V6 Scenario at 15 sec (a), 20 sec (b), 25 sec (c) and 30 sec (d) of the SPH simulation. The dashed lines show the maximum runout for V6 scenario; the main critical infrastructures are also showed (i.e. telecommunication and electrical power plant)

The thickness of the final lahar deposit (60 seconds of the simulation) has also been investigated for the two scenarios V3 and V6 (Fig. 7). The thickness of the final deposit at the base of the cone where houses are located is larger than the flow thickness at any previous time step (Figures 12, 13 and 14 in Appendix B). Four classes of deposit thickness have been selected in ArcGIS with user–defined classification according to the statistical distribution of the values for V3 scenario (<1.0 m; 1.0-2.5 m; 2.5-7.0 m and >7-0 m), in order to better compare V3 and V6 deposit thickness. For the V3 scenario, the maximum thickness values (>7 m) have been recorded at the base of the two watersheds of the Pietre Cotte lava flow and in a small portion at the W sector of the cone and inside the Forgia Vecchia old crater (Fig. 7a). For the V6 scenario, the maximum value of deposit thickness (2.5 m) has been recorded at the base of the left watershed of the Pietre Cotte lava flow, while the majority (95.2%) show thickness <1.0 m (Fig. 7b). The thickness distribution of final lahar deposits suggests that large thickness values (>7.0 m) can be considered as outliers (probability below 5%) or model artefacts due to DSM limitations, that could derive from unexpected holes or also small “dams” crossing drainages which are artifacts of the gridding.

Fig. 7

Final lahar deposit thickness at the end of simulations (60 seconds) for a) V3 scenario and b) V6 Scenario. The dashed lines show the maximum runout for V3 and V6 scenarios; the main critical infrastructures are also shown (i.e. telecommunication and electrical power plant)

Exposure analysis and physical vulnerability of buildings

58 residential buildings, 2 critical infrastructures (electrical power plant and telecommunication), 2 km of main roads and 3.7 km of secondary roads are located in the inundation areas associated with scenarios V3 and V6 (Fig. 8). Given the homogeneity of building characteristics on Vulcano (Biass et al. 2016), we combined results from both field and remote surveys. Out of the 60 buildings identified (including the 2 critical infrastructures), 20 buildings have, therefore, been analysed with direct field observations (from the ground-based survey carried out in 2017) and 40 based on remote survey using Google Earth and 3D Google street view (Fig. 8).

Fig. 8

Exposure map showing both buildings and infrastructures that are located in the lahar inundation areas for the three selected scenarios. The vulnerability of buildings was assessed either in situ (light violet) or remotely (with Google Earth; dark violet). The dashed lines show the maximum runouts for V1, V3 and V6 scenarios

20 of the residential buildings investigated in situ in October 2017 are inside the lahar V6 inundation area while 9 are inside the V3 inundation area. A total of 40 and 29 exposed buildings (two critical infrastructures included) have been investigated based with Google Earth and 3D street view for the V6 and V3 scenarios, respectively. The maximum, minimum and average score of the vulnerability index for both scenarios are 8.00, 3.75 and 5.45, respectively. Three classes of vulnerability have been selected in ArcGis through the Jenks natural breaks method (3.75-5=low; 5-6.25=medium; 6.25-8=high) (Fig. 9). A total of 25 buildings shows a low vulnerability index (42%), 25 buildings show a medium vulnerability index (42%) and 10 buildings have a high vulnerability index (16%). The height of all investigated buildings was obtained from the difference between the DSM and the DTM and ranges from 3.0 to 6.2 m, with an average value of 3.6 ± 0.15 m. The protection walls next to the buildings, below the Pietre Cotte lava flow and along the main road to the port area, have a minimum, maximum and medium height of 0.8 m, 1.9 cm and 1.3 cm, respectively; the wall around the drainage grids have a height of 2 m (Fig. 2). The protection walls around the two critical infrastructures (telecommunication and electrical power plant) are 0.8 m high. In particular, the electrical power plant is located next to the main road at a lower level with respect to the road (i.e. on a downhill slope) and the protection wall does not surround the whole structure (i.e. a metallic lattice gate is facing the main road and the volcanic slopes).

Fig. 9

Vulnerability index associated with the buildings located in the lahar inundation area for the V3 and V6 scenarios

Risk assessment

Semi-quantitative lahar risk assessment related to lahar flow velocity (scenarios V3 and V6)

The lahar risk assessment related to the impact of lahar-flow velocities has been carried out based on the SPH simulations (Figs. 5 and 6) considering their maximum values at first arrival on buildings and combined with the building vulnerability assessment described in the previous section (Fig. 9). Risk classes have been defined based on the Jenks natural breaks algorithm according to the distribution of values for V6 scenario; the same classes have been considered for V3 scenario.

For the V3 lahar scenario (Fig. 10a), a total of 37 buildings are concerned (Table 5). Most of the investigated buildings (78.4%) show a medium-risk value, while 13.5% and 8.1% show a low-risk and a high-risk value, respectively (Table 5). The main affected buildings associated with high risk are mostly located on the left side of the main watershed of the Pietre Cotte lava flow, where the highest lahar impact velocities are observed (Fig. 5). The telecommunication infrastructure is at medium risk because of high impact lahar velocity and for the inadequate protection walls around it (i.e. low wall height), while the electrical power plant is not reached by the associated lahar flow.

Fig. 10

Lahar risk map for the Vulcano Porto area based on impact velocities (m/s) for the V3 (a) and V6 (b) scenarios

Table 5 Affected buildings in relation to lahar impact velocity associated with scenarios V3 and V6

In the V6 lahar scenario (Fig. 10b), most buildings (38 out of 60; i.e. 63.3%) show a medium risk and they are mainly located at the base of the Pietre Cotte lava flow in both watersheds; 11.7% show a high risk (7 buildings), which are located mainly at the base of the La Forgia Vecchia and at the western side of the La Fossa Cone; and 25% (15 buildings) show a low risk (Table 5). The telecommunication infrastructure and the electrical power plant show a medium risk value (Fig. 10b).

Qualitative risk assessment of lahar related to the final deposit thickness (scenarios V3 and V6)

Three main classes of risk have been defined to classify the potential impact based on the thickness of the final lahar deposit (Fig. 11): a low-risk class, with a fraction of buried building height <25%; a medium-risk class, with a fraction of buried building height between 25-50%; and a high-risk class, with the fraction of buried building height >50%.

Fig. 11

Percentage of buildings buried by the final lahar deposit for the V3 and V6 Scenario (i.e. building height / lahar deposit thickness × 100). Low-risk class: burial fraction < 25%; medium-risk class: burial fraction 25-50%; high-risk class: burial fraction > 50%. The thickness of deposit is also indicated with colour points. The dashed lines show the maximum runout for V3 and V6 scenarios; the main critical infrastructures in this area are also showed (i.e. telecommunication and electrical power plant)

Half of the affected buildings in scenario V3 (18 buildings, i.e. 49%) are in the high-risk class, with a mean lahar-deposit thickness of 3.9 m and an average building height of 3.6 ± 0.15 m. 30% of the buildings (11) are in the medium-risk class, with a mean deposit thickness of 1.4 m and the 21% (8) are in the low-risk class with an average lahar-deposit thickness of 0.5 m (Table 6). A length of 850 meters of the main road and 1.6 km of the secondary roads are affected by lahar in this scenario, with the maximum deposit thickness recorded at the main road “Strada Provinciale 179” and at the secondary road “Via sotto Cratere” located at the foot of Pietre cotte lava flow with an average value of 2 m (Fig. 11).

Table 6 Affected building height (%) for lahar-deposit thickness associated with V3 and V6 scenario

For the V6 Scenario (Table 6), a total of 58 buildings (96.6%) with an average height of 3.6 ± 0.15 m are in the low-risk class, with an average deposit thickness of 0.2 m. Only one building located on the north side of the Pietre cotte lava flow is in the high-risk class with a mean lahar-deposit thickness of 1.8 m. A building with a height of 3.7 m is in the medium-risk class, with a deposit thickness of 1.3 m. A total of 1.1 km and 2.1 km of the main and secondary roads respectively are affected by lahar. The maximum lahar deposit thickness (up to 2.1 m) is recorded in one segment of the main road “Strada Provinciale 179”, at the secondary road of “Via Sotto Cratere” and below the Forgia Vecchia (Fig. 11).

Discussion

Lahar rheology and dynamics

Due to the intense rainy seasons and to the physical characteristics of the pyroclastic deposits on Vulcano, lahars represent a common process in the stratigraphy of La Fossa cone (e.g., Di Traglia et al. 2013). Nonetheless, lahar deposits on Vulcano are difficult to map and correlate as they commonly remobilize material of the same nature and are frequently buried or eroded by following events. The only lahars that can be mapped are those in individual channels. In fact, several lahar deposits covering the 1888-90 primary tephra fallout deposit have been observed in stratigraphic sections located in channels on the NW flank of the la Fossa cone (Baumann et al. 2019). However, these lahars were remobilized on a topography different from the existing one and cannot be used as a back analysis. Despite the difficulty to correlate the modelling results with the 1888-90 lahar deposit thickness, we obtained final lahar thicknesses (between 0 and 1 m) for the scenario V6 which are similar to the 1888-90 lahar deposits measured in the field (between 0.2 and 1 m) (Baumann et al. 2019).

As lahar models could not be accurately calibrated and validated, eight scenarios have been selected in this work to best represent the lahars associated with a potential Vulcanian eruption at La Fossa volcano of similar magnitude as the 1888-90 event (Table 1). These eight scenarios are characterized by end members of the model input parameters to be expected for Vulcano, especially in relation to the basal friction angle of the propagating flow, the relative water height, the ratio of pore water pressure to liquefaction pressure, the consolidation factor, and the empirical parameter for the bed entrainment law (Table 1). The main geotechnical parameters, in particular h_w^rel and p_w^rel, are related to the water pressure in the lahar source at the triggering stage, derived from the previous studies on La Fossa lahars (Baumann et al. 2019); the other parameters are taken from the studied lahars of the Vesuvius area (Sarno, Nocera and Bracigliano; Cuomo et al. 2016). In particular, V1 is representative of an old tephra deposit consolidated along the slopes due to wetting-drying cycles (high tanφ_o), V3 is representative of a fresh tephra deposit associated with high entrainment capability (high K), while V6 scenario is representative of a lahar generated by fresh tephra deposit associated to a low entrainment capability (low K) (Table 7 in Appendix A). Both V3 and V6 can, therefore, be considered as representative of syn-eruptive lahar scenarios. In particular, the volume at the end of V3 and V6 simulations are 481,570 m³ and 41,924 m³, respectively (with respect to an initial volume of 24,857 m³, Tables 1). Runout distances are similar for V3 and V6 (i.e., 694 m and 717 m, respectively), even though V3 is associated with thicker lahar deposits (i.e. average deposit thickness of about 3.0 m for V3 and 0.3 m for V6) and V6 is associated with higher flow velocities (Fig. 4, Table 4 and Appendix A). Given the associated large volumes, these flows are very different with respect to the volcanoclastic debris flows occurring every year during heavy thunderstorms that are characterised by low mobility and small volumes (20.5 m³ - 51.3 m³) with 60-100 m of runout on the La Fossa cone slopes (Ferrucci et al. 2005).

Risk assessment

The primary damage mechanism of lahars is the increased dynamic pressure which overcomes the structural design causing structures to fail; however, disruption can also occur when there is insufficient dynamic pressure to cause physical damage, but deposition occurs (e.g., Baxter et al. 2005; Rheinberger et al. 2013; Thouret et al. 2020; Wilson et al. 2014). In particular, flow density, velocity, depth, and the least angle between the flow direction and an obstacle represents the most meaningful predictors of relative damage and proportional loss (e.g., Rheinberger et al. 2013; Thouret et al. 2020). It is important to stress that the damage categories of Wilson et al. (2014) related to lahar velocities are based on the fragility functions of Zuccaro and De Gregorio (2013) and are in good agreement with damage categories associated with empirical analysis that indicate total building destruction for flow velocities >5 m/s (Zhang et al. 2018). In addition, according to the damage classification of Hu et al. (2012) for debris flows, a building (brick-concrete or with a reinforced-concrete frame) can be considered as completely damaged if two or more storeys are buried or seriously damaged if more than one storey is buried. Two strategies have, therefore, been adopted here to assess the risk associated with lahars on Vulcano. One associated with the maximum impact velocity of the first lahar arrival on the buildings (Fig. 10), and a second one associated with the potential burial of buildings (Fig. 11). Lahars associated with the V1 scenario do not reach the base of the La Fossa cone and, therefore, are ignored from the risk assessment. Based on our results, V3 represents the worst-case scenario in terms of impact associated with building burial by the final lahar deposit, with 49% of the investigated buildings falling in the high-risk class and the 30% in the medium-risk class. Nonetheless, V6 represents the worst-case scenario in terms of lahar-flow velocities, with 11.7 % of the buildings in the high-risk class and 63.3% in the medium-risk class.

The majority of the investigated buildings located at the base of the La Fossa cone in the Vulcano Porto area are associated with a medium-risk class in terms of lahar-flow velocity (i.e. 78.4% for the V3 scenario and 63.3% for the V6 scenario; Fig. 10). These are mostly located at the base of the Pietre cotte lava flow, at the WNW sector of the volcano and below the Forgia Vecchia old crater. They are affected by lahars after 20 seconds from simulation onset with a maximum value of 11 m/s for V3, and at 15 seconds of simulation with a maximum value of 21 m/s for V6. The buildings associated with a high-risk class for both scenarios represent 8.1% (3) and 11.7% (7) for V3 and V6, respectively, and are mainly located in the W sector for V3 (Fig. 10a) and on the N sector for V6 (Fig. 10b). 13.5 and the 25 % are in the low-risk class for the V3 and V6 scenario, respectively (Fig. 10; Table 5).

In terms of risk associated with the lahar deposit thickness (Fig. 11 and Table 6), 49% of the investigated buildings are in the high-risk class for V3 scenario with a mean deposit value of 3.9 m with an averaged building height of 3.6 ± 0.15 m. These are mostly located at the right side of the Pietre cotte lava flow and in the western sector of the La Fossa cone. In the V6 scenario, only a building at the north-side base of the Pietre cotte lava flow is in the high-risk class, with a mean deposit thickness of 1.8 m with a building height of 3.2 m. The buildings associated with the medium-risk class are the 30% (11) for the V3 scenario and only one for the V6. 66.7 % of the total (40 buildings) for the V6 scenario and 14% (5 buildings) for the V3 are in the low-risk class. 30% and 8% of the investigated buildings were not affected by the final lahar deposit for the V6 and V3 scenario, respectively.

Risk to infrastructures (telecommunication, power plant, road network)

Only two strategic infrastructures are exposed to V3 and V6 scenarios: the telecommunication service and the electrical power plant (the only two infrastructures of this type present on the island; Fig. 1). The telecommunication infrastructure is associated with a medium-risk class related to lahar-flow impact velocity for both V3 and V6 scenarios (Fig. 10), as the first arrival is associated with a value of 18 m/s at 20 seconds for V6 and with a value of 11 m/s at 30 seconds for V3. However, it is associated with a medium-risk class associated with the final lahar-deposit thickness in the V3 case (Fig. 11a), where the lahar is partially stopped by the protection walls on the left side of the structure and the deposit reaches a maximum value of 1.3 m. The same infrastructure falls in the low-risk class of the final lahar deposit thickness for the V6 scenario with a maximum thickness value of 0.1 m (Fig. 11b). Considering the high impact velocities (> 5 m/s) at the telecommunication infrastructures in both scenarios, a complete loss of functionality with significant structural damage is expected. According to Wilson et al. (2014), for these hazard classes a replacement or a financially expensive repair could occur due to the destruction of ground level component (e.g., lines, cabinets, exchanges).

The electrical power plant, on the other hand, is only affected by the V6 scenario and is associated with a medium-risk class, with a lahar-flow velocity of 11 m/s reached at 25 seconds of the simulation (Fig. 10b). In relation to the lahar-deposit thickness, this infrastructure is associated with a low-risk class (0.25 m deposit thickness) (Fig. 11b). In this case, being located on a downhill slope next to the main road with both the main building and the protection walls being lower with respect to the road level, the lahar could go through the lattice gate and reach a maximum thickness value of 0.25 m around the whole building, with a mean thickness of 0.2 m.

Lahars have also been found to affect a total road length of 850 m for V3 scenario (with mean thickness of 2 m) and 1.1 km of road length for the V6 scenario (with a mean thickness of 0.2 m) (Fig. 11). In particular, lahars associated with these two scenarios are expected to affect the main road where the telecommunication and electrical power plant are located and the secondary road at the base of the Pietre Cotte lava flow. The main road (Strada Provinciale 179) connects the Vulcano Porto area to Vulcano Piano and to the emergency harbour at Gelso (on the south of the Island) (Fig. 1), and, therefore, has a strategic role in rescue and emergency operations. The secondary roads at the base of La Fossa cone are affected in both scenarios for a total of 1.6 km and 2.1 km for V3 and V6 respectively, causing their potential loss of functionality (Fig. 11a-b).

Conclusions

Quantifying the hazard and the risk associated with lahar inundation is complex. Multiple generation mechanisms (e.g., deposit instability, sheet erosion), interactions with other phenomena (e.g., properties of the source tephra deposit, rainfall duration and intensity) and evolving rheological behaviour resulting from interactions with surfaces and topography during the flow represent as many parameters that require to be parameterized as input to numerical models. Additional technical (e.g., accuracy of topographic models) and engineering (e.g., choice of appropriate geotechnical and hydrological parameters) aspects combined with the range of potential impact mechanisms (e.g., inundation from hyper-concentrated flows to burial and destruction by debris flows) further complicate deterministic choices required when attempting to estimate lahar impacts. In addition, the hazard metrics used to assess the potential impact (e.g., lahar flow velocity, deposit thickness, dynamic pressure, static pressure) are often case specific and require more studies to be widely applied. As a result, both hazard (inundation areas, flow velocities and deposit thickness) and impact data (related both to lahar-flow velocity and lahar deposit thickness) presented in this study are associated with uncertainties and should thus treated with caution. However, this work represents the first attempt to quantify the potential impact of syn-eruptive lahars as a pipeline of interlinked processes, including the generation of the source material, the parameterization of conditions leading to instability to estimate source conditions, the modelling of the lahar flow, the vulnerability of the associated exposed buildings and infrastructures up to the estimation of physical impacts. Several deterministic choices were required to achieve this process (e.g., choice of eruption scenarios, rainfall properties, probability thresholds), which limit the direct generalisation of the proposed methodology to other case studies; nevertheless, our study well illustrates various critical aspects that need to be addressed in pre-event lahar impact assessments. The resulting semi-quantitative and qualitative risk analyses (for lahar flow velocity and deposit thickness, respectively) can help prioritize the buildings and/or infrastructures for which a detailed in-situ physical vulnerability assessment could be carried out in order to compile a quantitative risk assessment. In particular, we found that the level of potential impact depends on the selection of the hazard metrics. As an example, we found that the syn-eruptive lahar scenario with medium entrainment capability (V3) produces the highest impact associated with building burial by the final lahar deposit (nearly 80% of the total number of buildings analysed are in the medium- to high-risk class). In contrast, the syn-eruptive lahar scenario with low entrainment capacity (V6) results in the highest impact associated with lahar-flow velocities (45% of the total number of buildings analysed are in the medium- to high-risk class). In addition, a key critical infrastructure on Vulcano (telecommunication) is affected by both scenarios (V3 and V6) (Figs. 13 and 14), while the power plant is only affected by the syn-eruptive lahar scenario with low entrainment capacity (V6). The critical road that connects the two key inhabited areas of the island (Piano and Porto) would also be inundated in case of both scenarios. Clearly, locating the only telecommunication infrastructure and the main power plant present on the island at the base of La Fossa cone increases the general systemic vulnerability of the island (e.g., Galderisi et al. 2013). The presence of residential buildings in the inundation areas of syn-eruptive lahars increases the potential of loss of lives and/or significant structural damage with associated large economic cost. These results demonstrate that considering the potential impact of volcanic hazards in the urban and spatial planning of the island could represent a valuable strategy to reduce the risk in the long term. However, in order to develop an accurate and efficient land-use plan, further studies should also include the assessment of additional dimensions of vulnerability, such as systemic vulnerability related to both cascading effects associated with the loss of functionality of infrastructures (i.e., telecommunication and power plant) and to reduced accessibility due to the inundation of primary roads connecting the key areas on the island.

Appendices

Appendix A - Additional features of selected lahars scenarios

Table 7 Qualitative description of the simulated scenarios and their outputs (refer to Table 1 for associated rheological parameters). In bold the three scenarios considered for hazard assessment in this work

Table 8 Percentage of simulated SPH nodes in each lahar-flow velocity hazard class (Table 4) at each investigated time step for V3 and V6 scenario

Appendix B - Variation of lahar flow thickness with time at control points

Fig. 12

Location of the Control Points used to assess the variation of lahar-flow thickness at different times steps for V3 and V6 scenarios

Fig. 13

V3 lahar flow thickness between 20 and 60 seconds of simulation at Control Points (CP) shown in Fig. 12

Fig. 14

V6 lahar flow thickness between 15 and 60 seconds of simulation at Control Points (CP) shown in Fig. 12