Assessment of the 1783 Scilla landslide–tsunami's effects on the Calabrian and Sicilian coasts through numerical modeling

The 1783 Scilla landslide–tsunami (Calabria, southern Italy) is a well-studied event that caused more than 1500 fatalities on the beaches close to the town. This paper complements a previous work that was based on numerical simulations and was focused on the very local effects of the tsunami in Scilla. In this study we extend the computational domain to cover a wider portion of western Calabria and northeastern Sicily, including the western side of the Straits of Messina. This investigation focuses on Capo Peloro area (the easternmost cape of Sicily), where the highest tsunami effects outside Scilla were reported. Important tsunami observations, such as the wave height reaching 6 m at Torre degli Inglesi and flooding that reached over 600 m inland, have been successfully modeled but only by means of a high-resolution (10 m) topo-bathymetric grid, since coarser grids were inadequate for the purpose. Interestingly, the inundation of the small lake of Pantano Piccolo could not be reproduced by using today's coastal morphology, since a coastal dune now acts as a barrier against tsunamis. Historical analysis suggests that this dune was not in place at the time of the tsunami occurred and that a ground depression extending from the lake to the northern coast is a remnant of an ancient channel that was used as a pathway in Roman times. The removal of such an obstacle and the remodeling of the coeval morphology allows the simulations to reproduce the tsunami penetration up to the lake, thus supporting the hypothesis that the 1783 tsunami entered the lake following the Roman channel track. A further result of this study is that the computed regional tsunami propagation pattern provides a useful hint for assessing tsunami hazards in the Straits of Messina area, which is one of the most exposed areas to tsunami threats in Italy and in the Mediterranean Sea overall.

1 Introduction

The recent catastrophic tsunamis of Sumatra (2004), Japan (2011), and Sulawesi (2018) have raised the interest in such natural phenomena worldwide, including Europe. In the Mediterranean Sea, tsunamis are known to be of smaller magnitude than in the Pacific and in the Indian Oceans, but their effects can be as lethal, owing to the high coastal exposure and vulnerability, which has been constantly growing in recent decades as the result of an increased coastal occupation (see Papadopoulos et al., 2014). This creates the need for more detailed assessments of tsunami hazards and of the consequences of tsunami impact, which also implies the need for more accurate numerical simulation tools.

Figure 1 View of Scilla and its surroundings including the two main beaches of Marina Grande and Chianalea di Scilla (see the position on the map of Italy on the left and in the zoom in the bottom-left corner). The yellow contour delimits the landslide scar, measuring approximately 750 m × 400 m, about 1 km west of the town. The red contour marks the submarine deposit observed in the geo-marine surveys (Bozzano et al., 2011; Casalbore et al., 2014). Satellite view obtained by Google Earth ©.

In this framework, the reconstruction of historical events that significantly affected the coasts of the Mediterranean Sea is very important, since they can be used to test numerical simulation codes as well as for better estimates of the hazard and risk. One of the most relevant cases in the Mediterranean, owing to the disastrous impact and the unusual availability of plenty of very detailed observations, is the 6 February 1783 landslide-induced tsunami in Scilla. This tragic event, causing the deaths of more than 1500 people, was the subject of several coeval reports, providing details on the devastating effects along the littoral but also on the coastal landslide that was the cause of the tsunami (see Graziani et al., 2006; Tinti et al., 2007; Gerardi et al., 2008). In a previous paper (Zaniboni et al., 2016) we investigated the tsunami's impact close to the source area, which is on the beaches of the town of Scilla (Marina Grande and Chianalea di Scilla), where most of the tsunami fatalities were counted; see Fig. 1. This was accomplished through numerical simulations under the assumption that the causative landslide was purely subaerial.

Table 1 Summary of pieces of evidence of the 6 February 1783 tsunami in Calabria and Sicily, with the exclusion of local effects in the area of Scilla. For geographic locations, see Fig. 2.

As reported by eyewitnesses (descriptions and references are summarized in Table 1), tsunami waves affected at least as much as 40 km of the Calabrian coast, the cape of Sicily facing Scilla, named Capo Peloro, and the northeastern coast of Sicily, just south of Messina. In this paper, we carry out three tsunami simulations, one to cover coasts located farther away from the source and two to study the inundation in the area of Capo Peloro on the opposite side of the Straits of Messina. For this purpose, we have performed a preliminary reconstruction of the coeval coastal morphology and topography of Capo Peloro by using features and observations recovered from historical sources. Comparison between numerical results and observations will allow us to evaluate the quality of the simulations.

2 The 6 February 1783 landslide–tsunami

