Tsunamis From Submarine Collapses Along the Eastern Slope of the Gela Basin (Strait of Sicily)

Geophysical surveys in the eastern slope of the Gela Basin (Strait of Sicily, central Mediterranean) contributed to the identification of several episodes of sediment mass transport, recorded by scars and deposits of various dimensions within the Pleistocene succession. In addition to a huge failure called Gela Slide with volume exceeding 600 km3, the most studied events show volumes estimated between 0.5 and 1.5 km3, which is common to many other submarine landslide deposits in this region and that can therefore be considered as a characteristic value. In this work, the tsunamigenic potential of two of such landslides, the so-called Northern Twin Slide and South Gela Basin Slide located about 50 km apart along the eastern slope of the Gela Basin, are investigated using numerical codes that describe the onset and motion of the slide, as well as the ensuing tsunami generation and propagation. The results provide the wave height of these tsunami events on the coast of southern Sicily and Malta and can be taken as representative of the tsunamigenic potential of typical landslides occurring along the slope of the Gela Basin.


INTRODUCTION
Continental margins are one of the most favorable environments for the generation of relevant landslide-tsunamis (Masson et al., 2006;Tappin, 2010;Kawamura et al., 2014), due to many factors. Amongst these, one very relevant is the continuous supply of unconsolidated sediments from rivers, which may be activated in a submarine landslide by both seismic shaking and gravitational load. When the collapse starts from relatively shallow water the tsunami generation is particularly efficient: the perturbation is more easily transmitted to the whole water column, and the sliding mass soon attains high velocities, due to the steep slope typical of such environments, that can exceed 10°.
The most adopted approach in the description of mass transport deposits (MTD) along submarine slopes considers the size distribution and their frequency, providing in this way an assessment of the potential hazard connected to these occurrences. This has been repeatedly applied for the hazard analysis in several margins around the world, mainly in North America: the United States Atlantic margin (Chaytor et al., 2009;ten Brink et al., 2014); the United States Pacific coast (McAdoo et al., 2000;Greene et al., 2006); the Puerto Rico northern platform (ten Brink et al., 2006); the Gulf of Mexico (Pampell-Manis et al., 2016;Fan et al., 2020); Alaska (Sawyer et al., 2017). The Norway margin (North-East Atlantic Ocean) as well, has been object of investigation from this point of view (Solheim et al., 2005).
In the Mediterranean Sea, the tsunami hazard connected to continental slopes is still poorly constrained. The comprehensive study by Urgeles and Camerlenghi (2013) represents the first step toward the characterization of MTDs in the whole basin. In total, 696 events have been mapped and described, spanning wide area and volume ranges, (10 − 3 ÷10 5 km 2 ) and (10 − 4 ÷10 4 km 3 ) respectively. Among these, 28 events are reported to have generated tsunamis, both by direct observations and by deposits characterization. While it is not surprising the presence in this subset of huge masses (9 events in the volume range 10÷100 km 3 , 3 exceeding 100 km 3 ), it is significant the incidence in terms of tsunami generation of smaller occurrences: 9 between 1 and 10 km 3 , and 7 below 1 km 3 . Such typology of MTD is scarcely considered in the study of non-seismic tsunami hazard, since in general they generate considerable waves only when occurring in shallow water and in proximity of the coast.
One of the main characteristics of the Mediterranean Sea is the high recurrence of the combination of these two elements: mass wasting features (scars, headwalls, canyons) are recognizable along several margins close to populated coastal communities. Some of them are here recalled: • The Balearic Sea, where one of the most impressive underwater sliding bodies has been found along the Ebro margin, the so-called BIG'95 (Lastras et al., 2005;Lastras et al., 2007), whose tsunamigenic potential has been explored through numerical modeling (Iglesias et al., 2012;Zaniboni et al., 2014a;Løvholt et al., 2014). • The margin of the Ligurian Sea (French-Italian Riviera), with very steep slopes and relatively frequent seismicity that can mobilize sediments (Ioualalen et al., 2014), such as the case of the 1979 Nice tsunami (Assier-Rzadkiewicz et al., 2000). • The Tyrrhenian and Ionian margins, where many mass wasting processes covering different spatial scales have been mapped in the framework of the Italian project MaGIC (Chiocci and Ridente, 2011;Casalbore et al., 2014;Rovere et al., 2014;Casalbore et al., 2019). Among the many potential occurrences, numerical simulations for the study of the generated tsunamis have been performed for the 1977 Gioia Tauro event (Zaniboni et al., 2014b), on the Tyrrhenian Calabrian side, with an estimated volume of approximately 0.005 km 3 volume. On the Ionian side of Calabria, the Assi landslide (Ceramicola et al., 2014), has been object of investigation (see Table 1 for details). Also, the area of Crotone, Calabria, is worth of mention, with the (indeed still questioned) homonymous potential megalandslide involving a very thick sedimentary sequence (Zecchin et al., 2018). • The southern Adriatic Margin, where structures favoring mass transport such as the Bari Canyon  are found and pieces of evidence of vast movements exist, such as the large Gondola slide, a complex of events mobilizing deposits in the order of tens of km 3 (Ridente et al., 2008). The tsunamigenic potential of small landslides on the eastern margin has been investigated as well (see Table 1; Argnani et al., 2011). • The Hyblean-Malta Escarpment (Ionian coast of Sicily), where many canyons and scars are evident (Micallef et al., 2014), and the potential for tsunami generation has been examined . To mention also that the possibility of a 5 km 3 submarine landslide occurring in the occasion of the 1,693 earthquake, that might have enhanced considerably the effects of the earthquake-tsunami at a local scale (Argnani et al., 2012). • The margins close to the coasts of Crete and Cyprus (Papadopoulos et al., 2007b;Papadopoulos et al., 2014) and along the Corinth Gulf, where on 1963 a coastal slump generated relevant waves (Papadopoulos et al., 2007a;Tinti et al., 2007). • The African coast of the Mediterranean, that is still scarcely investigated. Some indications of collapses have been reported in Loncke et al. (2009), describing mass wasting processes covering a large volume range (from 1 km3 to around 500 km 3 ) along the Nile river submarine fan, offshore the town of Alexandria (Egypt). • The Levantine Basin is a place of large mass transport complexes, ranging from 35 to 94 km 3 , occurring along the continental margin off Israel and Lebanon (see Eruteya et al., 2016, and references therein).
The above list constitutes only a subset of the potential margins that are prone to sliding in the Mediterranean Sea and can potentially generate tsunamis. In this paper, we consider the Gela Basin eastern slope (GBES from now on) that is found in the northern part of the Strait of Sicily and extends from the shelf-edge at relatively shallow water (∼200 m) rapidly deepening to about 900 m depth. Studies based on morpho-bathymetric data show that the GBES has been extensively affected by submarine mass wasting during the Late Quaternary, involving volumes of sediments in the order of magnitude of 1 km 3 (e.g., Minisini et al., 2007;Kuhlmann et al., 2017). Among the recorded landslide events, we select two as representative of potential scenarios along the GBES. To evaluate their tsunamigenic potential, a comprehensive investigation through numerical codes is performed, mainly consisting into three phases: the simulation of the landslide dynamics; the computation of the tsunamigenic impulse, that is time-dependent; the simulation of wave propagation, with the assessment of the tsunami hazard for the neighboring coasts of Sicily and Malta. The analysis is limited to the effects of the tsunami on a regional scale, leaving the study of inland flooding for future publications.
The study of these two landslide scenarios along the GBES has been performed in the framework of the Italian project SPOT (Sismicità Potenzialmente innescabile Offshore e Tsunami; Antoncecchi et al., 2020), aiming at assessing the tsunamigenic potential of earthquakes and landslides possibly triggered by hydrocarbon production offshore the Italian coasts.

