Crustal and Upper Mantle Density Structure Beneath the Qinghai-Tibet Plateau and Surrounding Areas Derived from EGM2008 Geoid Anomalies

Abstract

As the most active plateau on the Earth, the Qinghai-Tibet Plateau (TP) has a complex crust–mantle structure. Knowledge of the distribution of such a structure provides information for understanding the underlying geodynamic processes. We obtain a three-dimensional model of the density of the crust and the upper mantle beneath the TP and surrounding areas from height anomalies using the Earth Gravitational Model 2008 (EGM2008). We refine the estimated density in the model iteratively using an initial density contrast model. We confirm that the EGM2008 products can be used to constrain the crust–mantle density structures. Our major findings are: (1) At a depth of 300–400 km, high-D(ensity) anomalies terminate around the Jinsha River Suture (JRS) in the central TP, which suggests that the Indian Plate has reached across the Bangong Nujiang Suture (BNS) and almost reaches the JRS. (2) On the eastern TP, low-D(ensity) anomalies at a depth of 0–300 km and with high-D anomalies at 400–670 km further verified the current eastward subduction of the Indian Plate. The ongoing subduction process provides force that results in frequent earthquakes and volcanoes. (3) At a depth of 600 km, low-D anomalies inside the TP illustrate the presence of hot weak material beneath it, which contribute to the inward thrusting of external material.

1. Introduction

The Qinghai-Tibet Plateau (TP) is the highest and largest plateau in the world. Its average height is approximately 4500 m, and its area is approximately 2,500,000 km² [1]. As a response to the subduction of the lithospheric mantle of the Indian Plate, large-scale motion of the crust occurs in this area. It has a complex crust–mantle structure [2]. Its average crust thickness is 65–70 km [3]. Knowledge of the structures of the crust and upper mantle beneath the TP is particularly informative for understanding the underlying formation. Different models of the formation of the TP have been developed [4,5,6,7]. Some geophysical constraints regarding the velocity structure [8,9,10,11], the electrical structure [12,13] and the thermal structure [14,15] of this region have been obtained. However, few results relating to the density structure of the whole plateau have been reported. Considering the fact that the Earth’s density is the direct product of its geodynamics, a detailed density structure can provide strong geophysical evidence for the underlying geodynamic processes.

Next to gravity, the geoid is another physical quantity that is sensitive to the Earth’s interior density. Both of them play important roles in studies of the Earth’s interior density structure. The gravity of one point is inversely proportional to the square of the distance between that point and the source of the structure. As the depth increases, the amplitude of gravity decreases rapidly, and the gravitational anomaly primarily reflects shallow short-wavelength density changes. However, another measure—the geoid—has a larger difference than gravity does. It is inversely proportional to the distance from the source. It reflects the long-wavelength field characteristics of the density changes deep in the earth. Therefore, the geoid anomaly is often used to study deep mantle density anomalies [16].

The absence of underground gravity observations in this region inhibits the progress of density research. However, as various observation methods have developed [17,18], sensitive, high-resolution, global satellite measurements of gravity have emerged to provide an effective means for the study of the earth’s interior density structure. The Earth Gravitational Model 2008 (EGM2008) [19] has been publicly released by the U.S. National Geospatial-Intelligence Agency (NGA). This gravitational model is complete in terms of spherical harmonics to degree and order 2159 and contains additional coefficients extending to degree 2190 and order 2159. It has an accuracy of 20 cm on the mainland, 12 cm in east-central China, 9 cm in northern China and 24 cm in western China at sea level [20]. In our study, we choose to use the geoid anomaly derived from the EGM2008 to invert the density.

A method of gravitational modeling based on a spherical prism was proposed to take the curvature of the earth into account. In this study, the objective function is determined using the Lagrange multiplier method and Tikhonov regularization theory [16]. The damped least-squares method is used during the inversion to obtain unbiased results with minimal variance [21]. Seismic velocities are imposed as constraints. Finally, a three-dimensional model of the density of the crust and the upper mantle beneath the TP and the surrounding areas is derived by inverting a set of geoid anomalies. This paper concludes with a geodynamic analysis.

2. Geology

