An up-to-date crustal deformation map of Iran using integrated campaign-mode and permanent GPS velocities

SUMMARY

We present the most extensive and up-to-date unified GPS velocity field for Iran. We processed the data collected during 10 yr (2006–2015) from the Iranian Permanent GNSS Network (IPGN) and combined them with previously published velocity solutions from the GPS survey measurements during 1997–2013. We analysed this velocity field using a continuum approach to compute a new strain-rate map for this region and we designed a block model based on the main geological, morphological and seismic structures. Comparison between both approaches suggests similar results and allows us to present the first comprehensive first-order fault-slip-rate estimates for the whole of Iran. Our results confirm most of the results from previous geodetic studies. However, we also show a trade-off between the coupling ratio of the Iranian Makran subduction interface and the kinematics of the faults north of the Makran in the Jazmurian depression. Indeed, although too scarce to accurately estimate a coupling ratio, we show that coupling higher than 0.4 on the plate interface down to a depth of 25 km will induce extension on the E–W faults in the Jazmurian region. However, the sites close to the shoreline suggest a low coupling ratio; hence, the coupling on this plate interface is probably more complicated than previously described and the Iranian Makran subduction interface mechanical behaviour might be similar to that on the Hellenic subduction zone.

1 INTRODUCTION

Iran encompasses a large part of the domain accommodating the Arabia–Eurasia convergence (e.g. Jackson & McKenzie 1984; Vernant et al.2004a). Most of the country is characterized by continental collision accommodated by several mountain ranges, including the Zagros, Alborz and Kopeh Dagh, but subduction takes place in the Makran. As a result of the present-day active tectonics, intensive seismicity occurs, and the long recorded history of the country provides a fairly well-documented record of historical seismicity (Berberian & Yeats 1999).

In order to study the present-day crustal deformation, several survey-GPS (sGPS) networks were established and measured in Iran from 1997 to 2013 (Tatar et al.2002; Nilforoushan et al.2003; Vernant et al.2004a,b; Masson et al.2005, 2006, 2007, 2014; Bayer et al.2006; Hessami et al.2006; Walpersdorf et al.2006, 2014; Tavakoli et al.2008; Peyret et al.2009; Djamour et al.2010, 2011; Mousavi et al.2013). Raeesi et al. (2017) combined all the aforementioned sGPS studies in a unified reference frame velocity field. However, following the Bam earthquake in 2003 (M_w 6.5, Talebian et al.2004), the National Cartographic Center of Iran (NCC) established the Iranian Permanent GNSS Network (IPGN) distributed throughout the whole country for continuous-GPS (cGPS) measurements. Parts of the IPGN measurements were published (e.g. Djamour et al.2010, 2011, Mousavi et al.2013; Walpersdorf et al.2014, Frohling & Szeliga 2016), but no velocity field for the whole network of 110 sites, with almost 10 yr of recording, has been reported. In this paper, we present the velocity field for the IPGN network based on the data collected from 2006 to the end of 2015. In order to provide the densest velocity field for Iran, we combine our IPGN velocity solution with the Raeesi et al. (2017) solution, and we extend its geographical coverage by also combining velocities from Frohling & Szeliga (2016) for the Pakistani Makran and Reilinger et al. (2006) for the Caucasus and eastern Turkey area.

Finally based on our unified velocity field we use both continuum and microplate approaches to analyse the velocity field and provide the most comprehensive strain-rate map of Iran and the slip-rates for the main Iranian faults.

2 IPGN

After the destructive 2003 Bam earthquake (Talebian et al.2004), which struck the Kerman province of southeastern Iran, the NCC established the IPGN for continuous monitoring of active crustal deformation in the whole country to allow more detailed studies of distributed and localized deformation (Djamour et al.2006). Presently, this network consists of about 130 cGPS stations, which were set up on geodetically designed concrete pillars, deeply rooted in stable ground. The network is measured with a mixture of Ashtech Z12 and Trimble NetR9 GPS receivers with choke-ring antennas. At the end of 2014, some new receivers such as CHC N71, COMNAV and SOUTH NETS8+ were added for other applications (e.g. DGPS and network RTK).

