Probabilistic tsunami hazard analysis for Tuzla test site using Monte Carlo simulations

In this study, time-dependent probabilistic tsunami hazard analysis (PTHA) is performed for Tuzla, Istanbul, in the Sea of Marmara, Turkey, using various earthquake scenarios of Prince Island Fault (PIF) within the next 50 and 100 years. The Monte Carlo (MC) simulation technique is used to generate a synthetic earthquake catalogue, which includes earthquakes having moment magnitudes between M_w6.5 and 7.1. This interval defines the minimum and maximum magnitudes for the fault in the case of an entire fault rupture, which depends on the characteristic fault model. Based on this catalogue, probability of occurrence and associated tsunami wave heights are calculated for each event. The study associates the probabilistic approach with tsunami numerical modeling. The tsunami numerical code NAMI DANCE was used for tsunami simulations. According to the results of the analysis, distribution of probability of occurrence corresponding to tsunami hydrodynamic parameters is represented. Maximum positive and negative wave amplitudes show that tsunami wave heights up to 1 m have 65 % probability of exceedance for the next 50 years and this value increases by 85 % in the Tuzla region for the next 100 years. Inundation depth also exceeds 1 m in the region with probabilities of occurrence of 60 % and 80 % for the next 50 and 100 years, respectively. Moreover, probabilistic inundation maps are generated to investigate inundated zones and the amount of water penetrated inland. Probability of exceedance of 0.3 m wave height ranges between 10 % and 75 % according to these probabilistic inundation maps, and the maximum inundation distance calculated in the entire earthquake catalogue is 60 m in this test site. Furthermore, synthetic gauge points are selected along the western coast of Istanbul by including Tuzla coasts. Tuzla is one of the areas that shows high probability exceedance of 0.3 m wave height, which is around 90 %, for the next 50 years while this probability reaches up to more than 95 % for the next 100 years.

1 Introduction

The Marmara region, especially highly populated cities along the coasts of the Marmara Sea, is the heart of the Turkish economy in terms of having a great number of industrial facilities with the largest capacity and potential, refineries, ports and harbors. The Marmara Sea and surrounding area is one of the most seismically active areas in Turkey. Main active faults of the region pass through the Marmara Sea. Thus, coastal cities in the Marmara region, especially Istanbul, which has significant importance in terms of the economy and historical and sociocultural heritage with a population of more than 15 million, are under the threat of high damage due to possible big earthquakes and also triggered tsunamis. Recent studies and evaluation of earthquake recurrence periods revealed that there is a high possibility of having an earthquake with a magnitude larger than M_w7.0 in the Prince Island Fault (PIF). According to Ambraseys (2002), the last earthquake on this fault system occurred in 1766 and since that time this fault has been accumulating a huge amount of energy. According to Parsons (2004), the probability of occurrence of a M>7 earthquake beneath the Marmara Sea was estimated to be 35 %–70 % in the following 30 years. The region has distinctive characteristics in terms of its complex tectonic structure and the high possibility of an earthquake occurrence with a magnitude larger than 7.0 offshore of Istanbul. Therefore, there has been a wide range of studies in the Marmara Sea region regarding the fault mechanisms, seismic activities, earthquakes, and triggered tsunamis (Armijo et al., 2002, 2005; Okay et al., 1999; Le Pichon et al., 2001; Yaltırak, 2002; McNeill et al., 2004; Aksu et al., 2000; Imren et al., 2001; Pondard et al., 2007; Yalçıner et al., 1999, 2000, 2002; Aytore et al., 2016; Hébert et al., 2005; Altınok et al., 2003, 2011; Guler et al., 2015; Cankaya et al., 2016; Tufekci et al., 2018; Latcharote et al., 2016).

The North Anatolian Fault Zone (NAFZ) controls a great part of the seismic activity in the Marmara Sea region. The fault zone sets apart Anatolia (Asian part of Turkey) and Eurasia due to the northward migration of the Arabian Plate in the east and southward rollback of the Hellenic subduction zone in the west as seen in Fig. 1 (Armijo et al., 1999; Flerit et al., 2004; Le Pichon et al., 2015).

Figure 1 Seismicity map of the Marmara region and general tectonic map of Turkey at the top left. In the seismicity map, the size of the circles changes with magnitude of the earthquakes, and the color of the circles defines the depth change of the earthquakes. Red lines show the known active faults (modified from Emre et al., 2013) in the region, and the white square is the area with the PIF. In the general tectonic map of Turkey, red arrows show the direction of the plate motion, black lines show the active faults in the region (modified from Emre et al., 2013), and the red rectangle shows the Marmara region (created using The Generic Mapping Tools, version 5.4.1). The duration of the catalog, used for the seismicity map, is 05 April 1905–31 December 2018 (KOERI, 2001).

The Marmara Sea region is a transition zone between the strike-slip regime of the NAFZ and the extension regime of the Aegean Sea area (top left of Fig 1). The northern branch of the NAFZ forms a major transtensional NW–SE right bend under the Sea of Marmara at the Çınarcık trough (Murru et al., 2016). The fault trace is attached to the complex Central Marmara and Tekirdağ pull-apart basins, before joining the NE–SW-striking Ganos fault on land by following the northern margin of the Marmara Sea. Finally, the fault exits into the Aegean Sea by way of Saros Gulf (Wong et al., 1995; Armijo et al., 1999, 2002; Okay et al., 1999; Le Pichon et al., 2001; Yaltırak, 2002; McNeill et al., 2004). The fault trace beneath the Marmara Sea is not directly observable. Therefore, making a segmentation model for the offshore parts of the NAFZ is quite challenging, which causes the fault dimensions, such as its length and width, to include a sum of error margin (Aksu et al., 2000; Imren et al., 2001; Le Pichon et al., 2001; Armijo et al., 2002, 2005; Pondard et al., 2007).