Geological Setting
The Gela Basin is a bathymetric depression of limited water depth (up to 936 m) located south of central Sicily. It represents the foreland basin of the Maghrebian thrust belt of Sicily Lickorish et al., 1999), and is filled by up to 2000 m of turbidites and pelagic sediments of Pliocene-Quaternary age. The sedimentation rate from Pliocene to Middle Pleistocene (800 ka) was 150 m/Myr and reached 900 m/Myr in the last 800 kyr (Gauchery et al., 2021). The upper part of the sedimentary fill is characterized by abundant mass transport deposits. The northern margin of the basin is partly shaped by the arcuate front of the Gela Nappe (Argnani, 1987), which represents the most recent Maghrebian thrust front ( Figure 1). A sedimentary prograding set, fed from the north, developed on top of the Gela Nappe. This prograding set extends eastward, away from the thrust front, fringing the Hyblean Plateau. The most recent clinoform of this prograding set represents the northern and eastern bathymetric slope of the Gela Basin.
The eastern slope, denoted as GBES here, has been the site of several mass transport events, as evidenced by the abundant slide scars which are visible on the morpho-bathymetry map (see Figures 1B,C), and as reported in detailed studies of selected sectors (e.g., Minisini et al., 2007). The morphological evidence indicates that complex and recurrent sediment failures affected the GBES during the Late Quaternary.

