Isostatic GOCE Moho model for Iran

Journal of Asian Earth Sciences 138 (2017) 12–24

Contents lists available at ScienceDirect

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

Full length Article

Isostatic GOCE Moho model for Iran Mehdi Eshagh a,⇑, Sahar Ebadi b, Robert Tenzer c a

Department of Engineering Science, University West, Sweden Department of Surveying Engineering, University of Tehran, Iran c New Technologies for the Information Society (NTIS), University of West Bohemia, Czech Republic b

a r t i c l e

i n f o

Article history: Received 7 October 2016 Received in revised form 23 January 2017 Accepted 26 January 2017 Available online 31 January 2017 Keywords: Crust Integral inversion Moho Satellite gradiometry Isostasy

a b s t r a c t One of the major issues associated with a regional Moho recovery from the gravity or gravity-gradient data is the optimal choice of the mean compensation depth (i.e., the mean Moho depth) for a certain area of study, typically for orogens characterised by large Moho depth variations. In case of selecting a small value of the mean compensation depth, the pattern of deep Moho structure might not be reproduced realistically. Moreover, the definition of the mean compensation depth in existing isostatic models affects only low-degrees of the Moho spectrum. To overcome this problem, in this study we reformulate the Sjöberg and Jeffrey’s methods of solving the Vening-Meinesz isostatic problem so that the mean compensation depth contributes to the whole Moho spectrum. Both solutions are then defined for the vertical gravity gradient, allowing estimating the Moho depth from the GOCE satellite gravity-gradiometry data. Moreover, gravimetric solutions provide realistic results only when a priori information on the crust and upper mantle structure is known (usually from seismic surveys) with a relatively good accuracy. To investigate this aspect, we formulate our gravimetric solutions for a variable Moho density contrast to account for variable density of the uppermost mantle below the Moho interface, while taking into consideration also density variations within the sediments and consolidated crust down to the Moho interface. The developed theoretical models are applied to estimate the Moho depth from GOCE data at the regional study area of the Iranian tectonic block, including also parts of surrounding tectonic features. Our results indicate that the regional Moho depth differences between Sjöberg and Jeffrey’s solutions, reaching up to about 3 km, are caused by a smoothing effect of Sjöberg’s method. The validation of our results further shows a relatively good agreement with regional seismic studies over most of the continental crust, but large discrepancies are detected under the Oman Sea and the Makran subduction zone. We explain these discrepancies by a low quality of seismic data offshore. Ó 2017 Published by Elsevier Ltd.

1. Introduction The three gravity-dedicated satellite missions, namely the CHAllenging Mini-satellite Payload (CHAMP), the GRavity field and Climate Experiment (GRACE) and the Gravity field and steady-state Ocean Circulation Explorer (GOCE), significantly improved our knowledge about the external gravitational field of the Earth at the long-to-medium wavelengths. The gravitygradiometry data from the latest satellite mission GOCE (ESA 2012) provide the information about the gravitational field with a spectral resolution approximately up to the spherical harmonic degree of about 280–300. Since the gravity gradients have a more pronounced regional support than the gravity, many of recent studies have focused on a regional gravity modelling and interpretations using gravity gradients. Moreover, the observed gravity ⇑ Corresponding author. E-mail address: [email protected] (M. Eshagh). http://dx.doi.org/10.1016/j.jseaes.2017.01.033 1367-9120/Ó 2017 Published by Elsevier Ltd.

gradients at satellite level are less-processed, giving more options for applying numerical schemes optimised for a particular purpose of the study. Among these studies we could mention work by Reed (1973), Xu (1992, 1998, 2009), Janák et al. (2009), Eshagh (2009a), and Eshagh and Sjöberg (2011). In studies of the lithospheric structure, GOCE data have been used extensively to recover the Moho interface. The Moho interface, otherwise known as the Mohorovicˇic´ discontinuity (e.g., Fowler, 1990), represents the boundary between the lowermost crust and the underlying uppermost mantle. Either seismic or gravity data are used to detect this interface. The seismic Moho surface is the interface at which a seismic wave velocity jump occurs, and the gravimetric one is the surface obtained from the gravity data inversion. Several hypotheses were proposed to explain the isostatic compensation mechanism, for instance, by Pratt (1854), Airy (1855) and Vening Meinesz (1931). Heiskanen (1931) presented a method to estimate the crustal thickness by considering a regional, instead of a local compensation scheme of topographic

13

M. Eshagh et al. / Journal of Asian Earth Sciences 138 (2017) 12–24

mass surplus and bathymetric mass deficiency. Parker (1972) modified the Vening Meinesz theory and presented an iterative method for finding the Moho depth, and Oldenburg (1974) stabilised this method by adding a low-pass filtering technique. Gomez-Oritz and Agarwal (2005) and Shin et al. (2007) generalised the Parker-Oldenburg method for a 3-D gravity inversion, and Kiamehr and Gomes-Ortiz (2009) applied this method to estimate the Moho depth in Iran by combining the terrestrial and satellite gravity data. Sünkel (1985) converted the Airy-Heiskanen Moho depth to the Vening Meinesz one based on minimising the global mean-square-error of differences between the disturbing and topographic-isostatic potentials. Moritz (1990) improved the Vening Meinesz hypothesis by adopting the spherical approximation of the Earth. Nevertheless, his method was not suitable for the gravity data inversion, because some theoretical deficiencies due to applied approximations. Sjöberg (2009) reformulated this problem by solving the Fredholm’s integral equation of the first kind. He called this method as the Vening Meinesz-Moritz (VMM) inverse problem of isostasy. Braitenberg et al. (2000) presented an iterative inversion method to obtain the Moho depth under the Tibet Plateau. Later, Braitenberg et al. (2006) presented a crustal model of the South China Sea obtained from combined processing of satellite-gravity, bathymetric, sediment, crustal thicknesses and isostatic flexure models. Tenzer et al. (2012b) studied the estimation of the Moho density contrast using the EGM2008 global gravitational model (Pavlis et al., 2008) and the CRUST2.0 global seismic crustal model (Bassin et al., 2000). Eshagh (2009b, 2010) investigated the effect of lateral density variations within the crustal and topographic masses on the GOCE data, and Eshagh et al. (2016) utilised these data for a gravimetric recovery of the Moho density contrast in central Eurasia. The crustal density-contrast stripping corrections were applied systematically to the topographically-corrected gravity field by Tenzer et al. (2009, 2012a, 2012c, 2015a, 2015b); see also Vajda et al. (2008) and Rexer et al. (2016). Braitenberg and Ebbing (2009) combined the GRACE and terrestrial gravity data to study the crustal structure. In a simulation study, Sampietro (2009) investigated a possibility of using the GOCE gravity-gradiometry data to recover the Moho depth globally. Sampietro (2011) considered a local inversion of the satellite gravity-gradiometry data by simulating a Moho surface and generating the data from it. Sampietro and Reguzzoni (2011) applied the collocation and Fast Fourier Transform (FFT) techniques for a Moho recovery from GOCE data, and Reguzzoni and Sampietro (2012) presented a global crustal model based GOCE data. Later, Reguzzoni et al. (2013) developed a new global Moho model based on combining seismic and GOCE data. Barzaghi et al. (2013) presented a collocation method for combining the global gravity-gradient information from GOCE with terrestrial gravity data. Sampietro et al. (2014) estimated the Moho depth beneath the Tibet Plateau and Himalaya. Reguzzoni and Sampietro (2014) applied stripping corrections due to sediments and consolidated crust in deriving the Moho model from GOCE data. Eshagh (2014a) presented a simple, linear method for a direct Moho recovery from the satellite gravity-gradiometry data by solving the VMM problem. Eshagh and Hussain (2015) investigated the relation between the geoid, gravity anomaly, deflection of the vertical and the Moho geometry. Tenzer and Chen (2014a, 2014b) derived the spectral and spatial expressions of solving the gravimetric Moho problem from gravity and crustal structure models. Eshagh et al. (2016) used GOCE data to determine the Moho density contrast in central Asia, and Eshagh and Hussain (2016) used these data to determine the Moho depth for the Indo-Pakistan region. In theoretical study, Eshagh (2016a,b) demonstrated that the VMM theory is a generalisation of the Airy-Heiskanen isostatic model and compared it to the flexural isostasy models.

