Intracontinental deformation of the Tianshan Orogen in response to India-Asia collision

How the continental lithosphere deforms far away from plate boundaries has been long debated. The Tianshan is a type-example of ongoing lithospheric deformation in an intracontinental setting. It formed during the Paleozoic accretion of the Altaids and was rejuvenated in the Cenozoic, which might be a far-field response to the India-Asia collision. Here we present seismic images of the lithosphere across the central Tianshan, which were constructed from receiver functions and Rayleigh wave dispersions along a N–S-trending linear seismic array. We observe an extensively deformed lithosphere in the Tianshan with inherited, structurally controlled brittle deformation in the shallow crust and plastic deformation near the Moho. We find that earlier multiple accretionary structures were preserved in the crust, which was deformed by pure-shear shortening in the south and thick-skinned tectonics in the north but was limitedly underthrusted by surrounding blocks. A balanced cross-section of Moho discontinuities supports the concept that intracontinental deformation in the Tianshan intensified synchronously with the direct contact between the underthrusting Indian slab and the Tarim Craton in the Late Miocene (~10 Ma). These findings provide a robust and unified seismic model for the Tianshan Orogen, and confirm that effective delivery of the India-Asia collision stress induced the rejuvenation of this intracontinental orogen.

O rogens, viewed as a result of rigid plate collision, are largely confined to plate boundaries 1 , whereas orogens inside continental plates are referred to as intracontinental orogens 2 . The Tianshan Orogen (Fig. 1), situated in the southern Altaids 3,4 , records long-lived accretion and continental growth in the Paleozoic [3][4][5] . In the Cenozoic, this orogen was rejuvenated and re-deformed by the India-Asia collision that was thousands of kilometers away [6][7][8] . Hence, the Tianshan Orogen is an ideal site to investigate ongoing intracontinental deformation.
Numerical models have suggested that rheological heterogeneity of the lithosphere controlled the intracontinental deformation of the Tianshan 9-11 , but it remains debatable how the deformation was accommodated across the orogen and what the tectonic trigger was. The Cenozoic deformation of the Tianshan has been related to the progressive growth of the Tibetan Plateau, which created the high gravitational potential energy necessary to transmit the pressure northwards 5,8,9 . Thermochronologic and stratigraphic studies indicate that the intracontinental deformation of the Tianshan commenced at a low strain in the Early Miocene and significantly increased at~10 Ma [12][13][14] , which was much later than the initial continental collision between India and Asia at~55 Ma 15,16 . Another model links the Tianshan to the transpressional system of the right-slip shears induced by the oblique convergence between Arabia and Asia 16 . Despite some of the estimated age of continental collision between Arabia and Asia being close to the uplift of the Tianshan 17 , it remains unclear whether the soft continental collision suggested there 18     the spatial and temporal development of the structures across the Tianshan could significantly improve the understanding of how and where the continental lithosphere responds to the far-field plate collisions. Geophysical investigations support the hypothesis that the upper mantle beneath the Tianshan is so weak [19][20][21] that it accumulated strain from the distant but ongoing India-Asia collision. As a prominent intracontinental deformation orogen, the Tianshan has accommodated~20-mm/yr N-S shortening, which is approximately half of the present-day India-Asia convergence rate 22,23 . Previous seismic images seem to suggest that the N-S shortening was accommodated by large-scale underthrusting of the Tarim and Kazakh lithospheres 20,[24][25][26][27] . Had such a scenario existed, most of the N-S shortening would have occurred near the northern and southern discrete boundaries, as is the case in the Himalayas where the Indian slab is underthrusting beneath the Tibetan Plateau (Fig. 1a). However, geodetic and neotectonics observations indicate that the shortening in the Tianshan is distributed throughout the orogen without a significant increase close to the boundaries 22,23,28 (Fig. 1a).
The Cenozoic N-S shortening of the Tianshan shows an alongstrike variation and decreases from west to east 22,23 . The central Tianshan, with the western boundary at the Talas-Ferghana Fault, is adjacent to the sub-horizontal underthrusting Indian slab beneath the Pamir Plateau and the West Kunlun [29][30][31][32][33] (Fig. 1a) and experienced the strongest Cenozoic N-S shortening along the entire orogen 23,34 . At the eastern, topographically low, termination of the Tianshan, the horizontal shortening in the Late Cenozoic was only 10-15 km and the deformation changed into the left-lateral transpressional type 8 .
Here, we studied the seismic structures of the crust and mantle lithosphere beneath the central Tianshan to provide a comprehensive image of the ongoing intracontinental deformation. The seismic cross-section of this study is close to the western end of the Tianshan Orogen, which is tectonically subdivided, from the north to the south, into the Issyk Kul Arc (IKA), Naryn Arc (NA), and South Tianshan Accretionary Complex (STAC), separated by ophiolitic mélanges and/or large-scale faults 4 ( Fig. 1). We collected data mainly from a N-S aligned linear broadband seismic array of the Middle AsiaN Active Source project (MANAS) 25 (Fig. 1), and applied the joint inversion scheme of receiver functions and Rayleigh wave dispersions 35 to construct a S-wave velocity (V S ) model, and the common conversion point (CCP) stacking technique 36 to image the geometry of the Moho and intracrustal interfaces. The results exhibit a spatial relationship between the surface deformation and the deep structures of the Tianshan, and further illustrate the mechanism and essential conditions for the intracontinental deformation.

