The Cenozoic Multiple-Stage Uplift of the Qiangtang Terrane, Tibetan Plateau

Cenozoic collision between the Indian and Asian continents is generally considered as the main driver forming the high Tibetan Plateau (TP). However, it remains hotly debated when and how the relatively flat and highly elevated TP was formed. Here, we present combined analyses of the apatite fission track (AFT) and apatite (U-Th)/He (AHe) of 18 granite samples along three steep topographic transects in the central part of the Qiangtang Terrane (QT), TP. The results indicate that the AFT ages of all samples are mostly between 130 Ma and 80 Ma, while the AHe ages range from 80 Ma to 40 Ma. Further thermal history modeling indicates that no significant cooling occurred after 40 Ma for most samples, except those lying close to the Reganpei Co extensional fault in the QT. The results are generally consistent with other low-temperature thermochronological data, as well as structural and sedimentologic data from the QT, suggesting that low relief and the relatively flat topography of the QT were almost completely formed before ∼40 Ma. As both megafossils and pollen had undergone a sharp change from subtropical- to psychro-species, indicating a relatively low elevation (∼2 km) at ∼40 Ma and >2 km uplift during the Oligocene. We use simple one-dimensional isostatic modeling to assess the contribution of convective removal of the lithospheric mantle to the present elevation of the QT. The results suggest that a combined effect of isostatic rebound (≥2 km) and thermal expansion related to asthenosphere upwelling and subsequent crustal base heating (∼0.4 km) led to the final uplift of the QT. Therefore, the QT experienced multiple-stage uplift processes which were controlled by crustal thickening before ∼40 Ma and lithospheric mantle delamination during the Oligocene, respectively.


INTRODUCTION
The Tibetan Plateau (TP) is a product of continent-continent collisions (Yin and Harrison, 2000). Until now, there has been a heated debate on growth models for the TP, which can be summarized into three end-member as "Proto-Tibetan Plateau Expansion" model (e.g., Wang C et al., 2008), "Valley" model (e.g., Ding et al., 2014) and "Two-Stage" model (e.g., Lu et al., 2018). Although the spatial-temporal patterns of uplift and the topographic development of the TP are largely controlled by crustal shortening and induced thickening patterns (e.g., Kapp et al., 2007), high elevation growth is clearly modulated by the negative buoyancy of denser lithospheric mantle because of isostasy (Hyndman and Currie, 2011;Cao and Paterson, 2016;Gvirtzman et al., 2016). Hence, the topographic evolution of the TP, the largest plateau on earth, is controlled by different factors and mechanisms, such as crustal shortening induced thickening and/or lithospheric mantle delamination caused by magmatic inflation and isostatic rebound (Molnar et al., 1993;Lee et al., 2015;Cao and Paterson, 2016;Chen et al., 2018). These mechanisms may act independently or in combination to produce punctuated episodes of uplift of the TP. In addition, the TP is an important site of weathering and erosion, and the spatial-temporal extent of its surface rise is accompanied by basin subsidence Li et al., 2012;Li et al., 2017), such as the eastern Himalayan syntaxis provides abundant sedimentary debris into the Bengal Fan carried by the Brahmaputra River (e.g., Gemignani et al., 2018). Thus, when the flat and high topography formed within the central TP, and its implications on plateau-building mechanisms, which are a combination of uplift and basin infilling, is still one of the most important focuses of geoscience research.
Based on structural geology, thermochronology, and sedimentology studies, the "Proto-Tibetan Plateau Expansion" model suggests that a significant portion of the central TP was >4 km in elevation prior to~40 Ma ( Figure 1A) (e.g., Murphy et al., 1997;Xu et al., 2013;Wang et al., 2014). This model suggests that the modern TP was probably formed by thrusting and related foreland basins that propagated from the central protoplateau outward both to the north and south during continued Indian plate indentation (Wang C et al., 2008). However, paleoelevation estimates by isotopic and paleobotanical methods in many locations have been used to argue that at least some regions were at low elevations in the central TP until the Miocene (e.g., Deng and Ding, 2015;Su et al., 2019;Spicer et al., 2020a). Based on these observations, they proposed the "Valley" model in which an E-W striking low-elevation valley (<2 km) is bound by high mountains (>4 km) to the north and south ( Figure 1B) (e.g., Ding et al., 2014;Wei et al., 2016;Su et al., 2019). The third "Two-Stage" model proposes that crustal thickening, induced by shortening, elevated the central TP to~2 km before~55 Ma, and~26 Ma lithospheric mantle delamination led to the crustal surface uplift to >5 km without obvious further crustal thickening ( Figure 1C) (Lu et al., 2018;Zhao et al., 2020;Zhao et al., 2021). Thus, there are three end-member models for the non-uniform spatialtemporal patterns of deformation and uplift in the central TP, which has led to debates about the timing and topographic evolution of the central TP (e.g., Wang C et al., 2008;Ding et al., 2014;Spicer et al., 2020a;Spicer et al., 2020b;Fang et al., 2020;Su et al., 2020;Spicer et al., 2021;Xiong et al., 2022;Zhang et al., 2022).
The growth of the modern TP has involved multi-stage collisions throughout the Mesozoic until the present, thus different pieces of the TP must possess distinct crustal thickening and relief generating processes that they have experienced different regional tectonic events (e.g., Spicer et al., 2020a). The Qiangtang Terrane (QT), occupying northern central Tibet, is an ideal place to investigate the topographic evolution of a key part of the TP and evaluate key components of above-mentioned models because it has exceptionally thick terrestrial sediments and abundant thermochronologic records since the Cretaceous (e.g., Zhao et al., 2017;Zhao et al., 2020). Moreover, because crustal thickness is internally coupled with its elevation by isostasy, therefore the internal relationships among crustal shortening/ thickening, lithospheric mantle delamination, and surface uplift can be used to reconstruct the topographic evolution of the QT. This paper publishes eighteen new apatite fission track (AFT) and apatite (U-Th)/He (AHe) ages, which are derived from four granite batholiths in the QT, and then reviews a large dataset of thermochronological ages to interpret the spatial-temporal pattern of crustal deformation and erosion across the width of the QT (Figure 2A). In addition, other well-documented geological records, such as paleoelevation and magmatic records, are used to evaluate QT elevation changes through the Cenozoic. By integrating the available geological and thermochronological data of the QT and applying simple one-dimensional isostatic modeling, we aim to address how the QT obtained its high elevation through the Cenozoic, especially after~26 Ma.
central QT, the Late Triassic Longmu Co-Shuanghu Suture subdivides the QT into the North and South QT . The Qiangtang Culmination is located in the South QT ( Figure 2C) (Kapp et al., 2003b), and exposes a Paleozoic metamorphic and sedimentary basement and Triassic subduction mélange, while the North and South QT are unconformably overlain by Late Triassic to Late Jurassic marine strata. The latter was again intruded by early-middle Cretaceous plutonic units (Kapp et al., 2003b). All these units are unconformably overlain by mostly flat-lying terrestrial Cretaceous to Tertiary volcano-sedimentary units with variable thicknesses in both subterranes (e.g., Sun et al., 2014;Li et al., 2017;Zhao et al., 2017).
The QT surface is currently relatively flat with elevations varying from 4,700 to 5,500 m ( Figure 2C) and is bounded to the north and south by positive topographic fronts with more than 500 m relief ( Figure 2C). The traces of the main thrust faults are roughly aligned with the east-west trending topographic features in both the South and North QT ( Figure 2B). These thrusts cause intense crustal shortening along foreland basins (e.g., Dai et al., 2012;Ma P et al., 2017;Dai et al., 2019;Zhao et al., 2020). The North QT is dominated by northvergent thrusts that control deposition and deformation in the Hoh Xil Basin ( Figure 2B) Dai et al., 2019). Additionally, southvergent faults in the South QT and the Gaize-Seling Co reverse faults in the northern Lhasa Terrane determine the deposition of a string of Cenozoic basins along the Bangong-Nujiang Suture, including the Gaize, Nima, and Lunpola Basins ( Figure 2B) Li et al., 2015). During Late Cenozoic times, an E-W tensile stress regime caused generalized extension in the QT and the development of wide NE-SWstriking intracontinental rifting basins ( Figure 2B) (e.g., Lu et al., 2018). Structural and Thermochronological data suggest that accelerated exhumation along the Qiangtang Culmination and positive topographic front occurred mostly during the Cretaceous-Paleogene, which may have had a profound impact on building of the modern TP (Dai et al., 2012;Zhao et al., 2017;Zhao et al., 2020).

