The Sensitivity of Tsunami Impact to Earthquake Source Parameters and Manning Friction in High-Resolution Inundation Simulations

In seismically active regions with variable dominant focal mechanisms, there is considerable tsunami inundation height uncertainty. Basic earthquake source parameters such as dip, strike, and rake affect significantly the tsunamigenic potential and the tsunami directivity. Tsunami inundation is also sensitive to other properties such as bottom friction. Despite their importance, sensitivity to these basic parameters is surprisingly sparsely studied in literature. We perform suites of systematic parameter searches to investigate the sensitivity of inundation at the towns of Catania and Siracusa on Sicily to changes both in the earthquake source parameters and the Manning friction. The inundation is modelled using the Tsunami-HySEA shallow water code on a system of nested topo-bathymetric grids with a finest spatial resolution of 10 m. This GPU-based model, with significant HPC resources, allows us to perform large numbers of high-resolution tsunami simulations. We analyze the variability of different hydrodynamic parameters due to large earthquakes with uniform slip at different locations, focal depth, and different source parameters. We consider sources both near the coastline, in which significant near-shore co-seismic deformation occurs, and offshore, where near-shore co-seismic deformation is negligible. For distant offshore earthquake sources, we see systematic and intuitive changes in the inundation with changes in strike, dip, rake, and depth. For near-shore sources, the dependency is far more complicated and co-determined by both the source mechanisms and the coastal morphology. The sensitivity studies provide directions on how to resolve the source discretization to optimize the number of sources in Probabilistic Tsunami Hazard Analysis, and they demonstrate a need for a far finer discretization of local sources than for more distant sources. For a small number of earthquake sources, we study systematically the inundation as a function of the Manning coefficient. The sensitivity of the inundation to this parameter varies greatly for different earthquake sources and topo-bathymetry at the coastline of interest. The friction greatly affects the velocities and momentum flux and to a lesser but still significant extent the inundation distance from the coastline. An understanding of all these dependencies is needed to better quantify the hazard when source complexity increases.

Introduction

Earthquake tsunamis pose a significant hazard to coastal communities, and account for approximately 80% of tsunami events globally (e.g. NCEI, 2021). A strong component of the global earthquake tsunami hazard is the one induced by large subduction zone events (e.g. Løvholt et al., 2012a; Løvholt et al., 2014; Davies et al., 2018) where sources are aligned along well known tectonic interfaces along the Pacific Ring of Fire and Sunda Trench in the Indian Ocean. On the other hand, in regions where the tectonic setting is more complex, the dominant tsunamigenic structures are less well known, and the tsunamigenesis and hazard is consequently subject to much larger uncertainty. Examples of such complex areas include for instance eastern Indonesia and the Philippines (Løvholt et al., 2012b; Horspool et al., 2014), the Caribbean (Parsons and Geist, 2008; Harbitz et al., 2012), and the Mediterranean Sea (Lorito et al., 2008; Selva et al., 2016; Basili et al., 2021; Lorito et al., 2021; Selva et al., 2021). This large source uncertainty poses a challenge for modelling the hazard, with accurate and efficient numerical modelling required for simulating tsunami generation, propagation, and inundation. Estimating the tsunami hazard that a given coastal region is exposed to requires an understanding of the likelihood and variability of virtually all possible seismic sources, the generation and propagation of the tsunami, and the coastal inundation, necessitating Probabilistic Tsunami Hazard Analysis (PTHA). A comprehensive review is provided by Grezio et al. (2017). We hence need to understand better the basic relationship between fault parameters and tsunamigenic strength to more efficiently sample the source variability and in turn correctly determine the hazard.

While there is increasing understanding that spatially variable slip on ruptures can be critical in determining the tsunami impact (Geist, 2002; McCloskey et al., 2007; Davies et al., 2015; Li et al., 2016; Murphy et al., 2016; Davies, 2019; Scala et al., 2020), it still remains essential to understand the source sensitivity from the simpler constant slip models on tsunami inundation (see also An et al., 2018). In addition to the heterogeneity of the slip, also the fault geometry may affect the tsunami modeling, as shown recently by a sensitivity analysis for subduction earthquakes in the same region of this study by Tonini et al. (2020). However, as stressed above, for many tsunamigenic earthquake scenarios, we may not have a good model of the fault geometry. Moreover, large earthquakes with complex rupture geometries are typically modelled as slip on multiple fault segments with simpler geometry. It is important to understand well the variability due to simpler uniform-slip models in its own right, also as a basis for more complex slip realizations. The model typically applied is that of Okada (1985) which computes analytically through a complex set of equations the seabed surface deformation resulting from a uniform dislocation (slip) occurring along a rectangular fault plane embedded into a homogeneous linear elastic half-space. In this approach the earthquake source, in addition to the position, is parameterised through its size (length and width of the rectangle), its orientation (strike and dip angles), and two kinematic rupture parameters (the slip and its direction in the fault plane, the rake angle). Several sensitivity studies using Okada’s source formulation have been performed. Geist (1998) reviews the basic tsunamigenic principles due to co-seismic slip. Gica et al. (2007) and Burbidge et al. (2015) perform systematic parameter sweeps to determine the sensitivity of the maximum offshore wave-height at a given location to the characteristics of the earthquake source as specified by the Okada parameters. Løvholt et al. (2012c) investigated sensitivity to slip distribution for different dip angles and depth for variable slip. However, none of the above studies performed a broad sensitivity analysis with focus on inundation; this is partly expected as this requires large computational resources.

Recent advances in efficient GPU (Graphical Processing Unit) shallow water-type tsunami models, combined with advances in High Performance Computing facilities, have now made it feasible to perform very large numbers of tsunami simulations with inundation modelled at high spatial resolution (down to a few meters). If further source down-sampling is still necessary due to limited computational resources, there are several techniques for limiting the number of simulations calculated (Lorito et al., 2015; Volpe et al., 2019; Williamson et al., 2020). Gibbons et al. (2020) performed a comprehensive PTHA for the town of Catania on the Eastern coast of Sicily comprising over 32,000 tsunami simulations for different earthquake sources in the Mediterranean Sea. The calculations were performed using the GPU-based Tsunami-HySEA code (de la Asunción et al., 2013; Macías et al., 2016, 2017). Each simulation in that study took approximately 30 min on an NVIDIA V-100 GPU, modelling the tsunami wave from source to inundation on a numerical Digital Elevation Model (DEM) with 10 m lateral resolution. However, a rather rough discretization of source parameters has been applied to reduce the computational effort, inherited from regional studies from which scenarios have been selected (see Selva et al., 2016; Basili et al., 2021). For such local studies, a finer discretization of earthquake sources, both spatially and with regards to source mechanism, can instead help improving the accuracy of the hazard analysis. To this end, a systematic sensitivity study on its impact of hazard is still required, to define efficient strategies for source description.

The accuracy of tsunami hazard estimates depends also upon the accuracy of the tsunami numerical models themselves for a fixed set of source parameters. Their accuracy will depend in turn both upon the fidelity of the DEM (e.g. Griffin et al., 2015) and on how well the resistance that surface friction, vegetation, and infrastructure provides is represented (e.g. Kaiser et al., 2011). In the numerical model used to perform the calculations described here, the combined friction effects are often described with the Manning equation where the bottom friction is proportional to the velocity squared through the Manning coefficient n. The significance of the choice of n for run-up and inundation was demonstrated by Gayer et al. (2010) who applied both spatially varying and uniform friction to tsunami simulations for Bali, Sumatra, and Java. Bricker et al. (2015) argue that values of n typically applied in tsunami simulations are too low, in particular for urban and highly vegetated areas. Here, for a limited number of earthquake scenarios, we systematically vary this parameter while keeping the earthquake source definition constant to estimate the influence of n on inundation extent and other properties of the flow. The intention is also to provide insight into the degree to which typical inundation estimates (e.g. those used in Gibbons et al., 2020) were determined by the choice of this parameter.

Here, we study inundation for the stretches of coastline surrounding the towns of Catania and Siracusa in Eastern Sicily, each represented numerically in the simulations with a system of nested grids. Based on Gibbons et al. (2020), we select a subset of the relevant parameters for which we investigate the sensitivity of model parameters on the inundation. Figure 1A displays the topo-bathymetric model of the region of the Mediterranean Sea in which the numerical tsunami simulations are performed. With the four earthquake locations in Figure 1, we can study earthquakes both offshore and with significant crustal deformation at the coast of interest. In section Numerical Simulations of Earthquake Generated Tsunamis off the Coast of Sicily we provide more comprehensive details about the sequences of simulations performed. In section Sensitivity of Inundation to the Manning Friction we discuss the sequence of calculations performed to investigate the sensitivity of inundation to the Manning friction coefficient and in section Sensitivity of Inundation to Earthquake Source Parameters we discuss the sequences of calculations performed to explore the sensitivity to earthquake source parameters. In section Conclusion and Discussion we draw conclusions.

