Inundation and evacuation of shoreline populations during landslide-triggered tsunamis: an integrated numerical and statistical hazard assessment

The volcanic island of Stromboli (southern Tyrrhenian Sea, Italy) is renowned for its persistent, periodic, low-intensity explosive activity, whose spectacular manifestations attract tens of thousands of tourists every year. However, sporadic more intense major explosive and effusive eruptions and paroxysms pose serious threats to the island. In addition to direct hazards, granular slides of volcanic debris and pyroclastic avalanches, which can rapidly reach the sea and potentially generate tsunamis, are often associated with such unpredictable eruptive activity. Due to the very fast propagation of the tsunami around the island and the consequent short tsunami warning time (ranging from less than a minute to only a few minutes), mitigation efforts and evacuation from the Strombolian coast must be carefully planned. In this paper, we describe a new GIS-assisted procedure that allows us to combine the outputs of an ensemble of 156 pre-computed landslide-generated tsunami hazard scenarios (with variable landslide volume, position, and density), statistical exposure data (i.e. the number of inhabitants and tourists), and digital geographic information to obtain a quantitative (scenario-based) risk analysis. By means of the analysis of the road network and coastal morphology, we develop a model with routes and times to reach a safe area from every pixel in the inundated area and an appraisal of the time needed to escape versus the wave arrival time. This allows us to evaluate and quantify the effectiveness of potential risk mitigation by means of evacuation. The creation of an impact score linking the predicted inundation extent and the tsunami warning signals is intended, in the long term, to be used to predict the intensity of future tsunamis and to adapt evacuation plans accordingly. The model, here applied to Stromboli, is general and can be applied to other volcanic islands. Evacuating an island hosting several thousand tourists every summer with very little warning time underlines the absolute necessity for such mitigation efforts, aimed at informing hazard planners and managers and all other stakeholders.

1 Introduction

NOAA's NCEI (National Centers for Environmental Information) estimates that about 81 % of recorded tsunamis originate from earthquakes, with 7 % originating from landslides and 5% from volcanoes, the remaining 5 % being from unknown sources (NCEI, 2023; Harbitz et al., 2014). Landslide-generated tsunamis vary greatly in their size and origin, with volcano flank collapse being a frequent source. The range of volumes and the position of the collapse on the volcano flank produce great variation in terms of coastal impact (Grezio et al., 2017). Large-volume collapse of the flank of a volcano island can, for example, generate local tsunamis with waves tens of metres high in proximal fields (Paris et al., 2018). However, tsunamis caused by volcanic landslides are characterised by shorter wavelengths compared to earthquake-induced tsunamis (Grezio et al., 2017), and consequently their far-field effect is often considerably reduced compared to their proximal impact, as observed, for example, during the 2018 Anak Krakatau tsunami (Muhari et al., 2019). However, even a 1 m high tsunami inundating a beach can present a high risk to populated coastlines and tourist seaside resorts.

At Stromboli, an active volcano of the Aeolian Islands (southern Tyrrhenian Sea, Italy; Rosi et al., 2000, 2013), tsunamis can be potentially generated by instability and collapse of volcano flanks (in particular, the steep north-western slope of the Sciara del Fuoco) or by pyroclastic currents generated by paroxysmal eruptions (Pistolesi et al., 2020; Giordano and De Astis, 2020). One example of a tsunami generated by a landslide is that of Stromboli in December 2002. In that case, two landslides generated two tsunamis 7 min apart with maximum run-up of 10.9 m (Tinti et al., 2006a). The first tsunami was due to a nearshore submarine landslide, probably involving up to 20×10⁶ m³ of material (Chiocci et al., 2003), and the second to a subaerial landslide of 4×10⁶–9×10⁶ m³ that broke off at about 500 m a.s.l. (Tinti et al., 2006b). The resulting tsunami caused significant damage along the island coast up to an elevation of about 10 m a.s.l. (Tinti et al., 2006a; Fornaciai et al., 2022), as well as at the island of Panarea 21 km distant (Maramai et al., 2005a).

The risk associated with tsunamis is particularly high on Stromboli due to the high density of population on the coasts, where travel time for the tsunami from the source is minutes (Bonilauri et al., 2021). This is enhanced at Stromboli, where, with the development of global tourism, vacation rentals, restaurants, shops, and hotels have been built close to the beaches. With the development of tourism and influx of seasonal workers, the collective memory linked to the tsunami risk also becomes reduced (Tulius, 2020; Riskianingrum and Yogaswara, 2022).

For these reasons, after the 2002 events, the scientific community and the Italian Civil Protection Department undertook a series of initiatives aimed at quantifying the hazards associated with landslide-generated tsunamis, providing an early-detection and alert system to mitigate the associated risks. In particular, the Laboratorio di Geofisica Sperimentale (LGS) of the University of Florence installed two tsunami detection beacons: the first in 2008 south-west of the foot of the Sciara del Fuoco and the second in 2017 in the north-east (Selva et al., 2021). To be able to generate an automatic alert and warn the relevant authorities and the population as quickly as possible, they have developed a tsunami detection algorithm by studying the STA $/$ LTA (short-time average over long-time average) ratio and the dispersion of surface waves (Lacanna and Ripepe, 2020; Selva et al., 2021). They were able to test the algorithm successfully during the two paroxysms in July and August 2019 (Lacanna and Ripepe, 2020; Selva et al., 2021) and during the avalanche caused by the partial collapse of the NE crater on 4 December 2022.

In parallel, a systematic study has been carried out to define tsunami hazard scenarios based on several complementary approaches. The first is the identification of past events in the sedimentological record (tsunami deposits) and historical archives in order to build a tsunami catalogue (Maramai et al., 2014; Pistolesi et al., 2020). The second is based on the interpretation of observations and instrumental data (Bonaccorso et al., 2003; Pino et al., 2004; Boldini et al., 2005; Calvari et al., 2005; Maramai et al., 2005a, b; Tinti et al., 2005, 2006a, b; Tommasi et al., 2005; Acocella et al., 2006; Chiocci et al., 2008a, b; Fornaciai et al., 2022). The third comprises model-based scenarios through data-driven simulations (Fornaciai et al., 2019; Bonilauri et al., 2021; Esposti Ongaro et al., 2021; Cerminara et al., 2024).

Finally, a series of initiatives have been carried out to enhance awareness of the tsunami risk and to enforce mitigation actions aimed at the timely evacuation of the risk zones in case of tsunami via the “Io non rischio” awareness campaign (https://www.iononrischio.gov.it/en/get-ready/volcanoes/stromboli/what-do/, last access: 8 January 2024) and via several work tasks as part of Task 4.3, which is an ongoing task of the Italian Civil Protection Department–Istituto Nazionale di Geofisica e Vulcanologia (DPC–INGV) 2022–2024 agreement for the implementation of the service activities: “Survey on risk perception of explosive paroxysms and tsunamis to better define a communication strategy and informative materials”.

In this paper, we describe a new procedure that allows us to combine the outputs of an ensemble of pre-computed tsunami hazard scenarios and exposure data and digital geographic information to obtain a quantitative (scenario-based) risk analysis and to quantify the effectiveness of potential risk mitigation by means of evacuation. In our procedure, the tsunami hazard is represented by a set of (static, georeferenced raster) maps of inundation (maximum wave height above ground) and wave arrival times, one for every individual scenario produced by numerical simulations. Every tsunami scenario is in turn defined, in a deterministic approach, by a corresponding landslide scenario, and it is associated with a trigger time (which approximately corresponds to the initial time at which an alert is issued). For every scenario, we are able to calculate an “impact score”, which classifies in a simple and intuitive way the scenarios in terms of their impact on the island shores. For exposure, we consider the population distribution on the island and the geometry of the road network, with which we can compute escape times from any pixel of the map towards a refuge area entry point (RAEP). By convolving hazard and exposure maps, we are then able to obtain maps expressing the level of risk of the different areas along the Stromboli shores in terms of the potential impact of a given scenario and potential mitigation in terms of evacuation capacity.

The procedure here applied to the island of Stromboli is a general one, and we discuss in the following sections its main components, assumptions, and criticalities.

2 Methods

2.1 Tsunami hazard scenarios

In this paper, we base our analysis on tsunami scenarios produced by numerical simulations with a coupled landslide–water multilayer, non-hydrostatic shallow-water model (Esposti Ongaro et al., 2021; Cerminara et al., 2024). The landslide is considered a granular fluid, having a given initial volume, position, and density, which is dynamically two-way-coupled with the water layers. The Multilayer Hyperbolic Systems and Efficient Algorithms (Multilayer-HySEA) numerical model (Fernández-Nieto et al., 2018; Macías et al., 2021a, b) adopted to generated the scenarios is particularly suited to the case of a volcanic island, since typical landslide-generated tsunamis have short wavelengths and develop over steep topo-bathymetry, making the usual hydrostatic approximation fail.

Simulations using Multilayer-HySEA were performed using 10 different initial landslide positions (positions 1 to 10, Fig. 1), five different volumes (5×10⁶, 8×10⁶, 14×10⁶, 21×10⁶, and 30×10⁶ m³) and three different densities (2.5, 2.0, and 1.7 kg m⁻³ or water–landslide density contrasts of 0.4, 0.5, and 0.6, respectively), based on the ranges hypothesised for the 2002 tsunami event at Stromboli (Fornaciai et al., 2019; Esposti Ongaro et al., 2021; Cerminara et al., 2024). Additional simulations from a higher subaerial position (position 0, Fig. 1) were run with two different volumes for a total of 156 simulations analysed. For each simulation, we consider (1) the coastal inundation (i.e. the estimated inundation depth for each onshore pixel), (2) the tsunami arrival time for each pixel in the inundated area, and (3) the offshore water surface elevation.

Figure 1 Location of the island of Stromboli and landslide position numbers simulated by the INGV of Pisa with their respective coordinates and altitudes and position number 0 corresponding to December 2002 characteristics. The authors made maps using Quantum GIS version 3.16.7 software (2021). Image data © 2021 Google Earth: SIO, NOAA, US Navy, NGA, GEBCO.

For each simulation, we also define an onset time T₀, when the landslide starts; a final time T₆₀₀, when the simulation is stopped; and “trigger time” T_g, when the generated wave overcomes a given threshold at one of the two gauges offshore the Sciara del Fuoco. These are located on the position of the actual beacons, installed close to Punta Labronzo (north-east beacon, BNE) and Punta dei Corvi (south-west beacon, BSW). In this work, we decided to adopt a threshold of +0.3 or −0.3 m for the wave detection at the gauges, to analyse only “significant” scenarios in terms of impact. The choice of such a threshold was made on a subjective basis and might be the object of further analysis in the future.

The trigger time T_g should be approximatively equivalent to the time at which the alert is issued, following the procedure described by Selva et al. (2021) and Lacanna and Ripepe (2020). Finally, the wave arrival time at a given pixel is computed as the difference between the actual tsunami arrival time and the trigger time T_g.

All simulation and post-processed data were integrated in QGIS (Quantum GIS) 3.16.7 with GRASS (Geographic Resources Analysis Support System) version 7.8.5. Topography was represented by a digital elevation model (DEM) of 31 635 pixels (20 m by 20 m) from the lidar (light detection and ranging) campaign carried out in July 2010 by the INGV, and bathymetry was from the Marine Geohazards along the Italian Coasts (MaGIC) project (for more details, see Favalli et al., 2009; Chiocci and Ridente, 2011; and Fornaciai et al., 2019).

2.2 Inundation impact score and link with the proximal wave height

For each landslide and tsunami scenario, we define an impact score S, equal to the number of on-land pixels (in the digital model) that are inundated at a given time during the numerical simulation (we might also consider our impact score in terms of inundated area in m²). The impact score allows us to classify the landslide–tsunami scenarios based on their coastal hazard and to link such a hazard to the features of the tsunami at the proximal gauge. In particular, we have used several statistical methods to try to establish a robust link between the impact score and the maximum wave height at the proximal gauge, which is, in principle, a measurable quantity. The classic method for this type of problem (that is, to find an unknown value from several knowns) is linear regression. This approach, however, lacks robustness when the number of explanatory variables (here, 40) is too large compared with the number of individuals (here, the 156 simulations). LASSO linear regression is thus a modification of traditional linear regression that identifies a subset of explanatory variables (in this case, times of interest for measuring wave heights) of sufficiently small size for the results to be robust.

In a real case, the volume, density, or position of the landslide is not immediately known; thus we only have the signals of the two gauges. As a result, the aim of LASSO is to find out whether there is any chance of detecting the impact of tsunamis on coastlines both quickly and accurately (without counterproductive false alarms) before the tsunami arrives and without knowledge of the landslide characteristics. We thus determine how we can use simulated wave signals to correctly determine the impact score. To set the impact score, we used the signals from the two beacons, i.e. 40 variables (one wave height every 2 s between 0 and 40 s for two beacons). The statistical procedure is described in detail in Appendix A.

2.3 Pedestrian horizontal evacuation model

For every simulated scenario, we applied an evacuation model to each inundated pixel using the approach of Bonilauri et al. (2021). This is a macroscopic model, i.e. a model for global evacuation and not an agent-based approach, which determines the fastest pedestrian evacuation paths from a danger point to a safe point and compares the escape time with the wave arrival time.

To apply the evacuation model, a grid with the same mesh size as the inundation model (20 m) was created for the village of Stromboli (Fig. 2), and a centroid, i.e. an evacuation starting point, was placed centrally in each pixel. A speed reduction coefficient was assigned to each road segment according to its width, slope, and surface type (Table 1). For each tsunami simulation, the quickest escape routes were projected from each inundated centroid, and times required to evacuate each inundated pixel were calculated, with evacuation time being the time needed to move from the inundated pixel to a “safe” zone defined by the limit of the area impacted by the tsunami.

Figure 2 Establishment of baseline data: starting centroids for the numerical simulation of a 30×10⁶ m³ tsunamigenic landslide that detached from position number 6 (294 m b.s.l.) with a density of volcanic materials of 2.5 kg m⁻³ and the Stromboli road network and refuge area entry points used to run evacuation models. The map was produced by the authors using Quantum GIS version 3.16.7 software (2021). Image data © 2014 LGS.

Table 1 Pedestrian evacuation speeds depending on land slope and road types (Péroche, 2016; Bonilauri et al., 2021).

To reach the safe zone, evacuees must pass through a refuge area entry point (RAEP). Two categories of RAEP were established:

1. A normal-event RAEP was assigned for tsunamis with run-ups less than or approximately equal to the run-up of the December 2002 event; i.e. the RAEP was placed at the intersection between the road network and 15 m contour line.

2. An extreme-event RAEP was set for tsunamis greater than the December 2002 run-up and placed at the intersection between the road network and 35–40 m contour line.

For each inundated pixel, the difference between the time needed for escape and the wave travel time to the pixel was calculated. In particular, we used two source-to-pixel travel times. Firstly, we used T₀, i.e. the onset time of the landslide trigger, and then T_g, the time of the first detection of the generated tsunami at one of the two beacons. It is worth highlighting that the use of T₀ gives a longer escape time, but using T_g is more realistic given the current beacon-based alert system.

Figure 3 Spatial distribution of population on the island of Stromboli based on Internet public data (accommodation, restaurant–bar–cafe, and beach–port). The authors produced maps using Quantum GIS version 3.16.7 software (2021). Image data © 2014 LGS.

Finally, a spatial-population-distribution layer (Fig. 3) was created using publicly available data, allowing the creation of three categories of population:

• The first category is populations in accommodation, i.e. those people resident in holiday rentals, hotels, and bed and breakfast (B&B) locations. We generally used websites associated with any given establishment, and other websites such as Booking.com were used if no direct booking website was available. Capacities for all active establishments were thus taken from as individual rental/hotel/B&B websites and, if unavailable, from Booking.com, Airbnb, Tripadvisor, Vrbo, Abritel, https://www.gites.fr/ (last access: 12 October 2023), or https://www.cybevasion.fr/ (last access: 12 October 2023).

• The second category is populations using restaurants, bars, and cafes, with capacities assessed from Tripadvisor and Google Images.

• The third category is pedestrian traffic in and out of the port and visitors to beaches. This was assessed from our own pictures and head counts made during surveys in September 2022 and in June and September 2023.

For locations where information was not available, we averaged the available capacity of the holiday rentals found and used a density of one person per 12.3 m². We then determined occupation scenarios by season and time of day. This distinguished between winter, when only the permanent population (≤150) is present, and the summer tourist season, when as many as 5000 visitors can be on the island. It also distinguished between morning, midday, afternoon, evening, and night, between which the distribution of workers and tourists across accommodation, restaurants–bars–cafes, and beaches–the port varies.

Combining the vulnerability determined from these surveys with the hazard output for each of the 156 tsunami simulations provided the number of pixels inundated and, for each inundated pixel,

• wave arrival time (T_arrive), where T_arrive is the moment when the generated tsunami has reached a pixel and inundation begins to be detected (1 cm threshold). T_arrive being calculated from T_g;

• inundation depth (i.e. water thickness above ground level);

• evacuation time (T_evac), where T_evac is the moment when evacuees reach a refuge area entry points (RAEPs);

• “evacuability” (T_arrive−T_evac);

• the number of people in need of evacuation.

We statistically convolved these metrics to allow a full and robust scenario-based risk assessment, which includes generation of the probability of inundation and an impact score.

Figure 4 Example of a numerical simulation for a 30×10⁶ m³ tsunamigenic landslide that detached from position number 6 (294 m b.s.l.) with a density of volcanic materials of 2.5 kg m⁻³. (a) Sea level variation around the island of Stromboli due to the propagation of the tsunami wave. (b) Time of travel of the tsunami wave around Stromboli. The authors produced maps using Quantum GIS version 3.16.7 software (2021). Image data © 2014 LGS.

3 Results

3.1 Analysis of individual scenarios: arrival time and coastal impact

The simulation outputs range from tsunamis of just 0.22 m in amplitude, with run-ups of 0.29 m and impacting only some parts of the beaches, to tsunamis with amplitudes of 48.1 m and run-ups of 24.2 m that inundate almost the entire village (Cerminara et al., 2024). For means of demonstration, we take a mid-range example that falls between these two extremes, this being scenario P6V30CD0.4. This scenario was selected as being slightly larger than the event of 2002 but within the same order of magnitude, as we have already conducted a detailed analysis of the 2002 events (Bonilauri et al., 2021). Scenario P6V30CD0.4 involves the simulation of a submarine landslide involving 30×10⁶ m³ of volcanic material with a density of 2.5 kg m⁻³ from position number 6 (centred at 294 m b.s.l.). Fornaciai et al. (2019) compared the model output with the extent of inundation and run-up recorded following the 2002 tsunamis and showed that the model favours subaqueous landslide volumes of 15×10⁶–20×10⁶ m³ and/or a subaerial landslide of 4×10⁶–6×10⁶ m³ on the Sciara del Fuoco. The P6V30CD0.4 scenario used here as the example thus produces a slightly larger tsunami, with a 12.5 m maximum run-up at Spiaggia Lunga, whereas a 10.9 m run-up for the 2002 events was measured by Tinti et al. (2006a). Such a scenario is likely to occur and could be mitigated for, unlike doomsday end-member scenarios.

The generated tsunami inundates an area of 544 pixels, i.e. a surface of approximately 217 600 m², with the wave trapped around the entire island (Fig. 4a). Refraction of the wave causes some variation in tsunami travel times over short distances, where arrival times (computed using the trigger time T₀ as an initial time) vary from 48 s at Spiaggia Lunga (north of Piscità) to 154 s at Ficogrande and to 188 and 242 s at Punta Lena and the port at Scari, respectively (Fig. 4b). These four vulnerable sites are about 2.4, 3.5, 4.1, and 4.9 km distant, respectively, from the Sciara del Fuoco. The arrival times computed using the landslide onset time T₀ as an initial time are 20 s longer.

Figure 5 Example of a numerical simulation for a 30×10⁶ m³ tsunamigenic landslide that detached from position 6 (294 m b.s.l.) with a density of volcanic materials of 2.5 kg m⁻³. (a) Extent of inundation zone and flow depths. (b) Arrival time of the landslide-induced tsunami on Stromboli coasts. The authors produced the maps using Quantum GIS version 3.16.7 software (2021). Image data © 2014 LGS.

The coastal inundation simulation for scenario P6V30CD0.4, i.e. the maximum height above ground of water reached by the tsunami in each inundated pixel, is represented in Fig. 5a. Inundation depths are highly variable due to the topography of the village, which is underlain by around 11 lava flow fields that erupted from eccentric vents just above the village between 15 and 2 ka (Calvari et al., 2011; Speranza et al., 2008). To the north of the village, while the Spiaggia Lunga is well-exposed, the 20–30 m high sea cliffs behind it protect the road and buildings from inundation. The same is true for Piscità, the first district of Stromboli to be reached by the tsunami, which is protected due to its location on the San Bartolo lava flow field. The same presence of the lava flow field along a coastal length of 1 km means that 10 m high cliffs protect the population as far as Ficogrande. In contrast, Ficogrande is a bay located at the SE edge of the San Bartolo lava flow field, which focuses the tsunami wave and produces an increase in run-up behind the bay. To the SE of Ficogrande, another 0.2 km stretch of coast is protected by the 15 m high sea cliffs associated with the San Vincenzo lava flow field. However, beyond this in the southern part of Punta Lena and the northern part of Scari, the relatively flat coastal topography for an area extending up to 150–200 m inland results in an increase in the extent of inundation. Refraction of the tsunami wave around the island produces a decrease in the wave height to the south, causing low levels of inundation south of Scari.

The peculiar behaviour of wave propagation around the island impacts the arrival times of tsunamis too. This is particularly true on the low promontory of Punta Lena. Here we observe a narrow, in-land, extending band of pixels on the north side of the promontory with an arrival time of 300–600 s, while pixels immediately to the south have an arrival time of between 120 and 300 s (Fig. 5b).

3.2 Impact score

Figure 6 shows the impact score function of landslide position and volume. After no change in the score between positions 0 and 3, the impact score systematically decreases as a function of landslide position between locations 4 and 10 (Fig. 6a). The level of impact increases with volume but with each curve having a similar shape. Subaerial positions (i.e. position numbers 0–3, Fig. 1) have higher impact levels but with little dependence on altitude. In the case of submarine positions (i.e. position numbers 4–10, Fig. 1), the impact decreases with depth (below sea level) of the landslide source. For all positions there is a positive, and broadly linear, relation impact and volume (Fig. 6b). In the case of the subaerial positions, position 3, which is closest to sea level (Fig. 1), shows a slightly higher impact score than the higher-elevation subaerial positions 0, 1, and 2. This effect is related to the fact that the sliding volume of volcanic material has not had time to deform before reaching the sea surface, so its entry is more focused.

Figure 6 Evolution of the impact score as a function of (a) the position for each volume of volcanic material landslide in the Sciara del Fuoco and (b) the volume for each landslide position (position 0 to position 3: subaerial landslide; position 4 to position 10: submarine landslide). The authors produced graphs using MATLAB version 2019b.

3.3 Analysis of the signals at virtual gauges and link with the impact score

The purpose of LASSO is to find the impact score based on the beacon data (without having to analyse the shape of the waves or the volume of the slip from the signal). Our LASSO method is robust in terms of its ability to adapt to landslides with volumes of between 5×10⁶ and 30×10⁶ m³, densities of between 1.7 and 2.5 kg m⁻³, and landslide source positions of between 500 and −584 m. That is, it is robust within the limits of our simulations. For our simulated beacon signals, we work with landslide models whose smallest volume is 5×10⁶ m³. We are thus not considering the same volume scale as scales that are associated with tsunamis generated by pyroclastic density currents (PDCs) at Stromboli, whose volumes were an order of magnitude smaller than our smallest-volume simulation: the 2019 volumes were estimated by Ripepe and Lacanna (2024) as 2.08×10⁵ m³ (July) and 1.05×10⁵ m³ (August), and the May 2021 event volume was estimated at around 8.4×10⁵ m³ by Calvari et al. (2022). Even so, Ripepe and Lacanna (2024) showed that the shape of the waveform did not change with slip volume and that it was possible to determine the inundation extent/run-up from the waveform amplitude. Ripepe and Lacanna (2024) also showed that real waveforms matched model-derived (NHWAVE) waveforms, and NHWAVE has been shown through benchmarking to produce results comparable to the model used here (Esposti Ongaro et al., 2021). We thus have confidence that out synthetic waveforms are valid.

In Fig. 7, we link the impact score for each scenario to the absolute value of the wave height. It is worth noting that in the evacuation analysis we used an initial time defined by T_g, i.e. the moment when a signal with an amplitude higher than a threshold of 0.3 m is detected. Please note that in the cases in which the signal never reaches this threshold, we assume that the landslide will not generate a tsunami large enough to require an evacuation and will therefore not be considered. Figure 8 shows the wave heights registered at the NE gauge (Fig. 7a) and the SW gauge (Fig. 7b) in the first 40 s (starting from T₀). The first peak is in most cases detected between 20 and 28 s at the NE beacon and between 8 and 16 s at the SW beacon. Not surprisingly, we observe that the highest impact scores are always associated with the highest wave heights.

Figure 7 Gauge signals with the (a) north-east gauge and (b) south-west gauge. The colour code corresponds to the impact score (in m² inundated). The bold red lines correspond to the 11 explanatory variables retained by the LASSO linear regression model.

Figure 8 Inundation expected in m² by tsunami gauge signals according to (a) a simple regression approach and (b) a LASSO linear regression approach. Blues points represent each simulation of tsunamigenic landslide, and the orange ellipse represents tsunamis with a high impact score that are underestimated by the approach. The authors produced graphs using MATLAB version 2019b.

To refine this relationship between wave amplitude and impact score, two approaches were adopted (see Methods and Appendix A). The first “rough” approach is a simple linear regression; i.e. we take the maximum amplitude value and relate it to the impact score (Fig. 8a). With this approach, some tsunamis with high impact scores are underestimated because the score predicted by the gauge is lower than the score determined by the inundation models (see orange ellipse in Fig. 8a). In the second method, the LASSO penalised linear regression algorithm automatically chooses which gauge data values are most likely to explain the impact scores. The method finds 11 “explanatory variables” that best explain the relationship between gauge data and impact scores (Fig. 7). As detailed mathematically in Appendix A, these variables are the 11 most critical times at which a point on the waveform better reproduces the inundation area. The 11 explanatory variables selected were times of 6, 8, 12, 24, 28, 38, and 40 s for the NE beacon (Punta Labronzo) and 6, 26, 34, and 40 s for the SW beacon (Punta dei Corvi). This is then used to best define and distinguish the specific shape of each given waveform and links it to its inundation capacity. The LASSO regression shows a closer relationship between the gauge data and impact scores than the simple regression and does not underestimate any tsunami with a high impact score (Fig. 8b). The LASSO regression also results in small errors, defining a linear relation with a low degree of scatter.

3.4 Evacuation capacity

We here define “warning time needed” and “real warning time”. The former is the time needed to move from any given point to a safe point, thus corresponding to the exit from the inundation zone refuge area entry point (RAEP), while the real warning time is the time available for escape prior to wave arrival. We thus have two types of point, (1) those for which the RAEP can be reached inside the threshold time and (2) those that cannot be reached inside the threshold time (i.e. outside of the threshold time). The threshold time is, for any given scenario, the wave arrival time minus the time needed to reach the RAEP. If this is negative, the point is outside of the threshold time.

Evacuation modelling gives the fastest routes from each inundated pixel to an RAEP, i.e. an entry point into a safe zone. Results for Piscità, the most proximal community to the source, are given in Fig. 9 by way of example. Results for the Ficogrande, Punta Lena, and Scari are given in Appendix B, Figs. B1–B3. As evacuation speeds consider the slope and nature of the escape route (path width, surface type, presence of steps, etc.), some routes are not used because, although they appear short, they are slow. This explains why RAEP4 and RAEP9 in Piscità (Fig. 9), as well as RAEP12 in Ficogrande (Fig. B1), RAEP15 in Punta Lena (Fig. B2), and RAEP22 in Scari (Fig. B3), are not chosen as viable refuge area entry points.

Figure 9 Example of an evacuation model at Piscità created for a numerical simulation with a tsunamigenic landslide of 30×10⁶ m³ having broken away from position 6 (294 m b.s.l.) with a density of volcanic materials of 2.5 kg m⁻³. A total of 64 pixels will be inundated (25 600 m²), and their associated flow depth and the fastest evacuation routes are shown. In green, points are inside the threshold time, and in red, points are outside of the threshold time. The authors produced maps using Quantum GIS version 3.16.7 software (2021). Image data © 2014 LGS.

For each pixel in the inundated zones, we used the time for the wave to arrive (T_arrive) minus the time needed to travel between the pixel centroid and the closest RAEP, in terms of time (evacuation time, T_evac), to set the threshold time A positive difference for T_arrive−T_evac thus means that the pixel can reach an RAEP inside the threshold time for the given scenario, these being the green centroids (shown as circles) in Fig. 9 and in Appendix B, Figs. B1–B3. In contrast, a negative T_arrive−T_evac value means that the pixel reaches an RAEP outside of the threshold time, these being the red centroids in Figs. 9 and B1–B3.

Figure 10 Minimum and maximum real warning time given by the tsunami gauges and warning time needed for each district to evacuate from the coastal inundated area.

In our selected simulation example (scenario P6V30CD0.4), out of 544 inundated pixels (217 600 m²), only 132 (24 %) pixels (52 800 m²) are theoretically inside the threshold time (i.e. the RAEP can be reached before the tsunami arrives). To understand if these pixels can reach a safe RAEP point before the arrival of the tsunami, we compared the real warning times (minimum and maximum for any given zone) available from the tsunami detection gauges and the warning times needed (minimum and maximum) to evacuate the inundated area (Fig. 10). In terms of number of inundated pixels, we observe that for Piscità, Punta Lena, and Scari there is time to evacuate some pixels but not the whole area (warning time needed > real warning time). To evacuate the whole area, a maximum of 102, 174, and 368 s of extra warning time is needed for Piscità, Punta Lena, and Scari, respectively. For Ficogrande, the two curves “warning time needed” and “real warning time” are very close (Fig. 10), which means that most pixels can reach an RAEP inside the threshold time (i.e. before the wave arrives). At most an extra 41 s is needed to evacuate the whole area.

Figure 11 Comparison between the real warning time given by the tsunami gauges and the warning time needed to evacuate 25 %, 50 %, 75 %, and 100 % of the Stromboli population from the four main inundated districts. The difference between the filled and empty bars is the number of cases for which we lack the time to carry out evacuation at a point for each evacuation scenario.

Figure 12 Table and bar chart of the global number of people who can reach an RAEP inside or outside the threshold time (i.e. before the arrival of the tsunami) in the case of a 30×10⁶ m³ tsunamigenic landslide at position 6 (294 m b.s.l.) with a density of volcanic materials of 2.5 kg m⁻³.

In Fig. 11 we assess the time needed to evacuate 25 %, 50 %, 75 %, and 100 % of the population from the inundated zone and thus the warning time needed to complete partial or full evacuation. In doing this, we need to distinguish between an evacuation in winter and one in summer. In summer when the coastline is highly populated, differences in the time needed to evacuate 25 %, 50 %, 75 %, and 100 % of the population can be observed depending on the time of day (Fig. 11).

We see in Fig. 11 that the time of day during which evacuation is most effective varies from district to district. For example, for Ficogrande the best time is daytime when tourists gather in locations with easy access to RAEP13. In contrast, for Scari and Punta Lena all times of day are equally bad due to the long distances to the RAEPs. For Piscità, the situation is slightly better during the night, when people leave the highly exposed beach and gather in the village closer to RAEPs. The situation is exacerbated in Piscità due to its proximal location and thus the very short wave arrival time.

To summarise the evacuation capacity of the zone inundated by scenario P6V30CD0.4, we created bar charts giving the distribution between the number of agents that can reach an RAEP inside and outside of the threshold time depending on the season and time of day (Fig. 12). For this scenario, it is extremely difficult to reach a safe point on time. Out of 544 inundated pixels (217 600 m²), in the best case (a tsunami occurring during a summer night), only 30 % of the population can reach an RAEP inside the threshold time (Fig. 12). We also prepared individual bar charts for each of the four main neighbourhoods of the village of Stromboli (see Appendix C, Figs. C1–C4). Depending on the district, we can see that during certain periods of the year and certain times of the day, some areas are not occupied (i.e. their bar charts have empty columns). This is particularly the case during the winter period, when a “winter desert” can be found in the districts of Piscità and Ficogrande, which are the locations of summer rentals, B&Bs, and hotels. Unfortunately, permanent residents are concentrated in the districts of Punta Lena and Scari, with Punta Lena being a zone with mainly (red) points that reach an RAEP outside of the threshold time (Fig. B2). This explains why for a winter night tsunami, when the Punta Lena population are at home and close to the shore, the evacuation capacity is 0 % (Fig. C3).

4 Discussion

4.1 Landslide source parameters

The choice of the source input parameters was driven by the need of the Italian Civil Protection Department to analyse tsunami impacts in the range of the 2002 event (Tinti et al., 2006b). Larger events are possible (with the upper bound of the 1.81 km³ mass movement for the 5 ka event that formed the Sciara del Fuoco; Kokelaar and Romagnoli, 1995), and they are documented in the geological history of the island (Rosi et al., 2019; Pistolesi et al., 2020).

For what concerns the landslide initial positions, the highest impact is always associated with landslides triggered at the lowest subaerial positions (position 3, which is just above sea level at 58 m a.s.l). This is due to the fact that granular flows starting at higher elevations deform and are diluted during the downward motion, reducing the front thickness at the impact with water and thus producing smaller waves. Finally, analysis of the whole database allows us to state that the bulk density of the landslide has a second-order effect on its tsunamigenic capability. Our analysis, carried out on simulation 6 with a density contrast of 0.4, is substantially unchanged for lower densities.

In Appendix D, Fig. D1, we show the maximum thickness of a 5×10⁶ m³ subaerial landslide (measured throughout the entire simulation) for three different initial positions along the Sciara del Fuoco (positions 0, 2, 3), thus demonstrating that the maximum landslide thickness at the impact with water is higher for the lowest-elevation landslide. Conversely, the highest initial position of landslide results in significant lateral spreading and a reduced thickness of the landslide.

We notice in our simulations that, even if the first peak is higher for low-elevation landslides at the proximal gauges, the later peaks are generally of lower amplitude (Fig. E1). This might support the idea that the tsunami energy is more focused on the first peak when the subaerial landslide is less dispersed in terms of lateral extent, but the total energy might still be correlated with the initial landslide potential energy. Verifying such a correlation for granular landslides is thus another of the objectives of our future studies.

In this study, landslide position and volume are considered the key parameters in determining the hazard score (Cerminara et al., 2024). The dependency on other parameters such as the landslide initial aspect ratio and horizontal position across the Sciara del Fuoco are of second order. Numerical simulations used in this work assume an initial position of the landslide along the central axis of the Sciara del Fuoco (consistently with the 2002 event). This assumption provides an average travel time for the tsunami to the village of Stromboli (Fig. 1). Off-axis landslides will give longer or shorter travel times (by ±30 s) for the direction of the village of Stromboli depending on whether the landslide is in the western (furthest from the village) or eastern (closest to the village) sector of the Sciara del Fuoco, thus affecting the results of the evacuation model. However, this measure can be very sensitive to the wave height, especially in areas characterised by small slopes or wall-bounded regions (e.g. topographic pools). Future numerical investigations will sample initial landslide positions considering potential off-axis sources. Finally, it is worth remarking that different assumptions about the granular flow rheology and water–landslide friction might significantly affect the results. Such effects will be analysed in a dedicated future work.

4.2 Human exposure

For the Aeolian Islands, demographic data are available from the official public censuses of 2021 (yearly survey by ISTAT, https://www.istat.it/, last access: 17 October 2023), Istat.it, 2023 and 2011 (https://www.citypopulation.de/, last access: 10 October 2023, Citypopulation.de, 2023). However, data are generally grouped by municipality for ISTAT and by island and sometimes by district for City Population. There are no data for the number of people living in each building or their demographics. Thus, we had to rely on publicly available data, such as hotel booking forms on the Internet, as well as a ground-based census completed by us for permanently occupied houses in January 2020 (Bonilauri et al., 2021). Some areas have therefore probably been overestimated in terms of population capacity, under the assumption of full capacity for all holiday facilities during the summer season. Ideally, precise population distribution and demographic data would allow us to adjust our evacuation model, but we consider the data used useful in providing a worst-case (full-capacity) scenario.

For our evacuation time calculations, we considered a “standard” person with an associated forced walking speed (see Table 1; Péroche, 2016; Bonilauri et al., 2021). These calculations were calibrated and checked as valid during escape simulations carried out by us in January 2020. Obtaining data on the age and physical abilities of the Stromboli inhabitants would allow construction of evacuation plans tailored to the individuals' capacities (see Bonilauri et al., 2021). Our escape times also involve no reaction time: the reaction is immediate, spontaneous, and correct, meaning that the agent escapes immediately using the fastest route.

Provitolo et al. (2015) presented two main kinds of behaviours in psychological reactions to imminent danger. The first was “instinctive” behaviours, where panic or disbelief prevail, and the second was “learned” behaviours, where reactions to danger are taken in a reflective manner, are adapted to the context, and are no longer instinctive. These learned behaviours can be found in preventive or spontaneous evacuations when the context is already known by the population, has been anticipated, and is likely the case for the resident populations (Bonilauri et al., 2021). However, visitors are likely to fall into the instinctive category. This poses a problem for risk management and mitigation. That is why risk awareness campaigns have been organised in volcanic contexts by the Italian Civil Protection Department for several years to reduce reaction times in the face of real danger among the visitor population.

Reaction time is hard to consider because it is specific to each person, depending, for example, on knowledge of, experience of, and familiarity with the hazard. Establishing the reaction time of both informed and uninformed agents present in vulnerable areas when the sirens sound is thus our next focus and will allow us to factor this uncertainty into our evacuation models. The same tests will allow us to assess, and refine, the escape times.

In winter, tourism on the island of Stromboli is greatly reduced, with visitors numbering below 5 per day and rising to 20 for two nights when a visiting football team was on the island (based on our January 2020 survey). Thus, generally, only permanent residents, who numbered 98 during our January 2020 count, are present in the winter months. Because most of the population live along the upper roads in the village, under winter conditions only around 20 residents are in vulnerable zones (mostly Punta Lena). However, these zones are at high risk, with only an 11 %–13 % escape percentage on a winter evening or night (Fig. 12), and thus require special attention in terms of evacuation needs.

The real warning time available to us (given by beacons) to evacuate as many people as possible before the tsunami reaches the coast is extremely short due to the short (<5 km) distance between the source and point of impact. This means that, even for the best-case scenario, only 30 % of the vulnerable population can reach a safe point inside the threshold time (Fig. 12). As the warning time cannot be reduced, only accelerating the evacuation can guarantee the safety of the people living on or visiting the island. Moreover, the greater the volume of the landslide is, the higher and faster the tsunami generated will be, reducing the warning time. Thus, this highlights the need to link tsunami waveforms to impact scores and tailor evacuation times accordingly.

To speed up horizontal evacuation, a number of RAEPs and ease of access to them are necessary, as are the distribution of evacuation plans and installation of signage. However, in cases where horizontal evacuation is clearly identified as impossible, vertical evacuation capable of taking the required capacity is the only solution (Bonilauri et al., 2021; Turchi et al., 2022). Since the response needs to be evacuation, either vertical or horizontal (Leone et al., 2013, 2018; Péroche, 2016; Solís and Gazmuri, 2017), the evacuation plans need to be developed well in advance. Our approach can be used to identify zones suitable for horizontal and vertical evacuation and be used to set capacity.

4.3 Anticipated inundation

The use of the LASSO penalised method (see Method, Fig. 8, Appendix A) to link a beacon signal to a predicted inundation extent is a major advance in our tsunami risk assessment and management. Testing this method on real tsunami detection beacon signals instead of simulated signals would make it possible to determine whether the method is valid and reliable in real conditions, and eventually integrating it directly into the beacon algorithms would make it possible to estimate the impact zone before the tsunami reaches the inhabited coast. Currently LASSO needs 40 s to recognise and classify the waveform. This time is the result of a payoff between precision in inundation area and the time needed to classify the waveform to an acceptable degree of accuracy. We could reduce this time, but that would make model output decreasingly accurate. This means that for proximal locations, maps will be delivered after the wave arrival, but the locations will have been alerted by the siren.

Having an estimation of how much coastline will be inundated and therefore needs to be evacuated could also aid in disaster planning and management, as well as guiding search-and-rescue follow-up. The refuge area entry points could be automatically adapted to each tsunami case by lowering them in elevation for small tsunamis and raising them in elevation for larger tsunamis. In our current study, we considered two types of RAEPs: 15 m a.s.l. for smaller tsunamis and 35–40 m a.s.l. for larger ones. This would modulate the population in need of evacuation and regulate the load placed on escape routes, which could be blocked for agents at risk lower in the zone by agents in the upper zone undergoing an unnecessary evacuation. Such an approach allows output of the hazard maps and risk assessments in real time, tailored to the case in hand. We are currently testing such problems using on-site escape tests with multiple agents faced with differing traffic flow scenarios, starting positions (beach, water, bed), and degrees of preparation (footwear, reaction time). This will enable us to better calibrate our escape times and to produce a distribution of potential times depending on traffic and route conditions.

5 Conclusion

At Stromboli, except the two tsunamis of December 2002, landslide-related tsunamis have few or no data available. In such cases, hazard scenarios must be supported by modelling that is validated by the well-constrained cases. The same can be argued for vulnerability and the modelling of mitigation efforts, such as evacuation, where gaps need to be identified and filled through extrapolation. Our integrated modelling approach, which involves a statistical combination of model-output hazard metrics (wave run-up, arrival time, and seafloor pressure sensor waveform) with seasonally and diurnally distributed populations as well as a model for evacuation times, allows risk mapping to assess whether or not horizontal evacuation is possible at a given point.

The simulated scenarios allow us to consider a wide range of source conditions, enabling a statistical approach to hazard assessment and providing scenario-tailored evacuation maps. Our approach involves 156 hazard (tsunami) scenarios and 5 diurnal (morning, midday, afternoon, evening, and night) population distributions divided by season (winter and summer). This results in $\mathrm{156}\times \mathrm{5}\times \mathrm{2}$ combinations of risk assessment scenarios. Each risk assessment scenario can be linked to the tsunamic type through the waveform recorded by a seafloor pressure gauge and is automatically selected accordingly. Given a warning time of less than 3–6 min on a small volcanic island where the tsunami source is just a few kilometres from the vulnerable shoreline population, such a pre-prepared and automatic event outcome procedure is indispensable.

This method is tested and applied here to Stromboli, but it can serve as a blueprint for any other volcanic islands where volcanic activity and landslides cause local tsunamis with little warning time. Key elements in this regard are to collect the two most important source terms for hazard modelling (landslide location and volume) as well as spatial and temporal assessments of population distributions. Given uncertainty in the source terms and population distributions our scenario-based approach considers all simulations to have the same probability of future occurrence, as well as a worst-case scenario for population distribution. This approach is thus consistent with the precautionary principle of disaster management while being open to refinement subject to further data collection from, for example, geological studies, expert elicitation, escape tests, and/or population surveys, which will be our next step.

Appendix A

Here we present our statistical analysis of tsunamigenic landslides and penalised linear regression LASSO.

Our aim is to define a tsunami impact score based on the coastal inundation and to understand the link between this score, the associated landslide characteristics, and the gauge signals.

For each tsunamigenic landslide i (i.e. individual i) we have the following information:

• the landslide characteristics given by (V_i, P_i, D_i), where V_i is the volume of volcanic material, P_i is the onset position of the landslide along the Sciara del Fuoco, and D_i is the density of the landslide material;

• the area of coastal inundation described by a vector ${\mathit{X}}_i=(x_{(i,\mathrm{1})}$, …, x_(i,j), …, x_(i,31 635)) of 31 635 values (i.e. the number of pixels of the DEM) of 0 or 1, where $x_{i,j}=\mathrm{1}$ if the jth pixel is inundated by the tsunami generated by the landslide i;

• the associated gauge signals BNE_i and BSW_i, which are two vectors describing the wave heights recorded every 2 s since the beginning of the landslide.

For each scenario, the inundation is given by a map of high-dimensional data (the dimension corresponds to 31 635 values). To reduce the data dimension, we tested three classical statistical methods (PCA, multidimensional scaling, and isomap), and we found that the inundation data are very close to a one-dimensional space characterised by the total number of inundated pixels.

Since we have data for 156 simulations, we used a regression model to reduce the number of scenarios by searching for a relationship between the explanatory variables and the variables to be explained from Y=f(X₁, …, X_k)+ε; we usually look for the “best” function f within a family of function parameters (e.g. for polynomials of degree less than d, and then d is the parameter). Due to the high dimensionality (relative to the number of experiments), we have evaluated two approaches: (1) a rough approach, i.e. retaining only the maximum amplitude of the first major wave detected by each gauge, and (2) a more classical approach, which consists in applying a LASSO penalised linear regression to select a smaller number of explanatory variables (Giraud, 2021).

The purpose of a multiple regression is to explain a variable y belonging to R using p explanatory variables x₁, …, x_p also in R. In our case, y is the hazard score and x₁, …, x_p denotes the signals emitted by the two tsunami detection beacons.

1. Linear regression. It is assumed that y can be determined as a function of x₁, …, x_p in a linear way, i.e. that there is an affine function f such as

   $$\begin{array}{}\text{(A1)}& y=f\left(x_{\mathrm{1}},\mathrm{\dots },x_p\right)=a_{\mathrm{1}}{x}_{\mathrm{1}}+\mathrm{\dots }+a_p{x}_p+b+\mathit{\epsilon }.\end{array}$$

   Here, a₁, …, a_p and b are the unknown coefficients that we will try to estimate from measurements made on individuals. In our case, an individual is a tsunamigenic landslide scenario. The random part ε synthesises the measurement noise and the experimental uncertainties.

   For each individual i we measure Y_i; the value of the variable y; and X_1,i, …, X_p,i, indicating the values of the variables x₁ to x_p.

   To carry out a linear regression, we look for coefficients ${\widehat{a}}_{\mathrm{1}}$, …, ${\widehat{a}}_p$ and $\widehat{b}$ which best allow us to estimate ${\widehat{a}}_{\mathrm{1}}$, …, ${\widehat{a}}_p$ and $\widehat{b}$ from the results of our N individuals.

   For this purpose, we define the following for each individual i:

   $$\begin{array}{}\text{(A2)}& {\widehat{Y}}_i\left({\mathit{\alpha }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\alpha }}_p,\mathit{\beta }\right)={\mathit{\alpha }}_{\mathrm{1}}{X}_{\mathrm{1},i}+\mathrm{\dots }+{\mathit{\alpha }}_p{X}_{p,i}+\mathit{\beta }.\end{array}$$

   If α₁=a₁, α_p=a_p, β=b, and ε=0, then $Y_i={\widehat{Y}}_i$ (α₁, …, α_p, β) and $|{\widehat{Y}}_i-Y_i|=\mathrm{0}$.

   We can show that there is equivalence between “the α_i values are all close to a_i and β is close to b” and “all ${\widehat{Y}}_i$ values are close to Y_i”.

   This is why the coefficients ${\widehat{a}}_{\mathrm{1}}$, …, ${\widehat{a}}_p$ and $\widehat{b}$ calculated in a linear regression are those that make the error criterion, $\sum \left({\widehat{Y}}_i\right({\mathit{\alpha }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\alpha }}_p,\mathit{\beta })-Y_i{)}^{\mathrm{2}}$, minimal.

   Linear regression does not fit the following two situations:

   • when N (total number of experiments done) ≤ p (total number of explanatory variables), i.e. when there are fewer experiments than the number of explanatory variables;

   • when explanatory variables are too correlated.

   Our problem lies in both cases, so we cannot use classical linear regression.

2. LASSO penalised linear regression. There are ways to get around the limitations of linear regression. One of them is the LASSO penalty. LASSO is described and studied in the book Introduction to High-Dimensional Statistics by Giraud (2021), in which it has been proven that this penalty works and is adequate to solve problems similar to ours.

   The LASSO penalty consists in adding to the error criterion of the linear regression, a penalty proportional to the sum of the absolute values of the coefficients. This has the effect of cancelling out many coefficients (which means working on fewer variables and therefore limits the negative effects of the above-mentioned situations).

   The criterion to be minimised becomes

   $$\begin{array}{}\text{(A3)}& \sum {\left({\widehat{Y}}_i\left({\mathit{\alpha }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\alpha }}_p,\mathit{\beta }\right)-Y_i\right)}^{\mathrm{2}}+\mathit{\lambda }\left(\left|\widehat{b}\right|+\sum \left|\widehat{a}i\right|\right).\end{array}$$

   For each λ value we have an associated LASSO regression. To choose a λ value, one can use classical “model choice” methods. For our part, we used the “leave-one-out” method thanks to the parameter selection program present in the MATLAB software.

   The leave-one-out method applied to our case involves, for a λ value and each subscript i (one of our tsunamigenic landslide scenarios), the following steps:

   • Virtually remove individual i values from our experiments (we pretend we have no information for individual i values).

   • Apply the LASSO method to calculate the coefficients (${\widehat{a}}_{\mathrm{1}}$, …, ${\widehat{a}}_p$ and $\widehat{b}$) for all other individuals.

   • Calculate ${\widehat{Y}}_i^{*}$, the estimated value of Y_i from the previously estimated coefficients, and then $E_i\left(\mathit{\lambda }\right)=({\widehat{Y}}_i^{*}-Y_i{)}^{\mathrm{2}}$.

   The error associated with lambda is defined by $E\left(\mathit{\lambda }\right)=\sum E_i\left(\mathit{\lambda }\right)$.

   Then we start again with several values of λ and choose the λ that makes E(λ) minimal.

