Lithospheric Structure of Eastern Tibetan Plateau from Terrestrial and Satellite Gravity Data Modeling: Implication for Asthenospheric Underplating


 The lithosphere of the eastern Tibetan plateau is underlain by a low-velocity zone at shallow depths which is interpreted as asthenospheric material in the upper-most mantle in various seismic tomography studies. The driving mechanism for the presence of asthenospheric material in the upper-most mantle is not well understood, and the spatial extent of the asthenospheric material is not well delineated. We use 2.5D gravity models to assess what drove the asthenospheric flow upwards in the past and determine the lateral extent of the asthenospheric material in the upper-most mantle. The models also allow us to determine the Indian slab configuration below the Tibetan plateau. The gravity models show that lithospheric thickness increases from ~120 km in the central and eastern parts of the plateau to ~150 km in the west, indicating that the lithosphere in the central and eastern parts of the plateau may have been delaminated. The ~30 km shallower Lithosphere-Asthenosphere Boundary in the central and eastern Tibetan plateau may indicate that asthenospheric flow could have been induced in the past by a combination of lithospheric delamination and a slab break-off event of the Greater Indian slab. The spatial extent of the asthenospheric material in the upper-most mantle beneath the Tibetan plateau is ~15,000 km2 (N−S length=500 km and thickness=30 km) between 85°E and 88°E, which could even extend east of 92°E. The Indian slab is dipping more steeply in the east. The slab dip along the Indian plate increases from ~10° in the west to ~18° in the central (~87°E) and ~25° in the eastern part (~91°E) of the plateau, indicating that the style of lithospheric deformation changes from underthrusting to slab roll-back from west to east.


