Introduction

Offshore Crete (Greece) at the Eastern Mediterranean Basin (EMB) was the site of a small tsunami following an Mw 6.6 earthquake on 2 May 2020 (Fig. 1). The earthquake occurred at 12:51:06 UTC with an epicenter located at 25.712 oE and 34.205 oN according to the United States Geological Survey (USGS). The USGS W-phase moment-tensor focal mechanism solution resulted in strike angle, 229°; dip angle 31°; rake angle, 46° and depth, 11.5 km (Fig. 1). Table 1 gives a summary of seismic parameters of the earthquake from different seismological agencies. The earthquake was felt on the south of Crete Island but made no damage or injuries according to media reports. The tsunami was small and made no damage although it was observed on the south coast of the Crete Island. Papadopoulos et al. (2020) studied this tsunami and earthquake and provided a discussion on the responses of tsunami service providers in the region.

Fig. 1
figure 1

The map of the East Mediterranean Basin showing the epicenter of the 2 May 2020 earthquake (blue and red stars), the tide gauge stations used in this research (orange triangles) and the Tsunami Travel Time (TTT) contours in hours with 0.25 h (= 15 min) intervals (dashed lines). Some other tsunamigenic earthquakes are also shown (pink stars). References: USGS, GCMT, Shaw et al. (2008), Papadopoulos et al. (2007), Yolsal-Çevikbilen and Taymaz (2012)

Full size image
Table 1 Seismic parameters of the 2 May 2020 Crete Island earthquake from different seismological agencies
Full size table

Although the May 2020 event was small, it occurred in an important tectonic setting within the Mediterranean Sea; thus providing important information on some aspects of the region’s seismotectonics. This thrust-dominant mechanism earthquake took place along the Hellenic Subduction Zone (HSZ) where the African Plate subducts beneath the Eurasian Plate (Fig. 1) at the rate of 25–60 mm/yr based on various estimates (Taymaz et al. 1990; Kassaras et al. 2020; Reilinger et al. 2010). The HSZ is the most active seismic zone in the EMB. The slip along the HSZ was long considered as occurring aseismically (e.g. Papazachos and Kiratzi 1996; Shaw et al. 2008); however, Laigle et al. (2002), demonstrated that a complete seismic coupling can be seen at the HSZ and that the western HSZ is the most active seismic zone in Europe. The largest known event in the HSZ occurred in July 365 at the west of Crete Island (Fig. 1). The study by Shaw et al. (2008) concluded that a subduction-zone thrust-fault earthquake (Mw 8.3 – 8.5) was responsible for the event of July 365. Another major earthquake and tsunami in the HSZ was reported in 1303 at its eastern end (Fig. 1) whose magnitude was estimated at ~ M7.8 – 8.0 (Ambraseys et al. 1994; El-Sayed et al. 2000; Hamouda 2006; Salamon et al. 2007; Papadopoulos et al. 2007; Papadimitriou et al., 2016). Recently, on 12 October 2013, an Mw 6.6–6.7 thrust earthquake occurred along the western HSZ (Fig. 1) although no noticeable tsunami was reported afterwards (Vallianatos et al. 2014; Papadimitriou et al. 2016). The EMB region also experienced two moderate tsunamis following an Mw6.6 earthquake in July 2017 and an Mw7.0 in October 2020 on the Turkey/Greece border (Fig. 1) which made some damage to coastal areas in both countries (Karasözen et al. 2018; Heidarzadeh et al. 2017). Although the July 2017 and October 2020 events are not directly linked to HSZ’s seismicity, they are related to the complex extensional pattern in the back-arc region (e.g. Ring et al. 2010).

