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

Reservoir-triggered seismicity has been observed near dams during construction, impoundment, and cyclic filling in many parts of the earth. In Turkey, the number of dams has increased substantially over the last decade, with Atatürk Dam being the largest dam in Turkey with a total water capacity of 48.7 billion m3. After the construction of the dam, the monitoring network has improved. Considering earthquakes above the long-term completeness magnitude of MC = 3.5, the local seismicity rate has substantially increased after the filling of the reservoir. Recently, two damaging earthquakes of Mw 5.5 and Mw 5.1 occurred in the town of Samsat near the Atatürk Reservoir in 2017 and 2018, respectively. In this study, we analyze the spatio-temporal evolution of seismicity and its source properties in relation to the temporal water-level variations and the stresses resulting from surface loading and pore-pressure diffusion. We find that water-level and seismicity rate are anti-correlated, which is explained by the stabilization effect of the gravitational induced stress imposed by water loading on the local faults. On the other hand, we find that the overall effective stress in the seismogenic zone increased over decades due to pore-pressure diffusion, explaining the enhanced background seismicity during recent years. Additionally, we observe a progressive decrease of the Gutenberg-Richter b-value. Our results indicate that the stressing rate finally focused on the region where the two damaging earthquakes occurred in 2017 and 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., 2013Dahm et al., , 2015Shapiro 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).
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, 1962M w 6.2, Kariba Zambia-Zimbabwe, 1963M w 6.2, Koyna Dam-India, 1967M w 6.3, and Zipingpu Reservoir-Wenchuan, 2008M w 7.9 (Gupta and Rastogi, 1976Gupta, 1992Gupta, , 2002Ge et al., 2009). Wilson et al. (2017) have recently constructed a database with 186 reported cases of RTS (The Human-Induced Earthquake Database HiQuake) 1 . 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) 2 . 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 3 ) and among the largest dam sites in terms of 1 https://inducedearthquakes.org; last accessed September 2020. 2 https://www.dsi.gov.tr/Sayfa/Detay/754; last accessed April 2021. electricity production (DSI) 3 . It is also the largest clay-cored rockfill 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) 3 . 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 rockfill 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).
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 Seyitoglu, 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 anḑ Saroglu, 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 . 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 Kadirioglu, 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ıg-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) 4 .  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ıg-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.
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 Eyidogan 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 Kadirioglu, 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, Bogaziç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 Minimum water level 513 m asl* Maximum water level 542 m asl* 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 Eyidogan et al. (2010). 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.

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 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).

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 5 , 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
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 . 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).
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.
Focal Depths of the M w 5.5 2017 and M w

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 6 https://pyrocko.org/grond 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-tonoise (SNR) ratio. We rely on the Abedeto algorithm 7 , 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.

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  Figures 6, 7).
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.
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 2 /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 2 /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. 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, 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.

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. 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 2 /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 2 /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 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.
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). 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.

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., Eyidogan et al., 2010;Kartal and Kadirioglu, 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 Eyidogan et al. (2010), by comparing earthquakes between 1992 and 2009 with the water level information obtained by DSI. However, Kartal and Kadirioglu (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 (Eyidogan et al., 2010;Kartal and Kadirioglu, 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.

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.