Numerous studies have been conducted to determine a regional Moho model for Iran from seismic data. However, seismic data over the country are distributed only sparsely and irregularly. Moreover, each study was mostly restricted to a specific area, like profiles between Shiraz-Mashhand, Tehran-Mashhad and Mashhad-Tabriz (Asudeh, 1982), the southern part of the Caspian Sea (Mangino and Priestley, 1998), Tehran region (Hatzfeld et al., 2003), Mashhad (Javan Doloei and Roberts, 2003; Javan Doloei, 2003), Central Alborz and the northern Iran (Sodoudi et al., 2009; Radjaee et al., 2010), central Zagros (Paul et al., 2006; Manaman et al., 2011), Kopeh Dagh (Nowrouzi et al., 2007), Naein (Nasrabadi et al., 2008), the northwest Iran (TaghizadehFarahmand et al., 2010), and Sanandaj-Sirjan zone (Sadidkhouy et al., 2012), while seismic data over whole Iran were used so far only by Taghizadeh-Farahmand and Afsari (2015). The terrestrial gravity data in Iran have, on the other hand, significantly better coverage and resolution than seismic data. The combined gravimetric and seismic data for Iran were used, for instance, by Dehghani and Makris (1984). Snyder and Barazangi (1986) determined the Moho depth under the Persian Gulf and Zagros fold and thrust belt. Jiménez-Munt et al. (2012) calculated the crustal and lithospheric mantle structure from the geoidal undulations and thermal analysis. Safari et al. (2014) estimated the Moho depth in Iran using the Euler deconvolution method, and validated their gravimetric result at 91 seismic stations. Eshagh (2014a) performed a simulation study to compare a numerical performance of two methods for a Moho recovery from GOCE data with a case study for Iran, and Eshagh (2014b) estimated the Moho depth for Iran from terrestrial gravity data. The aforementioned gravimetric methods, however, have some theoretical and numerical deficiencies related to a definition of the mean compensation attraction as well as the treatment of variable crustal density within the whole crust and the uppermost mantle. To address these issues, we reformulate Jeffrey and Sjöberg’s methods for a Moho determination from GOCE data so that the mean compensation depth is considered over the whole Moho spectrum. Moreover, we also take into consideration the variable Moho density contrast, and further incorporate the effects of sediments and consolidated crust into the solution. Developed theoretical and numerical models are then applied to determine a regional Moho model for Iran from (in-orbit) GOCE data. 2. Theory After reviewing Jeffrey and Sjöberg’s methods for solving the VMM isostatic problem, their modifications (as mentioned above) are presented in this section. 2.1. Jeffrey’s solution Let us begin with a definition of the isostatic gravity disturbance dgI. According to Eshagh and Hussain (2016), we write

dg I ¼ dg  dg TB  dg S  dg Crys þ dg C ;

ð1aÞ TB

where dg is the gravity disturbance, dg is the topographic/bathymetric (TB) effect on dg, dg C is the compensation attraction, dg S is the gravitational effect of sediments, and dg Crys is the gravitational effect of the consolidated crust. The formula in Eq. (1a) represents the fundamental condition used by Moritz (1990) and Sjöberg (2009) to solve the VMM problem, while assuming that dgI equals zero. Following their theoretical concept, we define the gravitational contribution dV (at any point outside the Earth’s masses) generated by an infinitesimal volumetric mass element in the following form

14

M. Eshagh et al. / Journal of Asian Earth Sciences 138 (2017) 12–24

dV ¼ Gq

r 02 dr 0 dr ; l

ð1bÞ

where G is the Newton’s gravitational constant, q is the mass density distribution function, r0 is the geocentric distance of the infinitesimal mass element which is defined by the radial and surface components dr0 and dr ¼ sin hdhdk respectively, and l is the Euclidean spatial distance. The 3-D position is defined by the geocentric distance r, spherical co-latitude h and longitude k. In spectral domain, the reciprocal spatial distance in Eq. (1b) reads (e.g., Heiskanen and Moritz, 1967, p. 33) 1 1 X r 0n ¼ P ðcos wÞ; nþ1 n l r n¼0

ð1cÞ

where Pn ðcos wÞ is the Legendre polynomial of degree n for the argument of the geocentric angle w. We further define the compensation potential VC for solving the VMM problem according to Jeffrey (1976). Instead of dividing the compensation potential into two parts, as done by Moritz (1990, p. 255) and later by Sjöberg (2009), we define it in the following form

ZZ Z VC ¼ G

RT 0 RT 0 t

r

Dqr 02 0 dr dr l

Finally, we apply the summation over spherical harmonics on both sides of Eq. (2g) and solve the solution for t. This procedure yields

t¼

ð2aÞ

2.2. Jeffrey’s solution for GOCE data Here we extend Jeffrey’s definition (given in Section 2.1) for the vertical gravitational-gradient component Vrr. From Eshagh et al. (2016) and Eshagh and Hussain (2016), the relation between the Laplace harmonics for Vrr and dg reads

dg n ¼

VC ¼ G

ZZ X 1 Dq 1 ½ðR  T 0 Þnþ3 nþ1 n þ 3 r n¼0 r

 ðR  T 0  tÞnþ3 P n ðcos wÞdr:

ð2bÞ

r2 1  r nþ1 V rr;n : R nþ2 R

1  1 X 2n þ 1 ðnþ2Þ  TB j dg n þ dg Sn þ dg Crys  W; n 4pGDq n¼0 n þ 1

ZZ X 1 Dq nþ1 r n¼0 r

nþ3

ðR  T 0 Þ  nþ3

Wn ¼

 r nþ1 1 r2 2n þ 1 jðnþ2Þ V rr;n : 4pGDq R ðn þ 1Þðn þ 2Þ R

"

 1 1

r

where

1 X nþ1

Dq

dg C jr¼R ¼ GR

j¼1

t R  T0

nþ3 # Pn ðcos wÞdr:

For a regional inversion of GOCE data, the solution of W is conveniently converted into the spatial form by means of solving integral equations. First, we simplify Eq. (3c) and solve it for rV rr;n . This procedure yields

rV rr;n ¼ 4pGDq

Wn ¼ ð2cÞ

n¼0

nþ3

"

nþ3

j

 1 1

t R  T0

nþ3 # P n ðcos wÞdr

ð2dÞ

1 GR X ¼ ðn þ 1Þjnþ3 R  T 0 n¼0

ZZ ðDqtÞPn ðcos wÞdr:

ð2eÞ

r

The definition of the surface integral convolution for Dqt in terms of the Laplace harmonics yields

dg C;n

n þ 1 nþ2 ¼ 4pG j ðDqtÞn : 2n þ 1

ð2fÞ

From Eqs. (1c), (2e) and (2f), we get

ðDqtÞn ¼ 

 ðnþ2Þ ðn þ 1Þðn þ 2Þ R  T 0 Wn: 2n þ 1 r

ð3dÞ

2n þ 1 4p

ZZ

W 0 Pn ðcos wÞdr:

ð3eÞ

r

Substituting from Eq. (3e) to Eq. (3d) and applying the summation over the Laplace harmonics, we arrive at

ZZ

rV rr ¼ GDq

Kðr; wÞW 0 dr;

ð3fÞ

r

with the kernel function K given by

T0 : R

Approximating the term [t/(R  T0)]n+3 by a binomial series (up to the linear term) and further simplifications, we arrive at

dg C;n

ð3cÞ

The Laplace harmonics W n in Eq. (3d) are defined by the following integral convolution (Heiskanen and Moritz, 1967, p. 33)

With reference to Heiskanen and Moritz (1967, p. 85), the harmonic expansion of dgC is given by

ZZ

ð3bÞ

where the parameter W represents the contribution of vertical gravity gradient to the Moho geometry. In spectral domain, it reads

We further rewrite Eq. (2b) into the following form

VC ¼ G

ð3aÞ

Inserting from Eq. (3a) to Eq. (2h) and some further simplifications, we arrive at

r

where t is the Moho depth correction (relative to the mean Moho depth T0). Solving the radial integral and subsequently substituting for integral limits, we get

ð2hÞ

For the estimated values of t, the Moho depth is found as follows: T ¼ T 0 þ t. As seen from these definitions, Jeffery’s method determines a Moho flexure from gravity data, not the total Moho depth.

t¼

Z ZZ X 1 Dq RT 0 0nþ2 0 ¼G r dr Pn ðcos wÞdr; rnþ1 RT 0 t n¼0

1  1 X 2n þ 1 ðnþ2Þ  S Crys j dg n  dg TB : n  dg n  dg n 4pGDq n¼0 n þ 1

 1 2n þ 1 ðnþ2Þ  S Crys j dg n  dg TB : n  dg n  dg n 4p G n þ 1

ð2gÞ

Kðr; wÞ ¼

 ðnþ2Þ 1 X R  T0 ðn þ 1Þðn þ 2Þ Pn ðcos wÞ: r n¼2

ð3gÞ

The integral kernel K in Eq. (3g) is a function of T 0 , meaning that the behaviour of this kernel is dependent on the value of the mean Moho depth chosen for a particular study area. The closed-form of the kernel K reads 2