FIGURE 1

FIGURE 1. Maps at different scales for the earthquake sources and tsunami inundation. (A) shows the fault centers for all earthquake scenarios considered. Points A, B, C, and D are all at latitude 37.15^o, with longitudes 14.70^o, 15.25^o, 16.50^o, and 17.50^o respectively. Yellow asterisks mark the locations of epicenters for scenarios where we examine the effect of the Manning friction coefficient. White circles are locations for which we conduct a fine scan of strike and dip parameters. The white square marks the epicenter of sources where we search the strike, dip, and depth parameter. Panels (B–D) show higher resolution maps of the regions surrounding the towns of Catania and Siracusa, where the tsunami inundation is to be evaluated. The black circle labelled S in panel (B) is the location 37.5368^oN, 15.1333^oE for which the sea height time-series are displayed later.

Numerical Simulations of Earthquake Generated Tsunamis off the Coast of Sicily

In this paper, we have carried out sensitivity tests by performing 1,516 separate tsunami simulations to understand how coastal inundation at high spatial resolution varies as a function of several key parameters. For these simulations, we employ the shallow water model Tsunami-HySEA (de la Asunción et al., 2013; Macías et al., 2016, 2017). The shallow water model assumes that the wave length is much longer than the water depth and hence neglects higher order characteristics such as frequency dispersion that is important for some types of tsunamis (e.g. Glimsdal et al., 2013; Løvholt et al., 2013). Manning friction terms, which assume that the bottom friction is proportional to the Manning number n and the squared current velocity, are added to the to the momentum terms in the shallow water equation system. Instantaneous earthquake deformations assuming planar uniform slip using the Okada (1985) model are provided as initial conditions to the simulations. Tsunami-HySEA solves the non-linear shallow water equations through a Finite Volume scheme by means of Riemann solvers, with shock-capturing properties that ensure that wave-breaking is incorporated and with wetting and drying processes necessary for capturing inundation. Further details of the specific Tsunami-HySEA implementation of friction are provided in de la Asunción et al. (2013). Finally, we note that Tsunami-HySEA includes a telescopic grid functionality that allow the modeller to couple coarsely resolved wave propagation over an ocean scale with high-resolution inundation simulations in a local domain. The simulations are performed with a system of four telescopic nested grids. The coarsest grid has the spatial extent displayed in Figure 1A and a spatial resolution of 640 m. The others have progressively finer resolutions of 160, 40, and 10 m. The inundation is modelled for the regions surrounding the towns of Catania and Siracusa on the East coast of Sicily (see also Tonini et al., 2021).

All tsunamis simulated are generated by hypothetical earthquakes whose epicenters are located at the positions A, B, C, and D in Figure 1A. These earthquake epicenters are based on a subset of the earthquake sources considered in the probabilistic earthquake model of the TSUMAPS-NEAM (NEAMTHM18) Tsunami Hazard Model (Basili et al., 2021). The NEAMTHM18 model consists of a discretization of tsunamigenic earthquake scenarios covering the entire NEAM (North East Atlantic and Mediterranean) region with different magnitudes, depths, fault orientations, and slip distributions. Here, we identified two crustal earthquake scenarios from the NEAMTHM18 model which had resulted in the greatest inundation in the bay of Catania in the tsunami simulations from Gibbons et al. (2020). Both were associated with hypothetical Mw 8.0 earthquakes with fault ruptures of length 183.5 km and width 38.1 km, and with slip of 7.5 m. Both had epicenters on the coast (point B in Figure 1), or inland (close to point A in Figure 1), and both earthquakes had resulted in significant crustal deformation both onshore and offshore. The coseismic deformations resulting from these two earthquakes are displayed in Figure 2 (panels A and C). These scenarios are in the low-probability and high-consequence end of the source discretization.

FIGURE 2

FIGURE 2. Sea floor deformation using the Tsunami-HySEA code with Okada sources with strike, dip, rake, and depth parameters as indicated. All sources have length 183.5 km, width 38.1 km, and slip 7.5 m. Red and blue indicate a raising and a sinking of the ocean floor respectively. The Depth quoted here is the depth of the fault plane centroid. The sea floor displacements in panels (A,C) correspond to scenarios based on definitions in the NEAMTHM18 model (subsequently labelled M1 and M2 respectively). The remaining panels (B,D–F) show non-NEAMTHM18 scenario definitions considered in this study.

Figure 1B displays the locations of the numerical grids with the highest spatial resolution (10 m) and Figures 1C,D detail the high-resolution topo-bathymetry with a colour scale chosen to emphasize spatial differences in elevation in the zero to 20-m range of greatest relevance for tsunami inundation and consequence for human life and infrastructure. Regions with elevation over sea-level in the range 0–5 m are coloured green, regions with elevation in the range 5–10 m are coloured brown, and regions with elevation exceeding 10 m coloured red with the colour fading as rising topography takes land and infrastructure above the tsunami threat. Details of the Digital Elevation Model (DEM) are provided in Gibbons et al. (2020) and the sources are listed in the Data Availability Statement. The DEM does not include buildings but, even at a 10-m resolution, some choices have been made regarding built infrastructure with, for example, raised road constructions visible in Figure 1. In the DEM and Tsunami-HySEA inundation modelling, these provide a barrier to the incoming water.

The first part of the analysis in the present study is dedicated to exploring systematically how the inundation varies with the choice of the Manning friction parameter. To this end, we perform a suite of tsunami simulations considering four different earthquake sources and 19 values of n (76 scenarios in total). A constant value of the Manning n value for the entire domain is assumed for simplicity. Such an approach is common practice in tsunami simulations even though this does not resolve spatially dependent friction properties of the Earth’s surface. We justify it here on the grounds that we want to examine inundation over varying topography as a function of this parameter. A spatially varying n here would represent a vast number of parametric choices and would complicate interpretation. Results from our uniform n parameter searches will hopefully provide guidance in experiment design for studies with spatially varying n. The second part of the study explores the dependency of the tsunami inundation on the variability of fault orientation (i.e. strike and dip angles), the slip direction (i.e. the rake angle), and the depth of the source; in that case we analyse a suite of 1,440 tsunami simulations. Figure 2 shows the sea floor deformation resulting from six different specifications of earthquake parameters selected from the total set of combinations. We adopt for each scenario a uniform slip value along the fault plane and all simulations were run for a total of 4 h following the earthquake.

Sensitivity of Inundation to the Manning Friction

Testing the sensitivity of the inundation to the Manning friction was performed for four hypothetical earthquake scenarios, M1, M2, M3, and M4. Two of them, labelled M1 and M2, are almost identical to scenarios from the NEAMTHM18 model, changed only in the exact coordinates of the earthquake fault centers. Scenarios M3 and M4 have the same Okada parameters as scenarios M1 and M2 respectively, except for the epicenters that are moved significantly offshore to point D to model the sensitivity of a more distal earthquake. These scenarios were selected since the corresponding models from NEAMTHM18 had demonstrated inundation at significant distances inland from the shoreline, such that changes in the degree of inundation due to changes in the friction would be easy to visualize. Out of the subset of scenarios simulated in Gibbons et al. (2020) which generated high inundation, the selection for the friction study here is otherwise quite arbitrary. The seafloor displacements for all four scenarios are displayed in Figure 3 together with the areas of inundations at Catania and Siracusa resulting from the generated tsunamis.

FIGURE 3

FIGURE 3. Cumulative inundation area plots for Catania and Siracusa for the scenarios M1, M2, M3, and M4. For a given value of h, the curve indicates the area of land that was dry prior to the event which has experienced a maximum flow depth exceeding h m. The 19 curves in each panel cover values of the Manning friction, n, between 0.05 and 0.95 in intervals of 0.005. The curves for n = 0.03 are highlighted and coloured black to mark that this was the value used for all calculations in Gibbons et al. (2020).