The objective of this research is to improve the current understanding of earthquake occurrence and tsunami behavior from thrust earthquakes originating in the HSZ through modelling the May 2020 event. The above short review reveals that the HSZ is comparatively less understood in terms of earthquake and tsunami generation, mostly due to the limited number of instrumentally-recorded events as well as long return periods of megathrust events. Therefore, every event, in particular those associated with tsunamis, is of great regional importance. We analyze tsunami records of this event, conduct tsunami inversion to obtain a slip model for the earthquake, and perform numerical modeling of the tsunami to better understand tsunami characteristics in the semi-enclosed and narrow EMB.

Data and methods

Sea level data

The data used in this study comprises sea level records from tide gauge stations in the region as well as bathymetry grids. The sea level data of three tide gauge stations at NOA-03 (Kasos port, Greece), NOA-04 (Ierapetra port, Greece) and Alexandria (Egypt) are used; all have sampling intervals of 60 s (Fig. 2). National Observatory of Athens (NOA), which is a tsunami warning center, supplied the data of NOA-03 and NOA-04 while Alexandria’s record was obtained through the sea level station monitoring facility of the UN Intergovernmental Oceanographic Commission (UN/IOC) (www.ioc-sealevelmonitoring.org). The two stations of NOA-04 and NOA-03 are located at distances ~ 100 km and ~ 180 km, respectively, from the epicenter while Alexandria is located ~ 500 km away from the epicenter. Maximum tsunami zero-to-crest amplitudes are measured 15.2 cm at NOA-04 and 6.5 cm at NOA-03. In Alexandria, the tsunami signal is hidden within the noise level and thus is undetectable (Fig. 2). It can be seen that the tsunami period (the peak-to-peak interval) is 4.0–4.5 min based on visual inspections of the waveforms (Fig. 2b).

Fig. 2
figure 2

Tide gauge records of the 2 May 2020 tsunami occurring offshore the Crete Island. a Original tide gauge records. b De-tided tsunami waveforms. The vertical red dashed lines indicate the origin time of the earthquake

Full size image

Bathymetry data

It is essential to provide high-resolution bathymetry data to accurately reconstruct the tsunami event and to develop a reliable numerical model. Therefore, we constructed high-resolution bathymetry grids with a spatial resolution of 0.1 arc-sec × 0.1 arc-sec (equivalent to 3 m × 3 m) for the port areas around the NOA-03 and NOA-04 tide gauge stations (Fig. 3) as such data were unavailable for the region. We note that interpolation of low-resolution bathymetry data is not helpful towards making a reliable model (e.g. Synolakis and Bernard 2006; Papadopoulos and Fokaefs 2005). We digitized existing nautical charts for both ports followed by converting the digitized bathymetry points into a uniform grid using the open-source software Generic Mapping Tools (GMT) (Wessel and Smith 1998). Figure 3 compares our newly-constructed high-resolution bathymetry grids with the existing low-resolution (15 arc-sec) grids of the General Bathymetric Chart of the Ocean (GEBCO) (Weatherall et al. 2015). It can be seen that our new high-resolution bathymetry grids appropriately incorporate all small-scale coastal features such as the 5-m-wide breakwaters and variations of water depths inside the two ports (Fig. 3).

Fig. 3
figure 3

Comparison of low-resolution (15 arc-sec) bathymetry data from GEBCO (left panels) and high-resolution (0.1 arc-sec; right panels) bathymetry data constructed in this study for areas around the two tide gauge stations: a NOA-04 (Ierapetra port), and b NOA-03 (Kasos port). The term “GEBCO bathy” indicates the 15 arc-sec bathymetry data from the General Bathymetric Chart of the Ocean (GEBCO) (Weatherall et al. 2015)

Full size image

Tsunami simulations

