Articles | Volume 14, issue 10
Research article
25 Oct 2023
Research article |  | 25 Oct 2023

The 2022 MW 6.0 Gölyaka–Düzce earthquake: an example of a medium-sized earthquake in a fault zone early in its seismic cycle

Patricia Martínez-Garzón, Dirk Becker, Jorge Jara, Xiang Chen, Grzegorz Kwiatek, and Marco Bohnhoff

On 23 November 2022, a MW 6.0 earthquake occurred in the direct vicinity of the MW 7.1 Düzce earthquake that ruptured a portion of the North Anatolian Fault in 1999. The Mw 6.0 event was attributed to a small portion of the Karadere fault off the main North Anatolian Fault that did not rupture during the 1999 sequence. We analyze the spatiotemporal evolution of the MW 6.0 Gölyaka–Düzce seismic sequence at various scales and resolve the source properties of the mainshock. Modeling the decade-long evolution of the background seismicity of the Karadere fault employing an Epistemic Type Aftershock Sequence model shows that this fault was almost seismically inactive before 1999, while a progressive increase in seismic activity is observed from 2000 onwards. A newly generated high-resolution seismicity catalog from 1 month before the mainshock until 6 d after was created using artificial-intelligence-aided techniques and shows only a few events occurring within the rupture area within the previous month, no spatiotemporal localization process and a lack of immediate foreshocks preceding the rupture. The aftershock hypocenter distribution suggests the activation of both the Karadere fault, which ruptured in this earthquake, and the Düzce fault that ruptured in 1999. First results on the source parameters and the duration of the first P-wave pulse from the mainshock suggest that the mainshock propagated eastwards, which is in agreement with predictions from a bimaterial interface model. The MW 6.0 Gölyaka–Düzce event represents a good example of an earthquake rupture with damaging potential within a fault zone that is in a relatively early stage of the seismic cycle.

1 Introduction

Large strike–slip fault zones such as the San Andreas Fault in California, USA, or the North Anatolian Fault in Türkiye, among others, host some of the largest shallow earthquakes (typically up to M  8; see, e.g., Wesnousky, 1988; Stirling et al., 1996; Martínez-Garzón et al., 2015) worldwide. Some of these hazardous faults run near urban areas; hence, they have an associated risk. This is the case for the western portion of the North Anatolian Fault, which is a seismic threat for the Istanbul metropolitan area and nearby population centers. There, the return period of M > 7 earthquakes rupturing the main fault trace has been estimated to be approximately 250 years (Murru et al., 2016). The most recent large earthquakes along the North Anatolian Fault zone (NAFZ) were the 17 August 1999 MW 7.4 İzmit event and the 11 November 1999 MW 7.1 Düzce event (Fig. 1). Together, they ruptured approximately 180 km and connected the Marmara segment in the west to the 1944 rupture in the east (Bürgmann et al., 2002; Sengör, 2005; Bohnhoff et al., 2013). On 23 November 2022, a MW 6.0 earthquake occurred around 6 and 10 km away from the cities of Gölyaka and Düzce, respectively, and about 200 km eastward of the Istanbul metropolitan area. In the following, we refer to this event as the Gölyaka–Düzce earthquake, as it was felt especially in the province of Düzce and its districts (Eyidoğan, 2022). Hypocenter locations provided by KOERI (Kandilli Observatory and Earthquake Research Institute;, last access: 18 October 2023) and AFAD (Disaster and Emergency Management Presidency;, last access: 18 October 2023) reported that the Gölyaka–Düzce earthquake occurred close to the intersection of the ruptures of the 1999 MW 7.4 İzmit and MW 7.1 Düzce earthquake ruptures (Bouin et al., 2004; Konca et al., 2010; Fig. 1). Such MW 6 earthquakes are relatively infrequent in the region (according to the ISC-GEM – International Seismological Centre Global Instrumental Earthquake Catalogue;, last access: 18 October 2023; similar-sized events could have occurred in 1926 and 1944), but they hold potential to damage key infrastructure and lead to casualties. Hence, analyzing this earthquake and its pre- and post-seismic deformation is important to illuminate the local fault architecture (e.g., Ross et al., 2020) and recover how moderate events nucleate in the region (e.g., Malin et al., 2018; Durand et al., 2020). In addition, analysis of the earthquake source properties can help us to decipher any preferential directions of seismic energy release, which is essential for seismic hazard estimation.

Figure 1Seismotectonic setting. (a) Location map with a red rectangle indicating the area enlarged in panel (b). (b) Study area with the location of main population centers along the North Anatolian Fault zone (NAFZ). Colored lines denote rupture extents of historical earthquakes along NAFZ, with their respective magnitudes and dates indicated in the legend. Arrows are an updated GPS velocity field for Türkiye (Kurt et al., 2022), considering the Eurasia-fixed reference frame. The magenta star indicates the Gölyaka earthquake epicenter, along with its focal mechanism from KOERI, and the focal mechanisms of the İzmit (green) and Düzce (purple) earthquakes in 1999 (data from the global centroid moment tensor (CMT) catalog; Ekström et al., 2012; Dziewonski et al., 1981). (c) View of the region struck by the Gölyaka earthquake.

The Gölyaka–Düzce earthquake was the first M  6 in the area after the 1999 M > 7 earthquakes and their aftershock sequences. As such, it represents the occurrence of an earthquake with damaging potential in a fault zone that is broadly still in the early stage of the seismic cycle, which is when the accumulated elastic strain is relatively low. However, at this location, the geometry of the fault zone is highly complex, with slip partitioning along two main fault traces bounding the Almacık Block and numerous secondary faults being obliquely oriented (Fig. 1c). The moment tensor estimation by different agencies (Table S1 in the Supplement) consistently reported on a strike–slip mechanism with a small normal faulting component for this earthquake. The hypocenter location and fault plane appeared consistent with a portion of the Karadere fault (Fig. 1c). Different lines of evidence suggest that at least a few kilometers of the Karadere fault did not rupture during the 1999 İzmit–Düzce sequence and thus remained loaded and then ruptured in the Gölyaka–Düzce earthquake (Bohnhoff et al., 2016; Özalp and Kürcer, 2022). These include (i) the magnitude of the latest mainshock, (ii) the spatial extension of aftershocks and (iii) the previous surface mapping of local faults. First field surveys immediately after the MW 6.0 Gölyaka–Düzce earthquake found no surface rupture, indicating that the slip did not extend to the surface (Özalp and Kürcer, 2022).

In this study, we analyze the spatiotemporal evolution of the Gölyaka–Düzce seismic sequence and the preceding seismicity in the area in detail, with a focus on the source mechanisms and in context of the local seismotectonic setting. Our primary goals are to determine how the earthquake initiated, what the ongoing deformation mechanisms in the region are and how the energy from this mainshock was radiated. Insights from this earthquake sequence are important for us to learn about the dynamics of potentially damaging earthquakes in complex transform fault settings and the rupture of highly stressed faults in the immediate vicinity of a major strike–slip transform with relatively low-stress conditions corresponding to the fault zone being early in its seismic cycle.

Figure 2Regional seismicity. (a) AFAD catalog from 1990 to 30 November 2022. The yellow star denotes the MW 6.0 Gölyaka–Düzce epicenter. Blue boxes indicate the target regions in which the seismicity is analyzed. The column on the right shows the catalog's probability distribution function of earthquake magnitude (PDF), where the vertical dashed line denotes the MC. (b) Same as panel (a) but for the KOERI catalog. See Figs. S3 and S4 for the spatial distribution of seismicity inside each region.

2 Data and methodology applied

2.1 Background seismicity evolution

To put the Gölyaka–Düzce seismic sequence in a regional and long-term context, we established a consistent regional seismicity catalog between 1990 and 2022 from the two national seismic agencies. The AFAD catalog reported 31 081 events with magnitudes M (1.3–7.6) for such a period, whereas the KOERI catalog reported 42 050 events M (1.4–7.4). We analyzed the decade-long evolution of the AFAD and KOERI regional seismicity catalogs through a declustering process (Fig. 2). Both catalogs changed the reported magnitude type at the end of 2011 and beginning of 2012 from the duration magnitude Md to local ML. For consistency, we homogenized both catalogs, converting uniformly to moment magnitude, MW, following the empirical relationships proposed for the study region by Kadirioğlu and Kartal (2016):


After both catalogs were homogenized to MW, we estimated the magnitude of completeness MC of each catalog, following a probabilistic approach to fit the frequency–magnitude curve (Ogata and Katsura, 1993; Daniel et al., 2008; Jara et al., 2017). In contrast to the maximum curvature method, this technique fits a function to explain the magnitude–frequency distribution, including all catalog events. The number of earthquakes is fit as a function of magnitude as follows:

(3) N ( m ) = A × 10 - b m × q ( m ) ,