In the first months of 1783, a catastrophic seismic sequence, with at least five main events, hit the Calabria region (southern Italy), provoking about 30 000 casualties and the destruction of over 200 villages (Guidoboni et al., 2007). Cascade effects were also reported, including tsunamis and diffuse landslide occurrences in the Apennines that form the backbone of the region (Tinti and Piatanesi, 1996; Graziani et al., 2006; Porfido et al., 2011).

The most tragic event took place during the night of 6 February, in the town of Scilla on the Tyrrhenian coast of Calabria (see map in Fig. 1). The day before, one of the strongest earthquakes of the sequence caused severe destruction in the village, forcing most of the population to find refuge on the wide beach of Marina Grande, west of the town. A strong aftershock in the night caused the collapse of the coastal flank of Mount Pacì in the sea, 1 km west of the Marina Grande beach. The ensuing tsunami killed more than 1500 people (Mercalli, 1906) and destroyed houses and churches (Tinti and Guidoboni, 1988; Graziani et al., 2006).

The huge death toll and the impression raised by the disaster induced the Bourbon government ruling the region at that time to finance and arrange several ad hoc surveys. Further, such a calamity was seen as a big challenge by the coeval international scientific community that was engaged in understanding the nature of earthquakes. Therefore, a number of additional investigations were carried out in the first few years after the tragedy by scholars and travelers from different countries (see notes about the mainshock in Guidoboni et al., 2007). The overall result of these observations and research activities that is relevant for our study is that a description of the tsunami source and of the tsunami consequences was made available for future generations (see the discussion in Zaniboni et al., 2016).

Figure 2 Toponyms reported in Table 1 (yellow circles). The area swept by the 1783 landslide is delimited in red.

In this work, our attention is focused on the tsunami impact outside the surroundings of Scilla. The most relevant reported effects are listed in Table 1 and are mainly gathered from recent papers that quote historical sources (i.e., Graziani et al., 2006; Gerardi et al., 2008; the Italian tsunami catalog by Tinti et al., 2007). The last column includes, among others, the main coeval references. The affected locations are shown in Fig. 2.

From Table 1, one can observe that the effects of the tsunami were much less severe in the neighboring area than in Scilla. This is not surprising since it is typical of landslide–tsunamis to attenuate rapidly with distance (see, for example, Masson et al., 2006 and Harbitz et al., 2013, and references therein). However, remarkably the waves were seen in Calabria from Nicotera (point 1, Fig. 2), about 40 km northeast of Scilla, to Reggio Calabria (point 5, Fig. 2), 20 km to the southwest. On the other side of the Straits of Messina, the place closest to the source is the easternmost part of Sicily, which is Capo Peloro (point 6, Fig. 2), next to the village of Torre Faro (point 7). Here the tsunami was strong and disastrous. In addition to historical sources, the impact of high-energy tsunami waves in this area was confirmed by recent geological investigations carried out in a site called Torre degli Inglesi (Capo Peloro, point 6), about 40 m away from the modern shoreline, where a 15 cm thick sand deposit in a trench layer sequence was attributed to the 1783 event (Pantosti et al., 2008).

3 Numerical methods

The numerical procedure and techniques used in this paper were developed by some of the authors and presented in previous applications (see, e.g., Tinti et al., 2011; Argnani et al., 2012; Zaniboni et al., 2013, 2014a, b, c; 2016; Ceramicola et al., 2014). Therefore, only a brief description will be given here, since the reader can refer to the quoted papers.

3.1 The landslide motion

The numerical code to simulate the landslide motion, UBO-BLOCK1, adopts a Lagrangian approach. The mass is partitioned into a chain of contiguous blocks, and their centers of mass (CoM) are taken as representative points of the landslide. Blocks are allowed to change shape (length, width, and height, the last one strongly influencing tsunami generation) depending on material rheology, while the total volume is conserved. Bulk (gravity, buoyancy) and surface (bottom friction, drag) forces acting on the blocks are accounted for, together with a block–block interaction governed by some parameters that regulate the sliding behavior during the descent. More details can be found in the aforementioned applications and in the original description given in Tinti et al. (1997).

3.2 Tsunami generation and propagation

The propagation of the waves induced by the landslide moving underwater is modeled through the Navier–Stokes equations in the nonlinear shallow-water approximation and is solved by the code UBO-TSUFD by means of a finite-difference scheme over a regular grid covering the domain of interest. An intermediate code, UBO-TSUIMP, provides a mapping between the Lagrangian and Eulerian grids (used, respectively, by the codes UBO-BLOCK1 and UBO-TSUFD). Furthermore, it computes the instantaneous tsunamigenic impulse transferring the sea bottom perturbation due to the landslide to the sea surface by means of an appropriate filter function. This term, which is time dependent, is included in the hydrodynamic equations as a source term. The tsunami simulation code accounts also for coastal inundation. This is accomplished by adopting the moving boundary technique (the instantaneous shoreline being identified as the dynamic boundary between wet and dry cells) and by using a two-flow grid nesting scheme where a higher resolution is needed. This is usually the case in near-shore regions for obtaining accurate inundation results. An extended description of the model UBO-TSUFD can be found in Tinti and Tonini (2013).