The current right-lateral slip rate along the NAFZ is about 25 mm yr⁻¹ (Meade et al., 2002; Reilinger et al., 2006). On the western side, the motion between the Anatolian and Eurasian plates is accommodated across the Marmara region by ∼19 mm yr⁻¹ of right-lateral slip and 8 mm yr⁻¹ of extension (Flerit et al., 2003, 2004). Slip rates of the main Marmara fault range between 17 and 28 mm yr⁻¹ (Le Pichon et al., 2003; Reilinger et al., 2006). On the other hand, Hergert and Heidbach (2010) suggest that the right-lateral slip rate on the main Marmara fault is between 12.8 and 17.8 mm yr⁻¹ due to slip partitioning and internal deformation. The right-lateral slip rate for the PIF and Çınarcık basin is 15±2 mm yr⁻¹ and in addition to this the fault has 6±2 mm yr⁻¹ of extension (Ergintav et al., 2014).

The main characteristic of the NAFZ is that the earthquakes systematically propagate westward, and historical records show that the northern strand of the NAFZ generates earthquakes with the recurrence interval of about 250 years beneath the Marmara Sea; the last event occurred in 1766 (Ambraseys, 2002; Bohnhoff et al., 2013). This event caused the rupture of the 58 km long northern part of the NAFZ from İzmit to Tekirdağ (Ambraseys and Finkel, 1995; Ambraseys and Jackson, 2000). However, the earthquake that happened on 2 September 1754 can be considered the last characteristic event for the PIF segment, and it caused the rupture of a 36 km long fault segment (Ambraseys and Jackson, 2000). The NAFZ has experienced two M>7 earthquakes: in August 1912 in Ganos and August 1999 at İzmit. After the 1999 İzmit event, seismic energy along the 150 km long northern part of the NAFZ has been accumulating continuously since the 22 May 1766 earthquake. This fault zone extends right to the south of Istanbul beneath the Marmara Sea, and this situation increases the rupture possibility of the PIF and the risk for Istanbul (Stein et al., 1997; Barka, 1999; Bohnhoff et al., 2013). Ergintav et al. (2014) also indicated that the PIF segment accumulates stress 15±2 mm yr⁻¹, and the 3.7 m slip deficit has been accumulating since the 1766 events. This makes the PIF most likely to generate the next M>7 earthquake along the Sea of Marmara segment of the NAF.

Besides these seismic activities in the region, studies of the historical tsunami records show that 35 tsunami events happened between 330 BCE and 1999 BCE in the Marmara Sea region, and the majority of them are earthquake-related tsunami events (Altinok et al., 2011; Yalçıner et al., 2002). The 1509 earthquake, with an estimated magnitude around 7.5 (Ambraseys, 2002), is one of the examples for these events. This earthquake triggered a tsunami and the tsunami waves inundated the Istanbul coast, reaching the city walls and killing around 4000–5000 people in the city (Ambraseys and Finkel, 1995). The 1894 earthquake is also one of the important events that happened in the Marmara Sea. The earthquake triggered a tsunami and the sea inundated 200 m of coast in Istanbul (Altinok et al., 2011). The most recent event happened after the 17 August 1999 İzmit earthquake, and after the earthquake E–W-trending tectonic deformation along the basin and submarine failures generated a tsunami. The International Tsunami Survey Team (Yalciner et al., 1999, 2000) investigated the region and they observed 2.66 m run-up along the coast from Tütünçiftlik to Hereke and 2.9 m run-up at Değirmendere (Yalçıner et al., 2002).

Several tsunami hazard estimation studies (Ozer Sozdinler et al., 2020; Hancilar, 2012; Aytore et al., 2016; Hébert et al., 2005) were also conducted in the region. These tsunami analyses were mostly performed in a deterministic manner using various earthquake scenarios depending on the combinations of different fault parameters without considering probability of occurrences. The 40 km long fault in the eastern basin of the Marmara Sea, with a significant normal component, may generate tsunami waves which can reach maximum 2 m heights along the Istanbul coast with considerable local inundation (Hébert et al., 2005). The rupture of the Yalova Fault, PIF, or Central Marmara Fault can also cause serious damage along the coast of Istanbul. Tsunami wave heights can reach 4.8 m and can penetrate 340 m inwards from the coast in Haydarpaşa Port (Aytore et al., 2016).

A few probabilistic seismic and tsunami hazard analyses (Murru et al., 2016; Erdik et al., 2004; Hancilar, 2012) were also performed in this region. Seismic hazard maps were prepared in the Marmara Sea region by describing fault segments and peak ground accelerations with the periods corresponding to 10 % and 2 % probabilities of exceedance in 50 years (Erdik et al., 2004). In addition, tsunami inundation maps are prepared based on probabilistic and deterministic analyses by depending on these segmentations (Hancilar, 2012). Time-dependent and time-independent earthquake ruptures are also estimated in the Marmara Sea region for the next 30 years (Murru et al., 2016). These previous studies have been conducted for entire Marmara Sea region and therefore they give general and rough information about probability of occurrence in the region without focusing on any specific region in high resolution. However, probabilistic tsunami hazard assessment is important to calculate the tsunami exposure and risk on human populations and infrastructures, since probability calculations consider all possible earthquakes in a fault even if they occur with very low probability (Løvholt et al., 2012, 2015; Grezio et al., 2017). The results of probabilistic studies should be considered when decision makers design coastal zones and structures, especially critical ones. Different from previous probabilistic approaches in the Marmara Sea, the probability of earthquake occurrences in one fault segment, PIF, are taken into account for the preparation of high-resolution tsunami inundation maps and distribution of hydrodynamic parameters due to the probability of occurrence of associated earthquakes on the PIF determined by Monte Carlo (MC) simulations.

This probabilistic tsunami hazard analysis (PTHA) study depends on the fully characteristic fault model, and the main purpose is to perform PTHA for selected test sites. The Tuzla test site is one of the coastal districts of Istanbul and located on the southernmost part of the city (Fig. 2). The region includes several residential areas, but the most critical point about the region is that Tuzla has the biggest shipyard area not only in the Marmara Sea but also in Turkey (Fig. 3). In this study we mainly focused on this region because it is about 20 km away from the PIF and therefore has a high risk of both earthquake and tsunami damage.