Thermochronology Strategy
We collected 18 granite samples from three steep topographic transects, except the RGL-1 and RGL-2 samples, for integrated apatite fission track (AFT) and apatite (U-Th)/He (AHe) analyses. The granite exposures are late Triassic to Jurassic in age (Geological Publishing House (GPH), 2019). The samples, extending east-west for approximately 400 km, are located in the central portion of the QT ( Figure 2B). The easternmost six samples (from G-1 to G-13) range from 4,839 to 5,313 m in elevation, with an average elevation interval of~100 m (Supplementary Table S1). The central five samples (from R-1 to R-10) have a wider elevation range (from 4,922 to 5,802 m), thus encompassing a larger average elevation interval of 176 m (Supplementary Table S1). The westernmost five samples have a similar elevation (from 4,894 to 5,720 m) and elevation interval (~165 m) to those from the central region (Supplementary  Table S1).
Apatite grains were concentrated using standard heavy liquid and magnetic separation procedures. The AHe and AFT analyses were conducted using standard procedures described in detail in Li et al. (2016) and Gleadow et al. (2015), and also detailed and illustrated in the supplementary file. The number of single grains analyzed per sample is shown in Supplementary Table S2. Weighted mean ages of AHe are used to discuss the cooling age of individual samples. All obtained mean AHe and AFT ages are plotted in Figure 3A. To evaluate the quality of thermochronological data, we also plotted the individual AHe age of each analyzed grain against the eU value, which could indicate the effect of radiation damage on the AHe age ( Figure 3B) (Guenthner and Ketcham, 2013). To avoid the effect of large age ranges, which are generally related to a lower average effective U content and smaller average grain size, we tested all grains with eU × Radius and chose the values >100 for further data processing. Dating of small, low eU apatite is problematic due to potential implanted He (Reiners and Farley, 2001). Average AHe ages were calculated only for samples with grains yielding a relative standard deviation (RSD) < 50%. The calculation method of average age was followed by Zhao et al. (2017). The average AHe ages of all samples are presented in Supplementary Table S2. The approximate cooling rate was determined from the correlation of cooling ages with elevations ( Figure 3C).
AFT analysis results are shown in Supplementary Table S3. To better constrain the thermal history of the samples, the AFT and AHe ages were modeled in combination where possible using the inverse HeFTy approach of Ketcham et al. (2007), using Dpar as a kinetic parameter; AHe data were modeled using the radiation damage accumulation and annealing model (Flowers et al., 2009). We use a temperature of 10 ± 5°C for the mean present-day surface (e.g., Zhao et al., 2020). Where applicable, broad temperature constraints 60-120°C were set around the time range covering the AFT ages allowing more than ±2σ analytical uncertainties to set an initial time-temperature constraint box. C-axis projected track length data were applied. These premodeling settings were always applied with large uncertainties so as to give the inversion algorithm sufficient freedom to search for a wide range of possible thermal histories. Inverse thermal history modeling was run until >1,000 acceptable paths or >100 good paths (Supplementary Figures S1, S2). All thermal modeling results are summarized in Figure 4.
We then present thermochronological data of the large-scale age pattern across the QT by focusing on a cross-orogen transect extending from the Hoh Xil Basin to the northern Lhasa Terrane (Supplementary Table S4

One-Dimensional Isostatic Modeling
Crustal elevation is controlled by the thickness of both the crust and underlying lithospheric mantle. The available evidence shows Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 818079 that most high (>2 km) mountain belts are isostatically supported by a largely thickened crust (e.g., Lee et al., 2015). To gain a firstorder estimate of QT paleoelevation changes and evaluate summarized geological observations, we modeled the relationships among crustal and lithospheric mantle thickness and paleoelevation according to the equation of isostasy (Gvirtzman et al., 2016). The lithospheric mantle is denser than the surrounding asthenospheric mantle and thus has a "pull-down" effect on the elevation, depending on its thickness (L ml ) and density (ρ ml ) (Cao and Paterson, 2016). Equation 1 shows that the elevation was determined both by crust-and lithospheric mantle-induced elevation, which depends on their thickness and density as well as density in contrast to asthenospheric mantle density (ρ a ) (Eqs. 2, 3) (Gvirtzman et al., 2016). To calculate the absolute elevation, which is mountain height from sea level, we use H 0 = 2.4 km, where H 0 is a constant that allows using sea level as a reference datum instead of a theoretical free asthenospheric surface (Eq. 4) (Gvirtzman et al., 2016).
Here, H is the elevation, which indicates the distance from the mean mountain height to the theoretical free asthenosphere surface, and H c and H ml are the crust-and lithospheric mantle-induced elevations, respectively (Figures 7, 8), while H A is the elevation from sea level. L c and L ml are crustal and lithospheric thicknesses, respectively. ρ c , ρ ml , and ρ a are the characteristic densities of crust, lithospheric mantle, and asthenospheric mantle, respectively ( Figure 8). For a firstorder approximation, we assume that the crust and mantle have their own uniform densities. We also assume that the isostatic adjustment is instantaneous due to a shorter Maxwell time of the mantle (~1 kyr) (Ranalli, 1995) compared to the time scale of calculation (>1 Myr).

THERMOCHRONOLOGY RESULTS
Sixteen samples yielded valid weighted AHe ages ranging from 40 Ma to 134 Ma ( Figure 3A). The AHe ages of the sample R1, R3, and R7 are older than its corresponding AFT ages ( Figure 3A) that may due to the apatite remaining in the partial annealing zone for a long time as a result of slow exhumation (e.g., Galbraith and Laslett, 1993;Brandon, 2002). When we plotted AHe data against the eU value, the central samples (from R-1 to R-10) recorded positive relationships between AHe ages and the eU value ( Figure 3B). This most likely reflects radiation damage (Guenthner and Ketcham, 2013), resulting in erroneous cooling ages. After exclusion of these ages, the remaining AHe ages were between 80 and 40 Ma ( Figure 3A FIGURE 4 | A plot displaying the modeled temperature-time paths for all samples that yielded sufficient confined track data for modeling purposes. Modeling was performed using HeFTy (Ketcham et al., 2007). Where available, apatite (U-Th)/He data were incorporated into its respective thermal history model. Individual sample histograms and temperature-time plots are available in Supplementary Figures S1, S2.
Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 818079 lines in Figure 3C are best least-squares fitting lines of the cooling rates for samples considered in vertical profile according to Glotzbach et al. (2011). In total, 16 samples yielded valid AFT ages, approximately ranging from 80 Ma to 130 Ma ( Figure 3A). The modeling results show that eight samples experienced rapid (~2.5°C/Ma) cooling during the 60-35 Ma period, which was followed by extremely low cooling rates (<0.5°C/Ma) after~35 Ma ( Figure 4). Six samples show a rapid cooling rate of~4.5°C/Ma at approximately 30-20 Ma ( Figure 4). Assuming a paleogeothermal gradient of 39°C/km (Mechie et al., 2004), these eight samples experienced fast exhumation rates of 0.06 mm/yr and then a slow cooling rate of~0.01 mm/yr after 35 Ma. The R series samples underwent a rapid exhumation rate of 0.11 mm/yr between 30 and 20 Ma.

Data Combination Results
In constructional orogens such as that of the TP, most exhumation is due to erosion rather than normal faulting or ductile thinning. Thus, to the extent that erosion can be related to contractional deformation, cooling dates can be used to elucidate spatial-temporal patterns of deformation (e.g., Zhao et al., 2020). To recognize the variability in deformation and exhumation, combined low-temperature thermochronological ages from the QT and surrounding basins, including AFT, AHe, zircon fission track (ZFT), and zircon (U-Th)/He (ZHe) ages, were collected and plotted in Figure 5. To constrain exhumation and thus probability of deformation in the QT, we plotted dates and topography of the QT and neighboring basins as a function of distance from the location of significant topographic relief, represented by two N-S striking lines and one E-W striking line (Supplementary Figures S3, S5).
Taking all the thermochronological ages together, there is a broad correlation in the difference between the minimum and maximum ZFT, ZHe, AFT, and AHe ages with distance across the QT ( Figure 5). The fact that most ZFT and ZHe dates across the QT are older than ca. 80 Ma precludes Cenozoic exhumation deeper than approximately 6-8 km in most places on the QT, because of that thermal sensitivity ranges of ZFT and ZHe methods span from~175 to 250°C to~20-200°C, respectively (e.g., Ault et al., 2019). Few younger ZHe ages than ca. 80 Ma are found on the eastern QT, which belongs to an external drainage region ( Figure 5C) (Kapp and DeCelles, 2019). AFT ages mainly range between 120 and 40 Ma, while the Qiangtang culmination records >130 Ma cooling ages and the Kunlun Range shows 20-13 Ma cooling ages ( Figures 5A,B). The youngest range Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 818079 6 FIGURE 7 | Plot of best fit 1-D isostatic simulation. The density of the crust decreased from 2.89 g/cm 3 to 2.87 g/cm 3 due to thermal expansion after lithospheric mantle delamination. The red dashed line shows the possible maximum lithospheric mantle thickness, for more details, see the main text.

Relative Flat Topography Formed Beforẽ 40 Ma in the QT
The termination of shallow marine deposition and the onset of exhumation of the Qiangtang Culmination suggest that the accumulation of crustal thickening in the QT most likely began during late Early Cretaceous based on the 150-130 Ma ZHe cooling ages . Our new AFT data are mainly between 130 and 75 Ma, and the corresponding AHe ages range from 80 to 50 Ma, which indicates continued exhumation and erosion of the QT between 130 and 50 Ma ( Figure 3A). Thermochronology modeling results of the AFT data show that the first-stage fast cooling was between 60 and 35 Ma in the QT, which shows a 0.06 mm/yr exhumation rate that is slower thañ 0.36-0.41 mm/yr deduced by the AHe ages vs. elevation diagram ( Figures 3C, 4). This cooling event probably continued until~40 Ma, as supported by the documented AFT and AHe ages ( Figure 5). In other words, the patterns in our summarized AHe and AFT data suggest that the QT almost stopped obvious exhumation-and erosion-induced cooling after 40 Ma ( Figure 5).
Regionally, large parts of the northern Lhasa and Qiangtang terranes experienced rapid exhumation between~70 and 45 Ma, followed by a rapid decline in the exhumation rate (Hetzel et al., 2011;Rohrmann et al., 2012). The peneplain on the north Lhasa Terrane was formed by the establishment of low-relief occurred at low elevation by~45 Ma (Haider et al., 2013). Moreover, it is understood that the main topography of the QT was built before ca. 40 Ma (Figure 5), which resulted in less than 1 km of erosion across the QT after ca. 40 Ma (e.g., Rohrmann et al., 2012). Therefore, based on new and published thermochronological ages are older than~40 Ma across the QT, we infer that erosion across the QT has been less than 1 km and the QT may have had a relative flat topography since~40 Ma. Younger dates for the apatite systems are found farther to the north in the Kunlun Range, which possesses large relief and is adjacent to the Qaidam Basin ( Figures 5A,B). These late Cenozoic dates for the Kunlun Range require rapid and deep erosion of the hanging wall of a thrust fault from depths below the AHe partial retention zone at 20-5 Ma . AHe ages younger than~40 Ma also appear in the eastern QT, which indicates late Cenozoic exhumation within the external drainage region in the QT ( Figure 5C) (Kapp and DeCelles, 2019). These minimum AHe dates may not reflect the predominance of contractional deformation of the QT, but instead the cooling of older continental crystalline rocks by river incision-induced erosion (Dai et al., 2013). We also observed a fast cooling event between 30 and 20 Ma by thermomodeling of AFT data (Figure 4) sampled near the Reganpei Co extensional fault (Taylor and Peltzer, 2006). Because these samples are collected from granite not far from a Cenozoic NE-SW striking rift (Supplementary Figure S3) and the timing of this cooling event is consistent with the~N-S rifting developing in the TP (e.g., Styron et al., 2013), we attributed this cooling stage to represent faster erosion near the NE-SW rift basin that was induced by E-W extension. However, <40 Ma is not recorded in any AHe age, which may be because erosion depths do not exceed~1-2 km.

Evidences of pre~40 Ma Crustal Shortening and Thickening in the QT
Geophysical results indicate that the current crust beneath the QT is approximately 62-67 km thick (Gao et al., 2013). The earliest stage of crustal thickening in the QT occurred~166 Ma , and this led to the entire QT being uplifted above sea level by ca. 118 Ma (Kapp et al., 2003a). From then to the Eocene, >50% of N-S crustal shortening has been recorded in the Nima Basin and north Lhasa Terrane (Kapp et al., 2003a;Kapp et al., 2007;Volkmer et al., 2014). In the northern front of the Tanggula Range (North QT), a parallel thrust system developed during 71-41 Ma and absorbed~50% of the N-S shortening Li et al., 2017). In the South QT, crustal shortening and thickening also occurred mainly before~40 Ma (Zhao et al., 2020). South of the QT,~25 km shortening (~28%) of Nima Basin strata occurred between 30 and 23 Ma (Kapp et al., 2003a;Kapp et al., 2007), which could raise basin floor by certain degree. In the Hoh Xil Basin, thermochronological data and balanced cross section determined that 24% of upper crustal shortening occurred during the mid-Eocene (Staisch et al., 2016). Less than 10% of upper crust shortening occurred after~40 Ma in the interior of the South QT (Zhao et al., 2020). Inside the QT, strata older than~28 Ma are gently folded (Xu et al., 2013), and undeformed flat-lying early Miocene strata overlie these strata (Chen et al., 2018), which indicates that little deformation occurred after the Oligocene (Wu et al., 2012). According to these results, we conclude that crustal shortening within the QT mostly occurred before~40 Ma due to that the Lhasa Terrane subducted beneath the QT Zhao et al., 2020) during the Cretaceous and the India-Asian collision during early Cenozoic (e.g., Yuan et al., 2013;Jackson et al., 2018;Jian et al., 2018;Jackson et al., 2020;Li et al., 2020;He et al., 2021), and most of the <40 Ma crustal shortening mainly occurred within the Cenozoic basins, which are located in the northern and southern fronts of the north and south QT, respectively.
These structural study results are consistent with palaeomagnetic data and geochemical conclusions. Palaeomagnetic data derived from the Qiangtang and Lhasa Terranes revealed that they two experienced~15.3°latitudinal crustal shortening between 64-60 Ma and 54-43 Ma (Tong et al., 2017). However, N-S convergence within the central TP was considerably reduced at approximately 26 Ma, corresponding to a transition from compression to extension within the TP . The geochemical characteristics of granitic suites derived from the northern QT suggest that parental magmas were generated at depths of 40-50 km in the crust during the Late Cretaceous (Zhang et al., 2015;Lu et al., 2019). Combined with~47-34 Ma adakitic porphyry dikes, these results indicate that a large part of the QT had already attained near-maximum crustal thickness by the end of the Eocene (Wang Q et al., 2008;Ou et al., 2017;Wang et al., 2010). These volcanic rocks lie flat and unconformably on top of strongly shortened Mesozoic strata at modern elevations of 5,000 m (Law and Allen, 2020), which demonstrates that shortening ceased and that low-relief topographic conditions were achieved in this region by~40 Ma.