4 The landslide simulation

4.1 Scenario reconstruction

The surveys carried out soon after the 1783 event described a wide scar along the Mount Pacì flank, 1 km west of Scilla, which is still visible today (see Fig. 1, yellow contour). Minasi (1785) reported a mass failure starting from 425 m a.s.l. and a landslide front penetrating 100 m into the sea. A later source (De Lorenzo, 1895) mentioned a bigger mass with a larger front (up to 2 km) but a similar offshore mass penetration.

Modern studies performed by means of geophysical surveys and numerical modeling seem to reinforce the hypothesis of a purely subaerial event, though some elements suggest the possibility of a subaerial–submarine collapse. Interpreting land and submarine survey data, Bosman et al. (2006) proved that the scar continues underwater, at least down to 100 m depth. This feature is reported also in Casalbore et al. (2014) and can be seen as evidence that the landslide also involved a submarine failure. However, it can alternatively be interpreted as the effect of the erosion of the sliding surface by the motion of a purely subaerial mass, which is the picture we adopted in the previous paper (Zaniboni et al., 2016) and we maintain herein.

Figure 3 Scilla 1783 slide simulation by means of the code UBO-BLOCK1. (a) Profile of the initial sliding body (in blue) and of the final simulated deposit at about 300 m depth (in red) over the undisturbed sliding surface (grey area). The sea level is shown in light blue. (b) Acceleration of the individual CoM (black circles) and of the CoM of the whole slide (green line) vs. distance along the sliding profile. (c) CoM velocities (black dots) and average velocity (green line) plotted on the sliding track distance.

Further, it is worth mentioning that Bozzano et al. (2011) applied a stress–strain numerical model to evaluate the response of the Mount Pacì slope to the shaking of the 5 and 6 February earthquakes and concluded that, though the slide scar admittedly extends underwater, the tsunamigenic source was a purely subaerial collapse, with the toe located at about 150 m a.s.l., involving a volume of 5.4 million cubic meters. They commented that the extension of the scar beyond the initial slide foot might be explained in two ways: either it was the effect of the already cited erosion, produced by the sliding material on the soft surface sediments nearshore, or it was due to a minor submarine failure that occurred before the main slide. Further, Mazzanti and Bozzano (2011), using an oversimplified single-block landslide model, simulated the subaerial tsunami scenario. Avolio et al. (2009) used cellular automata to simulate subaerial and submarine landslide dynamics, without paying attention to the generated tsunami.

As regards the landslide deposit offshore, it is sufficient to mention that various submarine surveys found an underwater blocky deposit at about 300 m b.s.l. (see Fig. 1, dashed red contour), located less than 2 km away from the coast, and associated it to the 1783 slide. The area covered by the deposit is around 700 000 m² and can be delimited with a certain degree of accuracy (see Mazzanti and Bozzano, 2011). The total volume is, however, difficult to assess. According to Bozzano et al. (2011), it is compatible with the subaerial scenario landslide; instead, Casalbore et al. (2014) estimated it in the range of 6–8 million cubic meters.

For the sake of completeness, we also mention that the 1783 Scilla landslide evolution has been the object of a recent application (Wang et al., 2019), where the same subaerial slide scenario has been investigated, mainly to test a new numerical method for computing the landslide motion based on a Eulerian fluid solver using a second-order central scheme with a minmod-like limiter: in that paper the final deposit was used as a constraint to test the accuracy of the landslide simulations, with no attention given to the generated tsunami.

4.2 Results and discussion

In this paper, we make use of the same landslide simulation illustrated by Zaniboni et al. (2016). It is, therefore, sufficient to synthesize only the main features of the landslide motion here, as reported in Fig. 3.

The sliding body (Fig 3a, in blue), reconstructed by considering that the crown is placed at about 400 m a.s.l. and the toe is close to the coastline, results in a subaerial volume of over 6 million cubic meters, with a maximum thickness of about 100 m in the central part. The main parameter governing the motion is the friction coefficient: in view of the comparison between the final simulated position and the observed deposit boundary, the back-analysis provided the values 0.23 and 0.06 for the subaerial and underwater portions of the slide motion, respectively. The red profile in Fig. 3a shows the simulated deposit at about 300 m depth.

