# Bulletin of the Seismological Society of America

## Abstract

The city of Istanbul is characterized by one of the highest levels of seismic risk in the Mediterranean region. An important source of such increased risk is the high probability of large earthquake occurrence during the coming years, which stands at about 65% likelihood owing to the existing seismic gap and the post‐1999 earthquake stress transfer at the western portion of the North Anatolian fault zone.

In this study, we simulated hybrid broadband time histories from selected earthquakes having magnitude *M*_{w}>7.0 in the Sea of Marmara within 10–20 km of Istanbul, the most probable scenarios for simulated generation of the devastating 1509 event in this region. Physics‐based rupture scenarios, which may be an indication of potential future events, are adopted to estimate the ground‐motion characteristics and its variability in the region. Two simulation techniques are used to compute a realistic time series, considering generic rock site conditions. The first is a full 3D wave propagation method used for generating low‐frequency seismograms, and the second is a stochastic finite‐fault model approach based on dynamic corner‐frequency high‐frequency seismograms. Dynamic rupture is generated and computed using a boundary integral equation method, and the propagation in the medium is realized through a finite‐difference approach. The results from the two simulation techniques are then merged by performing a weighted summation at intermediate frequencies to calculate a broadband synthetic time series.

The simulated hybrid broadband ground motions are validated by comparing peak ground acceleration, peak ground velocity (PGV), and spectral accelerations (5% damping) at different periods with the ground‐motion prediction equations in the region. Our simulations reveal strong rupture directivity and supershear rupture effects over a large spatial extent, which generate extremely high near‐fault motions exceeding the 250 cm/s PGV along the entire length of the ruptured fault.

## Introduction

Quantitative estimation of ground motion for expected earthquake scenarios is crucial, particularly in regions with high seismic potential for the near future (Parsons, 2004; Murru *et al.*, 2016). Istanbul and its surrounding areas are regions of concern for future high seismic hazards in the Euro‐Mediterranean region, the focus of which has increased following the 1999 İzmit and Düzce earthquakes. Recent studies confirm that present‐day tectonic loading, seismic slip deficit, and thus seismic hazards appear to be particularly high along segments of the Sea of Marmara, southwest of Istanbul (Armijo *et al.*, 2005; Pondard *et al.*, 2007; Bohnhoff *et al.*, 2013). In agreement with previous results of Parsons (2004), Murru *et al.* (2016) found that the probability of an earthquake (*M*_{w}≥7.0) close to Istanbul is based on a Poisson estimate of 21%–41% under the time‐dependent interaction model. Therefore, for future risk mitigation and city planning, a detailed estimate of seismic hazards in Istanbul is needed. The fundamental question for the region is that the probable rupture scenarios of the target earthquake are uncertain, even if a seismic gap is relatively well identified. This is a critical problem in ground‐motion estimation. Because near‐fault ground motion is sensitive to the faulting process, reasonable source models are required for simulation (Aochi and Ulrich, 2015).

Standard practice estimates scalar ground‐motion intensity measures using empirical ground‐motion prediction equations (GMPEs). However, modern earthquake engineering applications require the entire time series, which accounts for the effects of amplitude, phase, frequency content, and duration particularly for large‐magnitude earthquakes at closer distances. Because the ground‐motion databases have limited close source observations for large earthquakes, an accurate and realistic ground‐motion simulation that returns the seismic waveforms for future moderate‐to‐large earthquakes has increasingly been required by engineers.

Since the works of Hartzell (1978) and Irikura (1986), numerous methodologies for simulating realistic time series of ground motions have been proposed that range from deterministic physics‐based models to stochastic and hybrid models. Generations of earthquake scenarios may be obtained through deterministic methods, which have been used widely for engineering purposes (e.g., Irikura and Miyake, 2011). Strongly heterogeneous finite‐source models can be constructed stochastically as well (Guatteri *et al.*, 2003, 2004; Schmedes and Archuleta, 2008). Such methods are often known as pseudodynamic approaches and are based on dynamic simulations (Song *et al.*, 2014) or on the statistics of inverted kinematic source models (e.g., Mai and Beroza, 2003). The deterministic procedures are more reliable than the stochastic techniques for lower frequencies in heterogeneous models (Douglas and Aochi, 2008), owing in part to the difficulty in accurately representing the incoherency of wave propagation at higher frequencies. Stochastic models are generally simpler; with more rapid generation and estimation in ground motion. They consequently are more convenient to use than deterministic physics‐based models. Ground motions can also be simulated using hybrid methods that leverage the strengths of a stochastic approach at high frequencies and a deterministic approach at low frequencies (e.g., Berge *et al.*, 1998; Hartzell *et al.*, 1999, 2005; Mai and Beroza, 2003; Liu *et al.*, 2006; Graves and Pitarka, 2010; Mai *et al.*, 2010; Mena *et al.*, 2012).

