Modelling earthquake rates and associated uncertainties in the Marmara Region, Turkey

Modelling the seismic potential of active faults and the associated epistemic uncertainty is a fundamental step of probabilistic seismic hazard assessment (PSHA). We use SHERIFS (Seismic Hazard and Earthquake Rate In Fault Systems), an open-source code allowing us to build hazard models including earthquake ruptures involving several faults, to model the seismicity rates on the North Anatolian Fault (NAF) system in the Marmara Region. Through an iterative approach, SHERIFS converts the slip rate on the faults into earthquake rates that follow a magnitude frequency distribution (MFD) defined at the fault system level, allowing us to model complex multi-fault ruptures and off-fault seismicity while exploring the underlying epistemic uncertainties. In a logic tree, we explore uncertainties concerning the locking state of the NAF in the Sea of Marmara, the maximum possible rupture in the system, the shape of the MFD and the ratio of off-fault seismicity. The branches of the logic tree are weighted according to the match between the modelled earthquake rate and the earthquake rates calculated from the local data, earthquake catalogue and palaeoseismicity. In addition, we use the result of the physics-based earthquake simulator RSQSim to inform the logic tree and increase the weight on the hypotheses that are compatible with the result of the simulator. Using both the local data and the simulator to weight the logic tree branches, we are able to reduce the uncertainties affecting the earthquake rates in the Marmara Region. The weighted logic tree of models built in this study will be used in a following article to calculate the probability of collapse of a building in Istanbul.

1 Introduction

The North Anatolian Fault system (NAFS) runs through the north of Turkey for a distance of more than 1500 km (Fig. 1), accommodating the westward extrusion of the Anatolian plate with a right lateral motion of around 25 mm/yr (e.g. Reilinger et al., 2006). Along the western portion of the NAFS, this right lateral motion is partitioned between two fault branches, with the northern branch accommodating most of the motion (Reilinger et al., 2006). This northern branch of the NAFS crosses the Sea of Marmara 20 km south of the city of Istanbul (Fig. 2). The last earthquakes occurring in the Sea of Marmara and severely damaging the city of Istanbul date back to 1766 and 1894 (Ambraseys, 2002). Growing from a population of less than 1 million in 1950, Istanbul is now a city of more than 15 million people. The M_w = 7.6 1999 Izmit earthquake that ruptured the fault eastward of the Sea of Marmara has been a reminder of the hazard faced by Istanbul. The NAFS is expected to be the major source of seismic hazard looming over the megalopolis. Motivated by better constraining the seismic hazard, the two decades following the 1999 earthquake have been rich in geological, geophysical and geodetic studies of the Marmara Region. These studies have brought new light on the structure (Armijo et al., 1999; Le Pichon et al., 2001; Le Pichon et al., 2003) and active deformation (Hergert and Heidbach, 2010; Ergintav et al., 2014; Bohnhoff et al., 2017) of the NAFS and raised important questions such as the locking condition of the fault and the possibility of a large rupture in the Sea of Marmara.

Figure 1 Regional tectonic setting. Modified from Duman et al. (2018). The active faults are in black except for the North Anatolian Fault Zone (NAFZ) which is in red. The black and red arrows show the plate motion relative to a fixed Eurasia. The dotted box indicates the study area detailed in Fig. 2.

Figure 2 Fault system and the earthquake catalogue used in this study (Emre et al., 2018). (a) The fault sections are delimited by the purple lines, and the names of the fault sections are indicated in black. The purple triangles indicate the locations of palaeoseismic studies, and the numbers refer to the references of the specified fault section listed in Table 3. (b) Earthquakes of magnitude greater than magnitude 4.5 and during the complete period (see Table 2) in the Marmara Region. The years of occurrence of the largest earthquakes (M>7) are indicated.

Building on the findings of these studies, several earthquake rate forecast (ERF) models and hazard maps have been developed (Erdik et al., 2004; Gülerce et al., 2017; Murru et al., 2016; Sesetyan et al., 2018; Demircioğlu et al., 2018). While these studies improve our understanding of the hazard for Istanbul, none of them fully explore the epistemic uncertainties in the parametrization of the ERF parameters that have the potential to affect the seismic hazard in Istanbul. The NAFS being offshore in its portion closer to Istanbul, satellite-based geodesy (GPS and INSAR) is not able to resolve accurately the locking condition of the fault which remains debated (Klein et al., 2017). The possibility of creep in the Sea of Marmara might have an impact on the seismic hazard for Istanbul and should be explored in the seismic hazard assessment.

Erdik et al. (2004) and Demircioğlu et al. (2018) have used the classic approach to model faults and background seismicity in seismic hazard. The rate of earthquakes in the background is based on the local seismic catalogue up to a threshold magnitude (M5.5), and the rate of larger magnitude earthquakes is based on the slip rate of the faults. This approach can lead to several limitations in probabilistic seismic hazard assessment (PSHA) and eventually in the seismic risk assessment. First, this approach lacks the consideration that earthquakes with a magnitude larger than the threshold magnitude can occur in the background as has been observed in July 2019 in California where a magnitude 7.1 earthquake occurred on a fault that was not modelled as active in the UCERF3 model (Brandenberg et al., 2019; Field et al., 2014), but the possibility of such an earthquake occurring in the background was considered. Second, this approach when used with a spatially uniform distribution of the seismicity cannot consider that intermediate magnitude earthquakes (lower than the threshold magnitude) can be more likely to occur on the fault than on any given point of the background.

In this study, we use the recently developed the SHERIFS (Seismic Hazard and Earthquake Rate In Fault Systems) code (Chartier et al., 2019) that allows us to explore these uncertainties using a logic tree approach. Furthermore, SHERIFS considers a system level approach, in which different rupture scenarios are explored in an aleatory manner, relaxing fault segmentation and defining complex multi-fault ruptures. Through the SHERIFS approach we can also address two issues: (1) the uncertainty in the size of the largest rupture that may occur along the NAFS and (2) the shape of the magnitude frequency distribution (MFD).

