1. Introduction

The East Anatolian Fault Zone (EAFZ), one of the most important neotectonic and seismogenic fault systems in Turkey, is a sinistral strike-slip active fault zone extending from Karliova in the east of Erzincan to the Hatay region, where it reaches the Dead Sea Fault Zone further south [Arpat and Saroglu 1972; Westaway and Arger 1996]. The EAFZ has intense seismic activity depending on the northward relative motions of Arabian and African plates, the westward motion of the Anatolian plate, and continental collision occurring along the Bitlis–Zagros Fold and Thrust Belt [Saroglu and Yilmaz 1990].

In terms of seismicity, the Karasu Fault (KF; Amanous Fault) segment of the EAFZ and Toprakkale, Karatas, and Yumurtalık Faults are considerable fault segments within the complex structure of the northeast Mediterranean region. These faults were evaluated as active faults on the Active Fault Map of Turkey updated recently [Emre et al. 2018]. On the other hand, the fact that few geodetic studies conducted to investigate the kinematics of these active faults makes it difficult to understand the structure of these faults and analyze earthquake hazards [Aktuğ et al. 2016; Mahmoud et al. 2012; Nocquet 2012; Meghraoui et al. 2011; Reilinger et al. 2006]. To estimate the potential hazard of a probable earthquake, there is a need to study the slip rates along these active faults and also fault locking parameters using geodetic methods. However, the use of permanent Global Positioning System (GPS) stations together with campaign GPS observations is the most efficient method for this purpose, but the geometry and distinctive features of relevant faults are important criteria for geodetic networks.

This study was conceived with the aim of monitoring active faults with a high probability of generating earthquakes around the southern end of the EAFZ. Accordingly, we intend to estimate the recent slip rates along these faults and fault locking parameters with high accuracy and spatial resolution by modeling the observed GPS data. For this purpose, campaign observation sites were used along profiles perpendicular to the related faults. A wide geodetic network consisting of campaign observation sites together with permanent GPS stations was observed during the execution of the above proposed study. Raw data obtained from previous studies related to geodesy, paleoseismology, and active tectonics were analyzed together with the latest GPS observations to estimate current strain accumulations for the Eastern Mediterranean region. By analyzing all the obtained data, slip rates along the faults and fault locking depths were studied spatially and temporally in detail.

During the past decades, strain accumulations on the EAFZ have sometimes caused life-threatening earthquakes in eastern Turkey such as Elazig–Karakocan (Mw 6.1, 2010) and Bingol (Mw 6.4, 2003) earthquakes. The last destructive earthquake on the EAFZ occurred with an epicenter close to the town of Sivrice in Elazig province on the Puturge segment of the EAFZ in January 2020. The earthquake solution published by the United States Geological Survey (USGS) is a magnitude of Mw 6.7 and a depth of 10.0 km. However, seismic hazard analysis may prevent substantial loss of life and property associated with such relatively major earthquakes. In this context, the potential moment magnitude for a possible earthquake within the modeled segment of an active fault can be predicted using empirical relationships. These involve either total displacement considering historical and instrumental seismic data [Wells and Coppersmith 1994] or fault slip rate together with the length of fault rupture [Anderson et al. 1996]. Since the predicted moment magnitude is converted into average displacement with the help of the approach proposed by Wells and Coppersmith [1994], the recurrence interval of an earthquake can be calculated. From this perspective, the determination of tectonic deformations in active fault zones within the Adana, Osmaniye, Hatay, and Gaziantep provinces is very crucial in terms of drawing attention to seismic hazards in the Eastern Mediterranean region.

2. Kinematic setting

