Upper‐mantle velocity structures beneath the Tibetan Plateau and surrounding areas inferred from triplicated P waveforms

P‐wave waveforms in the distance range between 12° and 30° were analyzed to investigate upper‐mantle P velocity structures beneath the Tibetan Plateau and surrounding areas. The waveform data from 504 earthquakes with magnitudes larger than 5.0 between 1990 and 2005 that occurred within 30° from the center of the Plateau were modelled. We divided the study area into 6 regions and modeled upper‐mantle‐distance P waveforms with turning points beneath each region separately. The results show that the upper‐mantle P‐wave velocity structures beneath India, the Himalayas, and the Lhasa Terrane are similar and contain a high‐velocity lid about 250 km thick. The upper‐mantle velocities down to 200 km beneath the Qiangtang Terrane, the Tarim Basin, and especially the Songpan‐Garzê Terrane are lower than those in the south. The 410‐km discontinuity beneath these two terranes is elevated by about 20 km. High‐velocity anomalies are found in the transition zone below 500 km under the Lhasa and Qiangtang Terranes. The results suggest that the Tibetan Plateau was generated by thrusting of the Indian mantle lithosphere under the southern part of Tibet. Portions of the thickened Eurasian mantle lithosphere were delaminated; they are now sitting in the transition zone beneath southern Tibet and atop of the 410‐km discontinuity underneath northern Tibet.