Mass Failures Along the GBES
Mass transport events characterize the whole extent of the GBES, and in some cases, the reconstruction of the sliding mechanism and sequence is quite difficult, due to the superposition of different occurrences.
Starting from the north, the first occurrence is the Gela Slide (GS, see Figure 1A), a 630 km 3 landslide affecting an area of more than 1,500 km 2 , characterized by a few km downslope movement and occurred presumably in the Late Pleistocene. First described in the geological work by Trincardi and Argnani (1990), this collapse has been taken as one of the scenarios in the study of tsunami hazard on the coasts of the Malta archipelago by Mueller et al. (2020). According to their tsunami simulations, the flow depth exceeded 10 m on the island of Gozo, which was hit around 18 min after onset of the landslide.
Moving to the east of the GS, several other landslides have been recognized by Trincardi and Argnani (1990) and in more recent investigations (Minisini and Trincardi, 2009;Kuhlmann et al., 2017). The most interesting cases are the so-called Twin Slides located about 30 km far from the GS (see Figure 1B). These collapses, that occurred probably simultaneously in Late Holocene, are characterized by well-defined scars (see the headwalls, Figure 1B), deepening from 200 to 500 m water depth with well-recognizable deposits of comparable size at the toe (enclosed within the dashed-red and dashed-green boundaries respectively), down to 700 m b.s.l. The estimated volume is slightly less than 0.5 km 3 for both slides. The Northern and Southern Twin Slides (NTS and STS) have been interpreted as the final stage of a very complex sliding sequence, tentatively reconstructed in Kuhlmann et al. (2017), that started with a larger "Father Slide" (black line in Figure 1B) 87 ka ago and that repeats periodically about every 10 kyr.
Another submarine slide is placed southward, about 40 km north of the Malta archipelago ( Figure 1C). In this work, it is named South Gela Basin Slide (SGBS) according to Gauchery et al. (2021). Pieces of evidence of collapse are very clear: a large scar at the shelf-edge, more than 5 km long (blue line in Figure 1C), with a 5°average slope from 200 to 500 m water depth; a large deposit, easily recognizable at the toe of the slope (350-550 m water depth), covering an area of more than 50 km 2 (delimited by the dashed-blue contour in Figure 1C).
In this paper, we have opted to compute the tsunamis produced by two slides out of the many that have been identified along the GBES, namely the NTS and the SGBS. Although there is not any direct piece of evidence of the occurrence and the effects of such tsunamis, nonetheless exploring these scenarios is important mainly because 1) the involved volumes are the most frequently observed in the Mediterranean Sea (see the frequency distribution in Urgeles and Camerlenghi, 2013) and 2) landslides with comparable volumes are known to produce relevant (though local) waves, which stresses the need to systematically include also such sources in tsunami hazard evaluation. Some interesting examples of tsunami studies from this kind of landslides are listed in Table 1.

NUMERICAL METHODS
The study of the tsunamis associated with submarine mass movements along the GBES is performed through a well-tested numerical procedure that has been developed in-house and applied to several cases of landslide-generated tsunamis (Zaniboni et al., 2014b;Ceramicola et al., 2014;Zaniboni et al., 2016;Zaniboni et al., 2019;Gallotti et al., 2020;Triantafyllou et al., 2020). Under the assumption that the submarine slope will fail, the simulation sequence covers the whole process including 1) the dynamics of the sliding motion, 2) the tsunamigenic impulse caused by the movement of the mass along the sea bottom, and 3) the propagation of the tsunami over the computational domain.

Landslide Dynamics
When the sliding body reaches instability conditions, it starts moving along the slope. The dynamics of such motion is computed through the code UBO-BLOCK1, which implements a Lagrangian approach. The sliding body is partitioned into a set of interacting blocks, whose centers of mass (CoMs) motion is determined by the interaction with the surrounding environment, i. e., by the body forces (gravity, buoyancy), the surface stresses (basal friction, water drag on the exposed surfaces), and the internal interactions between blocks. The blocks conserve their mass and cannot penetrate nor superimpose with each other. This approach allows one to quantify how much the slide changes its shape during its descent, a crucial factor in tsunami generation. The simulation is stopped when the mass exits the computational domain, or when the mean velocity lowers a predefined threshold. The application of UBO-BLOCK1 requires as input: i. The undisturbed sliding surface; ii. The upper surface of the initial sliding mass; iii. The predefined CoM trajectory; iv. The lateral boundaries of the sliding surface.
Further details on the code can be found in Tinti et al. (1997).