Ground‐motion simulation in the Marmara region has been conducted by various researchers (Pulido *et al.*, 2004; Sesetyan, 2007; Sørensen *et al.*, 2007; Ansal *et al.*, 2009). All of whom supposed a scenario earthquake of magnitude larger than 7.0 along the main trace of the North Anatolian fault, running along the north side of the Sea of Marmara (Le Pichon *et al.*, 2003; Erdik *et al.*, 2004; Armijo *et al.*, 2005). Increased knowledge of the North Anatolian fault zone within the Sea of Marmara includes historical seismicity coupled with documented paleoearthquakes, instrumental recordings, faulting behavior, and physical parameters of the surface ruptures. Such research has provided an excellent opportunity to test enhanced broadband simulation methodology, as well as more robust estimates of the ground shaking expected in future earthquakes.

In this study, we introduce dynamically simulated earthquake scenarios (Aochi and Ulrich, 2015). In these simulations, the rupture scenarios are simulated along nonplanar structures of the North Anatolian fault, supposing pure strike‐slip faulting, under a tectonic loading. Moreover, the likelihood of macroscopic features—such as rupture area and directivity—are discussed. Each simulation provides an estimation of a particular ground motion. Assembling the different simulations enables discussion and exploration of the variability of ground motion at selected sites produced by specific ruptures. These estimates offer several advantages considering the effects of earthquake source, wave propagation path, and local site effects on ground motion compared with estimating the ground‐shaking levels in potential future earthquake based on empirical GMPEs. In this study, we simulate ground motions for five *M*_{w}>7 events selected according to their high likelihood determined by Aochi and Ulrich (2015) for the Marmara Sea faults. Two independent methodologies, deterministic physics‐based model by Aochi and Madariaga (2003) and Aochi and Ulrich (2015) and stochastic source‐based model by Boore (2009), are used to generate synthetic ground motion. Finally, to obtain the broadband time series, the low‐frequency (LF) and high‐frequency (HF) synthetics were combined to cover the entire frequency band of engineering interest.

## Hybrid Broadband Simulation Method

Broadband synthetic time histories were generated at 2574 virtual stations in the Marmara Sea region using the following hybrid broadband simulation technique. The LF portion of the synthetics was calculated using a full 3D wave propagation method, and a stochastic finite‐fault simulation model was employed to attain the HF portion. Then, the separately calculated HF and LF ground motions were combined at each station to produce more realistic broadband time histories associated with scenario earthquakes having magnitudes *M*_{w}>7.0 located within 30 km south of Istanbul. Rather than implementing kinematic earthquake scenarios, simulated dynamic rupture models along nonplanar faults were used, as reported by Aochi and Ulrich (2015). The wave propagation in the medium was calculated using a finite‐difference approach, as explained below. The HF radiation was computed using a stochastic finite‐fault model approach based on dynamic corner frequency (Motazedian and Atkinson, 2005; Boore, 2009). The results from the two simulation techniques were then merged in the frequency domain following Mai and Beroza (2003). First, LF component was selected to merge with the HF part using the consistency of the plateau level of acceleration in the Fourier space with the criterion around the 1 Hz frequency. The two horizontal LF components were rotated in increments of 1° to determine the best match with the HF part for the random horizontal component of motion. This choice was made to obtain a regular shape for the merged spectrum. Before the merge, a classical long‐ and short‐time average automatic‐picking algorithm was used for synchronizing the time for both signals. Therefore, the result of this process is hereafter referred to as the hybrid broadband horizontal component. Figure 1 depicts an example of the construction of a hybrid broadband horizontal seismogram from the summation of filtered deterministic LF and stochastic HF synthetics.

## Model Parameters

The hybrid broadband simulation method considers the physics‐based rupture process and takes into account different fault rupture scenarios using two different nucleation points and five different slip distributions to investigate their effects on the amplitude and frequency of simulated waveforms.