The study region is defined by coordinates 75°E and 115°E and 20°N and 50°N, as shown in Figure 1. Several studies regarding the overall structural features and physical characteristics of the upper mantle under the TP have been conducted in the past few years. The crust is approximately 70–75 km thick under southern (S) Tibet and 60–65 km thick under northern (N), northeastern (NE) and southeastern (SE) Tibet [22,23]. The structures of the crust and the upper mantle in Tibet have obvious zonal characteristics and inhomogeneous distributions of the lower velocity in the crust and the higher velocity in upper mantle [24]. Weak zones occur in the upper crust in the Himalayas, in the middle crust in southern Lhasa, in the middle and lower crust in northern Lhasa, and in the middle and lower crust in Qiangtang and Songpan-Ganzi [25,26,27]. The Indian crust collides with the Tarim Basin at 80°E and reaches the Bangong Nujiang Suture (BNS) belt at 88°E (in the west at 90°E) [11]. The northern TP is characterized by a sharp Moho offset beneath the Altun Tagh Fault (ATF), a steep plateau margin, and a major south-directed thrust fault. Similarly, the characteristics of the southeastern TP are a low-gradient topography change; a crustal lower velocity zone (LVZ) that extends for a long distance; and a gradual Moho change at the plateau’s margin [25]. The crust/mantle interface beneath Tibet is anisotropic [2]. However, there are few results that describe the dedicated structure, such as quantitative descriptions of the magnitude and range of anomalies. The objective of our study is to develop a detailed 3-Dimensional model of the density structure of the crust and upper mantle below the TP.

3. Methodology

3.1. Fundamental Equations

According to the famous Bruns formula [29], which describes the relationship between the geoid height and the perturbing potential, at any point outside the Earth, the induced geoid height is described by $N\left(P\right)=\frac{T\left(P\right)}{\gamma \left(P\right)}$, where $P\left(r,\phi ,\lambda \right)$ is the computation point; $r$, $\phi$ and $\lambda$ are the radius, latitude, and longitude, respectively; and $\gamma \left(P\right)$ is the normal gravity on the geodetic reference ellipsoid. The perturbing potential, $T\left(P\right)$, is given by numerical integration over a set of spherical prisms, the so-called tesseroids, which are bounded by two meridians, two parallels, and two concentric circles. Therefore, the approximate expression for the geoid height developed by Heck and Size [30] based on tesseroids is shown in Equation (1),

$$N\left(P\right)=\frac{G\rho \Delta r\Delta \phi \Delta \lambda }{\gamma \left(P\right)}\left[K_{000}+\frac{1}{24}\left(K_{200}\Delta r^2+K_{020}\Delta {\phi }^2+K_{002}\Delta {\lambda }^2\right)+O\left({\Delta }^4\right)\right]$$

(1)

where $G$ is Newton’s gravitational constant; $\rho$ is the density contrast, which is assumed to be constant within each tesseroid; $\Delta r=r_1-r_2$, $\Delta \phi ={\phi }_1-{\phi }_2$, and $\Delta \lambda ={\lambda }_1-{\lambda }_2$ are the tesseroid dimensions in the radial and horizontal directions, respectively; $K_{000}$, $K_{002}$, $K_{020}$, and $K_{200}$ are coefficients that depend on the relative positions of the computation point and the geometric center of the tesseroids, which can be obtained in the article; and $O\left({\Delta }^4\right)$ is the Landau symbol, which indicates omitted fourth-order terms in $\Delta r$, $\Delta \phi$, $\Delta \lambda$.

Because the derivative of Equation (1) is linear with respect to the density contrast, $\rho$, a set of $n$ calculated geoid anomalies caused by a model composed of $m$ tesseroids can be written as

$$N_{\rho (n\times 1)}=G_{\rho \left(n\times m\right)}{\rho }_{\left(m\times 1\right)}$$

(2)

where $\rho$ is the vector that represents the estimated density contrast; $N_{\rho }$ is the vector that represents the geoid anomaly; and $G_{\rho }$ is the matrix of the sensitivity function relating the density contrast and the geoid anomaly.

3.2. Density Inversion

To derive the density contrast from gravitational observations, we established the objective function in Equation (4) based on the theory of Lagrange multipliers, the method which has been used to invert the geoid anomalies in Yellowstone Province [16]. This theory considers both the error in the observations and the error in the ridge regression of the parameter. We need to minimize the objective function as follows:

$$L\left(\rho ,\alpha \right)={\left(\rho -{\rho }_{apr}\right)}^T{W}_{\rho }^T{W}_{\rho }\left(\rho -{\rho }_{apr}\right)+2{\alpha }^T{W}_d\left(G_{\rho }\rho -N_{\rho }\right)$$

