Middle East Moho Relief Estimation Using Spherical Prisms in Gravity Inversion

Study of Moho in Middle East and surrounding region is of great importance for geoscientists. The area contains parts of the Eurasian, Indian, African and Arabian plates as well as active tectonic structures which made it rich in geodynamic and tectonic interests. In this study we investigated Moho relief in the Middle East region using gravity and seismic data. Regarding the extent of the study area, spherical prism (tesseroid) modeling is used to calculate the gravity effect of the topography and crustal anomalies. Determining of Moho depth from gravity data is a nonlinear inverse problem. We inverted the gravity data using Uieda’s inversion method where the process was constrained by the available seismic data over the study region. The effect of topography and crustal sediments were excluded using global topography and crustal models. The resulting Moho relief is in accordance with plate boundaries and correlates well with the prominent tectonic features of the Middle East region. According to the results, the thinnest part of the crust found to be about 12 km in the Indian Ocean and the thickest part appeared in the west of Tibetan plateau with depths of about 54 km. In some parts of the study area discrepancies were seen between our results and Moho depths from seismological studies these differences are most probably caused by different approaches used in the different studies. Since we imposed smoothness by regularization on estimated Moho map, this can also be additional source of the discrepancies.


Introduction
The Middle East region extends from eastern Mediterranean to the Arabian Sea and has a rich geological history that stretches back to the Proterozoic [Mart, 2013]. This area is a collection of different geodynamic and tectonic structures, each one with its own complex nature, and of great importance for scientific studies. This region contains parts of the Eurasian, Indian, African and Arabian plates as the large plates and some small plates [Bird, 2003]. Figure 1 shows prominent geodynamics of the Middle East. In this region there are different types of boundary interactions like oceanic divergent boundary in the Red Sea, continental divergent boundary in East African rift, oceanic-continental subduction in Makran, Zagros thrust belt, north and east Anatolian fault, just to mention a few.
One of the most significant tectonic features for study is the Mohorovic discontinuity (Moho) which is the boundary surface between the crust and the mantle. Moho surface characteristics, including the depth, density contrast between crust and upper mantle, shear wave velocity, just to name a few, were evaluated all over the world by gravimetric, thermal, and seismic methods like CRUST1.0 [Laske et al., 2013], LITHO1.0 [Pasyanos et al., 2014] and GEMMA [Reguzzoni et al., 2015]. Over the study area there are vast amount of published studies such as: East African rift system and Afar Triangle related researches [e.g. Dugda et al., 2005;Tiberi et al., 2007;Hammond et al., 2011;Seber et al., 2001;Almadani, 2011;Globig et al., 2016]. Moho depth and crustal studies have been carried out in Red Sea and Gulf of Aden [e.g. Tiberi et al., 2005;Hansen et al., 2007;Salem et al., 2013]. There are also studies in Arabian plate [e.g. Al-Hashmi et al., 2011;Al-Damegh et al., 2005;Mechie et al., 2013].
There are also variety of studies in Anatolian, Aegean and Mediterranean region [e.g. Saunders et al., 1998;Li et al., 2003;Van der Meijde et al., 2003;Starostenke et al., 2004;Zhu et al., 2006;Tezel et al., 2013;Karabulut et al., 2013;Sampietro et al., 2018 andKarabulut et al., 2019]. Parastoo Jalooli et al. Presence of Global Geopotential Models in the satellite era, made the access to gravimetric regular grid data feasible all over the earth. In gravimetric studies over the large area (similar to the area of interest in this study) using rectangular prisms will not take into account the curvature of the earth. Since the study region is larger than 15° × 15°, tesseroids will produce more accurate results in comparison with prisms in an acceptable computation time [Wild-Pfeiffer, 2008]. In this paper we aim to determine continues model for Moho in the Middle East area by isolating the gravitational effect of Moho relief and then inverting data using tesseroids. Existing Moho depths from seismological studies over the study area will be used to validate the results.