Introduction
The Tibetan Plateau, bounded by the Tarim Basin and North China Blocks to the north and the Indian sub-continent to the south, has been created by collision between the Indian plate and the Eurasian plate (Figure 1). The collision started about 70-50 Ma ago (Molnar and Tapponnier, 1975). Since then, the India plate has continued to move northward into the interior of Eurasia. As a result, about 1,000 km of N-S convergence between India and Eurasia has taken place (Yin A and Harrison, 2000;Johnson, 2002). For collisions involving at least one oceanic plate, the convergence is always accommodated by subduction of one or more of the oceanic plates. It is still not clear how large convergence can be accommodated when both plates are continental. Various collision models have been proposed to explain the convergence between the Indian and Eurasian plates and the uplifting of the Tibetan Plateau, including 1) under-thrusting of the Indian plate beneath the Eurasian plate (Argand, 1924), 2) injection of the Indian crustal material into the weak Tibetan lower crust (Zhao WL and Morgan, 1985;Xu Q et al., 2015), 3) distributed shortening and vertical thickening of the Eurasian lithosphere (England and Houseman, 1986), 4) subduction of the Indian plate underneath the Eurasian plate in the middle of the Plateau (Owens and Zandt, 1997;Royden et al., 2008), and 5) multiple southward subductions of Eurasian mantle lithosphere underneath Tibet (Fan GW and Wallace, 1991;Tapponnier et al., 2001). Figure 2 shows simplified sketches of these models.
Each of these models predicts a unique crust and upper mantle structure. The under-thrusting model suggests a shield-like mantle lithosphere beneath the Plateau, similar to the one beneath India. The Eurasian mantle lithosphere should be removed to accommodate the advancing of Indian mantle lithosphere (Figure 2a). Injection of the Indian crust into the Tibetan lower crust would require a weak lower crust beneath the Tibetan Plateau ( Figure 2b). In the distributed shortening model, in order to raise the northern Plateau to reach an elevation of 5 km, the thickened Tibetan mantle lithosphere has to be removed by a thermal instability or small-scale mantle convection (England and Houseman, 1986) (Figure 2c). Subduction of the Indian plate in the middle of Tibet Plateau would result in a thicker crust in southern Tibet than in northern Tibet and a high-velocity subducting slab beneath the Bangong suture ( Figure 2d). Models of intracontinental subductions of both the Indian plate and Eurasian plate beneath the Tibetan Plateau are characterized by multiple slabs in the upper mantle beneath Tibet (Figure 2e). Therefore, the upper mantle structure beneath the Plateau is key to discriminating among different models.
A lot of progress has been made to reveal the seismic structures beneath the Tibetan Plateau. Seismological efforts include classical modeling waveforms in the upper-mantle distances (e.g. Zhao WJ et al., 1991;Chen WP and Tseng, 2007;Zhang RQ et al., 2011), surface-wave dispersion analyses (e.g. Yang YJ et al., 2010;Li HY et al., 2014), travel-time tomography (e.g. He RZ et al., 2010;Bao XW et al., 2013, Liang XF et al., 2016, and receiver function analysis (e.g. Wang CY et al., 2010;Zhao JM et al., 2010;Yue H et al., 2012;Pan SZ and Niu FL, 2011;Shen XZ et al., 2014). Most of these studies, however, focus on structures in the crust and uppermost mantle shallower than 100 km (e.g. Zhu LP et al., 1995;Zhu LP and Helmberger, 1998;Guo XY et al., 2013;Liu QY et al., 2014). Less at- tention has been given to deeper regions, which would also be expected to yield good indications of the tectonic history of the Tibetan Plateau. In this study, we used waveforms in the distance range of 12° and 30°, whose ray paths turn between depths of 100 km and 850 km, to obtain the upper-mantle velocity structures beneath the Tibetan Plateau and surrounding areas. The waveforms and differential travel times between different phases are very sensitive to upper mantle velocity structure. Such a study is now possible because large amounts of high-quality waveform data have accumulated over the last decade from the increasing number of seismic stations deployed in and around Tibet.

Method
Because waveforms return both arrival time and shape information from different seismic phases, waveform modeling imposes unique constraints on seismic structures (Chu RS et al., 2012a;Chu RS and Helmberger, 2014). Among these approaches, uppermantle triplications provide a powerful way to determine upper mantle velocity structures, as applied to USArray data of the United States (Chu RS et al., 2014, 2012b. As seismic velocities increase with depth, a seismic wave turns at a certain depth before being recorded by a receiver on the surface. If there is a rapid velocity increase or velocity jump, three seismic rays can arrive at the same station: 1) a ray turning above the discontinuity, 2) a ray reflected from the discontinuity, and 3) a ray turning below the discontinuity ( Figure 3a). Consequently, three phases appear on the seismogram. This phenomenon is often referred as upper-mantle triplication. In the upper mantle, there are two major velocity discontinuities, at depths of about 410 km and 660 km, both of which are due to mineral phase changes (Bina, 1991). Thus, two sets of triplication exist. For the global IASP91 velocity model and a shallow source, the 410-km discontinuity's triplication appears from 12° to 21°. The AB and CD branches cross at 18°. The triplication of the 660-km discontinuity begins at about 17° and ends at 28°. The CD and EF branches meet at around 23°( Figure 3b).
Upper-mantle triplications provide a unique way to determine upper mantle velocity structure. As shown in Figure 3a, the seismic wave of the AB branch travels nearly horizontally around its turning depth above the 410 discontinuity and is mostly sensitive to the velocity at this depth. The seismic waves of other branches turn in the transition zone and the lower mantle, where velocity structures are relatively more homogeneous laterally. Their arrivals can serve as reference time marks. Because all the rays share similar paths near the source and station, the differential travel times between the AB branch and the deeper branches are less influenced by the earthquake origin-time error and shallow velocity structural variations.
To demonstrate this, we placed a high-velocity anomaly between 200 km and 250 km in the velocity model used in Figure 3. The Pwave record section and the perturbed velocity model are shown in Figure 4a. It shows that this high-velocity anomaly changes the separations between AB and BC branches only in the distance range of 16° to 18°, with the largest separation change at about 17°. Waveforms at other distances are almost intact.
Since the differential arrival times between the AB branch and other branches are used, there is a trade-off between the velocity in the shallow upper-mantle and the deep part. For instance, a large separation between the AB and BC branches could also be explained by a low-velocity anomaly just above the 410 discontinuity or by a deeper 410 discontinuity. Figure 4b shows this trade-off by two sets of triplicated waveforms. A deeper 410 discontinuity (at 426 km) can fit the large separation in waveforms at distances between 16° and 18°. In this situation, it cannot be decided whether these waveforms are produced by a deeper 410 discontinuity or by high velocities in the top upper-mantle.
To reduce the trade-off in the waveform modeling, we used a "bottom-up" strategy to determine velocity structure in the deep part first by modeling waveforms at large distances, and then moving upward to the shallower depths. The advantage of this modeling technique is that it can reduce the trade-off and eliminate the effect of large velocity uncertainties of the shallow structures, especially the crust. We will use the phase with deep turning depths as the reference phase, on which the synthetic and the observed seismograms are aligned. Velocities near the turning depth of the shallow phase will be adjusted through trial and error to fit both the arrival time and amplitude of the shallow phase.