Results
The V S model and CCP stacked receiver function image constructed in this study are shown in Fig. 2. The cross-section passes through different tectonic units of the Tianshan Orogen 4 , i.e., IKA, NA, and STAC. The V S model reveals a widespread low-V S middle-lower crust (~3.4 km/s) beneath the STAC, and a series of elongate low-V S anomalies (~3.5 km/s) in the middle crust confined to several bands beneath the NA and IKA (Fig. 2d). These bands of low-V S anomalies correlate well with the negative Ps phases in the CCP image, the southernmost of which can be traced as a clear interface that dips to the south at~30°beneath the NA (Fig. 2c).
In the CCP image (Fig. 2c), the Moho can be traced by strong positive Ps phases. Along the cross-section, the Moho has been thrust-imbricated into several segments beneath the entire Tianshan Orogen and a doublet beneath the Nikolaev Line (NL).
The deep Moho segments are located at~70 km depth below the STAC and at~60 km depth below the IKA, whereas shallow segments are found at 40-50 km depth below the Tarim Craton, NA, and Kazakh Shield. The Moho thrusts appear to correlate with the sub-Moho Vs variations. Contrasting structures in the uppermost mantle are shown from south to north across the NL. A 40-km thick high-V S mantle lid (~4.5 km/s) lies above the low-V S anomaly (~4.1 km/s) beneath STAC and NA, whereas a thick high-V S anomaly (~4.5 km/s) spreads over the entire depth range beneath IKA and Kazakh Shield in the north and beneath Tarim Craton in the south (Fig. 2d).