Kðr; wÞ ¼ s

2

2 sðs þ 4kÞ 3s2 k þ 5  3 l l l

!

;

ð3hÞ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where l ¼ s2  2s cos w þ 1, and k ¼ s  cos w. The behaviour of the integral kernel K (in Eq. (3h)) for the mean Moho depth values: T0 = 10, 20 and 30 km is illustrated in Fig. 1. As seen, the maximum signal is in vicinity of the computation point. The signal further attenuates uniformly with the increasing spher-

15

M. Eshagh et al. / Journal of Asian Earth Sciences 138 (2017) 12–24

ical distance (i.e., the geocentric angle), and reaches approximately zero value at 5°. The truncation error beyond this spherical distance is almost negligible (cf. Eshagh, 2011). Moreover, the signal of K slightly magnifies with the decreasing mean Moho depth T0. The signature of a more detailed pattern of the Moho geometry is thus to some extent filtered out with an increasing value of the mean Moho depth. In practice, however, the Moho pattern under a thin oceanic crust is typically very smooth, so that the choice of a larger value of the mean Moho depth for a more realistic modelling of the Moho geometry under the thick continental crust will not affect significantly the Moho modelling under the oceanic crust. In contrast, the choice of a small value of the mean Moho depth will affect the accuracy of Moho modelling under the continental crust, and particularly under orogens with the largest Moho deepening.

Assuming that T = T0, the expression in Eq. (4c) is further rewritten as

dg C ¼ GRDq

  ZZ X 1 T T 1  ðn þ 2Þ Pn ðcos wÞdr ðn þ 1Þ R 2R n¼0 r

ZZ X 1 ðn þ 1Þ  1  jnþ3 Pn ðcos wÞdr  GRDq ðn þ 3Þ n¼0 r

The corresponding spectral form of Eq. (4d) reads

dg C;n ¼ 4pGDqbn

nþ1 4pGRDq Tn  ð1  jnþ3 Þdn0 ; 2n þ 1 3

S Crys dg n  dg TB þ 4pGDqbn n  dg n  dg n

V C ¼ GDq

Z R

ZZ X Z R 1 1 0nþ2 0 0nþ2 0 dr P ðcos wÞ r dr  r dr n r nþ1 RT RT 0 n¼0 r

¼ dV C þ V C0

ð4aÞ

The compensation potential in Eq. (4a) is divided into the residual and mean compensation potentials dV C and V C0 respectively. Solving the radial integrals in Eq. (4b) and simplifying the result of integration, we get

"  nþ3 # ZZ X 1 Rnþ3 1 T nþ3 Pn ðcos wÞdr: V C ¼ GDq j  1  R r nþ1 n þ 3 n¼0

ð4bÞ

r

For the gravity disturbance dg, the compensation attraction becomes

dg C ¼ GRDq

"  nþ3 # ZZ X 1 nþ1 T Pn ðcos wÞdr 1 1 nþ3 R n¼0 r

 GRDq

ZZ X 1 nþ1 r

n¼0

nþ3

½1  jnþ3 Pn ðcos wÞdr

ð4eÞ

T0 where bn ¼ 1  ðn þ 2Þ 2R . Inserting from Eq. (4e)–(1a), we get

2.3. Sjöberg’s solution A principal difference between Jeffery and Sjöberg’s definitions is attributed to the way they treat the compensation potential. Moritz (1990, p. 255) divided the compensation potential VC into two parts, and the same principle was used by Sjöberg (2009) for solving the VMM problem. From Moritz (1990, p. 255), we have

ð4dÞ

nþ1 4pGRDq Tn  ½1  jnþ3 dn0 ¼ 0: 2n þ 1 3 ð4fÞ

The solution of Eq. (4f) for Tn yields

Tn ¼

 R 1 1 2n þ 1 1  S Crys b dg n  dg TB þ bn ð1  jnþ3 Þdn0 : n  dg n  dg n 4pGDq n þ 1 n 3 ð4gÞ

Applying the summation (from degree 0 to 1) on both sides of Eq. (4g), we arrive at

T ¼ AC0 þ

1  1 X 2n þ 1 1  TB b dg n þ dg Sn þ dg Crys  dg n ; n 4pGDq n¼0 n þ 1 n

ð4hÞ

where

AC0 ¼

  R R : ð1  j3 Þ 3 R  T0

ð4iÞ

The resulting expression is now defined in the form where T0 contributes to all frequencies of the Moho spectrum. 2.4. Sjöberg’s solution for GOCE data

ð4cÞ

In order to incorporate the GOCE data into Sjöberg’s solution, we separate the contribution of dg n from the solution for W. To do so, we first write

W¼

1 1 X 2n þ 1 1 b dg n : 4pGDq n¼0 n þ 1 n

ð5aÞ

Substituting from Eq. (4a)–(5a) and applying a binomial series (up to the second-order term), the spectrum of Eq. (5a) becomes

Wn ¼

1 2n þ 1 1 1  r nþ1 2 r V rr;n : b 4pGRDq n þ 1 n n þ 2 R

ð5bÞ

Now we solve Eq. (5b) for r2 V rr;n . Hence, we get

4pGRDq

 nþ1 nþ1 R W n ¼ r 2 V rr;n : bn ðn þ 2Þ 2n þ 1 r

ð5cÞ

By analogy with Eq. (4e), we write

GRDqðn þ 1Þðn þ 2Þbn

 nþ1 ZZ R W 0 Pn ðcos wÞdr ¼ r 2 V rr;n r

ð5dÞ

r

Applying the summation on both sides of Eq. (5d) (from degree 2 to 1), the following integral equation is found

ZZ GRDq

Sðr; wÞPn ðcos wÞW 0 dr ¼ r 2 V rr ;

r Fig. 1. Behaviour of the integral kernel K (in Eq. (3h)).

where the kernel S reads

ð5eÞ

16

M. Eshagh et al. / Journal of Asian Earth Sciences 138 (2017) 12–24

Sðr; wÞ ¼ ðn þ 1Þðn þ 2Þbn

 nþ1 R Pn ðcos wÞ: r

ð5fÞ

The integral kernel S is a function T0, but comparing to the kernel K defined in Eq. (4g), it is functionally more complex. To find its closed-form analytical solution, we first formally separate the kernel into two parts:

Sðr; wÞ ¼

 nþ1 1 X R ðn þ 1Þðn þ 2Þ Pn ðcos wÞ r n¼2  nþ1 1 T0 X R ðn þ 1Þðn þ 2Þ2 Pn ðcos wÞ;  r 2R n¼2

ð5gÞ

where the closed-form formula of the first constituent on the righthand side of Eq. (5g) is equivalent to the expression given in Eq. (3h), while replacing the argument R by R  T0 in the definition of the parameter s. The closed-form formula for the second term reads

 nþ1 1 X R ðn þ 1Þðn þ 2Þ2 Pn ðcos wÞ r n¼2 2

¼

3

4s 3s4 k 8s2 ðs þ kÞ s3 k ð24 þ 564skÞ 15s4 k þ 2  þ  : 5 7 3 l l l l l

ð5iÞ

The behaviour of the integral kernel S (in Eq. (5g)) is illustrated in Fig. 2. By comparison, the integral kernels K (Fig. 1) and S (Fig. 2) have a very similar behaviour.

However, such detailed information about the geological structure is not required for a Moho modelling from GOCE data, due to a limited resolution of GOCE. Consequently, we determine the Moho depth with a spectral resolution complete to the spherical harmonic degree of 180, which corresponds (by means of the halfwavelength) to about 110 km on the equator. This choice of the spatial/spectral resolution is also supported by finding of Turcotte and Schubert (2014, p. 251), who mentioned that the loads having smaller wavelengths than 100 km are not isostatically compensated. The regional map also illustrates a tectonic boundary between the Arabian block and the Eurasian plate (in red). As seen, most of Iran comprises a rough topography, except for depressions of Central Iran, Lut block, Jazmourian and Makran basins, and southern border with the Caspian Sea. Zagros ranges extend from the south to northwest of Iran, and separate the Arabian block from the rest of Eurasian tectonic plate. The area includes a low amount of the Paleozoic outcrops modified by an ongoing sedimentation (Ghorbani, 2013) and without magmatic and metamorphic events. In the northeast edge of Zagros, the Sanandaj-Sirjan zone accommodates magmatic and metamorphic rocks. Margins between the Caspian Sea and Central Iran (central Alborz) are composed of different sedimentary rocks. The Kopeh Dagh Mountains and basin, in the northeast of Iran, are composed mainly by extrusive igneous rocks belong to Paleogenic volcanic areas. The Makran and Jazmourian depressions, in the southeast of Iran, are separated by a long range of ophiolites extending from the west to east (cf. Ghorbani, 2013).