Data
In order to study upper-mantle velocity structures beneath Tibet, seismograms at epicentral distances between 10° and 30° are needed. To ensure that their turning points are located in the study area, we restricted earthquakes to those with distances less than 30° from the center of the Tibetan Plateau (85°E, 33°N, Figure 1). Using the Preliminary Determinations of Earthquakes (PDE) catalog, we found 1,715 earthquakes between January 1990 and February 2005 with magnitudes larger than 5.0. We then verified the Harvard CMT solutions of all events using their teleseismic P waveforms. We found 504 earthquakes with correct focal mechanisms ( Figure 1). Earthquake depths and source time functions were re-determined using a method described in Chu RS et al. (2009).
Most waveforms recorded by these broadband stations were acquired from the IRIS Data Management Center. We used a time window of 30 s before the predicted P-wave arrival time and 60 s after to retrieve the desired P waveforms from the raw data. We removed instrument responses from each waveform, re-sampled the trace with a sampling rate of 10 points per second, and bandpass filtered the data between frequencies of 0.02 Hz and 1.0 Hz.
To ensure waveform quality, we visually inspected each waveform and discarded bad traces. In total, we used 11,306 seismograms in the distance range between 10° and 30°. Their turning points cover most parts of the western Tibetan Plateau, the Hindu-Kush area, and the Hindu-Burma region. Western Tibet, eastern Tibet, and the Tarim Basin have excellent coverage by waveforms with a distance range between 12° and 20°, corresponding to turning depths between 100 to 400 km, while central Tibet and southern Tibet have less-than-ideal coverage for this distance range ( Figure 5).

Travel-Time Residuals
We first handpicked the first P arrival times in all waveforms and computed their travel-time residuals with respect to the IASP91 velocity model. In order to remove the arrival time delays caused by inaccurate event origin times and near-source structural heterogeneities, we calculated the source delay for each event by minimizing the root-mean-square of travel-time residuals between observations and predictions by the IASP91 model for the EF branch.
We averaged the travel-time residuals using a bin size of 2.5° × 2.5° × 100km. Figure 6 shows the averaged residuals for traces turning at 800 km, 500 km, 300 km, and 100 km, respectively. We placed the travel-time residuals at the turning points because the residuals are most sensitive to the velocities near the turning points. Negative travel-time residuals indicate that the velocity around the turning point is faster than the IASP91, and positive residuals suggest slow velocity anomalies. Travel-time residuals below the 400 are small throughout the study area (Figures 6a and  6b). Their magnitudes are usually smaller than 1.0 s, which indicates that the velocities in the transition zone and lower mantle beneath the whole area are close to the IASP91. In contrast, negative residuals above the 400 dominate most of the area (Figures 6c  and 6d), including India, southern Tibet, central Tibet, and eastern Tibet. The difference between the observed arrival times and predictions can be as large as -7.0 s, which means that the velocities above the 410 are faster than the IASP91 upper-mantle velocities in these areas. Positive travel-time residuals, which suggest lower velocities, are mainly concentrated between 200 km and 300 km beneath the southern Tarim Basin, western Tibet, and southeastern Asia ( Figure 6).

Waveform Modeling
Based on the major fault and suture systems, we divided the area along the NS direction in Figure 7 into six regions: India, the Himalayas (HL), the Lhasa Terrane (LS), the Qiangtang Terrane (QT), the Songpan-Garzê Terrane (SG), and the Tarim Basin (TB) (Figure 7). Waveforms with turning points beneath each region were grouped for modeling to obtain the best 1-D upper-mantle velocity model for the region. To ensure good data quality, we inspected each waveform in the group and kept only those with a clear first P arrival. We then used the verified Harvard CMT solutions and revised focal depths to compute synthetic seismograms. This helped identify other late phases. Since the modeling technique relied on using a deep-turning phase, such as the CD and EF branches, as the reference phase to align the data and synthetic waveforms, only traces with at least one clear deep-turning phase were used for the modeling.
Waveforms in each group were from different earthquakes at various depths, which made the waveform modeling more difficult. Figure 8a shows synthetic P waveforms from two events, one at a depth of 30 km and the other at 10 km. If we align these two profiles at the first arrivals directly, the late phases are misaligned at a given epicentral distance. In order to construct a single epicentral distance waveform profile for modeling, it is necessary to correct for the source-depth effect because waveforms from a deep earthquake have larger "effective" epicentral distances on the distance profile than those from a shallow event. We chose 10 km as the reference source depth because most of the events used in this study were shallow earthquakes with focal depths at 4-30 km (Chu RS et al., 2009). Differences between the true and effective epicentral distances were estimated using the ray parameter of the first arrivals and the background velocity model. After correcting for the source-depth effect, all seismograms can be placed on a single epicentral distance profile for modeling (Figure 8b).

