Zhong-bao Zhao, Chao Li, Xu-xuan Ma

Institute of Geology, Chinese Academy of Geological Sciences, Beijing 100037, China

Keywords:Crustal thickness Paleoelevation Isostatic modeling Qiangtang terrane Lhasa terrane Tibetan Plateau China

ABSTRACT When and how the Tibetan Plateau formed and maintained its thick crust and high elevation on Earth is continuing debated. Specifically, the coupling relationship between crustal thickening and corresponding paleoelevation changing has not been well studied. The dominant factors in crustal thickness changing are crustal shortening, magmatic input and surface erosion rates. Crustal thickness change and corresponding paleoelevation variation with time were further linked by an isostatic equation in this study. Since 120 Ma crustal shortening, magmatic input and surface erosion rates data from the central Tibetan Plateau are took as input parameters. By using a one-dimensional isostasy model, the authors captured the first-order relationship between crustal thickening and historical elevation responses over the central Tibetan Plateau,including the Qiangtang and Lhasa terranes. Based on the modeling results, the authors primarily concluded that the Qiangtang terrane crust gradually thickened to ca. 63 km at ca. 40 Ma, mainly due to tectonic shortening and minor magmatic input combined with a slow erosion rate. However, the Lhasa terrane crust thickened by a combination of tectonic shortening, extensive magmatic input and probably Indian plate underthrusting, which thickened the Lhasa crust over 75 km since 25 Ma. Moreover, a longstanding elevation ＞4000 m was strongly coupled with a thickened crust since about 35 Ma in the central Tibetan Plateau.©2021 China Geology Editorial Office.

1. Introduction

Since the onset of the India-Asia continental collision at ca. 60 Ma (Hu XM et al., 2016), a continued India-Asian convergence has produced the highly elevated Tibetan Plateau(TP) (Harrison TM et al., 1992; Tapponnier P et al., 2001).These hypotheses include the following, potentially coexisting, scenarios: (1) The distributed shortening and thickening of the Tibetan Plateau crust as a “viscous sheet”(England P and Houseman G, 1989); (2) extensive thrusts thickened the Tibetan Plateau crust tectonically (Harrison TM et al., 1992), and (3) successive intra-continental subduction and sedimentary infilling of the intermontane basins(Tapponnier P et al., 2001). However, recent research results argued that the central Tibetan crust had already thickened and uplifted before the India-Asian collision (Kapp P et al.,2007; Murphy M et al., 1997). And, a high-standing Early Cenozoic paleoelevation existed in the central Tibetan Plateau(Ding L et al., 2014; Rowley DB et al., 2006; Wang CS et al.,2008, 2014), which suggested that a proto-plateau formed before ca. 45 Ma, including Qiangtang and Lhasa terranes where belongs to the core of the central Tibetan Plateau(Wang CS et al., 2008). As such, how and when the thick continental crust and high elevation formed in the central Tibetan Plateau remains an unanswered question, more specifically, the interaction between changes in crustal thickness and surface expressions during plateau formation are still vague.

Factors controlling crustal thickness include the interplay between tectonic shortening, magmatic addition and surface erosion (Cao WR et al., 2016; Lee C-TA et al., 2015). Both tectonic shortening and the influx of mantle-derived magma can thicken the crust and trigger surface uplift, which corresponds to faster erosional processes that thinning the crust (Jiang HH and Lee C-TA, 2017). In addition, crustal rooting or lithospheric mantle foundering can further raise the crust elevation (Cao WR et al., 2016; Chen M et al., 2017).Although there are an abundant number of studies on the Tibetan Plateau formation (Guillot S et al., 2013; Li YL et al.,2015; Wang CS et al., 2014), the causes of plateau crustal thickening and corresponding topographic evolution remain unknown to us. Most studies invoke the role of crustal shortening and thus high elevation that induces thickening(Kapp P et al., 2005), and others have invoked Indian plate underthrusting underneath the Tibetan Plateau (DeCelles PG et al., 2002), Asian continental subduction (Replumaz A et al.,2016) or lithospheric mantle remove (Chen M et al., 2017).Nonetheless, most of these hypotheses have not yet been quantitatively tested.

