3-D Geothermal Model of the Lurestan Sector of the Zagros Thrust Belt, Iran

Abstract

The Zagros thrust belt is a large orogenic zone located along the southwest region of Iran. To obtain a better knowledge of this important mountain chain, we elaborated the first 3-D model reproducing the thermal structure of its northwestern part, i.e., the Lurestan arc. This study is based on a 3-D structural model obtained using published geological sections and available information on the depth of the Moho discontinuity. The analytical calculation procedure took into account the temperature variation due to: (1) The re-equilibrated conductive state after thrusting, (2) frictional heating, (3) heat flow density data, and (4) a series of geologically derived constraints. Both geotherms and isotherms were obtained using this analytical methodology. The results pointed out the fundamental control exerted by the main basement fault of the region, i.e., the Main Frontal Thrust (MFT), in governing the thermal structure of the crust, the main parameter being represented by the amount of basement thickening produced by thrusting. This is manifested by more densely spaced isotherms moving from the southwestern foreland toward the inner parts of orogen, as well as in a lateral variation related with an along-strike change from a moderately dipping crustal ramp of the MFT to the NW to a gently dipping crustal ramp to the SE. The complex structural architecture, largely associated with late-stage (Pliocene) thick-skinned thrusting, results in a zone of relatively high geothermal gradient in the easternmost part of the study area. Our thermal model of a large crustal volume, besides providing new insights into the geodynamic processes affecting a major salient of the Zagros thrust belt, may have important implications for seismotectonic analysis in an area recently affected by a M_w = 7.3 earthquake, as well as for geothermal/hydrocarbon exploration in the highly perspective Lurestan region.

1. Introduction

The Zagros thrust belt forms part of the Alpine-Himalayan orogen. This mountain chain formed by the convergence between the Eurasian and Arabian plates, which occurred since the Late Cretaceous [1,2,3,4,5,6,7]. The orogenic process is still active at present, as confirmed by recent geodetic calculations which estimate a northward motion of the Arabian plate of ~20 mm yr⁻¹ in the fixed-Eurasia reference frame [8]. In the Lurestan arc, active plate convergence produces high-magnitude seismicity, as testified by the recent seismic events of 12 November 2017 (M_w = 7.3) and 25 November 2018 (M_w = 6.3) [9].

During the last years, several authors studied the thermal structure of the Lurestan arc, resulting in a series of 2-D models (i.e., [10,11,12]). In this study, an improved analytical calculation procedure was implemented to obtain for the first time a crustal 3-D thermal model of the Lurestan arc that took into account the temperature variation due to a re-equilibrated conductive state associated with thrusting, frictional heating, heat flow density data, and a series of geologically derived constraints. The computational procedure was applied to a 3-D geometric model built using published geological sections [4,12,13,14] integrated with the Moho depth obtained by the work of Jiménez-Munt et al. [15]. Obtaining a comprehensive 3-D picture of the thermal structure of the Lurestan arc is of pivotal importance for any further geodynamic modelling of this seismically active region (e.g., [13]), as well as for enhanced geothermal/hydrocarbon exploration and development in this highly perspective portion of the Zagros Mountains [16].

2. Geological Background

The NW-SE striking Zagros mountain belt forms part of the Alpine-Himalayan orogen (Figure 1). Plate convergence in the study area started during Late Jurassic–Early Cretaceous times as the Neo-Tethys oceanic lithosphere initiated to be subducted beneath the Eurasian plate. Subsequent continent–continent collision occurred during the Early Miocene [1,2,3,4,5,6,7].