Figure 2 The Marmara Sea region, Tuzla test site, and the location of the PIF segment which is used in the analysis like a straight line (created using ArcMap version 10.5).

Figure 3 Some important locations at the Tuzla domain. (a) Northern part of the Tuzla domain. (b) Southern part of the Tuzla domain. (c) Tuzla shipyard (created using ArcMap version 10.5).

2 Probabilistic analysis

Probabilistic tsunami hazard analysis (PTHA), as it is recently becoming a widely used procedure for coastal zones, is performed for the Tuzla region, Istanbul. This method has been applied for various tsunami sources, such as earthquakes, landslides, volcanic activities, etc., on various scales, local, regional, and global (Grezio et al., 2017). For earthquake-generated tsunamis, the method is generally adapted from seismic hazard assessment methods (González et al., 2009). Such kinds of studies consider the events that are generated by coseismic seafloor displacement, using seismic probabilistic hazard analysis (SPTHA), but numerous tsunami simulations are required to consider all expected combinations of seismic sources. This problem can be solved by applying a simplified event tree approach and a two-stage filtering procedure to reduce the number of required source scenarios without decreasing the quality and accuracy of inundation maps (Lorito et al., 2015). The earthquake source itself is very uncertain and the investigation of this uncertainty can be carried out by building an event tree instead of using a logic tree and hazard integrals (Selva et al., 2016). The logic tree approach can be applied to the generation of tsunami hazard curves to decrease the uncertainties by including branches, which are the combination of tsunami sources, magnitude distribution of characteristic tsunamigenic earthquakes, their recurrence interval, and the tsunami height estimation procedure based on a numerical simulation (Annaka et al., 2007). For regional studies, hazard curves can be generated by empirical analysis using available tsunami run-up data. However, if such data are not available, MC simulations, a computationally based method widely used in probabilistic seismic hazard analysis (PSHA), can be considered as a primary method to generate tsunami hazard curves (Geist and Parsons, 2006; Horspool et al., 2014). Submarine landslides, on the other hand, are the major tsunami source for passive margins, which are the transition zone between the oceanic and continental lithosphere that is not an active plate boundary, and they have been included in PTHA methodologies (Geist and Lynett, 2014). Probabilistic studies are also applied to develop multi-hazard loss estimation methodology for coastal regions that are exposed to cascading shaking-tsunami hazards due to offshore mega-thrust subduction earthquakes (Goda and De Risi, 2018).

In this study, a characteristic earthquake model is used to estimate the earthquake recurrence on the PIF. Paleoseismologic studies (Ryall et al., 1966; Allen, 1968; Schwartz and Coppersmith, 1984) suggest that an individual fault tends to generate characteristic earthquakes having a very narrow range of magnitudes. These individual faults have a different frequency distribution than the log linear Gutenberg–Richter frequency–magnitude relationship (Aki, 1984; Schwartz and Coppersmith, 1984; Youngs and Coppersmith, 1985). According to Aki (1984), a characteristic earthquake is generated as a result of constancy of barriers to rupture through repeated seismic cycles.

PIF is fully characteristic and a characteristic earthquake will rupture an entire fault as a whole and release all the energy. Therefore, while performing MC simulations, the area of the fault and fault parameters (strike, dip, and rake angles) are used as constants referring to the outcomes of EU Seventh Framework Programme project MARSITE (Ozer Sozdinler et al., 2020). One of the work packages of this project aimed to define the geometry of the possible tsunamigenic faults in the Marmara Sea and 30 different earthquake scenarios with the different rupture combinations of 32 possible fault segments. Based on these 30 different earthquake scenarios, tsunami numerical modeling is performed. The definition of fault segments depends on extensive review of the literature (Alpar and Yaltırak, 2002; Altınok and Alpar, 2006; Armijo et al., 2005; Ergintav et al., 2014; Gasperini et al., 2011; Hébert et al., 2005; Hergert et al., 2011; Hergert and Heidbach, 2010; Imren et al., 2001; Kaneko, 2009; Le Pichon et al., 2001, 2003, 2014; Oglesby and Mai, 2012; Şengör et al., 2014; Tinti et al., 2006; Utkucu et al., 2009). As a result of this review, each fault segment is defined as a rectangular area with hypothetical uniform slip. According to the results of the project, the fault parameters of the PIF are given in Table 1. The 3-D fault configuration given by Armijo et al. (2002), which explains fault segmentation in the region depending on morphology, geology, and long-term displacement fields, also fits with the PIF parameters that are used in the project. These parameters are used as constants in this study while assessing probability of occurrence of each earthquake to allow full fault rupture at different depths with different magnitudes.

Table 1 The area and the focal mechanism of the PIF zone. These are the constant parameters during the MC simulation application.

The MC simulation technique is generally applied to generate an earthquake catalogue of a given length of time. In this technique, a list of earthquakes can be generated using the frequency–magnitude relationship for each seismic source (Zolfaghari, 2015). Seismic zonation should be performed by considering regions that have relatively homogeneous earthquake activity and faulting regimes (Sørensen et al., 2012). In this study, the fault segment model proposed in Ozer Sozdinler et al. (2020) is used and PIF is the only segment that is a seismic source. After that, tsunami numerical modeling is performed for each event of this synthetic catalogue, and tsunami hydrodynamic parameters, mainly maximum wave heights, inundation depth, current velocities, and tsunami inundation zones, are estimated. Tsunami risk assessment will serve the needs of societies best when regional studies are associated with the local ones (Sørensen et al., 2012).

The MC simulation technique allows the generation of a list of earthquakes based on a frequency–magnitude relationship. This technique depends on a uniformly distributed source model and it provides an equal likelihood to each earthquake source. As a result, the synthetic earthquake catalogue will have uniformly randomly distributed earthquake sources (Zolfaghari, 2015).

