Next Article in Journal / Special Issue
Combined Forecasting Method of Landslide Deformation Based on MEEMD, Approximate Entropy, and WLS-SVM
Previous Article in Journal
Evaluating the Impact the Weekday Has on Near-Repeat Victimization: A Spatio-Temporal Analysis of Street Robberies in the City of Vienna, Austria
Previous Article in Special Issue
Dam Deformation Monitoring Data Analysis Using Space-Time Kalman Filter
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
State Key Laboratory of Geodesy and Earth’s Dynamics, Institute of Geodesy and Geophysics, Chinese Academy of Sciences, Wuhan 430077, China
2
University of Chinese Academy of Sciences, No. 19A Yuquan Road, Beijing 100049, China
*
Author to whom correspondence should be addressed.
ISPRS Int. J. Geo-Inf. 2017, 6(1), 4; https://0-doi-org.brum.beds.ac.uk/10.3390/ijgi6010004
Submission received: 14 October 2016 / Revised: 11 December 2016 / Accepted: 19 December 2016 / Published: 30 December 2016
(This article belongs to the Special Issue Recent Advances in Geodesy & Its Applications)

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 km2 [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 ( P ) = T ( P ) γ ( P ) , where P ( r , φ , λ ) is the computation point; r , φ and λ are the radius, latitude, and longitude, respectively; and γ ( P ) is the normal gravity on the geodetic reference ellipsoid. The perturbing potential, T ( P ) , 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 ( P ) = G ρ Δ r Δ φ Δ λ γ ( P ) [ K 000 + 1 24 ( K 200 Δ r 2 + K 020 Δ φ 2 + K 002 Δ λ 2 ) + O ( Δ 4 ) ]
where G is Newton’s gravitational constant; ρ is the density contrast, which is assumed to be constant within each tesseroid; Δ r = r 1 r 2 , Δ φ = φ 1 φ 2 , and Δ λ = λ 1 λ 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 ( Δ 4 ) is the Landau symbol, which indicates omitted fourth-order terms in Δ r , Δ φ , Δ λ .
Because the derivative of Equation (1) is linear with respect to the density contrast, ρ , a set of n calculated geoid anomalies caused by a model composed of m tesseroids can be written as
N ρ ( n × 1 ) = G ρ ( n × m ) ρ ( m × 1 )
where ρ is the vector that represents the estimated density contrast; N ρ is the vector that represents the geoid anomaly; and G ρ 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 ( ρ , α ) = ( ρ ρ a p r ) T W ρ T W ρ ( ρ ρ a p r ) + 2 α T W d ( G ρ ρ N ρ )
where the first term is the 2-norm of the fitting error, the condition term is the constraint, ρ a p r is an initial model, and W ρ 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 ρ , we use the matrix W ρ = d i a g ( 1 ( z + ε ) β ) to calculate the weight in the matrix G ρ , where z is the depth of each tesseroid, ε is a small constant to avoid singularities, β is a number that changes the weight of W ρ ; W d is the data weight matrix; and α is the vector of Lagrange multipliers. To minimize the objective function, we differentiate Equation (3) with respect to ρ and α ; by setting them equal to zero, we obtain the following solution:
{ ρ ^ = ρ a p r + ( W ρ T W ρ ) 1 G ρ T W d T ( W d G ρ ( W ρ T W ρ ) 1 G ρ T W d T ) 1 W d ( N ρ G ρ ρ a p r ) α = ( W d G ρ ( W ρ T W ρ ) 1 G ρ T W d T ) 1 W d ( N ρ G ρ ρ a p r )
where ρ ^ 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 ρ , we add the constant μ to Equation (4),
ρ ^ = ρ a p r + ( W ρ T W ρ ) 1 G T W d T ( W d G ( W ρ T W ρ ) 1 G T W d T + μ I ) 1 W d ( N ρ G ρ a p r )
where μ 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:
ρ ^ k + 1 = ρ ^ k + Δ ρ ^ k
where ρ ^ k is the density contrast after the kth iteration; for k = 0 , 1 , 2 , , n , where ρ ^ 0 is ρ a p r , which is an induced initial density model given by a 3-Dimension Vs model; and Δ ρ ^ k is as follows:
Δ ρ ^ k = ( W ρ T W ρ ) 1 G T W d T ( W d G ( W ρ T W ρ ) 1 G T W d T + μ I ) 1 W d ( N ρ G ρ ^ k )
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 ρ G ρ ρ ^ k must be less than or equal to a constant value, ε 1 . Another criterion for convergence is that E 2 = Δ ρ ^ k / ρ ^ k is smaller than a certain previously determined empirical value, ε 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 ρ = N ρ N ρ 1 N ρ 2 N ρ 3 N ρ 4 N ρ 5
where N ρ is the vector of the observed geoid anomaly; N ρ 1 is the vector of density changes below the upper mantle’s geoid effect; N ρ 2 is the vector of the topographic geoid effect; N ρ 3 is the vector of the geoid effect due to variations in the sedimentary interface; N ρ 4 is the vector of the Moho interface variation effect; N ρ 5 is the vector of the effect of lithosphere interface variations on the geoid; and n ρ 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/m3 and 1640 kg/m3 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/cm3 [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/cm3, the average depth of lithosphere is 100 km and the density contrast between the crust and the uppermost mantle is 10 kg/cm3 [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 ( Δ r = 50 km) with horizontal dimensions of approximately 1554 km and 2109 km ( Δ φ = 14, Δ λ = 19°) and a density contrast of −50 kg/m3. 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 ( Δ r = 450 km) and has horizontal dimensions of 333 km ( Δ φ = Δ λ = 3°) and a density contrast of −50 kg/m3. 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 β = 0.9; ε 1 = 10−4; and ε 2 = 10−2. The numbers of iterations, modeling errors and fitting errors are listed in Table 1.
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],
δ ρ = α δ V S 2
where δ ρ is the density contrast g/cm3; δ 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 α is the transfer coefficient, which has an empirical value of 3.13 × 105.
We invert the residual geoid anomaly based on Equation (7). The initial density contrast, ρ ^ 0 , used in the inversion is δ ρ . 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, μ , 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/m3 occur in blocks on the eastern TP from north to south. Moreover, high-D anomalies with amplitudes of up to 60 kg/m3 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 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/m3) 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/m3 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/m3 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.

Acknowledgments

We gratefully acknowledge the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB18010304), and the Major State Research Development Program of China (Grant No. 2016YFC0601101) for supporting the work. We appreciate the sharing of the S-wave velocity. The U.S. National Geospatial-Intelligence Agency (NGA) is thanked for the use of the EGM2008 data. We also acknowledge the use of the Generic Mapping Tools (GMT).

Author Contributions

Jian Fang provided the initial idea for this study; Honglei Li conceived and designed all the synthetic inversion experiments and applied them to the eastern Qinghai-Tibet Plateau; Honglei Li and Jian Fang analyzed the experimental results; and Honglei Li wrote the paper.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Kang, S.; Xu, Y.; You, Q.; Flügel, W.-A.; Pepin, N.; Yao, T. Review of climate and cryospheric change in the Tibetan Plateau. Environ. Res. Lett. 2010, 5, 015101. [Google Scholar] [CrossRef]
  2. Nábělek, J.; Hetényi, G.; Vergne, J.; Sapkota, S.; Kafle, B.; Jiang, M.; Su, H.; Chen, J.; Huang, B.-S. Underplating in the Himalaya-Tibet collision zone revealed by the Hi-CLIMB experiment. Science 2009, 325, 1371–1374. [Google Scholar] [CrossRef] [PubMed]
  3. Molnar, P. A review of geophysical constraints on the deep structure of the Tibetan Plateau, the Himalaya and the Karakoram, and their tectonic implications. Philos. Trans. R. Soc. Lond. A Math. Phys. Eng. Sci. 1988, 326, 33–88. [Google Scholar] [CrossRef]
  4. England, P.; Houseman, G. Finite strain calculations of continental deformation: 2. Comparison with the India-Asia collision zone. J. Geophys. Res. Solid Earth 1986, 91, 3664–3676. [Google Scholar] [CrossRef]
  5. Tapponnier, P.; Xu, Z.; Roger, F.; Meyer, B.; Arnaud, N.; Wittlinger, G.; Yang, J. Oblique stepwise rise and growth of the Tibet Plateau. Science 2001, 294, 1671–1677. [Google Scholar] [CrossRef] [PubMed]
  6. Tapponnier, P.; Peltzer, G.; Le Dain, A.; Armijo, R.; Cobbold, P. Propagating extrusion tectonics in Asia: New insights from simple experiments with plasticine. Geology 1982, 10, 611–616. [Google Scholar] [CrossRef]
  7. Clark, M.K.; Royden, L.H. Topographic ooze: Building the eastern margin of Tibet by lower crustal flow. Geology 2000, 28, 703–706. [Google Scholar] [CrossRef]
  8. Bourjot, L.; Romanowicz, B. Crust and upper mantle tomography in Tibet using surface waves. Geophys. Res. Lett. 1992, 19, 881–884. [Google Scholar] [CrossRef]
  9. Kind, R.; Yuan, X.; Saul, J.; Nelson, D.; Sobolev, S.; Mechie, J.; Zhao, W.; Kosarev, G.; Ni, J.; Achauer, U. Seismic images of crust and upper mantle beneath Tibet: Evidence for Eurasian plate subduction. Science 2002, 298, 1219–1221. [Google Scholar] [CrossRef] [PubMed]
  10. Li, S.; Mooney, W.D. Crustal structure of China from deep seismic sounding profiles. Tectonophysics 1998, 288, 105–113. [Google Scholar] [CrossRef]
  11. Zhang, Z.; Deng, Y.; Teng, J.; Wang, C.; Gao, R.; Chen, Y.; Fan, W. An overview of the crustal structure of the Tibetan plateau after 35 years of deep seismic soundings. J. Asian Earth Sci. 2011, 40, 977–989. [Google Scholar] [CrossRef]
  12. Bai, D.; Unsworth, M.J.; Meju, M.A.; Ma, X.; Teng, J.; Kong, X.; Sun, Y.; Sun, J.; Wang, L.; Jiang, C. Crustal deformation of the eastern Tibetan plateau revealed by magnetotelluric imaging. Nat. Geosci. 2010, 3, 358–362. [Google Scholar] [CrossRef]
  13. Sun, J.; Jin, G.; Bai, D.; Wang, L. Sounding of electrical structure of the crust and upper mantle along the eastern border of Qinghai-Tibet Plateau and its tectonic significance. Sci. China Ser. D Earth Sci. 2003, 46, 243–253. [Google Scholar]
  14. Mahéo, G.; Guillot, S.; Blichert-Toft, J.; Rolland, Y.; Pêcher, A. A slab breakoff model for the Neogene thermal evolution of South Karakorum and South Tibet. Earth Planet. Sci. Lett. 2002, 195, 45–58. [Google Scholar] [CrossRef]
  15. McKenzie, D.; Jackson, J.; Priestley, K. Thermal structure of oceanic and continental lithosphere. Earth Planet. Sci. Lett. 2005, 233, 337–349. [Google Scholar] [CrossRef]
  16. Chaves, C.A.M.; Ussami, N. Modeling 3-D density distribution in the mantle from inversion of geoid anomalies: Application to the Yellowstone Province. J. Geophys. Res. Solid Earth 2013, 118, 6328–6351. [Google Scholar] [CrossRef]
  17. LaCoste, L.J. Measurement of gravity at sea and in the air. Rev. Geophys. 1967, 5, 477–526. [Google Scholar] [CrossRef]
  18. Zhao, L.; Wu, M.; Forsberg, R.; Olesen, A.V.; Zhang, K.; Cao, J. Airborne gravity data denoising based on empirical mode decomposition: A Case study for SGA-WZ GREENLAND test data. ISPRS Int. J. Geo-Inf. 2015, 4, 2205–2218. [Google Scholar] [CrossRef] [Green Version]
  19. Pavlis, N.K.; Holmes, S.A.; Kenyon, S.C.; Factor, J.K. An earth gravitational model to degree 2160: EGM2008. In Proceedings of the EGU General Assembly, Vienna, Austria, 13–18 April 2008.
  20. Zhang, C.Y.; Guo, C.X.; Chen, J.Y.; Zhang, L.M.; Wang, B. EGM 2008 and its application analysis in Chinese Ma inland. Acta Geod. Cartogr. Sin. 2009, 38, 283–289. [Google Scholar]
  21. Koch, K.-R. Parameter Estimation and Hypothesis Testing in Linear Models; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2013. [Google Scholar]
  22. Unsworth, M.; Wenbo, W.; Jones, A.G.; Li, S.; Bedrosian, P.; Booker, J.; Sheng, J.; Ming, D.; Handong, T. Crustal and upper mantle structure of northern Tibet imaged with magnetotelluric data. J. Geophys. Res. Solid Earth 2004, 109. [Google Scholar] [CrossRef]
  23. Shin, Y.H.; Xu, H.; Braitenberg, C.; Fang, J.; Wang, Y. Moho undulations beneath Tibet from GRACE-integrated gravity data. Geophys. J. Int. 2007, 170, 971–985. [Google Scholar] [CrossRef]
  24. Teng, J.-W. My understand and ponder for few important kernel problems in and great achievements for first time and developmental direction of Geophysical research in Tibetan Plateau of China. Prog. Geophys. 2008, 23, 301–318. [Google Scholar]
  25. Klemperer, S.L. Crustal flow in Tibet: Geophysical evidence for the physical state of Tibetan lithosphere, and inferred patterns of active flow. Geol. Soc. Lond. Spec. Publ. 2006, 268, 39–70. [Google Scholar] [CrossRef]
  26. Yue, H.; Chen, Y.J.; Sandvol, E.; Ni, J.; Hearn, T.; Zhou, S.; Feng, Y.; Ge, Z.; Trujillo, A.; Wang, Y. Lithospheric and upper mantle structure of the northeastern Tibetan Plateau. J. Geophys. Res. Solid Earth 2012, 117. [Google Scholar] [CrossRef]
  27. Yang, Y.; Ritzwoller, M.H.; Zheng, Y.; Shen, W.; Levshin, A.L.; Xie, Z. A synoptic view of the distribution and connectivity of the mid-crustal low velocity zone beneath Tibet. J. Geophys. Res. Solid Earth 2012, 117. [Google Scholar] [CrossRef]
  28. Zhang, P.; Deng, Q.; Zhang, G.; Ma, J.; Gan, W.; Min, W.; Mao, F.; Wang, Q. Active tectonic blocks and strong earthquakes in the continent of China. Sci. China Ser. D Earth Sci. 2003, 46, 13–24. [Google Scholar]
  29. Heiskanen, W.A.; Moritz, H. Physical geodesy. Bull. Géod. (1946–1975) 1967, 86, 491–492. [Google Scholar] [CrossRef]
  30. Heck, B.; Seitz, K. A comparison of the tesseroid, prism and point-mass approaches for mass reductions in gravity field modelling. J. Geod. 2007, 81, 121–136. [Google Scholar] [CrossRef]
  31. Tikhonov, A.N.; Arsenin, V.Y. Solutions of Ill-Posed Problems; Winston & Sons: Washington, DC, USA, 1977. [Google Scholar]
  32. Bowin, C. Depth of principal mass anomalies contributing to the earth’s geoidal undulations and gravity anomalies. Mar. Geod. 1983, 7, 61–100. [Google Scholar] [CrossRef]
  33. Amante, C.; Eakins, B. ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis; NOAA Technical Memorandum NESDIS NGDC-24; Biblioteca Digital ILCE; National Geophysical Data Center: Boulder, CO, USA, 2009.
  34. Fang, J.; Xu, H. Three dimensional distribution of liothsperic density beneath the china and its anjacent region. Chin. J. Geophys. 2002, 45, 42–48. [Google Scholar]
  35. Laske, G.; Masters, G.; Ma, Z.; Pasyanos, M. CRUST1.0: An updated global model of earth’s crust. In Proceedings of the EGU General Assembly, Vienna, Austria, 22–27 April 2012.
  36. Pasyanos, M.E.; Masters, T.G.; Laske, G.; Ma, Z. LITHO1.0: An updated crust and lithospheric model of the Earth. J. Geophys. Res. Solid Earth 2014, 119, 2153–2173. [Google Scholar] [CrossRef]
  37. Woodhouse, J.H.; Dziewonski, A.M. Mapping the upper mantle: Three-dimensional modeling of Earth structure by inversion of seismic waveforms. J. Geophys. Res. Solid Earth 1984, 89, 5953–5986. [Google Scholar] [CrossRef]
  38. Simmons, N.A.; Forte, A.M.; Boschi, L.; Grand, S.P. GyPSuM: A joint tomographic model of mantle density and seismic wave speeds. J. Geophys. Res. Solid Earth 2010, 115. [Google Scholar] [CrossRef]
  39. Grand, S.P.; Helmberger, D.V. Upper mantle shear structure of North America. Geophys. J. Int. 1984, 76, 399–438. [Google Scholar] [CrossRef]
  40. Wang, X.; Fang, J.; Hsu, H. Three-dimensional density structure of the lithosphere beneath the North China Craton and the mechanisms of its destruction. Tectonophysics 2014, 610, 150–158. [Google Scholar] [CrossRef]
  41. Liang, C.; Song, X. A low velocity belt beneath northern and eastern Tibetan Plateau from Pn tomography. Geophys. Res. Lett. 2006, 33. [Google Scholar] [CrossRef]
  42. Lei, F.; Qu, Y.; Song, G. Species diversification and phylogeographical patterns of birds in response to the uplift of the Qinghai-Tibet Plateau and Quaternary glaciations. Curr. Zool. 2014, 60, 149–161. [Google Scholar] [CrossRef]
  43. Zhao, J.; Zhao, D.; Zhang, H.; Liu, H.; Huang, Y.; Cheng, H.; Wang, W. P-wave tomography and dynamics of the crust and upper mantle beneath western Tibet. Gondwana Res. 2014, 25, 1690–1699. [Google Scholar] [CrossRef]
  44. Zhao, J.; Yuan, X.; Liu, H.; Kumar, P.; Pei, S.; Kind, R.; Zhang, Z.; Teng, J.; Ding, L.; Gao, X. The boundary between the Indian and Asian tectonic plates below Tibet. Proc. Natl. Acad. Sci. USA 2010, 107, 11229–11233. [Google Scholar] [CrossRef] [PubMed]
  45. Lei, J.; Zhao, D. Teleseismic P-wave tomography and mantle dynamics beneath Eastern Tibet. Geochem. Geophys. Geosyst. 2016, 17, 1861–1884. [Google Scholar] [CrossRef]
  46. Sun, X.; Bao, X.; Xu, M.; Eaton, D.W.; Song, X.; Wang, L.; Ding, Z.; Mi, N.; Yu, D.; Li, H. Crustal structure beneath SE Tibet from joint analysis of receiver functions and Rayleigh wave dispersion. Geophys. Res. Lett. 2014, 41, 1479–1484. [Google Scholar] [CrossRef]
  47. Li, C.; van der Hilst, R.D.; Meltzer, A.S.; Engdahl, E.R. Subduction of the Indian lithosphere beneath the Tibetan Plateau and Burma. Earth Planet. Sci. Lett. 2008, 274, 157–168. [Google Scholar] [CrossRef]
  48. Li, Y.; Wu, Q.; Zhang, R.; Tian, X.; Zeng, R. The crust and upper mantle structure beneath Yunnan from joint inversion of receiver functions and Rayleigh wave dispersion data. Phys. Earth Planet. Inter. 2008, 170, 134–146. [Google Scholar] [CrossRef]
  49. Zha, X.H.; Lei, J.S. Crustal thickness and Poisson’s ratio beneath the Yunnan region. Sci. China Earth Sci. 2013, 43, 446–456. [Google Scholar] [CrossRef]
  50. Lei, J.; Li, Y.; Xie, F.; Teng, J.; Zhang, G.; Sun, C.; Zha, X. Pn anisotropic tomography and dynamics under eastern Tibetan plateau. J. Geophys. Res. Solid Earth 2014, 119, 2174–2198. [Google Scholar] [CrossRef]
  51. Lei, J.; Zhao, D.; Su, Y. Insight into the origin of the Tengchong intraplate volcano and seismotectonics in southwest China from local and teleseismic data. J. Geophys. Res. Solid Earth 2009, 114. [Google Scholar] [CrossRef]
  52. Lei, J.; Zhao, D. Structural heterogeneity of the Longmenshan fault zone and the mechanism of the 2008 Wenchuan earthquake (Ms 8.0). Geochem. Geophys. Geosyst. 2009, 10. [Google Scholar] [CrossRef]
  53. Zhao, D.; Pirajno, F.; Dobretsov, N.L.; Liu, L. Mantle structure and dynamics under East Russia and adjacent regions. Russ. Geol. Geophys. 2010, 51, 925–938. [Google Scholar] [CrossRef]
  54. Zhao, L.; Zheng, T.; Lu, G. Distinct upper mantle deformation of cratons in response to subduction: Constraints from SKS wave splitting measurements in eastern China. Gondwana Res. 2013, 23, 39–53. [Google Scholar] [CrossRef]
Figure 1. Topography and major tectonics of the density study area; in the figure: MBT—Main Boundary Thrust, MCT—Main Central Thrust, BNF—Bangong Nujiang Fault, JRS—Jinsha River Suture, ATF—Altun Tagh Fault, LMSF—Longmen Shan Fault, HM—Himalaya Block, LSB—Lhasa Block, QTB—Qiangtang Block, QB—Qaidam Block, SG—Songpan Ganzi Block, Tarim Basin—Tarim Basin, Ordos—Ordos Block, SCB—Sichuan Basin; Red stars mark the locations of earthquakes that have occurred in recent years and the volcano in the study area: WCH—epicenter of the 2008 8.0 M Wenchuan earthquake, LSH—epicenter of the 2013 7.0 M Lushan earthquake, LD—epicenter of the 2014 6.5 M Ludian earthquake, YJ—epicenter of the 2011 5.8 M Yingjiang earthquake; TC—Tengchong Volcano. (Note that, this is a new figure reusing only the fault data from [28].)
Figure 1. Topography and major tectonics of the density study area; in the figure: MBT—Main Boundary Thrust, MCT—Main Central Thrust, BNF—Bangong Nujiang Fault, JRS—Jinsha River Suture, ATF—Altun Tagh Fault, LMSF—Longmen Shan Fault, HM—Himalaya Block, LSB—Lhasa Block, QTB—Qiangtang Block, QB—Qaidam Block, SG—Songpan Ganzi Block, Tarim Basin—Tarim Basin, Ordos—Ordos Block, SCB—Sichuan Basin; Red stars mark the locations of earthquakes that have occurred in recent years and the volcano in the study area: WCH—epicenter of the 2008 8.0 M Wenchuan earthquake, LSH—epicenter of the 2013 7.0 M Lushan earthquake, LD—epicenter of the 2014 6.5 M Ludian earthquake, YJ—epicenter of the 2011 5.8 M Yingjiang earthquake; TC—Tengchong Volcano. (Note that, this is a new figure reusing only the fault data from [28].)
Ijgi 06 00004 g001
Figure 2. (a) Topography of the TP; (b) Moho depth of the TP; (c) sedimentary depth of the TP; (d) lithospheric depth of the TP; (e) effect of topographic masses; (f) effect of variations in the Moho interface; (g) effect of variations in the sedimentary interface; (h) effect of variations in the lithospheric interface; (i) geoid anomaly from EGM2008 to degree 180; (j) geoid anomalies from EGM2008 of degree 2–6; and (k) residual geoid anomaly after the above-motioned contributions are reduced.
Figure 2. (a) Topography of the TP; (b) Moho depth of the TP; (c) sedimentary depth of the TP; (d) lithospheric depth of the TP; (e) effect of topographic masses; (f) effect of variations in the Moho interface; (g) effect of variations in the sedimentary interface; (h) effect of variations in the lithospheric interface; (i) geoid anomaly from EGM2008 to degree 180; (j) geoid anomalies from EGM2008 of degree 2–6; and (k) residual geoid anomaly after the above-motioned contributions are reduced.
Ijgi 06 00004 g002
Figure 3. (a) Synthetic M1; (b) Estimated M1; (c) Synthetic M2; (d) Estimated M2; (e) Synthetic M3; and (f) Estimated M3.
Figure 3. (a) Synthetic M1; (b) Estimated M1; (c) Synthetic M2; (d) Estimated M2; (e) Synthetic M3; and (f) Estimated M3.
Ijgi 06 00004 g003aIjgi 06 00004 g003b
Figure 4. (a) L-curve; (b) the fitting error decreasing with the number of iterations; and (c) the final residual geoid anomaly (m).
Figure 4. (a) L-curve; (b) the fitting error decreasing with the number of iterations; and (c) the final residual geoid anomaly (m).
Ijgi 06 00004 g004
Figure 5. (a) Map of the horizontal density distribution at a depth of 50 km; (b) map of the horizontal density distribution at a depth of 150 km; (c) map of the horizontal density distribution at a depth of 250 km; (d) map of the horizontal density distribution at a depth of 350 km; (e) map of the horizontal density distribution at a depth of 450 km; and (f) map of the horizontal density distribution at a depth of 600 km.
Figure 5. (a) Map of the horizontal density distribution at a depth of 50 km; (b) map of the horizontal density distribution at a depth of 150 km; (c) map of the horizontal density distribution at a depth of 250 km; (d) map of the horizontal density distribution at a depth of 350 km; (e) map of the horizontal density distribution at a depth of 450 km; and (f) map of the horizontal density distribution at a depth of 600 km.
Ijgi 06 00004 g005aIjgi 06 00004 g005bIjgi 06 00004 g005c
Figure 6. (a) Vertical density sections at latitude 25°N from longitude 90°E to 105°E; the black solid line is the elevation across the profile; the red star shows the location of the Yingjiang earthquake; and the red triangle shows the location of the Tengchong volcano; (b) Vertical density sections at the latitude 31°N from longitude 100°E to 110°E; the black solid line is the elevation across the profile and the red stars show the locations of the Lushan and Wenchuan earthquakes.
Figure 6. (a) Vertical density sections at latitude 25°N from longitude 90°E to 105°E; the black solid line is the elevation across the profile; the red star shows the location of the Yingjiang earthquake; and the red triangle shows the location of the Tengchong volcano; (b) Vertical density sections at the latitude 31°N from longitude 100°E to 110°E; the black solid line is the elevation across the profile and the red stars show the locations of the Lushan and Wenchuan earthquakes.
Ijgi 06 00004 g006aIjgi 06 00004 g006b
Table 1. Conditions and results for three different synthetic inversion models.
Table 1. Conditions and results for three different synthetic inversion models.
Observations ( n ρ )Initial ( ρ ^ 0 )L-corner ( μ )Iterations ( k )MEFE
G1 + 5% noise022572.790.0134
G2 + 5% noise022552.630.0132
G3 + 5% noise026863.070.0821
G1 (or G2 or G3) is the geoid ( N ρ ) data for M1 (or M2 or M3); ME is the RMS modeling error ( ρ ^ k - ρ t u r e ); FE is the RMS fitting error ( G ρ ρ ^ k N ρ ); and ρ t u r e is the true model.
Table 2. Key parameters of the earthquakes involved in this study.
Table 2. Key parameters of the earthquakes involved in this study.
NameOrigin DateEpicenterDepth of HypocenterMagnitude
Wenchuan12 May 200830.99°N, 103.32°E10~20 km8.0 M
Lushan20 April 201330.3°N, 103°E13 km7.0 M
Ludian3 August 201427.1°N, 103.3°E12 km6.5 M
Yingjiang10 March 201125°N, 97.8°E12 km6.1 M

Share and Cite

MDPI and ACS Style

Li, H.; Fang, J. Crustal and Upper Mantle Density Structure Beneath the Qinghai-Tibet Plateau and Surrounding Areas Derived from EGM2008 Geoid Anomalies. ISPRS Int. J. Geo-Inf. 2017, 6, 4. https://0-doi-org.brum.beds.ac.uk/10.3390/ijgi6010004

AMA Style

Li H, Fang J. Crustal and Upper Mantle Density Structure Beneath the Qinghai-Tibet Plateau and Surrounding Areas Derived from EGM2008 Geoid Anomalies. ISPRS International Journal of Geo-Information. 2017; 6(1):4. https://0-doi-org.brum.beds.ac.uk/10.3390/ijgi6010004

Chicago/Turabian Style

Li, Honglei, and Jian Fang. 2017. "Crustal and Upper Mantle Density Structure Beneath the Qinghai-Tibet Plateau and Surrounding Areas Derived from EGM2008 Geoid Anomalies" ISPRS International Journal of Geo-Information 6, no. 1: 4. https://0-doi-org.brum.beds.ac.uk/10.3390/ijgi6010004

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop