Origins of the Tsunami Following the 2023 Turkey–Syria Earthquake

On 6 February 2023, a local tsunami was recorded in the southeastern Mediterranean Sea following the Mw 7.7 Turkey–Syria inland strike‐slip earthquake. Due to the lack of underwater observation, the tsunami generation mechanism remains mysterious. To understand the source mechanisms, we analyzed the tsunami waveforms of four nearby tide gauges and located possible sources using a backward tsunami ray tracing approach. We then conducted forward numerical modelings for a range of possible source parameters. We show that there were probably two tsunami sources, inside and outside Iskenderun Bay, which may be related to thick coastal sediments. A source inside the Bay with a characteristic length of 7 km produced dominant periods of 10–30 min with negative initial motion, possibly generated by a landslide. Another source of 6 km length outside the Bay produced dominant periods of 2–10 min with positive initial motions, possibly related with liquefaction.

form, suggesting that they were triggered by seismogenically induced submarine landslides/slumps (Salamon et al., 2007) (Figure 1a).The 2023 Turkey-Syria event broke the long tsunami quiescence in this region, which had persisted since the tsunami associated with the 1953 M L 6.2 Cyprus earthquake (Ambraseys & Adams, 1993).
The 2023 tsunami source is mysterious.The epicenter was ∼90 km inland from the nearest coastline at Iskenderun Bay.The fault rupture of the first mainshock did not extend into the ocean (Melgar et al., 2023).The tsunami was recorded at four tide gauges in southeastern Turkey in the eastern Mediterranean with a maximum height of 40 cm (Figure 1b).Currently, no submarine evidence has been released to identify the tsunamigenic source.The physical origin of the tsunami generation is still not understood due to its mystery and complexity.
In the study, we attempt to provide answers for the tsunamigenic origin.First, we analyze wave behaviors and wavelet energy distribution of the tsunami waveforms at the near-source tide gauges.Second, we use backward tsunami ray tracing to locate the tsunami sources.Finally, we assume Gaussian-shaped sources with dimensions derived from the tsunami waveforms, and then adopt forward numerical modeling to fit the recorded tsunami waveforms.In conjunction with available ground motion observations, remote sensing data, and field-survey investigations, we will demonstrate that the tsunami was generated by two sources, one involving landslide movement nearshore Iskenderun port, and another potentially related to liquefaction offshore the coastal plain of Antakya (Figure 1a).

Data and Method
We downloaded sea level records of the tsunami event from a sea level monitoring website funded by Intergovernmental Oceanographic Commission (IOC) (IOC, 2023).The tsunami waveforms at three near-field tide gauge stations, namely, Hatay, Erdemli, Tasucu in Turkey, are accessible with sampling intervals of 30 s.The epicentral distances of the tide gauges range between 132 and 301 km.The time series between 09:17:35 UTC on 5 February and 01:17:35 UTC on 8 February, or from 16 hr before to 48 hr after the first mainshock, were used for the analysis.The data quality is quite high without abnormal spikes within the time series.In data processing, we first filled short gaps by linear interpolation, and then applied a fourth-order Butterworth high-pass filter with a corner frequency of 3.5 × 10 −5 Hz (∼8 hr) to remove low-frequency tidal components (Heidarzadeh & Satake, 2013;Hu et al., 2022).We also note properties of the initial wave at the Gazimagusa station in Turkey for wave analysis from Yalcinera et al. (2023).
Continuous wavelet analysis was calculated for the tide gauge records based on the Morlet wavelet mother function (Kristeková et al., 2006).Wavelet analysis can detect the periodic variations in spectral peaks over time series.
We applied the backward tsunami ray tracing to locate the tsunami sources.First, we conducted forward tsunami ray tracing to calculate tsunami propagations along different paths using the tide gauge locations as input sources (Satake, 1988;Text S1 in Supporting Information S1).Second, we drew travel-time arcs corresponding to the observed tsunami travel times to the tide gauges.The bathymetric data in the computation was downloaded from the European Marine Observation and Data Network, whose resolution is ∼115 m (EMODnet Bathymetry Consortium, 2020).We adopted JAGURS software for the tsunami modeling (Baba et al., 2015).The simulation solves the tsunami equations with the staggered and leap-frog finite difference method.We applied a nonlinear long-wave model to simulate the tsunami propagation in a spherical coordinate system.We adopted three nested grids (grids 01-05) with the bathymetric and topographic grid resolutions of 18, 6, and 2 arc-second, respectively.The spatial coverages of the nested grids are shown in Figure S1 in Supporting Information S1.The bathymetric data used in the modeling is high-resolution EMODnet data set, which was resampled for the nested grid scheme.The simulations were performed for 3 hours with a time step of 0.2 s.
To quantify the similarity between two time series in the work, we adopted two algorithms, Dynamic Time Warping (DTW) (Berndt & Clifford, 1994) and Geometric Mean Ratio (GMR) (Aida, 1978).Algorithms details are elaborated in Text S2 in Supporting Information S1.DTW takes advantage of a warping technique to calculate the distance between time series, representing the differences and offset between them.Conversely, GMR algorithm evaluates how well they match.Thus, the two input time series with smaller DTW or larger GMR indicate more similarity between them.