where the b value represents the slope of the Gutenberg–Richter law, A is a normalization constant, and q(m) is the probability that one earthquake of magnitude m is listed in the catalog. Then, we modeled q as follows (Ogata and Katsura, 1993):

(4) q ( m ) = 1 2 + 1 2 erf m - μ ^ 2 σ ^ ,

where erf is the error function, and μ^ and σ^ correspond to the mean and standard deviation of the probability distribution function, respectively. We optimized [A,b,μ^,σ^] for each catalog, following a Bayesian approach to derive the parameters' posterior probability density function (PDF). We used the Markov chain Monte Carlo sampler of PyMC (Salvatier et al., 2016) to draw 500 000 samples from the posterior PDF. The inferred parameters and their associated uncertainties are given in Table S2. Then, the completeness magnitude MC is computed as follows:

(5) M C = μ ^ + 2 σ ^ ,

i.e., with a 97.7 % probability threshold, which yielded a MC=3.4 for the AFAD catalog, whereas, for the KOERI one, we obtained a MC=4.1 (see the insets in Fig. 2 and the obtained fitting in Figs. S1 and S2). Once MC was estimated for both catalogs, we declustered them using an epidemic-type aftershock sequence model (Marsan et al., 2017; Jara et al., 2017). Such an approach considers the total seismicity rates λ(x,y,t) to be the following sum:

(6) λ ( x , y , t ) = μ ( x , y , t ) + ν ( x , y , t ) ,

where ν(x,y,t) accounts for the aftershock productivity, and μ(x,y,t) is the background seismicity rate for earthquakes occurring at a given location (x,y) and time t. The aftershock rate was estimated following the Omori–Utsu law modified by a power spatial density, as follows:

(7) ν i ( x , y , t ) = i t i < t κ ( m i ) ( t + c - t i ) p × ( γ - 1 ) L ( m i ) γ - 1 2 π ( x - x i ) 2 + ( y - y i ) 2 + L i 2 γ + 1 2 ,

where c, γ and p are constants, and κ(m) is the productivity law with a constant α (Ogata, 1988). L(m)=L0×100.5(m-MC) is the characteristic length in kilometers (Utsu and Seki, 1955; van der Elst and Shaw, 2015). Here, we imposed realistic values for the following parameters: α=2, p=1, c=10-3d, γ=2, and L0= 1.78 km (Marsan et al., 2007; Jara et al., 2017; Karabulut et al., 2022). Parameters κ and μ(x,y,t) were inverted. The background seismicity was computed as follows:

(8) μ ( x , y , t ) = i μ ( x i , y i , t i ) λ ( x i , y i , t i ) e - ( x - x i ) 2 + ( y - y i ) 2 l × e - t - t i τ × 1 2 π l 2 a i ,

with l and τ as space and time and being smoothing parameters, and ai=2τ-τ(ets-tiτ-ete-tiτ), where ts and te are the temporal beginning and end of the catalog, respectively. κ is inferred as

(9) κ = i 1 - μ ( x i , y i , t i ) λ ( x i , y i , t i ) i e α m i ( ln ( t e + c - t i ) - ln ( c ) ) .

Figure 3Picks and detections from the AI-derived catalog. (a) Associated picks as a function of time per station. The vertical blue line marks the MW 6.0 Gölyaka–Düzce earthquake. (b) Venn diagram showing the earthquakes included in our catalog of detections vs. the events included in the KOERI catalog for the same spatiotemporal region.


Here, we used a smoothing length l= 100 km and a smoothing duration τ= 100 d. Such choices are able to preserve the potential accelerations or decelerations from the catalogs (Marsan et al., 2017; Jara et al., 2017). We declustered the catalogs using the obtained MC for each catalog and the fault regions in Fig. 2. When doing so, we observed an apparent increase in the seismicity from the AFAD declustered catalog at around 2012 (Fig. S1). Around that time, AFAD changed the reported magnitudes from Md to ML. Although we converted the corresponding magnitudes to MW, this change in the magnitude estimation might still produce spurious acceleration or deceleration in the background seismicity rate. We then tested higher MC values, finding that such behavior disappears at around MC= 4.1 (Fig. S1). Thus, we finally used MC=4.1 for both catalogs (Fig. 3a and b). The final parameters utilized for each catalog are provided in Table S1.

2.2 Enhanced seismicity catalog framing the mainshock

To generate an optimized enhanced seismicity catalog with the lowest possible magnitude detection threshold around the MW 6.0 Gölyaka–Düzce earthquake epicenter, we processed continuous waveform recordings from 16 local seismic stations and 9 local accelerometer stations. We covered a time period from 1 month before the mainshock up to almost 6 d after it (23 October 2022 at 00:00 LT up to 29 November 2022 at 00:00 LT). The employed stations belong to the AFAD and KOERI seismometer and strong motion networks.

We detected P- and S-wave onset times embedded in the continuous recordings by applying the supervised artificial intelligence (AI) method of PhaseNet (Zhu and Beroza, 2019) that is trained on the seismicity database from northern California. This method has proven to improve the detection process, especially for small earthquakes (e.g., Martinez-Garzon et al., 2023). With this method, 148 948 body wave onsets were detected, of which 78 410 were detections of P waves and 70 568 were detections of S waves.

The P and S picks were associated with seismic events using the unsupervised technique of Gaussian Mixture Model Association (GaMMA; Zhu et al., 2022). To classify an event to be an earthquake, a minimum of four necessary picks (either P and/or S) was set. The picks were spatiotemporally clustered using the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) method. About 19 % of the total number of picks were associated with earthquakes. In this way, we have obtained a catalog of detections containing 3361 possible seismic events (Fig. 3a).

As a comparison, KOERI reported a total of 505 events with magnitudes ML 0.5 for the same region and time interval analyzed here (Fig. 3b). Out of these, 440 correspond to common events from both catalogs.

In the next step, the waveforms from all events corresponding to the period before the MW 6.0 mainshock were visually inspected. About 343 detections from the time period before the mainshock were removed, as they showed signals in only one or two of the accelerometers, typically exhibiting tS-tP>5s, which is larger than what is expected for a small local event. Additionally, 35 events were identified as being regional events with locations outside the study region, and an additional 11 events were identified as being duplicates and removed. In the following, we refer to this catalog as the “catalog of detections”.

We calculated event locations by employing the probabilistic location software NonLinLoc (Non-Linear Location; Lomax et al., 2009). Here, only events with a minimum of six P and/or S picks were further processed, which implicitly removes possible false signal associations with fewer than six phases from the catalog of detections. The local 1-D velocity model from Bulut et al. (2007) was employed, assuming a constant vp/vs ratio of 1.73. The search area encompassed a 400 km× 200 km region centered around the mainshock epicenter. In the following, we refer to this refined catalog as the “catalog of absolute locations”. Further details on the refining of the catalog of absolute locations are provided in Sect. S1 in the Supplement. In this way, we obtained a catalog of 1290 events with absolute locations that contain 8927 P-wave picks and 7822 S-wave picks for further processing (Fig. S5). In this catalog, the median errors in the x, y and z directions are 2.3, 3.1 and 3.4 km, respectively.

In the next step, a relative event relocation was performed using HypoDD (Waldhauser and Ellsworth, 2000; Waldhauser et al., 2004). We utilized both catalog differential travel times derived from the automatic PhaseNet P and S picks and the cross-correlation time differences derived from the event waveforms. To estimate the waveform cross-correlations, we employed time windows covering 1 and 2 s centered at the P and S onset, respectively. The waveforms were filtered with a third-order Butterworth bandpass filter between 2 and 10 Hz. The retrieved correlations with a normalized cross-correlation coefficient of at least 0.7 were kept, and the square of the cross-correlation coefficient was used as weight in the relocation procedure. To look at the spatiotemporal evolution of the seismicity, we demanded a minimum of eight catalog time differences (either P and/or S phases) for each event combination, resulting in a catalog of 918 relocated events. In the following, we refer to this further-refined subset as the “relocated catalog”. The median formal relative relocation errors in the x, y and z directions are 11, 13 and 12 m, respectively.

2.3 Earthquake source parameters and directivity

The estimation of point source parameters of the MW 6.0 mainshock was performed using the spectral fitting method (Kwiatek et al., 2011, 2015). We used 249 high-gain seismometers of Kandilli Observatory (KO) and the National Seismic Network of Türkiye (TU). The combination of these networks offers a sufficient azimuthal coverage and signal-to-noise ratio for the analyzed catalog. From these networks, stations with epicentral distances of 200–800 km are used to derive source parameters of the mainshock.

