Long Tsunami Oscillations Following the 30 October 2020 Mw 7.0 Aegean Sea Earthquake: Observations and Modelling

Eastern Mediterranean Sea has experienced four tsunamigenic earthquakes since 2017, which delivered moderate damage to coastal communities in Turkey and Greece. The most recent of these tsunamis occurred on 30 October 2020 in the Aegean Sea, which was generated by an Mw 7.0 normal-faulting earthquake, offshore Izmir province (Turkey) and Samos Island (Greece). The earthquake was destructive and caused death tolls of 117 and 2 in Turkey and Greece, respectively. The tsunami produced moderate damage and killed one person in Turkey. Due to the semi-enclosed nature of the Aegean Sea basin, any tsunami perturbation in this sea is expected to trigger several basin oscillations. Here, we study the 2020 tsunami through sea level data analysis and numerical simulations with the aim of further understanding tsunami behavior in the Aegean Sea. Analysis of data from available tide gauges showed that the maximum zero-to-crest tsunami amplitude was 5.1–11.9 cm. The arrival times of the maximum tsunami wave were up to 14.9 h after the first tsunami arrivals at each station. The duration of tsunami oscillation was from 19.6 h to > 90 h at various tide gauges. Spectral analysis revealed several peak periods for the tsunami; we identified the tsunami source periods as 14.2–23.3 min. We attributed other peak periods (4.5 min, 5.7 min, 6.9 min, 7.8 min, 9.9 min, 10.2 min and 32.0 min) to non-source phenomena such as basin and sub-basin oscillations. By comparing surveyed run-up and coastal heights with simulated ones, we noticed the north-dipping fault model better reproduces the tsunami observations as compared to the south-dipping fault model. However, we are unable to choose a fault model because the surveyed run-up data are very limited and are sparsely distributed. Additional researches on this event using other types of geophysical data are required to determine the actual fault plane of the earthquake.


