Simulation of tsunami generation , propagation and coastal inundation in the Eastern Mediterranean

Introduction Conclusions References Tables Figures


Introduction
The 2004 tsunami in Southeast Asia and its devastating effects brought to the public's attention the long-neglected risk tsunamis pose for coastal areas.The issue had already alerted -to a certain extent -the scientific community (e.g. the review of Dawson et al., 2004 for Europe); however, it is evident that the 2004 event contributed significantly to the rise of awareness in public authorities and policy makers, resulting in a notable shift in related research as well.Accordingly, there has been a continuous effort post-2004 towards the improvement of the tools and methods used for the assessment of coastal vulnerability to tsunami-related hazards, with numerical modeling being the basis of all respective attempts.
Tsunami generation and propagation has been steadily studied since the late 1980s; nevertheless, the main gap in relevant knowledge can be identified as to what happens when tsunami waves approach the nearshore and run inland.The sequence of a tsunami hitting the coast comprises a series of processes: from the tsunami generation and propagation, to coastal-zone hydrodynamics (including surf and swash zone dynamics), coastal inundation and wave-structure interactions with the built environment.Regarding the modeling part -and focusing on coastal inundation -exemplary reference can be made to the work of Borrero et al. (2006), who used the MOST model (Titov and González, 1997) for tsunami generation and inundation in western Sumatra; Gayer et al. (2010), who used the MIKE21 Flow Model FM to simulate inundation based on roughness maps for Indonesia; Omira et al. (2010), who applied a modified version of the COMCOT model (Liu et al., 1998) to selected cases in Casablanca, Morocco; Apotsos et al. (2011), who used the Delft3D model to study inundation and sediment transport by the 2004 SE Asia tsunami in measured and idealized morphologies; and Løvholt et al. (2012), who used models based on the Boussinesq equations for tsunami propagation and nonlinear shallow-water wave equations for coastal inundation to simulate the 2011 Tohoku tsunami.Extending to coastal planning, vulnerability assessment and tsunami hazard mitigation, one may refer to the work of Bernard (2005), Published by Copernicus Publications on behalf of the European Geosciences Union.González et al. (2009), Post et al. (2009), Kumar et al. (2010), Sørensen et al. (2012) and González-Riancho et al. (2014).
Tsunamis in the Eastern Mediterranean have a long and significant history, and have attracted awareness due to the well-established geotectonic regime of the area (i.e.Papadopoulos and Chalkis, 1984;Papazachos and Papazachou, 1998;Soloviev et al., 2000;Papadopoulos, 2003;El-Sayed et al., 2004;Tinti et al., 2004;Papadopoulos and Fokaefs, 2005;Stefatos et al., 2006;Papadopoulos et al., 2014).The Aegean Sea and its surrounding areas, in particular, are not only the most active Mediterranean regions in terms of seismicity and tectonic movements, but their coastlines have also experienced numerous tsunami events in recent, historic and pre-historic times.
Earthquakes and submarine slides are the two principal tsunamigenic mechanisms in the aforementioned region, although volcanic eruption and collapse could not be ignored as a potential mechanism as well (e.g. the Late Minoan Thera event).The generation and propagation of tsunamis in the Eastern Mediterranean has been numerically studied by relatively few researchers, especially in comparison to the geotectonic regime of the area.One may refer to the work of Tinti et al. (2005), for scenarios of tsunamis of tectonic origin from the Algerian earthquake of 1980, the Eastern Sicily Arc and the Western/Eastern Hellenic Arc; Salamon et al. (2007), for tsunamis generated from landslide and/or earthquake scenarios impacting the coasts of Syria, Lebanon and Israel; Lorito et al. (2008) for earthquake-generated tsunamis from the Algeria-Tunisia, Southern Tyrrhenian, and Hellenic Arc source zones; as well as of Yolsal et al. (2007) and Periáñez and Abril (2014), covering all generation mechanisms (geological faults, landslides, entry of pyroclastic flows into the sea and the collapse of a volcano caldera).However, the adequate representation of nearshore dynamics and coastal inundation remains an issue in all relevant attempts for the area.
In the present work, an advanced tsunami generation, propagation and coastal inundation 2-DH model -developed by the authors -is applied to simulate representative earthquake-induced tsunami scenarios in the Eastern Mediterranean.Regarding the coastal hydrodynamics, the nonlinear wave transformation in the surf and swash zone is computed by a nonlinear breaking wave model based on the higher-order Boussinesq equations for breaking and nonbreaking waves (Karambas and Samaras, 2014).Tsunami generation is simulated through additional time derivative terms in the continuity and momentum equations in order to represent displacements at the sea bed or surface.Inundation is simulated based on the dry bed boundary condition (Karambas and Koutitas, 2002); the model's capability in representing swash zone hydrodynamics is validated through the comparison with both two-dimensional (crossshore) and three-dimensional experimental data by Synolakis (1987) and Briggs et al. (1995), respectively.After evaluating tsunamigenic zones and possible sources in the region, two areas of interest were selected for the applications: one at the southwest of the island of Crete in Greece and one at the East of the island of Sicily in Italy.Model results are presented in the form of extreme water elevation maps, sequences of snapshots of water elevation during the propagation of the earthquake-induced tsunamis, and inundation maps of the studied low-lying coastal areas.Regarding the inundation, in particular, this work marks one of the first successful applications of a fully nonlinear model based on the Boussinesq equations for the 2-DH simulation of tsunamiinduced coastal inundation, thus not resorting to estimates of the flooded area from simple superelevations of the water surface or from the spatial extension of cross-sectional runup results.