Three-component ground velocity waveforms from stations with a signal-to-noise ratio of at least 4 dB were filtered using a second-order 0.02 Hz high-pass Butterworth filter and then integrated in the time domain. We utilized window lengths of 25 s from the P- and S-wave ground displacements, allowing for an additional 4 s before the P- or S-wave onsets. The selected window length of 25 s, in combination with the minimum station distance of 200 km, prevents contamination of the P-wave window by S-wave energy, as the S–P time difference is larger than 25 s and assumes average velocities and vp/vs ratios along the travel path. The edges of the selected windows were smoothed using von Hann's taper. We utilized the three components from the far-field ground displacement spectra using the Fourier transform and then combined altogether. The observed ground displacement spectra were fit to Brune's point source model, which is expressed as

(10) u th ( f ) = R c 4 π ρ V c 3 R M 0 1 + f f c 2 exp - π f R V c Q c ,

where R is the source–receiver distance; M0 is the seismic moment; fc is the corner frequency, Qc is the quality factor, and Rc is the average radiation pattern correction coefficient, where c represents either the P or S wave. Following Boore and Boatwright (1984), we applied RP=0.65 and RS=0.7 for P and S waves, respectively, that are representative constants for the regional strike–slip faults. We used VP= 5680 m s−1 and VS= 3280 m s−1 (from Bulut et al., 2007; averaged at the depth interval where the earthquakes occurred), assuming VP/VS= 1.73 and a density ρ= 2700 kg m−3. We inverted for (M0, fc; Qc), optimizing the cost function

(11) | | log 10 u obs ( f ) - log 10 u th ( f ) | | L 1 = min ,

where uth(f) and uobs(f) are the theoretical and observed ground displacement amplitude spectrum for a given station and P or S wave. The starting model for M0 and fc was taken, using Snoke's integrals (Snoke, 1987), and we assumed initial values of Q=400 for both P- and S-wave trains. We utilized a grid search optimization (assuming a starting model), followed by simplex algorithms (starting from the best model of a grid search). Source parameters that deviated from the average by more than 3 standard deviations were eliminated from the calculation. The final source parameters (i.e., seismic moment, corner frequency and quality factor) were calculated as average values from all stations. An example of such calculated spectra for a station is provided in Fig. 4.

In the following, we calculated the static stress drop using the formula valid for a rectangular strike–slip fault (Shearer, 2019):

(12) Δ σ = 2 π M 0 W 2 L ,

where W= 8 km represents the fault width, which is assumed from the depth extent of the aftershocks (see Sect. 3.2.2). The rupture length L is estimated to be double the source radius, using Brune's source model constants and, for comparison, Haskell's rectangular source, assuming VR=0.9VS (see Savage, 1972; Table 1 for details) and using

(13) L = C C V c 2 π f c ,

in which CC is the geometrical correction coefficient (CP=CS=4.7 for Brune's model and CP=1.2 and CS=3.6 for Haskell's model; see Savage, 1972). fc and Vc represent the corner frequency and seismic velocity of either P or S waves, respectively.

For the earthquakes comprising the absolute locations catalog, we estimated only moment magnitudes MW using a simplified approach, following the Snoke (1987) integrals:


where u˙(f) and u(f) are ground velocity and displacement S-wave spectra corrected for attenuation and prepared from S-wave waveforms processed in the same way as for the mainshock. The original seismograms were filtered with 1 Hz high-pass second-order Butterworth filter, and we used a shorter 5 s window framing the first S wave arrival to limit the influence of low-frequency noise for predominantly small earthquakes (M < 4). The integrals in Eqs. (14) and (15) were corrected for the finite frequency band, following Di Bona and Rovelli (1988). The seismic moment has been estimated with the following (Snoke, 1987):

(16) M 0 = 8 π ρ V S 3 R K S 3 J S 0.25 ,

and the moment magnitude was calculated using the standard relation (Hanks and Kanamori, 1979).

(17) M W = ( log 10 M 0 - 9.1 ) 1.5

Similar to the mainshock, for each event, the final seismic moment and moment magnitude were calculated as the average values from all stations containing S-wave arrivals. Due to the limited number of S waves with a sufficient signal-to-noise ratio for the smallest earthquakes, the uncertainties were estimated using the mean absolute deviation.

Large earthquake ruptures potentially involve a propagation process along a fault plane. The rupture propagation direction could be deduced from the azimuthal variations in the amplitude and frequency content of the apparent source time functions (ASTFs) providing important information for seismic risk assessment. For a unilateral rupture, ideally this would lead to shorter ASTFs displaying larger amplitudes in the direction of rupture propagation and longer-duration and smaller-amplitude ASTFs in the opposite direction. To obtain the ASTFs, we initially tested the application of empirical Green's function (EGF) technique and tried four EGF candidates (Table S3) to recover the directivity of the mainshock (see Sect. S2 for all details) that led to inconclusive results.

Figure 4Source parameter analysis. (a) Station distribution employed for the source parameter estimation (upward-pointing red triangles). These stations lie within a source station distance between 200 and 800 km. The yellow star shows the 2022 MW 6.0 Gölyaka–Düzce mainshock. (b) Three-component displacement waveforms for the mainshock recorded by station AKO (with epicentral distance  530 km; see the white triangle in panel a). The red rectangle highlights the employed P-wave window. (c) Ground displacement spectrum of the P-wave signal (red) and the noise before the signal (black). The thick blue line indicates the modeled spectrum yielding the source parameters of M0=2.62×1018,fc=0.19Hz,QP=308.

We therefore tested the azimuthal variations in the duration and frequency content of the initial P wave arrivals for seismometers located at epicentral distances of 50–100 km from the mainshock. For comparison, we additionally included integrated signals from accelerometers located at much closer distances. We only used the unclipped first P-wave pulses that were rotated into the radial direction from three-component seismograms to enhance the signal-to-noise ratio of the initial portions of the P wave. The first P-wave pulses contain a combination of information, including the source time function and effects related to the wave propagation. However, comparing P-wave pulse characteristics for stations located at similar distances from the mainshock epicenter allows us to suppress the propagation effects. Therefore, the initial portion of the seismogram can be taken as a proxy for the ASTFs. The variation in the rise time and duration of the P pulses can then be used to infer whether the earthquake displays rupture directivity.

Figure 5Temporal evolution of the background seismicity at different segments of the NAFZ. (a) Complete (red) and declustered (black) AFAD catalog using MC=4.1. (b) Same as panel (a) but with KOERI catalog. (c) Cumulative background seismicity, color-coded by region. (d) Same as panel (c) but for the KOERI catalog. The vertical lines represent the time of occurrence of the MW 7.4 İzmit, MW 7.1 Düzce and MW 6.9 Saros earthquakes. (e) Normalized cumulative background seismicity for the AFAD catalog. (d) Same as panel (c) but for the KOERI catalog.


3 Results

3.1 Long-term evolution of background seismicity in the Gölyaka–Düzce region

We analyzed the evolution of the background seismicity along defined segments of the NAFZ, including the Marmara, İzmit, Düzce, Bolu and Karadere segments (Fig. 2). Both of the national Turkish catalogs introduced above show that the Bolu segment displays a low background seismicity rate when compared to, for example, the İzmit or Marmara segments (Fig. 5). Aseismic slip (surface creep) has been reported to occur along this segment and has occurred for at least 70 years (Ambraseys, 1970; Cakir et al., 2005; Cetin et al., 2014; Bilham et al., 2016; Jolivet et al., 2023). This might be a possible explanation for the low seismicity rate. The Marmara, İzmit and Düzce segments appear to host a constant background seismicity rate with time, especially after the 1999 İzmit and Düzce sequence (Fig. 5). Both catalogs report a deceleration of background seismicity after the 2014 MW 6.9 Saros earthquake (Bulut et al., 2018), supporting the idea that some significant deformation process not yet understood in detail was affecting the seismicity along the NAFZ (Karabulut et al., 2022).

The Karadere fault hosted a comparatively low background seismicity before the 1999 İzmit and Düzce earthquake sequence. A change in its seismic behavior is observed afterwards, which is when this segment experienced an increase in the seismic activity that included more than five events with M>MC= 4.1. The shape of the background rates is different for the AFAD and KOERI catalogs. This difference might be due to the different number of seismic stations operated by the agencies in this area; hence, this affects the monitoring capabilities and detection thresholds. Therefore, it is likely that the region was tectonically activated by the earthquake sequence in 1999 and progressively loaded since then, leading to the MW 6.0 Gölyaka–Düzce earthquake 23 years later. Interestingly, the region did not exhibit a lower background rate after the 2014 MW 6.9 Saros earthquake, which is different to the other NAFZ segments in the area.

3.2 Spatiotemporal seismicity distribution before and after the 2022 Gölyaka–Düzce earthquake

We obtained an enhanced seismicity catalog with 1290 refined hypocenter locations, as described in Sect. 2.2, covering the area of latitude (40–42 N) and longitude (30–32 E) for the time period 23 October 2022 at 00:00 LT up to 29 November 2022 at 00:00 LT and including moment magnitudes as low as MW 0.7. Out of these, a total of 222 and 1032 seismic events correspond to events preceding and following the 2022 Gölyaka–Düzce mainshock, respectively. For the same region and time interval, the seismicity catalog provided by the KOERI agency contained 529 events, out of which 23 and 506 corresponded to events preceding and following the mainshock, respectively (Fig. 3b).