Different Tsunami Properties Inside and Outside the Iskenderun Bay
To evaluate the source-related tsunami behaviors, we first obtain the tsunami period band and quantify the similarity between waves in different bands, and then analyze the recorded initial wave properties and the maximum tsunami heights.
The period band of the observed tsunami is distributed in 2-30 min, which is revealed by continuous wavelet energy analysis for three tide gauges in Figure 1c, where strong and clear tsunami energy appears after the first mainshock.At Hatay station, stronger energy following the first mainshock is distributed in 10-30 min.At Erdemli station, tsunami energy is divided into two parts, below and above 10 min, with stronger energy in 2-10 min.At Tasucu, the strongest background energy is above 30 min but tsunami energy is between 10 and 30 min.We compare the tsunami wave components distributed in the period band of 2-10 min and 2-30 min to explore the tsunami source effect on the observed waveforms.We apply DTW and GMR to quantify the similarity between wave components in both period bands for each station (Figures 2a and 2b). Figure 2c indicates that the Hatay record inside the Iskenderun Bay exhibits a larger DTW and a smaller GMR than the other stations outside the Bay, which means that the two wave components are similar outside the Bay, but more different inside the Bay.This suggests that the tsunami outside the Bay is dominated by the wave in 2-10 min, while the tsunami inside the Bay contains two different dominant period bands, 2-10 min and 10-30 min period.The difference implies that the wave in 10-30 min comes from inside the Bay.
Different arrival times and opposite initial wave motions are found between the tide gauges inside and outside the Iskenderun Bay.Figures 2a and 2b, and Figure S2 in Supporting Information S1 show that the tsunami first reached Hatay station at 25 min after the first earthquake with an initial downward motion.The tsunami waves then were recorded by stations at Gazimagusa, Erdemli, Tasucu with first arrival times of 36 min, 40 min, 50 min (Figure S2 in Supporting Information S1), respectively, whose initial wave motions were all upward, opposite to that at Hatay.The initial wave information at Gazimagusa is from Yalcinera et al. (2023).The arrival time and initial motion of the tsunami at Tasucu are read from waves in the 2-10 min period, which are more discernible than those in the 2-30 min period.
The maximum tsunami height (trough-to-crest) at Erdemli and Tasucu appeared at 1-2 hr after their first arrivals (Figure 1b).However, at Hatay station, the time lag of the maximum tsunami height from the initial arrival reaches 40 hr (Figure 1b).As the location of the station is in the narrow semi-closed Iskenderun Bay, we suppose the long delayed maximum waves may be due to tsunami resonance, a coupling of similar period ranges between the Bay's fundamental natural oscillation and source period band (Hu et al., 2022;Satake et al., 2020).We calculate the fundamental oscillation of the Bay based on a semi-parabolic depth profile with   = 2.22 × 2∕ √ (gH) (H is the water depth at the entrance of the Bay, L is the length of the Bay) (Wilson, 1972), in which we set H of 60-80 m and L of 60-70 km based on the bathymetry data set.Figure 2d shows that the fundamental oscillation period of the bay ranges between 159 and 213 min, whose period is of one-to-two orders of magnitude larger 10.1029/2023GL103997 5 of 9 than the tsunami period of 2-30 min, suggesting that tsunami resonance of the entire Bay did not occur.The long delayed maximum tsunami waves at Hatay may be amplified by resonance of smaller scale or another tsunami generation related to other shock in the earthquake sequence.In addition, tsunami energy in 2-4 min period at Hatay appears around the origin time of the second mainshock (∼9 hr after the first mainshock) (Figure 1c).Both unusual characteristics imply contributions from other tsunami generations inside the Iskenderun Bay.
In summary, the tsunami waves at Hatay inside the Iskenderun Bay is distinctly different from the waves outside the Bay, including different initial wave motion (negative) and smaller similarity between waves in the 2-10 min and 2-30 min period bands than those outside the Bay.The differences suggest two different tsunamigenic sources affect the wave behaviors inside and outside the Bay.

