Importance of earthquake rupture geometry on tsunami modelling: the Calabrian Arc subduction interface (Italy) case study

The behaviour of tsunami waves at any location depends on the local morphology of the coasts, the encountered bathymetric features, and the characteristics of the source. However, the importance of accurately modelling the geometric properties of the causative fault for simulations of seismically induced tsunamis is rarely addressed. In this work, we analyse the effects of using two different geometric models of the subduction interface of the Calabrian Arc (southern Italy, Ionian Sea) onto the simulated tsunamis: a detailed 3-D subduction interface obtained from the interpretation of a dense network of seismic reﬂection proﬁles, and a planar interface that roughly approximates the 3-D one. These models can be thought of as representing two end-members of the level of knowledge of fault geometry. We deﬁne three hypothetical earthquake ruptures of different magnitudes ( M w 7.5, 8.0, 8.5) on each geometry. The resulting tsunami impact is evaluated at the 50-m isobath in front of coastlines of the central and eastern Mediterranean. Our results show that the source geometry imprint is evident on the tsunami waveforms, as recorded at various distances and positions relative to the source. The absolute differences in maximum and minimum wave amplitudes locally exceed one metre, and the relative differences remain systematically above 20 per cent with peaks over 40 per cent. We also observe that tsunami energy directivity and focusing due to bathymetric waveguides take different paths depending on which fault is used. Although the differences increase with increasing earthquake magnitude, there is no simple rule to anticipate the different effects produced by these end-member models of the earthquake source. Our ﬁndings suggest that oversimpliﬁed source models may hinder our fundamental understanding of the tsunami impact and great care should be adopted when making simplistic assumptions regarding the appropriateness of the planar fault approximation in tsunami studies. We also remark that the geological and geophysical 3-D fault characterization remains a crucial and unavoidable step in tsunami hazard analyses.


I N T RO D U C T I O N
Tsunami generation, propagation and impact depend on several factors, such as the bathymetry (Satake 1988), the relative position between source and target (Geist 2009) and the local morphology of coastal areas (Tonini et al. 2014). For seismically generated tsunamis, their aspects depend on the parameters describing the seismic source (Geist 1998;Davies 2019), such as the focal mechanism (Tonini et al. 2011), the slip and the shear modulus distribution (Geist & Bilek 2001;Goda 2015;Murphy et al. 2016;Davies & Griffin 2019;Sallares & Ranero 2019;Scala et al. 2020) and the variability with the depth of seismic rupture mechanical features (Murphy et al. 2018;Scala et al. 2019).
Among all the seismic source features, the geometry of the fault interface is a less specifically explored contribution at defining the tsunami characteristics. Being the surface onto which the earthquake rupture nucleates and propagates, the fault geometry plays a fundamental role in determining the coseismic seafloor deformation and then in defining the corresponding sea surface displacement, which in turn is the initial condition of tsunami propagation.
Most of the tsunamis are caused by large to great earthquakes generated along the converging boundaries of the tectonic plates (Lay 2015;Bletery et al. 2016;Lorito et al. 2016;Wang & Tréhu 2016;Grezio et al. 2017;Romano et al. 2020). The rupture area of megathrust earthquakes in subduction zones extends for tens to hundreds of kilometres both downdip and along strike (Strasser et al. 2010), and deviates from the planar geometry (Hayes et al. 2018); hence, the larger the earthquake size, the more the curvature should matter.
Using a fault model that honours as much as possible the real 3-D geometry is, in principle, beneficial to reduce the epistemic uncertainty in the modelling when studying the characteristics of the seismic source. In turn, properly modelling the geophysical observables associated with an earthquake (e.g. geodetic, seismic, or tsunami data) is one of the crucial points when inverting data to estimate the seismic source properties Ryan et al. 2015;Lorito et al. 2016). Using a planar or a 3-D fault surface might result in a more accurate seismic rupture history model, as discussed for earthquakes spanning a wide range of magnitudes and different types of inverted data, such as for the M w 7.5, 1999 Chi-Chi earthquake (Lee et al. 2006), the M w 6.8, 2003 Boumerdes-Zemmouri earthquake (Belabbès et al. 2009), the M w 9.5, 1960 Valdivia earthquake (Moreno et al. 2009), or the M w 8.3, 2013 Santa Cruz Islands earthquake (Lay et al. 2013;Hayes et al. 2014;Romano et al. 2015). Accordingly, this aspect becomes particularly important when the estimated or tested seismic source is used to model the tsunami generation, propagation in the open sea and impact onto the coast.
It has been argued that like the other seismic source features, the fault shape or curvature count more in the near-field than in the far-field (Geist 2002). Hence, both in long-term Probabilistic Tsunami Hazard Analyses (PTHA) (González et al. 2009;Omira et al. 2014) and for quasi-real-time tsunami forecasting (Titov et al. 1999;Tang et al. 2009), planar geometries are sometimes adopted for modelling distant tsunamis while for local tsunamis a subdivision of the fault into subfaults of different strike and dip is usually adopted (González et al. 2009;Lorito et al. 2015).
Rather than making a distinction based on the relative position of the source and the target coastline, a different perspective was introduced in the PTHA methodology by Selva et al. (2016). This approach was subsequently adopted both for the NEAMTHM18 hazard model Basili et al. 2019) and for the site-specific PTHA approach proposed by Volpe et al. (2019). The modelling strategy is selected based on how well the fault geometry can be constrained. For some relatively well-constrained sources, mostly subduction interfaces, a 3-D geometry and discretization into small triangles with varying orientation is adopted; for others, mostly offshore crustal settings hosting background seismicity around major faults, a simpler planar model is adopted in the absence of substantial constraints on their geometry.
In this work, we focus on the subduction zone of the Calabrian Arc (Fig. 1a), located in the Ionian Sea (central Mediterranean Sea) to present an analysis aimed to explore the importance of including realistic geometric features of the subduction interface to model seismically induced tsunamis.