3. Numerical realisation Theoretical definitions (in Section 2) were applied in this section to determine a Moho depth at the regional study area comprising the Iran block. 3.1. Study area We specified the study area of Iran between the latitudes 25°N and 45°N and the longitudes 40°E to 65°E. This region is characterised by a rough topography with a complex geological structure. The topography (including the bathymetric depths offshore) of the study area, retrieved from SRTM dataset (Farr et al., 2007), is shown in Fig. 3. The geological classification of the study area is according to Stöcklin’s (1968). It is worth mentioning that more detailed geological classifications of the study area can be found in Nabavi (1976), Nogol-e-Sadat (1993) and Aghanabati (2004).

Fig. 2. Behaviour of the integral kernel S (in Eq. (5g)).

3.2. Data acquisition Inspecting theoretical definitions given in Eqs. (2h) and (4h), we could readily recognise that a Moho depth determination depends considerably on a choice of the Moho density contrast, because changes between the Moho depth and density contrast are related proportionally. Until recently, constant values 600 kg m3 or 480 kg m3 of the Moho density contrast were typically assumed for the continental and oceanic lithosphere respectively. Today, however, our knowledge about the Moho density contrast variations improved over some lithospheric settings. Moreover, the

Fig. 3. Topography and geological setting of the study area of Iran.

17

M. Eshagh et al. / Journal of Asian Earth Sciences 138 (2017) 12–24

CRUST1.0 seismic model (Laske et al., 2013) provides the information on the Moho density contrast globally, despite the real accuracy of this model for certain parts of the world might still be disputable. To utilise this information, we used the CRUST1.0 data to define the variable Moho density contrast by taking a weighted mean of the consolidated crustal density layers (with respects to their thickness) and subsequently subtracted the CRUST1.0 upper mantle density. The result of this procedure is shown in Fig. 4a. A regional map of the combined TB gravitational effect dg TB , computed using the solid topography coefficients up to degree/ order 180 generated from SRTM30 data (Farr et al., 2007), is shown in Fig. 4b. A constant density of 2670 kg m3 was adopted for the topographic masses, and 1000 kg m3 for the seawater density. It is worth mentioning here that a more accurate model of the seawater was developed and applied by Gladkikh and Tenzer (2011) and Tenzer et al. (2011, 2012c). The gravitational effect of sediments dg Sed (see Fig. 4c) was computed from the CRUST1.0 sediment data, and the CRUST1.0 consolidated crust data were used to compute the respective effect (see Fig. 4d). A statistical summary of these effects is given in Table 1. As seen the contribution of these effects on the input gravity data, and consequently on the Moho result is significant. 3.3. Moho models To validate results of a regional Moho inversion, we used three seismic Moho models, specifically CRUST1.0, Meier’s (Meier et al.,

2007) and M13_Eurasia (Stolk et al., 2013), as well as the gravimetric Moho model GEMMA computed from GOCE data. These models are shown in Fig. 5, and their statistical summary is given in Table 2. As seen in Fig. 5a, the CRUST1.0 Moho model is rather smooth with a maximum Moho deepening under the SanandajSirjan zone. The Moho depth minima are, under the Oman Sea, close to the Makan subduction zone. Some more localised features are also recognised, corresponding to the Kopeh Dagh sedimentary basin and the Lut and Jazmourian depressions. As seen in Fig. 4b, Meier’s Moho model is smooth, without capturing some more detailed Moho features (otherwise seen in Fig. 4a). The GEMMA GOCE gravimetric model in Fig. 4d, confirmed Moho deepening under the Sanandaj-Sirjan zone, including a more localised deepening under the Makran zone and the Lut block. The M13_Eurasia Moho model is shown in Fig. 4c. The Moho depth maxima are situated under the Sanandaj-Sirjan zone, Makran area and Lut block. However, the maximum Moho deepening to 66 km is likely unre-

Table 1 Statistics of the TB, sediment and consolidated crust corrections to gravity disturbances [mGal].

TB Sediments Crust Total

Max

Mean

Min

STD

227.3 183.5 117.4 195.1

35.0 66.0 13.5 230.6

371.1 3.9 169.4 486.8

96.7 34.5 42.4 123.2

Fig. 4. (a) CRUST1.0 density contrast [kg/m3]; (b) the TB gravitational effect [mGal]; (c) the gravitational effect of sediments [mGal]; and (d) the gravitational effect of consolidated crust [mGal].

18

M. Eshagh et al. / Journal of Asian Earth Sciences 138 (2017) 12–24

Fig. 5. Moho models of (a) CRUST1, (b) Meier, (c) M13_Euarisa, and (d) GEMMA [km].

Table 2 Statistics of solid topography, Moho depth (CRUST1.0, Meier, M13_Eurasia and GEMMA models) and CRUST1.0 Moho density contrast.

H [km] CRUST1 [km] Meier [km] M13_Eurasia [km] GEMMA [km] DC [kgm3]

Max

Mean

Min

STD

2.8 50.1 52.9 65.5 46.7 697.5

0.5 38.3 38.5 40.8 35.9 542.5

3.7 7.1 20.9 11.6 17.9 334.3

1.0 6.9 5.8 9.0 4.2 67.9

alistically large. In contrast, a large Moho deepening in the northwest Iran, including mountainous parts of Turkey, seen in CRUST1.0 and Meier’s models, are in this model absent. Based on a standard deviation of Moho depth variations (in Table 2), the GEMMA Moho surface is the smoothest, while the M13_Euarisa Moho surface has the largest variations. However, according to geologic setting and topography of the study area, the maximum regional Moho deepening is expected under the central Zagros Mountains and Sanandaj-Sirjan belt, therefore, locations of the maximum Moho deepening in the M13_Eurasia model (Fig. 5c) are likely unrealistic. Moreover, the maximum Moho depth of 66 km is also likely overestimated. Eshagh et al. (2016) demonstrated that this model is relatively good in central

Asia, but the same does not apply for Iran. The CRUST1.0 and Meier’s models show similar values of the maximum Moho deepening and the same mean values, but Meier’s model shows slightly smaller range of values. The GEMMA Moho depths under orogens are systematically underestimated, while overestimated under the oceanic crust and some parts of continental basins. Such pattern is explained by numerical techniques applied for a Moho recovery from GOCE data. Whereas the gravity-gradient corrections are first continued upward to the GOCE satellite altitude (250 km), the resulting gravity-gradient data are then continued downward in the process of gravity inversion. Since the inversion is ill-posed, a regularisation (or a frequency damping) is applied. As a result, this procedure removes some higher-frequency spectrum in inverted GOCE data, with consequent smoothing of the Moho geometry. To reduce this smoothing effect, Eshagh (2014a) proposed a method to continue downward only the GOCE data to a Moho level, then adding the corrections to recovered quantities in order to obtain the final solution. 3.4. GOCE data To reduce the spatial truncation error of the integral formula, we extended the study area by 5° on each side for the GOCE data inversion. The GOCE vertical gravity gradients were corrected for the GRS80 (Moritz, 2000) normal gravity-gradient component.

M. Eshagh et al. / Journal of Asian Earth Sciences 138 (2017) 12–24

The spherical harmonic coefficients of the GRS80 normal gravity field were inserted into the spherical harmonic expansion of GOCE data, and the computed gradients were then subtracted from the observed ones to obtain the second-order radial derivative of the disturbing potential. The GOCE vertical gravity gradients are shown in Fig. 6. Maxima are distributed over the Zagros and Sanandaj-Sirjan belts, and the corresponding minima over the Caspian Sea and the Persian Gulf. 3.5. GOCE data inversion Since the discretised integral equations are ill-conditioned, we applied Tikhonov’s (1963) regularisation with a quasi-optimal method (Hansen, 1998). The inversion was carried out over the data area (see Fig. 6) on 1°  1° grid according to both, Jeffery (Eq. (3f)) and Sjöberg’s methods (Eq. (5e)). Computed values of the W parameters are shown in Fig. 7. Both maps are very similar, showing maxima over mountains and minima in the Caspian Sea and the Persian Gulf, including some parts of the Arabian plate. The only apparent difference between these two solutions is a slightly smoother result from Sjöberg’s method.

19

