Numerical modeling of stresses and deformation in the Zagros–Iranian Plateau region

The Zagros orogenic system resulted due to collision of the Arabian plate with the Eurasian plate. The region is characterized by ocean–continent subduction and continent–continent collision, and the convergence velocity shows variations from east to west. Therefore, this region shows the complex tectonic stress and a wide range of diffuse or localized deformation between both plates. The in situ stress and GPS data are very limited and sparsely distributed in this region; therefore, we performed a numerical simulation of the stresses causing deformation in the Zagros–Iran region. The deviatoric stresses resulting from the variations in lithospheric density and thickness and those from shear tractions at the base of the lithosphere due to mantle convection were computed using thin-sheet approximation. Stresses associated with both sources can explain various surface observations of strain rates, $S_{{\mathrm{H}}_{\max }}$, and plate velocities, thus suggesting a good coupling between lithosphere and mantle in most parts of Zagros and Iran. As the magnitude of stresses due to shear tractions from density-driven mantle convection is higher than those from lithospheric density and topography variations in the Zagros–Iranian Plateau region, mantle convection appears to be the dominant driver of deformation in this area. However, the deformation in the east of Iran is caused primarily by lithospheric stresses. The plate velocity of the Arabian plate is found to vary along the Zagros belt from the north–northeast in the southeast of Zagros to the northwest in northwestern Zagros, similarly to observed GPS velocity vectors. The output of this study can be used in seismic hazards estimations.

1 Introduction

The Zagros Mountains are a part of the Alpine–Himalayan belts that originated due to the Arabian plate colliding with the southern boundary of the Eurasian plate. This collision resulted from the closing of the Neotethys Ocean and formed the Zagros fold-and-thrust belt (Agard et al., 2005, 2011; Alavi, 1980; Mouthereau et al., 2012). The Zagros Mountains extend from the eastern part of Anatolia for over 1500 km in the NW–SE direction until the Makran subduction zone, showing large-scale diffuse deformation. Despite the first-order characteristics of the active deformation and present-day kinematics of the Zagros orogen being relatively well understood (Allen et al., 2011; Le Dortz et al., 2009; Reilinger et al., 2010; Vernant et al., 2004; Walker, 2006), the timing of the collision is debated. The timing of the collision ranges from Cretaceous (Alavi, 1994; Mohajjel and Fergusson, 2000) to Miocene (Berberian and King, 1981) or Eocene (Allen and Armstrong, 2008; Jolivet and Faccenna, 2000). However, there has been an increasing consensus on Late Eocene to Oligocene for the onset of collision (Jolivet and Faccenna, 2000; Agard et al., 2005, 2011; Vincent et al., 2005; Ballato et al., 2011; Mouthereau et al., 2012; Koshnaw et al., 2019). The Arabia–Eurasia collision zone is a tectonically active region, where ongoing convergence is accommodated by distributed shortening across the Zagros Mountains and the northern and eastern margins of the Iranian Plateau and the southern Caspian Sea. The rate of convergence of Arabia relative to Eurasia also varies significantly, decreasing from 36 mm yr⁻¹ in the east to 16 mm yr⁻¹ in the west (Fig. 1). The diverse structures, tectonic history, and convergence velocity variations in the Zagros–Iran plateau region lead to variable tectonic stresses and deformations, thus making it the focus of various geophysical, geological, and geodesy studies (Engdahl et al., 2006; Hatzfeld et al., 2010; Khorrami et al., 2019; Masson et al., 2006; Tunini et al., 2016, 2017). Based on earthquake focal mechanisms, fault slips, and GPS velocities, the Zagros–Iran region has been categorized as a highly seismic region; thus, a better constraint on stresses and deformation in this region may be helpful in disaster mitigation studies.

Figure 1 Tectonic overview of central Eurasia. Abbreviations: CF – Chaman Fault, MSZ – Makran subduction zone, MZT – Main Zagros Thrust, HZF – High Zagros Fault, MFF – Mountain Front Fault, SSZ – Sanandaj–Sirjan Zone, UDMA – Urumieh–Dokhtar Magmatic Arc, and MRF – Main Recent Fault.

Generally, tectonic stress refers to the forces acting on the Earth's crust that cause it to deform or undergo changes, and it is classified by first, second, and third order on the spatial scale (Heidbach et al., 2007; Zoback, 1992). The first-order stresses originate due to plate boundary forces like ridge push, slab pull, and continental collision, and second-order stresses originate due to rifting, isostasy, and deglaciation. Moreover, third-order stresses are caused by local sources like interaction between fault systems, topography, and density heterogeneity. Therefore, to understand the origin of these stresses, in situ stress measurements are done using the focal-mechanism inversion, wellbore breakouts, hydraulic fracturing, and overcoring and are compiled under the World Stress Map project. However, in situ stress data are sparsely distributed and limited, so numerical modeling plays an important role in understanding the kinematics and dynamics of the Zagros–Iran region. Numerical modeling of tectonic stresses and deformation is generally conducted using two approaches: (1) using 2D and 3D geometrical structures (obtained from geophysical studies), plate boundary forces (like ridge push, slab pull, and continent collision forces), and rheological properties (like Young’s modulus, Poisson's ratio, viscosity, and density) (Coblentz and Sandiford, 1994; Dyksterhuis and Müller, 2008; Koptev and Ershov, 2010; Richardson et al., 1976; Yadav and Tiwari, 2018) or (2) considering gravitational potential energy and shear tractions from mantle convection with thin-sheet approximation (Bird, 1998; Flesch et al., 2001; Ghosh and Holt, 2012; Lithgow-Bertelloni and Guynn, 2004; Singh and Ghosh, 2020).

There are various studies that have tried to investigate present-day stresses and deformations of the Zagros–Iranian Plateau region using focal-mechanism inversions, GPS data, and numerical modeling. The stresses were computed through the inversion of focal mechanisms in areas like the Zagros fold-and-thrust belt (Nouri et al., 2023; Sarkarinejad et al., 2018; Yaghoubi et al., 2021), the Zagros–Makran transition zone (Ghorbani Rostam et al., 2018), western Zagros (Navabpour et al., 2008), northwestern Iran to southeastern Turkey (Mokhoori et al., 2021), the northeastern Lut block, eastern Iran (Rashidi et al., 2022; Raeesi et al., 2017), and the southern Caspian Sea (Jackson et al., 2002). The GPS studies also provided constraints on the present-day deformation in the Zagros–Makran transition zone (Bayer et al., 2006), Makran subduction zone (Frohling and Szeliga, 2016), Iran (Khorrami et al., 2019; Masson et al., 2006, 2007; Vernant et al., 2004; Walpersdorf et al., 2014), and the Nubia–Arabia–Eurasia plate system (Reilinger and McClusky, 2011). Sobouti and Arkani-Hamed (1996) studied the large-scale tectonic processes of the region and reproduced observed faulting patterns by considering highly rigid central Iran and the southern Caspian Sea using a viscous thin-sheet approximation. On the other hand, Md and Ryuichi (2010) used finite-element modeling (FEM) to analyze the neotectonic stress field of Zagros and the adjoining area and showed N–S- and NNE–SSW-oriented $S_{{\mathrm{H}}_{\max }}$ in Lurestan and the eastern Zagros simple folded belt, whereas $S_{{\mathrm{H}}_{\max }}$ was NW–SE-oriented around the Main Recent Fault (MRF) and in the northern High Zagros Fault (HZF). Further, the kinematic model by Khodaverdian et al. (2015) provided constraints on the fault slip rates, plate velocities, and seismicity of the Iranian Plateau. Most of the deformation studies done in this region focus on different tectonic fragments of the Arabia–Eurasia collision zone. Moreover, the previous studies do not include the role of shear tractions associated with mantle convection in affecting the deformation and stresses in the Zagros–Iran regions.