India and the Himalayas
There were fewer usable traces with turning points beneath India compared to other regions, mainly because there are fewer available stations and less seismicity in the stable continent. Therefore, we treated all of India as one unit. In order to minimize the effect of the upper-mantle structure beneath the collision front, we chose those paths with turning points at least 200 km away from the MBT. A total of 12 traces with distances ranging from 16° to 28° were selected. Figure 7 shows locations of event-station pairs and their turning points used for obtaining velocity structure beneath India.
We first analyzed 4 traces at distances of 22° to 28° to constrain the velocity structure between the 660 and 410 discontinuities. Depths of the turning points, which provide a good coverage of the transition zone, are shown in black circles. At this distance range, only CD and EF branches and their depth phases are clear, as indicated by dashed lines (Figure 9a). Since the EF branch turns below the 660 discontinuity, we choose it as our reference phase. The observed data and synthetic seismograms are aligned along the EF branch. Observed CD branch agrees with the IASP91 prediction, and synthetic seismograms from the IASP91 model match the observations well for these four traces (Figure 9a), which suggest that the velocity structure below the 410 is the same as the global average model.
At shorter distances between 16° and 21°, the EF branches have small amplitudes and are hard to identify. Fortunately, use of CD branches larger than 20° can resolve the velocities in the transition zone very well (Figure 9b). Therefore, we chose the CD branch as our reference phase to align the 8 remaining waveforms. The separations between the AB and CD branches in the observations are about 5.0 s larger than the IASP91 predictions ( Figure 9a). The early AB branch suggests high-velocity anomalies in the upper mantle above the 410 discontinuity. We estimated the turning depths of these waveforms using the IASP91 model and adjusted velocities around the turning depths until the synthetics matched the data. Figure 9a shows the final waveform fits. The best-fit velocity model has a thick high-velocity lid to a depth of 250 km with a P-wave velocity of 8.41 km/s (Figure 9b).
For the Himalayas, we show 11 waveforms with distances between 14° and 26°, whose ray paths sample depths from 100 km to 660 km (Figure 9). The same modeling strategy was used to obtain upper mantle velocities. Similar to the structure beneath India, a high-velocity layer is obtained at the depths of 70 km to 230 km beneath the Himalayas with waveform fits shown in Figure 9c and velocities shown in Figure 9d.

The Lhasa Terrane
In Figure 10a we show a profile of 7 seismograms used to constrain the upper-mantle velocity structure beneath the Lhasa Terrane. Three waveforms at distances between 23° and 26°s ample the lower part of the transition zone ( Figure 10b). The observed CD and DE branches are 1.0-1.5 s earlier than the IASP91 predictions and the crossover distance of CD and EF is larger than that of IASP91. Since the CD branch is the P waves turning above the 660 discontinuity and the DE branch consists of P waves reflected from the discontinuity, the observations suggest a highvelocity anomaly in the transition zone. Modeling of the waveforms yields a velocity increase of 2.5% between 550 km and 660 km and the 660 discontinuity is depressed about 15 km, as shown in Figure 10b. We verified the model using P waveforms from another earthquake at this distance range and compared them with synthetic seismograms calculated using IASP91 (Figure 10c) and the best-fit model (Figure 10d). It is clear that IASP91 synthetics do not fit the data while the new model predicts correctly the close separation between the EF and CD/DE branches.
The AB and CD branches of the next two seismograms at the distances of 22° and 18° sample above and beneath the 410 discontinuity, respectively (Figure 10b). Their separations from the EF branches are matched by the IASP91 predictions, yielding no velocity perturbation at this depth range (Figure 10a). The remaining two P waveforms at closer distance range, however, have the AB branch 5.0 s earlier than the IASP91 predictions, when the CD branch is used as the reference. This corresponds to a high-velocity anomaly of about 5% at depths 70-230 km beneath the Lhasa Terrane, similar to those beneath India and the Himalayas (Figure 10b).