Terrestrial Basin Deposition Corresponds to Topographic Flattening During the Eocene Within and Neighbor the QT
During crustal thickening by thrusts before~40 Ma, the QT probably possessed large relief, providing detritus northward into the nonmarine Hoh Xil Basin and southward into a linear depocenter along the Bangong-Nujiang Suture, i.e., the lowland of the valley model, probably since the Late Cretaceous (e.g., DeCelles et al., 2007;Dai et al., 2013;Staisch et al., 2016;Kapp and DeCelles, 2019). The upper crustal continued exhumation and erosion, and the basin architecture in and around the QT are consistent with a decrease in topographic relief, perhaps reflecting the establishment of a relatively flat topography during the Eocene (Kapp and DeCelles, 2019).
To the north of the QT, terrestrial deposition was initialized at 71.9 Ma in the Hoh Xil Basin at the northern front of the Tanggula Range, and an abrupt sediment accumulation rate increase after 53.9 Ma was likely a response to the India-Asia collision (Jin et al., 2018). The thrust system in the Tanggula Range stopped and migrated into the Hoh Xil Basin at~47 Ma, which led to the Hoh Xil Basin becoming a source region for the younger terrestrial sandstones . Within the Hoh Xil Basin, shortening probably occurred no later than~34 Ma and largely ceased by~27 Ma (Staisch et al., 2016;Dai et al., 2019). In the early Cenozoic, the Hoh Xil Basin and Qaidam Basin evolved as a unified basin (McRivette et al., 2019), but after fault activity stopped in the Tanggula Range and within Hoh Xil Basin, deformation transformed to the Kunlun Range, which led to the separation of the Hoh Xil Basin from the Qaidam Basin by a >2 km relief (Chen et al., 2018).
To the south of the QT, the latest Cretaceous to Eocene basins along the Bangong-Nujiang Suture consist of narrow, topographically partitioned fluvio-lacustrine depocenters, i.e., a valley as in the valley model ( Figure 2B) Kapp et al., 2007;Ding et al., 2014). The Jurassic-Cretaceous rocks are unconformably overlain by up to 4 km of Palaeocene to Eocene fluvial red-bed deposits and~500 m Oligocene-Pliocene lacustrine sediments in the Lunpola Basin of the eastern Bangong-Nujiang Suture ( Figure 2B) (Wei et al., 2017). Bỹ 29 Ma, basin deposits of the Lunpola Basin reached the maximum burial depth and began to enter the uplift and denudation stage Xiong et al., 2022). This is consistent with the AHe ages being older than~40 Ma on both the northern and southern sides of the Lunpola Basin ( Figure 5B), which indicates little erosion near the Lunpola Basin and probably a lower relief after the Oligocene. In summary, strong tectonic activity, mainly thrusts, induced large relief and accompanied fast basin deposition in front of these thrusts before~40 Ma, coeval with Qiangtang and Lhasa Terranes crust thickening and flattening.