The physical dimensions of the zone of interest are 200 km (east–west) × 120 km (north–south) with a depth of 40 km, covering the entire Sea of Marmara and Istanbul. The region includes the main North Anatolian fault running close to the northern coastline. The ground motions were calculated using the 3D finite‐difference method (Aochi and Madariaga, 2003), and the 3D structure model (Aochi and Ulrich, 2015) was adopted integrating the recent ocean‐bottom seismograph tomography (Bayrakci *et al.*, 2013) and a regional 1D model (H. Karabulut, personal comm., 2014; Table 1). The Sea of Marmara was considered by fixing the *P*‐ and *S*‐wave velocities at 1500 and 0 m/s, respectively. The receiver points were distributed every 3 km on the surface at a depth of 0 km on flat ground or at depth on the ocean bottom and at 47 distributed seismic stations of the Kandilli Observatory of Earthquake Research Institute (KOERI) network for further investigation. A regular grid was set with a spacing of Δ*s*=200 m and a time step of Δ*t*=0.01 s. The minimum velocity structure involved in the 3D simulations was about *V*_{min}=1000 m/s at the bottom of the Sea of Marmara so that the maximum frequency would be estimated as *V*_{min}/(5Δ*s*)=1 Hz, corresponding to the desired merging frequency.

The HF stochastic simulations were computed using well‐resolved region‐specific attenuation and source parameters derived by Akinci *et al.* (2006) through a set of regressions on weak‐ and strong-motion datasets of 2400 recordings from mainshock and aftershocks from the 17 August 1999 İzmit sequence (*M*_{w} 2.5– 7.2). Akinci *et al.* (2006) also provided a stress parameter of 10 MPa for the mainshock by verifying the predictions against the strong‐motion data recorded during the *M*_{w} 7.4 İzmit earthquake. The spectral source parameters used for the HF ground‐motion simulations are listed in Table 2. Using the input parameters given in the table together with the rupture model and the slip distribution resulting from the dynamic rupture process, as described subsequently, five random horizontal‐component time histories were simulated for each scenario earthquake at each virtual site. The results are given as the root mean square of those five simulations for each site. The time‐domain simulations were iteratively repeated five times following Boore (2005), who reported that the uncertainties in peak ground motions are lower than 10% when the number of iterations is increased from 10 to 640 (fig. 12 in Boore 2005).

## Selection of Possible Scenario Earthquakes in the Sea of Marmara

The likelihood of the probable scenarios in the Sea of Marmara has been previously discussed based on dynamic simulations for quantitative seismic‐hazard assessment (Aochi and Ulrich, 2015). These authors adopted a logic‐tree approach on input model parameters (such as fault geometry models, accumulated stress levels, and hypocenter positions) in which 27 scenarios in total were generated with varying probabilities between 0.5% and 7.5%. Figure 2 presents the statistical distribution of the earthquake scenarios for *M*_{w}>7.0 events, 17 out of 27, based on the logic‐tree approach. The probability was greater than one‐half, with 54% for a rupture beginning from a point located at the center of the model at about 28.5° that propagates bilaterally eastward or westward. Moreover, the probability was 41% for a rupture beginning from the eastern part at about 28.8° E that propagates westward.

In this study, the available statistical information (Fig. 2) was used as the basis for the reference earthquake, identified mainly by the magnitude, location, and resulting slip distribution over the fault. Five reference earthquakes were selected that are most likely to occur according to the probable scenario likelihood used by Aochi and Ulrich (2015) to estimate the level of ground motion around the Istanbul area. In this way, performing a large number of earthquake scenarios at the same time was avoided, reducing the uncertainties in ground motion and saving computational time. The three selected reference earthquakes (scenarios 1, 4, and 5 in Fig. 3a) have the largest probability occurrences with *M*_{w} 7.04, 7.17, and 7.27, respectively; their ruptures propagate mainly eastward from the central point. However, the other two scenarios present unilateral east–west rupture propagation over a longer distance, ending with *M*_{w} 7.19 and 7.26 (scenarios 2 and 3 in Fig. 3b). In the dynamic simulations, the rupture process was spontaneous following the stress field and a prescribed friction relation, including namely a slip‐weakening law, and all of the kinematic parameters such as ruptured area, rupture velocity, and slip‐rate functions were obtained through the simulations. Figure 3a,b also presents the final slip distributions for the five selected scenarios. The slip distribution associated with the LF simulation was smooth and inadequate for the HF part. Thus, the *k*^{2} slip distribution based on a fractal sum of asperities was applied following Ruiz *et al.* (2011) and Herrero and Bernard (1994). To link the HF part to the large wavelength slip model, it was used as a probability density function to draw the asperities on the fault. The result was a *k*^{2} slip distribution compatible at LF with the smooth model. Using the modified slip distribution indicated as the input model in Figure 3a,b, the LF ground motions were calculated at each site utilizing a full 3D wave propagation method for each scenario earthquake. HF stochastic simulations used the same *k*^{2} slip model, in agreement with the LF ground‐motion simulations.