Using MC simulation, a synthetic earthquake catalogue is generated by selecting earthquake magnitude and depth as uniformly distributed random numbers in a given interval and using area and directivity of the fault as a constant variable (Table 1). We performed MC simulations 100 times for 100 different earthquake scenarios. The number of earthquakes in the catalog is selected as a reasonable number that represents the number of iterations randomly performed in MC simulations for having a synthetic earthquake scenario. As mentioned earlier, NAFZ generates an earthquake with the recurrence interval of about 250 years beneath the Marmara Sea. Therefore, selecting 100 earthquake scenarios would cover a time period of 100×250 yr = 25 000 yr, which is considered as an adequate catalog duration in this study. However, because of having time-dependent probabilistic analyses, this catalog duration is not used for PTHA in this study.

Earthquake magnitude is one of the parameters randomly selected by the MC technique. Based on a characteristic earthquake model, individual faults tend to rupture the entire fault when a large earthquake occurs. This model assumes that a characteristic earthquake releases all of the seismic energy during the fault rupture, and the magnitude of the earthquake depends on the dimension of the fault (Abrahamson and Bommer, 2005).

As mentioned previously, only the PIF is considered an earthquake source approximately 34 km in length and 14 km in width (Ozer Sozdinler et al., 2020; Karabulut et al., 2002). This fault zone is assumed to have the potential to generate a characteristic earthquake and rupture the entire fault. According to the Wells and Coppersmith (1994) scaling relation between fault area and magnitude (Eq. 1), this fault can generate a characteristic earthquake with magnitude varying between M_w6.5 and 7.1.

$$\begin{array}{}\text{(1)}& M_{\mathrm{w}}=a+b\cdot \log (L\cdot W)\end{array}$$

In this equation, a and b are coefficients, which are 4.33 and 0.9, respectively, L is fault length, and W is the fault width.

Displacement on the fault surface calculations is carried out for each randomly selected magnitude using the formulation of Aki (1966),

$$\begin{array}{}\text{(2)}& D={\displaystyle \frac{M_{\mathrm{0}}}{\mathit{\mu }A}}={\displaystyle \frac{{\mathrm{10}}^{\left(M_{\mathrm{w}}+\mathrm{6.07}\right)\cdot \mathrm{1.5}}}{\mathit{\mu }A}},\end{array}$$

where D is displacement on the fault surface, M_w is moment magnitude, μ is the shear modulus (μ=30 GPa), and A is the fault area.

Seismogenic thickness and the location of the earthquake is another important parameter required for earthquake and tsunami source. At first, the PIF zone is accepted as fully characteristic and an earthquake should rupture the entire fault area. Therefore, it is assumed that if the rupture starts at the center of the fault and continues in both directions, the fault will rupture entirely. For this reason, the locations of the earthquakes are accepted as the midpoint of the PIF zone for each earthquake scenario (Ozer Sozdinler et al., 2020).

For the seismogenic thickness, the seismic activity of the northern segment of NAFZ starts at the depth of 5 km (Karabulut et al., 2003). The bottom of the seismogenic thickness can be determined based on the aftershock activity of the 17 August 1999 İzmit earthquake. The earthquakes on the northern scarp of the Çınarcık basin are observed between the depths of 5 and 14 km. The mechanism of events between the depth of 5 and 10 km shows the behavior of normal faulting. On the other hand, the strike-slip mechanism dominates the depths below 10 to 14 km. As a result, seismic activity can be observed between the depths of 5 and 14 km, and fault plane solutions show normal and strike-slip mechanisms in this area (Karabulut et al., 2002). Therefore, the depth of events varies between 5 and 14 km in MC simulations.

In time-independent earthquake occurrence models, probability of an event occurrence follows a Poisson distribution in a given period of time. Therefore, the result of this model does not vary in time. However, probability of an earthquake occurrence is based on the time that has passed since the occurrence of the last event and it follows a Brownian passage time (BPT), lognormal, or other probability distribution (Matthews et al., 2002; Ellsworth et al., 1999; Davis et al., 1989; Rikitake, 1974). In this model, in addition to the recurrence time of earthquakes, variability of the frequency of events and the elapsed time from the last characteristic event are the additional required information and the longer elapsed time causes an increase in probability of an event occurrence (Cramer et al., 2000; Petersen et al., 2007).

Calculation of probability in multi-segment ruptures and more complicated models includes the Gutenberg–Richter magnitude–frequency relationship (Gutenberg and Richter, 1944). The application of time-dependent models is based on a characteristic earthquake model, which assumes all large events occurring along a particular fault segment would have similar magnitudes, rupture area, and average displacements (Schwartz and Coppersmith, 1984). Therefore, this model is suitable for calculating the probability of occurrence of an earthquake on a single fault.

It should be noted that, in this study, PIF is considered to be the only source for the earthquake and tsunami. A time-dependent probabilistic model is followed for the probability calculations because this probabilistic model allows us to consider only one fault instead of using multi-segment rupture scenarios through a characteristic earthquake model.

In the time-dependent approach, the BPT probability model is used to obtain the recurrence time probability of the earthquake in the fault segment. This model does not show a significant difference with the lognormal distribution except for consideration of very long elapsed times from the last characteristic event (Petersen et al., 2007). A characteristic event occurs when the load-state process reaches the failure threshold; an earthquake releases all energy loaded on the fault and then starts the new failure cycle. The time interval between consecutive earthquakes shows a Brownian passage time distribution and that can be useful to forecast long-term seismic events by generating a time-dependent model (Matthews et al., 2002). The Working Group on California Earthquake Probabilities (1999) and the Earthquake Research Committee (2001) have already implemented this time-dependent approach in the San Francisco Bay area and Japan, respectively, for the prediction of long-term events (Petersen et al., 2007). This model depends on the time period passed since the last characteristic event and recurrence time of the earthquake. The probability density function for the BPT model (Matthews et al., 2002) is given by