Here, this work aims to integrate crustal thickening and reasonable surface uplift processes together, which could offer more constraints on how and when the central Tibetan Plateau attained its thick crust and high elevation. The authors first collected and summarized basic geological observations,including the crustal shortening rate, magmatic input flux and erosion rates, which dominated crustal thickness and paleoelevation changes since 120 Ma for the Qiangtang and Lhasa terranes (Fig. 1). A simple one-dimensional mass balance isostatic model was applied to link all of these processes together. Primary results showed that crustal thickening and elevation variation were strongly coupled and comparable with documented geological observations.

2. Methods

Due to crustal thickness is determined by tectonic shortening/extension, magmatic addition and surface erosion process (Cao WR et al., 2016; Lee C-TA et al., 2015) and could be estimated by empirical geochemical methods(Mantle GW and Collins WJ, 2008; Profeta L et al., 2015),therefore the authors firstly collected these data from documented literature (Fig. 2, and supplementary data file).The tectonic shortening/extension and surface erosion rates could be directly obtained from literature, but magma input flux and historical crust thickness have to be transformed from other parameters. Thus, here the authors briefly introduced that how the authors calculated magmatic addition rates and estimated historical crustal thickness in the following section.

2.1. Magmatic addition rate calculations

The authors choose longitude between 85°E and 90°E for the Qiangtang and Lhasa terranes to estimate magmatic addition rates (MARs) separately, which are relatively well studied in the field of magmatism. The MAR can be estimated from areas of different intrusive suites during special temporal intervals within the target area (detailed data in supplementary Table S1). The zircon U-Pb ages of plutons and volcanoes from the literature are contoured at intervals of the geological time scale on a geologic map (Fig. 3). Areas with plutons and volcanoes within each geological time interval measured from the map by ImageJ are plotted against the ages in Fig. 4. The total volume of magma can be calculated from the areal distribution of the plutons multiplied by the thickness of the plutons. The thickness is calculated by an empirical line,which plots the pluton distribution area against the pluton thickness with logarithmic coordinates (Brown M, 2013).Final magmatic thickening rates or magmatic flux rates(km3/km2/Ma) are then obtained by dividing the volumetric addition rate by the total area (Jiang HH and Lee C-TA,2017). Therefore, these calculated magmatic flux rates are averages of magmatic fluxes for chosen area Fig. 4.

Fig. 1. Elevation map of the Tibetan Plateau with main tectonic units and their boundaries. Elevation data were downloaded from http://www2.jpl.nasa.gov/srtm/. Locations of the collected data are shown, and more information is in their supplementary Tables S1-S5.JSS-Jinsha suture; LSS-Longmu Co-Shuanghu suture; BNS-Bangong-Nujiang suture; IYS-Indus-Yarlung suture.

Fig. 2. Collected and reprocessed data of geological observations are shown. The best fit modeled results are illustrated by curves. a-timing of two major episodes during tectonic deformation for the Lhasa and Qiangtang terranes (see supplementary Table S3). b, c-magmatic input flux rates estimates for the Lhasa and Qiangtang terranes, respectively (Fig. 4; supplementary Table S1). Blue and red lines are obtained from modeling results which are used to compare with calculated line. Green line is cited from Ma XX et al. (2021) which is statics of documented magma ages on the Lhasa terrane. d-geological constraints for crustal thickness in the Lhasa and Qiangtang terranes (see Zhu DC et al., 2017 and supplementary Table S2). e-geological constraints for paleoelevation. f-in this figure, red and blue crosses show the collected results for the Lhasa and Qiangtang terranes, respectively. The blue and red curves show the results of the reference simulation for the Qiangtang and Lhasa terranes,respectively. Crosses show the geological constraints for erosion rates.

2.2. Crustal thickness estimate

