Kinematics of the Kahramanmaraş triple junction and of Cyprus: evidence of shear partitioning

The recent tectonics of the Eastern Mediterranean result from the interaction of the Arabia, Africa, and Eurasia plates since Miocene (McKenzie, 1972; McKenzie et al., 1970; Biddle, 1985; Fig. 1a). During the middle Miocene, Arabia was separated from Africa along the left-lateral Levant (or Dead Sea) fault zone (e.g., Le Pichon and Francheteau, 1978). The subsequent collision of Arabia with Europe, resulting in gravitational potential build-up in Eastern Anatolia, combined with an acceleration of slab rollback in the Hellenic subduction where the northern African slab is subducted (Brun et al., 2016) and a possible contribution from underlying asthenospheric flow, have been driving the westward extrusion of the Anatolian Plate (Le Pichon and Kreemer, 2010; Özeren and Holt, 2010). Currently, the boundary between Arabia and Anatolia is a left lateral transform plate boundary, the East Anatolian Fault zone (EAF) (McKenzie, 1976). Africa has been divided in the Eastern Mediterranean into a Nubia plate (Le Pichon and Gaulier, 1988; McKenzie et al., 1970) and a Sinai sub-plate (Mahmoud et al., 2005). The Sinai plate thus subducts beneath Anatolia along the Cyprus Arc and moves southward concerning Arabia along the Levant Fault zone. The Kahramanmaraş triple junction is the junction of the Levant Fault zone, the East Anatolian Fault zone (EAF), and the Cyprus Arc Subduction (Karig and Kozlu, 1990; Sengör et al., 1980, 1985) (Fig. 1b). This triple junction also gives rise to a very complex pattern of deformation in the northeastern Mediterranean because of the convergence of the two continental plates, namely Arabia and Anatolia. Additionally, the marine part of the Sinai Plate is formed of thinned continental crust, comprising a thicker region, the Eratosthenes seamount, which is currently impinging the Cyprus Arc (Le Pichon, 2019). In this incipient collision, the relative motion of Sinai and Anatolia is distributed between the Cyprus Arc subduction, the Kyrenia fault running along the Northern Coast of Cyprus, and the Taurus range onshore further north. In fact, two triple junctions are currently active in the area, namely those of Hatay (see Fig. 1b) where the eastward prolongation of the Cyprus arc reaches the Levant Fault, and Kahramanmaraş where the NE prolongation of the Kyrenia fault reaches the EAF (Sengör et al., 2019; Özkan et al., 2023). However, the motion of Cyprus with respect to Anatolia has been ignored so far in regional kinematic models because of a lack of constraints from geodesy, due to insufficient data on the island (see for example Mahmoud et al., 2005; Reilinger et al., 2006; Gomez et al., 2020).

This paper presents a refined GNSS (Global Navigation Satellite Systems) velocity field encompassing southern Turkey and Cyprus. Leveraging these data, we construct regional kinematic models that not only illuminate the motion of the study area but also incorporate higher-resolution boundary conditions and internal deformation characteristics specific to Cyprus. Through these novel data and models, we foster a more comprehensive understanding of the neotectonics of the northeastern Mediterranean. In addition, we will discuss the kinematic context of the major earthquakes that occurred in the Kahramanmaraş triple junction area on February 6, 2023: Mw (moment magnitude) 7.8 Pazarck earthquake on the EAF and Mw 7.6 Elbistan earthquake on Çardak Fault (Barbot et al., 2023; Hussain et al., 2023).

Previous geodetic studies either measured Arabia and Anatolia plate motion and deformation (Bletery et al., 2020; Cavalié and Jónsson, 2014; England et al., 2016; Kurt et al., 2022; Özeren and Holt, 2010; Reilinger et al., 2006; Viltres et al., 2022; Weiss et al., 2020) or focused on the kinematics of the Levant fault zone (Al Tarazi et al., 2011; Alchalbi et al., 2010; Gomez et al., 2020; Gomez et al., 2007; Hamiel and Piatibratova, 2021; Le Beon et al., 2008; Sadeh et al., 2012) and of the East Anatolian fault zone (Bletery et al., 2020; Cavalié and Jónsson, 2014; Walters et al., 2014). However, only a few studies account for the active deformation of secondary faults such as Malatya and Çardak Faults (also called Sürgü Fault by some authors) (Westaway, 2003), or the Karatas-Osmaniye Fault (Mahmoud et al., 2013; Özkan et al., 2023). Ascertaining the kinematics of the Cyprus Arc Subduction using land-based GNSS observations has proven to be a difficult task since much of the subduction arc is expressed beneath the Mediterranean Sea. The Island of Cyprus is the only place within the Eastern Mediterranean Basin where the Cyprus Arc Subduction kinematics can be studied using onshore geodetic constraints. Still, the data have been too sparse to date. We thus conducted new GNSS surveys in Cyprus between 2019 and 2021, along with data from permanent GNSS sites, providing comprehensive spatial coverage for constructing a kinematic model of the eastern Mediterranean. Our velocity field includes follow-up surveys in the Turkish mainland and velocities published by previous studies. In Section 2 we present our GNSS processing strategy and integration workflow details.

We applied two kinematic inversion methods with our combined velocity field as input. On the one hand, we employed a continuum velocity field interpolation method to calculate the strain rate field of the study area (Beavan and Haines, 2001; Haines and Holt, 1993). This strain distribution may be compared with seismicity distribution and tectonic strain regimes indicated by fault maps. This information also contributes to the definition of boundaries for a block model. The block model calculates rigid block motions and coupling on the block boundaries defined as dislocation sources (McCaffrey et al., 2007). The output of this model may thus be interpreted in terms of long-term slip rates and seismic coupling on major faults. This model also allows for the internal deformation of the blocks by calculating a unique and uniform strain rate tensor for each block. In parallel,

