# Strain rate tensor in Iran from a new GPS velocity field

Frédéric Masson, Maximilien Lehujeur, Yann Ziegler, Cécile Doubre

► To cite this version:

Frédéric Masson, Maximilien Lehujeur, Yann Ziegler, Cécile Doubre. Strain rate tensor in Iran from a new GPS velocity field. Geophysical Journal International, 2014, 197, pp.10-21. 10.1093/gji/ggt509. hal-00967710

HAL Id: hal-00967710

https://hal.science/hal-00967710

Submitted on 25 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Strain rate tensor in Iran from a new GPS velocity field

Frédéric Masson, Maximilien Lehujeur, Yann Ziegler and Cécile Doubre

IPGS/EOST, Université de Strasbourg/CNRS, 5 rue Descartes, F-67084 Strasbourg Cedex, France. E-mail: frederic.masson@unistra.fr

Accepted 2013 December 16. Received 2013 December 3; in original form 2012 November 19

## SUMMARY

The aim of this paper is to determine the strain rate tensor (SRT) for the Iranian region. In this study, (1) we apply a method of computation of the SRT never used for the Iranian area and (2) we use a new GPS velocity field obtained from several previously published velocity fields. First, the method is described and tested on a synthetic case, which mimics the real Iranian case. The synthetic tests confirm that the method allows us to both retrieve high gradients of the strain rate field and reduce the effect of an erroneous velocity vector. Second, the method is applied to a real data set covering the Arabia-Eurasia collision zone in Iran. We particularly focus on the Zagros–Makran transition zone, the Central Iran region and the northernmost part of the Arabia-Eurasia collision zone (NW Iran–Caucasus–East Turkey). Whereas the main characteristics of the obtained SRT are consistent with known tectonic features, important new results are found in the Central Iran, with the strike-slip style along the Anar and Deshir faults, and the Zagros–Makran transition zone, with a north–south variation of the SRT along the Zendan–Minab–Palami fault system. We link these results to recent active tectonic studies.

**Key words:** Tomography; Satellite geodesy; Seismicity and tectonics; Dynamics: seismotectonics.

## 1 INTRODUCTION

The active tectonics in Iran results from the north–south convergence between the Arabian and Eurasian plates (Jackson & McKenzie 1984) at a rate of 22 mm yr⁻¹ (Vernant et al. 2004; Masson et al. 2007; Fig. 1). It involves a continental collision (Falcon 1974; Berberian & King 1981) except along the Makran, where a remnant part of the Tethys oceanic lithosphere subducts northward beneath southeastern Iran (Byrne et al. 1992). Several GPS networks have been installed and measured during the last 15 yr (Vernant et al. 2004; Bayer et al. 2006; Walpersdorf et al. 2006; Masson et al. 2007; Djamour et al. 2010). Using these recent data, we compute an up-to-date strain rate tensor (SRT) of the region that we compare with the seismicity.

Computing the SRT from a geodetic velocity field is a standard method used for studying regional deformation (e.g. Tarantola et al. 1979). The representation of the SRT is very convenient because (1) it gives a broad picture of the geodetic results, leading to an easy comparison with other data sets such as seismicity and (2) the SRT is independent of any reference frames, unlike the computation of the velocity field. Three main approaches are followed to calculate the SRT: the triangle approach (e.g. Frank 1967), the gridded approach using either a bi-cubic spline interpolation (e.g. Haines et al. 1998) or a spherical wavelet-based multiscale approach (Tape et al. 2009) and the inversion of baseline variations (Spakman & Nyst 2002). In the triangle approach, the area under study is divided into triangles (in general using Delaunay triangulation approach), whose summits correspond to sites of velocity measurements (GPS

sites) and the strain parameters are computed within each triangle from the three nearest velocity vectors. In the gridded approach, a functional relation between position and displacement with specific basis functions is defined to obtain an interpolated velocity grid and the SRT is computed from the interpolated velocity grid. The quality of the results obtained from both methods is therefore strongly dependent on the quality of the velocity field. This is especially true in the case of the triangle method, as shown by Masson et al. (2005): an error on the velocity at only one station induces a large error on the SRT calculated for the triangles adjacent to this station. Wu et al. (2011) give a review of the first two types of methods.

In order to avoid the strong dependency of the velocity field on the quality of the strain field, Spakman & Nyst (2002) have developed a method (hereafter called SNM for Spakman and Nyst method) of computation of the strain that uses all the information included in the velocity field and not merely the variation of velocity between two neighbour stations. To do so, they do not use the velocity field but the variations of distance between all pairs of stations. In addition, their method offers the possibility to include discontinuities in order to estimate the slip on active faults.

While this last method appears very powerful, it is too rarely implemented. Bos et al. (2003) used the SNM to describe the present-day surface deformation of Taiwan by computing the velocity gradient field and fault slip from GPS velocity vectors. Bos et al. (2004) applied the SNM to determine the coseismic surface deformation field associated with the 1999 Izmit earthquake, by using GPS velocity vectors and InSAR data as well. Pollitz & Nyst (2005)
<figure><img src="image_2.png" /><figcaption>**Figure 1.** GPS Iranian network and velocity field. In this study, we use the last homogeneous publicated velocity field of Iran (Nocquet, 2012). Deh: Dehshir fault, MRF: Main Recent Fault, MZT: Main Zagros Thrust, ZMP: Zendan Minab Palami fault zone, Nay: Nayband fault, NTF: North Tabriz Fault, KFZ: Kazerun Fault Zone, KF: Khazar fault, Hor. Str.: Hormuz Strait.</figcaption></figure>

and Bos & Spakman (2005) employed the SNM to study the San Francisco Bay Region and the southern California, respectively. Finally, Howe & Bird (2010) have developed a method inspired by the SNM, adding principal stress directions in the inverted data, in order to compute interpolated velocity fields and fault velocities. They do not focus on the strain tensor.