We further validated our results using available models. The statistical summary of Moho depth differences is given in Table 4. Since both, Jeffrey and Sjöberg’s methods are sensitive to the choice of T0, we tested models for different values of T0. As seen in Table 4, change of T0 affects the mean and RMS of differences, but not the STD of differences. The RMS of differences indicates that our results better agree with GEMMA. Since we used the same data and similar gravimetric-isostatic principle, this finding is not surprising. As mentioned above, Sjöberg’s method provides smoother solution. This is also evident when comparing results obtained after applying all gravity-gradient corrections. The RMS of differences between our solutions and GEMMA is 5.3 km (for Jeffery’s method) and 4.6 km (for Sjöberg’s method). When applying only the TB correction, the RMS of differences improved to 4.7 km (for Jeffery’s method) and 4.4 km (for Sjöberg’s method), indicating large uncertainties within the CRUST1.0 sediment and consolidated crust data. Tis is also evident from comparing the gravimetric Moho solutions, especially when applying the Jeffery’s method which strongly amplifies the signal. A principal difference between the Jeffery and Sjöberg’s methods is the influence of the mean Moho depth T0 on the whole Moho

3.6. Moho recovery The mean Moho depth of 30 km was chosen for a regional Moho recovery. The computation was realised individually for the corrections of the combined TB effect and all corrections including sediments and consolidated crust, while applying Jeffrey and Sjöberg’s methods. Results are shown in Fig. 8, and statistics of results are summarised in Table 3. For all corrections (see Figs. 8a and b), the Moho results are very similar, but Sjöberg’s method provides smoother solution. The maximum Moho deepening is detected under the Zagros Mountains and the Sanandaj-Sirjan belts, with extension under the Alborz Mountains and Kopeh Dagh. A relatively shallow Moho dominates central Iran and Lut block. The minimum Moho depth is detected under the Oman Sea and the Makran subduction zone. The results after applying all corrections are shown in Figs. 8c and d. In this case, Moho variations are considerably larger than those obtained after applying only the TB correction. A minimum Moho depth of 1.2 km from Jeffrey’s method under the Makran subduction zone is unrealistically small. Similarly, a minimum value of 5.9 km there according to Sjöberg’s solution seems to be too low. The application of sediment and consolidated crust stripping corrections does not change significantly results over continents.

Fig. 6. GOCE vertical gravity gradients [E] over the data area.

Fig. 7. W function recovered from integral equations (a) Jeffrey and (b) Sjöberg principles [km].

20

M. Eshagh et al. / Journal of Asian Earth Sciences 138 (2017) 12–24

Fig. 8. Moho model computed by (a) Jeffrey’s method for all corrections, (b) Sjöberg’s method for all corrections, (c) Jeffrey’s method for TB correction, and (d) Sjöberg’s method for TB correction. Table 3 Statistics of Moho depth computed according to Jeffrey and Sjöberg’s methods [km].

Jeffrey

All corrections Only TB correction All corrections Only TB corrections

Sjöberg

Max

Mean

Min

STD

51.2 50.3 48.1 47.8

30.8 32.9 30.7 32.8

1.2 15.9 5.9 18.1

7.0 4.7 6.3 4.2

Table 4 Statistics of differences between Moho models [km].

All corrections

Jeffrey

Sjöberg

Only TB correction

Jeffrey

Sjöberg

T0

Max

Mean

Min

STD

RMS

T0

Max

Mean

Min

STD

RMS

CRUST1.0 Meier M13 GEMMA CRUST1.0 Meier M13 GEMMA

30 30 30 30 30 30 30 30

18.3 14.3 11.6 10.4 17.8 12.6 9.6 8.1

7.5 8.3 10.3 5.1 7.5 8.3 10.3 5.1

30.1 30.0 29.7 20.2 24.5 27.3 28.0 19.4

5.9 6.2 6.3 5.2 5.5 5.6 6.0 4.5

9.6 10.4 12.1 7.2 9.3 10.1 11.9 6.8

38.3 39.3 41.3 36.1 38.3 39.3 41.3 36.1

26.6 22.6 22.9 16.5 26.1 21.9 20.9 14.2

0.8 0.0 1.0 1.0 0.7 1.0 1.0 1.0

21.9 21.7 18.4 14.1 16.3 17.9 16.7 13.3

5.9 6.2 6.3 5.2 5.5 5.6 6.0 4.5

6.0 6.2 6.4 5.3 5.5 5.7 6.1 4.6

CRUST1.0 Meier M13 GEMMA CRUST1.0 Meier M13 GEMMA

30 30 30 30 30 30 30 30

16.3 14.2 13.1 9.1 16.5 11.8 12.6 7.7

5.4 6.2 8.2 3.0 5.4 6.3 8.3 3.1

19.4 22.6 25.3 13.9 16.4 20.7 24.8 12.6

5.6 5.2 6.4 3.6 5.4 4.8 6.3 3.2

7.8 8.1 10.4 4.7 7.7 7.9 10.4 4.5

38.3 39.3 41.3 36.1 38.3 39.3 41.3 36.1

24.5 23.5 24.4 15.2 24.8 21.2 23.9 13.8

2.9 3.1 3.1 3.1 2.8 3.0 3.0 3.0

11.1 13.3 14.0 7.8 8.1 11.3 13.5 6.5

5.6 5.2 6.4 3.6 5.4 4.8 6.3 3.2

6.3 6.0 7.1 4.7 6.1 5.7 7.0 4.4

21

M. Eshagh et al. / Journal of Asian Earth Sciences 138 (2017) 12–24 Table 5 Statistics of the differences between Moho models computed by Jeffrey and Sjöberg’s methods [km].

All corrections Only TB correction

Max

Mean

Min

STD

RMS

3.6 3.2

0.0 0.0

5.6 3.5

0.8 0.6

0.8 0.6

Fig. 9. Moho depth differences [km] between results obtained by applying Jeffery and Sjöberg’s methods, for (a) all corrections, and only (b) TB correction.

spectrum. Jeffrey’s method amplifies more the signal at higher degrees than Sjöberg’s method. Eshagh (2016a) mentioned that the difference between the signal spectra globally reach up to about 3 km. In case of the regional study area of Iran, the differences between Moho solutions from applying Jeffery and Sjöberg’s methods are seen in Table 5 (see also Fig. 9). For all corrections, the RMS of differences is 0.8 km and reduces to 0.6 km when applying only the combined TB correction. Locally, however, these differences could reach large values; up to 5.6 km (in absolute sense) when applying all corrections, and 3.5 km when applying only the combined TB correction. 4. Discussion We utilised two methods for a Moho recovery from GOCE data, while applying the combined TB correction and all corrections (including sediments and consolidated crust). Generally, the Moho

solution computed according to Jeffery’s method is too shallow over the Oman Sea and Makran subduction zone. It reaches there only 1.4 km, which is unrealistic. Sjöberg’s method, on the other hand, provides smoother solution, with much less pronounced maxima and minima, for instance, the maximum Moho deepening is less than 50 km. In case of ignoring the effect of sediments and consolidated crust, results from Jeffrey and Sjöberg’s methods are more similar, with maximum differences between them up to 3.2 km. We further compared our gravimetric results with existing regional gravimetric and seismic studies for Iran. Safari et al. (2014) estimated the crustal thickness in Iran using the Euler deconvolution method. Their results revealed that the best RMS fit between the gravimetric and seismic models (at 91 seismic stations) was attained for the mean Moho depth of 45 km. According to Taghizadeh-Farahmand and Afsari (2015), however, the mean Moho depth is less than 40 km. Since we applied regional gravimetric methods for a Moho recovery, the optimal choice of the mean Moho depth is an essential parameter for the final result. According to our results, the Moho depth maxima of about 45 km are situated in central Zagros and the southern Iran. This finding fully agrees with results of Dehghani and Makris (1984) and Hatzfeld et al. (2003) who mentioned values around 46 km. A depth of 41 km was proposed by Taghizadeh-Farahmand and Afsari (2015). Sadidkhouy et al. (2012), by using receiver functions, concluded that the crustal thickness of the Sanandaj-Sirjan zone is about 50–55 km, which to some extent agrees with our finding, because our result shows the local Moho deepening to about 50 km in the southern part of this zone. However, Snyder and Barazangi (1986) reported the Moho depth of 65 km over this area, and Paul et al. (2006) mentioned that the Moho ranges from 50 to 90 km in the Sanandaj-Sirjan belt. Manaman et al. (2011) showed a significant crustal thickening under Zagros and the Sanandaj-Sirjan zone, with the maximum Moho depth reaching 65–70 km. The existing Moho models of CRUST1.0, M13_EURASIA, Meier and GEMMA indicate that the maximum Moho deepening to about 70 km is under Himalaya and Tibet Plateau. Therefore, it is unlikely that in Iran, the Moho could deepen to 70 or even 90 km. Meier’s Moho model has its deepest part in the southern Sanandaj-Sirjan zone and the Yazd block, which is in agreement with the result of Manaman et al. (2011); see Fig. 3. Our result also shows the localised Moho deepening in this region, characterised by large values of the vertical gravity gradients due to presence of heavy igneous rocks. Asudeh (1982) used seismic data over three profiles connecting Mashhad with Shiraz, Tehran with Mashhad, and Shiraz with Tabriz, and determined an average crustal thickness of 43, 45 and 46 km, respectively, along these three profiles. Javan Doloei (2003) and Javan Doloei and Ashtiany (2004) suggested the Moho depth of 52 km in Mashhad. Javan Doloei and Roberts (2003) estimated, based on processing seismic data, a crustal thickness of 46 km under Tehran region. Their result very closely agrees with the value reported by Asudeh (1982). Tehran regional is located in the southern part of central Alborz. Our result there indicates the Moho depth of about 40 km, thus more or less at the same range as in Mashhad, which are in agreement with the results of Asudeh (1982) and Javan Doloei and Roberts (2003). However, in Mashhad, the results of Asudeh (1982) and Javan Doloei (2003) and Javan Doloei and Ashtiany (2004) are different. A possible