The substantial strike-slip component associated with oblique plate convergence was accommodated by the composite lineament including the Main Recent Fault (MRF) and the Main Zagros Fault (MZT). This lineament separates the Arabian Plate to the SW from the Sanandaj-Sirjan Zone to the NE (Iran Block, Figure 2) [17,18,19,20,21]. The Sanandaj-Sirjan Zone was formed mainly by metamorphic rocks, by Jurassic to Early Eocene calc-alkaline magmatic rocks, and by the products of Middle Eocene gabbroic plutonism [1,22,23,24]. SW of the MRF-MZT composite lineament, the typical succession of the Arabian margin (Upper Triassic–Quaternary), largely folded and faulted due to the continental collision, sits on top of crystalline basement [25,26,27,28,29,30,31]. Here, the High Zagros Fault (HZF) separates the Simply Folded Belt to the SW from the Imbricate Zone to the NE (Figure 2). The latter zone includes the telescoped distal part of the original continental margin of the Arabian Plate (e.g., [31]). The mountain chain is bounded to the SW by the morphotectonic feature known as Mountain Front Flexure (MFF). This represents the topographic evidence of the crustal ramp of the major thrust fault of the study area, which we term here Main Frontal Thrust (MFT). The MFT is a blind crustal fault forming an upper flat at the basement-cover interface SW of the MFF (Figure 2) [7,13,17,26,32,33]. Recent geodetic measurements show that the northward relative motion of the Arabian Plate is still active today and occurs at ca. 20 mm yr⁻¹ with respect to fixed Eurasia [8]. The MFT represents the main thrust in the outer part of the orogen and it is marked by intense seismic activity occurring at depths between 10 and 20 km [34], including the main seismic event (M_w = 7.3) of 12 November 2017 [13].

3. Methods

In order to model the thermal structure of the Lurestan region of the Zagros thrust belt, a 3-D geological model was built for the area, shown in Figure 3. Using an analytical procedure, we calculated a geotherm for a series of vertical lines (pseudo-wells, shown in Figure 3). The resulting geotherms were then interpolated, producing a series of 3-D isotherms.

3.1. The 3-D Geological Model

The 3-D geological model (Figure 4a) was implemented based on of the geological sections by Basilici et al. [12], Tavani et al. [13], Tavani et al. [14], and Vergès et al. [4], located in Figure 3. It results in a volume of 6,144,000 (320 × 240 × 80) km³ of rock. The topography was integrated in the model using a 30-m resolution ASTER GDEM. The resulting volume was divided into three different parts, assumed as homogeneous taking into account the resolution of thermal models at the crustal scale (e.g., [12]): (1) Arabian sedimentary cover, (2) Arabian basement, and (3) Iran block (or Sanandaj–Sirjan zone). The geological model (Figure 4b) took into account the MRF and the MFT, while minor faults were not included because of their negligible effect on temperature variations at the crustal scale. The Moho discontinuity was included in the model using the Moho depth calculated by Jiménez-Munt et al. [15] (Figure 4c).

3.2. Constraints and Assumptions

The heat sources considered in the computation of the thermal structure were: (1) The heat rising from the mantle and (2) the radiogenic heat due to the presence of radioactive elements in the crust [35]. In thrust belts such as the Zagros, a further heat source is represented by frictional heating associated with thrust faulting. In the calculation of this latter heat source the shear stress was considered, which was obtained by Sibson’s [36] formulation for favorable oriented thrust faults at the chosen depths under hydrostatic pore-fluid conditions with a fixed friction coefficient.

To perform temperature calculations for the study area, we considered a series of pseudo-wells (located in Figure 3 and Figure 5).

Input parameters included data on heat sources and geologically derived constraints (thrust depth, lithology) based on Tavani et al. [13]. A basal heat flow coming from the mantle was assumed, corresponding to a shield-type heat flow at the Moho [37]. For each pseudo-well (Figure 6) we considered the following parameters and assumptions:

(1)

Altitude. We took into account the topography using the 30-m resolution ASTER GDEM.

(2)

Thickness of the Arabian cover ($h_C$). We considered a variable thickness of the sedimentary cover based on the 3-D geological model (Figure 4), taking into account present-day topography (the model did not take into account erosion). The density of the Arabian cover was fixed to ρ_c = 2.55 × 10 ³ [38].

(3)

Thickness of the Arabian basement ($h_B$). We considered a variable thickness of the basement based on the depth of the Moho discontinuity calculated by Jiménez-Munt et al. [15]. The density of the basement was fixed to ρ_b = 2.80 × 10 ³ [38].

(4)

Thrust depth and angle. The pseudo-wells crossed the MFT at different depths with a certain angle.

(5)

Amount of basement offset by the MFT (Δz). Taking into account the reconstructed geometry of the Arabian margin by Le Garzic et al. [39], we considered a total shortening of 20 km in the range of 20–10 Ma. The resulting slip rate ($v$) was in the range of 1–2 mm yr⁻¹, and, therefore, we used an average value of 1.5 mm yr⁻¹ [12]. The resulting basement thickening was in the range of 3.5–5.8 km. The friction coefficient (μ) was fixed to 0.6, according to Byerlee [40,41].

(6)

Constant heat production rate for the sedimentary cover ($H_C$) = 1.0 µW/m³ [11].

(7)

Constant heat production rate for the basement ($H_B$) = 0.4 µW/m³ [11].

(8)

Thermal conductivity for the sedimentary cover ($K_C$) = 2.0 Wm⁻¹K⁻¹ [11].

(9)

Thermal conductivity for the basement ($K_B$) = 2.2 Wm⁻¹K⁻¹ [11].

(10)

Heat flux at the Moho ($Q_m$) = 20 mW/m² [12,37].

3.3. Analytical Procedure

To obtain the geotherms we separated the calculation into two parts, one for the sedimentary cover ($0\le z<h_C$) and the other for the basement ($h_C\le z\le h_C+{h\prime }_B$), with ${h\prime }_B=h_B+\Delta z$.

Basement thrusting, besides generating a new source of heat due to friction along the fault, produces a perturbation in the propagation of the mantle heat and an increase in radiogenic heat from the basement. Therefore, the basement temperature depends on the frictional heat (σv) and by the new equilibrium states reached both for the mantle heat ($Q_m$) considered as constant and for crustal radiogenic heat ($H_B$).

For the sedimentary cover, we considered the related radiogenic heat ($H_C$) and added the three perturbed conductive states due to: (1) The increase in radiogenic heat from the basement, (2) the perturbation of the mantle heat transfer, and (3) the frictional heat. Furthermore, we calculated the new perturbed states as a function of time based on the analytical equations by Molnar et al. [35], appropriately modified for overthrusting within a layered lithosphere [12,42].

Finally, the temperatures for the two, basement and sedimentary cover, layers can be obtained by the following equations:

$$T\left(z,t\right)=T_{DH}\left(z,t\right)+T_{R{H}_c}\left(z\right)+T_{MH}\left(z,t\right)+T_{FH}\left(z,t\right)$$

(1)

$$0\le z<h_C$$

(2)

$$T\left(z,t\right)=T_{DH}\left(t\right)+T_{R{H}_C}+T_{R{H}_B}\left(z,t\right)+T_{MH}\left(z,t\right)+T_{FH}\left(z,t\right)$$

(3)

$$h_C\le z\le h_C+{h\prime }_B$$

(4)

where $T_{DH}\left(t\right) \mathrm{and} T_{R{H}_C}$ are not depth-dependent terms, with $z=h_C$. Furthermore, each time-dependent term has the following generic expression:

$$T_{XH}\left(z,t\right)=T_f\left(z\right)-{\displaystyle \sum }_0^{\infty }{A}_n \sin (a_{n}z)e^{-a_n^{2}Kt}$$

(5)

