Front. Earth Sci. Frontiers in Earth Science Front. Earth Sci. 2296-6463 Frontiers Media S.A. 1067239 10.3389/feart.2023.1067239 Earth Science Original Research Patched removal of the mantle lithosphere under orogens: A systematic numerical study Lei et al. 10.3389/feart.2023.1067239 Lei Tian 1 Liao Jie 1 2 3 * Li Zhong-Hai 4 Liu Mian 5 Li Lun 1 1 School of Earth Sciences and Engineering, Sun Yat-sen University, Guangzhou, China 2 Southern Marine Science and Engineering Guangdong Laboratory, Zhuhai, China 3 Guangdong Provincial Key Lab of Geodynamics and Geohazards, Guangzhou, China 4 Key Laboratory of Computational Geodynamics, College of Earth and Planetary Sciences, University of Chinese Academy of Sciences, Beijing, China 5 Department of Geological Sciences, University of Missouri, Columbia, MO, United States

Edited by: Weijia Sun, Institute of Geology and Geophysics (CAS), China

Reviewed by: Liang Qiu, China University of Geosciences, China

Eh Tan, Academia Sinica, Taiwan

Shaohua Yang, Chinese Academy of Geological Sciences (CAGS), China

Yajin Pang, China Earthquake Administration, China

 *Correspondence: Jie Liao, Liaojie5@mail.sysu.edu.cn

This article was submitted to Solid Earth Geophysics, a section of the journal Frontiers in Earth Science

 13 03 2023 2023 11 1067239 11 10 2022 31 01 2023 Copyright © 2023 Lei, Liao, Li, Liu and Li. 2023 Lei, Liao, Li, Liu and Li

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Delamination or convective thinning could cause large-scale and complete removal of the mantle lithosphere under orogens. However, geological and geophysical observations suggest that patched removal of the mantle lithosphere has occurred in some orogens, such as the northeastern Tibetan Plateau, the central Tianshan, and the central Andes. Dislocation-creep-induced strain localization cannot promote effective removal of the mantle lithosphere to the Moho on a small-scale. Recent rheological studies propose that dislocation-accommodated grain boundary sliding (DisGBS) may dominate upper mantle deformation. DisGBS could make the lower lithospheric mantle rheologically weaker than dry olivine. With 2-D high-resolution thermo-mechanical modeling, we systematically investigated the conditions for the initiation of small-scale lithospheric thinning under orogens and explored the minimum range of removal of the mantle lithosphere. The numerical results indicate that classic convective drip cannot effectively thin the mantle lithosphere to the Moho on a small-scale. In contrast, small-scale thinning can be induced by lithospheric heterogeneity with DisGBS and plasticity. The rheological heterogeneity can be verified by magmatism and metasomatism under the central Andes and orogens between terranes under the northeastern Tibetan Plateau or in Tianshan.

 small-scale thinning rheological heterogeneity dislocation-accommodated grain boundary sliding numerical modeling orogen 1 Introduction

Lithospheric mantle removal under orogens (e.g., the Himalayan–Tibetan Plateau, Central Anatolia, and Colorado) is often explained by delamination or convective thinning (Bird 1978; Bird 1979; Owens & Zandt 1997; Molnar et al., 1998; Gogus et al., 2017). Both delamination and convective thinning are driven by the negative buoyancy in the lithospheric mantle (and sometimes the eclogitized lower crust), even though the growth period of instability and the rheological control of them are different (Levander et al., 2011; Beall et al., 2017; Lei et al., 2019). The long wavelength of downwelling thickened lithosphere leads to mantle lithosphere thinning, and the range of the initial wavelength is commonly larger than the thickness of the mantle lithosphere (Conrad & Molnar 1997; Houseman & Molnar 1997). A laterally homogeneous rheological structure of the lithosphere with dislocation creep or diffusion creep is usually assumed in previous models of delamination or convective thinning, which predict large-scale peeling off of the mantle lithosphere (Bird, 1978; Bird 1979; Schott & Schmeling 1998; Morency & Doin 2004; Göğüş & Pysklywec 2008a; Göğüş & Pysklywec 2008b; Bajolet et al., 2012; Wang & Currie 2015; Beall et al., 2017; Lei et al., 2019). However, patched removal of the mantle lithosphere, with insignificant mafic magmatism, has been observed in the interior of collisional plates and the margin of subducting plates, such as the northeastern Tibetan Plateau, the Tianshan, and the central Andes (Beck & Zandt 2002; Schurr et al., 2006; Ducea 2011; Ducea et al., 2013; Li et al., 2022; Zhang et al., 2022). The width of detachment approaches the layer thickness of the mantle lithosphere under these orogens (∼100–200 km). Such patched removal of the mantle lithosphere under orogens cannot be readily explained by existing models of delamination or convective thinning because the convective dripping resulting from the short wavelength of perturbation tends to occur near the bottom of the mantle lithosphere and cannot remove the entire mantle lithosphere to the Moho (Conrad & Molnar 1997; Houseman & Molnar 1997; Currie et al., 2008; Lorinczi and Houseman, 2009; Ducea 2011; Beall et al., 2017; Liao et al., 2017).

Rheological weakening is the key factor for small-scale mantle lithosphere removal under orogens, that is, removal from the bottom of the mantle lithosphere to the Moho. However, dislocation creep and diffusion creep overestimate the upper mantle’s viscosity (van Hunen et al., 2005; Billen, 2008), which may not effectively promote the patched removal of the mantle lithosphere (Currie et al., 2008; Liao et al., 2017). Recent rheological studies indicate that dislocation-accommodated grain boundary sliding (DisGBS), which differs from dislocation creep or diffusion creep, may dominate the deformation of olivine in the Earth’s upper mantle (Tomohiro and Takaaki, 2015). The estimated viscosity of the upper mantle controlled by DisGBS is independent of depth and is lower than that of dislocation creep, indicating that the local mantle lithosphere may experience significant deformation via strain localization or hydration activity.

The mantle lithosphere in the northeastern Tibetan Plateau, the Tianshan, and the Central Andes manifests obvious rheological heterogeneities. Geological and geophysical observations indicate that the northeastern Tibetan Plateau consists of multiple terranes with contrasting lithologies, thermal states, and rheological structures. Patches of seismically fast and slow anomalies have been observed beneath the northeastern Tibetan Plateau (Deng et al., 2018). For example, the mantle lithosphere in the Songpan–Ganzi and southern Kunlun–Qaidam terranes has probably been removed, with traces of delaminated pieces. The entire mantle lithosphere in North Qilian has been removed, together with the base of the crust, leading to the sudden decrease of the Moho depth (Deng et al., 2018). This interpretation is consistent with the rapid uplift of the Qilian block at the Miocene-Quaternary from recent magnetostratigraphy and tectonosedimentology (Fang et al., 2013). Young volcanic rocks are exposed (<6 Myr) in the western North Qilian Mountains (Xia et al., 2011). Therefore, lithospheric heterogeneities may play a key role in the shortening of the northeastern Tibetan Plateau (Wang et al., 2008; Yin and Dang, 2008; Zhang and Wang, 2014; Deng et al., 2018). Compared with the Tarim Craton and Kazakh lithosphere (Li et al., 2022), the crust and mantle lithosphere of the sandwiched central Tianshan are assumed to be much rheologically weaker. The intracontinental deformation of the central Tianshan is controlled by rheological heterogeneity (Huangfu et al., 2021; Li et al., 2022). In contrast, the mantle lithosphere under the central Andes does not involve multiple terranes and sutures, thus differing from the northeastern TP and central Tianshan. The tectonic history of subduction has induced the migration of aqueous fluids from the subducting Nazca plate and may also lead to the weakening of the localized mantle lithosphere under the central Andes, which is related to variations in magmatism and metasomatism (Kay & Kay 1993; Beck & Zandt 2002; Jing et al., 2020; Contreras-Reyes and Diaz, 2021). Geochemical constraints show that the peridotite melting of mantle-derived magmatism may have occurred within a volcanic field over short time scales (1–5 Myr) in the Altiplano–Puna Plateau. The pattern of melting is consistent with the convective removal of the small-scale mantle lithosphere under the Altiplano–Puna Plateau (Ducea et al., 2013). Could such rheological heterogeneity in the mantle lithosphere account for the observed small-scale lithosphere thinning? Furthermore, with lithospheric heterogeneities, how does DisGBS regulate the removal of the whole lithospheric mantle? Can Rayleigh–Taylor instability developing in the perturbations of short wavelengths remove the mantle lithosphere with DisGBS to the Moho in the homogenous rheological model? In order to address these questions, we constructed a series of 2-D high-resolution thermo-mechanical models with homogeneous mantle lithosphere, including uniformly distributed perturbations, localized weak mantle lithosphere, and multiple terranes/blocks to systematically investigate the dynamics and rheological constraints of the small-scale removal of the mantle lithosphere under orogens.

 2 Numerical model 2.1 Governing equations