Our newly acquired data and model results help us to address several important questions regarding the kinematics and active tectonics of the Eastern Mediterranean region. We quantify the partitioning of deformation among the Cyprus Arc Subduction, Kyrenia fault, and Taurus range. We evaluate the distribution of deformation around the Kahramanmaraş and Hatay triple junctions, notably between the EAF and the Karatas-Osmaniye-Çardak-Malatya fault system and discuss implications for earthquake recurrence intervals. Regarding the Levant fault, the present study is based on fewer velocity vectors than, for instance, presented in (Gomez et al., 2020) and has little to add to their demonstration that slip on the main fault strand is decreasing northward toward the triple junction as part of the motion is diverted offshore.

We present a GNSS velocity field that unites newly derived with previously published velocities. We conducted GNSS surveys in Cyprus between 2019 and 2021 and revisited 18 points that had been previously measured in 1998 and 2001. We also incorporated the data from seven permanent GNSS sites in the southern part of the island. For the first time, we now have reasonable spatial coverage of space geodetic data in Cyprus. This enables us to construct a kinematic model of the easternmost Mediterranean that, in turn, provides us with a more detailed picture of the deformation within Cyprus. We also conducted some follow-up surveys at several GNSS survey-mode sites in the Turkish mainland to further constrain the kinematics of Anatolia with better coverage. We then integrated into our velocity field previously published velocities acquired over the vicinity of the Levant fault and EAF.

The raw data of the continuous stations were obtained from both international networks (International GNSS Service − IGS hereafter) and regional networks (Turkey Continuous GNSS Network and Cyprus Positioning System) (Fig. 1b). We primarily analyzed the dataset from 2009 to 2021 for continuous stations, although some of them had data gaps during that period. However, we have approximately 10 years of time series for almost every continuous site. The survey mode GNSS sites were selected from the Turkey Fundamental GNSS Network to utilize valuable existing observations and data resources. Each observation was carried out using dual-frequency receivers and a filtering cut-off angle of 10 degrees to minimize atmospheric noise. Each survey mode site has at least 7 different sessions, except for 3 sites in the northern part of Cyprus that were set up in 2019, and all sessions have at least 8 hours of observations. The readers can find further details of the data span of survey mode sites in the supplementary file. We evaluated the raw data of 137 GNSS stations (65 continuous − 72 survey modes) and estimated their velocities.

