Combined Gravimetric-Seismic Moho Model of Tibet

Substantial progress has been achieved over the last four decades to better understand a deep structure in the Himalayas and Tibet. Nevertheless, the remoteness of this part of the world still considerably limits the use of seismic data. A possible way to overcome this practical restriction partially is to use products from the Earth’s satellite observation systems. Global topographic data are provided by the Shuttle Radar Topography Mission (SRTM). Global gravitational models have been derived from observables delivered by the gravity-dedicated satellite missions, such as the Gravity Recovery and Climate Experiment (GRACE) and the Gravity field and steady-state Ocean Circulation Explorer (GOCE). Optimally, the topographic and gravity data should be combined with available results from tomographic surveys to interpret the lithospheric structure, including also a Moho relief. In this study, we use seismic, gravity, and topographic data to estimate the Moho depth under orogenic structures of the Himalayas and Tibet. The combined Moho model is computed based on solving the Vening Meinesz–Moritz (VMM) inverse problem of isostasy, while incorporating seismic data to constrain the gravimetric solution. The result of the combined gravimetric-seismic data analysis exhibits an anticipated more detailed structure of the Moho geometry when compared to the solution obtained merely from seismic data. This is especially evident over regions with sparse seismic data coverage. The newly-determined combined Moho model of Tibet shows a typical contrast between a thick crustal structure of orogenic formations compared to a thinner crust of continental basins. The Moho depth under most of the Himalayas and the Tibetan Plateau is typically within 60–70 km. The maximum Moho deepening of ~76 km occurs to the south of the Bangong-Nujiang suture under the Lhasa terrane. Local maxima of the Moho depth to ~74 km are also found beneath Taksha at the Karakoram fault. This Moho pattern generally agrees with the findings from existing gravimetric and seismic studies, but some inconsistencies are also identified and discussed in this study.


Introduction
Starting from the 1980s, numerous tomographic surveys have been carried out to investigate a deep structure in the Himalayas and Tibet [1][2][3][4][5][6]. For a more detailed overview of these studies, we refer readers to [7] and further in Section 2.2. Despite this effort, large parts of this region are still not yet sufficiently covered by high quality seismic data mainly due to its remoteness and extreme climate.
regional Moho models are briefly summarized in Section 2. Methods applied for this purpose are described in Section 3. Results are presented in Section 4 and discussed in Section 5. The study is concluded in Section 6.

Study Area
The gravity data analysis was conducted within the study area bounded by parallels 20 and 50 arc-deg of the northern latitude and meridians 60 and 110 arc-deg of the eastern longitude. It is noted that the seismic analysis was restricted by the availability of data (see Figure 1). The most prominent tectonic feature in this study area is an active continent-to-continent collision of the Indian plate with the Tibetan block, responsible for the uplift of the Himalayas. Tomographic evidence also suggests the southward subduction of the Eurasian lithosphere beneath the Tibetan block. A significant compressional tectonism due to these processes resulted in the uplift of the whole Tibetan Plateau. Major geological provinces include also continental basins (Indo-Ganges, Qaidam, Tarim, Sichuan) and orogens (Hindu Kush, Tien Shan). The geological setting of the study area is shown in Figure 1.

Seismic Data
For the Bengal lowland and the Shillong Plateau in the east of the Indian peninsula, we used receiver functions data from [37,38] and data from P and S waves prepared by [39]. According to these studies, the Moho depth under the Bengal lowland is 38-40 km, while it deepens to 48-50 km towards the collision zone between the Indian tectonic plate and the Tibetan block. For the western part of India, we used the receiver functions data from [40], showing that the Moho depth there is mostly within 36-44 km. For the central part of the Indian tectonic plate, we used the receiver functions data from [41][42][43][44][45]. The analysis of these data revealed a relatively complex Moho geometry with depths within a wide range from 35-48 km in the central part, deepening to 40-44 km towards

