Seismological and engineering characteristics of strong motion data from 24 and 26 September 2019 Marmara Sea earthquakes

An MW 4.5 earthquake took place on September 24, 2019 in the Marmara Sea. Two days after, on September 26, 2019, Marmara region was rattled by an MW5.7 earthquake. With the intention of compiling an ample strong ground motion data set of recordings, we have utilized the stations of Istanbul Earthquake Rapid Response and Early Warning System operated by the Department of Earthquake Engineering of Boğaziçi University and of the National Strong Motion Network operated by AFAD. Altogether 438 individual records are used to calculate the source parameters of events; namely, corner frequency, radius, rupture area, average source dislocation, source duration and stress drop. Some of these parameters are compared with empirical relationships and discussed extensively. Duration characteristics are analysed in two steps; first, by making use of the time difference between P-wave and S-wave onsets and then, by considering S-wave durations and significant durations. It is observed that they yield similar trends with global models. PGA, PGV and SA values are compared with three commonly used ground motion prediction models. At distances closer than about 60 km, observed intensity measures mostly conform with the GMPE predictions. Beyond 60 km, their attenuation is clearly faster than those of GMPEs. A frequency-dependent Q model is developed by using both events. Its consistency with existing regional models is confirmed.


Introduction
Istanbul was hit by a moderate size earthquake of M W 5.7 on September 26, 2019, which was preceded by a M W 4.5 event approximately two days earlier. It was the largest earthquake that took place on the central Marmara segment of the North Anatolian Fault (NAF) in recent years. The earthquake was widely felt across Istanbul and caused structural 1 3 damage in private and public buildings. Mobile communication was lost for about two hours. It has stirred many public and administrative discussions and evaluations, and reignited the long-standing question of whether a destructive earthquake with a larger magnitude will occur in the near future. In the past, two destructive earthquakes that occurred on the NAF, immediately east of Istanbul (M W 7.4 in August 1999 andM W 7.2 in November 1999), caused thousands of deaths along with excessive structural damage. It is crucial to advance the existing knowledge in understanding the potential hazard and risk that a major earthquake may cause in this region. Existing strong ground motion networks in Istanbul are rich sources of seismological information for studying the past and future potential earthquakes in this region.
In this study, we focus on estimating seismological and engineering parameters of these two events. Seismic activity on the segments of the NAF to the south of Istanbul has been significantly low. Between 1999 and 2020, there were only nine earthquakes with magnitudes 4.5 and above that could be associated with the central Marmara and Princes' Islands segments (KOERI-RETMC 2020). We used recordings from the strong motion stations in and around Istanbul. We estimated the following source parameters: corner frequency, source duration, source radius, rupture area, average source dislocation and stress drop. We have examined the regional distance dependence of the duration between P-wave and S-wave arrivals. Then, we have studied possible correlations between the hand-picked S-wave durations with hypocentral distance. The significant duration models have been prepared for 5-75% and 5-95% cases separately and the results are compared with the global models. Next, ground motion parameters (PGA, PGV, SAs at various periods) have been calculated and a discussion on engineering parameters has been presented. They have been compared with three commonly used ground motion prediction equations (GMPEs) for the region. In the last part of the paper, we have discussed the validity of frequencydependent anelastic attenuation models developed by comparing them with global and regional models proposed in the literature.

Earthquake information
An earthquake measuring 4.5 on moment magnitude scale was detected off the coast of Silivri (Sea of Marmara, northwestern Turkey) on 24th of September at 11:00:21 local time. Two days later, another earthquake with a moment magnitude of 5.7 struck the same region. Brief information on both earthquakes is summarized in Table 1. Their depths are estimated as 9.9 km and 13.3 km, respectively, by Kandilli Observatory and Earthquake Research Institute (KOERI) (KOERI-RETMC 2019a; b). The Ministry of Interior's Disaster and Emergency Management Presidency (AFAD) reported slightly larger moment magnitudes and smaller depths than KOERI for both events (AFAD 2020). Moreover, the United States Geological Survey (USGS) and European Mediterranean Seismological Centre (EMSC) developed similar body wave magnitudes (M b = 4.9) and depths (≈ 10 km) for the smaller event (USGS 2021;EMSC 2021). In the case of 26th September earthquake, the moment magnitudes stated by both organizations are identical to those of KOERI (KOERI-RETMC 2019a;. The focal mechanism solution of the 24th September earthquake is compatible with the fault mechanism of NAF's central Marmara segment and exhibits a strike-slip dominancy (KOERI-RETMC 2019a; Karabulut et al. 2021). However, there are different views about the faulting mechanism of the 26 September 2019 event. While the solutions by KOERI-RETMC (2019b), AFAD (2020), USGS (2021) and

3
Eyidoğan and Sevilgen (2019) indicate dominantly thrust faulting with dipping planes ranging from 35 to 55 degree, according to Karabulut et al. (2021) the event is dominantly strike-slip and took place on a 60-degree north-dipping fault.