Inundation scenarios are obtained for Catania and Siracusa for each of the four earthquakes considered and each of the 19 values of n, covering the range from 0.005 to 0.095 in increments of 0.005. The value n = 0.03 (which is widely used in operational tsunami modelling) is highlighted in Figure 3 as this was the value applied uniformly to all calculations in Gibbons et al. (2020). In all cases, a constant value of n for the entire domain is assumed for simplicity. The Tsunami-HySEA code calculates the tsunami propagation on the bathymetry and topography as deformed by the earthquake and records, for each grid location, the highest water level (relative to the undisturbed sea surface) attained over the entire simulation. The flow depth for a given location is obtained in post-processing by subtracting the deformed topography from this value. The flow-depth is therefore sensitive to the onshore co-seismic displacement which varies from scenario to scenario, in addition to the water displacement offshore. We consider only areas that were above sea-level prior to the earthquake. Figure 3 displays, for a given height, h, the area of land (in square kilometers) for which the maximum flow depth exceeds h. We note that both the horizontal and vertical scales vary from panel to panel in Figure 3. Each set of panels corresponds to a different earthquake scenario, and the curves for that scenario are reported for all the 19 n coefficient values. We compare directly only the shapes of the curves for the different values of n. All scales were chosen to optimize the resolution of the relationships between the different curves, even though this is at the expense of direct comparison between curves in different panels.

The space between the curves in the y-direction gives a measure of the sensitivity of inundation to the n parameter. If the lines are very close together, it means that the inundated area is not highly sensitive to the parameter for the given flow depth. If the lines are far apart, it means that the modelled inundation varies greatly for different values of the friction. The ratio between the value of a curve for a given h and the value of the same curve at h = 0 indicates the proportion of the inundated land experiencing the indicated flow depth or higher. That the curves diminish towards zero indicates that relatively small regions experience the maximum flow depths calculated for each scenario. For smaller flow depths, the lines diverge. Simulations with very low values of n model inundation that covers far wider regions than simulations with much higher values.

The results are evaluated for both Catania and Siracusa. We see, in Figure 1, three fundamental differences between the regions exposed to inundation at Catania and Siracusa. First, that there is far less low-lying land at Siracusa than at Catania. Second, that the Bay of Catania has almost 20 km of exposed coastline compared with an inlet less than 2 km wide at Siracusa. Thirdly, much of the Bay of Catania has several km of relatively shallow water (less than 100 m deep) separating the coast from the deep sea. At Siracusa, a large shallow bay is partially protected by a large sea wall and the sea floor falls to great depth within a far shorter distance than at Catania.

For the Bay of Catania, the inundation area vs n curves are similar for scenarios M1 and M2. While the epicenters and angle of strike for the two cases differ somewhat, both result in comparable co-seismic subsidence in the region where the inundation is measured. The curves for offshore earthquake scenarios M3 and M4 also display somewhat similar shapes, but the inundation is several times greater for scenario M4 than for scenario M3. Scenarios M3 and M4 differ in the choice of both the strike and dip angles. The strike angle for scenario M4 (157^o) makes the rupture more parallel to the coastline than is the case for scenario M3 (strike = 67^o), and the angle of dip for case M4 (50^o) is greater than for case M3 (30^o). Both these factors may contribute to the higher inundation for scenario M4. Scenarios M1 (strike = 67^o) and M2 (strike = 157^o) differ similarly in orientation to scenarios M3 and M4. However, being far closer to shore and with significant co-seismic displacement within the coastal zone, the tsunami generation relative to the shore is more complicated than for the more distant sources. There are many factors affecting the inundation area metrics for scenarios M1 and M2 and the similarity in the inundated extents may be quite coincidental.

For Siracusa, the two scenarios M2 and M4 (both with strike = 157^o) result in significantly greater inundation than either scenarios M1 or M3 (both with strike = 67^o). That the inundation for scenario M4 significantly exceeds the inundation for scenario M3 at Siracusa is, as for Catania, likely the result of the strike angle being more aligned with the coastline, resulting in an unfavourable directivity of the tsunami. The corresponding difference between scenarios M1 and M2 is likely exacerbated by the significant co-seismic subsidence at Siracusa for scenario M2 which does not happen in scenario M1 (see uppermost panel of Figure 3).

