Lahar Risk Assessment: the case study of Vulcano Island, Italy

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 absence of a hazard-based land-use planning. Their destructive power is mostly associated with their velocity (related to flow rheology and surrounding topography) and to their ability to bury buildings and structures (related to the deposit thickness). The distance reached by lahars depends on their volume, on sediments/water ratio, as well as on the overall characteristics of the path where they propagate. Here we present a novel strategy for the assessment of risk associated with lahar inundation related both to flow velocity and deposit thickness using Vulcano island (Italy) as a case study. First, a range of hazard scenarios has been identified that are related to the mobilization by intense rain events of tephra fallout deposited on the slopes of the La Fossa cone by a future Vulcanian eruption. Second, a numerical model has been used to identify the potential lahar impact areas on the northern sector of Vulcano, where both residential and touristic facilities are present. In this specific case we have used the Smoothed Particle Hydrodynamic (SPH) model that provides information on both flow velocity and deposit thickness. Finally, exposure and vulnerability surveys were carried out in order to compile risk maps for both lahar-flow velocity and final lahar-deposit thickness. Our analyses show the importance of carrying out accurate and detailed risk assessments exploring a variety of initial conditions in order to best quantify the potential damage and identify suitable mitigation strategies.


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 impact settled areas nearby active and dormant volcanoes at different spatial and temporal scales. In this framework, lahars represent one of the most important volcanic hazards that can potentially struck areas up to hundreds of kilometres from active volcanoes and kill a large number of people as demonstrated by the worst event of the 20 th century that occurred in 1985 at Nevado del Ruiz volcano (Colombia) (Lowe et al. 1986). Lahar is an Indonesian term used by modern volcanologists to describe mudflows or water-saturated debris flows.
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). The mixture is primarily formed by pyroclastic materials (tephra fallout, loose PDCs or older volcaniclastic sediments) and various sources of water (e.g. rain, melted ice caps and snow, lake water). Lahars can be produced both during (primary) and after (secondary) eruptive events 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. Jenkins et al. 2015;Wilson et al. 2014;Mead et al. 2017).
The deadly lahars associated with the 1985 eruption of Nevado del Ruiz represent an example of primary lahars generated by the rapid melting of snow and ice on volcanoes due to the interaction with eruptive products (Major and Newhall, 1989). In this case, a moderate explosive eruption melted the volcano ice cap generating a large lahar that inundated areas up to 100 km downstream destroying more than 5,000 houses and killing more than 23,000 people (Mileti et al. 1991). The Armero village, 73 km downstream of the volcano, was completely destroyed and 3/4 of inhabitants were killed.
As an example of secondary lahars, large flows were generated by heavy rainfall on freshly deposited tephra (fallout and PDC deposits) after the 1991 eruption of Mt. Pinatubo (Philippines) (Rodolfo 4 1989; Pierson et al. 1992;Arboleda and Martinez 1996;Rodolfo et al. 1996). Between June 12 th and September 10 th , about 140 lahars affected a surface of 400 km 2 around the volcano and caused the displacement of more than 50,000 people. The Pinatubo aftermath clearly demonstrated how lahars can create a persistent and chronic hazard in time, forcing communities affected to develop adaptation measures that include relocation and multi-level structures. A second example of secondary lahars is that of Casita volcano in the NW sector of Nicaragua, where, in 1998, heavy rainfall associated with the Hurricane Mitch contributed to the mobilisation of fractured and hydrothermally altered material resulting in a lahar that killed approximately 2,500 people (Kerle et al. 2001). More recently, the material deposited by the 2010 VEI 4 (Volcanic Explosivity Index) eruption at Merapi volcano Debris flows composed of volcaniclastic materials can also occur hundreds of years after volcanic activity due to heavy rain events. An example is given by the May 1998 events of Pizzo d'Alvano massif (Italy), located about 10 km east of Mt. Vesuvius volcano. In this case, after some days of intense rains, the upper slope of the Appenine belt, mantled with up to five meters of loose volcanoclastic deposits (ash and lapilli deposited during the 1800 BC and 472 AD activity of Vesuvius) was mobilized and impacted 4 municipalities, causing 159 fatalities and heavy damage to buildings and properties (Cascini 2004).
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 speeds exceeding 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). In medial to distal areas, lahars are typically confined in river valleys and flood plains, which favour their long runouts and make them a widespread hazard. Lahars become rapidly less energetic with increasing distance from vent due to a decrease of flow confinement and slope steepness Jenkins et al. 2015).

