Tectonic Origin Tsunami Scenario Database for the Marmara Region

This study presents a tsunami scenario database in Marmara Sea, Turkey referring to 30 different earthquake 10 scenarios obtained with the combinations of 32 possible fault segments, complemented with three additional worst-case scenarios based on catastrophic historical earthquakes. The main motivations of this study are to create a tsunami database which may assist the Regional Earthquake and Tsunami Monitoring Center (RETMC) of Kandilli Observatory and Earthquake Research Institute (KOERI) in situation assessment in case of a real event in the Marmara Sea and to investigate the nature and generation of historical tsunamis that were significantly affected Marmara region. Tsunami simulations have been 15 performed using tsunami numerical code NAMI DANCE which solves Nonlinear Shallow Water Equations (NLSWE) using leap-frog scheme. For each earthquake scenario, tsunami hydrodynamic parameters, mainly maximum water surface elevations, arrival time of first wave and maximum wave, and time histories of water level fluctuations were calculated at 1333 synthetic gauge points meticulously selected in shallow zone offshore the coasts of Marmara Sea. Among all analyses in this study the rupture of a segment along Main Marmara Fault generates the largest tsunami source and consequently the highest 20 wave amplitudes in the Sea of Marmara, reaching up to 2.2 m. According to the integrated maximum wave heights, the most significant effects are expected along Silivri at the north, Kursunlu, Bayramdere and Cinarcik at the south, and western Kadikoy, Princes’ Islands, Maltepe and Pendik at the east. Arrival times of first waves are quite short, almost within 5 minutes for the locations close to the earthquake origin. It should be noted that, however, the maximum wave height obtained due to this deterministic study is relatively lower in comparison to the available probabilistic studies published so far. The overview 25 of the results confirm that higher historical tsunami wave heights observed in Marmara Sea cannot be explained by only earthquake-generated tsunamis, as also indicated by recent related publications as there is strong agreement on considering submarine landslides as the primary component of tsunami hazard in the Marmara.