Paleoelevation Rise Before and After the Oligocene in and Around the QT
Paleoelevation of the central TP is a long-lasting and hotly debated topic (e.g., Deng and Ding, 2015;Spicer et al., 2020a;Spicer et al., 2020b;Su et al., 2020;Spicer et al., 2021). The paleoelevation of the QT and surrounding basins, including the Hoh Xi, Nangqian, Gaize, Nima, and Lunpola Basins, was close to sea level during the Cretaceous (e.g., Wu et al., 2012). The main controversy surrounds when the paleolevation of central TP reached its modern elevations after~60 Ma (Figure 1).
The "Proto-Tibetan Plateau Expansion" model proposed that the central TP, including the QT, was elevated to current elevations by~40 Ma (e.g., Wang C et al., 2008;Wang et al., 2014). >5 km paleoelevation during 51-28 Ma of the North QT was directly determined by isotopic methods, which is similar to the present elevation in the area (Xu et al., 2013). Nonetheless, 23 Ma time span of this paleoelevation estimate on the North QT is relative large, hence when the QT rise to >5,000 m has remained unknown until now. Moreover, >2 km thick red-bed strata (Xu et al., 2013) were difficult to deposit on top of the narrow QT highland which was bounded by >3 km significant relief on both north and south sides. In addition, the continental effect on the isotopic composition, which method was used in paleoelevation estimate in Xu et al. (2013), is complex in this inland setting (Spicer et al., 2020a). Thus, this paleoaltimetry need to be reappraised.
To the north of the QT, the Hoh Xil Basin attained its modern elevation by the mid-Eocene (Quade et al., 2011). Oxygen isotope-based estimates of the paleoelevation indicate that the Lunpola basin has been at an elevation of more than 4 km sincẽ 35 Ma (Rowley and Currie, 2006). However, these paleoaltimetry estimates derived from isotopic methods need be reappraised due to that complex isotopic fractionation effects in continental interiors and complex terrains (e.g., Spicer et al., 2020a). Recently, new estimated paleoelevation of the Bangor Basin, located South of the QT along the Bangong-Nujiang Suture, show the basin floor to have been at~1.5 km at 47 Ma (Su et al., 2020). This and other paleoaltimetric measurements (e.g., Su et al., 2019) indicate low elevation along the Bangong-Nujiang Suture before the Oligocene, which was also indicated by isotopic and paleoentological constrains in the Gaize Basin (Wei et al., 2016). The fossil finds from basins along the Bangong-Nujiang Suture indicate elevations below 2.5 km until sometime in the Oligocene (Spicer et al., 2020a;Spicer et al., 2020b;Spicer et al., 2021). This is supported by the conclusion that the elevation of the Lunpola Basin was <2.0 km at 50-38 Ma and reached >4 km by~29 Ma (Xiong et al., 2022). This is also consistent with the observation that relatively high paleoelevation (4.5-5 km) is obtained bỹ 26 Ma in the Nima Basin Huntington et al., 2015).
Relief between the QT and basins distributed along the Bangong-Nujiang Suture was probably not larger than 1 km since~40 Ma due to following reasons. Firstly, almost all AHe ages derived from the QT and surrounding basins and the northern Lhasa Terrane are older than~40 Ma, which means less than 1-2 km erosion have occurred after 40 Ma in these regions if we adopt 30°C/km geothermal gradient, due to that thermal sensitivity ranges of AHe method span from~40-60°C ( Figure 5) (e.g., Ault et al., 2019). Because of significant relief will inevitably lead to rapid denudation and can be recorded by low temperature thermochronological methods, therefore >40 Ma AHe ages could indicate relative low relief existed within and around the QT since then. Secondly, main crustal thickening of the QT and surrounding basins has terminated before 40 Ma according to tectonic and igneous evidences (Figure 6), and there was only~500 m thick sediments deposited in these basins along the Bangong-Nujiang Suture since the Early Oligocene (e.g., Wei et al., 2017;Xiong et al., 2022). Hence, relative flat topography probably already formed between the QT and these basins as nowadays since~35 Ma Li et al., 2022;Xue et al., 2022). Otherwise it's difficult to explain how large relief between the QT and these basins has been reduced after the early Oligocene.
In summary, during the Eocene, the relief between the QT and the surrounding basins gradually decreased to present-day topography, which was reached in the late Eocene by shedding sediments northward into the Hoh Xil Basin (e.g., Dai et al., 2012) and southward into basins along the Bangong-Nujiang Suture (e.g., Kapp et al., 2007). Later on, the paleoelevation of this low relief region gained >2 km rise in the central TP during the Oligocene ( Figure 6) Su et al., 2019;Spicer et al., 2020a;Spicer et al., 2020b;Xiong et al., 2022), which did not involve strong crustal thickening hence no intensive erosion. This is consistent with palaeoaltimetric studies suggesting that the Hoh-Xil Basin underwent 1,700-2,600 m of surface uplift during the late Eocene and early Miocene (Polissar et al., 2009;Quade et al., 2011). Evidence from igneous geochemistry and geophysics indicates that the QT may have experienced lithospheric mantle delamination and subsequently hot asthenospheric upwelling during~36-26 Ma (e.g., Lu et al., 2018). The~30 Ma K-rich mafic rocks and crust-derived trachytes in the central QT were attributed to lithospheric delamination and asthenoshperic upwelling (Ou et al., 2018;Ou et al., 2019). On the South QT, the widely distributed Nading Co volcanic rocks are~36-28 Ma and sit flat at relatively high modern elevations (generally~5 km) atop a basement of shortened Paleozoic-Mesozoic rocks . >25 Ma High-Mg alkaline basanite occurs within the Hoh Xil Basin (Yakovlev et al., 2019). This was attributed to underplating of mantle-derived magma and correlated with the sharp uplift of the QT and Hoh Xil Basins compared to the stable Qaidam Basin at~25 Ma (Chen et al., 2018). The mantle delamination model is also supported by the thin lithosphere depth of ca. 100 km in the north-central TP (Jiménez-Munt et al., 2008). Thus, the spatial correlation with ultrapotassic and adakitic magmatism supports the hypothesis of convective removal of thickened Qiangtang lithosphere causing major uplift of the QT during 37-26 Ma (Kapp and DeCelles, 2019).