The network design and the optimal site locations were determined by the Geological Survey of Iran (GSI) and the International Institute of Earthquake Engineering and Seismology (IIEES), and also by French scientists from Montpellier, Grenoble-Alpes and Strasbourg Universities based on tectonics and population density. Areas with high population density and high earthquake potential were specifically targeted with a higher density of GPS sites, including Central Alborz and Tehran, northwest Iran, northeast Iran and Zagros.

3 GPS DATA PROCESSING AND STRAIN-RATE CALCULATION

3.1 GPS data and processing of IPGN

Parts of the IPGN data have already been used for previous regional studies (Djamour et al. 2010, 2011; Mousavi et al.2013; Walpersdorf et al.2014). However, this is the first time that a comprehensive velocity field of the IGPN network obtained from all the measurements between 2006 and the end of 2015 is presented. The horizontal velocities were determined for 123 Iranian cGPS stations but 8 sites were excluded due to their high subsidence rates, up to 30 cm yr⁻¹, related to underground water pumping and also affecting the horizontal components (e.g. Argus et al.2005). In addition, the velocities of 5 stations were removed because of suspicious transient behaviour. Finally, we present the velocities of 110 IPGN stations. All the velocities are based on longer time-series than previously published, and hence reduced uncertainties, and the velocities of 37 IPGN sites are presented for the first time (listed in Supporting Information Table S1).

We processed the GPS data using the GAMIT/GLOBK 10.6 (Herring et al.2015) software. To tie the IPGN network to the latest version of the International Terrestrial Reference Frame, ITRF2014 (Altamimi et al.2016), 20 nearby IGS stations were included in the daily solution processing. In the first step, GPS data were processed in daily loose network solutions following the strategy of Mousavi et al. (2013). In the second step the daily solutions were combined with global quasi-observations of about 100 IGS stations around the world using a Kalman filter strategy. Finally, in a third step, we estimated a consistent set of coordinates and velocities by minimizing the velocity and coordinate residuals for a subset of sites between our solution and the ITRF2014 solution (Altamimi et al.2016). To estimate the velocities, we removed outliers and corrected for offsets in the time-series. A unique noise model for each permanent station was calculated to account for correlated errors in the time-series analysis. The algorithm used to model the data noise spectrum assumes that each time-series can be adequately modelled using a first-order Gauss Markov (FOGM) process noise (Gelb 1974; Reilinger et al.2006; Djamour et al.2010, 2011).

The FOGM model is estimated from individual time-series of stations by averaging the residuals over increasingly longer intervals that range from a minimum of 7 d to a maximum of 1/10th of the total time-series span. This estimated FOGM model is then used to predict the site velocity uncertainty based on the time span of the time-series (Reilinger et al.2006; Djamour et al.2010, 2011). The Eurasia-fixed reference frame coordinate is based on the Eurasia plate as defined by Altamimi et al. (2013).

3.2 Combination of the GPS velocity fields

We combine our newly estimated velocity field obtained for the IPGN stations with other available velocity solutions. We use the combination done by Raeesi et al. (2017) of the GPS velocity solutions of Iran. We also include some GPS velocities reported by Reilinger et al. (2006) in the Caucasus and eastern Turkey regions to better characterize the eastern Turkey–NW Iran–Caucasus collision zone deformation style. And finally, to add velocities in the eastern Makran region, we use the Frohling & Szeliga (2016) velocity solution.

Reilinger et al. (2006) studied the deformation of the Africa–Arabia–Eurasia convergence zone using GPS measurements of sGPS and cGPS sites from 1988 to 2005. Frohling & Szeliga (2016) processed 20 GPS stations covering the whole Makran subduction zone. Raeesi et al. (2017) combined the velocity field (1999–2011) reported by Zarifi et al. (2014) with published velocities (1997–2008) for eastern Iran (Walpersdorf et al.2014) in order to have an updated unified velocity field covering Iran.