Introduction
The Aegean Sea coasts of Turkey and Greece were struck by a moderate tsunami on 30 October 2020, which was generated by an M w 7.0 normalfaulting earthquake (Fig. 1). With an epicenter close to the city of Izmir (Turkey) and the Samos Island (Greece), the earthquake was destructive and left 117 and 2 deaths in Turkey and Greece, respectively, along with significant damage to residential buildings and other structures in both countries. According to the United States Geological Survey (USGS), the earthquake occurred at 11:51:27 UTC (13:51:27 local time in Greece), its epicenter was located at 26.784°E and 37.897°N, and the focal depth was 21.0 km. The W-phase Moment Tensor solution of the USGS revealed a normal-faulting mechanism for the earthquake with strike angle: 93°; dip angle: 61°a nd rake angle: -91°for the south-dipping fault plane. The corresponding USGS focal mechanism parameters for the north-dipping fault plane are 276°( strike), 29°(dip), and -88°(rake). The USGS chose the south-dipping fault plane as the actual plane of the earthquake based on the agreement between observed and synthetic seismic waveforms, 1 whereas Kiratzi et al. (2020) selected the north-dipping fault plane as the actual fault plane based on regional tectonic setting and aftershocks analysis.
The earthquake produced co-seismic uplift of 15-20 cm along part of the coast of the Samos Island (Triantafyllou et al., 2021).
The consequent tsunami impacted the coasts of Seferihisar district in Izmir province (Turkey) (Fig. 2) and Samos Island (Greece) (Fig. 3), where it caused one fatality in Sigacik (Turkey); no tsunami death was reported from Greece. Dogan et al. (2021) measured a maximum tsunami runup height of 3.8 m in the central Aegean coast of Turkey (Fig. 1). For Greece, the maximum runup was 3.4 m, measured on the north coast of Samos Island (Triantafyllou et al., 2021). The October 2020 tsunami was an important event and gained high media attention as it was a fresh call for strengthening tsunami preparation, mitigation, and awareness in the highly seismic zone of the Eastern Mediterranean and the Aegean Sea regions (Heidarzadeh & Gusman, 2021;Heidarzadeh et al., 2017Heidarzadeh et al., , 2019Ozel et al., 2011;Papadopoulos et al., 2020;Triantafyllou et al., 2021).
From a plate tectonic viewpoint, the Aegean Sea microplate is bounded from the north by the North Figure 1 Epicentral area of the 30 October 2020 M w 7.0 earthquake and the tectonic setting of the Aegean Sea region. The epicentres and mechanisms of two other earthquakes in June-July 2017 are also illustrated. Green triangles show the locations of tide gauges. Data of epicentres and focal mechanisms belong to the United States Geological Survey (USGS) earthquake catalogue. The blue boxes show the geographical areas of Grid-2 (spatial resolution = 10.8 arc-sec) and Grid-3 (spatial resolution = 3.6 arc-sec), which form the nested grid system that we used for numerical modelling of the tsunami c Figure 2 Photos showing tsunami damage due to the 30 October 2020 tsunami in three locations along the Turkish coast. Photos are from authors' field surveys. Zeytineli, Sigacik, and Akarca are the most impacted areas along the Turkish coast. In Zeytineli, all of the boats and coastal structures in the fishery shelter were damaged by the tsunami requiring a total reconstruction. In the Kaleici region of Sigacik Bay, many cafes and shops were destroyed, with a local maximum tsunami height of 2.3 m. In Akarca, the tsunami damaged seafront facilities and summerhouses, where a 1.9 m tsunami splash height was measured on the walls of the damaged house shown in the photo 1532 M. Heidarzadeh et al. Pure Appl. Geophys. mechanism earthquakes have occurred in the vicinity of the recent event on 12 June 2017 (M w 6.3) and 20 July 2017 (M w 6.6) (Fig. 1). The latter event was followed by a damaging tsunami (Dogan et al., 2019;Heidarzadeh et al., 2017;Karasözen et al., 2018). Some unusual observations were reported following the 30 October 2020 tsunami including longlasting tsunami oscillations, which were over a day, and an unusual tsunami run-up height of 3.8 m which are not usually expected from an M w 7.0 normalfaulting earthquake. Long tsunami oscillations are considered as one of the cascading risks of tsunami and earthquake events which further intensify the damaging effects of tsunamis. The other question is that which nodal plane of the earthquake was the actual fault plane as seismic data are unclear about this important question; the USGS and Kiratzi et al. (2020) picked opposite fault models (i.e., north-dipping vs. south-dipping planes) as the actual fault plane of the earthquake. This research is an attempt to address these questions and provides some answers while we leave detailed analyses to future works on this important event. Our methodology is based on tsunami waveform analysis and numerical modeling of tsunami and comparison with actual tsunami tide gauge and run-up data.