## Simulation Results

In this section, the performances of the hybrid ground‐motion simulations in the Marmara region are discussed. First, the hybrid broadband peak ground acceleration (PGA) and peak ground velocity (PGV) values are assessed to gain insight into the spatial distribution and extent of the strong ground shaking. Then, the hybrid broadband horizontal‐component time histories, and Fourier and response spectra are presented for the five selected sites close to Istanbul. Because no strong‐motion data have been recorded for large earthquakes at these sites thus far, the attenuation trend and the ground‐motion level of the synthetics as a function of distance and frequency with the recent GMPEs were obtained for comparison with the actual data.

### Spatial Distribution of Ground Motion from Hybrid Broadband Simulation

To better capture the ground‐motion variability caused by the source complexities, hybrid broadband synthetic records were generated for 2574 virtual stations every 3 km with the rock‐site classification. The occurrence of site‐amplification effects was not included so that the predicted ground‐motion levels refer to bed/generic rock conditions corresponding to the National Earthquake Hazards Reduction Program (NEHRP) site class A. In Figure 4, the distribution of ground‐motion parameters, predicted PGA, and PGV for the five different scenario earthquakes are examined.

In the scenario 1, 4, and 5 (*M*_{w} 7.04, 7.17, and 7.27) earthquakes (Fig. 4a,d,e), the rupture began at the western and central part of the ruptured fault and propagated mainly to the east toward the city of Istanbul. This generated a long spatial extent of strong near‐fault motion along nearly the entire length of the fault with PGV exceeding 150 cm/s. The PGA distributions in Figure 4a,d,e show that the strongest ground shaking occurred around the location of the large‐slip asperity located in the fault plane. The onshore position of the maximum shaking area is located along the southwestern coast of the European side of Istanbul close to the Bakırköy and Avcılar districts, with PGVs reaching up to 80–100 cm/s. In the case of scenario 4, the rupture propagates bilaterally producing larger ground motions and having slightly larger magnitude, *M*_{w} 7.27, with respect to scenarios 1 and 5.

In the scenario 2 and 3 (*M*_{w} 7.19 and 7.26) earthquakes, the rupture was nucleated in the eastern part of the ruptured fault, close to the easternmost point of the central Marmara fault (Fig. 3b). These two scenarios produced maximum PGA and PGV values higher than those obtained for scenarios 1, 4, and 5. In fact, the PGA and PGV values exceeded 1.0*g* and 250 cm/s, respectively (Fig. 4b,c), close to the fault. The maximum PGA pattern close to the fault was always controlled by the asperity locations in the slip model (e.g., Fig. 3b), varying from 1.5*g* to 0.5–0.6*g* from the eastern to western sectors of the fault (Fig. 4b,c). This tendency is similar to the PGV distribution, especially in the case of scenario 2 (Fig. 4b), with higher levels observed in the eastern sector of the fault. The reason for the lack of significant PGV along the western sector, west of about 28.3° E in the case of scenario 2, is that the rupture in the dynamic rupture simulation was decelerated owing to a small fault jog and subsequent changes in the fault strike based on the complexity of the fault geometry (for the details on the faults geometry and complexity, see Aochi and Ulrich, 2015). The PGV distribution presents more complex behavior owing to the heterogeneous slip distribution and the rupture directivity together with the high‐rupture velocity. In fact, the PGV increased again to about 150 cm/s at the western Marmara Sea region, clearly reflecting the westward rupture‐directivity effect. Another notable phenomenon for these scenario earthquakes is that strong seismic radiation was generated by the supershear rupture along the fault segment in the form of a Mach cone around the rupture front. For this reason, high peak velocities become important even far from the fault because the waves radiated by a supershear rupture hardly attenuate with distance (Bernard and Baumont, 2005).