Because these individual velocity solutions were obtained by different processing strategies, or by choosing different reference frames, the GPS velocity fields could not be directly merged. Therefore, we combined all the previously published velocity solutions by minimizing the residuals between the velocities of the common sites. We used the VELROT program version 1.01 (Herring et al.2015) to combine various velocity fields into a single homogeneous Eurasia-fixed reference frame based on ITRF2014 (Altamimi et al.2016). We used our IPGN solution as the reference velocity field and aligned all the previously published velocity fields (i.e. Reilinger et al.2006; Raeesi et al.2017) on it. The computed residual horizontal velocities have RMS values of 0.89 mm yr⁻¹ for the 64 common sites with the Raeesi et al. (2017) solution and 0.33 mm yr⁻¹ for the 9 common sites with the Reilinger et al. (2006) solution (Table 1). Finally, we have a combined velocity field including 110 permanent and 184 campaign GPS stations, plus 93 stations from Reilinger et al. (2006) and 7 sites from Frohling & Szeliga (2016) (Fig. 1 and Supporting Information Table S1).

Our denser velocity field (Fig. 1) confirms previous results on the deformation pattern and active tectonics of the region (e.g. Vernant et al.2004a; Masson et al.2014; Raeesi et al.2017). The Arabia–Eurasia convergence rate increases towards the east (from ∼19 mm yr⁻¹ at longitude 48° to ∼27 mm yr⁻¹ at longitude 58°, Fig. 1) as predicted from the Arabia–Eurasia Euler pole located in North Africa at 19.4 ± 1.4°E and 27.9 ± 0.5°N (Vernant et al.2004a, corroborating Euler pole locations of Sella et al.2002; McClusky et al.2000, 2003). As described by Vernant et al. (2004a), the Arabia–Eurasia convergence accommodation changes drastically at 58˚E of longitude. East of 58°E, most of the shortening is accommodated by the Makran subduction and the remaining part is absorbed by the Kopeh Dagh. The differential between the velocities of the sites located on the Makran coast (jasc, chbr) with the sites on the Oman gulf (musc) suggests that subduction rates increase from west to east (Supportign Information Table S1). West of 58°E, the shortening is distributed between the Zagros, Alborz and Caucasus mountain ranges. This difference induces shear motion between Central Iran and the Helmand block, which is accommodated in eastern Iran along the Lut block boundaries (Vernant et al.2004a; Walpersdorf et al.2014).

The velocities of sites located east of the Iranian border are very small and suggest that the Helmand block is part of the Eurasian plate (Jackson & McKenzie 1984; Vernant et al.2004a). West of Iran, along the Arabian promontory, the velocities are large: the Anatolia plate rotates and moves westwards towards the Hellenic Trench (Sengor et al.1985; Reilinger et al.2006) while the northern part of the Iranian plateau and Lesser Caucasus moves towards the northeast.

3.3 Continuum description—strain rates

Earlier studies presented strain rates derived from GPS velocity fields for various parts of Iran and in some cases the whole country (e.g. Walpersdorf et al.2006; Masson et al.2005, 2007, 2014; Zarifi et al. 2014; Raeesi et al.2017) using different approaches and velocity fields. In this study, we choose to use an approach based on the inversion of baselines first proposed by Spakman & Nyst (2002) and further developed by Masson et al. (2014) and named STIB (Strain Tensor from Inversion of Baselines). Details can be found in Masson et al. (2014), but to sum it up, the method is similar to that used for seismic tomography, the seismic rays being replaced by the GPS baselines. Therefore, this method avoids the strong dependency of the strain-rate field on the quality of the velocity field by integrating all the information included in the GPS data linked by their baselines crossing the nodes of the grid where the strain-rate tensor are computed. Masson et al. (2014) demonstrated the STIB method efficiency by applying it to the GPS velocity field published by Nocquet (2012).