Since the motions of tectonic plates in the Earth’s lithosphere are defined in accordance with plate tectonics theory, it is believed that the Arabian, Anatolian, and Sinai plates interact around the region of the Eastern Mediterranean. According to this hypothesis, the Arabian plate moves to the north at an approximate rate of 20 mm/yr and the Sinai plate converges to a subduction zone along the Cyprus Arc at approximately 10 mm/yr. Furthermore, the westward motion of the Anatolian plate is accommodated by a sinistral strike slip along the EAFZ and a dextral strike slip along the North Anatolian Fault Zone at rates of 10 mm/yr and 20–25 mm/yr, respectively [Reilinger et al. 1997; McClusky et al. 2000; Reilinger et al. 2006; Mahmoud et al. 2012; Aktuğ et al. 2016].

Considering the role of geodetic techniques in active tectonics, it is now possible to estimate the slip rates more precisely along active faults. The slip rates for EAFZ segments are estimated in different geodetic studies as 15 ± 3 mm∕yr [Reilinger et al. 1997] for the whole EAFZ and 10.0 ± 0.3 mm∕yr for the Palu–Hazar Lake segment and 9.9 ± 0.2 mm∕yr for the Hazar Lake–Turkoglu segment along the EAFZ [Reilinger et al. 2006]. The slip rates are calculated to be 9.7 ± 0.9 mm∕yr for EAFZ and 5.5 ± 1.5 mm∕yr for the Karatas–Osmaniye Fault (KOF) segment [Bertrand 2006].

More recently, one of the GPS block modeling studies about slip rate estimations has reported the results as 8.8 ± 0.3 mm∕yr for the EAFZ and 3.6 ± 0.6 mm∕yr for the KOF segment [Meghraoui et al. 2011]. However, Aktuğ et al. [2016] obtained slip rates ranging from 10.3 ± 0.7 mm∕yr to 13.5 ± 1.3 mm∕yr for the segments between Turkoglu and Karliova. Additionally, the rates were indicated to be 4.5 ± 1.1 mm∕yr for the KF and 2.7 ± 1.4 mm∕yr along the KOF in the same study.

Tectonic processes in the Eastern Mediterranean region trigger seismic activities around the area. This is why active fault zones generating seismic hazards need to be monitored continuously using geodetic techniques (Figure 1). Some well-documented destructive historical seismic activities indicate potential risks in seismic zones, namely between the Turkoglu–Golbasi segment of the EAFZ and the KOF (1114 AD earthquake, Ms > 7.8), northern segment of the Dead Sea Fault (DSF) (1408 AD earthquake, M 7.4), KOF (1513 AD earthquake, M ∼ 7.4), and KF (1822 AD earthquake, M 7.0 and 1872 AD earthquake, M 7.2) [Ambraseys and Jackson 1998; Sbeinati et al. 2005; Ambraseys 2004].

3. GPS data set and analysis

Using the GPS network presented in Section 2, campaign GPS observations were carried out at 24 campaign sites for four campaign sets in 2009, 2010, 2011, and 2019. Since 4 of the 24 campaign sites were destroyed, the last campaign in 2019 was achieved with the remaining 20 sites. The strategy adopted during campaign GPS observations is based on data logging in two sessions on consecutive days. In addition, each session covers an observation span of at least 8 hours. All campaign observations were performed within the same time periods of years to mitigate seasonal effects on GPS observations.

The GPS data analysis was performed in the GAMIT/GLOBK software [King et al. 2009] using validated strategies [Reilinger et al. 2006] for both campaign GPS observations and for also data obtained from municipalities, Turkish National Permanent GNSS Network-Active, and IGS permanent GNSS stations. First, in the GAMIT stage, station coordinates, zenith delay of the atmosphere at each station, and orbital and Earth orientation parameters (EOPs) were estimated from daily GPS phase observations using a weighted least squares algorithm. Second, in the GLOBK stage, Kalman filtering was applied to loosely constrained estimates of station coordinates, orbits, and EOPs and their covariances for a reference frame definition with a consistent set of coordinates and velocities. Eighteen stations selected from the IGS global networks were incorporated both to tie our local network with a global network and to estimate the orbital and earth rotation parameters with high accuracy. Apart from IGS stations, a total of 30 permanent GPS stations from local municipalities and TUSAGA-Aktif networks were introduced into the daily GAMIT processing. Loosely constrained daily solutions from GAMIT processing were investigated on expected normalized root mean square (nrms) values between 0.15 and 0.25. After discarding outliers, coordinate repeatabilities were generated by controlling weighted root mean square (wrms) parameters approximately 2–4 mm for horizontal components and 10–15 mm for vertical components. Since the monumentation of GPS stations fluctuates [Langbein and Johnson 1997], a random walk of 2mm∕√yr in horizontal positions of the stations was applied to obtain more realistic uncertainties for GPS velocities [McClusky et al. 2000].