Characteristics of strong ground motion database
The strong ground motion recordings from the networks of the Ministry of Interior's Disaster and Emergency Management Presidency (AFAD) and the Istanbul Earthquake Rapid Response and Early Warning System (IERREWS) of Boğaziçi University's Kandilli Observatory and Earthquake Research Institute have been utilized in this study. For both events, the number of selected stations and the epicentral distances (R epi ) of closest and furthest recordings are summarized in Table 2.
A total of 70 and 76 stations located in and around Istanbul were selected to analyse the M W 4.5 and M W 5.7 earthquakes, respectively. The spatial distribution of selected stations and earthquakes' epicentres including their focal mechanism solutions are illustrated in Fig. 1. 1 3

Data processing
All selected recordings are processed and analysed individually. Firstly, to provide physically meaningful velocities and displacements, all raw acceleration time histories are uniformly detrended. P-wave arrivals, S-wave windows and noise windows are hand-picked following Kishida et al. (2016) as shown in Fig. 2. The majority of recordings are sufficiently long enough to select the pre-event noise window, similar to the length of the hand-picked S-wave window. Fourier spectra of S-wave and noise windows are estimated to check the signal-to noise ratios. S-wave spectra and noise spectra multiplied by the factor of three are excessively smoothed and checked for any intersections between them. The frequencies of intersection points (if they exist, most of the time the signal is way above the noise level) are used as the low-pass corner of the filters. It has been noted that the recordings in the database are mostly of high quality with low noise levels and a wide usable frequency band-reaching up to 50 Hz. Each acceleration trace is examined one by one and individual high-pass corner frequencies are determined for each recording separately via visual inspection of the velocity and displacement waveforms. As the final step of data preparation, band-pass filtering is applied to the recordings using the individually selected low-pass and highpass filter limits in order to eliminate the noise contamination. However, it should be noted that the low pass filter frequency limit is kept constant at 20 Hz only in the calculation of engineering parameters (PGA, PGV, and SA values) by considering the total duration of the time histories in order to obtain comparable information in the short period range of response spectrum.

Source parameters
To estimate the corner frequency of the source models, the S-wave spectra only with very high usable bandwidths are used to obtain the source spectra. For this estimation, the spectral decomposition of S-wave Fourier amplitude spectra has been applied by correcting the spectra for geometrical spreading, site response and high-frequency attenuation effects. We consider the geometrical spreading effect by a simple 1/R model. The site response is accounted for by frequency-dependent amplification factors estimated by Boore and Joyner (1997). The high-frequency attenuation has been measured for each individual recording, following the model proposed by Anderson and Hough (1984). The corner frequency is calculated using Eq. (1), introduced by Andrews (1986), using the source spectra for the 0.1-10 Hz range. This equation is selected due to the equation's objectivity in calculating the corner frequency: where V(f) and D(f) are the Fourier amplitude spectra of velocity and displacement, respectively. All individual f c calculations are averaged over all components to obtain the corner frequencies of these events using the equation: where f c,i is the corner frequency of the ith recording and N refers to total number of recordings. To provide further accuracy, we have plotted the mean f c over the displacement Fourier amplitude spectra of all components for both events (Fig. 3). The source radii of the seismic sources are estimated through Brune's model (1970) as, where is the crustal shear wave velocity in km/s and it is assumed as 3.42 km/s for the Marmara region (RETMC, personal communication 2018). Once the radii are known, a circular rupture plane area can be estimated with the following equation: The seismic moment (M 0 ) of the M W 4.5 event is reported as 6.494 × 10 15 Nm by KOERI-RETMC (2019a). The conversion from M W to M 0 for M W 5.7 event is applied by inverting the well-known relationship as defined by Hanks and Kanamori (1979), The shear modulus can be calculated by using the equation of = 2 . Assuming the density of earth crust as ρ = 2.7 g/cm 2 , the shear modulus is calculated as μ = 3.16 × 10 10 Pa . Thus, from the equation of scalar seismic moment, the average source dislocation, D, can simply be estimated by formulating the seismic moment equation as, Several source duration models are proposed in the past (Boatwright and Choy 1992;Hanks and McGuire 1981;Beresnev 2002). All of these models assume that source duration is inversely proportional to f c , while many of them use different coefficients in their formulations. The approximation presented in Boatwright and Choy (1992) and in Boore (1983) is used in this study, The static stress drop, ∆σ, is the difference in shear stress on the fault plane before and after the earthquake. Based on Eshelby (1957) and Keilis-Borok (1959), the following relationship holds for a circular fault rupture with a homogeneous stress drop: The unit of stress drop in this equation is in Pascal (Pa = 10 -5 bar). The results obtained for the M W 4.5 and M W 5.7 earthquakes regarding the source parameters are listed in Table 3.
The relatively high corner frequency (≈ 0.5 Hz) estimated for M w 5.7 event is directly affecting the values of succeeding parameters. However, Irmak et al. (2021) observed a spectral decay at 0.4 Hz that is observed in coda-wave fitting process for a different dataset than ours. Approximated rupture areas (A) are calculated also with two empirical source scaling relationships developed by Wells and Coppersmith (1994) and Thingbaijam et al. (2017), referred as WC1994 and T2017, respectively (see Table 3). Rupture areas calculated through these relationships for the M W 4.5 event are slightly underestimated, while those for M W 5.7 are significantly different than the ones presented in Table 3. However, our circular rupture plane area estimation is matching well with the slip area proposed in Karabulut et al. (2021).
The estimation of 0.71 m as the average source dislocation of the M W 5.7 event is in good agreement with the slip range considered in Fig. 7 of Karabulut et al. (2021). Average source dislocations calculated through T2017 are shown in Table 3 as well.
The calculated stress drop value of 122.2 bar for the M W 5.7 event is relatively high when compared to stress drop estimations globally. Nevertheless, some high values are observed for earthquakes that occurred in the Sea of Marmara (Irmak et al. 2020), while high stress drop values also exist globally in shallow strike-slip and thrust fault earthquakes (Allmann and Shearer 2009). We have also checked the consistency of this value with a Δσ model suitable for elliptic ruptures. For an elliptic rupture, we assume the major-axis, denoted as a, ranging from 3.5 to 5 km and the minor-axis, denoted as b, ranging from 1.3 to 1.8 km -each a and b combination resulting in approximately the same area to that calculated by assuming a circular rupture (Table 3). a and b are considered to have approximate proportionality with the length and width estimations for the rectangular rupture assumptions for these events in WC1994 and T2017. We calculated Δσ by using the equation B12 in Denolle and Shearer (2016), arranged for elliptic slip distribution based on Aki (1972): where C = 16∕(3π 3∕2 )∼0.95. For a set of axes combinations mentioned above, we obtain Δσ values ranging from 117.7 to 158 bar. The Δσ estimated for a circular rupture (122.2 bar), falls within this range-showing consistency between the results. 6 Duration of strong ground motions