$$\begin{array}{}\text{(3)}& f\left(t,T_{\mathrm{r}},\propto \right)={\left({\displaystyle \frac{T_{\mathrm{r}}}{\mathrm{2}\mathit{\pi }{\mathit{\alpha }}^{\mathrm{2}}{t}^{\mathrm{3}}}}\right)}^{\mathrm{1}/\mathrm{2}}\exp \left({\displaystyle \frac{{\left(t-T_{\mathrm{r}}\right)}^{\mathrm{1}/\mathrm{2}}}{\mathrm{2}{T}_{\mathrm{r}}{\mathit{\alpha }}^{\mathrm{2}}t}}\right),\end{array}$$

where t is the elapsed time from the last characteristic event and α is the aperiodicity (also known as the coefficient of variation). Aperiodicity defines the regularity of the expected characteristic earthquakes on the fault and varies between 0.3 and 0.7. This parameter, which is known as the parameter defining how much an expected characteristic earthquake occurs regularly or irregularly on any fault segment (Murru et al., 2016), was taken as 0.5 in this study (Parsons, 2004). The mean recurrence interval of earthquakes (T_r) can be defined as the ratio between the mean moment of repeating earthquakes (seismic moment) and the long-term moment accumulation rate on the fault (moment rate) (Ren and Zhang, 2013). Seismic moment can be obtained using the formulation of Kanamori (2004), and the moment rate of the fault is calculated from fault area and long-term slip rate of the fault (WGCEP, 2003).

$$\begin{array}{}\text{(4)}& T_{\mathrm{r}}={\displaystyle \frac{M_{\mathrm{0}}}{\dot{M_{\mathrm{0}}}}}={\displaystyle \frac{{\mathrm{10}}^{\left(M_{\mathrm{w}}+\mathrm{6.07}\right)\cdot \mathrm{1.5}}}{\mathit{\mu }VA}}\end{array}$$

In this equation, M_w is moment magnitude, μ is the shear modulus, V is long-term slip rate in millimeters per year, and A is the fault area. The moment magnitude value in Eq. (4) was selected randomly using MC simulations. Thus, seismic moment (M₀) and the mean recurrence time (T_r) were calculated for each earthquake scenario. Long-term slip rate is also selected as 17 mm yr⁻¹ for this equation (Ergintav et al., 2014).

Probability of the earthquake occurrence on the fault is calculated based on the probability density function approach. The probability of occurrence of an event in the next ΔT years, given that it has not occurred in the last t years, is given by (Erdik et al., 2004)

$$\begin{array}{}\text{(5)}& P(t,\mathrm{\Delta }T)={\displaystyle \frac{\underset{t}{\overset{t+\mathrm{\Delta }T}{\int }}f\left(t\right)\mathrm{d}t}{\underset{t}{\overset{t+\mathrm{\infty }}{\int }}f\left(t\right)\mathrm{d}t}}.\end{array}$$

In this case, probability of a characteristic earthquake was calculated using ΔT as 50 and 100 years.

Figure 4 Nested domains for tsunami numerical modeling. Red rectangles show the limits of these domains. Grid size of these domains have a certain 1:3 ratio between each other (created using ArcMap version 10.5).

3 Tsunami numerical modeling

Tsunami simulations are performed for each earthquake in the synthetic catalogue using the tsunami numerical model NAMI DANCE (NAMI DANCE, 2011). The code is the user-friendly version of TUNAMI-N2 (Imamura et al., 2001) developed in C$++$ language, which computes all fundamental parameters of tsunami motion in shallow water and in the inundation zone. It uses an explicit numerical solution of shallow water wave equations with the finite-difference technique and allows for better understanding of the effect of the tsunami waves (Shuto et al., 1990; Imamura, 1989). NAMI DANCE can solve both linear and nonlinear shallow water (NSW) equations with a selected coordinate system (Cartesian or spherical) and calculates the tsunami motion. Linear shallow water (LSW) equations are preferable in deep water because of reasonable computer time and memory, and they calculate the results at an acceptable error limit (Insel, 2009). NAMI DANCE is validated and verified using NOAA standards and criteria for tsunami currents and inundation (Synolakis et al., 2007, 2008). The numerical solutions of NAMI DANCE are also tested, validated and verified against analytical solutions, laboratory measurements, and field observations (NTHMP, 2015; Lynett et al., 2017; Velioglu, 2009).

NAMI DANCE calculates tsunami generation using Okada (1985) equations. In this study, water surface distribution of tsunami source (initial wave amplitude) is calculated with this method for 100 earthquakes of the synthetic earthquake catalogue prepared by MC simulations. As an example, Fig. 5 shows the initial water surface calculated due to one of 100 tsunami sources generated by MC simulations (Fig. 5).

Before starting tsunami simulations, the necessary inputs should be prepared precisely in order to obtain reliable results. Bathymetry–topography data are one of the most important inputs in NAMI DANCE that significantly effects the reliability of results, especially in the shallow water zone due to the nature of the NSW equations. NAMI DANCE can perform nested analyses under the condition that the grid sizes of the study domains have a certain 1:3 ratio between each other. Therefore, we generated four nested domains having the coarsest grid size as 81 m and the finest grid size as 3 m with a 1:3 ratio in the GIS environment. Bathymetric data for the biggest domain are the combination of the 30 arcsec resolution General Bathymetric Chart of the Oceans (GEBCO) and data produced by navigational charts in shallow zones. Topographic data, on the other hand, are at a high resolution, which is obtained from the Department of Housing and Urban Development of Istanbul Metropolitan Municipality digital elevation model (DEM) and vector data with resolution of 5 and 1 m, respectively. The bathymetry–topography data in the smaller domains are the downscaled version of the 81 m grid bathymetry–topography data; however high-resolution digitized coastline and sea and land structures are also included in the data to generate the smallest grid domain of 3 m (Fig. 4).

Figure 5 Initial water surface distribution of one of the 100 tsunami sources. The red frame indicates the boundary of the largest nested domain with an 81 m grid size (plotted using ArcMap version 10.5).