Figure 3b and c show the acceleration and velocity of the CoM of the landslide blocks (black dots) as a function of the distance along the profile, together with the average values (green line). Observe that the initial positive acceleration phase, with values around 3 m s⁻², is limited to the subaerial part of the motion (the first 700 m of the profile approximately), while the acceleration adjusts around slightly negative constant values once underwater, causing a slow uniform deceleration phase. This is confirmed by the velocity profile: the average velocity peak is reached just before the entrance of the mass into the water. Notice also that all the CoM have a similar velocity evolution with distance. The mass has an initial length of around 700 m. For most blocks, the runout is about 2 km, with the exception of the frontal one that travels more than 2.2 km, which results in a lengthening of the slide by about 200 m.

5 Tsunami propagation and impact on the coast

The first step of tsunami simulations is the assemblage of the computational grids; the area covered and the adopted grid resolution depend deeply on the goal of the simulation. In general, the higher the grid resolution, the more accurate the results are but also the heavier the related computational effort is. The tsunami simulations in Zaniboni et al. (2016) were concentrated in the town of Scilla, and the extent of the very local regular mesh was a few kilometers in length and width with 10 m spaced nodes (see Grid 1, in green, Fig. 4).

Figure 4 Computational grids adopted for numerical simulations of tsunami propagation: Grid 1 (green) was used for the tsunami simulations by Zaniboni et al. (2016); Grid 2 (in red) was used in the first simulation of this paper; Grids 3 and 4 (in blue) are an extension of Grid 1 covering the target area of Capo Peloro. Coordinates are in UTM-WGS84 projection, zone 33^∘ N.

In this work, our purpose is twofold: (i) to compute tsunami propagation over a wider area, involving tens of kilometers of the Calabrian and Sicilian coasts; and (ii) to compute detailed inundation in the specific area of Capo Peloro, located in front of Scilla. For the first goal, we have built a lower-resolution grid (50 m space step), denoted in Fig. 4 in red as Grid 2, by using three different databases: SRTM for topography (about 90 m spaced, Jarvis et al., 2008); EMODnet DTM, in its 2018 version, providing bathymetry every ∼160 m; and the nautical chart N.138 from Istituto Idrografico della Marina. For this latter case, we have digitized isolines with the aim of better characterizing the shallow water area, which was not adequately described by previous datasets. The resulting computational mesh covers an area of 30×27 km². As concerns the Calabrian coast, it runs from Bagnara Calabra, 10 km west of Scilla to Reggio Calabria, about 20 km southwestward. Regarding Sicily, the grid includes the city of Messina and 7–8 km of coast southward, as well as a piece of the northern coast facing the Tyrrhenian Sea via an extension of at least 10 km.

The choice of Capo Peloro as the target area for the grid resolution enhancement will be justified later, in the discussion of the simulation results obtained on the coarser domain. Figure 4 illustrates the extent of Grids 3 and 4 (how they differ from one another will be explained later as well) in blue that are actually a westward extension of Grid 1 with the same spatial resolution (10 m).

Before discussing the simulation results, it is worth providing some terminological specifications: with “wave height” or “wave elevation” offshore, we mean the difference between the offshore water top and the mean sea level. The “flow depth” is the height of the water with respect to the ground, onshore. Consequently, the onshore “wave height” is the sum of flow depth and local topographic level. Finally, with “run-up” and “inundation distance” we denote the maximum onshore ground elevation and the maximum onshore penetration (approximately normal to the shore) reached by the tsunami, respectively.

Figure 5 Maximum tsunami elevation, computed for each point of Grid 2. Numbers represent the positions of the available observations (see Table 1). The respective toponyms are reported in the legend, together with the observed inundation distance (I) and run-up (R). The dashed blue rectangle marks the limits of the 10 m resolution Grids 3 and 4, zooming on the area of Capo Peloro, that differ from one another only in a small part of their topography. The purple boundary delimits the sliding surface of the 1783 Scilla landslide. For geographical coordinates, see Fig. 4.

5.1 Simulation over the 50 m spaced grid

Figure 5 shows the maximum water elevation, obtained by considering the maximum value of the sea surface height reached during the simulation time span in each node of Grid 2. This category of plots provides a general spatial pattern of the tsunami energy distribution that, in the intermediate- and far-field, depends more on the bathymetry than on the source type. Therefore, it allows one to distinguish the areas that are more prone to tsunami occurrence.

The area close to the slide shows, as expected, the highest water elevations with more than 8 m in Scilla and waves at least 4 m high (in red) affecting the surrounding coast for about 6 km. The wave height distribution illustrated in Fig. 5 also provides indications on the extension of the coast affected by 1.5 to 2 m (in yellow) and 1 m maximum wave (in green), resulting in approximately 10 and 20 km, respectively, along the Calabrian coast. Southwestward, in the Messina Strait, the maximum height reaches 0.5 m.