For tsunami simulations, a nested grid system was employed to obtain a balance between model precision and computational cost, in which the spatial resolutions of the grids increase as we move from the source region towards the coastal areas (Fig. 4). The nested grid system comprised four levels of A, B, C, and D with respective spatial resolutions of 30, 10, 3.33, and 1.11 arc-sec (Fig. 4). The first three layers (A, B, and C) were resampled from the 15 arc-sec bathymetry data of GEBCO. The high-resolution grids (0.1 arc-sec), constructed in this study, were resampled to 1.11 arc-sec and were used as the fourth layer (i.e. D) to model the computational domains around NOA-03 and NOA-04 tide gauges. Linear and nonlinear Shallow Water equations were applied for tsunami inversion and forward simulations, respectively. Tsunami travel time (TTT) analysis was carried out using the TTT software of Geoware (2011).

Fig. 4
figure 4

The nested grid system used in this study around the two tide gauge stations NOA-03 and NOA-4 comprising four grids: A (resolution: 30 arc-sec), B (resolution: 10 arc-sec), C (resolution: 3.33 arc-sec) and D (resolution: 1.11 arc-sec). The four bathymetry grids of the nested system for areas around tide gauge NOA-04 are marked as Grid A, NOA-04B, NOA-04C, and NOA-04D

Full size image

Spectral analysis (Fourier and Wavelet)

Sea level data underwent quality control which helped to remove spikes, gaps and outliers from the original data. The tidal analysis package TIDALFIT (Grinsted, 2008) was applied to calculate tidal signals which were then removed from the original records (Fig. 2, left) to obtain the de-tided tsunami waveforms (Fig. 2, right) (e.g. Heidarzadeh et al. 2018, 2020). For Fourier analysis of the tsunami waveforms, the Welch (1967) algorithm was applied by considering Hanning windows and 30% overlaps following the methodology previously adopted in a number of studies (e.g., Rabinovich 1997; Heidarzadeh and Satake, 2013, 2015, 2017). The length of data used for spectral analysis was 450 min for both tsunami and background signals, which corresponds to 450 data points as our data sampling interval was 1 min. Wavelet analysis was performed following the algorithm by Torrence and Compo (1998). Several authors have applied Wavelet analysis to tsunami waves (e.g. Rabinovich and Eblé 2015; Allgeyer et al. 2013; Heidarzadeh et al. 2015).

We calculate tsunami spectral ratio which is the result of dividing tsunami spectra to that of the background sea level signals at each station. It was shown that spectral ratio is an efficient way to separate tsunami source periods from those generated by bathymetric features during tsunami propagation (i.e. non-tsunami periods) (Rabinovich 1997; Vich and Monserrat 2009; Heidarzadeh and Satake 2014a, 2017). Large spectral ratios are expected for the tsunami dominant periods, while a ratio of approximately one is achieved for non-tsunami periods (Heidarzadeh et al. 2018). Based on the theoretical modeling of Rabinovich (1997), tsunami spectral ratios are independent of instrument location and represent only tsunami forcing. As tsunamis may trigger some local modes, the spectral ratio plots likely include some local modes as well. The peak periods of the spectral ratio plots are considered as the dominant tsunami source periods.

Tsunami inversion and validation by checkerboard tests

For obtaining a source model for the earthquake, we employed the tsunami inversion method of Satake (1987, 2014) which has been applied to numerous tsunamis worldwide (e.g. Satake et al. 2013; Gusman et al. 2015, 2016a, b). Tsunami inversion is based on splitting the source of an earthquake into a number of sub-faults (\(j=1, 2, 3,\dots , N\)) and calculating tsunami Green’s functions \(\left\{{G}_{ij}(t)\right\}\) at the location of observation points (e.g. tide gauge stations) (\(i=1, 2, 3, \dots , M\)) by assuming a unit slip in a sub-fault. For example, \({G}_{31}(t)\) represents the Green’s function at tide gauge number three resulting from a unit slip at sub-fault number one. The simulated waveform at each tide gauge station \(\left\{{Z}_{sim\_i}(t)\right\}\) can be considered as a linear combination of all Green’s functions at that particular station with an amplifying factor (\({d}_{j}\)) for each Green’s function; in fact, \({d}_{j}\) represents the slip in each sub-fault “\(j\)”:

$${Z}_{sim\_i}\left(t\right)= \sum_{j=1}^{N}{G}_{ij}\left(t\right)\times {d}_{j}$$
(1)

The final value of slip in each sub-fault (\({d}_{j}\)) can be obtained by minimizing the difference between simulated \(\left\{{Z}_{sim\_i}(t)\right\}\) and observed \(\left\{{Z}_{obs\_i}(t)\right\}\) waveforms in all stations by solving the following minimization problem:

$$min: \sum_{i=1}^{M}\sum_{j=1}^{N}{\left[{Z}_{sim\_i}\left(t\right) - {Z}_{obs\_i}\left(t\right)\right]}^{2}$$
(2)

To solve Eq. (2), it is evident that the coefficients \({d}_{j}\) cannot take negative values as they represent slip values on the sub-faults. We use the non-negative least-squares solver “lsqnonneg” imbedded in the optimization toolbox of MATLAB software (Mathworks 2020) to solve Eq. (2).

The Green’s functions were built using a tsunami numerical model that solves the linear shallow water equations in a spherical coordinate system (Gusman et al. 2010). The time step for the numerical tsunami simulation using the nested grid system was 0.25 s to satisfy the Courant-Friedrichs-Lewy condition (Courant et al. 1928) in the finest grid (i.e. grid level D). The initial seafloor displacement field due to unit-slip sources in each sub-fault was calculated using the Okada’s (1985) formula and considering fault parameters: strike, dip, rake, depth, sub-fault dimensions, and slip amount. Initial tsunami simulation results using single fault models show that the fault geometry from the GCMT solution gives waveforms that fit the tsunami data better than those from USGS solution. The simulation results from the GFZ solution are similar to those from the GCMT solution. Therefore, we applied the GCMT focal mechanism for our tsunami inversion (i.e. strike, 257°; dip, 24°; rake 71°). A rectangle fault with dimensions of 40 km (along strike) × 30 km (along dip) was subdivided into twelve sub-faults with each having a size of 10 km × 10 km. The depth of the shallowest edge of the fault was taken as 8 km.

Checkerboard tests were conducted to validate our tsunami inversion process. Slip amounts of 0.5 and 0.25 m were used to generate a target model with a checkerboard pattern (Fig. 5a). We evaluated the resolvability of the target model from synthetic tsunami waveforms at the NOA-04 and NOA-3 stations. The length of synthetic tsunami waveforms used in the inversion was 15 min. The synthetic tsunami waveforms were corrupted by adding Gaussian random noise. Our tsunami inversion code was then applied along with the target synthetic waveforms from the target model to produce an inverted slip model. Our results show that the inverted model moderately well reproduces the target model (Fig. 5b).

Fig. 5
figure 5

Checkerboard test result. a Target slip model. b Inverted slip model using waveforms generated by the target model. Slip amounts of 0.25 and 0.5 m are used for the target model

Full size image

Spectral analyses and estimating source dimension

Spectra for the tsunami waveforms (Fig. 6b, blue lines) at two stations of NOA-03 and NOA-04 are compared with the background spectra (black lines). Here, “background” refers to part of the waveform before tsunami arrival at an observation station. A clear increase of spectral energy is seen at the period band of approximately 3–5 min (Fig. 6b) which is attributed to the presence of tsunami signals in this period band. Wavelet analysis (Fig. 7) reveals two major period bands of 3–5 min and 5–10 min with the former band (i.e. 3–5 min) dominating during the first hour since tsunami generation and thus can be attributed to the tsunami source. The later period band of 5–10 min is most likely a local/regional mode triggered by the tsunami. Both tsunami and background spectra match at very low frequencies (period > 20 min; Fig. 6b) which implies there is no tsunami energy at this band. Results of spectral ratio analysis (Fig. 6c) reveal that the peaks of the spectral ratio plots are 3.8 min and 4.7 min for NOA-04 and 4.2 min for NOA-03. This suggests that the dominant tsunami periods must be in the range of 3.8–4.7 min.