with (

Table 1. Terms of the equations used for the analytical procedure.

		$0\le \mathit{z}\le {\mathit{h}}_{\mathit{C}}$	${\mathit{h}}_{\mathit{C}}<\mathit{z}\le {\mathit{z}}_0$	${\mathit{z}}_0<\mathit{z}\le {\mathit{h}}_{\mathit{C}}+\mathit{h}{’}_{\mathit{B}}$
$T\left(z,t\right)$		$T_{DH}\left(z,t\right)+T_{R{H}_C}\left(z\right)+T_{MH}\left(z,t\right)+T_{FH}\left(z,t\right)$	$T_{DH}\left(t\right)+T_{R{H}_C}+T_{R{H}_B}\left(z,t\right)+T_{MH}\left(z,t\right)+T_{FH}\left(z,t\right)$	$T_{DH}\left(t\right)+T_{R{H}_C}+T_{R{H}_B}\left(z,t\right)+T_{MH}\left(z,t\right)+T_{FH}\left(z,t\right)$
	$T_f$	$\frac{H_B{h}_B^{\prime }}{K_C}z$	$\frac{H_B{h}_B^{\prime }}{K_C}{h}_C$	$\frac{H_B{h}_B^{\prime }}{K_C}{h}_C$
$T_{DH}$	$a_n$	$\frac{\left(2n+1\right)\pi }{2h_C}$	$\frac{\left(2n+1\right)\pi }{2h_C}$	$\frac{\left(2n+1\right)\pi }{2h_C}$
	$A_n$	$\frac{2H_B\Delta z}{K_C{h}_C{a}_n^2}\sin (a_n{h}_C)$	$\frac{2H_B\Delta z}{K_C{h}_C{a}_n^2}\sin (a_n{h}_C)$	$\frac{2H_B\Delta z}{K_C{h}_C{a}_n^2}\sin (a_n{h}_C)$
$T_{R{H}_C}$		$\frac{H_C{h}_C}{K_C}z-\frac{H_C}{2K_C}{z}^2$	$\frac{H_C}{2K_C}{h_C}^2$	$\frac{H_C}{2K_C}{h_C}^2$
	$T_f$		$\frac{H_B{h}_B^{\prime }}{K_B}z-\frac{H_B}{2K_B}{z}^2$	$\frac{H_B{h}_B^{\prime }}{K_B}z-\frac{H_B}{2K_B}{z}^2$
$T_{R{H}_B}$	$a_n$		$\frac{\left(2n+1\right)\pi }{2h_B^{\prime }}$	$\frac{\left(2n+1\right)\pi }{2h_B^{\prime }}$
	$A_n$		$\frac{2H_B\Delta z}{K_B{h}_B^{\prime }{a}_n^2} \left[\sin (a_n{z}_0)+a_n\left(h_B^{\prime }-z_0-\frac{\Delta z}{2}\right)\cos \left(a_n{z}_0\right)\right]$	$\frac{2H_B\Delta z}{K_B{h}_B^{\prime }{a}_n^2} \left[\sin (a_n{z}_0)+a_n\left(h_B^{\prime }-z_0-\frac{\Delta z}{2}\right)\cos \left(a_n{z}_0\right)\right]$
	$T_f$	$\frac{Q_m}{K_C}z$	$\frac{Q_m}{K_B}z$	$\frac{Q_m}{K_B}z$
$T_{MH}$	$a_n$	$\frac{\left(2n+1\right)\pi }{2\left(h_C+h_B^{\prime }\right)}$	$\frac{\left(2n+1\right)\pi }{2\left(h_C+h_B^{\prime }\right)}$	$\frac{\left(2n+1\right)\pi }{2\left(h_C+h_B^{\prime }\right)}$
	$A_n$	$\frac{2Q_m\Delta z}{K_C\left(h_C+h_B^{\prime }\right)a_n}\cos (a_n{z}_0)$	$\frac{2Q_m\Delta z}{K_B\left(h_C+h_B^{\prime }\right)a_n}\cos (a_n{z}_0)$	$\frac{2Q_m\Delta z}{K_B\left(h_C+h_B^{\prime }\right)a_n}\cos (a_n{z}_0)$
$T_{FH}$	$T_f$	$\frac{\sigma v}{K_C}z$	$\frac{\sigma v}{K_B}z$	$\frac{\sigma v}{K_B}{z}_0$
	$a_n$	$\frac{\left(2n+1\right)\pi }{2\left(h_C+h_B^{\prime }\right)}$	$\frac{\left(2n+1\right)\pi }{2\left(h_C+h_B^{\prime }\right)}$	$\frac{\left(2n+1\right)\pi }{2\left(h_C+h_B^{\prime }\right)}$
	$A_n$	$\frac{2\sigma v}{K_C\left(h_C+h_B^{\prime }\right)a_n^2} \sin (a_n{z}_0)$	$\frac{2\sigma v}{K_B\left(h_C+h_B^{\prime }\right)a_n^2} \sin (a_n{z}_0)$	$\frac{2\sigma v}{K_B\left(h_C+h_B^{\prime }\right)a_n^2} \sin (a_n{z}_0)$

)

$$\mathrm{X}\equiv \mathrm{D}, \mathrm{R}, \mathrm{M}, \mathrm{F}\hspace{1em}{a}_n, A_n$$