Following the procedures for the definition of the reference frame for both coordinates and velocities [Aktuğ et al. 2009], transformation parameters to ITRF2008 coordinates of selected IGS stations were estimated for the stabilization frame. The stabilization was realized for the Eurasia-fixed reference frame by minimizing the adjustments to the horizontal velocities of the 17 stations where these velocities were obtained with a post root mean square (rms) value of 0.31 mm/yr.

The GPS velocity field obtained from this study was combined individually with published velocities while introducing the velocities by Reilinger et al. [2006] as the reference solution. The transformation accuracy of the velocity field solutions was expressed by the wrms parameter depending on the common observation sites (Table 1). The combined velocity field obtained from the transformation model was referenced to the ITRF2000 Eurasia-fixed reference frame (Figure 2).

4. Kinematic modeling

Kinematic modeling relying on GPS velocity fields instead of interpreting site velocities individually strengthens tectonic analysis in the region. Blocks bounded by faults have different characteristics inside the mass and along their boundaries. Although the block motions are systematic, the resistance due to the friction of rocks prevents free movement at their boundaries. This velocity difference, known as “slip deficit”, causes strain accumulations generally ending in an earthquake. Slip deficit describes the divergence of surface velocities from expected velocities of the blocks depending on geological factors. The block modeling approach on the basis of inversion of GPS velocities provides the estimation of block rotation rates, fault locking parameters, and internal strain rates in an elastic and homogeneous half-space [Okada 1985]. The components of the strain-rate tensor aim to reveal internal deformations within the blocks that are caused by unsettled faults. Apart from GPS site velocities, block and fault geometries formed in a three-dimensional (3D) elastic half-space, euler poles for blocks and locking depths are considered as model inputs for parameter estimation in the inversion. However, defects in block modeling emerge from poor density of observation sites. In this case, slip rate estimations become more susceptible to introduced block geometries and a priori locking depths. The TDEFNODE software [McCaffrey 2009] based on the simulated annealing method [Press et al. 1986] was used for block modeling by minimizing the residuals between observed and model velocities. According to the approach in McCaffrey [2002], minimizing the misfit of the velocities allows the simultaneous estimation of block rotations using Euler poles and coupling fractions along block boundaries. The data misfit minimization is achieved by reduced chi-square statistic tests.

The locking depth in a seismogenic zone roughly reflects the characteristic mechanism of the fault. The state of being locked or creeping for a fault segment is controlled by coupling fraction. This parameter denoted by phi (𝜑) depends on the short-term slip rate (V_c) and the plate velocity (V ), where 𝜑 = 0 corresponds to creeping and 𝜑 = 1 to a fully locked fault while partially locked faults have values between 0 and 1 as per (1).