In this study we similarly apply the STIB method to our denser velocity field. For studying the strain-rate pattern at the scale of Iran, a grid spacing of 1° in latitude and longitude was used. Fig. 2 shows the strain-rate field deduced from the velocity field shown on Fig. 1. The first-order features observed on this large-scale map are: (1) the significant shortening between Arabia and southern Iran across the Makran Trench; (2) the transition zone between the Zagros collision and Makran subduction (54°–58°E) with a right-lateral strike-slip motion related to Minab–Zendan–Palami (MZP) fault system; (3) the decreasing compressional strain rate from eastern Zagros to western Zagros; (4) the right-lateral strike-slip pattern in NW Iran–eastern Turkey along the North Tabriz Fault (NTF), the Gailatu-Siah Cheshmeh-Khoy (GSK) fault, the Chalderan fault, and continuing along the North Anatolian fault (NAF) in Turkey; (5) the large shortening accommodated across the northern Kura Basin and the southern Greater Caucasus in Azerbaijan; (6) the pure N–S compressive regime of the Apsheron Balkan sill (ABS). Most of these features have been described in previous studies (e.g. Masson et al.2005, 2007, 2014; Reilinger et al.2006; Kadirov et al.2012; Raeesi et al.2017). However, the almost pure compressive regime focused on the ABS had been already proposed by Djamour et al. (2010) using a kinematic model, but was debated with proposed motion of the South Caspian block implying a more transpressive motion on the ABS (e.g. Jackson et al.2002; Hollingsworth et al.2006). The ABS is the offshore eastward continuation of the Caucasus within the Caspian Sea, separating the South Caspian Basin (SCB) from the North Caspian Basin (NCB). The SCB is composed of oceanic crust while the NCB is of continental crust (Mangino & Priestley. 1998). Even without GPS data close to the ABS, the large strain is well located on the ABS and does not spread to the whole SCB. This unexpected sharp location of strain shows the specificity of the STIB method.

3.4 Block model description

The best way to describe the surface deformation is still debated. Thatcher (1995) considers the issue in terms of two end members: the micro plate and the continuum descriptions. Currently no consensus has been reached. Accordingly, we also investigate a block model derived from previously proposed models (Vernant et al.2004a; Bayer et al.2006; Djamour et al.2010, 2011; Walpersdorf et al.2014) and the main active tectonic structures of Iran. The models include block rotations on a sphere with elastic strain accumulation on block-bounding faults following the formulation of Okada (1985). All faults must be located on block boundaries and no ‘unconnected’ fault segments can exist. No permanent deformation exists within the blocks. Relative block motions (relative Euler vectors) are estimated using the DEFNODE code (McCaffrey 2002) by minimizing the GPS residual motions within the blocks in a least squares sense. We report the fault normal slip rate (normal and thrust, positive compression) and fault parallel slip rate (strike-slip, positive right-lateral) component for each fault segment. Fig. 3 shows the block geometry and the fault slip rates obtained. Fault slip rates are not only GPS velocity dependant but also block geometry dependant (see discussion in Karakhanyan et al.2013), therefore we do not provide fault slip rate uncertainties, but looking at previously published block models (Djamour et al.2010, 2011; Walpersdorf et al.2014), it is clear that these uncertainties can reach 2–3 mm yr⁻¹. Indeed, for the Alborz region, fault slip rates can vary by 2–3 mm yr⁻¹ with respect to Djamour et al. (2010) even though the geometries are quite close. We choose to use a simpler block geometry for the Sanandaj-Sirjan and Central Iran area than Walpersdorf et al. (2014) which appears to be supported by the very low residuals with this geometry. However, Walpersdorf et al. (2014) had several more blocks and fault slip rates of 2–5 mm yr⁻¹ where it does not seem our model requires block boundaries. Plotting the block geometry on the maximum shear strain-rate magnitude map (Fig. 4a) shows consistent results between the STIB method and the block model since the highest amplitude maximum shear strain rates are consistent with the highest fault slip rates on the block boundaries. Moreover, following Allmendinger et al. (2007), the rotation rates derived from the continuum description are consistent with strike-slip motions along block boundaries (Fig. 4b). The right-lateral displacement along the Lut block, the western SCB boundary, the Kazerun fault, the transition between the Zagros and the Makran or the Chalderan-Tabriz faults are clearly expressed as clockwise rotations. On the other hand, the left-lateral displacement along the eastern Alborz range and SW of the Zagros corresponds to counter clockwise rotations. These rotations are a result of strain accumulation indicating a high likelihood of earthquakes along these strike slip fault segments.

4 DISCUSSION

4.1 NW Iran–Caucasus–eastern Turkey