The long-term seismicity catalogue between 1905 and 2019 represented in this figure was compiled from Kandilli Observatory and Earthquake Research Institute’s seismicity catalogue between 1905 and 2019 (Kandilli et al., 2001; http://www.koeri.boun.edu.tr). The magnitude of completeness of the catalogue is ∼Mc = 4. The mean horizontal location uncertainty is less than 5 km in N-S and E-W directions. The mean of the depth uncertainty is ∼3.5 km, varying between 2 and 8 km. We filtered the original catalogue based on quality factors such as horizontal location uncertainty <5 km and RMS <0.5 s. This dataset does not include the February 2023 earthquakes and their aftershocks. In the East Anatolian shear zone, this data set emphasizes seismic activity on the Malatya Fault (MF), between the Çardak Fault (CF) and Karatas-Osmaniye Fault (KOF), and south of Hatay but displays relatively little activity along and on the recently ruptured segments of Çardak fault and the EAF. A concentration of earthquakes observed south of Cyprus is largely associated with the subduction plane of the Cyprus Arc and will here be used to constrain its geometry.

We performed data processing using a combination of GAMIT/GLOBK software (Herring et al., 2018a, and the extensive literature cited therein) and a stochastic approach. GAMIT/GLOBK integrates the double differences method and carrier phase combinations to eliminate geometric and non-dispersive delays in the solution. We used IGS final orbit and clock products as orbit parameters and the VMF1 (Vienna Mapping Function 1) mapping function to minimize the effect of the tropospheric delay (Boehm et al., 2006). In addition, we incorporated over 20 IGS stations into our network to define a well-constrained global network. We processed daily data from 2009 to 2022, while for the period between 1998 and 2009, we only evaluated days with observations for our survey-mode sites. We verified our daily solutions by following the steps outlined in (Herring et al., 1990). Once daily solutions were obtained for each station, combinations were carried out using a Global Kalman filter approach (Herring et al., 2018a). This approach facilitates the sequential estimation of parameters, providing an advantage over other estimation methods due to its ability to define the state vector and its stochastic model for future measurement times (Herring et al., 2018b). Global network solutions from various institutions were integrated into the evaluation to create a comprehensive network encompassing the existing regional one. The time series for all sites were generated relative to the International Terrestrial Reference Frame 2014 (ITRF14) (Altamimi et al., 2016).

The time series analysis was carried out in three main steps. Firstly, outliers were detected and removed, with a particular emphasis on permanent sites due to their extensive data. Secondly, efforts were made to acquire more realistic sigma values for each site. To achieve this, a first-order Gauss-Markov Extrapolation was implemented. Given the substantial differences in data quantity between continuous and survey-mode sites, these steps were executed with multiple approaches. Random walk noise of 0.02 mm/yr was added to all permanent sites’ horizontal and vertical components. Conversely, the random walk noise added to survey-mode sites was five times greater. The subsequent step involved generating velocities for each site. The velocity field was determined with respect to the Arabian plate fixed reference frame, utilizing Euler pole parameters from (Altamimi et al., 2017).

Station velocities and uncertainties for each component (North, East) and time span The supplementary document provides the corresponding time series in Tables S1–S3.

After generating an initial velocity field, we combined our resulting GNSS velocity field with published velocities from previous studies (Gomez et al., 2020; Hamiel and Piatibratova, 2021; Kurt et al., 2022; Viltres et al., 2022; Özkan et al., 2023). To minimize the effect of some well-known sources of noise such as those stemming from different data evaluation strategies and pre-defined reference frames from different studies, we rotated all velocity fields individually with respect to our dataset.

The rotation is based on a least-square approach that aims to optimize the transformation matrix of common stations for each velocity field pair. Though this approach has been applied in several studies, we made some critical changes to the weight matrix of the objective function (see Ozbey et al., 2021, (Eq. (4)). The weight parameter r_i has been constructed as the function of both distance D_i between the ith common site pair, and the number of observations n_x and n_y for the two related stations x and y (Eq. (1)).

Stations with a distance closer than 5 km are considered to be co-located while stations with a distance closer than 1 km are considered to be the same points. In addition, our approach takes into account the plate and block boundaries during its decision-making process. If the related station pair is located on different blocks or plates, the algorithm rejects it. The second parameter that is taken into account is the number of epochs for each site. Here it is important to note that the number of observations of permanent sites has been postulated as 365 for a year. Each velocity field has been rotated separately by taking the velocity field obtained by this study as the reference system. The statistical outcomes of these processes are listed in Table 1.

Figure 2 shows the final velocity field leveraging in the kinematic models. The unified velocities with respect to both the Arabian-fixed and ITRF reference frames, defined from (Altamimi et al. 2017), can be found in the supplementary material as Table S4.

We present two different modelling approaches to reveal the present-day kinematics around the Kahramanmaraş triple junction and Cyprus. We first introduce a continuum kinematic model and generate a continuum velocity field to monitor the deformation of the region. We then introduce a block model that describes the rates of interseismic block motions occurring along the block boundaries. We suggest a block geometry for the study area and testing this geometry with previously published models including our own (Klein et al., 2022) that we developed prior to the two devastating earthquakes of February 2023. Once the best-fitting geometry is determined, we present rotation poles for each block, slip rates, and coupling ratios for the faults. It is important to note that the block model geometry presented herein was defined before 2023 (Klein et al., 2022).

Here, we generate a contemporary strain rate field to characterize deformation styles in the region that comprises northeast Nubia, the eastern Mediterranean Sea, the Levant fault, Cyprus, Adana/Cilicia basin, and neighbouring southern Turkey. We aim to shed light on the kinematics of the region therein. This can be useful for future dynamical models as kinematic constraints or for future seismic hazard models as geodetically inferred moment constraints. Our kinematic continuum model is based on the method described by (Beavan and Haines, 2001; Haines and Holt, 1993).

The method is essentially a least-squares fit to the GNSS data. The horizontal velocity field in the interpolation domain is derived from a vector function W(r):

The function Wr, in equation (2), is defined at the knotpoints of a quadrilaterals mesh on the spherical earth surface and interpolated with bicubic spline functions (Beavan and Haines, 2001; Haines and Holt, 1993). The W(r) values at knotpoints are inverted in order to minimize a penalty function, which in our application case is of the form:

where V_α,β represent the data variance-covariance matrices for the geodetic velocity measurements v^obs with subscripts αβ ranging over longitude φ and latitude θ. The are strain rate tensor components for each cell, S corresponds to the surface area of the related cell.The first double summation of the penalty function is the misfit to the observed GNSS velocity field subject to observational errors. The second double summation represents an a priori constraint to minimize strain. The weighting factor ν determines the relative weight of velocity data and the minimal strain assumption in the penalty function, and thus the amount of smoothing in the interpolation.

Here, we solely utilized the GNSS velocity field without imposing any plate motion boundary conditions. We also assigned a uniform ν value to achieve the objective function. Experiments were made using lower ν values in the grid cells in SW Cyprus and further offshore, where seismic catalogues show the clustering of earthquakes. If these are “damage” zones, their bulk deformability might be higher than other zones. However, these experiments did not significantly improve the fit to the GNSS velocities even in the near field sites in southern Cyprus justifying our decision not to include laterally varying ν values.

Figure 3 shows the second invariant of the strain rate tensor (obeying the formula where are the tensor components) overlain by the seismicity of the region. The solution indicates, in general, low strain rates within Cyprus, where the second invariant rarely exceeds 20 nanostrain/yr. However, in Figure 4, deformation styles indicate a clear spatial variability of deformation in the island. A transpressive strike-slip regime is found along the Kyrenia range, while roughly N-S compression dominates in the southern part of the island (Fig. 4). On a larger scale, this strain partitioning within Cyprus seems to act as a diffuse transition that rotates the predominant compression from NW-SE in the Sinai block onto NE-SW in the Cilicia basin immediately to the north of Cyprus. Between Cyprus and Turkey, the shortening integrated along the principal strain rate axis between the coasts of Cyprus and Turkey amounts to a maximum of 0.8 mm/yr. The lack of GNSS data probably leads us to a strain rate field that is much smoother than reality in the Cilicia basin. Despite this, the solution shows a progressive transition from compression to strike-slip toward the NE, associated with a rotation of principle axes to N-S compression and E-W extension, eventually matching the dominant strain regime found on land in the Adana Basin.

A swath of higher strain rate (more than 30 nanostrain/yr) over a width of 50–60 km is found east of the Adana basin and extends NE along the northern side of the EAF, thus defining a broader East Anatolian shear zone (Fig. 3). The principal strain directions (E-W extension and N-S compression) are consistent with left-lateral strike-slip motion. The areal strain (⁠⁠) is positive (see Fig. 4) except for at a few locations, indicating transtensive to extensional deformation, consistent with focal mechanisms (Fig. 3). Principal strain rate orientations retain the same orientation over a broader area, with lower strain rates, that include the Adana basin and the mountains north of it (Aladağlar, see Fig. 5). Areal strain indicates extension is dominant in these mountains while both mildly transpressive and transtensive styles are found in the Adana basin. This may suggest that gravity influences strain distribution between topographic highs and the basin. The westward limit of this zone of E-W extension coincides with the Ecemiş Fault (see Fig. 4). This fault zone has taken up 60 km of left lateral slip since late Eocene and has been under transtension since Miocene (Sarikaya et al., 2015; Jaffey and Robertson, 2005; Umhoefer et al., 2020; Yildirim et al., 2016). The principal strain rate axes rotate to NE-SW compression and NW-SE extension west of the Ecemiş fault, which are respectively parallel and perpendicular to this part of the Taurus mountain range. The Ecemiş fault thus appears to bound a zone of east-west extension related to the escape of the Anatolian plate. Strain rates within the Taurus range, west of the Ecemiş fault, are low (less than 10 nanostrain/year), and areal strain there is dominantly positive, but changes sign toward the coast in the south-western part of the range. This is the only part of Taurus where compressive strain is currently observed.

Along the Levant fault, the dominating principal strain orientations (see Fig. 4) are consistent with left-lateral shear on the fault. The GNSS coverage, however, provides poor kinematic constraints in the region from south of Turkey along the coast towards Israel and, in particular, along the east side of the Levant fault where station coverage is sparse. Some short wavelength variability of the strain rate field, with compressional axes trending largely NW-SE, is evident within Israel, where the GNSS coverage is dense but mostly located on the western side of the Levant fault. The strain rates become less coherent towards the southern tip of Israel. North of Israel the zones of higher extension (positive areal strain) and compression (negative areal strain) do not match the location of the Lebanon restraining bend. This puzzling observation has been reported previously (cf. Fig. 8 in Gomez et al., 2020) and was possibly explained by the transfer of the compression onto offshore faults. Overall, zones of higher strain rate appear to roughly correlate with zones of higher seismic activity in the 0–30 km depth range corresponding to crustal seismicity (Fig. 3). We already mentioned that the recently ruptured segments of Çardak Fault and the EAF have been relatively silent before the earthquakes. On the other hand, the zone of E-W extension north of Karatas-Osmaniye Fault (KOF) and east of Çardak Fault (CF) displays relatively high seismic activity. A cluster of seismic activity between Mersin and Bolkardağ also appears as a zone of relatively high strain rate in the GNSS interpolation with positive areal strain indicating extension (see Figs. 3 and 4). However, the orientation and style of the strain tensors determined in this area lack consistency and velocities present relatively higher interpolation residuals (Fig. 5). On the other hand, distributed compressive to transpressive deformation between northern Cyprus and the Anatolian coast may explain seismicity beneath the Cilicia basin. Clusters of seismic activity are also present within the Arabian plate and offshore Lebanon in zones of apparently low strain rates, but these areas on the edges of GNSS data coverage are not well constrained in the continuum model interpolation.

Our second modelling effort aims to determine the kinematic behaviour of the Nubia-Cyprus-Anatolia tectonic system in the context of an elastic block-based approach (McCaffrey et al., 2007). Such a block model involves solving an inverse problem where the unknowns are Euler vectors for individual blocks, uniform strain rate tensor for each block, and coupling ratios on the fault node points. The model can, in principle, be constrained by GNSS velocities, geological fault slip rates and the azimuth of these rates, and earthquake focal mechanisms. Here we only used GNSS data.

The block geometry follows the main active fault zones (Levant fault, Cyprus Arc subduction, Kyrenia fault, East Anatolian fault), which are thought to be critical in shaping the regional tectonics (Fig. 6a). The distribution of seismicity is another key feature defining some block boundaries. The main differences in block architecture with previous studies result from the definition of a Cyprus block. On one hand, large-scale studies considered Cyprus to be a part of Anatolia and did not feature the Kyrenia range as a block boundary (Gomez et al., 2020; Reilinger et al., 2006). This is justifiable as it would not have been possible to constrain the motion of a Cyprus block with the very limited GNSS data available on the island at that time. On the other hand, a detailed kinematic study of the Hatay triple junction proposed connecting the Kyrenia arc to the EAF along the KOF (Özkan et al., 2023) and this implies considering the area between Hatay and KOF as part of the Cyprus block. Inversion results obtained with this block geometry are presented in the appendix and we will show that this assumption results in large misfits. An alternate solution may be drawn taking into account a concentration of mostly extensional focal mechanisms along an N-S trend to the east of the Adana/Cilicia basin, 30 km west of the Levant fault. We assume this trend delineates a block boundary, separating Anatolia and the Cyprus block on its western side from a zone of complex deformation along the Levant and East Anatolian faults (the East Anatolian Shear Zone, Fig. 3). As internal block deformation is taken into account in the model (approximated as a uniform strain rate field in each block) this deforming zone can be defined as a block that stretches eastward along the East Anatolian fault toward Malatya (Klein et al., 2022). Defining part of the northern boundary of this block along Çardak fault appears as an obvious hypothesis after the February 6 earthquakes. Moreover, Çardak (sometimes also called Surgu Fault) and Malatya faults (see Fig. 6a) are known to be active and have been included in previous tectonic models of the triple junction (Acarel et al., 2019; Sançar et al., 2019; Westaway, 2003).

To the south of Cyprus lies the Sinai block which is largely offshore. Its motion is crucial for the kinematics of the Cyprus Arc subduction but can only be constrained by velocities along the Levant coast and in the Sinai Peninsula. However, GNSS velocity fields in the Levant and in Southern Sinai do not fit in the same rigid block reference frame. Previous studies of Levant Fault kinematics proposed that the NE part of the Sinai microplate is fragmented in order to account for geodetic slip rates decreasing northward on the Levant fault (Gomez et al., 2020). We thus consider an additional block, referred to as the “Latakia” block, to extend along the Levant coast (dashed blue line in Fig. 6a) north of Israel and compare solutions with and without this additional block. Note that (Gomez et al., 2020) considers several blocks west of the Levant Fault. Hence, our Latakia block should be considered, like the Malatya block, as a deformable block. Moreover, south Sinai may be affected by extension around the Gulf of Aqaba but this extension cannot be well accounted for in our model, which simplifies the prolongation of the Levant fault into the Gulf of Aqaba as a vertical fault. We thus removed all GNSS rates south of 30° N for the south Sinai.

Along the block boundaries, we defined 5 main dislocation sources, capable of accumulating elastic deformation, on which the coupling ratio will be calculated by inversion (Fig. 6a). The Levant and East Anatolian faults, which obey nearly pure strike-slip motion, are modelled as vertical planar sources. The boundary between the Anatolia and Malatya blocks is simplified as a vertical fault. However, to generate the geometry of dipping faults, such as the Cyprus subduction and Kyrenia fault, we followed published interpretations of seismic profiles (Aksu et al., 2014a; Aksu et al., 2005; Aksu et al., 2021; Burton-Ferguson et al., 2005; Calon et al., 2005; Feld et al., 2017; Hall et al., (2005); Welford et al., 2015). In addition, we utilized the seismicity of the region to compile earthquakes greater than Mw = 2.8. The earthquake locations validate the geometry of the main subduction seismogenic zone between Sinai and Cyprus down to a depth of 40 km but do not help define the geometry at the depth of Kyrenia fault between Cyprus and Anatolia (Fig. 6b).

We conducted a series of synthetic tests to determine the optimal spatial resolution of node points for the slip rate distribution, which was determined by the spatial coverage of the GNSS data along both strike and dip directions. To achieve this, we utilized a checkerboard test, in which we divided the main thrust interface south of Cyprus into planar cells, and monitored the level of recovery of the given slip rate boundary conditions using synthetic GNSS velocities at the same geographic locations as our data. We tested two different average cell sizes for the thrust interface: one with an average cell size of 35 km², and the other with an average cell size of 10 km². The tests were carried out without any synthetic observation noises. Our results indicated that the test conducted with larger patches had slightly better misfits than the one consisting of finer patch resolution. Figure S1 in the supplementary displays the checkerboard test solutions.

We evaluate inversion results using three different block geometries. It is crucial to emphasize that the models were inverted employing the same velocity field. In Model 1, we employed the geometry introduced by (Özkan et al., 2023) around the triple junction. This model does not include a Malatya block between Anatolia and the East Anatolian Fault (EAF). Additionally, the Karata s-Osmaniye fault (see Fig. 7a) connects northeastward with the EAF, extending the boundary between Anatolia and Cyprus, and adding part of Hatay to the Cyprus block. Model 3 is the block model illustrated in Figure 6a. It includes a Latakia block in an attempt to account for the fragmentation of the Sinai plate (Gomez et al., 2020). Model 2 is based on Model 3, without a Latakia block, the extent of which is here considered part of the Sinai plate. As indicated in Table 2, the residuals exhibit a substantial improvement in Models 2 and 3, which consider a Malatya block. This notable improvement is primarily attributed to larger misfits in Model 1 for stations in the Hatay region and Cyprus when these stations are forced to be part of the same block.

However, the improvement associated with the assumption of a Latakia block is subtle. Model 3 has a slightly lower misfit than Model 2 but with a small reduction of the chi-2 over the degree of freedom ratio (Tab. 2). An F-test, which is a powerful statistical method that allows us to assess whether the variability of the variances is significantly different, indicates this improvement is not significant (Tab. 2). Figures S2 and S3 in the supplementary show the first 2 model solutions.

Determining robust Euler pole parameters is of utmost importance for addressing the horizontal motion of the region. This problem is complicated by the fact that a large part of the Sinai plate is underwater and that the prolongation of the Sinai plate along the Levant coast (where GNSS stations are located) may be deforming. Euler poles of Sinai relative to Anatolia and to Arabia were obtained by the block model inversion. In order to estimate the Euler pole of the Nubia plate with respect to Anatolia we combined our determination of the Arabia-Anatolia pole with a Nubia-Arabia pole calculated in the ITRF No Net Rotation reference frame by (Altamimi et al., 2017). The Nubia poles from this study, and (Reilinger et al., 2006) are similar to each other, albeit with a slight difference in the rotation rates. The location of the Euler pole of (Reilinger et al., 2006) for Sinai relative to Anatolia and ours are also close to each other. The rotation rate of the (Reilinger et al., 2006) pole, however, is markedly faster resulting in a slower subduction velocity in our model. The data we used to constrain the motion of the Sinai block are essentially the same as in the previous studies, but we excluded data from South Sinai (below 30 N) as these cannot be fit in the same rigid reference frame, and this explains in large part the differences. The pole we determined provides a better fit of GNSS data along the Levant coast, but a worse fit of the GNSS data in the southern part of Sinai. We believe that this pole provides a better description of the motion of the Mediterranean seafloor as it subducts beneath Cyprus (Tab. 3).

After an Euler pole and a uniform strain rate tensor are estimated for each block, a series of nonlinear inversions (i.e., grid search and simulated annealing (Press et al., 2007)) has been run iteratively to solve the coupling ratio on each node point. The Green’s function that coincides with the location of the GNSS stations on the surface is determined with a rigorous approach to discretize the planar fault into rectangular patches (Okada, 1992). For the parametrization of fault coupling, we express the coupling ratio as constant between the surface and depth z1 (an inversion parameter at each node) and decaying exponentially below z2 where the fault starts fully slipping, as proposed by (Wang et al., 2003).

Figure 8a represents the block motions along the boundaries defined as dislocation sources in Model 3. As the residuals we obtained are relatively small, we present histogram of both the north and east residuals in Figure 8b. On the Levant fault, the slip rate decreases steadily from south to north. It accommodates a 4.7 ± 0.6 mm/yr slip rate from the Gulf of Aqaba to the Dead Sea with an almost purely left lateral strike-slip regime. At the Lebanese restraining bend, slip remains dominantly strike-slip at a 3.0 ± 6 mm/yr rate, with a poorly constrained 0.8 ± 0.7 mm/yr compressional rate. The Cyprus arc accommodates 3.5–6.2 mm/yr convergent rates reducing progressively from west to east. The motion of the Kyrenia fault, on the other hand, is mainly left lateral strike-slip from the northwesternmost tip to the east of the island with rates of 3.2–4.2 mm/yr. Although the slip rate along the western prolongation of the Kyrenia fault also obeys a left lateral strike-slip behaviour, the spatial distribution of our dataset may not be considered capable of resolving this particular region. Along the boundary between the Anatolian and Cyprus blocks on one side and the so-called Malatya block on the other side, the inversion indicates a significant ∼1.3 mm/yr extensional motion, which is in agreement with the predominantly extensional style of the earthquakes (see Fig. 6a). Meanwhile, there are relatively larger velocity residuals on the left side of the boundary, roughly coinciding with Adana Basin, that may represent active deformation not properly modelled with the assumed block geometry. Where the block boundary changes its azimuth from N-S to E-W and along the Çardak fault, the slip rate is an almost left lateral strike-slip behaviour with a 2.0 ± 1.0 mm/yr rate. The Malatya fault, which extends from the eastern tip of the Çardak fault towards the north, has also left-lateral strike-slip motion accounting for 1.7 ± 1.2 mm/yr. The motion on the East Anatolian fault zone, on the other hand, decreases from 6.2  ± 1.2 to 5.2 ± 1 mm/yr from the northernmost tip of the EAF to Kahramanmaraş triple junction where it connects with the Levant fault. Some additional shear is taken up by internal deformation of the Malatya block. Furthermore, the block motions for the two alternative model scenarios are provided in the supplementary document (Figs. S2 and S3).

Most block boundaries inverted for dislocation sources appear fully locked down the depth z1. The only notable exception is the western part of the Kyrenia fault in Cyprus, but this result may not be reliable because of the distribution of velocity data, all located south of the fault with few points near the fault zones. The estimated z1 value for the Levant fault is consistently around 10 km. However, the inversion result indicates a depth of approximately 15 km (following the convention of (Wang et al., 2003)) for z2 at the segments north of the Dead Sea, while it is 20 km at the segments south of it. For the EAF, the z1 value decreases from 20 to 13 km from south to north, while z2 remains relatively constant at 25 km. The locking behaviour of the Çardak fault is homogeneous along strike with z1 = 8 km and z2 = 20 km (see Fig. 9a). On the Cyprus arc subduction the inversion estimates a z1 value around 20–25 km depth, but due to the lack of GNSS data for the offshore part, the uncertainties of the coupling coefficient increase with the distance between the node points and Cyprus (Fig. 9b-c).

We discuss our solutions and describe some implications of them including (1) the shear partitioning between the Cyprus subduction and Kyrenia fault; (2) the effect of this partitioning system around the East Anatolian shear zone toward the Kahramanmaraş triple junction; (3) a comparison of our solution with the previous studies around the Levant fault and the prolongation of distributed deformation for the northern part of the Levant fault; and (4) the results of our interseismic fault coupling model and coseismic slip behaviour of the 2023 Kahramanmaraş earthquakes.

The slip characteristics on both the Cyprus arc and Kyrenia range indicate that ongoing shear partitioning is the dominant regime for the area. The subduction is still active despite an incipient collision with the Eratosthenes seamount southwest of Cyprus, with near frontal convergence at a 6.0 mm/yr rate. The continuum deformation field (Fig. 4) obtained in and around the Taurus mountain range shows little shortening in the upper crust and, offshore, convergence between Taurus and Cyprus across the western part of the Cilicia basin is less than 1 mm/yr. The dominantly left-lateral strike-slip motion we find on the Kyrenia fault with a rate of 3.5–4.2 mm/yr, suggests a nearly perfect shear partitioning between the subduction and Kyrenia range. However, the distribution of seismicity in the area suggests the Kyrenia Fault may not be the only fault system involved. The Anatolia-Cyprus pole we calculated would also predict pure strike-slip relative motion on a fault running along the coast (the coast nearly follows a small circle for this pole) in the prolongation of The Kozan fault. This fault has been proposed to move at rates of 4–8 mm/yr based on sediment depocenter migration (Aksu et al., 2014b). However, the GNSS residuals (although consistent in orientation with left-lateral motion along the coast, see Fig. 5) are barely above noise level at about 1 mm/yr. Most probably, the strike-slip motion between Cyprus and Anatolia is dominantly taken up along the Kyrenia fault and its NE prolongation (the KOF).

The Euler poles published by previous studies (Gomez et al., 2020; Reilinger et al., 2006) only predict moderate obliquity on the Cyprus subduction, at about 20° south of Cyprus, which corresponds to the critical obliquity threshold for the onset shear partitioning above a subduction zone (McCaffrey, 1992). The rotation motion of Anatolia vs. Sinai and the arc shape of the subduction cause a lateral variation of obliquity so that 20° (McCaffrey, 1992). It is thus possible that forces applied on the E and W boundaries of the Cyprus block play a role. The Cyprus block interacts with the Malatya block at its NE end and kinematic conditions on this boundary are extensional. It is thus possible that the forces applied in this zone near the triple junction influence the motion of the Cyprus block and particularly the amount of strike-slip taken up by the Kyrenia fault.

The fault coupling model of the Cyprus arc subduction is characterized by full locking from the surface to 20 km, which is consistent with (Feld et al., 2017; Welford et al., 2015), transitioning to partial locking between 20 and 30 km. The seismic activity during the instrumental period aligns with our coupling model, with the majority of earthquakes occurring within the 12–20 km depth range of the Cyprus subduction zone (see Fig. 6b). Conversely, the Kyrenia fault exhibits a distribution of locking extending from 0 to 7 km, transitioning to freely slipping behaviour beyond that depth. However, seismic activity on this fault is too low to provide an independent constraint on the depth range of the seismogenic zone.

In the vicinity of the triple junction the shear between Anatolia and Arabia is distributed over a zone that we defined as a block (the Malatya block). On the East Anatolian fault, we obtain 5.1–6.2 mm/yr strike-slip rates that are lower than most previous studies (Reilinger et al., 2006). However, analyses of a high-resolution velocity field obtained by combining GNSS and InSAR data (Weiss et al., 2020) found a laterally varying interseismic loading rate on the EAF (Güvercin et al., 2022; Li et al., 2023). The strike-slip rate on the Malatya fault we obtained (1.7 mm/yr) is consistent with previous studies (2013). The strike-slip rate on the Çardak fault we obtained from the block model is 1.8 mm/yr and is comparable with a 2 mm/yr slip rate calculated from geomorphological offsets (Westaway, 2003). In addition, the principal strain rate orientations between these two faults indicate left-lateral shear co-linear with shear along the East Anatolian fault zone (see Figs. 8 and 4) with an average rate of 35 nanostrain per year and a small extensional component (Fig. 3). Assuming simple shear, the internal deformation of the Malatya block over an average width of 50 km amounts to about 1.7 mm/yr. The strike-slip motion between Anatolia and Arabia thus appears to be distributed between the East Anatolian fault, the faults defining the northern boundary of the Malatya block and the internal deformation of the block. These three components add up to 8.5–9.7 mm/yr, which is consistent with previous estimations of Arabia/Anatolia plate motion (Bletery et al., 2020; Reilinger et al., 2006). The block boundaries we propose differ from the geometry proposed by (Özkan et al., 2023) in that they connect the Karata s-Osmaniye fault to the EAF near Kahramanmaraş and thus do not allow shear partitioning east of the triple junction. Their geometry cannot account for the loading of Çardak Fault which was ruptured during the second earthquake of the February 6, 2023, earthquake doublet (Barbot et al., 2023; Toda et al., 2023). Our solution with an N-S transtensional boundary connecting to Çardak fault has a lower misfit and also better represents mapped active faults (Emre et al., 2018; Fig. 1b). We conclude that part of the EAF motion is partitioned from the main fault in a broader zone around the triple junction, resulting in decreased interseismic loading rates on the main strand of the EAF.

Block-based model inversion indicates that the Levant fault accommodates a 3.5–4.7 mm/yr slip rate, and it decreases slightly but steadily from south to north. This result is consistent with previous studies such as (Al Tarazi et al., 2011; Gomez et al., 2020; Hamiel and Piatibratova, 2021; Li et al., 2024; Sadeh et al., 2012). Gomez et al. (2020) found a decreasing motion northward from 5.0 mm/yr to 2.2 ± 0.5 mm/yr as they consider two block boundaries transferring part of the Levant fault motion to hypothetical offshore structures. Our inversion results find less transfer of motion to offshore faults. For instance, only 0.7 ± 0.7 mm/yr are transferred to the western boundary of the Latakia block compared to 1.5 ± 0.5 mm/yr for the corresponding block boundary in (Gomez et al., 2020). However, uncertainties are large so that inversion results remain compatible. Several factors may contribute to a lower velocity on the offshore block boundary and a higher uncertainty in our results. Our data set is different and possibly with a higher noise level as it combines several different studies, a different inversion code is used that allows internal block deformation, and block boundaries also differ.

During the Feb 6, 2023 sequence, ruptures occurred on the southern and northern boundaries of Malatya block. The main shock (Mw 7.8) occurred on the East Anatolian and Levant faults while the large Mw 7.6 aftershock occurred on Çardak Fault (CF), corresponding to a moment magnitude about half of that of the main shock (Hussain et al., 2023; Toda et al., 2023). The occurrence of these events shows that both boundaries are seismically active and present a high seismic hazard. Our interseismic coupling inversion results find both faults are fully and homogeneously locked down to at least 10 km depth, but we cannot exclude that heterogeneities may be present but not resolvable with the GNSS data set we used. Based on InSAR and seismicity distributions, previous studies found shallow locking depth and aseismic creep on the EAF east of E038.5° longitude (Bletery et al., 2020; Cakir et al., 2023; Konca et al., 2021). Moreover, the depth of locking near the bend between EAF and Levant Fault (Fig. 9a) may be overestimated because of the simplified geometry of the block boundary where two faults branches are in fact present (the Nurdai-Pazarck Fault and the Pazarck segment of the EAF). Several research groups worked on the coseismic slip distributions on both EAF and CF using seismological records, SAR interferometry, optical (Sentinel-2) images, and coseismic GNSS data and various inversions and joint inversions have been published using these data (Barbot et al., 2023; Chen and Zhou, 2024; Jia et al., 2023; Li et al., 2023; Melgar et al., 2023; Toda and Stein, 2024). The average distribution of coseismic slip with depth is consistent with a locking depth of 10–15 km for both faults with a progressive transition down to 20 km (Jia et al., 2023). These coseismic slip models also find that the smaller earthquake has, in fact, larger displacements but over a smaller rupture length. This implies, taking into account interseismic loading rates, that the recurrence interval is very different on these two faults. Estimates of surface displacement from Sentinel-2 image correlation are 4 m on average over a large part of the main shock rupture with a local maximum of about 7 m and 6 m on average for the aftershock with a maximum of 8–9 m (Barbot et al., 2023; Hussain et al., 2023). Using our interseismic backslip estimates, recurrence intervals of 750 to 1500 years are inferred for the East Anatolian Fault and from 3000 to 5000 years on the Çardak Fault. It thus appears that triggering of the Çardak Fault does not occur each time a large earthquake occurs on the East Anatolian Fault, and probably does so less than once every three cycles on the EAF. Moreover, the Mw 7.8 earthquake is a multi-segment rupture (Barbot et al., 2023), and such events have longer recurrence intervals than those estimated for characteristic earthquakes on individual segments. Estimates from single segments range from 100 years to about 900 years along the East Anatolian Fault with the longest intervals and largest maximum magnitude (Mw 7.4) in the Kahramanmaraş triple junction area (Güvercin et al., 2022). It is important to note that that study took into account a westward decrease of interseismic loading rates along the EAF segments. The duration of seismic cycles involving a multi-segment rupture is even longer. Their estimation depends on complex scenarios that combine multi-segment and single-segment ruptures and ranges from 700 to 2500 years (Karabulut et al., 2023). We conclude that shear partitioning in the Kahramanmaraş triple junction is one of the factors contributing to very long earthquake cycles on the EAF. Moreover, the earthquake hazard on secondary faults must not be ignored even though large events on these slower faults may have even longer recurrence intervals.

Acquisition of new GNSS data on Cyprus and southern Turkey brings new insight into the deformation of the region around the Kahramanmaraş triple junction. It shows that the present-day deformation of Cyprus may be understood as a shear partitioning system between the Cyprus arc subduction and the Kyrenia fault, which appears to be a dominantly strike-slip boundary. Thus, the incipient collision with the Eratosthenes Seamount may not have yet perturbed much the kinematics of the Cyprus subduction. On the other hand, the northeast continuation of the shear partitioning system toward the Anatolia/Arabia collision zone brings further complexity to the Kahramanmaraş triple junction.

In the Anatolia-Arabia plate boundary, our study demonstrates that part of the motion on the East Anatolian Fault is distributed away from the main fault in the vicinity of the triple junction, as was previously shown for the Levant fault in the Arabia-Nubia plate boundary (Gomez et al., 2020). Thus, the EAF is not the only deformation source in the East Anatolian shear zone as deformation is also distributed on secondary faults such as Çardak and Malatya faults. The earthquake sequence that occurred on February 6 2023 emphasizes the earthquake hazard presented by secondary faults in complex plate boundary zones. In the study area, the kinematics of faults offshore of the Levant and Hatay coasts still need consideration. On land, the slip rates calculated on the faults bounding the Malatya block result in very long earthquake cycles of 750-2000 years for the part of the EAF that ruptured recently and probably 3000-5000 years for Çardak Fault. How such rare events may be considered in hazard assessment poses questions.

We thank our editors and reviewers whose constructive feedback and thoughtful suggestions have significantly improved this manuscript. We want to thank the T.C. Ministry of Defence General Directorate of Mapping for their contribution during the field studies and for sharing previous GNSS observations they measured. We are grateful to Ali Fahri Özten and Sebat Proje ve Mühendislik for the instrument supply and accommodation in Cyprus. We also want to thank Fatih Taskran, Mert Topal, Mustafa Ozan Güldoğan and Dr Ali Ĭhsan Kurt for their priceless effort during the field observations. Dr. Robert McCaffrey’s remarks, especially in establishing the block model, guided us. We also thank Xavier Le Pichon and Solène Antoine for their comments and discussions on understanding the recent earthquakes. This study was funded by the Istanbul Technical University Scientific Research Projects Coordination Unit with an MGA-2020-42584 ID Number research project. Moreover, it is part of the Ph.D. thesis of the corresponding author. The modeling part was carried out mostly at Centre Européen de Recherche et d’Enseignement des Géosciences de l’Environnement (CEREGE) in the scope of the TUBITAK 2214A International Research Scholarship during Ph.D. for Ph.D. candidates program with the project number 1059B142000638. The calculation of daily coordinates of GNSS sites reported in this paper was performed at TUBITAK ULAKBIM, High Performance and Grid Computing Center (TRUBA resources).