Moving to the coast of Sicily, especially north of the Messina Strait, some significant features are observed. First, a strong tsunami energy concentration can be noted towards the easternmost end of Sicily, named Capo Peloro, where waves higher than 6 m hit the coastline. Along the Tyrrhenian coast of Sicily other tsunami beams are clearly visible (in green, meaning at least 1 m height), the most interesting of which is the western one, affecting the coastal area close to the small village of San Saba. Overall, the simulation confirms the generally known feature that for landslide–tsunamis the wave height decays rapidly with distance.

When comparing the simulation results with the observed effects a general underestimation can be noticed: for example, in point 3 (Cannitello) of the map of Fig. 5 a run-up of 2.9 m has been reported, higher than the 1 m maximum height obtained in the simulation. The same holds for Messina (point 8) and Reggio Calabria (point 5), the most important towns in the Messina Strait, where the difference is by far bigger. Only the area of Capo Peloro (points 6 and 7) seems to fit the observations.

Such discrepancies can be ascribed to the low resolution of the computational domain (50 m) not allowing one to describe the coastal zone and all the nonlinear effects prevailing in shallow water properly, such as the inland inundation. Further, low-resolution simulations often cannot account for tsunami wave amplification caused by resonances induced inside coastal basins, such as harbors (see, for example, Dong at al., 2010; Vela et al., 2010). The latter can explain, for example, the underestimations of the observations reported in Messina and Reggio Calabria.

These considerations push for the use of a more detailed computational grid for the sites of major interest, suggesting that the optimal approach is (1) to assess the general tsunami energy distribution in the first run of simulations carried out in a wider low-resolution domain and (2) to pick up the areas most exposed to tsunami occurrence (e.g., areas hit by energy beams) to perform higher-resolution investigations. With this approach in mind, the focus has been moved to the area enclosed in the dashed blue rectangle of Fig. 5, including Capo Peloro, where most of the devastating effects of the tsunami outside the source zone were reported (see Table 1, nos. 6 and 7).

Figure 6 Area covered by Grids 3 and 4 employed to simulate the impact of the 1783 tsunami at Capo Peloro, the easternmost point of Sicily. The magenta contour delimits the strip swept by the landslide. Coordinates are in UTM-WGS84 projection, zone 33^∘ N.

5.2 First tsunami simulation in the Capo Peloro area

An additional computational grid (Grid 3) has then been built to simulate the tsunami propagation and inundation in this zone.

Actually, in this application, due to the proximity of the landslide to the target area, the choice was to build Grid  3 by extending the high-resolution grid (Grid 1, 10 m node spacing) used by Zaniboni et al. (2016) westward for the simulation of the local tsunami effects in Scilla. The relevant data for the Sicilian region was provided by the Civil Protection Department of Sicily and retrieved from the regional cartography service, SITR (Sistema Informativo Territoriale Regionale), in the form of a DTM describing the present topography of the coastal zone in high detail. As for bathymetric data, they have been retrieved from the same sources used for Grid 2. Grid 3 covers an area of 12.5 km (east–west) by 4 km (north–south), for a total of more than 500 000 nodes. As can be noticed from the map in Fig. 6, Capo Peloro presents a morphology that is radically different from the Calabrian coasts. In Calabria, coasts are steep, rapidly ascending to 400 m a.s.l., while in the area of Capo Peloro a wide lowland is found, extending for about 2 km×1 km with elevation on the order of 1 to 5 m. This area is now densely inhabited, with many houses facing the sea, especially along the southern coast. Another interesting feature is the presence of two brackish lakes, called Lago di Ganzirri and Pantano Piccolo (also known as Faro Lake), characterized by 6 and 29 m water depth, respectively. They are connected to the sea (and to each other) by narrow channels, built under British authority at the beginning of 19th century (Leonardi et al., 2009; Manganaro et al., 2011; Ferrarin et al., 2013). One of the most relevant observations reported in the historical reconstruction of the tsunami effects is that the basin of Pantano Piccolo was reached by sea water in 1783 (Minasi, 1785). This feature should be considered in the evaluation of the simulated tsunami accuracy.

Figure 7 Propagation frames for the 1783 Scilla landslide–tsunami in the target area of Capo Peloro, eastern Sicily. Geographical coordinates are given in Fig. 6.

Figure 7 shows different frames of the tsunami propagation simulated by means of the code UBO-TSUFD, starting from the snapshot taken at 80 s. It is worth pointing out that the origin time coincides with the landslide motion initiation and recalling that, according to our simulation, the slide enters the sea between 10 and 30 s and settles at the final position around 110 s (see Zaniboni et al., 2016, for the accurate description of these results).