In this region, the deformation is partitioned between almost pure shortening across the southern Greater Caucasus, presumably on the Main Caucasus thrust fault (MCTF), and right-lateral strike-slip along the NT, GSK and Chalderan fault system. This partitioning is due to the motion of the Lesser Caucasus area towards the NNE in a Eurasia-fixed reference frame (Fig. 1). It induces E-W compression between the Talesh and the SCB, since the Caspian basin is moving NNW. We do not have the resolution to define if the shortening along the southern limit of the Greater Caucasus is accommodated only on the MCTF. Based on Sokhadze et al. (2018) it seems that even though most of the deformation probably occurs on the MCTF, some might be accommodated south on the Lesser Caucasus thrust fault (LCTF) and in the Kura Fold and Thrust Belt (Forte et al.2010, 2013; Kadirov et al.2012). Based on our block model (Fig. 3) and maximum shear strain-rate map (Fig. 4a), the shortening on the Lesser Caucasus thrust fault has to be low (∼2 mm yr⁻¹) in agreement with Sokhadze et al. (2018) results. The total shortening accommodated by the MCTF, LCTF and Kura fold and thrust belt decrease from ∼11 mm yr⁻¹ in the east to ∼6 mm yr⁻¹ in the west, as reported by Reilinger et al. (2006).

The Lesser Caucasus has been divided in several blocks by Karakhanyan et al. (2013), but the overall deformation was low, mainly below 1 mm yr⁻¹. Our results are consistent with Karakhanyan et al. (2013) findings and show that as suggested by Reilinger et al. (2006) and Djamour et al. (2011), the Lesser Caucasus is a low deformation area bounded by relatively high slip rate faults. South of the Lesser Caucasus, the NT, GSK and Chalderan fault system has an average slip rate of 8 mm yr⁻¹. The fairly well located strain along the Tabriz fault broadens towards the west in the Chalderan fault area (Fig. 4a). Broadening of the strain-rate field is consistent with crustal shortening and associate thrust faulting in the Van Lake area, both indicated in Fig. 2, in agreement with the focal mechanism of the 2011 M_w 7.2 Van earthquake (Dogan & Karakas, 2013 ). If the long-term seismicity level of the Chalderan fault is low as suggested by Hisarli et al. (2016), then creep must be occurring in this area. East of Tabriz, the deformation seems to be distributed in the mountainous area north of the Tabriz fault (Fig. 4a), from the Mount Bozgush to the Ahar fault, consistent with Djamour et al. (2011) observations and the 2012 twin Ahar earthquakes (Copley et al.2013).

East of the Talesh range, EW shortening is observed indicating compression between the SCB and the Talesh consistent with the seismicity of the NS faults extending from the Alborz to the Kura basin along the Caspian shoreline (Jackson et al.2002). This pattern was previously proposed by Masson et al. (2014) but without the large spatial consistency and amplitude observed here. The clockwise rotation indicates right lateral motion along theses faults in agreement with the block model results and geological observations (i.e. Philip et al.1989). This motion of the Lesser Caucasus area appears to be inconsistent with the Arabia–Eurasia collision and as yet to be explained but is most likely related to deep processes (Vernant & Chéry 2006a) and the complex lithosphere structure of the region (e.g. Mumladze et al.2015).

4.2 Alborz and SCB

Given the high GPS site density we use a denser grid spacing of 0.4° for the Alborz region. The strain-rate map shows that the deformation within the Alborz range is changing in style and intensity across and along the belt (Fig. 5). The eastern part accommodates mainly left-lateral motion across the range and shortening on the southern margin of the belt. This is consistent with the counter clockwise rotation observed for the eastern Alborz (Fig. 4b) and the large left-lateral slip rates obtained from the block model, while fault normal rates are marginal (Fig. 3). In the central Alborz (from 51°E to 52.5°E) where the high counter clockwise rotation ends (Fig. 4b) and the highest maximum shear strain rate occurs, compression is mainly perpendicular to the strike of the belt (Fig. 5). In the block model, it is expressed as a high shortening rate on the Khazar thrust north of the Alborz (Fig. 3). On the faults south of the central Alborz, the block model suggests very low slip rates (∼1 mm yr⁻¹), given the uncertainties on the slip rates, it is consistent with geomorphologically derived fault slip rates for North Tehran fault (Ritz et al.2006) or the Pardisan fault (Talebian et al. 2016). In spite of denser GPS coverage, it is impossible to decipher more precisely the slip rates of all the faults bordering the Alborz range. However, the deep locking depth found by Djamour et al. (2010) for the Khazar thrust, in agreement with the depth of the M_w6.2 Baladeh earthquake in 2004 (Tatar et al.2007), suggests that some deformation may occur offshore, north of the range.