Discussion
Inherited properties control shallow crustal shortening. Along the profile, a number of active thrust faults can be identified by large slip rates (>2.0 mm/yr) 28 (Fig. 2b) and significantly increased strain rates 37 (Fig. 2a). In our crustal model, we observe several south-dipping low-V S bands, which extend continuously from the shallow to middle-lower crust (Fig. 2c, d) and correlate well with active faults at the surface. These bands represent crustal shear zones with weak rheology at depth 38 . Given the position and geometry of these crustal structures, most of the shallow crustal shortening in the NA and IKA was likely accommodated by thick-skinned tectonics, which is consistent with the intense seismicity near the shear zones with high-angle thrust focal mechanisms 39 . The high and flat topography of the STAC and the strongly depressed Moho below it are indications of extensive crustal shortening. However, the fault slip rate observed in the STAC is extremely small (~0.3 mm/yr) (Fig. 2b), implying that the STAC underwent a pure-shear shortening with diffusive deformation accommodated by the rheologically weak crust. This model is supported by the stable low stain rate (~20 × 10 −9 /yr, Fig. 2a) and lack of seismic activity (Fig. 2d).
These differences in shallow crustal deformation structures correlate well with the locations of tectonic terranes accreted in the Paleozoic, i.e., the pure-shear shortening in the STAC contrasts with the thick-skinned tectonics in the NA and IKA ( Fig. 2a-d). Magnetotelluric studies also show that the crust of the central Tianshan consists of isolated high resistivity bodies separated by conductive zones in the north, but low resistivities in the southern part adjacent to the Tarim Craton 40 . The STAC was generated by the long-lived consumption of the South Tianshan Ocean in the Paleozoic. It contains high/ultrahigh-pressure rocks, ophiolitic mélanges, and imbricated sedimentary rocks, which are interpreted as off-scraped fragments 4 (Fig. 2b). Such an accretionary wedge was too weak to deliver the stresses and in consequence experienced pure-shear shortening in response to the India-Asia collision, whereas the crust of the Paleozoic volcanic arcs (NA and IKA) to the north was so strong 4,28 that active faults were able to penetrate into the deep crust to accommodate the contractional stress ( Fig. 2b-d).
Moho deformation responses to deep crustal shortening. Previous models have proposed large-scale underthrusting of the Tarim Craton and/or the Kazakh Shield beneath the Tianshan 24-27 . Our results demonstrate that the underthrusting of the lower crust is only limited in the Tarim Craton and even negligible in the Kazakh Shield, instead, prominent Moho steps are visible near the southern and northern boundaries of the Tianshan that were accompanied by the fade of seismic activities in the lower crust ( Fig. 2c, d, see details in the Supplementary Text 1). Another accordance of the underthrusting model is the high-velocity anomalies in the upper mantle of the Tianshan observed in the teleseismic body wave tomographic images 20 , which were interpreted as a deeply subducted lithospheric mantle.   However, a large amount of shortening (>500 km, according to the scale of anomalies in Lei & Zhao 20 ) expected in the interpreted model is difficult to reconcile with the short rejuvenation history of the Tianshan (~10 Ma) [12][13][14] , even with the current high contractional rate (~20 mm/yr) 22,23 .
Based on the almost linear decrease of GPS velocities across the Tianshan, England and Molnar 41 argued that the Tianshan lithosphere has a spatially invariant strength controlled by its low-temperature plasticity, which resulted in plastic deformation around the Moho that has accommodated the stress in the ductile lithosphere. Interestingly, our refined crustal model shows a significant discontinuous Moho with thrust structures across the Tianshan (Fig. 2c, d). Similar Moho discontinuities were imaged in the eastern Tianshan 42,43 , although the lateral variation along the strike of the orogen is unclear. It should be noted that the peneplanation in the Tianshan is generally inferred to have developed during the Mesozoic tectonic quiescence 44 after the Paleozoic accretion of the Altaids 4,5 . Accordingly, the Moho can be used as a marker horizon for estimating the amount of N-S shortening across the Tianshan in response to the India-Asia collision in the Cenozoic. As shown in Fig. 2e, we estimate that the Cenozoic N-S shortening was~85 km by balancing the Moho length and~130 km by balancing the crustal volume across the Tianshan, respectively. It can be inferred that the significant N-S shortening started at~11-6 Ma by combining the present-day N-S shortening rates of~15 mm/yr from GPS observations (Fig. 2a) and~12 mm/yr from Late Quaternary fault slip rates (Fig. 2b) across the estimated section. Such an inferred time is similar to the onset of rapid uplift of the Tianshan (~10 Ma) [12][13][14] (Fig. 1a), which confirms that the deep deformation of the central Tianshan is performed as the Moho segments are discontinued by localized thrusts. Therefore, we propose that the plastic deformation close to the Moho mainly accommodated the N-S shortening in the deep lithosphere of the central Tianshan, and that the underthrusting of the Tarim Craton or the Kazakh Shield should be on a limited scale. This deformation model correlates well with the shallow crustal structures revealed in this study and with the comprehensive observations from geodesy and neotectonics, and supports the concept that intracontinental deformation in the Tianshan intensified in the Late Miocene (~10 Ma).
Rejuvenation of the Tianshan Orogen. Our images reveal that the lithosphere of the Tianshan was extensively deformed in response to the Cenozoic India-Asia collision. The question remains as why the amount of the India-Asia convergence accommodated in the Tianshan with the N-S shortening of 20 mm/yr is nearly equivalent to that in the Tibetan Plateau 22,23 . The former should be much lower considering the lithospheric strength of the Tianshan is an order of magnitude greater than that of the Tibetan Plateau 41 . After the break-off of Greater India (i.e., the northern passive continental margin of the Indian plate, as suggested by plate reconstructions and seismic tomography studies 45 ), the India-Asia convergence changed from continental subduction to continental collision due to the buoyancy of the Indian Craton with its thick and depleted lithospheric mantle 46,47 . Although the Tibetan Plateau and the Indian Craton are now abutted along the Indus-Yarlung ophiolite belt at the surface, seismic images suggest that the cratonic Indian slab has underthrust sub-horizontally farther northwards for several hundred kilometers beneath the western Tibetan Plateau 31,48 . Previous receiver function studies confirm that the front of the cratonic Indian slab has reached the Tarim Craton in the West Kunlun [29][30][31] and thrust under the Pamir Plateau northwards to the intermediate-depth seismic zone 32,33 (Fig. 1a). The synchronous foundering of the deep crust beneath the Pamir Plateau at~1 1 Ma, documented by the xenoliths in the southeastern Pamir and related to the indentation process of the cratonic Indian slab beneath the Pamir Plateau 49 , provides an age record for the initial contact between the Indian slab and the Tarim Craton. The balanced cross-section analysis of the Moho discontinuities in this study (Fig. 2e) and the onset of rapid exhumation in the Tianshan 12-14 jointly support the concept that the intracontinental deformation in the Tianshan intensified after the Late Miocene (~10 Ma) synchronously with the India-Tarim collision. Therefore, we propose that the direct contact between the cratonic Indian slab and the Tarim Craton enabled the compressional stresses to reach the distant intracontinental setting 50 and intensified the rejuvenation of the Tianshan Orogen. Such a stress transformation in the upper mantle was a significant contributor to the slowing-down of the India-Asia convergence 51 and to the clockwise rotation of the Tarim Craton 34 , but also can well explain why the Tianshan was rejuvenated much later than the initial India-Asia collision 11 .
Previous studies have also reported low-velocity anomalies in the upper mantle beneath the Tianshan 19-21 , implying a high temperature and weak rheology there, which were likely caused by lithospheric detachment 52,53 . Small-volume Cretaceous-Paleogene basalts scattered in the Tianshan 54 and contemporaneous highheat flux (∼80-85 mW/m 2 ) inferred from xenoliths 55 suggest that the uppermost mantle in the Tianshan should have been hot and weak before the Cenozoic shortening. Although our Vs model has a limited resolution in the mantle depth, it does image upper mantle low-V S anomalies overlain by a thin mantle lid extending from the Moho to 80-90-km depth beneath the STAC and the NA (Fig. 2d). A similar mantle lid structure is also indicated in a previous S-wave receiver function study 56 . Such a pre-weakened lithosphere of the Tianshan weaker than its surroundings is plausible to accommodate the vast contractional stresses 2,10,57 . The thin mantle lid also might be beneficial for the Moho thrust deformation, as indicated by the Moho doublet developed in the northern edge beneath the NL (Fig. 2c, d).
Therefore, with the direct contact between the cratonic Indian slab and the Tarim Craton, the India-Asia collisional stress could effectively reach the pre-weakened Tianshan and contribute to its rejuvenation (Fig. 3). Our refined crustal model well demonstrates that inherited properties in the crust controlled shallow brittle deformation as indicated by the pure-shear shortening of the accretionary wedge (STAC) and the thick-skinned deformation of the Paleozoic volcanic arcs (NA and IKA). The plastic deformation accommodated N-S shortening in the deep and thrust-imbricated the Moho into several segments across the central Tianshan. The balanced cross-section of Moho discontinuities supports the idea that intracontinental deformation in the Tianshan intensified synchronously at the time of the direct contact between the cratonic Indian slab and the Tarim Craton in the Late Miocene (~10 Ma). These findings provide a robust and unified rejuvenation model for the Tianshan Orogen, and lead to promising implications for depth-dependent deformation in intracontinental orogens worldwide.