Tsunami Generation and Propagation
The motion of the sliding body on the sea bottom changes the sea depth and mobilizes the whole water column, generating a perturbation that propagates throughout the water body. The tsunami impulse provided by the slide is not instantaneous, since the time scale of the two phenomena (landslide motionwave propagation) is comparable, in contrast to the process of earthquake-generated waves, where the source can be considered instantaneous.
The tsunamigenic impulse is computed as the time history of the seabed change due to the passage of the mass, over each node of the tsunami computational grid. This perturbation is filtered with the water depth through a function cutting higher frequencies. These tasks are performed by a specific code, named UBO-TSUIMP, which is described in full detail in .
The wave propagation is modeled by the application of the classic non-linear shallow-water equations (SWE), that are solved by a finite difference approach (leap-frog numerical scheme) implementing the staggered grids technique. When the computational grid boundary is the open sea, a pure transparency condition is imposed, while the interaction with the coast is handled in two possible ways: in case of land inundation, the model implements the moving boundary technique, considering the flooded inland cells as part of the bathymetry; when the no-inundation condition is selected, the coast is considered as a vertical wall and the wave is reflected seaward. The choice between the two approaches depends on the aims of the investigation. If one wants to evaluate the tsunami hazard at a regional scale over a wide domain, the second is preferable. If the interest is on the effect on coastal communities and buildings, one should select the first option. The two alternatives usually require domains with different characteristics: high resolution to compute inundation, low resolution to simulate propagation in oceanic regions. The code, named UBO-TSUFD, includes also the possibility to manage domains with different spatial steps implementing the nested-grid technique, useful if a heterogeneous resolution is convenient for the simulation.
The input datasets needed to run the code are: i. The tsunami computational grid, or set of grids; ii. The tsunami initial condition.
This code is more extensively described in Tinti and Tonini (2013). Its applications are reported in the same references cited in the previous section. Besides, one can find an application to a case of an earthquake-tsunami in Loreto et al. (2017).

BUILDING THE LANDSLIDE SCENARIOS AND THE TSUNAMI COMPUTATIONAL GRID
In this paper we study the NTS and the SGBS since we consider their volume as typical of tsunamigenic mass movements in this region, and, more in general, because events of this size require

Landslide Scenario for NTS
The simulation of the landslide motion requires the definition of the four elements listed in Landslide Dynamics Section. These have been devised mainly on morphological considerations, starting from the present seabed bathymetry given in Figure 2. In the NTS scenario, the sliding surface (green line in Figure 2B) coincides with the slide scar in the steeper part of the margin, uncovered after the sliding event down to 500 m depth. The remaining portion, now hidden by the slide deposit, has been inferred by extending the outer isobaths inside the sliding boundary (green boundary, Figure 2A) since these are supposed to represent the undisturbed surface under the slide deposit. After the sliding surface has been reconstructed, it is straightforward to obtain the sliding deposit, simply by difference with the present bathymetry (red line, Figure 2B). Though not essential for the simulation itself, the observed final distribution of the mass is very useful as a constraint for the parameters governing the sliding model. The initial mass has been obtained simply by filling the scar, again extending the isobaths inside the initial landslide contour (blue line, Figure 2A). The result is a body with volume of 0.46 km 3 , consistent with the deposit at the toe of the slope, covering an area of almost 7 km 2 . The initial slide mass distribution, obtained as the difference with the gliding surface, evidences a thicker central part, reaching 150 m ( Figure 2A). The CoM path follows the direction of local maximum steepness, while the slide lateral boundaries (Figure 2A, in green) follow the observed slide deposit (red line, Figure 2A).

