Slip rate of the Kazerun Fault and Main Recent Fault (Zagros, Iran) from 3D mechanical modeling

Journal of Asian Earth Sciences 41 (2011) 89–98

Contents lists available at ScienceDirect

Journal of Asian Earth Sciences journal homepage: www.elsevier.com/locate/jseaes

Slip rate of the Kazerun Fault and Main Recent Fault (Zagros, Iran) from 3D mechanical modeling H.R. Nankali ⇑ National Cartographic Center of Iran, Geodesy and Geodynamics Department Tehran, Iran

a r t i c l e

i n f o

Article history: Received 13 January 2010 Received in revised form 14 October 2010 Accepted 14 December 2010 Available online 26 January 2011 Keywords: Zagros Finite-element methods Shortening GPS Iran

a b s t r a c t The aim of this work is, using GPS velocities and 3D mechanical modeling including preexisting faults with frictional properties that obey Coulomb-failure process to explore the long-term slip rates of the MRF fault and Kazerun Faults in the Zagros mountain of Iran. A model was constructed using spatially varying crustal thickness, geothermal gradient, and two major faults. Structural boundaries are derived from the several deep seismic soundings carried out in the area. Far field GPS velocities and late Quaternary fault slip rates are used to compare the model results. Rheological tests show that effective fault friction lower than 0.05 leads to high slip rates that fit with geologically and geodetically determined slip rates of MRF and Kazerun Fault. These results suggest that MRF and Kazerun Fault are weak faults like San Andreas and North Anatolian fault. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction Among the present active belts, the Zagros mountain fold-and thrust belt, located in Southwest Iran, is one of the youngest and most seismically active in the world. It results from the collision of the Arabian plate with the continental blocks of Central Iran, beginning in Miocene time and continuing today (e.g. Stocklin, 1968). The main feature of the Zagros is a linear, asymmetrical folded belt, which forms a 200–300 km wide series of ranges extending for about 1500 km from eastern Turkey to the Strait of Hormoz. The Zagros contains a thick and almost continuous sequence of shelf sediments deposited on the 1–2 km thick infraCambrian Hormoz Salt formation. The sediments, of Paleozoic to Late Tertiary age, are believed to be separated from the Precambrian metamorphic basement by this Hormoz salt (Berberian and King, 1981). The Main Zagros Reverse Fault (MZRF) separates the Zagros mountain belt from Central Iran, which is a major structural discontinuity (Stocklin, 1968). The NW trending Zagros fold-and-thrust belt is affected by two major dextral faults: (1) the NW trending Main Recent Fault that accommodates partitioning of oblique convergence at the rear of the western Zagros. The Main Recent Fault (MRF) trends NW–SE and forms the northeastern border of the Zagros Mountains (Tchalenko and Braud, 1974) and (2) the north trending Kazerun Fault located in the central Zagros. The assessment of the slip rates of these faults is important to understand the deformation behavior ⇑ Fax: +98 2166071000. E-mail address: [email protected] 1367-9120/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.jseaes.2010.12.009

of the crust and mantel beneath this area and because fault slip rate is directly related to the repeat time of earthquakes rupturing the fault, its knowledge is critical for the assessment of seismic hazard. Both geodetic and geologic methods are used to determine slip rates. Recent GPS studies in Iran (Tavakoli et al., 2008; Walpersdorf et al., 2006; Hessami et al., 2006; Vernant et al., 2004) have produced new data on the kinematics of continental deformation in the Iran and Zagros region, see (Fig. 1). The maximum slip rate along the MRF deduced from GPS measurements would be of 4 ± 2.5 mm/yr if the fault achieves complete partitioning of the 4–7 ± 2 mm/yr of north trending shortening across the western Zagros (Vernant et al., 2004). This is not in agreement with a Pliocene (3–5 Ma) initiation of the MRF (Talebian and Jackson, 2002) and with a cumulative lateral slip of 50 km along that fault, which will lead to a long-term slip rate of 10–17 mm/yr (Talebian and Jackson, 2002) and also with the estimate of 10 mm/yr of (Bachmanov, 2004) based on the offset of a river valley incised into a surface of likely postglacial age. The Zagros is divided into two distinct along-strike blocks by the north–south-trending Kazerun Fault. The Kazerun Fault is one of the basement structures inherited from a Neo-Proterozoic tectonic phase (Talbot and Alavi, 1996). Present-day activity of the Kazerun Fault is evidenced by recent earthquakes with right-lateral mechanisms located on the fault (Baker et al., 1993; Yamini-Fard et al., 2007). It extends from the southeastern termination of the MRF, in the north, to the Persian Gulf to the south. GPS measurements provide a precise estimate of the presentday dextral strike–slip motion on the Dena (3.7 ± 0.8 mm/yr) and Kazerun (3.6 ± 0.6 mm/yr) segments (Tavakoli et al., 2008). These

90

H.R. Nankali / Journal of Asian Earth Sciences 41 (2011) 89–98

Fig. 1. Framework of the Arabia–Eurasia collision in Iran, arrows displays the GPS velocity field as computed by (Vernant et al., 2004), thick solid lines represent active faults, and thin solid lines match political boundaries.

