The Role of Lower Crustal Rheology in Lithospheric Delamination During Orogeny

The continental lower crust is an important composition- and strength-jump layer in the lithosphere. Laboratory studies show its strength varies greatly due to a wide variety of composition. How the lower crust rheology influences the collisional orogeny remains poorly understood. Here I investigate the role of the lower crust rheology in the evolution of an orogen subject to horizontal shortening using 2D numerical models. A range of lower crustal flow laws from laboratory studies are tested to examine their effects on the styles of the accommodation of convergence. Three distinct styles are observed: 1) downwelling and subsequent delamination of orogen lithosphere mantle as a coherent slab; 2) localized thickening of orogen lithosphere; and 3) underthrusting of peripheral strong lithospheres below the orogen. Delamination occurs only if the orogen lower crust rheology is represented by the weak end-member of flow laws. The delamination is followed by partial melting of the lower crust and punctuated surface uplift confined to the orogen central region. For a moderately or extremely strong orogen lower crust, topography highs only develop on both sides of the orogen. In the Tibetan plateau, the crust has been doubly thickened but the underlying mantle lithosphere is highly heterogeneous. I suggest that the subvertical high-velocity mantle structures, as observed in southern and western Tibet, may exemplify localized delamination of the mantle lithosphere due to rheological weakening of the Tibetan lower crust.


INTRODUCTION
It is widely thought that the lithosphere mantle thickens in concert with the crust during continental collision (England and Houseman, 1986). Owing to its positive buoyancy, thickening of the continental crust by horizontal shortening often gives rise to orogens or orogenic plateaus. However, increasing evidence shows the thick orogenic crust is not necessarily underlain by an equivalently thick lithosphere mantle. For example, seismological studies indicate that the India-Asia collision results in the double normal thickness crust of Tibet (e.g., Hirn et al., 1984), but the Tibetan lithosphere mantle does not thicken proportionally. In the central Tibetan plateau, the lithosphere mantle is thinner than the normal continental lithosphere or even missing (Owens and Zandt, 1997;Nabelek et al., 2009). Similar phenomena are also observed in the Altiplano plateau. Seismic studies show low velocities and high attenuation in the shallow mantle in the Central Andes, implying the absence of the lithospheric mantle beneath much of the orogen, even though the crust thickness ranges from 50 to 75 km (Yuan et al., 2000;Beck and Zandt, 2002;Schurr et al., 2006). The discrepancy between the crustal and lithospheric thicknesses under orogens raise an important question: What is the fate of the mantle portion of the lithosphere during orogeny?
Several hypotheses have been proposed to explain the lithospheric mantle loss during orogeny (e.g., Gogus and Pysklywec, 2008). They can be categorized into two types. One FIGURE 1 | (A) Initial geotherm used for calculating the lower crust viscosity. The continental lithosphere is composed of a 40-km-thick crust (a 20-km-thick upper crust and a 20-km-thick lower crust) and a 110-km-thick mantle. The lithosphere-asthenosphere boundary (LAB) is determined by the intersection of the geotherm (black) with the mantle adiabat (red). (B) Density (dark green) and pressure (sky blue) profiles for the lithosphere based on the geotherm presented in (A). (C) Effective viscosity of the lower crust as a function of depth for a variety of flow laws assuming a strain rate of 10 −16 1/s. Diorite and plagioclase flow laws are from Ranalli (1995), granulite from Wang et al. (2012), anorthite from Rybacki and Dresen (2000), dry clinopyroxene (cpx) from Moghadam et al. (2010), and dry diabase from Mackwell et al. (1998). (D) Relative viscosity profiles compared to the diorite-derived viscosity in (C).
is convective removal through gravitational instability (e.g., Houseman et al., 1981;Conrad and Molnar, 1999;Gorczyk et al., 2012). This hypothesis requires no pre-existing structural weakness, but requires that the lithosphere mantle is denser than the underlying asthenosphere or a thickened region (sufficient perturbation with a suitable wavelength) in the lithosphere (Schott and Schmeling, 1998;Elkins-Tanton, 2007). Recent studies indicate that the continental lithosphere is largely neutrally buoyant (e.g., Lee et al., 2011), which deviates from the requirement. The other is commonly called delamination in literature, which originally indicates that the cold lithosphere mantle peels away as a coherent slice from the crust via an elongated conduit connecting the base of the crust with the asthenosphere (Bird, 1979). The delaminated cold mantle is replaced by hot asthenosphere and would cause surface uplift, elevated heat flow, reduced seismic velocities, and magmatism (Bird, 1979;Kay and Kay, 1993;Ducea and Saleeby, 1998). This mechanism requires a pre-existing weakness to initiate delamination, which is not necessarily met during collisional orogeny. Previous modeling studies on lithospheric mantle removal have highlighted the importance of viscous rheology of the mantle lithosphere (Buck and Toksoz, 1983;Conrad and Molnar, 1997;Houseman and Molnar, 1997), thermal diffusion (Conrad and Molnar, 1999), mechanical thickening (Molnar et al., 1998;Conrad, 2000), compositional density variations of the lower crust (Jull and Kelemen, 2001;Elkins-Tanton, 2005;Krystopowicz and Currie, 2013), Moho temperature (Morency and Doin, 2004), lower crust radioactive heating (Pysklywec and Beaumont, 2004), and prescribed weak zones or lithosphere perturbation root (Schott and Schmeling, 1998;Gogus and Pysklywec, 2008;Bajolet et al., 2012) in the lithosphere mantle removal. However, few studies paid attention to the lower crust rheology (e.g., Meissner and Mooney, 1998).
Laboratory results show that the continental lower crust rheology differs strongly as a function of composition, temperature, stress, and fluids (Burgmann and Dresen, 2008). There is great variety in composition and chemical and physical properties of the continental lower crust (Lee, 2014;Hacker et al., 2015). This results in large uncertainty about the lower crust rheology. Figure 1 shows the comparison among viscosities calculated from six laboratory-derived lower crust flow laws using the same geotherm and pressure curves. They encompass a wide variation of the lower crust minerals or rocks, ranging from felsic diorite to mafic diabase. It is evident that the diabase viscosity is at least three orders of magnitude greater than that of diorite at the same lower crust condition. Thus, an unresolved problem is how the uncertainty in the lower crust rheology influences crustal deformation and mantle dynamics during collisional orogeny.
In this study, I use 2D high-resolution thermal-mechanical to examine the role of the lower crust rheology in the deformation of the continental lithosphere subject to horizontal shortening. I systematically investigate variations in the lower crust rheology derived from laboratory studies to understand the conditions required for lithospheric mantle removal and its surface expressions. In contrast to some previous studies, I do not impose initial lithosphere thickness perturbation or local weaknesses.
Instead, lithospheric thickening and/or delamination develop self-consistently during model evolution. This study focuses on the variations in the laboratory-derived rheologies of the lower crust and fate of the lithosphere mantle during continental compression.