Tsunami Sources Inside and Outside the Iskenderun Bay
To gain more knowledge of the tsunamigenic source, we conduct forward tsunami modeling for a range of source scenarios to fit the recorded tsunami waveforms.Conventional inversion approach is inappropriate for the case, because (a) the source of the tsunami event is unknown, thus pre-described fault geometry setting required by conventional inversion approach cannot be satisfied; (b) As the tracked tsunami locations (described below) are nearshore coastal areas with shallow water where tsunami presents strong nonlinear effects, linear relationship between fault slip and observation in conventional inversion is inapplicable.In this study, we first use backward ray tracing to locate tsunami sources (Figure 3).Then, we set Gaussian-shaped sea surface displacements as sources (details in Text S3 in Supporting Information S1).Dimensions of the sources are derived from empirical relationship between tsunami period and typical source length.Finally, we set a range of maximum source amplitudes in Gaussian basis function to fit waveforms at the three tide gauges shown in Figure 1.The best solution is the one with optimal DTW and GMR results (Figures 4a and 4b).
To estimate the source locations using backward tsunami ray tracing, the origin time of the tsunami should be set first.The clear tsunami energy begins to appear around the time of the first mainshock in Figure 1c.Besides the first mainshock, another candidate is the largest aftershock of Mw 6.6 at 11 min after the first mainshock.We assume the origin times of the two tsunami scenarios to be the same with both candidates.The travel-time arcs of the rays to corresponding stations are marked with different colors (Figure 3).The potential source location for initial wave is at the converging zone of the rays.For the scenario following the first mainshock (Figure 3a), the travel-time rays converge on the shelf offshore the coastal plain of Antakya, where source2 (source outside the Iskenderun Bay) is located.Adding the ray from Tasucu station does not change the results, whose tsunami arrival time is blur in the source period of 2-30 min but clear in 2-10 min (Figures 2a and 2b).For another scenario following the largest aftershock (Figure 3b), the rays fail to converge on a zone, even for those outside the Bay, indicating that tsunami was not generated by the largest aftershock.The scenario following the first mainshock suggests that only one source exists.Tsunami simulations in Figures 4c and 4d suggest source2 can only explain the observed waveforms at Erdemli and Tasucu outside the Bay but can not reproduce those at Hatay inside the Bay.The unreproduced waveforms inside the Bay imply another source inside the Bay exists (source1), which is consistent with waveform property analysis in Section 3. The location of the source1 inside the Bay is very likely nearshore the Iskenderun Fishery port (Figure 1a), where considerable water inundation has been observed due to coastal area subsidence most probably caused by the liquefaction-induced landslide according to field survey (Maria et al., 2023).Therefore, the exact location of source1 should be the nearest point of the red line (Figure 3a) to the Iskenderun Fishery port.Both sources are distributed in the proximity to the largest PGAs of the first mainshock (Figures 4e-4g) (Erdik et al., 2023), suggesting the sources are closely associated with the strong ground shaking generated by the earthquake.Additionally, sensitivity tests suggest that the resultant ray fronts from coarser bathymetric data are still similar (Figure S3 in Supporting Information S1).
Typical dimensions of the Gaussian-shaped sources are calculated by  =  √ (gH)∕2 (H is the mean water depth at source, T is dominant tsunami period).According the EMODnet bathymetry, we set H of 15 m and T of 10-30 min for the source1.We set H of 130 m and T of 2-10 min for the source2, thus the typical dimensions range 4-10 km and 2-10 km, respectively.The median of the dimension ranges are chosen for both sources, that is 7 and 6 km, which we set as characteristic lengths (see Text S3 in Supporting Information S1).
Finally, we set a series of combinations of maximum initial displacements for both sources to fit the tide-gauge waveforms.The optimal maximum displacements are obtained iteratively by finding the solution that yields the best fit.Figures 4a and 4b show that the source1 and source2 with respective maximum displacements of −0.8 and 1.1 m can well account for the observed waveforms inside and outside the Iskenderun Bay.If we increase or decrease the maximum displacements by 0.3 m for both sources, the agreements get worse quantitively.If we only use source2 for the tsunami simulation, the waveforms outside the Bay can be well explained but those inside the Bay can not (Figure 4c).When adding source1 inside the Bay, the waveforms inside the Bay can be well explained as well (Figure 4d).The results manifest that both sources together can reproduce the observed waveforms.We also did sensitivity tests to the source parameters, which indicates varying locations and characteristic length (2, 4, and 8 km) for the source2 perform not as good as original parameters in reproducing either tsunami periods or amplitudes of later phase at Erdemlt.Therefore, the tsunami simulations demonstrate that two sources are necessary to jointly account for the tsunami event.
To explore the physical mechanisms of both sources, we conduct forward tsunami modelings for the first mainshock based on two finite fault models from USGS and Melgar et al. (2023) (Figures S4a and S4b in Supporting Information S1).The underwater parts of their initial surface elevations are very limited with the maximum values less than 5 cm (Figures S4c and S4d in Supporting Information S1), whose modeled waveforms are too small to match the observations.According the post-event investigations, the largest PGAs of the first mainshock were recorded near both source zones (Figures 4f-4h).The coastal areas near both source zones are covered by dense liquefaction sites based on remote sensing (Taftsoglou et al., 2023).According to filed survey, at Iskenderun port, near source1, considerable seawater invaded and inundated the subsided area during earthquake shaking most probably caused by the liquefaction-induced landslide (Maria et al., 2023) or subfault rupture triggered by the first mainshock.For source2, current evidences support it is likely caused by a gravity flow of liquefaction-induced collapse near coastal land or landslides due to the first mainshock, where dense liquefactions are found in the coastal plains of Antakya (Taftsoglou et al., 2023).Maria et al. (2023) show that the two sources are located in areas where rich liquefactions are found distributed near coastal/fluvial areas covered by thick Holocene loose sediments.Similar cascading hazards occurred during the 2018 Sulawesi earthquake event.Most of the landslides and liquefied gravity flows contributing to the 2018 Sulawesi tsunami are located near estuaries with dense sediments (Liu et al., 2020;Sassa & Takagawa, 2019).