22

M. Eshagh et al. / Journal of Asian Earth Sciences 138 (2017) 12–24

explanation is that Asudeh (1982) estimated the Moho depth along a profile from Tehran to Mashhad and not in Mashhad. However, our result does not show a deep Moho in Mashhad, where the isostatic model cannot reach such large values by means of compensating masses. Dehghani and Makris (1984) reported the Moho deepening to 35 km under the Alborz Mountains and 54 km in central Alborz. Sodoudi et al. (2009) used receiver functions to investigate the crustal and upper mantle structure beneath Central Alborz in the northern Iran. They analysed data from 290 teleseismic events at 12 short-period stations of the Tehran telemetric network. They estimated the average Moho depth between 44 and 46 km, which closely corresponds to our results. On the other hand, Manaman et al. (2011) reported the Moho depth within 55–60 km in the central Alborz, and Radjaee et al. (2010) found a crustal thickening from 48 km beneath the northern part of the Central Iranian Plateau to 55–58 km below the central part of the Alborz Mountains. Taghizadeh-Farahmand and Afsari (2015) reported the Moho depth of 54 km in the central Alborz. From these studies, the Moho depth estimates in the central Alborz are more or less similar. In contrast, our result indicates the Moho depth of about 45 km in that area. Paul et al. (2006) computed the Moho depth in UrumiehDokhtar magmatic area and reported the Moho deepening to 42 km. Nowrouzi et al. (2007) estimated there values around 44– 50 km using teleseismic data. Nasrabadi et al. (2008) reported an average crustal thickness of 40 km in the northwest Iran in Maku. Taghizadeh-Farahmand et al. (2010) used P and S seismic waves to investigate the crustal structure beneath the northwest Iran. According to their results, the average Moho depth there is 48 km, with variations between 38.5 and 53 km. TaghizadehFarahmand and Afsari (2015) reported an average Moho depth of 41 km in the northwest Iran, and the Moho variations from 45 km to 49 km in the northeast. Jiménez-Munt et al. (2012) reported a crustal thicknesses around 50 km under the Alborz and Kopet Dagh Mountains. According to our results, the Moho depth there reaches only about 40 km. Mangino and Priestley (1998) estimated the Moho depth of 33 km under the South Caspian and 30–33 km under the southwest and southeast Caspian basin, showing a good agreement with our results for this area. Manaman et al. (2011) reported similar values of the Moho depth between 30 and 33 km for the South Caspian Basin. Manaman et al. (2011) estimated that across the Makran subduction zone, the Moho depth increases from the Oman seafloor and Makran to 25–30 km. Taghizadeh-Farahmand and Afsari (2015) reported the Moho depth in Makran to be about 35 km. Our results show similar depths in Makran region, but along the subduction zone and tectonic margin in the Oman Sea the Moho should be about 17–18 km based on our estimates when disregarding the sediment and crustal density structures. According to Dehghani and Makris (1984), the Moho depth in the eastern Iran is between 45 and 48 km, while our models give values of about 40 km, thus closer to the value of 43 km estimated by Paul et al. (2006). Paul et al. (2006) also mentioned that Moho deepens more than 40 km under the Lut block, but our results indicate deepening not more than 30 km. According to Jiménez-Munt et al. (2012), the Moho depth is less than 36 km under the Lut block, and less than 38 km under the Jazmurian basin. Manaman et al. (2011) estimated the Moho depth between 35 and 40 km in the central Iran and under the Lut block. Nasrabadi et al. (2008) indicated that the Moho deepens to about 56 km beneath the Naein station in the central Iran. Sadidkhouy et al. (2012) inferred the Moho depth between 38.5 and 43 km in Isfahan area. From our results, the Moho depth in the central Iran is 35–40 km, thus closer to result of Sadidkhouy et al. (2012) for Isfahan, which is located in this area. Paul et al. (2006) suggested depths 36–52 km for the cen-

tral Iran, while our models show values only about 35–40 km. Paul et al. (2006) also estimated the Moho depth of 25 km for coastal regions along the Persian Gulf, which is in a good agreement with our results. These comparisons indicate a relatively good agreement between our results and existing seismic studies. This is true especially for most of the continental crust, where our results obtained after taking into consideration effects of sediments and consolidated crust have better agreement with local seismic studies. Significant differences in the estimated Moho depth were, however, found under the Oman Sea and the Makran subduction zone. In case of applying the sediment and consolidated crust corrections, our results there are unrealistically small. These findings indicate a low quality of the CRUST1.0 seismic and consolidated crust data offshore. Generally, our gravimetric Moho models are consistent with seismic studies in regions of Zagros, the northwest Iran, central Alborz, the Caspian Sea and Persian Gulf, including some parts of the eastern Iran. Therefore, we expect that our Moho model is reliable over these regions. In contrast, some discrepancies with respect to seismic studies were identified in Mashhad and Kopeh Dag, except for a relatively close agreement of our results with findings of Dehghani and Makris (1984). Our results also do not confirm likely the unrealistically large Moho deepening under the Sanandaj-Sirjan. 5. Summary and concluding remarks The regional differences between the Moho models computed according to Sjöberg and Jeffrey’s theoretical principles reach maxima of about 3 km after applying the combined TB correction. The reason is that Sjöberg’s method smoothes the Moho geometry more substantially than Jeffrey’s method. This was evident from the comparison of extreme values of the Moho depth obtained after applying both methods. We could see that the Moho maxima and minima are considerably reduced (in absolute sense) when applying Sjöberg’s method. In case of considering the sediment and consolidated crust corrections, our gravimetric results relatively closely agree with seismic studies over most of continental areas, but differ significantly under the Oman Sea and the Makran subduction zone. This could be due to a low quality of the CRUST1.0 sediment and consolidated crust data in these regions. Otherwise, our Moho models agree with the seismic Moho depth reported in most of the continental regions of Zagros, Alborz, the central Iran, the Caspian Sea, the Persian Gulf, but some discrepancies were found under the SanandajSirjan zone and Kopeh-Dag mountain ranges. The Moho depth computed by applying Sjöberg’s method does not exceed 50 km in Iran, and only slightly exceeds this value when applying Jeffrey’s method. Such finding holds also for the CRUST1.0, M13_Eurasia and Meier’s seismic models. However, some regional seismic studies show areas in Iran where Moho depth reaches 70 km or even 90 km. Such large Moho deepening is, however, likely unrealistic. The validation of our results revealed that the RMS of the Moho depth differences between our and CRUST1.0 (as well as Meier’s) is 5.1 km, whilst 5.6 km with respect to M13_Eurasia, and 4 km with GEMMA. In case of considering only the combined TB correction, the RMS improved to 2.9 km relative to GEMMA. Our models thus better agree with the GEMMA Moho model. This is due to using GOCE data and similar isostatic principles in both computations. Acknowledgment The first author is thankful to Mehrdad Eshagh for scientific discussions and helps about the geological structure of Iran.

M. Eshagh et al. / Journal of Asian Earth Sciences 138 (2017) 12–24