The basic principle for estimating historical crustal thickness is that the Ce/Y ratio of arc basalts (Mantle GW and Collins WJ, 2008) and the La/Yb and Sr/Y ratios of arc intermediate rocks (Profeta L et al., 2015) intrinsically reflect the presence of mineral assemblages (amphibole + plagioclase± garnet) in the magma source region when slab-derived and fractionated samples are excluded. These ratios could indicate crustal thickness changes (Chapman JB et al., 2015; Profeta L et al., 2015). Methods for translating (La/Yb)n[where n denotes the normalized chondritic values (Mcdonough WF and Sun SS, 1995)] to crustal thickness are based on the method of Profeta F et al. (2015), which selects the intermediate intrusive rocks first (55%-68% SiO2). The empirical fit defined by the (La/Yb)nratios of global intermediate rocks with a crustal thickness (Profeta F et al.,2015) is used to track the crustal thickness of the Gangdese arc in the southern Tibetan Plateau (Zhu DC et al. 2017) (Fig.2d, red cross). For the southern Qiangtang terrane, the authors follow the methods of Profeta F et al. (2015) and Zhu DC et al. (2017), which filtered the geochemical data first to eliminate the effect of differentiation on primary magma composition. The detailed results are available in supplementary Table S2, Figs. 2d, 5 (blue cross).

3. Geological observations and constraints

Fig. 3. Magmatism contoured by geological time, scale from 85°E to 90°E in the Qiangtang and Lhasa terranes. Abbreviations are same as Fig. 1.

Fig. 4. Magmatic input flux estimates for the Lhasa and Qiangtang terranes, respectively (data are shown in the supplementary Table S1). Gray histograms show the pluton-exposed area, and the black squares show magmatic influx rates.

Fig. 5. The crustal thickness variation of the southern Qiangtang terrane from 170 Ma to 55 Ma. The crustal thickness variation of the southern Qiangtang terrane was based on data from supplementary Table S2. Different colors means collected from different literatures.

The Tibetan Plateau (Fig. 1) is a typical continentcontinent collisional region that was built on the Mesozoic-Cenozoic convergence between the Indian and Eurasian plates(Yin A and Harrison TM, 2000; Meng J et al., 2020). The Qiangtang terrane collided with Eurasia during the Late Triassic, while the Lhasa-Qiangtang collision occurred along the Bangong-Nujiang suture during the Early Cretaceous,which formed basic framework prior to the modern central Tibetan Plateau (Fig. 1; Bian WW et al., 2017; Zhao ZB et al.,2017). This is consistent with the paleomagnetic data, which indicates that the Lhasa-Qiangtang collision occurred at ca.120 Ma (Li C et al., 2019; Meng J et al., 2018). The Indian continent finally collided with Eurasia at ca. 60 Ma along the Indus-Yarlung suture leading to the Tibetan Plateau finally raised (Hu XM et al., 2016; Wang CS et al., 2014). Hence, the Tibetan Plateau formations resulted from a combination of several geological processes and were inherited from an early geological basis.

For building up logical model the authors first summarized geological observations from the central Tibetan Plateau since 120 Ma, which determined the historical crustal thickness and relative paleoelevation of the Qiangtang and Lhasa terranes (Figs. 1, 2, 5). Detailed methods of the data collection and calculations are included in the methods section and summarized results are plotted in the Fig. 2.

Both crustal shortening and mantle-derived magmatic additions are important processes that induce crustal thickening (Jiang HH and Lee C-TA, 2017). A large amount of crustal shortening and thickening occurred during 120-50 Ma, with an average shortening rate of 50% and 55% for the Qiangtang and Lhasa terranes, respectively (Fig. 2a,supplementary Table S3) (Murphy M et al., 1997; Kapp P et al., 2005, 2007; Volkmer JE et al., 2007, 2014; Zhao ZB et al., 2020). While, after 50 Ma there are only ca. 10% crustal shortening occurred within the south Qiangtang terrane (Zhao ZB et al., 2020) and very weak crustal deformation on the Lhasa terrane (Fig. 2a; Table S3) (Volkmer JE et al., 2014).These geological observations are in agreement with paleomagnetic results (Meng J et al., 2017).