We simulated lithospheric thinning in collisional orogens by numerically solving the governing equations of mass, momentum, and energy conservation in two-dimensional (2-D) finite difference models with a marker-in-cell technique (Gerya & Yuen 2003). The governing equations are as follows: ∂ v i ∂ x i = 0 , ∂ σ i j ′ ∂ x j − ∂ P ∂ x i + ρ g i = 0 , ρ c p ∂ T ∂ t + v i ∂ T ∂ x i = ∂ ∂ x i k ∂ T ∂ x i + H r + H a + H s , where v is velocity, σ ′ is the deviatoric stress tensor, P is pressure, ρ is density, g is gravity acceleration, c p is heat capacity, T is temperature, k is thermal conductivity, H r is the radioactive heating rate, H a is the adiabatic heating rate ( H a = T α D P / D T , α is thermal expansivity, and H s is the shear heating rate ( H s = σ i j ′ ε ˙ i j , ε ˙ strain rate tensor).

For the density of a specific rock type, ρ depends on pressure (P) and temperature (T): ρ P , T = ρ 0 1 − α T − T 0 1 + β P − P 0 , where ρ 0 is the reference density under the conditions of P 0 = 0.1 MPa and T 0 = 298 K (i.e., pressure and temperature at Earth’s surface); α and β are the thermal expansion coefficient and the compressibility coefficient.

 2.2 Viscoplastic rheology 1 Viscous rheology

η d u c t i l e = ε ˙ I I 1 − n n A D − 1 G p f H 2 O − r 1 n ⁡ exp E + P V n R T , where ε ˙ I I = ε ˙ i j ε ˙ i j 1 / 2 is the second invariant of the strain rate tensor; A D , the pre-exponential viscous factor; E, the activation energy; V, the activation volume; G, the grain size; f H 2 O , the water fugacity; and n (creep exponent) are experimentally determined flow law parameters (Table 1). The parameter R is the gas constant, p is the grain size exponent, and r is the water fugacity exponent.

Viscous flow laws used in the numerical experiments.

Id symbol Flow law E ( k J m o l − 1 ) V ( J M P a − 1 m o l − 1 ) n log 10 A D m p MPa − n − r s − 1
RA Disl (wet quartzite) 154 8 2.3 − 3.49
RB Disl (plagioclase A n 75 ) 238 8 3.2 − 3.48
RC Disl (dry olivine) 532 8 3.5 4.39
RD DisGBS (wet olivine) 423 17.6 3.0 − 4.89

Table based on Ranalli (1995), which compiled original data from Kirby (1983), Kirby and Kronenberg (1987), and Ranalli and Murphy (1987). Data on dislocation-accommodated grain boundary sliding (DisGBS) are from Tomohiro and Takaaki (2015).

Dislocation-accommodated grain boundary sliding is used in our numerical model. Recent studies have proposed that DisGBS may dominate upper mantle deformation (Hansen et al., 2011; Tomohiro and Takaaki, 2015). DisGBS is more effective at a larger strain rate and a smaller grain size, which is consistent with the upper mantle of orogens under convergence. We have applied the same rheological profile for DisGBS to the weak mantle lithosphere and asthenosphere, indicating that the mantle lithosphere of weak terrane is weakened from potential fluid/melt during previous oceanic subduction or terrane accretion. Consequently, the relatively low viscosity based on wet olivine rheology and low plastic strength are applied for the lithospheric mantle of the weak terrane (Tables 1, 2). We used the effective viscosity of “wet quartzite” for both the lower and upper crust of weak terrane (Tables 1, 2). In contrast, the flow law of “wet quartzite” is used for the upper continental crust and “Plagioclase A n 75 ” for the lower continental crust in the strong terrane (Tables 1, 2). “Dry olivine” for the lithospheric mantle and the asthenosphere (Tables 1, 2) is generally used for the strong continental lithosphere of the strong terrane, where the fluid/melt weakening effects are neglected; thus, a high plastic effective friction coefficient ( sin φ d r y λ = 0.6 ∼ 0.3 ) is used for the dry and strong lithospheric mantle.

Material properties used in the numerical experiments. a

Section Material b ρ 0 k g   m − 3 k c W m − 1 K − 1 Viscous flow law d Plastic e C 0 M P a Plastic e sin ϕ d r y λ
All plates Asthenosphere (7) 3300 K 3 RD 10–1 0.6–0.3
Strong terrane Upper crust (1) 2700 K 1 RA 10–1 0.15–0.075
Lower crust (2) 3000 K 2 RB 10–1 0.15–0.075
Mantle lithosphere (5) 3300 K 3 RC 10–1 0.6–0.3
Weak terrane Upper crust (3) 2700 K 1 RA 10–1 0.15–0.075
Lower crust (4) 3000 K 2 RA 10–1 0.15–0.075
Mantle lithosphere (6) 3300 K 3 RD 10–1 Varied
References f 1,2 3 4,5 6 6

Isobaric heat capacity C p = 1000 J   k g − 1 K − 1 , thermal expansion coefficient α = 3 × 10 − 5 K − 1 , and compressibility coefficient β = 1 × 10 − 5 M P a − 1 are used for all rock types.

Numbers of materials correspond to Figure 1.

K 1 = 0.64 + 807 T k + 77 exp 0.00004 P M P a ; K 2 = 1.18 + 474 T k + 77 exp ⁡ ⁡ 0.00004 P M P a ; K 3 = 0.73 + 1293 T k + 77 exp 0.00004 P M P a

Parameters of viscous flow laws are shown in Table 1.

strain weakening effect is applied for plastic rheology, in which both cohesion (C 0) and the effective friction coefficient sin φ d r y λ decrease with strain increase. The cut-off values shown in the table correspond to the strain ‘0–1’.

References 1–6 are Turcotte and Schubert (2002), Bittner and Schmeling (1995), Clauser and Huenges (1995), Ranalli (1995), Tomohiro and Takaaki (2015), and Vogt et al. (2012), respectively.

 2 Drucker–Prager plasticity

The extended Drucker–Prager yield criterion (e.g., Ranalli, 1995) is adopted in our model to simulate the viscoplastic behavior of the lithosphere: η p l a s t i c = σ y i e l d 2 ε ˙ I I , σ y i e l d = C 0 + P ⁡ sin ϕ d r y λ , where η p l a s t i c is the viscosity of the Drucker–Prager plasticity, σ y i e l d is the yield stress, P is the dynamic pressure, C 0 is the residual rock strength at p=0, and φ d r y is the internal frictional angle of dry rocks. λ is the pore fluid/melt coefficient that controls the brittle strength of fluid/melt containing porous or fractured media: λ = 1 − P f l u i d / m e l t P . Because the fluid/melt parameter P f l u i d / m e l t is not directly calculated in the model, we use a range of “ sin φ d r y λ ” values to represent the effective friction coefficient, based on previous systematic investigations (e.g., Gerya & Meilick 2011; Vogt et al., 2012; Li et al., 2016). The strain weakening effect is included in the plastic rheology, in which both the cohesion C 0 and effective friction coefficient sin φ d r y λ decrease with increased strain, as shown in Table 2.

The minimum value of the viscous or plastic viscosity defines the effective viscosity in the model (Ranalli 1995): η e f f = min η d u c t i l e , η p l a s t i c .

2.3 Initial model configuration and boundary conditions

Large-scale models ( 4000 × 400   k m ) were built to study the dynamics of the small-scale removal of the mantle lithosphere. Using a non-uniform rectangular numerical grid, the collision zone is represented by 1 × 1 km high-resolution grids, while 5 × 1 km grids are used for the rest of the model domain. More than 15 million active Lagrangian markers are used to trace and mark internal lithological boundaries, material properties, and temperature. Although the model length is 4,000 km, the deformation localizes in a relatively narrow (ca. 250 km and 1,000 km) region of interest. The effects of the Earth’s curvature are, therefore, neglected in this simplified Cartesian model.

The lithosphere includes a 20-km thick upper crust, a 15-km thick lower crust, and a 105-km thick lithospheric mantle, underlain by the asthenosphere. Some orogens, such as the central Andes, northeastern Tibetan Plateau, and Tianshan, could have experienced a series of subduction and lithospheric detachments (Ducea et al., 2013; Deng et al., 2018; Jing et al., 2020; Huangfu et al., 2021; Li et al., 2022) which could have reduced the viscosity of the mantle (Lei et al., 2019). Moreover, the lithospheric mantle under continental orogens that are adjacent to a subduction zone, such as South America and the Gibraltar Arc, may be modified by subducting plates, leading to the development of secondary downwelling in the continental interior (Levander and Berzada, 2014). Consequently, three groups of models were constructed to simulate the small-scale removing mantle under orogens, and a series of initial perturbations are integrated into the models of Type I∼Type III. In the model of Type I, the lithosphere is set to be weak and homogeneous, and a series of initial perturbations with a constant wavelength (u) and perturbation ( Z 0 ) are designed at the bottom of the mantle lithosphere (Figure 1A). In the model with heterogeneous rheology (Type II and Type III), the lithosphere is divided into strong and weak terranes, with the initial perturbation imposed on the bottom of the mantle lithosphere in the weak terrane (Figures 1B, C). The properties of various rock types are summarized in Tables 1 and 2. We assumed the same reference density ρ 0 for the mantle lithosphere and asthenosphere (Table 2). Thus, the gravitational instability of the mantle lithosphere arises from the thickening of the weak lithospheric mantle, which becomes denser than the asthenosphere because of the lower temperature. The initial thermal structure of the lithosphere (the white lines in Figure 1A) is laterally uniform with a linear gradient from 0 ℃ at the surface to 1350 ℃ at the bottom of the lithosphere (Turcotte & Schubert 2002). The initial adiabatic thermal gradient in the asthenosphere is 0.5 ℃ • k m − 1 .

Reference model and boundary conditions. (A) Composition field of the homogenous model with initial wavelength (u=150 km) and perturbation ( Z 0 = 10 k m ). Enlargement (2000×400 km) of the numerical box (4,000×400 km) is shown. White lines are isotherms in °C with intervals of 200 °C. Initial convergence rates ( V x = ± 4.5 + 4.5 c m / y r ) are imposed on the small internal domains in the left and right parts, respectively (yellow arrow). (B) Composition field of heterogeneity model with a single weak block. (C) Composition field of the heterogeneity model with two weak blocks. Colors indicate the different rock types; 1, 2: continental upper and lower crust of the strong terrane; 3, 4: continental upper and lower crust of the weak terrane; 5: lithospheric mantle of the strong terrane; 6: lithospheric mantle of the weak terrane; 7: asthenosphere.

The velocity boundary conditions in the model of heterogeneity rheology include free slip for the left, right, and top boundaries and a permeable boundary for the lower boundary (Burg & Gerya 2005; Li et al., 2016). This infinity-like external free-slip condition along the lower boundary implies a free slip condition to be satisfied at about 100 km below the base of the model. The external free slip allows global conservation of mass in the computational domain and is implemented by using the following limitation for velocity components at the lower boundary: ∂ v x / ∂ y = 0 and ∂ v y / ∂ y = − v y / Δ y e x t e r n a l , where Δ y e x t e r n a l is the vertical distance from the lower boundary to the external boundary where free slip ( ∂ v x / ∂ y = 0 , v y = 0 ) is satisfied. The lithosphere is pushed from both sides with a constant convergence velocity ( V x ) imposed on the left and the right parts of the models (Figure 1A).

The thermal boundary conditions contain a fixed temperature ( 0   ℃ ) at the upper boundary and zero horizontal heat flux across the vertical boundaries. For the lower thermal boundary, a constant temperature condition is imposed at a great depth (1,000 km) below the bottom of the model to allow both temperature and vertical heat flux to vary along the permeable lower boundary of the model domain and to be adjusted dynamically during the model’s evolution (Li et al., 2016).

 3 Model results 3.1 Type I model with uniformly distributed perturbations

In the case of the homogeneous mantle lithosphere, models of classical delamination predict a large-scale peeling of the mantle lithosphere (Bird 1978, 1979; Lei et al., 2019), which may not adequately explain the small-scale lithospheric mantle removal under some orogens. Rayleigh–Taylor instability could occur in small-scale drips when viscosity is relatively low (Beall et al., 2017; Lei et al., 2019). However, will short wavelength drips remove the mantle lithosphere from its base to the Moho? In order to answer this question, we have performed a series of 2-D numerical experiments to systematically investigate the effects of DisGBS and plasticity on the removal of the mantle lithosphere with uniformly distributed perturbations.

3.1.1 Effect of perturbation wavelength on model evolution

In this case, we used dislocation-accommodated grain boundary sliding (DisGBS) and plasticity ( λ = 0.001 ) to study the effects of perturbation wavelength on convective thinning. Figure 2 shows the model results of different perturbation wavelengths from 150 km to 800 km with constant amplitude ( Z 0 = 20 k m ) at 3.4 Myr, where convergent velocity ν x = ± 4.5 + 4.5 c m / y r . When the initial perturbation wavelength is small ( u = 150 k m − 300 k m ), the drips cannot effectively remove the mantle lithosphere (Figures 2A, B). In contrast, with a larger wavelength initial perturbation ( u = 500 k m − 800 k m ), the lower mantle lithosphere can be easily removed by the drips (Figures 2C, D).

Effects of various wavelengths on the homogenous model. The model integrates the dislocation-accommodated grain boundary sliding (n=3, r=1.25, p=1), with a constant initial perturbation of Z 0 = 20 k m and a convergence velocity of v x = ± 4.5 + 4.5 c m / y r . (A–D) Constant initial perturbation ( Z 0 = 20 k m ) and different wavelengths (u=150 km–1000 km). White lines are isotherms with 200 ℃ intervals. (E–H) Initial model at t = 0 Myr corresponding to the red boxes in (A–D).

 3.1.2 Effect of perturbation amplitude on model evolution

In addition to various perturbation wavelengths, different perturbation amplitudes can also play an important role in convective thinning. Figure 3 shows the model results of different amplitudes from 1 km to 30 km with a constant of wavelength ( u = 500 k m ) at 3.4 Myr. With a small perturbation amplitude (Figure 3A, Z 0 = 1 k m ), no convective dripping occurs. When the initial perturbation amplitude Z 0 = 10 k m , the drips developed on the bottom of the mantle lithosphere; however, this process cannot effectively remove the lithospheric mantle to the Moho (Figure 3B). With a large amplitude ( Z 0 = 20 k m ), obvious convective dripping developed on the bottom of the mantle lithosphere, leading to lower mantle lithosphere removal from ∼150 km to ∼100 km (Figure 3C). Furthermore, in the case of perturbation amplitude Z 0 = 30 k m , the drips effectively induced a removal of the mantle lithosphere from ∼150 km to ∼50 km (Figure 3D).

Effects of various amplitudes on the homogenous model. The model integrates the dislocation-accommodated grain boundary sliding (n=3, r=1.25, p=1) with a constant initial wavelength of u = 500 k m and a convergence velocity of v x = ± 4.5 + 4.5 c m / y r . (A–D) Constant initial wavelength (u = 500 km) and different perturbation ranges ( Z 0 = 1 k m − 30 k m ). White lines are isotherms with 200 ℃ intervals. (E–H) Initial model at t = 0 Myr corresponding to the red boxes in (A–D).

 3.2 Type II model with a weak block

Natural orogens often consist of numerous terranes with contrasting lithologies, thermal states, and rheological structures, or have experienced localized weakening in the mantle lithosphere through metasomatism induced by subduction or collision (Ducea et al., 2013; Deng et al., 2018; Jing et al., 2020). Could such rheological heterogeneity explain the observed small-scale lithospheric thinning? In order to solve these problems, we performed a series of experiments with localized weak lithospheric mantle to systematically investigate the effects of DisGBS and plasticity on the removal of the mantle lithosphere with a weak block.

3.2.1 Effect of weak crust and mantle

In the first set of models, we studied the effects of the weak crust and mantle lithosphere on the removal of the mantle lithosphere. Figure 4 shows the model results with DisGBS and significant plastic yielding ( λ = 0.001 ). At the beginning of the experiment, the localized weak mantle lithosphere is thickened by small scale drips (Figure 4B). With further convergence, the obvious fragmentary or small-scale thinning of the mantle lithosphere occurs in the section of weak lithosphere with weak crust (Figures 4C, D). Figures 4E–H show the corresponding effective viscosity and the second invariant of the strain rate. It shows that the strain rate localization goes through the entire mantle lithosphere by conjugate slip and has reduced viscosity in the same way.

Model evolution with heterogeneity. Dislocation-accommodated grain boundary sliding (n=3, r=1.25, and p=1) and plastic yielding ( λ = 0.001 ) are used for the weak mantle. (A–D) Results of the composition field (L = 300 km). (E,F) Effective viscosity of the model. (G,H) Second invariants of the strain rate in the model. Colors of the composition field indicate the different rock types (Figure 1). White lines are isotherms with 200°C intervals. Time (Myr) of evolution is shown in each panel.

The model in Figure 4 assumes a constant width of L = 300 km for the weak terrane. Furthermore, we tested the model with different widths of weak terrane (Figure 5). When the width of the weak block is small (L = 100 km), converging continents lead to pure shear thickening under weak terrane (Figure 5A). In contrast, the obvious removal of lithospheric mantle under weak terrane, with thickened crust, occurs at a width of L = 200 km (Figure 5B). The decreasing of plastic equivalent viscosity with high strain rates explains the quick development of instability in the lithospheric mantle on both sides of the weak terrane. This leads to the upwelling of the asthenosphere to the bottom of the crust to trigger delamination under weak terrane. With the increasing width of weak terrane, convective dripping can effectively remove the lithospheric mantle to the Moho without apparently thickened crust (Figures 5C–E).

Range of effective thinning in the weak terrane. Dislocation-accommodated grain boundary sliding (n=3, r=1.25, and p=1) and plastic yielding ( λ = 0.001 ) are used for the weak mantle. (A–E) Composition field with different widths L = 100 km–500 km in the weak mantle lithosphere. Time (Myr) of the models (A–E) is, respectively, 0.6 Myr, 1.6 Myr, 2.1 Myr, 2.6 Myr, and 2.3 Myr. Colors of the composition field indicate the different rock types (Figure 1). White lines are isotherms with 200°C intervals.

 3.2.2 Effect of localized weak mantle lithosphere

In this set of models, the localized weak mantle lithosphere, with a strong crust, is added into model setup to simulate the effects of the local weak mantle lithosphere on the small-scale thinning. Figure 6 shows the evolution results of the composition field with localized weak mantle, where the plastic yield i s   0.001 . When the width of the weak mantle lithosphere is small (L = 100–200 km), the drips could not cause effective removing under the weak mantle lithosphere (Figures 6A, B). In comparison with models of weak terrane (Figure 5), the removal of lithospheric mantle occurs in the critical width of the weak region (L = 300 km), which is larger than the weak terrane model (Figure 5). With the increasing width of weak regions, convective dripping can effectively remove the lithospheric mantle to the Moho without apparently thickened crust (Figures 6D, E).

Range of effective thinning in the weak mantle lithosphere. Dislocation-accommodated grain boundary sliding (n=3, r=1.25, and p=1) and plastic yielding ( λ = 0.001 ) are used for the weak mantle. (A–E) Composition field with different widths L = 100 km–500 km in the weak mantle lithosphere. Time (Myr) of the models (A–E) is, respectively, 0.7 Myr, 2.0 Myr, 1.4 Myr, 1.7 Myr, and 2.1 Myr. Colors of the composition field indicate the different rock types (Figure 1). White lines are isotherms with 200°C intervals.

 3.3 Type III model with multiple rigid and weak blocks

We have further considered the effects of multiple weak terranes (Deng et al., 2018). Figure 7 shows the corresponding composition, effective viscosity, and the second invariant of the strain rate field. These results illustrate the feedback between strain rates and effective viscosity when DisGBS and plastic yielding are considered. Developed drips occur in the weak terrane and promote the mantle lithosphere peeling off to the crust (Figures 7A–D). As with the model results of a single section of weak terrane (Figure 4), strain weakening reduces both viscous and plastic viscosity, causing localized drips to develop in the bottom of the mantle lithosphere and the further removal of crust (Figures 7E–J). Like the model results of a single weak block (Figure 5), the effective thinning of the mantle lithosphere occurs in the weak terrane, when its width falls into the range of 200 km to 400 km (Figures 8B–D), except when L = 100 km (Figure 8A).

Model evolution with the heterogeneity of multiple weak terranes. Dislocation-accommodated grain boundary sliding (n=3, r=1.25, and p=1) and plastic yielding ( λ = 0.001 ) are used for the weak mantle. (A–D) Results of the composition field (L = 200 km). (E–J) Effective viscosity and the second invariants of the strain rate of the model.

Model evolution with the heterogeneity of multiple weak terranes with different widths. Dislocation-accommodated grain boundary sliding (n=3, r=1.25, and p=1) and plastic yielding ( λ = 0.001 ) are used for the weak mantle. (A–D) Composition field with different widths L = 100 km–400 km in the weak mantle lithosphere. Time (Myr) of the model (A–D) is 1.6 Myr, 2.8 Myr, 2.7 Myr, and 3.2 Myr, respectively. Colors of the composition field indicate the different rock types (Figure 1). White lines are isotherms with 200°C intervals.

 4 Discussion 4.1 Comparisons of the three types of models

Our study indicates that small-scale mantle lithosphere removal can develop under orogens, depending on the rheology of the mantle lithosphere. Here, we compare the effects of the rheological structure on the patched or small-scale removal of mantle lithosphere.

4.1.1 Homogeneous model (type I)

Given that dislocation-accommodated grain boundary sliding and plasticity is integrated into the mantle lithosphere of weak terranes, convective removing could occur with dripping; however, in that case it tends to occur in the model with a large perturbation wavelength and amplitude (Figures 2, 3). To quantify the effect of the initial perturbation wavelength and amplitude on the magnitude of the removal of the mantle lithosphere, we conducted an “available buoyancy” scaling analysis. The degree of convective instability depends on the thickness of a potentially unstable layer, which can be caused by the mechanical thickening of the layer. In this study, based on Conrad and Molnar (1999), the gravitational instability of the mantle lithosphere with thickness h can be represented by a Rayleigh number, R a n : R a n = ρ m g α T 0 n B m n λ 2 h Z 0 n − 1 κ F n , where n = 3 is the creep exponent of the accommodated grain boundary sliding (Table 1), T 0 is the temperature difference across the layer, g is the gravitational acceleration, α is the thermal expansion coefficient, κ is the thermal diffusivity, ρ m is the density of mantle lithosphere, B m is the rheological strength parameter of the mantle lithosphere, h is the layer thickness of the mantle lithosphere, λ is the initial wavelength, and Z 0 is the initial perturbation amplitude. This number is similar to the standard Rayleigh number with the additional factor F n , which is the “available buoyancy” of the layer. F n is calculated by the integral of the negative buoyancy divided by viscosity (Conrad & Molnar 1999). This integrates the variation of density, viscosity, and temperature with depth in an unstable layer; the dimensionless growth rate should depend only on the wavelength of the initial perturbation. Numerical experiments show that the layer will be unstable when R a 3 is greater than 100 (Conrad & Molnar 1999; Conrad 2000).

In our models, the mantle lithosphere density profile depends on pressure and temperature, while the effective viscosity profile is given by both plasticity and the background strain rate associated with convergence. Figure 9 shows the calculated R a 3 as a function of layer thickness at different initial perturbation wavelengths and amplitudes, respectively. At the small amplitude ( Z 0 = 1   k m ), the analysis indicates that the whole mantle lithosphere remains stable, with R a 3 < 100 for all thicknesses (Figure 9A). This agrees with the model result that indicates that the mantle lithosphere is stable and thickened during convergence when the small amplitude is applied (Figure 3A). When the initial perturbation amplitude Z 0 falls in the range of 10 km to 30 km, the calculations show that the large wavelength of the initial perturbation becomes easier to super-exponentially grow and further remove part of the mantle lithosphere (Figures 9B–D). The thickness h c r i t represents the minimum amount of mantle lithosphere that may be removed by gravitational instability, which indicates that the critical R a 3 of 100 is obtained for layer thickness. In the models with the larger amplitude ( Z 0 = 20   k m ), the calculated critical R a 3 of 100 is obtained for the layer thickness of ∼45 km, ∼18 km, ∼8 km, and ∼5 km, corresponding to perturbation wavelengths of 300 km, 500 km, 800 km, and 1,000 km (Figure 9C). Because the initial amplitude ( Z 0 = 20   k m ) is larger than the critical thickness h c r i t , corresponding to a wavelength range of 500 km–1000 km, the gravitational instability will develop at the base of mantle lithosphere. The results agree with the models, which indicate that the mantle lithosphere is unstable with dripping at a large wavelength of the initial perturbation (Figures 2C, D). For the same reason, in the models with constant wavelength (u=500 km), the calculated critical R a 3 of 100 is obtained for a layer thickness of ∞ , ∼65 km, ∼18 km, and ∼10 km, corresponding to the respective perturbation amplitudes of 1 km, 10 km, 20 km, and 30 km (Figures 9A–D), where ∞ indicates that the mantle lithosphere is stable during convergence (Figure 3A). When the initial perturbation amplitude is 10 km, the minimum amount of removal of the mantle lithosphere h r (∼65 km) is larger than the initial amplitude Z 0 (∼10 km) with a 500 km initial wavelength (Figure 9B); this indicates that the evolution of the model with a 10 km initial amplitude will have difficulty reaching super-exponential growth and the removal of the mantle lithosphere (Figure 3B). In contrast, with the initial amplitude Z 0 = 20   k m , super-exponential growth develops at the base of the mantle lithosphere (Figure 3C) because the initial amplitude ( Z 0 = 20   k m ) is larger than the critical thickness h c r i t (∼18 km), corresponding to a wavelength of 500 km (Figure 9C).

(A) Predicted stability ( R a 3 ) as a function of layer thickness (measured from the base of the lithosphere) for the constant initial perturbation Δ L 0 = 1 k m with variations in the initial wavelength. The critical layer thickness ( h c r i t ) occurs at R a 3 = 100 , indicating the minimum lithospheric mantle thickness that may be removed by Rayleigh–Taylor instability. (B) Predicted stability ( R a 3 ) as a function of layer thickness for the constant initial perturbation Δ L 0 = 10 k m with variations in the initial wavelength. (C) Predicted stability ( R a 3 ) as a function of layer thickness for the constant initial perturbation Δ L 0 = 20 k m with variations in the initial wavelength. (D) Predicted stability ( R a 3 ) as a function of layer thickness for the constant initial perturbation Δ L 0 = 30 k m with variations in the initial wavelength. Dislocation-accommodated grain boundary sliding (Table 1) and plastic yielding ( λ = 0.001 ) are used for the semi-analytical model (Eq. 6).

The thickness h c r i t only estimates the minimum amount of the layer of removal by instability. Larger thicknesses are also potentially unstable ( R a 3 >100 in Figures 9B, C). The rapid growth triggered by larger instability also depends on the convergence and weakening of strain rate localization. Thickening of the mantle lithosphere by horizontal shortening can lead to dripping to remove it in several obvious ways, which include the lithosphere weakening with increasing strain rate—as is expected for the mantle lithosphere with non-Newtonian viscosity (Molnar et al., 1998; Lei et al., 2019). Moreover, the effect of plasticity promotes dripping to remove the larger quantity of the mantle lithosphere than the model without plasticity. Figure 10 shows the amount of removed mantle in the models with an initial wavelength of 800 km and amplitude of 20 km, as well as a convergence velocity of 4.5 cm/yr. First, the models show that super-exponential growth occurs with a thickening of ∼30 km, which is calculated from the depth of the initial base of the thickened layer ( 1350 ℃ ) to the depth at which the initial thickening develops super-exponential growth (Figures 10A, B). Second, the model with DisGBS (Figure 10A) shows that the thickness of the removed layer ( R a 3 > 100 ) is ∼30 km. In comparison, the larger amount of the removed layer (∼110 km) of R a 3 > 100 occurs in the model with plasticity (Figure 10B), which indicates that plasticity promotes instability to remove the larger amount of potentially unstable mantle lithosphere ( R a 3 > 100 ).

Evolution of the average depth of the 1350 ℃ isotherm for the model of Type I with the constant initial perturbation Z 0 = 20 k m and initial wavelength u = 800 km. (A) Dislocation-accommodated grain boundary sliding (Table 1) is integrated into the model; (B) dislocation-accommodated grain boundary sliding (Table 1) and plastic yielding ( λ = 0.001 ) are integrated into the model.

 4.1.2 Heterogeneity model (type II and type III): Strain rate localization

Classical delamination and convective thinning are used to account for the removal of the mantle lithosphere with a homogeneous rheological structure, driven by dripping behaviors (Beall et al., 2017; Lei et al., 2019). Delamination would trigger large-scale peeling-off of the mantle lithosphere, which cannot explain the patched removal of the lithosphere under the Central Andes and the northeastern Tibetan Plateau (Ducea 2011; Lei et al., 2019). The feedback between the strain rate localization and reduced viscosity, based on dislocation creep, is concentrated at the base of the homogeneous lithosphere.

For orogens related to ocean-to-continent subduction, strong hydration processes can lead to local weakening of the mantle lithosphere and cause rheological heterogeneity. The feedback between the strain rate and effective viscosity can reduce the rheological strength of weak mantle lithosphere, causing patched removal of the mantle lithosphere under orogens (Figure 6). On the other hand, the model with multiple weak terranes and weak crust displays obviously patched removal of the mantle lithosphere, indicating a larger parameter range for lithospheric removal than the model of a local weakening mantle lithosphere (Figure 4 and Figure 5). However, the strain rate localization caused by dislocation creep (Ranalli 1995) cannot lead to conjugate slip, which instead promotes the small scale removal of mantle lithosphere (Lei et al., 2020). In contrast, dislocation-accommodated grain boundary sliding can promote strain rate localization in weak terranes, which causes the local weak mantle lithosphere to sink into the asthenosphere by dripping (Figure 11).

(A) Second invariants of the strain rate in the weak mantle lithosphere with dislocation creep (modified from Lei et al., 2020). (B) Second invariants of the strain rate in the weak mantle lithosphere with dislocation-accommodated grain boundary sliding (Table 1).

Systematic numerical modeling suggests that the rheological range of the small-scale removal of the mantle lithosphere and the associated critical width depends on the heterogeneity in the rheological structure. The mantle lithosphere removal to the Moho in small-scale regions requires low plastic yielding stress ( λ = 0.001 ) (Figure 5). In comparison, the model results with localized weak mantle lithosphere and strong crust indicate that obvious removal can occur with the condition of low plastic yielding stress ( λ = 0.001 ) and larger critical width (∼300 km) than is the case with weak block/terranes (Figure 6).

 4.2 Geological applications

The history of surface uplift and magmatism is key to discerning delamination or convective thinning under orogens (Göğüş & Pysklywec 2008a; Göğüş & Pysklywec 2008b; Li et al., 2016; Huangfu et al., 2018; Lei et al., 2019). However, it is difficult to identify the patched or fragmental removal of mantle under orogens by means of the two geological indicators because patched removal causes neither symmetric pattern of convective thinning nor migration of delamination.

Adiabatic upwelling of asthenospheric mantle triggered by dripping has been regarded as the most significant expected geological response (Kay & Kay 1993; Ducea & Saleeby 1998). However, it is difficult for small drip with short wavelength to remove the whole mantle lithosphere, which further restricts adiabatic upwelling of the asthenosphere in a few million years (Beall et al., 2017). Thus, it cannot account for the patched removal of lithospheric mantle (to the Moho) and rapid magmatism in the northeastern Tibetan Plateau, Tianshan, and central Andes (Ducea 2011; Ducea et al., 2013; Deng et al., 2018). The results of numerical modeling shed light on the tectonic evolution of orogens with localized weak mantle of metasomatism or the suture between terranes.

4.2.1 Northeastern Tibetan Plateau and Tianshan

The growth of the northeastern TP has been accompanied by continuous lithosphere-scale shortening since the Neogene, during which the inherited lithospheric heterogeneities have played a key role in this shortening (Wang et al., 2008; Yin and Dang, 2008; Zhang and Wang, 2014). The joint inversion of receiver functions reveals patches of seismically fast and slow anomalies in the mantle (Deng et al., 2018). This removal of the mantle lithosphere with patches under the northeastern Tibetan Plateau is reconciled in the models with weak lithospheric blocks (Figures 7A–D).

The central Tianshan together with the adjacent Kazakh Shield to the north and Tarim Craton to the south manifest obvious rheological heterogeneity (Bing et al., 2022; Li et al., 2022; Zhang et al., 2022). Previous studies attributed large-scale intracontinental subduction to the rapid uplift of the central Tianshan since ∼11 Myr, including the northward subduction of the Tarim Craton or the southward subduction of the Kazakh Shield (Gilligan et al., 2014; Zhang et al., 2022). However, the major shortening of the lithosphere is distributed throughout the Tianshan orogen (Thompson et al., 2002; Zubovich et al., 2010). Geophysical investigations reveal a weak block in the central Tianshan and a range of low-velocity anomalies in the lithosphere beneath the Tianshan at a depth of ∼180 km, indicating a local high temperature and weak rheology there (Lei and Zhao, 2007; Li et al., 2009). The weak mantle lithosphere and crust beneath the Tianshan accommodates strain induced by the distant India–Asia collision, which may have promoted the small-scale removal of the mantle lithosphere in the Tianshan since the Miocene. Our numerical model with localized weak crust and the mantle lithosphere shed light on the mantle lithosphere of the weak block, which is consistent with a high temperature mantle lithosphere presently under the Tianshan and can explain the onset of rapid uplift of the Tianshan since ∼11 Myr (Figure 12).

Patched/small-scale lithospheric removal and rheological heterogeneity under central Tianshan. (A) Model result with a weak block, based on dislocation-accommodated grain boundary sliding and plastic yielding ( λ = 0.001 ). (B) Weaker lithosphere under the central Tianshan than the surroundings (Kazakh Shield and Tarim Craton), modified from Li et al., 2022.

 4.2.2 Central Andes

Another case with small-scale removal of the mantle lithosphere is the central Andes. In contrast to the northeastern TP, the mantle lithosphere under the central Andes does not include multiple terranes. The orogeny of the central Andes commenced in the Cenozoic (particularly in the past 30 Myrs). However, the Nazca Plate has been subducting along the western margin of the South American Plate for more than 200 Myrs (Isacks 1988; Allmendinger and Gubbels, 1996; Sobolev and Babeyko, 2005). This long history of subduction has caused a migration of aqueous fluids from the subducting Nazca plate, which may weaken the localized mantle lithosphere under the central Andes (Jing et al., 2020; Contreras-Reyes and Diaz, 2021; Wu et al., 2022). Geochemical constraints and geophysical data indicate a small-scale removal of the mantle lithosphere under the central Andes, with a pattern of melting in accordance with the process of small-scale foundering/dripping (<50 km diameter) of the thickened mantle lithosphere in the Altiplano–Puna Plateau, where mafic volcanic regions on the plateau manifest individual dripping (Drew et al., 2009; Ducea et al., 2013). This may be caused by the localized removal of the weak mantle lithosphere under the Altiplano–Puna Plateau. Geophysical observations also indicate a similar patched removal of lithospheric mantle under the southern Puna Plateau. This process is related to the activity of water released from the Nazca slab. Figure 13 shows the patched/small-scale removal of the mantle lithosphere under the southern Puna Plateau and that the low-velocity body beneath Cerro Galan (L2) reaches the greatest depth as it extends up to ∼50 km deep into the crust. This indicates that the localized weak mantle lithosphere may be removed, followed by adiabatic upwelling of the asthenosphere. Meanwhile, a high velocity body (H2) is located beneath the Moho, indicating patched or small-scale removal of the lithosphere (Jing et al., 2020). These processes may be well explained by the results of numerical modeling with a localized weak lithosphere (Figure 13C), in which the local weakness of the mantle lithosphere is dominated by DisGBS and strong plastic yielding ( λ = 0.001 ).

Patched/small-scale lithospheric removal under the southern Puna Plateau. (A) Vp perturbations under the southern Puna Plateau. The high-velocity body “H2” is interpreted as the small-scale removed lithosphere (modified from Jing et al., 2020), which causes the upwelling of fluids and melts into the crust, triggering the volcanic eruption of Cerro Galan (L2). (B) Interpretation cartoon under the southern Puna Plateau (modified from Jing et al., 2020). Water was released from the slab due to the dehydration reaction, which promotes the localized weak mantle lithosphere and further causes small-scale thinning under the southern Puna Plateau. (C) Model result with a weak mantle lithosphere, based on dislocation-accommodated grain boundary sliding and plastic yielding ( λ = 0.001 ).

 5 Conclusion

Using high-resolution thermomechanical modeling, we have numerically investigated the small-scale removal of the mantle lithosphere under orogens. Based on both localized weak mantle lithosphere and numerous terranes, we explored the water fugacity and plastic yielding for the model of the small-scale removal of the mantle lithosphere. The main conclusions from this study are as follows:

(1) Classical convective dripping cannot effectively thin the mantle lithosphere on a small scale because it could occur with a large initial wavelength and an amplitude of perturbation that even comprises dislocation-accommodated grain boundary sliding (DisGBS) and plasticity.

(2) Patched removal of the mantle lithosphere is induced by a localized weak block (weak crust and mantle lithosphere) with DisGBS and low plastic yield stress, which can occur in the model with the width of the weak block larger than 200 km, similar to the central Tianshan. The strain rate localization extends through the whole mantle lithosphere by conjugate strike-slip, which causes patched removal of the mantle lithosphere. In comparison, when a low rheological strength is only set to the mantle lithosphere, dripping can effectively remove the mantle lithosphere with the larger width of the weak mantle (∼300 km) than the weak block. This result can better account for the patched removal of the mantle lithosphere under the Altiplano–Puna Plateau in the central Andes, which may have been weakened by previous terrane accretion or oceanic subduction.

(3) Lithospheric heterogeneities due to terrane accretion can promote patched thinning beneath weak terranes. Such segmental removal of the mantle lithosphere differs from models with a homogeneous rheological structure. These results can better explain the patched thinning of the mantle lithosphere under the northeastern Tibetan Plateau, where the lithosphere consists of multiple terranes.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding authors.

 Author contributions

TL and JL designed the study. TL wrote the manuscript. Z-HL, ML, and LL helped improve the manuscript with some suggestions. All authors contributed to the article and approved the submitted version.

 Funding

This research is financially supported by NSFC projects (U1901214, 41974104, and 91855208) and the Guangdong project 2017ZT07Z066. Numerical simulations were run with the clusters of the National Supercomputer Center in Guangzhou (Tianhe-II).

 Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 References Allmendinger R. W. Gubbels T. (1996). Pure and simple shear plateau uplift, Altiplano-Puna, Argentina and Bolivia. Tectonophysics 259, 1–13. 10.1016/0040-1951(96)00024-8 Bajolet F. Galeano J. Funiciello F. Moroni M. Negredo A. M. Faccenna C. (2012). Continental delamination: Insights from laboratory models. Geochem. Geophys. Geosystems 13, 1525–2027. 10.1029/2011gc003896 Beall A. P. Moresi L. Stern T. (2017). Dripping or delamination? A range of mechanisms for removing the lower crust or lithosphere. Geophys. J. Int. 210, 671–692. 10.1093/gji/ggx202 Beck S. L. Zandt G. (2002). The nature of orogenic crust in the central Andes. J. Geophys. Res. Solid Earth 107 (10), ESE7-1–ESE7-16. 10.1029/2000jb000124 Billen M. I. (2008). Modeling the dynamics of subducting slabs. Annu. Rev. Earth Planet 36, 325–356. 10.1146/annurev.earth.36.031207.124129 Bing F. Z. Xue W. B. Yi X. X. (2022). Seismic anisotropy in the central Tien Shan unveils rheology-controlled deformation during intracontinental orogenesis. Geology 50, 812–816. 10.1130/g49633.1 Bird P. (1979). Continental delamination and the Colorado plateau. J. Geophys. Res. Solid Earth 84 (13), 7561–7571. 10.1029/jb084ib13p07561 Bird P. (1978). Initiation of intra continental subduction in the Himalaya. J. Geophys. Res. Solid Earth 83 (10), 4975–4987. 10.1029/jb083ib10p04975 Bittner D. Schmeling H. (1995). Numerical modelling of melting processes and induced diapirism in the lower crust. Geophys. J. Int. 123, 59–70. 10.1111/j.1365-246x.1995.tb06661.x Burg J. P. Gerya T. V. (2005). The role of viscous heating in barrovian metamorphism of collisional orogens: Thermomechanical models and application to the lepontine dome in the central alps. J. Metamorph. Geol. 23, 75–95. 10.1111/j.1525-1314.2005.00563.x Clauser C. Huenges E. (1995). “Thermal conductivity of rocks and minerals,” in Rock Physics and phase relations. Editor Ahrens T. J. (Washington D. C.: AGU), 105–126. Reference Shelf 3. Conrad C. P. (2000). Convective instability of thickening mantle lithosphere. Geophys. J. Int. 143, 52–70. 10.1046/j.1365-246x.2000.00214.x Conrad C. P. Molnar P. (1999). Convective instability of a boundary layer with temperature and strain-rate-dependent viscosity in terms of ‘available buoyancy. Geophys. J. Int. 139, 51–68. 10.1046/j.1365-246x.1999.00896.x Conrad C. P. Molnar P. (1997). The growth of Rayleigh–Taylor type instabilities in the lithosphere for various rheological and density structures. Geophys. J. Int. 129 (1), 95–112. 10.1111/j.1365-246x.1997.tb00939.x Contreras-Reyes E. Diaz D. Bello-Gonzalez J. P. Slezak K. Potin B. Comte D. (2021). Subduction zone fluids and arc magmas conducted by lithospheric deformed regions beneath the central Andes. Sci. Rep. 11, 23078. 10.1038/s41598-021-02430-9 Currie C. A. Huismans R. S. Beaumont C. (2008). Thinning of continental backarc lithosphere by flow-induced gravitational instability. Earth Planet. Sci. Lett. 269 (3-4), 436–447. 10.1016/j.epsl.2008.02.037 Deng Y. F. Li J. T. Song X. D. Zhu L. (2018). Joint inversion for lithospheric structures: Implications for the growth and deformation in northeastern Tibetan plateau. Geophys. Res. Lett. 45 (9), 3951–3958. 10.1029/2018gl077486 Drew S. T. Ducea M. N. Schoenbohm L. M. (2009). Mafic volcanism on the Puna Plateau, NW Argentina: Implications for lithospheric composition and evolution with an emphasis on lithospheric foundering. Lithosphere 1, 305–318. 10.1130/l54.1 Ducea M. N. (2011). Fingerprinting orogenic delamination. Geology 39 (2), 191–192. 10.1130/focus022011.1 Ducea M. N. Seclaman A. C. Murray K. E. Jianu D. Schoenbohm L. (2013). Mantle-drip magmatism beneath the Altiplano-Puna plateau, central Andes. Geology 41 (8), 915–918. 10.1130/g34509.1 Ducea M. Saleeby J. B. (1998). A case for delamination of the deep batholithic crust beneath the Sierra Nevada, California. Int. Geol. Rev. 40, 78–93. 10.1080/00206819809465199 Fang X. Liu D. Song C. Dai S. Meng Q. (2013). Oligocene slow and Miocene-quaternary rapid deformation and uplift of the Yumu Shan and north Qilian Shan: Evidence from high-resolution magnetostratigraphy and tectonosedimentology. Geol. Soc. Lond. Spec. Publ. 373 (1), 149–171. 10.1144/SP373.5 Gerya T. V. Meilick F. I. (2011). Geodynamic regimes of subduction under an active margin: Effects of rheological weakening by fluids and melts. J. Metamorph. Geol. 29 (1), 7–31. 10.1111/j.1525-1314.2010.00904.x Gerya T. V. Yuen D. A. (2003). Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties. Phys. Earth Planet. Interiors 140 (4), 293–318. 10.1016/j.pepi.2003.09.006 Gilligan A. Roecker S. W. Priestley K. F. Nunn C. (2014). Shear velocity model for the Kyrgyz Tien Shan from joint inversion of receiver function and surface wave data. Geophys. J. Int. 199, 480–498. 10.1093/gji/ggu225 Gogus O. H. Pysklywec R. N. (2008a). Mantle lithosphere delamination driving plateau uplift and synconvergent extension in eastern Anatolia. Geology 36 (9), 723–726. 10.1130/g24982a.1 Gogus O. H. Pysklywec R. N. (2008b). Near-surface diagnostics of dripping or delaminating lithosphere. J. Geophys. Reserch 113, B11404. 10.1029/2007jb005123 Gogus O. H. Pysklywec R. N. Şengör A. Gün E. (2017). Drip tectonics and the enigmatic uplift of the central anatolian plateau. Nat. Commun. 8, 1538. 10.1038/s41467-017-01611-3 Hansen L. N. Zimmerman M. E. Kohlstedt D. L. (2011). Grain boundary sliding in San Carlos olivine: Flow law parameters and crystallographic-preferred orientation. J. Geophys. Reserch 116, B08201. 10.1029/2011jb008220 Houseman G. A. Molnar P. (1997). Gravitational (Rayleigh-Taylor) instability of a layer with non-linear viscosity and convective thinning of continental lithosphere. Geophys. J. R. Astronomical Soc. 128, 125–150. 10.1111/j.1365-246x.1997.tb04075.x Huangfu P. P. Li Z. H. Gerya T. Fan W. Zhang K. J. Zhang H. (2018). Multi-terrane structure controls the contrasting lithospheric evolution beneath the Western and central–eastern Tibetan plateau. Nat. Commun. 9, 3780. 10.1038/s41467-018-06233-x Huangfu P. P. Li Z. H. Zhang K. J. Fan W. Zhao J. Shi Y. (2021). India-Tarim lithospheric mantle collision beneath Western Tibet controls the Cenozoic building of Tian Shan. Geophys. Res. Lett. 48, e2021GL094561. 10.1029/2021GL094561 Isacks B. L. (1988). Uplift of the Central Andean plateau and bending of the Bolivian orocline. J. Geophys. Res. 93, 3211–3231. 10.1029/jb093ib04p03211 Jing C. Sofia K. K. Yuan X. Heit B. Wu H. Yang D. Schurr B. (2020). Lithospheric delamination beneath the southern Puna plateau Resolved by local earthquake tomography. J. Geophys. Res. 125, e2019JB019040. 10.1029/2019JB019040 Kay R. W. Kay S. M. (1993). Delamination and delamination magmatism. Tectonophysics 219, 177–189. 10.1016/0040-1951(93)90295-u Kirby S. H. Kronenberg A. K. (1987). Rheology of the lithosphere: Selected topics. Rev. Geophys 25, 1219–1244. 10.1029/rg025i006p01219 Kirby S. H. (1983). Rheology of the lithosphere. Rev. Geophys. Space Phys. 21, 1458–1487. 10.1029/rg021i006p01458 Lei J. Zhao D. (2007). Teleseismic P-wave tomography and the upper mantle structure of the central Tien Shan orogenic belt. Phys. Earth Planet. Interiors 162, 165–185. 10.1016/j.pepi.2007.04.010 Lei T. Li Z. H. Liu M. (2020). Dynamics of small-scale removing of mantle lithosphere under orogens. Chin. J. Geophys. 63 (10), 3727–3739. Lei T. Li Z. H. Liu M. (2019). Removing mantle lithosphere under orogens: Delamination versus convective thinning. Geophys. J. Int. 219 (2), 877–896. 10.1093/gji/ggz329 Levander A. Berzada M. J. Niu F. Humphreys E. D. Palomeras I. Thurner S. M. (2014). Subduction-driven recycling of continental margin lithosphere. Nature 515, 253–256. 10.1038/nature13878 Levander A. Schmandt B. Miller M. Liu K. Karlstrom K. Crow R. (2011). Continuing Colorado plateau uplift by delamination-style convective lithospheric downwelling. Nature 472 (7344), 461–465. 10.1038/nature10001 Li W. Chen Y. Yuan X. H. Xiao W. Windley B. F. (2022). Intracontinental deformation of the Tianshan Orogen in response to India-Asia collision[J]. Nat. Commun. 13, 3738. 10.1038/s41467-022-30795-6 Li Z. H. Liu M. Gerya T. V. (2016). Lithosphere delamination in continental collisional orogens: A systematic numerical study. J. Geophys. Res. Solid Earth 121, 5186–5211. 10.1002/2016jb013106 Li Z. W. Roecker S. Zhihai L. Bin W. Haitao W. Schelochkov G. (2009). Tomographic image of the crust and upper mantle beneath the Western Tien Shan from the MANAS broadband deployment: Possible evidence for lithospheric delamination. Tectonophysics 477, 49–57. 10.1016/j.tecto.2009.05.007 Liao J. Wang Q. Taras G. Ballmer M. D. (2017). Modeling Craton destruction by hydration-induced weakening of the upper mantle. J. Geophys. Res. Solid Earth 122, 7449–7466. 10.1002/2017jb014157 Lorinczi P. Houseman G. (2009). Lithospheric gravitational instability beneath the southeast Carpathians. Tectonophysics 474 (1), 322–336. 10.1016/j.tecto.2008.05.024 Molnar P. Houseman G. A. Conrad C. P. (1998). Rayleigh—taylor instability and convective thinning of mechanically thickened lithosphere: Effects of non-linear viscosity decreasing exponentially with depth and of horizontal shortening of the layer. Geophys. J. Int. 133 (3), 568–584. 10.1046/j.1365-246x.1998.00510.x Morency C. Doin M. P. (2004). Numerical simulations of the mantle lithosphere delamination. J. Geophys. Res. Solid Earth 109 (3), B03410. 10.1029/2003jb002414 Owens T. J. Zandt G. (1997). Implications of crustal property variations for models of Tibetan plateau evolution. Nature 387, 37–43. 10.1038/387037a0 Ranalli G. Murphy D. C. (1987). Rheological stratification of the lithosphere. Tectonophysics 132 (4), 281–295. 10.1016/0040-1951(87)90348-9 Ranalli G. (1995). Rheology of the Earth, deformation and flow process in Geophysics and geodynamics . 2nd ed. London: Chapman & Hall, 413. Schott B. Schmeling H. (1998). Delamination and detachment of a lithospheric root. Tectonophysics 296 (3-4), 225–247. 10.1016/s0040-1951(98)00154-1 Schurr B. Rietbrock A. Asch G. Kind R. Oncken O. (2006). Evidence for lithospheric detachment in the central Andes from local earthquake tomography. Tectonophysics 415 (1-4), 203–223. 10.1016/j.tecto.2005.12.007 Sobolev S. V. Babeyko A. Y. (2005). What drives orogeny in the Andes? Geology 33, 617–620. 10.1130/g21557ar.1 Thompson S. C. Weldon R. J. Rubin C. M. Abdrakhmatov K. Molnar P. Berger G. W. (2002). Late quaternary slip rates across the central Tien Shan, Kyrgyzstan, central Asia. J. Geophys. Res. Solid Earth 107, ETG 7-1–ETG 7-32. 10.1029/2001jb000596 Tomohiro O. Takaaki K. Higo Y. Funakoshi K. I. Suzuki A. Kikegawa T. (2015). Dislocation-accommodated grain boundary sliding as the major deformation mechanism of olivine in the Earth’s upper mantle. Sci. Adv. 1 (9), e1500360. 10.1126/sciadv.1500360 Turcotte D. L. Schubert G. (2002). Geodynamics. 2nd ed. New York: Cambridge Univ, 456. van Hunen J. Zhong S. Shapiro N. M. Ritzwoller M. H. (2005). New evidence for dislocation creep from 3-D geodynamic modeling of the Pacific upper mantle structure. Earth Planet, Lett. 238, 146–155. 10.1016/j.epsl.2005.07.006 Vogt K. Gerya T. Castro A. (2012). Crustal growth at active continental margins: Numerical modeling. Phys. Earth Planet. Interiors 192-193, 1–20. 10.1016/j.pepi.2011.12.003 Wang H. L. Currie C. A. (2015). Magmatic expressions of continental lithosphere removal. J. Geophys. Res. Solid Earth 120, 7239–7260. 10.1002/2015jb012112 Wang Q. Wyman D. A. Xu J. Dong Y. Vasconcelos P. M. Pearson N. (2008). Eocene melting of subducting continental crust and early uplifting of central Tibet: Evidence from central-Western Qiangtang high-K calc-alkaline andesites, dacites and rhyolites. Earth Planet. Sci. Lett. 272 (1-2), 158–171. 10.1016/j.epsl.2008.04.034 Wu Y. Liao J. Guo F. Wang X. C. Shen Y. (2022). Styles of trench-parallel mid-ocean ridge subduction affect Cenozoic geological evolution in circum-Pacific continental margins. Geophys. Res. Lett. 49. 10.1029/2022GL098428 Xia L. Q. Li X. M. Ma Z. P. Xu X. Xia Z. (2011). Cenozoic volcanism and tectonic evolution of the Tibetan Plateau. Gondwana Res. 19 (4), 850–866. 10.1016/j.gr.2010.09.005 Yin A. Dang Y. Q. Zhang M. Chen X. H. McRivette M. W. (2008). Cenozoic tectonic evolution of the Qaidam basin and its surrounding regions (Part 3): Structural geology, sedimentation, and regional tectonic reconstruction. Geol. Soc. Am. Bull. 120 (7-8), 847–876. 10.1130/B26232.1 Zhang B. F. Bao X. W. Xu Y. X. (2022). Seismic anisotropy in the central Tien Shan unveils rheology-controlled deformation during intracontinental orogenesis. Geology 50, 812–816. 10.1130/g49633.1 Zhang Z. Wang Y. Houseman G. A. Xu T. Wu Z. Yuan X. (2014). The Moho beneath Western Tibet: Shear zones and eclogitization in the lower crust. Earth Planet. Sci. Lett. 408, 370–377. 10.1016/j.epsl.2014.10.022 Zubovich A. V. Wang X. q. Scherba Y. G. Schelochkov G. G. Reilinger R. Reigber C. (2010). GPS velocity field for the Tien Shan and surrounding regions. Tectonics 29, TC6014. 10.1029/2010tc002772