References Aghanabati, A., 2004. Iran Geology. Geological Survey of Iran. Airy, G.B., 1855. On the computations of the effect of the attraction of the mountain masses as disturbing the apparent astronomical latitude of stations in geodetic surveys. Trans. Roy. Soc. Lond. 145 (B), 101–104. Asudeh, I., 1982. Seismic structure of Iran from surface and body wave data. Geophys. J. R. Astron. Soc. 71, 715–730. Barzaghi, R., Reguzzoni, M., Borghi, A., De Gaetani, C., Sampietro, D., Marotta, A.M., 2013. Global to local Moho estimate based on GOCE geopotential model and local gravity data. In: Proceedings of IAG 2013 Symposium. Bassin, C., Laske, G., Masters, T.G., 2000. The current limits of resolution for surface wave tomography in North America. EOS Trans AGU 81 (F897), 2000. Braitenberg, C., Ebbing, J., 2009. New insights into the basement structure of the westSiberian basin from forward and inverse modelling of Grace satellite gravity data. J. Geophys. Res. 114, B06402. http://dx.doi.org/10.1029/2008JB005799. Braitenberg, C., Wienecke, S., Wang, Y., 2006. Basement structures from satellitederived gravity field: South China Sea ridge. J. Geophys. Res. 111, B05407. http://dx.doi.org/10.1029/2005JB003938. Braitenberg, C., Zadro, M., Fang, J., Wang, Y., Hsu, H.T., 2000. The gravity and isostatic Moho undulations in Qinghai-Tibet plateau. J. Geodyn. 30 (5), 489– 505. http://dx.doi.org/10.1016/S0264-3707(00)00004-1. Dehghani, G., Makris, J., 1984. The gravity field and crustal structure of Iran. N. Jahrb. Geol. Palaont. Abh. Stuttgart 168, 215–229. Eshagh, M., 2009a. On satellite gravity gradiometry Doctoral dissertation in geodesy. Royal Institute of Technology (KTH), Stockholm, Sweden. Eshagh, M., 2009b. The effect of lateral density variations of crustal and topographic masses on GOCE gradiometric data: a study in Iran and Fennoscandia. Acta Geodaetica et Geophysica Hungarica 44 (4), 399–418. http://dx.doi.org/ 10.1556/AGeod.44.2009.4.3. Eshagh, M., 2011. The effect of spatial truncation error on integral inversion of satellite gravity gradiometry data. Adv. Space Res. 47, 1238–1247. http://dx.doi. org/10.1016/j.asr.2010.11.035. Eshagh, M., 2014a. Determination of Moho discontinuity from satellite gradiometry data: linear approach. Geodyn. Res. Int. Bull. 1 (2), 1–13. Eshagh, M., 2014b. From satellite gradiometry data to subcrustal stress due to mantle convection. Pure Appl. Geophys. 171, 2391–2406. http://dx.doi.org/ 10.1007/s00024-014-0847-2. Eshagh, M., 2016a. A theoretical discussion on Vening-Meinesz-Moritz inverse problem of isostatsy. Geophys. J. Int. 207, 1420–1431. Eshagh, 2016b. On Vening Meinesz-Moritz and flexural theories of isostasy and their comparison over Tibet Plateau. J. Geodetic Sci. 6, 139–151. Eshagh, M., Sjöberg, L.E., 2011. Determination of gravity anomaly at sea level from inversion of satellite gravity gradiometric data. J. Geodyn. 51 (5), 366–377. http://dx.doi.org/10.1016/j.jog.2010.11.00. Eshagh, M., Hussain, M., 2015. Relationship amongst gravity gradients, deflection of vertical, Moho deflection and the stresses derived by mantle convection: a case study over Indo-Pak and surroundings. Geodyn. Res. Int. Bull. 3 (4), 2345–4997. ISSN 2345-4997. Eshagh, M., Hussain, M., 2016. An approach to Moho discontinuity recovery from on-orbit GOCE data with application over Indo-Pak region. Tectonophysics 690 (B), 253–262. Eshagh, M., Hussain, M., Tenzer, R., Romeshkani, M., 2016. Moho density contrast in central Eurasia from GOCE gravity gradients. Remote Sens. 8 (5), 418. http://dx. doi.org/10.3390/rs8050418. Farr, T.G., Paul, A.R., Edward, C., Robert, C., Riley, D., Scott, H., Michael, K., Mimi, P., Ernesto Rodriguez, Ladislav, R., David, S., Shaffer, S., Joanne, S., 2007. The Shuttle Radar Topography Mission, Rev. Geophys. 45, RG2004, http://dx.doi.org/10. 1029/2005RG000183. Fowler, C.M.R., 1990. The Solid Earth: An Introduction to Global Geophysics. Cambridge University Press, Cambridge, UK. Ghorbani, M., 2013. The Economic Geology of Iran. Mineral Deposits and Natural Resources, Springer. Gladkikh, V., Tenzer, R., 2011. A mathematical model of the global ocean saltwater density distribution. Pure Appl. Geophys. 169 (1–2), 249–257. Gomez-Oritz, D., Agarwal, B.N.P., 2005. 3DINVER.M: a MATLAB program to invert the gravity anomaly over a 3D horizontal density interface by ParkerOldenburg’s algorithm. Comput. Geosci. 31, 13–520. Hatzfeld, D., Tatar, M., Priestley, K., Ghafory Ashtiany, M., 2003. Seismological constraints on the crustal structure beneath the Zagros Mountain belt (Iran). Geophys. J. Int. 155, 403–410. http://dx.doi.org/10.1046/j.1365-246X.2003. 02045.x. Hansen, P.C., 1998. Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion. SIAM, Philadelphia. Heiskanen, W., Moritz, H., 1967. Physical Geodesy. W.H. Freeman and company, San Francisco and London. Heiskanen, W.A., 1931. New isostatic tables for the reduction of the gravity values calculated on the basis of Airy’s hypothesis, Publ. Isos. Inst. IAG (Helsinki), No. 2, 1931. Janák, J., Fukuda, Y., Xu, P., 2009. Application of GOCE data for regional gravity field modelling. Earth Planet. Sci. 61, 835–843. Javan Doloei, G., 2003. Transfer/Receiver functions and its application on crustal structure of Uppsala, Tehran and Mashhad area Ph.D Thesis. International Institute of Earthquake Engineering and Seismology (IIEES).

23

Javan Doloei, J., Roberts, R., 2003. Crust and uppermost mantle structure of Tehran region from analysis of teleseismic P-waveform receiver functions. Tectonophysics 364, 115–133. Javan Doloei, G., Ashtiany, M.G., 2004. Crustal structure of Mashhad Area from time domain receiver functions analysis of teleseismic earthquake. Res. Bull. Seismol. Earthq. Eng. 27, 30–38. Jeffrey, H., 1976. The Earth: Its Origin, History and Physical Constitution. Cambridge University Press. Jiménez-Munt, I., Fernàndez, M., Saura, E., Vergés, J., Garcia-Castellanos, D., 2012. 3D lithospheric structure and regional/residual Bouguer anomalies in the ArabiaEurasia collision (Iran). Geophys. J. Int. 190, 1311–1324. http://dx.doi.org/ 10.1111/j.1365-246X.2012.05580.x. Kiamehr, R., Gomes-Ortiz, D., 2009. A new 3D Moho depth model for Iran based on the terrestrial gravity data and EGM2008 model. EGU General Assembly, Vol. 11, EGU2009-321-1, Vienna, Austria (Abstract). Laske, G., Masters, G., Ma, Z., Pasyanos, M., 2013. Update on CRUST1.0 – A 1-degree Global Model of Earth’s Crust. In EGU General Assembly Conference Abstracts 15, 2658. Manaman, N.S., Shomali, H., Koyi, H., 2011. New constraints on upper mantle S-velocity structure and crustal thickness of the Iranian plateau using partitioned waveform inversion. Geophys. J. Int. 184, 247–267. http://dx.doi. org/10.1111/j.1365-246X.2010.04822.x. Mangino, S., Priestley, K., 1998. The crustal structure of the southern Caspian region. Geophys. J. Int. 133, 630–648. Meier, U., Curtis, A., Trampert, J., 2007. Global crustal thickness from neural network inversion of surface wave data. Geophys. J. Int. 169, 706–722. http://dx.doi.org/ 10.1111/j.1365-246X.2007.03373.x. Moritz, H., 1990. The Figure of the Earth. H Wichmann, Karlsruhe. Moritz, H., 2000. Geodetic reference system 1980. J. Geodesy 74, 128–133. Nabavi, M.H., 1976. History of Iran Geology, Geological Survey and Mineral Exploration Country. 109 p. Nasrabadi, A., Tatar, M., Priestley, K. & Sepahvand, M.R., 2008. Continental lithosphere structure be neath the Iranian Plateau, from analysis of receiver functions and surfaces waves dispersion, The 14th World Conference on Earthquake Engineering, 2008 October 12–17, Beijing, China. Nowrouzi, G., Priestley, K., Ghafory-Ashtiany, M., Doloei, J., Rham, D., 2007. Crustal velocity structure in the Kopet-Dagh, from analysis of P-waveform receiver functions. J. Seismol. Earthq. Eng. 8, 187–194. Oldenburg, D.W., 1974. The inversion and interpretation of gravity anomalies. Geophysics 39, 526–536. Parker, R.L., 1972. The rapid calculation of potential anomalies. Geophys. J. Roy. Astron. Soc. 31, 447–455. 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. 166, 227–237. Pavlis, N., Holmes, S.A., Kenyon, S.C., Factor, J.K., 2008. An Earth Gravitational model to degree 2160: EGM2008. In: Presented at the 2008 General Assembly of the European Geosciences Union, Vienna, Austria, April 13–18, 2008. Pratt, J.H., 1854. On the attraction of Himalayan mountain and of the elevated regions beyond them, upon the plumb line in India. Philos. Trans. Roy. Soc. Lond. 146, 53–100. Radjaee, A., Rham, D., Mokhtari, M., Tatar, M., Priestley, K., Hatzfeld, D., 2010. Variation of Moho depth in the central part of the Alborz Mountains, northern Iran. Geophys. J. Int. 181, 173–184. http://dx.doi.org/10.1111/j.1365246X.2010.04518.x. Reed, G.B., 1973. Application of kinematical geodesy for determining the short wave length components of the gravity field by satellite gradiometry (No. DGS-201). Ohio State Univ Columbus Dept of Geodetic Science. Reguzzoni, M., Sampietro, D., Sanso, F., 2013. Global Moho from the combination of the CRUST2.0 model and GOCE data. Geophys. J. Int. 195 (1), 222–237. http://dx. doi.org/10.1093/gji/ggt247. Reguzzoni, M., Sampietro, D., 2012. A new global crustal model based on GOCE data grids, Presented at first international GOCE solid Earth workshop, Enschede, The Netherlands. Reguzzoni, M., Sampietro, D., 2014. GEMMA: an earth crustal model based on GOCE satellite data. Int. J. Appl. Earth Obs. Geoinf. 35, 31–43. http://dx.doi.org/ 10.1016/j.jag.2014.04.002. Rexer, M., Hirt, C., Claessens, S., Tenzer, R., 2016. Layer-based modelling of the Earth’s gravitational potential up to 10-km scale in spherical harmonics in spherical and ellipsoidal approximation. Surv. Geophys. http://dx.doi.org/ 10.1007/s10712-016-9382-2. Sadidkhouy, A., Alikhani, Z., Sodoudi, F., 2012. The variation of Moho depth and VP/ VS ratio beneath Isfahan region, using analysis of teleseismic receiver function. J. Earth Space, Phys. 38 (3), 15–24. Safari, A., Sharifi, M., Bahroudi, A., Parang, S., 2014. Determination of the best window size and structural index in estimating Moho depth through Euler deconvolution method (Case Study: The Zagros Zone). J. Geomatics, Sci. Technol. 4 (2), 123–137. Sampietro, D., 2011. GOCE exploitation for Moho modelling and applications. In: Proceedings of the 4th international GOCE user workshop, Munich, Germany, vol. 31. Sampietro, D., Reguzzoni, M., 2011. A study on the Austrian Moho from GOCE data. Poster Presented at European Geosciences Union General Assembly 2011 Vienna, Austria, 08 April 2011.

