Evaluation of different gravimetric methods to Moho recovery in Iran

The complexity of geological units in Iran because of several unique events like tectonics and orogenic activities in this region led to extensive investigations for Moho recovery by seismic methods therein. In this research, three gravimetric methods have been evaluated by some point-wise seismic data. We applied collocation method as an iterative process as well as modified forms of Sjöberg and Jeffrey’s theory of isostasy for local Moho depth recovery. The gravity data has been generated by GOCO03S model reduced by topography/bathymetry, sediment and consolidated crust effects. Although the iteration process in collocation approach only slightly changed the estimated depths, this method led to a better agreement with seismic data rather than others. Differences between collocation, Jeffrey and Sjöberg’s solutions with seismic studies are similar but Jeffrey and Sjöberg’s methods displayed a systematic bias. The standard deviations of the residuals among seismic data and gravimetric solutions are around 6 km. Overall, the evaluation of these approaches indicated that Moho from gravimetric approaches reduced only slightly the standard deviation of seismic Moho estimates. Significant discrepancies with seismic data have been detected in Makran subduction zone, Oman Sea, Persian Gulf and Caspian Sea. The explanation of such inconsistency can be partially due to the poor quality of CRUST1.0 data in these areas, as this model has been used to correct the gravity values that were input in the inversion procedures.