4.3 Zagros fold-and-thrust belt

The Zagros mountain belt extends from eastern Turkey to the Strait of Hormuz. The belt is divided into two parts by the Dena-Kazerun-Borzajan (DKB) NS fault system at ∼51°E. The structures northwest of the DKB (western Zagros) are oblique to the Arabia–Eurasia shortening direction. Those southeast of the DKB (eastern Zagros) are more perpendicular to the direction of Arabia–Eurasia shortening (Talebian & Jackson 2004). Walpersdorf et al. (2006) presented the first dense GPS velocity field covering the total Zagros mountain belt suggesting that the deformation in the western Zagros is partitioned between shortening perpendicular to the axis of the mountain belt and right-lateral strike-slip motion on several northwest–southeast trending faults, in particular the Main Recent fault (MRF). This interpretation was previously proposed by Talebian & Jackson (2004) but is still debated since numerical modelling of oblique convergence suggests that partitioning cannot be complete in the Zagros (Vernant & Chéry 2006b). On the other hand, in the south eastern Zagros, the deformation is closer to pure shortening perpendicular to the fold and thrust belt and restricted to the Persian Gulf shore (Walpersdorf et al.2006; Hessami et al.2006). Our block model quantifies the differences in the accommodation of the deformation east and west of the DKB fault system (Fig. 3). The maximum shear strain rate and the rotation maps suggest also that the highest deformations are in the southeastern Zagros and point a counter clockwise rotation along the eastern part of the southeastern Zagros front in agreement with the left-lateral component suggested by the block model. From the Zagros–Makran transition zone towards the southeast Zagros, we observe a rotation of the compressional axes of the strain-rate tensors (Fig. 2), roughly following the direction perpendicular to the shoreline. The direction of the shortening is in good agreement with the focal mechanisms. According to analogue models, this rotation of the main compressional axes of the strain-rate tensors would be partly due to the existence of the Hormuz salt formation acting as a viscous decollement beneath the sedimentary cover (Nilforoushan & Koyi 2007).

In the western Zagros, the block model suggests a partitioning between the shortening along the southern front of the Zagros and strike slip along the northern boundary. However, the extensional component along the MRF between the Zagros and the Sanandaj-Sirjan zone is unexpected. It might be related to a kinematic issue of the block model, the block geometry forcing this extension, or it might be a geological deformation. Talebian & Jackson (2002) reported evidences of normal fault scarp along the MRF, but they interpreted them as strike-slip induced extension. This remains to be studied in detail, but the strain-rate tensors also suggest extension along the MRF (Fig. 2). The DKB fault transfers the fault slip and active deformation along the MRF to the southern front of the eastern Zagros preventing the northern part of the eastern Zagros from having a high seismic activity as suggested by previous authors (Talebian & Jackson 2004; Tavakoli et al.2008).

North of the Zagros, the Sanandaj-Sirjan zone shows very low strain rate (Figs 2–4) and can be considered to the first order as a block as previously proposed based on geodetic observations (Vernant et al.2004b; Masson et al.2005, 2007), and in agreement with its low seismic activity. If the northern part of the Dehshir fault seems to be active with mainly a right lateral component of ∼2 mm yr⁻¹ in agreement with geomorphological evidences (Le Dortz et al.2009), it is hard to see an evidence of this fault actively offsetting the Sanandaj-Sirjan zone in its southern part as suggested by Le Dortz et al. (2011).

4.4 Zagros–Makran transition zone and Makran subduction