(6)

where $T_f\left(z\right)$ is the final temperature of the new state of equilibrium ($t\to \infty )$, and $T_{DH}\left(t\right)$ and $T_{R{H}_C}$ are the respective temperatures calculated at the top of the basement. In particular, $T_{R{H}_C}\left(z\right)$, the only not-time-dependent term, is the temperature due to the radiogenic heat in the sedimentary cover, while $T_{DH}\left(z,t\right)$ is the temperature due to the propagation of the radiogenic heat in the lower (i.e., basement) layer considering the heat production due to overthrusting, $T_{R{H}_B}\left(z,t\right)$ is the temperature associated with the perturbed radiogenic heat in the basement, and $T_{MH}\left(z,t\right)$ and $T_{FH}\left(z,t\right)$ are the temperatures associated with the new perturbed state of the mantle and of the new heat source due to overthrusting.

To calculate the isotherms along the considered sections (AB, KH, WJ, and DC of Figure 5) we interpolated model data by a second-order power-law polynomial equation, obtaining the best-fit analytical curves approximating the geothermal values.

4. Results

Using the analytical procedure described in the previous paragraph, we calculated the geotherm for each pseudo-well shown in Figure 5. The result was a series of 20 geotherms, shown in Figure 6.

The geotherms relative to pseudo-wells 5, 10, 15, and 20 (yellow lines in Figure 7) showed a constant increment of temperature with depth due to the absence of basement offset by the MFT. On the other hand, all other pseudo-wells showed an increase of temperature variation with depth in correspondence to the MFT. In particular, the different trend of these geotherms depended on both the basement thickness modified by thrusting ($\Delta z$) and on the thrust depth. A weaker temperature increase occurred for the geotherms of pseudo-wells 4, 9, 14, and 19 (green lines in Figure 7), where the thrust was shallower and implied a smaller amount (3.5 km) of basement offset. All other geotherms showed higher temperature variations due to larger basement offsets ($\Delta z$ in the range of 5.3–5.8 km) associated with thrusting occurring at progressively greater depths (from 18 to 36 km) moving towards the hinterland (i.e., from SW to NE).

The sensitivity of this type of analytical procedure was checked by Megna et al. [43]. A variation of about 10% of thrust depth or slip rate produced a change of frictional heat that modified the final temperature by values never exceeding 2%.

In order to provide an effective overview of the temperature trend along each cross section (AB, KH, WJ, and DC of Figure 5), we interpolated the temperature data from 50 °C to 500 °C using a 50 °C step. The isotherms obtained by interpolating the data using a second-order power-law polynomial equation consisted of best-fit analytical curves with an rms (root mean square) error comprised between 0.6 and 3.75 (for the 100 °C and 500 °C isotherms, respectively).

The isotherms shown in Figure 8 (where a temperature step of 100 °C was used to allow a better visualization) define a general trend of increasing temperature moving from the foreland to the hinterland (i.e., from SW to NE) of the chain along each of the four cross-sections. Only the isotherm pattern along cross-section AB defined an upward curvature, thus confirming the previous 2-D study of Basilici et al. [12]. The isotherm pattern along cross-sections HK and JW show a gentle downward curvature, while more straight, SW-ward dipping isotherms characterize cross-section CD.

5. Discussion

Based on the analytical results, isothermal surfaces were constructed using the minimum curvature spline interpolation technique. This allowed us to produce a series of maps for the 100 °C to 500 °C isotherms (Figure 9), as well as a 3-D model of the thermal structure of the Lurestan arc subsurface (Figure 10). A comparison of the latter with the thermal structure obtained by second-order power-law polynomial interpolation (Figure 11) shows a good fit between the two methods, thus confirming the consistency and robustness of the thermal structure depicted in Figure 10.