Based on the simple statistic results, there was one strong magmatic addition stage with a magmatic flux rate of ca.0.01 km3/km2/Ma for the Qiangtang terrane during the Cenozoic Era (Figs. 2b, 4; supplementary Table S1).Nonetheless, three high-magmatic addition events have occurred over the Lhasa terrane since 120 Ma (Figs. 2c, 4;supplementary Table S1). The largest magmatic input flux reached ca. 0.034 km3/km2/Ma during the Paleocene on the southern Lhasa terrane (Figs. 2c, 4). However, the other two peaks flux numbers were relative ). Due to erosion and a lack of exposures, the estimated magmatic flux was at its minimum. The authors also compared the calculated magmatic flux rate curve with magmatic age probability distribution curve from ompared the calculated magmatic flux rate curve with magmatic age probability distribution curve from Ma XX et al. (2021) for the Lhasa terrane. Although peaks of these two lines are shifted a bit, the basic trending fit quite well (Figs. 2b, c).

These magmatic and deformational events changed the crustal thickness and further affected the elevation and surface erosion of the central Tibetan Plateau (Figs. 2d, f). The erosion rate decreased from 0.25 km/Ma at ca. 110 Ma to a very slow rate of ＜0.1 km/Ma at ca. 75 Ma on the southern Qiangtang terrane (Fig. 2f; supplementary Table S4; Lu L et al., 2015; Zhao ZB et al., 2017). While the erosion rate decreased from 0.3 km/Ma at 65 Ma to 0.05 km/Ma at 45 Ma on th;hasa terrane (Fig. 2f; Hetzel R et al., 2011; Rohrmann A et al., 2012), during the Cenozoic, the erosion rate was obviously small in both the Lhasa and Qiangtang terranes(＜0.1 km/Ma) because the central Tibetan Plateau became plateau-like by ca. 45 Ma (Gourbet L et al., 2016; Rohrmann A et al., 2012).

Using above mentioned inputs, the authors can calculate crust thickness changing with time and further compare them with crustal thickness estimated by geochemical methods(Fig. 2d). For the Qiangtang terrane, the authors collected and filtered intermediate rock data following the methods of Zhu DC et al. (2017) to calculate historical crustal thickness.Detailed results are in Figs. 2d, 5 (blue crosses) and supplementary Table S2. The Qiangtang terrane crust increased from 48 km to ＞60 km between 120 Ma and 60 Ma(Fig. 2d). Afterwards, the Qiangtang crustal thickness was nearly maintained at ca. 60 km. Crustal thickness of the south Lhasa terrane was obtained by the (La/Yb)nmethods (Fig. 2d,red crosses) (Hu FY et al., 2020; Zhu DC et al., 2017). The Lhasa crust thickened from ca. 37 km at 120 Ma to ca. 50 km at 70 Ma, then thickened sharply from ca. 50 km to ＞60 km during 70-50 Ma;Fig. 2d;Hu FY et al., 2020; Zhu DC et al.,2017). However, the northern and central Lhasa terrane crustal thickness was suddenly dropped to ca. 40 km and thus elevation to ca. 2.5 km during 50-70 Ma in Hu FY et al.(2020) by delamination Geological which is a bit wired with that delamination normally lead to crust uplifting (Jiang HH and Lee C-TA, 2017). Finally, crustal thickness of the Lhasa terrane increased from ＞60 km to ＞75 km after 50 Ma (Zhu DC et al., 2017).

Geological records suggested that the plateau elevation and erosion rate fluctuated during the plateau evolution (Figs.2e, f; supplementary Tables S4, S5). Based on the paleoelevations of the Qiangtang and Lhasa terranes (Fig. 2e),the authors found that both terranes stood ＞4 km high since ca. 50 Ma (Ding L et al., 2014; Liu X et al., 2016; Xu Q et al.,2013). The result has benefited from the development of a number of paleoaltimetry analyses and their application to studies of the TP (Ding L et al., 2014; Polissar PJ et al., 2009;Liu X et al., 2016). The original ca. 2 km elevation of the Qiangtang terrane at ca. 120 Ma was documented by Zhang KJ et al. (2014) and Hu FY et al. (2020). A special exception is the Lunpola Basin, where located along the Bangong-Nujiang Suture (Fig. 1) and recorded ＜2.3 km elevation until 40 Ma (Fang XM et al., 2020).