The major role of the MZP fault system consists of accommodating the transition from subduction in the Makran and continental collision in the Zagros. The high values of the maximum shear strain-rate magnitudes (Fig. 4a) and the high clockwise rotation rates (Fig. 4b) highlight its high right-lateral slip-rate. This is confirmed by the high overall right-lateral component obtained by the block model (∼17 mm yr⁻¹), this value is consistent with former geodetic estimates (Bayer et al.2006; Peyret et al.2009) and its geomorphologically determined slip rate (Regard et al.2005). To the east of this fault zone, subduction along the Makran accommodates most of the shortening with a rate of ∼20 mm yr⁻¹ (Fig. 3) confirming earlier results (Vernant et al.2004a). Most of the shortening takes place offshore and due to its scarcity, our network is not well designed to study the coupling of the Makran subduction zone interface. However, following Vernant et al. (2014), we ran forward models where the subduction interface dips at 10˚ to the north according to the seismicity and the location of the volcanic arc. We used different coupling ratios ranging from 0 (i.e. no coupling) to 1 (i.e. full coupling) down to a depth of 25 km. The results are reported on Fig. 6(a). Because most of the coupled part of the modelled subduction interface is offshore, only the sites close to the shoreline are highly affected by the coupling. Based on these limited observations, we suggest that the Pakistanis part of the Makran may be more highly coupled than the Iranian part segment. This is consistent with the seismicity of the Makran (Byrne et al.1992). In addition, if the Iranian Makran subduction interface is coupled, the most coupled part must be offshore with a swallower locking depth than 25 km. However, we cannot rule out a low coupling ratio below a depth of 25 km, in this case large earthquakes could still strike the Iranian Makran but with longer repeat times as has been proposed for the Hellenic subduction interface (Shaw & Jackson 2010; Vernant et al. 2014). Our result differs from previous studies suggesting a rather highly coupled subduction interface for the whole Makran (Frohling & Szeglia 2016; Penney et al.2017) and show that further studies with denser GPS networks are needed to constrain better the coupling ratio of the Makran subduction zone.

Based on the strain-rate tensor and the block model (Figs 2 and 3), the deformation between the Makran prism and the Lut block seems to be low but with an NS shortening. Depending on the coupling ratio of the Iranian Makran subduction interface, the EW faults between the Makran and Jazmurian can be either in extension or in compression (Fig. 6b). If they are in extension as suggested by Penney et al. (2017), then the coupling ratio should be higher than 0.4, but in this case larger residuals will be observed for the sites along the shore. Without a denser network, it will be hard to answer these questions more precisely.

4.5 East and northeast of Iran

The large difference between the shortening accommodated by the Zagros and the Makran induces a significant NS shear between Central Iran and Afghanistan (e.g. Berberian & Yeats 1999; Berberian et al.1999; Wellman 1966; Jackson & McKenzie 1984; Vernant et al.2004a). Most of this NS shear is accommodated in eastern Iran by NS oriented faults; from east to west, the Neh, Nayband, Kuhbanan, Anar and Dehshir right-lateral strike-slip faults. The strain-rate tensors of Fig. 2 indicate that most of this NS shear between central Iran and the Helmand block is localized on the fault systems bounding the Lut block to the east (the Neh fault systems) and to the west (Kuhbanan fault and Nayband–Gawk–Bam fault systems). The maximum amplitude shear strain rate and the rotation rate (Fig. 4) show that the most eastern fault is the most active, consistent with the block model that indicates decreasing right lateral slip rates from east to west (Fig. 3). The high maximum amplitude shear strain rates north of the Neh fault system (NE Lut) and south of the Gowk–Kuhbanan–Nayband fault system in the Bam area (SW Lut) are correlated with high seismic activity. For the Bam area we also have large residuals from the block model (Fig. 4). These observations for the Bam region are most likely related to post-seismic activity of the 2003 Bam earthquake (Talebian et al.2004). Consequently, the high maximum shear strain amplitude north of the Neh fault system might be also related to post-seismic processes. The obtained fault slip rates are in agreement with previously published work (Vernant et al.2004a; Walpersdorf et al.2014). The Neh fault system has a right-lateral strike-slip rate of ∼7 mm yr⁻¹, the Nayband fault system a rate of ∼5 mm yr⁻¹ and the northern Dehshir fault a right-lateral motion of ∼2 mm yr⁻¹.

