The incompressible continuity equation, which accounts for mass conservation using the Boussineq approximation, is as follows: ∂ ν i ∂ x j = 0 , (1) where i and j are spatial directions according to the Einstein summation convention, ν i is the component of the velocity vector, and x j is the spatial coordinate.

Momentum conservation is governed by the three-dimensional Stokes equations, which are given as follows: ∂ σ ′ i j ∂ x j = ∂ P ∂ x i − g i ρ ( T , P , C , M ) , (2) where σ ’ i j is the deviatoric stress tensor, g i is the gravitational acceleration, and ρ is the density, which depends on temperature (T), pressure (P), composition (C), and melt fraction (M).

The heat conduction equation for the purpose of energy conservation is ρ C P D T D t = ∂ ∂ x i ( k ( T , P ) ∂ T ∂ x i ) + H r + H a + H s + H L , H a = T α D P D t , H S = σ ′ i j ε ˙ i j , (3) where t is time; D / D t is the substantive time derivative; C P is the isobaric heat capacity; k is the thermal conductivity, which is temperature and composition dependent; α is the thermal expansion; ε ˙ i j is the strain rate tensor; H r and H s are radiogenic and shear heating; and H a and H L are adiabatic and latent heat production.

The deviatoric stress tensor σ ′ i j is related to the effective viscosity η e f f and the strain rate tensor ε ˙ i j via σ ′ i j = 2 η e f f ε ˙ i j , (4) ε ˙ i j = 1 2 ( ∂ v i ∂ x j + ∂ v j ∂ x i ) . (5)

The code considers the viscous and plastic rheologies, as well as the thermomechanical properties of different rocks. The lithospheric strength is determined by a combination of ductile and brittle deformation mechanisms at a mountain building time scale. The brittle deformation follows the Drucker–Prager yield criterion, which describes the linear relationship of the material resistance on the total pressure (Ranalli, 1995): σ y i e l d = C 0 cos ( φ e f f ) + P s i n ( φ e f f ) , s i n ( φ e f f ) = sin ( φ d r y ) ( 1 − λ ) , (6) η p l a s t i c = σ y i e l d 2 ε ˙ II , (7) where σ y i e l d is the yield stress, P is the dynamic pressure, C 0 is the cohesion at p = 0, φ is the internal frictional angle ( φ d r y ) (stands for dry rocks), λ is the pore fluid pressure factor and ε ˙ II is the second invariant of the strain rate, and η p l a s t i c is the viscosity for plastic rheology.

The viscosity for ductile creep takes the form (Beaumont et al., 2004)   η d u c t i l e = B ∗ ( ε ˙ II ) ( 1 n − 1 ) ⁡ exp ( E a + P V a n R T ) , B ∗ = f ∗ A D ( − 1 / n ) , (8) where E a is the activation energy, V a is the activation volume, n is the stress exponent, R is the gas constant, T is the absolute temperature, A D is the material constant, and f is a pre-exponential scaling factor. The factor f is applied to scale the effective ductile viscosity calculated from the reference flow laws.

The effective viscosity ( η e f f ) is defined as the minimum value between the plastic and ductile viscosities (Ranalli, 1995): η e f f = min { η p l a s t i c ,   η d u c t i l e } . (9)

Laboratory-determined flow laws of ‘wet quartzite,’ ‘plagioclase An₇₅,’ and ‘dry olivine’ are used for the continental upper crust, lower crust, and asthenospheric mantle, respectively (Ranalli, 1995). Note that a modulated ‘plagioclase An₇₅’ is used for the lower crust of the hinterland basin. The detailed rheological parameters used in this study are shown in Table 1. We set a lower cutoff viscosity of 10¹⁸ Pa s and an upper cutoff viscosity of 10²⁶ Pa s for all rocks.