Using a goodness-of-fit method (Wiemer and Wyss, 2000), the magnitude of completeness of the derived catalog within selected region is MWC=1.5. Calculating the b value for events above MC, from both the periods before and after the mainshock, we find a value of b=0.95±0.05 (Fig. S6). This could be related to the fact that we utilized MW, while many other estimates use ML, possibly leading to a larger b value (see, e.g., Raub et al., 2017). Alternatively, the relatively low b value may suggest that the fault did not yet release all of its accumulated strain (e.g., Gulia and Wiemer, 2019). Given the magnitude of the mainshock and the spatial extent of the rupture, we consider the latter option to be rather unlikely.

3.2.1 Seismic activity preceding the Gölyaka–Düzce earthquake

The MW 6.0 Gölyaka–Düzce earthquake hypocenter is located at the northeastern portion of the Karadere fault that remained unbroken during the 1999 MW 7.4 İzmit and 1999 MW 7.1 Düzce events.

Figure 6Seismicity located during the preceding month. Seismicity distribution included in the absolute location catalog (colored circles) during the month preceding the Gölyaka–Düzce earthquake (red star). Symbol size is encoded with magnitude. Surface ruptures of the 1999 MW 7.4 İzmit and MW 7.1 Düzce earthquakes are shown with dashed green and pink lines, respectively. For comparison, the seismic activity for three different time periods around the 1999 İzmit and Düzce mainshock is shown in cyan (from Bohnhoff et al., 2016; Bulut et al., 2007). Fault traces are from Emre et al. (2018).

The area that ruptured in the MW 6.0 Gölyaka–Düzce earthquake and its surroundings only displayed a small number of seismic events during the 30 d preceding the mainshock. The catalog of absolute locations reported 222 seismic events during this time, out of which 55 could be successfully relocated. Most of the relocated seismic activity occurred away from the future MW 6.0 earthquake rupture, extending up to 50 km to the east (Fig. 6). The locations of these seismic events show a good correspondence with the mapped local faults (Emre et al., 2018). A small cluster of events is visible at the eastern edge of the analyzed region, coinciding with the termination of a local fault, near a quarry area (see Fig. 6a for the location). The presence of a quarry in the area suggests that some of these events could be quarry blasts. However, these events appear to be regular seismic events based on the following: (i) these detections display regular P- and S-wave trains, (ii) their hypocentral depth is deeper than 8 km and (iii) these events occur randomly in time. Within a 25 km radius from the epicenter of the MW 6.0 Gölyaka–Düzce earthquake, 23 events were included in the catalog of absolute locations. The most active time period was between 6 and 11 November, where a small spatially clustered seismic sequence with magnitudes up to MW 2.2 occurred about 7.5 km to the north of the mainshock epicenter (Figs. 6 and S7). The location of this cluster coincides with the deepest part of the fault activated during the aftershock sequences.

Figure 7Temporal evolution of seismic activity and seismic moment. (a) Bars show the histogram of seismicity rates, where every bar represents a time period of 12 h. Gray and blue colors represent the seismicity included in the catalog of detections and absolute locations, respectively. Lines show the cumulative number of seismic events as a function of time. Lighter and darker colors represent the time periods before and after the mainshock. (b) Evolution of cumulative seismic moment release from the catalog of absolute locations.


Both the catalog of detections and the catalog of absolute locations show that seismicity rates were time-invariant, with a transient increase in seismic activity around 10 November, reflecting the transient cluster north of the future mainshock. This increase in the seismicity rates quickly decayed back to the level before the occurrence of the cluster and remained constant until the occurrence of the mainshock (Fig. 7). The regional seismicity did not display any significant acceleration at the scale of days to hours before the mainshock.

Figure 8Seismicity distribution after the Gölyaka–Düzce earthquake (colored dots). Cyan dots in the background reflect 1999 İzmit and Düzce aftershocks (from Bohnhoff et al., 2016; Bulut et al., 2007).

Figure 9Magnification of the spatiotemporal distribution of the seismicity during 6 d following the Gölyaka–Düzce earthquake. (a) Map view. Depth profiles along (b) A-A' (approximately perpendicular to the Karadere fault strike), and (c) B-B' (approximately perpendicular to the strike of the Düzce fault). The symbol size and color are encoded with magnitude and date, respectively.

3.2.2 The aftershock sequence following Gölyaka–Düzce earthquake

After the MW 6.0 Gölyaka–Düzce earthquake, vigorous seismic activity struck the region during the following days. Compared to the scattered seismicity in a much larger region, most of this early aftershock activity occurs within an area extending 15 km to the east and west, as well as 8 km to the north and south, of the mainshock epicenter, respectively (Figs. 8 and 9). Generally, aftershocks typically occur around the mainshock rupture area, and they may also activate nearby faults, due to stress changes induced by the mainshock. In a first-order approximation, the relocated aftershock activity delineates a planar structure trending SW–NE that is dipping towards the NNW, which is consistent with the geometry of the Karadere fault (Fig. 9). The plane best fitting to the seismicity (contained within 1 km distance) has a strike of φ=257 and a shallow dip of approximately δ=45 (Fig. 9). The strike of this plane is thus in good agreement with the moment tensor solutions for the MW 6.0 mainshock (Table S2). However, the dip of our plane is shallower than the δfm= 72–82 reported by the moment tensor solutions (Fig. 9c). The depth of the seismicity along the strike of the fault segment is not uniform, with the southwestern portion of the fault generally displaying shallower seismicity from 5 to 13 km depth and the northeastern portion of the activated fault between 9 and 16 km depth (Fig. 9b). Along the strike, the hypocentral location of the mainshock coincides with this depth change, suggesting the presence of a fault jog or a heterogeneity that could have promoted a stress concentration.

The mainshock triggered an aftershock sequence that within the first 6 d can be fitted with an Omori law of the shape N(t)=kt-p, with p=0.90 and k=2.5 (see Fig. S8). Typical values observed for the p value representing the decay rate oscillate around 1.0, suggesting that the aftershock decay associated with this sequence is fairly standard, including three to four MW > 4 earthquakes occurring within the analyzed time window.

3.3 Source parameters and directivity of the Gölyaka–Düzce mainshock

Earthquake source parameters for the 2022 MW 6.0 Gölyaka–Düzce mainshock are provided in Table 1, with the average values and multiplicative error factors calculated in log10 domain (García-García et al., 2004). The averaged seismic moment is 8.80 × 1017, leading to a moment magnitude of MW 5.9, which is equal to the moment magnitude given by AFAD. The average corner frequency fc values obtained for P and S waves are 0.23 Hz and 0.24 Hz, respectively, with a ratio of fcPfcS=0.96. The obtained ratio of corner frequencies from P and S waves is lower than the VPVS=1.73, which holds for a stationary source and can be decreased due to the rupture propagation effects (Sato and Hirasawa, 1973; Kwiatek and Ben-Zion, 2013). In general, the fcP>fcS arises for roughly equidimensional source models (length equal to width). For long and thin faults, lower fcPfcS ratios are to be expected; for example, fcPfcS=0.77, assuming rupture velocity VR=0.9VS (Savage, 1972). Nearly equal fcP and fcS values are given in a dislocation model with a unilateral rupture propagation (Haskell, 1964). The small fcP/fcS ratio might imply that the fault width W could be overestimated from the aftershock distribution and could be smaller than 8 km, which is also supported by the narrower fault width estimated from the early aftershock distribution (Fig. 9).

Attenuation greatly affects the amplitude and frequencies included in the seismograms. Commonly, S waves tend to have larger attenuation than P waves. The ratio between the P and S quality factors is

(18) Q P Q S = 3 4 V P 2 V S 2 .

For a Poisson solid, VP=3VS, thus resulting in QPQS= 2.25. Our observations provide a considerably lower QPQS= 1.2. Such lower ratios are not uncommon and have been interpreted as being the attenuating effect of pore fluids (Olsen et al., 2003; Hauksson and Shearer, 2006; Kwiatek and Ben-Zion, 2013; Kwiatek et al., 2015).

Utilizing the average source size and seismic moment from both P and S waves, the static stress drop of the mainshock is estimated as 0.61 and 1.48 MPa, while using a Brune model (Eq. 12) and a Haskell model, respectively (Table 1). The estimated rupture length varies at around 14 and 6 km for the Brune and Haskell models, which yields a rectangular source with a small L/W ratio. A relatively small aspect ratio was also observed for the MW 7.1 in 1999 in the direct vicinity of this area (Bürgmann et al., 2002).