### Synthetic Broadband Waveforms, Fourier Amplitude, and Response Spectra at Selected Sites

Figure 5 presents an ensemble of the generated broadband acceleration and velocity time histories for the five scenario earthquakes at the five selected sites (BOTS, CTKS, SINB, ISK, and BUYA), which correspond to seismic stations operated by KOERI, Turkey. For the five reference earthquakes, Figure 5a–e shows acceleration time histories of the hybrid broadband synthetics (left panels) and velocity time histories (right panels). Above each signal, we reported the maximum registered PGA (cm/s^{2}) and PGV values (cm/s). Similarly, as an example, Figure 6 shows both Fourier response (5% damping) spectra calculated on both acceleration (left panels) and velocity (right panels) time histories (hybrid horizontal component) from two scenarios (1 and 2) for the five selected sites.

For the scenario 1, 4, and 5 earthquakes (Fig. 5a–c), the Joyner–Boore distance (*R*_{JB}) from the five selected stations to the fault was 13.4 (SINB), 13.7 (BUYA), 22.3 (ISK), 36.1 (BOTS), and 40.4 km (CTKS). The largest PGA values of 0.23*g*, 0.44*g*, and 0.36*g* were observed at one of the closest stations, SINB, located north of the ruptured fault from three scenarios, respectively. Station BUYA is located at about the same distance but showed lower PGA values of 0.11*g*, 0.20*g*, and 0.15*g* for scenarios 1, 4, and 5, respectively (Fig. 5a–c, left panels). The high PGA value of station SINB can be explained by the closeness to the nucleation point and to the position of the main asperity of the rupture process (e.g., Fig. 4a,d,e). On the contrary, in scenarios 1 and 5, BUYA presented a larger ground velocity, at 19.8 and 25.8 cm/s than that of SINB, at 12.5 and 15.8 cm/s (Fig. 5a,e, right panels). This discrepancy is attributed to rupture directivity, where the former station is located in the forward‐directivity area, and the latter is perpendicular to the ruptured fault. This effect was more pronounced in the high periods (>2.0 s) presented at SINB and BUYA acceleration spectra (Fig. 6a, top right panel). Similar phenomena, 10.9 cm/s, were observed at station ISK located in the forward‐directivity area compared with those at station BOTS at 3.4 cm/s. The latter shows lower ground velocity and is located in the back forward‐directivity area (Fig. 6a, bottom left panel).

For the scenario 2 and 3 earthquakes (Fig. 5d,e), the *R*_{JB} from the five selected stations to the fault was 13.8 (SINB), 16.8 (BUYA), 29.1 (ISK), 20.1 (BOTS), and 40.9 km (CTKS). The largest PGA of 0.33*g* was recorded at station SINB, which is the closest to the ruptured fault and to the maximum fault‐slip area for the two scenarios. The directivity effects for scenarios 2 and 3 were consistent with those of scenarios 1, 4, and 5. Stations BUYA and BOTS correspond to comparable distances and showed similar PGAs of about 0.11*g* (Fig. 5d,e, left panels). Nevertheless, these two stations presented quite different PGVs depending on their azimuths with respect to the fault strike; station BOTS, located in the forward‐directivity area, presented a larger ground velocity, at 28.8 and 37.6 cm/s, than that at station BUYA, at 7.6 and 12.8 cm/s for both scenarios, respectively, which is located in the backward‐directivity direction (Fig. 5d,e, right panel). At the same stations, large differences were observed in the Fourier velocity spectra at lower frequencies (<0.5 Hz) depending on their azimuths (Fig. 6b, top panels) for scenario 2. Similar behavior was observed at larger distances. For example, although stations ISK and CTKS are located at distances of *R*_{JB}=29.1 and 40.9 km, respectively, station CTKS showed more than twice the PGV, at 16.0 cm/s, compared with that in station ISK at 6.6 cm/s, located in the backward‐directivity area (Fig. 6b, bottom panels). In general, it is obvious that the stations located east of the epicenter in the backward‐directivity direction showed lower ground‐motion velocities than those located west of the epicenter in the forward‐directivity direction. This effect was quite pronounced in their spectra. Stations BOTS and CTKS, located in the forward‐directivity direction, showed larger amplitudes at long periods of about 3 s owing to the cumulative effect of the radiation from the fault. In contrast, stations BUYA and ISK in the backward‐directivity region showed very low amplitude and no clearly dominant frequencies (Fig. 6b).