Duration between P-wave and S-wave arrivals
The first part of this section deals with the duration between hand-picked P-wave and S-wave onsets. These durations are detected in each recording individually. The P-wave and S-wave onsets are selected by plotting all three components of the recordings one under the other simultaneously for a single visual selection. Figure 4 clearly shows a robust linearity between the durations from P-wave arrival to S-wave onsets and hypocentral distance (R hyp ), which yields a regression coefficient of R 2 = 99%. The fitted model intersects the R hyp -axis at 0.34 s. Due to its insignificant weight in the equation, it can be ignored and only the path contribution can be considered as: It can be said that the estimated linear model is in good agreement with the global linear model that is R hyp /8, used in previous studies (Ancheta et al. 2013;Kishida et al. 2016).

Duration of S-waves
Following the S-wave windowing in the Data Processing section, the individual S-wave duration of each recording for the two events are plotted against hypocentral distance in Fig. 5. The data points are slightly scattered because of the subjectivity in the selection of the end of the S-wave windows, yet, they are following a linear trend supported by R 2 = 49% and R 2 = 46% for the M W 5.7 and M W 4.5 events, respectively. However, given the subjectivity in the selection of S-wave window edges and the lack of the data at close distances, any regression model applied to this dataset will lack accuracy.
(10) D P−wave arrival−S−wave arrival = 0.118 × R hyp Fig. 4 Duration between hand-picked P-wave and S-wave arrivals versus hypocentral distance of both events; the solid line is a first-degree polynomial fitting, while the dashed lines represent the 95% confidence intervals; the fitting equation is shown in the legend with ± 95% confidence intervals

Significant duration
Duration of strong shaking defined through various methods is considered as the combinations of source, path and site contributions (Kempton and Stewart 2006;Bommer et al. 2009;Ghofrani and Atkinson 2015). In this section, we only focus on the path contributions estimated through significant duration approach (Arias 1970). Significant duration is defined as the time interval across which a certain amount of energy is dissipated. Arias (1970) utilize the integral of the square of the ground acceleration to represent energy, which is commonly known as Arias Intensity (I A ), where g is the gravitational acceleration, T is total duration of strong ground motion and a(t) represents acceleration time history. We have considered the two widely used measures of significant durations, which are expressed as the time interval needed to accumulate 5-75% and 5-95% of I A . They will be denoted as SD 5-75 and SD 5-95 hereafter.
The significant durations for the average horizontal and vertical components are presented against hypocentral distance in Figs. 6 and 7. The linearity of these distributions and fitted first-degree polynomial lines are checked as well. The R-squared values of the distributions of horizontal components are relatively lower than that of the vertical component. Therefore, the 95% confidence intervals of the fitted first-degree polynomials are higher in horizontal components than in the vertical component (Table 4). In this regard, the highest R-squared values with tightest confidence intervals are of SD 5−75 of the vertical component, showing the best compatibility with linear regression. The path components of the regressed models are 0.13 ± 0.03 and 0.12 ± 0.03 for M W 4.5 and M W 5.7 events, respectively. This term is 0.1 for events with M3.1-6.7 in Southern California (Raoof et al. 1999) (11)  Kempton and Stewart (2006) for SD 5-75 and SD 5-95 , respectively. The path component of S-wave durations is modelled linearly as 0.1 × R hyp for California, in Kishida et al (2016). These estimations are in good agreement with our results. However, the same term is 0.05 for moderate magnitude earthquakes in eastern North America (Atkinson, 1993). From visual check of individual time series and from duration distributions against hypocentral distance, we infer that the SD 5-75 appears to be a favourable representation of the strong shaking window for these earthquakes, whereas the SD 5-95 is a better estimate regarding the total duration of the motion.   The blue and purple data points show the durations needed to accumulate 5-75% and 5-95% of the energy, respectively; whereas, the error bars represent µ ± 1σ for each range of R hyp every 10 km  Figure 8 feature the spatial distribution of stations considering the geometric means of peak ground accelerations (PGAs) for both events. The figure explicitly unveils the decrease in PGAs with the rise in epicentral distances while the stations at the immediate north coast of Marmara Sea, which are the closest mainland to the epicentres, present relatively larger PGAs. It is striking that beyond 60 km of R epi , there exists no PGA greater than 4 cm/s 2 for M W 4.5 event and similarly, only 4% of the stations exhibit PGAs greater than 30 cm/s 2 for M W 5.7 event. Also, it should be stated that although the stations of 1 3 FBZMD and 3411 are positioned slightly further than 60 km of R epi , they produce marginally higher values than 30 cm/s 2 .