velocities indicate that the Dena, Kazerun and Kareh Bas faults accommodate most of the differential motion between North and Central Zagros, with 3–6 mm/yr and 10 mm/yr shortening across North and Central Zagros, respectively. Effectively, the 3–6 mm/ yr NS shortening across North Zagros are increasing from 3 mm/ yr in the west to 6 mm/yr in the east, close to the Kazerun Fault System. Horizontal slip rates have been estimated on these two faults over the last 140 ka from lateral offsets of streams and fans and in situ cosmogenic 36Cl exposure dating of cobbles sampled on the surface of these geomorphic features. Compared to GPS data, the obtained minimum slip rate of 3.5–12.5 mm/yr on the Main Recent Fault implies strike-slip partitioning of the convergence along this fault. Minimum slip rate of the Kazerun Fault is 2.5– 4 mm/yr for its northern strand, 1.5–3.5 mm/yr for its central segment and is negligible for its southern segment (Authemayou et al., 2009). These results are consistent with southward distribution of the slip from along the Main Recent Fault to the longitudinal thrusts and folds of the fold-and-thrust belt through the Kazerun Fault, with a decrease of slip from the southeastern tip of the Main Recent Fault towards the southern termination of the Kazerun Fault. Because fault slip rate is directly related to the repeat time of earthquakes rupturing the fault, its knowledge is critical for the assessment of seismic hazard. The idea of evaluating the geological reliability of the different slip rate estimates of the two faults is important. This study presents another proxy that allows us to better constrain fault slip rates, which help us to increase our knowl-

edge and understanding for assessment of seismic hazard and distributed deformation in this part of Iran. Using GPS, heat flow data, seismic and faults (MRF and Kazerun) data available in the Zagros, I use mechanical modeling to investigate the mechanical relation between far field loading of Arabian and Eurasia and slip rate along the MRF fault and Kazerun Fault in Zagros. 2. Previous images of deep structures in the Zagros Published data on the Zagros crustal structure including moho depth estimates are very scarce. The only available profiles of crustal thickness variations have been computed from Bouguer anomaly modeling by Dehghani and Makris, (1984) for whole of Iran and Snyder and Barazangi (1986) and for Zagros. Hatzfeld et al. (2003) estimated a crustal thickness of 46 km from receiver function computed at a single station close to the town of Ghir in central Zagros. Based on microearthquake survey and receiver functions for the eastern part of the central Zagros Hatzfeld et al. (2003) suggested, a crustal structure consisting of an 11 km thick sedimentary layer overlying an 35 km thick crystalline crust. The crystalline crust consists of an upper layer extending from 11 to 20 km depth and a lower layer extending from 20 km depth to a crust mantel interface at 46 km depth. Moreover, in that region most well-located microearthquakes occurred between 10 and 14 km depth, with none occurring deeper than 20 km (Tatar et al., 2004; Engdahl et al., 2006), see (Fig. 2). Using the waveform-modeled depths in the Zagros, Talebian and Jackson (2004) found no earthquakes

H.R. Nankali / Journal of Asian Earth Sciences 41 (2011) 89–98

Fig. 2. Earthquake depth distribution in Zagros (Engdahl et al., 2006) depths determined by waveform modeling are not shaded. Shaded are EHB focal depths.

deeper than 20 km anywhere except near the Oman Line in the extreme SE Zagros i.e., beneath the Zagros there is no evidence in the form of mantle earthquakes for present-day active subduction of continental crust, with shortening apparently accommodated entirely by crustal thickening and distributed deformation. Analysis of the seismological network data across central Zagros by (Paul et al., 2006) show Moho depth variations with a lateral resolution of a few km along a 620-km-long profile transverse to the range. The average crustal thickness is 45 km beneath most of the Zagros fold-and-thrust belt and 42 km beneath the Urumieh–Dokhtar volcanic zone and the southern part of the Iranian microcontinent. The region between the suture zone in the High Zagros to the Sanandaj–Sirjan metamorphic zone is characterized by a marked crustal thickening from 45 to 70 km in a narrow region. Bird (1978) reached similar conclusions using 2-D finite element modeling of lithospheric deformation in the Zagros assuming pure anelastic rheology. He showed that the shear stress concentrates on the suture zone for all slab–pull boundary conditions whatever the densities and flow laws considered. The thin viscous sheet model has been applied to Iran by Jackson (1995) and Sobouti and Arkani-Hamed (1996). Thin viscous sheet (TVS) models of lithospheric deformation were introduced to the study of continental tectonics about two decades ago by England and McKenzie (1983). In these models the lithosphere is permitted to thicken or thin, or undergo horizontal displacement, but horizontal displacement vectors are assumed invariant with depth in the lithosphere. Assuming also that tractions on vertical planes only act in the horizontal directions and those on horizontal planes only act in the vertical direction allows the force-balance equation to be reduced to a 2D viscous flow problem. Both of these studies concluded that deformation of the Iranian lithosphere is determined by the shape of the rigid boundaries and the disposition of the rigid blocks within the collision zone. Vernant and Chery (2006) use a simple 2.5 D mechanical modeling for the Zagros and shows that for the case of the NW Zagros, Main Recent Fault can accommodate only 25% of the whole tangential motion. As the GPS inferred range parallel motion across the NW Zagros is 5 ± 2 mm/yr, the present day slip rate of this fault could be as low as 1.2 mm/yr. Their model suggests that no strikeslip activity takes place on a vertical fault when the direction of convergence is more frontal. Vernant and Chery (2006) using a simple 2 layer numerical modeling with linear Maxwell rheology

91