Fig. 6
figure 6

Spectral analysis of the sea-level records of the 2 May 2020 tsunami occurring offshore the Crete Island. a Sea level records. b Spectra of the tsunami parts of the sea level records and the spectra of the background parts of the records before the tsunami arrival in each station. c The spectral ratio which is the division of the tsunami spectrum by the background one

Full size image
Fig. 7
figure 7

Wavelet analysis of the sea level records of the 2 May 2020 tsunami occurring offshore Crete Island at station NOA-03

Full size image

The dominant period band of 3.8–4.7 min for the May 2020 event appears to be shorter than that for other tsunamis; for example, the tsunami dominant period was 7–13 min for the 2017 Bodrum/Kos event from an Mw 6.6 earthquake (Heidarzadeh et al., 2017) and was 10–22 min for the 2013 Santa Cruz tsunami from an Mw 8.0 earthquake (Heidarzadeh et al. 2016). Although the May 2020 offshore Crete and the 2017 Bodrum/Kos tsunamis were both generated by a same-size earthquake (Mw 6.6), the May 2020 event has a shorter wave period.

Dimension of a tsunami source is correlated with its dominant period: the larger the source dimension, the longer the period of the tsunami will be. Earthquake fault length can be estimated from the dominant period of the tsunami (\(T\)) through an inverse analysis. In a forward analysis, an earthquake ruptures a length \(L\) of the seafloor followed by a tsunami with a dominant period of \(T\). Tide gauges at the coastal areas record the tsunami and reveal its dominant period. In an inverse analysis, we have \(T\) through Fourier analysis of the tide gauge records and we conduct an inverse analysis to estimate the length of the earthquake rupture zone (\(L\)) around the source area. Heidarzadeh and Satake (2014b) developed an equation for estimating the tsunami characteristic length (\(L\)) from its dominant period (\(T\)):

$$L= \frac{T}{2} \sqrt{gh}$$
(3)

where, \(g\) is gravitational acceleration (i.e. 9.81 m/s2) and \(h\) is water depth around the source region (i.e. epicenter). Application of Eq. (3) to the May 2020 event gives an earthquake fault length of 18–27 km by assuming a water depth of 2500–3500 m around the epicenter (Fig. 8b) and using a tsunami period of 3.8–4.7 min. We note that it would be useful having tsunami measurements around the source region, which is unavailable for this event. However, our tsunami spectral analysis (Figs. 6–7) is able to reveal tsunami source periods using coastal observations. The tsunami characteristic length (\(L\) ) estimated by Eq. (3) is a rough estimate of the source length and cannot be interpreted as a definite value. In addition, the tsunami period dictated by the fault length and fault width could be different as reported by Heidarzadeh and Satake (2013). For this event, such an effect appears to be weak due to the relatively small size of the event.

Fig. 8
figure 8

Estimating the source dimension of the 2 May 2020 tsunami at offshore the Crete Island. a Bathymetry contours around the source region. The focal mechanism is based on the USGS. b Variations of the tsunami source dimension (\(L\)) based on tsunami dominant period (\(T\)) and the water depth (\(h\))

Full size image

Alternatively, the earthquake fault length can be estimated given the moment magnitude of the earthquake (\({M}_{w}\)) using the empirical relationship proposed by Wells and Coppersmith (1994) for reverse (i.e. thrust) earthquakes:

$${\mathrm{log}}_{10}L= -2.8+0.63 {M}_{w}$$
(4)

where \({\mathrm{log}}_{10}\) is logarithm with base 10. Considering a moment magnitude between 6.6 (USGS) and 6.7 (our analysis; next Section) for the May 2020 earthquake, Eq. (4) yields a source length of 20 – 26 km. It can be seen that the source length estimates made using the two aforementioned methods are close to each other. It is noteworthy that the source length estimated here for the Crete event (\(L\) = 18–27 km) represents the main slip area (or main deformation area) rather than the total length of the source.