The landslide-generated wave heads towards Capo Peloro with a strong positive front 5 m high, directly caused by the slide entering the sea and characterized by an almost circular shape when moving in deep water. At t=100 s the wave begins changing form, due to the interaction with the platform characterizing the seabed east of Capo Peloro. When meeting shallow water, the tsunami is intensely decelerated and begins to deform. Frames at 100, 110, and 120 s display the subdivision of the wave into two fronts, north and south of the Capo Peloro end, hitting the coast obliquely. At 140 s the southern coast of the cape is reached by the tsunami, which floods the mainland by some tens of meters and then inundates the other parts. In the eastern extreme the water penetration is maximum, reaching 600–700 m distance (as visible from the 200 s frame). At the same time, the northern branch of the tsunami reaches the coastline, with a positive 2 m high front that tends to align to a direction parallel to the shoreline. The following oscillations constituting the train of waves (already visible behind the tsunami front in the first snapshots) do not produce relevant effects on the coast, apart from the cape, where a second relevant positive signal, exceeding 2 m, can be noticed at t=240 s.

Figure 8 Maximum water inundation at Capo Peloro for the high-resolution Grid 3 (western portion). The black squares mark the cells inundated in the simulation with the low-resolution Grid 2. The purple symbols represent the positions of historical buildings: Torre degli Inglesi, the site where tsunami deposits were recognized and associated with the 1783 event (Pantosti et al., 2008) and Torre Bianca and which is partially buried by sand deposits (Bottari and Carveni, 2009). For geographical coordinates, see Fig. 6.

The maximum water elevation for each node of the western portion of Grid 3 is illustrated in Fig. 8. Some of the features already observed in the 50 m grid (Grid 2; see Fig. 5) are confirmed, such as the maximum wave height concentration in the easternmost zone of Capo Peloro. Here, in correspondence with the edifice named Torre degli Inglesi (“English Tower”, referring to a building realized during British rule, marked in black in Fig. 8), some tsunami deposits were recognized and associated with the 1783 event (Pantosti et al., 2008). The simulation shows the maximum water elevation exactly for this spot, exceeding 6 m, and the maximum inland penetration, more than 700 m. Furthermore, here one can notice the main discrepancy between the low-resolution Grid 2 and the high-resolution Grid 3 simulations: in the former case, only a few inland cells (purple in Fig. 8) are inundated (corresponding to an inundation distance of 150 m), while in the latter the inland penetration is larger by far. For the remaining coastal stretch, the water flooding is limited to the scale of tens of meters or less, showing little difference between the two simulations. It is very relevant to stress that, concerning the small lake of Pantano Piccolo, both simulations exclude the possibility that it is reached by the tsunami, contrary to what reported in the historical reconstructions.

6 Capo Peloro: morphology reconsideration

The transition from Grid 2 to Grid 3 produced interesting results and precious insights on the tsunami effects in the target area. Yet, as already noted, one of the most relevant features (i.e., the inundation of the lake of Pantano Piccolo) was not reproduced, nor can the uncertainties naturally connected with the simulation of natural events like landslides and tsunamis be invoked. Indeed, to generate a tsunami able to penetrate further the flat area beyond Torre degli Inglesi and to cover the distance of more than 2 km separating this from the lake, a much larger landslide would be needed, which goes against the geological evidence. Another possibility would be that the tsunami finds a penetration route from another direction, for example from the north. In fact, the strip of land between the basin and the sea is narrow (around 150 m) and is characterized by 3 to 5 m topography, which is again really unlikely to be overrun by the simulated tsunami.

A possible alternative solution comes from the considerations of the coastal morphology of this area. In a historical study of this zone, Randazzo et al. (2014) gave evidence of the strong variability of the coastal stretch around Capo Peloro, due to the action of atmospheric agents. The dataset used to build Grid 3 refers to the present ground level that could have experienced great changes since the 1783 event. Based on the historical investigations by Bottari et al. (2006) and Bottari and Carveni (2009), we know that in Pantano Piccolo hundreds of ships took shelter during the early Roman period (5th century BCE) and that a small harbor probably also existed, meaning the presence of direct access from the sea. In the following centuries the channel was filled with sediments, but a look at the topography clearly shows that the lower ground area starts from the northeastern corner of the lake (as also evidenced by the modern 2 m isoline taken from Bottari et al., 2006) and is the best candidate for the channel that was connecting the lake to the sea more than 2000 years ago. Currently, the access is closed by a 3–4 m sand dune on the coast that can be reasonably associated with progressive sedimentation due to wind and storms and removed, in order to reconstruct a tentative 1783 morphology of the area.

Figure 9 (a) Shaded-relief map of the current topography of Capo Peloro, with a 2 m isoline (dashed red line) obtained from Bottari et al. (2006) showing the most likely connection between the lake and the sea, just east of Torre Bianca. The corrections are shown in green and blue when negative (meaning “digging” with respect to the present ground level) and in yellow and red when positive, meaning increased ground elevation. (b) Contour plot of Grid 4, with the low topography from the northeastern corner of the lake to the sea. In green, the new 2 m isoline is evidenced. With the corrected topo-bathymetry we have built, Grid 4 has the same areal coverage as Grid 3. For geographical reference, see Fig. 6.