and elastic constant show that the mechanical behavior of the Iranian lithosphere seems to be partly controlled by the large strike-slip faults. However, some deformation is still taken up by the continuum medium, suggesting a compromise between the microplate and continuum descriptions. Liu and Bird (2007) use a kinematic finite-element method to derive a self-consistent long-term-average velocity field by combining space geodesy, geological slip rates and principal ‘stress directions’ (e.g. seismic moment tensor orientations) in PersiaTibet and Burma Orogen. They show that Arabia–Eurasia relative plate motion in eastern Iran (east of 58°, E) is mainly accommodated by the Makran subduction zone and Kopet Dagh. Despite little constraint from input geological slip rates, their model predicts the correct sense (compared to earthquake focal mechanisms) for strike-slip faults on the west and east sides of the Lut block. The central Iran and Lut blocks appear to be relatively rigid, with strain localization on their boundary faults. West of 58°, E, Arabia–Eurasia convergence is accommodated by the Zagros and Alborz mountain ranges. Output slip rates of 16 mm/yr on the Main Recent Fault are compatible with input slip rates of 10–17 mm/yr (Talebian and Jackson, 2002) but much larger than geodetic slip rates of 3 ± 2 mm/yr (Vernant et al., 2004). They suggest that a transition from strong to weak coupling between the underthrusting Arabian plate and the Zagros in western Iran and Makran subduction in the east may affect surface deformation patterns in the Zagros and the Makran subduction zone. Previous numerical models (Vernant and Chery, 2006; Bird et al., 1975) of Zagros deformation have been constructed using the assumption of a continuum deformation, but without direct constraints on the velocity field and seismic data structure in Zagros. Using a combination of viscoelastic rheology for the lithosphere continuum and coulomb friction condition for the fault motion, I study the mechanical processes controlling fault slip rate and deformation of bounding crust.

3. Modeling approach I use ANSYS (finite element program) to solve mechanical equilibrium. ANYSYS employs the Newton–Raphson approach to solve non-linear problems. In this method a load is subdivided into a series of increments applied over several steps. Before each solution the Newton–Raphson method evaluates the out-of-balance load vector, which is the difference between the restoring forces and the loads corresponding to the element stresses and the applied loads (such as displacement). A linear solution is performed using the out of balance loads and is checked for convergence. If convergence criteria are not satisfied the out-of balance load vector is reevaluated, the stiffness matrix updated and a through direct elimination of equations until the problem converges. 3.1. Model formulation The three-dimensional model covers a wide rage surrounding the MRF and Kazerun Faults. The model incorporates four key features: (1) The Arabian plate moves towards the Eurasian plate at a rate of 22 mm/yr at Bahrain position (26.209°, 50.608°) on a N10°_E oriented vector (e.g., Sella et al., 2002; Vernant et al., 2004), the GPS velocity field at large distance from the fault is not influence by the seismic cycle and is therefore a proxy of continental long term surface motion. This hypothesis seems reasonable regarding the fair agreement between geological and corresponding geodetic plate motion (Sella et al., 2002).

92

H.R. Nankali / Journal of Asian Earth Sciences 41 (2011) 89–98

(2) The Zagros region is crossed by near-vertical strike-slip faults that cut through the crust (Vernant and Chery, 2006). (3) the crustal velocity structure is known (Hatzfeld et al., 2003; Tatar et al., 2004). (4) The crustal geotherm is derived from initial thermal condition of the model. The finite element model of the Zagros region was built in three compositional layers inferred from measured crustal velocity structure (e.g., Hatzfeld et al., 2003). Horizontal boundaries divide the crust into the upper and lower part. This division is approximately based on the P-wave velocity isolines for the boundary between the upper and the lower crust, and the Moho boundary for division between the crust and the mantle. We have tried to keep the model structurally as simple as possible and therefore the finest details found in the seismic velocity sections are excluded. Rock composition determines material elastic constants and, consequently, the strain response to an imposed stress at a given rate and temperature. According to Hatzfeld et al. (2003) the upper 20 km of the model represents the upper crust. Elastic properties of these rocks were approximated by granite and are given in Table 2. The lower crust is 28 km thick in the model and has elastic properties representative of basalt-diabase composition. The upper mantle layer has properties associated with dunite samples (Table 2) and is 22 km thick. All elements of the Zagros finite element model strained through a combination of linear elasticity (Hook law) and rate dependent creep behavior. Time-independent elastic strain occurred in the model according to

e ¼ r=E;

ð1Þ

where E is Young’s modulus and r is differential stress and e is the elastic strain. Modeled time-dependent inelastic strain rate was controlled by the creep equation (e.g., Kirby, 1983)

e_ ¼ A expðQc=RTÞrn ;

ð2Þ

where A (elastic constant), Qc (activation energy), and n are experimentally derived elastic constants, R is the universal gas constant, T is temperature, and r is differential stress (Table 2). Eq. (2) adds an important contribution to modeled strain with increasing temperature, and simulated a smooth transition from elastic to ductile strain with depth. At low temperatures the exponential term of Eq. (2) approached unity, causing the creep term to become vanishingly small. With increasing temperature, the contribution of the exponential term became more important. The finite element model has cuts in it that represent the two major strike-slip faults of the Zagros i.e. Main Recent Fault and Kazerun; The faults are deformable and are constructed from contact elements that obey the Coulomb failure relation.

Table 1 Thermal Properties of the model. Surface temperature (K)

Ts

273

Lithospheric basal flow (W m2)

q

0.019

Upper crust Thermal conductivity (W m1 k1) Radiogenic heat Production (W m3)

K H

3 1.5  106

Lower crust Thermal conductivity (W m1 k1) Radiogenic heat Production (W m3)

K H

2.5 0.5  106

Upper mantle Thermal conductivity (W m1 k 1) Radiogenic heat production (W m3)

K H

3 0.02  106

CF ¼ s þ lðrn Þ;

ð3Þ

where s is the shear stress acting on a fault surface, l is the friction coefficient, and rn is the stress component normal to the fault surface. Contact elements have zero thickness and are welded to the sides of elements. Faults are permitted to extend down to the base of the model. Our mesh covers a rectangular area from south east to northwest Iran with horizontal dimensions of 1500 km  600 km (Fig. 4). It consists of 28,125 elements with 119,260 active nodes. The mesh is finest where volumes change shapes the most, and in regions of greatest complexity such as fault terminations or intersections.