Source model based on tsunami inversion

Initial tsunami simulations of different single fault models showed that the synthetic waveform at NOA-03 always arrives earlier than the observation. This type of time mismatch is often found for tide gauge records and thought to be the result of an inaccurate bathymetric model (Romano et al. 2016). An optimum time alignment procedure in tsunami waveform inversion was proposed by Romano et al. (2016) to reduce waveform misfit which can be obtained simultaneously with the source model. In this study, we manually find an optimum time alignment of 1 min for station NOA-03. The overall waveform misfit using the optimum time shift is only slightly better than the one without time shift because the tsunami amplitude at NOA-04 is much larger than the one at NOA-03. Moreover, the amplitudes of the first two wave cycles at NOA-03 are below the noise level. These suggest that the inversion result is mostly controlled by the tsunami waveform at NOA-04.

Guided by the results of spectral analysis and estimation of the source dimension, an initial fault with a size of 40 km × 30 km was considered comprising 12 sub-faults (Fig. 9a). The slip distribution estimated from the tsunami waveform inversion has a main slip region at the depth of about 12 km and its maximum slip amount is 0.64 m (Table 2 and Fig. 9a). The average slip, considering only sub-faults with non-zero slip values, is 0.28 m. By assuming an earth rigidity of 4 × 1010 N/m2, the seismic moment calculated from our slip distribution is 1.13 × 1019 Nm which is equivalent to a moment magnitude of Mw 6.67. Our calculated seismic moment is close to the GCMT seismic moment of 1.01 × 1019 Nm. Several inversion attempts were made with different sub-fault locations to obtain the best-estimated slip distribution for the earthquake. We initially used the USGS epicenter as the center of the fault which resulted in a major-slip area located on the southern sub-faults (Figure S1). The GCMT centroid location (latitude = 34.06° E, longitude = 25.63° N; Table 1) was also assumed to be the center of the fault; this time, the major-slip area we moved to the eastern sub-faults (Figure S2). From this experiment, we concluded that the center of the fault should be located between the USGS epicenter and the GCMT centroid.

Fig. 9
figure 9

a Source model (slip distribution) of the 2 May 2020 tsunami offshore Crete Island from tsunami waveform inversion. Blue star represents the epicenter; open circles indicate the aftershocks from 2 to 27 May 2020 based on the USGS catalogue. b The calculated seafloor displacement from the source model. c Tsunami waveform comparison between observation (black) and simulation (red) at NOA-04 and NOA-03 tide gauges. The blue lines indicate the observed waveforms that were used for the inversion

Full size image
Table 2 Fault parameters for the
Full size table

The estimated slip distribution (Fig. 9a) gives calculated maximum uplift of ~ 20 cm and maximum subsidence of ~ 4 cm on the seafloor (Fig. 9b). The simulated tsunami waveform at NOA-04 matches the observation very well (Fig. 9c). Although the simulated first tsunami wave peak slightly overestimates the observation in NOA-04, the second and the third wave peaks match the observation very well (Fig. 9c). The amplitudes of the later waves are also very close to the observation. We note that for our final source model, we used the first two wave cycles of the recorded tsunami at NOA-04 for the inversion process (Fig. 9c, blue wave), because, by inverting only the first tsunami wave cycle (Additional file 1: Figure S3), the simulated later wave amplitudes were underestimated and were out of phase (Additional file 1Figure S3). Therefore, it was necessary to use the first two wave cycles in the inversion process to obtain the best results. The time windows used for the inversion of NOA-04 and NOA-03 records were 13 min and 10 min, respectively.

Coastal hazards from tsunami reflections and basin oscillations in the EMB