(3)

where the first term is the 2-norm of the fitting error, the condition term is the constraint, ${\rho }_{apr}$ is an initial model, and $W_{\rho }$ is the parameter-weighting matrix because the geoid anomaly is inversely proportional to the distance from the source. To ensure that each tesseroid has the same probability according to the matrix $G_{\rho }$, we use the matrix $W_{\rho }=diag\left(\frac{1}{{\left(z+\epsilon \right)}^{\beta }}\right)$ to calculate the weight in the matrix $G_{\rho }$, where $z$ is the depth of each tesseroid, $\epsilon$ is a small constant to avoid singularities, $\beta$ is a number that changes the weight of $W_{\rho }$; $W_d$ is the data weight matrix; and $\alpha$ is the vector of Lagrange multipliers. To minimize the objective function, we differentiate Equation (3) with respect to $\rho$ and $\alpha$; by setting them equal to zero, we obtain the following solution:

$$\{\begin{array}{l}\widehat{\rho }={\rho }_{apr}+{\left(W_{\rho }^T{W}_{\rho }\right)}^{-1}{G}_{\rho }^T{W}_d^T{\left(W_d{G}_{\rho }{\left(W_{\rho }^T{W}_{\rho }\right)}^{-1}{G}_{\rho }^T{W}_d^T\right)}^{-1}{W}_d\left(N^{\rho }-G_{\rho }{\rho }_{apr}\right)\\ \alpha =-{\left(W_d{G}_{\rho }{\left(W_{\rho }^T{W}_{\rho }\right)}^{-1}{G}_{\rho }^T{W}_d^T\right)}^{-1}{W}_d\left(N^{\rho }-G_{\rho }{\rho }_{apr}\right)\end{array}$$

(4)

where $\widehat{\rho }$ is the estimated density contrast. To solve the inversion problem, which is ill-posed because of the small singular values in the sensitivity matrix, $G_{\rho }$, we add the constant $\mu$ to Equation (4),

$$\widehat{\rho }={\rho }_{apr}+{\left(W_{\rho }^T{W}_{\rho }\right)}^{-1}{G}^T{W}_d^T{\left(W_{d}G{\left(W_{\rho }^T{W}_{\rho }\right)}^{-1}{G}^T{W}_d^T+\mu I\right)}^{-1}{W}_d\left(N_{\rho }-G{\rho }_{apr}\right)$$

(5)

where $\mu$ is the regularization parameter, which can be determined using Tikhonov regularization theory [31].

To obtain the best solution to Equation (5), we consider solving it iteratively as follows:

$${\widehat{\rho }}_{k+1}={\widehat{\rho }}_k+\Delta {\widehat{\rho }}_k$$

(6)

where ${\widehat{\rho }}_k$ is the density contrast after the kth iteration; for $k=0,1,2,\cdots ,n$, where ${\widehat{\rho }}_0$ is ${\rho }_{apr}$, which is an induced initial density model given by a 3-Dimension Vs model; and $\Delta {\widehat{\rho }}_k$ is as follows:

$$\Delta {\widehat{\rho }}_k={\left(W_{\rho }^T{W}_{\rho }\right)}^{-1}{G}^T{W}_d^T{\left(W_{d}G{\left(W_{\rho }^T{W}_{\rho }\right)}^{-1}{G}^T{W}_d^T+\mu I\right)}^{-1}{W}_d\left(N_{\rho }-G{\widehat{\rho }}_k\right)$$

(7)

Then, the process is to search for the solution that fits the observed data best. As an acceptance criterion, we use the root-mean-square (RMS) error of the fitting error function, i.e., $E_1=N_{\rho }-G_{\rho }{\widehat{\rho }}_k$ must be less than or equal to a constant value, ${\epsilon }_1$. Another criterion for convergence is that $E_2=\Delta {\widehat{\rho }}_k/{\widehat{\rho }}_k$ is smaller than a certain previously determined empirical value, ${\epsilon }_2$. The inversion may also be terminated after a pre-specified number of iterations.

4. Data Processing

The three-dimensional density model of the crust and the uppermost mantle is derived by inverting the geoid anomalies from the EGM2008 and expanding them to degree 180 using a starting density model derived from a seismic S-wave model.

We calculate the geoid height data at sea level using the spherical harmonic expansion to degree and order N = 180 from the EGM2008 (Figure 2i). The geoid model, which has an RMS accuracy of 0.15 m worldwide, is computed with respect to the WGS-84 [19]. It provides a uniform and global field, which enables interpretation on both regional and global scales.

However, the anomaly is generated by all sources within the earth. It represents the responses to all the interface undulations and density inhomogeneities below the earth’s surface. With the aim of delineating the density variations in the crust and upper mantle density, the signals of topographic masses above sea level, the effects of variations in density interfaces (the lithosphere, Moho, and sedimentary interfaces), and the influence of density changes below the upper mantle should be reduced before the inversion.

$$n_{\rho }=N_{\rho }-N_{\rho }^1-N_{\rho }^2-N_{\rho }^3-N_{\rho }^4-N_{\rho }^5$$

(8)

where $N_{\rho }$ is the vector of the observed geoid anomaly; $N_{\rho }^1$ is the vector of density changes below the upper mantle’s geoid effect; $N_{\rho }^2$ is the vector of the topographic geoid effect; $N_{\rho }^3$ is the vector of the geoid effect due to variations in the sedimentary interface; $N_{\rho }^4$ is the vector of the Moho interface variation effect; $N_{\rho }^5$ is the vector of the effect of lithosphere interface variations on the geoid; and $n_{\rho }$ is the vector of the residual geoid anomalies caused by variations in the upper mantle density.

The influence of density changes below the upper mantle (Figure 2j) is given by order and degree 2–6 of EGM2008 based on the point mass source theory developed by Bowin [32], i.e., $z_n=R/n-1$, where $z_n$ is the depth of the source; $R$ is the radius of the Earth; and $n$ is the order and degree of the spherical harmonic coefficient. The topographic effects (Figure 2e) are estimated using the 1 arc min global relief model of the earth’s surface, ETOPO1 [33] (Figure 2a), with a set of spherical prism discretizations based on Equation (2). The layer-specific density values are 2670 kg/m³ and 1640 kg/m³ for rocks above and below the surface, respectively. The same strategy is applied to the collected 0.5° × 0.5° sedimentary undulations (Figure 2c), where the average sediment depth is 4 km and the density contrast between sediment and crust is 200 kg/cm³ [34].

We then obtain the 1° × 1° Moho interface (Figure 2b) included in the global crustal model, crust1.0 [35], and the 1° × 1° lithospheric interface (Figure 2d) from the global lithospheric model, litho1.0 [36], respectively. The average depth of the Moho interface is 35 km, the density contrast between the crust and the uppermost mantle is 420 kg/cm³, the average depth of lithosphere is 100 km and the density contrast between the crust and the uppermost mantle is 10 kg/cm³ [34]. We use the strategy described above.

The effects are illustrated in Figure 2e–h. The contribution of the topographic masses is the largest with an amplitude of approximately 650 m. The effect of the amplitude of the basins is smaller (approximately 570 m) than that of the topography. The largest lithospheric effect is in the southwestern part of the study area and nearly 30 m. The sedimentary effects in the Sichuan, Qaidam and Indian Basins are obvious. Their amplitudes can be up to 25 m.

The residual geoid anomaly is shown in Figure 2k. It is computed after the above-motioned contributions have been reduced, and its value is ±55 m. The residual geoid exhibits a “positive-negative-positive” trend from southwest to northeast.

5. Results

5.1. Synthetic Tests

We constructed some synthetic models based on geoid-induced terms in our designed density model for testing our above-proposed method. We simulated three synthetic models (see Figure 3a,c,e). M1 is a synthetic model of a negative geoid anomaly caused by one low density (low-D) rectangular prism ranging from 150 km to 200 km in depth ($\Delta r$ = 50 km) with horizontal dimensions of approximately 1554 km and 2109 km ($\Delta \phi$ = 14, $\Delta \lambda$ = 19°) and a density contrast of −50 kg/m³. M2 is a synthetic model of a contoured geoid anomaly caused by one low-D rectangular prism on the right and one high density (high-D) rectangular prism on the left that ranges from 600 km to 150 km in depth ($\Delta r$ = 450 km) and has horizontal dimensions of 333 km ($\Delta \phi$ = $\Delta \lambda$ = 3°) and a density contrast of −50 kg/m³. M3 is similar to M2 but corresponds to two low-D rectangular prisms with different depths. We added 5% Gaussian random noise to the synthetic data, which are referred to as the “observed geoid data”. We implemented two tests based on two different types of models whose inversion parameters are $\beta$ = 0.9; ${\epsilon }_1$ = 10⁻⁴; and ${\epsilon }_2$ = 10⁻². The numbers of iterations, modeling errors and fitting errors are listed in

Table 1. Conditions and results for three different synthetic inversion models.

Observations ($\stackrel{\mathbf{\rightharpoonup }}{{\mathit{n}}_{\mathit{\rho }}}$)	Initial (${\widehat{\mathit{\rho }}}_{\mathbf{0}}$)	L-corner ($\mathit{\mu }$)	Iterations ($\mathit{k}$)	ME	FE
G1 + 5% noise	0	225	7	2.79	0.0134
G2 + 5% noise	0	225	5	2.63	0.0132
G3 + 5% noise	0	268	6	3.07	0.0821

G1 (or G2 or G3) is the geoid ($N_{\rho }$) data for M1 (or M2 or M3); ME is the RMS modeling error (${\widehat{\rho }}_k-{\rho }_{ture}$); FE is the RMS fitting error ($G_{\rho }{\widehat{\rho }}_k-N_{\rho }$); and ${\rho }_{ture}$ is the true model.

.

As suggested by Figure 3b (or Figure 3d,f), which shows the estimated counterpart of Figure 3a (or Figure 3c,e), we always obtain acceptable results, particularly for the first two models, M1 and M2, which indicates the validity and correctness of our inversion method.

5.2. The TP Density Inversion

The model of the crust and upper mantle that is used in the forward model of the geoid is discretized with tesseroid blocks [30], each of which has a uniform density. The model is composed of six layers that image the density structures of the TP at depths of up to 670 km over a horizontal range of approximately 222 km (2°) and a variable vertical thickness. The first five layers are 100 km thick, and the last layer is 170 km thick. The transition depths of the model are 50 km, 150 km, 250 km, 350 km, 450 km, and 600 km.

The initial density model used in the inversion is based on seismic velocities from the tomography using the empirical relationship between velocity and density derived by Woodhouse and Dziewonski [37],

$$\delta \rho =\alpha \delta V_S^2$$

(9)

where $\delta \rho$ is the density contrast g/cm³; $\delta V_S$ is the velocity anomaly in km/s based on the seismic Vs model derived from seismic tomography developed by Simmons et al. [38] and the starting average Vs model TNA-SNA [39]; and $\alpha$ is the transfer coefficient, which has an empirical value of 3.13 × 10⁵.

We invert the residual geoid anomaly based on Equation (7). The initial density contrast, ${\widehat{\rho }}_0$, used in the inversion is $\delta \rho$. We refine the density model iteratively to reduce the fitting error. We terminate the iterations as soon as the difference between the fitting errors of two consecutive iterations is below a pre-specified value (here 0.01 m). Figure 4a, which shows the L-curve, demonstrates that the regularization parameter (namely, $\mu$, which is given in Equation (5)) is equal to 300 in this case. As illustrated by Figure 4b, the fitting error decreases with the number of iterations, and after approximately 50 iterations, the inversion is almost stable. The final residual geoid anomaly is shown in Figure 4c. Most of the final residuals are high-frequency errors.

Our results reveal some large-scale, structurally controlled density variations at depths illustrated in Figure 5a–f.

Across the first layer, density anomalies at an average depth of 50 km from the crust vary intensely. As shown in Figure 5a, in the eastern TP, low-D anomalies in the crust terminate near the Longmen Shan fault (LMSF). In the southern TP, low-D anomalies in the crust extend over a long distance. In the northern TP, crustal low-D anomalies stop beneath the ATF. These results agree well with the LVZ found by previous studies [11,27]. Low-D anomalies are also found in the Qilian Fold Belt. These results are in line with the high-precision underground gravity anomalies [40]. Additionally, low-D anomalies with the amplitudes of up to −20 kg/m³ occur in blocks on the eastern TP from north to south. Moreover, high-D anomalies with amplitudes of up to 60 kg/m³ are clustered on the western TP. From west to east in the study area, crustal high-D anomalies appear in Tianshan, the Tarim Basin, the Qiadam Basin, and the Sichuan Basin; they could be due to the hard basements in those regions. High-D anomalies can also be found in the southwestern TP (the Himalaya, Lhasa, and Qiangtang Blocks) extending downward to 300 km. An analysis of recent earthquakes (see

Table 2. Key parameters of the earthquakes involved in this study.

Name	Origin Date	Epicenter	Depth of Hypocenter	Magnitude
Wenchuan	12 May 2008	30.99°N, 103.32°E	10~20 km	8.0 M
Lushan	20 April 2013	30.3°N, 103°E	13 km	7.0 M
Ludian	3 August 2014	27.1°N, 103.3°E	12 km	6.5 M
Yingjiang	10 March 2011	25°N, 97.8°E	12 km	6.1 M

for their parameters) in this region showed that all of them occurred in transition zones between high-D and low-D density anomalies in the crust.

In the second layer, the characteristics of the density anomalies at an average depth of 150 km, which are in the upper mantle (Figure 5b), agree well with the results for the first layer, except that the higher densities in the center of the Songpan-Ganzi Block decrease in amplitude. In southeastern Tibet, particularly below Yunan Province, blocked low-D anomalies form a closed loop. Furthermore, earthquakes occurred in the transition zones between high-D and low-D anomalies in the uppermost mantle.

The density anomalies in the third layer (Figure 5c), which have an average depth of 250 km, are analogous to the results shown in Figure 5b, except that Figure 5c shows the amplitudes of the anomalies. Prominent low-D anomalies occur near the higher-density periphery of the TP. High-D anomalies below southern Qiangtang and northern Lhasa are always present from 0 km to 300 km and stretch to levels deeper than 300 km.

In the fourth layer, the density anomalies, which have an average depth of 350 km and are in the 300–400 km mantle, vary significantly. They exhibit a pattern of apparently weak density zones across the TP, with an EW strike in central Tibet, where they rotate until they are aligned in an approximately NW-SE direction toward eastern Tibet, which is similar to the low-Pn velocity belt extending from west to east beneath Qiangtang and Songpan-Ganzi, then turn southward beneath the eastern TP and extend farther south [41,42]. Obvious high-D anomalies (approximately 40 kg/m³) are present in the southern TP and on the northwestern TP. In the southern TP, blocked high density anomalies have already reached over the BNS and almost reach the Jinsha River Suture (JRS). In the northwestern TP, the high density anomaly should represent the south-dipping subduction of the Asian lithospheric plate.

In the fifth layer, the density anomalies, which have an average depth of 450 km, which is part of the 400–500 km upper mantle, exhibit a striking change. Low density anomalies are present throughout the TP. In the southeastern TP, three high-density anomaly zones (HDZs) with amplitudes of up to 30 kg/m³ are obviously present in the mantle transition zone along the Bruma arc. On the eastern side of those three HDZs, a weak HDZ with a value of approximately −5 kg/m³ emerged.

In the sixth layer, the density anomalies, which have an average depth of 600 km, which is in the 500–670 km mantle transition zone, are different from the anomalies discussed above. The trend here is for low-D anomalies to occur in the western region while high-D anomalies occur in the eastern region. As shown in Figure 5f, in the inner TP, low-D anomalies are blocked, whereas high-D anomalies occur outside the TP. The low-D anomalies in TP illustrate the existence of hot weak material beneath the TP, which is beneficial to the inward-thrusting external material. It is interesting to note that on the eastern TP, beneath the epicenters of all the recent earthquakes (WCH—Wenchuan, LSH—Lushan, LD—Ludian, YJ—Yingjiang), the density anomalies at this depth are high.

We conclude this section with Figure 6, which consists of two subfigures that show the vertical density at latitude 25°N from longitude 90°E to 105°E and at latitude 31°N from longitude 100°E to 110°E. According to Figure 6a, there is an evident low-D anomaly in the upper mantle in the south-eastern TP. Moreover, in the mantle transition zone, it is found that the high-D anomaly experienced by the subducted Indian Plate is more prominent than that experienced by the Bruma Plate, which implies that the Indian Plate may have dived into the mantle transition zone and consequently, formed a mantle wedge. Note that the presence of the mantle wedge mainly accounts for the deep origin of the Tengchong volcano and for the crustal earthquake. As for Figure 6b, the resistance originating from the hard, high-D Sichuan basin is responsible for the strain accumulated in the locking region along the northeastern boundary of the TP. Therefore, this effect leads to the potential for earthquakes to occur along the edges of the plates.

6. Discussion

Studies of deep seismic soundings showed that the northern edge of the Indian crust is approximately 50–100 km south of the BNS (east of 90°E) [11]. Receiver function results signified that the Indian crust can be traced to 31°N [2]. However, P wave tomography derived by Zhao et al. showed that the Indian Plate is currently sub-horizontal and under-thrusting to the JRS at depths from approximately 100 to 250 km [43]. The high velocity zone (HVZ) estimated by P wave tomography inferred that Indian Plate underlies only the southwestern plateau [44]. Similar to the HVZ, the high-D anomalies also trace the subduction of the Indian Plate and can be used to determine how far the Indian lithospheric plate is subducted beneath the Tibetan Plateau. Our results show that high-D anomalies are predominant in the southwestern TP (the Himalaya, Lhasa, and Qiangtang Blocks) at depths from 0 to 300 km and terminate near the JRS in the south. Based on the range of high-D anomalies inside the TP, we can infer that the under-thrusting Indian Plate has reached over the BNS and almost reaches the JRS at a depth greater than 300 km.

The P wave tomography conducted by Lei and Zhao found that broad high-V anomalies are visible from the Burma arc northward in the mantle transition zone (MTZ) beneath eastern Tibet, which is also the case in several other studies [45,46,47], and it has been inferred that they may represent the eastward subduction of the Indian slab along the Burmese arc. That our high-D anomalies are at a depth of 400–500 km (the depth of the MTZ) provides additional evidence for eastward subduction.

Teleseismic P-wave tomography shows that the velocities are lower in both the lower crust and the uppermost mantle under southwestern Yunnan based on regional waveform inversion [48]. Beneath the Yunnan region, the Poisson’s ratio is high and the crust is thin [49]. In our results, low-D anomalies emerge in both the crust and the lithospheric mantle (at depths between 0 and 300 km). Analogous to the velocity anomalies, low-D anomalies in the crust are the result of the extrusion of the Indian Plate and the obstruction by the surrounding hard blocks [50]. This results in a higher temperature and lower Poisson’s ratio in the lower crust. The low-D anomalies in the lithospheric mantle may be due to deep slab dehydration. The low-D anomalies in the MTZ are due to the subducted heavy plate. This provides forces and environments that are conducive to frequent earthquakes [50,51,52] (such as the Wenchuan, Lushan, Ludian, and Yingjiang earthquakes) and the formation of the Tengchong volcano in this area [51].

In the south, the Indian Plate underthrusts Tibet. In central and northern Tibet, there is a separate, thin Tibetan Plate that is underthrust from the north by the Asian Plate [53,54]. Hot thermal anomalies are present [15]. At an average depth of 600 km, low-D anomalies occur inside the TP in combination with high-D anomalies outside the TP, which result in incidental inward-thrusting from the cold, hard external TP to the hot, weak internal TP.

7. Conclusions

In this paper, we inverted the geoid in both simulations and real cases. The geoid anomaly can be used to see long-wavelength, detailed density features underground. We obtained three-dimensional structures of the crust and upper mantle beneath the eastern margin of the Tibetan Plateau. We summarize our major findings as follows:

(1)

A certain degree of consistency in the density anomalies at depths from 0 km to 300 km, except in the middle Songpan Ganzi, indicates a coupling between the crust and the upper mantle.

(2)

At depths between 0 km and 300 km depth, high density (high-D) anomalies are apparent in the southwestern TP; these terminate near the Jinsha River Suture (JRS) in the south. This illustrates that the under-thrusting Indian Plate has reached over the Bangong Nujiang Suture (BNS) and almost reaches the JRS at a depth greater than 300 km.

(3)

On the eastern TP, high density anomalies occur at depths between 400 km and 500 km (the depth of mantle transition zone, MTZ) while lower density (low-D) anomalies occur at depths between 100 km and 300 km (the depth of upper mantle). Low velocities in the crust and uppermost mantle and high velocities in the MTZ beneath this region further verify the current eastward subduction of the Indian Plate along the Burma arc. The ongoing subduction provides forces and environments that are conducive to frequent earthquakes and the appearance of the Tengchong volcano in this area.

(4)

At a depth of 600 km, low-D anomalies appear inside the TP and high-D anomalies outside the TP; these account for the incidental inward-thrusting from the cold, hard external TP to the hot, weak internal TP.