Appendix B

Here we present examples of the evacuation model at Ficogrande (Fig. B1), Punta Lena (Fig. B2), and Scari (Fig. B3) created for a numerical simulation with a tsunamigenic landslide of 30×10⁶ m³ released from position 6 (294 m b.s.l.) with a density of volcanic materials of 2.5 kg m⁻³. Here we show the pixels that will be inundated and their associated flow depth, as well as the fastest evacuation routes. In green (circles) are points from which an RAEP safe point can be reached inside the threshold time (i.e. before the tsunami arrives), and in red are points from which the RAEP is reached outside of the threshold time. The authors produced maps using Quantum GIS version 3.16.7 software (2021) (image data © 2014 LGS).

Figure B1 Ficogrande – 79 inundated pixels (31 600 m²). Image data © 2014 LGS.

Figure B2 Punta Lena – 178 inundated pixels (71 200 m²). Image data © 2014 LGS.

Figure B3 Scari – 223 inundated pixels (89 200 m²). Image data © 2014 LGS.

Appendix C

Tables and bar charts showing the number of people who can reach an RAEP point inside or outside of the threshold time (before the arrival of the tsunami) in the case of a 30×10⁶ m³ tsunamigenic landslide (initiated from position 6 at 294 m b.s.l. and with a density of 2.5 kg m⁻³) for four different locations (Piscità – Fig. C1; Ficogrande – Fig. C2; Punta Lena – Fig. C3; Scari – Fig. C4).

Figure C1 Piscità – 64 inundated pixels (25 600 m²).

Figure C2 Ficogrande – 79 inundated pixels (31 600 m²).

Figure C3 Punta Lena – 178 inundated pixels (71 200 m²).

Figure C4 Scari – 223 inundated pixels (89 200 m²).

Appendix D

Figure D1 Maximum thickness of a 5×10⁶ m³ subaerial landslide from positions 0 (top panel), 2 (middle panel), and 3 (bottom panel) (Cerminara et al., 2024).

Appendix E

Figure E1 Waveforms at the proximal gauge (Punta dei Corvi) for a 5×10⁶ m³ subaerial landslide initiated from positions 0 (purple), 2 (green), and 3 (blue).