### Comparison with GMPEs

The hybrid broadband accelerations and velocities, as well as the spectral accelerations (SAs) simulated from the five scenario earthquakes, were compared with the three different GMPEs for the active shallow crustal regions. These include the GMPE developed within the context of the Next Generation Attenuation (NGA) models developed by Boore *et al.* (2014; hereafter, BA14); the Turkish empirical GMPEs of Akkar and Cagnan (2010; hereafter, AC10); and the Europe and Middle East GMPE model of Akkar and Bommer (2010; hereafter, AB10).

In Figures 7 and 8, the simulated hybrid broadband PGAs and PGVs together with their standard deviations for five scenario earthquakes are plotted up to 100 km as a function of *R*_{JB} and are compared with the selected GMPEs by BA14, AB10 (with standard deviations), and AC10, respectively. All GMPEs were derived from the strike‐slip‐faulting style and for hard‐rock conditions (NEHRP site class A).

To assess the level of fit between our simulations and GMPEs, we calculated residuals *R* at each *R*_{JB}, *j* as (1)in which *Y* is the considered parameter (PGA, PGV, or SA) from the GMPE model. indicates, for the distance *j*, the mean values from the present hybrid horizontal broadband simulations. The model bias *B*(*T*_{i}) for each period *T*_{i} was also computed by summing the residuals over all distances: (2)A perfect match between the empirical model and the broadband simulation would have *R*_{j}=0, whereas a positive residual shows an overprediction of the simulations with respect to the empirical model from GMPEs. Figures 9 and 10 show the plots of PGA, PGV, and SA residuals at 10 different periods versus the distance to evaluate the agreement with the selected GMPEs.

For all scenarios, the simulated ground motions in terms of PGA values together with their standard deviations (±*σ*) fell almost well within ±1*σ* of the empirical model of BA14 (Fig. 7 and Fig. 9, left panels) up to 50–60 km distance. Among the GMPEs used in this study, it was observed that the mean of the simulated PGAs showed better agreement with BA14 and AB10 than that in AC10 which lies close to the lower limit (−1*σ*) of the simulated values at short distances *R*_{JB}<20 km. However, the PGA curve of the AC10 matches better with the trend of the simulated PGAs at larger distances (*R*_{JB}>70 km). The mean of the simulated PGAs was slightly underestimated compared with that of BA14 and AB10 for distances larger than 20 km, whereas the simulations were overestimated with respect to AC10 for distances <70 km. In scenarios 2 and 5, at closer distances the PGAs were more scattered, whereas at higher distances, they also tended to decrease quickly with respect to GMPEs.

As shown in Figure 8, the simulated PGVs were more scattered than the PGAs, which presented challenges in finding good agreement with the empirical models at closer distances (Fig. 8b,c,e). At intermediate and longer distances (>30 km), the PGVs from the simulations were less scattered but decreased faster than those predicted by GMPEs. Only for scenario 1, most of the PGV residuals fell within the range ±1.0*σ* of the BA14 empirical equation at larger distances of ∼70 km. The agreement of the BA14 and AB10 models with the simulations was better than AC10 in the case of PGV. Again AC10 underestimated the simulated peak velocities up to 60–70 km. The standard deviations of the simulated PGVs were significantly larger than those of AB14 at very close distances and decreased at larger distances. As clearly shown in Figure 9 (right panel), the PGV residuals were within almost constant standard deviation ±1.0*σ* of the considered GMPE of BA14 at distances greater than 20–30 km. Hence, the simulated PGVs significantly overestimated all selected, especially that of the AC10 GMPE at short distances (*R*_{JB}<20 km). The sites that fell out of the ±*σ* range were those subjected to full backward‐ or strong forward‐directivity effects and both heterogeneity of the fault slip and rupture velocity at the very close distances to the fault (*R*_{JB}<20 km).

Figure 10 compares the SA estimates of the ground motions simulated in the present study with the predictive equations of BA14 at different periods as a function of distance. The ground‐motion variability was higher than the predicted values of BA14 at long periods >1.0 s and with longer distances, which contrasts with the lower than predicted values at short periods ≤1.0 s at all distances. The variability resembled the pattern observed for all simulations. The simulated SA values were higher than the predicted values of the BA14 equation presenting positive bias for longer periods >1.0 s and were underpredicted compared with the BA14 predictive SA values at short periods <1.0 s over all distances.