The degree of the partial melting of rocks is calculated by P-T–dependent solidus and liquidus curves (Table 1). The volumetric fraction M of the melt is assumed to linearly increase between the solidus and liquidus temperatures at a given pressure (Burg and Gerya, 2005): M = { 0               w h e n   T ≤ T s o l i d u s   T − T s o l i d u s T l i q u i d u s − T s o l i d u s   w h e n   T s o l i d u s < T < T l i q u i d u s 1               w h e n   T ≥ T l i q u i d u s , (10) where T s o l i d u s and T l i q u i d u s are the solidus and liquidus temperatures of the rocks, respectively.

The effective density of partially molten rocks changes with the amount of melt fraction and P-T conditions: ρ e f f = ρ s o l i d − M ( ρ s o l i d − ρ m o l t e n ) , (11) ρ P , T = ρ 0 [ 1 − α ( T − T 0 ) ] [ 1 + β ( P − P 0 ) ] , (12) where ρ s o l i d and ρ m o l t e n are the densities of the solid and molten rock, respectively; ρ 0 is the density at P 0 = 0.1   MPa and T 0 = 298 K; and α ( 3 × 10 − 5 K − 1 ) and β ( 1 × 10 − 5 MPa − 1 ) are the thermal expansion and compressibility coefficients, respectively.

The effects of latent heating generated by melting or crystallization are accounted for by an increased effective heat capacity ( C p e ) and thermal expansion ( α e ) (Burg and Gerya, 2005): C p e = C p + Q L [ ( ∂ M / ∂ T ) P = c o n s t a n t ] , (13) α e = α + ρ Q L [ ( ∂ M / ∂ T ) T = c o n s t a n t ] , (14) where C p is the heat capacity of the solid rock and Q L is the latent heat of the melting rock (Burg and Gerya, 2005).

We investigate the topographic response of a hinterland basin imbedded in the upper plate under horizontal compression condition. The computational domain is a 4000 × 200 × 4000 km (in the order of x, y, and z) Cartesian box, and is resolved by 501 × 101×501 grid points with a uniform resolution of 8 × 2 × 8 km (Figure 3). There are approximately 200 million Lagrangian markers randomly distributed inside the box, which are used to advect the physical properties. As Chen L. et al. (2017) demonstrated that the lower crust rheology plays an important role in the deformation of the continental lithosphere. We vary the scaling factor f (Eq. (8)) of the lower crust flow law to represent the different strengths of the basin inside the plateau. Previous studies indicate that the composition of the continental lower crust ranges from basaltic to andesitic (Hacker et al., 2015) and the basaltic lower crust is ∼1 order of a magnitude more viscous than the andesitic lower crust (Shinevar et al., 2015). Therefore, in the reference model, we set a rectangle basin with the scaling factor f = 10, which means that the viscosity of its lower crust is 10 times higher than that of the surrounding regions. The basin is 500 × 800 km in size, and its location is varied in the study.

In the model, back and side walls are set to be free slip. We take 55 Ma as the initial continental collision time (Molnar and Stock, 2009; Hu et al., 2016), which corresponds to the onset of the model. A starting convergence rate of 10 cm/yr is imposed at the central region of the front wall (1000 < X < 3000). It gradually decreases to 5 cm/yr after 5 Ma, and is fixed at 5 cm/yr henceforth. The time-dependent convergence rate used here is generally consistent with the estimates from plate reconstructions or paleomagnetic data (Molnar and Stock, 2009; Copley et al., 2010; and Cande and Stegman, 2011) (Figure 3A). It should be noted that when the India–Asia collision began is still controversial. Recently, Xiao et al. (2017) proposed that the terminal India–Asia collision occurred after 14 Ma based on the anatomy of composition and tectonic nature of the Himalayas, which supports a multi-stage collision process between India and Asia. The material influx at the front boundary is limited in the central region to simulate the indentation of the Indian continent, beside which (x = 0–1000 km and x = 3000–4000 km) the free slip condition applies. To keep mass balance within the layer, a constant outflow velocity is applied at the upper boundary, which is determined by V = ( W ′ x W x   t a W z ) ×   V i , where t a is the thickness of the sticky air; W x and W z are the model widths in the x and z directions, respectively; and W x ' is the horizontal range on the front wall that the convergence is imposed. At the bottom boundary, an infinity-like external free slip boundary condition is imposed (Burg and Gerya, 2005). In order to simulate the quasi-free surface and hence allow for the topography development, a 20 km-thick “sticky air” layer, which is characterized by a density of 1 kg / m 3 and a viscosity of 10 19 Pa s, is placed above the rocky model (Schmeling et al., 2008).

The initial temperature first increases linearly from 0°C at the surface to the Moho with 500°C for the indenter and 600°C for the upper plate, and then continues to increase to 1,360°C at the lithosphere base. The thermal gradient gradually changes in the transition zone between the indenter and its surrounding upper plate. The underlying asthenospheric mantle has an initial temperature gradient of 0.5°C/km. The top thermal boundary maintains a constant temperature of 0°C, and all the vertical thermal boundaries are insulating (no horizontal heat flow). At the bottom thermal boundary, an infinity-like external constant temperature condition is applied (Burg and Gerya, 2005). This implies that a constant temperature condition is satisfied at ∼200 km below the model base, allowing a spontaneous adjustment of the temperature and heat flux at the bottom of the model.

Taking into account that multiple ocean closure and intercontinental suturing events occurred at the southern Asian continental margin prior to the India–Asia collision (Yin and Harrison, 2000; Kapp and DeCelles, 2019), a relatively weak plate is set to simulate the pre-collisional Asian continental margin. The Asian plate has a 40 km thick continental crust composed of a 17 km thick upper crust and a 23 km thick lower crust, and a 80 km thick lithospheric mantle. The indenter has a 35 km thick continental crust composed of a 15 km thick upper crust and a 20 km thick lower crust, and a 105 km thick lithospheric mantle, to simulate the northern moving Indian craton. In consideration of oceanic subduction before continental indentation, we set a weak zone with a dip angle of 30 in front of the indenter cutting through the entire lithosphere (Figure 3B). Our focus is on the topography evolution in the hinterland of the Tibetan Plateau, where the erosion and sedimentation is not as significant as that at the plateau margins. Here, we do not consider the effects of the surface processes on the topographic evolution, but we will discuss the possible limitations for simplification.

Figure 4 shows the topographic evolution of the reference model (model-1 in Table 2). The deformation propagates northward quickly and localizes around the borders of the strong basin soon after the convergence starts. From the cross-sections, the existence of the strong basin causes the surrounding crusts to intensely thicken and buckle, resulting in the development of a high topography around the basin. In contrast, the strong basin itself experiences little deformation and underthrusts beneath the surrounding thickened crust. As a consequence, a lowland area develops in the interior of the plateau, characterized by a notable topographic relief. As the convergence continues, the lowland gradually shrinks and migrates northward. The lowland exists in the interior of the plateau until ∼46 Myr, after which it merges into the surrounding high plateau.

The lowland area displays significant topographic relief in its interior. At the early stage of the collision, topographic depressions with negative elevation developed along the margins of the strong basin and uplifted belts developed adjacently (Figure 4). In terms of the uplift process, the elevation of a lowland does not rise monotonously. During the convergence, the distributed width of the strong lower crust changes slightly. Here, we assume that the length of the basin’s lower crust remains unchanged. This allows us to track the elevation variations of the basin’s different regions (points A1, B1, C, B2, and A2) with time (Figure 5). The surface of the lowland uplifts progressively from the margins to the center of the strong basin. The uplift histories of different regions vary significantly, but all undergo a considerable elevation drop followed by a rapid uplift. It is noteworthy that the lowland center (Point C) maintains a constant elevation for a long period and then experiences a noticeable elevation drop, followed by a quick uplift to the elevation of the surrounding mountains (Figure 5).

In this section, we design a group of models (models two to three in Table 2) to investigate the impacts of a strong basin’s position on the topographic evolution of the plateau. Model-2 and model-3 are identical to the reference model, except that the distances of the strong basin to the convergent boundary are 0 and 1,000 km, respectively, instead of 500 km in the reference model. Figure 6 shows the topographic evolution of model-2 and model-3. The strong basins in model-2 and model-3 represent the rigid basins located at the south and north Tibetan Plateau, respectively. The plateau growth patterns in model-2 and model-3 are similar to the reference model, with lowlands developing in the interior of strong basins and plateaus forming in the surrounding region. Differently, the formation of the hinterland basin in model-2 occurs roughly 10 Myr earlier than that in model-3 (Figure 6). In addition, the time of the lowland merging into the plateau in model-2 also occurs earlier than that in model-3, at approximately 5 Myr. At 38 Myr of convergence, the lowland in model-2 has shrunk into a small residual basin, whereas the lowland in model-3 is still broad (Figure 6).

The topographic features in the interior of the lowland region are significantly different in the model-1, model-2, and model-3 (Figure 7). In model-2, the center region of the lowland has been elevated from the early stage and displays a locally positive topography. Specifically, the lowland center (Point C) in model-2 uplifts nearly 2,000 m after 20 Myr of convergence, whereas the height of the lowland center in model-1 remains unchanged until 32 Myr of convergence. In addition, in model-2, the depression formed at the south margin of the lowland is more visible and deeper than that in the reference model. In detail, the subsidence of the south depression in model-2 exceeds 3,000 m, compared to less than 2,000 m in the reference model. By contrast, the lowland in model-3 has similar topographic features to the reference model.

Considering the fact that the hinterland basins inside the Tibetan Plateau vary in size, we designed a model (model-4, Table 2), in which the width of the strong basin is half that in the reference model, to investigate the effects of the strong basin’s width on the topographic evolution of the plateau. Figure 8 shows the topographic evolution and the lowland uplift history of model-4. The topographic growth pattern of model-4 is similar to that of the reference model, except that the lowland in model-4 vanishes nearly 10 Myr earlier than that in the reference model (Figure 8A–D). As a result, it takes only 40 Myr to form a large uniform plateau in model-4 (Figure 8E), instead of >45 Myr in the reference model.

The topographic features and uplift history of the lowland in model-4 are dramatically different from those observed in the reference model. The lowland center region in model-4 was elevated soon after the convergence started, rather than maintaining constant in the reference model (Figure 8F). To be specific, after 10 Myr of convergence, the lowland center (Point C) in model-4 had increased to a height of 3,000 m and then kept unchanged until 20 Myr. Then, the height dropped by more than 2,500 m, close to sea level after 30 Myr of convergence (Figure 8F). In contrast, the elevation of the lowland center in the reference model was nearly constant during the first 30 Myr. The elevation fall of the lowland center in model-4 exceeds 2,500 m, whereas this value in the reference model is less than 1,000 m.

In this section, we designed a model with a less strong basin (f = 2, model-5 in Table 2) and a weak basin (f = 0.5, model-6 in Table 2) to test the effects of the basin’s strength on the topographic expressions of the plateau. The model results showed that the less strong basin had similar effects on the topographic growth of the plateau to that produced in the reference model (Figure 9A–E). The lowland basin also survives in the interior of the plateau for approximately 30 Myr. However, the basin’s uplift history is quite different to that in the reference model. The lowland center does not experience an elevation drop at a late stage, but keeps rising until reaching the height of the surrounding plateau (Figure 9F).

The topographic response to a weak basin differs significantly from that of a strong basin. In model-6, the weak basin undergoes intense deformation in the early stage of the collision and uplifts earlier than the surrounding areas, leading to the formation of a topography high in the interior of the upper plate (Figure 10). The highland continues to rise and maintains its size until 20 Myr. After that, it propagates southward and northward. After 30 Myr of convergence, the plateau develops a flat-surfaced morphology. This set of models illustrates that the strength of the hinterland basin is a key factor in controlling the uplift pattern of the plateau. A strong hinterland basin favors a later uplift of the plateau’s interior, whereas a weak hinterland basin facilitates the formation of a topography high in the early stage of collision.

Different to the aforementioned models that test the inherent properties of the basin, we perform an additional model (model-7 in Table 2) to test the influence of the convergence rate in the topographic response of a hinterland basin. Here, we choose the convergence rate reconstructed by Cande and Stegman (2011), in which the velocity is 13 cm/yr in the first 5 Myr and is then gradually reduces to a fixed value of 5 cm/yr after 15 Myr, keeping the other parameters the same as the reference model. The model results show that the growth pattern of the plateau and the topographic features in the interior of the lowland are similar to that in the reference model (Figure 11A–E). The lowland basin develops soon after the initial collision and survives for ∼35 Myrs in the interior of the plateau. The basin center also experiences an evident elevation drop and rapid uplift before merging into the plateau (Figure 11F). It is necessary to mention that the lowland basin migrates northward and shrinks more quickly due to the increased convergence amount. Also, it is worth noting that the lowland center was elevated about 1,000 m at an early stage of the basin’s lifespan, which is not noticeable in the reference model (Figure 11F). By comparing the faster convergence model with the reference model, we can conclude that the initial convergence velocity has little influence on the major characteristics of the basin’s topographic evolution and plateau growth.

Our models demonstrate that the strength of the hinterland basin exerts a first-order control over the topographic evolution. The basin with a rigid basement develops into a lowland surrounded by high mountains soon after the initial collision and survives in the interior of the plateau for a long period of time (30–40 Myrs). In contrast, the basin with a soft basement uplifts quickly and reaches its current elevation after 20 Myr of the convergence. This provides an explanation to why some basins inside the Tibetan Plateau have long-term sedimentary records and keep a low elevation until the Late Oligocene, such as the Lunpola and Hoh Xil basins (Wang et al., 2008; Sun et al., 2014; Sun et al., 2015; Wu et al., 2017; and Fang et al., 2020), while some other basins ceased to deposit and have uplifted during the Eocene, such as the Linzhou and Shuanghu basins (Ding et al., 2014; Currie et al., 2016). Depending on the locations and the widths, the rigid hinterland basins display quite different lifetimes in the interior of the plateau. The time required for the formation of a lowland basin in the north Tibetan Plateau (model-3, Dis = 1000 km) delays by 10 Myr in comparison to the basin at central Tibet (model-1, Dis = 500 km). The time of merging into the plateau for the narrow basin (model-4, width = 250 km) is approximately 10 Myrs earlier than a vast basin (model-1, width = 500 km).

The hinterland basins also exhibit a variety of topographic features and uplift histories, depending on their locations, widths, and strengths. When the strong basin is close to the collisional zone (model-2) or is narrow (model-4), the central region rises at the early stage and then experiences an evident elevation drop at the late stage (Figure 12). However, when located further north (e.g., model-3), the lowland basin displays low-relief topography in the interior and the elevation keeps unchanged for a long period of time (Figure 7D). The topographic feature of the basin is closely related to the deep tectonic structure. Depressions and uplifted belts formed owing to the underthrusting and bending of the strong lower crust. Therefore, due to a weak surface topographic response to the bending, the less strong basin (model-5, f = 2) displays a less topographic relief in the interior (Figure 12). In model-2 and mode-4, the uplifted belts occur in the lowland center owing to the more intense bending and the narrowness of the strong basin, and the elevation drop is actually the result of the marginal depressions merging as the basin gradually shrinks (Figure 12).

The Lunpola Basin is a narrow basin and is located in the northern Lhasa terrane (Figure 1). Geophysical observation and geochemical data show that it is characterized by fast velocity anomalies in the crust (Yang et al., 2012) and negative Hf isotopic values of felsic igneous rocks (Hou et al., 2015; Hou et al., 2020) (Figure 2), which indicate that the Luopala Basin has a strong and old basement. Therefore, the uplift history of the Lunpola Basin can be compared with that of the lowland in model-4, in which the strong basin is 250 km in width and located 500 km north of the convergent boundary. A topographic analysis of model-4 shows that the basin was elevated to a medium elevation after 10 Myr of convergence and experienced a noticeable elevation drop of ∼3,000 m from 20 to 30 Myr of convergence, then followed by a rapid uplift (Figure 8F). We compiled the paleoelevation data published over the past 2 decades in the Lunpola Basin, including oxygen isotope, hydrogen isotope, and paleontology (Rowley and Currie, 2006; Polissar et al., 2009; Sun et al., 2014; Jia et al., 2015; Wu et al., 2017; Farnsworth et al., 2018; and Fang et al., 2020) (Figure 13A). The estimates based on isotope data tend to predict a higher elevation, whereas fossil-based data tend to reflect lowland elevations (Spicer et al., 2021). Fossil-based data indicate that the Lunpola Basin was elevated to a medium elevation during the Eocene (50–38 Ma) and suffered an evident elevation drop during the Oligocene (38–25 Ma), and then followed by a rapid uplift in the Early Miocene (23 Ma) (Figure 11). The elevation fluctuation of the Lunpola Basin, as revealed by fossils, is in agreement with the uplift history of the lowland center predicted by our model-4 (Figure 13A). In addition, the isotopic data also suggest that the Lunpola Basin had a high elevation during the Eocene and suffered a rapid uplift during the Early Miocene. Therefore, the Lunpola Basin was likely to be elevated to 2,500–3,000 m in the Middle Eocene and suffered an evident elevation fall of ∼2,000 m during the Oligocene, and then experienced a rapid uplift to achieve its present height in the Early Miocene.

The Hoh Xil Basin is a large basin far from the Indus–Yarlung suture zone (IYS). Weak crustal deformation (Staisch et al., 2016) and the negative Hf isotopic feature of felsic igneous rocks (Hou et al., 2020) imply a mechanically strong basement beneath the Hoh Xil Basin. Therefore, the topographic evolution of the Hoh Xil Basin can be compared with the lowland in model-3, where the strong basin is wide and far from the collisional zone. Our results show that the lowland center in model-3 keeps a constant elevation for ∼30 Myr and then experiences a slight elevation drop before the rapid uplift (Figure 7D). In the Hoh Xil Basin, a compilation of paleoelevation data indicates that the elevation did not change significantly throughout the Eocene (Cyr et al., 2005; Polissar et al., 2009; Miao et al., 2016; and Song et al., 2021) (Figure 13B). This is in accordance with the long-term unchanged elevation of the lowland center in model-3. As expected by the uplift histories of two sub-basins within the Hoh Xil Basin (Liu et al., 2016) (Figure 13B), the Hoh Xil Basin experienced a slight elevation drop during Oligocene followed by a fast uplift in the Early Miocene. This is consistent with the late-stage elevation variation of the lowland center in model-3.

How the Tibetan Plateau uplifted is a subject of intense debate. Several hypotheses for the growth pattern of the Tibetan Plateau have been proposed, including: 1) northward stepwise uplift model (Tapponnier et al., 2001; Mulch and Chamberlain, 2006), 2) a central proto-Tibetan Plateau model (Wang et al., 2008; Wang et al., 2014), and 3) the “Central Tibetan Valley” model (Deng and Ding, 2015; Farnsworth et al., 2018; and Xiong et al., 2022). According to the stepwise uplift model, the Tibetan Plateau rises gradually from the south to the north, with the Lhasa rising by Eocene, the Qiangtang rising during the Oligocene, the Songpan–Ganzi rising during the Miocene, and the Qilian rising during the Pliocene to Quaternary (Mulch and Chamberlain, 2006). The proto-Tibetan Plateau model argues that the Lhasa and southern Qiangtang terranes were elevated into a central Tibetan highland during the Middle Eocene, and then the highland expanded to the Himalayas during the mid-Miocene, and to the Hoh Xil Basin during the Early Miocene (Wang et al., 2008; Wang et al., 2014). However, the “Central Tibetan Valley” model supports a lowland central Tibetan along the Bangong–Nujiang suture zone (BNS) sandwiched between the Gangdese and the Tanggula mountains (or central watershed) during the Early Eocene to the Late Oligocene, and then uplifted in the early Miocene (Farnsworth et al., 2018; Fang et al., 2020; Xiong et al., 2022).

Based on the predicted uplift histories of the Lunpola and Hoh Xil basins in this study and previous paleoelevation data, we propose a revised Cenozoic uplift model for the Tibetan Plateau south of the Kunlun Mountain (Figure 14). We subdivide the Tibetan Plateau’s growth into four stages according to the surface topographic evolution. At the first stage (50–38 Ma), the Tibetan Plateau is characterized by two highlands separated by three lowlands (Figure 14A). Paleoelevation data from the Linzhou and Namling–Qiyug basins suggest that a high Andean-type Gangdese mountain range (∼5000 m) developed in the southern Lhasa during the Early Eocene (Ding et al., 2014; Ingalls et al., 2017). Also, the central watershed highlands were contemporaneously elevated to 4,000–5,000 m in central Qiangtang, according to the paleoelevation data of the Heihuling and Gonjo basins (Xu et al., 2013; Xiong et al., 2020). The fossil evidence indicates the existence of a lowland with an elevation of <2,000 m south of the Gangdese Mountains (Ding et al., 2017), a lowland region with an elevation of 2,500–3,000 m between the Gangdese and the central watershed highlands (Wei et al., 2016; Fang et al., 2020; Su et al., 2020), and a broad lowland region with an elevation of ∼2,000 m north of the central watershed highlands during the Eocene (Miao et al., 2016; Song et al., 2021).

At the second stage (35–25 Ma), the central Tibetan lowland along the Bangong–Nujiang suture (BNS) experienced an elevation drop of ∼2,000 m, forming a central Tibetan deep valley (Figure 14B). Plant and foraminifera fossils indicate that the elevation of the Lunpola, Nima, and Gerze basins were no more than 1,000 m during the Late Oligocene (Wei et al., 2016; Wu et al., 2017). Sedimentary records show a noticeable facie transition from lacustrine mudstones to fluvial and alluvial conglomerates and a threefold increase of the sedimentary rate in the Lunpola Basin during the Late Oligocene (Fang et al., 2020). This consolidates our prediction that the central Tibetan lowland suffered an evident subsidence during the Oligocene. However, paleoelevation data reveal that the lowland south of the Gangdese Mountains and the lowland north of the central watershed mountains remained at an elevation of ∼2,000 m at this stage (Ding et al., 2017; Dai et al., 2019). The elevation drop of the central Tibetan lowland occurred soon after the reduction of the India–Asia convergence rate (Molnar and Stock, 2009; Copley et al., 2010; and Cande and Stegman, 2011). Thus, we suggest that the slow convergence rate may result in a relaxation of the lithospheric compression bending, thereby sinking the surface.

At the third stage (25–18 Ma), the central Tibetan valley suffered a rapid uplift and merged into the plateau (Figure 14C). Fossil evidences from plants, mammals, and fish indicate that the ecosystem of the Lunpola and Nima basins experienced a significant turnover from a warm tropical environment in the Late Oligocene to a cool climate since the Early Miocene (Deng et al., 2019; Deng et al., 2021), whereas the lowland south of the Gangdese Mountains and the lowland north of the central watershed mountains remained at a low elevation according to carbon isotopic and palaeobotanical evidences (Wang et al., 2006; Sun et al., 2015; and Li et al., 2020). Seismic evidence suggests the convective removal of thickened lithosphere under the central Tibet during the Early Miocene (Chen M. et al., 2017), implying that the lithospheric delamination and mantle upwelling may have played important roles in the uplifting of the central Tibetan valley (Xiong et al., 2022). A recent study reveals an Early Miocene transition from a cold and strong middle-lower crust to a hot and weak middle-lower crust in the Qiangtang Block (Zhang et al., 2022). This supports that the Early Miocene uplift of the central Tibetan valley was driven by the crustal flow from surrounding highlands. Thus, we suggest that the lateral crustal flow and the vertical lithospheric delamination and mantle upwelling are the primary drivers of the central Tibetan valley’s rapid uplift.

At the last stage (18–0 Ma), the lowland south of the Gangdese Mountains and the lowland north of the central watershed mountains were uplifted and merged into the plateau since the Middle Miocene, forming the topography similar to the present Tibetan Plateau (Figure 14D). Oxygen isotopic data from the Gyirong Basin and the Thakkhola graben indicate that the Tethyan Himalaya was elevated to its current elevation by the Middle Miocene (Garzione et al., 2000; Rowley et al., 2001). Oxygen and hydrogen isotopes-based paleoelevations show that the Hoh Xil Basin was elevated to a height of >4,000 m after the Early Miocene (Polissar et al., 2009). The small degree of shortening and the Late Oligocene cessation of deformation in the Hoh Xil Basin imply that the north Tibetan lowland may have been uplifted as a result of the crustal flow and lithospheric removal (Staisch et al., 2016; Li et al., 2017).