Patched removal of the mantle lithosphere under orogens: A systematic numerical study

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.


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 mantlederived magmatism may have occurred within a volcanic field over short time scales (1)(2)(3)(4)(5) 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 thermomechanical 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 smallscale removal of the mantle lithosphere under orogens.
2 Numerical model

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: 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α(DP/DT), α is thermal expansivity, and H s is the shear heating rate (H s σ ij ′ _ ε ij , _ ε strain rate tensor). For the density of a specific rock type, ρ depends on pressure (P) and temperature (T): 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.

Viscoplastic rheology
1 Viscous rheology where _ ε II (_ ε ij _ ε ij ) 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 H2O , 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.
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 An 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(φ dry )λ 0.6~0.3) is used for the dry and strong lithospheric mantle.

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: where η plastic is the viscosity of the Drucker-Prager plasticity, σ yield is the yield stress, P is the dynamic pressure, C 0 is the residual rock strength at p=0, and φ dry 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 fluid/melt P .  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).  Table 1. e strain weakening effect is applied for plastic rheology, in which both cohesion (C 0 ) and the effective friction coefficient sin(φ dry )λ decrease with strain increase. The cut-off values shown in the table correspond to the strain '0-1'. f 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.

Frontiers in Earth Science frontiersin.org
Because the fluid/melt parameter P f luid/melt is not directly calculated in the model, we use a range of "sin(φ dry )λ" 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(φ dry )λ 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):

Initial model configuration and boundary conditions
Large-scale models (4000 × 400 km) 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,

FIGURE 1
Reference model and boundary conditions. (A) Composition field of the homogenous model with initial wavelength (u=150 km) and perturbation (Z 0 10km). 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.5cm/yr) 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.

Frontiers in Earth Science
frontiersin.org 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℃•km −1 . 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: zv x /zy 0 and zv y /zy −v y /Δy external , where Δy external is the vertical distance from the lower boundary to the external boundary where free slip (zv x /zy 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).

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 Frontiers in Earth Science frontiersin.org mantle lithosphere (Bird 1978(Bird , 1979Lei 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.

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 20km) at 3.4 Myr, where convergent velocity ] x ± 4.5 + 4.5cm/yr. When the initial perturbation wavelength is small (u 150km − 300km), the drips cannot effectively remove the mantle lithosphere (Figures 2A, B). In contrast, with a larger wavelength initial perturbation (u 500km − 800km), the lower mantle lithosphere can be easily removed by the drips (Figures 2C, D).

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 500km) at 3.4 Myr. With a small perturbation amplitude ( Figure 3A, Z 0 1km), no convective dripping occurs. When the initial perturbation amplitude Z 0 10km, 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 20km), 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 30km, the drips effectively induced a removal of the mantle lithosphere from 150 km to~50 km ( Figure 3D).

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  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.

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.
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

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 is 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). Frontiers in Earth Science frontiersin.org

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).

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.

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, Ra 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 Ra 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 Ra 3 as a function of layer thickness at different initial perturbation wavelengths and amplitudes, respectively. At the small amplitude (Z 0 1 km), the analysis indicates that the whole mantle lithosphere remains stable, with Ra 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 Frontiers in Earth Science frontiersin.org ( 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 superexponentially grow and further remove part of the mantle lithosphere ( Figures 9B-D). The thickness h crit represents the minimum amount of mantle lithosphere that may be removed by gravitational instability, which indicates that the critical Ra 3 of 100 is obtained for layer thickness. In the models with the larger amplitude (Z 0 20 km), the calculated critical Ra 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 km) is larger than the critical thickness h crit , 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 Ra 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 km, super-exponential growth develops at the base of the mantle lithosphere ( Figure 3C) because the initial  (Table 1) and plastic yielding (λ 0.001) are used for the semianalytical model (Eq. 6).
Frontiers in Earth Science frontiersin.org amplitude (Z 0 20 km) is larger than the critical thickness h crit (~18 km), corresponding to a wavelength of 500 km ( Figure 9C). The thickness h crit only estimates the minimum amount of the layer of removal by instability. Larger thicknesses are also potentially unstable (Ra 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

FIGURE 10
Evolution of the average depth of the 1350℃ isotherm for the model of Type I with the constant initial perturbation Z 0 20km 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.  Lei et al., 2020). (B) Second invariants of the strain rate in the weak mantle lithosphere with dislocation-accommodated grain boundary sliding (Table 1).

Frontiers in Earth Science
frontiersin.org thickening develops super-exponential growth ( Figures 10A, B). Second, the model with DisGBS ( Figure 10A) shows that the thickness of the removed layer (Ra 3 > 100) is~30 km. In comparison, the larger amount of the removed layer (~110 km) of Ra 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 (Ra 3 > 100).

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 largescale 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).
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).

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.

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

Frontiers in Earth Science
frontiersin.org 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).

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).

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.