ISOSTATIC MODELING RESULTS AND ITS IMPLICATIONS
Our new thermochronological data indicate that crustal shortening and thickening had already occurred beforẽ 40 Ma in the QT and relatively flat topography formed in and around the QT. Moreover, the QT probably has experienced fast paleoelevation rise by~26 Ma, similar with basins located to the south of the QT. In addition, lithospheric mantle had delaminated beneath the QT before~26 Ma. Because topographic growth is clearly controlled by crustal thickening and beneath lithospheric thinning because of isostasy (Cao and Paterson, 2016;Gvirtzman et al., 2016), lithospheric mantle delamination is an important driving force for fast QT uplift during~36-26 Ma. Hereafter, we use 1-D isostatic modeling, by which we aim to describe that the crust of the QT has not thickened but largely elevated by~26 Ma (e.g., Lu et al., 2018).

Model Constraints
From the above summary, there are three main time-duration constraints for our calculation ( Figure 6). The first age range is 120-40 Ma, during which both crust and lithospheric mantle experienced~50% shortening ( Figure 6). During 40-30 Ma, both crust and lithospheric mantle shortened by another~10%. Original crustal and lithospheric mantle thicknesses use global average numbers for sea level at 30 and 75 km, respectively (Figures 7, 8A) (Gvirtzman et al., 2016). This indicates that the Qiangtang crust was close to sea level at~120 Ma (Wu et al., 2015). The modern constraint is that the Qiangtang crust reaches >62 km thick and the surface elevation is close to~5 km (Figures 7, 8D) (e.g., Zhao et al., 2020).
We assume plane strain deformation without shortening or stretching parallel to the orogenic belt. The crustal thickness is thickened from 30 to 60 km from 120 Ma to 40 Ma and further thickened to 66 km by another~10% shortening during 40-30 Ma, but it is reduced to~62 km thick due to continued erosion and probably extensional rifts induced by lower crustal flow (e.g., Ryder et al., 2014). Erosion-induced thinning should be considered because of the exposure of the <40 Ma AHe sample and Late Cenozoic rift in the eastern QT (e.g., Dai et al., 2013).
The lithospheric mantle also thickened during crustal shortening with the same shortening ratio because there was no evidence showing any decoupling between crust and lithospheric mantle. Therefore, lithospheric mantle is also thickened from 75 km to >165 km thick during upper crustal shortening . Nonetheless, the lithospheric mantle delaminated during 36-26 Ma beneath the QT according to magmatic records ( Figure 6) (e.g., Guo and Wilson, 2019). Therefore, the lithospheric mantle probably has never reached >165 km thickness (e.g., Wu et al., 2021). We changed crustal density from 2.7 to 2.9 g/cm 3 and lithospheric mantle density from 3.28 to 3.3 g/cm 3 to find the best fit model according to other similar studies (e.g., Cao and Paterson, 2016;Lamb et al., 2020). The asthenospheric mantle density is constant at 3.25 g/ cm 3 (Gvirtzman et al., 2016;Lamb et al., 2020).