Engineering parameters
The East-West (EW) and North-South (NS) acceleration traces of both events are presented for two spatially close stations in Fig. 9. Station 3408 is the closest station to both epicentre while stations 5906 and SINB, which are deployed on softer soil site (V s,30, 5906 ≈ 224 m/s; V s,30, SINB ≈ 313 m/s) and slightly further away from 3408 (V s,30, 3408 ≈ 639 m/s), yield the largest PGAs for both M W 4.5 and M W 5.7 earthquakes, respectively. A more detailed discussion on the peak ordinates of ground motions will be provided in the further paragraphs.
The appraisal of the geometric mean of PGAs and peak ground velocities (PGVs) has been performed by considering the epicentral distances and average shear wave velocities for the upper 30 m depth (V s,30 ) of each recorded station (Figs. 10 and 11). The preliminary findings demonstrate that soft soils and proximity to the epicentre increase the PGAs and PGVs when compared with the further-distance and stiffer soil records. However, a detailed examination unveils some exceptions which will be detailed in what follows.
For M W 4.5 earthquake, one of the closest stations of AFAD, station 5906 (R epi = 26 km; V s,30, 5906 ≈ 224 m/s), deployed on relatively soft soil sites provides the largest geometric mean PGA as about 12.2 cm/s 2 . At a similar distance, another AFAD station (station 3408: Fig. 9 EW and NS components of acceleration traces for a the stations of 5906 and 3408 for M W 4.5 event, b the stations of SINB and 3408 for M W 5.7 event R epi = 21 km; V s,30, 3408 ≈ 639 m/s), the closest one to the epicentre, gives a lower geometric mean (PGA g,3408 ≈ 5.9 cm/s 2 ) than the station 5906. This may be associated with relatively stiffer soil conditions of station 3408 compared to station 5906. However, component-based peak values reveal that PGA (PGA EW,3408 ≈ 10.3 cm/s 2 ) of the EW component of the closest station is identical to that of station 5906 (PGA EW,5906 ≈ 12.7 cm/s 2 ) while its NS component exhibits remarkably lower PGA (PGA NS,3408 ≈ 3.4 cm/s 2 < PGA NS,5906 ≈ 11.8 cm/s 2 )  Figs. 9a and 10a). Another striking point is that the EW components at stations having epicentral distances less than 50 km seem to be larger than their NS counterparts except for the AVCLI station. Moreover, the PGV of the EW component of station 3408 is slightly greater than that of 5906, while the difference between their PGVs of NS components is significant (PGV NS,3408 ≈ 0.14 cm/s < PGV NS,5906 ≈ 0.40 cm/s) (Fig. 10b). It should also be noted that the stations associated with higher PGA and PGV values, are at locations very close to the Fig. 11 a Variation of PGAs and elastic response spectra of single degree of freedom systems (SDOF) for selected records and b variation of PGVs as a function of epicentral distance and V s,30 for the M W 5.7 earthquake northern coast of the Marmara Sea. Another feature that comes forward from Figs. 10a and 10b is the seeming existence of two zones that are marked with strikingly different average PGA and PGV levels separated at about 60 km. We will elaborate more on this in the next section.
Similar to PGAs and PGVs, spectral accelerations (SAs) reach their maxima in EW components (Fig. 10a). For closer stations with V s,30 values lower than 250 m/s (3412 and 5906), elastic acceleration response spectra show a bump at periods near 0.1 s. However, the response spectra of recordings from the other three stations (3408, AVCLI, SINB; V s,30 > 250 m/s) with largest PGAs and PGVs give maximum accelerations at the period of around 0.2 s.
From the PGA distribution of the M W 5.7 event, it is seen that recordings from four stations (5906, 3408, SINB, HVHR) generate relatively high PGAs (Fig. 9b and Fig. 11a). SINB of IEEWS, the third closest station to the epicentre (R epi ≈ 31 km) registers the largest PGA as 145.8 cm/s 2 in the EW component. Although 5906 is the second closest station and located on a site with the lowest V s,30 value (V s,30 ≈ 225 m/s) among the closest stations, the PGA of its EW component is lower than that of SINB, which possesses a slightly larger V s,30 and is at a similar epicentral distance. However, their NS component-PGAs are about identical. Furthermore, the closest station (R epi ≈ 22 km), 3408, which is deployed on relatively stiff soil (V s,30 ≈ 639 m/s), produces lower PGAs (PGA EW,3408 ≈ 51.3 cm/s 2 ) in the EW direction than the aforementioned stations at equivalent distances. However, its PGA in the NS component (PGA NS,3408 ≈ 83.7 cm/s 2 ) appears to be slightly higher than those of the other two stations. The PGV distribution is shown in Fig. 11b and includes interesting features worth discussing. The largest PGV is registered in the EW component of station SINB similar to its PGA; however, the maximum of geometric mean PGV is associated with station HVHR deployed at an intermediate epicentral distance (R epi ≈ 54 km). As a general observation the largest PGVs are recorded at stations at epicentral distances less than 60 km and that have V s,30 values less than 350 m/s. The two zones separated at about 60 km epicentral distance and characterized by two distinctly different average PGA and PGV levels realized for the case of the M w 4.5 event, exists for this larger event as well. As mentioned earlier they will be discussed further with the help of GMPEs.
For the M W 5.7 earthquake, response spectral amplitudes follow the trend observed for the PGAs (Fig. 11a). The EW component of station SINB generates a very high salient peak SA (≈ 558 cm/s 2 ) at 0.23 s. However, its NS component produces compatible SAs with other three stations (3408, HVHR and 5906). NS components of the closest stations, 3408 and 5906, also have high SAs at the periods around 0.1 s. The EW components of these stations have SA peaks around 0.2 and 0.25 s, respectively. It is worth noting that stations 3408 (V s,30 ≈ 639 m/s) and 5906 (V s,30 ≈ 225 m/s) produce similar levels of SAs at similar frequencies although they have different soil conditions. Furthermore, while the EW component of HVHR station experiences two high SA peaks between 0.3 and 0.6 s periods, the maximum SA (≈ 300 cm/s 2 ) arises in NS component around 0.65 s, which corresponds to relatively longer periods.
Discussions above motivate the comparison of recorded ground motion peaks with established GMPEs, which is the subject of the next section.