Introduction
The knowledge of the Mohorovičić discontinuity (Moho) can provide valuable information to understand some topical issues in solid Earth sciences (Sampietro, 2016). Its knowledge in the Iran block is one of the crucial issues in this region. Several unique events like tectonics and orogenic activities in Iran led to a complex geological structure of the area. Thus, it is important to study in deep the peculiar structure of the Iranian Moho applying different methods and different observations in order to have a comprehensive definition of its main features. Moho interface is commonly estimated by either seismic or gravimetric methods.
Although seismic Moho estimates have a significant accuracy (at the level of 1-2 km), their coverage over entire of the Earth is quite poor. Thus, in regions where seismic data are sparse or missing, results of gravimetric studies can be profitably used. At global scale, this has been made possible after dedicated gravity-satellite missions, namely the Gravity Recovery and Climate Experiment (GRACE) (Tapley et al., 2004a, Tapley et al., 2004b and the Gravity field and steady-state Ocean Circulation Explorer (GOCE) (Floberghagen et al., 2011). These two missions have provided global gravity field with an accuracy and a resolution which is suitable for investigating the Moho structure. Among successful seismic estimates, the early results dated back to Beloussov et al. (1980). The most widely known crustal model based on seismic refraction is CRUST5.1 model (Mooney et al., 1998). Bassin et al. (2000) further upgraded it and called the new version, CRUST2.0. The CRUST1.0 (Laske et al., 2013) including ice layers, water, sediments and consolidated crustal layers, is the most recent version compiled with a 1×1 arc-deg spatial resolution. As already pointed out, because of insufficient seismic data coverage over large areas, gravimetric or combined gravimetric/seismic solutions have been utilized. Based on some isostasy hypotheses on compensating the Earth's topographic masses, gravity data can be used to determine the Moho depth. Several basic theories were suggested to explain the mechanism of isostasy like those proposed by Pratt-Hayford (Pratt, 1855, Hayford, 1909 and Airy-Heiskanen (Airy, 1855, Heiskanen, 1931. Vening Meinesz (1931) modified the Airy-Heiskanen theory by considering a regional instead of a local compensation based on a thin plate lithospheric flexure model (Watts, 2001). Vening Meinesz theory had been modified by Parker (1973) in an iterative approach for Moho determination and Oldenburg (1974) made an attempt to stabilise this method by applying a low-pass filtering technique. The combination of these two methods was known as a Parker-Oldenburg method and it has been generalized for the 3-D gravity inversion by Gómez-Ortiz and Agarwal (2005) and Shin et al. (2007). Braitenberg et al. (2000) developed a similar method with integration of seismic data; also they estimated variation of the Moho under the Tibet plateau by an iterative inversion method. Moritz (1990) improved the Vening Meinesz inverse problem for a global compensation by adopting the spherical approximation model of the Earth. Since there were some theoretical deficiencies in this isostatic method, Sjöberg (2009) reformulated Moritz's theory and called it as the Vening Meinesz Moritz (VMM) problem. In this approach he solved a non-linear Fredholm's integral equation of the first kind. Bagherbandi and Sjöberg (2012) made a comparison between the gravimetric VMM and local Airy-Heiskanen methods in the determination of Moho depth.
Another approach based on the inversion of gravity data to determining the Moho depth is collocation method (Krarup, 1969, Tscherning, 1985, Moritz, 1990. Barzaghi et al. (1992) proposed an approach based on collocation principle, which propagate the covariance structure of the depth interface to the covariance function of the observed gravity field. Barzaghi et al. (2015) applied collocation method by combining the global gravity-gradient information of GOCE and local gravity data to Moho recovery. Barzaghi and Biagi (2014) implemented then a further version of the collocation method by including seismic Moho depths as input data. Braitenberg and Ebbing (2009) studied the structure of the crust by combination of GRACE and terrestrial gravity data. Some other scientists made attempt to estimate Moho depth by applying GOCE gravity gradient data (Braitenberg et al., 2010, Sampietro, 2011, Sampietro et al., 2014. Reguzzoni et al. (2013) combined seismic and GOCE data to obtain a new global Moho model. Tenzer et al. (2009Tenzer et al. ( , 2011Tenzer et al. ( , 2012aTenzer et al. ( , 2012b proved that applying the crustal density-contrast stripping corrections is appropriate for a gravimetric Moho recovery. They showed that gravitational contributions of topography and major known crustal density structures have a large spatial correlation with the Moho geometry. In the area under investigation in this paper, several studies based on seismic data have been performed for the determination of the regional Moho model. Most of these studies concentrated on a specific area, for instance profiles between Shiraz-Mashhad, 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 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, Shad Manaman et al., 2011), Kopeh-Dagh (Nowrouzi et al., 2007, Naein (Nasrabadi et al., 2008), the northwest Iran (Taghizadeh-Farahmand et al., 2015), and Sanandaj-Sirjan zone (Sadidkhouy et al., 2012). Furthermore, based on terrestrial gravity data in this area with a quite homogenously coverage, Dehghani and Makris (1984) combined gravimetric and seismic data to determination of the crustal structure of Iran. Abbaszadeh et al. (2013) compared the effective elastic thickness of the lithosphere estimated by terrestrial and satellite data in Iran. Eshagh et al. (2017) reformulated two isostatic methods for Moho recovery by considering contributions of mean Moho over the whole Moho spectrum.
In this paper, we adopt collocation method as well as two isostatic approaches to determine the Moho depths inverting gravity data over Iran. The collocation inversion method for a twolayer model devised by Barzaghi and Biagi (2014) is applied. In addition, the generalized form of Jeffrey and Sjöberg's method proposed by Eshagh et al. (2017) have been utilized in this study. These methods have been applied to gravity from GOCO03S gravitational model (Mayer-Gürr et al., 2012) reduced for the SRTM30_PLUS topographic/bathymetric data (Becker et al., 2009), the sediment and the crystalline data of CRUST1.0 (Laske et al., 2013).
Comparisons of these Moho depth estimates have been finally performed with the seismic estimates available from literature.

The theoretical background
In this section, we review three gravimetric inversion methods that can be successfully applied to estimate the Moho depth. At first, we describe the procedure based on the collocation principle. Subsequently, we also discuss two alternative isostatic approaches, which are reformulations of Jeffrey and Sjöberg's theories.

The collocation solution
The collocation method for Moho estimation is based on a stochastic approach derived from the Wiener filtering and prediction theory that allows estimating the signal correlated component based on the covariance structure of the data (Barzaghi and Biagi, 2014). In the years, it has been applied to determine the gravimetric estimate of the Moho (Barzaghi et al., 1992, Barzaghi and Biagi, 2014, Barzaghi et al., 2015. In this context, the covariance structure of the Moho depth is propagated to the covariance of the observed gravity in a simple two-layer model. In order to apply collocation for estimating the Moho depth, we should consider the following linear relationship which holds in planar approximation (Barzaghi and Biagi, 2014 Also, the following conditions must be fulfilled (Barzaghi and Biagi, 2014): 1) is a weak stationary stochastic process, ergodic in the mean and in the covariance 2) The noises in gravity and depth, g n and n are spatially uncorrelated zero mean signals 3) The cross-correlations between signals and noises are zero To perform the computation, the gravity auto-covariance and the cross-covariances between gravity and depth are needed, i.e.: Assuming that the observed gravity values contain a noise component, we can write: The collocation estimate of is (Moritz, 1980, Barzaghi et al., 1992, Barzaghi and Biagi, 2014: As the first step to use this approach, the empirical covariance function of Dg OBS should be estimated and modelled with appropriate positive definite model functions (Moritz, 1980).
Under the hypotheses previously stated, the empirical covariance can be estimated as (Barzaghi et al., 1992): Also, auto and cross-covariance models must be defined in order to derive the estimator of for Moho determination. One possible model for the auto covariance of Dg (see Barzaghi and Biagi (2014)) is: Where 1 (. ) is the first order Bessel function.
If this auto-covariance is considered, one can prove (Barzaghi et al., 1992) that:  Jeffrey (1976) has solved the problem of isostasy for Moho modelling in a very similar way to the VMM method (Sjöberg, 2009). Eshagh and Hussain (2016) proposed a different modelling of the isostatic gravity disturbance I g by considering the gravitational effects of topography and bathymetry, sediments and consolidated crystalline basement as follows:

The Jeffrey's solution
(2-9) In this equation, g denotes the observed gravity disturbance which is the difference between measured and normal gravity at a computation point on the Earth, C g is the compensation attraction, TB g is the topographic/bathymetric effect on g , S g and Crys g represent the gravitational effects of sediment and consolidated crust, respectively. Moritz (1990) and Sjöberg (2009) applied this equation as a fundamental condition to solve the VMM problem.
In order to define the compensation potential in equation (2-10), Eshagh et al. (2017) arranged the formula in Jeffrey (1976) in the following way: As stated before, is the variation of the Moho depths relative to the mean value indicated by T0. In this equation, R is the radius of the Earth, σ is the unit sphere, r´ denotes the geocentric distance of the infinitesimal mass element while dr´ and d sin d d σ θ θ λ are the radial and the surface integration elements, cos n P ψ is the Legendre polynomial of degree n for the argument of the geocentric angle ψ . θ and λ represent the spherical co-latitude and longitude of the element, l is the Euclidean spatial distance which is a function of the r´ (see. Heiskanen and Moritz (1967)).
Further, Eshagh et al. (2017) solved the radial integral and expanded C g in a spectral form according to the Heiskanen and Moritz (1967) Also Eshagh et al. (2017) approximated the term [ / (R-T0)] n+3 by a binomial series and after simplification, (2-11) reduces to: Since is the undulation of Moho depths with respect to 0 T , the total Moho depth is given as

The Sjöberg's solution
The Vening Meinesz Moritz (VMM) inverse problem of isostasy (Sjöberg, 2009) has been developed and applied successfully over different areas of the Earth. The main difference between Jeffrey and Sjöberg's method is how to write the compensation potential. Sjöberg (2009) used the compensation potential VC derived by Moritz (1990): Eshagh et al. (2017) solved the integrations by assuming T = T0 and obtained the compensation attraction of the gravity disturbance as follows: Furthermore, they considered the spectral form of C δg and got: They proved by inserting Eq. (2-18) into (2-9) that the Moho depth can be obtained as: As it can be seen in Eq. (2-19), T0, besides the zero-degree term, affects all frequencies.

The Iran case study
In this section, we present the Iran case study and the application of the three methods previously described to the gravimetric estimate of the Moho in Iran. We divide this section into four parts. In section 3.1 we present a brief geological structure of the area. The used gravity data are described in section 3.2 while in section 3.3 the local gravity Moho estimates over the study area by collocation, the Jeffrey and Sjöberg methods are presented. Finally, in section 3.4 comparisons among different gravimetric and seismic derived Moho values are shown.

The study area
All methodologies have been applied in Iran over an area limited by the parallels 20° and 45° North and the meridians 40° and 65° East (see Fig. 1). Several unique events like tectonics and orogenic activities led to complicated structural units in this area. By considering the geological features of Iran, some models and interpretations have been proposed for this region (Nabavi, 1976, Eftekharnezhad, 1980, Alavi-Naini, 1993, Aghanabati, 2004, Ghorbani, 2013. According to these studies, a geological setting of this area showing a geological classification of various structural zones of Iran has been devised and is presented in Fig. 1, superimposed on the regional topography. Topography/bathymetric heights were generated by the SRTM30_PLUS model to degree and order 2160 with a resolution of 5'×5' over the study area (Becker et al., 2009). As seen in Fig. 1 King, 1981, Berberian et al., 1982).

The gravity data set
In this study, the considered gravity data have been synthetized from the GOCO03S gravitational model up to degree and order 180 and were computed on a 0.5×0.5 arc-deg surface grid. In order to obtain the topography/bathymetry (TB) corrections, we applied the spherical harmonics expansion of digital elevation model from SRTM30_PLUS to degree and order 180 consistent with the resolution of CRUST1.0 model. Moreover, for computing the gravity corrections due to sediment and consolidated crustal layers, we used the Earth's crustal model CRUST1.0, which in this area gives mean densities values for the sediments and the consolidated crustal layers that are, respectively, 2060 kg m -3 and 2670 kg m -3 (see Fig. 2b and 2c). By applying all these corrections we obtained the gravity data that were used in the inversion procedures. Fig. 2d represents the map of the refined Bouguer gravity, in unit of mGal, reduced for the gravitational contribution of crustal density heterogeneities over the study area. This illustrates the largest and positive value at Oman Sea and the negative one over the Zagros and Sanandaj-Sirjan belts and in the North-East part of Iran around its border to Azerbaijan and Turkey and in the south part in Bam. All the mentioned corrections are displayed in Fig. 2 and the related statistics are given in Table. 1.

The gravimetric Moho estimates
The three methods described in section 2 for determination of Moho were applied in the study area. As stated above, the collocation approach in planar approximation has been applied to GOCE derived reduced Bouguer gravity, i.e. Bouguer gravity anomalies reduced for sediment and consolidated crust effects. Moho depths have been determined with respect to mean value 0 T while the values have been obtained according to the scheme described in section 2. In this computation, the constant density contrast between crust and mantle has been set to 600 kg m -3 . The gravity data, obtained by applying combined topography/bathymetry gravitational effect from SRTM30_PLUS data and corrections for the sediment and consolidated crust data generated from the CRUST1.0, has been regridded on a regular (x,y) grid in the investigation area (see Fig. 3):

Fig. 3 Bouguer gravity anomaly reduced by topography/bathymetry, sediment and consolidated crust corrections on a (x,y) grid [mGal]
The collocation procedure has been then applied iteratively. The behaviour of the empirical covariance function of gravity data and the best-fit model with the related parameters can be seen in Fig. 4 for the three steps that have been performed to get the final estimate.  (Barzaghi et al., 1992). In the first step, the model function which best described the empirical values is J1 Bessel function divided by its argument, with A=14000 mGal 2 and parameter ˆ =0.0055 km -1 . A value is sufficiently close to initial one of empirical function, which is 14815 mGal 2 . The difference between the two values in the origin, i.e. the noise variance, has been fixed to 815 mGal 2 . In the second step, residuals of gravity from the first inversion step have been used as input data. As can be seen in Fig. 4 the variance of these residuals, i.e. the value in the origin of the empirical function, decreased drastically to 3490 mGal 2 and the signal variance, i.e. the A value, has been set to 3200 mGal 2 . The ˆ and 2 n quantities have been fixed at 0.0135 km -1 and 290 mGal 2 , respectively. Furthermore, in third step, based on the residual gravity from the second iteration, model covariance parameters were fixed at ˆg A 900 mGal 2 , ˆ 0.0317 km -1 , 2 n 192 mGal 2 . This iterative process has been stopped at the third iteration since the covariance function of the third step residuals has a correlation length (the distance at which the covariance function is half of its value in the origin) that is comparable with the grid step.
Results of this procedure are shown in Fig. 5 and related statistics are summarized in Table 2  The final estimate (see Fig. 5c In first step it can be observed that the Moho depth varies between 35.5 km and 49.0 km and in most areas, depth of Moho is below 40 km. In second step, statistics are similar to previous with a slight increase of the standard deviation of the estimated Moho depths. As seen in Fig.   5c the Moho depths coming from the third step are more high frequency (standard deviation increases to 3.4 km) and ranges between 31.6 km and 54.6 km.
Based on gravity disturbance data, coming from the GOCO03S, completed to degree and order 180, we then estimated the Moho depth over Iran with the gravimetric approaches devised by Sjöberg and Jeffrey. As done in the collocation solution, the mean Moho depth has been set to 44 km. Also, the constant density contrast between crust and mantle has been set to 600 kg m -3 and corrections including topography/bathymetry, sediment and consolidated crust has been accounted as well. The results are plotted in Fig. 6 and the summary of statistics is given in Table.  As seen in Fig. 6, Sjöberg method gives a smoother solution than the one based on Jeffrey method.
Minimum Moho depths in both these methods are lower than minimum depths estimated using the collocation method. Statistics in Table. 3 show a minimum Moho depth of 17.6 km and 33.9 km for Jeffrey and Sjöberg methods, respectively. Also, the maximum Moho depth according to Jeffrey solution seems to be overestimated since this method has given a depth around 67.4 km near the northwest of Iran. On the contrary, by applying the Sjöberg method, the maximum Moho depth in Iran has been estimated to 59.6 km.

Comparing gravimetric and seismic Moho estimates
To validate the gravimetric Moho solutions, we compared them with existing regional seismic studies for Iran (Mangino and Priestley, 1998, Paul et al., 2006, Taghizadeh-Farahmand et al., 2010, Radjaee et al., 2010, Tatar and Nasrabadi, 2013, Taghizadeh-Farahmand et al., 2015, Motaghi et al., 2015, Abdollahi et al., 2018. A compilation of these numerous seismic datasets has been prepared and checks for their consistency were performed. In this way, we defined a selected collection of seismic Moho values in the Iran area, which consists of 277 points. These models are shown in Fig. 7 and their statistics are given in Table 4. The maximum Moho depth is located under the Sanandaj-Sirjan zone and surrounding mountains.
By seismic estimates, the minimum depth of Moho is under the Oman Sea and Makran subduction zone. As it can be seen, Moho depths derived from seismic studies varies between 18.5 km and 66 km.  These seismic values were then compared with the gravity derived estimates. In order to perform the comparison, we interpolated the gridded gravity estimates on the sparse seismic point using linear interpolation.
The differences are plotted in Fig. 8 and the statistics of the differences are summarized in Table. 5. As one can see in Fig. 8a, the differences between the collocation solution and seismic data are between -17 km and 21.9 km. The differences between Jeffrey's and Sjöberg's solutions and seismic data range from -11.1 km to 22.2 km and -10.8 km to 23.4 km, respectively (see Fig. 8b and 8c).  As it can be seen in Table 5, the STDs of the differences with respect to seismic data are around 6 km, the smaller STD values being obtained by the Sjöberg estimate. Since the STD of seismic data 8.2 km (see This overall analysis can be further specified for different sub-areas in the Iran region where existing regional seismic solutions are available.

Seismic
Central  Sadidkhouy et al. (2012) showed that depth of Moho in Isfahan area is variable between 38.5 and 43 km. Our solutions by gravimetric methods give the Moho depth ranging from 47 km in central Iran to 44 km in Lut block, which is in a good agreement with Paul et al. (2006) and Sadidkhouy et al. (2012). In the coast of the Persian Gulf, a Moho depth of about 25 km has been suggested by Paul et al. (2006) where we have different values (here our estimates are around 35 km). All these comparisons between seismic estimates and regional Moho depths obtained from the different gravimetric inversion methods for each of the different sub-areas have been summarized in Table. 6.
All in all, these comparisons show that our results are in most cases in the same range with seismic estimations in literature. By considering the effects of sediment and consolidated crust we provided satisfactory results for most of the continental crust where our estimates have a relatively good agreement with local seismic studies.
Larger discrepancies are present between seismic and gravimetric estimates in the offshore area around 22 degrees of latitude. This problem has to be further investigated also following the discussion presented in Eshagh et al. (2017). Table. 6. The regional seismic Moho depths in sub-areas of Iran and differences with gravimetric solutions in this area [km] 48 59 48 6.0 6.1 6.1 6.8 6.3 6.5

Summary and concluding remarks
In this study, we applied collocation method as well as two approaches based on isostasy principle presented by Sjöberg and Jeffrey for estimating the regional Moho in Iran using gravity observations. For this purpose, we considered the data of the GOCO03S satellite only global geopotential model that were subsequently reduced by topography/bathymetry, sediment and crystalline crust data effect by using the SRTM30_PLUS DTM and the CRUST1.0 model. The three different gravimetric approaches gave coherent estimates. The Sea. This can be the effect of the poor quality of the CRUST1.0 data in this region (see Eshagh et al. (2017) ).
The comparisons performed in this paper prove that further analyses are needed to come to a better consistency between gravity and seismic derived Moho depths in Iran before computing any joint seismic/gravimetric estimate.