Reservoir-Triggered Earthquakes Around the Atatürk Dam (Southeastern Turkey)

Introduction

Discriminating induced or triggered seismicity related to industrial activities from natural seismicity has been a highly debated subject. Since the beginning of the last century, many earthquakes associated with anthropogenic activities have been reported, and the number of cases has been increasing due to the expanding man-made operations, such as gas and oil production, wastewater injection, mining, geothermal operations, and water impoundment (Dahm et al., 2010; Grigoli et al., 2017; Rinaldi et al., 2020). The most recent, outstanding cases debating potential induced or triggered seismicity, attracting societal interest, include the 2011 M_w 5.7 and 2016 M_w 5.8 Oklahoma earthquake sequences (Ellsworth, 2013; Keranen et al., 2013, 2014; Walsh and Zoback, 2015; Manga et al., 2016; Yeck et al., 2016, 2017), the 2012 M_w 6.1 and 5.9 Emilia, Italy, earthquakes (Cesca et al., 2013a; Dahm et al., 2015; Juanes et al., 2016), the 2017 M_w 5.5 Pohang, South Korea, earthquake (Grigoli et al., 2018; Kim et al., 2018), the 2011 M_w 5.1 Lorca, Spain, earthquake (González et al., 2012; Martínez-Díaz et al., 2012), the 2013 M_w 4.3 Castor, Spain, earthquake sequence (Cesca et al., 2014; Gaite et al., 2016; Villaseñor et al., 2020), and the 2012 M_L 3.6 Huizinge earthquake at the Groningen gas field (Richter et al., 2020). Many more reported cases caused by human-related activities have been compiled in different studies (McGarr et al., 2002; Davies et al., 2013; Ellsworth, 2013; Foulger et al., 2018).

In recent years, some attempts have been made to differentiate induced and triggered seismicity (McGarr and Simpson, 1997; McGarr et al., 2002; Dahm et al., 2013, 2015; Shapiro et al., 2013). In the case of induced earthquakes, the nucleation, growth, and rupture process are determined by human-related stress perturbations (Dahm et al., 2013). In the case of triggered seismicity, the background stress field plays a more important role, and human activities are only responsible for the earthquake nucleation, while the rupture evolution is controlled by the background stresses (Dahm et al., 2013). This latter case may include large earthquakes, which could be triggered by small perturbations near their nucleation point, but then grow considerably, with the final size and magnitude not being controlled by the original anthropogenic stress changes but depending on fault dimensions and strain (Dahm et al., 2013; Grigoli et al., 2017). In our study, we use the term reservoir triggered seismicity (RTS) for the earthquakes that occurred close to the Atatürk Dam, as this region is located between tectonically active faults, and the background tectonic stresses presumably play a role in the size and magnitude of observed seismicity (Figure 1).

FIGURE 1