Table 1Source parameters for the 2022 MW 6.0 Gölyaka–Düzce earthquake. fc is the corner frequency. M0 is the seismic moment. L is the source rupture length. Q is the quality factor. Δσ is the stress drop.

Download Print Version | Download XLSX

Figure 10First-peak duration recorded at seismic stations between 50 and 100 km from the mainshock epicenter. (a) Station distribution near the epicenter (yellow star). Colored triangles highlight the stations used in panel (b). (b) Normalized ground displacement recordings on radial components. The waveforms are aligned relative to the P-wave arrival (0 s) in the time axis and are ordered according to the azimuthal angles relative to the mainshock. The time duration of the colored segments is shown color-coded for in the station symbols in panel (a).

Figure 10 shows P-wave arrivals highlighting the initial portion of the ground displacement record Δt. Longer Δt rise times and durations of first P-wave displacement pulses are observed for western stations with azimuth angles of 196–293 (i.e., stations SUSU, GEYV, KAYN and KAND). At the same time, eastern stations at comparable distances and azimuth angles ranging from 32 to 130 display shorter rise times and visibly higher-frequency content (see discussion in Douglas et al., 1988; Fig. 10b), especially for stations RUZG and BCAM near 90 azimuth. These observations suggest an eastward rupture propagation, while assuming a unilateral rupture. However, in the case of a more complicated rupture process, the shorter rise time could also be promoted by a closer large local slip asperity in the easterly direction. We also estimated the azimuthal variations of fc for the stations between 200 km and 800 km from the mainshock (Fig. S9). Larger fc values are observed at approximately 100, hence being roughly consistent with the eastward rupture propagation. However, we note that scattered large fc values were also observed at other azimuths.

4 Discussion

The various spatiotemporal scales covered by the different methodologies applied in this study provide insights into the processes leading to and involved in the rupture of the MW 6.0 Gölyaka–Düzce earthquake. In the following, we discuss the most important patterns that emerged from the obtained results and their relation with the rheology of the region, the development of previous large earthquakes (i.e., the M > 7 İzmit and Düzce earthquakes) and its stage in the seismic cycle.

4.1 The 1999 M > 7 İzmit and Düzce earthquakes promoted the seismic activation of the Karadere fault

The Karadere fault connects the Akyazı and Düzce basins, which are both pull-apart structures in response to the regional transtensional tectonic setting (Pucci et al., 2006; Ickrath et al., 2015; Bohnhoff et al., 2016). The spatiotemporal evolution of seismicity along different portions of the broader Marmara region since 1990 has shown that the Karadere fault was primarily quiet until the occurrence of the 1999 M> 7 İzmit and Düzce events (Fig. 5). Most of the Karadere fault was activated during the 17 August 1999 MW 7.4 İzmit earthquake, while its northeastern portion likely hosted fewer aftershocks (Bohnhoff et al., 2016). The 1999 İzmit rupture was then extended further eastwards 87 d later, with the 11 November 1999 MW 7.1 Düzce earthquake onto the east–west trending Düzce fault splaying off the Karadere fault, and also dipping towards the north, with a dip of around 55 (Bürgmann et al., 2002).

The 23 November 2022 MW 6.0 Gölyaka–Düzce rupture likely occurred on the northeastern portion of the Karadere fault that remained inactive in 1999, marking the western flank of the Düzce Basin as a topographic depression north of the Düzce fault as a releasing bend. The fact that the İzmit rupture stopped on the Karadere fault and redirected onto the Düzce fault indicates that the northeastern Karadere fault acted as a barrier in 1999. This is supported by the observation of a lower seismic velocity contrast in the Karadere fault with respect to the fault regions west of it (e.g., the Mudurnu fault; see Najdahmadi et al., 2016). Nevertheless, our results show increased background seismic activity from 1999 onwards in the Karadere segment, with a visible increase in 2004–2005. One hypothesis is that the stress redistribution from the 1999 İzmit and Düzce earthquake sequences brought the Karadere segment closer to failure by stress transfer, leading to a progressive activation of this segment over the years. In that way, after 23 years of additional continuous tectonic loading, it was finally activated with a MW 6 event within a region of the fault zone that is still in a relatively early phase of the seismic cycle. Some support for this scenario comes from a reported change in stress regime, together with a rotation of the Shmin orientation in the Karadere segment before and after the 1999 İzmit and Düzce sequences (Ickrath et al., 2015). Before the earthquakes, a predominantly normal faulting stress regime was observed, while the strike–slip regime was observed after the Düzce earthquake. As the magnitude of SV at a certain depth is mostly given by the weight of the overburden, it is expected to remain approximately constant during the earthquake cycle. This suggests that the horizontal shear stresses on the fault increased after the 1999 sequence. We additionally note that the average recurrence period of M > 7 earthquakes in the area is around 250 years (Murru et al., 2016). Therefore, the recurrence time of a M > 6 earthquake should be about 25 years, which roughly fits with the occurrence of the last M > 7 earthquakes 23 years before the Gölyaka–Düzce event.

The observed changes in the background seismicity rates could also be related to a change in the seismic coupling of the region (e.g., Marsan et al., 2017; Jara et al., 2017). In particular, the occurrence of the 1999 M > 7 İzmit and Düzce earthquakes and their post-seismic deformation could have resulted in the promotion of the occurrence of an aseismic slip at depth, hence leading to a progressive decoupling of the fault. The buildup of stresses from the occurrence of enhanced aseismic slip can increase the background seismicity rates over the region with distributed deformation over a large area. Indeed, an additional proposed mechanism for the 1999 Düzce rupture was a viscoelastic post-seismic relaxation at depth, which affected a broad area from the 1999 İzmit rupture (e.g., Bürgmann et al., 2002; Ergintav et al., 2009). A detailed study on the microseismicity from this area also suggested that this possibility could account for the larger seismicity rates at depth (Beaucé et al., 2022).

4.2 How did the mainshock start?

Our catalog of absolute locations revealed at least 23 seismic events with an epicentral location less than 25 km from the MW 6.0 Gölyaka–Düzce event during the month before its occurrence. Out of these, only two are located in the northeastern segment of the Karadere fault, as the main fault segment that ruptured in the Gölyaka–Düzce event. The spatiotemporal evolution of these events does not suggest clustering but rather a scattered activation of the area (Fig. S7).

Likewise, the foreshocks do not generally resemble a spatial or temporal localization of the seismicity prior to the mainshock. This is of relevance, since a number of moderate to large earthquakes in this region displayed systematic foreshock activity (Bouchon et al., 2011; Ellsworth and Bulut, 2018; Malin et al., 2018; Durand et al., 2020). A similar spatiotemporally scattered precursory activity pattern as seen for the mainshock was also found for the 1999 MW 7.1 Düzce earthquake, where the largest event in the region of the earthquake rupture in the preceding 65 h was a M 2.6 event (Wu et al., 2013). Additional small events detected around the future Düzce 1999 rupture did not show any clear signatures of acceleration. The few seismic events preceding the 2022 MW 6.0 event, together with their lack of spatiotemporal localization, suggest the existence of relatively homogenous local stress conditions along this fault segment or, alternatively, homogeneous fault strength that would allow a progressive fault loading without rupturing small heterogeneities in the medium reflecting foreshock activity. In this respect, laboratory rock deformation experiments have shown that seismic precursors are more frequent on rough fault surfaces, while seismic foreshocks are much less frequent on polished fault surfaces (Dresen et al., 2020). This is consistent with the linear and relatively simple geometry of the eastern portion of the Karadere segment. In fact, the decade-long seismicity along the Karadere fault shows that it is notoriously more localized within the fault trace than in other fault areas (see e.g., Wu et al., 2013).

The fault area that was activated in the 1999 M > 7 İzmit and Düzce earthquakes is documented to continue displaying post-seismic deformation almost 20 years after (Ergintav et al., 2009; Aslan et al., 2019), which is mainly related to afterslip and viscoelastic relaxation. In this respect, one possibility is that the initiation of the mainshock was also promoted by the occurrence of distributed aseismic slip in the region at depth over a broad area (e.g., Beaucé et al., 2022; Karabulut et al., 2022). This is supported by the observation of a small number of seismic events around 11 November, at the bottom of the Düzce fault, near the place where the 1999 Düzce earthquake nucleated (Fig. 6). Another hypothesis is that a regional or local stress perturbation could have destabilized the northeastern Karadere fault that was close to failure. Some examples for such a potential stress perturbation may include tidal effects or seasonal effects, such as the effect of precipitation (e.g., Hainzl et al., 2013) or barometric pressure changes (Martínez-Garzón et al., 2023). Regarding seasonal or semi-periodic stress perturbations, it is worth mentioning that the MW 6.0 Gölyaka earthquake, a MW 4.9 event in 2021, as the largest and most recent event in this area, and the 1999 MW 7.1 occurred within the second-half of November. Further statistical analysis is not conducted in the frame of this study but may give a further indication with respect to whether earthquakes in this region show any significant temporal pattern.