In total, two slab foundering events were predicted by documented data and adopted into the model. The first involved the Neotethyan oceanic lithosphere that foundered before 45 Ma, and the second consisted of the Indian lithosphere that delaminated at ca. 30 Ma (Chen M et al.,2017; Razi AS et al., 2016). After this, the long-term low velocity of convergence rate supports a continuous hard collision between India and Asia (van Hinsbergen DJJ et al.,2011). Hence, the final elevation was also controlled by lithospheric/crustal root foundering at ca. 45 Ma and 30 Ma,respectively. This process deeply affected paleoelevation of the central TP. Asthenospheric upwelling associated with break-off events may explain the patterns of Cenozoic volcanism in the TP (DeCelles PG et al., 2002).

4. Isostatic modeling

Aiming to combine all of these observations together and attempting to find their internal relations, the authors applied the isostatic modeling method of Cao WR and Paterson S(2016). Here, the model used deformational, magmatic and erosion rates parameterized from observations as inputs. The outputs were historical crustal thickness and elevation values from 120 Ma to present, which were tested by a comparison with independent geological observations (Fig. 2). Lithospheric/crustal root foundering was taken into account as well.

To calculate the elevation as a function of time [h(t)], the authors used equation (8) and equation (9) of Cao WR and Paterson S (2016) and their Matlab scripts to solve the following equations:

4.1. Isostatic modeling assumption

For a first-order approximation, the authors assume that the crust and mantle have uniform densities, respectively. The authors also assume that the isostatic adjustment is instantaneous due to a shorter Maxwell time of the mantle (ca.1 ka) (Burov EB, 2011) compared to the time scale for the simulation (＞10 Ma). A plane strain deformation without shortening or stretching parallel to the orogenic belt is assumed. The erosion rate is proportional to the elevation divided by the erosion response time (τE), which becomes a function of bedrock erodibility and precipitation rate[equation 6 and equation 7 in Cao WR and Paterson S(2016)]. The continental root, garnet-rich pyroxenites or eclogites in the lower crust or upper mantle are denser than the surrounding mantle, which creates a “pull-down” effect on the elevation depending on its thickness (R) and density (pr)(Saleeby J, 2003). All rate curves follow a Gaussian distribution, and the area beneath the Gaussian curve is scaled to the total amount of variation Cao WR and Paterson S(2016). The Matlab code is available in Cao WR and Paterson S (2016).

Because the Lhasa and Qiangtang terranes are separated in the central Tibetan Plateau and probably possess historical disparate uplifting (Hu FY et al., 2020), the authors modeled the two terranes separately. The original crustal thickness of the Qiangtang and Lhasa terranes were 48 km (supplementary Table S2) (Hu FY et al., 2020) and 37 km (Zhu DC et al.,2017), respectively. The initial elevation of the Qiangtang terrane was approximately 2 km (Hu FY et al., 2020);however, the Lhasa terrane was close to sea level (Kapp P et al., 2005; Hu FY et al., 2020). The authors assigned three separate evolution stages: 120-65 Ma, 65-45 Ma and 45-0 Ma. Deformational and magmatic input parameters are shown in Table 1, which derived from the collected data. The erosion rate, which is a function of elevation, was tested for several ranges to fit the observed data. To reduce the crustal root foundering effect, the authors used a root growing related number ofγ=1 (Cao WR and Paterson S, 2016). More detailed parameters for the model setup are shown in Table 1. The representative modeling results are shown in Figs. 6, 7.

4.2. Modeling results compared with geological observations

The best fit models are identical with the geologicalobservations and reveal several important crustal thickening processes for the Lhasa and Qiangtang terranes (Figs. 2, 6, 7).For a range of reasonable input conditions, the authors found that observed tectonic shortening alone can achieve modern crustal thickness for the Qiangtang terrane (Fig. 6a-1). Hence,tectonic shortening is the dominant factor in crustal thickening from the Cretaceous to the Paleocene to form＞60 km thick crust for the Qiangtang terrane (Figs. 2d, e).Modern Qiangtang crustal thickness is ca. 62 km (Gao R et al., 2013). The crustal thickening curve also coincides well with the (La/Yb)nresults (Fig. 2d).