4. Heat flow/power law rheology Conductive heat transfer, understood to be the main mechanism of the heat transport in the lithosphere, is described with Fourier’s law of heat conduction where the heat flow is directly proportional to the temperature gradient

rðr  TÞ ¼ A=k;

Q ¼ k

@T ; @z

ð4Þ

where Q is the near surface heat flow (W m2), T is the absolute temperature (K), k = k(z) is the coefficient of thermal conductivity (W m1 K1)as a function of depth and temperature (Artemieva and Mooney, 2001) and z is the depth (m) and A = A(z) is the heat production as a function of depth (z). This equation can be used to determine the heat flow density (HFD) with temperature measurements in boreholes as a function of a depth, together with the laboratory estimates of thermal conductivity. If we assume steady state conditions, we end up with an equation that can be used to calculate the stable geotherms of the lithosphere. Temperature has a major importance in rheological and the tectonic modeling as, for instance, ductile flow is dependent on temperature in a non-linear way. Thus, the temperature in the lithosphere must be known in order to resolve the rheological structure properly, so it’s necessary to simulate the transition from elastic, seismic deformation into viscous, aseismic deformation. The thermal field in Zagros is difficult to determine, because the progressive crustal thickening increases the radiogenic production inside the belt, second, the onset of the Arabian mantle subduction below the Zagros could cause a temperature decrease in the Zagros crust above (Vernant and Chery, 2006). Bird et al. (1975) compute a 2D simple thermal model for continental areas such as Zagros and show the temperature distribution in this area. Estimated heat flow is between 37 mW m2 and 45 mW m2. This result is consistence with the obtained surface heat flow measured in 3–5 km deep oil well drilled by National Iranian Oil Company (NIOC) Bird et al. (1975). They did not talking about the correlation between surface heat flow and seismicity depth, but this thermal filed appear compatible with a temperature of 350 °C occurring between 10 to 20 km. I compute 1D steady state thermal models according to the present-day Zagros crustal structure. The radiogenic heat production is 0.02 lW m3 taken for the mantle, and a basal heat flux of 0.019 mW m2 coming from the mantle is assumed (Artemieva and Mooney, 2001). We use radiogenic productions for upper and lower crust of 1.5 and 0.5 lW m3, respectively. In the Zagros finite element model, the proportion of inelastic to elastic behavior of a given element node is governed by a crustal geotherm derived from surface heat flow analysis (Figs. 5 and 6).

93

H.R. Nankali / Journal of Asian Earth Sciences 41 (2011) 89–98

Table 2 Material constants used in the three layers of the finite element model elements in the model are all viscoelastic, with viscous behavior controlled by the temperature gradient. Material properties of the model Parameter

E Young’s modulus (MPa) A physical constant (MPan S1) n Power law exponent Q activation energy (kJ/mo1) Poisson’s ratio (m) Density103 (kg/m3 q)

Layer 1

Layer 2

Layer3

0–20 km

Source

20–48 km

Source

>48 km

Source

104  6 104  2 1.9 140.6 0.25 2.775

1 2 2 2 3 4

104  9.4 102  6.3 3.1 276 0.26 2.905

1 5 5 5 3 4

105  1.6 103  5 3.8 492 0.28 3.200

1 6 6 6 3 4

Fig. 3. Geology and tectonic setting of Zagros, thick solid line shows Area of this study.

Fig. 4. The finite element mesh, model boundaries, and fault structures. Elements are all viscoelastic. The faults are discontinuities in the model and have deformable contact elements on their surfaces that obey a Coulomb failure criterion.

94

H.R. Nankali / Journal of Asian Earth Sciences 41 (2011) 89–98

parameter variation has been kept to a minimum in horizontal directions. Elastic material parameters, i.e., Young’s Modulus and Poisson’s ratio (Table 1) are taken from the ratio of P- and S-wave velocities along seismic profiles (Snyder and Barazangi, 1986) were used to confine the Poisson’s ratio. Values of densities (Table 1) are taken from the velocity models of seismic profiles and from the gravity model of (Snyder and Barazangi, 1986) for the Zagros profile. A simplification is made, as densities are kept constant in each horizontal layer, i.e., in the upper crust 2775 kg m3, in the lower crust 2950 kg m3, in the mantle 3200 kg m3. Thermal conductivity follows the temperature dependence with different starting values in the crust and in the mantle (Table 1). Value of the heat production in the upper crust is 1.5 lW m3 and in the lower crust is 0.5 lW m3. In the mantle values 0.02 lW m3 is used. 6. Numerical experiments

Fig. 5. The geothermal gradient used in the model input to the creep equation that controls material behavior. Model elements behave according to Hook’s law where temperatures are low and become transitionally more viscoelastic at higher temperatures. Horizontal axis represents the temperature in Kelvin and the vertical axis represents depth in km.