NUMERICAL MODELLING OF COUPLED LITHOSPHERIC DEFORMATION AND TOPOGRAPHY EVOLUTION
Numerical Method 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 nonuniform 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.

Model Setup
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-kmthick 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 3 ) and low viscosity (10 18 Pa s). The large density and viscosity contrast between the sticky air and rocky domain ensures small shear stresses (<10 4 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.

Rheology
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): where σ yield 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 fluidcontaining porous or fractured rocks; _ ε II (0.5_ ε ij _ ε ij ) 1/2 is the second invariant of the strain rate tensor; η plastic is the viscosity for plastic rheology. φ eff 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): Model setup illustrating the initial model configuration and thermal-mechanical boundary conditions. The red vertical bars with lateral arrows show the locations where convergence rates are imposed. White thin lines are isotherms with an interval of 300℃. The composition codes are: 0, sticky air; 1, water; 2 and 3, orogen upper crust; 4 and 5, orogen middle crust; 6, orogen lower crust; 7, craton upper crust; 8, craton lower crust; 9, orogen lithosphere mantle; 10, craton lithosphere mantle; 11, asthenosphere; 12, partially molten lower crust. (B) Lithosphere yield stress profiles for orogen lithosphere. The profiles are calculated with a constant strain rate of 10 −16 1/s and initial temperature shown in (A), both of which are variable in modeling. See Table 1 for all the rheological parameters. LC, lower crust.
where A D (material constant), E a (activation energy), V a (activation volume), n (stress exponent) are experimentally determined flow law parameters. C OH 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): 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 18 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.

MODELING RESULTS
I conducted ten numerical simulations by varying the lower crust flow law and convergence rate. Three distinct modes of model deformation have been identified: 1) lithosphere mantle delamination, 2) localized thickening, and 3) lithospheric mantle underthrusting. Table 2 shows a summary of all the simulations. 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  Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 755519 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).

Intermediately Strong Lower Crust and Localized Thickening
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.

Strong Lower Crust and Underthrusting
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.

Topographic Evolution
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 FIGURE 4 | Time evolution of the model "delc", in which the orogen lower crust rheology is represented by "plagioclase". (Ai-C) composition evolution for selected model times (colors as in Figure 2). (Aii-Cii) viscosity distribution for selected model times (colors as in Figure 3).
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.

Effects of Convergence Rate
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).

Comparisons With Previous Delamination Modeling
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 FIGURE 9 | Schematic illustration of three types of lithosphere mantle removal during collisional orogeny. (A) Conventional delamination model (Bird, 1979). (B) Convective removal via Rayleigh-Taylor instability (Houseman et al., 1981). (C) Weakened lower crust-assisted delamination.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 755519 convergence proceeds. A comprehensive investigation of the surface expression resulted from delamination will be presented elsewhere.

Implications for Tibetan Mantle Dynamics
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 Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 755519 ( 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).

CONCLUSION
In this study, I investigate the role of the lower crust rheology and convergence velocity based on a range of laboratory-derived flow laws in collisional orogeny using 2-D thermomechanical modeling method. I have drawn the following conclusions: 1) The lower crust strength and velocity are critical factors controlling how the orogen lithosphere mantle is accommodated during horizontal shortening. If the lower crust is extremely weak, the upper portion of the orogen crust decouples from the mantle lithosphere. Crustal folding and buckling are accompanied by localized downwelling of the mantle lithosphere, followed by detachment as a coherent slab. Conversely, if the lower crust rheology is represented by strong end-member of flow laws (e.g., diabase), the convergence is accommodated by underthrusting of the strong block lithospheres on both sides of the orogen. For an intermediately strong lower crust (e.g., plagioclase), concurrent thickening of the orogen crust and mantle lithosphere at the middle and sides of the orogen accommodates most of convergence, and thickened mantle lithosphere remains below the orogen. 2) The surface expression of an orogen subjected to horizontal shortening depends on the strength of the orogen itself. For a weak orogen, topography highs first develop at orogenic margins. The rise of the central region, resulted from lithosphere delamination, postdates the uplift of the margins. In contrast, topography highs only form at the orogenic margins for a strong orogen, and no punctuated surface uplift occurs in the orogen's interior. 3) When compared to the Tibetan plateau, I suggest that the subvertical high-velocity structures, as revealed by seismic tomography studies in southern Tibet and western syntax, possibly exemplify localized delaminating of the thickened mantle lithosphere owing to rheological weakening of the lower crust during the India-Asia convergence.