Table 1. Parameters used in the simulation experiments.

However, the elevation of the Qiangtang terrane is approximately 4 km high in the Set 1 model, which does not fit the modern elevation (ca. 5 km) (Xu Q et al., 2013; Fig. 6a-3). One possible cause for this is that due to the asthenosphere upwelling beneath the Qiangtang terrane probably decrease crust density by thermal expansion (Hyndman RD and Currie CA, 2011). In this model, the authors decrease the crustal density from 2.8 kg/m3to 2.6 kg/m3and maintain the other parameters, which raises the modeled Qiangtang elevation close to the modern elevation (Figs. 6b-3, 2e). Early crustal shortening enhances the erosion rate, but the erosion rate decreases during the late crustal shortening stage, which is probably due to intermountain sedimentation (Fig. 2f;Tapponnier P et al., 2001). In this model, the authors reduce the erosion rate and set it to be similar to the documented numbers (Fig. 6b-4).

According to geological observations (Fig. 2a), the simulations suggest that the Lhasa terrane increased crustal thickness together with the Qiangtang terrane due to tectonic shortening (Fig. 7a-1) from 120 Ma to 65 Ma. However, the second significant phase of crustal thickening and elevation increased due to magmatic input from 65 Ma to 45 Ma (Fig. 2c),during which the crustal thickness increased from ca. 53 km to ca. 65 km and the elevation increased from ca. 2 km to ＞4 km during ca. 20 Ma (Figs. 7a-1, a-3). Here, because our observed magmatic input rates were at a minimum, the authors used a larger magmatic flux rate than the observed magmatic flux rate to obtain a reasonable crustal thickness(Fig. 2c; Table 1). During the last stage, the Set 3 simulation predicted that the thickness of the Lhasa crust was ca. 65 km(Fig. 7a-2) which is much thinner than it’s real thickness(Fig. 2d; Gao R et al., 2013). For achieving a ＞75 km crustal thickness, the model need additional ca. 10 km thickness(Fig. 7b-1). The authors prefer that this additional crust thickening was caused by the Indian plate underthrust beneath the Lhasa terrane (Gao R et al., 2017). In this model, the authors add additional crustal thickness by allowing a specific amount of magma input for achieving the Indian crust underthrusting (Figs. 7b-1, b-4).

The ca. 4 km elevation at the Paleocene from the modeling results was similar with the published elevation estimate (Fig. 2e; Ding L et al., 2014). Paleoelevation estimated for the Eocene to Pliocene Oiyug Basin, Penbo Basin, Eocene Nianbo Formation, all on Lhasa terrane, is ＞4.1 km by at least ca. 30 Ma (Ingalls M et al., 2018). Finally, the modeled elevation rose to ＞5.5 km, but induced a relative faster exhumation rate than the documented data (＜0.1 km/Ma)(Fig. 2f).

5. Crustal thickening and elevation rising in the central Tibetan Plateau

Crustal thickening and a rise in topography are mainly due to tectonic shortening, magmatic addition, large scale underthrusting, crustal root delamination or a combination of these processes (Cao WR and Paterson S, 2016; Chen M et al., 2017; DeCelles PG et al., 2002; Jiang HH and Lee C-TA,2017).

5.1. Crustal thickening and elevation rising in the Qiangtang terrane

Isostatic modeling results show that the Qiangtang terrane obtained ＞60 km thick crust and probably ＞4 km elevation during 120-50 Ma (Figs. 2d, e). That was because of the crust of the Qiangtang terrane tectonically shorten ca. 50% during 120-50 Ma which was induced by the Lhasa-Qiangtang collision (Kapp P et al., 2007; Zhao ZB et al., 2020). The Qiangtang terrane crust continuously shortened and elevated at a small rate due to the successive convergence of Lhasa-Qiangtang and a minor magmatic input after 80 Ma (Figs. 2d,2e; Zhao ZB et al., 2017, 2020).