The kinematic model consists of creeping, fully or partially locked segments [Wang et al. 2003], while the segments in the tectonic structure accommodate purely strike-slip or along with dip-slip components in the Eastern Mediterranean region [Yavaşoğlu 2009; Tiryakioğlu 2012]. The kinematic model in our study is assumed to be locked from the surface down to a locking depth fixed at 7 km in the crust. This approach constrains the fault from creeping since a thicker zone greater than 7 km in depth is recommended when exploring continental transform faults [Genrich et al. 2000; Le Béon et al. 2008; Walters et al. 2011]. In a deeper layer below the locked zone, the fault is assumed to be locked partially down to the depth below where full creeping begins. In our model, the faults were introduced as not creeping with a constraint of 𝜑 = 1 from 0 to 7 km in the zone close to the surface. Following this, 𝜑 parameters were estimated between 0 and 1 in the inversion for depths ranging from 7 km to the fully creeping limit. Below the partially locked zone, the faults creep fully with 𝜑 = 0. Between locked and fully creeping zones, which is called the “effective transition zone” by Wang et al. [2003], the down dip transition was estimated using a factor of 𝜆 = 0.2 based on an exponential function, which indicates the shape of the slip distribution. During the inversion process, along-strike smoothing was introduced using a penalty function with a smoothing factor referred to above. No special data weighting strategy was followed since all the velocity solutions in modeling were equally weighted. However, the velocity field published by Mahmoud et al. [2012] was not incorporated into the kinematic model since the site velocity uncertainties on average are greater than those in each independent solution.

This study was initiated from a point of view comprising several blocks and main faults, namely Arabia (arab-reference block), Anatolia (anat), and Sinai (sina) blocks and East Anatolian Fault (EAF), KF, Cyprus Arc (CA), and DSF. The east of the Amanous (E-AMNS) block boundary equates to the KF, while the western boundary of the Amanous block was assumed to be free slipping (nonlocked) without any elastic strain accumulated. As discussed in Mahmoud [2012] and Mahmoud et al. [2012], block and fault geometries in the kinematic model vary from the simplest to the most complex related to the tectonic structures in the region. The different kinematic models from Mahmoud [2012] were considered. However, evaluating the statistical results of the inversion, the geometry of the proposed kinematic model fits better with observations. The simplest model with the blocks of arab, sina, and anat comprises EAF, KF, CA, and DSF. The reduced chi-square parameter of this simplest model was estimated as 4.816 with a wrms of 1.476 mm and an nrms of 2.124. The more complex model consists of one additional block and one fault segment compared to the simplest model. The addition of the Iskenderun Block (iskd) and the KOF improves the reduced chi-square parameter to 3.869 with a wrms of 1.122 and an nrms of 1.943. On the other hand, the Amanous block (amns) was introduced additionally into the model in the last stage. Eventually, the GPS observations fit well with the kinematic model formed with five blocks and five fault segments (Figure 3). The significant nonsystematic residuals are noted with 1.35 mm overall data rms between the model and observations (Figure 4). Depending on the block and fault geometries introduced, the reduced chi-square parameter of the model converges to 1.814.

In light of the findings from this research, a left-lateral strike slip at a rate of 7.5 mm/yr on the EAF allows significant stress transfer to the KF at 4.4–5.4 mm/yr sinistral strike-slip rates (Figure 5). Moreover, it contributes to a sinistral strike slip at a rate of 3.4 mm/yr and a reverse slip at a rate of 3.1 mm/yr along the KOF. Reverse faulting on the KOF indicates contraction in agreement with the result by Aktuğ et al. [2016]. Normal faulting on the KF has a slip at a rate approximately 3 mm/yr indicating extension, which is consistent with the result by Reilinger et al. [2006]. The internal block strain-rate tensors, which are shown by blue double-headed arrows in the same figure, clearly point out a compression of approximately 30 nanostrain/year within the Amanous Block. However, the greatest strain accumulation was loaded on the Iskenderun Block at approximately 40 nanostrain/year, inducing extension parallel to the main faults and compression perpendicular to them. It is also possible to see a similar pattern on the Sina Block with a strain rate almost half of the Iskenderun Block. The slowest straining region seems to be the Anatolian Block with an internal deformation of approximately 10 nanostrain/year. Furthermore, a locking depth of less than 10 km was estimated for the EAF. However, it is partially locked down to a 25 km depth. In contrast, just like the KOF, the KF has a locking depth of 20 km (Figure 6).