The observed tsunami waveform at NOA-04 reveals at least three distinct wave trains (Fig. 10). After each of the first and the second wave trains, an approximately 10–15 min time lapses of quiescence can be seen. The first and the second wave trains carry waves with visual periods of 3–5 min, which is the source period of the tsunami, while the third train comes with a significantly-different period of 5–10 min (Fig. 10). These three wave trains and their timings are also clear in the Wavelet plot (Fig. 7). The first wave train is directly from the tsunami source region (i.e. epicentral area) while the second train appears to be a reflected wave (Fig. 11) which still carries the tsunami source period (i.e. 3–5 min). The reason for considering that the second wave train is most likely a reflected one and not from harbor resonance are twofold: I) the harbor resonance waves usually have different periods and are made of shorter-period waves (e.g. < 1 min); II) harbor resonance usually form monochromatic and constantly-amplifying waves which last at least a few hours; both of the aforementioned conditions cannot be seen for the second wave train in the NOA-04 record (Fig. 10).

Fig. 10
figure 10

Different wave trains and associated visual wave periods (\(T\)) observed on the tide gauge record of NOA-04 (Ierapetra port) during the 2 May 2020 offshore Crete tsunami

Full size image
Fig. 11
figure 11

Snapshots of the propagation of the 2 May 2020 offshore Crete tsunami at different times showing reflected waves from the Crete Island (black dashed lines) and from the North African coast (blue dashed lines). The times indicated at the bottom-left corners of the panels indicate the times after the earthquake origin time

Full size image

It is most likely that the third wave train, whose wave period is double the tsunami source period, is a component of basin oscillation. The durations of the first and the second wave trains are approximately 10 min and 20 min, respectively, whereas the third train lasts for hours (Figs. 2, 10). A basin such as the semi-enclosed region between the Crete Island and the North African coast, with a number of islands in between, normally possesses several oscillation modes; some of which could be excited by a tsunami or a storm depending on the size and location of the event (Raichlen and Lee 1991; Rabinovich 1997; Synolakis 2003; Yalciner and Pelinovsky 2007; Heidarzadeh and Satake 2014a). It is beyond the scope of this study to determine various oscillation modes of this basin as it requires separate investigations. However, at least we can say that the third train is most likely a basin mode because: a) it lasts for longer time, and b) it was initiated approximately 90 min after the earthquake origin time (Fig. 10) which implies that the tsunami had enough time to complete at least one return journey between Crete and North African coast (see TTT in Fig. 1) and set all of the sea surface in motion which is a prerequisite for triggering a basin mode.

Effects of bathymetry resolution and nonlinearity on tsunami hazard analysis

To evaluate the effects of the spatial resolution of the bathymetry grids on the accuracy of tsunami simulations, we compare tsunami simulation results using low-resolution GEBCO data (15 arc-sec) with those using our high-resolution grids (i.e. the nested grid system including grid D with resolution of 1.11 arc-sec; Fig. 4). The source model for both simulations was the one obtained by our tsunami waveform inversion (Fig. 9a). On the low-resolution grid, the NOA-03 tide gauge is located on a dry cell; hence, we shifted its location to the nearest wet cell (Fig. 12b and d). Figure 12f demonstrates that the simulations using a low-resolution grid are significantly underestimating the tsunami at the NOA-04 station (Fig. 12f, green) while that using high-resolution bathymetry reproduces the tsunami very well at this station (Fig. 12f, red). At NOA-03, the difference in tsunami simulations between low- and high-resolution bathymetries is not significant for the first few cycles; however, after that, the low-resolution bathymetry significantly underestimates the tsunami. The maximum tsunami amplitude simulated using the low-resolution grid (Fig. 12c, d) is approximately twice smaller than the one obtained using the high-resolution grids (Fig. 12a, b). This experiment suggests that it is essential to employ high-resolution bathymetry data in order to obtain accurate tsunami simulation results and to deliver reliable tsunami hazard analysis.