This paper is dedicated to a reassessment of the SRT in Iran using the SNM in its simple configuration: we only use the SNM to evaluate the SRT, without fault discontinuities. Unlike the classical SNM, this simplified version is particularly well suited to many parts of the Earth, such as the regions affected by the post-glacial rebound where the deformation is continuous (no fault) or in very low tectonic activity regions (most of the continental areas) where the velocity of the faults is far below the resolution of geodetic data. It is also adapted when the active faults are not well known in terms of location and interseismic slip rate, such as Central Iran. In this paper, first, we describe the simplified SNM in details, second we

apply the method to synthetic data and third we apply the method to the GPS data collected in Iran.

## 2 METHOD

Data used for the strain computation are the temporal variations of the baseline components. It is assumed that the deformation is constant over time. Thus, the slopes of the baseline components are directly related to the deformation that occurs between each pair of stations. As shown in Fig. 2, the main idea is to distribute along the baseline the deformation observed between the two stations. Using the full set of baselines, it is possible to get a map of the deformation on the surface covered by the whole network. This method is directly derived from delay-time seismic tomography (e.g. Aki et al. 1977) with the baseline variations in the role of the time delay (data) and the strain in the role of the velocity (model).
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Sketch of the STIB method. The space is divided in rectangular cells. The strain rate tensor in each cell is obtained from the variation of the length of all the baselines crossing the cell. The computation of the strain rate tensor is obtained by inversion of the variation of the length of all the baselines obtained from the geodetic network.</figcaption></figure>

The following development is valid for any 2-D geometry. Here, we are in the spherical case with the baselines being arcs of great circles and the velocities being considered as angular velocities. For the discretization (see eq. 2.13), we use a grid whose squares are limited by parallels and meridians.

We define the length $L_{AB}$ of the baseline between two stations A and B, that is, the length of the geodesic that links the two stations. This baseline can be characterized by its two components along two coordinate lines in the north-south and east-west directions, $L_N$ and $L_E$, respectively. The time derivatives of $L_N$ and $L_E$ equal the difference of the corresponding components of the velocity $v_N$ and $v_E$ between each stations, namely