Methods
Quality control and harmonic analysis of receiver functions. We collected data from the MANAS network (40 stations) 25 and three other stations close to the MANAS network to construct a~475-km long cross-section with a~10-km station space across the central Tianshan (Fig. 1a, b). We selected events with Mw ≥ 5.0 at an epicentral distance of 30-90° (Fig. 1c) to calculate P-wave receiver functions (RFs) with Gaussian coefficients of 1.0, 2.0, and 3.5 by using the time-domain iterative deconvolution method 58 . After normalizing both time and amplitude of the RFs to a reference slowness of 6.4 s/°(moveout correction 59 ), all corrected RFs for each station were stacked in bins of 30°back-azimuth within a delay time window of −5-15 s that include Ps phases converted at the Moho (Pms) but exclude the subsequent multiples (PpPms and PsPms + PpSms). We selected reliable RFs according to the quality control procedure proposed by Shen et al. 60  (3) that were significantly different from the corresponding bin-averaged RF (cross-correlation coefficient <0.5). After the quality control process, a reliable, successive RF dataset was established ( Supplementary Fig. 1a-c).
We applied the harmonic analysis method 60 to extract the azimuthally independent RFs with minimal potential influences of tilting interface and azimuthal anisotropy ( Supplementary Fig. 1d-h). The piercing points of the Ps phases at a depth of 60 km sample a circular area with a diameter of 0.5°centered on each station (Fig. 1b), similar to the horizontal grid scale of the Rayleigh wave dispersion (DPs) (0.5°× 0.5°grids) used in this study 61 . However, the harmonic analysis method may be less effective when the azimuthal coverage of RFs is poor, due to the limits of the distribution of events and data quality. To minimize the bias in estimating the average RF for stations with uneven azimuthal coverage, we calculated bin-averaged RFs within a bin of 30°back-azimuth, and then stacked them to estimate the average RFs and corresponding uncertainties. Subsequently, the average RFs were used in the joint inversion with the DPs. Supplementary  Fig. 2b-d present the average RFs with Gaussian coefficients of 1.0, 2.0, and 3.5 for each station, and the Pms phases are traceable continuously along the cross-section.
Joint inversion of receiver function and surface wave dispersion. Based on the results of the ambient noise tomography (0.5°× 0.5°grids) in the Tianshan and the Pamir Plateau 61 , we interpolated group (6-40 s) and phase velocity (6-42 s) DPs at each station ( Supplementary Fig. 2e, f). For the central Tianshan, checkerboard tests indicate that group and phase velocity DPs can well resolve anomalies with a size of 0.5°× 0.5°in the shorter periods and 1.0°× 1.0°in the longer periods ( Supplementary Fig. 3a), and posterior errors are <2% for all periods (Supplementary Fig. 3b). In this study, we performed a joint inversion of the interpolated DPs and average RFs to obtain a one-dimensional (1-D) V S model at each station using the Computer Programs in Seismology (CPS, version 3.30) software package 62 that is based on an iteratively damped least-squares algorithm 35 . The V S was set to a constant value of 4.5 km/s in the initial model. We used the smoothed crustal V P /V S ratio for each station by averaging the recent H-κ-c stacking results 27 within a distance of 50 km from each station, and set the V P /V S ratio in the upper mantle corresponding to the AK135 model 63 . The density (ρ) was calculated as a simple approximate from the Nafe-Drake empirical relationship: 64 ρ = 1.73 + 0.2V P . The recovery tests using a similar dataset suggest that our datasets and inversion schemes can obtain a robust result down to 120-km depth 61 .
As shown in Supplementary Fig. 4, the final 1-D V S models are well consistent with the interfaces, especially the Moho, as indicated along raypaths in the individual RFs plots and by the single station CCP stacking. The synthetic RFs and DPs calculated from the final 1-D V S models of the joint inversion fit well with the observed ones. To evaluate the errors of the final model, we used resampled DPs and RFs as input to the same inverse procedure described above for 200 times. The resampled DPs and RFs were generated by adding Gaussian random noise with a standard deviation corresponding to the posterior error (Supplementary Fig. 3b) and the difference between the raw RFs and harmonic RFs ( Supplementary   Fig. 1d), respectively. The average values of the 200 V S models are similar to the final 1-D V S model with small uncertainties (standard deviation <1%) at depths shallower than 120 km (Supplementary Fig. 4).
The 1-D V S models at all stations inverted by the above-mentioned scheme ( Supplementary Fig. 2g) were assembled to construct the final two-dimensional (2-D) V S model along the cross-section. A most prominent feature resolved by our V S model is the contrasting crustal structure across the Tianshan orogen, i.e., an extensive low-V S (~3.4 km/s) middle-lower crust beneath the STAC, whereas low-V S anomalies (~3.5 km/s) in the middle-lower crust confined within several bands beneath the NA and the IKA (Fig. 2d). Three elongate low-V S anomalies, which can be easily recognized in the 2-D V S model beneath stations KARD, BESM, and SOUR (Fig. 2d), are located at depths of~20-50 km as seen in the 1-D V S models, respectively ( Supplementary Figs. 4b-d). To test the robustness of the low-V S anomalies, we modified the 1-D V S models of stations KARD, BESM, and SOUR, by excluding the crustal low-V S anomalies, and calculated synthetic RFs and DPs. We find significant differences between the synthetic and observed DPs, and the synthetic RFs from the modified models fail to fit the negative Ps phases at~3-4 s ( Supplementary Fig. 5a-c). These large mismatches between the observed and calculated RFs and DPs indicate the reliability of the low-V S anomalies in our V S model. We also performed robustness tests of the mantle lid structure imaged beneath the STAC and the NA (Fig. 2d). Both synthetic DPs from the final and the modified 1-D V S models are similar to the observed ones, but the synthetic RFs from the modified models, by excluding the high-V S mantle lid at depths of~50-80 km, fail to fit the Ps phases at 6-8 s ( Supplementary Fig. 5d, e). These features indicate that the mantle lid structure is also reliable although DPs used in this study have a limited resolution of the mantle structure.
Common conversion points stacking of receiver functions. Referring to the final 1-D V S model at each station, we projected the reliable RFs into the depth domain, and use the CCP stacking technique 36 to image the geometry of the Moho and intracrustal interfaces. Due to the V P /V S ratio being fixed during the joint inversion, the V P /V S ratios at each station used in the CCP stacking are the same as we set in the initial model of the joint inversion, which are smoothed from previous Hκ-c stacking results in the crust 27 and correspond to the AK135 model in the upper mantle 63 . We constructed CCP stacked RF images with different Gaussian parameters (i.e., 1.0, 2.0, and 3.5) to distinguish a true interface, which should be consistent in different frequency bands. The Moho is visible as a strong positive Ps phase in all three CCP images with depths between~45 km and~70 km (Supplementary Fig. 6). Referring to the same velocity model at each station, we also constructed CCP images using multiples (i.e., PpPs phases) that reverberated between the free surface and interfaces in the crust and upper mantle. Given the fact that the PpPs phases exhibit weaker sensitivities to the V P /V S ratio than the Ps phases, the consistency in Moho depths of the CCP images constructed with the PpPs phases and the Ps phases confirms the reliability of the CCP images and indicates the suitability of the V P /V S ratio used in the inversion of the V S model ( Supplementary Fig. 7).
Besides the Moho, our CCP images constructed with Ps phases also reveal clear but intermittent intracrustal signals (Supplementary Fig. 6). The marked intracrustal N P P P P P P P P Pa a a a a a a a am m m m m m m e e e e e ea a a a a au u u u u u  t teau teau a a au u u u u u u u  at  a a a a a a at t t t t t  at at at  a a a a a a a a a a P P P P P P P P P P P P e e e e e e e e e t t t t t t t t t t a a a a a a a a a a a l l l l l l l l l l l l l P P P P P P P P P P P a a a a a a a a a a e e e e e  m m m m m ma a a a a a a a an n n n n n n n n n n n nt t t t t t t t t t t t p p p p p p p p p p p p n n n n n n n nt t t t t t t t t t t p p p p p p l l l l l l l an n n n n n n n n n n n n n n n n n n n n n n n m m m m ma a a a a a a a a a  Ma. An earlier detachment of the lithospheric root led to a weaker lithosphere in the Tianshan than in the surroundings. When the cratonic Indian slab directly impinged on the Tarim Craton, the contractional stress from the India-Asia convergence was effectively delivered to the Tianshan. Consequently, the lithosphere of the Tianshan underwent intense deformation, which was accommodated by inherited, structurally controlled brittle deformation in the shallow crust and by plastic deformation close to the Moho. IKA Issyk Kul Arc, NA Naryn Arc, STAC South Tianshan Accretionary Complex.