Figure 1. Map of the Atatürk Dam and its vicinity. Thick red lines indicate active faults in the region obtained from the European Database of Seismogenic Faults (EDSF; Basili et al., 2013), and thin red lines illustrate the local faults; SF: Samsat Fault, KF: Kalecik Fault: LF: Lice Fault (General Directorate of Mineral Research and Exploration of Turkey; Perinçek et al., 1987; MTA, 2020). Yellow circles show historical earthquakes (Soysal et al., 1981; Ambraseys, 1989; Ambraseys and Jackson, 1998; the improved historical earthquake catalog of Turkey and its surroundings (https://deprem.afad.gov.tr/tarihseldepremler; AFAD, 2020; last accessed June 2020). The green star shows the location of the recent 2020 Elazığ-Sivrice earthquake caused by the reactivation of EAF, which had been silent for more than one century. Black arrows show GPS velocity vectors in the area (McClusky et al., 2000). The inset panel shows tectonic plates and boundaries surrounding Turkey, where red lines indicate main plate boundaries (Bird, 2003), and black arrows show the relative motion of Arabian and Anatolian Plates roughly. The blue rectangle shows the study area, which is enlarged on the map.

An influence of water reservoir loading on earthquake activity was first proposed by Carder (1945) at Lake Mead, United States. Many case studies of RTS have been reported since that time; most known RTS cases (M_w > 6) were observed at Xinfengjiang Dam–China, 1962 M_w 6.2, Kariba Zambia–Zimbabwe, 1963 M_w 6.2, Koyna Dam–India, 1967 M_w 6.3, and Zipingpu Reservoir–Wenchuan, 2008 M_w 7.9 (Gupta and Rastogi, 1976; Gupta, 1992, 2002; Ge et al., 2009). Wilson et al. (2017) have recently constructed a database with 186 reported cases of RTS (The Human-Induced Earthquake Database HiQuake)¹. In Turkey, although one of the richest countries in geothermal, mining, and water resources potentials, only a few case studies of induced/triggered earthquakes have been reported in the literature so far.

Around 860 active dam sites are currently existing in Turkey, and the number is expected to increase (The General Directorate of State Hydraulic Works of Turkey; DSI)². Most of these dams are located in southeast Turkey since the Southeastern Anatolia Project (GAP) was launched in 1977. Today, the Atatürk Dam is the fifth largest dam on Earth in terms of water storage capacity (48.7 billion m³) and among the largest dam sites in terms of electricity production (DSI)³. It is also the largest clay-cored rock-fill dam in Turkey, with 169 m height. The construction of the Atatürk Dam was initiated in 1983, the water impoundment started in 1990, and the dam became operational in 1992 (Tosun et al., 2007; Tosun, 2012). The annual variation of the water level in the Atatürk Dam is in the range of 30 m, between 513 m and 542 m above sea level (DSI)³. Table 1 summarizes the characteristics of the Atatürk Dam and its reservoir. After the water impoundment in the reservoir, field and laboratory experiments showed damages along the crest (Çetin et al., 2000) attributed to the rising amount of water. Consequently, the rock-fill part of the dam started to slake by May 1992; the upper part of the dam was then reconstructed to its original height (549 m), and the dam was maintained operational by keeping a 7 m freeboard (Çetin et al., 2000).

TABLE 1

Table 1. The characteristics of Atatürk Dam and its reservoir (The General Directorate of State Hydraulic Works of Turkey; DSI; http://www.ataturkbaraji.com; last accessed April 2021; asl* = above sea level).

The Atatürk Reservoir (AR) is located on the Euphrates River between the Adıyaman-Samsat region and the Şanlıurfa province, in southeast Turkey. This region is tectonically influenced by the relative motion of the African, Arabian, and Eurasian Plates resulting in the movement of the Anatolian Plate to the west (McKenzie, 1972; Şengör and Yılmaz, 1981; Şengör et al., 1985; Reilinger et al., 2006) as shown in Figure 1. AR is situated between two fault systems: the Bozova Fault (BF), which is an NW-SE right-lateral strike-slip structure with 60 km length passing through the southwest of AR (Şahbaz and Seyitoğlu, 2018), and the East Anatolian Fault (EAF), in a distance of about 60 km in the northeast of AR, which is a ∼580 km left-lateral strike-slip active fault striking NE-SW direction and one of the most prominent structural elements in the region (Arpat and Şaroğlu, 1972; Duman and Emre, 2013). A recent destructive earthquake (M_w 6.8) occurred on January 24, 2020, along the EAF, within approximately 100 km distance from AR, and caused serious damages not only in the epicentral area but also in the neighboring regions (Jamalreyhani et al., 2020). Figure 1 also shows the existence of local faults dominated by the regional tectonics in the study area. The Samsat Fault (SF) and Kalecik Fault (KF) showing parallel alignment to the BF and crossing the AR in the NW-SE direction with a right-lateral strike-slip mechanism. On the other hand, The Lice Fault (LF) indicates a left-lateral strike-slip mechanism in the NE-SW direction with respect to the BF, SF, and KF (Perinçek et al., 1987; Kartal and Kadirioğlu, 2019; Irmak et al., 2020).

Considering the historical seismicity (B.C. 1800–A.D. 1905), no strong (M ≥ 6.0) and damaging earthquakes were reported near AR (Soysal et al., 1981). Historical earthquakes occurred mainly in the EAF to the north of the dam in 1866, 1893, 1905 (Figure 1). The fault zone was remarkably inactive during the 20^th century (Ambraseys, 1989; Ambraseys and Jackson, 1998), facing earthquakes with magnitude up to M_w 6.6–6.8, until the most recent destructive M_w 6.8 Elazığ-Sivrice earthquake on January 24, 2020. On the other hand, some historical earthquakes occurred south of the dam, e.g., in 718, 1003, 1037 (Figure 1; data from the historical earthquake catalog of Turkey and its surroundings; AFAD)⁴.

On March 2, 2017, and April 24, 2018, two moderate earthquakes (M_w 5.5 and M_w 5.1, respectively) struck Samsat town near AR (Figure 1). These earthquakes were responsible for dozens of injuries and significant damages to buildings. The occurrence of the 2017 earthquake, which is the largest event in this region, and its potential anthropogenic source triggered the interest of seismologists. Reservoir-triggered seismicity around AR was first hypothesized by Eyidoğan et al. (2010) after the occurrence of the M_L 5.2 earthquake on September 3, 2008, which was the largest earthquake prior to the 2017–2018 earthquakes. They pointed out that small-magnitude earthquakes started to occur in the vicinity of the dam soon after the water level reached its first maximum in 1994. Furthermore, they depicted a clear anti-correlation between water-level change and seismicity in the region and suggested that the September 3, 2008 earthquake (M_w 5.0) might have been triggered upon a drastic decrease in the water level, accompanying the low rainfall in the summer of 2008. On the other hand, Kartal and Kadirioğlu, 2019 listed several earthquakes from the catalog of DSI local network mostly after the dam construction and claimed that the seismicity in the region occurs independently of changes in the water loading. Thus, it has remained a matter of debate whether or not the local seismicity near the AR is correlated with water level changes and its related effective stress changes.

In this study, we use new satellite altimetry open access data that has not been analyzed in previous studies (Database for Hydrological Time Series of Inland Waters, DAHITI; Schwatke et al., 2015). This database provides nowadays reliable and accurate water level data (average uncertainty <0.01 m) for AR over the long period 2002–2020, allowing us to compare water level and seismicity rate for 18 years in the region. Furthermore, we use sophisticated methods to obtain reliable earthquake characteristics (e.g., moment tensor solutions, focal depth estimations, and Gutenberg-Richter b-values) and decluster the earthquake catalog to remove aftershocks that might distort the signatures of RTS. Finally, we calculate the time-dependent Coulomb-stress changes at seismogenic depth resulting from water loading and pore-pressure diffusion based on the observed water-level evolution and extended reservoir geometry. Based on this analysis, the comparison of stress and seismicity pattern indicates a causal relationship between the reservoir and local seismicity.

Seismic Networks and Seismicity Distribution

Local seismicity around the AR is monitored by permanent networks of the Disaster and Emergency Management Authority (AFAD) (1990) and Kandilli Observatory and Earthquake Research Institute, Boğaziçi University (KOERI) (1971). The distribution of active broadband seismic stations around the dam is marked in Figure 2 for various periods. The network has been densified with time particularly after the March 2, 2017 M_w 5.5 earthquake when additional AFAD stations were installed. To get a complete earthquake catalog for our analysis, we combine the seismic catalogs from AFAD and KOERI networks between 37.3°–37.8°E and 38.1°–39.0°N (Figure 2). Both catalogs are compiled with the SeisAn - earthquake analysis software (Havskov and Ottemöller, 1999). Hypocenters and magnitudes are updated according to the travel time residuals (RMS) and the number of stations used in the location. The magnitude of completeness (M_C), being ∼3.5 in the beginning, decreased with the densification of the network to an ultimate value of about 1.95 in 2017 (Figure 3 and Supplementary Figure 1). Taking into account only the largest events with M ≥ 3.5, which are homogeneously detected over the whole period, the earthquake rate has significantly increased in the AR region between 1990 and 2020 (Figure 3B) as previously pointed out in the study of Eyidoğan et al. (2010).

FIGURE 2

Figure 2. (A–D) Long term (1905–2020) spatio-temporal evolution of recorded earthquakes with their magnitudes and depth cross-sections based on a joint KOERI-AFAD catalog near Atatürk Reservoir. The spatial earthquake distribution of historical and instrumental seismic activity around AR for (A) the period 1905–1990 before the water filling, (B) 1990–2004 after the impoundment of the dam – depth versus longitude in panels (A–C) 2004–2017 and two orthogonal cross-sections crossing the dam (AA′,BB′). The brown dots here show the events during 2004–2008 to highlight the background seismicity, the M_w 5.0 2008 earthquake, and its aftershocks (D) 2017–2020, same cross sections as in panel (C). The existing broadband seismic stations during the corresponding periods are marked by the green (AFAD stations) and orange (KOERI stations) triangles. Yellow and pink stars refer to the hypocenters of March 2, 2017, M_w 5.5 and April 24, 2018, M_w 5.1 earthquakes, respectively. The white dashed rectangle marks the region where the seismic catalogs are combined (37.3°-37.8°E, 38.1°-39.0°N). The red lines represent the Bozova, Samsat, Kalecik, and Lice Faults (see Figure 1 to detail).

FIGURE 3

Figure 3. Frequency-magnitude distribution of the observed earthquakes: (A) Number of seismic stations with a distance less than 50, 100, or 150 km from the 2017 M_w 5.5 epicenter as a function of time; (B) cumulative number of M ≥ 3.5 events and (C) magnitudes of all recorded events versus time, where colors refer to the different periods analyzed separately in panel (D) and the horizontal dashed line indicates the completeness magnitude M_C = 2.75; (D) histograms and cumulative distributions with GR-fits (dashed lines, the corresponding b-values are provided in the legend). Additional plots of the frequency-magnitude distributions and the b-value calculation are presented in the Supplementary Figure 1.

The spatial pattern of the seismicity around the AR is demonstrated in Figures 2A–D. The lack of recorded seismic activity before the start of the dam impoundment in 1990 is shown in Figure 2A. After the dam operation begins and the reservoir is filled, the seismicity rate increases, especially in the vicinity of the AR (Figures 2B,C). The densification of the seismic stations after the March 2, 2017 M_w 5.5 has led to a significant improvement in the detection of weak events, contributing to the observed increase in detected events (Figure 2D). Most earthquakes are shallow, predominantly occurring at less than 11 km depth. Such shallow seismicity around the dam hints that the water filling and leakage with the consequent induced changes in pore-pressure and effective stress can have affected the local seismicity.

Statistical Analysis of the Earthquake Catalog

Frequency-Magnitude Distribution

The Gutenberg-Richter (GR) relation states that the number of events (N) above a certain magnitude (M) follows the simple relation log(N) = a – bM. While the a-value, which defines the earthquake production rate above M = 0, varies largely between different tectonic regions, the b-value is usually found to be rather universal and scatters around 1.

Figure 3D shows the frequency-magnitude distributions of the earthquakes together with GR-fits for three successive periods with progressively decreasing M_C, from M_C = 2.75 in the period 2004–2012, M_C = 2.45 for 2012–2017, to M_C = 1.95 for t > 2017. At the same time, the b-value of GR distribution is estimated as 1.8, 1.4, and 0.9 for the periods of 2004–2012, 2012–2017, and 2017–2021, respectively (see also Supplementary Figure 1). Thus the b-value is found to systematically decrease with increasing time after the impoundment of the dam, which might be a result of increased stresses in the crust as previously hypothesized (Schorlemmer et al., 2005; Scholz, 2015).

Anti-Correlation Between Water-Level and Declustered Seismicity

For a detailed analysis of the spatial and temporal characteristics of the earthquake activity with the reservoir, aftershocks which are triggered by preceding earthquakes should be removed from the catalog. The remaining declustered earthquakes are the so-called background events, which are related to tectonic stressing or transient aseismic forcing such as reservoir-induced stress changes. To separate aftershocks and background events, we use an established scheme based on nearest-neighbor distances (Baiesi and Paczuski, 2004, 2005; Zaliapin et al., 2008; Zaliapin and Ben-Zion, 2013). This is a purely statistical method that does not rely on any particular aftershock-triggering mechanism, such as static/dynamic coseismic stress changes or afterslip. The method is described in more detail in Appendix A.

After declustering, we compare the temporal evolution of earthquake activity of the complete part of the catalog (M ≥ 2.8 for t > 2004) with the water level variations. Figure 4 shows the monthly seismicity rate in comparison to the reversed water-level variations. A significant anti-correlation between earthquake rate and water-level is observed in the period between 2004 and 2014 with a correlation coefficient of r = −0.48. Afterward, the anti-correlation becomes weak, particularly because of the significantly reduced activity of M ≥ 2.8 background events (see also Supplementary Figure 2).

FIGURE 4

Figure 4. Time evolution of the water level and the earthquake activity (declustered, M ≥ 2.8). Here the y-scale is reversed with values referring to the changes relative to the start of the impoundment. A similar plot with a normal y-scale is provided in Supplementary Figure 2. The corresponding correlation coefficients (r) for the periods 2004–2014 and 2014–2020 are provided in the title line together with the corresponding p-value (significance for p < 0.05).

Analysis of Source Properties

Construction of a Local Velocity Model

Firstly, we estimate a local 1D velocity model for our further analysis of focal depths and source mechanisms. For that purpose, we use PyVelest⁵, based on the travel time inversion VELEST program (Kissling et al., 1994). We carefully select 470 well-located earthquakes in the region with root mean squared (RMS) misfit of the solution ≤0.5 s and azimuthal gap ≤180°, respectively. The recent model of Acarel et al. (2019) is used as the initial model since it is the closest model to the reservoir area. By perturbing velocities in the ±0.3 km/s range, randomly 500 synthetic velocity models are generated and inverted for the study region. Consequently, a well-defined velocity model from the mean of inverted models is constructed. The final velocity model has a good agreement with the reference model, with velocity changes not exceeding 0.2 km/s for each corresponding layer. Details of the selected events, ray coverage and final 1D velocity model are presented in the Supplementary Figures 3–5 and Table 1.

Moment Tensor Inversion

The inversion and further decomposition of regional moment tensor solutions can be used to discuss cases of natural or anthropogenic seismicity (Cesca et al., 2013b). Here, moment tensor inversion has been performed using a probabilistic inversion method, provided by the software Grond⁶ (Heimann et al., 2017, 2018). This method has been successfully applied in different studies (e.g., Dahm et al., 2018; Jamalreyhani et al., 2019, 2021; Cesca et al., 2020; Dost et al., 2020; Kühn et al., 2020; López-Comino et al., 2021) and described in some of them (e.g., Dahm et al., 2018; Dost et al., 2020; Kühn et al., 2020).

We obtain full moment tensor solutions for the two largest events, namely the M_w 5.5 2017 and M_w 5.1 2018 earthquakes, using 3-components waveform inversion in the time domain and the frequency band of 0.02–0.05 Hz. A prior data quality assessment is applied to all stations to prevent systematic errors in the moment tensor solutions due to sensor misorientations (Büyükakpınar et al., 2021). An example of waveform fits and MT solution is illustrated in Figure 5. Additionally, the waveform fits and MT solution are shown for the M_w 5.1 2018 earthquake in the Supplementary Figures 6, 7. The full moment tensor decomposition (ISO-CLVD-DC) components show a relatively large CLVD component, 17 and 41% for 2017 and 2018 earthquakes, respectively. Furthermore, moment tensors for 66 weaker events down to M 2.8 were calculated, applying a simplified double couple representation and frequency band of 0.06–0.11 Hz. This is possible due to the presence of the densified networks in the area after 2017. All obtained source parameters are listed in Supplementary Table 2. Figure 6 shows the result of all focal mechanisms. Most of the solutions are characterized by strike-slip mechanisms with relatively shallow centroid depths (<6 km and uncertainty <1 km, see Supplementary Table 2).

FIGURE 5

Figure 5. Waveform fits in the time domain for the 2017-08-31 (UTC) 16:11:09 M 3.1 earthquake. Red and gray waveforms represent synthetic and observed records, respectively. The fuzzy MT shows the solution with its uncertainties (see also Supplementary Figures 6, 7).

FIGURE 6

Figure 6. (A) Focal mechanism solutions for 68 events that occurred between 2017 and 2020. The recorded epicenters are color-coded in time. The red lines show the faults with their slip direction. (B) The cross-section of the (A-B) profile in the study area showing the centroid depths mostly less than 6 km.

For the largest earthquake, the 2017 M_w 5.5 event, we study the distribution of the direct aftershocks to determine the rupture plane. The aftershocks are found to extend in the NW direction in agreement with the nodal plane with a strike of 313.5° (Figure 7). Only minor activity is found in the strike direction of the second nodal plane with a strike of 46.5°, which is interpreted as the auxiliary plane. Both the focal mechanism strike and the trend of aftershocks are in general agreement with the strike direction of the Bozova, Samsat, and Kalecik faults (see Figure 6). Thus we conclude that the rupture mechanism of the largest event is best described by strike = 313.5°, dip = 64.3°, and rake = 173.1°. These values are used to estimate stresses on receiver faults in the following chapter.

FIGURE 7

Figure 7. Aftershocks of March 2, 2017, M_w 5.5 earthquake: (A) Spatial distribution of the first ten-day aftershocks, which were directly triggered by the mainshock according to the declustering method (green points), and all aftershocks occurred within the first day (blue points). The epicenter of the mainshock is marked by the star. (B) Rose diagram of the directions of the same events relative to the mainshock location, where the preferred strike value of the mainshock focal mechanism is marked as a dashed black line.

Focal Depths of the M_w 5.5 2017 and M_w 5.1 2018 Events

Depth values of seismic catalogs may have large uncertainties and often suffer from trade-offs between origin time and depth. A better depth constraint is substantial to discriminate induced seismicity from natural ones. We adopt here a method for an accurate source depth estimation based on the delay between direct P phases and surface reflected (pP and sP) phases at teleseismic distances. Since the waveforms of these moderate events at large distances are weak, we investigate these delays at the location of seismic arrays, where waveforms of many stations can be shifted and stacked to improve the signal-to-noise (SNR) ratio. We rely on the Abedeto algorithm⁷, which has been previously used in similar studies (Negi et al., 2017; Braun et al., 2018; Gaebler et al., 2019). The qualitative comparison of observed and synthetic beams for different source depths (Figure 8), which are built assuming source and receiver specific crustal models and a global model for the propagation of the seismic waves in between, allows to estimate accurate focal depths. Here, the global crustal velocity model Crust 2.0 (Bassin et al., 2000) and the estimated earthquake source models are used to calculate the focal depths (see Figure 6 and Supplementary Table 2). In this way, we estimate the focal depths of 11 and 5 km for the 2017 and 2018 earthquakes, respectively. We find similar results for different arrays (see Supplementary Figures 8, 9). The obtained focal depths are compatible with the hypocentral depths given in the seismic catalogs.

FIGURE 8

Figure 8. (A,B) Modeled teleseismic depth phases for the BCA array for different assumed source depth and fixed source mechanism (black lines). For the modeling of the source side crust, we use the velocity model (see Supplementary Figure 5 and Table 1) which we have obtained using PyVelest. Stacked array beams (blue lines) are seen consistent with synthetics for a depth of 11 and 5 km in the case of the M_w 5.5 2017 earthquake (A) and the M_w 5.1 2018 earthquake (B), respectively.

Reservoir-Induced Stress Changes

We calculate the reservoir-related variations of the Coulomb-Failure Stress (CFS) relative to the initial stress state at the beginning of the water impoundment. For that purpose, we assume a uniform and isotropic half-space and represent the Atatürk Reservoir by 246-point sources covering the reservoir surface (see Supplementary Figure 10). In particular, we calculate the stress induced by the water load using Boussinesq-Cerruti solutions (see, e.g., Deng et al., 2010), and pressure changes related to pore-pressure diffusion by convolution of the observed reservoir water level (extended to include the filling phase, see Supplementary Figure 10A) with the Green’s function (Gahalaut and Hassoup, 2012; Hainzl et al., 2015). Details of the stress calculation are provided in Appendix B. The main model parameter, which influences the stress evolution, is the hydraulic diffusivity (D), while the other parameters such as friction coefficient and Skempton coefficient only have a minor impact. We use a friction coefficient of μ = 0.8 and calculate CFS for the obtained rupture mechanism of the 2017 M_w 5.5 mainshock, namely strike = 313.5°, dip = 64.3°, and rake = 173.1° (see the previous section). This mechanism is assumed to be representative of the wider region due to the strong similarities of estimated focal mechanisms and the correspondence to the strike of the main regional fault.

At first, we calculate the induced Cauchy stress resulting from the water load alone, which is independent of the choice of the diffusivity (D). In particular, we determine the induced stresses for the estimated mainshock mechanism in the case of the reservoir with a 60 m water column ignoring any fluid diffusion. The result is shown in Figure 9 indicating that the whole region is unloaded at all depth layers due to water load. This means that all faults with the mainshock mechanism are firstly stabilized by the reservoir impoundment. As a counterpart, the pore-pressure diffusion leads to an increase of stress with time, but with some delay depending on the distance to the reservoir and the value of diffusivity (D). The stabilization effect of the water loading can explain the observed anti-correlation between water-level and seismicity rate. A sudden increase of the water-level leads to an immediate reduction of the CFS-value on the faults with the predominant mechanism because the related increase of the pore-pressure diffusion is delayed. This can also explain the decreased seismicity rate at peaks of the water levels.

FIGURE 9

Figure 9. Stress changes due to reservoir loading (60 m water column) ignoring fluid diffusion and calculated for receiver faults with the orientation of the 2017 M_w 5.5 mainshock mechanism (strike = 313.5°, dip = 64.3°, rake = 173.1°). Contour lines refer to stress values in units of MPa. Points refer to M ≥ 2.8 earthquakes (black = background, gray = aftershocks) in the depth range indicated in each title line. The colored crosses mark the three locations for which Figure 10 shows the total stress history, including pore pressure changes.

While the stabilization effect is found to dominate in the short-term, pore-pressure diffusion leads to an increase of the effective stress with time. This is demonstrated in Figure 10A, where the total Coulomb stress is calculated for three different values of the hydraulic diffusivity (D = 0.05, 0.1, and 1.0 m²/s) at 5 km depth in the three locations marked in Figure 9. All curves firstly show negative values, but later a crossover to positive values is observed, after which induced seismicity might be expected. The time of this crossover is strongly dependent on the assumed D-value and the distance to the reservoir. For the location beneath the reservoir (blue cross in Figure 9 and blue lines in Figure 10), it already occurs between 1993 and 2000 for D-values in the range between 0.05 and 1 m²/s. At the farthest distance (green cross) and smallest D-value (dashed line), it occurs only in 2016. Overall, the general increase due to the pore-pressure diffusion is modulated by the instantaneous stress changes induced by changes in the water level.

FIGURE 10

Figure 10. (A) Time evolution of the total Coulomb stress calculated at the three locations indicated in Figure 9 at 5 km depth, assuming the mainshock mechanism as a receiver. The colors of the lines refer to the location, while the line style refers to different diffusivity values (see legend); (B) Temporal increase of the maximum CFS value relative to the year 2004 (for the same cases). Note that the number of triggered events should be proportional to these curves at the given location according to the simple CFS-model. For comparison, the cumulative events of the recorded M ≥ 2.8 background events in the whole area are shown by the bold gray line (with a scale on the right).

The simplest seismicity model, which builds on CFS-values, assumes that the number of earthquakes is proportional to the stress change, if it is positive, while no triggering is expected for negative changes. Furthermore, considering stress shadowing (Kaiser effect), the model only assumes triggering if the absolute stress exceeds all precursory values. Thus the seismicity rate R(t) is proportional to the stressing rate, R(t)∼d/dt CFS(t), if CFS(t) > max (CFS (time < t)), otherwise R(t) = 0.

The cumulative number of events becomes simply N(t) ∼ max (CFS (time ≤ t)).

Based on this model, we calculate the expected number of events in the three locations marked in Figure 9. In Figure 10B, the result is shown for the cumulative number of earthquakes N(t) in the case of three different D-values (0.05, 0.1, and 1 m²/s), where we arbitrarily set the proportionality factor to 1. We cannot compare this model prediction pointwise with observations due to the limited number of recorded events. Thus, the shape of the resulting curves is compared to the evolution of the cumulative number of the homogeneously recorded M ≥ 2.8 background events in the whole region after 2004 (gray line). Despite some variations depending on the location and chosen D-value, the shapes have a similar tendency, with a steep increase in the first period and a flattening in the later period.

For the case of D = 0.1 m²/s and a depth of 5 km, we also calculate the CFS-values on a spatial grid at different time points. Figure 11 (left column) shows the total CFS-stress in our study area at the beginning of 2010, 2015, and 2020, respectively. Southwest of the dam, including the Bozova fault, the total CFS values remain negative for the whole period, and no earthquakes with the mainshock mechanism are expected; only a few events are observed in this region. On the other hand, the highest stresses of about 0.7 bar occur just beneath the centers of the two arms of the lake, while the absolute value is slightly smaller (approximately 0.4 bar) in the Samsat region, where the two largest earthquakes occurred in 2017 and 2018. The earthquake probability depends not on the absolute value, but the stress increases according to the simple CFS-model, which is shown in the right column of Figure 11. The contour lines show that this increase just focused on the Samsat region in the period of the mainshocks.

FIGURE 11

Figure 11. Left column: contour lines of the total CFS-stress [kPa] at the year 2010 (top), 2015 (middle), and 2020 (bottom), calculated for the mainshock mechanism at a depth of 5 km and D = 0.1 m²/s. Right column: increase of the maximum CFS value in the period 2005–2010 (top), 2010–2015 (middle), and 2015–2020 (bottom), which is proportional to the number of triggered earthquakes according to the simple CFS-model. For comparison, the epicenters of earthquakes (M ≥ 2.8, black = declustered, gray = aftershocks) recorded in the corresponding time intervals are plotted. In the bottom row, the star refers to the 2017 M_w 5.5 mainshock.

Discussion

Based on the combined catalog from regional networks in this study, the seismicity rate increases in the AR region between 1990 and 2020. However, it should be noted that the seismic network was sparse before the impoundment and only improved with time. On the other hand, the estimated magnitude of completeness (M_C) decreased from 3.5 in the beginning to approximately 1.95 after 2017 because of the seismic network’s densification in the reservoir area. Nevertheless, considering only the complete part of the catalog (M ≥ 3.5) indicates a significant activation of seismicity after the reservoir filling. Besides the changes in the seismicity rate, the shape of the magnitude distribution also changed with time. The estimated b-value successively decreased from a high value of 1.8 in the period of 2004–2012, to 1.4 in 2012–2017, and finally to 0.9 in 2017–2021, respectively. According to field and lab experiments, this significant b-value decrease may be explained by an increase in stress in the AR region (Schorlemmer et al., 2005; Scholz, 2015).

Our moment tensor solutions for 68 earthquakes, including the largest earthquakes and aftershocks down to M 2.8, indicate a clear dominance of strike-slip mechanisms that are in agreement with the results of previous studies (e.g., Eyidoğan et al., 2010; Kartal and Kadirioğlu, 2019; Irmak et al., 2020). The estimated mechanisms are consistent with the regional stress field (Figures 1, 6). The spatial distribution of early aftershocks suggests that the NW-SE fault orientation is the causative fault plane, which is in general agreement with the strike of the faults in the region (e.g., the Bozova and Samsat faults). In particular, the Samsat fault is most probably responsible for the 2017 and 2018 earthquake sources, as previously indicated by Irmak et al. (2020). Furthermore, the MT solutions in this study reveal shallow centroid depths mostly below 6 km. Shallow source depths are also estimated in the study of Irmak et al. (2020). Our full moment tensor results also show relatively large CLVD components, 17 and 41% for 2017 and 2018 earthquakes, respectively. Shallow source depths and the percentage of non-double components also support the triggering mechanism in the study area.

The focal depths of the largest 2017 and 2018 events (11 and 5 km, respectively) are independently confirmed using an accurate array beam technique. Shallow focal depths (<10 km) are also revealed in the study of Irmak et al. (2020). Shallow hypocenters are often found for seismicity induced or triggered by reservoir loading in other regions (Simpson et al., 1988; Lizurek et al., 2019; Ruiz-Barajas et al., 2019) and they can be linked to the significant damage in the vicinity of the AR in 2017 and 2018. As previously indicated in the study of Irmak et al. (2020) most of the epicenters are located close to the Samsat town and mapped along the Samsat fault, which has been known but not considered as an active fault in Turkey’s current faults map (Emre et al., 2018).

Induced gravitational stress due to surface loading/unloading can favor or inhibit normal faulting and thrust faulting (Simpson, 1986) but is expected to have only a small direct influence on vertical strike-slip faults. However, dipping fault planes increase the impact of the water loading also for strike-slip events. The dip of the 2017 M_w 5.5 Samsat earthquake is estimated to be 64.3°. Furthermore, pore pressure transients in response to water loading favor the occurrence of any faulting type with some delay depending on distance and hydraulic diffusivity.

The joint interpretation of the temporal evolution of water level and seismicity rate is crucial to identifying seismicity simulated by water reservoir operations and understanding its spatio-temporal evolution. Water level and seismicity rate have been analyzed here through a declustering method, resolving a clear anti-correlation between water level and seismicity by using recent additional datasets. This pattern was firstly observed by Eyidoğan et al. (2010), by comparing earthquakes between 1992 and 2009 with the water level information obtained by DSI. However, Kartal and Kadirioğlu (2019) found no evidence of any correlation between water load and earthquake activity pointing to induced/triggered seismicity; this might result from the fact that they analyzed and correlated seismicity and water level only for short periods, such as three months sequences in the intervals of August-October, 2008 and February–April, 2017. Analyzing too short periods does not allow to resolve correlations whenever the triggering mechanism requires a considerable temporal delay. Our study reveals a correlation between reservoir impoundment and triggered seismicity, which is attempted to discriminate from the induced seismicity here, in a more extended data period. We can particularly explain this anti-correlation by the immediate stabilization effect of the surface water load on the dominant rupture mechanism of regional crustal faults.

On the other hand, we show that pore-pressure diffusion increased the effective stress with time, leading to positive total stresses and fault destabilization, explaining the observed seismicity in the proximity of the Samsat fault where the estimated total stress is high, and the two largest events occurred in 2017 and 2018. Interestingly, the Bozova fault in the SW of the reservoir area remains so far in the stress shadow, which agrees with the low seismic activity observed in this region.

Conclusion

Two recent, damaging earthquakes, with magnitude M_w 5.5 and M_w 5.1, struck in close vicinity to the Atatürk Reservoir in 2017 and 2018, raising the question of whether they have been induced or triggered by water loading operations. In this study, we analyzed the spatio-temporal evolution of seismicity and earthquake source characteristics in relation to the stresses induced by the Atatürk Reservoir, one of the largest dam reservoirs on Earth.

The local seismicity rate has substantially increased after constructing the dam and its impoundment, which began in 1990. Despite the overall seismicity increase, our analysis confirms a clear anti-correlation among seismicity rate and water level on shorter time scales, which was so far debated (Eyidoğan et al., 2010; Kartal and Kadirioğlu, 2019). Our stress calculations show that the anti-correlation can be explained by the stabilization effect of the water load, while the overall seismicity activation is attributed to pore pressure diffusion. We also observe a significant b-value decrease with time after the impoundment operations with a reduction from 1.8 in 2004 to 0.9 in 2020, suggesting a progressive increase of the effective stresses due to increased pore pressure. Furthermore, moment tensor solutions show that the NW-SE oriented strike-slip mechanism, which is compatible with the general trend of the existing tectonic regime, is dominant in the dam area.

Our analysis provides strong indications that the observed seismicity is partly triggered by the impounding of the Atatürk Reservoir based on the data consisting of seismicity, source mechanisms, and long-term water level information (2002–2020) which is not previously taken into account in other studies.

This work shows how combining accurate seismicity analysis, stress estimations, and statistical approaches helps to better discriminate and understand reservoir-triggered seismicity. Our results provide a solid base to assess the seismic hazard near the Atatürk Dam in Turkey.

Data Availability Statement

All seismic and satellite altimetry data used in this study are free and open access to users via the International Federation of Digital Seismograph Networks (FDSN) client code (KO; https://doi.org/10.7914/SN/KO) in the databases of the Boğaziçi University Kandilli Observatory and Earthquake Research Institute (KOERI; http://eida.koeri.boun.edu.tr; last accessed September 2020), the Disaster and Emergency Management Authority (AFAD; https://tdvms.afad.gov.tr; last accessed September 2020), and the Database for Hydrological Time Series of Inland Waters (DAHITI; https://dahiti.dgfi.tum.de/en/; last accessed September 2020). Active faults in the region were taken from the European Database of Seismogenic Faults (EDSF; http://diss.rm.ingv.it/share-edsf; last accessed September 2020). Some of the maps were prepared using the Pyrocko toolbox (https://pyrocko.org) and GMT software (http://gmt.soest.hawaii.edu). The topography data were obtained from https://topex.ucsd.edu/WWW_html/srtm15_plus.html (last accessed May 2020). Bathymetry data were taken from https://www.gebco.net/data_and_products/gridded_bathymetry_data (last accessed May 2020).

Author Contributions

All the authors contributed to the work and discussed the results presented in this manuscript. All the authors performed the investigation and research.

Funding

This study was supported by the International Training Course “Seismology and Seismic Hazard Assessment” which has been funded by the GeoForschungsZentrum Potsdam (GFZ) and the German Federal Foreign Office through the German Humanitarian Assistance program, grant S08-60 321.50 ALL 03/19.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor declared a past co-authorship with one of the authors, SHa and TD.

Acknowledgments

We would like to thank Haluk Eyidoğan for his inspiration in this study. We are very grateful to Mustafa Aktar and Mohammad Mohseni Aref for constructive discussion and comments on this work. We also thank Claus Milkereit and Dorina Kroll. We also thank the editor AR and the reviewers SP and LV for their valuable comments and suggestions that helped us to improve the article.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2021.663385/full#supplementary-material

Footnotes

1. ^ https://inducedearthquakes.org; last accessed September 2020.
2. ^ https://www.dsi.gov.tr/Sayfa/Detay/754; last accessed April 2021.
3. ^ http://www.ataturkbaraji.com; last accessed April 2021.
4. ^ https://deprem.afad.gov.tr/tarihseldepremler; last accessed June 2020.
5. ^ https://github.com/saeedsltm/PyVelest; last accessed September 2020.
6. ^ https://pyrocko.org/grond
7. ^ https://github.com/HerrMuellerluedenscheid/ArrayBeamDepthTool; last accessed September 2020.

References

Acarel, D., Cambaz, M. D., Turhan, F., Mutlu, A. K., and Polat, R. (2019). Seismotectonics of Malatya Fault, Eastern Turkey. Open Geosci. 11, 1098–1111. doi: 10.1515/geo-2019-0085

CrossRef Full Text | Google Scholar

Ambraseys, N. N. (1989). Temporary seismic quiescence: SE Turkey. Geophys. J. Int. 96, 311–331. doi: 10.1111/j.1365-246X.1989.tb04453.x

CrossRef Full Text | Google Scholar

Ambraseys, N. N., and Jackson, J. A. (1998). Faulting associated with historical and recent earthquakes in the Eastern Mediterranean region. Geophys. J. Int. 133, 390–406. doi: 10.1046/j.1365-246X.1998.00508.x

CrossRef Full Text | Google Scholar

Arpat, E., and Şaroğlu, F. (1972). The east anatolian fault system: thoughts on its development. Bull. Miner. Res. Explor. Inst. Turkey 78, 33–39.

Google Scholar

Baiesi, M., and Paczuski, M. (2004). Scale-free networks of earthquakes and aftershocks. Phys. Revi. E 69:066106. doi: 10.1103/PhysRevE.69.066106

PubMed Abstract | CrossRef Full Text | Google Scholar

Baiesi, M., and Paczuski, M. (2005). Complex networks of earthquakes and aftershocks. Nonlinear Process. Geophy. 12, 1–11. doi: 10.5194/npg-12-1-2005

CrossRef Full Text | Google Scholar

Basili, R., Kastelic, V., Demircioglu, M. B., Garcia Moreno, D., Nemser, E. S., Petricca, P., et al. (2013). The European Database of Seismogenic Faults (EDSF) Compiled in the Framework of the Project SHARE. Available online at: http://diss.rm.ingv.it/share-edsf/ (accessed September 2020). doi: 10.6092/INGV.IT-SHARE-EDSF

PubMed Abstract | CrossRef Full Text | Google Scholar

Bassin, C., Laske, G., and Masters, G. (2000). The current limits of resolution for surface wave tomography in North America. EOS Trans. AGU 81:F897.

Google Scholar

Bird, P. (2003). An updated digital model of plate boundaries. Geochem. Geophys. Geosyst. 4:1027. doi: 10.1029/2001GC000252

CrossRef Full Text | Google Scholar

Braun, T., Caciagli, M., Carapezza, M. L., Famiani, D., Gattuso, A., Lisi, A., et al. (2018). The seismic sequence of 30th May–9th June 2016 in the geothermal site of Torre Alfina (central Italy) and related variations in soil gas emissions. J. Volcanol. Geothermal Res. 359, 21–36. doi: 10.1016/j.jvolgeores.2018.06.005

CrossRef Full Text | Google Scholar

Büyükakpınar, P., Aktar, M., Petersen, G. M., and Köseoğlu, A. (2021). Orientations of broadband stations of the KOERI seismic network (Turkey) from two independent methods: P- and rayleigh-wave polarization. Seismol. Res. Lett. 92, 1512–1521. doi: 10.1785/0220200362

CrossRef Full Text | Google Scholar

Carder, D. S. (1945). Seismic investigations in the Boulder Dam area, 1940-1944, and the influence of reservoir loading on local earthquake activity. Bull. Seismol. Soc. Am. 35, 175–192.

Google Scholar

Cesca, S., Braun, T., Maccaferri, F., Passarelli, L., Rivalta, E., and Dahm, T. (2013a). Source modelling of the M5–6 Emilia-Romagna, Italy, earthquakes (2012 May 20–29). Geophys. J. Int. 193, 1658–1672. doi: 10.1093/gji/ggt069

CrossRef Full Text | Google Scholar

Cesca, S., Grigoli, F., Heimann, S., González, Á, Buforn, E., Maghsoudi, S., et al. (2014). The 2013 September–October seismic sequence offshore Spain: a case of seismicity triggered by gas injection? Geophys. J. Int. 198, 941–953. doi: 10.1093/gji/ggu172

CrossRef Full Text | Google Scholar

Cesca, S., Letort, J., Razafindrakoto, H. N. T., Heimann, S., Rivalta, E., Isken, M. P., et al. (2020). Drainage of a deep magma reservoir near Mayotte inferred from seismicity and deformation. Nat. Geosci. 13, 87–93. doi: 10.1038/s41561-019-0505-5

CrossRef Full Text | Google Scholar

Cesca, S., Rohr, A., and Dahm, T. (2013b). Discrimination of induced seismicity by full moment tensor inversion and decomposition. J. Seismol. 17, 147–163. doi: 10.1007/s10950-012-9305-8

CrossRef Full Text | Google Scholar

Çetin, H., Laman, M., and Ertunç, A. (2000). Settlement and slaking problems in the world’s fourth largest rock-fill dam, the Ataturk Dam in Turkey. Eng. Geol. 56, 225–242. doi: 10.1016/S0013-7952(99)00049-6

CrossRef Full Text | Google Scholar

Dahm, T., Becker, D., Bischoff, M., Cesca, S., Dost, B., Fritschen, R., et al. (2013). Recommendation for the discrimination of human-related and natural seismicity. J. Seismol. 17, 197–202. doi: 10.1007/s10950-012-9295-6

CrossRef Full Text | Google Scholar

Dahm, T., Cesca, S., Hainzl, S., Braun, T., and Krüger, F. (2015). Discrimination between induced, triggered, and natural earthquakes close to hydrocarbon reservoirs: a probabilistic approach based on the modeling of depletion-induced stress changes and seismological source parameters. J. Geophys. Res. Solid Earth 120, 2491–2509. doi: 10.1002/2014JB011778

CrossRef Full Text | Google Scholar

Dahm, T., Hainzl, S., Becker, D., Bisschoff, M., Cesca, S., Dost, B., et al. (2010). “How to discriminate induced, triggered, and natural seismicity,” in Proceedings of the Workshop Induced seismicity: November 15 - 17, 2010 (Walferdange: Centre Européen de Géodynamique et de Séismologie), 69–76.

Google Scholar

Dahm, T., Heimann, S., Funke, S., Wendt, S., Rappsilber, I., Bindi, D., et al. (2018). Seismicity in the block mountains between Halle and Leipzig, Central Germany: centroid moment tensors, ground motion simulation, and felt intensities of two M ≈3 earthquakes in 2015 and 2017. J. Seismol. 22, 985–1003. doi: 10.1007/s10950-018-9746-9

CrossRef Full Text | Google Scholar

Davies, R., Foulger, G., Bindley, A., and Styles, P. (2013). Induced seismicity and hydraulic fracturing for the recovery of hydrocarbons. Mar. Pet. Geol. 45, 171–185. doi: 10.1016/j.marpetgeo.2013.03.016

CrossRef Full Text | Google Scholar

Deng, K., Zhou, S., Wang, R., Robinson, R., Zhao, C., and Cheng, W. (2010). Evidence that the 2008 Mw 7.9 Wenchuan earthquake could not have been induced by the Zipingpu reservoir. Bull. Seismol. Soc. Am. 100, 2805–2814. doi: 10.1785/0120090222

CrossRef Full Text | Google Scholar

Disaster and Emergency Management Authority (AFAD). (1990). Turkish National Seismic Network [Data set]. Department of Earthquake, Disaster and Emergency Management Authority. Available online at: https://doi.org/10.7914/SN/TU (accessed September 2020).

Google Scholar

Dost, B., Stiphout, A., van Kühn, D., Kortekaas, M., Ruigrok, E., and Heimann, S. (2020). Probabilistic moment tensor inversion for hydrocarbon-induced seismicity in the groningen gas field, the netherlands, Part 2: application. Bull. Seismol. Soc. Am. 110, 2112–2123. doi: 10.1785/0120200076

CrossRef Full Text | Google Scholar

Duman, T. Y., and Emre, Ö. (2013). The east anatolian fault: geometry, segmentation and jog characteristics. Geol. Soc. Lond. Special Publ. 372, 495–529. doi: 10.1144/SP372.14

CrossRef Full Text | Google Scholar

Ellsworth, W. L. (2013). Injection-induced earthquakes. Science 341, 1225942. doi: 10.1126/science.1225942

PubMed Abstract | CrossRef Full Text | Google Scholar

Emre, Ö, Duman, T. Y., Özalp, S., Şaroğlu, F., Olgun, Ş, Elmacı, H., et al. (2018). Active fault database of Turkey. Bull. Earthq. Eng. 16, 3229–3275. doi: 10.1007/s10518-016-0041-2

CrossRef Full Text | Google Scholar

Eyidoğan, H., Geçgel, V., and Pabuçcu, Z. (2010). “Correlation between water level decrease in Atatürk Dam, Turkey and Mw5.0 earthquake on September 3, 2008,” in Proceedings ESC 2010, 6-10 September 2010, Montpellier, 61.

Google Scholar

Foulger, G. R., Wilson, M. P., Gluyas, J. G., Julian, B. R., and Davies, R. J. (2018). Global review of human-induced earthquakes. Earth Sci. Rev. 178, 438–514. doi: 10.1016/j.earscirev.2017.07.008

CrossRef Full Text | Google Scholar

Gaebler, P., Ceranna, L., Nooshiri, N., Barth, A., Cesca, S., Frei, M., et al. (2019). A multi-technology analysis of the 2017 North Korean nuclear test. Solid Earth 10, 59–78. doi: 10.5194/se-10-59-2019

CrossRef Full Text | Google Scholar

Gahalaut, K., and Hassoup, A. (2012). Role of fluids in the earthquake occurrence around Aswan reservoir, Egypt. J. Geophys. Res. 117, 2303. doi: 10.1029/2011JB008796

CrossRef Full Text | Google Scholar

Gaite, B., Ugalde, A., Villaseñor, A., and Blanch, E. (2016). Improving the location of induced earthquakes associated with an underground gas storage in the Gulf of Valencia (Spain). Phys. Earth Planet. Interiors 254, 46–59. doi: 10.1016/j.pepi.2016.03.006

CrossRef Full Text | Google Scholar

Ge, S., Liu, M., Lu, N., Godt, J. W., and Luo, G. (2009). Did the zipingpu reservoir trigger the 2008 wenchuan earthquake? Geophys. Res. Lett. 36:L20315. doi: 10.1029/2009GL040349

CrossRef Full Text | Google Scholar

González, P. J., Tiampo, K. F., Palano, M., Cannavó, F., and Fernández, J. (2012). The 2011 Lorca earthquake slip distribution controlled by groundwater crustal unloading. Nat. Geosci. 5, 821–825. doi: 10.1038/ngeo1610

CrossRef Full Text | Google Scholar

Grigoli, F., Cesca, S., Priolo, E., Rinaldi, A. P., Clinton, J. F., Stabile, T. A., et al. (2017). Current challenges in monitoring, discrimination, and management of induced seismicity related to underground industrial activities: a European perspective. Rev. Geophys. 55, 310–340. doi: 10.1002/2016RG000542

CrossRef Full Text | Google Scholar

Grigoli, F., Cesca, S., Rinaldi, A. P., Manconi, A., López-Comino, J. A., Clinton, J. F., et al. (2018). The November 2017 Mw 5.5 Pohang earthquake: a possible case of induced seismicity in South Korea. Science 360, 1003–1006. doi: 10.1126/science.aat2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Gupta, H. K. (1992). Reservoir-Induced Earthquakes. Amsterdam: Elsevier, 364.

Google Scholar

Gupta, H. K. (2002). A review of recent studies of triggered earthquakes by artificial water reservoirs with special emphasis on earthquakes in Koyna, India. Earth Sci. Rev. 58, 279–310. doi: 10.1016/S0012-8252(02)00063-6

CrossRef Full Text | Google Scholar

Gupta, H. K., and Rastogi, B. K. (1976). Dams and Earthquakes. Amsterdam: Elsevier, 229.

Google Scholar

Hainzl, S., Aggarwal, S. K., Khan, P. K., and Rastogi, B. K. (2015). Monsoon-induced earthquake activity in Talala, Gujarat, India. Geophys. J. Int. 200, 627–637. doi: 10.1093/gji/ggu421

CrossRef Full Text | Google Scholar

Havskov, J., and Ottemöller, L. (1999). SeisAn earthquake analysis software. Seismol. Res. Lett. 70, 532–534. doi: 10.1785/gssrl.70.5.532

CrossRef Full Text | Google Scholar

Heimann, S., Isken, M., Kühn, D., Sudhaus, H., Steinberg, A., Daout, S., et al. (2018). Grond—A Probabilistic Earthquake Source Inversion Framework. V. 1.0. Potsdam: GFZ Data Services, doi: 10.5880/GFZ.2.1.2018.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Heimann, S., Kriegerowski, M., Isken, M., Cesca, S., Daout, S., Grigoli, F., et al. (2017). Pyrocko—An open-source seismology toolbox and library. V. 0.3. Potsdam: GFZ Data Services, doi: 10.5880/GFZ.2.1.2017.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Irmak, T. S., Bulut, I., and Doğan, B. (2020). Regional seismicity in Adıyaman – Samsat (SE Turkey). Boll. Geofisica Teorica ed. J. 61, 451–468.

Google Scholar

Jamalreyhani, M., Büyükakpınar, P., Cesca, S., Dahm, T., Sudhaus, H., Rezapour, M., et al. (2020). Seismicity related to the eastern sector of anatolian escape tectonic: the example of the 24 January 2020 Mw 6.77 Elazığ-Sivrice earthquake. Solid Earth Discussions 1–22. doi: 10.5194/se-2020-55

CrossRef Full Text | Google Scholar

Jamalreyhani, M., Pousse-Beltran, L., Büyükakpınar, P., Cesca, S., Nissen, E., Ghods, A., et al. (2021). The 2019-2020 Khalili (Iran) Earthquake Sequence—Anthropogenic Seismicity in the Zagros Simply Folded Belt? [Preprint]. Earth and Space Science Open Archive. Available online at: http://www.essoar.org/doi/10.1002/essoar.10506454.1 (accessed April 2021).

Google Scholar

Jamalreyhani, M., Rezapour, M., Cesca, S., Heimann, S., Vasyura-Bathke, H., Sudhaus, H., et al. (2019). The 2017 November 12 Mw 7.3 Sarpol-Zahab (Iran-Iraq border region) Earthquake: Source Model, Aftershock Sequence and Earthquakes Triggering, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-759. Available online at: https://doi.org/10.5194/egusphere-egu2020-759 (accessed September 2020).

Google Scholar

Juanes, R., Jha, B., Hager, B. H., Shaw, J. H., Plesch, A., Astiz, L., et al. (2016). Were the May 2012 emilia-romagna earthquakes induced? a coupled flow-geomechanics modeling assessment. Geophys. Res. Lett. 43, 6891–6897. doi: 10.1002/2016GL069284

CrossRef Full Text | Google Scholar

Kandilli Observatory and Earthquake Research Institute, Boğaziçi University (KOERI) (1971). Bogazici University Kandilli Observatory and Earthquake Research Institute [Data set]. Istanbul: International Federation of Digital Seismograph Networks, doi: 10.7914/SN/KO

PubMed Abstract | CrossRef Full Text | Google Scholar

Kartal, R. F., and Kadirioğlu, F. T. (2019). Impact of regional tectonic and water stress on the seismicity in Ataturk Dam Basin: Southeast of Turkey. J. Seismol. 23, 699–714. doi: 10.1007/s10950-019-09830-5

CrossRef Full Text | Google Scholar

Keranen, K. M., Savage, H. M., Abers, G. A., and Cochran, E. S. (2013). Potentially induced earthquakes in Oklahoma, USA: Links between wastewater injection and the 2011 Mw 5.7 earthquake sequence. Geology 41, 699–702. doi: 10.1130/G34045.1

CrossRef Full Text | Google Scholar

Keranen, K. M., Weingarten, M., Abers, G. A., Bekins, B. A., and Ge, S. (2014). Sharp increase in central Oklahoma seismicity since 2008 induced by massive wastewater injection. Science 345, 448–451. doi: 10.1126/science.1255802

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, K.-H., Ree, J.-H., Kim, Y., Kim, S., Kang, S. Y., and Seo, W. (2018). Assessing whether the 2017 Mw 5.4 Pohang earthquake in South Korea was an induced event. Science 360, 1007–1009. doi: 10.1126/science.aat6081

PubMed Abstract | CrossRef Full Text | Google Scholar

Kissling, E., Ellsworth, W. L., Eberhart-Phillips, D., and Kradolfer, U. (1994). Initial reference models in local earthquake tomography. J. Geophys. Res. Solid Earth 99, 19635–19646. doi: 10.1029/93JB03138

CrossRef Full Text | Google Scholar

Kühn, D., Heimann, S., Isken, M. P., Ruigrok, E., and Dost, B. (2020). Probabilistic moment tensor inversion for hydrocarbon-induced seismicity in the groningen gas field, the Netherlands, Part 1: testing. Bull. Seismol. Soc. Am. 110, 2095–2111. doi: 10.1785/0120200099

CrossRef Full Text | Google Scholar

Lizurek, G., Wiszniowski, J., Giang, N. V., Van, D. Q., Dung, L. V., Tung, V. D., et al. (2019). Background seismicity and seismic monitoring in the Lai Chau reservoir area. J. Seismol. 23, 1373–1390. doi: 10.1007/s10950-019-09875-6

CrossRef Full Text | Google Scholar

López-Comino, J. Á, Braun, T., Dahm, T., Cesca, S., and Danesi, S. (2021). On the source parameters and genesis of the 2017, Mw 4 montesano earthquake in the outer border of the Val d’Agri Oilfield (Italy). Front. Earth Sci. 8:617794. doi: 10.3389/feart.2020.617794

CrossRef Full Text | Google Scholar

Manga, M., Wang, C.-Y., and Shirzaei, M. (2016). Increased stream discharge after the 3 September 2016 Mw 5.8 Pawnee, Oklahoma earthquake. Geophys. Res. Lett. 43, 11588–11594. doi: 10.1002/2016GL071268

CrossRef Full Text | Google Scholar

Martínez-Díaz, J. J., Bejar-Pizarro, M., Álvarez-Gómez, J. A., Mancilla, F., de, L., Stich, D., et al. (2012). Tectonic and seismic implications of an intersegment rupture: the damaging May 11th 2011 Mw 5.2 Lorca, Spain, earthquake. Tectonophysics 546–547, 28–37. doi: 10.1016/j.tecto.2012.04.010

CrossRef Full Text | Google Scholar

McClusky, S., Balassanian, S., Barka, A., Demir, C., Ergintav, S., Georgiev, I., et al. (2000). Global positioning system constraints on plate kinematics and dynamics in the eastern mediterranean and caucasus. J. Geophys. Res. Solid Earth 105, 5695–5719. doi: 10.1029/1999JB900351

CrossRef Full Text | Google Scholar

McGarr, A., and Simpson, D. (1997). “Keynote lecture: a broad look at induced and triggered seismicity, “Rockbursts and seismicity in mines”,” in Proceeding of 4th international symposium on rockbursts and seismicity in mines, Poland, 11–14 August 1997, eds S. J. Gibowicz and S. Lasocki (Rotterdam: A.A. Balkema), 385–396.

Google Scholar

McGarr, A., Simpson, D., and Seeber, L. (2002). 40—Case histories of induced and triggered seismicity. Int. Geophys. 81A, 647–661.

Google Scholar

McKenzie, D. (1972). Active tectonics of the mediterranean region. Geophys. J. Int. 30, 109–185. doi: 10.1111/j.1365-246X.1972.tb02351.x

CrossRef Full Text | Google Scholar

MTA (2020). General Directorate of Mineral Research and Exploration. 1:100 000 scaled Geological Maps of Turkey, sheet M40-M41, Geological Research Department. Ankara: MTA.

Google Scholar

Negi, S. S., Paul, A., Cesca, S., Kamal Kriegerowski, M., Mahesh, P., and Gupta, S. (2017). Crustal velocity structure and earthquake processes of Garhwal-Kumaun Himalaya: Constraints from regional waveform inversion and array beam modeling. Tectonophysics 712–713, 45–63.

Google Scholar

Perinçek, D., Günay, Y., and Kozlu, H. (1987). “New observations on strike-slip faults in east and southeast Anatolia”, in Proceedings 7th Biannual Petroleum Congress of Turkey, UCTEA Chamber of Petroleum Engineers, Turkish Association of Petroleum Geologists, Ankara, 89–103.

Google Scholar

Reilinger, R., McClusky, S., Vernant, P., Lawrence, S., Ergintav, S., Cakmak, R., et al. (2006). GPS constraints on continental deformation in the Africa-Arabia-Eurasia continental collision zone and implications for the dynamics of plate interactions. J. Geophys. Res. 111:B05411. doi: 10.1029/2005JB004051

CrossRef Full Text | Google Scholar

Richter, G., Hainzl, S., Dahm, T., and Zöller, G. (2020). Stress-based, statistical modeling of the induced seismicity at the Groningen gas field, The Netherlands. Environ. Earth Sci. 79, 252. doi: 10.1007/s12665-020-08941-4

CrossRef Full Text | Google Scholar

Rinaldi, A. P., Improta, L., Hainzl, S., Catalli, F., Urpi, L., and Wiemer, S. (2020). Combined approach of poroelastic and earthquake nucleation applied to the reservoir-induced seismic activity in the Val d’Agri area, Italy. J. Rock Mechan. Geotech. Eng. 12, 802–810. doi: 10.1016/j.jrmge.2020.04.003

CrossRef Full Text | Google Scholar

Roeloffs, E. A. (1988). Fault stability changes induced beneath a reservoir with cyclic variations in water level. J. Geophys. Res. Solid Earth 93, 2107–2124. doi: 10.1029/JB093iB03p02107

CrossRef Full Text | Google Scholar

Ruiz-Barajas, S., Santoyo, M. A., Benito Oterino, M. B., Alvarado, G. E., and Climent, A. (2019). Stress transfer patterns and local seismicity related to reservoir water-level variations. A case study in central Costa Rica. Sci. Rep. 9:5600. doi: 10.1038/s41598-019-41890-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Şahbaz, N., and Seyitoğlu, G. (2018). The neotectonics of NE Gaziantep: the Bozova and Halfeti strike-slip faults and their relationships with blind thrusts, Turkey. Bull. Miner. Res. Explor. 156, 17–40.

Google Scholar

Scholz, C. H. (2015). On the stress dependence of the earthquake b value. Geophys. Res. Lett. 42, 1399–1402. doi: 10.1002/2014GL062863

CrossRef Full Text | Google Scholar

Schorlemmer, D., Wiemer, S., and Wyss, M. (2005). Variations in earthquake-size distribution across different stress regimes. Nature 437, 539–542. doi: 10.1038/nature04094

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwatke, C., Dettmering, D., Bosch, W., and Seitz, F. (2015). DAHITI – an innovative approach for estimating water level time series over inland waters using multi-mission satellite altimetry. Hydrol. Earth Syst. Sci. 19, 4345–4364. doi: 10.5194/hess-19-4345-2015

CrossRef Full Text | Google Scholar

Şengör, A. M. C., Görür, N., and Şaroğlu, F. (1985). “Strike-slip faulting and related basin formation in zones of tectonic escape: Turkey as a case study,” in Strike-slip Deformation, Basin Formation and Sedimentation, Vol. 37, eds K. T. Biddle and N. Christie- Blick (Society of Economic Mineralogist and Paleontologists Special), 227–264. doi: 10.2110/pec.85.37.0211

PubMed Abstract | CrossRef Full Text | Google Scholar

Şengör, A. M. C., and Yılmaz, Y. (1981). Tethyan evolution of Turkey a plate tectonic approach. Tectonophysics 75, 181–241.

Google Scholar

Shapiro, S. A., Krüger, O. S., and Dinske, C. (2013). Probability of inducing given-magnitude earthquakes by perturbing finite volumes of rocks. J. Geophys. Res. 118, 3557–3575. doi: 10.1002/jgrb.50264

CrossRef Full Text | Google Scholar

Simpson, D. W. (1986). Triggered Earthquakes. Annu. Rev. Earth Planet. Sci. 14, 21–42. doi: 10.1146/annurev.ea.14.050186.000321

CrossRef Full Text | Google Scholar

Simpson, D. W., Leith, W. S., and Scholz, C. H. (1988). Two types of reservoir-induced seismicity. Bull. Seismol. Soc. Am. 78, 2025–2040.

Google Scholar

Soysal, H., Sipahioğlu, S., Kolçak, D., and Altınok, Y. (1981). Earthquake Catalog for Turkey and surrounding area (2100 B.C.–1900 A.D.). The Scientific and Technological Research Council of Turkey, Report Nr: TBAG-341. Ankara: The Scientific and Technical Research Council of Turkey (TUBİTAK), 86.

Google Scholar

Tosun, H. (2012). “Re-analysis of ataturk dam under ground shaking by finite element models,” in Proceedings of the CDA Annual Conference, Saskatoon, SK.

Google Scholar

Tosun, H., Zorluer, İ, Orhan, A., Seyrek, E., Türköz, M., and Savaş, H. (2007). Seismic hazard and total risk analyses for large dams in Euphrates Basin in Turkey. Eng. Geol. 89, 155–170.

Google Scholar

Villaseñor, A., Herrmann, R. B., Gaite, B., and Ugalde, A. (2020). Fault reactivation by gas injection at an underground gas storage off the east coast of Spain. Solid Earth 11, 63–74. doi: 10.5194/se-11-63-2020

CrossRef Full Text | Google Scholar

Walsh, F. R., and Zoback, M. D. (2015). Oklahoma’s recent earthquakes and saltwater disposal. Sci. Adv. 1:e1500195. doi: 10.1126/sciadv.1500195

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilson, M. P., Foulger, G. R., Gluyas, J. G., Davies, R. J., and Julian, B. R. (2017). HiQuake: the human-induced earthquake database. Seismol. Res. Lett. 88, 1560–1565. doi: 10.1785/0220170112

CrossRef Full Text | Google Scholar

Yeck, W. L., Hayes, G. P., McNamara, D. E., Rubinstein, J. L., Barnhart, W. D., Earle, P. S., et al. (2017). Oklahoma experiences largest earthquake during ongoing regional wastewater injection hazard mitigation efforts. Geophys. Res. Lett. 44, 711–717. doi: 10.1002/2016GL071685

CrossRef Full Text | Google Scholar

Yeck, W. L., Weingarten, M., Benz, H. M., McNamara, D. E., Bergman, E. A., Herrmann, R. B., et al. (2016). Far-field pressurization likely caused one of the largest injection induced earthquakes by reactivating a large preexisting basement fault structure. Geophys. Res. Lett. 43, 10198–10207. doi: 10.1002/2016GL070861

CrossRef Full Text | Google Scholar

Zaliapin, I., and Ben-Zion, Y. (2013). Earthquake clusters in southern California I: identification and stability. J. Geophys. Res. Solid Earth 118, 2847–2864. doi: 10.1002/jgrb.50179

CrossRef Full Text | Google Scholar

Zaliapin, I., Gabrielov, A., Keilis-Borok, V., and Wong, H. (2008). Clustering analysis of seismicity and aftershock identification. Phys. Rev. Lett. 101, 018501. doi: 10.1103/PhysRevLett.101.018501

PubMed Abstract | CrossRef Full Text | Google Scholar

Appendix A: Declustering

The method quantifies the correlation between an event i and a preceding event j by its magnitude-weighted space-time distance $n_{ij}=(t_i-t_j){|\overrightarrow{x_i}-\overrightarrow{x_j}|}^d{10}^{-b{M}_j}$with $t,\overrightarrow{x},M$ being the time, location, and magnitude of the events. b is the Gutenberg-Richter b-value, and d is the fractal dimension of the hypocenter distribution, which would be 2 for a planar distribution and 3 for a homogeneous three-dimensional distribution. The distance can be written as n_ij = (T_ijR_ij) with rescaled time T_ij = (t_j−t_i)10^−0.5bM_i and rescaled distance $R_{ij}={|\overrightarrow{x_j}-\overrightarrow{x_i}|}^d{10}^{-0.5b{M}_i}$. Among all events j preceding i, the identification of the (most likely) trigger of i results from selecting that event with the lowest n_ij-value. To distinguish between triggered and background activity, a threshold value of n_c is set, and only events with n_ij≤n_c are considered as plausible mainshock-aftershock pairs. By means of epidemic-type aftershock sequence (ETAS) simulations, the applied detection method has been previously demonstrated to be robust with respect to (1) changes of the involved parameters of the method, (2) catalog incompleteness, and (3) location errors (Zaliapin and Ben-Zion, 2013). Here we used standard parameters d = 2.3,b = 1, and a threshold value of n_c = 10^−3.5.

Appendix B: Reservoir-Induced Stress Changes

To calculate the reservoir-related stress variations, we follow the calculations of Gahalaut and Hassoup (2012), assuming a uniform and isotropic half-space. The total pressure changes p related to a reservoir is the sum of p_c and p_d, which are the change in pore pressure due to the instant compression caused by the reservoir load, and the change in pore pressure due to the diffusion of reservoir water load, respectively (Roeloffs, 1988). The instant effect p_c can be calculated by −B(σ₁₁ + σ₂₂ + σ₃₃)/3, where B is the Skempton coefficient. Additionally, we consider for a given receiver mechanism (strike, dip, rake), the induced shear stress τ and compressional normal stress σ_n related to the loading and finally calculate the corresponding Coulomb Failure Stress CFS = τ−μ(σ_n−p) with friction coefficient μ. For that, we calculate the loading induced stress tensor σ_ij using 3-D Boussinesq-Cerruti solutions for a point force acting on the surface of an infinite half-space (see, e.g., Deng et al., 2010). The pressure change related to diffusion is calculated by the convolution of the observed reservoir water level L with the Green’s function G (Gahalaut and Hassoup, 2012; Hainzl et al., 2015):

$\Delta p_d(x,y,z,t)=D\underset{0}{\overset{t}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}L(\overline{x},\overline{y},\overline{t})G(x-\overline{x},y-\overline{y},z,t-\overline{t})d\overline{t}d\overline{x}d\overline{y}\left(1\right)$

with $G(x-\overline{x},y-\overline{y},z,t-\overline{t})=\frac{z}{8{\mathrm{\pi }}^{\mathrm{\hspace{0.33em}1.5}}\left[D(t-\overline{t})\right]{}^{2.5}}\mathit{exp}(-\frac{(x-\overline{x}){}^2+(y-\overline{y}){}^2+z^2}{4D\left(t-\overline{t}\right)})$

where x,y and $\overline{x},\overline{y}$ refer to horizontal coordinates of the observation and source points, respectively, and z is the depth of the observation. For our calculations, we use the values B = 0.5, μ = 0.8.