$$ \left\{ \begin{aligned} \frac{dL_E}{dt} &= v_E^B - v_E^A, \\ \frac{dL_N}{dt} &= v_N^B - v_N^A, \end{aligned} \right. \qquad (2.1) $$

where the exponents A and B refer to the corresponding station.

The variation rates $\Delta L_E$ and $\Delta L_N$ for each component are defined as the derivatives that appeared in eq. (2.1) multiplied by a sign function that depends on the chosen convention on the sign of the deformation. Here, we chose a positive variation when the stations move away one from each other in longitude or latitude. These signed variation rates are the formal data of the problem :

$$ \left\{ \begin{aligned} \Delta L_E &= \operatorname{sgn}(L_E) \frac{dL_E}{dt}, \\ \Delta L_N &= \operatorname{sgn}(L_N) \frac{dL_N}{dt}, \end{aligned} \right. \qquad (2.2) $$

with $\operatorname{sgn}(x)$ the function giving the sign of $x$ (1, 0 or -1 for $x$ being positive, null or negative, respectively).

In a first time, we linearize the problem by assuming that the stations A and B are infinitely close together. If the station A is in $r_A$, the velocity at station B, located in $r_B = r_A + dr$, is then given at first order as a function of the velocity gradient in A, $\nabla v|_A$:

$$ \mathbf{v}(\mathbf{r}_A + \mathbf{dr}) \approx \mathbf{v}(\mathbf{r}_A) + \nabla \mathbf{v}|_A \cdot \mathbf{dr}. \qquad (2.3) $$

Using the so-called midpoint method reduces the error made on eq. (2.3) of order 1. This method consists in using the gradient not in A but in the middle of the $[AB]$ segment, $\nabla v|_{[AB]/2}$. A better approximation is thus

$$ \mathbf{v}(\mathbf{r}_B) \approx \mathbf{v}(\mathbf{r}_A) + \nabla \mathbf{v}|_{[AB]/2} \cdot \mathbf{dr}. \qquad (2.4) $$

Combining the eqs (2.1) and (2.4), we deduce

$$ \frac{d\mathbf{L}_{AB}}{dt} \approx \nabla \mathbf{v}|_{[AB]/2} \cdot \mathbf{dr}. \qquad (2.5) $$

By using eq. (2.5) along with eq. (2.2) and defining $\mathbf{dr} = dr_1\mathbf{e}_1 + dr_2\mathbf{e}_2$ ($\mathbf{e}_1, \mathbf{e}_2$ the basis vectors of the chosen coordinate system), we can write

$$ \begin{aligned} \Delta \mathbf{L}_{AB} = & \begin{pmatrix} \operatorname{sgn}(dr_1)dr_1 & \operatorname{sgn}(dr_1)dr_2 & 0 & 0 \\ 0 & 0 & \operatorname{sgn}(dr_2)dr_1 & \operatorname{sgn}(dr_2)dr_2 \end{pmatrix} \\ & \times \nabla \mathbf{v}|_{4 \times 1}, \end{aligned} \qquad (2.6) $$

where $\Delta \mathbf{L}_{AB}$ is the column vector of the variation rates between stations A and B:

$$ \Delta \mathbf{L}_{AB} = \begin{pmatrix} \Delta L_E \\ \Delta L_N \end{pmatrix}_{AB}, \qquad (2.7) $$

and $\nabla \mathbf{v}|_{4 \times 1}$ is the gradient reshaped into a column matrix of four lines for the sake of the inversion:

$$ \nabla \mathbf{v}|_{4 \times 1} = \begin{pmatrix} \nabla \mathbf{v}_{11} \\ \nabla \mathbf{v}_{12} \\ \nabla \mathbf{v}_{21} \\ \nabla \mathbf{v}_{22} \end{pmatrix}. \qquad (2.8) $$

We can slightly simplify eq. (2.6) by letting

$$ s = \operatorname{sgn}(dr_1) \operatorname{sgn}(dr_2), \qquad (2.9) $$
then by defining

$$ \mathbf{g} = \begin{pmatrix} s|dr_1| & s|dr_2| & 0 & 0 \\ 0 & 0 & s|dr_1| & s|dr_2| \end{pmatrix}. \qquad (2.10) $$

In this way and by using the handier notation $\nabla v$ for $\nabla v|_{4 \times 1}$, eq. (2.6) finally reads

$$ \Delta \mathbf{L}_{AB} = \mathbf{g} \cdot \nabla v, \qquad (2.11) $$

and we can already identify the data matrix $\Delta \mathbf{L}_{AB}$, the theory matrix $\mathbf{g}$ and the model matrix $\nabla v$ for two stations infinitely close together.

The ultimate step consists in extending the previous results to any network whose stations are separated by finite distances. To do this, we just have to integrate eq. (2.11) along a path $AB_i$ between the stations A and B of the ith couple, for all i in [1, $n(n-1)/2$], with n being the number of stations. Thus, we have

$$ \Delta \mathbf{L}_{AB_i} = \int_A^B \mathbf{g}_i(\mathbf{r}) \cdot \nabla \mathbf{v}|_r, \qquad (2.12) $$

with the index i denoting the ith couple. In the discrete case, with a grid made of m numbered squares, if C is the set of the indices of the squares crossed by the path $AB_i$, we have

$$ \Delta \mathbf{L}_{AB_i} = \sum_{k \in C} \mathbf{g}_i k \cdot \nabla \mathbf{v}|_k, \qquad (2.13) $$

which is strictly equivalent to the matrix notation

$$ \Delta \mathbf{L}_{AB_i} = (\mathbf{g}_1 \quad \cdots \quad \mathbf{g}_m)_i \cdot \begin{pmatrix} \nabla \mathbf{v}|_1 \\ \vdots \\ \nabla \mathbf{v}|_m \end{pmatrix}, \qquad (2.14) $$

if $k \notin C$ implies that $\mathbf{g}_k$ is the $2 \times 4$ null matrix. Eq. (2.14) is valid for an unique couple of stations; we just have to vertically concatenate the matrices $\Delta \mathbf{L}_{AB_i}$ and $(\mathbf{g}_1 \quad \cdots \quad \mathbf{g}_m)_i$, for all couples i in the network. Finally, with $\ell = n(n-1)/2$

$$ \begin{pmatrix} \Delta \mathbf{L}_{AB_1} \\ \vdots \\ \Delta \mathbf{L}_{AB_\ell} \end{pmatrix} = \begin{pmatrix} \mathbf{g}_{11} & \cdots & \mathbf{g}_{1m} \\ \vdots & \ddots & \vdots \\ \mathbf{g}_{\ell 1} & \cdots & \mathbf{g}_{\ell m} \end{pmatrix} \cdot \begin{pmatrix} \nabla \mathbf{v}|_1 \\ \vdots \\ \nabla \mathbf{v}|_m \end{pmatrix}. \qquad (2.15) $$

By making profit of the matrix notation, we can write eq. (2.15) in a very compact form. In fact, by letting, respectively, the data, theory and model matrices be

$$ \mathbf{D} = \begin{pmatrix} \Delta \mathbf{L}_{AB_1} \\ \vdots \\ \Delta \mathbf{L}_{AB_\ell} \end{pmatrix}, \quad \mathbf{G} = \begin{pmatrix} \mathbf{g}_{11} & \cdots & \mathbf{g}_{1m} \\ \vdots & \ddots & \vdots \\ \mathbf{g}_{\ell 1} & \cdots & \mathbf{g}_{\ell m} \end{pmatrix} \quad \text{and} \quad \mathbf{M} = \begin{pmatrix} \nabla \mathbf{v}|_1 \\ \vdots \\ \nabla \mathbf{v}|_m \end{pmatrix}, \qquad (2.16) $$

eq. (2.15) reads

$$ \mathbf{D} = \mathbf{G} \mathbf{M}. \qquad (2.17) $$

Due to the 1-D indexation described above, $\mathbf{D}$ and $\mathbf{M}$ are column matrices, so it is possible to introduce covariance matrices on the data and a priori model. We use a null a priori model, that is, at any point, velocity gradients are zero. The covariance matrix of the data is assumed diagonal (independent of all the measurements) and each diagonal element is the square of the uncertainty of the corresponding measure. Thus, the inversion of the direct problem

given by eq. (2.17) follows the method of Tarantola (2005) for a linear problem that leads to the solution

$$ \mathbf{M} = \mathbf{M}_0 + (\mathbf{G}'\mathbf{C}_D^{-1}\mathbf{G} + \mathbf{C}_M^{-1})^{-1}\mathbf{G}'\mathbf{C}_D^{-1}(\mathbf{D} - \mathbf{G}\mathbf{M}_0), \qquad (2.18) $$

with $\mathbf{M}_0$ the a priori model and $\mathbf{C}_D$ and $\mathbf{C}_M$ the covariance matrices on the data and the a priori model, respectively.

The matrix $\mathbf{C}_M$ is used to regularize the inversion through a damping parameter $d_c$ that defines the smoothing of the inversion over the grid. This matrix is based upon the distance matrix that gives the distance between each pair of cells in the grid. We assume that the farther a cell will be from another one, the less the strain rate will be correlated between these two cells. To achieve this goal, the coefficient $(i, j)$ of the covariance matrix on the a priori model is defined as

$$ (C_M)_{ij} = \sigma_0^2 \exp \left( \frac{d_{ij}}{d_c} \right), \qquad (2.19) $$

with $d_{ij}$ the coefficient $(i, j)$ of the distance matrix, that is, the distance between the ith and jth cell of the grid, $d_c$ being the damping parameter defined above and $\sigma_0$ the uncertainty on the a priori model.

In the next section, the method presented here is called STIB method for strain tensor from inversion of baselines. The STIB method is a simplified version of the SNM. We call the strain rate tensor the SRT.

# 3 APPLICATION TO IRAN

## 3.1 GPS network and data

Here, we apply the STIB method to the GPS velocity field covering the Arabia-Eurasia collision published by Nocquet (2012). He combines recently published GPS results in order to derive a geodetic horizontal velocity field consistent at the scale of the Mediterranean basin and its surrounding Alpine belts (Fig. 1). Nocquet (2012) gives the velocity field as supplementary material. In this study, we use the data acquired from 1999 to 2010 from longitude 40° to 63°E and latitude 20° to 45°N. Most of the GPS velocities result from campaign measurements (Vernant et al. 2004; Bayer et al. 2006; Djamour et al. 2010), preventing any computation of the vertical velocities. Therefore, we do not include the vertical component of the velocity and we only estimate the horizontal components of the strain. The distribution of velocity data is not homogeneous over the region under study, with a high density in NW Iran and the Caucasus, in the central part of the Alborz (around Tehran) and at the transition from the Zagros to the Makran. In order to apply the STIB method, we use the baselines of all the paired stations of the network, resulting to a baseline density distribution plotted in Fig. 3(a). The grid used for the inversion of the SRT corresponds to cells of 1° in latitude and longitude.

## 3.2 Synthetic test

First, we apply the STIB method to a synthetic model. Our synthetic test reproduces the Iranian network configuration. The a priori diagonal covariance matrix for the model is taken with a constant value of $1.10^{-08}$ yr$^{-1}$, which corresponds to the expected amplitude of the deformation. The diagonal element of the resolution matrix of the inverse problem is shown in Fig. 3(b). At first order, the resolution is clearly dependent on the ray density (Fig. 3a). However, a more detailed analysis shows that the well-resolved part of the model,
<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a) Normalized ray density. The ray density consists of the sum of the length of the baselines crossing the cell. The value is normalized between 0 (no ray) and 1 (maximum coverage). (b) Diagonal term of the resolution matrix.</figcaption></figure>

roughly corresponding to a triangle covering all of Iran except the eastern regions, is wider than the region of high ray density. This underlines that the resolution depends not only on the baseline density but on the baseline-crossing pattern as well. For two cells with the same number of baselines, the resolution is higher when the number of baseline crossings is high. This is obviously the case for the region around the Straight of Hormuz, where the baseline density is high but where the resolution is lower than in the northwestern part of Iran, due to the fact that most of the baselines are parallel, more or less northwest-southeast.

Inverted SRT has been estimated first without noise on the synthetic velocities and second adding a Gaussian noise at 1 mm yr⁻¹ on the synthetic velocities. Synthetic and inverted SRTs are plotted in Fig. 4(a), drawing strain amplitude along the main axes, using the results of the eigenvalues and eigenvectors computation. Fig. 4(b) shows the histogram of the difference of the azimuth of the principal axis of the strain tensor in the two synthetic inversions (with and without noise). Both the synthetic and the noiseless inverted models are quite similar in terms of shape and amplitude, that is, positive and negative values are correctly located. Without noise, 75 per cent of the differences of the azimuth are smaller than 3° (respectively, 10° with noise).

Small discrepancies are observed along the boundaries of the models where the resolution is poor. Introducing the noise in the velocity field, we observe some differences. The major discrepancy is related to the amplitude. If the maximum of the amplitude is broadly retrieved at the correct location, the amplitude is generally underestimated. This is a typical problem for inverse problem resolution: the inversion process drives to a minimum norm solution. This problem is well known, for example, for most of the synthetic tests made in seismic tomography inversion (Evans & Achauer 1993; Lévque & Masson 1999). However, in the NW of the network, where the ray density is high and consequently where the resolution is high, the SRT is particularly well retrieved in tectonic style (strike-slip, extension or compression), direction and amplitude.

### 3.3 Application to real data

We now apply the STIB method to the real GPS velocity field in Iran. In Fig. 5, the SRT is obtained at each degree, corresponding to the cell size used in the inversion.

At first glance, the main features observed using the STIB method are:

(i) the significant shortening related to the Makran subduction,

(ii) the strike-slip pattern in front of the straight of Hormuz related to the NS-trending faults separating the Central Iran Plateau and the Lut Block,

(iii) the shortening within the Zagros and the Alborz mountain chains,

(iv) the rigid behaviour of the western Central Iran Plateau, which contrasts with its eastern part,

(v) a complex pattern of the strain field in NW Iran.

The three first features have been roughly described by Masson *et al.* (2007) from a SRT representation based on a triangle method applied in the same area. However, it is compelling to note that the STIB method gives a finer description of the strain than the triangle method, since the STIB method allows us to describe the SRT at a scale which is significantly smaller than the interstation distance. Therefore, we are now able to describe the SRT with short
<figure><img src="image_5.png" /><figcaption>**Figure 4.** (a) Synthetic tests. The middle map is the synthetic strain rate field that must be recovered through the inversion of the corresponding synthetic velocity field (we put it at the middle for easier comparison with the two other maps). The stars are the GPS stations. On the left is the result of the inversion for noiseless velocities (the data are the original synthetic velocities). On the right is the result of the inversion after a Gaussian noise at 1 mm yr⁻¹ has been added to the synthetic velocity field. The colour scale gives the trace of the strain rate tensor (strain yr⁻¹). (b) Histograms of the difference of the azimuth of the principal axis of the strain tensor in the two synthetic inversions (with and without noise). The three red lines indicate, respectively, the 50, 75 and 90 per cent percentiles.</figcaption></figure>

wavelength variability, for example, in NW Iran and in Central Iran, which were intrinsically included in the data but not visible using the triangle method.

## 3.4 Robustness of the inversion

As underlined above, the main limitation of the triangle or gridded methods comes from the significant effect of an erroneous velocity vector on the final results. One of the objectives of the STIB method is to reduce this effect using the information provided by the full set of baselines, including baselines between two stations with correct velocities passing close to the site with the wrong velocity. To illustrate this point, we introduce an inaccurate velocity vector in the Iranian velocity field and we compare the SRT obtained using the triangle method and the STIB method. Two tests are shown in Fig. 6. In the first, we increase the northward component of velocity at ARD1 of 5 mm yr⁻¹. This introduces compression north of ARD1 and extension south of ARD1. However, the effect is limited

to the cells surrounding ARD1, independently of the location of the closest stations. In the second, we do the same thing but with the eastward component of velocity. Using the triangle or gridded methods, the false velocity at one site strongly impacts the SRT in the area around this site. If the network is sparse, the perturbed zone can be very large. In our example, due to the fact that ARD1 is an isolated site in the centre of Iran, this zone could correspond to a circle of 300 km of radius. Using the STIB method, the impact of the velocity perturbation in the area surrounding ARD1 is strongly reduced. However, the impact of the erroneous input velocity is not confined to the area as in the triangle or gridded method but the strain amplitude variations out of this area remain very small and can be considered negligible. This is true in both cases.

## 4 DISCUSSION

The active tectonics in Iran results from the north-south convergence of the Arabian and Eurasian plates (Jackson & McKenzie
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Eigenvectors of the inverted strain rate tensor over the Iranian domain using the STIB method. The colour scale gives the trace of the strain rate tensor (strain/yr).</figcaption></figure>

1984). Within Iran, most of the convergence is accommodated by compressive tectonic structures forming the major mountain belts (Zagros, Alborz) and by large strike-slip faults surrounding rigid blocks (Central Iran, Lut and the southern Caspian sea). East of longitude 058°E, most of the shortening occurs in the Makran subduction zone, with an horizontal velocity related to the fixed Eurasia plate of 19.5 mm yr⁻¹ decreasing to 6 mm yr⁻¹ in NE Iran. West of 058°E, the deformation is distributed over several folds and thrust belts. At the longitude of Tehran, the Zagros and the Alborz mountain ranges accommodate 7.5 and 6 mm of shortening per year, respectively (Fig. 1). From the results using the STIB method, we have here the opportunity to give a fine description of the SRT variations over Iran, with a high spatial resolution, larger than the distance between the stations. In Fig. 7, the STR is plotted together with the focal mechanisms from 1909 to 2012. At first glance, it appears that the style of the seismicity is, in general, consistent with the deformation defined by the SRT, with compression in the mountain ranges (Zagros, Alborz, Kopet-Dag and Caucasus) and

in the Makran subduction zone and strike-slip around the stable areas (border faults of the Lut block or Main Recent Fault along the northern edge of the Zagros). However, this general sketch is frequently too simple to explain both the SRT and the seismicity. The SRT is more complex than expected from previous studies. For instance, Masson *et al.* (2007) proposed a very simple sketch of the SRT, which is not fully in agreement with the SRT proposed here. Also, the seismicity appears heterogeneous. This can be examined in details for the main regional structures of Iran.

*The Makran subduction zone:* The Makran belt is the subaerial section of an extensive accretionary prism. Its frontal thrust is located about 70 km seaward. The velocity field indicates a large shortening between the Omani and the Iranian coasts, related to the subduction of oceanic crust of the Gulf of Oman below the Makran. The shortening is mainly located within the Gulf of Oman, with an abrupt decrease of the shortening onland, that is, north of the coastline. This suggests that the Makran massif concentrates only a small part of the compressive deformation associated with
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Robustness tests. The middle map is a zoom of Fig. 5, that is, the inverted strain rate field from the original data in the central part of Iran (we put it at the middle for easier comparison with the two other maps). On each side are the computed strain rate fields after a perturbation in the velocity at station ard1 (red arrow). On the left, the perturbation is +5 mm yr⁻¹ on the north component of velocity, while on the right, the perturbation is +5 mm yr⁻¹ on the east component. The colour scale gives the trace of the strain rate tensor (strain yr⁻¹).</figcaption></figure>

the Arabia-Eurasia convergence. However, these results should be analysed with caution since this area located on the edge of the geodetic network certainly suffers from a poor resolution. Unlike the eastern Pakistan side of the Makran subduction zone which was stroked by large earthquakes in 1945 (M 8.1) and 1947 (M 7.3) (Smith *et al.* 2013), the western Iranian Makran subduction zone is only affected by small-magnitude seismicity over the last century. Particularly, no reverse fault seismicity is related to the large shortening within the Gulf of Oman. The seismogenic potential of this area is therefore not well known and is the purpose of current researches (Smith *et al.* 2013). The recent 2013 April 16 Mw 7.8 Khash earthquake is a deep normal-faulting earthquake within the subducting plate, far from the trench, which cannot be directly related to the surface processes and the SRT.

A abrupt change of the style and amplitude of the deformation occurs west of 57° of longitude, which corresponds to the location of the transition from the oceanic crust of the Gulf of Oman to the continental crust.

**The Zagros–Makran transition zone:** The Strait of Hormuz area (Fig. 8) corresponds to the transition from the continental collision of the Zagros (west) to the oceanic subduction of beneath the Makran wedge (east). Analysing a regional GPS velocity field, Bayer *et al.* (2006) and Peyret *et al.* (2009) show that the difference of shortening rate is mainly accommodated along the NNW–SSE trending reverse right-lateral Zendan–Minab–Palami (ZMP) fault system (Fig. 8). This system is extended northward by large faults, which mark the western border of the Lut block. We specifically compute a model for this zone with a grid spacing of 0.5°. Around the Straight of Hormuz, the strain tensor shows an interesting evolution from east to west: the direction of the maximum strain switches from NNE–SSW along the Oman Gulf to NE–SW along the coast of the Musandam peninsula. This change is very abrupt, in less than 50 km. In this region, according to Bayer *et al.* (2006), rotations are very large due to the accommodation of a large velocity change within a very short distance, in front of the Musandam Peninsula which indents the Iranian coast. Rotations, plotted in Fig. 8(b), are large east of the ZMP fault system and in front of the Musandam Peninsula.

The SRT shows NS variations along the ZMP fault system. From south to north, we observe that the right-lateral strike-slip component switches to a more compressive NE-SW fault system. The dense seismicity north of the ZMP is consistent with shortening, while the sparse seismicity along the ZMP does not allow us to conclude about the consistency. A recent seismic crisis (May–June 2013) should help to answer this question.

**Zagros:** The Zagros thrust and fold belt extends for more than 1500 km in NW–SE direction from the ZMP fault system in southern Iran to eastern Turkey. This belt results from the closure of the Neo-Tethyan ocean involving a northeast-dipping subduction below the Iranian microcontinent. From a regional GPS network, Walpersdorf *et al.* (2006) propose that the deformation is partitioned in northwest Zagros into a 3–6 mm yr⁻¹ shortening perpendicular to the axis of the mountain belt (i.e. by purely reverse faults), and a 4–6 mm yr⁻¹ lateral motion along NW–SE-trending dextral faults. In southeast Zagros, the shortening of 8 ± 2 mm yr⁻¹ occurs along perpendicular to the main features of the mountain belt. The SRT shows a pure shortening in the south-eastern part of the belt (Fig. 7), close to the Straight of Hormuz. This pattern is consistent with the seismicity, since most of the focal mechanisms show pure compression. In the other part of the belt, the SRT is more variable with both compression and strike-slip components. Interestingly, the dextral strike-slip north–south Kazerun fault is consistent with the results of the STIB method.

**Central Iranian block:** Mostly from seismological observations, Jackson & McKenzie (1984) suggest that the central Iranian block, the elongated plateau parallel to the northern limit of the Zagros range, can be considered as rigid. Using five GPS stations, Vernant *et al.* (2004) estimated an Euler pole for the central Iranian block. The residual velocities for those five sites are below the 1σ uncertainty, indicating that the rigid description of the central Iranian block is appropriate, since an internal deformation occurs at less than 2 mm yr⁻¹. The SRT obtained from the STIB method, due to the inversion process which does not takes only into account the velocity of the nearest GPS stations, does not indicate a null strain within the Central Iranian block. Even if the amplitude of the strain remains very small and the interpretation has to be cautious, we
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Strain tensor and focal mechanisms from 1909 to 2012 (grey: seismicity from Jackson *et al.* 1995, black: Harvard catalogue). The bold black line indicates the region where the diagonal term of the resolution matrix is greater than 0.5.</figcaption></figure>

note a slight right-lateral strike-slip component in the strain tensor at the vicinity of both the Deshir and Anar faults (Fig. 9). This result, supported by the high resolution (Fig. 3b) obtained in this area, is consistent with previous geomorphologic and palaeoseismologic studies performed along the Anar and Deshir faults which estimate on these two faults displacements of about 1 mm yr⁻¹ and attest of a existence in the past of large earthquakes (magnitude 7). Meyer *et al.* (2006), using geological and geomorphologic offsets to constrain the total displacement, estimate the slip rate of the right-lateral Deshir fault close to 2 mm yr⁻¹. This rate has been slightly reduced by Le Dortz *et al.* (2011) in the range 0.9–1.5 mm yr⁻¹ using 10Be cosmic ray exposure and optically stimulated luminescence. Meyer & Le Dortz (2007) estimate also a slip rate averaged over the Holcene of 0.5–0.7 mm yr⁻¹ along the Anar fault. To do this, they used high-resolution SPOT5 (pixel size 2.5 m) images to measure recent cumulative offsets. Le Dortz *et al.* (2009) confirm this result finding a rate of 0.8 mm yr⁻¹ along the same fault, once again using 10Be cosmic ray exposure and optically stimulated luminescence. These rates along the Deshir and Anar faults have been confirmed by several palaeoseismological studies. Foroutan *et al.* (2012) show that the Anar fault experiences three large ($M_w \approx 7$)

earthquakes during the Holocene, while Nazari *et al.* (2009) and Fattahi *et al.* (2010) show the occurrence of large earthquakes on the Deshir fault. These results are in agreement with the sparse seismicity observed east of the Anar fault.

*NW Iran–Caucasus–East Turkey:* For this region as well, we compute a dedicated model using a grid spacing of 0.5° (Fig. 10). This region is the northernmost part of the Arabia-Eurasia collision zone. It is characterized by complex active tectonics, including both N–S compressive structures (E–W thrusts and folds) and E–W extensional structures (N–S normal faults and dykes) together with important Neogene to Quaternary volcanism, and NE–SW left-lateral and NW–SE right-lateral strike-slip faults (Rebai *et al.* 1993). These structures are all imbricated and can be observed at different scales. The western part of the zone, the Armenian plateau centred on the Van lake in eastern Turkey, is a large zone of compression, which gradually gives way to the rigid Anatolian plateau limited by the North and East Anatolia Faults. The northern border of the zone corresponds to the Caucasus. The south-eastern part of the zone, from the Ourmia lake to the Caspina sea in NW Iran, is mainly characterized by the large right-lateral strike-slip Tabriz fault. In spite of a global N–S shortening from Arabia to
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** (a) Same as Fig. 7 for the Zagros–Makran transition zone. Hor. Str.: Hormuz Strait, Mus. Pen.: Musandam Peninsula, ZMP: Zendan–Minab–Palami fault system. (b) Rotations obtained for the Straight of Hormuz region. The differential movement east and west of the Musandam peninsula induces large block rotation. This pattern of rotation is very similar to the results obtained by Bayer *et al.* (2006).</figcaption></figure>

Eurasia, the SRT underlines the complexity of the tectonics of the region. It is strongly variable, in terms of both style and amplitude. The example of the Caucasus illustrates perfectly this observation. At large scale, a decreasing shortening rate is observed from east to west. However, at small scale, we observe pure shortening south of the eastern Caucasus in the Kura basin in Azerbaijan, strike-slip south of the central part of the chain (from 44° to 46°) and once again pure compression in the western inner part of the chain. Another example is given by the Talesh block, which appears rigid between the Tabriz fault and the Kura basin.

A SRT covering the eastern Turkey and western Caucasus has been published by Özeren & Holt (2010) using another velocity field and a gridded approach. It extends eastward up to the Ourmia Lake, although they used only few Iranian GPS sites. At large scale, especially west of the Van Lake and around the Karliova triple junction, both SRT have significant similarities.

When the seismicity is dense, the deformation pattern is generally consistent with the seismicity. Most of the earthquakes have focal mechanisms in agreement with the SRT. This is the case in the western inner part of the Caucasus or in Eastern Turkey (North and East Anatolian faults zone). However, this is not true when the seismicity is sparse. In this case, the seismicity appears more heterogeneous in terms of mechanisms and not always in agreement with the SRT. This is observed in the Kura basin (while a high strain rate is observed) or in Northern Armenia.

The SRT allows a fine description of the tectonics based on geodetic data, which complements the seismic data.

## 5 CONCLUSION

In order to avoid the strong dependency of the strain field on the quality of the velocity field, we have used a method of computation of the strain that takes into account all the information included in the GPS data. To do so, we use the baseline time-series, that is, the variations of distance between all pairs of stations.

The STIB method applied to an Iranian GPS velocity field offers the opportunity to describe the strain distribution with high spatial precision. Therefore, we are able to identify regions where strong deformation occurs (about 1.5 cm for 1500 km, i.e. 10⁻⁸) together with regions where the deformation is very small. New maps of the SRT are produced, drawing the pattern of deformation at global scale for the whole Iran and at regional scale for the NW Iran–Caucasus–East Turkey region, the Zagros–Makran transition zone and the Central Iran region.

The three main results are

(1) a new illustration of the very sharp transition zone from the Makran (subduction) to the Zagros (collision) underlined by the ZMP fault system. This system generates recently a seismic swarm (two earthquakes of Mw 5.5 in May 2013).

(2) the strongly heterogeneous behaviour of the NW Iran–Caucasus–East Turkey region, which presents a large variability of strain patterns, clearly explaining the large heterogeneity of observed structures at the surface (reverse, normal, strike-slip faults, volcanoes).

(3) the unexpected possibility of deformation in the eastern Central Iran consistent with previous geomorphological studies of the Anar and Deshir faults, not proposed from the original velocity field proposed by any previous geodetic studies.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Same as Fig. 7 for the Central Iranian zone.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Same as Fig. 7 for the northern part of the Arabia-Eurasia collision zone including the Caucasus and the eastern Anatolian plateau. NTF: North Tabriz Fault.</figcaption></figure>

## REFERENCES

Aki, K., Christoffersson, A. & Husebye, E., 1977. Determination of the three-dimensional seismic structure of the lithosphere, *J. geophys. Res.*, **82**, 277–296.

Bayer, R. *et al.*, 2006. Active deformation in Zagros–Makran transition zone inferred from GPS measurements, *Geophys. J. Int.*, **165**, 373–381.

Berberian, M. & King, G.C.P., 1981. Towards a paleogeography and tectonic evolution of Iran, *Can. J. Earth Sci.*, **18**(2), 210–285.

Bos, A. & Spakman, W., 2005. Kinematics of the southwestern U.S. deformation zone inferred from GPS motion data, *J. geophys. Res.*, **110**, doi:10.1029/2003JB002742.

Bos, A., Spakman, W. & Nyst, M., 2003. Surface deformation and tectonic setting of Taiwan inferred from a GPS velocity field, *J. geophys. Res.*, **108**, doi:10.1029/2002JB002336.

Bos, A., Usai, S. & Spakman, W., 2004. A joint analysis of GPS motions and InSAR to infer coseismic surface deformation of the Izmit, Turkey earthquake, *Geophys. J. Int.*, **158**, 849–863.
Byrne, D.E., Sykes, L.R. & Davis, D.M., 1992. Great thrust earthquakes and aseismic slip along the plate boundary of the Makran subduction zone, *J. geophys. Res.*, **97**, 449–478.

Djamour, Y. et al., 2010. GPS and gravity constraints on continental deformation in the Alborz mountain range, Iran, *Geophys. J. Int.*, **183**, 1287–1301.

Evans, J. & Achauer, U., 1993. Teleseismic velocity tomography using the ACH method: theory and application to continental-scale studies, in *Seismic Tomography: Theory and Practice*, pp. 319–360, eds Iyer, H.M. & Hirahara, K., Chapman & Hall.

Falcon, N.L., 1974, Southern Iran: Zagros Mountains, in *Mesozoic-Cenozoic Orogenic Belts, Data for Orogenic Studies*, pp. 199–211, ed., Spencer, A.M., Geological Society of London Special Publication.

Fattahi, M. et al., 2010. Refining the OSL age of the last earthquake on the Dheshir fault, Central Iran, *Quat. Geochronol.*, **5**, 286–292.

Foroutan, M., Sébrier, M., Nazari, H., Meyer, B., Fattahi, M., Rashidi, A., Le Dortz, K. & Bateman, M.D., 2012. New evidence for large earthquakes on the Central Iran plateau: palaeoseismology of the Anar fault, *Geophys. J. Int.*, **189**, doi:10.1111/j.1365-246X.2012.05365.x.

Frank, F.C., 1967. Deduction of earth strains from Survey data, *Bull. seism. Soc. Am.*, **56**, 35–42.

Haines, A., Jackson, J., Holt, W. & Agnew, D., 1998. *Representing distributed deformation by continuous velocity fields, Sci. Rep.* **98/5**, Institute of Geological and Nuclear Sciences, Wellington, New Zealand.

Howe, T. & Bird, P., 2010. Exploratory models of long-term crustal flow abbd resulting seismicity across the Alpine-Aegean orogen, *Tectonic*, **29**, TC4023, doi:10.1029/2009TC002565.

Jackson, J. & McKenzie, D., 1984. Active tectonics of the Alpine-Himalayan Belt between western Turkey and Pakistan, *Geophys. J. R. astr. Soc.*, **77**, 185–264.

Jackson, J.A., Haines, J. & Holt, W., 1995. The accomodation of the Arabia-Eurasia plate convergence in Iran, *J. geophys. Res.*, **100**, 15 205–15 219.

Le Dortz, K. et al., 2009. Holocene right-slip rate determined by cosmogenic and OSL dating on the Anar fault, Central Iran, *Geophys. J. Int.*, **179**, 700–710.

Le Dortz, K. et al., 2011. *Dating inset terraces and offset fans along the Dehshir Fault (Iran) combining cosmogenic and OSL methods. Geophys. J. Int.*, **185**, 1147–1174. doi:10.1111/j.1365-246X.2011.05010.x.

Lévque, J.-J. & Masson, F., 1999. From ACH tomographic models to absolute velocity models, *Geophys. J. Int.*, **137**, 621–629.

Masson, F., Chery, J., Hatzfeld, D., Martinod, J., Vernant, P., Tavakoli, F. & Ghafory-Ashtiani, M., 2005. Seismic versus aseismic deformation in Iran inferred from earthquakes and geodetic data, *Geophys. J. Int.*, **160**, 217–226.

Masson, F., Anvari, M., Djamour, Y., Walpersdorf, A., Tavakoli, F., Daignieres, M., Nankali, H. & Van Gorp, S., 2007. Large-scale velocity field and strain tensor in Iran inferred from GPS measurements: new insight for the present-day deformation pattern within NE Iran, *Geophys. J. Int.*, **170**, 436–440.

Meyer, B. & Le Dortz, K., 2007. Strike-slip kinematics in Central and Eastern Iran: estimating fault slip-rates averaged over the Holocene, *Tectonics*, **26**, TC5009, doi:10.1029/2006TC002073.

Meyer, B., Mouthereau, F., Lacombe, O. & Agard, P., 2006. Evidence of quaternary activity along the Deshir Fault, *Geophys. J. Int.*, **164**, 192–201.

Nazari, H. et al., 2009. First evidence for large earthquakes on the Deshir Fault, Central Iran Plateau, *Terra Nova*, **21**, 417–426.

Peyret, M. et al., 2009. Present-day strain distribution across the Minab-Zendan-Palami fault system from dense GPS transects, *Geophys. J. Int.*, **179**, 751–762.

Özeren, M.S. & Holt, W.E., 2010. The dynamics of the eastern Mediterranean and eastern Turkey, *Geophys. J. Int.*, **183**, 1165–1184.

Pollitz, F. & Nyst, M., 2005. A physical model for strain accumulation in the San Francisco Bay Region, *Geophys. J. Int.*, **160**, 302–317.

Rebai, S., Philip, H., Dorbath, L., Borisoff, B., Haessler, H. & Cisternas, A., 1993. Active tectonics in the Lesser Caucasus: coexistence of compressive and extensional structures, *Tectonics*, **12**(5), 1089–1114.

Smith, G.L., McNeill, L.C., Wang, K., He, J. & Henstock, T.J., 2013. Thermal structure and megathrust seismogenic potential of the Makran subduction zone, *Geophys. Res. Lett.*, **40**, 1528–1533.

Spakman, W. & Nyst, M., 2002. Inversion of relative motion data for estimates of the velocity gradient field and fault slip, *Earth planet. Sci. Lett.*, **203**, 577–591.

Tape, C., Musé, P., Simons, M., Dong, D. & Webb, F., 2009. Multiscale estimation of GPS velocity fields, *Geophys. J. Int.*, **179**, 945–971.

Tarantola, A., 2005. *Inverse Problem Theory and Methods for Model Parameter Estimation*, pp. 163–169, SIAM.

Tarantola, A., Ruegg, J.-C. & Lépine, J.-C., 1979. Geodetic evidence for rifting in Afar: a brittle-elastic model of the behaviour of the lithosphere, *Earth. planet. Sci. Lett.*, **45**, 435–444.

Vernant, P. et al., 2004. Present-day crustal deformation and plate kinematics in the Middle East constrained by GPS measurements in Iran and northern Oman, *Geophys. J. Int.*, **157**, 381–398.

Walpersdorf, A. et al., 2006. Difference in the GPS deformation pattern of north and central zagros (Iran), *Geophys. J. Int.*, **167**, 1077–1088.

Wu, Y., Jiang, Z., Yang, G., Wei, W. & Liu, X., 2011. Comparison of GPS strain rate computing methods and their reliability, *Geophys. J. Int.*, **185**, 703–717.