Selected ground motion prediction equations
PGAs, PGVs and SAs are compared with three commonly used GMPEs in the region to understand the attenuation characteristics of observed data as a function of epicentral distance. Selected equations are as in the following: Boore et al. 2014 (BSSA14), Akkar et al. 2014 (ASB14) and Kale et al. 2015 (KAAH15). Information on the input parameters of these expressions and the values that we used as well as our assumptions in employing them are given below.
• Moment magnitude (M W ) Moment magnitudes are taken as 4.5 and 5.7 for 24th September and 26th September earthquakes respectively, as given by KOERI. • Distance (R) BSSA14 and KAAH15 adopt Joyner-Boore distance (R JB ), while ASB14 uses R epi . For the M W 4.5 event, due to lack of a fault plane solution, we have assumed that R JB equals to epicentral distance (R epi ). As for the M W 5.7 event, Karabulut et al. (2021) suggested an elliptic finite-fault model, whose rupture zone has a length of 5 km on the east and 3 km on west of the epicentre. There is about 8% and 3% decrement in distance when R JB with respect to its surface projection is considered instead of R epi for stations 3408 (closest) and YLVH (farthest) respectively. This is a difference which will have a relatively insignificant effect on our general evaluations. Therefore, we have preferred to use R epi for the M W 5.7 event as well. Eventually R epi has been the distance definition that we adopted for the three GMPEs. According to NEHRP site classification, this value corresponds to "Soil Class C: very dense soil and soft rock", whose V s,30 values are between 360 -760 m/s. • Region BSSA14 and KAAH15 offer regional options. Parameter values corresponding to Turkey have been utilized in the generation of GMPE curves. ASB14 does not specify any option for any region. However, since ASB14 were derived from strong ground motion records obtained in the Mediterranean region and the Middle East, they could directly be used without any regional modification.

Variation of peak ground motion parameters
Median level PGAs and PGVs estimated by the three previously mentioned GMPEs, as well as corresponding bounds for standard deviations (σ) (median ± σ and median ± 2σ lines) are shown in Fig. 12. There are no records at distances less than 10 km. Therefore, the considered distance range is 10-100 km in Fig. 12. The recorded PGAs and PGVs plotted correspond to the geometric mean of two horizontal components (EW and NS). When PGAs are compared with the estimations of BSSA14 for both earthquakes, observed values up to about 60 km mainly remain within the median ± σ interval. However, with the increase in epicentral distance, observed PGAs start to decay faster than BSSA14. While a few observed PGAs approach to the median-2σ line (2.3rd percentile) for the smaller event, for the larger event, there is an observed PGA cluster below − 2σ line after about 70 km. For M W 4.5 earthquake, ASB14 and KAAH15 appear to have a better agreement with observed PGAs varying between 16 and 84 th percentile (median ± σ) prediction lines with a few exceptions. PGAs of M W 5.7 event seems to be underestimated by ASB14 and KAAH15 up to about 70 km of distance, beyond which they begin to approach the median lines and at further distances, they arrive at 16th percentile line. To sum up, for both earthquakes, faster attenuation is unveiled in the observed PGAs especially above the distances around 60-70 km compared with GMPE predictions. Herein, it is noteworthy to mention that Sect. 8.4 scrutinizes the estimation tendencies of these GMPEs more quantitatively through the evaluation of residuals for all regarded ground motion parameters.
The PGVs of both events are more dispersed than the PGAs. For M W 4.5 earthquake, observed PGVs are within the median ± σ band of three selected GMPEs up to epicentral distances around 60-70 km. Beyond this distance, PGVs display a sharp reduction down to the median-2σ line. The dispersion of PGVs is larger in the case of the M W 5.7 event. Similar to the observations we made about the PGAs and the PGVs of the M W 4.5, two distance ranges are easily identified, in which the PGV attenuation has different trends. They seem to separate at about 50 km. At distances smaller than approximately 50 km, the PGVs follow the trend set by the GMPEs, although the majority of them are above the median. After 50 km, they start to attenuate faster than the GMPE estimates, on a general path that quickly approaches the median-2σ line.