Fig. 6. Results of Set 1-2 simulations. a-1 and b-1 represent crustal thickness; a-2 and b-2 represent total thickness, here, total thickness is equal to crustal thickness plus crustal root thickness; a-3 and b-3 represent simulated elevation; and a-4 and b-4 represent the magmatic thickening rates and erosion rates. The Set 1 simulation indicates the lowest and final paleoelevation of the Qiangtang terrane, with a crustal density equal to 2.8×103 kg/m3. Set 2 is the best fit simulation for the Qiangtang terrane when the crustal density is decreasing to 2.6×103 kg/m3.

East-west calc-alkaline and potassium-rich volcanic rocks developed from 45 Ma to 28 Ma in the Qiangtang terrane,which were related to lithospheric thinning (Ding L et al.,2007; Chung SL et al., 2005), which might be accompanied by a large amount of mafic intrusion that will lead to thermal expansion of the Qiangtang crust, and further increasing elevation (Fig. 6b; Hyndman RD and Currie CA, 2011).Actually, the negative P-wave velocity beneath the Qiangtang terrane was usually interpreted as a thin and delaminated lithospheric mantle that was replaced by a hot asthenosphere since 35 Ma, which lead to the final uplift of the central Tibetan Plateau (Chen M et al., 2017; Molnar P et al, 1993).

Fig. 7. Results of Set 3-4 simulations. a-1 and b-1 represent crustal thickness; a-2 and b-2 represent total thickness, here, total thickness is equal to crustal thickness plus crustal root thickness; a-3 and b-3 represent simulated elevation; and a-4 and b-4 represent the magmatic thickening rates and erosion rates. The Set 3 model does not consider the underthrusting of the Indian Plate beneath the Lhasa terrane. Set 4 is the best fit simulation for the Lhasa terrane after adding an additional ca. 10 km crustal thickness.

5.2. Crustal thickening and elevation rising in the Lhasa terrane

In this model, the Lhasa terrane crust was thickened during 120-65 Ma, 65-45 Ma and 45-0 Ma, thus corresponding elevation raised from sea level to ca. 2 km at 65 Ma, to ca. 4 km at 45 Ma and ＞5 km since 25 Ma (Figs.2d, e). Using geological constrains, the authors discussed these stages separately.

Stage 1: The Late Cretaceous crustal thickening of the Lhasa terrane was due to the Lhasa-Qiangtang collision induced by tectonic shortening (Murphy M et al., 1997), or this crustal thickening stage was likely related to the oceanic subduction of the Neotethys (Ding L et al., 2014; Lai W et al.,2019). Basin inversion and regional compression initiated during deposition of the uppermost Shexing strata (ca. 96 Ma), as indicated by active thrust faults and widespread accumulation of syntectonic conglomerates in the western part of the Lhasa block (Wang JG et al., 2020). Detrital zircons separated from modern river sands in the Gangdese belt,Lhasa terrane, reveal that the Lhasa crust has thicken to 60-70 km since the Cretaceous (Tang M et al., 2020).

Stage 2: The Lhasa terrane experienced strong magmatic input during 65-45 Ma, which led to the crustal thickening coeval with the rising elevation (Ma XX et al., 2021; Zhu DC et al., 2017). The 65-45 Ma intensive magmatism was caused by a rollback and consequential break-off of the subducted Neotethyan lithosphere, which represented the transition from oceanic to continental subduction (Huang et al., 2017). During 65-45 Ma, the crust of the Lhasa terrane thickened rapidly from 50 km to 65 km due to magmatic input (Fig. 2d; Ma XX et al., 2021; Tang M et al., 2020; Zhu DC et al., 2017). The ca. 40 Ma adakitic magmatism on the south Lhasa terrane suggested that the Lhasa crust was already thickened to ca.65 km (Guan Q et al., 2012).