Figure 4 illustrates what the inundation curves mean in terms of the spatial pattern of the inundation for scenario M1 at Catania (as shown in the top left panel of Figure 3). We focus on the pre-earthquake topography between zero and 10 m with a greyscale indication of elevation. A sequence of coloured lines marks the maximum extent of inundation for each of the specified levels of flow depth as a function of the Manning friction. The Southernmost part of the Bay of Catania is low-lying land bounded by a steep rise in topography. Panel A shows that the inundation of the low-lying region is stopped only by topography for all but the highest (and least physical) values of n. In the Northern part of the bay, the topography rises somewhat more gradually. The extent of the inundation increases steadily further inland for lower values of n. In this northern part, the extent of inundation with the high friction parameter of n = 0.09 is approximately 50% of the extent of inundation with the respective low value of n = 0.01 with intermediate (and more realistic) values resulting in inundation extents between these extremes. Greater flow depths (panels B through F) are increasingly restricted to the near-shore region. In the more sloping northern part of the bay, the limit of inland reach approaches the coast steadily for each increase in flow depth considered. In the far flatter southern part of the bay, the area of inundation is less sensitive to n. For flow depths of 0–2 m (panels A–C, only the very highest values of n prevent almost total inundation of the low-lying region. For a flow depth of 5 m (panel E), only at the lowest values of n is there inundation significantly inland. As expected from the inundation curves in Figure 3, a flow depth of 10 m is recorded only close to the shoreline to the north of the harbour wall (panel F).

FIGURE 4

FIGURE 4. Inundation regions for Catania given earthquake scenario M1 for flow depth as indicated for different values of the Manning friction, n. We note that n applies both in onshore and offshore regions and that high n in very shallow offshore water may dissipate some energy in the tsunami before it reaches the shore. Panels (A–F) show respectively the maximum extent of inundation for a flow depth of 0, 1, 1.0, 2.0, 3.0, 5.0, and 10.0 m.

Figure 5 displays the maximum momentum flux attained as a function of location for each of the scenarios and for three chosen values of n: 0.01, 0.03, and 0.05. This measure of the tsunami impact changes more significantly with n than the spatial extent of the inundation. The wave period of the incident wave is in the order of 10 min; this gives rise to a positive wave forcing on the coastline that lasts for several minutes. This will feed the inland wave over a long time before the wave can retreat. In the lower lying regions, the inundation will continue until stopped by an ultimate topographic barrier. The flow depth, and consequently the momentum flux, will decrease more with distance from the sea for higher values of n. For areas with gently increasing topography, the reduction in flow depth with distance resulting from the higher n will lead to a corresponding reduction in the inundation distance.

FIGURE 5

FIGURE 5. Maximum momentum flux recorded as a function of position in the Bay of Catania for earthquake scenarios M1, M2, M3, and M4 in panels (A–D) respectively. Each panel shows the results of the simulations with Manning friction values 0.01, 0.03, and 0.05 as indicated.

Figure 6 displays corresponding patterns of maximum momentum flux but for the region surrounding Siracusa. Whereas the inundation in the Bay of Catania is very similar for scenarios M1 and M2 (Figures 5A,B), the inundation at Siracusa for scenario M1 (Figure 6A) is significantly lower than that for scenario M2 (Figure 6B). This is easy to understand in relation to the fault specifications (Figure 3). Both scenarios have significant downward displacement along the Catania coastline but only scenario M2 results in ground displacement directly at Siracusa. Even though the inundation for scenario M1 is significantly less at Siracusa, there are clear differences in both the extent of inundation and in the maximum momentum flux for the different values of n. In all cases, the inundation is greatest for the town of Siracusa (at the northern end of the bay) although inundation penetrates several km inland along the lowest-lying ground both for scenarios M2 and M4.

FIGURE 6

FIGURE 6. Maximum momentum flux recorded as a function of position around Siracusa for earthquake scenarios M1, M2, M3, and M4 in panels (A–D) respectively. Each panel shows the results of the simulations with Manning friction values 0.01, 0.03, and 0.05 as indicated.

Sensitivity of Inundation to Earthquake Source Parameters

We perform four different sweeps of parameters, with one being described in each of the following four sections. The first of these sweeps considers sensitivity to the angle of rake. The second and third sweeps scan jointly the parameters strike and dip, considering separately the far offshore and the near-shore earthquake sources. The fourth and final sweep investigates the dependency on strike, dip, and depth.

Sensitivity to the Angle of Rake

The angle of rake gives the relative direction in which the hanging wall moves over the foot wall. A rake angle of zero or 180° is a strike-slip fault, in which there is no vertical slip. A rake angle of +90^o is a reverse thrusting fault, where the hanging wall moves upwards relative to the foot wall. A rake angle of 270° or −90° is a normal fault in which the hanging wall moves downwards relative to the foot wall. All other angles of rake indicate some form of oblique faulting. The eastern coast of Sicily is approximately North-South aligned such that a tsunami generated by an earthquake to the East of Sicily will likely have greatest impact if the rupture has a North-South strike angle. We also anticipate that an dip angle of around 45–50°, accompanied by dip-slip faulting (rake around 90° or −90°), will likely have the greatest effect (Kajiura, 1980; Geist, 1998). We ran 48 simulations, with strike angles of both 0° and 180° and the full range of rake angles in increments of 15^o. The epicenter for all the earthquake sources selected was set to the point labelled D in Figure 1; this ensured that no co-seismic displacement occurred on or close to land at the coastline of interest.

In addition to the inundation maps, we also save time-series of tsunami surface elevations at several offshore points of interest from the numerical tsunami simulations. Figure 7 displays these time-series for each of the 48 simulations, evaluated at a single location offshore Catania: synthetic marigrams. Figure 7B is for strike = 0^o, where the fault dips to the East and away from Sicily, and panel (B) is for strike = 180^o, where the fault dips to the West and towards Sicily. As expected, the greatest amplitudes are observed for the pure reverse faulting (rake = 90^o) and the pure normal faulting (rake = −90^o) whereas the strike-slip faulting results in very small wave-heights. As is quite well established, the oblique faulting scenarios result in very similar waveforms to the purely reverse and normal faulting scenarios, just at lower amplitudes, as the faulting style controls the proportion between horizontal and vertical coseismic sea floor displacement; dip-slip faulting provokes a larger amount of the more tsunamigenic vertical displacement.

FIGURE 7

FIGURE 7. Varying the angle of rake. (A,B) Offshore surface elevation time-series generated by the Tsunami-HySEA code for the location 37.5368^oN, 15.1333^oE, offshore of the Bay of Catania on Sicily as a function of minutes after the earthquake origin time. All channels are displayed to the same vertical scale with the red scale bar indicating a vertical displacement of 4 m (C,D) Area of land inundated with a flow depth exceeding the value indicated by the colour for the tsunamis generated by the earthquakes in the scenarios as indicated. The areas measured for Catania are within the region displayed in Figure 1C and the areas measured for Siracusa are within the region displayed in Figure 1D. The curves are displayed for the following values of the height, h: 0.01, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, and 4.0 m.

The choice between strike 0^o and strike 180^o has a second order effect on the amplitudes, affecting the shape of the waves more than their amplitude. For the strike = 0^o simulations (fault dipping away from Sicily: Figure 7A), reverse faulting (rake > 0^o) results in a small initial drop in the offshore sea elevation, around 15 min after the earthquake, followed by a rapid increase, for which the global maximum is achieved about 20 min after the origin time. For the normal faulting earthquakes (rake < 0^o) the opposite occurs; a small increase is followed by a significant decrease with maximum inundation happening some 20 min after the initial wave arrival. For the strike = 180^o simulations (fault dipping towards Sicily: Figure 7B), with reverse faulting, the maximum sea elevation still happens on the first upward wave but with a somewhat slower rise. For the normal faulting sources, the offshore elevation is somewhat lower than both the reverse faulting sources and the normal faulting sources for the strike = 0^o case.

Figures 7C,D shows the dependence of the inundation features to the faulting parameters. We choose several discrete target flow depth values, h, and for each we plot, as a function of the rake angle, a curve of the area of land inundated with maximum flow depth greater than h. This is done separately for Catania and Siracusa, and for strike = 0^o and strike = 180^o, as indicated. For both Catania and Siracusa, and both strike = 0^o and strike = 180^o, there is an approximate order of magnitude scaling factor between the inundated area given a normal or reverse thrust fault movement (rake = 90^o or -90^o) and the inundated area given a purely strike-slip fault. For Catania, the reverse thrusting faults generate a 20–25% greater inundation than the corresponding normal faults, which reflects the ratio between the positive and negative offshore amplitudes (Figures 7A,B). For Siracusa, there is less of a difference. As expected there is a quite simple relation between the rake angle and the area of inundation, supporting the assumption that we can limit our parameter searches in PTHA earthquake scenarios to the cases 180^o, 90^o, 0^o, and −90^o (as with e.g. Selva et al., 2016; Basili et al., 2021) and estimate from calibration curves what the inundation would be for oblique faulting. We appreciate that there may be situations in more complicated topo-bathymetry where the angle of rake and severity of inundation could have a more complicated relationship.

Sensitivity to the Fault Strike and Dip (Far Offshore Case)

For the second sweep of earthquake parameters, we stay at the further offshore source location (Figure 1D) and vary the angles of strike and dip while assuming rake angles of +90^o or −90^o. We cover the full range of strike angles in intervals of 10^o and consider dip angles in the range 10^o–70^o, also in intervals of 10^o. Figure 8 displays the inundated area in the Bay of Catania (panels A and C) and at Siracusa (panels B and D) as a function of strike and dip using intensity of colour in a wheel; the sources with shallow dip angles (dip = 10^o) are displayed in the outer ring and sources with increasingly steep dip angles displayed progressively closer to the center. The upper panels (A and B) are for reverse thrust faulting earthquakes and the lower panels (C and D) are for normal earthquake sources.

FIGURE 8

FIGURE 8. Area of tsunami inundation in km² for Catania (A,C) and Siracusa (B,D) for earthquake scenarios with fault centers at point D (significantly offshore). Panels (A,B) display inundation for rake = 90^o (reverse faulting) and panels (C,D) for rake = −90^o (normal faulting). The colour scales differ for the two locations, reflecting the areas subject to inundation, but are identical for the normal and reverse sources at both locations.

Strike and dip seem to influence the inundation area almost independently, with the peaks induced by the dip at intermediate angles that modulate the trend imposed by the strike angle. The strongest dependence is with the strike angle with by far the greatest inundation for sources with strike angles in an approximately North-South direction. The strike angles resulting in the most significant inundation differ slightly for the two coastal stretches, matching almost perfectly the corresponding coastline orientations at the two locations (see Figure 1B). A strike angle differing from this direction by 30° or more results in a significant decrease in the inundation area.

The dip angle dependence is weaker than the strike angle dependence. Maximum inundation typically occurs for a dip angle around 45^o degrees although there are exceptions. The steeper dip angle sources result in greater inundation for some strike angles (at Siracusa in particular) and only for the shallowest dip angles (e.g. 10^o) is the inundation extent significantly lower than for the steeper angles. There is significant symmetry between the inundation areas resulting from the reverse fault earthquakes (Figures 8A,B) and the corresponding normal counterparts (Figures 8C,D). Inundation varies relatively smoothly with both strike and dip and it appears that the chosen discretization is sufficiently fine that no qualitatively different inundation patterns have been omitted. Each of the rose diagrams in Figure 8 contains the results from 252 calculations such that, with this sampling density, 504 scenarios are required to cover the strike, dip, and rake parameters alone for a single location, depth, dimensions, and slip specifications. In a PTHA, where many fault centers, magnitudes, and slip specifications need to be covered, we would ideally like to reduce the sampling density for these angles. Given the relatively small range of the strike angles resulting in the greatest inundation, it seems necessary to have 30-degree intervals or finer for the angle of strike for resolving this parameter space in hazard analysis. The smaller sensitivity to the dip angle would suggest that the sampling interval in dip could be reduced without fear of introducing significant bias in inundation quantification.

Figures 9A,C show maximum flow depth maps at both Catania and Siracusa for four of the scenarios displayed in Figure 8 resulting in the greatest inundation areas. All scenarios have a dip angle of 40^o and we consider strike = 0^o and 180^o for both normal and reverse earthquake mechanisms. The patterns of maximum inundation are relatively consistent for all scenarios with no qualitative differences observed. Panels B and D shows synthetic marigrams at the same offshore point as displayed in Figure 7 for each of the calculations with dip = 50^o. Whereas the wave arrival time was almost identical for all the traces in Figure 7, with a fixed strike angle and varying rake, there is an almost sinusoidally varying arrival time for the traces in Figure 9. The East-West striking sources result in weaker inundation but an earlier arriving wave, as part of the seafloor displacement occurs closer to the sensor. The maximum surface elevations achieved for strike = 0^o and strike = 180^o are similar although there are subtle differences in the shape of the waveform. The gradual change in the shape of the wavelet with the 10^o increment in the strike angle is again a good indication that the applied parameter sampling density is sufficient for a qualitative understanding of the behaviour over the full range of strike angles. We note that these single-point offshore time-series do not provide any information about the directionality of the wave. The offshore wave amplitudes are lower for the East-West striking sources than for the North-South striking sources and the oblique directivity of the incoming wave will likely reduce the inundation further.

FIGURE 9

FIGURE 9. Inundation maps for Catania and Siracusa for the relatively distant offshore earthquakes as indicated together with surface elevation time-series as simulated at the location 37.5368^oN, 15.1333^oE. Panels (A,B) display inundation maps and time-series for rake = 90^o, panels (C,D) for rake = −90^o. The time-series are all for the dip = 50^o scenarios. Vertical surface elevation scales and colour scales are identical for all panels.

The simulated offshore marigrams in Figures 9B,D indicate the most significant differences between the normal and reverse faulting sources. For the reverse mechanisms (Figure 9B), the maximum surface elevation is attained on the first peak; the corresponding times in the traces for the normal mechanism earthquakes (Figure 9D) has a negative peak with similar amplitude. Along the wavetrain, there is considerable symmetry between each trace for the reverse thrusting earthquake in Figure 9B and the negative of the corresponding trace for the corresponding normal source in Figure 9D. As we observed for the normal faulting sources in Figure 7, the first major peak in amplitude is downwards and the maximum upward offshore wave-height arrives some 20 min later.

Sensitivity to the Fault Strike and Dip (Near-Shore Case)

We repeat an identical set of calculations to those in the previous section, with the significant difference that all sources here have an epicenter at the point B: the one closest to the shore of interest. Figure 10 displays the exact analogies of the quantities displayed in the panels of Figure 8 but for the near-shore earthquakes. Almost every symmetrical feature in Figure 8 is replaced by an asymmetric (or antisymmetric) feature in Figure 10. For the significantly offshore earthquake sources (Figure 8), there is a relatively simple relationship between the inundation area and the strike angle, for a given dip angle, and similarly between the inundation area and the dip, for a given strike angle. (The functions of strike have maxima near 0^o and 180^o and minima near 90^o and 270^o; the functions of dip have maxima close to 45^o.) For the close-to-shore earthquakes, there are also trends indicating significant dependencies of inundation on the specific combination of strike and dip.

FIGURE 10

FIGURE 10. Area of tsunami inundation in km² for Catania and Siracusa for the earthquake scenarios with fault centers at point B (near-shore). Panels (A,B) display inundation for rake = 90^o, panels (C,D) for rake = −90^o. The colour scales differ for the two locations, reflecting the areas subject to inundation, but are identical for the normal and reverse sources at both locations.

The patterns of inundation area for the distal reverse faulting sources (Figures 8A,B) are almost identical to the corresponding patterns for the distal normal faulting sources (Figures 8C,D). All panels of Figure 10 show complicated patterns of dependence upon the strike and dip angles with non-linear dependencies on the two different parameters. The inundation rose plots for the near-shore reverse faulting sources (Figures 10A,B) are almost inverse images of the corresponding plots for the near-shore normal faulting sources (Figures 10C,D). The relationship between the land orientation and the source orientation varies greatly from case to case for the near-shore earthquakes, with the exact spatial form of the uplift or subsidence depending strongly on both angles of strike and slip. From the present analysis, we cannot determine the causative reason for this parametric control on the inundated area, but we note the following aspects that can provide a substantial influence: First, the rake angle controls the displaced volume of water and the depth of generation, this provides two different proxies for the tsunamigenic strength. The displaced volume and distribution of the initial wave naturally influences this strength and the depth influences the relative contribution of shoaling. Second, the source orientation influences the relative orientation towards the shoreline, influencing the wave directivity. Third, as the angle of rake determines either upward or downward displacement, the complementary nature of panels (A, B) and (C, D) in Figure 10 indicate that the co-seismic displacement on land can influence the spatial extent of the tsunami inundation for these sources. These trends are difficult to predict and to be generalized. Therefore, for a PTHA, a rather fine discretization of the parameters seems to be required for near-field sources.

Figure 11 shows corresponding inundation areas for Catania and Siracusa for four selected earthquake scenarios. All these sources have strike angle 20^o but differ in dip and rake angles and show very different inundation patterns on both stretches of coastline. Figure 11A shows inundation maps for two reverse thrusting earthquakes with dip angles 10^o and 70^o. The shallow dip reverse fault earthquake results in moderate inundation at both Catania and Siracusa whereas the steep dip earthquake results in severe inundation at Catania and exceptionally little at Siracusa. The maps in Figure 11A provide an interpretation of the relative colour scales in Figures 10A,B. For the normal fault earthquakes with the same dip angles (Figure 11C), the situation is the exact reverse. There is severe inundation at Catania for the shallow dip source and less for the steep dip source. The opposite is also seen at Siracusa; there is far greater inundation for the steep dip earthquake than for the shallow dip source.

FIGURE 11

FIGURE 11. Inundation maps for Catania and Siracusa for the relatively distant offshore earthquakes as indicated together with surface elevation time-series as simulated at the location 37.5368^oN, 15.1333^oE. Panels (A,B) display inundation maps and wave height for rake = 90^o, panels (C,D) for rake = −90^o. The time-series are all for the dip = 50^o scenarios. Vertical surface elevation scales and colour scales are identical for all panels.

The surface elevation time-series (Figures 11B,D) are also distinctly different from the corresponding traces in Figure 9. For the offshore sources (Figure 9), there is a delay of between 10 and 20 min between the time of the earthquake and the water wave reaching the sensor; the exact time of arrival varies near-sinusoidally with the angle of strike. For the near-shore sources (Figure 11), the displacement of water at the location of the sensor is immediate or almost immediate after the earthquake.

Sensitivity to the Earthquake Depth

Our final parameter sweep explores the sensitivity of the inundation to the depth of the earthquake, in addition to the angles of strike and dip. The fault centers of all earthquakes are located at point C in Figure 1. As it is primarily depth we are focusing on, we present a coarser discretization of the angles of strike and dip. The earthquakes are located such that little or no co-seismic displacement is to be expected at the coast. However, the distance between the extremes of the fault plane and the coast becomes relatively small for values of the strike close to 90 or 270°. Figure 12 displays the inundation rose plots for the bay of Catania for the reverse thrusting scenarios at six different depths. The depth displayed is that of the top of the fault, such that one adds half the product of the fault width and the sine of the angle of dip to get the depth of the fault center.

FIGURE 12

FIGURE 12. Area of inundation in the bay of Catania as a function of strike, dip, and depth for tsunamis with the earthquake sources with rake = 90^o. Panels (A–F) indicate respectively inundation resulting from earthquakes with depths 0, 1, 5, 10, 15, and 20 km.

Figures 12A–C, for the shallowest events, all display quite similar rose plots. The inundation is strongest for the strike = 0^o and strike = 180^o scenarios, although the inundation for the 45^o and 225^o scenarios is more significant than for the more distant earthquakes displayed in Figure 10. The dependence upon the dip angle is quite strong with significant inundation for dip = 50^o and very little inundation for dip = 10^o. However, this dependency is significantly influenced by the depth of the source; the deeper the source, the weaker the dip-dependence becomes. Indeed, Figures 12D–F show the dip and strike dependence for progressively deeper earthquakes and, while the strike-dependence remains quite constant, the variability with the angle of dip diminishes. For the deepest earthquakes (panel F), the dependence is almost purely a function of strike alone.

This is interpreted as an effect of the deeper earthquakes providing more smoothing of the initial seabed displacement, which reduces the tsunami variability. Figure 13 shows the inundation area at Siracusa for the same set of scenarios displayed in Figure 12 for Catania. The severity of tsunami inundation varies with strike angle somewhat differently for the two locations. However, the first order observations are consistent; the inundation is most significant for North-South striking earthquake sources and the dependence upon the angle of dip diminishes with greater depth. For PTHA, the trend is quite smooth. Thus, the consideration of different depth is very important, but a rough discretization may be sufficient to cover the natural variability of tsunami inundation.

FIGURE 13

FIGURE 13. Area of inundation in the bay of Siracusa as a function of strike, dip, and depth for tsunamis with the earthquake sources with rake = 90^o. Note that the colour scales in Figures 12, 13 are different with lower areas of inundation for the smaller bay at Siracusa. Panels (A–F) indicate respectively inundation resulting from earthquakes with depths 0, 1, 5, 10, 15, and 20 km.

Figure 14 displays offshore surface elevation time-series with the panels (A), (B), (C), and (D) showing the time-series as a function of depths for the angles of dip for the panels as indicated. For the sources closest to the surface, the surface elevation is very sensitive to the angle of dip with the uppermost traces in panels (A), (B), (C), and (D) increasing in amplitude. The lowest amplitude wave-heights for the dip = 10^o case is consistent with the minimum area of inundation displayed in Figures 12A, 13A. Like with the inundation areas in Figures 12, 13, the surface elevation time-series displayed in Figure 14 become more and more similar for the different angles of dip as the depth decreases.

FIGURE 14

FIGURE 14. Offshore time-series of tsunami surface elevation as a function of dip angle and depth for tsunamis from earthquake sources with rake = 90^oand strike = 0^o. The scenario with the top of the fault at 1 km depth is displayed in red and was included first and foremost as a check that there was no fundamental or qualitative difference between having the top of the fault plane at zero or close to zero. Panels (A–D) show the results using dip angles of 10°, 30°, 50°, and 70° respectively.

The maximum flow depths recorded in the simulations for which the time-series are displayed in Figure 14D (uppermost and lowermost traces) are mapped out in Figure 15. The shallow source results in a higher amplitude wavetrain with shorter period waves than the deeper source. The inundation maps for the Bay of Catania (Figures 15A,B) for the two scenarios are very similar, the greatest differences being observed to the North where the near-shore water is deepest. The difference in the maximum flow depth for Siracusa is much greater with far higher flow depths recorded for the shallow source than the deep source.

FIGURE 15

FIGURE 15. Maximum flow-depth in meters recorded for simulations with strike = 0^o, rake = 90^o, and dip = 70^o with depth 0 km (panels (A,C)) and depth 20 km (panels (B,D)).

Conclusion and Discussion

We have performed multiple numerical simulations modelling earthquake-generated tsunamis in the Mediterranean Sea, where the inundation at the towns of Catania and Siracusa on the island of Sicily is estimated on 10 m by 10 m resolution grids. The non-linear shallow water Tsunami-HySEA model was used with a system of nested grids, with a scale of four increase in cell dimensions for each level. All earthquake sources were modelled using uniform slip on a single rectangular fault segment, modelling large earthquakes centered at four different locations. We performed parameter sweeps studying how the inundation at Catania and Siracusa varied with the Manning friction, n, and different combinations of the geometrical Okada fault parameters. We have performed both an analysis of the dependence of the inundation characteristic to the input parameter variation (seismic parameters) and to the model parameterization (Manning coefficient), but we also briefly discussed the sensitivity to the adopted discretization, that is to say how densely the allowed parameter ranges are sampled. Just as Burbidge et al. (2015) demonstrated that offshore tsunami height is a complex function of earthquake parameters, we demonstrate significant differences in the intensity and extent of onshore inundation both as a function of source and friction parameters.

We confirm that the Manning friction, n, is a significant parameter in determining the extent of inundation in numerical tsunami simulations. The sensitivity of the inundation to the choice of this parameter is dependent upon the amplitude of the incoming wave and on the topography of the coastal region. In our study of the inundation of the Bay of Catania in Sicily, we showed that the sensitivity of the inundation area to the Manning coefficient is likely greatest when the topography increases gently. While increasing n reduces the geographical extent of the inundation, the reduction in the momentum flux is likely far larger. This metric may have greater significance as a hazard indicator than the tsunami flow depth. When PTHA is focused on quantifying inundation area, the results of our test show that there is a significant dependency of the results on the selected Manning coefficient. This means that more attention should be devoted in future local PTHA analyses to constraining the arising epistemic uncertainty, which may significantly impact the overall results.

For PTHA, we need a meaningful discretization of all potential earthquake tsunami sources. To make a PTHA computationally feasible, and calculable in a time-frame of interest, we want the minimum number of scenario definitions possible, while adequately covering the expected aleatory variability and related parameter space, and eventually different parameterizations or ranges stemming from the description of the epistemic uncertainty. In this study, we have considered only uniform slip on faults of a single dimension and only a single value of the slip. In addition, we have greatly restricted the number of earthquake locations. These limitations have allowed us to cover the angles of strike, dip, and rake using a fine discretization. Knowing how finely the source parameter space needs to be discretized is not possible, or is highly subjective, without a comprehensive sensitivity study. Or, seen the other way around, the risk is to adopt a too fine discretization to avoid undersampling. We wish for a discretization of each parameter such that no scenarios with significantly different tsunami impact are missed. In other words, the tsunami impact for a scenario with a given value of one parameter can be reasonably well predicted from the tsunami impact resulting from scenarios with neighbouring values of that parameter. Our goal is that the behaviour at arbitrarily specified points in our parameter space can be well interpolated from the behaviour at the points in our discretization.

We see examples of relatively smoothly varying behaviour as a function of the angle of rake (Figure 7) and the angles of strike and dip (Figure 8) for significantly offshore earthquakes. The continuity of the inundation metrics and the shapes of the offshore time-series with changes in these parameters give us confidence that the number of simulations performed can be limited to relatively few points in this parameter space for a given fault size and epicenter. In practice, this means for offshore earthquakes that we can probably consider only rake = 90^o and rake = −90^o, angles of strike in 30-degree increments, and 3 or 4 different values of the dip angle. For scenarios similar to those studied here, the tsunami impact from earthquakes with intermediate parameters can likely be well estimated using a relatively sparse set of scenarios. When the earthquake sources are far closer to the shore where inundation occurs, a far finer discretization of the hazard may be needed. The spatial distribution and intensity of the inundation is also greatly influenced by the co-seismic deformation (c.f. Volpe et al., 2019). We found that the impact of far offshore earthquakes was very similar for normal and reverse faulting, despite significant differences in the wave height time-series and implied consequences for the time before maximum inundation (Figure 9). For near-shore earthquakes, the impact metrics for normal faulting (rake = −90^o) are almost an inverse of the impact metrics for reverse faulting (rake = 90^o).

Finally, in a more limited sweep of strike and dip angles for offshore earthquakes, we demonstrate that the dependency upon the angle of dip decreases for deeper sources. For earthquakes close to the sea floor, shorter period waves are generated and there are significant changes in the tsunami intensity as we change both angles of strike and dip. For deeper earthquakes, only the longer period waves are excited and the tsunami intensity becomes less affected by the angle of dip. The strike angle is still significant for deeper earthquakes.

Despite the large number of simulations performed in this study, only a tiny fraction of parameter space has been covered. It is first noted that only two small stretches of coast have been studied, with a relatively simple relation between orientation of the coast to the locations of the prescribed earthquakes. Stretches of coastline with other geometries will likely experience tsunami impact which varies with source parameters in different ways. Some of the asymmetries displayed in this study are likely related to specific characteristics of the topo-bathymetry combined with the source orientations. These complications will increase as more complicated source geometries and more realistic models of heterogeneous slip are employed. Serra et al. (2021), for example, recently studied the sensitivity of tsunami height along the Iberian coast as a function of complex fault geometry and variable slip and find significant variability in the tsunami impact for different realizations. Such offshore studies illustrate how the relationship between slip distribution and bathymetric features influences near-shore wave height for extensive coastal stretches. However, the inundation problem represents an enormous increase in the computational cost and yet remains essential for assessing long term local tsunami hazard.

There may be several strategies for predicting onshore tsunami impact from offshore models. Davies et al. (2021) apply efficient Monte-Carlo sampling to estimate onshore PTHA from offshore. Liu et al. (2021) apply different Machine Learning (ML) techniques to predict tsunami amplitude at forecast points based upon short-term time-series at more offshore observation points; it is conceivable that an extension of such ideas will allow realistic inundation forecasts with high spatial resolution based upon offshore models and ML with a sufficient set of training data. For both ideas, simulated high-resolution inundation as a function of finely resolved parameter spaces will provide good validation of model forecasts. As PTHA with high-resolution inundation simulations becomes increasingly feasible for vast numbers of tsunamigenic scenarios, we seek a tsunami hazard quantification that is converged with respect to the granularity of the source discretization and modelling parameters. This work is a step in that direction.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: The Tsunami-HySEA code can be obtained from https://edanya.uma.es/hysea/index.php/models/tsunami-hysea (last accessed January 2022). The Digital Elevation Model used was constructed at INGV using data obtained from the following sources: https://www.gebco.net/ (GEBCO, accessed January 2022) https://www.ingv.it/it/monitoraggio-e-infrastrutture-per-la-ricerca/laboratori/laboratorio-geologia-e-geotecnologia (The Geological and Geotechnical Laboratory, INGV, last accessed January 2022) https://sextant.ifremer.fr/record/18ff0d48-b203-4a65-94a9-5fd8b0ec35f6/ (EMODnet Digital Bathymetry, DTM2018, last accessed January 2022) http://dati.protezionecivile.it/geoportalDPC/catalog/main/home.page (Dipartimento della Protezione Civile, Italy, last accessed January 2022).

Author Contributions

SG designed the experiment, performed the calculations, generated figures, and prepared the manuscript. SL, MVolpe, JS, BB, RT, and FR designed the calculation framework and computational grids, and contributed to interpretation of results. MdlA, JM, and CS-L extensively rewrote the Tsunami-HySEA code for this set of calculations. MVöge and SGli designed and tested the procedure for generating the input for the simulations and automating the procedure on HPC facilities. PL assisted with installing and implementing the procedure on Marconi-100 at CINECA. JCM assisted with installing and implementing the procedure on SAGA at the Norwegian University of Science and Technology. FL supervised the project, contributed to interpretation of the results, and prepared the manuscript. All authors contributed to review, revision, and quality control of the manuscript.

Funding

This work is partially funded by the European Union’s Horizon 2020 Research and Innovation Program under grant agreement No 823844 (ChEESE Center of Excellence, www.cheese-coe.eu). Additional funding was received from the Research Council of Norway under grant 322236 (INITIATOR). Computational resources made available through Sigma2/UNINETT on Saga at NTNU, Trondheim, Norway (in project nn5008 k) and through PRACE on Marconi-100 at CINECA, Rome, Italy (through a PRACE grant with project number 2020225386: Pra21_5,386/TsuHazAP). Some authors also received funding from the project “Assessment of Cascading Events triggered by the Interaction of Natural Hazards and Technological Scenarios involving the release of Hazardous Substances,” funded by the Italian Ministry MIUR PRIN (Progetti di Ricerca di Rilevante Interesse Nazionale) 2017—Grant 2017CEYPS8. This work benefited of the agreement between Istituto Nazionale di Geofisica e Vulcanologia and the Italian Presidenza del Consiglio dei Ministri, Dipartimento della Protezione Civile (DPC). This paper does not necessarily represent DPC official opinion and policies.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Maps and tsunami intensity plots were created using GMT software (Wessel et al., 2019). We gratefully acknowledge the contribution of three reviewers whose comments made significant improvements to the manuscript.

References

An, C., Liu, H., Ren, Z., and Yuan, Y. (2018). Prediction of Tsunami Waves by Uniform Slip Models. J. Geophys. Res. Oceans 123, 8366–8382. doi:10.1029/2018JC014363

CrossRef Full Text | Google Scholar

Basili, R., Brizuela, B., Herrero, A., Iqbal, S., Lorito, S., Maesano, F. E., et al. (2021). The Making of the NEAM Tsunami Hazard Model 2018 (NEAMTHM18). Front. Earth Sci. 8, 1–29. doi:10.3389/feart.2020.616594

CrossRef Full Text | Google Scholar

Bricker, J. D., Gibson, S., Takagi, H., and Imamura, F. (2015). On the Need for Larger Manning's Roughness Coefficients in Depth-Integrated Tsunami Inundation Models. Coastal Eng. J. 57, 1550005-1–1550005–13. doi:10.1142/S0578563415500059

CrossRef Full Text | Google Scholar

Burbidge, D., Mueller, C., and Power, W. (2015). The Effect of Uncertainty in Earthquake Fault Parameters on the Maximum Wave Height from a Tsunami Propagation Model. Nat. Hazards Earth Syst. Sci. 15, 2299–2312. doi:10.5194/nhess-15-2299-2015

CrossRef Full Text | Google Scholar

Davies, G., Griffin, J., Løvholt, F., Glimsdal, S., Harbitz, C., Thio, H. K., et al. (2018). A Global Probabilistic Tsunami hazard Assessment from Earthquake Sources. Geol. Soc. Lond. Spec. Publications 456, 219–244. doi:10.1144/SP456.5

CrossRef Full Text | Google Scholar

Davies, G., Horspool, N., and Miller, V. (2015). Tsunami Inundation from Heterogeneous Earthquake Slip Distributions: Evaluation of Synthetic Source Models. J. Geophys. Res. Solid Earth 120, 6431–6451. doi:10.1002/2015JB012272

CrossRef Full Text | Google Scholar

Davies, G. (2019). Tsunami Variability from Uncalibrated Stochastic Earthquake Models: Tests against Deep Ocean Observations 2006-2016. Geophys. J. Int. 218, 1939–1960. doi:10.1093/gji/ggz260

CrossRef Full Text | Google Scholar

Davies, G., Weber, R., Wilson, K., and Cummins, P. (2021). From Offshore to Onshore Probabilistic Tsunami hazard Assessment via Efficient Monte-Carlo Sampling. earthArXiv. doi:10.31223/X5J04W

CrossRef Full Text | Google Scholar

de la Asunción, M., Castro, M. J., Fernández-Nieto, E. D., Mantas, J. M., Acosta, S. O., and González-Vida, J. M. (2013). Efficient GPU Implementation of a Two Waves TVD-WAF Method for the Two-Dimensional One Layer Shallow Water System on Structured Meshes. Comput. Fluids 80, 441–452. doi:10.1016/j.compfluid.2012.01.012

CrossRef Full Text | Google Scholar

Gayer, G., Leschka, S., Nöhren, I., Larsen, O., and Günther, H. (2010). Tsunami Inundation Modelling Based on Detailed Roughness Maps of Densely Populated Areas. Nat. Hazards Earth Syst. Sci. 10, 1679–1687. doi:10.5194/nhess-10-1679-2010

CrossRef Full Text | Google Scholar

Geist, E. L. (2002). Complex Earthquake Rupture and Local Tsunamis. J. Geophys. Res. 107, 2086. doi:10.1029/2000JB000139

CrossRef Full Text | Google Scholar

Geist, E. L. (1998). Local Tsunamis and Earthquake Source Parameters. Adv. Geophys. 39, 117–209. doi:10.1016/S0065-2687(08)60276-9

CrossRef Full Text | Google Scholar

Gibbons, S. J., Lorito, S., Macías, J., Løvholt, F., Selva, J., Volpe, M., et al. (2020). Probabilistic Tsunami Hazard Analysis: High Performance Computing for Massive Scale Inundation Simulations. Front. Earth Sci. 8, 1–20. doi:10.3389/feart.2020.591549

CrossRef Full Text | Google Scholar

Gica, E., Teng, M. H., Liu, P. L.-F., Titov, V., and Zhou, H. (2007). Sensitivity Analysis of Source Parameters for Earthquake-Generated Distant Tsunamis. J. Waterway, Port, Coastal, Ocean Eng. 133, 429–441. doi:10.1061/(asce)0733-950x(2007)133:6(429)

CrossRef Full Text | Google Scholar

Glimsdal, S., Pedersen, G. K., Harbitz, C. B., and Løvholt, F. (2013). Dispersion of Tsunamis: Does it Really Matter? Nat. Hazards Earth Syst. Sci. 13, 1507–1526. doi:10.5194/nhess-13-1507-2013

CrossRef Full Text | Google Scholar

Grezio, A., Babeyko, A., Baptista, M. A., Behrens, J., Costa, A., Davies, G., et al. (2017). Probabilistic Tsunami Hazard Analysis: Multiple Sources and Global Applications. Rev. Geophys. 55, 1158–1198. doi:10.1002/2017RG000579

CrossRef Full Text | Google Scholar

Griffin, J., Latief, H., Kongko, W., Harig, S., Horspool, N., Hanung, R., et al. (2015). An Evaluation of Onshore Digital Elevation Models for Modeling Tsunami Inundation Zones. Front. Earth Sci. 3. doi:10.3389/feart.2015.00032

CrossRef Full Text | Google Scholar

Harbitz, C. B., Glimsdal, S., Bazin, S., Zamora, N., Løvholt, F., Bungum, H., et al. (2012). Tsunami hazard in the Caribbean: Regional Exposure Derived from Credible Worst Case Scenarios. Continental Shelf Res. 38, 1–23. doi:10.1016/j.csr.2012.02.006

CrossRef Full Text | Google Scholar

Horspool, N., Pranantyo, I., Griffin, J., Latief, H., Natawidjaja, D. H., Kongko, W., et al. (2014). A Probabilistic Tsunami hazard Assessment for Indonesia. Nat. Hazards Earth Syst. Sci. 14, 3105–3122. doi:10.5194/nhess-14-3105-2014

CrossRef Full Text | Google Scholar

Kaiser, G., Scheele, L., Kortenhaus, A., Løvholt, F., Römer, H., and Leschka, S. (2011). The Influence of Land Cover Roughness on the Results of High Resolution Tsunami Inundation Modeling. Nat. Hazards Earth Syst. Sci. 11, 2521–2540. doi:10.5194/nhess-11-2521-2011

CrossRef Full Text | Google Scholar

Kajiura, K. (1980). “Some Statistics Related to Observed Tsunami Heights along the Coast of Japan,” in Tsunamis: Their Science and Engineering (Dordrecht: Springer Netherlands), 131–145. doi:10.1007/978-94-009-7172-1_11

CrossRef Full Text | Google Scholar

Li, L., Switzer, A. D., Chan, C.-H., Wang, Y., Weiss, R., and Qiu, Q. (2016). How Heterogeneous Coseismic Slip Affects Regional Probabilistic Tsunami hazard Assessment: A Case Study in the South China Sea. J. Geophys. Res. Solid Earth 121, 6250–6272. doi:10.1002/2016JB013111

CrossRef Full Text | Google Scholar

Liu, C. M., Rim, D., Baraldi, R., and LeVeque, R. J. (2021). Comparison of Machine Learning Approaches for Tsunami Forecasting from Sparse Observations. Pure Appl. Geophys. 178, 5129–5153. doi:10.1007/s00024-021-02841-9

CrossRef Full Text | Google Scholar

Lorito, S., Alessandro, A., Cugliari, L., Romano, F., Tonini, R., Valbonesi, C., et al. (2021). Tsunami hazard, Warning, and Risk Reduction in Italy and the Mediterranean Sea: State of the Art, Gaps, and Future Solutions. Turkish J. Earth Sci. 30, 882–897. doi:10.3906/yer-2110-7

CrossRef Full Text | Google Scholar

Lorito, S., Selva, J., Basili, R., Romano, F., Tiberti, M. M., and Piatanesi, A. (2015). Probabilistic Hazard For Seismically Induced Tsunamis: Accuracy And Feasibility Of Inundation Maps. Geophys. J. Int. 200, 574–588. doi:10.1093/gji/ggu408

CrossRef Full Text | Google Scholar

Lorito, S., Tiberti, M. M., Basili, R., Piatanesi, A., and Valensise, G. (2008). Earthquake-generated Tsunamis in the Mediterranean Sea: Scenarios of Potential Threats to Southern Italy. J. Geophys. Res. 113, B01301. doi:10.1029/2007JB004943

CrossRef Full Text | Google Scholar

Løvholt, F., Glimsdal, S., Harbitz, C. B., Horspool, N., Smebye, H., de Bono, A., et al. (2014). Global Tsunami hazard and Exposure Due to Large Co-seismic Slip. Int. J. Disaster Risk Reduction 10, 406–418. doi:10.1016/j.ijdrr.2014.04.003

CrossRef Full Text | Google Scholar

Løvholt, F., Glimsdal, S., Harbitz, C. B., Zamora, N., Nadim, F., Peduzzi, P., et al. (2012a). Tsunami hazard and Exposure on the Global Scale. Earth-Science Rev. 110, 58–73. doi:10.1016/j.earscirev.2011.10.002

CrossRef Full Text | Google Scholar

Løvholt, F., Kühn, D., Bungum, H., Harbitz, C. B., and Glimsdal, S. (2012b). Historical Tsunamis and Present Tsunami hazard in Eastern Indonesia and the Southern Philippines. J. Geophys. Res. 117, 1–19. doi:10.1029/2012JB009425

CrossRef Full Text | Google Scholar

Løvholt, F., Lynett, P., and Pedersen, G. (2013). Simulating Run-Up on Steep Slopes with Operational Boussinesq Models; Capabilities, Spurious Effects and Instabilities. Nonlin. Process. Geophys. 20, 379–395. doi:10.5194/npg-20-379-2013

CrossRef Full Text | Google Scholar

Løvholt, F., Pedersen, G., Bazin, S., Kühn, D., Bredesen, R. E., and Harbitz, C. (2012c). Stochastic Analysis of Tsunami Runup Due to Heterogeneous Coseismic Slip and Dispersion. J. Geophys. Res. Ocean. 117, n/a. doi:10.1029/2011JC007616

CrossRef Full Text | Google Scholar

Macías, J., Castro, M. J., Ortega, S., Escalante, C., and González-Vida, J. M. (2017). Performance Benchmarking of Tsunami-HySEA Model for NTHMP's Inundation Mapping Activities. Pure Appl. Geophys. 174, 3147–3183. doi:10.1007/s00024-017-1583-1

CrossRef Full Text | Google Scholar

Macías, J., Mercado, A., González-Vida, J. M., Ortega, S., and Castro, M. J. (2016). Comparison and Computational Performance of Tsunami-HySEA and MOST Models for LANTEX 2013 Scenario: Impact Assessment on Puerto Rico Coasts. Pure Appl. Geophys. 173, 3973–3997. doi:10.1007/s00024-016-1387-8

CrossRef Full Text | Google Scholar

McCloskey, J., Antonioli, A., Piatanesi, A., Sieh, K., Steacy, S., Nalbant, S. S., et al. (2007). Near-field Propagation of Tsunamis from Megathrust Earthquakes. Geophys. Res. Lett. 34, L14316. doi:10.1029/2007GL030494

CrossRef Full Text | Google Scholar

Murphy, S., Scala, A., Herrero, A., Lorito, S., Festa, G., Trasatti, E., et al. (2016). Shallow Slip Amplification and Enhanced Tsunami hazard Unravelled by Dynamic Simulations of Mega-Thrust Earthquakes. Sci. Rep. 6, 35007. doi:10.1038/srep35007

PubMed Abstract | CrossRef Full Text | Google Scholar

NCEI (2021). National Geophysical Data Center/World Data Service: NCEI/WDS Global Historical Tsunami Database. doi:10.7289/V5PN93H7

CrossRef Full Text | Google Scholar

Okada, Y. (1985). Surface Deformation Due to Shear and Tensile Faults in a Half-Space. Bull. Seismol. Soc. Am. 75, 1135–1154. doi:10.1785/bssa0750041135

CrossRef Full Text | Google Scholar

Parsons, T., and Geist, E. L. (2008). Tsunami Probability in the Caribbean Region. Pure Appl. Geophys. 165, 2089–2116. doi:10.1007/s00024-008-0416-7

CrossRef Full Text | Google Scholar

Scala, A., Lorito, S., Romano, F., Murphy, S., Selva, J., Basili, R., et al. (2020). Effect of Shallow Slip Amplification Uncertainty on Probabilistic Tsunami Hazard Analysis in Subduction Zones: Use of Long-Term Balanced Stochastic Slip Models. Pure Appl. Geophys. 177, 1497–1520. doi:10.1007/s00024-019-02260-x

CrossRef Full Text | Google Scholar

Selva, J., Amato, A., Armigliato, A., Basili, R., Bernardi, F., Brizuela, B., et al. (2021). Tsunami Risk Management for Crustal Earthquakes and Non-seismic Sources in Italy. Riv. Nuovo Cim. 44, 69–144. doi:10.1007/s40766-021-00016-9

CrossRef Full Text | Google Scholar

Selva, J., Tonini, R., Molinari, I., Tiberti, M. M., Romano, F., Grezio, A., et al. (2016). Quantification of Source Uncertainties in Seismic Probabilistic Tsunami Hazard Analysis (SPTHA). Geophys. J. Int. 205, 1780–1803. doi:10.1093/gji/ggw107

CrossRef Full Text | Google Scholar

Serra, C. S., Martínez‐Loriente, S., Gràcia, E., Urgeles, R., Gómez de la Peña, L., Maesano, F. E., et al. (2021). Sensitivity of Tsunami Scenarios to Complex Fault Geometry and Heterogeneous Slip Distribution: Case‐Studies for SW Iberia and NW Morocco. J. Geophys. Res. Solid Earth 126. doi:10.1029/2021JB022127

CrossRef Full Text | Google Scholar

Tonini, R., Basili, R., Maesano, F. E., Tiberti, M. M., Lorito, S., Romano, F., et al. (2020). Importance of Earthquake Rupture Geometry on Tsunami Modelling: the Calabrian Arc Subduction Interface (Italy) Case Study. Geophys. J. Int. 223, 1805–1819. doi:10.1093/gji/ggaa409

CrossRef Full Text | Google Scholar

Tonini, R., Di Manna, P., Lorito, S., Selva, J., Volpe, M., Romano, F., et al. (2021). Testing Tsunami Inundation Maps for Evacuation Planning in Italy. Front. Earth Sci. 9. doi:10.3389/feart.2021.628061

CrossRef Full Text | Google Scholar

Volpe, M., Lorito, S., Selva, J., Tonini, R., Romano, F., and Brizuela, B. (2019). From Regional to Local SPTHA: Efficient Computation of Probabilistic Tsunami Inundation Maps Addressing Near-Field Sources. Nat. Hazards Earth Syst. Sci. 19, 455–469. doi:10.5194/nhess-19-455-2019

CrossRef Full Text | Google Scholar

Wessel, P., Luis, J. F., Uieda, L., Scharroo, R., Wobbe, F., Smith, W. H. F., et al. (2019). The Generic Mapping Tools Version 6. Geochem. Geophys. Geosyst. 20, 5556–5564. doi:10.1029/2019GC008515

CrossRef Full Text | Google Scholar

Williamson, A. L., Rim, D., Adams, L. M., LeVeque, R. J., Melgar, D., González, F. I., et al. (2020). A Source Clustering Approach for Efficient Inundation Modeling and Regional Scale Probabilistic Tsunami Hazard Assessment. Front. Earth Sci. 8, 1–22. doi:10.3389/feart.2020.591663

CrossRef Full Text | Google Scholar