24

M. Eshagh et al. / Journal of Asian Earth Sciences 138 (2017) 12–24

Sampietro, D., 2009. An inverse gravimetric problem with GOCE data. Ph.D. thesis. Politecnico di Milano Polo regionale di Como. Sampietro, D., Reguzzoni, M., Braitenberg, C., 2014. The GOCE Estimated Moho Beneath the Tibetan Plateau and Himalaya. C. Rizos and P. Willis (Eds.), Earth on the Edge: Science for a Sustainable Planet, International Association of Geodesy Symposia 139, http://dx.doi.org/10.1007/978-3-642-37222-3-52. Shin, Y.H., Xu, H., Braitenberg, C., Fang, J., Wang, Y., 2007. Moho undulations beneath Tibet from GRACE-integrated gravity data. Geophys. J. Int. 170, 971– 985. http://dx.doi.org/10.1111/j.1365-246X.2007.03457.x. Sjöberg, L.E., 2009. Solving Vening Meinesz-Moritz inverse problem in isostasy. Geophys. J. Int. 179, 1527–1536. http://dx.doi.org/10.1111/j.1365246X.2009.04397.x. 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 (3), 361–373. Sodoudi, F., Yuan, X., Kind, R., Heit, B., Sadidkhouy, A., 2009. Evidence for a missing crustal root and a thin lithosphere beneath the Central Alborz by receiver function studies. Geophys. J. Int. 177, 733–742. http://dx.doi.org/10.1111/ j.1365-246X.2009.04115.x. Stöcklin, J., 1968. Structural history and tectonics of Iran, A review. AAPG Bull. 52 (7), 1229–1258. Stolk, W., Kaban, M., Beekman, F., Tesauro, M., Mooney, W.D., Cloetingh, S., 2013. High resolution regional crustal models from irregularly distributed data: application to Asia and adjacent areas. Tectonophysics 602, 55–68. http://dx. doi.org/10.1016/j.tecto.2013.01.022. Sünkel, H., 1985. An isostatic Earth model, Report 367. Department of Geodetic Science and Surveying. Ohio State University, Columbus. Ohio. Taghizadeh-Farahmand, F., sodoudi, F., Afsari, N., Ghassemi, M.R., 2010. Lithospheric structure of NW Iran from P and S receiver functions. J. Seismol. 14, 823–836. Taghizadeh-Farahmand, F., Afsari, N., 2015. Determination of Crustal thickness beneath INSN stations in Iran Plateau by using of the modelling of the P receiver function, Publication of new research in earthquake, Volume 10, Issue 35.65-72. Tenzer, R., Hamayun, K., Vajda, P., 2009. Global maps of the CRUST2.0 components stripped gravity disturbance. J. Geophys. Res. (Solid Earth) 114, B05408. http:// dx.doi.org/10.1029/2008JB006016. Tenzer, R., Novák, P., Gladkikh, V., 2011. On the accuracy of the bathymetrygenerated gravitational field quantities for a depth-dependent seawater density distribution. Stud. Geophys. Geod. 55 (4), 609–626 (SCI).

Tenzer, R., Novák, P., Gladkikh, V., 2012a. The bathymetric stripping corrections to gravity field quantities for a depth-dependent model of the seawater density. Mar. Geodesy 35, 198–220. Tenzer, R., Hamayun, Novak, P., Gladkikh, V., Vajda, P., 2012b. Global crust-mantle density contrast estimated from EGM2008, DTM2008, CRUST2.0 and ICE-5G, Pure and Applied Geophysics 169 (9), pp. 1663–1678, http://dx.doi.org/10. 1007/s00024-011-0410-3. Tenzer, R., Gladkikh, V., Vajda, P., Novák, P., 2012c. Spatial and spectral analysis of refined gravity data for modelling the crust-mantle interface and mantlelithosphere structure. Surv. Geophys. 33 (5), 817–839. Tenzer, R., Chen, W., 2014a. Regional gravity inversion of crustal thickness beneath the Tibetan plateau. Earth Sci. Inf. 7, 265–276. Tenzer, R., Chen, W., 2014b. Expressions for the global gravimetric Moho modeling in spectral domain. Pure Appl. Geophys. 171 (8), 1877–1896. Tenzer, R., Chen, W., Jin, S., 2015a. Effect of the upper mantle density structure on the Moho geometry. Pure Appl. Geophys. 172 (6), 1563–1583. Tenzer, R., Chen, W., Tsoulis, D., Bagherbandi, M., Sjöberg, L.E., Novák, P., Jin, S., 2015b. Analysis of the refined CRUST1.0 crustal model and its gravity field. Surv. Geophys. 36 (1), 139–165. Tikhonov, A.N., 1963. Solution of incorrectly formulated problems and regularization method, Soviet Math. Dokl., 4: 1035–1038, English translation of Dokl. Akad. Nauk. SSSR 151, 501–504. Turcotte, D., Schubert, G., 2014. Geodynamics. Cambridge University Press. Vajda, P., Ellmann, A., Meurers, B., Vanícek, P., Novák, P., Tenzer, R., 2008. Gravity disturbances in regions of negative heights: a reference quasi-ellipsoid approach. Stud. Geophys. Geod. 52, 35–52. Vening Meinesz, F.A., 1931. Une nouvelle methode pour la reduction isostatique regionale de l’intensite de la pesanteur. Bull. Geod. 29 (1), 33–51. Xu, P., 1992. Determination of surface gravity anomalies using gradiometric observables. Geophys. J. Int. 110, 321–332. http://dx.doi.org/10.1111/j.1365246X.1992.tb00877.x. Xu, P., 1998. Truncated SVD methods for discrete linear ill-posed problems. Geophys. J. Int. 135, 505–514. http://dx.doi.org/10.1046/j.1365246X.1998.00652.x. Xu, P., 2009. Iterative generalized cross-validation for fusing heteroscedastic data of inverse ill-posed problems. Geophys. J. Int. 179, 182–200. http://dx.doi.org/ 10.1111/j.1365-246X.2009.04280.x.