The Qiangtang Terrane
Because of its location and earthquake-station distribution (Figure 1), the Qiangtang Terrane has the best path coverage. The Terrane's CD-EF separation at distances larger than 23° is smaller than the prediction of IASP91 (Figure 11a and 11c). Modeling waveforms at these distances therefore also resulted in high-velocity anomaly in the transition zone, similar to beneath the Lhasa Terrane. Because the CD-EF separation between 20° and 23° agrees with the IASP91 prediction, we constrained the high-velocity anomaly to the bottom half of the transition zone.
The 410 triplication location is further than the IASP91 prediction and the B-cusp of the AB and BC branches extends to 25° as opposed to 21° predicted by IASP91 (Figures 11a and 11c). The broader BC branch requires a larger velocity jump across the 410 discontinuity. Since the velocity value below the 410 discontinuity has been fitted by modeling of the CD branch in waveforms between 20° and 23°, we had to reduce the velocity above the 410 discontinuity by about 2.0%. This velocity reduction will significantly delay arrival times of AB branches at larger distances (20°a nd 26°). The observed data, however, do not show this velocity reduction (Figure 11a). Therefore, we move the 410 discontinuity 20 km shallower, to 390 km (Figure 11b). Using the CD branch as the reference phase, the AB branch is about 3.0 s earlier than the IASP91 predictions. This requires a high-velocity layer in the upper mantle down to 270 km beneath the Qiangtang Terrane. The magnitude of the velocity perturbation, however, is smaller than those beneath India, the Himalayas, and the Lhasa Terrane. Between 13° and 18°, the AB branch has an extra set of arrivals, as indicated by the green lines (Figures. 11a). A small discontinuity of 2.5% velocity jump at the depth of 270 km is needed to explain this minor triplication of the AB branch.

The Songpan-Garzê Terrane
The separation of CD and EF branches agrees with that of the IASP91 for waveforms at distances larger than 19° (Figure 12a), which suggests that the transition-zone velocity structure beneath the Songpan-Garzê Terrane fits the IASP91 model. No transition zone high-velocity anomaly was detected. The Songpan-Garzê Terrane has a 410 discontinuity velocity structure similar to that of the Qiangtang Terrane, with a shallow 410 and a reduced velocity above it. This is also due to a broader BC branch in the observed waveforms and early AB branch in the distance range of 20° to 22° shown in Figures 12c and 12d for a 21.0 km deep earthquake in the Hindu-Burma region.
For waveforms between 12° and 17°, arrivals of the AB branch are about 2.0 s earlier than the IASP91 predictions. These early arrivals suggest a high-velocity layer in the uppermost mantle, similar to that beneath the Qiangtang Terrane. In order to fit the relatively later AB arrival near 13°, a low-velocity layer between depths of 100 to 180 km is needed (Figure 12b).

Tarim Basin
The Tarim Basin upper-mantle velocity structure was obtained by modeling 10 waveforms with turning points beneath the region (Figure 7). The waveform fits of large distances beyond 19° indicate that the transition-zone velocity structure is the same as the global average model. However, the AB branch is about 1.5 s earlier than predicted by the IASP91 model. This requires a high-velo-    Figure 13). This high-velocity layer is similar to those beneath the Qiangtang Terrane and Songpan-Garzê Terrane.