Stage 3: In this model, the Lhasa terrane needs an additional ca. 10 km thickness after 45 Ma to obtain ＞75 km crust (Fig. 7b). Due to there were weak deformation and magma records (Figs. 2a, b), here, the Indian Plate underthrusting was proposed to explain additional crustal thickening and the elevation rising (Figs. 7b-1, b-3). The crustal thickness increase after 45 Ma was attributed to underthrusting of the Indian Plate (Styron R et al., 2015). The ca. 15 km thick Indian crust possibly retained its strength as it underthrusted beneath the Tibetan Plateau (Copley A et al,2011; Gao R et al., 2017; Freymueller JT, 2011). In addition,the India-Asia convergence rate dramatically decreased from ca. 16 cm/a at approximately 50 Ma to ca. 5 cm/a at approximately 35 Ma (van Hinsbergen DJJ et al., 2011),which has normally been suggested to be related to the resistance of the Indian Plate underthrusting (Yakovlev PV and Clark MK, 2014). After the underthrusting, the crustal thickness and paleoelevation of the Lhasa terrane exceeded the Qiangtang terrane (Figs. 2d, e). The underthrusting was followed by the foundering of thickened crustal root in the Lhasa and Qiangtang terranes (Chen M et al., 2017; Razi AS et al., 2016).

After 30 Ma, both the Lhasa and Qiangtang terranes continued to thicken with a small tectonic shortening rate(＜10%) (Fig. 2a; Yin A et al., 1999). Meanwhile, the Himalayas began to increase in elevation, and the plateau propagated towards both northern and southern Tibetan Plateau (Ding L et al., 2017).

Though constraints on elevation and erosion rate were too scattered to reconstruct continuous histories of the Qiangtang and Lhasa terranes, the collected data indicate that the high elevation did not enhance the erosion rate in the central Tibetan Plateau, which was probably due to intermontane drainage (Figs. 2e, f). This may have a profound influence on climate change during the Tibetan Plateau formation (Chen J et al., 2013; Garzione CN, 2008; Molnar P et al., 2010), or vice versa (van der Beek P, 2016).

6. Conclusion

Using the available published and newly estimated data,the authors affirmed a strong coupling between crustal thickening and elevation obtained in the central Tibetan Plateau since 120 Ma. Hence, the Tibetan Plateau was most likely built as a consequence of isostatic equilibrium from its over-thickened crust.

Based on the one-dimensional isostatic modeling results,the authors concluded that the Qiangtang Plateau was formed mainly by tectonic shortening induced by crustal thickening and a slow erosion rate. However, several processes have collaborated to reconstruct the Lhasa Plateau under different geological backgrounds, including the Qiangtang-Lhasa continental collision during 120-65 Ma, the large magmatic input during 65-45 Ma, and underthrusting of the Indian crust beneath the Lhasa terrane after 30 Ma. The authors also showed that significant crustal thickening triggered central Tibetan Plateau uplift and reduced long-term erosion after the India-Asia collision.CRediT authorship contribution statement

Zhong-bao Zhao and Xu-xuan Ma developed the theoretical formalism, performed the analytic calculations and performed the numerical simulations. Chao Li mainly collected data. All authors contributed to the final version of the manuscript.Declaration of competing interest

The authors declare no conflicts of interest.

Acknowledgement

This work is financially supported by the Second Tibetan Plateau Scientific Expedition and Research Program (STEP)(2019QZKK0901), the Natural Science Foundation of China(41672211), the China Geological Survey (DD20190059),and the Fundamental Research Funds for institute of geology,CAGS (JYYWF201810). The authors thank Dr. Jing-chao Li for helping with the data collection. The authors also thank three reviewers’ comments which helped this paper to improve a lot. Associated editor Dr. Xi-jie Chen and Dr. Cong Zhang was also appreciated by their efficient organizing and carefully-worded edits.

Supplementary dataset: The supplementary dataset contain the collected and reprocessed data: (1) the estimated magmatic input flux rate for the two terranes is shown in Table S1; (2) the historical crustal thickness of the southern Qiangtang terrane is shown in Table S2; (3) the published crustal shortening amount and time duration of the TP is shown in Table S3; (4) the published erosion rate of the Qiangtang and Lhasa terranes is shown in Table S4and (5) the paleoelevation of the Qiangtang and Lhasa terranes is shown in Table S5. References for supplementary materials are included in the ‘Supplementary data’ file. Supplementary data(Table S1, Table S2, Table S3, Table S4and Table S5) to this article can be found online at doi: 10.31035/cg2021013.