Another element that needs to be considered is the existence, on the narrow coastal stretch separating Pantano Piccolo from the Tyrrhenian Sea, of a partially buried tower named Torre Bianca (“White Tower”, known also as “Mazzone”; see Figs. 8 and 9 for the location), which now is filled with sand up to the first floor (Bottari and Carveni, 2009). Yet, in the 18th century, this building served as a storehouse, so it was probably fully available and unburied, meaning that the previously cited assumption of removing 3–4 m of sediments (corresponding to the average height of one floor) in this area is perfectly reasonable.

Keeping all these considerations in mind, a rearrangement of Grid 3 has been done resulting in a new mesh, the Grid 4, with the same spatial extension. Fig. 9a shows a shaded-relief map of the present topography and bathymetry and the corrections that have been done to reconstruct the likely morphology of 1783. Notice again that in Fig. 9 only the western portions of Grids 3 and 4, covering the target area, are reported, since the source area remained unchanged. Such reconstruction is based on hypotheses and some evidence. For example, the two small channels now connecting Pantano Piccolo to the sea did not exist in 1783, since they were excavated during the following British rule. For this reason, they have been “filled” in Grid 4. The most relevant changes, regarding the area between the northeastern corner of Pantano Piccolo and the Torre Bianca site, is the removal of 1 to 5 m of the surface ground layer. The highest amount of “digging” is located close to the coast, where the sand dune currently closing the possible access of water (and presumably filling the first floor of Torre Bianca) has been leveled. Morphological changes in other areas of the cape are much less realistic since they would entail lowering of high topography that is less intensely affected by natural agents such as storms, winds, or sea waves.

Figure 10 Propagation frames of the 1783 landslide-generated tsunami in the western part of Grid 4. Notice the water propagation inside the lake. Geographical coordinates are given in Fig. 6.

The tsunami simulation in Grid 4 coincides with the simulation in Grid 3 for the first 200 s, as is shown in Fig. 10. The main differences start from the 220 s frame, when the coastal area east of Torre Bianca begins to be inundated. The following frames show that the water reaches more than 1 m height, channeling through the modified area and reaching the lake between 240 and 260 s. After that, a small positive wave propagating inside the basin can be noticed, 10 to 20 cm high, with some reflections. The adopted correction to the topography, then, facilitates the tsunami flooding of the Pantano Piccolo lake.

Figure 11 (a) Maximum water elevation for the western part of Grid 4. Notice that tsunami penetrates also into the lake. (b) Maximum flow depth for each topographic node. For geographical coordinates, see Fig. 6.

The map with maximum water elevation (Fig. 11a) is analogous to the one of Fig. 8. The main difference comes from the modified zone, where the water climbs up to 2.5 m close to the coast and 1 m just before entering the lake. Inside the basin, the maximum wave elevation reaches 40 cm in the western part but less in the central area, which is deeper. Figure 11b displays the height of the water column characterizing the flood, known usually as flow depth. The plot shows that the northern coast is affected by higher flow depth, reaching 2 m almost everywhere , with the maximum located at the Torre degli Inglesi (more than 5 m). The flood, just before entering Pantano, is characterized by a water column of about half a meter.

7 General discussion

In this paper we have performed a deep analysis of the 1783 tsunami, extending the results of a previous study where attention was restricted to the area of the highest impact, i.e., Scilla and surroundings. Here, we have simulated the tsunami in a broader region and we can draw the following remarks. An interesting, yet unexpected, feature that we found is the strong tsunami beam reaching the village of San Saba, located along the northern coast of Sicily, 10 km to the west (see Fig. 5). There is no historical account of a tsunami impact here. But this can be easily explained by considering that the tsunami probably only affected the beach and passed unobserved. As for the Calabrian coast, we observe that the simulated wave height rapidly decreases with distance from the tsunami origin, as expected for landslide-induced tsunami events. In addition to the positive results highlighted above, we also have to point out as a drawback that in the region covered by the 50 m resolution grid (Grid 2) a general underestimation of the observed wave height is noted, which has to be ascribed to an insufficient description of the local, shallow bathymetry and of the coastal topography features. The conclusion is that 50 m computational domain is suitable to assess the general pattern of the tsunami energy distribution, highlighting the areas where the highest waves concentrate, although this is inadequate to provide good coastal tsunami data for the 1783 case.