Modeling Results and Limitations
The outputs of our isostatic modeling are the histories of elevation, crustal thickness, and lithospheric thickness from 120 Ma to present, which can be tested by comparison with independent geological observations ( Figure 6). The best fit model could explain >2 km of fast crustal uplift at~26 Ma and other geological observations (Figures 6, 7). By our modeling results, we determined the following: 1) the crust and lithospheric mantle density should be~2.89 g/cm 3 and 3.3 g/cm 3 , respectively, which could account for the initial sea level height and QT elevation before~26 Ma (Figures 8B,C); 2) the elevation remains at~2.3 km when the crustal thickness is already over 65 km ( Figure 7); and 3) the lithospheric mantle delamination time determines elevation changes but not the final crustal elevation and topography. Here, we decrease crustal thickness from 66 km at 26 Ma to 61 km at present because QT crust thickness probably decreased as a result of extension, crustal flow, and/or erosional loss in externally drained parts of the orogeny ( Figures 8C,D) (e.g., Kapp and DeCelles, 2019). In addition, differences between delamination-induced magmatic records during~37-26 Ma and fast uplift at~26 Ma are probably due to isostatic rebound being delayed compared to mantle delamination-induced magmatic responses.
However, the best model requires a~166 km thick lithosphere, which is at odds with the fact that the maximum lithosphere thickness is generally <130 km (Rychert and Shearer, 2009;Lamb et al., 2020). When the lithospheric mantle reaches a critical thickness, it may start to locally delaminate, preventing it from becoming extremely thick (e.g., Lamb et al., 2020). This may have an effect on the paleaoelevation of the QT. Moreover, when we kept these modeling settings, the predicted current elevation was less than 4.9 km after lithosphere delamination, which conflicts with the >5 km high elevation in the QT (Zhao et al., 2020). One possible solution is thermal isostasy (Hyndman and Currie, 2011). Due to lithospheric mantle delamination, Qiangtang crust became hot, reaching 1,000°C at the Moho (Hacker et al., 2014), which may cause thermal expansion of the overlying crust and further contribute to crustal uplift (Hasterok and Chapman, 2007;Hyndman and Currie, 2011).