Variation of spectral accelerations
SAs have been examined for four different periods (T), namely, T = 0.4; 1.0; 2.0 and 4.0 s. In order to pinpoint explicitly the extent to which observed spectral amplitudes in different periods conform to GMPEs, the observed and estimated SAs are depicted in Figs. 13 and 14 for M W 4.5 and M W 5.7 events respectively.
At first glance, the notable feature that can be detected in Fig. 13 is the over-prediction of SAs by three selected GMPEs after about 60 km of distance for the M W 4.5 earthquake. On the other hand, for short periods, i.e. T = 0.4 s and 1.0 s, the observed SAs are within the range of 2.3rd and 98th percentiles but the majority of the data are very close to the latter. With the increase in the period, the overestimations become clearer and observed values are descended below the 98th percentiles. It should be noted that these judgements are more consistent particularly for the IERRS data, which contains the most data. The residual analysis section goes into considerable detail concerning network-based assessment. Before 60 km, SAs are generally within one standard deviation range of three selected GMPEs for all periods.
Beneath 60 km, the majority of SAs of the M W 5.7 event take larger values than GMPEs' estimations except for SAs at T = 4.0 s, where a relatively better agreement with median lines is evident (Fig. 14). Beyond this distance threshold, SAs start to attenuate faster. While this is more explicit in SAs with T = 0.4 s, SAs at greater periods (T = 1.0 and 2.0 s) are assembled around the 16 th percentile line and the 2.3 rd percentile for T = 4.0 s. At the furthest distances, wide dispersion of observed SAs which extend beyond the 2.3rd and 98th percentile predictions is evident.

Evaluation of residuals
In Kuehn et al. (2016), the comparison of various regionalized models unveils the different attenuation features of Turkey particularly with magnitude and distance and also emphasizes the smaller median estimations of the region-specific models than other regions. A recent study by Kotha et al. (2020)  regions consisting of eight regions from Turkey and particularly highlights that the geological differences of Turkey affect regional ground motion predictions. The results of the study points out the need for regional ground motion models since residual variances increase with the rising contribution of data from different regions when compared the previously proposed pan-European ground motion models. These evaluations have urged us to the more detailed examination of the observed variability of the events. Therefore, an analysis of residuals is employed to assess the variation of misfit between estimated and actual ground motion parameters with the epicentral distances, soil conditions, and source-to-site azimuths. The total residuals R GMP,ij refer to differences between the natural logarithm of the estimated ground motion parameter (GMP) and its observed value as in Eq. (12) and can be decomposed into between-event (event term or inter-event) E i and within-event (intra-event) W ij components (Al- Atik et al. 2010).