The SA modeling bias for all scenario earthquakes was also derived using equation (2) for sites *R*_{JB}<60 km (Fig. 11). The simulations in the present study demonstrate systematic and significant positive bias, showing overprediction with respect to the BA14 predictive equation that increases at higher periods. The SA modeling bias for the scenario 2 and 3 earthquakes for sites with *R*_{JB}<60 km again demonstrated a systematic and significant positive bias, showing maximum overprediction with respect to BA14, reaching its maximum value at 3.0 s (Fig. 11).

A comparison with empirical ground‐motion prediction illustrated that all of the simulated ground motions (PGAs, PGVs, and SA) decayed faster than the global empirical ground‐motion equations, at both moderate and larger distances. This feature was captured by the AC10 model derived from the Turkish earthquake database. The faster attenuation of ground motion owing to the high attenuation of seismic waves in the Marmara region was in agreement with the present horizontal hybrid broadband simulations. This highlights the importance of retrieving specific regional seismic parameters for GMPEs.

Moreover, a slightly larger spread in the PGV was observed, which might be associated with rupture propagation along the fault (i.e., rupture‐directivity and supershear‐rupture effects). These effects were considered in our simulations, although most of the empirical GMPEs tend to smooth them (Ripperger *et al.*, 2008). For the scenario 2 and 3 earthquakes, the maximum PGVs reached 350 cm/s near the fault, which is almost three times larger than the +1*σ* of the empirical model BA14, at ∼100 cm/s. These extremely high PGV values in the very near field directly reflect the slip rate simulated in the dynamic rupture process. Although these two scenarios are considered to be less probable than scenarios 1, 4, and 5, it is stressed that the dynamic rupture approach under a given mechanical setting may produce an extreme scenario not recorded thus far in the area near the fault and not considered in empirical models. Most of the conventional GMPEs do not explicitly account for the characteristics of near‐fault ground motions. Nevertheless, recent studies (i.e., Oglesby and Mai, 2012) demonstrated that the inclusion of such types of records containing high‐velocity pulses in the databases used for developing GMPEs can strongly influence the resulting predictions (Spudich *et al.*, 2014).

## Discussion and Conclusions

In this study, a set of simulations was performed for selected earthquakes having magnitudes *M*_{w}>7 in the Sea of Marmara to explore their variability and level of ground motion in the Sea of Marmara close to the city of Istanbul. Hybrid broadband ground motions were generated and the experiment was made to use full 3D simulations for frequencies lower than 1 Hz and stochastic simulations with dynamic corner frequency for higher frequencies. In particular, the dynamically simulated earthquake scenarios were integrated with fractal complexity directly from the beginning of the process (e.g., *k*^{1} stress distribution).

Ground‐motion parameters including PGA, PGV, and SA (5% damping) at different periods computed at hard‐rock sites were displayed to highlight the azimuthal variations that depend on the fault geometry and rupture model. As expected, the near‐field results were governed by source effects such as fault geometry, the distribution of asperities on the fault plane, and rupture directivity. However, in this study, we did not attempt to investigate which of those parameters provide potentially important insight into the ground motions of the Marmara earthquakes.

Comparison with empirical GMPE demonstrated that the simulated ground motions were consistent with the prediction equation for short (PGA) and medium (PGV) periods for almost all of the scenarios but were mostly underestimated, especially at lower frequencies by the AC10 model derived from the Turkish earthquake database. It is evident that the simulated PGAs decayed faster than the global empirical ground‐motion equations at both moderate and larger distances. The faster attenuation of ground motion can most likely be attributed to the region‐specific low inelastic *Q* parameter or high attenuation of seismic waves in the Marmara region (Horasan *et al.*, 1998; Akinci *et al.*, 2006) used for the HF part of the ground‐motion simulations. This again highlights the importance of retrieving specific regional seismic parameters for GMPEs (Akinci and Antonioli, 2013). The large dissipation of the simulated PGV depends on the position of the nucleation point and corresponds to more than one standard deviation of the considered GMPE. The standard deviation of the simulated PGVs and the SAs (>1–2 s) was larger than that of the BA14 model at very close distances but decreased at greater distances.

