Introduction
Plate boundary forces caused by relative movements between Africa, Arabia, and Eurasia resulted in complex tectonic structures in the western Anatolia-Aegean region (Fig. 1). Plate movements involving the convergence of Arabia and Africa with Eurasia in the north; 1) give rise to the westward extrusion of Anatolia along the strike-slip North Anatolian Fault (NAF) and East Anatolian Fault (EAF), 2) cause crustal extension in western Anatolia-Aegean, and 3) are responsible for complex deformation in the Eastern Mediterranean region (McKenzie, 1972; Dewey and Şengör, 1979; Le Pichon and Angelier, 1979; Angelier et al., 1981; Jackson and McKenzie, 1984, 1988; Taymaz et al., 1990; 1991; Over et al., 2010; Özden et al., 2018). The complex tectonics of western Anatolia-Aegean involve different extensional regimes: N-S, NE-SW, and NW-SE (McKenzie, 1972; Jolivet and Brun, 2010; Alçiçek et al., 2006; Shah, 2015; Över et al., 2016; Ocakoğlu et al., 2018).
Gökova Gulf, which runs east-west, is an asymmetrical basin that is growing towards the sea. The Gökova Graben formed on the Lycian Nappes and filled with Plio-Quaternary units known as the Gökova Formation (Görür et al., 1995; Gürer and Yilmaz, 2002; Gürer et al., 2013; Tur et al., 2015). Gürer et al. (2013) suggested that the northern edge of the gulf is bounded by the south-dipping E-W Gökova Fault, one of the most active structures in the southern part of west Anatolia-Aegean, based on geomorphological and geological features. The southern border of the basin is not as clear and uninterrupted as the northern border. The active fault map also shows that the Datça Fault Zone, which represents the western part of the southern boundary of the basin, is more active (Fig. 2). The Datça Fault Zone is divided into two parts: north and south. The northern segment extends to the east of Kos (Karasözen et al., 2018), while the southern segment is listric (Kurt et al., 1999). Analyzing both bathymetric and seismic reflection data, Datça et al. (2013) revealed the existence of active strike-slip faults with various directions in the Gulf of Gökova. From analysis of the focal mechanisms, Shah (2015) asserted that the area is dominated by normal events with few strike-slip events at the western margin of Gökova Gulf. The interpretation of new multichannel seismic profiles by Ocakoğlu et al. (2018) also shows the existence of active strike-slip and normal faults in different directions in the Gulf of Gökova. They pointed out southern-dipping normal faults on the north border (i.e. Gökova Fault Zone) and northern-dipping normal faults on the southern border of the basin (i.e. İşcan Fault Zone). Normal faults were also implied in the Gökova Basin by several studies (e.g., Kurt et al., 1999; Uluğ et al., 2005; içcan et al., 2013; Tur et al., 2015; Ocakoğlu et al., 2018).
Following the earthquake on July 20, 2017, researchers obtained contrary findings about the source fault and the fault's dip direction: (1) from geophysical and geodetic data the Bodrum-Kos earthquake ruptured a south-dipping E-W normal fault (Saltogianni et al., 2017; Tiryakioğlu et al., 2018; Ocakoğlu et al., 2018) or (2) a north-dipping E-W normal fault (Karasözen et al., 2018; Ganas et al., 2019; Konca et al., 2019). Karasözen et al. (2018) asserted that the source fault was the northern Datça Fault. A section of Gökova Bay extending from the Bodrum Peninsula to the Datça Peninsula shows the presence of E-W striking faults dipping both south and north (Fig. 3).
In the last 30 years, GPS studies about the tectonics of the Eastern Mediterranean have yielded important results, including: 1) the velocity and direction of movement of Africa and Arabia relative to Eurasia (McClusky et al., 2000; Reilinger et al., 2006), 2) the rate and direction of movement of Anatolia towards the west relative to Eurasia along the dextral NAF (North Anatolian Fault) and sinistral EAF (East Anatolian Fault) (McClusky et al., 2000; Reilinger et al., 2006), 3) rate and direction of motion in west Anatolia-Aegean toward the Hellenic Trench (McClusky et al., 2000; Aktuğ et al., 2009; Howell et al., 2017), and 4) the crustal extension in the west Anatolia-Aegean (LePichon et al., 1995; McClusky et al., 2000; Reilinger et al., 2010). This area is currently undergoing a N-S extensional process at a rate of 30-40 mm/yr, which is well documented (Oral et al., 1995; LePichon et al., 1995), 6) the determination of both earthquake parameters and the source fault geometry (Saltogianni et al., 2017; Tiryakioğlu et al., 2018; Karasözen et al., 2018), and finally 7) revealed the coseismic deformation of the 20th July Bodrum-Kos earthquake (Tiryakioğlu et al., 2018). The GPS velocity vectors calculated relative to Anatolia (fixed) indicate an increase in the velocity of the southwestward extrusion of Anatolia in western Turkey and the Aegean region compared to other areas of Turkey accompanied by anticlockwise rotation (McClusky et al., 2000; Howell et al., 2017).
To define the stress regimes and their tectonic implications, we used a method proposed by Carey-Gahillardis and Mercier (1987). We introduced the focal mechanisms of both the main shock and its aftershocks into the inversion algorithm. We also compiled the published focal mechanisms of earthquakes that occurred in Gökova and its environs between 1933 and 2017 in order to determine the tectonic regime responsible for the great Bodrum-Kos Earthquake in terms of regional significance.
Seismicity
The historical catalogues show that Gökova Gulf and its surrounding settlements have been hit many times by a number of moderate to major destructive earthquakes, such as in 24 BC, 412 BC (Ambraseys and White, 1997), 227 BC, 199-198 BC, AD 141, 144, 174, 344, 1493, 1851, 1863, and 1869 (Ergin et al., 1967; Guidoboni et al., 1994; Yolsal et al., 2007). Luttrell (1999) reported that the city of Bodrum was totally demolished by the 1493 earthquake. The activity in the region is also confirmed by earthquakes occurring during the instrumental period, such as 23 April 1933 (Mw: 6.2), 23 May 1941 (Mw: 5.4), 13 December 1941 (Mw: 6.3), 25 April 1959 (Mw: 6.1), 19 February 1989 (Mw: 5.6), 5 October 1999 (Mw: 5.2), 4 August 2004 (Mw: 5.4), and 10 January 2005 (Mw: 5.2) (Eyidoğan and Barka, 1996; Çaroglu et al., 1992; Uluğ et al., 2005; Yolsal and Taymaz, 2010). Kalafat and Horasan (2012) observed swarm type activity in the western Gökova Gulf during 2004 and 2005. Historical and instrumental era earthquakes show that the basin contains significant faults that enable its development. On July 20, 2017, one of the faults triggered a 6.6 (Mw) earthquake. Following the Bodrum-Kos Earthquake, more than 3000 aftershocks ranging in magnitude from 1.0 to 4.8 were recorded across Gökova Gulf according to KOERI (Kandilli Observatory and Earthquake Research Institute) data (Fig. 4). Aftershocks were scattered across a wide region, mainly in the west of the basin, and did not follow a linear distribution (Fig. 4). This event created a tsunami with a water wave reaching almost 2 meters height, which damaged the southern coast of Bodrum and the northern coast of Kos Island (Heiderzadeh et al., 2017; Yalçiner et al., 2017). Several studies reported tsunamis in relation to historical and instrumental period events (i.e. in 365, 554, 1303, 1481, 1822 and 1948) in the Eastern Mediterranean (Altinok and Ersoy, 2000; Yolsal et al., 2007; Çevikbilen and Taymaz, 2012, Çevikbilen et al., 2014).
Methodology
Moment Tensor Inversion
In this study we used Moment Tensor (MT) inversion algorithms that model the waveforms recorded at one or more 3-component broadband seismic stations (e.g., Kuge, 2003; Sokos and Zahradnik, 2008). These techniques perform waveform inversions to find source parameters of small to moderate sized earthquakes. One of the most user-friendly MT inversion routines recently adopted by Seiscmop3 is ISOLA, which is based on FORTRAN codes and offers a MATLAB graphic interface. ISOLA allows for both single and multiple point source iterative deconvolution (Kikuchi and Kanamori, 1991) and inversion of complete regional and local waveforms. In addition, we also used the technique developed by Kuge (2003). In this method, waveform fitting between the observed and synthetic displacement seismograms from one or more stations at local distances is achieved by searching for a MT point on a grid scheme for the best fit between the observed and the synthetic displacement seismograms. During the inversion process uniform weight is given to all seismograms. The broadband seismic station networks operated by KOERI, Disaster and Emergency Management Presidency (AFAD) and National Observatory of Athens (NOA) were utilized. For data selection to be used for the inversion process, firstly the signal quality check was completed for three-component broadband stations. Seismograms with gaps and signals with signal-to-noise ratio lower than 4.0 were not considered. Good azimuthal distribution of stations plays significant role in constraining the MT inversion results. The Green's function (GF) calculation is performed using the frequency-wavenumber method (Bouchon, 1981). GFs were computed using the crustal structure from Akyol et al. (2006). When GFs are calculated, the positions of epicenters are fixed, and the depth is allowed to deviate from close to the surface to a depth of up to 30 km in steps of 2 km. The quality of fit between the observed and predicted seismograms is measured by variance reduction (VR); the larger the value of VR, the better the fit. The variance reduction is calculated for various depths for each time shift, and the faulting mechanism was selected with maximum VR for the analyzed event (Fig. 5).
Inversion of seismic slip-vector data sets to determine the present-day stress state
To compute the state of stress responsible for present-day faulting from the population of focal mechanisms of earthquakes that occurred in the Gökova Gulf, we used an inversion method proposed by Carey-Gailhardis and Mercier (1987); this is one of several existing algorithms (Vasseur et al., 1983; Gephart and Forsyth, 1984). The inversion is a statistical method which allows calculation of the best mean fitting stress state from a population of focal mechanisms by selecting one of two nodal planes as the fault plane.
According to this inversion process, the slip (s, defined by a slip vector corresponding to a striation for geological data or a rake for seismological focal mechanisms) on each fault plane occurs in the direction of the resolved shear stress (t), the fault plane being a pre-existing fracture. In the case of a seismic event, the inversion computes a mean best-fitting deviatoric stress tensor by minimizing the angular deviation between an expected slip vector (maximum shear, t) and the observed slip vector (s) deduced from the focal mechanism (Carey and Brunier 1974; Carey 1979). All inversion results include the orientation (azimuth and plunge) ofthe principal stress axes of a mean deviatoric stress tensor, as well as a 'stress ratio' [R = (σ2 - σ 1)/( σ 3 - σ 1)], a linear quantity describing relative stress magnitudes, where the principal stress axes, σ1, σ2, and o3 correspond to compressional, intermediate, and extensional deviatoric stress axes, respectively. It is necessary to know the seismic slip vector, and consequently to select the preferred seismic fault plane for each pair of nodal planes to compute the stress state from earthquake focal mechanisms. For major earthquakes, the selection can be made if there is a co-seismic rupture, or from the spatial epicenter distribution of the aftershock sequence. There is another option for earthquake populations with low magnitude and no surface rupture-computation. It is possible to calculate the fault slip vector using Bott's (1959) model since only one of the two slip vectors of a focal mechanism solution is in accordance with the principal stress axes. For this slip vector, the R ratio, defined [R = (σ2 - σ1)/( σ3 - σ1)], is such that 0<R<1 (Carey-Gailhardis and Mercier, 1987). Furthermore, unless both nodal planes of a focal mechanism converge along a principal stress axis, if one of the nodal planes meets this condition, the other does not (Carey-Gailhardis and Mercier, 1987). We used the calculation method described above to select the nodal plane of each focal mechanism in this analysis. In general, a collection of seismic event focal mechanisms leads to a well-defined assessment of the regional or local stress state that is consistent with the geologically determined stress state: that is, the stress state resulting from inversions of striae measured on fault planes (e.g., Sebrier et al., 1988; Mercier et al., 1991, 1992; Bellier and Zoback, 1995; Bellier et al., 1997; Över et al., 2010). As mentioned above, fault slip inversion schemes assume that the slip vector on each plane corresponds to the direction of the maximum resolved shear stress on that plane. The inversion involves at least four different fault sets in this case since there are four unknowns (three defining the orientation of the principal axes and one defining the stress ratio R). Faults with variable dip angles and distinct strike directions, rather than a continuum of strikes around a single mean direction, are included in ideal data sets. When the deviation angle between the measured slip vector 't' and the observed slip vector 's' is less than 20°, the slip vector determined from a focal mechanism is usually considered to be mechanically explained by the computed stress deviator. Results of stress inversions are considered reliable if 80 percent of the deviation angles between t and s are less than 20° and if the computed solution is stable; that is, the inversion tends towards the same solution regardless of the initial parameter values (Carey, 1979; Carey-Gailhardis and Mercier, 1987; Mercier et al., 1989, 1991; Bellier and Zoback, 1995).
Result and Discussion
Focal Mechanisms of the Main Earthquake and Aftershocks
An earthquake (M=6.6) hit Bodrum, Kos, and the surrounding area in the west of the Gulf of Gökova on July 20, 2017. The main shock and the majority of the aftershocks show predominantly normal faulting character. The focal mechanism analysis was made using Kuge (2003)'s method and ISOLA software. The main shock (Mw: 6.6) with a N-S (N4°E) T-axis must have occurred on a relatively large E-W fault (Table 1). The MT analysis yields normal faulting mechanisms for 34 of the events occurring after the Bodrum-Kos earthquake (Mw 6.6). The focal mechanisms of the aftershocks show two distinct nodal plane directions. They are approximately E-W (Fig. 8 and Table 1) and NE-SW to ENE-WSW (Fig. 9 and Table 2), respectively.
The E-W directional nodal planes indicate N-S extension, while the ones with NE-SW to ENE-WSW orientation indicate approximately NW-SE extension. Several marine geophysical and geodetic studies revealed a large number of faults, mostly normal and a few strike-slip, in the basin with varying directions (Kurt et al., 1999; Uluğ et al., 2005; içcan et al., 2013; Tur et al., 2015; Ocakoğlu et al., 2018; Tiryakioğlu et al., 2018). Ocakoğlu et al. (2018) described NE-SW normal faults at the western end of the bay off Kos Island (their profile 3 in Fig. 11).
N-S extension
In this study 58 events, comprising new and published focal mechanisms (Tables 1, 2 and 3), are included in the Carey-Gailhardis and Mercier (1987) inversion method in order to obtain present-day stress tensors and their tectonic implications. These are the focal mechanisms for earthquakes that occurred in Gökova and its environs before July 20, 2017 (references in Table 3), as well as the major shock and its aftershocks (see Tables 2 and 3). Figure 10 illustrates the available focal mechanisms for earthquakes with magnitudes ranging from 4.0 to 6.6 between 1933 and 2017 including the 20th July 2017 main event (Mw: 6.6). All focal mechanisms that give normal to oblique slip faulting are introduced into the inversion in order to determine the stress state acting in an area wider than Gökova Gulf. The inversion of 24 selected nodal plane sets (Table 3) yields a normal faulting stress regime (e.g., σ1 is vertical) characterized by approximately N-S (N356°E) o3 axis, i.e. minimum horizontal stress axis (Fig. 11). The computed R value is 0.571 indicating a triaxial extensional stress state and clearly differs from a radial extension stress state, which requires a high R value (R≥0.85) as defined in Bellier and Zoback (1995).
Figure 8 shows the focal mechanisms of 12 aftershocks with magnitude ranging from 3.6 to 4.5 as well as that of the main shock. The inversion of the 13 nodal plane sets (Table 1) gives a normal faulting stress regime with approximately N-S (N5°E) minimum horizontal stress axis (o3) (Fig. 11). The determined R value of 0.367 implies a triaxial stress state rather than a radial extension which needs a low R value (R<0.15) as described by Bellier and Zoback (1995). The inversion of the moment tensor of the 20th July 2017 earthquake (Mw 6.6) yields fault parameters from both nodal planes (NP1: strike, 280°, dip. 50° and rake, -79°; NP2: 83°, 41°, -103°; see Table 1 for details). The inversion of the focal mechanism algorithm chose NP1 as the suitable fault plane (i.e. the plane that dips to the north) for calculating the convenient stress tensor. Similar to other inversion methods, the analysis we used in this study is an effective and successful method for determining the stress tensor, or stress field orientations, from a population of nodal planes, but it is not capable of determining stress magnitude. Furthermore, even if the algorithm selects one of the nodal planes as the fault plane based on the computed stress tensor, surface rupture and/or aftershock distribution are required to confirm this. Saltogianni et al. (2017) proposed that the Bodrum-Kos earthquake caused a nearly E-to ESE striking rupture in the upper crust from the surface (sea bed) to a depth of 12 km based on seismological and geodetic data. However, no surface ruptures were observed on land. As a result, it is difficult to tell which of the nodal planes is the true fault (i.e. whether the seismic fault dips north or south) based on the available evidence.
The dip and location of the source fault are also debated; those who argue that the fault dips north indicated that the fault is located in the center of the gulf, such as the North Datça Fault extending to the east of Kos (Karasözen et al., 2018) or a fault with 10 km E-W trend following the Gökova Rift (Ganas et al., 2019). The seismic fault is situated in the northwest of the basin, offshore of Kos and Bodrum according to those who suggested the source fault dips south (Saltogianni et al., 2017; Tiryakioğlu et al. 2018; Ocakoğlu et al., 2018). The fault planes measured near Ören town in the middle part of the Gökova Fault as shown in Figure 12, are E-W striking normal faults with high dip angle (80o). Our MT parameters for the 20 July 2017 earthquake, on the other hand, show E-W nodal planes with relatively low dip (Table 1). The current results indicate that the fault that ruptured during the 20 July 2017 earthquake is not the same as the Gökova Fault or its western continuation in the gulf. The disagreements about the consequences of the 20 July 2017 earthquake (i.e. the inconsistency of the results and interpretations regarding the fault that ruptured during the earthquake and its dip) show that the basin is quite complex and remains to be further explored.
NW-SE extension
The 21 focal mechanisms of the aftershocks, which range in magnitude from 3.8 to 4.8, are given in Figure 9. The focal mechanisms of the earthquakes show normal to oblique slip faulting (Table 2 and Fig. 9). The inversion of 21 nodal plane sets yields a normal faulting stress regime (σ1 is vertical) characterized by approximately NW-SE ((N330°E) o3 axis, i.e. minimum horizontal stress axis (Table 2 and Fig. 11). The computed R value is 0.45 indicating a triaxial extensional stress state. The approximately NW-SE extension was obtained from the inversion of the aftershock focal mechanisms in the western part of the Gulf of Gökova (Fig. 4). The presence of these normal faults is revealed by seismic-reflection profiles with GPS, land DEM and multibeam bathymetry in the marine area (Tur et al. 2015; Ocakoğlu et al., 2018; Tiryakioğlu et al., 2018). Some of the mapped faults were suggested to have strike-slip character as a result of previous geophysical, multibeam bathymetry, and geodetic studies (Uluğ et al., 2005; içcan et al., 2013; Ocakoğlu et al., 2018; Tiryakioğlu et al., 2018). None of the aftershocks in this analysis are strike-slip faulting; instead, they all are normal faulting to oblique slip (Figs. 8 and 9, and Table 1-2). Based on combined geophysical and geodetic data from land and sea, Tur et al. (2015) demonstrated that the faults mapped as strike-slip faults are normal type faults. They also argued that all faults mapped in the gulf which are associated with the tectonic development of the basin have normal character with distinct directions: 1) NW-SE faults controlling submarine valleys, 2) E-W and 3) approximately NE-SW faults associated with the development of Gökova Basin. A seismic profile taken from the eastern coast of Kos Island in the northwest of Gökova Gulf in the NW-SE direction reveals the presence of listric and synthetic faults (see profile no:3 in Ocakoğlu et al., 2018). These faults are in NE-SW direction, suggesting the existence of NW-SE extension, according to this seismic profile. According to Ocakoğlu et al. (2018), there is decelerated tectonism in this area, and minor extension may be enough to create these structures.
Shah (2015) also found the normal extensional regime in western Anatolia, from N-S to NNE-SSW, suddenly changes direction to NW-SE at the exit of Gökova Gulf, using stress tensor analysis of the earthquakes. Shah (2015) inverted 67 nodal plane sets and found a normal faulting stress regime (o1 is vertical) with a NW-SE (N313°E) o3 axis, and a R value of 0.39, which supports our findings (Fig. 11). Shah (2015) suggested also that the Gökova area is distinguished by a clockwise rotation of stress directions from N-S to NW-SE. The anticlockwise rotation in the Aegean was also noted in previous geodetic studies (McClusky et al., 2000; Howell et al., 2017). GPS slip vectors (Saltogianni et al., 2017; Tiryakioğlu et al., 2017; Ganas et al., 2019), which show displacements (i.e. horizontal co-seismic displacement) changing direction from N-S to NW-SE, support this faulting process and agree with our findings.
Conclusions
The extensive work of processing the data generated by the main shock and its aftershock activity is presented here to contribute to understanding of the stress regime which is responsible for the deformation in the study area. The results are as follows:
1) The MT solution of the main shock yields approximately E-W striking nodal planes which are compatible with N-S extension. However, the MT solution analysis of aftershock events gives two distinct directions for the nodal planes: approximately E-W and NE-SW to NNE-SSW. The inversion of the focal mechanisms yields N-S extension from E-W nodal planes and NW-SE extension from NE-SW to NNE-SSW nodal planes. The Bodrum-Kos earthquake (Mw: 6.6) and its aftershocks with E-W focal planes occurred under N-S extension, while the aftershocks related to the approximately NE-SW trending nodal planes were generated under NW-SE extension. The NW-SE extension appears to be local and probably contributes to the growth of the western part of the asymmetric Gökova Basin.
2) The inversion of the focal mechanisms of earthquakes occurring in Gökova Basin and its surrounding areas for the period between 1933 and 2017 yields approximately N-S extension. The N-S extension is also supported by field data. The N-S extension which is responsible for the 20 July 2017 Bodrum-Kos earthquake (6.6 in Mw) seems to be related to African subduction beneath the Aegean.