Seismic Data
For the Bengal lowland and the Shillong Plateau in the east of the Indian peninsula, we used receiver functions data from [37,38] and data from P and S waves prepared by [39]. According to these studies, the Moho depth under the Bengal lowland is 38-40 km, while it deepens to 48-50 km towards the collision zone between the Indian tectonic plate and the Tibetan block. For the western part of India, we used the receiver functions data from [40], showing that the Moho depth there is mostly within 36-44 km. For the central part of the Indian tectonic plate, we used the receiver functions data from [41][42][43][44][45]. The analysis of these data revealed a relatively complex Moho geometry with depths within a wide range from 35-48 km in the central part, deepening to 40-44 km towards north, and reaching 50-58 km at the beginning of the collision zone. Across the collision zone characterized also by a rising topography, the Moho deepens from 46 km to 56-58 km according to the analysis of receiver functions data from [5]. In Ladakh, Moho locally deepens to~80 km [46]. To the east, Moho deepens from 38 km for the Indian plate to 44-48 km under the Himalayan foothills [47]. The receiver functions data from [47] also revealed that further north, Moho deepens under the greater Himalayas and southern Tibet to 65-70 km. Beneath the Lhasa terrane, Moho deepens locally tõ 80 km. The analysis of receiver functions data from Pakistan also revealed a complex Moho geometry with depths ranging from 52 km in the Peshawar basin to 68-72 km in the higher Himalayas [47].
The seismic dataset for China comprises mainly the receiver functions data and the seismic profiles collected by [48][49][50][51]. Receiver functions calculated for data collected by the INDEPTH-IV seismic array [52] indicate an average Moho depth 65-70 km beneath north-eastern Tibet, 50-55 km for the Qaidam basin with an abrupt change (from 50-70 km) towards the Kunlun Mountains. The Moho topography detected from receiver functions data beneath western Tibet (ANTILOPE-1 profile) is rather uneven [53]. Moho there deepens from~50 km to 75-78 km under margins of the Himalayas and then shallowing to 66-70 km under the central Lhasa terrane. For the Qiangtang terrane, the Moho depth is typically~70 km, but abruptly shallows to 43-45 km at the beginning of the Tarim basin. The deep seismic sounding profiles from [54] exposed a sharp Moho depth change from 60 km under the east Kunlun to 48-50 km under the Zaige basin. For the north-eastern margin of the Tibetan Plateau, the receiver functions data from the ASCENT project [55] exhibited large Moho undulations within the range of 40-60 km. In the Qinling-Qilian block, Moho becomes shallower gradually from the west to the east. To the east of 105 • E, the average Moho depth is~45 km, whereas it exceeds 50 km to the west. Seismic data from [56][57][58][59][60][61], comprising the receiver functions data at 252 stations and the deep seismic sounding profiles, provide a detailed Moho map covering a broad region of the eastern part of the Tibetan Plateau and including also the Sichuan basin, the Ordos block, and the Yangtze Craton. According to these data, the Moho depth decreases from the eastern margin of the Tibetan Plateau (68-72 km) to the southern part of Yangtze Craton (28-30 km) and the Ordos block (38-42 km). Other regions are characterized by a complex Moho topography with depths: 52-56 km for the Qilian orogen, 40-44 km for the Qinling-Dabie orogen, 50-54 km for the Songpan-Ganzi terrane, 40-44 km for the Sichuan basin, 48-52 km for the Chuandian Plateau, and 28-36 km for the Yangtze Craton. According to results obtained from processing the receiver functions data by [62], the Moho depth under the north-east part of the Tibetan Plateau is typically within 50-52 km. Beneath the Liupan Mountains, Moho deepens to 52-54 km. Moho in the south-western part of the Ordos block is typicallỹ 50 km. To the south, the deep seismic sounding profiles from [63][64][65][66] indicate that Moho changes from 34-44 km.

Gravity Data
The isostatic gravity disturbances used in this study for a gravimetric Moho inversion were computed according to numerical steps, explained next.

Free-Air Gravity Data
The free-air gravity disturbances δg FA were computed from the EIGEN-6C4 [67] gravitational coefficients T n,m corrected for the GRS80 (Geodetic Reference System 1980) [68] normal gravity component using the following expression [69]: where GM is the geocentric gravitational constant, R is the Earth's mean radius, Y n,m are the surface spherical functions of degree n and order m, and n is the upper summation index of spherical harmonics. The 3D position in Equation (1) and thereafter is defined in the spherical coordinate system (r, Ω); where r is the radius and Ω = (ϕ, λ) is the spherical direction with the spherical latitude ϕ and longitude λ.

Bouguer Gravity Data
The Bouguer gravity disturbances δg B were obtained from the free-air gravity disturbances δg FA after applying the topographic g T and stripping gravity corrections due to density contrasts of the ocean (i.e., bathymetry) g B , the ice g I , sediments g S , and the consolidated crust g C . Hence, we write: The gravity corrections applied in Equation (2) were computed using the following generalized expression [70][71][72]: The potential coefficients V n,m are defined by: where ρ Earth is the Earth's mean density, and the coefficients {Fl The coefficients {L  (5) describe the geometry and density (or density contrast) distribution within a particular volumetric mass density layer.
Topographic and stripping gravity corrections due to the bathymetry and the ice were computed using the Earth2014 datasets [73]. A depth-dependent seawater density model was adopted in the definition of the ocean density contrast. For the reference crustal density of 2900 kg m −3 and the surface seawater density of 1027.91 kg m −3 [74], the nominal ocean density contrast (at zero depth) equals 1872.09 kg m −3 . The depth-density parameters (up to the second-order density term) were given in [75,76]. The glacial density of 917 kg m −3 [77] was adopted to compute the ice-stripping gravity correction. The stripping gravity corrections attributed to sediments and the consolidated crust were computed from the CRUST1.0 global seismic crustal model [78]. It is noted that the atmospheric gravity correction is completely negligible in the context of the gravimetric Moho modelling, having maxima less than 1 mGal [79].

Isostatic Gravity Data
Finally, the isostatic gravity disturbances δg I were computed from the Bouguer gravity disturbances δg B by applying the compensation attraction g VMM , so that: The compensation attraction g VMM in Equation (6) reads [34]: where G is the Newton's gravitational constant, ∆ρ c/m is the Moho density contrast, and the values of the Moho depth D were used from a new seismic model (Section 3.1).

Gravity Maps
Regional maps of the free-air, Bouguer, and isostatic gravity disturbances are shown in Figure 2 (for statistics, see Table 1). All gravity computations were realized on a 1 × 1 arc-deg grid with a spectral resolution complete to the spherical harmonic degree of 180. Values of the free-air gravity disturbances were mostly within ±50 mGal, except for gravity highs over the Himalayas, Tien Shan, and Hindu Kush. Gravity lows were along the Indo-Ganges basin (Figure 2a). Negative values of the free-air gravity disturbances also prevailed over the Tarim, Qaidam, and Sichuan basins. The Bouguer gravity map reveals the isostatic signature of major orogens, marked by large negative values (Figure 2b). The isostatic gravity disturbances were everywhere negative, indicating a systematic bias mainly due to unmolded mantle density heterogeneities (Figure 2c). where G is the Newton's gravitational constant, Δρ c/m is the Moho density contrast, and the values of the Moho depth D were used from a new seismic model (Section 3.1).

Gravity Maps
Regional maps of the free-air, Bouguer, and isostatic gravity disturbances are shown in Figure 2 (for statistics, see Table 1). All gravity computations were realized on a 1 × 1 arc-deg grid with a spectral resolution complete to the spherical harmonic degree of 180. Values of the free-air gravity disturbances were mostly within ±50 mGal, except for gravity highs over the Himalayas, Tien Shan, and Hindu Kush. Gravity lows were along the Indo-Ganges basin (Figure 2a). Negative values of the free-air gravity disturbances also prevailed over the Tarim, Qaidam, and Sichuan basins. The Bouguer gravity map reveals the isostatic signature of major orogens, marked by large negative values (Figure 2b). The isostatic gravity disturbances were everywhere negative, indicating a systematic bias mainly due to unmolded mantle density heterogeneities (Figure 2c).

Method
The methodology applied to compute the seismic and gravimetric Moho models and their combination is described in this section.

Seismic Moho Model
We used the seismic data to compile the Moho model by applying a processing strategy similar to that used for the construction of recent continental-scale crustal models, for instance, by [31]. For this purpose, we first inspected the quality of seismic data and then used them to generate the Moho contours by applying a standard kriging technique. This method was chosen from among others as a widely-applied geo-statistical technique. The idea behind this geo-statistical method is to reproduce trends that were estimated from combining the seismic data and the relief. The kriging parameters were chosen as follows: the interpolation area was 60 • -110 • E, 20 • -50 • N, grid step 1 degree, no anisotropy, with a linear variogram, while setting a scale factor equal to one. The linear variogram was intended to find a local vicinity of the observed point and for weighting the observed points used in the function interpolation at a given grid point.

Isostatic Moho Model
We used isostatic gravity disturbances δg I (Figure 2c) to determine the Moho depth by solving the VMM inverse problem of isostasy defined in the following generic form [35]: where dσ is the surface integration element and σ is the solid angle. The kernel K in Equation (8) is a function of the spherical angle ψ, and the parameter s = 1 − D/R is a function of the Moho depth D.
Its spectral form reads: where the Legendre polynomials P n are defined for the argument t = cos ψ. The expression in Equation (8) is the (non-linear) Fredholm integral equation of the first kind. Its direct solution (up to a second-order term) was derived in the following form [34]: The Moho term D 1 in Equation (10) was computed from the isostatic gravity coefficients δg i n,m as follows: The singularity for ψ → 0 in the third constituent on the right-hand side of Equation (10) was solved by applying the surface integration over the inner zone [34].

Combined Moho Model
We applied the method of [36] to constrain the gravimetric solution by seismic data. The principle of this method is to compute the non-isostatic gravity correction in order to account for differences between the gravimetric and seismic models. The non-isostatic gravity correction was then applied to the isostatic gravity disturbances. The resulting isostatic gravity disturbances obtained after applying the non-isostatic correction were then used to compute the combined Moho model according to the VMM isostatic model (Equations (10) and (11)).

Results
The seismic and gravimetric Moho models are shown in Figure 3 (for statistics, see Table 2). For the comparison, we also plotted the CRUST1.0 model. The final combined gravimetric-seismic Moho model is presented in Figure 4. Moho depth differences are plotted in Figure 5 (with the statistical summary given in Table 3).  Table 3. Statistics of the regional Moho depth differences.

Results
The seismic and gravimetric Moho models are shown in Figure 3 (for statistics, see Table 2). For the comparison, we also plotted the CRUST1.0 model. The final combined gravimetric-seismic Moho model is presented in Figure 4. Moho depth differences are plotted in Figure 5 (with the statistical summary given in Table 3).  Table 3. Statistics of the regional Moho depth differences.