5. Discussion and conclusion

In this study, we explored the strain accumulations on major tectonic features around the Eastern Mediterranean region by the inversion of GPS velocity estimations to obtain block rotation rates, fault locking parameters, and internal strain rates in a 3D elastic dislocation model. For this purpose, a region-specific GPS network, composed of campaign observation sites and permanent stations, was formed to clarify kinematic characteristics of the tectonically active zones around the Hatay Triple Junction in the region. Initially, former campaign data sets from 2009, 2010, and 2011 acquired in joint campaign observation sites from previous studies were collected. The campaign data sets from previous years provide a historical kinematic background for the region. Moreover, it was required to determine recent tectonic deformations. Ongoing strain accumulations lead this research to review the estimations for slip distributions and locking depths using new GPS observations across major faults. Therefore, GPS observations performed in 2019 on the same network, which was the fourth time after 2009, 2010, and 2011, have great importance for understanding the present kinematics of the Eastern Mediterranean region. The campaign GPS observations in 2019 helped improve the data quality in estimating velocity uncertainties. In addition to the new campaign data set acquired in 2019, permanent GPS stations from local networks operated by municipalities were merged into our regional network to densify observation sites spatially and temporally. Both the new campaign data set and permanent GPS stations, regarded as the true value of this study, constitute unique data sets that were not available in the past.

Accordingly, slip rates for the faults mentioned above were published in several studies recently. Earlier, Aktuğ et al. [2016] and Mahmoud et al. [2012] estimated 10.5 mm/yr and 9.0 mm/yr strike-slip rates, respectively, while Reilinger et al. [2006] reported 10.0 mm/yr slip rates for the left-lateral EAFZ. However, there seems to be no consensus on normal/reverse fault components of the EAF since Reilinger et al. [2006] argued for extension at a 5.1 mm/yr dip-slip rate, while Aktuğ et al. [2016] proposed contraction at 2.4–6.3 mm/yr reverse-slip rates. On the other hand, Mahmoud et al. [2012] and our study assert no significant normal/reverse fault components for the EAF. Similarly, the KF was also studied many times recently so that the estimations for strike-slip rates are 4.5 mm/yr by Aktuğ et al. [2016], 4.0 mm/yr by Mahmoud et al. [2012], 6.8 mm/yr by Reilinger et al. [2006], and 4.4–5.4 mm/yr from this study. There are also some debates about normal/reverse faulting on the KF. In contrast, the KOF needs to be studied further in detail.

In conclusion, a sinistral slip rate of 7.5 mm/yr estimated from our model was initialized to accumulate a strain of approximately 6.8 m in total after the 1114 AD earthquake (Ms > 7.8) between the Turkoglu–Golbasi segment of the EAFZ and the KOF. The total strain accumulation within the Golbasi–Turkoglu segment supports the prediction of the next probable major earthquake with a magnitude of 7.2–7.6 if it ruptures entirely over its full 90 km length. Similarly, a 3.4 mm/yr dextral slip rate for the KOF corresponds to a 1.7 m total strain since the 1513 AD earthquake (M ∼ 7.4) and makes M 7.1 possible for the next probable major earthquake. However, Mw 6.8–7.0 is predicted in the case of two segments with a rupture length of approximately 40 km in each. Considering the 1822 AD (M 7.0) and 1872 AD (M 7.2) earthquakes on the KF, a total strain accumulation of 0.65–0.80 m at a 4.4–5.4 mm/yr sinistral slip rate confirms the prediction of a possible M 6.8–6.9 earthquake along 150 km based on Wells and Coppersmith [1994]. As per the approach by Anderson et al. [1996], another prediction is a magnitude of Mw 7.2 if one of the 75 km segments ruptures along the KF.