hosting mainly the Istanbul megacity with a population of more than 15 million and huge capacity of trading and tourism, various industrial facilities, ports, airports as well as other densely populated cities.
Marmara Sea is seismically very active as part of the North Anatolian Fault Zone (NAF) with a history of large earthquakes and possibility of experiencing earthquakes with magnitudes larger than 7.0 in the future (Ambraseys, 2002). There is a large number of references in literature proving the generation of historical tsunamis in the Marmara Sea caused by these earthquakes 5 (Mihailovic 1927;Gundogdu 1986;Oztin and Bayülke 1991;Oztin 1994;Ambraseys andFinkel, 1987, 1995;Altınok and Ersoy 2000;Arel and Kiper, 2000;Altınok et al. 2001Altınok et al. , 2003Altınok et al. , 2011Yalciner et al. 2002;Ambraseys 2002;Cetin et al. 2004;Rothaus et al. 2004;Tinti et al. 2006;Özel et al. 2011). The most significant ones are the Istanbul earthquakes on 10 September 1509, on 22 May 1766, and on 10 July 1894; the Sarkoy-Murefte earthquake on 9 August 1912; and the 17 August 1999 Kocaeli earthquake (Necmioğlu, 2016). During the 1509 earthquake, with a magnitude close to 7.5 (Bulut et.al, 2019), the sea 10 flooded the shores along Istanbul coasts, waves crashed against city walls and around 4000-5000 people died in the city (Ambraseys and Finkel 1995). Tsunami waves with probably more than 6.0 m height overtopped the city walls and caused flooding (Oztin and Bayülke 1991;Yalciner et al, 2002). 1766 Istanbul earthquake, on the other hand, triggered considerable tsunami in Besiktas, Istanbul and inner parts of Bosphorus strait, Gemlik Bay and Mudanya in Eastern Marmara (Ambraseys and Finkel 1995;Altınok et al. 2003Altınok et al. , 2011 (The estimated origin of these earthquakes as well as observed runup values 15 according to historical records are shown in Figure 3 in the following section).
The main motivation for this study is investigating the nature of historical tsunamis in Marmara Sea, namely whether they are generated solely due to those significant earthquakes or additional sources, such as submarine landslides triggered by the earthquakes. Due to short arrival times of first waves in Marmara coasts, having prepared tsunami scenarios covering various possible earthquakes is considered as vital, and a comprehensive tsunami database may also assist the Regional Earthquake For this purpose, we have identified a set of credible worst-case earthquake scenarios for the whole Marmara Region obtained 25 by the compilation of historical records, past studies in literature and empirical results. These scenarios constitute the basis for tsunami numerical modelling conducted to obtain tsunami scenario database in Marmara Sea.

Identification of Earthquake Scenarios
The main structural element controlling the morphological and structural features in Marmara Sea region is the northern strand 30 of the NAF, which is considered as a predominantly strike-slip displacement zone (Alpar and Yaltırak, 2002). The NAF, which is a major right-lateral strike-slip fault, extends more than 1200 km from eastern Turkey to the north Aegean Sea (Sengör et al., 2005). It accommodates the relative right-lateral motion between the Anatolian region and Eurasia at a geodetic rate of ~25 mm/yr (Meade et al., 2002;Reilinger et al., 2006). Along its westernmost segment, the fault bifurcates into northern and southern branches, the northern branch following Izmit Bay and entering the Sea of Marmara southeast of Istanbul. By far the majority of long-term fault slip occurs on the northern fault branch following the northwest striking Princes' Islands Fault (PIF) and joining the east-west striking Central Marmara Fault (CMF) immediately south of Istanbul 5 Armijo et al., 2005). After traversing much of the Sea of Marmara, the CMF merges with the Ganos Fault, exiting the Sea along the Ganos Peninsula. Ergintav et al. (2014) concluded that the Princes' Islands segment is most likely to generate the next M > 7 earthquake along the Sea of Marmara segment of the NAF. Armijo et al. (2005) stated that a zone of maximum loading with at least 4-5 m of slip deficit encompassing the 70 km long strike-slip segment between the Cinarcik and Central Basins would alone be capable of generating an earthquake with Mw 7.2. Hergert et al. (2011) argues that the Main Marmara 10 Fault can be interpreted as a through-going fault that slips almost purely in a strike-slip sense, but they also point out that, not contradictory to the previous statement, there is significant dip-slip motion at some sections of the Main Marmara Fault. The South Marmara Fault lies between the highly active northern branch and the weakly active (but still capable of generating magnitude 7 earthquakes) southern branch .
Based on the databases and literature review, 32 faults segments were simplified in order to be able to use them as input for tsunami modelling, where each segment correspond to a rectangular area with an associated hypothetical uniform slip. The locations of each segment are given in Figure 1 where the segments lying on the same fault system (i.e. MMF, SMF, CMF, PIF) are shown in corresponding boxes. All parameters required for the identification of the segments, such as geographical 25 coordinates for the start-and end-points of the segments, hypocentre, type of fault, strike, dip, rake, length and width of the segment, focal depth (where the top of the fault has been set to 0.5 km depth) and corresponding displacements according to empirical relations provided by Wells and Coppersmith (1994) are also presented in Table S.1. The formula provided by Wells and Coppersmith (1994) and associated parameters are given in Table 1. Standard errors defined by Wells and Coppersmith (1994) have been also considered in the Mw calculations to determine Mw (min) and Mw (max) values. Corresponding 30 Moment values have been calculated from the Mw= 2/3 logM0 −10.7 (Hanks and Kanamori, 1979). Corresponding displacement has been obtained from M=μAD, where A is the rupture area, D is the displacement in m and μ is the rigidity modulus taken as 3.25x1011 dyn/cm2.  This was followed by the definition of different hypothetical rupture scenarios reaching a total number of 30 scenarios. The total earthquake moment for each scenario is derived from the summation of the moments associated to the individual segments 15 considered to be ruptured in a given scenario. Slip values have been assigned randomly without any prior assumption, so that heterogeneous earthquake rupture scenarios can be represented. It should be noted, that while the individual segments are assigned with a homogeneous slip, scenarios (except SN06 and SN25) are composed of a combination of segments which represents a heterogeneous slip distribution to the extent possible. These 30 earthquake scenarios derived by this approach are shown in Figure 2 and the details of each scenario in terms of slip values and moment magnitudes are described in Table S.2  20 in Supplementary Material. It is arguable that the maximum earthquake scenarios with Mw 7.4 obtained by Wells and Coppersmith (1994)  Associated slip values for these earthquakes have been derived from Hanks and Kanamori (1979), considering the fault length 15 (L), fault width (W), Mw (Bulut et al., 2019) and rigidity modulus (3.25x1011 dyn/cm2), as shown in Table 2. It's noteworthy, however, that the empirical relationships proposed by Wells and Coppersmith (1994) results with lower Mw and associated slip values for these earthquakes (shown as WC94 in Table 2). The reason for that is mainly due to the fact that Mw values proposed by Bulut et al. (2019) are based on a mean slip deficit rate of 9.9 mm/year derived from a summation of the seismic moment released by historical earthquakes for a period of 1500 years. 20

Tsunami Numerical Modeling of 30 earthquake scenarios
Tsunami numerical modeling has been performed based on 30 earthquake scenarios defined in Section 2.1 using the numerical 5 code NAMI DANCE (NAMIDANCE, 2011). NAMIDANCE solves Nonlinear Shallow Water Equations (NLSWE) with leap-frog scheme both in Cartesian and Spherical coordinate system, which was tested, validated and verified against analytical solutions, laboratory measurements and field observations in several scientific articles (NTHMP, 2017;Lynett et al., 2017;Velioğlu, 2017;Ozer Sozdinler et al, 2015).
Tsunami numerical modelling is performed using 90m grid sized bathymetry -topography data as a single study domain. It 10 was prepared by compiling various data as multi-beam bathymetric measurements, 900m grid sized GEBCO data in the sea as well as 30m grid sized ASTER data on land. Besides, coastline and coastal defence structures i.e. breakwaters, groins and large docks in the ports were also digitized in GIS environment and added to bathymetry -topography data for increasing the resolution and precision in coastal zones. Higher-resolution ASTER data has an important role in data compilation process as it is denser compared with the bathymetry data. In that way, interpolation between less sensitive bathymetry data and much 15 denser topography data provides more reliable coastline in 90m grid sized study domain. The precision in coastline supports the process of selecting synthetic gauge points in shallow zone very close to shoreline, which is described in coming sections below.
The initial sea surface at the time of fault rupture for each segment has been calculated using Okada (1985) formula. In each scenario, it was assumed that all designated fault segments are ruptured simultaneously. NAMIDANCE calculates the sea 20 surface after the rupture of each segment and combines each calculation to output the resultant sea surface as the tsunami source of each scenario. For instance, for Earthquake Scenario #1 (SN01), segment-1, segment-2, segment-3 and segment-4 are the designated fault components for this scenario. The sea surface for each fault segment is calculated using Okada (1985) formula. NAMIDANCE then outputs a resultant sea surface as the combination of these simultaneously ruptured four segments as an overall tsunami source for SN01. We applied the same procedure for all earthquake scenarios accordingly. The maximum 25 and minimum water surface elevations corresponding to the tsunami initial condition computed for each scenario are given in Table S.3 in Supplementary Material. As seen from the table, the initial sea surface disturbances for all scenarios are less than 1 m. The highest sea surface was calculated for SN23, which includes the rupture of segments 17, 18 and 19 located at the center of Marmara Sea striking in NW-SE direction.
The synthetic gauge points along the coasts of Marmara Sea were first selected with very sensitive analysis so as to locate them in shallow zone at water depths less than 50 m. We basically considered the locations of industrial facilities, residential areas, harbors, marinas, factories and six Tsunami Forecast Points (TFPs) while selecting those gauge points (TFPs are synthetic gauge points located at Marmara Eregli, Haydarpasa, Yalova, Mudanya, Erdek and Degirmencik, where the arrival time of first wave and tsunami alert level are calculated and included in national tsunami alert messages disseminated from 5 Regional Tsunami and Earthquake Monitoring Center in KOERI; see Figure 3. There are additional 42 TFPs along the coasts of Turkey in Black Sea, Aegean Sea and Mediterranean Sea). After the selection of synthetic gauge points, test runs were performed in order to identify the water depth where NAMIDANCE located each gauge point as the software assigns each synthetic point at the nearest grid node in bathymetric and topographic data. In other words, although gauge points were selected in the sea within the shallow zone less than 50 m water depth they may be relocated on land or at locations deeper 10 than expected due to the input principles of NAMIDANCE. For that reason, test analyses are critical to validate that synthetic gauge points are located in shallow zone at the possible shallowest location. After these validation analyses, the total number of 1333 gauge points were defined, most of which were located at the water depths of less than 10 m (water depths at some of the gauge points are higher than 10m due to steep topographic conditions at some regions). It is noted that the northern part of the area has much more critical locations than the southern part; therefore, gauge points in that region are denser than the 15 southern part of the Marmara.
Tsunami simulations were conducted for each scenario for 2 hours using the corresponding tsunami sources. The simulation time used for 30 earthquake scenarios was defined after trial simulations performed for the selection of synthetic gauge points as described previously above. As a result of these analyses, it was decided as 2 hours wave propagation in the Sea of Marmara would be adequate to obtain maximum wave amplitudes as well as arrival time of first waves at gauge points. The water level 20 fluctuation was still observed after 2 hours; however, the wave heights were not so significant and less than the maximum values. It should be also noted that such wave motion was observed mostly at the gauge points located in enclosed basins where wave reflections are significant.
Tsunami hydrodynamic parameters as maximum and minimum wave amplitudes, arrival times of first and maximum wave, flow depths and current velocities were calculated throughout Marmara basin and at 1333 synthetic gauge points. We discuss The maximum wave amplitudes less than 75 cm were coloured with green as a representative of relatively safer coastal zones of Marmara Sea. Several experimental and numerical studies in literature prove that 30 cm inundation depth can be accepted as threshold for critical water level that has potential to cause a person to fall down (Jonkman and Penning-Rowsell, 2008;Takagi et al., 2016). So, according to the distribution and relation between the maximum wave amplitudes calculated at gauge points in the sea and inundation depth values calculated on land, the value of 75 cm wave amplitude was decided as to define relatively safer coastal zones.

Summary of Results
As described in previous section, the simulation results are presented here as an integrated distribution of maximum wave amplitudes for overall tsunami scenario database in Marmara Sea. The maximum wave amplitudes were calculated at each 15 synthetic gauge point for 30 earthquake scenarios. The calculated results of 30 scenarios at each gauge point were sorted from higher to smaller and the highest value was stored as the representative maximum wave amplitude at this gauge point. After defining all maximum wave amplitude values at 1333 synthetic gauge points, their integrated distribution was plotted for the entire Marmara Sea (Figure 4). As described in previous section, the coastal zones with green color represents relatively safer locations according to the earthquake scenarios in this database. 20 Following the same procedure, the integrated distribution of arrival time of maximum waves (that is the time of occurrence of maximum wave amplitude at each gauge point in entire earthquake scenarios) was plotted in Figure 5. The results show that the arrival of maximum waves is expected at Princes' Islands, Yalova coasts, some parts of Kadikoy and Silivri coasts within Summary of the simulation results are given in Table 3

Simulation of 1509, May 1766 and August 1766 earthquakes
Three catastrophic historical earthquakes, one in 1509 and two in 1766 (Bulut et al, 2019), were simulated as additional credible 5 worst-case scenarios referring to input data given in Table 2. The earthquake of 10 September 1509 Mw 7.5 is similar to scenario SN13, whereas earthquakes of 22 May 1766 Mw 7.3 (1766a) and 5 August 1766 Mw 7.4 (1766b) are close to scenarios SN4 and SN1, respectively, but with considerably higher slip values as described in earlier sections. The simulation time for the analyses of these historical earthquakes were used as 4 hours this time to observe the change of tsunami hydrodynamic parameters in longer wave propagation as the slip values here are much higher than previously simulated 30 scenarios. 10 Maximum and minimum initial water surface elevations of tsunami sources for these three scenarios are given in Supplementary Material (Table S.

Discussion
This study presents a deterministic tsunami scenario database, based on 30 different earthquake scenarios obtained with the combinations of 32 possible fault segments, complemented with three additional credible worst-case scenarios based on 10 historical earthquakes occurred in 1509, May 1766 and August 1766 referring to Bulut et al. (2019). The results of this deterministic scenario database were presented as arrival times of first and maximum waves, maximum wave heights along Marmara coasts as well as integrated maps, water level fluctuations at selected tsunami forecast points and critical important locations in the Marmara Sea. The overall simulation results reveal that the maximum wave height calculated in sensitively selected synthetic gauge points in shallow zone within entire tsunami scenario database does not exceed 2.2 m. Among all 15 scenarios, SN06 (rupturing of Segment-5 located along Main Marmara Fault) generates the largest tsunami source and consequently the highest wave amplitudes in the Sea of Marmara. The highest wave amplitude calculated among three historical earthquake scenarios is 2 m for 1509 earthquake. However, this value is quite lower than the historical wave records (more than 6.0 m wave height overtopped the city walls and caused flooding (Oztin and Bayülke 1991;Yalciner et al, 2002)). This outcome shows that higher historical tsunami wave heights observed in Marmara Sea cannot be explained by only earthquake generated tsunamis. Submarine landslides should be inevitably considered as the primary tsunami hazard 5 component in the Marmara Sea. As proved by several previous studies, possible tsunamis from submarine landslides in the Marmara Sea could be significantly higher than those from earthquakes depending on the landslide volume. In addition, waveforms from all the coasts around the Marmara Sea indicate that other coastal residential areas in the Marmara region might be subject to high risk of tsunami hazards from submarine landslides, which can generate higher tsunami amplitudes and shorter arrival times, compared to Istanbul (Latcharote et al., 2016). Latcharote et al (2016) argued that the maximum 10 tsunami height could reach 4.0 m along Istanbul shores for a full submarine rupture of the NAF, with a fault slip of 5.0 m in the eastern and western basins of the Marmara Sea, which would correspond to an earthquake with Mw 7.6. However, the maximum tsunami height for landslide-generated tsunamis from small, medium, and large of initial landslide volumes (0.15, 0.6, and 1.5 km3, respectively) could reach 3.5, 6.0, and 8.0 m, respectively, along Istanbul shores and therefore possible tsunamis from submarine landslides could be significantly higher than those from earthquakes, depending on the landslide 15 volume significantly.
Moreover it should be noted, however, that the maximum wave height calculated in this study is relatively lower in comparison to the available probabilistic studies published so far, such as Hancilar (2012), which provide inundation maps resulting from probabilistic tsunami hazard analysis for a 10% probability of exceedance in 50 yr including the building numbers and types, lifeline systems and demographic data in Istanbul, reaching run-up height of 5-6 m. Hancilar (2012) also highlights that the 20 residential buildings at risk are mainly located in Kadikoy, Tuzla, Bakirkoy and Princes' Islands where our study points out significant wave heights as well.
As indicated earlier, slip values for the scenarios considered in this study has been assigned randomly, reaching to a maximum of 4.5 m on the Main Marmara Fault, but comparable to the slip values used in many previous studies. Latcharote et al. (2016) proposes that the Mw 7.2 eastern rupture, the Mw 7.5 western rupture, and the Mw 7.6 full rupture scenarios considered by 25 Hebert et al. (2005) generated a mean displacement on the fault of 4.3, 4.0, and 4.1 m, respectively. Tinti et al. (2006) suggested the generating mechanisms of the 1999 Izmit bay tsunami, in which the Marmara fault in the eastern rupture had the slip of 1.0 m. Oglesby and Mai (2012) suggested the fault geometry and stress parameters along the NAF, including the final slip patterns of less than 3.5 m. Ergintav et al. (2014) suggested the slip deficits accumulated, since the last major earthquake which is less than 2.0 m for the western rupture and 3.7 m for the eastern rupture. Last but not least, Armijo et al. (2005)  It should be also noted that the recent literature addresses the possibility of a tsunami within a strike-slip fault system due to direct earthquake-induced uplift and subsidence, as it was the case in Palu Bay after the September 2018, Mw 7.5 Sulawesi earthquake (Ulrich, et al. 2019). Nevertheless, as also indicated by Goda et al. (2019), there are mixed opinions in the current literature, with regard to the devastating tsunami damage caused by this earthquake, and the possibility of landslide component cannot be excluded (Heidarzadeh et al, 2019;Pakoksung et al, 2019;Mikami et al., 2019;Yolsal-Çevikbilen and Taymaz, 5 2019). Future sensitivity studies complemented with probabilistic methods considering these phenomena could definitely provide a better insight to tsunami hazard and risk in the Marmara region.
The tsunami database considered in this study addresses only possible earthquake generated tsunamis and calculated maximum wave height does not reach the level of the values historically observed. This outcome points out that a tsunami warning system for the Marmara region should be decoupled from the earthquake parameters. To address this issue, Necmioglu (2016)  10 proposed a tsunami warning system in the Marmara region coupled with the existing earthquake early warning system, which could work without waiting for any focal mechanism parameter determination that may lead to underestimate tsunami hazards in the case of a strike-slip fault earthquake, due to the fact that submarine landslides could generate large tsunamis in the Marmara Sea. 15