The geometry of the network could potentially host larger ruptures than the one observed in historical time. For example, the change in azimuth in the geometry of the fault in front of Princes' Island (Fig. 2) that was not crossed by the 1999 Izmit earthquake could be crossed by a future earthquake as has been shown by dynamic rupture modelling (Oglesby et al., 2008; Aochi and Ulrich, 2015). Furthermore, ruptures that run through more complex fault geometries have been observed in fault systems (e.g. the 2016 Kaikoura earthquake; Klinger et al., 2018).

It has long been observed that the number of earthquakes in a region decreases with the magnitude of the earthquake. Gutenberg and Richter (1944) established the logarithmic decrease in the number of earthquakes with the magnitude (noted as GR law hereafter). However, the possibility for the seismicity on individual faults to follow a different type of MFD has been discussed, notably in California where several studies have discussed this aspect of earthquake statistics along the San Andreas strike slip fault systems, arguing either for the GR law (Page and Felzer, 2015) or for a discrepancy (Schwartz and Coppersmith, 1984). Recently, Stirling and Gerstenberger (2018) have analysed several fault zones in New Zealand and have argued for the systematic exploration of the uncertainty in the shape of the MFD when modelling faults in seismic hazard assessment. We explore this uncertainty in this study.

The earthquake rates modelled with SHERIFS using each combination of uncertainties will be compared to the earthquake rates calculated from the earthquake catalogue and the palaeoseismic records in order to give a score to each model. However, for some hypotheses, the comparison with the data is not sufficient for rating one hypothesis against another. We tackle this issue by modelling a synthetic catalogue of the fault system using the earthquake simulator RSQSim (Richards-Dinger and Dieterich, 2012). By analysing the statistics of the synthetic catalogue, we are able to discuss the physical validity of the different hypotheses explored and give more weight to the branches of the logic tree that are more compatible with physics-based models.

The final goal of this study is to enhance the understanding of the seismic risk in Istanbul. The earthquake occurrence models developed in this study will be used in a later study to calculate to probability of collapse of a theoretical building in Istanbul. The impact of each uncertainty affecting the earthquake rates on the uncertainty in the estimates of the probability of collapse will be quantified and discussed.

2 Model parameters for the western North Anatolian Fault system