Awareness of Atypical Tsunami Sources Triggered by Coastal Strike-Slip Earthquake
The 2023 Turkey-Syria tsunami demonstrates the significance of the tsunami hazards from atypical sources, especially due to coastal strike-slip earthquake.Typical tsunami sources generally refer to megathrust earthquakes in subduction zones, like the 2004 Mw 9.1 Sumatra-Andaman earthquake and the 2011 Mw 9.0 Tohoku earthquake.Recently several tsunamis from atypical sources present their underestimated and unexpected impacts.For example, the maximum runup of 3.8 m was produced by an Mw 6.9 normal-faulting earthquake in the eastern Aegean Sea (Doğan et al., 2021;Hu et al., 2022).The 2022 Hunga Tonga-Hunga Ha'apai volcanic eruption caused global ocean disturbance (Carvajal et al., 2022;Hu et al., 2023;Kubota et al., 2022;Lynett et al., 2022).Although the 2023 Turkey-Syria tsunami only caused limited impact, the event presents a variety of the atypical sources and the uncertainty of their spatial occurrences caused by inland/coastal strike-slip earthquake.A stacking of the tsunami waves distributed in different period bands bring additional complexity and unpredictability for tsunami impact.
Previous examples also demonstrate the dangers from seismically induced landslides/slumps and liquefactions triggered by strike-slip earthquakes (Li et al., 2022).On 26 December 1939, an earthquake of Mw 7.8 occurred along the right-lateral strike-slip North Anatolian Fault.A tsunami similar with the 2023 Turkey-Syria event occurred in the Black Sea.The rupture zone of the earthquake was also inland with epicentral distance of 160 km from the closest coast of Black Sea.Although all tide gauges recorded limited negative leading waves ranging ∼0.5-1 m, extremely large initial sea recessions were observed at some specific sites through eyewitnesses, such as 50 m at Fatsa and 100 m at Unye (Altinok & Ersoy, 2000;Papadopoulos et al., 2011).Induced submarine landslide or a secondary fault is favored mechanism for the tsunami event (Papadopoulos et al., 2011;Yalçiner et al., 2004).Another similar tsunami in the region was excited by 17 August 1999 dextral strike-slip earthquake of Mw 7.4.The large runup of 15 m at Degirmendere was due to coseismic subsidence and dominant contribution of at least a local coastal landslide (Tinti et al., 2006).For recent events, the tsunami of the 2018 Mw 7.5 Sulawesi strike-slip earthquake with the contribution of triggered landslides led to 11 m runup in the southern part of the Palu bay (e.g., Carvajal, 2019;Liu et al., 2020).An extensive liquefactions generated by the 2010 Mw 7.0 Haiti earthquake in the river delta area resulted in the maximum inundation height greater than 3 m at the coast (Hornbach et al., 2010).A difference between the 2018 Sulawesi and the 2010 Haiti cases is that the induced collapses caused tsunamis at least at nine locations throughout the Palu bay near coastal areas, not at one site, as in the Haiti event.