Landslide Scenario for SGBS
The second landslide scenario selected is placed around 50 km south of the NTS, closer to the Malta archipelago. The same procedure followed for the NTS has been used for the preparation of the SGBS simulation input. The sliding surface follows the uncovered steeper part of the slope and under the present deposit, is inferred by continuity with the outer isobaths (green profile in Figure 3B). The initial mass is obtained by filling the scar inside the respective boundary (blue boundary in Figure 3A), providing an initial sediment body with a volume of 1.48 km 3 (three times bigger than NTS), over an area of more than 26 km 2 (four times larger than NTS), implying a smaller thickness (maximum less than 100 m, Figure 3A). The final deposit is obtained by the difference between the present morphology and the sliding surface. As in the previous case, the obtained volume has been used as a further constraint on the reconstruction of the initial mass, since the two amounts have to be compatible.

Tsunami Computational Grids
The simulations of the landslide-generated tsunamis require a regularly spaced computational grid. The investigated area is shown in Figure 4: the larger grid (G1, in green) covers the southern corner of Sicily, more specifically the SE coast watered by the Strait of Sicily to the south and by the Ionian Sea to the east. This grid has been built with a spatial step of 250 m. The source data come from the EMODnet public database, covering this area with a resolution of about 115 m. Due to the limited extent of the tsunami sources, that would have been described by too few nodes with the resolution of grid G1, it was necessary to make use of finer domains, covering the two tsunami generation areas. Grid G2, marked in blue in Figure 4, accounts for the SGBS scenario and includes the Malta archipelago as well; grid G3 (in red) covers the NTS case and the coast of the Gulf of Gela. Both computational grids have been built with a 50 m spatial step by combining a swath bathymetry dataset acquired with 30-100 kHz multibeam systems and an overall 50 m resolution along the slope (Gauchery et al., 2021), with the already cited EMODnet bathymetry used to fill the data gaps. This grids configuration (a larger 250-m mesh, G1, and two smaller 50-m domains, G2 and G3) allows an acceptable and detailed reconstruction of the landslide dynamics and their tsunamigenic impulse in the source areas, and sufficient coverage of the tsunami propagation in the Strait of Sicily. The simulation with this grid set is possible through the nested-grid technique implemented in the model UBO-TSUFD, allowing to account for the passage of the tsunami wave across boundaries of contiguous computational domains with different resolutions.

LANDSLIDES SIMULATIONS
The code UBO-BLOCK1 provides the complete time-history of the motion of each block composing the landslide mass. As a consequence, the simulation accounts for the mass shape variation during the descent, a factor that deeply influences the perturbation produced on the water column, necessary to evaluate the time-dependent tsunamigenic impulses associated with the mass moving along the seabed.
The comparison with the observed deposit provides important constraints on the simulation parameters governing the slide motion. Concerning the friction coefficient, a value of 0.03 provided the best fit: this is a typical value for submarine landslides. The drag coefficients have been selected basing on values coming from previous applications, simulating similar failures. For the superficial stress, the value chosen for the drag parameter is C d 0.02; as for the frontal drag, C f 0.5 (for a more detailed description of these coefficients and of their range of values, refer to Tinti et al., 1997). Both sliding bodies are affected by a considerable elongation, as inferable also from the profiles of Figures 2,3: NTS passes from 4 km to more than 8 km at the end of the motion; SGBS from 6 km to more than 12 km. The values adopted for the internal interaction parameters, governing the mass deformation, have been tuned to account for this behavior. Figure 5 reports some of the motion features of the two scenarios. Panels A and B report the velocity evolution with time: it can be noticed that the two curves are similar, with an initial abrupt acceleration phase followed by a slower deceleration, typical of masses moving along steep margins, and then reaching the flat area at their toe. The NTS reaches the mean velocity peak (15 m/s) after around 100 s, while for SGBS the peak is slightly smaller, and attained at around 200 s. The black dots mark the individual CoM velocity record: notice that there is a generalized spread around the mean values (continuous line, red for NTS, and blue for SGBS), which is a natural consequence of the mass elongation during the descent. The mass comes to rest after almost 700 s for NTS and more than 850 s for SGBS, but it can be noticed that some blocks stop much earlier (already after 300 s for NTS; at 500 s for SGBS), while other have still residual velocity when the simulation is stopped.
The block thickness evolution is represented in Panels C and D of Figure 5: reflecting the elongation and spreading at the end of the motion, the sliding mass gets thinner, passing from an average of 75 m to around 20 m for NTS, and from 58 m to 31 m for SGBS. It is noticeable that in correspondence with the velocity peak, the block thickness increases considerably, especially for the NTS where some blocks reach height values of more than 120 m, with obvious influence on tsunami generation.
A typical indicator of the tsunamigenic efficiency is the Froude number, Fr. This is computed as the ratio between the horizontal component of the slide velocity and the wave phase speed ( gh , with g gravity acceleration and h local sea depth). When the two quantities are equal, the energy transfer from the slide to the water is maximum, and the resonance condition is attained (Fr 1). For values lower than 1 (subcritical condition, typical of deep slides) the wave moves faster than the slide; for Fr > 1 (supercritical condition, typical of fast subaerial slides entering the water, usually occurring close to the generation area) the mass runs away from the generated wave. In the present cases, one can notice that the motion is always subcritical (panels E and F, Figure 5). The maximum