The synthetic gauge point file is another required input of NAMI DANCE. In addition to the calculation of principal tsunami hydrodynamic parameters, the program can also calculate the change of water level, current velocity, and flow depth over time in every gauge point. Therefore, various gauge points are selected along the coast of nested domains, nearshore and offshore and close to some critical structures on land.

During the inundation of tsunami waves, current velocity is an important tsunami parameter in land and sea, especially in ports and bays. Strong current velocities may cause sea vessels to be dragged offshore by undertow or to ground inland. This parameter as well as tsunami wave amplitude, inundation depth, and Froude number can be calculated by NAMI DANCE. However, in this study, the results are represented based on only the probability of exceedance of threshold values for water surface elevation and inundation depth.

4 Results and discussion

In this study, tsunami hydrodynamic parameters are calculated in both the coarsest domain (whole Marmara Sea) and finest domain (Tuzla region). The main parameters focused in this study are the tsunami wave heights and inundation depths, and the results are shown in terms of probability of exceedance of threshold wave height and inundation depth values within the next 50 and 100 years. The situation for the next 500 years is not considered because the return period of the fault rupture is about 250 years, which means this fault generates at least one earthquake within the next 500 years. In other words, probability of exceedance for the next 500 years will be about 99 %.

We present the results of the PTHA for the Tuzla test site in terms of three different visualization categories for the next 50 and 100 years. First, distribution of probability of occurrence of the tsunami hydrodynamic parameters, which are minimum and maximum water surface elevation and inundation depth, is shown. Second, tsunami inundation maps that show the probability of exceedance of 0.3 m inundation depth for different time periods are generated for the Tuzla region in order to observe flooded areas and their probabilities clearly. Finally, the probability map of exceedance of 0.3 m wave heights at synthetic gauge points is represented as a bar chart.

4.1 Probability of exceedance for the entire synthetic earthquake catalogue

The graphics are generated to demonstrate the probabilities of occurrences corresponding to the minimum and maximum water surface elevations and inundation depth calculated from tsunami sources of each earthquake in the synthetic earthquake catalogue.

It should be noted that in the case of having same magnitude of earthquakes in two different earthquake scenarios of the catalogue, the probability of occurrences of these scenarios would be the same. However, since they would have different focal depths, the tsunami initial wave height calculated by Okada (1985) will be different, which results in the calculation of different hydrodynamic parameters. As a result, the graphs show different maximum water surface elevations having the same probability of occurrences.

Figure 6 Probabilities of exceedance corresponding to maximum water surface elevation, minimum water surface elevation, and inundation depth for the next 50 years. Black dots represent the probability of exceedance of the tsunami hydrodynamic parameter for each event in the catalog. The blue line is the best fit curve to the data and the dashed blue line is the 95 % confidence boundary of the fitted curve. The residual of the fit is represented for each probability curve.

In Fig. 6, graphics of probabilities of occurrences according to maximum and minimum water surface elevation (maximum water withdraw) and inundation depth for the next 50 years are represented. According to these graphs, tsunami wave heights up to 1 m and withdrawal of the waves around 1 m have approximately 65 % ± 15 % probability of occurrence. The Tuzla region includes various shipyards, ports, and other important facilities. Therefore, the probability of the withdrawal of the water is as important as maximum water surface elevation. The 1 m height of wave withdrawal may cause the ships to be stranded at the ports and results in extreme financial losses as observed in the 20 July 2017 Bodrum–Kos earthquake and tsunami (Yalçıner et al., 2017). The probability of having 1 m inundation depth, on the other hand, can be predicted as about 60 % ± 10 %. The residual of probability with respect to the fitted curve for each data point is demonstrated right after the percentage of probability with the ± sign.

The situation for the next 100 years (Fig. 7) obviously shows that probability of occurrences would increase with time. The probability of exceedance of 1 m water surface elevation and 1 m wave withdrawal reaches up to 85 % ± 10 %. Probability of exceedance of inundation depth also changes significantly. The probability of exceedance of 1 m inundation depth is found to be around 80 % ± 10 %.

Figure 7 Probabilities of exceedance corresponding to maximum water surface elevation, minimum water surface elevation, and inundation depth for the next 100 years. Black dots represent the probability of exceedance of the tsunami hydrodynamic parameter for each event in the catalog. The blue line is the best fit curve to the data and the dashed blue line is the 95 % confidence boundary of the fitted curve. The residual of the fit is represented for each probability curve.

Considering the results of the whole simulation, the worst-case earthquake scenario generated tsunami waves with maximum water surface elevation equal to 1.8 m, minimum water surface elevation (maximum withdraw) equal to −2.1 m, and inundation depth equal to 1.6 m. The probability of occurrence of this event is 35 % for the next 50 years and 60 % for the next 100 years.

4.2 Probabilistic tsunami inundation maps for the Tuzla test site

Inundation maps of the Tuzla domain are also prepared for the next 50 and 100 years in the GIS environment. Even if inundation depth is on the order of a few centimeters, it can lead to people being dragged by undertow in coastal regions due to the high current velocities of the waves (Jonkman and Penning-Rowsell, 2008). Therefore, these inundation maps have a great significance for understanding the flooded areas in the study domain and the amount of water penetrated inland.

Generation of inundation maps is based on the probability of exceedance of 0.3 m inundation depth. There are several studies in the literature proving both experimentally and numerically that tsunami waves with an order of 0.3 m height have the potential to crush a human body (Jonkman and Penning-Rowsell, 2008; Takagi et al., 2016). For this reason, only the earthquake scenarios that generated inundation depths larger than or equal to the 0.3 m threshold value are considered.

Inundation depth files, which are one of the outputs of the NAMI DANCE, are used for the calculation.

The inundation depth values at each grid node are replaced with the probability of occurrence of the respective earthquake scenario. We repeated this procedure for all earthquake scenarios, which has inundation depths larger than or equal to the 0.3 m threshold.

The mean (average) probability of occurrence is calculated at each grid node. Thus, the spatial distribution of probability of exceedance of 0.3 m inundation depth in the inundation zone is obtained for a specific time interval (Fig. 8).