Thermal Expansion After Lithospheric Delamination Contributed to QT Crustal Rise
Elevation is related to the thermal regime, as defined by heat flow data, after correction for crustal thickness and density variations (Hasterok and Chapman, 2007). Within the QT, the temperature gradient is constrained by several documented studies. The temperature at a depth of 18 km is~700°C, where an α-β quartz transition occurs (Mechie et al., 2004). At the Moho, the temperature is~1,000°C, inferred from the~7.9-8.0 km/s Pn speed reported by Liang and Song (2006). Xenoliths provide an important additional thermal constraint. Mineral compositions derived from xenoliths within the Cenozoic basalt indicate equilibration temperatures of~900-1,000°C at 30-40 km depth prior to 28 Ma and subsequent heating to 1,100-1,200°C by 3 Ma beneath the central QT crust (Hacker et al., 2000;Ding et al., 2007). This is similar to the temperature gradient assumed for the QT (Craig et al., 2012).
If the QT underlain by a uniformly hot asthenosphere is accepted, we expect that thermal expansion will contribute a nearly constant amount to surface elevation (e.g., Hyndman and Currie, 2011). We assume that the~60 km depth Moho temperatures before and after lithospheric delamination arẽ 550 and 1,000°C, respectively. Using these summary temperatures and a coefficient of thermal expansion of 3.2 × 10 −5°C−1 (Hasterok and Chapman, 2007), if the original density of the QT crust is 2.89 g/cm 3 , the average density change is calculated to be~0.02 g/cm 3 . Therefore, the average crustal density decreased from 2.89 g/cm 3 to 2.87 g/cm 3 ; if we retained other parameters, the QT crust would raise another 0.4 km height, and the predicted current elevation would bẽ 5.3 km (Figures 7, 8D) at 66 km crustal thickness. There may be an elevation effect due to systematic mantle composition differences, but this would be much smaller than that for temperature (e.g., Kaban et al., 2003).
Although our modeling results cannot predict detailed crustal thickness or elevation variation over time, the most preliminary conclusions show that pure isostatic rebound by lithospheric delamination could not account for enough elevation change at 26 Ma. Hence, the~26 Ma sudden elevation change in and surrounding the QT is interpreted to be primarily due to not only isostatic rebound after lithospheric delamination but also the thermal isostasy effect of nearly constant high temperatures at Moho depths.