Data and Methods
Tsunami data comprises seven tide gauge records (Figs. 1,4) in the Aegean Sea provided by the Sea Level Station Monitoring Facility of the Intergovernmental Oceanographic Commission (http://www. ioc-sealevelmonitoring.org/map.php) and the National Observatory of Athens (NOA) (http://hlntwc.gein.noa.gr/en/). The sampling interval for all tsunami records is 1 min except for Bodrum, which has a sampling interval of 0.5 min. The tidal analysis package TIDALFIT (Grinsted, 2008) was applied to predict the tidal signals and remove them from the tide gauge records (Fig. 4). TIDALFIT applies a Least Square approach to fit tidal components to the observation data. The mean root square amplitude of tsunami (MRSA) is used in this study as a measure of tsunami duration at tide gauge stations. Tsunami duration is defined here as the time interval that tsunami MRSA is above that of the background signal (i.e. before tsunami arrival at each station). Two types of spectral analyses of tsunami waves are conducted here: Fourier and Wavelet (frequencytime) analyses. Our Fourier analysis is based on an updated version of the Welch's (1967) algorithm embedded in MATLAB package (Mathworks, 2021) by considering Hanning windows and 40% overlaps (Heidarzadeh & Satake, 2015). The length of the waveforms inputted in the Fourier analyses was 450 min which corresponds to 450 sea level data points, as the sampling intervals of sea level data was 1 min; except for the Bodrum station, which entailed 900 data points due to its sampling interval of 0.5 min. We also calculated 95% confidence bounds for the spectra. The method of Torrence and Compo (1998) was applied in this study for Wavelet analysis using the Morlet wavelet mother function; this wavelet analysis is similar to the time-frequency analysis of Rabinovich et al. (2021).
The earthquake source model is constructed based on the finite-fault inversion method developed by Shimizu et al. (2020), which has been an efficient tool for robustly estimating the complex rupture evolution of earthquakes (e.g., Okuwaki et al., 2020;Shimizu et al., 2020;Tadapansawut et al., 2021). This method is based on reducing the modelling errors associated with uncertainties in Green's function calculations (Yagi & Fukahata, 2011) and fault geometries . We use the vertical components of 64 teleseismic P waveforms in our analyses (Fig. 5). The Green's functions are calculated based on the method of Kikuchi and Kanamori (1991). We use the ak135 model (Kennette et al. 1995) to calculate travel times, ray parameters, and geometric spreading factors. CRUST 1.0 model is used for the onedimensional layered medium near the source region to calculate Haskel propagator (Kikuchi & Kanamori, 1991) in Green's functions (Table 1). We use two possible model-plane geometries for the finite-fault models, based on the Global Centroid Moment Tensor (GCMT) solution (Dziewonski et al. 1981;Ekström et al. 2012): 96°strike and 53°dip angles (south-dipping model, Fig. 5a), and 270°strike and 37°dip angles (north-dipping model, Fig. 5b). For both model faults, the spatial intervals of sub-faults are 5 km (strike-wise) 9 5 km (dip-wise). The rectangular fault dimension is 105 km (length) 9 35 km (width) (Fig. 5). The total source duration is set as 30 s. The maximum rupture velocity is set at 3.5 km/ s based on the shear-wave velocity near the source region (Table 1). We use the hypocentre at 37.888°N and 26.834°E, and 10-km depth determined by KOERI-RETMC (2020) for the initial rupture point.
Both finite-fault models with the two possible model planes show that the dominant normal faulting is concentrated at the shallow-most, western part of the model faults, with the largest slips of 2.60 m and 1.87 m for the south-and north-dipping models, respectively (Fig. 5). The resultant seismic moments are 5.7 9 10 19 Nm (M w 7.1) and 6.3 9 10 19 Nm (M w 7.1) for the south-and north-dipping models, respectively (Fig. 5). The waveform fitting qualities between observed and synthetic waveforms were calculated using the following equation : where u obs j t ð Þ and u syn j t ð Þ are the observed and synthetic waveforms, respectively, at the jth station and at time t. Calculations of r 2 (Eq. 1) for both southand north-dipping fault models showed less than * 1% difference between the two models: r 2 is 0.372 for the south-dipping model and is 0.370 for the north-dipping model. Therefore, teleseismic data alone is not conclusive about the actual fault plane of the 30 October 2020 Aegean Sea earthquake.
Tsunami simulations are performed using the state-of-the-art numerical package JAGURS (Baba et al., 2015), employing a nested grid system and nonlinear shallow water equations. Bathymetry data are based on the European Marine Observation and Data Network (EMODnet), which comes with a spatial resolution of 3.75 arc-sec (approximately 115 m) (EMODnet Bathymetry Consortium, 2020). A three-level nested grids, named Grid-1 (geographical area 19°E-30°E and 33°N-42°N), Grid-2, and Grid-3 (blue boxes in Fig. 1), are used for tsunami simulations whose spatial resolutions are 32.4 arcsec, 10.8 arc-sec and 3.6 arc-sec, respectively. The time step for Finite Difference calculations was 0.5 s to satisfy the convergence condition of the numerical scheme. Tsunami simulations were conducted for a total duration of 48 h. We excluded inundation modeling but have recorded maximum coastal tsunami amplitudes which are considered as reasonable estimates of tsunami runups (e.g. Tinti et al., 2006). The reasons for excluding inundation modeling are: (i) The EMODnet bathymetry data lacks sufficient spatial resolution for conducting inundation modeling; (ii) The maximum coastal tsunami amplitude provides a good approximation of coastal run-up based on several past studies (e.g., Satake et al., 2013). Simulated tsunami waveforms and maximum coastal amplitudes are compared with observed tide gauge records and surveyed tsunami heights/run-ups, Figure 5 Slip models for the 30 October 2020 Aegean Sea earthquake based on the south-dipping fault plane (strike/dip angles = 96°/53°) (a) and the north-dipping fault plane (strike/dip angles = 270°/37°) (b). The figure also shows the source-time function (also known as moment-rate function) for each case (the gray-shaded plots) as well as the worldwide distribution of the teleseismic data used in this study (the earth's map with red triangles). The star represents the epicentre respectively. Tsunami survey data are reported by Dogan et al. (2021) and Triantafyllou et al. (2021). Quality of fit between observations and simulations is calculated using the root-mean-square misfit index (e) as given by the equation below: where N is the total number of data points of a station used for misfit calculation: i ¼ 1; 2; 3; . . .; N; the two parameters X obs i and X sim i are the observed and simulated values at data point i. Equation (2) gives e for one station; in case several stations are involved, the overall e is obtained by averaging the e of all stations.

Tsunami Amplitudes and Durations
Among the tide gauge records examined here, tsunami amplitudes are hardly above the maximum high tide levels with an exceedance height of 2.1-6.2 cm (Fig. 4a). The maximum zero-to-crest tsunami amplitude ranges from 5.1 cm (Bodrum) to 11.9 cm (Kos) (Fig. 4b). The maximum tsunami waves occur from 1.3 h (Plomari) to 14.9 h (Heraklion) after the first tsunami arrivals at tide gauge stations; this parameter is [ 3.1 h for all stations, except for Plomari (Fig. 6). Such long arrival times for the largest tsunami wave are relatively unusual as compared to other tsunamis. For example, for the case of the 2011 Japan mega-tsunami, Heidarzadeh and Satake (2013) demonstrated that the largest waves arrived within approximately 1-2 h after the first tsunami wave in most tide gauge stations across the Pacific Ocean; although some stations (e.g., Tosashimizu, Japan) received the largest wave approximately 7-8 h after the first arrival. The late arrival of the largest tsunami wave in the Aegean Sea can be potentially attributed to several wave reflections and long tsunami oscillations in the Aegean Sea due to its semi-enclosed nature and the presence of numerous islands (see Fig. 1).
We measured total tsunami duration at each tide gauge station through the tsunami mean root square amplitude (MRSA) index, which reveals how long the tide gauge amplitudes are above the normal oscillation level before the tsunami arrival in each station. Figure 7 illustrates that the tsunami duration is ranging from 19.6 h (Plomari) to [ 90 h (NOA-03). For five stations, the tsunami duration is [ 31.1 h. These tsunami duration times are considered as relatively long and are not usually expected from a tsunami generated by an M w 7.0 normalfaulting earthquake.
In general, it is expected that tsunamis would trigger various oscillation modes of the Aegean Sea basin and associated sub-basins. A semi-enclosed basin such as the Aegean Sea, where numerous small and large islands exist, normally possesses several oscillation modes with periods in the domain of 1-60 min depending on the dimensions and water depths of the basin and sub-basins. Therefore, we expect that the tsunami source period to be mixed with those basin and sub-basin oscillation modes. The long tsunami oscillations observed during the 30 October 2020 tsunami could be attributed to the excitation of such oscillation modes by the tsunami rather than directly by the tsunami source modes. In the next section, Fourier and Wavelet analyses are performed to identify those oscillation mods. Long tsunami oscillations and late arrival of the largest waves can significantly increase tsunami damage and fatalities as the coastal assets would be under the attack of large and energetic waves for longer times; also, local people may return to the coastal areas before the arrival of the largest waves and may sustain injuries or fatalities. Figure 8 presents the results of Fourier analyses for the tsunami waves (blue and pink spectra) and the background signals (black and gray spectra) at each station. The term ''background signal'' here implies part of the tide gauge record before tsunami arrival at each station. The large gap between the tsunami and background spectra is attributed to the tsunami energy (Fig. 8). The periods of the tsunami spectral peaks inform the period of the waves generated either by the tsunami source or by basin/sub-basin oscillations. Based on Fig. 8 (Fig. 8). According to Rabinovich (1997Rabinovich ( , 2009, tsunami source periods are those that appear at the spectra of most of the tide gauge stations. For the October 2020 tsunami, the peak period band of 14.2-23.3 min is present on the spectra of all tide gauge stations and is likely the tsunami source period band. Other peak periods, such as 4.5 min, are absent in some stations. As a way to estimate tsunami source period, we refer to the equation proposed by Heidarzadeh and Satake (2015), which provides an approximation of peak tsunami period (T) using earthquake rupture length (L) and water depth around the epicenter (d):

Tsunami Spectra and Basin Oscillation Modes
where g is the gravitational acceleration (9.81 m/s 2 ). Considering L = 25-30 km (width of the fault plane) and 90 km (length of the fault plane) (Fig. 5), and d = 200-400 m (Fig. 1), Eq. (3) gives tsunami source periods of 13.3-22.6 min (dictated by fault width) and 47.9-67.7 min (dictated by the fault length). The latter period band of 47.9-67.7 min appears to be absent in all tide gauge stations, which is natural because the tsunami period is mostly dictated by the fault width (25-30 km). Therefore, it can be seen that the tsunami source period of 14.2-23.3 min, revealed by spectral analyses of tide gauge data (Fig. 8), is fairly consistent with the estimate made by Eq.
(3). The other peak periods (4.5 min, 5.7 min, 6.9 min, 7.8 min, 9.9 min, 10.2 min, and 32.0 min) are attributed to non-source phenomena such as basin and sub-basin oscillations. Wavelet analysis (Fig. 9) discloses the variations of dominant tsunami peak periods over time. We conducted Wavelet analysis for four stations that are the closest to the epicenter (Fig. 1). The first signals at each station are commonly attributed to the tsunami source periods (Rabinovich, 1997). In other words, non-source periods (i.e. those belonging to basin and sub-basin oscillations) arrive later at tide gauge stations because basin/sub-basin oscillations require longer time to be generated. It can be seen that the first tsunami waves at all stations have dominant periods in the range of 14.2-23.3 min (Fig. 8). This is another confirmation that the period band of 14.2-23.3 min must have been the tsunami source period. According to wavelet plots, tsunami energy is channeled over multiple period bands over time. In Kos, two period channels of 6.9 min and 19.7 min are visible; the tsunami oscillations in the 6.9-min channel lasts approximately 5 h while tsunami energy persists for more than 14 h in the 19.7-min channel. The Wavelet plot of Plomari reveals a long-lasting (* 9 h) and persisting wave at period of 32.0 min which begins its oscillations approximately 2 h after the earthquake origin time and 1 h after the first tsunami arrival at this station. This wave is most likely a basin mode, triggered by the tsunami source waves, because it is generated 1 h after the first tsunami arrival and its period is beyond the tsunami source period of 14.2-23.3 min. A prerequisite for the generation of basin mode oscillations is that the water body must be set in motion by a triggering mechanism (i.e., tsunami source waves), and thus they commence sometime (e.g., one hour or more) after the tsunami arrival time (Heidarzadeh & Gusman, 2021).

Numerical Modelling and the Actual Fault Plane
As the results of teleseismic inversion were not conclusive regarding the identification of the actual fault plane of the earthquake, here we performed tsunami simulations for both fault models (southdipping and north-dipping models) (Figs. 10,11,12). Simulated tsunami waveforms at the location of tide gauges are compared with the observations for both nodal planes (Fig. 10). The average misfit (e) between observed and simulated waveforms is 1.46 for the south-dipping fault model whereas it is 1.10 for the north-dipping plane. Based on Fig. 10, simulations from both south-dipping and north-dipping fault models produce similar waveforms at the location of tide gauges; thus, it is not possible to favour a fault model over another one. So far, two types of geophysical data associated with the October 2020 event (i.e. teleseismic waveforms and tide gauge records) failed in clearly favoring a nodal plane of this normal-faulting earthquake over the other one.
This behavior has been seen for other M B 7 earthquakes in the past which can be partly attributed to the relatively small size of the tsunamigenic earthquake, and partly to the sparse tsunami observations (e.g. Lay et al., 2014). The third geophysical data was the tsunami runup and coastal heights measurements which were collected through post-tsunami field surveys (Dogan et al., 2021;Triantafyllou et al., 2021). We compared such measurements with simulated maximum coastal amplitudes in Fig. 11. It can be seen that the north-dipping fault model (green, Fig. 11) produces smaller misfits (e) relative to the south-dipping fault model (red, Fig. 11), with an average misfit (e) of 0.49 versus 1.1. This is consistent with the results of Figure 11 Comparison of surveyed tsunami heights and run-up data (yellow circles and black rectangles) with computed maximum nearshore tsunami amplitude (peak coastal tsunami amplitude) using the south-dipping fault model (red bars) and north-dipping fault model (green bars) along the coasts of Turkey and Greece for the 30 October 2020 tsunami. The surveyed data are from Dogan et al. (2021) and Triantafyllou et al. (2021). The colormap for the offshore area shows maximum tsunami amplitudes (in meters) at each computational grid point during the entire tsunami simulations. The red and green numbers are tsunami misfits using Eq. (2) for the south-and north-dipping models, respectively 1544 M. Heidarzadeh et al. Pure Appl. Geophys. Kiratzi et al. (2020), who also chose the north-dipping plane as the actual fault plane of the earthquake based on aftershocks analysis and considering regional tectonic setting. Despite this, it is not possible to favour a fault model over another one based on Fig. 11 because the surveyed run-up data are very limited and are sparsely distributed. Another type of geophysical data that could help towards identification of the earthquake's actual nodal plane is the coseismic uplift/subsidence data. Triantafyllou et al. (2021) reported a few such data points; however, the published data are insufficient. The long tsunami oscillations during the 30 October 2020 tsunami are reproduced through numerical simulations with varying degrees of success in different stations (Fig. 7, pink plots). Although the quality of simulated and observed duration times are good in most of the tide gauge stations, some discrepancies are seen in two stations of Bodrum and Heraklion (Fig. 7), which can be partly attributed to the insufficient quality of bathymetry data. Snapshots of tsunami propagation reveal that most of the tsunami is trapped in the Bay area bounded by Izmir and Samos Island (Fig. 12). Regarding coastal tsunami heights, our modeling was able to fairly reproduce the surveyed run-up heights (Fig. 11). Therefore, we may conclude that the relatively large coastal run-ups following this event can be potentially attributed to the confined nature of the source region, which is located in a bay, and to the short distance between the epicenter and the coast (Fig. 1).

Conclusions
We studied the 30 October 2020 M w 7.0 Aegean Sea earthquake and tsunami through analysis of sea level data and numerical modelling. Main findings are: • The maximum zero-to-crest tsunami amplitude was 5.1-11.9 cm on the tide records examined, and the maximum tsunami amplitude occurred 1.3-14.9 h after the first tsunami arrivals. We attribute the late arrivals of the largest tsunami waves to potential several wave reflections and long oscillations in the Aegean Sea as it is a semi-enclosed basin encompassing numerous islands. • Duration of tsunami oscillations on tide gauges varied from 19.6 h to [ 90 h. These long tsunami oscillations are fairly reproduced through numerical modeling. • We identified the period band of 14.2-23.3 min as belonging to the tsunami source. Spectral and Wavelet analyses revealed other peak periods such as 4.5 min, 5.7 min, 6.9 min, 7.8 min, 9.9 min, 10.2 min and 32.0 min, which we attribute to nonsource phenomena such as basin and sub-basin oscillations. • Teleseismic inversion and tsunami simulations and comparison with tide gauge records were inconclusive regarding identification of the actual fault plane of the earthquake; both south-dipping and north-dipping fault models give similar simulation results.
• Comparison of simulated maximum coastal amplitudes with surveyed run-up/heights shows that the north-dipping fault model gives smaller misfit; however, it is not possible to determine the actual fault plane of the earthquake because the surveyed run-up data are very limited and are sparsely distributed.