considers 46 territories from European and Mediterranean
where the related parameters above refer to values at the j th station of i th event. The standard deviation of the total residual (σ) can also be separated into within-event (ϕ) and between-event (τ) terms. In the case of single event, the between-event term E i , which reflects the event-to-event variability, become constant and the record basis variation of total or within-event residuals may provide us comprehension of the the potential effects of distances, soil conditions and azimuths.
For the M W 4.5 and M W 5.7 events, the total residuals of PGAs, PGVs, and SAs for four different periods (T = 0.4; 1.0; 2.0; and 4.0 s) with epicentral distance (R epi ) are illustrated in Fig. 15a and b, respectively. It should be noted that the station-specific V s,30 value is considered in the calculation of ground motion parameters estimation instead of an average value of V s,30 . At the first glance, the dominancy of negative residuals displays the overestimation tendency of selected GMPEs for both earthquake and it becomes more explicit in the SAs especially with the increase in periods. For the M W 4.5 event, although the total residual values of PGAs for all three GMPEs appear to be more or less steady through all distances, after 60 km, they are resulted in a reduction in residual values, with BSSA14 having the lowest residual (R PGA,BSSA14 ≈ − 1.64 at R epi ≈ 74 km). Moreover, the ASB14 and KAAH15 models yield mean PGA residuals (R PGA,ASB14 ≈ − 0.24 and R PGA,KAAH15 ≈ − 0.09) closer to zero than BSSA14 (R PGA,BSSA14 ≈ − 0.55), implying that the formers are more in line with the observational data. Similar considerations may apply to PGV and SAs. However, the previously mentioned decrease beyond 60 km is more pronounced in these parameters; even in SA (T = 4 s), the beginning of this residual reduction shift to closer distances (R epi ≈ 40 km). In the case of M W 5.7 event, all ground motion parameters also verify the decrease beyond 60 km and beneath this distance, the residuals appears to be more dispersed for all ground motion parameters.
The accelerated rise of residuals in the negative direction after roughly 60 km validates the previously stated argument that actual ground motion parameter values attenuate faster. Furthermore, highlighting the discrepancy in attenuation characteristics of Turkey in several studies, which are summarized in the beginning of the section, strengthens this argument. Nevertheless, in order to back up this argument, the tendency of the residuals variation with distance on a network basis has been investigated by performing simple linear regression ( Residual = a + bR epi ; herein a and b regression coefficients) for three GMPEs and all ground motion parameters. Since the overall trend is similar for all three GMPEs, for only KAAH15, this evaluation is presented as an example in Fig. 16, but also only for the 4 s period SA, in which the decrease is most noticeable. It should be stated that the numerical values of the regression coefficients are not supplied herein because the primary priority is to determine the general trend. Each network's residuals precisely lessen with distance for M W 4.5 event; however, for the larger event (M W 5.7), residuals of the IEEWS and AFAD data showcase a more stable attitude with distance for SA at T = 4 s. However, for PGA, PGV and SAs at the other regarded periods, similar results are obtained for both seismic events. Based on the residual analyses of M W 4.5 event, the previously highlighted argument, which point out the different attenuation features beneath and beyond 60 km, may be asserted for whole data regardless of the network; however, the dispersed feature and scarcity of AFAD and IEEWS data hinder accurate individual evaluations on the basis of these networks while the faster attenuation tendency of the ground motion parameters is clear in the recordings from IERRS network.
Moreover, the average of total residuals, which refer to within-event residual in the case of single event, and their standard deviations have been computed with 10 km interval of epicentral distances to increase the comprehensibility of the variation of residuals and also to compare with the within-event standard deviations of selected GMPEs. The paucity of data at a close distance makes it difficult for a reasonable judgment of standard deviation variation. However, one may clearly say that spectral accelerations at higher periods exhibit larger variability of residuals for both events. For M W 4.5 event, the largest variability of residuals is mostly seen in the data recorded in distant stations while the greatest standard deviations are detected in the closest distances for M W 5.7 for all ground motion parameters. These largest variabilities are also compared with the within-event standard deviations yielded by the regarded GMPEs. For that purpose, the ratios of the observed standard deviation of each bin to the within-event standard deviations of selected GMPEs are computed and the largest ratios for PGAs, PGVs, and SAs at four different periods are summarized in Table 5. The bold font styles in the table underline the larger values of residual standard deviations than within-standard deviations of the relevant GMPEs. Although observed variabilities only surpass the within-event standard deviations of SAs at higher periods for M W 4.5 event, the maximum standard deviations of calculated bin residuals, which are majorly resulted in close distances, are greater than the within-event standard deviations given by considered GMPEs for all ground motion parameters. The variation of total residuals is also investigated with respect to soil conditions (Fig. 17). Also, the mean ± standard deviations are seen with 200 m/s of V s30 intervals, as in the previous assessments with respect to epicentral distances. The ground motions recorded in stiff soil sites, which are majorly composed of IERRS data, yield relatively smaller residual fluctuations. It should be stated that considered stations of AFAD and IEEWS data are majorly deployed on soft and medium stiff soil. This could explain the reason that we cannot identify discern faster attenuation at great distances as clearly at AFAD stations while we can identify it in IERRS data.
In the scope of the study, the azimuthal variations are also investigated to understand whether it is responsible for the unique attenuation attributes of the region. The source-to-station azimuth angles, on the other hand, majorly range from 60 to 100 degrees (see Fig. 1 for the spatial distribution of stations) making it difficult to distinguish azimuthal alterations precisely. Another obstacle that makes discerning the azimuth impact harder may be the scarcity of nearfield recordings in our database. Nevertheless, the variability of path-corrected residuals, which computed by subtracting the mean of distance bins from total residuals, with azimuth has been explored for both events, but no definitive conclusion could be achieved as envisaged. For the sake of the brevity, the related plots are not presented herein. Nonetheless, Türker et al. (2022), scrutinize the spatial variability of strong ground motion in three seismic events, whose dataset includes the M W 5.7 and M W 4.5 Silivri earthquakes, using only AFAD stations that encircled the Sea of Marmara 360 degrees. Higher observed PGAs and SAs toward the eastern Marmara, Armutlu peninsula, and Kocaeli are attributed with the region's distinguishable azimuthal patterns. The study is state that the rupture propagation directivity of these medium-sized seismic events may be primarily responsible for the azimuthal variation in the region.