However, the simulated PGAs and PGVs better constrained ground‐motion variability and therefore provide more insight than the GMPEs into wave propagation in the region. In fact, none of the selected GMPEs explicitly model directivity effects. However, Spudich and Chiou (2008) developed empirical models and applied the corrective factor for directivity correction to four GMPEs of the NGA project. Moreover, neither the AC10 nor the AB10 relations have corrective factors and thus do not include a directivity term. In this study, only the GMPEs commonly used to assess the probabilistic seismic‐hazard assessments (PSHA) in Turkey are contemplated, without considering the modified GMPEs that include directivity effect (Spudich and Chiou, 2008). Finally, we compare our simulations in terms of PGA, PGV, and SA with those predicted using the not‐modified GMPEs. However, Spagnuolo *et al.* (2016) use the corrective factor proposed by Spudich and Chiou (2008) to improve the PSHA considering fault‐rupture‐related parameters in the Marmara Sea Region in Turkey. Their results demonstrate that rupture directivity tends to increase seismic‐hazard estimations up to 25% and can predict a maximum reduction of 15% along the two considered segments, the central Marmara fault and the Çınarcık fault of the North Anatolian faults system at 2 s. The effectiveness of the correction for directivity depends on the period of the SA, and its contribution is higher at larger periods. In all scenario earthquakes, hybrid broadband simulations, large PGV and SA at larger periods are mostly confined in the near field within 30 and 20 km from the fault, respectively. Nevertheless, directivity effects were prevalent, particularly for scenarios 2 and 3 earthquakes. They could be one of the responsible of the overestimation of long‐period ground motions (∼3 s) compared with standard GMPEs in which the directivity effect is not taken into account.

A strong‐directivity or supershear‐rupture propagation still constitutes an extreme scenario in which the impact on ground‐motion simulations requires improvement at both high and low frequencies. However, the present results again highlight the importance of retrieving specific regional seismic parameters and the consideration of rupture‐directivity effects in GMPEs (Shahi and Baker, 2011; Spudich *et al.*, 2014; Spagnuolo *et al.*, 2016). Moreover, advanced source models that consider dynamic effects of rupture propagation may lead to improvement in simulation results and may be used to minimize the uncertainties owing to the physical parameters of source rupture and wave propagation. These results have important implications for estimating seismic hazards in the Marmara Sea region and emphasize the need for improved understanding of earthquake rupture processes.

In all simulations, the average PGA and PGV distributions showed that the maximum onshore‐shaking area was located at the southeastern part of the city along the coast near Bakırköy and Avcılar districts, 15–20 km from the fault rupture, and the values of PGA and PGV were about 0.5–0.6*g* and 60–70 cm/s, respectively. This area is also densely populated and hosts several critical facilities such as Ataturk International Airport and industrial installations. Hence, identifying the upper and lower limits of the expected ground motions may be useful in risk mitigation strategies to provide important indications for decision makers in the study area.

In the present study, local site effects in this area were not considered. However, previous studies, following the 1999 İzmit earthquake, demonstrated the effects of the soft sediments located in the Bakırköy and Avcılar districts, which experienced significant damage caused mainly by the strong amplification owing to the site effects during the *M*_{w} 7.4 İzmit earthquake (Ozel *et al.*, 2002; Birgören *et al.*, 2004; Sørensen *et al.*, 2007). Therefore, understanding linear and particularly nonlinear site responses during large earthquakes is crucial for better modeling of ground motions to assess future hazards in the city of Istanbul and the surrounding areas.

## Data and Resources

Model information used for hazard calculations including Tables 1 and 2 came from many published sources listed in the References. Many of the plots are made using the Generic Mapping Tools v. 4.2.1 (www.soest.hawaii.edu/gmt, last accessed December 2008; Wessel and Smith, 1998).

## Acknowledgments

This study is supported by the Marmara Supersite (MARSite) “New Directions in Seismic Hazard Assessment through Focused Earth Observation in the Marmara Supersite,” European Integrated Project, THEME‐ENV.2012.6.4‐2 (Long‐term monitoring experiment in geologically active regions of Europe prone to natural hazards: the Supersite concept), Grant Agreement Number 308417. Numerical simulations at low frequencies were carried out at the French National Supercomputing Center, Grand Equipement National de Calcul Intensif, under Grant Number 46700. We thank Martin Mai for very helpful comments and suggestions; and an anonymous reviewer for thorough reviews that largely contributed to organizing and improving the article. We also thank Matteo Taroni of Istituto Nazionale di Geofisica e Vulcanologia for insightful discussions.

- Manuscript received 24 March 2016.