TSUNAMI PROPAGATION
Tsunami simulations have been run with the linear version of the code UBO-TSUFD, without considering land inundation. Concerning the Malta archipelago where coasts have complex morphology, this task would require more refined grids, to describe better the many small gulfs and inlets characterizing especially the northern coast.
The propagation of the NTS-generated tsunami is shown in Figure 6. In the first frames (4 and 8 min) the typical circular radiation of landslide-tsunamis can be noticed, with a positive front propagating in the same direction as the slide motion (offshore, toward south-west) and a negative signal on the opposite side. This entails that the tsunami firstly manifests with a sea withdrawal at the Sicily coast, a factor that can be important in terms of alert management. The wave hits the coasts after about 12 min, with the negative front, which is soon followed by a positive wave, meaning sea-level increase. Within 20 min the whole Gulf of Gela (see Figure 2 for geographic location) is affected by the wave, the same happens for Gozo (the north-westernmost island of the Malta archipelago), that in contrast to Sicily is hit by a positive wavefront. Notice that the travel time is very similar to the one obtained in Mueller et al. (2020) for a much bigger mass, the Gela Slide. The 36-min sketch shows that at this time almost the entire Malta islands are affected by the tsunami. Finally, notice the strong deceleration effect of the Malta Plateau, the shallow-water area between Sicily and Malta (whose location is shown in Figure 4 as well), on the wave propagation due to bathymetry and shoaling mechanisms.
The SGBS-tsunami propagation sketches (Figure 7) show some similarities with the NTS case, especially in the first frames, i. e., the circular radiation from the source and the polarity of the wavefront (positive westward, negative eastward). Conversely, the SGBS tsunami is generally higher than the NTS one. Gozo island is attacked by a positive front between 8 and 12 min, and the whole Malta archipelago within 28 min. Considering the Sicily coasts, one can notice that already at the 12 min frame the wavefront tends to deform, due to the interaction with the bathymetry. The first affected coastal stretch is the southern extreme of the Gulf of Gela, between 20 and 28 min. At this time a positive signal reaches the coast. Here the tsunami signal deceleration in the Malta Plateau area is evident as well, with the negative front (in blue) taking several minutes to cross this sea region.
Interesting insights on the tsunami characteristics come from Figure 8, representing the maximum sea surface elevation computed for each node of the tsunami grid, during the whole tsunami propagation history. This plot provides a useful glance at the spatial distribution of the tsunami energy. One can notice that most of the tsunami energy in the NTS scenario is captured within the Gulf of Gela. Some tsunami beams are visible, evidencing preferential directions for water maximum elevation. Noticeable are two rays hitting the central part of the gulf, two more affecting its eastern end, and another one moving south-eastward. The Malta archipelago coasts seem scarcely affected by waves generated from the NTS scenario.
Concerning the SGBS (Figure 8right panel), the pattern is quite different. The Malta islands are hit by relevant maximum waves, reaching also 2 m. A strong beam heads towards the Malta Plateau, south-east, but the most noticeable feature is the high concentration of energy directing towards the coast of Sicily, east of the Gulf of Gela. Here the water elevation exceeds 3 m. Conversely, the coasts of the gulf itself are moderately protected, since most of the tsunami energy is attracted to the east. The SGBS scenario, though more distant, produces more relevant and diffused effects on the Sicily coasts than the NTS scenario.
These observations are confirmed by Figure 9, representing the maximum water elevation along the coast of Sicily. The water height is computed along the 5-m isobath since the linear version of the simulation code with fixed coastal boundary has been run and no inundation has been computed. Therefore, the results described here can be considered as underestimations of the effective run-up heights in terms of hazard, also considering that non-linear terms would amplify the waves when approaching the coast. Figure 9 confirms that the NTS tsunami is mainly confined within the Gela Gulf (between points #1 and #4), where a maximum of 3 m is reached close to Marina di Acate (node #3) and another 2 m peak can be observed around node #4 (Santa Croce Camerina). Out of this coastal stretch, about 80 km long, the wave height rapidly drops below 1 m and is almost negligible west of Licata (#1) and beyond Capo Passero (point #6).
The curve representing SGBS coastal water height, on the other hand, shows the maximum elevations reached between points #4 and #5, i. e., the area towards which the tsunami beam mentioned above (Figure 9) is directed. Here waves reach 3.5 m and remain higher than 2 m for at least 30 km along the coast. Within the Gulf of Gela, the water elevation oscillates around the value of 1.5 m, comparable then to the NTS case. West of Licata (#1) the wave height remains in the range of 1-1.5 m for almost the entire coast. On the east, similar behavior is observed, with oscillations exceeding 1.5 m towards Capo Passero (node #6) and also farther, along the eastern coast of Sicily. In general, the SGBS affects the coast of Sicily with waves exceeding 1 m for almost the whole coastal stretch covered by the simulation, 300 km long, dropping down to 1 m elevation only at the plot borders. Focusing the attention on the Malta archipelago, we consider the maximum water elevation along the 20-m isobath rather than along the 5-m isoline ( Figure 10). This choice is motivated by the fact that the bathymetry in the shallow coastal zone here is quite complex, being characterized by numberless small inlets, bays, promontories that could be well described only by high-resolution grids, which in turn would require heavier computational costs. The results we will show on the 20-m isoline are expected to be underestimates of the maximum tsunami waves since they do not account for possible resonance amplification nor tsunami energy focusing. Nonetheless, they provide a good basis to evaluate the overall impact of the tsunami on the Malta islands.
The NTS tsunami produces limited effects, mainly on the northern part of the Gozo island, where 0.4 m height is reached.
The inlets along the northern coasts are affected by waves at most of some tens of centimeters, that can be considered negligible in terms of human hazard but can be amplified by resonance effects, producing heavy damage on boats and harbor facilities.
The SGBS source area is closer to the Malta archipelago, producing relevant effects especially along the northern coasts. Gozo island is the most affected, with the north-western coast hit by waves exceeding 2 m. Also, the northern half of the main Malta island is impacted by a wave at least 1.5 m high, with relevant waves entering St. Paul's Bay (point #1 in Figure 10), the most populated town of the island, with a substantial increase of population in the tourist season. Also, the inlet of La Valletta (point #2 in Figure 10) is affected by waves higher than 1 m. All the southern coast, on the contrary, seems protected with maximum waves barely reaching half a meter.

CONCLUSIONS
In this paper, we have studied two cases of tsunamis produced by mass failures along the GBES, selected mainly for two reasons. First, high-quality seafloor geomorphological data are available, accurately describing the scar and the slide deposit, which allows a suitable reconstruction of the landslide sources and provides useful constraints for the simulation model. Moreover, these slides can be considered representative of typical failure episodes along the slope.
In the framework of the SPOT project, aimed at assessing the influence of local earthquakes on tsunami generation and mobilization of sediments, the stability analysis of the two landslide scenarios have been performed (see results in Supplementary Appendix A). They showed that failures in the northern sector of the GBES, where the NTS is located, can be activated with a return period of few thousands of years. More in the south, farther away from the seismic faults, at least according to the present knowledge of the offshore tectonics, failures like the SGBS cannot be explained by invoking seismic load, and other destabilizing causes have to be found, which implies further research efforts.
Mass movements along the GBES require particular attention since they occur in relatively shallow water. Moreover, in the initial stage, they soon attain large velocity, due to the high steepness of the continental slope (5°-10°). The combination of these two elements (large velocity and shallow water) enhances the tsunamigenic potential in a considerable way.
The tsunami simulations show that masses such as the NTS and the SGBS can produce relevant waves impacting coastal stretches from tens to hundreds of km long. For the NTS, the arrival of the tsunami on the Sicily coasts manifests as a sea withdrawal. In this case, a useful precursor could be the earthquake, destabilizing the underwater body since it anticipates the tsunami arrival by some tens of minutes. This latter kind of phenomenon sometimes called "surprise tsunami" (Ward, 2001) needs increasing attention and continental margins should therefore be more extensively mapped, investigated and possibly monitored, especially off those coasts of the Mediterranean Sea that are densely populated. The continuous supply of sediments by rivers, the diffused seismicity in the whole basin, the presence of other destabilizing phenomena (such as strong submarine bottom currents and volcanoclastic material) are elements contributing to increasing the potential hazard.

Further Improvements and Perspectives
The GBES is particularly interesting for many reasons. First of all, it is placed in relatively shallow water, connecting a wide shelf at 200 m depth to the deeper sea (800-900 m). This permits the characterization of the morphology in detail, both for the reconstruction of past events and for the recognition of possible future occurrences. Many submarine landslide events can be mapped by high-resolution bathymetric data and described with sufficient accuracy, covering a large spectrum of volumes and return periods, making it possible to apply probabilistic approaches. The study of the tsunami effects on the coast should also consider land inundation and the impact on coastal communities. This issue was the scope of the Italian project SPOT (Antoncecchi et al., 2020) for the coast of Sicily, but not for the Malta archipelago. It was addressed by employing empirical laws for tsunami flow over coastal 1D transects and the results are left to further publications. As regards the Malta coasts, the inlets along the northern coasts, in particular, favor wave amplifications that should be investigated through more detailed computational grids. The two scenarios explored in this work, together with the ones reported in Mueller et al. (2020), could be integrated into a more comprehensive study of tsunami hazard of the Malta archipelago.
An exhaustive study of the tsunami hazard related to collapses along the GBES would require more scenarios, covering a larger spectrum of volumes, also using a probabilistic approach. Nonetheless, the authors are convinced that the two cases presented here provide interesting insights for the evaluation of the tsunami hazard along the coast of Sicily and Malta and that will stimulate the interest on these phenomena, which should require increasing consideration.
Validation of results from numerical simulations can come from the analysis of tsunami deposits in the sedimentary sequences on land and also offshore. This was the case, for example, of the geological investigations by Pantosti et al. (2008), De Martini et al. (2012 and Smedile et al. (2020), that were able to associate sediment layers found in trenches or cores in eastern Sicily to historical-as well as paleo-tsunami cases. The GBES morphology suggests that mass collapses along the slope repeat cyclically, and the generated tsunami can reach the coasts of Sicily and Malta with potentially relevant waves. The events investigated here are prehistorical. Finding their coastal signature is a hard, but not an impossible task and future research by tsunami geologists could fill this gap.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
FZ realized the datasets describing the landslide scenarios; contributed to develop and maintain the landslide simulation code; performed the landslide dynamics simulations; wrote the first draft of the paper; produced most of the figures in the paper; collected and harmonized all the contributions from coauthors; managed the manuscript progress. GP contributed to develop and maintain the tsunami simulation code; built the tsunami computational grids; run the simulations of the tsunami propagation. MAP contributed to the definition of the landslide scenarios; contributed to develop and maintain the stability analysis code; performed the slope stability analysis. TG contributed to the definition of the SGBS scenario and the discussion of the collapse triggering factors; performed the literature review on landslide-tsunamis; revised the paper. MR provided bathymetric data for landslide scenarios and tsunami computational grids; contributed to the definition of the landslide scenarios; supervised the geological part of the work; produced Figure 1; contributed to the manuscript revision. AArg: contributed to the definition of the landslide scenarios; supervised the geological part of the work; contributed to the manuscript revision; produced Figure 1. AArm contributed to develop and maintain the tsunami simulation code; supervised the modeling part of the work; revised the paper. ST developed the stability, landslide, and tsunami numerical codes; supervised the whole work; revised the paper.

ACKNOWLEDGMENTS
The bathymetric metadata and Digital Terrain Model data products have been derived from the EMODnet Bathymetry portal -http://www.emodnet-bathymetry.eu. The PGAs in the landslide areas originated by the selected faults were provided by Barbara Borzi and her research team (Eucentre), in the framework of the SPOT project. The SPOT project was funded by the Italian Ministry of Economic Development, Directorate General for Safety-National Mining Office for Hydrocarbons and Georesources (DGS-UNMIG) under the umbrella of the Offshore safety network "Clypea," which is in force since 2014, and was established with bilateral agreements between the Ministry and Research Centres, Governmental Bodies and Universities to increase the safety, also in terms of environmental protection, of offshore plants operations. The Ph.D. of Tugdual Gauchery is funded by the European Union's Horizon 2020 research and innovation program under the Marie-Skłodowska-Curie grant via project ITN SLATE (Grant Agreement No 721403). The authors are indebted to the editor and the reviewers, Roger Urgeles and Juan Horrillo, whose comments and suggestions contributed to considerably improve the manuscript.