Proposed Multiple-Stage Elevation Evolution Model for the QT
Basins along the Bangong-Nujiang Suture, central TP, did not obtain modern elevation at~40 Ma (e.g., Su et al., 2019), which is inconsistent with the "Proto-Tibetan Plateau Expansion" model that proposed >4 km elevation of the central TP at 40 Ma (Wang C et al., 2008). In the "Valley" model, large relief between the Lunpola Basin and QT until~29 Ma was proposed (Su et al., 2019;Xiong et al., 2022), but other studies suggested that low relief topography already existed since at least~35 Ma Li et al., 2022;Xue et al., 2022). Thermochronological data ( Figure 5) indicate <1 km of relief between the QT and surrounding basins by~40 Ma. In addition, the sedimentary basins of the Bangong-Nujiang Suture, which represent the valley remnants, accumulated up to 2-3 km of Cenozoic sediments, most of which are primarily Eocene (Wei et al., 2017;Xiong et al., 2022). Therefore, significant relief did not exist anymore between the QT and the depocenter along the Bangong-Nujiang Suture since the late Eocene.
Here, we prefer the "Two-Stage" model proposed by Lu et al. (2018). In this model, the QT crust first shortened and thickened before~40 Ma, and shortening by thrusts induced large relief and fast erosion, which was recorded by thermochronological data (Figure 5). Relatively flat topography in and around the QT was formed during the termination of this stage, which was indicated by almost all AHe data older than~40 Ma within the Qiangtang internal drainage region ( Figure 5). Fast erosion results in 1) relatively flat topography within the QT and 2) coeval thick terrestrial basins developed around the QT with a relative high deposition rate. But, the elevation of the QT and surrounding basins were still lower than nowadays until 39-38 Ma (e.g., Su et al., 2019;Fang et al., 2020;Xiong et al., 2022), probably finally rise to >4,500 m by~29-26 Ma (DeCelles et al., 2007;Huntington et al., 2015;Xiong et al., 2022).
During the second uplift stage, only a small amount of crustal shortening occurred with and around the QT, which was supported by flat-lying terrestrial strata and lava Ou et al., 2019). However, elevation changed sharply to modern elevation during~29-26 Ma (e.g., Hoke et al., 2014;Jia et al., 2015) and climate change records in the Oligocene-Early Miocene on the northern Tibetan Plateau (Cheng et al., 2019). The only reasonable cause for this sudden uplift of the QT is lithospheric mantle delamination, which was well recorded by late Cenozoic magmatism within the central Tibetan Plateau (e.g., Guo and Wilson, 2019). Moreover, lithospheric mantle delamination has twofold effects on the QT crust uplift: isostatic rebound and thermal isostasy.

CONCLUSION
Our new AFT and AHe ages show that the QT experienced a slow cooling rate of~0.07 mm/yr during~130-70 Ma, followed by a fast cooling rate of~0.36-0.41 mm/yr during~70-50 Ma. There are no cooling event records younger than 40 Ma in our new thermochronological ages. Moreover, low-temperature thermochronological dates along three transects across the QT display narrow trends of minimum AHe cooling ages at~40 Ma, indicating low relief and rapidly decreasing erosion rates after 40 Ma. This means that relatively flat topography was already formed in and surrounding the QT. Moreover, crustal thickening of the QT finished before~40 Ma, and in many places on the QT, erosion did not exceed 1-2 km after~40 Ma, although these regions are now 4.7-6 km above sea level.
Summarized geological data indicate that the lithospheric mantle delaminated beneath the QT during~37-26 Ma and accompanied asthenospheric mantle had upwelled, which accounted for >2.5 km of the sudden crustal uplift during the Oligocene. By using 1-D simple isotactic modeling, we primarily concluded that 1) the~2.3 km elevation of the QT was induced by crustal thickening during~120-40 Ma, which was associated with enhanced intermontane basin subsidence and flat topography evolution (stage 1); and 2) around~26 Ma, the elevation of the QT and surrounding basins increased from~2.3 km to >5 km without obvious crustal thickening involvement but instead was due to lithospheric mantle delamination, which led to not onlỹ 2 km isostatic rebound crustal uplift but also~0.4 km thermal isostatic rising (stage 2).
The timing and dynamics of the Tibetan Plateau surface rise have been proven difficult to reconcile with the independent crustal thickening process, especially where has experienced complex mantle processes, such as lithospheric mantle delamination. Our model certainly has some limitations due to that our explanation was mainly based on current existing data which are still too less to build a more sophisticated plateau evolution model. In the future, more detailed crustal thickening process, lithospheric mantle evolution, crustal density change and precise paleoaltimetry estimate during the Cenozoic across the Tibetan Plateau are needed for cross-validation via isostasy.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.