Comparison of Results
Despite overall similarities between the Moho solutions presented in Figures 3 and 4, large localized differences existed (see Figure 5). As already stated, the seismic Moho model (Figure 3a The interpolation of the Moho geometry in regions with a sparse seismic data coverage by using the gravity information modified the Moho solution only slightly within ±2 km for most of Tibet (Figure 5a). The corresponding modifications of the Moho geometry up to 8 km were, however, seen in the western and eastern parts of the Himalayas. Interestingly, the combined Moho model had a slightly better RMS fit with the CRUST1.0 than our gravimetric and seismic models (cf. Table 3). The interpolation of the Moho geometry in regions with a sparse seismic data coverage by using the gravity information modified the Moho solution only slightly within ±2 km for most of Tibet (Figure 5a). The corresponding modifications of the Moho geometry up to 8 km were, however, seen in the western and eastern parts of the Himalayas. Interestingly, the combined Moho model had a slightly better RMS fit with the CRUST1.0 than our gravimetric and seismic models (cf. Table 3).

Discussion
Our result based on the combined processing of gravity and seismic data (see Figure 4) revealed that the maximum Moho deepening of ~76 km occurred to the south of the Bangong-Nujiang suture under the Lhasa terrane. Deep Moho (~74 km) was also detected under Taksha at the Karakoram fault. This finding closely agrees with the results from teleseismic receiver function analysis of seismograms recorded on a ~700 km-long profile of 17 broadband seismographs traversing the north-west Himalayas conducted by [5]. They reported a progressive northward Moho deepening from ~40 km beneath Delhi south of the Himalayan foredeep to ~75 km beneath Taksha at the Karakoram fault. They also proposed that the Indian lithosphere might penetrate as far as the Bangong-Nujiang suture and possibly as far north as the Altyn Tagh. Alternative theories facilitated the hypothesis of crustal shortening and consequent crustal thickening attributed to the extrusion or escape tectonics mechanism [80]). According to these theories, the motion of the Indian plate pressed the Indochina block, and a proposed mechanism is that a large part of the crustal shortening was accommodated by thrusting and folding of the sediments of the passive Indian margin together with the deformation of the Tibetan crust [81].
A Moho depth largely between 60-70 km beneath most of the Himalayas and Tibet was reported in [7]. These values closely agree with our findings. Their result, however, indicated that the maximum Moho deepening of ~79 km occurred in the northern Tibet and the Himalayas. A maximum Moho deepening in northern Tibet along margins with the Tarim basin was also reported by [18]. According to their estimates, the Moho was more than 65 km deep under most of Tibet with maxima of ~82 km in its western part. A different Moho topography was presented by [21]. According to their result, the Moho depth in the west was deeper than that in the east, and the deepest Moho of ~77 km occurred beneath the western Qiangtang block. They also identified a Moho offset of approximately 5 km beneath the Yarlung-Zangbo suture at the juncture between the Himalayan and Lhasa blocks. Our combined Moho model also differed from the result presented by [14]. According to their result, the Moho depth under most of Tibet was largely within 70-75 km with the maximum deepening to ~80 km along margins of the plateau and a shallower depth of 65 km under the Bangong-Nujiang suture in central Tibet.

Accuracy Assessment
As follows from the discussion in the preceding paragraph, large inconsistencies exist between Moho models. Despite the assessment of accuracy of these models being generally not simple, two principal factors should be taken into consideration in this respect. These involve the quality of seismic and gravity data, as well as available information about the crustal density structure.

Discussion
Our result based on the combined processing of gravity and seismic data (see Figure 4) revealed that the maximum Moho deepening of~76 km occurred to the south of the Bangong-Nujiang suture under the Lhasa terrane. Deep Moho (~74 km) was also detected under Taksha at the Karakoram fault. This finding closely agrees with the results from teleseismic receiver function analysis of seismograms recorded on a~700 km-long profile of 17 broadband seismographs traversing the north-west Himalayas conducted by [5]. They reported a progressive northward Moho deepening from~40 km beneath Delhi south of the Himalayan foredeep to~75 km beneath Taksha at the Karakoram fault. They also proposed that the Indian lithosphere might penetrate as far as the Bangong-Nujiang suture and possibly as far north as the Altyn Tagh. Alternative theories facilitated the hypothesis of crustal shortening and consequent crustal thickening attributed to the extrusion or escape tectonics mechanism [80]). According to these theories, the motion of the Indian plate pressed the Indochina block, and a proposed mechanism is that a large part of the crustal shortening was accommodated by thrusting and folding of the sediments of the passive Indian margin together with the deformation of the Tibetan crust [81].
A Moho depth largely between 60-70 km beneath most of the Himalayas and Tibet was reported in [7]. These values closely agree with our findings. Their result, however, indicated that the maximum Moho deepening of~79 km occurred in the northern Tibet and the Himalayas. A maximum Moho deepening in northern Tibet along margins with the Tarim basin was also reported by [18]. According to their estimates, the Moho was more than 65 km deep under most of Tibet with maxima of~82 km in its western part. A different Moho topography was presented by [21]. According to their result, the Moho depth in the west was deeper than that in the east, and the deepest Moho of~77 km occurred beneath the western Qiangtang block. They also identified a Moho offset of approximately 5 km beneath the Yarlung-Zangbo suture at the juncture between the Himalayan and Lhasa blocks. Our combined Moho model also differed from the result presented by [14]. According to their result, the Moho depth under most of Tibet was largely within 70-75 km with the maximum deepening to~80 km along margins of the plateau and a shallower depth of 65 km under the Bangong-Nujiang suture in central Tibet.