The area that is of the highest interest is Capo Peloro, the easternmost cape of Sicily, where the most damaging effects besides Scilla were recorded, including a deep tsunami penetration and the death of 26 people, which is the only death toll of the tsunami away from the source. Thus, this zone has been selected as a good target for a grid refinement. A higher-resolution computational domain (Grid 3, 10 m spaced) has then been assembled, based on the present morphology. The tsunami simulation on Grid 3 shows that the inundation distance in the area of maximum impact changes radically, passing from 150 m found for Grid 2 to more than 600 m (Fig. 8), in full agreement with the observations. This feature clearly shows the improvement brought by the computational grid adjustment.

However, the tsunami does not reach the Pantano Piccolo lake, which contrasts sound historical accounts (see Minasi, 1785). This misfit calls for simulation changes that can be done via three possible options: (i)  producing a stronger tsunami; (ii) using a better resolution, or (iii) changing the basic topo-bathymetry dataset, especially in the target area. We have excluded option 1, which is to use Grid 3 and to produce a stronger tsunami by using a larger volume landslide, since the landslide geometry is well constrained by geological considerations of the scar morphology, the identified offshore deposit, and the good agreement with observations of the tsunami run-up height in the Scilla beaches computed through Grid 1 (see Zaniboni et al., 2016). We have also further ruled out option (ii), which is to use grid spacing finer than 10 m, since this resolution was enough for the target area, as mentioned before, and also because it was enough to fit the 600 m observed inundation of Capo Peloro with Grid 3. So the only option left was to keep the same resolution and to build another set, Grid 4, which is based on opening a way to the tsunami access through an ancient channel that was a pathway to a well-protected inner harbor in Roman times, and which is presently filled with sediments and obstructed by a coastal sand dune (Fig. 9). The signature of this channel is still recognizable today in a linear topographic depression (only 2 m a.s.l.) going from the northeastern corner of the lake to the sea. We believe that Grid 4, allowing the tsunami to reach the lake, represents reasonably well the Capo Peloro topography at the time of the 1783 tsunami occurrence, also based on the existence of Torre Bianca building, whose present status reinforces the plausibility of sand deposit removal along the coast. The tsunami reached the lake through the channel either because there was no coastal dune, as assumed in Grid 4, or because the dune was lower than today and was overcome and/or demolished by the incoming tsunami.

8 Conclusions

The work presented here is the natural extension of the Zaniboni et al. (2016) paper, where the initial sliding mass (subaerial) was reconstructed, its descent along the Mount Pacì flank and its immersion in the sea were simulated, and the effects of the generated wave on Marina Grande beach of Scilla (where most of damages and casualties were caused) were assessed. The tsunami simulations were performed in a high-resolution (10 m) grid, covering the town of Scilla (as depicted in Fig. 4). Encouraged by the very satisfactory results obtained, which fit the detailed historical reports well, in this work we have applied the same numerical codes to extend the study of the tsunami to a wider scale.

This has been done by means of additional simulations: one in a wider domain, covering a larger portion of the Calabrian and the northeastern Sicilian coasts, characterized by a lower resolution (50 m node spacing) and two over a high-resolution domain (10 m spacing), embracing Capo Peloro, where the present topo-bathymetry (Grid 3) and a reconstructed topo-bathymetry (Grid 4) have been used to investigate tsunami flooding.

The simulation over a wider domain provides interesting hints on the tsunami energy distribution on the regional scale, heading mainly towards Capo Peloro, where tsunami effects were strong, according to historical records and studies of tsunami deposits. The focus on some areas of interest, realized with grid refinement, allows one to obtain a better fit with the observations. In this particular case, the adjustment of Capo Peloro morphology was necessary to account for the inundation of the small basin of Pantano Piccolo.

This paper together with the previous paper by Zaniboni et al. (2016) provides a good reconstruction of the 1783 landslide-induced event, respecting the main available geological constraints and fitting the tsunami observations in those places where the tsunami was most lethal, namely the areas of Scilla and Capo Peloro (see simulations in Grids 1 and 4).

Still, some level for improvement remains for fitting observations in places covered by Grid 2, which probably means that local high-resolution grids have to be nested, to handle areas where a better agreement has to be found, as in the harbors of Messina and Reggio Calabria, for example. Here local effects such as eigenmode resonance probably strengthened the tsunami waves. But this will be left to further analyses.

A final consideration regards the area of Capo Peloro, which is densely populated with a flat inland topography and in a region with many potential tsunami sources, either due to earthquakes or coastal landslides. Hence, it needs particular attention in terms of hazard vulnerability and risk assessments. This paper, through the simulation carried out on Grid 4, has shown the effect of coastal channels as elements favoring tsunami penetration. Since today there are active channels connecting the lake of Pantano Piccolo and Ganzirri to the waters of the Straits of Messina, their role should be explored and taken into account when planning and designing civil protection measures.