To show whether the contemporary far field motion of the Arabia–Eurasia convergence could result in the long-term fault slip rate average over several seismic cycle fitting the geological rate, we run series of models with effective fault friction coefficients from 0.02 to 0.3. In all the following numerical experiments the effective fault friction is the only parameter that varies. We present five experiments; for five of them the effective fault friction is the same in the crust and mantel (case#1: l = 0.02, case#2: l = 0.05, case#3: l = 0.1, case#4: l = 0.2, case#5: l = 0.3). 7. Boundary conditions, loads

Fig. 6. Initial surface heat flow density in the Zagros. The surface heat flow is 0.063 mW m2 Vertical axes represent depth in km and horizontal axes represent heat flow.

5. Material parameters Material parameters are chosen so that they are representative of typical rock types occurring in the continent. Felsic composition for the upper crust, mafic for the lower crust and ultra mafic for the mantle. I have approximated these compositions with creep parameters from (Carter and Tsenn, 1987) see (Table 1). Material

The geographical boundaries are superimposed on the topographic map of (Fig. 3). Our model boundaries are parallel and perpendicular to the Zagros belt trend, i.e. extends from North Western Iran (West) near Urumieh to Southeast (East) end of the Zagros (transition zone). It is approximately 1600 km  600 km. The first step in each experiment is that running the model with gravity, which compressed the model and establish initial stress state. The bottom of the model was constrained to zero displacement in the vertical direction, and the model sides were not permitted to move laterally. The model elements had non-linear time dependence; thus the rate of gravity implementation was important. If the model was subjected to gravity too rapidly, then deformation rates higher than could be supported by the defined material time dependence ensued. In this study, gravitational loading was ramped over 20 k yr, and then an additional settling period of 5 k yr was applied to ensure that the model was fully compressed. Time steps were initially 1 yr and gradually increased by a factor of 1.5 per iteration. Next continue with tectonic loading which, simulated by moving the western edge of the model according to the convergence rate between Arabia and central Iran according to (Vernant et al., 2004). The eastern model edge in X and Y direction is held fixed due to the central Iran and is thus not free to move laterally. The model base is freely slipping laterally but cannot move vertically. The northern and southern model faces are deformable only parallel to their surfaces. The model free surface is fully deformable. All constraints are imposed on the model edges as described above; no constraints are imposed on elements within the model. The results shown below represent the steady state values. 8. Numerical results-fault slip rates The finite element model was constrained to nearly match observed long-term fault slip rates within the given error ranges. Long-term fault slip rates derived from geological observations

H.R. Nankali / Journal of Asian Earth Sciences 41 (2011) 89–98

in the Zagros area were compiled by the several researches (e.g. Berberian and King, 1981; Talebian and Jackson, 2002; Authemayou et al., 2009). The fault slip rates estimated by model are listed in Table 3. Strictly speaking GPS and modeled velocities are not directly comparable, because GPS measurements are made during an interseismic period when faults are locked. However, it is usually assumed that the shape of interseismic motions around continental faults is controlled by a locking depth of about 12– 15 km. In this case, the interseismic motion at some distance from the fault (i.e., more than 50 km) is close to the geological slip rate (Reilinger et al., 2006; Savage and Burford, 1973). The variation of fault slip rate resulting from the change in fault friction only weakly affects the GPS site velocities, and fault friction variation does not result in significant velocity changes at the GPS sites, it affects the fault slip rate distribution throughout the model. Therefore, I use the geologically estimated fault slip rate (Authemayou et al., 2009) as a control parameter for model results. The estimated slip rates of the each fault in numerical experiments are shown in Table 3.

9. Effects of model assumptions Results of numerical rheological modeling are subjected to many uncertainties and thus a brief review of the most important error sources is given next. Fernàndez and Ranalli (1997) have divided the uncertainties related to rheological calculations into two groups: operational and methodological errors. Operational errors include inaccurate estimation of lithospheric structure and composition, errors in temperature, improper material and rheological parameters, pore fluid pressure and friction coefficient, etc. Methodological errors originate from the basic assumptions used in the rheological calculations. They include, for example, inaccurate rheological flow laws both in the brittle and ductile regimes and also assumption of the uniform strain and constant strain rate. The measured surface HFD value is the parameter that indicates what the deep geothermal structure is like. The surface HFD is quite often used as a constraint when geotherms are derived in a numerical or analytical way. The surface HFD values are evaluated from temperature measurements taken in the boreholes and large uncertainties can be connected to these estimations. Quantitatively speaking variation in the surface HFD with an amount of 10 mW m2 can result in a temperature change of 100 °C at the depth of 50 km Moho depth. Thermal conductivity and heat production values are also required in lithospheric thermal modeling. These values are derived from laboratory measurements of rock samples but they usually only correspond to the conditions in the uppermost crust. For the rest of the lithosphere they are usually indirectly estimated, e.g., from relationships between the heat flow - heat production, the seismic P-wave velocity. From the methodological errors, the inaccuracy of rheological flow laws must be kept in mind as those expressions are derived from experimental data, which are furthermore extrapolated to

Table 3 Fault slip rates estimate from model with different friction coefficient. Experiment

Fault friction in the crust

Fault friction in the mantel

MRF slip rate from model (mm/yr)

Dena slip rate from model (mm/yr)

Kazerun slip rate from model (mm/ yr)

Case#1 Case#2 Case#3 Case#4 Case#5

0.02 0.05 0.1 0.2 0.3

0.02 0.05 0.1 0.2 0.3

2.3 2 1.5 1.2 1

4.8 3.7 2.5 1.3 0.9

2.1 1.6 1.3 1.1 1

95

lithospheric conditions. The brittle failure criterion of Byerlee’s law has been derived from results that are confirmed only to the depths of the upper crust. Therefore, using it to much larger depths results in strengths that are most likely somewhat unrealistic. In summary, an attempt was made to include as much detailed information on crustal structure and rheology, as possible, but calculated results should be interpreted with the caveat of considerable uncertainty.

10. Discussion During the last two decades, the debate of the long-term behavior of the Zagros and the kinematic of large faults has been obscured by the lack of a common mechanical viewpoint. As summarized by Thatcher (2007), the lithosphere can be viewed as two end-member models. One behavior consists of rigid blocks separated by major faults as in the global tectonic plate model (block model). This model allows in principle to accumulate all the motion provided by the boundary conditions on a few faults. Although this model has a kinematical origin, it corresponds to a mechanical view in which weak faults are surrounded by strong lithospheric blocks. Alternatively, continuous strain may also occur in the lithosphere according to a fluid-like behavior of a viscous material (continuum model). In this case, the strain in the brittle crust is distributed on many faults and is therefore close to a continuous behavior. In this view, most of the lithospheric strength is taken by ductile layers in the lower crust and the uppermost mantle. The model I use in this paper has rheological properties of both block model and continuum model. First, the use of faults with an a priori fault strength using an effective friction allows to create a discontinuous motion if the friction is low (block model behavior). Conversely, a high friction will lead the fault to remain inactive, the strain being kept inside the continuum medium (continuum model behavior). The 3D visco-elastic finite element model presented here provides a useful tool for simulating long-term lithospheric deformation with faults. I applied this technique to investigate deformation associated with of major strike-slip faults. As mentioned before, fault friction variation does not result in significant velocity changes at the GPS sites, it affects the fault slip rate distribution throughout the model. Therefore, I use the geologically estimated fault slip rates (Authemayou et al., 2009) as a control parameter to determine the best model. Case#1gives the highest slip rates on the MRF and Kazerun strike–slip faults. On the contrary, case#5 indicates all the faults to be locked. Cases#2 and #3, with intermediate effective fault friction, show intermediate faults slip rates. The comparison between numerical experiments and the geological estimates of the fault slip rates suggest that case#1 is the optimal model. In this experiment, the effective fault friction is so low, frictional fault slip occurs in the lower crust and the upper mantle (Fig. 7). Also, I use the case#6 to check if we can obtain fault slip rates in agreement with the geological estimations without having a discontinuity of the velocity field in the mantle under the crustal faults. To do so, I apply an effective fault friction l = 0.02 in the crust and l = 0.3 for the extent of the fault down to the mantle. The results of case#6 show much lower velocities than case#1 (l = 0.02) on the whole fault plane, and are not in agreement with geological estimates. Therefore, a low fault friction is also needed in the mantle. This may imply that slip may occur in the mantle or that a localized strain zone develops in response to a low strength of the fault zone. The most appropriate model to describe the tectonic deformation is the numerical experiment #1, since it provides the best agreement with the geological fault slip rates (Table 4). The strain rate map and the velocity field obtained for this experiment #1 show quasi-rigid blocks (i.e., low strain rate areas) in the Central

96

H.R. Nankali / Journal of Asian Earth Sciences 41 (2011) 89–98

tures and mechanisms accommodating the internal deformation remain unclear, but they are certainly correlated with the fact that most of the upper crust in the Zagros consists of limestone. Therefore twinning and promoting crystal–plastic mechanisms are more likely to occur at low temperature, creating internal deformation rather than slip on the faults. It seems that based on these results, the appropriate description of deformation in Zagros seems to be a compromise of high strain localization on faults, quasi-rigid blocks and continuum deformation.

10.1. Comparison with geological displacement rate

Fig. 7. Slip rate estimated from model, with fault friction (l = 0.02).

Iran (CIB). The high strain rate zones are well correlated with the seismicity (Fig. 8). A main weakness of this modeling is the lack of dipping faults to accommodate shortening. Slip on basement faults does not reach the Earth’s surface as it is decoupled from folding within sedimentary layers by thick evaporite deposits that occur at several levels within the Phanerozoic section, most importantly at the base, in the infra-Cambrian Hormuz Salt Formation (e.g. Berberian and King, 1981). However, if dipping faults could be added to the model, the conclusions would remain the same. Indeed, in the case#1, the deformation could be fully accommodated by a strike–slip fault, but some is still taken up by the continuum medium. This is the case for example of the Main Recent Fault (MRF) which accommodates only a major part of the Zagros range parallel component. Moreover deformation in Iran is known to be partly aseismic (Masson, 2006), the best example being the Zagros where less than 5% of the whole deformation seems to be seismic. The struc-

The identification of total horizontal fault displacements and dating of offset alluvial fans by (Authemayou et al., 2009) enable us to compare the model results with long term average slip rates. For the Dena and the Kazerun Fault segments, geological displacement rates are available and evaluated to 2.5–4 mm/yr and 1.5– 3.5 mm/yr, respectively. The comparison of these long-term slip rates with the model slip rate shows that the geological and the model results are consistence (within the uncertainties) and less than the 14 mm/yr suggested by Berberian (1995). Assuming a constant slip rate since the onset of the faults and considering the total fault offsets of 13 and 8 km on the Dena and the Kazerun Fault segments, respectively, the coherent long-term slip rates yield comparable onset ages for Dena 2.7 Ma by model with respect to 3.2–5.2 Ma by geology and Kazerun 3.9 Ma by model with respect to 2.3–5.4 Ma by geology. We do not have a geological onset time for the Borazjan fault, but the model results shows that 0.3–1 mm/yr slip rate for this fault as consistence by GPS (Tavakoli et al., 2008). Long-term Main Recent Fault slip rates derived from geological observations in the Zagros area were compiled by several researches (e.g. Berberian and King, 1981; Berberian, 1995; Talebian and Jackson, 2002; Authemayou et al., 2009). Instantaneous slip rate along the MRF deduced from GPS measurements (Vernant et al., 2004) would be of 5 ± 1.5 mm/yr if the fault is considered to achieve total partitioning of the 7 ± 2mm/yr of N–S shortening measured around 50°E. Based on local GPS network surveys, Walpersdorf et al. (2006) suggest 4–6 mm/yr of cumulated right-lateral strike-slip across the whole eastern Zagros. The range of published estimates for the slip rate along the MRF

Fig. 8. Strain rate map (s1).

H.R. Nankali / Journal of Asian Earth Sciences 41 (2011) 89–98 Table 4 Total horizontal fault slips and age of strike-slip onsets from Authemayou et al. (2009), comparison with the model results case#1 (l = 0.02). Fault

Total slip (km)

Age of fault onset (Ma)

Average geological slip rate (Authemayou, 2009) (mm/yr)

Finite element model (this study) (mm/ yr)

Model (this study) fault onset (Ma)

Dena

13

2.5–4

4.8

2.7

Kazerun

8.2

1.5–3.5

2.1

3.9

Borazjan MRF

0 50

3.2– 5.2 2.3– 5.4 – 4–14.2

– 2.1–5.9 mm/yr 1.5–6.5 mm/yr

0.3–1 2.3

– 21.7

therefore varies from 1.2 to 17 mm/yr, calling for refinement based on Quaternary, precisely dated offset along the MRF (Talebian and Jackson 2002; Vernant and Chery, 2006). I use the geologically estimated fault slip rate by (Authemayou et al., 2009), as a control parameter for the model results. For the MRF, geological displacement rates evaluated as 3.5–12.5 mm/yr (Authemayou et al., 2009). The slip rate estimated for the MRF is total slip rate, which add the Farsan (2.1–5.9 mm/yr) and Ardal (1.5–6.5 mm/yr) segments. The comparison of this long-term slip rate with the model slip rate (for this segment), 2.3 mm/yr shows that the geological and the model results are nearly consistence (close to minimum band) and less than the 10–17 mm/yr suggested by Talebian and Jackson, 2002, which may be due to poorly constrained age of slip initiation on the fault. 11. Conclusion I have developed a 3D visco-elastic finite element model for simulating lithospheric deformation with faults. In this study I adapted this model to simulate long-term, quasi-steady state slip rates of the MRF and Dena-Kazerun and Borazjon faults. For simplicity purpose, I was neglecting two other less important Sabzpushan and Sarvestan faults in the model. The thermal, rheological and structural models investigated in this case are only a scratch on the surface of tectonic modeling. However, they provide important information in describing possible areas of deformation, and also for developing more sophisticated and complex models. These results are consistent with southward distribution of the slip from along the Main Recent Fault to the longitudinal thrusts and folds of the fold-and-thrust belt through the Kazerun Fault and associated faults, with a decrease of slip from the southeastern tip of the Main Recent Fault towards the southeastern termination of the Kazerun Fault. Combined with GPS measurement results (Vernant et al., 2004), our results suggest complete strike-slip partitioning of plate convergence in the Zagros on the Main Recent Fault. Loaded with the present day far field velocity the model for Zagros predicts maximum dextral slip rates of the Main Recent Faults and the Kazerun Faults on their central segments of about 2.3 mm/yr and 2.1 mm/yr, respectively. This suggest that geological slip rate of these faults are compatible with the far field velocity field of the Arabia -Eurasia collision zone giving our current understanding of the rheology of the Zagros. Fault slip rate estimates from this model leads to a very low fault friction in order to fit geological slip rate values. This finding is in good agreement with the low fault friction proposed for some other intracontinental strike–slip faults (e.g., NAF and San Andreas fault). Acknowledgments I wish to thank anonymous reviewers for their useful comments and suggestions. Ansys is a trademark of ANSYS Inc.

97

References Artemieva, I.M., Mooney, W., 2001. Thermal thickness and evolution of Precambrian lithosphere, a global study. JGR 106 (B8), 16387–16414. Authemayou, C., Bellier, O., Chardon, D., Benedetti, L., Malekzade, Z., Claude, C., Angeletti, B., Shabanian, E., Abbassi, M.R., 2009. Quarternary slip rates of the Kazerun and Main recent faults: active strike-slip partitioning in the Zagros fold and thrust belt. Geophys. J. Int.. doi:10.1111/j.1365246X.2009.04191. Bachmanov, D.M., 2004. Active faults in the Zagros and central Iran. Tectonophysics 380 (3–4), 221–241. Baker, C., Jackson, J., Priestley, K., 1993. Earthquakes on the Kazerun Line in the Zagros mountains of Iran: strike-slip faulting within a fold and-thrust belt. Geophys. J. Int. 115, 41–61. Berberian, M., 1995. Master blind thrust faults hidden under the Zagros folds: active basement tectonics and surface morphotectonics. Tectonophysics 241, 193– 224. Berberian, M., King., 1981. Active Faulting and Tectonics of Iran, in Zagros-HinduKush-Himalaya Geodynamic Evolution, vol. 3, pp. 33–69. Bird, P., 1978. Finite element modeling of lithosphere deformation: the Zagros collision orogeny. Tectonophysics 50, 307–336. Bird, P., Toksoz, M.N., Sleep, N.H., 1975. Thermal and mechanical models of continent–continent convergence zones. J. Geophys. Res. 32, 4405–4416. Carter, N.L., Tsenn, M.C., 1987. Flow properties of continental lithosphere. Tectonophysics 136, 27–63. Dehghani, G.A., Makris, J., 1984. The gravity field and crustal structure of Iran. N. Jb. Geol. Paläont. Abh. 168, 215–229. Engdahl, E.R., Jackson, J., Myers, S.C., Bergman, E.A., Priestley, K., 2006. Relocation and assessment of seismicity in the Iran region. Geophys. J. Int. 167, 761–778. England, E., McKenzie, D., 1983. Thin viscous sheet model for continental deformation. Geophys. J. R. Astron. Soc. 73, 523–532. Fernàndez, M., Ranalli, G., 1997. The role of rheology in extensional basin formation modelling. Tectonophysics 282, 129–145. Hatzfeld, D., Tatar, M., Priestley, K., Ghafory-Ashtyany, M., 2003. Seismological constraints on the crustal structure beneath the Zagros mountain belt (Iran). Geophys. J. Int. 155, 403–410. Hessami, K., Nilforoushan, F., Talbot, C.J., 2006. Active deformation within the Zagros mountains deduced from GPS measurements. J. Geol. Soc. Lond. 163, 143–148. Jackson, J., 1995. The accommodation of the Arabia–Eurasia plate convergence in Iran. J. Geophys. Res. 100, 15205–15219. Kirby, S., 1983. Rheology of the lithosphere. Rev. Geophys. Space Phys. 21, 1458– 1487. Liu, Z., Bird, P., 2007. Kinematic modelling of neotectonics in the Persia–Tibet– Burma orogen. Geophys. J. Int.. doi:10.1111/j.1365-246X.2007.03640. Masson, F., 2006. Seismic versus aseismic deformation in Iran inferred from earthquakes and geodetic data. Geophys. J. Int. 160 (1), 217–226. Paul, A., Kaviani, A., Hatzfeld, D., Vergne, J., Mokhtari, M., 2006. Seismological evidence for crustal-scale thrusting in the Zagros mountain belt (Iran). Geophys. J. Int.. doi:10.1111/j.1365-246X.2006.02920. Reilinger, R., McClusky, S., Vernant, P., Laurence, S., Ergintav, S., Cakmak, R., Ozener, H., Kadirov, F., Guliev, I., Stepanyan, R., Nadariya, M., Hahubia, G., Mahmoud, S., ArRajehi, A., Abdulaziz, K., Paradissis, D., Al-Aydrus, A., Prilepin, M., Guseva, T., Evren, E., Dmitrotsa, A., Filikov, S.V., Gomez, F., Al-Ghazzi, R., Karam, G., 2006. GPS constraints on continental deformation in the Africa–Arabia–Eurasia continental collision zone and implications for the dynamics of plate interactions. J. Geophys. Res. 111 (B5). doi:10.1029/2005JB004051. Savage, J., Burford, F., 1973. Geodetic determination of relative plate motion in Central California. J. Geophys. Res. 95, 4873–4879. Sella, G.F., Dix, T.H., Mao, A., 2002. REVEL: a model for recent plate velocities from space geodesy. J. Geophys. Res. 107 (B4) (ETG 11-1, 11-32). Snyder, D.B., Barazangi, M., 1986. Deep crustal structure and flexure of the Arabian plate beneath the Zagros collisional mountain belt as inferred from gravity observations. Tectonics 5, 361–373. Sobouti, F., Arkani-Hamed, J., 1996. Numerical modelling of the deformation of the Iranian plateau. Geophys. J. Int. 126, 805–818. Stocklin, J., 1968. Structural history and tectonics of Iran: a review. Am. Associat. Petrol. Geolog. Bull. 52, 1229–1258. Talbot, C.J., Alavi, M., 1996. The past of a future syntaxis across the Zagros. Geol. Soc. Spec. Publ. Lond. 100, 89–110. Talebian, M., Jackson, J., 2002. Offset on the main recent fault of the NW Iran and implications on the late Cenozoic tectonics of the Arabia–Eurasia collision zone. Geophys. J. Int. 150, 422–439. Talebian, M., Jackson, J., 2004. A reappraisal of earthquake focal mechanisms and active shortening in the Zagros mountains of Iran. Geophys. J. Int. 156, 506–526. Tatar, M., Hatzfeld, D., Ghafory-Ashtiany, M., 2004. Tectonics of the Central Zagros (Iran) deduced from microearthquake seismicity. Geophys. J. Int. 156, 255–266. Tavakoli, F., Walpersdorf, A., Authemayou, C., Nankali, H.R., Hatzfeld, D., Tatar, M., Djamour, Y., Nilforoushan, F., Cotte, N., 2008. Distribution of the right-lateral strike‑‘slip motion from the Main Recent Fault to the Kazerun Fault System Zagros, Iran): Evidence from present-day GPS velocities. Earth Planet. Sci. Lett. 275 (3–4), 342–347. Tchalenko, J.S., Braud, J., 1974. Seismicity and structure of the Zagros (Iran)—the Main Recent Fault between 33 and 35.N. Phil. Trans. Roy. Soc. Lond. 277, 1–25. Thatcher, W., 2007. Microplate model for the present-day deformation of Tibet. J. Geophys. Res. 112, B01401. doi:10.1029/2005JB004244.

98

H.R. Nankali / Journal of Asian Earth Sciences 41 (2011) 89–98

Vernant, P., Chery, J., 2006. Mechanical modelling of oblique convergence in the Zagros, Iran. Geophys. J. Int. 165, 991–1002. Vernant, P., Nilforoushan, F., Hatzfeld, D., Abbassi, M., Vigny, C., Masson, F., Nankali, H., Martinod, J., Ashtiani, A., Bayer, R., Tavakoli, F., Chéry, J., 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., Hatzfeld, D., Nankali, H.R., Tavakoli, F., Nilforoushan, F., Tatar, M., Vernant, P., Ch´ery, J., Masson, F., 2006. Difference in the GPS deformation pattern of North and Central Zagros (Iran). Geophys. J. Int. 167, 1077–1088. Yamini-Fard, F., Hatzfeld, D., Tatar, M., 2007. Microearthquake seismicity at the intersection between the Kazerun fault and the Main Recent Fault (Zagros, Iran). Geophys. J. Int. 166, 186–196. doi:10.1111/j.1365-246X.2006.02891.x.