Accuracy Assessment
As follows from the discussion in the preceding paragraph, large inconsistencies exist between Moho models. Despite the assessment of accuracy of these models being generally not simple, two principal factors should be taken into consideration in this respect. These involve the quality of seismic and gravity data, as well as available information about the crustal density structure.
The best horizontal Moho resolution is typically inferred from reflection profiles [82], but such a technique is expensive, thus not widely used. The Moho detection from a two-way travel time might also be affected by a weak reflectivity. An intermediate spatial resolution could be obtained from a deep seismic sounding based on using refracted and wide-angle reflected waves. Uncertainties of detecting the Moho depth from the wide-angle reflection and refraction methods are typically 1-2 km. Another technique, which has become quite common during the last two decades, is based on inverting the P-or S-wave receiver functions [83,84]. The estimated uncertainty of this method is about 3 km. Intermediate-period surface waves are quite sensitive to crustal thickness, but are not able to discriminate it from the mantle velocity structure, thus they cannot be inverted uniquely [85,86].
To interpolate the Moho information over regions where seismic data are sparse or missing, we used the gravity data and information about the crustal density structure. Moho depth uncertainties attributed to errors in gravity data are relatively small, because the accuracy of the latest global gravitational models of about ±10 mGal or better is expected globally at a resolution of about 100 km [87,88]. Errors in topographic models are completely negligible in the context of a gravimetric Moho recovery. However, relative errors in the estimated Moho depth of 10% or more could be expected due to uncertainties within the CRUST1.0 sediment and consolidated crustal density layers [31].

Conclusions
We have used seismic, gravity, and topographic data to estimate the Moho depth of the orogenic formations of the Himalayas, Tibet, Hindu Kush, and Tien Shan, including also the surrounding continental basins of Indo-Ganges, Tarim, Sichaun, and Quadam. According to our combined gravimetric-seismic model, the Moho depth under the Himalayas and the Tibetan Plateau is mostly within 60-70 km. The maximum Moho deepening of~76 km occurs to the south of the Bangong-Nujiang suture under the Lhasa terrane. The additional significant Moho deepening tõ 74 km was detected beneath Taksha at the Karakoram fault. A large Moho depth of 60-65 km was also confirmed beneath the Hindu Kush and Tien Shan. In contrast to a large continental crustal thickness of these orogenic formations, the Tarim, Sichuan, and Indo-Ganges basins are characterized by a Moho depth typically within 38-47 km. Moho in the Quadam basin slightly exceeds the depth of 50 km. The mean Moho depth within the whole study area is estimated to be 46.5 km.
Reviewing seismic and gravimetric studies, we identified some significant discrepancies between existing Moho models of Tibet. Based on our assumptions, we expect that errors in estimated values of the Moho depth from the seismic data analysis are typically less than 5 km, but could reach as much as 10-15 km when interpolated using the gravity information. These large errors are mainly due to the fact that a Moho recovery from the gravity data requires the knowledge of crustal density structure. The accuracy of currently available global and regional crustal/lithospheric models is obviously restricted, with expected relative uncertainties of 10% or more, especially in geologically-complex regions. Another problem is that isostatic theories could not realistically reproduce a Moho geometry, especially along active tectonic margins. All these factors obviously affect the accuracy of the combined Moho model.

Conflicts of Interest:
The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; nor in the decision to publish the results.