Further north, we observe a very low strain on the Doruneh fault (Fig. 2), but the derived rotation suggests a clockwise rotation (Fig. 4b) that is in agreement with the right-lateral strike slip suggested by the block model north of the Lut (Fig. 3). Five mm yr⁻¹ of northward motion with respect to stable Eurasia (Fig. 1), are transferred to the north of the Doruneh fault system and gradually decrease to zero at the Turkmenistan border. The shortening pattern of ∼4 mm yr⁻¹ in the NW-striking Binalud Mountain ranges (Figs 2 and 3) changes to a left-lateral strike-slip pattern in the Eastern Alborz. This pattern is related to the north–westward motion of the SCB (Djamour et al.2010; Mousavi et al.2013, 2015), inducing in the same time a right-lateral motion along the Ashkabad fault bordering the Kopeh Dagh (Walters et al.2013).

5 CONCLUSION

We presented the first velocity field from the IPGN (110 stations). Adding the velocities obtained from the GPS survey measurements in Iran as well as from Turkey and the Caucasus, we produced a dense velocity field (402 stations) for Iran and adjacent areas. Using a continuum approach by computing strain-rate tensors we compare it to a block model that we designed following the main geological, morphological and seismic structures. Comparison between both approaches suggests similar results and allows us to present the first comprehensive, first-order fault-slip-rate estimates for Iran. Our results confirm a partitioning of deformation north of the Arabian promontory with right-lateral strike-slip along the Chalderan-Tabriz fault system south of the Lesser Caucasus and shortening north of the Lesser Caucasus along the Greater Caucasus southern front. We also confirm a first-order partitioning in the western Zagros with shortening along the southern front and right-lateral strike-slip along the Main Recent fault. This fault is the border between the western Zagros and the quasi-rigid Sanandaj-Sirjan area to the north. An unexpected low extension rate is reported normal to the MRF by both approaches, this is surprising, but more observations are definitely needed before drawing any conclusion. The strike-slip motion of the MRF is transferred to the southern front of the eastern Zagros by the Kazerun fault. In the eastern Zagros, most of the deformation is concentrated on the southern front and the transition between the Zagros collision and the Makran subduction; this transition is accommodated by the Zendan–Minab–Palami fault system with an overall right-lateral strike-slip rate of ∼17 mm yr⁻¹. The Makran subduction rate estimate depends in part on the coupling of the subduction interface, unfortunately the network is too scarce to accurately estimate the coupling directly. However, a trade-off exists between the coupling ratio and the kinematics of the E–W faults between Makran and the Jazmurian depression. With a low coupling (i.e. <0.4) on the Iranian Makran subduction interface, these faults should be in compression, while with a high coupling ratio, they should be in extension. Penney et al. (2017) suggest normal faults in this area, but the GPS measurements along shore appear to require for a low coupling below 25 km depth to the contrary of the measurements in Pakistan. In light of these uncertainties, further detailed studies are needed to better assess the seismic behaviour of the Makran subduction. Finally, the continuum approach shows a localized almost pure shortening on the ABS, which is consistent with the estimated rotation of the SCB from the GPS sites surrounding it.

SUPPORTING INFORMATION

Table S1. Combined GPS velocities, their locations, velocities (mm yr⁻¹) and their 1σ uncertainties with respect to the Eurasia fixed frame (EURA08). Cor is the correlation coefficient between E and N velocities. IPGN: Iranian Permanent GNSS Network; Raeesi: Raeesi et al. (2017); Reilinger: Reilinger et al. (2006); Frohling: Frohling & Szeliga (2016). Bold rows indicate the Iranian permanent GPS sites presented for the first time. The rows marked with star represent the IGS stations which have been processed with the Iranian Permanent GNSS Network.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

ACKNOWLEDGEMENTS

We would like to thank the National Cartographic Center (NCC) in Tehran for the Iranian Permanent GNSS Network (IPGN), all Geodynamic and DGPS experts, and in particular Y. Djamour, F. Tavakoli and M. Sedighi. We thank Morteza Talebian and an anonymous reviewer for their comments that have greatly improved this work. We acknowledge the support of the University of Strasbourg for EOST's visit awarded to FK and FN. This research was supported in part by the ‘Center National d'Etudes Spatiales’ (CNES-TOSCA). We also thank Yann Ziegler for providing STIB program for strain calculations. The maps in this paper were produced using the Generic Mapping Tools (GMT) software.

REFERENCES