Fig. 12
figure 12

Effects of the spatial resolution of bathymetry grids on the propagation and coastal amplification of the May 2020 offshore Crete tsunami in two locations of NOA-04 (Ierapetra port) and NOA-03 (Kasos port). The simulated maximum tsunami amplitude in the smallest modeling domain around NOA-04 using: a high-resolution, and c low-resolution bathymetry data. The simulated maximum tsunami amplitude in the smallest modeling domain around NOA-03 using: b high-resolution, and d low-resolution bathymetry data. e) The simulated maximum tsunami amplitude in the largest computational domain. F Simulated (red and green) and observed (black) tsunami waveforms using the high- and low-resolution bathymetry data at NOA-04 and NOA-03 tide gauges. Red triangles indicate positions of tide gauges NOA-04 and NOA-03

Full size image

We also simulated the tsunami from the estimated source model by solving nonlinear shallow water equations. Figure 13 shows that the first two wave cycles at NOA-04 from linear and nonlinear simulations are almost the same. This reconfirms that the use of the two wave cycles at NOA-04 in the tsunami waveforms inversion is indeed accurate. The difference between the two simulation results is more noticeable at NOA-03 although it is negligible. However, this does not change the reliability of our final source model, because the tsunami amplitude at NOA-03 is much smaller than that at NOA-04 and the tsunami wave used for the inversion at NOA-03 is smaller than the noise level. Another effect that might be important in tsunami source inversion is the dispersion effect. It has long been known that dispersion is significant for tsunamis with short wavelengths (500 m – 2000 m) such as those generated by underwater landslides (e.g. Synolakis, 2003; Heidarzadeh et al. 2020). As the source of the 2020 Crete tectonic tsunami is in the range of 30 – 40 km (Fig. 9), it is believed that the dispersion effect is negligible.

Fig. 13
figure 13

Effects of nonlinearity on tsunami propagation for the May 2020 offshore Crete tsunami. Tsunami waveforms simulated by solving the linear (red) and nonlinear (blue) shallow water equations

Full size image

Conclusions

We studied the 2 May 2020 offshore Crete tsunami generated by an Mw 6.6 earthquake using sea level records and applying spectral analysis, numerical modeling and tsunami source inversion. For accurate tsunami simulations, we constructed high-resolution (0.1 arc-sec) bathymetry grids for the Ierapetra and Kasos ports in which the two tide gauge stations NOA-04 and NOA-03 were installed, respectively. The main results are:

  1. (i)

    The tsunami signals were clear at the NOA-04 and NOA-03 tide gauge stations with respective maximum zero-to-crest amplitudes of 15.2 cm and 6.5 cm. Tsunami was not detectable on the Alexandria tide gauge station located on the North African coast at the distance of 510 km from the epicenter.

  2. (ii)

    Spectral analysis revealed that the dominant period band of the May 2020 tsunami was 3.8 – 4.7 min which is comparatively shorter than the period band of other tsunamis generated by the same-size earthquakes (i.e. Mw 6.6). This is attributed to the relatively deep water depth (approximately 3000 m) around the source area.

  3. (iii)

    Our tsunami source inversion resulted in a heterogeneous fault model with maximum slip of 0.64 m and an average slip of 0.28 m. This source model gives a seismic moment of 1.13 × 1019 Nm for the earthquake which is equivalent to a moment magnitude of Mw 6.67.

  4. (iv)

    At least three distinct wave trains were seen on the wave record of NOA-04, which is the nearest tide gauge station to the epicenter. The first and the second wave trains carry waves with periods of 4 – 5 min, which is the source period of the tsunami, while the third train comes with a significantly-different period of 5 – 10 min. The first wave train is directly from the tsunami source; and the second train most likely is a reflected wave which still carries the tsunami source period. It is possible that the third wave train is the result of basin oscillations because it starts around 90 min after the earthquake and lasts for hours.