Figure 8 Probabilistic tsunami inundation maps for the Tuzla study domain representing the probability of exceedance of 0.3 m inundation depth within the next 50 and 100 years. Change of colors from green to red represents the increasing probability of exceedance (created using ArcMap version 10.5).

Figure 8 shows the inundation maps of Tuzla shipyard for the next 50 and 100 years. Most of the area in the Tuzla shipyard region has a probability of exceedance between 10 % and 20 % for the next 50 and 100 years. However, some places in the northern and southern parts of the area and inside the bay show larger than 75 % probability of inundation within the next 100 years. Maximum inundation distance is observed at around 60 m at the test site.

In Fig. 9, probabilistic inundation maps of one of the most important facilities in the study region are represented for the next 50 and 100 years. The area has high potential to be exposed to tsunami waves with a probability larger than 50 % for the next 50 years. In 100 years, this probability increases and varies between 75 % and 90 %. No significant inundation zone is observed along the coast of the seawall and the peninsula. This may be due to the high ground elevation of these zones. Tsunami waves are inundated up to 45 m inside the small bay. This inundation distance could cause severe damage to the shipyard and other constructions if corresponding current velocities are also significant.

Figure 9 Probabilistic tsunami inundation maps of the northern part of the Tuzla study domain representing the probability of exceedance of 0.3 m inundation depth for the next 50 and 100 years. Change of colors from green to red represents the increasing probability of exceedance (created using ArcMap version 10.5).

In the next figure (Fig. 10), the southern part of the Tuzla shipyard is seen according to probabilities of inundation for the next 50 and 100 years. Very limited area in the coastal zone is inundated with the probability between 30 % and 50 % within the next 50 years. The probability decreases up to 10 % at some inner locations from the coastline. For 100-year recurrence time, the situation is almost the same. Only minor parts of the region in the south approach the 75 %–90 % probability of exceedance of 0.3 m inundation depth threshold. The maximum inundation distance is calculated at about 60 m. The inundated region does not include any important facility or structure, and the effect of the tsunami will be minimal. The inundation distance decreases to 10 m in the other parts of the region.

Figure 10 Probabilistic tsunami inundation maps for the southern part of the Tuzla study domain representing the probability of exceedance of 0.3 m inundation depth for the next 50 and 100 years. Change of colors from green to red represents the increasing probability of exceedance (created using ArcMap version 10.5).

The region indicated in Fig. 11 is located inside the bay and includes a large part of the shipyard area. This area includes lots of large and small piers and ship construction facilities. The situation is more or less the same as the previous region (Fig. 9). The probability of having larger than 0.3 m inundation depth changes between 30 % and 50 % within the next 50 years, while only a few places show 75 %–90 % probability for the next 100 years along the coast. Moreover, the maximum inundation distance is calculated as 25 m for this zone. Even if the probability of inundation is low, these zones should be taken into consideration before constructing a new structure.

Figure 11 Probabilistic tsunami inundation maps of the shipyard area in the Tuzla study domain representing the probability of exceedance of 0.3 m inundation depth for the next 50 and 100 years. Change of colors from green to red represents the increasing probability of exceedance (created using ArcMap version 10.5).

4.3 Synthetic gauges

Finally, the probability of exceedance of 0.3 m wave heights at synthetic gauge points is presented by bar charts to consider the nearshore effect of tsunami waves along the western coast of Istanbul. Because of the closeness to the fault zone, the southeast coasts of the city are under threat of significant tsunami damage. Similar to the method applied during the preparation of probabilistic inundation maps, the earthquake scenarios with wave heights at synthetic gauge points larger than or equal to 0.3 m are selected and replaced with the probability of each scenario according to wave heights, and after that the average probabilities at each synthetic gauge point are obtained accordingly.

Figure 12 demonstrates the probability of exceedance of 0.3 m wave height at synthetic gauge points, which are about 350 m apart from each other, along the western coast of Istanbul within the next 50 and 100 years. The probability increases while color scale changes from green to purple. According to this figure, minimum probability of exceedance is shown as 75 % at some points. Except for a few of the 228 synthetic gauge points, all points have larger than 90 % probability of exceedance of 0.3 m wave height within the next 50 years.

This condition is very serious since there are so many residential areas and important spots such as ports and recreational facilities in this region. The minimum probability of occurrence, which can generate tsunami waves with at least 0.3 m wave heights, reaches up to 90 % for the next 100-year time period. However, 95 % probability of exceedance of 0.3 m wave height dominates the region for this timescale.

Figure 12 Probability of exceedance of the 0.3 m tsunami wave height within the next 50 and 100 years at synthetic gauge points. The yellow rectangles show the Tuzla study domain; change of colors from purple to green on the bars represents the decreasing probability of exceedance (created using ArcMap version 10.5).

4.4 Uncertainties

PTHA studies include some uncertainties because of the rare occurrence of the large events. Quantification of these uncertainties generally includes the mixture of empirical analyses and subjective judgment.

Uncertainties of PTHA can be divided into two: as aleatory and epistemic variability. Aleatoric uncertainty is the natural randomness of the physical process. Including more data in the analyses does not contribute to the reduction of the aleatoric uncertainty. However, knowledge about the modeling process may decrease this unpredictability. The occurrence time of the earthquake is one of the most fundamental aleatory variables in PTHA. This parameter is generally assumed to be a time-independent variable. However, in this study we used a time-dependent probability model, which reduces the uncertainty on this parameter. The mechanism of the source is considered to be another aleatory variable for PTHA studies. The majority of earthquakes around the world occur at well-defined plate boundaries. However, some unidentified low-activity intraplate faults exist, which were recently included in PTHA studies (Selva et al., 2016). Moreover, the fault volume, which is used in scaling relations to calculate the source magnitude, is another aleatory term. Although homogenous slip distribution is a common implementation in PTHA, slip distribution of large events does not show homogenous behavior. Therefore, definition of asperities on the fault is another aleatoric variable which should be considered. Tsunami numerical modeling, itself, is also another aleatoric variable since they do not show correlation with real observations, which are more variable than earthquake scenarios incorporated in PTHA (Grezio et al., 2017). The aleatory variable affects the results because it is incorporated directly into the hazard calculations (Abrahamson and Bommer, 2005).

