2.1. Data
The first arrival times of the P- and S-waves of the 2500 earthquake records of the magnitude ≥ 2.3 between latitudes of 36.5° and 37.5° and longitudes of 26.75° and 29° that occurred during 2017 were determined using zSACwin software (Yılmazer 2003) with manual picking. The earthquake records were chosen from the 30 stations, including 18 strong-motion accelerometers and 4 weak-motion seismometer stations of the Disaster and Emergency Management Presidency (AFAD), 5 weak-motion seismometer stations of the Boğaziçi University Kandilli Observatory and Earthquake Research Institute (KOERI), and 3 weak-motion seismometer stations of the Geodynamic Institute and the National Observatory of Athens (NOA). (Fig. 2, Table A.1).
2.2. Minimum one-dimensional (1-D) starting velocity model
In the study area, 11-layer model (Fig. 3, Table A.2) obtained by deriving intermediate velocity values using linear interpolation method from the 7-layer 1-D model proposed by Doğanay-Özkan (2020) was used as the initial velocity model (Özdemirli-Esat 2020). The 1488 earthquakes, whose first arrival times were determined from the routine catalog data used in the 1-D inverse solution in the region, were selected according to an azimuthal gap smaller than 200° and 10 P-wave and 5 S-wave arrivals (Doğanay-Özkan 2020) in terms of data quality. After selecting the initial model and the earthquakes to be used, it was determined in which areas in the region to provide better resolution by generating a synthetic dataset.
2.3. Double-difference (DD) tomography
DD tomography is the generalized version of the DD location. It solves simultaneously for 3-D velocity structure and seismic event locations. DD tomography uses a combination of absolute and accurate differential arrival times and estimates the velocity structure from large to a smaller scale hierarchically. This method can generate velocity structure by identifying incident locations closer to the source region, and then the standard tomography uses only absolute arrival times (Zhang and Thurber 2006). It is based on the hypoDD (Waldhauser and Ellsworth 2000) code and uses both absolute and relative arrival time data. The method determines a 3-D velocity model with absolute and relative event locations. The advantage of this approach is to include relative arrival times in quality values along with absolute arrival times, thereby simplifying assumptions about ray path geometries or path anomalies without losing valuable information, using absolute selections, but also determining absolute positions (Zhang and Thurber 2003).
The velocity model obtained by DD tomography is better than standard tomography. With standard tomography, the event location is uncertain and slightly scattered due to errors. The use of differential arrival times in DD tomography eliminates most of these errors, which removes some uncertainties from the velocity model (Zhang and Thurber 2003).
2.4. Inversion with synthetic and real data
We applied the checkerboard analysis as a solution analysis. In the synthetic data, forward modeling was created for the medium where low and high velocities are 4 km/s and 6 km/s, respectively and Vp/Vs ratio is 1.73. In the inversion with the initial model, where velocities are defined as 5 km/s in all directions, the aim is to determine the regions where the resolution is high or low in the tomography to performed with real data. Synthetic arrival times were generated using the same station geometry, focal positions, weightings, and ray propagation, and then used as input data in the 3-D inversion. While generating synthetic data, noise with zero mean and 0.1s standard deviation suitable for Gaussian distribution was added. Several models with different grid spacings were produced, such as 20 km x 20 km (Fig. A.1), but the model having 10 km x 10 km grid spacing was chosen to be used in this study (Fig. A.2). Synthetic P-wave velocity perturbation model was obtained using synthetic data inversion. Because of the inversion of P-wave synthetic travel times, the checkerboard model with 10 km x 10 km node interval was generated (Fig. A.2). In this checkerboard model, high resolution is seen in the western part of the Gulf of Gökova where the 2017 main shock and its aftershocks are located. The model shows high resolution at 4 km depth in the Datça Peninsula south of the gulf and between 6 and 15 km depths in the Bodrum Peninsula and its surroundings. Apart from these specific areas, high resolutions are also seen in various parts of the gulf.
Accurate determination of damping and smoothing values used in the inverse solution affects velocity model results and resolution calculations (Kissling et al 2001). An appropriate damping value ensures smooth variation of the velocity model and minimization of data variance (Eberhart-Philips 1986). The trade-off curve was obtained by performing a single iteration in the inverse solution, with values between 200 and 1000 for the damping value and 0.5-5 for the smoothing value (Fig. 4). The most appropriate damping factor was chosen as 600, the smoothing value was chosen as 2 in the x, y, and z directions, and the curve formed by the damping values for the most appropriate smoothing value is given in Fig. 4. As seen in the trade-off curve, the most suitable damping value is the value that makes the condition number (CND), defined as the ratio of the largest to smallest eigenvalues, 60. To solve the damped least squares algorithm, different damping values should be chosen to make the CNDs similar (Zhang and Thurber 2006). The CND has no meaning during the inversion steps (DeShon et al 2007; Hansen et al 2013). Approximately 28,000 P- and 22,000 S-waves, absolute data from the first arrival times, 1,080,967 differential catalog P, and 688,769 differential catalog S data produced in ph2dt program (Waldhauser and Ellsworth 2000) were used as inputs in the tomoDD (Zhang and Thurber 2003). We did not use cross-correlation data.
The differential and absolute data files, hierarchical weighting scheme, and initial model were determined, and inversion was started in the tomoDD. Weightings, damping, threshold, and iteration numbers have been found through trial and error by varying until obtaining the best resolution.
While determining the hierarchical weighting scheme (Table A.3) used in the tomoDD, the relative weightings of the absolute and differential data were controlled. First, by starting the inversion with higher weighting of the absolute data, broader-scale velocity heterogeneity was obtained, and then higher weighting of the differential data resulted in better resolution near the sources.
According to Zhang and Thurber (2006), when a higher weighting (0.01*abs + diff) is given to the differential data, a better solution is obtained using high resolution values in earthquake locations. For the slowness parameters at the model edges, the system has high model resolution with higher weighting (abs + 0.1*diff) applied to the absolute data (the “abs” and “diff” refer to absolute and differential data, respectively).
When the inversion process is terminated in the tomoDD, derivative weighted sum (DWS) values are generated and written to a file along with the velocity models created after each iteration. DWS is the ray sampling density used to qualitatively define whether the model is well resolved. Since it is difficult to calculate the model resolution matrix for large seismic tomography problems, it is used to characterize the model resolution qualitatively (Zhang and Thurber 2007). DWS ≥ 40 was chosen in the inversion where the absolute data were higher weighting, and DWS ≥ 100 was selected where the differential data were higher weighting (Table A.3). The parameterization used in the synthetic data is also used in the inversion of the real data. 1442 earthquakes were relocated (Fig. 2, Fig. 5, Fig. 6). The 3-D velocity models of P- and S-waves (Fig. A.3, Fig. A.4) and Vp/Vs ratio (Fig. A.5) were generated up to 30 km. Vp/Vs was visualized by directly dividing the absolute P-wave velocities by absolute S-wave velocities obtained from the 3-D inversion. The models were illustrated with the eight north-south (Fig. A.6) and five east-west profiles (Fig. A.7). The north-south and east-west profiles of the DWS values are also given in Fig. A.8.
We observed that the resolution of up to 30 km depth is good around the Bodrum Peninsula where the earthquakes clustered and west of the gulf. However, information cannot be obtained from depths below 30 km because the earthquakes are mostly shallow.
2.5. Interpretation of the tomographic results
There are two different anomaly regions for the seismic velocity structure in the volume with the best resolution (Fig. A.6). The first velocity structure is a sedimentary basin, which is a water-saturated, fractured, and porous structure where both P- and S-wave velocities are quite low and Vp/Vs ratio is high. The basin is observed up to approximately 4–5 km depth. The second velocity structure beneath the basin area shows that the P-wave velocity is low, the S-wave velocity is relatively high, and the Vp/Vs ratio is low. If the system is steam-dominant, it exhibits a behavior similar to the dry state, and a low Vp/Vs anomaly is observed because the gas can be compressed higher than salt water. A sudden decrease in the P-wave velocity while the S-wave velocity decreases slowly causes a low Vp/Vs anomaly. In the case of high pore pressure, the gas has a negligible bulk modulus and acts as a saturated fluid (Mavko et al 2003). Tomographic studies conducted in geothermal fields show that vapor-filled volumes are associated with low Vp/Vs anomalies (Julian et al 1996). Accordingly, it can be interpreted that the active fracture zones in the study area may contain gas. Additionally, the regions out of the fractured zones and showing such anomaly were interpreted as hot magma intrusions (Fig. 7, Fig. 8).
Considering the visual tomographic data, we concluded that two different faults are active in the west of the Gulf of Gökova and the Bodrum Peninsula, considering the main shock of Mw6.5 and the earthquake data of 2017 used in the study. Tectonic interpretation of the north-south tomographic profile#2 (Fig. 7) demonstrates the DKBF and the younger GFZ. GFZ can be shown in two parts considering its morphological features at the surface. Karaada Islet is located between these two parts. The faults observed in the north and south of the islet are most probably merged at depth. The north-dipping fault immediately north of the GFZ is also seen in the profile. Only the main shock and nearby earthquakes were projected onto the profile. In other words, the earthquake cluster in the eastern part of the gulf not projected to avoid confusion in the interpretation. The distribution of the earthquakes on the profile is also compatible with the faults.
In profile#2 (Fig. 7), where the geometry of the DKBF is clearly observed, the basement rock and sedimentary basin separation can be seen. A young fault in the northernmost part cuts the DKBF again. In the southernmost part, we think that the low velocity zone up to 15 km depth is a magmatic intrusion.
In the north-south tomographic profile#3 (Fig. 8), we interpret that the younger south-dipping normal fault zone GFZ cuts the older low-angle breakaway fault DKBF and reaches a depth of around 20 km. The geometry of DKBF and positions of the antithetic and synthetic faults in the Gulf of Gökova with the basin geometry shown in the seismic reflection profiles in Kurt (2000) support the results of our study. Additionally, the northern continuation of DKBF, which outcropped at the north of Milas (Fig. 1), can be seen in the tomographic profile. Immediately north of the GFZ, two significant north-dipping and south-dipping normal faults are observed in the profiles. These faults are also geomorphologically evident. The earthquake activity is seen on the north-dipping fault extending to a depth of around 10 km. With the GFZ, these two young faults cut the older DKBF.
The east-west tomographic profile#11 (Fig. 9) is nearly parallel to the strike of the GFZ, especially in the western part of the gulf, where the earthquakes are concentrated. We interpret the area indicated by the dashed lines in the low-velocity zone, where the profile#11 cuts the zone, as the fault plane of the GFZ (Fig. 9).