i i i i i i i i i it t t t t t t t t t t th h h h h h h h h ho o o o o o o o os s s s s s s sp p p p p eri p p p p p p p p p p p ph h h h h h h h h h h h h h h he e e e e e e e er p p p p p p p p p p p p p p p p p p ri er h h h h h h h h h h h h h he e e e e e c L L L L L L L L L L L L Li i i i i i i i i it t t t t t t t t t th h h h h h h h h h p p p p p o o o o e e e e i i i i le le e t t t t t t t t tl le e
signals are traceable in the image of individual RFs traces and in the CCP images with different stacked-bin parameters (2.5, 5.0, and 7.5 km), which largely excludes the possibility that the intracrustal signals are derived from noise or scattering ( Supplementary Fig. 8). Moreover, the CCP images produced by a fraction (20% and 60%) of reliable RFs are consistent with those made by the entire dataset, implying that the effect caused by an uneven data distribution or anisotropy structure is negligible ( Supplementary Fig. 9). We undertook some synthetic tests to check the reliability of the intracrustal Ps phases. We constructed the synthetic model at each station by overlaying the final 1-D model of V S and V P /V S ratio at depths shallower than 10 and 20 km, respectively, on a homogeneous half-space model ( Supplementary Fig. 10b, d). Then, we synthesized RFs with the same data coverage as the real dataset, and constructed synthetic CCP images by using the same processing parameters as used for the real dataset ( Supplementary Fig. 10c, e). The significant negative Ps phases at depths of 10-20 km beneath the Tarim Basin and Kazakh Shield are clearly generated by the effect of the sedimentary layer, but not by real velocity discontinuities. Except for these artifact signals in the northern and southern ends of the cross-section, synthetic CCP images keep their purities across the Tianshan. These tests strongly suggest that the marked intracrustal Ps phases (Fig. 2c, d; Supplementary Figs. [6][7][8][9] in the CCP images are robust. Thus, the discussion and interpretation in this study focused on these reliable features mostly in the Tianshan, but not at the ends of the cross-section. Balanced cross-section estimates of N-S shortening across the Tianshan. The V S model and CCP images collectively reveal a thrusting Moho geometry along the cross-section (Fig. 2c, d). Because the peneplanation in the Tianshan is generally inferred to have developed during the Mesozoic tectonic quiescence 44 , the Moho segments and thrusts imaged in this study can be used as a marker for estimating the amount of N-S shortening across the Tianshan in response to the India-Asia collision in the Cenozoic. Assuming the Moho beneath Tianshan was flat before the Cenozoic intracontinental deformation, we estimate that the N-S width before the deformation of the imaged region was~560 km by summing the lengths of each Moho segment, which indicates~85-km N-S shortening when considering thẽ 475-km width at present (Fig. 2e). It should be noted that a balanced cross-section of the Moho length should represent the minimum estimate because of the limited resolution of seismic images on some small-scale thrusts or doublets. Furthermore, we assume that the crustal volume remained constant during the Cenozoic deformation, and that the Tianshan had a ∼40-km thick crust before the Cenozoic, as indicated by the P-T equilibration conditions of xenoliths 49 and by the presentday Moho depth beneath the Tarim Craton and the Kazakh Shield (Fig. 2c, d). The present-day crustal area of~2.42 × 10 4 km 2 imaged along the 2-D cross-section implies a crustal dimension of~605-km N-S width and~40-km thickness before the Cenozoic deformation, and thus indicates~130-km N-S shortening across the Tianshan in the Cenozoic. Given the variation of previous crustal thicknesses and the erosion of the uplifted area, the N-S shortening estimated by balancing the crustal thickening also has some uncertainties. We consider these balanced crosssection estimates as a first-order indication of the N-S shortening of a hundred kilometers across the Tianshan. Consequently, considering the present-day N-S shortening rates of~12-15 mm/yr suggested by Late Quaternary fault slip rates and GPS observations across the estimated section (Fig. 2a, b), the N-S shortening of 85-130 km implies a start time of~11-6 Ma for the significant N-S shortening across the Tianshan.

Code availability
The Computer Programs in Seismology (CPS, version 3.30) software package is available from the Saint Louis University Earthquake Center (http://www.eas.slu.edu/eqc/eqccps. html). The common conversion points (CCP) stacking code is available from the author, Prof. Lupei Zhu, at his homepage http://www.eas.slu.edu/People/LZhu/home.html.