Anelastic attenuation of the S-waves
The frequency-dependent attenuation of the direct S-waves of these events is analysed by observing the spectral decays at selected central frequencies on acceleration Fourier amplitude spectra. To eliminate the source and site effects on the observed spectra, we used the coda normalization method introduced by Aki (1980), which uses coda waves to normalize the observed S-wave amplitudes. The frequency-dependent model is created assuming a power-law function, in the following form,   where Q 0 is the value of Q S at 1 Hz and n is the frequency parameter. Q S varies within the crust in regard to inelasticity or scattering, or both at the same time and it is inversely proportional to attenuation. We make use of all the recordings with adequate bandwidths, covering the selected frequency range for the designated model. The overall considered frequency range is between 0.5 and 21 Hz. Selected central frequencies and their corresponding ± ranges that form the bins are presented in Table 6. We utilized only the vertical components of the recordings and those having very high signal-to-noise ratios.
The coda waves were windowed for pre-defined lapse times, which is increasing as the travel time of S-waves increase. As the closest recording's lapse time is selected as 124 s (27 km), the furthest one is 160 s (101 km). When lapse times are greater than approximately twice the S-wave travel time, the spectral amplitudes of coda waves are considered as distance-independent, consisting of scattered S-waves due to crustal heterogeneities. It should also be noted that coda wave amplitude is proportional to the S-wave source spectral amplitude (Tsumura 1967).
To establish the attenuation model, first, the average amplitudes are calculated by taking the mean of the amplitudes in each predefined frequency bin, both for S-waves and coda waves. The calculated mean amplitudes for the given bins are denoted as A s (f , r) and A c f , t c for the S-waves and coda waves, respectively. The ratio of A s (f , r) and A c f , t c is taken and multiplied by a distance-dependent geometrical spreading model, G R hyp , expressed in the form of R − hyp . In our study, we used the γ exponents as used in Akinci et al. (2013), that is a tested model for the west of Turkey, The calculation process can be described as, where β is the average S-wave velocity. The estimated values are plotted against hypocentral distance, and linearly regressed as shown in Eq. (15). Q S is calculated from the estimated slope (Eq. 16). The regression analysis and the following Q S estimation are repeated for each central frequency. Therefore, we obtain a Q S value for each central frequency. Figure 18 shows four of such analyses and Q S estimations at 1, 3, 6 and 9 Hz, whereas the rest of the estimations and the functional Q S (f) model are presented in Table 6. Horasan et al. (1998) studied anelastic attenuation for the Marmara Sea (1.5-24 Hz), while Horasan and Boztepe-Güney (2004) showed how Q S estimations may differ in the same region when datasets are grouped for smaller subregions (1.5-12 Hz). However, both generalised models in Horasan et al. (1998) and in Horasan and Boztepe-Güney (2004) for the Marmara Fig. 19 Frequency-dependent Q S models prepared in this study, compared with past studies Sea show similar results. Their results indicate relatively high attenuation at lower frequencies. Akinci et al. (2006) and Akinci et al. (2013) studied parametrization of the ground motions, along with an anelastic attenuation models for the Marmara region (0.5-16 Hz) and for the west of Turkey (0.25-9 Hz), respectively. Their models show lower attenuation at low frequencies, compared to the previous two studies mentioned here. Our estimations for the lower frequencies follow a trend between these models, whereas we observe high attenuation at higher frequencies. The graphical comparison of these models is presented in Fig. 19.

Conclusion
In this study, we examined the seismological and engineering parameters of two earthquakes that occurred at an exceptional location in the Marmara Sea, where not many events were reported in the recent past. We had the opportunity to use the majority of the existing recordings associated with these events that allowed us to study several aspects of these events. The overall conclusions can be summarized as follows: • Estimated source parameters of the M W 4.5 event are in reasonable agreement with those from global models for moderate magnitude, strike-slip events. However, the corner frequency and the stress drop of the M W 5.7 event are relatively high considering the global estimations. • The time difference between hand-picked P-wave and S-wave onsets is modelled as 0.118 ×R hyp . It is correlated only with distance, as the regressed model intersects almost at zero. This model is slightly different than the globally used relation of R hyp ∕8 (Ancheta et al. 2013;Kishida et al. 2016). • The 5%-75% significant duration estimates for the vertical component have a clear linear trend with hypocentral distance, yielding 2.90(± 1.85) + 0.13(± 0.03)× R hyp and 2.74(± 2.08) + 0.12(± 0.03)× R hyp for M W 4.5 and M W 5.7 events, respectively. Besides, the path components of these models match well with the global models. • Maxima of peak ground motion parameters and SAs are registered at closer stations, which are located on the north Coast of Marmara Sea, on softer soils.
• The PGAs, PGVs and SAs of the M W 5.7 earthquake are more scattered when compared to those of the smaller event. As a result, the data spread between median ± 2σ lines of the three GMPEs employed. • Our results suggest that beyond 60 km from the earthquake source zone the peak and spectral ground motion values start to attenuate distinctively faster especially for IERRS data. This observation needs to be explored further for underlying reasons, which could be related to regional geometrical spreading or intrinsic attenuation properties. Investigation of a larger data set that includes earthquakes originating from the same and other NAF segments in the Marmara Sea for that matter will be the next step. • We analysed the anelastic attenuation for the northeast of Marmara, by using the coda normalization method. We observed considerable differences exist between our model and the models prepared in the past. This comes due to the differences in the data used in all these studies (i.e., differences in magnitudes, station locations, distances etc.). Overall, our model shows a good agreement with the models prepared in Horasan et al. (1998) and in Horasan and Boztepe-Güney (2004). However, compared to all other models in Fig. 19, the model prepared in this study show higher attenuation at high frequencies.