Methodology
To estimating the relief of the Moho from gravity data, first one must obtain the gravitational effect of the target anomalous density distribution attributed to the Moho relief, and then apply the appropriate inverse method. This requires eliminating all gravity effects other than that of the target anomalous density from observed data. Since real Moho undulates around a reference Moho (Moho of the Normal Earth, an ellipsoidal reference Earth) the first correction is to remove the gravity of an ellipsoidal reference Earth (γ). The gravity disturbance is determined by subtracting the normal gravity γ(P) of reference ellipsoid WGS84 from the observed gravity at the same point "P" on the earth: The obtained disturbance , as shown in Figure 2, contains all masses attributed to the topography as well as the mass deficiency corresponding to the oceans and sedimentary basins, crustal sources, heterogeneities below the upper mantle, and the effect of the difference between the real Moho topography and the Moho of the Normal Earth . Next, total Bouguer effect at point "P" is obtained by eliminating the gravity effect of topography and oceanic masses from gravity disturbance ( ), using ETOPO1 digital train models [Amante and Ekins, 2009]: To remove the gravitational effect of known sedimentary basins (the effects of other crustal and mantle sources, shown in (

=sediments (3)
Gravitational effects of above mentioned reductions and corrections are calculated by forward modeling in a spherical Earth approximation using tesseroids. Each utilized tesseroid is a spherical mass element bounded by two meridians, two parallels and two co-center spheres. There are two approaches for calculating gravitational effect of tesseroids. One is based on Taylor series expansion [Hech and Seitz, 2007], and the other is based on Gauss-Legendre integration [Asgharzadeh et al., 2007]. In this paper we used the modified version of the second approach [Uieda et al., 2015 andUieda, 2013]. The residuals of all applied reduction and correction steps represents the gravity effect of the anomalous Moho relief (Figure 2-e) that will serve as an input data for inversion procedure.

Discretization of problem and inversion
In order to perform forward modeling to calculate the gravity attraction of the anomalous Moho relief during the inversion procedure, one needs to discretize the real Moho into N × M spherical prisms (Figure 2-f).
Here we used the inversion method introduced by Uieda and Barbosa [2016]. They merged the effective Bott's [Bott, 1960] method with smoothness regularization and a discretization of the anomalous Moho into tesseroids (spherical prisms). Since in this method the absolute value of the density-contrasts of the tesseroids is a fixed parameter, then Moho gravity anomaly will be a nonlinear function of tesseroids depth: where d i is the ith element of N-dimensional predicted data vector, k is M-dimensional parameter vector including M Moho depths (h p ), f i is the ith nonlinear function that maps parameters onto data. These functions are radial component of the gravitational attraction of the tesseroids.
We estimate the parameter vector p from a set of observed gravity data d o in a least-squares sense which minimizes the data-misfit function or fidelity term (p): Since estimation of the Moho depths from a set of observed gravity data is an ill-posed non-linear inversion problem (Silva et al. 2001), however under some assumptions the uniqueness of the solution is guaranteed [Sampietro and Sanso, 2012]. We needed to perform a regularized least-squares non-linear inversion in this study.
First-order Tikhonov regularization [Tikhonov and Arsenin, 1977] is used to impose smoothness on the solution.
Where is regularization parameter and θ is the smoothness penalty functional which is given by: Here R is an l × m finite-difference matrix representing l first order differences between the depths of adjacent tesseroids. The goal function (p) is also nonlinear with respect to p and can be minimized using the Gauss-Newton method.
Apart from Moho depths, we need to estimate three hyper parameters  namely the regularization parameter ( ), Moho reference level (h n ) and density contrast (Δ ). They have been estimated in two steps during the inversion by holdout-cross validation method [Kim, 2009]. The holdout approach first divides the observed data into two groups of training and testing sets. Training set is used for inversion procedure while testing set is used for evaluating the quality of above mentioned parameters. Mean Square Error (MSE) is also used to select inversion parameters optimally. The inversion was performed using hyperparameters: reference level (h n = 32.5 km), density contrast (Δ = 500 / ³ ) and regularization parameter ( = 1 -10), determined during cross validation step.

Data Sources
The gravity data used in this study are generated from satellite only spherical harmonic model EIGEN-6S4-V2 [Förste and Bruinsma, 2016] up to degree and order 300 (the model combines data from GOCE, GRACE and LAGOS satellites). The data were downloaded from the International Center for Global Earth Models (ICGEM) [Barthelmes and Kohler, 2012] in the form of the complete gravity field on a regular grid with 0.2° grid spacing at ellipsoidal height of 50 km.
The gravitational effect of topography and sediment basin was removed using spherical Earth approximation via utilizing spherical prisms. The effects are calculated using ETOPO1 digital terrain model [Amante and Eakins, 2009], CRUST1.0 [Laske et al., 2013] and standard densities of 2670 / ³ for continents and -1630 / ³ for oceans. The sediment-free Bouguer disturbance is obtained after subtracting the gravitational effect of topography and sediment basin from gravity disturbance by ignoring all other crustal sources.
The seismic information extracted from different seismological studies over the Middle East region including 448 Moho depths are shown in Table 1. The estimated Moho depths from our proposed method later on are validated by these seismic data. Figure 3 shows the sediment-free Bouguer disturbance and location of seismic Moho depths.

Results and Discussion
As mentioned earlier the study area, Middle East region, contains significant plate boundaries such as The African-Arabian divergent boundary located in the red sea, Arabian-Eurasian convergent boundary forming Zagros Mountains belt, African-Anatolian and Aegean Sea subduction zone in the Mediterranean Sea, African-Somalian divergent boundary, Arabian-Indian convergent and Indian-Somalian divergent boundaries in the Indian ocean. In this Section we concentrate on the observed correlations between above mentioned features and the Moho discontinuity map obtained from current study. The estimated Middle East Moho map is shown in Figure 4, where we adopted Peter Bird's digital model [Bird, 2003] Radjaee et al. [2010], that is, if there is any crustal root for Alborz Mountain, it is inadequate to say this region is in isostatic compensation, but still needs more research to be clarified. Parastoo Jalooli et al.
8 Figure 5. Differences between Moho depths from seismological studies and current study.

Conclusion
In this study the Middle East Moho map was determined with an efficient inversion method and considering spherical earth by replacing rectangular prisms with tesseroids. The resulted gravity residuals and correlation of the Moho depths with geological and tectonic structures of the study region show the performance of the inversion procedure that adapt in this study. The differences between Moho depths estimated in this investigation and those determined by seismological studies are mostly due to smoothness we imposed by regularization on estimated Moho map and different approaches used to determine the Moho depths. The realistic Moho depth model from employing of detailed sediment description will fluctuate around the Moho depth model obtained in current study.
The main cause of such fluctuation arises from employing of CRUST 1.0 model and imposing smoothness by Tikhonov regularization on solution of the inversion in this study.
In addition, we ignored crustal anomalies and considered certain density values for the oceans and continents.
The smallest discrepancy happened in Arabian plate, Anatolian and Aegean Sea plates and the maximum discrepancy was seen in some parts of Zagros and Alborz mountains, as well as in some parts of the African rift.
The estimated Moho map also showed a good correlation with topography of the region, this indicates the acceptable performance of the adapted inversion procedure and utilizing of spherical prisms. This correlation could be interpreted as the fact that almost the whole region is in an isostatic compensation, but there are areas that do not show acceptable correlation like Alborz Mountain. The estimated Moho depth for the study area ranges from 12 km in the deepest parts of Indian Ocean up to 54 km in profound parts of the Tibetan plateau.