I use I2VIS, which combines a conservative finite difference method and a marker-in-cell technique, to solve the governing equations of momentum, mass, and energy in incompressible flow (Gerya and Yuen, 2003). Rock properties are advected on the moving Lagrangian markers. A composite rheological model described below is employed to account for plastic deformation and creeping flow of rocks. The method considers the effects of adiabatic, shear, latent and radioactive heating. A full description of the method, including the implementation of partial melting, can be found in the textbook by Gerya (2010). All the related parameters used in the study follow the Table 17.2 in Gerya (2010).

Figure 2 shows the initial configuration and thermomechanical state of our experiments. The whole computational domain is 4,000 km wide and 820 km deep, and is resolved with a non-uniform 2041 × 481 rectangular grid. In the horizontal direction, the model domain between 1,000 and 3,000 km has a high resolution of 1 km. The two transitional zones extend 100 km away from the high-resolution region, where the resolution reduces from 1 to 5 km. The grid resolution for the rest model domain where no obvious deformation occurs is 5 km. In the vertical direction, the resolution for the upper 370 km is 1 km, for the depth range 370–470 km reduces from 1 to 5 km, and for the rest depth is 5 km. Over 32 million Lagrangian markers randomly distribute in the whole model domain. These markers are moved at each time step using the computed velocity field.

The initial model consists of a 2000-km-wide weak block in the middle and two strong blocks on either side (Figure 2A). The weak block is composed of a 55-km-thick crust (a 15-km-thick upper crust, a 20-km-thick middle crust, and a 20-km-thick lower crust) and a 65-km-thick mantle lithosphere, representing a juvenile orogen. The strong block is composed of a 40-km-thick crust (a 20-km-thick upper crust and a 20-km-thick lower crust) and a 100-km-thick mantle lithosphere, representing a confining craton. The design of a hot orogen with a relatively thick crust allows that it is sufficiently weak and more deformable, compared to the confining cratons. A small inverted-triangular seed is imposed at the middle of the orogen crust bottom, which allows for the deformation initiating at the center of the model. The initial geotherm linearly increases from 273 K at the model surface to 1623 K at the lithosphere base. A mantle adiabat with potential temperature 1573 K and adiabatic gradient of 0.5°C/km is used for the sub-lithosphere mantle. The mechanical boundary conditions are free slip for all the boundaries. The thermal boundary conditions are constant temperature (273 K) on the top, remote fixed temperature on the bottom (Gerya et al., 2006), and insulating (no horizontal heat flow) on the two sides. Horizontal shortening is implemented by imposing varying convergence rates ( V x ) at X = 500 km and X = 3,500 km (Figure 2A).

To account for surface topographic evolution, a 20-km-thick sticky air layer is added on the top of the rocky portion of the model. It is characterized with low density (1 kg/m³) and low viscosity (10¹⁸ Pa s). The large density and viscosity contrast between the sticky air and rocky domain ensures small shear stresses (<10⁴ Pa) along the interface and makes it an internal free surface (Schmeling et al., 2008). For simplicity, our models do not consider surface erosion and sedimentation processes.

The strength of the lithosphere is controlled, on the geologic timescale, by the combination of both brittle and ductile deformation mechanisms. The brittle part of the lithosphere can be described by the Drucker-Prager yield criterion which expresses the linear dependence of the materials resistance on the total pressure (Gerya, 2010): σ y i e l d = C + P s i n ( φ e f f ) (1) sin ( φ e f f ) = sin ( φ ) ( 1 − λ fluid ) (2) λ fluid = P fluid P total (3) η p l a s t i c = σ y i e l d 2 ε ˙ I I (4) where σ y i e l d is the yield stress; P is the total dynamic pressure; C is the cohesion at P = 0 ; φ is the internal frictional angle; λ fluid is the pore fluid pressure factor that reduces the yield strength of fluid-containing porous or fractured rocks; ε ˙ I I = ( 0.5 ε ˙ i j ε ˙ i j ) 1 / 2 is the second invariant of the strain rate tensor; η p l a s t i c is the viscosity for plastic rheology. φ e f f can be illustrated as the effective internal frictional angle that integrates the effects of internal frictional angle and pore fluid coefficient. In the modeling, a cohesion of 1 MPa is used for all the material. sin ( φ ) = 0 is used for the air, water and partially molten rocks, sin ( φ ) = 0.15 for all the solid crust rocks, and sin ( φ ) = 0.6 for all the mantle rocks.

The viscosity for ductile creep is given by (Hirth and Kohlstedt, 2003): η d u c t i l e = ( ε ˙ I I ) ( 1 − n ) / n A D − 1 / n C O H − r / n e x p ( E a + P V a n R T ) (5) where A D (material constant), E a (activation energy), V a (activation volume), n (stress exponent) are experimentally determined flow law parameters. C O H is the water content in ppm H/Si, r is the water content exponent (Hirth and Kohlstedt, 2003), R is the gas constant, and T is the absolute temperature. The rheological parameters ( A D ,   E a ,   V a ,   n ) are derived from laboratory studies of deformation of rocks.

The effective viscosity is determined by comparing brittle/plastic and creep viscosities as a function of depth (Ranalli, 1995): η e f f = m i n { η p l a s t i c ,   η d u c t i l e } (6) The application of laboratory-derived rheological parameters to deformation at geological time and space scales needs to extrapolate the laboratory strain rates over many orders of magnitudes. Burgmann and Dresen (2008) show that laboratory-based rheologies are well consistent with geodetic measurement and field observations, indicating that laboratory-derived rheologies can be used as a good description of geological deformation. However, there is significant uncertainty in the laboratory-derived rheologies of the continental lower crust. For example, the lower crust yield stress calculated with the laboratory-derived flow laws of “diorite” and “plagioclase” (Ranalli, 1995) is less than 100 MPa under the model condition (Figure 2B). When I use “wet clinopyroxene” (Zhang et al., 2006) or “diabase” (Mackwell et al., 1998) as the lower crust flow law, the yield stress reaches several hundreds of MPa.

In the modeling, the rock rheologies are based on the laboratory-established flow laws. In all models, I use a “wet quartzite” flow law for the continental upper crust, “plagioclase An75” for the craton lower crust and orogen middle crust, and “dry olivine” for the lithospheric and asthenospheric mantle (Ranalli, 1995). I examine changes in the orogen lower crust rheology by varying the laboratory-derived flow laws. All partially molten rocks are assigned to a reduced effective viscosity of 10¹⁸ Pa s (Pinkerton and Stevenson, 1992) which provides large viscosity contrast with the surrounding solid rocks. Table 1 lists the rheological parameters tested in this study.

Figure 3 shows the evolution of the model “dela” (see Table 2). This model uses a lower crust flow law of “diorite” and a half convergence rate of 2.5 cm/yr. The convergence of the strong cratonic lithosphere causes the orogen to undergo pure shear shortening. After 5 Myr of convergence (250 km of shortening), lithospheric thickening preferentially occurs at the center of the model (Figure 3Ai). Owing to the notably low viscosity of the orogen lower crust relative to its upper and lower layer, the upper crust is decoupled from the lithosphere mantle (Figure 3Aii). By 14 Myr of shortening (700 km), the orogen crust experiences marked folding and buckling, forming a number of V-shape strain localizations, and has a thickness of ∼70 km (Figure 3B). However, the orogen lithosphere mantle thickening does not concur with the crust. Convergence forces the orogen lithosphere mantle to slide along the weak orogen lower crust and sink downward at the middle of the model (Figure 3B), resulting in a large piece of downwelling similar to the two-sided subduction shown by Gerya et al. (2008). The downwelling undergoes necking in a short time span (<1 Myr) due to its negative buoyancy, and finally breaks off from the orogen crust (Figure 3C), which I refer to as delamination. The delamination of the lithosphere mantle makes the orogen crust expose directly to the asthenosphere, leading to partial melting of the lower crust (Figure 3D). The replacement of the lower crust flow law with “mafic granulite” (Wang et al., 2012) does not change the overall behavior (see Effects of Convergence Rate). This is because both “diorite” and “granulite” fall into the weak end of the lower crust rheology (Figure 1). Note the flow law “mafic granulite” presented by Wang et al. (2012) is significantly weaker than the early one presented by Wilks and Carter (1990) due to different experimental samples and conditions (see Wen et al., 2021 for details).

This model (“delc” in Table 2) differs from the reference model only in that the flow law of the orogen lower crust is changed to “plagioclase”, representing an intermediately strong lower crust. Because the viscosity of the orogen lower crust increases relative to that in the reference model, lithospheric thickening is no longer concentrated only at the center of the model (Figure 4A). Instead, localized thickening occurs at the orogen edges and center, forming three separate lithospheric bulges (Figure 4B). After 20 Myr of convergence, the strong craton lithosphere indents into the orogen lithosphere (Figure 4C). In comparison to the reference model, crustal folding and buckling are less pronounced. Localized thickening accommodates most of the horizontal shortening, but thickened lithosphere mantle remains below the orogen.

In this model (“delf” in Table 2), “diabase” is used as the flow law of the orogen lower crust and the other parameters are identical to those in the reference model. As mentioned above, the “diabase” flow law represents the strong end member of the lower crust rheology according to the available laboratory results (Figure 1). By 5 Myr of convergence (250 km), horizontal shortening first causes the whole orogen lithosphere to fold in a wavy manner with a wavelength of ∼250 km (Figure 5A). With continued convergence, orogen crust folding becomes more pronounced and the strong craton lithosphere starts to underthrust below the orogen (Figure 5B). In contrast to the reference model, no localized lithospheric thickening occurs in this model owing to enhanced resistance of the orogen to shortening. By 20.4 Myr of convergence (∼1,000 km), the wavelength of folding decreases to ∼200 km and the strong craton lithosphere subducts beneath the orogen at both sides. In this case, horizontal shortening is mainly accommodated by underthrusting of craton mantle lithosphere below the orogen, and orogen crustal thickening only plays a secondary role due to its high strength.

The surface response of lithosphere to convergence is associated with the strength of the crust, in the models. Figure 6 shows the comparison of the topography evolutions predicted by the above models. In the weak lower crust case (Figure 6A), the central region of the orogen first subsides below the sea level with an amplitude of ∼4 km. This is resulted from the loading on the surface of the denscending mantle lithosphere beneath the orogen center. Topography highs at sides of the orogen arise from the compression between the craton and orogen. By 14.7 Myr, the downwelling mantle lithosphere breaks off and sinks into the deep mantle (see Figure 3C). This results in an abrupt rebound of the central depression, rapidly attaining an elevation of ∼2 km (Figure 6aii). The negative surface defection neighboring the topography highs at the sides of the orogen reaches ∼4 km, forming two foreland deeps. After 20 Myr of convergence, the orogen is characterized by three marked topography highs with two at boundaries and one at the central region, reaching ∼4,000 m, while a large area remains flat with no elevation. The surface expression of the intermediately strong orogen to horizontal shortening is different from the weak case (Figure 6B). Because neither wholesale mantle lithosphere descending nor delamination occurs in this case, the central region of the deformed area does not change significantly, and does not develop topography. By 20 Myr, two topography highs, which has a width of ∼200 km and an elevation of ∼4 km, develop at both sides of the orogen. In contrast, the orogen interior is relatively flat and remains low. In the strong orogen case (Figure 6C), two notable topography highs also develop at the orogen sides. They are very high in comparison to the weak and intermediately strong cases, which I ascribe to the absence of erosion in the models. Interestingly, the region bounded by the two topography highs alternates between negative and positive surface deflection, forming a couple of evenly-spaced topographic ribbons. Most of the undulated surface is close to or below the sea level, thus they have potential to be infilled with the material eroded from the two topography highs.

In summary, the topography expression during collisional orogeny depends on the strength of the orogen itself. The above-mentioned topographic features can be used as surface diagnostics to distinguish the nature of orogeny.

The convergence rate has an important control on the deformation of the lithosphere. With increasing convergence velocity, the response of the weak lithosphere shifts towards the response of a stiff lithosphere, although delamination occurs nevertheless.

Figure 7 compares the evolution of the orogen shortened by different convergence rates, showing that, for the same rheology, the evolution stages are the same, yet occurs at different times. In these two models (“delb” and “delj” in Table 2), “mafic granulite” is used as the flow law of the orogen lower crust and the other parameters keep the same as those in the reference model. At the early stage of convergence, all experiments show a similar deformation behavior, which are comparable to the weak crust model presented before. The differences become remarkable when the amount of convergence exceeds ∼500 km. For a fast convergence rate (V_x = 2.5 cm/yr; see Figure 7A), the amount of convergence required for lithospheric mantle delamination is ∼1,080 km, which occurs after 22 Myr from the model start. When a slower convergence rate is used (V_x = 1.0 cm/yr; see Figure 7B), the amount of convergence required decreases to ∼820 km, and delamination occurs at 41 Myr from the model start. Thus, the slower the convergence rate, the later the delamination occurs, and the less convergence is required to reach similar geodynamic outcomes. This second-order difference is due to the considerably greater volume of the lower crustal melting in the slow convergence, enhanced by larger radiogenic heating. The slower convergence enhances the crustal flow and inhibit folding and faulting (Figure 7). Therefore, a faster convergence generates a more rugged topography throughout the evolution of the models, while a slower convergence produces a relatively smoother topography (Figure 8).

The models predict that the orogen lithosphere mantle is removed as a coherent slice in response to horizontal shortening when the lower crust is sufficiently weak (e.g., Figures 3, 7). This differs both from the conventional delamination model and convective removal model (Figure 9). The former argues that the lithosphere mantle peels away from the crust via an elongated conduit (Bird, 1979; see Figure 9A), and the latter suggests the cold and dense lithosphere mantle drips off in an un-plate-like manner (e.g., Houseman et al., 1981; see Figure 9B). In the new delamination model presented here, it is the weak lower crust that plays a critical role in decoupling the deformation of the upper crust and lithosphere mantle and detaching of the thickened lithosphere mantle as a coherent slab (Figure 9C). The weak lower crust acts like a lubricant which decouples the deformation of the orogen’s upper crust from its lithosphere mantle, allows sliding of the lithosphere mantle along the base of the crust, and promotes delamination of the lithosphere mantle as a coherent slice. Once a large scale lithospheric downwelling is formed, the delamination process occurs rapidly and can be completed in 2 Myr, which is much more efficient than the manner of viscous dripping (∼25 Myr; see Gogus and Pysklywec, 2008). In addition, the punctuated surface uplift due to lithosphere delamination is limited to a localized area (∼200-km wide) above the delaminated lithosphere, not in a broad region as previously demonstrated (e.g., Gogus and Pysklywec, 2008). Because the convergence is imposed at both sides of the orogen, the elevated region is not mobile as shown by Krystopowicz and Currie (2013), but becomes narrow as the convergence proceeds. A comprehensive investigation of the surface expression resulted from delamination will be presented elsewhere.

In this section, I apply the model results to understanding the mantle dynamics of the Tibetan Plateau. Because the models are not specially designed for the Tibetan plateau, the comparisons are only made in an intuitive and qualitative way.

While the Tibetan crust has been uniformly thickened to double normal thickness throughout the India-Asia convergence, the Tibetan mantle lithosphere underwent a more complex evolution. According to tectonic reconstructions, ∼1,000 km of north-south shortening at least, including the eastward extrusion, is accommodated by the Tibetan crust (Guillot et al., 2003; Replumaz and Tapponnier, 2003). This implies an equal amount of mantle lithosphere to be accommodated beneath the Tibetan plateau, which is proposed in the frame of two different scenarios. One suggests that the shortened lithosphere mantle is found beneath the plateau, i.e., it has not detached (Priestley et al., 2006). Evidence for this comes from a multimode surface waveform tomography study, showing high-velocity and cold material underlying most of Tibet to a depth of 225–250 km (Priestley et al., 2006). The other scenario invokes convective removal of Tibet’s thickened mantle and its replacement by a hot asthenosphere (Molnar et al., 1993; Turner et al., 1993). This inference is mainly based on the inefficient Sn propagation and low Pn velocities observed beneath northern Tibet (Barazangi and Ni, 1982; McNamara et al., 1995; McNamara et al., 1997).