Epistemic uncertainty, on the other hand, consists of the lack of knowledge of the physical process and data. Segmentation of the fault system is one of the epistemic variables since it is not certain where the rupture will be generated and which segments will be triggered. In addition, there are many different scaling relations, which cause another epistemic uncertainty, between the fault area and magnitude. It is also important for tsunami generation whether the fault rupture reaches the surface or not. Thus, updip and downdip limits of the fault rupture can be considered another epistemic variable (Grezio et al., 2017). Accurate probability distributions of input cannot be known. For example, one cannot assume that probability of occurrence of an event follows Poisson distribution. However, return periods of events do not simply fit this distribution (Gonzalez et al., 2013). Unlike aleatoric uncertainty, epistemic uncertainty can be decreased when more information is available (Godinho, 2007). Different techniques, such as logic tree, the Bayesian method, etc., have been developed to reduce these uncertainties.

In this study, a probabilistic model is established based on the characteristic fault model of PIF, which is a segment of NAF, one of the best studied fault zones in the world. It is also assumed that the entire fault area is ruptured, reaching the surface and generating a homogenous slip for each event. The maximum magnitude range of the fault is calculated with the Wells and Coppersmith (1994) scaling relation. All these assumptions naturally include uncertainties which are naturally reflected in this PTHA study. In addition, MC simulation itself also includes uncertainty as being performed 100 times to create synthetic earthquake scenarios. The effect of uncertainty in the aperiodicity parameter also exists and can be reduced by including different parameters for MC simulation. Therefore, the tsunami hydrodynamic parameters associated with the probability of occurrence of the corresponding scenario preserve the same uncertainty.

5 Conclusion

In this study, time-dependent PTHA is performed in the Tuzla region of Istanbul for the purpose of understanding the probability of having tsunami inundation after the PIF rupture. The study combines tsunami numerical modeling with a probabilistic approach, which is modified by probabilistic seismic hazard analysis. Probability calculations have been done based on the time-dependent BPT model, which depends on the time period passed since the last characteristic event and the recurrence time of the earthquake. After that, the synthetic earthquake catalogue is generated using the MC simulation technique, and tsunami numerical modeling was performed depending on this earthquake catalogue using NAMI DANCE code in a GPU environment.

Results of this PTHA study were presented in three different ways for the next 50 and 100 years. The first one was the graphs showing the change of probability with the maximum and minimum water surface elevation and inundation depth for different time intervals. Secondly, the probabilistic tsunami inundation maps are generated for the Tuzla region. Finally, the probability maps of exceedance of 0.3 m wave heights at synthetic gauge points are represented with bar charts.

The main results of this study can be summarized as follows.

• According to the distribution of probability with respect to tsunami hydrodynamic parameters, the probability of exceedance of 1 m maximum positive and negative water surface elevation is 65 % within the next 50 years. The probability for 1 m inundation depth is 60 %.

• Considering probabilities for the next 100 years, 85 % probability of exceedance of 1 m was calculated. For 1 m inundation depth, probability of exceedance of about 80 % is obtained.

• As a result of the whole simulation, 1.8, −2.1, and 1.6 m were calculated for maximum and minimum water surface elevation and inundation depth, with the probability of 35 % for the next 50 years and 60 % for the next 100 years.

• Inundation maps indicate that inundation of tsunami waves that are equal to or larger than 0.3 m have probability mostly higher than 10 % and 20 % for the next 50 years and 100 years, respectively. The probability of occurrence of 0.3 m inundation depth was calculated as a maximum of 75 % for the next 100 years. Maximum inundation distance is calculated as 60 m and observed in the southern part of the finest 3 m grid-sized study area.

• Probabilistic results for the exceedance of 0.3 m wave height at synthetic gauge points demonstrate that only a few of them have a probability between 75 % and 85 %; however several points have more than 90 % probability for the next 50 years. Probability of exceedance increases by more than 95 % for the next 100 years.

The tsunami impact of the PIF rupture along the Tuzla coast is very important as proposed by the results of this study. However, as further steps of this study, PTHA can be done for the other critical test sites along the Marmara Sea that are close to the PIF segment. In addition, it is also advantageous to consider the other fault segments, with their various rupture combinations and complex rupture probabilities in Marmara Sea, as further studies. Previously in the framework of the MARSITE project, tsunami arrival times and maximum wave amplitudes were calculated along the coast of the Marmara Sea using different earthquake scenarios, and a tsunami scenario database was obtained with a deterministic approach (Ozer Sozdinler et al., 2020). Results of this study show that arrival time of tsunami waves is very short in Marmara Sea for most of the scenarios, which complicates the tsunami early warning operations and evacuation actions. However, due to the short arrival times of the first tsunami waves along the Marmara coast, the tsunami inundation scenario databases would be of great importance in such conditions. It would be the best option for the decision makers and civil protection authorities to also have the inundation maps prepared with a probabilistic approach in order to realize the possibility of exceedance of selected threshold inundation depth for certain critical coastal locations.

This study shows a methodology for PTHA with a time-dependent probabilistic model using only one fault (PIF) as the earthquake and tsunami source. Furthermore, this study can be developed including some faults connected to the PIF in both time-dependent and time-independent probability calculations, and BPT probability can be combined with static Coulomb stress changes on the faults. The BPT model can also be improved by including different aperiodicity parameters. The probability of occurrence of earthquakes is the main focus of this study to perform tsunami hazard analyses. However, submarine landslides are other critically important sources for tsunami generation in the Marmara Sea. Probabilities of sliding areas and the sliding volumes can be considered in the analyses. Submarine landslide-generated tsunamis can be coupled with the earthquake-triggered tsunamis in order to obtain integrated PTHA in the Marmara Sea.