4.3 Fault segments potentially activated during the mainshock and aftershock sequence

Based on the estimated rupture length from the mainshock source parameters, the event activated a  12 km long segment of the Karadere fault, terminating just east of the Düzce Basin (Fig. 9). Although we tested the application of EGF methods to recover the directivity more accurately, the analysis did not yield clear results (see Sect. S1 for details). The reasons for this are not clear. It may be that the events used for the EGF deconvolution did not fulfill all necessary criteria (e.g., occurring on the same location, with a similar focal mechanism and at least a unit of magnitude difference). Alternatively, it could be that the mainshock rupture did not activate a single fault segment, resulting in some complexities obscuring the directivity pattern.

The rupture complexity is also somewhat consistent with the spatial distribution of aftershock seismicity, which shows a heterogeneous event distribution, possibly also illuminating fault structures that were not directly involved in the mainshock rupture. On the western section, the spatial distribution of aftershock seismicity is oriented SW–NE, and part of the distribution suggests the activation of a NW-dipping fault plane of the mainshock, in accordance with fault–plane solutions of the event (Table S2) and with the size of the mainshock rupture estimated from the source parameters.

However, the eastern part of the aftershock distribution is also compatible with the fault geometry of the main Düzce fault activated in 1999. Indeed, the main cluster of events is located at approximately 10 km distance from the mapped surface trace of the Düzce fault. As the deepest aftershock seismicity is located at about 15 km depth, the distribution is also consistent with a fault dipping at about 55, as we previously reported in Sect. 3.2.2 (Fig. 9). Indeed, this dip is more consistent with the fault geometry reported for the Düzce fault (Bürgmann et al., 2002) than with the dip of the Karadere fault extracted from the focal mechanism of the MW 6.0 Gölyaka–Düzce earthquake.

Therefore, we suggest that the aftershock distribution that we obtained is likely reflecting the activation of both faults, with the Karadere segment displaying a steeper dip towards the northwest, as observed from the focal mechanisms, and the main Düzce fault at depth dipping more gently (around 55) towards the north.

4.4 A proxy for rupture directivity suggests a larger radiation of energy towards the east

Higher-frequency P-wave pulses with shorter rise times were identified at seismic stations found eastward of the rupture, suggesting that the mainshock propagated towards the east. This is consistent with the rupture directivity derived for this earthquake from the analysis of ground motions (Türker et al., 2023). A statistically preferred rupture propagation towards the east was also resolved in the Karadere fault segment below 5 km depth, based on the analysis of fault zone head waves (FZHWs) and fault zone reflected waves (FZRWs; Najdahmadi et al., 2016). At depth, the authors identified the faster side as being the elevated crustal Almacık block to the SE. Together with models of bimaterial ruptures, these results suggest that earthquakes on the Karadere segment nucleating at > 5 km depth have a physically explainable preferred propagation direction to the east. However, at a shallower depth, the fault core was detected to host even slower material between both blocks to either side (Najdahmadi et al., 2016). This led the authors to decide on a narrow wedge-shaped structure of the fault rather than a simple first-order impedance contrast of the fault. A preferred rupture propagation towards the east was also resolved in the Mudurnu fault segment (about 70 km west of the mainshock epicenter; see Fig. 6) from detection of fault zone head waves (Bulut et al., 2012). From the moveout of the fault zone head waves, a velocity contrast of about 6 % was estimated, with a slower seismic velocity for the northern side of the fault. An eastward propagation of the rupture was also reported for the 1999 MW 7.1 Düzce earthquake rupture from the joint analysis of geodetic, seismic and strong motion data (Konca et al., 2010). We conclude that based on our observations of an eastward-directed rupture during the MW 6.0 Gölyaka–Düzce earthquake, the observations of fault zone head waves in the region and the existence of bimaterial faults in the area should be considered to be an important ingredient for refined seismic risk studies in the area, especially for the Istanbul metropolitan region further to the west. However, this only applies if earthquakes are located along the bimaterial interface. For earthquakes located on secondary splay faults or in the rock volume, their rupture directivities may be related to other factors. A future possible analysis of the source parameters from the smaller events of the sequence may reveal whether the eastward directivity is a persistent feature in the region.

5 Conclusions

We investigated the source parameters of the 2022 MW 6.0 Gölyaka earthquake in NW Türkiye and the evolution of the seismicity framing this mainshock at various spatial and temporal scales. This earthquake mainly ruptured the Karadere fault, a small fault segment located in the direct vicinity of the 1999 MW 7.1 Düzce earthquake. Hence, this case is an example of a medium-sized earthquake which ruptured a critically stressed fault embedded in a fault zone that is overall in a relatively early stage of the seismic cycle. Our primary goal was to determine how the earthquake initiated, what the ongoing deformation mechanisms in the region are and how the energy from this mainshock was radiated. The main conclusions extracted from our analysis are the following.

The decade-long evolution of background seismicity in the Karadere segment shows that the segment was mostly silent before the 1999 M > 7 İzmit and Düzce earthquakes. From the year 2000 to the present, the segment has been comparatively more seismically active, supporting the hypothesis of a progressive approach to critical stress level of the fault segment.

The high-resolution seismicity catalogs derived in this study report on 23 locatable events during the previous month within a 25 km radius from the 2022 Gölyaka–Düzce earthquake. Only few of them occurred close to the future earthquake rupture, suggesting relatively homogenous fault stress conditions, and no signatures of foreshock localization were observed.

The early aftershocks of the sequence (i.e., first 6 d) suggested the activation of the Karadere fault segment dipping steeply towards the NW, as reported by the moment tensor, and the Düzce fault in the southern part dipping shallower directly towards the north. This suggests that the mainshock rupture, located along the Karadere fault, was able to trigger abundant aftershocks in the neighboring fault segment.

An analysis of the mainshock rupture directivity patterns, including an attempt to employ empirical Green's function analysis, did not yield clear results. However, a shorter rise time and the higher-frequency content of the P-wave pulses is observed at seismic stations located east of the mainshock hypocenter. If the mainshock rupture did indeed show promoted directivity towards the east, then the observation is consistent with predictions from models of bimaterial interfaces and observations from fault zone head waves at this fault.

Code and data availability

The here-derived enhanced seismicity catalog with absolute locations derived with NonLinLoc and moment magnitudes estimated is provided in the supplementary text file “catalog_absolute_locations.txt”.

The seismicity catalog from AFAD can be accessed through the link provided in the reference Turkish Disaster Emergency Authority (, AFAD, 2023). The seismicity catalog from KOERI can be accessed through the link provided in the reference KOERI observatory (, KOERI Observatory, 2023).

The here-generated AFAD and KOERI catalogs correspond to the time period from 23 October 2022 at 00:00 LT up to 29 November 2022 at 00:00 LT. Latitude and longitude ranges of 40–41, 30–32, respectively, are used.


The supplement related to this article is available online at:

Author contributions

PMG and DB generated the enhanced seismicity catalog. JJ analyzed the long-term seismicity rate variations. XC estimated the stress drop, and EGF conducted the directivity estimations, using methodologies developed previously by GK. All authors discussed the results in terms of the seismotectonic setting. PMG drafted the paper, and all co-authors contributed to its finalization.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


We thank Elif Turker and Fabrice Cotton for discussions and comparison with their directivity results from ground motion.

Financial support

This research has been supported by the Helmholtz Association (grant no. VH-NG-1232; SAIDAN). Xiang Chen received funding from the Deutsche Forschungs Gemeinschaft (DFG) in the frame of the proposal “Earthquake source characterization and directivity effects near Istanbul: Implications for seismic hazard”.

The article processing charges for this open-access publication were covered by the Helmholtz Centre Potsdam – GFZ German Research Centre for Geosciences.

Review statement

This paper was edited by Caroline Beghein and reviewed by two anonymous referees.


AFAD – Turkish Disaster Emergency Authority: Seismicity catalog, (last access: 20 October 2023), 2023. 

Ambraseys, N. N.: Some characteristic features of the Anatolian fault zone, Tectonophysics, 9, 143–165,, 1970. 

Aslan, G., Lasserre, C., Cakir, Z., Ergintav, S., Özarpaci, S., Dogan, U., Bilham, R., and Renard, F.: Shallow Creep Along the 1999 Izmit Earthquake Rupture (Turkey) From GPS and High Temporal Resolution Interferometric Synthetic Aperture Radar Data (2011–2017), J. Geophys. Res.-Sol. Ea., 124, 2218–2236,, 2019. 

Beaucé, E., van der Hilst, R. D., and Campillo, M.: Microseismic Constraints on the Mechanical State of the North Anatolian Fault Zone 13 Years After the 1999 M7.4 Izmit Earthquake, J. Geophys. Res.-Sol. Ea., 127, e2022JB024416,, 2022. 