Discussion
In this study, waveforms in the upper-mantle distances were fitted using a band-pass filter between 0.02 Hz and 1.0 Hz; the main frequency content is about 0.5 Hz. Thus, the wavelength of the signals is about 20 km. Using the quarter-wavelength rule, we estimated that the waveform data had a vertical resolution power of less than 5 km.
The horizontal resolution of our results depended on station and earthquake distributions. In regions such as the Lhasa and Qiangtang Terranes, which were sampled by abundant seismic rays, waveforms of turning points separated by less than 4° were grouped together for modeling. The corresponding horizontal resolution was 400-500 km. In other regions like India and the Tarim Basin where data were scarce, the whole tectonic block was treated as one unit.
Uncertainties of the velocity models come from several sources, including failure to consider lateral velocity structural variations in the crust, earthquake epicenter mis-locations, and earthquake depth errors. Although all the focal depths of the earthquakes used were determined using the differential travel times between the direct P and the depth phases at teleseismic distances, a depth uncertainty of 2-3 km exists since a 1-D velocity model was used. To estimate the effect of source depth uncertainties on the modeling results, we calculated travel-time triplications from a 10km-deep earthquake and a 15-km-deep earthquake. After removing the arrival-time delay due to different source depths, the two triplications are almost identical. This indicates that a difference of 5 km in event depth causes less than 0.1 s in differential arrival times. Such a small error can be safely ignored.
The waveform modeling results show that a high-velocity lid exists beneath India, the Himalayas, and the Lhasa Terrane. The lid velocity beneath the Tarim Basin is relatively slower than these Terranes, but still faster than those beneath the Qiangtang and the Songpan-Garzê Terranes. Both India and Tarim are Archean cratons and are believed to lie over a cold and strong mantle lithosphere. Travel-time tomography studies suggested a high-velocity anomaly down to 250-300 km beneath India (Huang JL and Zhao DP, 2006;Li C et al., 2008;Zhang H et al., 2012). A thick highvelocity upper mantle lid has also been suggested by Lyon-Caen (1986)  Existing descriptions of the upper-mantle velocity structure beneath the Tarim Basin are controversial. Surface-wave dispersion and some travel-time tomographic results have indicated that upper-mantle velocities beneath the Tarim Basin are slow (Li C et al., 2008;Lebedev and Van Der Hilst, 2008), while other travel-time tomography studies have suggested high velocities (Liu M et al., 2004). Yao HJ et al. (2005) used inter-station Rayleigh-wave dispersions and found lateral velocity variations beneath the basin. The northeastern Tarim had higher upper-mantle velocities than the southwestern Tarim. Our waveform modeling did not reveal such a lateral variation. The 8.22 km/s P velocity at the top of the mantle agrees with the Pn velocity from the travel-time tomographic result (Liang CT and Song XD, 2006).
The upper-mantle velocities above 200 km beneath the Qiangtang and Songpan-Garzê Terranes are lower than in its north and south, especially beneath the Songpan-Garzê Terrane where a low-velocity zone exists ( Figure 14). These results agree with previous Sn and Pn studies that showed a region of inefficient Sn propagation and low Pn velocities in northern-central Tibet (Lü Y et al., 2011). Cenozoic volcanism has also been discovered in this region (Gansser, 1980).
Since the upper-mantle velocity structures beneath the Himalayas and the Lhasa Terrane are similar to those beneath India, we suggest that the Indian mantle lithosphere has thrusted beneath the Tibetan Plateau as far as to the Bangong-Nujiang Suture. Based on these results, an upper-mantle structural cross-section along a N-S profile across the Tibetan Plateau can be sketched ( Figure 15). Prior to the collision between India and Eurasia about 50 Ma ago, the Eurasian crust had been shortened from 2000 km to 1000 km, implying a stretching factor of 50%. The uppermantle portion did not change because the collision involved oceanic plates whose upper mantle was hot. With the marching of India northward, the Eurasian upper mantle underwent compression between two continents. Its length changed from 1000 km to 650 km, a stretching factor of 65% . In short, the 1000 km N-S convergence between India and Eurasia is accommodated by horizontal shortening of the Eurasian crust and mantle lithosphere with stretching factors of 50% and 65%, respectively ( Figure 15).
The high velocity anomalies in the transition zone beneath the Lhasa and Qiangtang Terranes are interpreted as the remnants of the delamination of the Eurasian mantle lithosphere at about 15 Ma ago, as suggested by Chen WP and Tseng (2007). This rem-

Conclusions
By fitting P waveforms in the distance range between 12°and 30°, we obtained upper-mantle velocity models beneath six regions along a N-S profile across the Tibetan Plateau. The results show that the upper-mantle velocity structures beneath India, the Himalayas, and the Lhasa Terrane are similar and contain a high-velocity lid about 250 km thick. The Tarim Basin also lies above a highvelocity upper-mantle lid. The upper-mantle velocities deeper than 200 km beneath the Qiangtang and Songpan-Garzê Terranes are lower than those in the north and south, especially those beneath the Songpan-Garzê Terrane. A high-velocity anomaly can be seen beneath these two terranes between 270 km and 410 km, and the 410 discontinuity is elevated by 20 km. We also found high-velocity anomalies in the transition zone below 500 km under the Lhasa and Qiangtang Terranes. The results suggest that the Tibetan Plateau was generated by the thrusting of the Indian mantle lithosphere under the southern part of Tibet. The India-Eurasian convergence is a two-phase process. Two portions of the thickened Eurasian lithosphere were delaminated in southern and northern Tibet. They are now sitting atop the 660 and 410 discontinuities, respectively.