In this study, we investigate the stress and deformation in the entire Zagros–Iranian Plateau region to constrain the forces acting in this region with gravitational potential energy (GPE) and shear traction of mantle tractions. We will use a thin viscous sheet model based on Flesch et al. (2001) to compute various deformation parameters such as deviatoric stresses, strain rates, most compressive principal stress ($S_{{\mathrm{H}}_{\max }}$), and plate velocities within the Zagros–Iran region.

2 Tectonics and geology

The evolution of the Zagros mountain belt is a direct consequence of the continental collision between the Arabian and Eurasian plates. Zagros is located on the northeastern margin of the Arabian plate, trending in the southwest direction (Fig. 1). It is bounded by the Main Zagros Thrust (MZT) in the northeast, while it joins the Tauras Mountains in southern Turkey in the northwest. In the southeast, the N–S-trending Minab–Zandan fault zone separates Zagros from the Makran range. In outer Zagros are the young folded mountains in the southwest of the orogeny (Falcon, 1974; Sattarzadeh et al., 2002). The High Zagros Fault (HZF) separates the highly deformed metamorphic rocks of inner Zagros from the simply folded mountains of outer Zagros (Hatzfeld and Molnar, 2010; Hatzfeld et al., 2010). Inner Zagros is bounded by the MZT in the northeast and is dominated by thrust faulting, possibly due to compression during the Late Cretaceous (Alavi, 1980). Northwestern Zagros is separated from central Zagros by a north–south-trending strike–slip zone of deformation known as the Kazerun fault system (KFS) (Authemayou et al., 2005).

The Zagros Mountains were formed between ∼35 and ∼23 Ma due to the convergence of the Arabian platform beneath the central Iranian crust (Agard et al., 2005; Ballato et al., 2011; Mouthereau et al., 2012). The Arabian plate moved towards Eurasia with a plate velocity of 22–35 mm yr⁻¹ (DeMets et al., 1990; McClusky et al., 2000; Jackson et al., 2002; McQuarrie et al., 2003; Reilinger et al., 2006) in N–S to NNE direction. The convergence of the two rigid plates of Arabia and Eurasia lead to a zone of widespread deformation in the form of the high plateaus of Iran. The Iranian Plateau extends from the Caspian Sea and the Kopet Dagh range in the north to the Zagros Mountains in the west. The Iranian Plateau is bounded by the Persian Gulf and Hormuz Strait in the south and the country's political borders on the eastern side. Several tectonic processes such as intracontinental collisions, subduction along the Makran, and the transition from the Zagros fold-and-thrust belt to the Makran subduction zone contribute to the complex tectonics of the Iranian plateau.

During the last few decades, various geophysical studies (receiver function, deep seismic, GPS, and tomographic) have been carried out in the Zagros–Iran region to investigate the structure and deformation in this region. Southeastern Zagros accommodates the convergence between Arabia and Eurasia by pure shortening that occurs through high-angle (30–60^∘) reverse faults that are perpendicular to the belt (Hessami et al., 2006; Irandoust et al., 2022; Walpersdorf et al., 2006). On the other hand, oblique convergence in central and northern Zagros is partitioned into a strike–slip component that is accommodated by the MRF and shortening occurring across the belt (Jackson et al., 2002; Talebian and Jackson, 2002). Zagros is separated from the Makran subduction zone (MSZ) by the Minab–Zendan–Palami (MZP) fault (54–58^∘ E), which is a right-lateral strike–slip fault (Bayer et al., 2006). East of MZP shows significant shortening that is accommodated through the subduction in the MSZ. Due to the difference between convergence rates, a shearing occurs in eastern Iran, which is accommodated by the N–S-trending faults bounding the Lut block. In northern Iran, the fold-and-thrust belt of Alborz accommodates a quarter of the Arabia–Eurasia convergence (Irandoust et al., 2022). The oblique convergence in eastern Alborz is also partitioned into shortening at the southern boundary and a left-lateral component across the mountain belt (Irandoust et al., 2022; Khorrami et al., 2019; Tatar and Hatzfeld, 2009). The Alborz Mountains extend into Talesh in the west, which shows thrust faulting on nearly flat faults. The Kopet Dagh range in the northeast accommodates the Arabia–Eurasia convergence through N–S shortening on major thrust faults in the south.

3 Modeling

3.1 Equations

To model the present-day stresses causing deformation in the Zagros–Iranian Plateau due to the Arabia–Eurasia collision, we solve the force balance equation considering the thin sheet approximation.

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial {\mathit{\sigma }}_{ij}}{\partial x_j}}+\mathit{\rho }{g}_i=\mathrm{0}.\end{array}$$

Here σ_ij, x_j, ρ, and g_i indicate the ijth component of the total stress tensor, the jth coordinate axis, the density, and the acceleration due to gravity, respectively (England and Molnar, 1997; Ghosh et al., 2013b).

In the above equation, total stress σ_ij is substituted by deviatoric stress using the following relation:

$$\begin{array}{}\text{(2)}& {\mathit{\tau }}_{ij}={\mathit{\sigma }}_{ij}-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\mathit{\sigma }}_{kk}{\mathit{\delta }}_{ij}.\end{array}$$

In the above equation, the Kronecker delta and mean stress are denoted by δ_ij and $\frac{\mathrm{1}}{\mathrm{3}}{\mathit{\sigma }}_{kk}$, respectively. The force balance equation (Eq. 1) is integrated up to the base of the lithospheric sheet (L), resulting in the following full horizontal force balance equations:

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{\mathit{\tau }}}_{xx}}{\partial x}}-{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{\mathit{\tau }}}_{zz}}{\partial x}}+{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{\mathit{\tau }}}_{xy}}{\partial y}}=-{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{zz}}{\partial x}}+{\mathit{\tau }}_{xz}\left(L\right),\text{(4)}& {\displaystyle }& {\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{\mathit{\tau }}}_{yx}}{\partial x}}+{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{\mathit{\tau }}}_{yy}}{\partial y}}-{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{\mathit{\tau }}}_{zz}}{\partial y}}=-{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{zz}}{\partial y}}+{\mathit{\tau }}_{yz}\left(L\right).\end{array}$$

In Eqs. (3) and (4), the over bars indicate integration over depth. Both equations (Eqs. 3 and 4) contain the first term representing horizontal gradients of GPE per unit area on the right-hand side. On the other hand, the shear tractions at the lithosphere base (L) arising due to mantle convection are denoted by τ_xz(L) and τ_yz(L) (Ghosh et al., 2009).

Both of the force balance equations (Eqs. 3 and 4) were solved using the finite-element technique (Flesch et al., 2001; Ghosh et al., 2009, 2013b; Singh and Ghosh, 2019, 2020) for a 100 km thick lithosphere of varying strength (Fig. S1a in the Supplement). The laterally varying viscosities for the lithosphere were assigned from Singh and Ghosh (2020). After solving these equations, we obtained the horizontal deviatoric stresses $S_{{\mathrm{H}}_{\max }}$, strain rates, and plate velocities and compared them with observations.

The quantitative comparison between predicted and observed $S_{{\mathrm{H}}_{\max }}$ axes (Fig. 3a) was performed by computing the misfit given by sin θ(1+R) (Ghosh et al., 2013a; Singh and Ghosh, 2019, 2020), where R represents the quantitative difference between stress regimes of observed and predicted $S_{{\mathrm{H}}_{\max }}$, while θ denotes the angular difference between both. Hence, this misfit accounts for both the angular and regime misfits with values lying between 0 and 3.

The correlation between predicted deviatoric stresses and the Global Strain Rate Model (GSRM) strain rates (Fig. 3b) (Flesch et al., 2007; Ghosh et al., 2013b; Singh and Ghosh, 2019, 2020) is given by the following equation:

$$\begin{array}{}\text{(5)}& -\mathrm{1}\le \sum _{\mathrm{areas}}(\mathit{ϵ}\cdot \mathit{\tau })\mathrm{\Delta }S/\left(\sqrt{\sum _{\mathrm{areas}}\left(E^{\mathrm{2}}\right)\mathrm{\Delta }S}\cdot \sqrt{\sum _{\mathrm{areas}}\left(T^{\mathrm{2}}\right)\mathrm{\Delta }S}\right)\le \mathrm{1},\end{array}$$

where

$$\begin{array}{ll}{\displaystyle }E& {\displaystyle }=\sqrt{{\dot{\mathit{ϵ}}}_{\mathit{\varphi }\mathit{\varphi }}^{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{\mathit{\theta }\mathit{\theta }}^{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{rr}^{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{\mathit{\varphi }\mathit{\theta }}^{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{\mathit{\theta }\mathit{\varphi }}^{\mathrm{2}}}\\ {\displaystyle }& {\displaystyle }=\sqrt{\mathrm{2}{\dot{\mathit{ϵ}}}_{\mathit{\varphi }\mathit{\varphi }}^{\mathrm{2}}+\mathrm{2}{\dot{\mathit{ϵ}}}_{\mathit{\varphi }\mathit{\varphi }}{\dot{\mathit{ϵ}}}_{\mathit{\theta }\mathit{\theta }}+\mathrm{2}{\dot{\mathit{ϵ}}}_{\mathit{\theta }\mathit{\theta }}^{\mathrm{2}}+\mathrm{2}{\dot{\mathit{ϵ}}}_{\mathit{\varphi }\mathit{\theta }}^{\mathrm{2}}},\end{array}$$

$$\begin{array}{ll}{\displaystyle }T& {\displaystyle }=\sqrt{{\mathit{\tau }}_{\mathit{\varphi }\mathit{\varphi }}^{\mathrm{2}}+{\mathit{\tau }}_{\mathit{\theta }\mathit{\theta }}^{\mathrm{2}}+{\mathit{\tau }}_{rr}^{\mathrm{2}}+{\mathit{\tau }}_{\mathit{\varphi }\mathit{\theta }}^{\mathrm{2}}+{\mathit{\tau }}_{\mathit{\theta }\mathit{\varphi }}^{\mathrm{2}}}\\ {\displaystyle }& {\displaystyle }=\sqrt{\mathrm{2}{\mathit{\tau }}_{\mathit{\varphi }\mathit{\varphi }}^{\mathrm{2}}+\mathrm{2}{\mathit{\tau }}_{\mathit{\varphi }\mathit{\varphi }}{\mathit{\tau }}_{\mathit{\theta }\mathit{\theta }}+\mathrm{2}{\mathit{\tau }}_{\mathit{\theta }\mathit{\theta }}^{\mathrm{2}}+\mathrm{2}{\mathit{\tau }}_{\mathit{\varphi }\mathit{\theta }}^{\mathrm{2}}},\end{array}$$

and

$${\displaystyle }\mathit{ϵ}\cdot \mathit{\tau }=\mathrm{2}{\dot{\mathit{ϵ}}}_{\mathit{\varphi }\mathit{\varphi }}{\mathit{\tau }}_{\mathit{\varphi }\mathit{\varphi }}+{\dot{\mathit{ϵ}}}_{\mathit{\varphi }\mathit{\varphi }}{\mathit{\tau }}_{\mathit{\theta }\mathit{\theta }}+{\dot{\mathit{ϵ}}}_{\mathit{\theta }\mathit{\theta }}{\mathit{\tau }}_{\mathit{\varphi }\mathit{\varphi }}+\mathrm{2}{\dot{\mathit{ϵ}}}_{\mathit{\theta }\mathit{\theta }}{\mathit{\tau }}_{\mathit{\theta }\mathit{\theta }}+\mathrm{2}{\dot{\mathit{ϵ}}}_{\mathit{\varphi }\mathit{\theta }}{\mathit{\tau }}_{\mathit{\varphi }\mathit{\theta }}.$$

In the above equations, the second invariants of the strain rate and stress tensors are denoted by E and T. GSRM strain rates, area, and predicted deviatoric stresses are represented by ${\dot{\mathit{ϵ}}}_{ij}$, ΔS, and τ_ij, respectively. We also get the relative plate velocities and strain rates as output from models. However, to calculate the absolute plate velocities and strain rates, we require absolute viscosity values. We compute the scaling factor for relative viscosities by placing the predicted velocities in a no-net-rotation (NNR) frame, such that $\int (v\times r)dS=\mathrm{0}$, minimizing the misfit between the predicted dynamic velocities and those from Kreemer et al. (2014). Here, v denotes the horizontal surface velocity at position r, and S is the area over the Earth's surface (see Ghosh et al., 2013b, for details).

3.2 Crustal models

On the right-hand side of Eqs. (3) and (4), the first term represents the vertically integrated vertical stress. It is computed and integrated from the top of the variable topography up to depth L (100 km) (England and Molnar, 1997; Flesch et al., 2001; Ghosh et al., 2013b; Singh and Ghosh, 2019, 2020) using the following relation:

$$\begin{array}{}\text{(6)}& {\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{zz}=-\underset{-h}{\overset{L}{\int }}\left[\underset{-h}{\overset{z}{\int }}\mathit{\rho }\left(z^{\prime }\right)g\mathrm{d}{z}^{\prime }\right]\mathrm{d}z=-\underset{-h}{\overset{L}{\int }}(L-z)\mathit{\rho }\left(z\right)g\mathrm{d}z,\end{array}$$

where ρ(z), L, and h denote density, the depth to the lithosphere base (100 km), and topographic elevation, respectively. z and z^′ are variables of integration, and g represents the acceleration due to gravity. We also calculated the stresses for a thicker lithosphere (L=150 and L=200 km) as studies have shown a much thicker lithosphere in the region (Robert et al., 2017; Tunini et al., 2017) (Fig. S2a–f).

The right-hand side of Eq. (6) is given by the negative GPE per unit area. To calculate GPE and the stresses associated with it, we used three global crustal models, namely CRUST2.0 (Bassin et al., 2000), CRUST1.0 (Laske et al., 2013), and LITHO1.0 (Pasyanos et al., 2014). The upper-crust thickness lies within 15–20 km in the Zagros–Iran region for the CRUST2 model (Fig. 2a). However, the thickness of the upper crust in the Zagros–Iranian region is much higher for CRUST1 and LITHO1 (>25 km) (Fig. 2b and c). The Zagros–Iran region has a thicker middle crust (>20 km) in the case of both the CRUST2 and LITHO1 models (Fig. 2d and f), while CRUST1 shows a much thinner middle crust (<12 km) in this region (Fig. 2e). The lower crust in the Zagros–Iran region is found to be very thin (<10 km) for all three models (Fig. 2g–i).

Figure 2 Thickness and density variations of different layers in all three crustal and lithospheric models: CRUST2 (a, d, j, m, p), CRUST1 (b, e, h, k, n, q), and LITHO1 (c, f, i, j, o, r).

Figure 3 (a) Most compressive horizontal principal axes ($S_{{\mathrm{H}}_{\max }}$) from WSM (Heidbach et al., 2016). Red indicates normal fault regime, and blue indicates thrust regime, whereas green denotes strike–slip regime. (b) Observed plate velocities in a no-net-rotation frame of reference from Kreemer et al. (2014) plotted on top of the second invariant of strain rate tensors obtained from Kreemer et al. (2014) plotted on $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid.

The density variations in the study area are minimal for the CRUST2 model. CRUST2 also shows the highest average density in all three layers (>2.7 g cm⁻³) (Fig. 2j, m, and p). CRUST1 also indicates an average density of ∼2.72 g cm⁻³ in the Zagros–Iran region for the upper crust (Fig. 2k). The middle and lower crustal layers of CRUST1 show average densities of 2.80 and ∼2.85 g cm⁻³, respectively (Fig. 2n and q). The LITHO1 model shows the lowest average density in the study area for all three layers (Fig. 2l, o, and r). The upper crust of the LITHO model shows an average density of ∼2.65 g cm⁻³. The central Iran block has a relatively denser upper crust (∼2.75 g cm⁻³), while the density decreases to ∼2.62–2.64 g cm⁻³ near the Zagros region. Similar patterns of density variations are observed in the middle and lower crusts of the LITHO1 model (Fig. 2o and r). Such differences in thickness and density data lead to varying GPE values and hence different stresses.

3.3 Mantle convection

We ran mantle convection models using the HC solver code (Hager and O'Connell, 1981). HC is a semi-analytical mantle convection code that uses density anomalies derived from seismic tomography models and radial viscosity as inputs. Here, we considered four global mantle convection models, namely S40RTS (Ritsema et al., 2011), SAW642AN (Mégnin and Romanowicz, 2000), 3D2018_S40RTS, and S2.9_S362 to infer the mantle density anomalies. 3D2018_S40RTS is a merged model of SV wave upper-mantle tomography model, 3D2018_Sv given by Debayle et al. (2016), and S40RTS. S2.9 is a global tomography model of the upper mantle with higher resolution, which is given by Kustowski et al. (2008b). We merged this model with the global shear wave velocity model, S362ANI (Kustowski et al., 2008a), to obtain the merged tomography model of S2.9_S362. We used two different radial viscosity structures, namely GHW13, which is the best viscosity model from Ghosh et al. (2013b), and SH08, given by Steinberger and Holme (2008). GHW13 is a four-layered viscosity structure, with a highly viscous lithosphere (∼10²³ Pa s⁻¹). The viscosity drops to ∼10²⁰ Pa s⁻¹ in the asthenosphere, which again increases to ∼10²¹ Pa s⁻¹ in the upper mantle and ∼10²² Pa s⁻¹ in the lower mantle (Fig. S1b). On the other hand, the viscosity in the SH08 model increases gradually with depth, and it has a slightly weaker lithosphere as compared to GHW13. It has the highest viscosity value of 10²³ Pa s⁻¹ around 2000–2300 km depth and significantly lower viscosity for the D” layer (Fig. S2b). The GHW13 viscosity model performed slightly better than SH08 in fitting the observed parameter; thus, we have shown results from the same throughout this paper. However, we have also included the predicted results and their fit to the observables in the Supplement (Table S1).

3.4 Data

To have better constraints on this study's models, we also estimated $S_{{\mathrm{H}}_{\max }}$ (most compressive horizontal principal axes) orientations, as well as plate velocities. Various deformation indicators such as $S_{{\mathrm{H}}_{\max }}$ orientations from the World Stress Map (WSM) (Heidbach et al., 2016), strain rates, and plate velocities from the Global Strain Rate Model (Kreemer et al., 2014) were used to perform a quantitative comparison with the predicted results of this study (Fig. 3).

The WSM project provides a global database of the crustal stress field obtained from various sources such as focal mechanisms, geophysical logs of borehole breakouts and drilled-induced fractures, engineering methods such as hydraulic fractures and overcoring, and geological indicators that are obtained from fault slip analysis and volcanic alignments. These data have been assigned quality ranks from A to E based on the accuracy range. A-type data suggest that the standard deviations of $S_{{\mathrm{H}}_{\max }}$ orientations are within ±15^∘ range; this range is ±20^∘ for B-type data, ±25^∘ for C-type data, and ±40^∘ for D-type data. However, E-type data indicate that the data records are either incomplete or from non-reliable sources or that the accuracy is $>\pm \mathrm{40}$^∘. This study uses A–C quality stress data records (Fig. 3a). Observed $S_{{\mathrm{H}}_{\max }}$ axes are aligned in a NNE–SSW direction in Zagros, with dominant thrust faulting. NW and central Iran show some strike–slip mode of deformation with NE–SW compressional directions.

The strain rates and plate velocities are taken from the GSRM v2.1 model (Kreemer et al., 2014) (Fig. 3b). GSRM v2.1 provides a global dataset of strain rates and plate motions that are determined using ∼22 500 geodetic plate velocities. Higher strain rates are observed along the simply folded mountains (∼40–$\mathrm{100}\times {\mathrm{10}}^{-\mathrm{9}}$ yr). Most of Iran shows strain rates between 4–$\mathrm{10}\times {\mathrm{10}}^{-\mathrm{9}}$ yr. The plate motions used in our study for comparison with predicted velocities are given in a no-net-rotation (NNR) frame interpolated on a $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid. The velocity vectors show an eastward motion in the study area, which becomes nearly E–W in the Afghan block (Fig. 3b).

Figure 4 (a, c, e) Deviatoric stresses predicted from GPE variations, plotted on top of their second invariants. The compressional stresses are denoted by solid black arrows, while white arrows show tensional stresses. $S_{{\mathrm{H}}_{\max }}$ axes predicted from GPE variations are plotted in (b), (d), and (f). Red denotes the tensional regime, blue is for thrust, and green is for the strike–slip regime.

4 Results

4.1 Stress and deformation due to GPE

Three crustal models (CRUST1.0, CRUST2.0, and LITHO1.0) were used to compute GPE within the study region. The second invariant of stress computed using GPE lies within ∼10–12 MPa along Zagros for the CRUST2 and CRUST1 models (Fig. 4a and c). The LITHO1 model predicts larger stress magnitudes along Zagros (Fig. 4e). NE–SW compressional stresses are observed along the frontal faults of Zagros (MFF) (Fig. 1a and c). The central part of the Zagros thrust faults (MZT) shows the strike–slip mode of faulting for nearly all three models (Fig. 4b, d and f). The strike–slip regime further extends into the Sanandaj–Sirjan Zone (SSZ), while it lies north of the MZT for the CRUST2 and LITHO1 models (Fig. 4b and f). It transitions to a thrust type of deformation in the north of the MZT for CRUST1 (Fig. 4d). The Urmia–Dokhtar magmatic arc (UDMA) and central Iran also show the strike–slip mode of faulting for CRUST2 and LITHO1. The north of the MRF shows tension for the CRUST2 model, while CRUST1 predicts this area to be predominantly strike–slip. On the other hand, the entire region shows significant compression for LITHO1 model.

We compared predicted $S_{{\mathrm{H}}_{\max }}$ from our three GPE-only models to observed $S_{{\mathrm{H}}_{\max }}$ orientations and type obtained from the WSM (Heidbach et al., 2016) by computing the regime misfit (Fig. 5, left panel). The average misfit is lowest for the LITHO1 model with a value of 0.59 (Fig. 5g), while the CRUST2 model shows the highest average misfit of 0.77 (Fig. 5a). High misfits (2–3) are observed to the north of the MRF and Tehran for CRUST2, while the lowest (<1) are observed in the case of LITHO1, suggesting that the dominant mode of faulting in this area is possibly thrust as opposed to normal deformation predicted by CRUST2. In central Iran, $S_{{\mathrm{H}}_{\max }}$ misfit is low (<1) when the dominant mode of deformation is strike–slip, as predicted by the LITHO1 model.

Figure 5 (a, d, g) Total misfit between observed and predicted S_H−max from GPE variations. Correlation coefficients between strain rate tensors obtained from Kreemer et al. (2014) and deviatoric stresses predicted using GPE variations are shown in (b, e, h), with average regional correlation coefficients in each figure's bottom right. (c, f, i) Observed velocities (black) and predicted plate velocities(white) from GPE variations in the NNR frame, plotted on the top of the angular misfit between both.

On calculating the correlation between the predicted deviatoric stresses and GSRM strain rates, the LITHO1 model shows the highest average correlation (0.92) (Fig. 5, middle panel). The correlation is found to be extremely poor ($\sim -\mathrm{1}$) for the CRUST2 model in the north of the MRF (Fig. 5b). Such poor correlation suggests that the predicted stresses differ entirely from those causing deformation. For example, anti-correlation north of the MRF suggests that the dominant mode of deformation in this area might be thrust rather than normal faulting. Again, the correlation coefficient is less than 0.2 in the central Iranian block for the CRUST2 and CRUST1 models (Fig. 5b and e), while the LITHO1 model shows a better correlation, suggesting that the strike–slip type of deformation is more prominent in central Iran (Fig. 5h).

We predicted the plate velocities for all three models in the NNR frame and compared them with observed plate velocities obtained from Kreemer et al. (2014) (Fig. 5, right panel). CRUST2 gives the lowest rms error (7.32 mm yr⁻¹) and the lowest angular misfit (5.5^∘) (Fig. 5c). The LITHO1 model shows high misfits (>20^∘) between observed and predicted velocities in the east of central Iran (i.e., Afghan block) (Fig. 5i). Both the CRUST2 and LITHO1 models predict the plate velocities very closely to observed ones in the Zagros mountains, as shown by nearly zero angular misfits along Zagros (Fig. 5c and i). CRUST1 performs at an average level in predicting the plate velocities in the study area (Fig. 5f).

Interestingly, the use of a thicker lithosphere to calculate GPE leads to the introduction of more compressional stresses in the region (Fig. S2a–f). The average misfit between predicted and observed $S_{{\mathrm{H}}_{\max }}$ is found to be lowest for the 200 km thick lithosphere (Table S2). Similarly, the correlations between strain rate tensor and predicted stresses and the rms errors between observed and predicted NNR velocities show significant improvement for a thicker lithosphere. However, the improvement in the fit is better for CRUST2 as opposed to the other two models, CRUST1 and LITHO1, where the misfit between observed and predicted velocities shows an increase. Thus, we can say that, while considering lithospheric contributions only, the thicker lithosphere does a better job of explaining the observed deformation indicators (Table S2).

Figure 6 (a, c, e, g) Deviatoric stresses predicted using mantle tractions derived from various tomography models for GHW13 viscosity structure, plotted on the second invariant of deviatoric stresses. The white arrows denote tensional stresses, and the black arrows indicate compressional stresses. $S_{{\mathrm{H}}_{\max }}$ predicted from mantle tractions is shown in (b), (d), (f), and (h). Red denotes the tensional regime, blue is for thrust, and green is for the strike–slip regime.

4.2 Stress and deformation due to mantle convection

The deviatoric stresses predicted using all four mantle convection models are found to be mostly compressional along MFF (Fig. 6). All models, except for SAW642AN, predict the strike–slip mode of faulting in the northwestern parts of Zagros with nearly E–W-oriented extensional axes and N–S compressional axes (Fig. 6a, e and g). On the other hand, SAW642AN shows predominant compression within this area (Fig. 6e). S40RTS, 3D2018_Sv, and S2.9_S362 show strike–slip deformation in the northwestern parts of SSZ, UDMA, and northwestern Iran. Central Iran is predicted to have mostly compressional stresses by all models, except for S40RTS. Thrust type of deformation is predicted in the Afghan block by all models with some intermittent strike–slip deformation. The SINGH_SAW model predicts the whole Afghan block in the strike–slip regime (Fig. 6g and h). S40RTS and S2.9_S362 predict higher stress magnitude in northwestern parts of the Zagros orogeny system and central Iran compared to other models.

Figure 7 Parameters predicted from mantle tractions and their comparisons with observables. (a, d, g, j) Total misfit between $S_{{\mathrm{H}}_{\max }}$ obtained from WSM (Heidbach et al., 2016) and those predicted using mantle tractions derived from various tomography models using GHW13 viscosity structure. Correlation coefficients between strain rate tensors obtained from Kreemer et al. (2014) and deviatoric stresses predicted using basal tractions are shown in (b), (e), (h), and (k) with average regional correlation coefficients in each figure's bottom right. (c, f, i, l) Observed velocities (black) and plate velocities predicted using mantle tractions (white) in the NNR frame plotted on the top of the angular deviation between both.

The misfit between observed and predicted $S_{{\mathrm{H}}_{\max }}$ is found to be much lower for mantle convection models (0.54–0.57) (Fig. 7, left panel) than for those of GPE-only models (Fig. 5, left panel), evidently showing the importance of mantle flow. The lowest average misfit is observed for SAW642AN (0.54) (Fig. 7d), though the misfit increases in the east, the Lut block, and near the MSZ. The correlation of predicted deviatoric stresses with GSRM strain rates improves over GPE-only models (Fig. 7, middle panel), with SAW642AN yielding the highest correlation coefficient (0.91) (Fig. 7e). Correlation drops below 0.4 in parts of central Iran. S40RTS predicts the plate velocities most closely to the observed one out of all models, with the lowest rms error (∼6.20 mm yr⁻¹) between predicted and observed plate velocities (Fig. 7c). On the other hand, SAW642AN and 3D2018_S40RTS models show high misfits (rms error ∼10 mm yr⁻¹) as they are unable to match observed plate velocities in the Zagros–Iran plateau, both in orientation and magnitude (Fig. 7f and i).

As discussed above, mantle convection models perform better in predicting deviatoric stresses in the study area, which is made evident by the high correlation between predicted stresses and observed strain rates and by the low misfits between observed and predicted $S_{{\mathrm{H}}_{\max }}$. However, the error in predicting plate velocities is higher for mantle convection models than for GPE-only models. As there are still significant misfits in fitting the observables, we added the deviatoric stresses predicted from GPE differences and mantle convection models to constrain the total stress field in the Zagros–Iranian plateau, which may account for both forces.

We also ran the S40RTS model with LAB (lithosphere–asthenosphere boundary) at 150 and 200 km (Fig. S2g and h). Similarly to GPE models, the fit to observed data shows an improvement when LAB is at 200 km, though the stress patterns do not change significantly (Table S2).

Figure 8 (a, c, e, g) Deviatoric stresses predicted using the combined effects of GPE computed from CRUST2 and mantle tractions derived from various tomography models plotted on top of their second invariants. The white arrows denote tensional stresses, and the black arrows indicate compressional stresses. Panels (b, d, f, h) show $S_{{\mathrm{H}}_{\max }}$ predicted from these models. The red lines denote the tensional regime, blue is for thrust, and green is for the strike–slip regime.

Figure 9 (a, c, e, g) Deviatoric stresses predicted using the combined effects of GPE computed from CRUST1 and mantle tractions derived from various tomography models plotted on top of their second invariants. The white arrows denote tensional stresses, and black arrows indicate compressional stresses. Panels (b, d, f, h) show $S_{{\mathrm{H}}_{\max }}$ predicted from these models. The red lines denote tensional regime, blue is for thrust, and green is for strike–slip regime.

4.3 Stress and deformation by GPE and mantle convection

Adding mantle contributions to GPE-only models led to significant changes in total deviatoric stresses for all models (Figs. 8–10). There is a significant increase in the total stress magnitude of the entire study area, except north of the MRF and SE of central Iran, which show slightly lower stresses (<16 MPa) for combined models of CRUST2 and mantle convection (Fig. 8). These models show predominant compression in most of Zagros, SSZ, UDMA, and NW and central Iran, except for the strike–slip type of deformation in the northwestern parts. The joint models of CRUST1 and mantle convection predict higher stresses (>25 MPa) in NW Iran and at MFF (Fig. 9). Interestingly, the stresses drop below 20 MPa towards the north of the HZF and the MRF till the southern Caspian. The combined models of CRUST1 and mantle convection show that compressional stresses are dominant in the study area, with occasional strike–slip faulting in the northwest (Fig. 9, right panel). The stresses predicted by combined models of LITHO1 and mantle convection models are higher in magnitude than other models in the study area (>25 MPa) (Fig. 10). S40RTS+LITHO and S2.9_S362+LITHO models show high stresses in Zagros (>50 MPa) (Fig. 10a and g).

Figure 10 (a, c, e, g) Deviatoric stresses predicted using the combined effects of GPE computed from LITHO and mantle tractions derived from various tomography models plotted on top of their second invariants. The white arrows denote tensional stresses, and the black arrows indicate compressional stresses. Panels (b, d, f, h) show $S_{{\mathrm{H}}_{\max }}$ predicted from these models. The red lines denote tensional regime, blue is for thrust, and green is for the strike–slip regime.

The combined models show a lower misfit between observed and predicted $S_{{\mathrm{H}}_{\max }}$ (Fig. 11), especially when compared to GPE-only models (Fig. 5, left panel). SAW642AN+LITHO showed the lowest average misfit of 0.47 (Fig. 11f). Interestingly, SAW642ANcr2 and 3D2018_S40RTScr2 show low misfits in the Zagros–Iranian plateau region despite not having the lowest average misfit (Fig. 11d and g). The higher misfits in NW Iran and in the SE of the central Iran block observed for GPE-only models get reduced significantly due to the addition of mantle-derived stresses, reflecting the importance of mantle convection in these areas.

Figure 11 Total misfit between observed $S_{{\mathrm{H}}_{\max }}$ from WSM (Heidbach et al., 2016) and $S_{{\mathrm{H}}_{\max }}$ predicted using the combined effects of GPE computed from different crustal models and mantle tractions derived from various tomography models.

As we look at the correlation between predicted stress tensors and GSRM strain rate tensors, the overall correlation is better for combined models (Fig. 12), especially for combined models of LITHO1 and mantle convection (Fig. 12, right panel). A high average correlation coefficient of 0.94 is observed for SAW642AN+LITHO, 3D2018_S40RTS+LITHO, and S2.9_S362+LITHO1 (Fig. 12f, i, and l). Despite an overall improvement in correlation between observed strain tensors and predicted deviatoric stresses, the correlation is found to be much poorer in areas such as the northwestern parts of Zagros and east of the central Iranian block for combined models of mantle convection and the GPE-only models of CRUST2 and CRUST1 (Fig. 12, left and middle panels). In NW Zagros, mantle-only models are found to perform much better as they show better correlation (Fig. 7, middle panel), thus suggesting that mantle-derived stresses should be much higher than those from GPE to explain the observed deformation in these areas.

Figure 12 Correlation coefficients between strain rate tensors from Kreemer et al. (2014) and deviatoric stress tensors predicted using the combined effects of GPE computed from different crustal models and mantle tractions derived from various tomography models. The figure shows the average correlation coefficient in the lower-right corner of the figure.

Again the combined models of GPE and mantle tractions give lower rms errors when predicted plate velocities are compared to the observed ones. S40RTScr2 shows the lowest rms error (3.28 mm yr⁻¹) and the least average angular misfit (3.0^∘) between predicted and observed plate velocities (Fig. 13a). Relatively, the combined models of S40RTS and S2.9_S362 and GPE perform much better than other models in predicting the orientation and magnitude of plate velocities. Significant misfits are observed for the SAW642ANcr1 and 3D2018_S40RTScr1 models. The joint models of S40RTS and GPE for a thicker lithosphere do not offer any significant changes in terms of stresses and their fit to observed data (Table S2) (Ghosh et al., 2009; Jay et al., 2018; Hirschberg et al., 2018). Thus, considering the lithosphere base at 100 km appears to be a satisfactory approach.

Figure 13 Plate velocities predicted using the combined effects of GPE computed from different crustal models and mantle tractions derived from various tomography models plotted on top of angular misfit (θ). Black arrows represent observed NNR velocities (Kreemer et al., 2014), and white ones denote predicted velocities.

5 Discussion

The Zagros–Iranian plateau region is formed due to the convergence of the Arabian plate towards the Eurasian plate. The Zagros mountain belt demarcates the southwestern boundary of the deformation zone, whereas it is bounded by the Makran subduction zone in the southeast and by the Afghan block in the east. Kopet Dagh and Arborz act as this region's northeastern and northern boundaries (Irandoust et al., 2022). We modeled the stresses and deformation parameters in the study area by solving the force balance equation using the finite-element method for a global grid of $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ resolution considering two primary sources of stresses, namely GPE and mantle tractions. GPE was calculated using the thickness and density variation from the different global models CRUST1.0, CRUST2.0, and LITHO1.0. The shear tractions were computed from a density-derived mantle convection model.

The magnitude of stresses due to GPE variations was below 15 MPa in the Iranian plateau for the CRUST2 and CRUST1 models (Fig. 4a and c). However, the LITHO1 model predicted higher stresses (>30 MPa) with predominant compression in parts of the Zagros–Iran region and Afghan block. Most of the convergence of the Arabian and Eurasian plates has been accommodated by shortening across Zagros (Irandoust et al., 2022; Khodaverdian et al., 2015). Walpersdorf et al. (2006) and Hessami et al. (2006) suggested a nearly pure N–S shortening of 8±2 mm yr⁻¹ in southeastern Zagros. The convergence occurs perpendicularly to the simply folded mountains and is restricted to the shore of the Persian Gulf. Earthquake focal mechanisms also show reverse faulting within this area (Berberian, 1995; Hatzfeld et al., 2010; Hatzfeld and Molnar, 2010; Irandoust et al., 2022). In our study, the LITHO1 model predicted the thrust mode of faulting within Zagros, which is consistent with these results. In NW Zagros, Hatzfeld et al. (2010); Hatzfeld and Molnar (2010), Jackson and McKenzie (1984), Khorrami et al. (2019), Talebian and Jackson (2002), and various others have suggested partitioning of deformation. The oblique shortening is partitioned into strike–slip faulting that is accommodated by the MRF, while shortening occurs perpendicularly to the mountain belt (Hatzfeld et al., 2010; Hatzfeld and Molnar, 2010; Jackson and McKenzie, 1984; Khorrami et al., 2019; Talebian and Jackson, 2002). In considering lithospheric models only, we predicted the normal mode of faulting to be dominant in this area for CRUST2. On the other hand, the CRUST1 model predicted strike–slip components in the northern segment of the MRF, while LITHO1 showed thrust type of deformation in this area. Interestingly, the misfits of predicted parameters with various observations of $S_{{\mathrm{H}}_{\max }}$, strain rates, and plate velocities were found to be lowest for the LITHO1 model, thus arguing for thrust type of deformation in this area. SSZ in the north of the MZT consists of various thrust systems (Alavi, 1994). CRUST1 predicted the thrust mode of faulting in this region, while the CRUST2 and LITHO1 models showed an intermittent strike–slip type of faulting. Alborz and Kopet Dagh in the north have also been subjected to reverse faulting (Allen et al., 2003; Hatzfeld and Molnar, 2010; Hollingsworth et al., 2010; Irandoust et al., 2022; Khodaverdian et al., 2015), which has also been shown by the CRUST1 and LITHO1 models. Models predicting thrust in the Talesh mountains showed low misfits to observations, suggesting thrusting of the mountain range over the basin, with slip vectors directed towards the southern Caspian Sea (Irandoust et al., 2022). The N–S convergence in the Kopet Dagh range was predicted by the LITHO1 model considering the contribution from lithospheric density and topographic variations only. The shearing between central Iran and the Afghan block caused due to varying rates of shortening across Zagros, Alborz, and Caucasus is accommodated by strike–slip faults near the Lut block boundaries (Khorrami et al., 2019; Vernant et al., 2004; Walpersdorf et al., 2014). Again, the LITHO1 model predicted similar strike–slip deformation in these areas; however, CRUST2 and CRUST1 failed to do so.

Table 1 Summary of quantitative comparison of predicted results of various models with observed data.

The stresses predicted using basal tractions were mostly compressional in southeastern Zagros owing to the convergence of Arabia–Eurasia (Fig. 6). However, all models except SAW642AN predicted strike–slip type of deformation in northwestern Zagros (MRF), which concurs with the results from various studies (Hatzfeld et al., 2010; Hatzfeld and Molnar, 2010; Jackson and McKenzie, 1984; Khorrami et al., 2019; Talebian and Jackson, 2002). The mantle-derived stress parameters showed a better fit to observables than those from GPE variations (Fig. 7, left and middle panels), though the correlation dropped below 0.5 in central Iran. Here, mantle convection models found a compressional type of deformation, while Baniadam et al. (2019) and Khorrami et al. (2019) suggested strike–slip faulting along the fault system bounding the Lut block. The velocity misfits were very high for all models except S40RTS (Fig. 7, right panel). Although we used four tomography models to compute the mantle-derived stresses, the stress regimes for all models are found to be similar, with varying magnitudes. Such results suggest that nearly all four seismic tomography models are relatively consistent in predicting the stresses in this region.

Adding the GPE-derived stresses to those from the mantle to obtain the total lithospheric stress field showed a notable improvement in constraining the observed deformation parameters. The final stress regimes also varied significantly depending on particular combinations of GPE and mantle convection models. All joint models of CRUST2 and mantle tractions showed lower magnitudes of stresses (<15 MPa) in the north of the MRF, Tehran, and the southern Lut block. The stresses showed an obvious increase in these areas for other models. Significantly higher stresses (>30 MPa) were also observed near the collisional front (MFF) for all models. On comparison with observations, combined models of CRUST2 and mantle tractions showed significant improvement in fit, except in areas north of the MRF and Tehran. The CRUST1 model, when added to mantle contribution, predicted thrust faulting along the faults bounding the Lut block, leading to poor correlation (<0.5). On the other hand, combined LITHO1 and mantle convection models gave a much better fit in this area as they predicted strike–slip faulting. The use of different mantle convection models is much less sensitive in the Iran–Zagros region as most models can match various surface observables reasonably well.

On running various models and comparing the stresses in Zagros–Iran, we try to explain the relative roles of GPE and mantle tractions in causing observed deformation. The contributions from both sources vary significantly among different models. However, these variations arise mainly from GPE-only models, which may be due to uncertainties in the crustal models of this area. Another interesting observation from this study is that the role of GPE in the study region may not be that significant as mantle-derived stresses were able to explain many of the deformation indicators. To get a quantitative constraint on the best model, we computed a total error as given below:

$$\begin{array}{}\text{(7)}& \mathrm{Total}\mathrm{error}=S_{{\mathrm{H}}_{\max }}\mathrm{error}+\mathrm{1}-C_{\mathrm{strain}}+V_{\mathrm{rms}}.\end{array}$$

The $S_{{\mathrm{H}}_{\max }}$ error in the above equation is calculated as described in Sect. 3.4, while C_strain is the correlation computed using Eq. (6). V_rms is the rms error between predicted and observed velocities. The total errors calculated using Eq. (7) have been tabulated in Table 1. S40RTScr2 is found to have the lowest error.

Figure 14 Predicted parameters of the best-fit model, S40RTScr2. (a) GPS (blue) and predicted (red) plate velocities with respect to a fixed Eurasian plate, (b) fast polarization directions (FPDs) (blue), and $S_{{\mathrm{H}}_{\max }}$ (red) are plotted for the best-fit model. (c) Correlation between deviatoric stresses predicted from GPE and mantle convection models and (d) ratio ($T_{\mathrm{1}}/T_{\mathrm{2}}$) of the second invariant of deviatoric stresses from GPE (T₁) to those from mantle tractions (T₂).

We also calculated plate velocities with respect to the Eurasian plate (Fig. 14) and compared them with observed GPS velocities relative to Eurasia. The GPS velocities were obtained from various studies conducted in this area (ArRajehi et al., 2010; Bayer et al., 2006; Frohling and Szeliga, 2016; Khorrami et al., 2019; Masson et al., 2006, 2007; Raeesi et al., 2017; Reilinger and McClusky, 2011; Vernant et al., 2004). GPS measurements show a northward convergence rate of ∼22 mm yr⁻¹ for Arabia relative to Eurasia (Reilinger et al., 2006; Vernant et al., 2004); however, it varies significantly along Zagros. Southeastern Zagros shows the highest convergence rates of ∼25 mm yr⁻¹, oriented in the north–northeast direction. GPS vectors are oriented northward in central Zagros, which transitions to become oriented north–northwest in northwestern parts of Zagros with the lowest convergence rates of ∼18 mm yr⁻¹ (Hatzfeld and Molnar, 2010; Hatzfeld et al., 2010; Khorrami et al., 2019). Vernant et al. (2004) suggested that the MSZ accommodates most of the shortening (19.5±2 mm yr⁻¹) east of 58^∘ E, while fold-and-thrust belts of Zagros, Alborz, and Caucasus collectively accommodate the shortening west of 58^∘ E. GPS velocities in the east of Iran (Afghan block) are very small in magnitude. To the west, velocities increase, showing a westward rotation of Anatolia (Khorrami et al., 2019; Reilinger et al., 2006). The northern part of Iran shows that GPS vectors are aligned towards the northeast. We found that the combined model of S40RTS and CRUST2 can approximately match the GPS velocities (Fig. 14a). Predicted plate velocities with respect to the fixed Eurasian plate show a northward movement of ∼2–3 cm yr⁻¹ in southeastern Zagros. The plate moves in an NNE direction east of central Zagros (53^∘ E). On the other hand, west of 53^∘ E shows a movement in an NNW direction, becoming much more prominent in the north. However, the convergence rates in the east of Iran, i.e., the Lut block and Afghan block, are predicted to be much higher (∼1–2 cm yr⁻¹) than those suggested by various observations. Plate velocities predicted by joint models, S40RTScr1 and S40RTS+LITHO1, show nearly N–S contraction of very high magnitudes (4–5 cm yr⁻¹) throughout the region (Fig. S3), which suggests much higher rates of deformation than those suggested by the above-mentioned studies.

We also used shear-wave-splitting measurements to further study the deformation in the Zagros–Iran region by comparing them with $S_{{\mathrm{H}}_{\max }}$ (Fig. 14b). The fast polarization directions (FPDs) are the indicators of seismic anisotropy. We consider two primary causes of seismic anisotropy, namely induction by stress and that due to the structure of the region (Yang et al., 2018). If the FPDs are parallel to $S_{{\mathrm{H}}_{\max }}$ orientations, it suggests that anisotropy is associated with stress. On the other hand, the latter kind of anisotropy is related with the alignment of fault, fast axes of minerals that may cause polarization, and sedimentary bedding planes. The FPDs in our study were obtained from Sadeghi-Bagherabadi et al. (2018) and Kaviani et al. (2009, 2021). The FPDs are subparallel to $S_{{\mathrm{H}}_{\max }}$ orientations in NW Zagros, the Arabian plate, northern Iran, and the MSZ. Such a correlation between both indicates that anisotropy in this region may be stress induced. Additionally, the correlation of $S_{{\mathrm{H}}_{\max }}$ orientations and FPDs argues for a good coupling between the lithosphere and mantle in those areas. In contrast, Sadeghi-Bagherabadi et al. (2018) showed FPDs parallel to the strike of the fault (subparallel to $S_{{\mathrm{H}}_{\max }}$ directions of CRUST2), suggesting that seismic anisotropy mainly reflects the deformation in the lithospheric mantle. Again, FPDs are subparallel to the strike of range in northeastern Iran, eastern Kopet Dagh, and central Alborz, indicating structure-induced anisotropy caused by strong shearing along the strike–slip faults (Gao et al., 2022; Kaviani et al., 2021).

To explore the relative roles of lithospheric and mantle-derived stresses, we compared the deviatoric stresses from CRUST2 to those from S40RTS. We performed a correlation between both stresses by using Eq. (5) and found a high correlation (>0.5) near the MSZ and central Zagros (Fig. 14c). The correlation degrades north of the simply folded mountains and NW Iran. The stresses are anti-correlated in northwestern parts of higher Zagros and north of the MRF and Tehran as CRUST2 predicted NNE–SSW tension (Fig. 4b) as opposed to the strike–slip faulting predicted by S40RTS (Fig. 6b). The Lut block also shows a slight anticorrelation between stresses ($\sim -\mathrm{0.5}$) as the stresses predicted by CRUST2 are very low. The log of the ratio of second invariants of deviatoric stresses from GPE variations (T₁) to those of mantle tractions (T₂) is plotted in Fig. 14d. Positive values of the logarithmic ratio suggest the dominance of GPE-derived stresses over mantle ones, as observed in the south of the collisional boundary (MFF). The ratio is negative in most parts of the Iranian plateau and Zagros, indicating that the magnitude of mantle-derived stresses is higher than that of GPE, especially in higher Zagros and central Iran (Fig. 14d).

The deformation in the Zagros–Iran plateau region has been found to exhibit various similarities to another similar complex collision zone, i.e., the Himalaya–Tibetan plateau region, as both continental collisions went through many of the same processes. The high topography in both collisions reflects ongoing crustal deformation through crustal thickening and shortening. However, there are differences in the convergence rates, total amounts of convergence, and various stages of development of the Zagros–Iran and Himalaya–Tibet regions (Hatzfeld and Molnar, 2010). Singh and Ghosh (2020) studied the deformation in the Himalaya–Tibet region by joint modeling of lithosphere and mantle. They showed that GPE plays a crucial role in the ongoing deformation of the India–Eurasia collision zone as it leads to the observed E–W extension in the Tibetan plateau. In contrast, we found that GPE has a much lesser role in the Zagros–Iran plateau region (Fig. 14d), and no normal mode of faulting is observed in this area. In the Zagros–Iran plateau region, mantle convection appears to be the primary driver of deformation in most parts, as discussed above. Despite these differences, numerical models argue for a good coupling between the lithosphere and mantle in both collision zones, which is also supported by seismic anisotropy studies in both regions (Kaviani et al., 2021; Singh et al., 2016; Sol et al., 2007).

6 Conclusion

The Zagros–Iranian plateau region has large deformations along and across the collision zones. Therefore, we conducted numerical simulation studies for stress and deformations. The stresses predicted in this region were primarily compressional, with magnitudes lower than 30 MPa. The southeastern boundary of Zagros was found to be under high stress, which is also reflected by higher convergence rates. Mantle convection models were able to constrain most observations in the Iranian plateau. However, the misfits with observations were much larger in the east of Iran, when only mantle contributions were considered. The combined models of lithosphere and mantle-derived stresses can explain surface observables in most of the area, suggesting a good lithosphere–mantle coupling, except east of Iran. The shearing in those areas was predicted by lithosphere models, though variation in the lithospheric and density structures given by these models lead to varying degrees of misfits. Hence, there is a need for better constraint on lithospheric structure in this area.

The mantle-derived stresses were found to be much higher than lithospheric stresses; thus, the overall stress regimes predicted by combined models were more biased towards the compressional type of stresses. This caused our combined models to predict the thrust mode of faulting in most cases, especially when lithospheric stresses were computed from the CRUST1 and LITHO1 models. The CRUST2 model predicted more extensional stress in the Iranian plateau, which in turn balanced the effect of compressional stresses predicted by mantle convection models, hence leading to the prominence of the strike–slip mode of faulting in the northwestern parts of the study region. The rate of convergence of Arabia relative to a fixed Eurasia was found to vary along the Zagros orogeny in a similar way to that found with GPS measurements.