DATA A N D M E T H O D S
We set up a test case by defining and modelling tsunami scenarios generated by earthquake ruptures located in the Calabrian Arc ( Fig. 1a), using two subduction interface geometry models: an accurate 3-D geometry (Fig. 1b) and a simplified planar geometry (Fig. 1c). We evaluated the effects of the 3-D complexities in the seismic source geometry in terms of the resulting tsunami waveforms, recorded along the 50-m isobath of the central and eastern Mediterranean and, in particular, the associated maximum (A max ) and minimum (A min ) tsunami wave amplitude (positive and negative water elevation relative to mean sea level, respectively). To isolate the relative contributions of the two models to the tsunami modelling, we kept constant all the earthquake parameters which do not directly depend on the geometry (magnitude, epicentre, area and slip). In the next paragraphs, we provide more details of the case study.

Tectonic setting of the case study
The Calabrian Arc is one of the smallest documented subduction worldwide (Hayes et al. 2018) and involves one of the possibly oldest oceanic crust in the subduction process (Müller et al. 2008;Speranza et al. 2012).
The definition of the geometric model of a tectonic structure depends on the type and amount of available data and can be derived from several types of methodologies (e.g. geology, geophysics and seismology). Here, the subducting slab is well imaged by the intraslab seismicity, which highlights a steeply dipping portion below 80 km, and the presence of a slab window in the northern part (Faccenna 2005) evolving in a horizontal incipient tear in the central sector (Neri et al. 2009;. The interpretation of seismic reflection profiles constrains the shallow (<20 km) portion of the subduction zone, showing an almost flat detachment level between the tectonically stacked Messinian evaporites and the relatively undeformed Mesozoic to Palaeocene succession overlying the Ionian oceanic crust. The transition toward the inner part of the subduction zone is characterized by a flat-ramp geometry (Fig. 1d). The detachment level of the subduction interface should be located at the top of the Mesozoic succession of the Ionian basin, which is in turn involved in the deformation within the accretionary wedge. The subduction terminates laterally with a gradual transition to a continental collision in the northern part (Ferranti et al. 2014;Volpi et al. 2017;Teofilo et al. 2018) and to a major subvertical lithospheric tear on the Alfeo Fault System (Polonia et al. 2011;Gallais et al. 2013;Gutscher et al. 2016;Maesano et al. 2020).
The seismogenic behaviour of the Calabrian Arc is still far from being well understood.
Although it accommodates the convergence between the African and European plates with rates ranging between 1.5 and 5 mm yr −1 (Devoti et al. 2008;Carafa et al. 2018), based on the ISIDe Working Group (2016) earthquake catalogue, it shows only intraslab significant seismic activity below 40 km depth and no relevant seismicity in the shallower portion of the interface. Nevertheless, the subduction interface is often considered capable of generating large earthquakes and suspected to be responsible for some past earthquakes and tsunamis (Gutscher et al. 2006;Polonia et al. 2013). This situation is somewhat similar, though on different spatial and temporal scales, to the Cascadia subduction zone where, despite the relatively low seismicity, evidence was found that it generated significant earthquakes and tsunamis in the past. Wang & Tréhu (2016) invoke the presence of heterogeneous locking to explain those observations. In the Calabrian Arc, however, recent geodetic observations, coupled with geodynamic modelling, led to the hypothesis that the southern part of the Calabrian Arc could either be locked or partly locked (Carafa et al. 2018) or negligibly active (Nijholt et al. 2018).

Subduction interface geometric models
We considered two end-member models of the slab interface: a complex and realistic 3-D geometry and a simple planar geometry.
The first is a curved 3-D surface, the more realistic geometric reconstruction of the subduction interface (Fig. 1b), which is the outcome of the analysis and interpretation of several seismic reflection profiles and the recorded seismicity of the area . The 3-D surface derives from the interpretation of 54 seismic reflection lines for a total of 8658 km of surveys exploring about 15 km depth below the seafloor within an area of ∼100 000 km 2 . The seismic data set grid has a spacing between profile pairs ranging from 10 to 50 km (see Fig. S1 in the Supporting Information) and has 438 intersections among the profiles which helped to check the consistency of the interpretation. A 3-D surface was then obtained through Delaunay triangulation of the interpreted profiles. This interpretation in the time domain was converted to metric units using the 3-D time-to-depth conversion algorithm by Maesano & D'Ambrogi (2017) and then connected with another surface based on the recorded seismicity of the area. The coupled data set was finally resampled using a 15-km-spaced grid .
The second is a planar fault (Fig. 1c) which is an extreme approximation of the 3-D geometry (see the profiles in Fig. 1e for comparison in the vertical dimension), assumed to represent the common case of limited knowledge of tectonic structures or their typical approximation.
In both cases (3-D and planar), we modelled the interface by a mesh of triangular elements using the Cubit meshing software (https://cubit.sandia.gov/, Casarotti et al. 2008). Although the triangular mesh is not strictly necessary in the planar case, we adopted this strategy to ensure the equivalence of the procedure in most aspects and to test alternative parametrizations not reported here. We constrained the triangles to have approximately the same area in both meshes to minimize possible dependence of the results from the discretization refinement. For the 3-D fault, each triangular element represents a subfault characterized by a specific orientation (strike and dip) and slip direction (rake, Fig. S3 in the Supporting Information). These angular values are based on the fault surface geometry and the available geological and tectonic knowledge. In particular, the rake values were imposed to reflect the average direction of convergence (Fig. 1c,and Fig. S2 in the Supporting Inormation) of the Calabrian Arc (Devoti et al. 2011). For the planar fault, all the triangles have the same strike and dip as the containing plane (constant values). The rake was fixed at 90 • which given the strike of the planar fault corresponds to the average direction of convergence (Fig. 1d,and Fig. S2 in the Supporting Information) of the Calabrian Arc (Devoti et al. 2011) and reflects the dominantly thrusting mechanisms of the subduction interface. This configuration is commonly used when only poor kinematic constraints are available because it maximizes the vertical component of seafloor displacement, which is more efficient than the horizontal component in generating tsunamis and avoid underestimation of its impact. Such configuration is also consistent with negligible strain partitioning for both fault geometries and roughly symmetric rake dispersion of the 3-D geometry in comparison with the planar geometry ( Fig. S3 in the Supporting Information).

Earthquake and tsunami scenarios
We defined a set of earthquake scenarios and modelled the associated coseismic seafloor deformation for the two slab interface models (Fig. 2a), varying the magnitude (M w 7.5, 8.0 and 8.5) to explore the possible dependence between the geometry and the size of the earthquake. The extent of the rupture area of each earthquake was determined from the earthquake scaling relations by Strasser et al. (2010) for subduction interfaces. The triangular elements of the meshes representing the rupture areas were identified by starting from an assigned geometrical centre of the earthquake rupture (longitude 16.75 • E, latitude 37.94 • N; Fig. 2a) and by iteratively adding neighbouring mesh elements, as in the procedure described by Scala et al. (2020). This procedure was arrested when the area of the selected triangles exceeded the expected area predicted by the scaling relation. The influence of heterogeneous slip distributions on the tsunami modelling (Geist 2002) and the tsunami hazard assessment (Li et al. 2016) has been widely investigated. In our analysis, we used uniform slip distributions to isolate, at the first order, the effect of different rupture geometries. Therefore, the slip value predicted by the scaling relations was assigned to each selected triangular element. To this end, we used the classical empirical relations between moment magnitude and seismic moment (Hanks & Kanamori 1979), and the equivalence of seismic moment to the product between rupture area, average slip and shear modulus (here fixed to 33 GPa).
Finally, we computed the tsunami initial conditions (Fig. 2b) of the earthquake scenarios by superposing the seafloor deformation due to each triangular element contained in the modelled rupture area, using the Volterra's formulation of elastic dislocation theory applied to a triangular source buried in an elastic half-space (Meade 2007). Since the water column acts as a low-pass filter between the sea-bottom and the sea-surface displacements (resulting in shortwavelength attenuation), we applied a 2-D filter to the static vertical seafloor deformation field (Kajiura 1963).
Tsunami simulations, for the equivalent of 8 h of propagation, were performed with the finite-volume Tsunami-HySEA GPU code (de la Asunción et al. 2013) that solves the nonlinear shallow water equations and was extensively benchmarked (Macías et al. 2016;Macías et al. 2017;Macías et al. 2020). The computational grid has a spatial resolution of 30 arcsec and covers the central and eastern Mediterranean basin. Seismically generated tsunamis are usually treated as non-frequency-dispersive long waves whose behaviour is essentially linear as long as their amplitude is much smaller than the local sea depth. Therefore, to minimize the influence of the nonlinear processes related to the inundation phase and keep the focus on the effect of the source, we recorded the tsunami waveforms at 394 points of interest (POIs), with an average spacing of ∼2 km, located along the 50-m isobath (Fig. 3). The analysis of the values along the POIs on an isobath is a commonly used procedure to make a rough (relative to detailed inundation analysis) but consistent evaluation of tsunami impact or tsunami hazard at a regional scale for the corresponding coastal target areas (Kamigaichi 2009;Sørensen et al. 2012;Basili et al. 2019). The two triangular meshes of the subduction interface, the selection of mesh elements for each earthquake rupture, and the maximum (A max ) and minimum (A min ) tsunami wave amplitudes are described in the Supporting Information together with the corresponding data files in comma-separated values format.

R E S U LT S
Starting from the two geometry models of the slab interface and the fixed earthquake rupture centre, we simulated two tsunami scenarios for each of the three magnitudes. The parameters of the six rupture models are listed in Table 1. The selection of the triangular elements is reflected in the rupture areas, which are different for the two geometries. These differences, however, are very small and, coupled with the fixed average slip, produce a relatively small uncertainty on the actual seismic moment. The recalculated moment magnitudes have maximum deviations of 0.02 (Table 1), which is well below the typical uncertainty affecting their estimations. More significant are the differences in rupture depth (Fig. 2a). Indeed, the average rupture depth for the planar approximation is 10 km deeper than the average rupture depth for the 3-D geometry. Moreover, the differences in the strike, dip and rake angles are significant (Fig. 4), with their variability increasing with increasing earthquake magnitude, because of the progressively larger occupancy of the rupture area on the mesh. In particular, the distribution of strike and rake angles for the 3-D geometry is asymmetrical in comparison with the fixed values of the planar geometry in the case of M w 7.5 and 8.0 earthquake ruptures and rather symmetrical in the case of the M w 8.5 earthquake rupture. The M w 7.5 and the 8.0 earthquake ruptures on the 3-D geometry result to be oblique reverse, that is with a reduced dip-slip component relative to the planar geometry.
The calculated sea water displacement (i.e. the tsunami initial condition) corresponding to each rupture model is shown in maps and profiles in Fig. 2(b). The maps show how the differences in the source geometry are immediately reflected on the initial condition of the generated tsunamis, being all other source parameters approximately the same. The most evident difference is the arched shape of both the upthrown and the downthrown areas of the 3-D geometry, compared to the straighter shape of the planar geometry. Nonetheless, it is worth noting that, in terms of maximum initial water displacement, the two geometry models produce similar values (Fig. 2b); this effect is due to the trade-off between fault depth and dip. Moreover, the initial water displacement profiles along the A-B transect (Fig. 2b, second and fourth columns) for the M w 8.5 scenarios highlight how a leading negative wave towards the Calabrian coasts for the 3-D geometry is missing because of the interface curvature (maximum subsidence inland).
The maximum wave amplitudes (A max ) are computed for all earthquake scenarios in the whole domain (Fig. 5) of the central and eastern Mediterranean Sea, which includes the maxima at the 394 selected POIs. Comparing the maps for the planar sources with the maps of the 3-D sources (Fig. 5), we notice how the geometric differences of the seismic source are further reflected in different propagation patterns within the entire domain. It is particularly evident that the pattern of the A max generated by planar ruptures is very focused in the direction orthogonal to the fault strike, whereas the pattern generated by the 3-D ruptures has jagged maxima in diverging directions.
The different aspects of the tsunami impact onto the coast can be compared from the profiles of the A max and A min retained at the POIs at three selected coastal segments in Italy, Libya and Greece (Fig. 6). The first is in the very near-field of the tsunami source, where the highest impact is expected (1000 km of coast in southern Italy, Fig. 6, left-hand column). Following the terminology used in Geist (2009), the near-field can be separated into two distinct regimes: broadside and oblique (grey and white backgrounds, respectively, in Fig. 6, left-hand column). The second segment is in the far-field of the source on the axis of predominant directivity (orthogonal to fault strike), where the expected A max and A min values are also among the highest (1300 km of coast in northern Africa, Fig. 6, middle column). The third is at an intermediate distance from the source and in the direction orthogonal to the axis of predominant directivity, where the expected A max and A min values are lower than in the other two cases (500 km of coast in western Greece, Fig. 6, right-hand column).
In deep waters, the tsunami wave propagation is strongly influenced by the most relevant features of the bathymetry, such as seamounts, escarpments, ridges, valleys, canyons, or any relief sharp enough to produce a variation of the water depth. Indeed, the offshore motion of tsunami waves is mainly governed by linear shallow-water equations, and the wave speed is proportional to the square root of the seafloor depth. Consequently, waves move slower in shallower waters, and most of the tsunami energy follows this refraction effect as it emerges in the pattern of the A max for the scenarios shown in Fig. 5. Note that the colour scale in Fig. 5 is different in the three-panel pairs for the different earthquake magnitudes to highlight the water depth-dependent wave-propagation features. The effect of the two different geometries on the wave directivity is evident. Although both geometric models have the same predominant directivity of the tsunami energy in the NW-SE direction, the simulations based on the 3-D geometry present additional propagation paths, especially toward the coasts of Albania and northern Greece, due to the curvature of the subduction interface (Fig. 2). This effect is particularly evident for the highest magnitude, and it means that a tsunami scenario that uses an oversimplified earthquake rupture geometry, that is the planar fault, may significantly misrepresent the wave amplitude at many locations. However, focusing only on bathymetric effects or only on source directivity associated to average source strike and curvature is just a descriptive simplification, since there is an interplay between the two as predicted by basic refraction theory.
Nevertheless, the systematic wave amplification toward the coasts between eastern Tunisia and western Libya (see also Fig. 6, middle column) is clearly controlled by the bathymetry. This pattern of wave propagation is present in a similar fashion for all the earthquake scenarios for both geometries. Its origin can be attributed to the presence of the Malta Escarpment (an east-facing step 300 km long, 1000-1500 m high), which is known to enhance the wave amplitude of any tsunami propagating westward from the eastern side of the basin (Lorito et al. 2008;Basili et al. 2013). Because of the seafloor morphology in front of the Gulf of Gabes (Tunisia), the reference 50 m bathymetric depth is several tens of kilometres away from the coast and the A max and A min at those POIs cannot be safely considered as representative of the actual impact on the coast in comparison to all the other POIs.
In the three selected coastal stretches (Fig. 6), the general pattern of both A max and A min values is very similar for both geometric models; nonetheless, some of the highest peaks occur at slightly different locations. In the near-field, there is no clear dominance of one model over the other for the M w 7.5 and 8.0 scenarios. For the largest earthquake magnitude (M w 8.5), however, the highest A max value for the planar geometry overtakes by almost a metre the highest A max value of the 3-D geometry. Moreover, the A max distributions along the coast confirm the different phenomenological behaviour between oblique and broadside regimes, the latter focusing the higher amplitudes. In the far field, the planar geometry seems to dominate in most cases, and especially for the M w 8.5 earthquake scenario. As mentioned above, the combined effect of bathymetry and, prevalently, directivity is particularly evident in this stretch of coast, where larger amplitudes are focused between POIs numers 50 and 70. The opposite occurs in the Greek coastal stretch where the A max values for the 3-D geometry significantly exceed those for the planar geometry in most POIs for all scenarios. It is also worth noting that the relative difference of A max values very often exceeds 40 per cent, whereas it remains systematically below 20 per cent only in the far-field case for the M w 7.5 earthquake. Importantly, the pattern of the relative difference is uncorrelated with that of the absolute values, showing that significant differences between the two source models arise in all cases, and even where A max or A min values are relatively small. This observation suggests that using an oversimplified source model can hinder even our basic understanding of the tsunami impact. This consideration is further confirmed by full tsunami waveforms recorded at the same selected POIs (Fig. 7). It can be observed that even secondary maximum and minimum amplitude peaks can record significant differences, depending on which of the two geometric models is used (Fig. 7b).
Moreover, a deeper look at the waveform plots suggests a possible change in the main frequency components of tsunami waves. The frequency spectrum of the signals (Fig. 7c) confirms this hypothesis, showing some differences in the main spectral components (see as an illustrative example the bottom right subplot in Fig. 7c), which is a very relevant feature to study resonance phenomena in semi-closed natural bays, basins or harbours (Rabinovich & Thomson 2007). In this specific example, it seems that the planar sources may have the effect of systematically enhance the spectral peak at ∼37 min, and perhaps in some cases of a second peak at ∼26 min, which might represent the spectral signature of the source. The 3-D sources enlighten a peak at 22 min, in particular for M w 8.5. In the framework of the linear theory, the tsunami wave period approximately increases with increasing characteristic source size (earthquake magnitude) and decreases with the square root of increasing water depth in the source area. For the considered sources, both parameters are controlled by the curvature (Fig. 2), suggesting this variation of the source period.
The scatter plots (Fig. 8) show the maximum wave amplitudes A max obtained for the 3-D geometry versus those obtained for the planar geometry for the selection of points in the central and eastern Mediterranean Sea (Fig. 3). Results are separated for the different earthquake magnitude (M w 7.5, 8.0 and 8.5) in the first three panels, and grouped for all the magnitudes together in the fourth panel. These plots show that the scatter of the A max values (rmse), around the ideal line of the 1:1 ratio, increases with increasing earthquake magnitude. For the largest magnitude, the orthogonal regression line is well below the 1:1 line, and the scatter of points is rather asymmetrical; it is thus evident that in this case the planar geometry slightly overestimates the A max values in the whole domain in comparison with the 3-D geometry. For the smaller magnitudes, the orthogonal regression line almost overlaps the 1:1 line, and the scatter is rather symmetrical and in general less pronounced in absolute values. The scatter of all earthquake scenarios together is asymmetrical and clearly dominated by the largest earthquake-magnitude scenario.

D I S C U S S I O N A N D C O N C L U S I O N
In this work, we compared the effects on modelled tsunami amplitudes of two end-member geometric reconstructions of a subduction interface. The first is a realistic 3-D geometry reconstructed from a dense grid of 2-D high-penetration seismic reflection profiles ; the second approximates the first one through a uniformly dipping plane. The two subduction interface models were discretized using triangular elements and the tsunamis generated by three earthquakes of increasing magnitude (M w 7.5, 8.0 and 8.5) were then simulated (Figs 2, 3 and 5), keeping fixed the rupture centre.
Downloaded from https://academic.oup.com/gji/article/223/3/1805/5898673 by guest on 03 March 2022 Figure 6. Maximum (A max ) and minimum (A max ) wave amplitudes calculated for both the 3-D (red) and the planar (black) models of the slab and corresponding relative differences (blue positive and orange negative) on three segments of coast in southern Italy (left), northern Africa (middle) and western Greece (right). Broadside and oblique near-field regimes (Geist 2009) are graphically shown (left-hand column).
Minimizing the differences in the rupture area and keeping the seismic slip uniform, allowed us to explain the subsequent results of the modelling in terms of the geometric construction of the ruptures, which are defined by the depth, strike, dip and rake values of the mesh triangles. On the one hand, this approach has the benefit of isolating the minimum set of parameters onto which the results depend. On the other hand, the main limitation is that the planar fault used to model the ruptures on the subduction interface is, on average, deeper than the complex 3-D fault (Fig. 1e). This discrepancy is due to the difficulty of fitting the complex reality with an ideal simple plane. Given the centre of the three ruptures, the depth range of the planar ruptures results systematically deeper than the 3-D ruptures (Fig. 2). However, a planar shallower fault would likely violate other geological constraints, such as the position of the upper tip within the accretionary wedge or the position of the lower tip relative to the crust bottom (Moho) of the overriding plate. The subjectivity on how the planar geometry is constructed to approximate the 3-D geometry cannot be avoided. In other words, the situation described here epitomizes the general case in which the planar geometry is chosen because of the lack of knowledge and the fact that using such a crude approximation has evident drawbacks. In more specific terms, the smaller ruptures sample fewer mesh triangles, and the strike and, especially, the dip deviations between 3-D and planar geometries are smaller (Fig. 4). The situation is different for the rake and the depth. As regards the rake, sampling fewer triangles for the smaller ruptures also introduces a larger amount of oblique slip for this specific geometrical centre of the ruptures, thereby decreasing the dip-slip component of the rake-parallel seismic slip (Table 1). As regards the depth, the profiles in Fig. 2(b) show that the planar ruptures systematically sample mesh triangles at a deeper depth than the 3-D ruptures. These two aspects, rake and depth, represent a trade-off in terms of seafloor displacement and consequent tsunami initial conditions.
In summary, the modelled planar ruptures are more effective because of the larger dip-slip component, but less effective because of the deeper depth at which the slip takes place (implying a larger volume to deform above the rupture with the same amount of slip). However, a very similar planar geometry of the subduction interface in the Calabrian Arc (dip 5 • , depth 5-20 km) was already adopted to model the source of the 1693 Catania earthquake and tsunami (Gutscher et al. 2006), which is one of the largest documented earthquake (M w 7.3) in this region (Rovida et al. 2020). The comparison with the historical data, therefore, suffers by the same approximations highlighted here.
The results obtained in our case study cannot be quantitatively generalized, possibly because several alternative planar approximations of a real interface could be explored. Each different planar solution would likely better preserve some of the real features at the expense of others. Although the modelled tsunamis strongly depend on the encountered bathymetric features, we could still recognize the imprint of the different geometric models in all of the considered scenarios. Indeed, the purpose of our work is mainly to raise the awareness that great care should be taken when using a planar fault approximation in all types of forward-modelling studies. For example, the modelling of slip complexities in tsunami hazard studies (e.g. Scala et al. 2020, who used the same 3-D geometry we used here) could provide ineffectual results if adopted in conjunction with simplistic planar fault geometries. Among the various alternative approaches that can be used to model fault elastic dislocation (Rudnicki & Wu 1995;Wendt et al. 2009;Geist & Oglesby 2014a;Geist & Oglesby 2014b;Murphy et al. 2016Murphy et al. , 2018Scala et al. 2020), whenever a planar fault is used, we inevitably introduce a bias in the results of our modelling. However, there is no simple rule to anticipate the source and amount of this bias and how the bias propagates in any specific case. For example, uniform slip may create artificially high vertical displacement near the upper edge of the dislocation, which can be avoided by using crack models (Geist & Dmowska 1999). However, this situation occurs at specific critical conditions (shallow source depth and very low dip angles, Geist & Oglesby 2014a) which is not the case for the earthquake scenarios considered in this study. Ultimately, alternative physicsbased approaches also regarding the tsunami generation phase in conjunction with the other geometric and slip distribution factors should be considered (Nosov & Kolesov 2011;Lotto et al. 2017;Lotto et al. 2018).
At the same time, it is important to highlight that it is not always possible to have a large amount of data to model the fault geometry properly; therefore, in these cases using a planar fault model parametrized by using geological constraints remains a viable and maybe necessary procedure. It should be considered though, particularly for hazard studies, to use several alternative representations to explore this epistemic uncertainty.
The variability of the tsunami wave impact along the coasts can also depend, in some cases, on the shape of the coseismically displaced surface that is directly related to the fault geometry. Paradigmatic in this sense might be the M w 7.8, 2006, Java, tsunami earthquake (i.e. an earthquake that generates a regional or teleseismic tsunami greater in amplitude than it would be typically expected from its seismic moment (Polet & Kanamori 2016). In this case, the average value of runup measurements (i.e. the maximum topographic height reached by the tsunami during the inundation stage) collected during post-event surveys was ∼5-7 m, with a significant peak of over 20 m on the south coast of Nusa Kambangan (Fritz et al. 2007). This localized and very large value could be attributed to the convex shape of the subduction zone nearby Java (Kânoglu et al. 2013) that, by analogy with optics, acts as a convex lens that makes the rays (the tsunami waves) converge towards a focal point. A similar situation occurs with the tsunami scenarios presented here for the Calabrian Arc. In our case study, the analysis of the impact was performed onto the offshore points, which are useful to illustrate the different effects obtained by the planar or the 3-D geometry. For the M w 7.5 and 8.0 earthquakes the A max distribution along the coast of southern Italy peaks around POI number (Fig. 6, left-hand panels) for both geometries in the broadside near-field, but larger values (about + 20 per cent) are obtained with the 3-D geometry. In these cases, since the maxima of the tsunami initial conditions are located approximately in the same position for both geometries (Figs 2a and b), the convex shape of the Calabrian Arc is likely the main feature responsible for the tsunami amplification toward POI number 28. The situation is different for the M w 8.5 earthquakes, which show two very different spatial distributions of the initial water displacement (Figs 2a and b), being the maxima of the tsunami initial conditions for the planar geometry closer to the coastline than those for the 3-D geometry. This difference influences the A max distribution determining a peak at POI number 28 for the 3-D geometry but an even larger peak (about + 30 per cent) at POI number 30 for the planar geometry.
To reduce this source of uncertainty, one needs to know beforehand the geometric 3-D reconstruction of any offshore crustal fault and subduction interface. In subduction zones, this include the accretionary wedge internal complexities (Strasser et al. 2009;Tsuji et al. 2014), the faults at their frontal and lateral terminations (Hananto et al. 2020;Maesano et al. 2020), as well as the possible splay faults that if activated can produce significantly different tsunami wave pattern than the subduction interface alone (Wendt et al. 2009). Pursuing this objective would imply a massive use of high-penetration seismic reflection data coupled with a great deal of geological data (e.g. well logs). 3-D seismic surveys would be the ideal tool to capture the details of the complexity of a subduction interface, but they are currently available only for limited portions of subduction interfaces worldwide, such as the Nankai subduction zone (Moore et al. 2007;Gulick et al. 2010). In comparison, 2-D seismic surveys provide less information in the space between the profiles, but they often have a high enough resolution to capture the first order geometrical complexities of megathrusts. They are also less expensive, more diffuse and more easily made available by the Oil and Gas industry, thereby allowing geologists to obtain realistic reconstructions of the subduction interface, as in the case presented here. The fault geometry reconstruction is, therefore, a key element that adds up to the various types of epistemic uncertainty encountered in the seismic source characterization (Basili et al. 2013) for tsunami hazard analyses. However, these types of data are not available everywhere and can be very expensive to collect. Our results, however, show that to overcome these difficulties, the use of alternative planar faults may not suffice. A recommendation for future studies is to adopt geometries derived from well-studied cases in tectonically analogue contexts.
The uncertainty related to the way the fault geometry affects the tsunami impact as well as other observables like ground motion, or geodetic deformation can be, in fact, relevant for seismic and tsunami hazards. This uncertainty also has implications for implementing more effective tsunami warning systems since early evaluations of tsunami alerts need to be supported by forecasting methods which must be as reliable as possible, by catching all the identified sources of uncertainty. In our case study, we explored only the case of a subduction interface. However, the differences we have found for the M w 7.5 scenario, suggest that future studies should be directed at examining the effects of fault geometry reconstructions also for offshore crustal faults, which were modelled as planar for example in the NEAMTHM18 tsunami hazard model (Basili et al. 2019), in the absence of better information regarding their geometry. Where relatively simple major subduction zones dominate the tectonic landscape, crustal faults should generate only smaller earthquakes and tsunamis than the megathrusts, with some exceptions though related, for example, to outer-rise earthquakes (Fujii & Satake 2008), or to supposedly first-order tsunamigenic contribution by splay faults (Waldhauser et al. 2012). However, crustal faults may locally dominate tsunami hazard at places, particularly in relatively complex tectonic contexts, as shown, for example, by Selva et al. (2016), or very likely demonstrated by the 2018 Palu, Sulawesi, tsunami which was indeed initiated in a strike-slip context (Ulrich et al. 2019). The effects of relatively complex crustal fault geometries are also important in ground shaking scenarios (Passone & Mai 2017). Developing studies adopting the 3-D characterization of fault geometries then, possibly including also crustal faults, may prove relevant. Doing so in both fields of ground shaking and tsunami effects can open the floor to more reliable multihazard risk analyses which would require a great deal of effort, and possibly a paradigmatic change on how the earthquake sources are characterized.

A C K N O W L E D G E M E N T S
This work was carried out in the framework of the TSUMAPS-NEAM project co-financed by the European-Union Civil Protection Mechanism (DG ECHO) through the agreement number: ECHO/SUB/2015/718568/PREV26 and benefitted from funding provided by the Italian Presidenza del Consiglio dei Ministri-Dipartimento della Protezione Civile (DPC). This paper does not necessarily represent the ECHO or DPC official opinion and policies. We thank Eric Geist and João Duarte for their precious reviews which substantially improved the quality and clarity of the paper.