Fault traces and some of the slip-rate estimates along the North Anatolian Fault system have been recently updated in Emre et al. (2018) and references therein. In this study, we rely on their map (available online at http://www.mta.gov.tr/v3.0/hizmetler/yenilenmis-diri-fay-haritalari, last access: 20 August 2019​​​​​​​) to digitize the fault traces (Fig. 2).

The parameters of the main faults used in this study are presented in Table 1. Slip rates are taken from Hergert and Heidbach (2010) who inverted the GPS velocity field in a geomechanical model in order to calculate the slip rate of the fault network. These slip rates are in general agreement with the slip rates calculated by Flerit et al. (2004) and Ergintav et al. (2014). Since they agree with the long-term geological slip-rates of the faults, these values of slip rates have been preferred to those estimated by Reilinger et al. (2006) calculated at a larger scale because they tend to overestimate the slip rates compared to the geological estimates (Hergert and Heidbach, 2010). Uncertainties in slip rates (Table 1) are taken into account by a random exploration in a uniform distribution between the minimal slip-rate value and the maximal slip-rate value set from the literature (Hergert and Heidbach, 2010; Klein et al., 2017) and weighted in order to give the most weight to the mean value. It can be noted that the range of slip rates explored in this study excludes values that vastly differ from the geological estimate (such as the 23 mm/yr according to, for example, Le Pichon et al., 2003) since they integrate the deformation on a wider area than the fault itself.

Table 1 Model parameters of the faults of the Marmara Region, closer to Istanbul. The full parameter table for all the faults in the model is available in the electronic Supplement. See the text for details on slip-rate setting and definition of partial and deep creep.

The scientific community has been debating over the possibility of the NAF creeping in the western Marmara Region, along the Terkidag, central basin, Kumburgaz and Avcilar sections of the fault (Table 1, Fig. 2). Creep has been observed on other strike-slip faults in the world like for San Andreas Fault (Nason, 1973) and the Longitudinal Valley Fault in Taiwan (Hsu and Bürgmann, 2006) and is expected to release an important part of the tectonic load of the fault without producing major earthquakes, hence possibly reducing the seismic hazard in the region. The large scientific efforts of the past two decades have led us to conclude that the eastern part of the fault in the Sea of Marmara is locked (Diao et al., 2016). However, the distance from the fault to the shore in the western part of the sea remains a challenge for traditional land-based or satellite-based geodetic instruments. While novel sea-floor geodesy experiments seem to suggest that the fault is locked at the surface at least at some places (Lange et al., 2019), these conclusions are still under discussion (Yamamoto et al., 2019). Relocated microseismicity and the identification of repeater earthquakes suggest that the lower part of the fault might be creeping (Schmittbuhl et al., 2016a, b). Based on the latest studies of the fault in this area, we propose exploring four models of locking condition for the NAF that represent the current state of knowledge. The “creep” hypothesis considers the fault as fully creeping, thus reducing its contribution to the seismic slip-rate budget to 0. The “partial creep” hypothesis considers the fault releasing half of its slip rate as creep and half as seismic moment, reducing its slip rate by half. The “deep creep” hypothesis considers that the fault is fully locked for the first 5 km and creeping below, and the “fully locked” hypothesis considers that 100 % of the fault slip rate can be released as earthquakes.

For the faults within the Armutlu peninsula (Fig. 2), we were not able to find slip-rate estimates. However, the velocity field shows very little deformation (Ergintav et al., 2014) within the Armutlu peninsula. Based on this observation, we assume the slip rate of these faults to be less than 1 mm/yr. It is worth noting that these faults are relatively far from Istanbul in comparison to the NAF, and their participation in the seismic hazard affecting the city is expected to be negligible.

In this study, we combined two catalogues: the earthquake catalogue from Kadirioğlu et al. (2018) homogenized in M_w for the period 1900–2012 and the catalogue SHEEC (Stucchi et al., 2013) for the period 1000–1899. The completeness period used in this study is based on these two studies and is presented in Table 2.

Woessner et al., 2015Kadirioğlu et al. (2018)

Table 2 Completeness time as a function of magnitudes used in this study.

Based on the analysis of the depth distribution of earthquakes in the instrumental catalogue (Kadirioğlu et al., 2018), 80 % of the earthquakes are in the first 15 km of the crust in the region of interest (Fig. 3). Furthermore, slip inversion of the Izmit earthquake show that most of the slip during the earthquake was contained in the first 15 km of the crust (Reilinger et al., 2000). Thus in this study we have limited the seismogenic thickness of the fault system to a depth of 15 km with the exception of the Adalar section of the fault stopping at 10 km deep in order to avoid crossing the Cinarcik fault at depth.

Figure 3 (a) Magnitude frequency distribution for the study area calculated from the catalogue and the completeness periods presented in Table 2. The red dashes are individual Monte Carlo samples of the earthquake magnitude uncertainties and the completeness period. The mean rate value is indicated by grey squares, and the 16th and 84th percentiles are indicated by the grey area. (b) Depth distribution of earthquakes of magnitude 4.5 and above in the background zone (Kadirioğlu et al., 2018) since 1992.

Dikbaş et al. (2018) present a review of the palaeoseismic studies that have been carried out on the Izmit, Adapazari, Karadere and Duzce segments of the fault. They propose an interpretation of the rupture history of these segments. Since we aim to compare the rate of observed palaeo-earthquakes to the modelled rates that are considered Poissonian in PSHA, we decided to consider the number of observed events during the observation time rather than the inter-event time for the earthquake rate calculation. Based on the number of ruptures observed in one location and the length of the observation period, we estimate the annual rate of earthquakes with a magnitude larger than or equal to 7.2±0.1 for different sections of these segments. The locations of these segments are presented in Fig. 2, and the estimated rates are presented in Table 3. At the west of the Sea of Marmara, for the Ganos segment, the palaeo-earthquake rates have been calculated from Rockwell et al. (2009) and Meghraoui et al. (2012). Both studies lead to similar rates for earthquakes larger than 7.2 (Table 3). The palaeo-earthquake rates are calculated by dividing the number of events by the observation time, and the uncertainty reflects the statistical uncertainty due to the small number of observation assuming a Poissonian process.

Rockwell et al. (2009)Meghraoui et al. (2012)Klinger et al. (2003)Dikbaş et al. (2018)Dikbaş et al. (2018)Dikbaş et al. (2018)Dikbaş and Akyüz (2014)Pantosti et al. (2008)

Table 3 Annual rate of M 7.2+ earthquakes on sections of the NAF deduced from palaeoseismic studies.

3 Earthquake rate modelling

In this section, we describe two approaches for calculating earthquake rates in fault systems. The first approach is SHERIFS (Chartier et al., 2019), a statistical approach that converts the slip rate of the faults into earthquake rates built to follow a given shape of MFD at the system level. The second approach is RSQSim that generates a long catalogue of synthetic ruptures using the loading rates and the rate and state equation (Richards-Dinger and Dieterich, 2012). SHERIFS allows us to explore a wide range of epistemic uncertainties in the input hypotheses and RSQSim allows us to discuss which of these hypotheses are physically plausible.

3.1 The statistical approach: SHERIFS

3.1.1 Core principle and main input hypothesis

SHERIFS uses an iterative budget spending approach of the slip rate of the fault to calculate the annual rate of occurrence of each rupture of a predefined set of ruptures. In an iterative manner, SHERIFS randomly selects user-defined rupture scenarios for which faults involved have a slip-rate budget to spend. The random selection is done in order to ensure that the resulting system level MFD has the shape imposed as input (b value in the case of a Gutenberg–Richter MFD, for example). It is important to recall that SHERIFS does not simulate earthquakes but only converts slip rates into earthquake rates. For this reason, SHERIFS is computationally light which allows aleatory (combination of ruptures) and epistemic (logic tree) uncertainties concerning the fault system to be easily explored.

SHERIFS takes as input the geometry and slip rate of the faults, the set of multi-fault ruptures that can be expected in the fault network, and the shape of the MFD defined at the fault system level. Before the calculation, the actual value of the MFD and the shape of the MFD of each individual fault are not known. They will be deduced from the fault slip-rate budget and the other hypotheses. Depending on the combination of input hypotheses and fault parameters, SHERIFS can consider part of the slip-rate budget of some faults as non-main-shock (NMS) slip in order to respect the target MFD shape. A NMS of more than 30 % is most likely an indication that the combination of input hypotheses used does not agree with the fault parameters in the SHERIFS framework and that they should be reconsidered.

3.1.2 Background seismicity

One uncertainty that SHERIFS allows us to explore is the proportion of seismicity that can occur in the background on faults that are unknown or not considered as active in the model. In most PSHAs, this is taken into account by a background zone with a GR MFD truncated at a given magnitude (Woessner et al., 2015).

In SHERIFS, it is possible to define a priori the proportion of earthquakes that can be expected on the faults and the proportion in the background for each range of magnitude. In order to assess these proportions, we analyse the spatial distribution of earthquakes compared to the fault traces (Fig. 2). Considering the poor knowledge on the epicentral location of the historical earthquakes, we only consider the instrumental catalogue after 1970. Since we are interested in the hazard in Istanbul, we only consider the spatial distribution of earthquakes around the Sea of Marmara (the central zone in Fig. 2).

Most of the observed seismicity of the Marmara occurs close to the known faults; therefore, we chose to consider that a large proportion of the seismicity is on the known faults for a wide range of magnitudes (Table 4, Fig. 4). However, in order to represent the epistemic uncertainty associated with the proportion of seismicity considered for the faults versus that for the background, we set up three branches of the logic tree corresponding to three different background hypotheses (Fig. 4, Table 4). The proportion of earthquakes considered to occur on the faults for each branch is presented in Fig. 4. Studies of the NAFS off-fault deformation might shed some light on these value in the future (Şengör and Zabcı, 2019).

Table 4 Ratio of earthquakes assumed to be on the faults over the total number of earthquakes in the system for each background hypothesis.

Figure 4 Distribution of earthquakes between the background seismicity and the faults for each background hypothesis. In black is the earthquake rate for the entire fault system (faults + background), the solid colour line is the faults only, and the dashed line is the background only.

3.1.3 MFD

The annual rates of earthquakes obtained by analysis of the Kadirioğlu et al. (2018) catalogue in the central zone (Fig. 2) indicate an MFD diverging from a GR MFD for magnitudes larger than 6.5 (Fig. 3). The rates are obtained while exploring the uncertainties in the magnitude of earthquakes, in the completeness period length using a Monte Carlo approach and in those linked to the low number of earthquakes for the larger magnitudes. Despite this large number of uncertainties explored, the resulting MFD significantly differs from the GR MFD shape. We can observe a larger number of earthquakes of magnitude 7.0 and above than predicted by the GR law. Like any earthquake catalogue, the Turkish catalogue is short in comparison with the seismic cycle. Only six earthquakes larger than magnitude 7.0 are in the catalogue, and it can be argued that the observed MFD is a result of the incomplete sample of the phenomenon and that the long-term MFD should follow a GR distribution. This debate on the catalogue MFD concerns the Californian earthquake catalogue as well (Page and Felzer, 2015; Parsons et al., 2018).

We explored two alternative hypotheses for the target MFD: one in which the target MFD follows a GR truncated between a minimum magnitude and a maximum magnitude, and in which the target MFD follows a shape tuned to that of the rates deduced from the catalogue. The tuned shape (TS) MFD is described by the following equation and is composed of two parts, both defined by a double truncated GR with a b value deduced from the lower-magnitude part of the catalogue ($\mathrm{4.5}<M<\mathrm{6.0}$) that follows an exponential decrease. This formulation of the MFD has been developed especially to closely resemble the MFD observed in the earthquake catalogue of our study area (Fig. 3). If other formulations could also have been explored, this set of equations was found appropriate for the expected use in this study: matching the rates observed in the earthquake catalogue while exploring the uncertainty in the set of possible rupture scenarios. While we do not argue for the universal nature of this formulation, variations around this formulation could be useful for other study regions.

$$P_i\left(M\right)=\left\{\begin{array}{ll}\frac{e^{-\mathit{\beta }\times (M-\mathrm{0.05}-M_{\min })}}{\mathrm{1}-e^{-\mathit{\beta }\times (\mathrm{6.8}-M_{\min })}}-\frac{e^{-\mathit{\beta }\times (M+\mathrm{0.05}-M_{\min })}}{\mathrm{1}-e^{-\mathit{\beta }\times (\mathrm{6.8}-M_{\min })}},& \text{if}M\ge \mathrm{6.6}\\ P_i\left(\mathrm{6.6}\right)/\mathrm{3},& \text{if}M=\mathrm{6.7}\\ P_i\left(\mathrm{6.3}\right),& \text{if}M=\mathrm{6.8}\\ P_i\left(\mathrm{6.8}\right)\times \left({\mathrm{10}}^{(-b\times (M-\mathrm{6.8}\left)\right)}\right),& \text{if}M>\mathrm{6.8}\end{array},\right.$$

where P_i is the density function depending on M, b describes the linear decrease with M, $\mathit{\beta }=b\times \ln \mathrm{10}$, and M_min is the minimal magnitude considered in the calculation.

3.1.4 Set of rupture scenario

The historical catalogue does not contain any earthquake with a magnitude larger than 7.5 since 1700. Based on the statistics of the earthquake catalogue, Bohnhoff et al. (2016) consider this value of 7.5 to be the largest possible magnitude in the Marmara Region, while earthquakes can reach up to magnitude 8.0 in the eastern part of the NAFS. A recent review of palaeoseismological studies and historical earthquake studies (Dikbaş et al., 2018) suggests that the largest earthquake occurring in the region could have been a rupture including both the 1999 Izmit earthquake rupture area and the 1999 Duzce earthquake rupture area resulting in a magnitude 7.7 earthquake. We define a set (Set 1) of more than 300 fault to fault (FtF) ruptures allowing a large diversity of ruptures to be possible in the system with the largest magnitude of 7.7 (Fig. 5a).

Figure 5 Example of the six largest ruptures included in the models using Set 1 (a) and the three largest ruptures included in the models using Set 2 (b). Ruptures are stacked on top of each other for readability. The full list of ruptures is available in the electronic Supplement.

The resolution of the earthquake records diminishes as we consider older events, and it is possible that larger earthquakes occurred in the NAFS but may not have been observed either in historical times or in the palaeoseismic records. Furthermore, since most bends in the system could be crossed by a rupture without jumping a large distance or by large changes in azimuth, we need to imagine that larger earthquakes might be possible. This hypothesis has been considered in previous studies (Murru et al., 2016; Mignan et al., 2015). We thus build a set (Set 2) of FtF ruptures to explore this possibility in which the largest possible ruptures correspond to a magnitude 8.0 earthquake (Fig. 5b).

Since there are several hundred ruptures considered in each hypothesis, we chose to only illustrate the larger ones for each set in Fig. 5. Although not illustrated in Fig. 5, in each set, all combinations of fault sections smaller than these ruptures are considered even if not illustrated in the figure.

3.1.5 Locking condition of the NAFS in the western Marmara Region

As exposed during the presentation of the NAFS earlier, there is uncertainty concerning the locking condition of the NAF in the western Sea of Marmara. Four hypotheses of locking conditions ranging from fully creeping to fully locked are explored in a logic tree to represent the current state of knowledge. The values of slip rates used for each hypothesis are presented in Table 1.

3.1.6 Results: earthquake rates modelled by SHERIFS

SHERIFS is run with the hypotheses of each branch of the logic tree (Fig. 6). We then compare the modelled annual earthquake rates with the seismicity rates from the earthquake catalogue and the palaeoseismic records.

Figure 6 Logic tree explored in this study. For each branch, the scaling law parameters and the slip-rate uncertainties are explored through 10 random samples.

In this study, we only modelled part of the NAF (inside the dotted box in Fig. 1); therefore, some boundary effects might occur on the faults that would rupture with faults outside of the modelled zone. With the aim of modelling the seismic risk in Istanbul, far from the system boundaries, we compare only the rates for a region smaller than the whole modelled system centred on Istanbul (Fig. 2) that is not affected by boundary effects. In this zone, we extract the earthquakes from the catalogue and the modelled rate of ruptures. For ruptures that are only partly located in the zone, the rate is corrected by taking into account the proportion of the fault in the zone.

The rates modelled using SHERIFS are slightly higher than the ones of the catalogue in the central zone (Fig. 7). However, some combinations of hypotheses better fit the rate from the catalogue than others. The branches using the GR target shape lead to a rate of earthquakes of magnitude between 4.5 and 6.5 that is larger than the observed rate for all models (Fig. 7). The tuned shape (TS) hypothesis is in better agreement with the catalogue for this range of magnitudes. Both MFD hypotheses reproduce the rate of large earthquakes equally well. The branches using Set 2 of rupture scenarios, allowing ruptures up to magnitude 8, have lower earthquake rates than the branches using Set 1, and they obtain a better fit with the rates from the catalogue. The four different creep hypotheses lead to similar earthquake rates at the level of the central zone (Fig. 7). In conclusion, the branches using the combination of the TS MFD and Set 2 of ruptures best fit the rates calculated from the earthquake catalogue.

Figure 7 Comparison between the rates calculated from the earthquake catalogue (in red with uncertainties in grey) and the model rates (in green) for the central zone of the fault system, close to Istanbul (Fig. 2). The dashed green lines are the MFD for each individual model of the logic tree. The dotted green lines are the 16th and 84th percentiles of the distribution, and the continuous green line is the mean value of the distribution. Each figure shows the match between model and data for each branch of the logic tree, organized as a table. A fit of 1 is a perfect fit, and a value close to 0 expresses a very poor fit.

In Fig. 8, we compare the participation rate of fault sections where an estimate of the rate of earthquakes is available based on palaeoseismic studies. The participation rate is the rate of all the ruptures in which the given fault section is involved. For most sites, the rate of earthquakes modelled is in general agreement with the rate deduced from the palaeoseismicity for all branches of the logic tree.

Figure 8 Comparison of the modelled rupture rates with the rates calculated from the palaeo-earthquake record at each palaeo-earthquake site (Fig. 2). The green curves are the modelled rates of ruptures rupturing the given section. The dashed lines are each individual model of the logic tree. The darker dashed lines are the 16th and 84th percentiles of the distribution, and the continuous line is the mean value of the distribution. Palaeo-earthquake rates are in purple (see Table 3).

At the Ganos 1 site (Fig. 9), the branch “creep”, assuming a complete creep in the west of the Sea of Marmara, leads to rupture rates significantly lower than the palaeoseismic rates. For this branch of the logic tree, the Ganos fault spends around 46 % of its slip-rate budget as non-main shock (NMS) slip. The proportion of NMS slip is the proportion of slip-rate budget that could not be spent in seismic moment rates in SHERIFS. A large NMS (>20 %–30 %) indicates the inability in SHERIFS to satisfy the hypotheses of MFD, slip-rate budget and the set of ruptures. The introduction of a fault completely creeping in the western part of the Sea of Marmara limits the number of large earthquakes along the Ganos fault. In the SHERIFS framework, the ability of the Ganos fault to release seismic moment jointly with the faults of the Sea of Marmara is required in order for it to spend its full budget.

Figure 9 Comparison of the modelled rupture rates with the rates calculated from the palaeo-earthquake record at the Ganos 1 site (Fig. 2) for different branches of the logic tree. The green curves are the modelled rate of ruptures rupturing the Ganos segment. The dashed lines are each individual model of the logic tree. The dotted lines are the 16th and 84th percentiles of the distribution, and the continuous line is the mean value of the distribution. Palaeo-earthquake rates are in purple (Table 3).

3.2 Insight from the physics-based approach: RSQSim

The comparison between the modelled earthquake rates and the data allows us to reduce part of the uncertainties explored in the logic tree, but some uncertainties still remain, notably the uncertainties concerning the MFD shape and the maximum rupture size.

In the hope of reducing these uncertainties, we implemented our fault system in the physics-based simulator RSQSim (Richards-Dinger and Dieterich, 2012) and analysed the resulting synthetic earthquake catalogues. In this study, we do not aim to reproduce perfectly the earthquake process using RSQSim but rather to bring additional information to the seismogenic source model generated using SHERIFS.

3.2.1 Core principle and main hypotheses

RSQSim is a boundary element model that applies the rate and state equation (Dieterich, 1978) in which each element of a fault system can be in three different states: loading, nucleating the earthquake rupture or sliding. Element interactions are considered to be quasi-static, therefore neglecting the dynamic influence of seismic waves that would be generated by sliding elements. Where the elements are sliding, they are doing so at a constant slip rate. These simplifications of the earthquake physics and the rate and state law allow RSQSim to be computationally efficient and to generate long earthquake histories for a large fault network.

RSQSim takes as input the same information as SHERIFS concerning the fault parameters but does not require a target MFD shape or a set of possible FtF ruptures. The MFD of the fault system and the ruptures will be deduced from the synthetic earthquake catalogue which results from the combination of the fault loading rates and the rate and state friction law as implemented in RSQSim.

In this study, we will use the same friction parameters a and b for all the elements of the fault system (0.01 and 0.015 respectively). These a and b value are based on empirical measurements from laboratory earthquake experiments and have been widely used for modelling earthquake cycles with the rate and state equation (Marone, 1998). The stress loading rate of the faults in the system is defined using back-slip loading. Back-slip loading calculates the stress rate for each element of the fault that corresponds to the long-term slip rate of the fault.

We are using triangular elements of 1 km size. For this reason, we will only discuss earthquakes larger than magnitude 6 that rupture a number of elements large enough (around 100) to be representative of an earthquake rupture.

3.2.2 Earthquake rates in RSQSim

Using the fault geometry and slip rate with the fully locked hypothesis (Table 1, Fig. 2) as input for RSQSim, we simulate a 10 000-year-long earthquake catalogue from which we remove the first 1000 years of seismicity to let the fault system run through several earthquake cycles before starting the statistical analysis.

Figure 10 (a) MFD of the synthetic catalogues generated from the RSQSim simulations exploring a range of input parameters. (b) Normalized MFD of the earthquake catalogue, the MFD of the faults modelled using SHERIFS with a GR MFD or a TS MFD (the background seismicity has been removed), and the RSQSim synthetic catalogue (a=0.01 b=0.015).

Richards-Dinger and Dieterich (2012) compared RSQSim with fully dynamic earthquake simulations and have shown that both modelling approaches lead to similar earthquake ruptures in terms of slip history and final slip on the fault. We observe that the scaling of earthquakes in RSQSim is in good agreement with the scaling laws (Tullis et al., 2012) in terms of geometrical scaling and of stress drop. Interestingly, RSQSim reproduces quite well the rupture extent and magnitude of the 1999 Izmit and Duzce earthquakes. Based on these comparisons we conclude that even if the physics of earthquake is simplified to allow for cost efficient calculations, the simulated catalogue is representative of the natural system at least to the first order. In order to not over-interpret the results we limit our use of the RSQSim simulated catalogue to the discussion of two first order questions: (1) does the long-term MFD of the fault system follow a GR law and (2) are magnitude 8.0 earthquakes possible in this fault system?

Considering the MFD of the synthetic catalogue, we can reach two conclusions: RSQSim cannot reproduce a GR MFD with the given fault system, and the maximum magnitude that can be generated is closer to magnitude 7.7 than to magnitude 8.0. These two conclusions are not affected by uncertainties in a and b values (Fig. 10).

We also explored the impact of having a region of the fault creeping in the western Sea of Marmara. In a way similar to the deep creep model used in SHERIFS, in RSQSim, elements that are below 5 km depth were attributed a b value much lower than the a value. The resulting MFD (Fig. 10) leads to the same first order conclusions as for the fully locked hypothesis. This is to be expected. Indeed, the region where creep is suspected is of a negligible size compared to the entire fault system and does not drastically affect the shape at the system level MFD. The creeping condition of the fault in the western Sea of Marmara most likely affects the participation rates of neighbouring faults; however, here we only consider the first order results of the calculation.

The results are also compared with the hybrid loading method used in Shaw et al. (2018) with the same input parameters as in the Californian application. The hybrid loading method allows us to remove the edge effects on the stress field of the fault that occurs with the standard back-slip loading (Shaw et al., 2018). Different smoothness parameters are tested. The synthetic catalogues do not follow GR MFD, strengthening our previous observation.

4 Discussion: the weight of the epistemic uncertainties in the logic tree

The flexibility of SHERIFS for modelling the earthquake rates in the NAFS allowed us to explore a wide range of uncertainties in a logic tree framework. The SHERIFS approach does not contain any physical constraints; it can therefore accommodate the different input hypotheses that are being discussed by the scientific community but also can lead to a large range of uncertainties in the earthquake rates. Since SHERIFS only uses the fault data as input, it is possible to compare the modelled rates with rates calculated from the earthquake catalogue and palaeo-earthquakes which can be considered as independent data from the SHERIFS input.

Branch weights in a logic tree are usually based on the scientific value of each hypothesis explored in the logic tree but not on the capacity of the modelled earthquake rates in each branch to reproduce the data. The weight of an individual branch of the logic tree is simply the product of the weights of each hypothesis used for the branch. In this discussion, we propose a novel approach to set the weight of each branch of the logic tree accounting for both the input hypotheses used and the ability to reproduce the independent data.

We set up a quantitative scoring system in order to set the weight of each individual branch of the logic tree. For each branch, four scores are calculated (Fig. 11). One score judges the ability of the combination of input parameters and hypotheses to spend the slip-rate budget of the faults as earthquake rate and not generate a large amount of NMS (S1). Two scores compare the model to the data, judging the capacity of the modelled earthquake rates to reproduce the earthquake rates deduced from the catalogue (S2) and to reproduce the earthquake rates deduced from the palaeoseismic records (S3). The last score rates a priory the input hypotheses based on the analysis of the RSQSim synthetic catalogue (S4). The final weight of each branch of the logic tree is given by the combination of the four scores (Fig. 11). The code allowing the computation of the scores and the excel file showing the computation of the weights are provided in the electronic Supplement.

Figure 11 Scoring system allowing the individual weighting of each branch of the logic tree according to both its capacity to spend the fault slip-rate budget in earthquake rates (S1), its capacity to reproduce the earthquake rates observed in the data, and the agreement between the input hypotheses (S2 and S3) and the results of the discussion on the physics-based synthetic catalogue (S4).

4.1 Uncertainty in the creep

4.1.1 Scoring based on the ratio of NMS

We calculate the NMS value for each model and each fault section of the central zone. Since a large NMS value is likely linked to incompatibilities between input hypotheses, given the SHERIFS framework, a low score is attributed to models with high NMS slip value. For a given model, if the mean NMS value for the fault sections of the central zone is greater than 40 % or if the NMS value of one of the faults of the central region is greater than 50 %, the score is 0. The score is 1 if the average NMS of the sections of the central region is less than 20 %. Between the average values of 20 % and 40 %, the score linearly decreases from 1 to 0 as the NMS value increases. As a result, the average score of the fully creeping models is 0.27, the average score of the partly creeping models is 0.6, and the average scores of both the deep creep and fully locked models are 0.96.

The modelling of the creep conditions assumed on the Terkidag section play a predominant role in the modelling of earthquake rates on the adjacent Ganos section and hence the scoring based on NMS. The models considering the Terkidag fault as completely creeping have difficulties reproducing the palaeo-earthquake rates estimated on the Ganos fault (Fig. 9) and lead to a high level of NMS on this fault. In the fully creeping model, the Ganos section is only able to convert 54 % of its slip rate into seismic moment rate, and 46 % is considered NMS slip. In the partly creeping model, the Ganos section spends 71 % of its slip rate as moment rate. In the deep creep model and the fully locked model, the Ganos section spends 90 % of its slip rate as seismic moment rate.

4.1.2 Scoring based on the fit to the rate of the palaeo-earthquake

Figure 9 exposes the underestimation of the rate of earthquakes on the Ganos section compared to data when considering the creep branch of the logic tree. The comparison with the palaeoseismicity is done at all sites where palaeoseismicity is available for each model of the logic tree (see above). At each site, the rate of palaeo-earthquake is extracted from the literature (Table 3) and represented as a 2D Gaussian distribution in order to capture the uncertainty both in terms of magnitude and annual rate of occurrence. At each site, the modelled rates are compared to this distribution. The closer the participation rate for a given section is from the maximum of the palaeo-earthquake rate distribution, the larger the score. The score of the branch is the average score for each palaeo-site. The average scores on the palaeo-earthquake fit for the creep, partial creep, deep creep and fully locked models are 0.65, 0.73, 0.82 and 0.81 respectively.

4.1.3 Scoring based on the fit to the earthquake catalogue

The average scores for models imposing a TS MFD are only slightly higher than those with a GR MFD as a target shape for the MFD of the fault system (0.52 versus 0.49). While both shapes show a good agreement with the rate of large earthquakes calculated from both the catalogue and the palaeoseismicity, the models with a GR MFD overestimate the rate of small to intermediate earthquakes, as seen in Fig. 7. Given the uncertainty in the annual rates calculated from the earthquake catalogue, the general better fit obtained when using the TS MFD is only leading to a slightly stronger score for this branch compared to the GR MFD.

4.2 Scoring based on RSQSim results – MFD

It can be argued that the deviation of the MFD shape from the GR shape is an artefact of the observation period, which is too short for accurately capturing the full earthquake cycle. Such arguments have been brought forward in California in support of a GR MFD target shape even when the apparent shape of the earthquake catalogue close to the faults might differ from the GR shape (Page and Felzer, 2015). Stirling and Gerstenberger (2018) have shown that for four fault systems in New Zealand, while the data do not show a good agreement with the GR MFD hypothesis, ETAS simulations suggest that the observation time in the data is too short to accurately calculate the long-term rate of a clustered earthquake catalogue.

However, the 10 000-year-long synthetic catalogue generated by RSQSim, complete and representative of the seismic cycle by definition, also diverges from the GR MFD shape (Fig. 10). In the RSQSim catalogue, there are almost no earthquakes of magnitudes between 6.5 and 7.0. We explored a number of combinations of a and b values as input for the RSQSim equation with b–a ranging from 0.003 to 0.007 and obtained minor changes in the resulting MFD, its shape remaining different from a GR MFD. A 20 % variation around the normal stress parameter set at 100 MPa had an impact on the number of simulated earthquakes but not on the shape of the MFD. It is therefore likely that the shape of the MFD is controlled by the geometry of this particular fault system.

According to the scaling laws (Thingbaijam et al., 2017), 6.5 is the magnitude of an earthquake rupturing the whole fault width. We speculate that in RSQSim, an earthquake rupturing the entire width of the fault is more likely to rupture neighbouring faults unless stopped by geometrical complexities or the presence of faults that are still in the early part of their cycle. The average length of the straight sections of the faults is around 100 km (Fig. 1), which corresponds to an earthquake of magnitude around 7.2. Ruptures are able to overcome bends and asperities if the neighbouring fault is close to rupture leading to larger magnitude earthquakes. Since they require several sections to be close to rupture simultaneously, these ruptures are less frequent. We obtain therefore a bi-modal distribution with the earthquakes up to magnitude 6.3 following a GR, rupturing different portions of the fault plane, very few earthquakes between magnitude 6.3 and 7.0, and a population of earthquakes larger than magnitude 7.0 with the number of earthquakes decreasing with the magnitude.

Based on the comparison between the modelled rates with the rates of the earthquake catalogue and the results of RSQSim, and in consideration of the scientific debate around the MFD, we suggest a stronger score (s_mfd in Fig. 11) on the tuned shape MFD branch. The branch exploring the GR MFD should remain part of the logic tree, however. Consequently, we set the score of the branch tuned shape MFD to 0.8 and the score of the branch GR MFD to 0.2.

4.3 Scoring based on RSQSIM results – Mmax

Based on the comparison between the modelled rates from SHERIFS and the rate calculated from the earthquake catalogue and the palaeoseismicity studies, it is not possible to weight differently the two branches of the logic tree exploring the uncertainty in the rupture scenarios. While the fit between the modelled rates and the catalogue rates using Set 2 of ruptures, allowing larger ruptures, leads to a slightly better fit than the models using Set 1 of ruptures, both fits can be considered as satisfying (Fig. 7). Whether ruptures up to magnitude 8.0 are considered or not, the rates in the models are matching the rates from the data when using the TS MFD as a target.

Figure 12 Weighted logic tree established in this study. For each branch, the scaling law parameters and the slip-rate uncertainties are explored through 10 random samples. The weight of each branch is indicated by the bold number.

In the several 10 000-year-long catalogues simulated by RSQSim, we do not observe earthquakes larger than 7.7 (Fig. 10). We explored various smoothing factors of the faults and input parameters in RSQSim leading to small variations in the maximum magnitude, but we never obtained magnitudes as large as 8.0 in the synthetic catalogue (Fig. 10). In the SHERIFS models, earthquakes with magnitude larger than or equal to 8.0 have an annual frequency of 5 × 10⁻⁴. The likelihood of observing at least one earthquake in a 10 000-year-long catalogue is 99.3 %. While the bends and other geometrical complexities of the fault trace are crossed by ruptures in the model, no rupture crosses enough complexities to generate such large earthquakes. We can assume that the fault system is very chaotic and that having a large number of neighbouring fault sections close to rupture is very unlikely. However, it is important to recall that RSQSim does not fully model the dynamic effects of the earthquake rupture. It is possible that this simplification in the modelling leads to underestimating the likelihood of a rupture propagating through fault complexities. In consideration of these arguments, we set the s_set score of the branch “Set 1” to 0.7 and the s_set score of the branch “Set 2” to 0.3.

4.4 Discussion on the use of RSQSim

In this study, we used RSQSim as an auxiliary tool to bring additional information for the weighting of the logic tree. As SHERIFS does not model the physics of earthquake but relies on a statistic approach, a physics-based model is complementary. We recognize that further development could be made to the RSQSim model and the exploration of the uncertainties in the input parameters done in this study is not exhaustive. However, we believe that the results exposed above bring sufficient information to affect the balance of the weights of the logic tree, in particular in the cases when the results corroborate the observations made in the data. On another hand, we do not believe these physics-based calculations are sufficiently exhaustive to bring the weight of a logic tree to zero and remove completely a hypothesis. Hypotheses could be removed from the logic tree only if additional modelling with alternative physics-based approaches are performed and lead to similar conclusions.

4.5 The weighted logic tree

In the discussion, we established four types of scores for the logic tree branches: the weights established from the analysis of the RSQSim synthetic catalogues and those established from the comparison of the earthquake rates modelled using SHERIFS with the rates calculated from the earthquake catalogue and the palaeo-earthquake record. For each individual branch of the logic tree, these scores are convolved into a final weight unique to the branch (see equation in Fig. 11).

The final weight of each hypothesis is calculated by summing all the branches using a given hypothesis (Fig. 12). Due to the high NMS value of the models using the creeping hypothesis and their underestimation of the palaeo-earthquake rates, this hypothesis has the smallest weight in the logic tree.

The weights of the MFD hypothesis branches and the set of rupture scenario branches are strongly influenced by the scores imposed after the analysis of the RSQSim synthetic catalogue. While the fit to the catalogue is better with the TS MFD hypothesis, the fit to the palaeo-earthquake rates and the ratio of NMS is similar for both branches. Therefore, the comparison with the data affects the weight only marginally. The same can be said about the set of rupture scenario branches.

The weights of the background hypothesis branches are only affected by the scores depending on the comparison with the data. Overall, the background 1 hypothesis reproduces the earthquake rates better than the other two background hypotheses. In Fig. 7, we can observe that the modelled MFD tends to be slightly higher than the catalogue MFD for the majority of models. Therefore, background 3, adding more seismicity to the background than the two other hypotheses (Table 4), necessarily leads to earthquake rates that are in general slightly higher than the observed rates. It is worth noting that we are discussing small differences in trends that are correctly captured by the weights in the logic tree, with the three hypotheses having similar weights with a slightly stronger weight for background 1.

The impact of the weights on the annual earthquake rates is presented in Fig. 13. The weighting leads to lower earthquake rates for the small and intermediate magnitudes and larger rates for the magnitudes greater than 7.0. This is due to the reduced weight of the branches using the GR hypothesis.

Figure 13 Density function distribution of the cumulative earthquake rates in the logic tree for the central zone for different magnitudes. In black is the density function distribution when only the scoring according to the NMS is applied. In green is the density function distribution weighted according to the fit with the data (catalogue and palaeoseismicity). In purple is the final density function distribution of earthquake rates weighted according to the match with the data and the discussion based on RSQSim. The vertical bar shows the mean annual earthquake rate of each distribution. The horizontal bar extends between the 16th percentile and the 84th percentile of each distribution.

5 Conclusions

In this study, we calculated earthquake rates in the Marmara Region relying on the fault slip rate and geometry as primary information. We combined two innovative approaches: the SHERIFS approach that relies on statistical rules and the RSQSim approach that relies on physical rules.

With SHERIFS, we explored an extensive logic tree of uncertainties concerning the locking condition of the NAF in the Marmara Region, the shape of the MFD, the ratio of seismicity between the background and the faults, the largest possible rupture, and uncertainties in the slip rates and the maximum magnitude predicted by the scaling law. Rather than basing the weights of the branches of the logic tree only on expert judgement of each hypothesis, we rather take advantage of model performance by comparing results with data (earthquake catalogue and palaeo-earthquake) and weighting the branch accordingly.

In addition, the analysis of the synthetic catalogue simulated by RSQSim showed an MFD diverging from a GR and no earthquake of magnitude larger than 7.7. We explore different hypotheses on the input parameters of RSQSim and found these conclusions stable. While we do not consider that the findings coming from our models are final and that additional modelling using alternative physics-based models is necessary, we can nonetheless use these results as an additional source of information to weight the logic tree.

This allowed us to attribute a stronger weight to the branches of the logic tree showing similar features, leading to a final weighted distribution of modelled annual rates that properly represents the state of knowledge of the NAFS in the vicinity of Istanbul (Fig. 13).

The logic tree of earthquake source models developed in this study will be later used in a seismic risk assessment study to evaluate the risk of collapse of a theoretical building in Istanbul.