Introduction
The present day Tibetan plateau is the result of an~270 Myr long convergence [1][2][3] between the Indian and Asian plates (Figure 1(a); after [4]). The dynamic process included different subduction and suturing episodes during the closure of the Tethys Ocean, leading to multiple accretion events on the southern border of Asia [5][6][7][8][9]. The final collision between India and Asia occurred at~50-55 Ma [10][11][12], and they coupled with the Himalayan-Tibetan orogenic system.
Throughout the Tibetan plateau, the primary contributing factor for plateau growth is the subduction or underthrusting of the continental lithosphere. Geophysical studies have found evidence of a northward subduction of the Indian plate beneath the southern Tibetan plateau [13][14][15] and a southward subduction of the Eurasian plate beneath the northern Tibetan plateau ( [16]; W. [17]). The continental subduction of the Indian and Asian plates together contributes greatly to the plateau growth in southern and northern Tibet [18]. The lithospheric scale models of the central and eastern parts of the plateau are more complex and provide us with several plausible coexisting geodynamic scenarios, such as the lateral block extrusion model which suggests underthrusting of the Indian plate below the Tibet plateau [19][20][21][22], underthrusting accompanied by lower crustal flow beneath the central Tibetan plateau with weak viscous layers (e.g., [23][24][25][26]), and whole lithospheric deformation coupled with crustal flow towards the east ( [27,28]; Chen et al. [29]. These models Western 0−500 500−1,000 1,000−1,500 1,500−2,000 2,000−2,500 2,500−3,000 3,000−3,500 3,500−4,000 4,000−4,500 4,500−5,000 5,000−5,500 5,500−6,000 6,000−6,500 6,500−7,000 7,000−7,500 7,500−8,000 8,000−8,500 8,500−9,000 2 Lithosphere explain the accommodation of convergence between the Indian and Asian plates. But in contrast to the wellstudied crustal shortening and thickening in the region, there are several questions related to lithospheric structure and asthenospheric underplating which have been debated for a long time and remain unanswered. Some of these questions are related to the lithospheric structure and asthenospheric flow below the central Tibetan plateau. The central part of the Tibetan plateau is associated with weak crustal rheology marked by high temperature, poor S n propagation, low P n velocities, and large seismic anisotropy in the Qiangtang and Songpan-Ganze terranes (Figure 1), indicating that the crust of the eastern Tibetan plateau was hot during the last two subduction events from Early (~55 Ma) to Mid-Tertiary time (~25 Ma) ( [30][31][32][33][34][35][36][37][38][39][40]). Neogene magmatism in the southern Tibetan plateau and Karakorum is evidence for the most recent slab break-off of the Greater Indian continental slab [18,41,42]. Intrusion of ultrapotassic mafic magmatic material between 25 and 8 Ma in the southern Tibetan plateau (Lhasa terrane) may have been caused by rapid removal of the mantle lithosphere beneath the central Tibetan plateau after the last slab break off event [43][44][45][46]. This is also consistent with shear wave anomalies seen in recent work from Chen et al. [29] which suggests recent asthenospheric upwelling and thermal erosion of thickened lithosphere responsible for the observed potassic volcanism. Seismic tomography images a low velocity anomaly in the uppermost mantle which is interpreted as asthenospheric material [47][48][49]. Hence, a thinner lithospheric mantle underlain by a low-velocity asthenospheric material is suggested below the eastern Tibetan plateau. The driving mechanism for this asthenospheric flow in the past has not been well understood. There are two mechanisms that can explain what may have driven the asthenospheric flow below Tibetan plateau in the past. The last slab break-off event may have been responsible for inducing asthenospheric flow towards the uppermost mantle below the Tibetan plateau. An equally convincing hypothesis is that lithospheric delamination along with the slab break-off event may have induced asthenospheric flow in the uppermost mantle (e.g., [8]). The lateral extent of the asthenospheric material (low-velocity zone) in the uppermost mantle is not well delineated either. It is different in various tomography studies depending on the resolution and method used (cf., e.g., [15,22,50,51]). The motivation of this work is to test which amongst these hypothesis applies to this region and can explain what drove asthenospheric flow upwards in the past.
The western part of the Indian-Asian collision zone is at a more advanced stage of collision than the eastern part [52]. Recent tomographic studies indicate that slab dip becomes steeper from west to east (~8°to~20°, [49]). The mechanism of deformation changes from underthrusting/flat slab in the west to slab roll-back/steeper slab in the east, which may have been caused by eastward subduction and westward retreat of the Burmese arc during clockwise rotation of the Sunda block ( Figure 1; [13,53]). A numerical model from Li et al. [54] also suggests a causal relationship between a steeper slab angle and the removal of lithospheric mantle in the eastern Tibetan plateau. Petrologic data from Guo et al. [43] indicate postcollisional ultrapotassic mafic magmatism before 8 Ma in the west, but not east of 87.5°E, suggesting a steeper slab angle in the east. From west to east, the ultrapotassic mafic material is younger in age indicating eastward migration of the most recent slab detachment [43]. Shi et al. [49] suggest that the angle change occurs around 91.5°E. More recent work from Chen et al. [29] also suggests a steeper slab angle of the Indian lithosphere below the Tibetan plateau (>10°) as compared to the west. The exact change in slab angle and location is still much debated.
The south-north trending rift system on the plateau has been explained as an extension in the Tibetan plateau associated with eastward extrusion [55]. The central Tibetan crust has been spreading eastward via distributed eastward extrusion of small (<150 km wide) wedge-shaped blocks [56] due to two shear zones in ENE (left lateral) and WNW (right lateral) direction leading to V-shaped conjugate strike-slip faults [28]. Some of the rifts generated in the presence of these opposite sense shears, such as the Xainza-Dingjye Rift (XDR), the Yadon-Gulu Rift (YGR), and the Nyima-Tingri Rift (NTR), are perpendicular to the collisional boundary, which rules out the possibility that the slab roll-back in the east is a possible reason for their extension, as slab roll-back cannot provide extensional forces perpendicular to a collisional boundary. Several models have been proposed to explain the extensional forces at play for the widening of these rifts, such as gravitational collapse [57], lithospheric mantle delamination below the rifts [48], and slab tearing within the Indian slab because of a change in slab angle [58]. The relation between slab-tear, asthenospheric flow, slab angle change, and formation of rift zones remains not well understood.
In this paper, we analyze what drove asthenospheric flow upwards in the past and determine the spatial extent of asthenospheric material beneath the Tibetan plateau using terrestrial and satellite gravity data modeling. There are two main goals: (1) to construct a gravity model of the Tibetan plateau, which includes the density structure of the continental crusts, the lithosphere, the asthenosphere, and the dipping slab, and (2) to use the model to assess the configuration of the Indian slab below the Tibetan plateau. We present five representative gravity models of the Tibetan plateau and discuss their primary features and geodynamic implications.

Geological Setting
The Tibetan plateau is composed of a collage of continental fragments (terranes) that were added successively to the Eurasian plate during the Paleozoic and Mesozoic Eras [1,3,[59][60][61]. In the Tibetan plateau, five terranes exist (Figure 1(b)): Lhasa, Qiangtang, Songpan-Ganze, Kunlun-Qaidam, and Qilian Shan [3]. The Indus-Tsangpo suture marks the contact between India and these terranes ( Figure 1). The Himalayan thrust belt is south of the suture. The five terranes form an essential part of the gravity models presented in this study. The origin and geological evolution of the terranes help us in constraining the orientation, depth, and density ranges of geologic units in the gravity models.

Lithosphere
On the north-eastern boundary of the Tibetan plateau lies Qilian Shan. Qilian Shan formed gradually after the collision of the Tarim-North China cratons with the Kunlun-Qaidam continent. The timing of collision is not well constrained, but~445-440 Ma plutonism was followed by relatively little magmatism [62]. Prolonged convergence led to a shortening and exhumation of the early Paleozoic granitoids, UHP rocks, and ophiolite suites in the Devonian [63][64][65][66]. Evidence from the past subduction event has been extensively studied along what now forms the Qilian suture. It marks the northern boundary with the Qaidam Basin. The Qaidam Basin is the largest depression (1:2 × 105 km 2 ) inside the Tibetan plateau with an average elevation of 2.8 km. The basin preserves a 10 to 15 km thick complete record of the Cenozoic detrital sedimentation [67]. The southern boundary of Qaidam Basin is the Permian-Triassic Kunlun suture, which corresponds to the closure of the Songpan-Ganze oceanic domain [7]. The Kunlun suture forms the northern boundary of the Songpan-Ganze terrane. Songpan-Ganze is bounded by the Longmen Shan [7] in the east, and by the left-lateral Altyn-Tagh fault in the west [68,69]. Evidence for a collision between India and Asia can be seen here through the presence of the Kunlun volcanic group in northern Songpan-Ganze ( Figure 1). It is the last remaining volcanic expression of the collision, which is mostly dormant after the 1951 Ashi eruption. Due to the 5-15 km thick turbidite sequence, an oceanic lower crust is inferred for the Songpan-Ganze terrane [7]. The Jinsha suture lies at the southern boundary of the Songpan-Ganze terrane, which marks the closure of Paleo-Tethys.
The Qiangtang terrane lies south of the Late Triassic-Early Jurassic Jinsha suture, with the Bangong-Nujiang suture (BNS) as the southern boundary [70]. This terrane is approximately 500-600 km wide and its central part is dominated by Early Mesozoic blueschist rocks and Carboniferous to Cretaceous marine to volcaniclastic strata [3,11,71]. The Jurassic suture zone rocks are unconformably overlain by Mid-Cretaceous nonmarine red beds and volcanic rocks [71]. The next fragment attached to the southern border of Qiangtang is Lhasa. The Lhasa terrane is 300 km wide in its central part and extends further west across the active Karakorum fault [72][73][74][75]. The basement of Lhasa has a Gondwana affinity and is Mesoproterozoic to Early Cambrian in age [3,76]. Different metamorphic events are related to heat advection caused by potassic and ultrapotassic magmatism, and the magmatism may be a result of a break-off of the Indian continental slab before 25 Ma (e.g., [41,77,78]). The southern boundary of the Lhasa terrane corresponds to the Indus-Tsangpo suture zone (ITSZ) which formed around Early Eocene (52 Ma), which can be traced for more than 3000 km and marks the collision boundary of the Indian continent and Greater India with the Lhasa terrane.
Within the Lhasa terrane, the formation of N-S trending rifts ( Figure 1) in the southern Tibetan plateau began in the Late Neogene as a consequence of extensional tectonics (e.g., [55]). These rifts spread throughout the Tibetan plateau at regular intervals of 200 km, with dimensions within a few hundred kilometers running in a north-south direction lying perpendicular to suture zones [58]. The most prominent rifts lie south of the BNS [55] and oriented at an angle of~10°-15°w ith respect to longitudes. Many of these rifts cross through ITSZ with the Yadon-Gulu Rift (YGR) at~88-89°E and the North Lunggar Rift (NLR) at ∼82-83°E, extending into the Himalaya. Gravitational collapse, asthenospheric upwelling and slab tearing, lateral extrusion, and the geometry of the collision zone are some of the reasons cited for the presence of these rifts. Within the central part of the plateau, crustal segments escaped away from the colliding Indian plate along prominent strike-slip faults that originate out of the Tibetan plateau. Large portions of what is now Southeast Asia are escaping eastwards along major strike-slip faults to make space for India as it continues its motion in the northern direction ( Figure 1; [1,79]). Geodetic data indicate that crustal materials currently experience a significant clockwise rotation around the Eastern Himalaya Syntaxes along several N-S trending fault systems [80][81][82][83][84].

Data and Method
3.1. Gravity Data. The gravity data are based on the European Improved Gravity Model of the Earth (EIGEN-6C4; [85]). The EIGEN-6C4 geopotential model includes terrestrial and satellite-derived gravity data from the GRACE (Gravity Recovery and Climate Experiment), LAGEOS (Laser Geodynamics Satellites), and GOCE (Gravity field and steady-state Ocean Circulation Explorer) missions. The gravity data were obtained from the data portal of the International Centre of Global Earth Models (ICGEM) [86,87]. A standard gravity correction was applied to the gravity data to determine the complete Bouguer gravity anomaly of the Tibetan plateau ( Figure 2). The gravity data correction was based on the global topographic gravity field model (dV_ELL_RET2014_plusGRS80, [88]) and includes curvature and terrain corrections. A standard density (2670 kg/m 3 for rock and 1645 kg/m 3 for water) was used for the calculation of the Bouguer anomaly from the Earth Gravitational Model. The degree and order of the EIGEN-6C4 geopotential model is 2190, which gives us a spatial resolution of~9 km. Gravity disturbance values from the topographic model were calculated using isGrafLab [89]. The EIGEN-6C4 model is suitable for modeling of lithospheric blocks and can be used to analyze geologic structures whose lateral dimensions lie above the resolution of the geopotential models [90,91].

Data Constraints.
To determine the crust and upper mantle structures of the Tibetan plateau, five 2.5D gravity models have been developed down to a depth of 400 km. 2.5D models include the effect of lateral structures on each side of the models up to 100 m perpendicular to strike. These structures help in reducing the error attached with lateral end effects. The same thing is done along the length; overall, this provides the 2.5D models with a rough 3D effect around the profile. To reduce the ambiguity inherent in potential field interpretation, the gravity model is constrained by all available geological and geophysical data ( Figure 3 after [52]) in the Tibetan region (e.g., [92][93][94][95][96][97][98][99]). Recent seismic tomography and receiver function data have been used to constrain the P-wave velocities of the crust and upper mantle structure 4 Lithosphere of the Tibetan plateau [14,52,100]. The P-wave velocities are converted to densities using the methodology of Sobolev and Babeyko [101], which accounts for temperature, pressure, and compositional variation of density in the Tibetan plateau. This method works well with velocities up to~7.8 km/sec. However, the velocities in the upper mantle in the Tibetan plateau reach up to 8.5 km/sec. For the upper mantle, we have used the Nafe and Drake [102] empirical relation to convert velocity to density. The temperatures are obtained from the thermal models of Brewer et al. [103], Hetenyi et al. [104], and Tunini et al. [52], and the pressures were determined assuming a pressure-gradient of 30 MPa/km. Table 1 shows the P-wave velocities of major tectonic elements in the eastern Tibetan plateau, along with the corresponding densities. The crustal thickness for the Indian plate was obtained from seismic reflection [94] and receiver function analysis [97], and the values range from 40 km to 80-88 km close to the collisional boundary [97,99]. Wide-angle seismic reflection studies from INDEPTH II projects [94] helped to determine the thickness of the crust and density variations in the Lhasa [98] and southern Qiangtang terranes [93,96,98,99,105]. The Tibetan crust is~74 km thick in the southern Lhasa terrane and decreases to~65 km towards the north [36,93,99,106]. The crustal structure between the Qiangtang and Songpan-Ganze terranes has been determined using wideangle seismic profiles [95] and magnetotelluric data [107]. The crust thins northward to 60-65 km below the Qiangtang terrane [36,99,106], and to 55-60 km in the Songpan-Ganze terrane and Qaidam Basin [16,36,99,106].
The Lithosphere-Asthenosphere Boundary (LAB) beneath the Tibetan plateau was determined using seismic tomography. The LAB depth for the Indian slab in the western part of the Tibetan plateau gradually deepens to 200 km below the southern Qiangtang terrane and suddenly thins to 150 km below the northern Qiangtang, Jinsha, and Songpan-Ganze terranes [17,108,109]. In the east, seismic studies suggest an Indian LAB at a depth of 200-240 km below the central to eastern Tibetan plateau with the northern end of the Indian plate lying below the Bangong-Nujiang suture [17,18,[109][110][111]. The LAB decreases to 120-140 km below the Lhasa, Qiangtang, and Songpan-Ganze terranes. On the northern end of the profiles, the LAB depth increases to 160-180 km below Qaidam Basin due to the presence of a downgoing Asian slab [47,52,109]. The densities of the asthenosphere were constrained using an eastern and western Tibetan plateau tomographic model [52]. It is important to note that the densities of the asthenosphere have less constraints than the densities of the blocks in the lithosphere because of the lower resolution of tomographic studies at greater depths. The uncertainties in the final model were calculated by varying the densities without introducing a significant change in the fitted models. The positions of the slabs were constrained based on earthquake hypocenters from the USGS database for the last 30 years along with seismic tomography studies in the region (Li et al., 2012; [49,52,54]; Chen et al. [29].  5 Lithosphere which uses an algorithm from Talwani et al. [112] to calculate the gravity anomaly of 2D and 3D geologic structures. The 2.5D modeling profiles are small in comparison with the modeling space and mass of the entire Earth. This can create an edge effect on the forward models and alter anomaly values. To take into consideration the effect of the surround-ing mass on the forward gravity models, the dimension of the modeling space was extended to 1000 km. A finite strike length of 100 m has been used on each side (along width) of the gravity profile. To improve the fit of the observed gravity to predicted gravity, density values were inverted using a ridge-regression algorithm. The errors in the misfit were   7 Lithosphere minimized in a least-square sense. It is important to note here that 2.5D modeling is a first order approximation to 3D modeling. A 3D model can be generated with enough 2.5D models (with extended lateral space) within a region. The 2.5D models presented in this study have sufficient spatial resolution and constraints to answer the research problems. To ensure that the gravity model is isostatically balanced, we determined the isostatic Moho depth of the eastern Tibetan plateau based on the Airy isostatic model and adjusted the gravity model (Supplementary Material, Table S1). The Airy isostatic model is justified, because effective elastic thickness is small, which is equivalent to the Airy type isostatic compensation [113][114][115][116]. It is important to note that the Airy isostatic Moho model considers low flexural rigidity for a region, while the Vening Meinesz-Mortiz (VMM) Moho includes flexural rigidity for calculation. Hence, our Moho calculation is a first order approximation compared to other Moho studies based on the VMM model (cf., e.g., [117,118]).

Results and Discussion
4.1. Gravity Anomalies. The complete Bouguer anomaly map of the Tibetan plateau shows gravity variations along different suture zones and terranes (Figure 2). The values for the gravity anomalies are the lowest, -518 to -450 mGal, in the Qiangtang terrane where there is a low-density zone beneath the highest elevation of the plateau. The highest values, -150 to 60 mGal, are observed on the Indian plate, Tarim Basin, and along the Kunlun fault ( Figure 2). This can be indirectly related to a decrease in elevation from an average of 4800 m on the plateau to 1000-1500 m for the Indian plate. As we move from the central part of the plateau towards the south, the gravity anomalies change from low (~-450 mGal) to high (~-250 mGal) values along the Indus-Tsangpo suture zone, indicating a transition from low-density crustal roots under the Tibetan plateau to a higher density under the Indian plate. Gravity data from the complete Bouguer anomaly map was used to model structure along the profiles shown in Figure 2.
To determine the regional gravity anomalies in the Tibetan plateau, we use upward continuation and wavelength filtering. At~25 km above the surface of the Earth, most of the short wavelengths of the Bouguer gravity anomalies attenuate, and the anomalies show a regional component. Furthermore, we filter the Bouguer gravity anomalies using a low-pass filter for various cut-off values. Based on the similarity of amplitudes of regional gravity anomalies obtained from upward continuation and wavelength filtering, a cutoff wavelength of 300 km is selected. The resulting regional gravity anomaly map ( Figure 4) shows a strong negative anomaly (~-500 mGal) in the middle of the plateau correlating with the Qiantang terrane. This agrees with previous studies suggesting decreased crustal velocities (densities) below the Qiantang terrane or equivalently interpreted as increased crustal thickness [113][114][115][116]119].

Lithospheric Structure and Spatial Extent of
Asthenospheric Material. The gravity models (Figures 5-9) of the continental collision zone between the Indian and Eurasian plates show the crust and upper mantle structure of the Tibetan plateau along longitudes 82°E, 85°E, 88°E, 90°E, and 92°E, with latitudes spanning from 26°N to 38°N. The depth to the LAB beneath the Indian plate ranges from~160 to 200 km in the western part of the Tibetan plateau ( Figures 5  and 6), which agrees with the receiver function study of Zhao  8 Lithosphere et al. [109]. The LAB depth decreases across the Bangong-Nujiang suture zone (BNS) to~150 km, which is in agreement with previous studies of LAB for the BNS (e.g., [120]). The deep LAB in the south is attributed to the underthrusting of the Indian plate below the Tibetan plateau [17,110]. In comparison to the western part of the Tibetan plateau, the central and eastern parts of the plateau exhibit a drastic change in LAB depth. The LAB below the Indian plate lies at a depth of~180-240 km (Figures 7-9), and it decreases to 110-120 km below the Lhasa, Qiantang, and Songapan-Ganze terranes. Comparing laterally, the lithospheric thick-ness below the central part of the plateau decreases from 150 km to 110-120 km (Figures 6-9), indicating~30 km lithospheric thinning below the eastern part of Tibetan plateau. The thin lithosphere is underlain by a low-density asthenospheric material in the upper-most mantle (Figures 7-9), which is in agreement with receiver function analysis [48]. Thus, part of the central and eastern Tibetan lithosphere may have been delaminated, and the delaminated lithosphere may have induced asthenospheric flow in the upper most mantle in the past, leading to the asthenospheric underplating that we see today.   The law of fluid mechanics states that a large body sinking in a fluid can direct the fluid to flow upwards to fill the position that the body is leaving behind. This analogy is valid for the sinking delaminated lithospheric mantle, which could be inducing an asthenospheric flow upwards to fill the space previously filled by the lithospheric mantle. In these gravity models, we see asthenospheric material occupying the place previously occupied by the delaminated Tibetan lithospheric mantle. The spatial extent of the asthenospheric material in the upper-most mantle beneath the central Tibetan plateau is~15,000 km 2 (N − S length = 500 km and thickness = 30 km) between 85°E and 88°E, which could even extend east of 92°E. Various numerical models place the timing of this delamination at~8 Ma ( [121]; [122]). This recent delamination may have been one of the driving forces for the upward asthenospheric flow below the Tibetan plateau in the past. However, it was not the only driving force for asthenospheric flow in the upper-most mantle. Based on stacked   10 Lithosphere teleseismic P-wave receiver function data, Duan et al. [123,124] also proposed a cold mantle lithosphere at the mantle transition zone beneath the Qiantang terrane, implying that large-scale lithospheric detachment has most likely occurred. An alternative idea for producing asthenospheric flow is slab break-off. Ultrapotassic mafic magmatic material between 25 and 15 Ma from the last slab break-off event of the Greater-Indian slab indicate that slab break-off had a major impact on the asthenospheric flow beneath the Tibetan plateau. It involved partial melting of pyroxenites induced by the upwelling of hot asthenospheric mantle in response to an Indian slab break-off, resulting in the generation of the ultrapotassic magmas [43]. The low rigidity, hotter uppermost mantle (~300°C more), and the presence of widespread fluid in the Tibetan crust suggest a weak rheology in the central part of the Tibetan plateau and can be considered as evidence for this event. Present day tomographic studies place the detached slab from the Greater Indian slab break-off event in the mantle transition zone below the Tibetan plateau [8,18]. Even though the slab is at a greater depth today, it is still sinking. Similar to the delaminated lithospheric mantle, the sinking slab may have induced upward asthenospheric flow, as shown by a numerical model from Huangfu et al. [47]. Thus, the past asthenospheric flow beneath the central and eastern Tibetan plateau might have been induced by the combination of a delaminated lithospheric mantle and the last slab breakoff event.
There have been several models which agree with our gravity models. These models suggest that the Tibetan plateau cannot be sustained isostatically solely by a thickened crust. An underlying thin and hot lithosphere is needed to account for the present day high topography [125]. The configuration of the lithospheric mantle and asthenospheric material depicted in the gravity models ( Figures 5-9) are in favor of studies suggesting underplating in the western part of Tibet [22,52,126] with a steeper slab in the eastern part of the plateau ( [15,51,127,128]). The lithospheric thicknesses discussed in this study are in direct support of distributed crustal thickening together with the subsequent convective removal of the lithospheric mantle below the central and eastern Tibetan plateau [47,48,52]. More recent research [45] suggests rapid removal below the Qiantang and Songpan-Ganze terranes but not below the Lhasa terrane due to the absence of the Tibetan lithospheric mantle material below Lhasa's mantle transition zone [123,124,129]. There is also a fair share of disagreement between our models and other studies which suggest a completely underthrust Indian plate in the east resulting in a thicker lithospheric mantle below the Tibetan plateau (e.g.; Chen et al. [29], [45]). Another view also suggests an underthrusting of the Asian plate below the Tibetan plateau [17,109,110]. These configurations lead to a presence of a double Moho below the Tibetan plateau which has not been witnessed in several recent studies discussed above.
The presence of asthenospheric material at shallow depths below the Tibetan plateau is reminiscent of the Anatolian plateau, where crustal thinning is much more pronounced comparatively. But, similar to the Tibetan plateau, lithospheric delamination and slab break-off are major drivers for asthenospheric flow in the upper-most mantle [91]. Even though the Anatolian plateau is part of the Himalayan-Alpine collisional belt, the Early Miocene Arabian-Eurasian collision is much younger, and the aftermath presents a very different picture as compared to the Tibetan plateau.
The crustal structure of the Tibetan plateau is very well resolved. The crustal thickness below the Indian plate 11 Lithosphere increases to 80 km near the Indus-Tsangpo suture zone (ITSZ; Figures 5 and 6), north of which it decreases to 65-70 km below the Lhasa terrane. The Moho depth remains in the similar range (60-70 km) below the Qiantang and Songpan-Ganze terranes. A common trend can be noticed in the crustal thickness for the central and eastern parts of the Tibetan plateau. The crustal thickness ranges from 45 to 70 km for most parts of the dipping Indian plate. At the collisional boundary, the Moho goes to depths of 80 km as constrained by earthquake data. The thickness changes sharply to~75 km and~65-70 km below the ITSZ and the Lhasa terrane, respectively (Figures 6-8). We get similar values for the western part of the Tibetan plateau in our gravity models. All the values match the crustal thicknesses suggested by Li et al. [130]. Apart from slight differences in Moho depth around the ITSZ, we do not see crustal thinning close to the collision zone. The absence of crustal thinning could mean that delamination may be limited to the lithospheric mantle, and the crust is still intact. These values support the fact that crustal shortening beneath the southern Tibetan plateau is accommodated in the form of crustal thickening [131]. The average crustal thickness in this region (~65-70 km) is much higher than the global average. The evidence for rheologically weak layers in the gravity model is not obvious in lower crustal densities beneath central Tibet as the crustal densities are moderately high (3.05-3.15 g/cm 3 ). Due to the static nature of the gravity model, it is difficult for us to comment on the scientific studies supporting mid to low crustal flow in the north-eastern Tibetan plateau. It is possible that the flows suggested in these studies are more localized and do not fully lie under the modeled area.
Comparing the crustal structure of this study to that of previous gravity models [117,125,132], we have a very similar depth for the Moho (75 − 60 ± 5 km). It is important to note that none of the previous studies consider the crust of the steeper Indian plate to be reaching deeper than 75 km at the junction between the Greater Indian and Lhasa terranes, when evidence from earthquake data suggests otherwise in our models (Figure 9). The angle of collision within the crust in our gravity models matches the angles discussed in Bai et al. [132] except for the structure in the west. Our gravity models indicate steeper angles for crustal collision in the east (90°and 92°, Figures 8 and 9) as compared to profiles from Bai et al. [132] which can be attributed to a difference in crustal densities for the lower crust. There is also a good agreement between our model and Bai et al. [132] in the upper-middle-lower crust configuration suggested for the central Tibetan plateau (Figures 5 and 9). A 3D density structure from large-scale gravity modeling [133] shows a very similar picture as given in our western most profile, with a thicker lithospheric mantle in the west and a high-density lithospheric mantle in the middle. Comparison of our models with a 3D density structure from Kaban et al. [133] shows a slightly higher density value at 200 km depth in the west as compared to the east, which matches with our study; however, the overall values of density are slightly higher for Kaban et al. [133], with differences ranging from 100 to 150 kg/m 3 . Overall, the advantage of the gravity models presented in this study lies in the details of the collision zone and slab angle at greater depths along with the configuration of the Asian slab depicted in the gravity models in the north.

Slab Geometry and the Formation of Rift
System. The slab dip along the Indian plate, as determined from the gravity model, increases from~8°in the west to~17°in the east (Figures 6-8). This change in angle confirms a slightly steeper Indian slab east of 87°E, as suggested by Guo et al. [43]. This is also in agreement with body and surface wave tomography [13,15,22,126,127,[134][135][136][137][138], receiver function analysis [17,110,[139][140][141], and controlled-source seismic investigations [142], which suggest a lateral variation of dip angle and slab movement of the Indian lithosphere under the Tibetan plateau from west to east [13,22,143]. The Indian slab dips at~25°at 92°E (Figure 8). This change in slab angle could be attributed to a change in the style of lithospheric deformation from underthrusting to slab rollback in the central and eastern parts of the Tibetan plateau (e.g., [139,144]).
While the slab configurations in our gravity models agree with an increase in slab dip from west to east, the modern view on this subject gives us a wide range of work to compare with. There are mainly three slab configurations which are discussed in this area. The first one, which our models represent closely, resembles the slab-tearing model presented by Chen et al. [50], where the slab consistently steepens moving east [13,43,49,51,109]. The second view on the topic suggests that the slab angle flattens out after steepening in the central Tibetan plateau [123,145]. The last configuration is completely opposite to the first two suggested models with the steeper slab angle in the west and a flat subduction in the east (Xiao et al. [146]) and does not match with other results from previous studies. It is important to note that even though our models match the first configuration, the last profile discussed in this paper is at 92°E. Thus, it is difficult to comment on the slab configuration east of this profile. It is possible that the slab angle eventually flattens out after steepening. This is a potential research topic that future work can be focused upon.
The gravity model shows a sudden change in the slab dip at two distinct locations (Figures 6 and 7 between 85°E and 88°E and Figures 8 and 9 between 90°E and 92°E). The changes can be attributed to the geometry of the collision zone. East of the eastern Tibetan plateau, the geometry of the collision zone changes dramatically, as the Indian plate is subducting eastward and sinks into the mantle transition zone along the Burmese arc ( Figure 1; [13,147]), which is very different from the main Himalayan-Tibetan collision zone. Hence, on a large scale, this transition may be explained by eastward subduction and westward trench retreat of the Burmese arc during clockwise rotation of the Sunda block (e.g., [13,53]). On a local scale, Chen et al. [50] suggest that the physical properties of colliding Indian lithosphere are far from homogeneous causing kinematic differences within the plate (e.g., a change in slab dip from west to east), while Shi et al. [49] try to explain it through eclogitization of the lower Indian crust, acting as a hinge for slab roll-back, hence a steeper slab angle. Our gravity models do not show a highdensity block (>3.4 g cm -3 ) in the Indian lower crust close to 12 Lithosphere the collision zone (Supplemental material, Figure S1). This is in agreement with previous crustal scale gravity modeling in the Tibetan plateau [132]. Thus, an eclogitized lower Indian crust may not be the reason for a steeper Indian slab in the central and eastern parts of the Tibetan plateau as suggested in Shi et al. [49]. The change in slab dip of the Indian plate may be mainly attributed to eastward subduction and westward trench retreat of the Burmese arc, which could be supported by far-field forces in the region. The regions with a lateral variation (west-east) of slab angle may also represent locations of strain localization within the Tibetan plateau. Slab angle change may have resulted in segmentations along internal weaknesses [50,148]. The presence of N-S trending rifts on the surface could be related to the locations of slab dip change. The change could be accommodated through slab tearing, and the presence of these slab tears may provide a window for asthenospheric material to go up and add to the extensional forces in the region which are responsible for opening the rifts on the surface. This relation between slab angle change, rifts, and slab tearing could help in explaining how the lateral variation in slab dip is accommodated from west to east. We see the presence of the Xainza-Dingjye Rift (XDR) and the Nyima-Tingri Rift (NTR) between models at 85°E and 88°E, where the first slab angle change is observed, and the Yadon-Gulu Rift (YGR) between models at 90°E and 92°E, where the second slab angle change is seen. This reasserts the hypothesis mentioned above, where the rifts can be considered as surface manifestations of the slab tears. YGR and NTR may mark places for a change in angle from gentle to steeper, where asthenospheric material could be penetrating. Xiao et al. [146] proposes a slab-tearing model in which the Indian lithosphere has been split into two parts separated by YGR, which was developed further to apply similar reasoning along other rift zones (e.g., XDR and NTR) [13,48,53]. Because of the presence of these rifts at many places within the plateau, slab tearing and asthenospheric flow may be considered as some of the major mechanisms for the formation of these rifts zones (Xiao et al. [146,13,48,50,53]).
It is important to note that there is always an inherent ambiguity in potential field interpretations. The densities and geometric structures of the sediments, crust, lithospheric mantle, and asthenosphere in this study are well constrained using velocity models from receiver function analysis, seismic tomography, and other geophysical data available in the region. However, the interpretations of gravity models discussed in this paper can, of course, be further improved with better constraining data in the future. Adding better terrestrial gravity data and constraining models for the uppermost mantle would definitely improve the resolution of gravity models and allow us to model the deeper parts of the plateau in better detail.

Conclusions
The lithosphere of the eastern Tibetan plateau is underlain by a low-velocity zone in the upper-most mantle which has been linked with asthenospheric material in the upper-most man-tle. Our understanding of what drove this asthenospheric flow in the past towards the upper-most mantle and its spatial extent has remained limited before. To get a better understanding of the lithospheric structure and asthenospheric flow, we developed 2.5D gravity models along five N-S trending profiles. The models show lateral variations in lithospheric thickness from 150 km in the west to 110-120 km in the central and eastern Tibetan plateau. Thẽ 30 km lithospheric thinning beneath the central and eastern parts of the Tibetan plateau indicates lithospheric delamination. Based on the area covered in the profiles, we see~15,000 km 2 (N − S length = 500 km and thickness = 30 km) of asthenospheric material in the upper-most mantle between 85°E and 88°E, which could even extend east of 92°E. Based on the gravity models, this study shows that the upward asthenospheric flow in the past and hence asthenospheric underplating below the central and eastern Tibetan plateau may have been induced by lithospheric delamination and the slab break-off event of the Greater Indian slab.
Utilizing the gravity models, we found that slab angle change occurs at two distinct places and the angle changes from~8°in the west to 17°in the central part and a further 25°in the eastern part of the Tibetan plateau, suggesting that the style of lithospheric deformation changes from underthrusting to slab roll-back from west to east. Slab tears within the Indian plate could accommodate an extension caused by a slab angle change from west to east. The presence of rifts at locations of slab angle change suggests that the genesis of rifts in the eastern Tibetan plateau may be related to slab tearing and asthenospheric flow induction. Thus, the rifts can be considered as the surface manifestation of slab tears.
This study concludes that eclogitization of the lower crust is not the reason for a steeper Indian slab in the central and eastern parts of the Tibetan plateau as the densities of the lower-crustal rocks are not ≥3.4 g cm -3 . This is in contrast to studies such as those of Le Picheon et al. (Le Picheon et al., 1992) and Shi et al. [49] which support eclogitization as a reason for a steeper slab angle. On a larger scale, the change in the slab angle of the Indian plate may be attributed to eastward subduction and westward trench retreat of the Burmese arc, along with a contribution from far-field forces like slab pull and ridge push.

Data Availability
The presented data are available via FTP transfer by contacting the authors.

Additional Points
Plain Language Summary. The presence of low-velocity asthenospheric material within the upper-most mantle beneath the eastern Tibetan plateau has been known for some time now, but what has driven this flow in the past has remained unknown. Here, we assess what drove the asthenospheric flow material upwards in the upper-most mantle in the past, and delineate its spatial extent based on satellite and terrestrial gravity data modeling. The results of modeling show that the lithospheric mantle below the central and eastern Tibetan plateau is delaminated. The low-velocity asthenospheric material in the upper-most mantle may have been induced by lithospheric removal processes and consequently underplated below the eastern Tibetan plateau. Key Points. (i) Delaminated lithospheric mantle is underplated by low-velocity asthenospheric material in the upper-most mantle. (ii) Asthenospheric flow could have been induced by a combination of lithospheric delamination and a slab break-off event of the Indian slab. (iii) The Indian slab becomes steeper laterally from west to east in the Tibetan plateau.

Disclosure
This paper was presented at the AGU Fall Meeting 2019 in San Francisco.

Conflicts of Interest
There is no conflict of interest regarding the publication of this article. Submitting authors are responsible for coauthors declaring their interests.

Acknowledgments
Harshpal Singh's work is supported by the Department of Geological Sciences, The University of Alabama. The authors thank An Yin, Zuzana Alasonati Tasarova, and Sarah Roeske for their thoughtful and constructive comments, which were of great help in the preparation of the final version of this paper.

Supplementary Materials
Supplementary 1. Figure S1 End member gravity model at 92°E showing misfit in the data with eclogitized crust. The two arrows in the fitted model on the far right corner are anchors used in modeling. Black dots connected through a black line in the model marks isostatic Moho depth calculated using the Airy isostatic model. Different terranes are marked with different colors, darkening with depth. Abbreviations are as in Figure 1. Table S1. The table presents the digital data for the isostatic Moho depth (Figures 5-9). The Airy Moho model assumes a low flexural rigidity and hence presents a first order approximation to Vening Meinesz-Moritz Moho estimates. Lithosphere Lithosphere Lithosphere