The isotherm depth contour maps of Figure 8 are characterized by a general SW-ward deepening of the isotherms, which clearly outline a higher temperature axial zone of the mountain belt and a “colder” outer zone including the foothills of the Simply Folded Belt and the adjacent foreland. The NE-ward temperature increment was mainly due to the crustal thickening associated with basement-involved faulting, as the NE dipping MFT produces a progressively increasing amount of basement offset by the fault ($\Delta z$). The SW-ward deepening of the isotherms involved depth changes in the range of ca. 7 km to ca. 20 km for the 100 °C and 500 °C, respectively. The resulting thermal structure tended to mimic the general arc-shaped geometry of the mountain belt, as readily observable by comparing the isotherm depth contours with the trace of the MFF in the maps of Figure 9. As a matter of fact, the junction between the Lurestan arc (forming a salient of the Zagros fold and thrust belt) and the Kirkuk embayment (representing a recess) to the NW was clearly reflected by the crustal thermal structure observed in map view (Figure 9). The thermal structure shows a clear correlation with the geometry of the main thrust fault of the region, i.e., the MFT. According to Tavani et al. [14], the MFT developed in the necking domain of the Jurassic rift system ahead of an array of inverted Jurassic extensional faults, the related structural architecture resembling that of a crustal-scale footwall shortcut. As the main crustal ramp of the MFT underlies at depth the MFF, the sinusoidal shape of the latter in the Lurestan region would derive from the re-use of the originally segmented, inverted Jurassic rift system [14]. Despite this long-lived paleo-tectonic control and the fundamental role of structural inheritance, the development of the MFF, and thus the activity of the underlying thrust (MFT) controlling it, are recent processes. These are dated to ca. 3 Ma based on stratigraphic information [44], and to ca. 5 Ma based on low-temperature thermochronology data [7]. According to the latter studies, the major morphotectonic feature constituting the MFF developed during the Pliocene by basement-involved thrusting as deformation migrated downward at deeper crustal levels. Therefore, late-stage, crustal ramp-dominated thrusting [45] not only led to the development of a prominent geomorphological boundary between the high Zagros Mountains and the low foothills to the SW [17,32,33] but also appears to have exerted a major control on the thermal structure of the Lurestan salient.

The first-order, general pattern of foreland-ward deepening of the isotherms was rendered more articulated by along-strike variations, mainly consisting of a SE-ward climbing of the isotherms (Figure 9). Also, this along-strike, SE-ward temperature increase was essentially controlled by the 3-D geometry of the MFT, which is characterized by a low-angle crustal ramp in its SE portion, changing to a steep ramp rapidly reaching lower crustal depths to the NW (Figure 10). This is essentially because the gently dipping MFT segment produced a larger basement offset in the SE sector, where it joins the MRF along a branch line located within crustal depths (Figure 4 and Figure 10). This larger basement offset resulted in a more marked temperature increase and related shallower isotherms in the SE sector of the Lurestan arc (Figure 9). The along-strike variation in the thermal structure is readily observable in Figure 11. There, the widely spaced isotherm pattern of the NW “colder” edge of the model (i.e., at the junction with the Kirkuk embayment) gives way to a pattern characterized by progressively more densely spaced isotherms, which become closer to each other moving to the SE portion of the Lurestan arc. As a result, the easternmost sector of the study area represents a zone of relatively high geothermal gradient, particularly concerning temperatures <300 °C in the shallower part of the upper crust (Figure 9). This could have interesting implications in terms of geothermal energy potential in the region.

6. Conclusions

This study provides the first 3-D model of the thermal structure of the Lurestan arc of the Zagros thrust belt. Our results emphasize the major role played by a crustal thrust fault (i.e., the MFT) in controlling the thermal setting of the region. This is testified by a very good correlation between high temperature gradients in the crust and the amount of basement offset by the MFT. Besides a general temperature increase towards the hinterland, manifested by shallower and more densely spaced isotherms moving to the NE, there is also a marked along-strike variation associated with a lateral transition from a gently to a moderately dipping crustal ramp of the MFT. This results in a zone of relatively high geothermal gradient in the easternmost part of the study area. The marked effect of recent (i.e., Pliocene) thick-skinned inversion [14] and associated basement thrusting in perturbing the thermal structure, particularly at upper crustal levels, represents an interesting outcome to be taken into account in future geodynamic studies, as well as in geothermal/hydrocarbon exploration in the Lurestan region. It may also have important implications for active tectonic studies and seismotectonic modelling in an area recently affected by a M_w > 7 earthquake.