Boussinesq equations for breaking/non-breaking waves and tsunami generation
Boussinesq-type equations are widely used for the description of the non-linear breaking and non-breaking wave propagation in the nearshore or long wave propagation in the open sea (Gobbi and Kirby, 1999;Gobbi et al., 2000;Ataie-Ashtiani and Najafi Jilani, 2007;Fuhrman and Madsen, 2009;Zhou and Teng, 2009;Zhou et al., 2011).Over the years, the classical Boussinesq equations have been extended so as to be able to include higher-order nonlinear terms, which can describe better the propagation of highly nonlinear waves in the shoaling zone.The linear dispersion characteristics of the equations have been improved as well, in order to describe nonlinear wave propagation from deeper waters (Zou, 1999).Antuono et al. (2009) and Antuono and Brocchini (2013) provide significant improvements with respect to typical Boussinesq-type models for both numerical solution features (Grosso et al., 2010) and the overall flow structures; a thorough overview on Boussinesq-type models can be found in Brocchini (2013).The higher-order Boussinesq-type equations for breaking and non-breaking waves used in this work are the following (Zou, 1999;Karambas and Koutitas, 2002;Karambas and Karathanassi, 2004;Karambas and Samaras, 2014): where M u is defined as: and G as: In Eqs. ( 1)-( 4) the subscript t denotes differentiation with respect to time, d is still the water depth, U is the horizontal velocity vector U = (U, V ) with U and V being the depth-averaged horizontal velocities along the x and y directions, respectively; ζ is the surface elevation, h the total depth (h = d + ζ ), g is the gravitational acceleration, τ b = (τ bx , τ by ) is the bottom friction term (shear stress components approximated by the use of the quadratic law according to Ribberink (1998), δ is the roller thickness (determined geometrically according to Schäffer et al., 1993).E is the eddy viscosity term (according to Chen et al., 2000), and u o is the bottom velocity vector u o = (u o , v o ) with u o and v o being the instantaneous bottom velocities along the x and y directions respectively.In the above, M u is the excess momentum term introduced to account for energy dissipation due to wave breaking; the process itself is based on a specific characteristic of the breaker: the presence of the surface roller, i.e. the passive bulk of water transported with the wave celerity.
Regarding the effects of unresolved small-scale motions, they are parametrized applying the philosophy of the large eddy simulation.The effects of subgrid turbulent processes are taken into account by using the Smagorinsky-type subgrid model (Chen et al., 2000;Zhan et al., 2003).The components of the eddy viscosity term E in Eq. ( 2) are defined as: with the eddy viscosity coefficient ν e estimated from (Zhan et al., 2003): Tsunami generation is simulated through additional terms in the continuity and momentum equations, Eqs. ( 1) and ( 2) respectively.The time derivative term ζ b,t is added to Eq. ( 1) to represent bed level changes (Mitsotakis, 2009), thus transforming it to: where ζ b is the bottom displacement; accordingly the term d 2 ∇ζ b,t is added to the right-hand side of Eq. ( 2).For the bottom displacement function the model follows the approach of Hammack (1973), considering two types of bed movements: an exponential and a half-sine one.The above method is called active tsunami generation.The model additionally includes the option for a passive tsunami generation (i.e. the introduction of the aforementioned displacement directly on the free surface), which is the one used in the present work as well.

Numerical scheme and boundary conditions
The numerical solution of the Boussinesq-type equations (Eqs. 1 and 2) is based on the accurate higher-order numerical scheme of Wei and Kirby (1995), who proposed a fourthorder predictor-corrector scheme for time stepping, discretizing the first-order spatial derivatives to fourth-order accuracy.
The specific discretization has the advantage -over lower order schemes -of automatically eliminating error terms that would be of the same form as the dispersive terms and would, therefore, need to be corrected.The scheme consists of the third-order in time explicit Adams-Bashford predictor step and fourth-order in time implicit Adams-Bashford corrector step (Press et al., 1992;Wei and Kirby, 1995).Energy absorption at the open boundaries is accounted for through the introduction of artificial damping terms in the momentum equation (Eq.2).In particular, terms F and G are added to the right-hand sides of the momentum equation expressions along the x and y directions, respectively; the terms are defined as (Wei and Kirby, 1995): where α r is the constant to be determined for the specific run, and r is the relaxation parameter that varies from 0 to 1 within the specified damping zone (r =1 at the outer edges of the zones and decreasing down to 0 at the edges facing the model domain) according to: with "NN" being the number of grid elements in the damping zone.
The above-described damping layer is applied along with a radiation boundary condition, which for principal wave propagation direction close to the x axis is expressed by the following (Wei and Kirby, 1995): where c 1 = (gd) 1/2 is the phase speed specified by the longwave limit.
The coast in the model can be considered either as a solid (fully or partially reflecting) boundary, or as a boundary allowing sea mass inland penetration and inundation.The first case for a fully reflective boundary derives from the conservative assumption expressed by: where n is the unit landward normal vector; for a partially reflective boundary, it is simulated by properly adjusting the value of the eddy viscosity coefficient ν e (see Eq. 7) in front of the coast.The second case is simulated based on the dry bed boundary condition for the simulation of run-up, as described in detail by Karambas and Koutitas (2002).

Model validation
The model's capability in representing swash zone hydrodynamics was validated through the comparison with both two-dimensional (cross-shore) and three-dimensional experimental data by Synolakis (1987) and Briggs et al. (1995), respectively.Synolakis (1987) studied the run-up and run-down of breaking and non-breaking solitary waves on a plane beach.The experiments were carried out in the wave tank facility of the W. M. Keck Laboratories of the California Institute of Technology.The glass-walled tank's dimensions were 37.73 m × 0.61 m × 0.39 m (length × width × depth); the sloping beach was constructed at a 1:19.85slope (tan a = 1:19.85);the still water depth in the constant depth region was set to 0.2 m.The profile of the solitary wave reproduced, centered at x = X 1 , is given by where H = solitary wave height and γ is defined as Figure 1 shows the comparison of the normalized surface elevation between model results and the measurements of Synolakis (1987) for a solitary wave of H /d = 0.28 amplitude ratio, as a series of snapshots at consecutive nondimensional time instances.Model predictions are in close agreement with the experimental data in both surf and swash zones.Run-up and run-down are simulated well, with the collapse of the bore identified in Fig. 1d, e and f, and the moment of maximum run-up in Fig. 1g.Briggs et al. (1995) studied three-dimensional tsunami run-up on a circular island.The experiments were carried out in the facilities of the US Army Corps of Engineers Waterways Experiment Station (WES) in Vicksburg, Mississippi.The physical model of a conical island was constructed in the center of a 30 m wide 25 m long flat bottom basin, shaped as a truncated right circular cone with a 7.2 m diameter at its toe and a 2.2 m diameter at its crest.The height of the cone was approximately 62.5 cm, with a beach face slope of β = 14 • .Tsunami waves were simulated using solitary waves, their surface profiles given by Eq. ( 14), following the rationale of Synolakis (1987).Experiments for symmetric source lengths and a depth of d = 32 cm in the basin were reproduced, for two different ratios of initial amplitude to depth (i.e.ε = H /d), namely ε = 0.10 and ε = 0.20.Tables 1, 2 and Fig. 2 show the comparison of the normalized maximum run-up height (i.e.run-up height to initial wave height is R/H ) distribution around the circular island between model results and the measurements of Briggs et al. (1995); the angle α = 0 corresponds to the front direction of the wave approaching the island and α = π to the back direction.Figure 3 shows snapshots of the free surface at different time instances for the experiment with ε = 0.20.Model predictions are in close agreement with experimental data for this test as well.For ε = 0.10, R/H values are practically overlapping around the entire island, with a relatively Table 2. Tsunami run-up on a circular island: comparison of the normalized maximum run-up height (R/H ) distribution around the island between model results and the measurements of Briggs et al. (1995) for ε = 0.20; the angle α = 0 corresponds to the front direction of the wave approaching the island and α = π to the back direction.

Tsunamigenic zones and sources
Figure 4 shows a map of the known tsunamigenic zones in the Mediterranean Sea region along with a relative scale of their potential for tsunami generation, calculated as a convolution of the frequency of occurrence and the intensity of tsunami events (Papadopoulos and Fokaefs, 2005;Papadopoulos, 2009).Sakellariou et al. (2007) summarized the possible tsunamigenic sources in the Eastern Mediterranean based on existing marine geological, bathymetric and seismic data; the results are presented in the map of Fig. 5.

Model applications
In the present work, the model for tsunami generation, propagation and coastal inundation presented in Sect. 2 was applied to two areas of interest: one at the southwest of the island of Crete in Greece and one at the East of the island of Sicily in Italy.The areas, indicated in Fig. 5, comprise the sources of the earthquake-induced tsunami scenarios and the low-lying coastal areas where inundation phenomena were studied (see also Fig. 6).
Regarding the earthquake-induced tsunami scenarios, earthquakes that would generate a normalized wave amplitude of ζ 0 = 1 m were considered; the lengths of the major and minor axes of the elliptical aftershock areas were estimated based on the empirical equations proposed by Karakaisis (1984), Papazachos et al. (1986) and Demetracopoulos et al. (1994): where M is the magnitude of the mainshock.It should be noted, of course, that due to the non-linearity of the studied phenomena the presented model results should not be used to estimate the absolute propagation/inundation characteristics of tsunamis with multiple or sub-multiple wave amplitudes at generation for the specific sources/areas of interest.The essence of the presented applications lies in testing the capabilities of the developed model and methodology for real case scenarios of operational interest in the Mediterranean; not in replicating single tsunami events, for which, furthermore, accurate inundation data would not be available.Bathymetric information was extracted by the EMODnet Bathymetry Portal (EMODnet, 2015); shorelines by the GSHHS database (Global Self-consistent Hierarchical Highresolution Geography database; NOAA/NGDC, 2015).The topographic information for the coastal areas of interest were extracted by Digital Elevations Models of the NASA Shuttle Radar Topography Mission, at the best resolution available for the areas of interest (3 arc seconds for Crete and 1 arc second for Sicily; USGS/GDE, 2015).Figure 6 shows the loca-  tion, the elevation and selected topographic contours of the low-lying coastal areas of interest where inundation phenomena were studied, at south-southwest Crete (Fig. 6a) and eastsoutheast Sicily (Fig. 6b).

Results and discussion
Figure 7 shows the simulated extreme water elevation (ζ /ζ 0 ) for the earthquake-induced tsunami scenarios at the southwest of Crete (Fig. 7a) and at the East of Sicily (Fig. 7b).Sequences of snapshots of water elevation during tsunami propagation are presented in Figs. 8 and 9 for the two aforementioned scenarios, respectively.A tsunami generated at the southwest of Crete (see Figs. 7a and 8) would impact most severely the adjacent coasts (as expected due to the source proximity), with calculated extreme water elevations reaching the normalized amplitude of the tsunami wave at generation.The impact is expected to be significant to the East part of the Libyan coast as well (approx.250 km away from the tsunami source; at x ≈ 650 km ÷ 750 km), with extreme elevations locally exceeding ζ /ζ 0 = 0.4.A tsunami generated at the East of Sicily (see Figs. 7b and 9) would have a similar Figure 10 shows the inundation maps of the studied lowlying coastal areas at (a) south-southwest Crete and (b) eastsoutheast Sicily.The inundated areas shown in Fig. 10a and  b, where the inundation extent was more easily represented at the scale of interest, cover 3.429 and 0.641 km 2 , respectively.Although for both cases inundation heights are comparable, the relatively steeper slopes at the studied coasts of Sicily result in an overall narrower inundation zone (with the inevitable added scale effect of the representation).Again, it should be noted that these areas are indicative of the ones to be affected by eventual larger events.Finally, regarding the simulation of inundation itself, Fig. 11 shows snapshots of the evolution of the phenomenon at the coasts of south-  southwest Crete for an exemplary exaggerated tsunami scenario (normalized wave amplitude of ζ 0 = 6 m at generation); it should be underlined that these results serve only to demonstrate more clearly the performance of the model at the scale of the representation, and do not relate to the scenarios presented in the previous.

Conclusions
This work presents an advanced tsunami generation, propagation and coastal inundation 2-DH model (developed by the authors) and its applications for two representative earthquake-induced tsunami scenarios in the Eastern Mediterranean.The model is based on the higher-order Boussinesq equations, and its capability in representing swash zone hydrodynamics is validated through the comparison with both two-dimensional (cross-shore) and threedimensional experimental data by Synolakis (1987) and Briggs et al. (1995), respectively.The model is applied to two areas of interest: one at the southwest of the island of Crete in Greece and one at the East of the island of Sicily in Italy.Model results, presented in the form of extreme water elevation maps, sequences of snapshots of water elevation during the propagation of the earthquake-induced tsunamis, and inundation maps of the studied low-lying coastal areas, highlight the model's capabilities and are indicative of how areas in the region would be affected by eventual larger events.It should be noted that this work marks one of the first suc-cessful applications of a fully nonlinear model based on the Boussinesq equations for the 2-DH simulation of tsunamiinduced coastal inundation, thus not resorting to estimates of the flooded area from simple superelevations of the water surface or from the spatial extension of cross-sectional run-up results.Similar attempts can constitute the basis of a more detailed coastal flooding risk assessment and mitigation along the coasts of the Eastern Mediterranean.

Figure 1 .
Figure 1.Run-up and run-down of a solitary wave of H /d = 0.25 on a 1:19.85plane sloping beach; comparison of the normalized surface elevation (ζ /d) between model results and the measurements of Synolakis (1987), at consecutive non-dimensional time instances (t = t (g/d) 1/2 ).

Figure 2 .
Figure 2. Tsunami run-up on a circular island: comparison of the normalized maximum run-up height (R/H ) distribution around the island between model results and the measurements of Briggs et al. (1995) for (a) ε = 0.10 and (b) ε = 0.20; the angle α = 0 corresponds to the front direction of the wave approaching the island and α = π to the back direction.

Figure 3 .
Figure 3. Tsunami run-up on a circular island: snapshots of the free surface for the experiment with ε = 0.20 at: (a) t = 6.0 s, (b) t = 7.5 s, (c) t = 9.0 s, and (d) t = 12.0 s.

Figure 4 .
Figure 4.The tsunamigenic zones of the Mediterranean Sea and their respective tsunami potential (adopted from Papadopoulos and Fokaefs, 2005).

Figure 5 .Figure 6 .
Figure 5. Possible tsunamigenic sources in the Eastern Mediterranean; the black rectangles outline the two areas of interest in the present work, comprising the sources of the earthquake-induced tsunami scenarios and the low-lying coastal areas where inundation phenomena were studied (adopted from Sakellariou et al., 2007; privately processed).

Figure 7 .
Figure 7. Simulated extreme water elevation (ζ /ζ 0 ) for the earthquake-induced tsunami scenarios at: (a) the southwest of Crete, and (b) the East of Sicily.

Figure 11 .
Figure 11.Snapshots of the evolution of inundation at the coasts of southwest Crete (see also Fig. 6a) for an exemplary exaggerated tsunami scenario (normalized wave amplitude of ζ 0 = 6 m at generation), presenting: (a) the propagation of the tsunami wave to the nearshore; (b) the wave breaking at the lower part and the run-up at the upper part of the figure; (c) the run-up at the lower part and run-down at the upper part of the figure; and (d) the run-down at the lower part of the figure.

Table 1 .
Briggs et al. (1995)circular island: comparison of the normalized maximum run-up height (R/H ) distribution around the island between model results and the measurements ofBriggs et al. (1995)for ε = 0.10; the angle α = 0 corresponds to the front direction of the wave approaching the island and α = π to the back direction.