Recent seismic studies reveal that the mantle lithosphere under the Tibetan plateau is strongly heterogeneous (Figure 10). P-wave travel time tomography reveals low speeds in the shallow mantle beneath much of central and eastern Tibet and decrease of the horizontal distance over which the Indian lithosphere slides northward beneath the plateau from west to east (Li et al., 2008). Surface wave studies show that shear-wave speeds within the mantle lithosphere are high beneath western and southwestern Tibet and low beneath central and eastern Tibet (Agius and Lebedev, 2013; Schaeffer and Lebedev, 2013), suggesting the Tibetan lithosphere is cold and thick in western Tibet and warm and thin in central and eastern Tibet. More recently, Chen et al. (2017) showed a T-shaped high wave speed structure beneath southern to central Tibet using adjoint seismic tomography method and interpreted it as an upper mantle remnant from earlier lithospheric foundering. More detailed P-wave travel time tomography from INDEPTH II and III data reveals a subvertical high-velocity structure extending to ∼400 km beneath southern Tibet, north of the Yarlung-Zangbo suture, which was interpreted as downwelling Indian mantle lithosphere (Tilmann and Ni, 2003). Such a high-velocity anomaly has been confirmed by the global P wave tomographic model (Li et al., 2008; also see Figure 10). A similar subvertical high-velocity structure was also observed beneath the western syntaxis of the India-Asia collision zone (Negredo et al., 2007; also see Figure 10). In addition, recent S receiver function results reveal a thin Tibetan lithosphere (∼100 km depth) overriding a northward subducting Indian lithosphere (∼250 km) and a southward subducting Asian lithosphere (∼200 km) (Zhao et al., 2011). This indicates strong heterogeneities and rules out a uniform lithosphere, as Priestley et al. (2006) suggested.

According to the model results, how the mantle lithosphere is accommodated during horizontal shortening is associated with the lower crust rheology. At the early stage of the India-Asia convergence, the Tibetan crust underwent minor shortening and was not thickened significantly. The widespread Cretaceous and Palaeocene red bed sediments covering much of Tibet (e.g., Yin, 2010) provides evidence for the Tibetan crust having a normal thickness at that time, although pre-collision crustal thickening possibly occurred in the southern Lhasa terrane (Murphy et al., 1997; Kapp et al., 2005). Thus, the crust beneath Tibet remained relatively cold and strong, allowing for transmitting the Indian horizontal forces to the far north. The accommodation of the lithosphere mantle shortening mainly takes the form of strong lithosphere subducting on both sides of Tibet, just like the case showed in the model “indf” (Figure 5). This is consistent with the P or S receiver function results showing northward subduction of India beneath southern Tibet and southward subduction of Asia beneath northeastern (Kind et al., 2002; Kumar et al., 2006; Zhao et al., 2011. The model results show that the deformation of the northern margins of Tibet can be explained with a relatively strong Asian crust at the early stage of convergence.

As the Tibetan crust was thickened close to the present values, the lower crust becomes hotter and weaker owing to shear and radioactive heating (e.g., Chen et al., 2019). The high conductivity zones imaged by magnetotelluric studies in northern and eastern Tibet have been interpreted as widespread partial melt with the lower crust (Wei, 2001; Bai et al., 2010), which provides relevant support for its rheologically weak state. A rheologically weak lower crust can facilitates decoupling of the crust and lithosphere mantle and local downwelling of the lithosphere mantle, as demonstrated by the model “dela” (Figure 3). Therefore, I speculate that the progressive heating of the thickened crust would lead to crustal weakening and induce the delamination of the lithospheric mantle at the later stage of the Tibetan Plateau evolution. The subvertical high-velocity structures seen in southern Tibet and western Himalayan syntax (Figure 9) may exemplify localized downwelling of orogen mantle lithosphere steered by a rheologically weak lower crust. As shown in Figure 6, the delamination of the downwelling lithosphere mantle would result in a local, not broad, surface uplift. This process possibly explains a proto-Tibetan plateau first uplifted in central Tibet by 40 Ma ago (Wang et al., 2008).