Conclusion
Unexpected tsunami was recorded in the gulf of southeast Turkey in the eastern Mediterranean Sea following the 6 February 2023 Mw 7.7 and Mw 7.6 Turkey-Syria inland strike-slip earthquakes.As both mainshocks' fault ruptures did not extend to the ocean, the physical generation mechanism of the tsunami remains unknown.We analyzed the tsunami records from four tide gauges and performed tsunami modeling for a range of possible sources to understand the tsunami mechanisms.Wavelet analysis reveals that the dominant tsunami period is distributed in 2-30 min, but the tsunami inside the Iskenderun Bay was dominated by wave components in 10-30 min with negative initial motion.Outside the Bay, tsunami waves mainly came from components in 2-10 min with positive initial motion.Backward tsunami ray tracing shows that the tsunami source (source1) inside the Bay was likely generated by the landslide movements at the Iskenderun port with initial negative elevation of 7 km in diameter.The tsunami source (source2) outside the Bay was related with liquefaction in the coastal plains of Antakya with initial positive elevation of 6 km in diameter.Although observed maximum tsunami height of the event was only 40 cm, remarkably different wave properties and generation mechanisms raise our awareness of the underestimated tsunami hazards due to the disaster chain caused by coastal strike-slip earthquake.

Figure 1 .
Figure 1. Background of the 2023 Turkey-Syria tsunami event.(a) Tectonic background.The black bold lines indicate the plate boundaries and general plate motion direction in the eastern Mediterranean.The black thin lines are the main faults (Styron & Pagani, 2020).The moment tensor solutions for the first mainshock, the largest aftershock, and the second mainshock in the earthquake sequence from the GCMT catalog.The yellow hexagons represent the locations of historical tsunamis between 2100 BCE and 2023 CE (NOAA, 2023).Triangles with different colors represent the locations of different tide gauge stations.(b) Detided tsunami waves at the near-source stations shown in (a).(c) Wavelet analysis of the tsunami waves in (b).

Figure 2 .
Figure 2. The filtered tsunami waveforms.Similarity of the waveforms in 2-30 min (a) with those in the period band of 2-10 min (b) are evaluated with DTW and GR methods (c).Panel (d) is the fundamental natural oscillation period in the semi-closed Iskenderun Bay.

Figure 3 .
Figure 3. Tsunami source locations.Backward tsunami ray tracing to locate the tsunami source area for two scenarios following the first mainshock (a) and the largest aftershock (b).The arrows in (a) mark the potential locations of the sources.Solid triangles represent the locations of tide gauges whose colors are the same as the ones of their traveling time contours.

Figure 4 .
Figure 4. Tsunami source estimations.Goodness of fit (GOF) between modeled and observed tsunami waves according to different source combinations are evaluated with average DTW (a) and mean GMR (b) methods.Panels (c) and (d) represent the data fits of the models between modeled (red lines) and observed (black lines) tsunami according to two scenarios, that is, source2 only, and source1 and source2 together.Panels (e)-(g) are the distributions of the maximum peak ground accelerations (PGAs) of north-south and east-west horizontal components, and up-down vertical components, respectively.The PGAs data sets are from Erdik et al. (2023).