Bilham, R., Ozener, H., Mencin, D., Dogru, A., Ergintav, S., Cakir, Z., Aytun, A., Aktug, B., Yilmaz, O., Johnson, W., and Mattioli, G.: Surface creep on the North Anatolian Fault at Ismetpasa, Turkey, 1944–2016, J. Geophys. Res.-Sol. Ea., 121, 7409–7431,, 2016. 

Bohnhoff, M., Bulut, F., Dresen, G., Malin, P. E., Eken, T., and Aktar, M.: An earthquake gap south of Istanbul, Nat. Commun., 4, 1999,, 2013. 

Bohnhoff, M., Ickrath, M., and Dresen, G.: Seismicity distribution in conjunction with spatiotemporal variations of coseismic slip and postseismic creep along the combined 1999 Izmit-Düzce rupture, Tectonophysics, 686, 132–145,, 2016. 

Boore, D. M. and Boatwright, J.: Average body-wave correction coefficients, B. Seismol. Soc. Am., 74, 1615–1621, 1984. 

Bouchon, M., Karabulut, H., Aktar, M., Özalaybey, S., Schmittbuhl, J., and Bouin, M.-P.: Extended Nucleation of the 1999 Mw 7.6 Izmit Earthquake, Science, 331, 877–880,, 2011. 

Bouin, M. P., Bouchon, M., Karabulut, H., and Aktar, M.: Rupture process of the 1999 November 12 Düzce (Turkey) earthquake deduced from strong motion and Global Positioning System measurements, Geophys. J. Int., 159, 207–211,, 2004.  

Bulut, F., Bohnhoff, M., Aktar, M., and Dresen, G.: Characterization of aftershock-fault plane orientations of the 1999 İzmit (Turkey) earthquake using high-resolution aftershock locations, Geophys. Res. Lett., 34, L20306,, 2007. 

Bulut, F., Ben-Zion, Y., and Bohnhoff, M.: Evidence for a bimaterial interface along the Mudurnu segment of the North Anatolian Fault Zone from polarization analysis of P waves, Earth Planet. Sc. Lett., 327–328, 17–22,, 2012. 

Bulut, F., Özener, H., Doğru, A., Aktuğ, B., and Yaltırak, C.: Structural setting along the Western North Anatolian Fault and its influence on the 2014 North Aegean Earthquake (Mw 6.9). Tectonophysics, 745, 382–394,, 2018. 

Bürgmann, R., Ayhan, M. E., Fielding, E. J., Wright, T. J., McClusky, S., Aktug, B., Demir, C., Lenk, O., Türkezer, A.: Deformation during the 12 November 1999 Düzce, Turkey, Earthquake, from GPS and InSAR Data. B. Seismol. Soc. Am., 92, 161–171,, 2002. 

Cakir, Z., Akoglu, A., Belabbes, S., Ergintav, S., and Meghraoui, M.: Creeping along the Ismetpasa section of the North Anatolian fault (Western Turkey): Rate and extent from InSAR. Earth Planet. Sc. Lett., 238, 225–234,, 2005. 

Cetin, E., Cakir, Z., Meghraoui, M., Ergintav, S., and Akoglu, A. M.: Extent and distribution of aseismic slip on the Ismetpaşa segment of the North Anatolian Fault (Turkey) from Persistent Scatterer InSAR, Geochem. Geophy. Geosy., 15, 2883–2894,, 2014. 

Daniel, G., Marsan, D., and Bouchon, M.: Earthquake triggering in southern Iceland following the June 2000 Ms6.6 doublet, J. Geophys. Res.-Sol. Ea., 113,,, 2008. 

Di Bona, M. and Rovelli, A.: Effects of the bandwidth limitation of stress drops estimated from integrals of the ground motion, B. Seismol. Soc. Am., 78, 1818–1825, 1988. 

Douglas, A., Hudson, J. A., and Pearce, R. G.: Directivity and the Doppler effect, B. Seismol. Soc. Am., 78, 1367–1372, 1988. 

Dresen, G., Kwiatek, G., Goebel, T., and Ben-Zion, Y.: Seismic and Aseismic Preparatory Processes Before Large Stick–Slip Failure, Pure Appl. Geophys., 177, 5741–5760,, 2020. 

Durand, V., Bentz, S., Kwiatek, G., Dresen, G., Wollin, C., Heidbach, O., Martínez-Garzón, P., Cotton, F., Nurlu, M., Bohnhoff., M.: A Two-Scale Preparation Phase Preceded an Mw 5.8 Earthquake in the Sea of Marmara Offshore Istanbul, Turkey, Seismol. Res. Lett., 91, 3139–3147,, 2020. 

Dziewonski, A. M., Chou, T.- A., and Woodhouse, J. H.: Determination of earthquake source parameters from waveform data for studies of global and regional seismicity, J. Geophys. Res., 86, 2825–2852,, 1981. 

Ekström, G., Nettles, M., and Dziewonski, A. M.: The global CMT project 2004-2010: Centroid-moment tensors for 13,017 earthquakes, Phys. Earth Planet. In., 200–201, 1–9,, 2012. 

Ellsworth, W. L. and Bulut, F.: Nucleation of the 1999 Izmit earthquake by a triggered cascade of foreshocks, Nat. Geosci., 11, 531,, 2018. 

Emre, Ö., Duman, T. Y., Özalp, S., Şaroğlu, F., Olgun, Ş., Elmacı, H., and Çan, T.: Active fault database of Turkey, Bull. Earthq. Eng., 16, 3229–3275,, 2018. 

Ergintav, S., McClusky, S., Hearn, E., Reilinger, R., Cakmak, R., Herring, T., Özener, H., Lenk, O., and Tari, E.: Seven years of postseismic deformation following the 1999, M = 7.4 and M = 7.2, Izmit-Düzce, Turkey earthquake sequence, J. Geophys. Res.-Sol. Ea., 114,,, 2009. 

Eyidoğan, H.: What happened in Düzce? A preliminary assessment of the November 23 quake in Western Turkey, Temblor,, 2022. 

García-García, J., Romacho, M., and Jiménez, A.: Determination of near-surface attenuation, with κ parameter, to obtain the seismic moment, stress drop, source dimension and seismic energy for microearthquakes in the Granada Basin (Southern Spain), Phys. Earth Planet. In., 141, 9–26,, 2004. 

Gulia, L. and Wiemer, S.: Real-time discrimination of earthquake foreshocks and aftershocks, Nature, 574, 193–199,, 2019. 

Hainzl, S., Ben-Zion, Y., Cattania, C., and Wassermann, J.: Testing atmospheric and tidal earthquake triggering at Mt. Hochstaufen, Germany. J. Geophys. Res.-Sol. Ea., 118, 5442–5452,, 2013. 

Hanks, T. C. and Kanamori, H.: A moment magnitude scale, J. Geophys. Res.-Sol. Ea., 84, 2348–2350, 1979. 

Haskell, N. A.: Total energy and energy density of elastic wave radiation from propagating faults, B. Seismol. Soc. Am., 54, 1811–1841, 1964. 

Hauksson, E. and Shearer, P. M.: Attenuation models (QP and QS) in three dimensions of the southern California crust: Inferred fluid saturation at seismogenic depths, J. Geophys. Res., 111, B05302,, 2006. 

Ickrath, M., Bohnhoff, M., Dresen, G., Martínez-Garzón, P., Bulut, F., Kwiatek, G., and Germer, O.: Detailed analysis of spatiotemporal variations of the stress field orientation along the Izmit-Düzce rupture in NW Turkey from inversion of first-motion polarity data, Geophys. J. Int., 202, 2120–2132,, 2015. 

Jara, J., Socquet, A., Marsan, D., and Bouchon, M.: Long-Term Interactions Between Intermediate Depth and Shallow Seismicity in North Chile Subduction Zone, Geophys. Res. Lett., 44, 9283–9292,, 2017. 

Jolivet, R., Jara, J., Dalaison, M., Rouet-Leduc, B., Özdemir, A., Dogan, U., Çakir, Z., Ergintav, S., and Dubernet, P.: Daily to Centennial Behavior of Aseismic Slip Along the Central Section of the North Anatolian Fault, J. Geophys. Res.-Sol. Ea., 128, 1–35,, 2023. 

Kadirioğlu, F. T. and Kartal, R. F.: The new empirical magnitude conversion relations using an improved earthquake catalogue for Turkey and its near vicinity (1900–2012), Turk. J. Earth Sci., 25, 300–310,, 2016. 

Karabulut, H., Bouchon, M., and Schmittbuhl, J.: Synchronization of small-scale seismic clusters reveals large-scale plate deformation, Earth Planets Space, 74, 158,, 2022. 

KOERI Observatory: Seismicity catalog, (last access: 20 October 2023), 2023. 