5
Their associated impacts transition 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;Daga et al. 2018). Buildings can be damaged through a number of mechanisms including: (i) direct damage resulting from static and dynamic forces imposed by the flow, (ii) damage to foundations through erosion and scour, (iii) buoyancy effects of the flow causing structures to float, (iv) direct damage from large debris (missiles) within the flow, (v) burial of buildings and other type of infrastructures under high and widespread thickness of deposits and (vi) indirect damage caused by chemical and biological processes such as seeping induced weakness of mortar (Kelman and Spence 2004;Mead et al. 2017). It is important to note also that the production of numerous lahars can cause the soil to gradually rise over time with many buildings and infrastructures remaining below the existing ground level (e.g. areas around Pinatubo, Merapi, and east of Vesuvius volcanoes). One of the main physical parameters typically considered to assess the potential impact of lahars are the flow depth (i.e. hydrostatic pressure) and the dynamic pressure (derived from lahar velocity and flow density). 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). Wilson et al. (2014) examined the observed impacts and hazard intensity relationships, showing that most infrastructures can be disrupted by the impact of lahars.
Various modelling approaches have been used in literature to simulate the runout path and the inundation area over a Digital Elevation Model (DEM) and implemented by the use of a Geographical Information System (GIS). Some of these models are of straightforward application but can only identify inundation areas with no information on flow velocity and thickness (e.g. LAHARZ; Schilling 1998). More recently, models have been developed that can also compute 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 (Williams et al. 2008, LaharFlow 6 (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). The model was able to satisfactorily reproduce all of the major lahar characteristics including trajectory, inundation areas, flow velocity, and deposit depth (Haddad et al., 2010(Haddad et al., , 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 Regardless of the many models developed to describe lahar inundation and the variety of applications to assess the associated hazard, only few studies have investigated the actual vulnerability (e.g. Jenkins et al. 2015;Mead et al. 2017). Nonetheless, the assessment of potential impacts is key to sustainable development plans around volcanic areas. In this paper, we present a new methodology for the assessment of potential lahar impact using Vulcano island (Italy) as a case study. Because the island is periodically affected by intense rainy events, lahar remobilization of tephra deposits is still nowadays an active process (e.g. Di Traglia et al. 2013, Baumann et al. 2019. Galderisi et al. (2013) had already shown how the north of Vulcano island is the most exposed to lahars 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 as the hazard assessment was carried out using the model LAHAR-Z. Here we describe all the steps to assess risk associated with lahar inundation (i.e. future potential impact) based on the assessment of hazard (flow velocity and thickness, here 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 7 southernmost of the seven emerged Aeolian Islands. It is part of the Lipari Municipality and has a surface of about 20 km 2 with a total of 1,282 permanent residents with almost 70% of the population located in the Porto-Vulcanello area and 30% in Piano. In the month of August (peak of the high season), the population reaches a maximum of about 28,000 people due to tourists and seasonal workers. The main economic activities and tourist attractions are located in the northern sector of the Island. Most critical infrastructures and facilities of the island 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). These roads have a total length approximately of 13 km. The secondary roads instead, reach a total length of about 54.8 km.
Vulcano Island is a composite volcanic edifice located at the southern extremity of the central segment of the Aeolian Islands, along the Tindari Letojanni (TL) fault. It has built through four main eruptive stages during which the activity shifted from SE (Primordial Vulcano stratovolcano) to NW (Lentia -Mastro Minico and La Fossa cone), and finally towards North with Vulcanello cones. Two multi-collapsed calderic structures formed during the Vulcano history: the "Caldera del Piano" and "Caldera de La Fossa", where the actual active La Fossa cone is located. La Fossa cone is the current active edifice (a 391 m-high cone with a basal diameter of 1 km) that was built during since 5.5 ka (up to historical times) in the middle of La Fossa caldera through recurrent hydromagmatic to Vulcanian explosive phases of activity, which produced multiple dilute PDCs and tephra fallout deposits, alternating with a few 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 and morphological variations in the cone and its surrounding area was investigated through a stratigraphic reconstruction by Frazzetta et al. (1983), Ferrucci et al. (2005) 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 8 and deposition of small lahars. Lahars occurred at La Fossa cone involving the tephra layers emplaced by the Vulcanian and phreatic explosive events during the last hundreds of years. Stratigraphic sections from exposure along the main creeks around La Fossa cone show several lahars events.
Lahars bed thickness vary 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 lahars generation (Ferrucci et al. 2005). According to meteorological data of the National Institute of Geophysics and Volcanology (INGV, section of Palermo) obtained from a Davis meteorological station located in the Lentia caldera rim, 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 h 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) (Ferrucci et al., 2005;Di Traglia, 2011), and the steep slopes without a significant vegetation cover favour the remobilization of volcanic deposits as lahars (Di Traglia et al. 2013). At Vulcano Island and, in particular in the Porto Levante area, flood-control engineering structures are present that are designed for rainwater collection (Fig. 2). At the base of the La Fossa Cone, along the main watershed on the left 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 engineering 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).

Methods
The main objective of this paper is the compilation of a new strategy for the assessment of risk associated with lahar inundation. In order to demonstrate the applicability of such a strategy, lahar risk assessment maps were compiled that are relative to a future Vulcanian eruption of La Fossa cone 9 of similar magnitude of that of 1888-90. The risk assessment was focused on the northern sector of the island, which is characterized by the highest exposure of residential buildings and critical infrastructures and facilities (Galderisi et al. 2011) (Fig. 1).
First, the inundation area of lahars was determined together with the associated flow velocity and final deposit thickness (i.e. hazard assessment). The lahar inundation area was determined for the Vulcanian scenario associated with the highest potential volume of remobilised material, which and applied in order to define a physical vulnerability score for buildings (i.e. vulnerability assessment). Buildings were analysed based both on a field survey carried out in 2017 and on a remote survey using 3D Google Earth images. Finally, two risk maps were compiled (i.e. risk assessment).
The first risk map related to flow velocity was compiled based on a combination of the assessment of flow velocity at the first arrival at buildings and the vulnerability score of each exposed element (buildings). The second risk map was compiled by comparing the final deposit thickness and the height of buildings to assess the percentage of building burial. The lahar flow velocities, together with the final deposit thickness outputs, were exported into a Geographical Information System (GIS) using a dedicated Matlab script created as part of this effort and available for future use at https://github.com/jrarrowsmith/Process-SPH-output.

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. In particular, Biass et al. (2016) have provided the potential accumulation of tephra on the La Fossa cone in the case of a variety of eruptive scenarios based on TEPHRA2 simulations (i.e. Vulcanian cycle and subplinian eruptions of VEI (Volcanic Explosivity Index) 2 and 3. Baumann et al. (2019) have used the associated probabilistic isomass maps to assess the scenario that would produce the most unstable deposit combining field characterization, geotechnical analysis, and numerical modelling with TRIGRS (Baum et al. 2002). For the Vulcanian scenario, the largest unstable volume was reached for an eruption duration of 18 months, while the subplinian deposits, which are characterised by a smaller cohesion, become more stable with thickness increase depending on rainfall intensity and slope angle (Baumann et al. 2019). Given that a Vulcanian cycle has the highest probability of occurrence at this volcanic system (Selva et al. 2020), in this study we consider the worst-case scenario identified by Baumann et al. (2019) for the Vulcanian cycle (i.e. 18 months of eruption duration with the characteristics of the last cycle occurred on Vulcano in 1888-90). 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). This scenario covers an area of 119,048 m 2 with a maximum tephra-fallout deposit thickness of 0.26 meters resulting in a potential initial lahar volume of 24,857 m 3.

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. (2014Cuomo et al. ( , 2016Cuomo et al. ( , 2018 for modelling the propagation of flows such as water saturated sediments and other mixture of multiphase materials (see also Braun et al. 2017).
The governing equations are reported in Pastor et al. (2009). Here we only highlight that 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 ( ), 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 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: where is the excess of pore water pressure to the hydrostatic value, h is the propagation thickness and cv is the consolidation coefficient.
The consolidation process reduces the pore water pressure while increases the normal stress and the shear resistance at the flow base. The initial pore water pressure is taken into account through two input parameters: ℎ , which is the relative water height resulting from the ratio of water table height to deposit thickness and , 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 landslide propagation has been demonstrated in the literature for many shallow landslide events (Pastor et al. 2009;Xu et al. 2012).
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: where b is the basal shear stress of the flow, '  is the submerged density defined as , 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 ̅ is the depth-averaged flow velocity.
In the equation of mass conservation of the flow, the capability of the flow to entrain further material from the ground surface is also accounted for. The entrainment rate (er) depends on the thickness (h) and 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: 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 (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 Elevation Model (DEM) of the area
The DEM that has been used in this study is derived from a 2-m resolution lidar point cloud acquired in 2017 by Ministero dell'Ambiente (MATTM; www.minambiente.it). For the numerical analysis 13 and in order to minimize the computational cost of SPH, the DEM has been cut in the northern sector of the island with a total of 295,693 cells while, in the source area, 11,832 initial computational points (tephra thickness) were considered with 3 m-spacing, corresponding to an initial volume of 24,857 m 3 , over an area of 119,048 m 2 .
The modelling was 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 (a rhyolitic flow which reached the foot of the slope but did not flow due to high viscosity) and two coalescent craters (Forgia Vecchia) related to phreatic explosions occurred in the XV-XVI centuries.

Selection of SPH Scenarios
In order to reproduce the frictional rheological behaviour of potential lahars associated with a Vulcanian eruptive cycle on Vulcano, we first performed a parametric analysis. The tephra porosity ranges from 0.5 to 0.7, and the internal friction angle spans from 30° to 40°, as measured in 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 gr cm -3 and fluid density as 1 gr cm -3 , we obtained an apparent friction angle (Eq. 2) varying from 11° to 22º, i.e. tan from 0.2 to 0.4. These values are consistent with the range 0.344-0.40 obtained from the back-analysis of past flow-like landslides occurred in the pyroclastic deposits at different sites of the Campania Region (Italy), a few kilometers far from Vesuvius (Cuomo et al. 2016). Furthermore, the relative pore water pressure ( ; i.e. ratio of pore water pressure to liquefaction pressure) was set to 0.4 and relative water height (ℎ ; 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 hw rel 14 and pw rel are also similar to those given by Cuomo et al. (2016). Finally, the consolidation coefficient (cv) and the erosion coefficient (K) were varied within ranges taken from literature (10 -3 to 10 -2 for cv 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 (Tables 1 and 2). The cases of fresh deposits (loose material with low friction angle) are considered in the scenario 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 scenarios representative, respectively, of old (dense materials with high friction, deposited some time ago and consolidated along the slope due to wetting-drying cycles) or intermediate age deposits (Table 2). It is important to note that each scenario is associated with one deposit remobilisation event which results in individual lahars in different catchements and channels.
Each scenario was ran for 60 seconds, which represent the duration for the lahars to reach the maximum runout whatever the scenario; the time step was 0.1 s with the output results plotted at each 5s.
Step 2: Exposure analysis The elements at risk considered in this study are residential buildings and critical infrastructures and facilities in the northern sector of the island (Fig. 1). The distribution of exposed elements depends on the propagation patterns of the lahars associated with the three selected scenarios (V1, V3 and V6) (Tables 1 and 2). Various propagation outcomes are expected for these scenarios because the lahar runout varies depending on the associated rheology (e.g., basal friction angle, consolidation coefficient or erosion parameter) and topography. Thus, the three investigated scenarios provide an overall description of the full range of expected propagation patterns given the selected hazard scenario (i.e. Vulcanian cycle of 18-months duration).

15
Step 3: Physical vulnerability of residential buildings and infrastructures The physical vulnerability of residential buildings has been estimated using the parameters proposed by Papathoma-Kohle et al. (2012), each of which was assigned with a vulnerability score (Table 3; low = 0.25, medium = 0.50, high = 0.75, very high = 1.0). Such a score is then used to calculate the total physical vulnerability index for each building obtained by summing the single score of each parameter ( Table 3) Step 4: Risk assessment The lahar risk has been assessed based on two separate strategies. First, the risk associated with the lahar-flow velocity is assessed following a semi-quantitative approach based on the following equation:

Risk = Hazard × Exposure × Vulnerability
A hazard score was assigned to the hazardous velocity categories identified by Wilson et al. (2014) 16 (Table 4). It is important to note that Wilson et al. (2014) assigned these velocity classes directly to damage of buildings. The first class (< 1 m/s) is associated with no damage; the second class (1-3 m/s) is associated with damage to openings and non-structural elements; the third class (3-5 m/s) is associated with moderate structural damage to walls with some partial collapse; and the fourth class (> 5 m/s) is associated with complete damage to buildings with few structural elements remaining and/or swept off foundations. In addition, the second to fourth class are associated with infiltration of debris and building material into buildings making them inhabitable with the fourth class reaching damages beyond economic repair. These considerations have helped the interpretation of our results.
Nonetheless, in this work, the hazard assessment has been decoupled from the damage assessment, and, therefore, the velocity classes have been assigned a different score (1 to 4) with respect to the damage levels of Wilson et al. (2014) (0 to 3). Exposure represents the location of buildings inside a given hazard zone, which is, therefore, considered with a value of 1. The vulnerability value is obtained at each building inside the hazard zone from the total vulnerability index ( Table 3). The risk results from the multiplication between the hazard scoring, the vulnerability index and the exposure.
Second, the risk associated with the thickness of the final lahar deposit is expressed as the proportion of the building height that could be affected: with hb being the building height and lth the lahar deposit thickness. In this case, the exposure remains, whereas the only vulnerability parameter considered is the height of the building.

Lahar inundation areas and deposit characteristics
Out of all the numerical results, those obtained for the scenarios V1, V3 and V6 ( Table 2) are discussed in detail. In particular, V1 has the minimum porosity (0.5) combined with the maximum internal friction angle (40°), which results in tanb=0.4, thus being representative for the case of old materials and high entrainment capability. V3 has the maximum porosity (0.7) with the minimum internal friction angle (30°), which gives tanb=0.2 to consider the freshly deposited materials, and also associated with high entrainment capability. V6 represents a cross-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 5).The remobilized lahar volume, at the end of the 60 seconds of simulation is 260,122 m 3 , one order of magnitude higher than the initial volume (24,857 m 3 ) with an averaged lahar-deposit thickness of 2.28 meters ( Fig. 4A and Table 5). The V3 scenario is characterized by the highest value of remobilised material (481,570 m 3 ), with an average lahar deposit thickness of about 2.97 m ( Table 5). The maximum runout distance of 694 meters 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 meters, with a final remobilized lahar volume of 41,924 m 3 and an average lahar deposit thickness of 0.29 m ( Fig. 4C and Table 5).

Lahar-flow velocity distribution (scenarios V3 and V6)
The spatial distribution of the lahar-flow velocity was investigated for scenarios V3 and V6 only, as the lahar runout associated with scenario V1 does not reach the built environment (Fig. 4).  Table 6). At 25 seconds the lahar has reached the base of the La Fossa cone, mainly at the N and WNW sectors, with velocity ranging between 6 and 11 m/s and 90.9% of the simulated SPH nodes falling in the highest velocity class (i.e. > 5 m/s; Fig. 5B and Table 6). At 30 seconds simulation, the lahar reached the N sector at the La Forgia Vecchia old crater with velocity ranging between 1.8 and 8.6 m/s, and in the W sector, with velocities between 2.5 m/s and 17.7 m/s.
The lahar flow velocity decreases down to 1.7 m/s at the foot of the Pietre Cotte Lava flow. At this time step, 63% of the simulated nodes is in the highest class ( Fig. 5c and Table 6). At 35 seconds, the lahar decreased its velocity mainly at the N sector of the foot of the La Fossa cone, while in the W sector remains high up to 14.8 m/s next to Palazzi valley. At this time step, most of the SPH simulated nodes are in the lowest velocity class (i.e. 45% in the < 1 m/s class; Fig. 5d and Table 6).
The V6 scenario is associated with higher lahar-flow velocities with respect to the V3 scenario ( Fig.   6 and Table 6). At 15 seconds, 98.6 % of the simulated SPH nodes are in the highest velocity class (i.e. velocities between 5.01 m/s and 27.44 m/s; Fig. 6a and Table 6). At 20 seconds the lahar has reached the base of the La Fossa cone with velocity ranging between 1.7 and 15 m/s in the N and WNW sector, and up to 25 m/s in the W sector; 92.9% of the total simulated SPH nodes fall in the highest velocity class ( Fig. 6b and Table 6). At 25 seconds, the lahar has propagated both in the N sector at the La Forgia Vecchia old crater, with velocity ranging between 2 and 8 m/s, and in W sector with velocities ranging between 4.2 and 19.2 m/s. Lahar velocities decrease at the foot of the Pietre Cotte lava flow with values between 0.03 m/s and 0.88 m/s; 68% of the simulated SPH nodes fall in the highest velocity class ( Fig. 6c and Table 6). At 30 seconds, the lahar decreased its velocity at the base of the La Fossa cone, both in the N and W sector, while inside the watersheds of the La Forgia Vecchia the velocity remains high up to 10 m/s; most of the simulated SPH nodes (45%) fall in the first velocity class (Fig. 6d and Table 6).

Lahar-deposit thickness distribution (scenarios V3 and V6)
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 (Appendix A  (Fig. 7a). Only 6.3% and 3.5 % show a thickness of < 1.02 m and > 7m, respectively. 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, 95.2% of the total nodes show thickness < 1.02 m, while 4.8 % have a thickness between 1.02 and 2.46 m (Fig. 7b). The maximum value of deposit thickness (2.46 m) has been recorded at the base of the left watershed of the Pietre Cotte lava flow. The thickness distribution of final lahar deposits suggests that large thickness values (>7 m) can be considered as outliers (probability below 5%) or model artefacts due to DEM limitations.

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). Out of the 60 buildings identified (including the 2 infrastructure buildings), 20 buildings have been analysed with direct field observations (from the ground-based survey carried out in 2017) and 40 based on remote survey (based on 3D Google Earth) (Fig. 8) (see following section for details).
A total of 20 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 3D Google Earth for the V6 20 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 (3.75-5=low; 5-6.25=medium; 6.25-8=high) (Fig. 9). A total of 25 buildings show 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 ranges from 3.0 to 6.2 m, with an average value of 3.6 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 80 cm, 190 cm and 133 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 85 cm 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).

Risk assessment
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 risk map (Fig. 10a), a total of 37 buildings are concerned ( Table 7). Most 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 7). 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). The electrical power plant is hardly reached by the lahar flow.
In the V6 lahar risk map (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 7). The telecommunication infrastructure and the electrical power plant show a medium risk value (Fig. 10b).

Lahar risk assessment 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 percentage of buried building height < 25%; a medium-risk class, with a portion of buried building height between 25-50%; and a high-risk class, with the portion of buried building height > 50%.
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.5 m. 30% of the buildings (11)  Cratere" located at the foot of Pietre cotte lava flow with an average value of 2 m (Fig. 11). Provinciale 179", at the secondary road of "Via Sotto Cratere" and below the Forgia Vecchia (Fig.11).

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 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-1890 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-1890 lahar deposit thickness, we obtained final lahar thicknesses (between 0 and 1 m) for the scenario V6 which are similar to the 1888-1890 lahar deposits measured in the field (between 0.2 and 1 m) (Baumann et al. 2019).
As lahar models couldn't 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 (Tables 1 and 2). These eight scenarios are characterised by end members of the SPH parameters to be expected on Vulcano, especially in relation to the basal friction angle of the propagating flow (o), the relative water height (hw rel ), the ratio of pore water pressure to liquefaction pressure (pw rel ); the consolidation factor (cv), and the empirical parameter for the bed entrainment law (K) ( Table 1). The main geotechnical parameters, in particular hw rel and pw rel , are related to the water pressure in the lahar source at the triggering stage, 23 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). Out of these 8 scenarios, V1, V3 and V6 scenarios are representative of specific conditions: high tano and K and low cv (V1), low tano and cv and high K (V3), and low tano, cv and K (V6).
V1 is therefore representative of an old tephra deposit consolidated along the slopes due to wettingdrying cycles (high tano; Table 2); V3 is representative of a fresh tephra deposit associated with high entrainment capability (high K; Table 2) while V6 scenario is representative of a lahar generated by fresh tephra deposit associated to a low entrainment capability (low K; Table 2). In particular, the volume at the end of V3 and V6 simulations are 481,570 m 3 and 41,924 m 3 , respectively (with respect to an initial volume of 24,857 m 3 , Tables 1 and 2). 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 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 (e.g. Wilson et al. 2014). However, disruption can also occur when there is insufficient dynamic pressure to cause physical damage, but deposition occurs (e.g. Baxter et al. 2005;Wilson et al. 2014). 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 (Figs. 10 and 11). Based on our results, V3 scenario represents the worst case 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 scenario represents the worst case 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, 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   7).
In terms of risk associated with the lahar deposit thickness (Fig. 11), 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.5 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 25 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 the V6 scenario with a maximum thickness value of 0.05 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 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.19 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 "Strada Provinciale 179" where the telecommunication and electrical power plant are located and the secondary road "Via sotto cratere" 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).

Land-use planning and mitigation measures
The risk maps elaborated in this paper both in terms of lahar-flow velocity (Fig. 10) and final lahardeposit thickness (Fig. 11) represent a useful tool for a sustainable land-use planning in active volcanic areas prone to lahar hazard. In fact, some conclusions can be drawn both for the positioning of residential buildings and infrastructures and on the implementation of efficient mitigation measures. 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 lahars increases the potential of loss of lives and/or significant structural damage with associated large economic cost. A relocation of both these two infrastructures and of the residential buildings currently located in the V3 and V6 inundation areas (total of 37 and 60 buildings, respectively) is, nonetheless, very costly. The implementation of mitigation measures could, therefore, represent a valuable strategy to reduce the risk. Based on our results associated with V3 and V6 scenarios, the existing mitigation measures located at the base of the Pietre Cotte lava flow (Fig. 2) do not seem to reduce the lahar impact. In fact, the large thickness of the final lahar-deposit and high lahar-flow velocities overflow both the channels and the protection walls around the buildings (Figs. 5, 6, 7). In addition, also the protection walls around the critical infrastructures and the exposed buildings (mainly at the base of Pietre Cotte lava flow), have been found to be potentially ineffective for both V3 and V6 scenarios due to their relative low height and for their scarce maintenance status.
In light of our results, the area at the base of La Fossa cone could be affected in case of lahars associated with a Vulcanian cycle, in particular the area at the base of the Pietre Cotte lava flow and at the W side of La Fossa cone. New drainage channels as well as more effective protection walls around residence buildings and critical infrastructures higher than the existing ones with reinforced concrete could intercept the lahar flows. In contrast, an early-warning system due to lahar would not be effective in Vulcano, due to the short time of propagation of the lahars from the top of the cone to