Konca, A. O., Leprince, S., Avouac, J.-P., and Helmberger, D. V.: Rupture Process of the 1999 Mw 7.1 Duzce Earthquake from Joint Analysis of SPOT, GPS, InSAR, Strong-Motion, and Teleseismic Data: A Supershear Rupture with Variable Rupture Velocity, B. Seismol. Soc. Am., 100, 267–288,, 2010. 

Kurt, A. I., Özbakir, A. D., Cingöz, A., Ergintav, S., Doğan, U., and Özarpaci, S.: Contemporary Velocity Field for Turkey Inferred from Combination of a Dense Network of Long Term GNSS Observations, Turk. J. Earth Sci., 32, 4,, 2022. 

Kwiatek, G. and Ben-Zion, Y.: Assessment of P and S wave energy radiated from very small shear-tensile seismic events in a deep South African mine, J. Geophys. Res.-Sol. Ea., 118, 3630–3641,, 2013. 

Kwiatek, G., Plenkers, K., Dresen, G., and JAGUARS Research Group: Source parameters of picoseismicity recorded at Mponeng deep gold mine, South Africa: Implications for scaling relations, B. Seismol. Soc. Am., 101, 2592–2608, 2011. 

Kwiatek, G., Martínez-Garzón, P., Dresen, G., Bohnhoff, M., Sone, H., and Hartline, C.: Effects of long-term fluid injection on induced seismicity parameters and maximum magnitude in northwestern part of The Geysers geothermal field, J. Geophys. Res.-Sol. Ea., 120, 7085–7101, 2015. 

Lomax, A., Michelini, A., and Curtis, A.: Earthquake Location, Direct, Global-Search Methods, in: Complexity In Encyclopedia of Complexity and System Science, Part 5, Springer, New York, 2449–2473,, 2009. 

Malin, P. E., Bohnhoff, M., Blümle, F., Dresen, G., Martínez-Garzón, P., Nurlu, M., Ceken, U, Kadirioglu, F. T., Kartal, R., Kilic, T., and Yanick, K.: Microearthquakes preceding a M4.2 Earthquake Offshore Istanbul, Sci. Rep.-UK, 8, 16176,, 2018. 

Marsan, D., Bouchon, M., Gardonio, B., Perfettini, H., Socquet, A., and Enescu, B.: Change in seismicity along the Japan trench, 1990–2011, and its relationship with seismic coupling, J. Geophys. Res.-Solid, 122, 4645–4659,, 2017. 

Martínez-Garzón, P., Bohnhoff, M., Ben-Zion, Y., and Dresen, G.: Scaling of maximum observed magnitudes with geometrical and stress properties of strike-slip faults, Geophys. Res. Lett., 42, 2015GL066478,, 2015. 

Martínez-Garzón, P., Beroza, G. C., Bocchini, G. M., and Bohnhoff, M.: Sea level changes affect seismicity rates in a hydrothermal system near Istanbul, Geophys. Res. Lett., 50, e2022GL101258,, 2023. 

Murru, M., Akinci, A., Falcone, G., Pucci, S., Console, R., and Parsons, T.: M 7 earthquake rupture forecast and time-dependent probability for the sea of Marmara region, Turkey, J. Geophys. Res.-Sol. Ea., 2015JB012595,, 2016. 

Najdahmadi, B., Bohnhoff, M., and Ben-Zion, Y.: Bimaterial interfaces at the Karadere segment of the North Anatolian Fault, northwestern Turkey, J. Geophys. Res.-Sol. Ea., 121, 2015JB012601,, 2016. 

Ogata, Y.: Statistical Models for Earthquake Occurrences and Residual Analysis for Point Processes, J. Am. Stat. Assoc., 83, 9–27,, 1988. 

Ogata, Y. and Katsura, K.: Analysis of temporal and spatial heterogeneity of magnitude frequency distribution inferred from earthquake catalogues. Geophys. J. Int., 113, 727–738,, 1993. 

Olsen, K. B., Day, S. M., and Bradley, C. R.: Estimation of Q for long-period (> 2 sec) waves in the Los Angeles basin, B. Seismol. Soc. Am., 93, 627–638, 2003. 

Özalp, S. and Kürçer, A.: 23 Kasım 2022 Gölyaka (Düzce) Depremi (Mw 6,0) Saha Gözlemleri ve Değerlendirme Raporu, MTA Jeoloji Etütleri Dairesi Başkanlığı, Ankara, Türkiye, 30 s, 2022 (in Turkish). 

Pucci, S., Pantosti, D., Barchi, M., Palyvos, N.: Evolution and complexity of the seismogenic Düzce fault zone (Turkey) depicted by coseismic ruptures, Plio-Quaternary structural pattern and geomorphology, Geophys. Res. Abstr. 808339, 2006. 

Raub, C., Martínez-Garzón, P., Kwiatek, G., Bohnhoff, M., and Dresen, G.: Variations of seismic b-value at different stages of the seismic cycle along the North Anatolian Fault Zone in northwestern Turkey, Tectonophysics, 712–713, 232–248,, 2017. 

Ross, Z. E., Cochran, E. S., Trugman, D. T., and Smith, J. D.: 3D fault architecture controls the dynamism of earthquake swarms, Science, 368, 1357–1361,, 2020. 

Salvatier, J., Wiecki, T. V., and Fonnesbeck, C.: Probabilistic programming in Python using PyMC3, PeerJ Computer Science, 2, e55,, 2016. 

Sato, T. and Hirasawa, T.: Body wave spectra from propagating shear cracks, J. Phys. Earth, 21, 415–432, 1973. 

Savage, J. C.: Relation of corner frequency to fault dimensions, J. Geophys. Res., 77, 3788–3795, 1972. 

Sengör, A. M. C.: The North Anatolian Fault: a new look, Annu. Rev. Earth Pl. Sc., 33, 37–112, 2005. 

Shearer, P.: Introduction to Seismology, in: 3rd Edn., Cambridge University Press, Cambridge,, 2019.  

Snoke, J. A.: Stable determination of (Brune) stress drops, B. Seismol. Soc. Am., 77, 530–538, 1987. 

Stirling, M. W., Wesnousky, S. G., and Shimazaki, K.: Fault trace complexity, cumulative slip, and the shape of the magnitude-frequency distribution for strike-slip faults: a global survey, Geophys. J. Int., 124, 833–868,, 1996. 

Türker, E., Yen, M., Pilz, M., and Cotton, F.: Significance of Pulse‐Like Ground Motions and Directivity Effects in Moderate Earthquakes: The Example of the Mw 6.1 Gölyaka–Düzce Earthquake on 23 November 2022, Bull. Seismol. Soc. Am.,, in press, 2023. 

Utsu, T. and Seki, A.: Relation between the area of aftershock region and the energy of the mainshock, Zisin, 7, 233–240, 1955. 

van der Elst, N. J. and Shaw, B. E.: Larger aftershocks happen farther away: Nonseparability of magnitude and spatial distributions of aftershocks, Geophys. Res. Lett., 42, 5771–5778, 2015. 

Waldhauser, F. and Ellsworth, W. L.: A double-difference earthquake location algorithm: Method and application to the northern Hayward fault, B. Seismol. Soc. Am., 90, 1353–1368, 2000. 

Waldhauser, F., Ellsworth, W. L., Schaff, D. P., and Cole, A.: Streaks, multiplets, and holes: high-resolution spatiotemporal behavior of Parkfield seismicity, Geophys. Res. Lett., 31, L18608,, 2004. 

Wesnousky, S. G.: Seismological and structural evolution of strike-slip faults, Nature, 335, 340–343,, 1988. 

Wiemer, S. and Wyss, M.: Minimum magnitude of completeness in earthquake catalogs: Examples from Alaska, the Western United States & Japan, Bull. Seismol. Soc. Am., 90, 859–869, 2000. 

Wu, C., Meng, X., Peng, Z., and Ben-Zion, Y.: Lack of Spatiotemporal Localization of Foreshocks before the 1999 Mw 7.1 Düzce, Turkey, Earthquake, B. Seismol. Soc. Am., 104, 560–566,, 2013. 

Zhu, W. and Beroza, G. C.: PhaseNet: a deep-neural-network-based seismic arrival-time picking method, Geophys. J. Int., 216, 261–273,, 2019. 

Zhu, W., McBrearty, I. W., Mousavi, S. M., Ellsworth, W. L., and Beroza, G. C.: Earthquake Phase Association Using a Bayesian Gaussian Mixture Model, J. Geophys. Res.-Solid, 127, e2021JB023249,, 2022. 

Short summary
We analyze the 2022 MW 6.0 Gölyaka sequence. A high-resolution seismicity catalog revealed no spatiotemporal localization and lack of immediate foreshocks. Aftershock distribution suggests the activation of the Karadere and Düzce faults. The preferential energy propagation suggests that the mainshock propagated eastwards, which is in agreement with predictions from models, where the velocity in the two sides of the fault is different.