Early Miocene Exhumation of High-Pressure Rocks in the Himalaya: A Response to Reduced India-Asia Convergence Velocity

Low-viscosity channel flow, originating from a melt-weakened mid-crustal layer, is one of the most popular tectonic models to explain the exhumation of deep-seated rocks in the Greater Himalayan Sequence (GHS). The driving mechanism of such channel flow, generally attributed to focused erosion in the mountain front, is still debated, and yet to be resolved. Moreover, the channel flow model cannot explain eclogites in the GHS. In this study, we present a new two-dimensional thermo-mechanical numerical model, based on lubrication dynamics to demonstrate the exhumation process of deep crustal rocks in GHS. The model suggests that a dynamic-pressure drop in the Himalayan wedge, following a large reduction in the India-Asia convergence velocity (15 cm/yr at 50 Ma to nearly 5 cm/yr at ∼22 Ma) localized a fully developed extrusion zone, which controlled the pressure-temperature-time (P-T-t) path of GHS rocks. We show that the wedge extrusion, originated in the lower crust (∼60 km), was initially bounded by two oppositely directed ductile shear zones: the South Tibetan Detachment systems (STDS) at the top and the Higher Himalayan Discontinuity (HHD) at the bottom. With time the bottom shear boundary of the extrusion zone underwent a southward migration, forming the Main Central Thrust (MCT) at ∼17 Ma. Our model successfully reproduces two apparently major paradoxical observations in the Himalaya: syn-convergence extension and inverted metamorphic isograds. Model peak P (10–17 kb) and T (700–820°C) and the exhumation P-T-t path estimated from several Lagrangian points, traveling through the extrusion zone, are largely compatible with the petrological observations in the GHS. The model results account for the observed asymmetric P-T distribution between the MCT and STDS, showing peak P-T values close to the MCT. The lubrication dynamics proposed in this article sheds light on the fast exhumation event (>1 cm/yr) in the most active phase of crustal extrusion (22-17 Ma), followed by a slowed-down event. Finally, our model explains why the extrusion zone became weak in the last 8-10 Ma in the history of India-Asia collision.


INTRODUCTION
Understanding the exhumation mechanism of high pressure (HP) and ultra-high pressure (UHP) rocks, such as eclogites is the key to modeling the crustal recycling processes in convergent plate boundaries (Ernst et al., 1997;Agard et al., 2009). The Himalaya represents an archetype collisional orogen, displaying HP rocks within the high-grade metamorphic terrain of the Greater (or Higher) Himalayan Sequence (GHS) ( Figure 1A). There are two contrasting types of eclogites reported from the GHS (Lombardo and Rolfo, 2000;Corrie et al., 2010;Grujic et al., 2011;O'Brien, 2019). The eclogites in northwestern Himalaya (Kaghan and Tso Morari) are of Early Tertiary age, and their peak pressure (P 26-28 kb) temperatures (T 580-600°C) suggest UHP metamorphism (Lombardo and Rolfo, 2000;O'Brien, 2019). On the other hand, those in the eastern Himalaya (Ama Drime/Arun valley/Kharta, north Sikkim, NW Bhutan) are much younger in age , and the calculated peak P ( 10-15 kb) and T ( ∼780°C) indicate their exhumation depths of 50 km or more (Lombardo and Rolfo, 2000;Corrie et al., 2010;Kohn, 2014;O'Brien, 2019). Extrusion mechanism of HP rocks from such depths has emerged as a stimulating topic of research in geosciences, giving way to a remarkable conceptual development in convergent tectonics (Mancktelow, 1995;Chemenda et al., 2000;Burov and Yamato, 2008;Warren et al., 2008;Li and Gerya, 2009). Diverse types of mechanical and thermo-mechanical models have been proposed to explain the extrusion tectonics of high-grade rocks in the GHS  and the development of fold-thrust belts in the Himalaya. These include the coulomb wedge model (Dahlen et al., 1984), the wedge extrusion model (Burchfiel and Royden, 1985;Hodges et al., 1992), the slab break-off model of the Indian crustal lithosphere during 25-20 Ma (Chemenda et al., 2000), the channel flow model (Grujic et al., 1996;Beaumont et al., 2001), the critical taper model (Kohn, 2008), and the insequence shearing model (Carosi et al., 2016;Carosi et al., 2018). However, the channel flow model, which explains the formation of GHS as a focused exhumation of low-viscosity mid-crustal materials through a channel-like passage between the Main Central Thrust (MCT) on the south and the South Tibetan Detachment System (STDS) on the north (Beaumont et al., 2001;Grujic et al., 2002) ( Figure 1B), revolutionized the idea about the Himalayan tectonics. According to this model, the pressure gradient between the thick elevated Tibetan plateau and a relatively thin foreland resulted in a Poiseuille flow in the large hot orogens, coupled with active erosion at the mountain fronts to facilitate the exhumation of mid-crustal rocks. Similarly, Clark and Royden, 2000 conceptualized channel flow dynamics of low-viscosity mid-crustal materials (∼10 18 Pa s), driven by lateral pressure gradients, to interpret the variations in topography building at the northeastern, southeastern, and eastern margins of the Tibetan plateau. Several studies predicted the occurrence of similar narrow extrusion zones in subduction settings, which mediate the return flows and exhumation of ultrahigh pressure rocks (Warren et al., 2008;Li and Gerya, 2009;Malusà et al., 2015) in convergent plate margins.
The channel flow model of Beaumont et al. (2001), although it gained enormous popularity in the context of Himalayan tectonics, has triggered counter-arguments in the light of various geological issues, as presented in the special volume on the Channel Flow conference (Godin et al., 2006;Harrison, 2006) and discussed in later studies (Kohn 2008;Carosi et al., 2018). For example, its validity essentially demands the existence of a mechanically weak middle crust, with viscosity in the order of 10 19 Pa s. Such a mechanical stratification is expected to cause decoupling of the Tibetan upper crust from the stronger lower part (lower crust and mantle) of the lithosphere (Copley et al., 2011). In fact, several workers invoked a diametrically opposite view, considering strong mechanical coupling between the upper crust and the underlying lithosphere, and thereby a coherent deformation for the whole Tibetan lithosphere (Zhang et al., 2004;Flesch et al., 2005;Copley et al., 2011). Another important condition required to sustain the channel flow is focused surface-erosion at the mountain front. Geochronological and geological evidence suggests that the extrusion process in the Himalaya was most active during the Oligocene-Miocene time. Despite the availability of data for high erosion rates in this time interval (Thiede et al., 2005;Clift, 2017), there is an uncertainty with the spatial correlations between the high orogenic precipitation and the exact location of focused exhumation (Corrie et al., 2010). Another line of modeling suggests the possibility of focused extrusion in subduction wedges, which excludes the surface erosion as a pre-requisite factor for the deep-crustal exhumation process (Brun and Faccenna, 2008;Malusà et al., 2015;Piccolo et al., 2018). It is thus necessary to readdress the question-can channel-like extrusion occur in the Himalaya without any significant contribution of focused surface erosion at the mountain front? Also, the channel flow model does not account for the occurrence of highpressure eclogites in the GHS, where their peak metamorphic conditions indicate pressures in excess of 1.5 GPa and exhumation depths >50 km, as pointed out by Grujic et al. 2011 (Figure 1C). According to the channel flow model, the GHS exhumation occurred between two synchronous shear zones, the STDS on the north and the MCT on the south. However, many recent field studies suggest that GHS were exhumed between 27 to 17 Ma, prior to the activation of MCT (Wang et al., 2016;Carosi et al., 2018). These studies recognized a more than 1000 km long new shear zone, called High Himalayan Discontinuity (HHD) within the GHS ( Figures 1A,B) Cottle et al., 2015;Khanal et al., 2015;Zeiger et al., 2015;Agustsson et al., 2016;Wang et al., 2016;Iaccarino et al., 2017;Walters and Kohn, 2017;Carosi et al., 2018). The HHD holds a cogenetic relation with the STDS and GHS exhumation event (Wang et al., 2016). All these open-ended issues motivate us to re-visit the dynamics of deep-crustal extrusion and the pressure-temperature evolution of GHS rocks in the Himalaya. Geological records suggest that the phase of most rapid GHS exhumation (Ganguly et al. 2000;White et al., 2002;Najman et al., 2009;Thiede et al., 2009;Grujic et al., 2011;Carrapa et al., 2014;Clift, 2017) occurred when India-Asia convergence velocity dropped to ∼5 cm/yr at ∼22 Ma from ∼15 cm/yr at 50 Ma (Molnar and Stock, 2009;Copley et al., 2010). Interestingly, this timing marks the onset of gravitational collapse (extensional deformations, manifested in widely distributed N-S and E-W normal faulting) in the Tibetan plateau, following a prolonged period of compressional deformation since 50 Ma (Wang et al., 2014;Li et al., 2015). The two events: gravitational collapse in the Tibetan plateau and exhumation of GHS rocks in the Himalayan wedge then appear to be broadly synchronous, indicating a dynamic connection between them.
In this study, we develop a thermo-mechanical geodynamic model within the premises of lubrication theory (Hamrock, 1991;Mancktelow, 1995;Maiti and Mandal, 2020) to provide an alternative explanation for the exhumation dynamics of highpressure rocks in the Himalaya, and constrain their P-T paths during the extrusion. The model shows the drop in dynamic pressure in the Himalayan wedge as the key factor in driving the extrusion of deep-crustal materials from a depth of 55-60 km, that accounts for the occurrence of high-pressure rocks (>1.5 GPa) in the GHS (Grujic et al., 2011). Moreover, many petrological studies deciphered a fast exhumation (>1 cm/yr) event at 22-17 Ma, followed by a phase of slowed-down rates (Najman et al., 2009). We use the thermo-mechanical model to explore the cause of such an unsteady exhumation history in the GHS. We also present an analysis of the pressure (P) and temperature (T) profiles across the GHS in the light of the available P-T estimates from different transects in the Higher Himalaya. Our model estimates yield a peak P∼17 kb and T∼820°C, fairly in agreement with the petrological data reported from exhumed rocks.

CONCEPTUALIZATION OF THE HIMALAYA-TIBET DYNAMIC COUPLING
We conceptualize the Himalaya-Tibet coupled dynamics in the framework of a subduction model, where the subducting plate motion builds a dynamic pressure in the tectonic wedge to support the load in the overriding crust ( Figure 2A). Many petrological studies reported that subduction wedges (the interface between two colliding plates) consist of fluid-rich relatively low viscosity and density crustal materials (Duarte et al., 2015;Bebout and Penniston-Dorland, 2016;Bayet et al., 2018). Moreover, spontaneously evolved subduction models (Burov and Yamato, 2008;Warren et al., 2008;Li and Gerya, 2009) show the development of low-viscosity downward-tapered wedges, filled with crushed materials of upper and mid crustal rocks. Based on these evidences, we consider the Himalayan wedge as a mechanically weak region with viscosity (∼10 21 to 10 20 Pa s) and density (2600-3050 kg/m 3 ) lower than the colliding plates (∼10 22 -10 21 Pa s and 2700-3200 kg/m 3 ). Such a system is mechanically analogues to engineering thrust or pad bearing devices (Hamrock, 1991), in which two inclined plates remain separated by a relatively low-density viscous wedge (Figures 2A,B). The upper load-bearing plate is held stationary, whereas the lower plate moves at a tangential velocity and produces a flow-induced dynamic pressure (P Dy ) in the intervening viscous wedge. The P Dy supports the load of the upper plate and keeps the two plates separated from each other. P Dy is calculated using lubrication approximation (Hamrock, 1991). Maiti and Mandal (2020) have shown the applicability of lubrication approximation in calculating P Dy in subduction wedges with large convergence angle (α) up to 15°. The lubrication dynamics suggests that the subducting plate motion can develop a dynamic pressure (P Dy ) in the subduction wedge (Hamrock, 1991;Mancktelow, 1995;Maiti and Mandal, 2020), which supports the excess gravitational load of the overriding plate (P Ob ) to maintain a stable state (P Dy P Ob ) (Figure 2Ci). In case P Dy > P Ob , the system tends FIGURE 2 | (A) Schematic presentation of an orogenic wedge, forming the interface between the subducting and overriding plates in a subduction system. (B) Geometrical parameters of an orogenic wedge considered for the analysis of dynamic pressure within the wedge, where h 0 and (h 0 + S h ) are the opening widths at the two converging ends of the wedge, l is the along-plate length of the wedge, and α is convergence angle. The overriding plate is kept stationary (u o 0) and the subducting plate moves at velocity u s (cm/yr). (C) Three dynamic states of a subduction system. I: Wedge dynamic pressure (P Dy ) equals to overburden pressure (P Ob ), acting on the wedge-overriding plate interface; the system is in a balanced state. II: P Dy >P Ob ; the overriding plate (OP) undergoes contractional deformation. III. P Dy <P Ob ; the OP undergoes rotational collapse. (D) Graphical plots of P Dy as a function of I) subduction plate velocity, u s and II) wedge viscosity, µ.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 632806 to achieve the equilibrium either by increasing the inter-plate distance and/or by transferring compressive stresses into the overriding plate ( Figure 2Cii). In contrast, when P Dy < P Ob , the system would reduce the inter-plate distance to raise the dynamic pressure level in order to restore a stable state ( Figure 2Ciii). The latter condition triggers a collapse of the overriding plate onto the wedge and forces the buoyant wedge rocks to extrude out of the plate confinement zone. For a given tapered geometry of the viscous wedge ( Figure 2B), the dynamic pressure (P Dy ) is obtained from Reynold's equation as a function of subducting plate velocity (u s ), and wedge viscosity (µ)  as, where h 0 and (h 0 + S h ) are the opening widths at the two converging ends of the wedge, l is the along-plate wedge length. Eq. 1 shows that P Dy holds a positive linear relation with u s and μ. A decrease either in u s or μ would reduce the wedge dynamic pressure. For example, in a subduction wedge (l 500 km; h o 10 km, S h 110 km) with μ 4 × 10 20 Pa s, a drop of u s from 8 to 5 cm/yr lowers the maximum dynamic pressure (P max ), localizing near the narrow convergence ends of the two plates, from 5.7 to 3.6 kb ( Figure 2Di). For u s 5 cm/yr, The density profiles in the reference model. Note that we have considered materials of relatively lower density to constitute the model Himalayan wedge. This is based on the results of earlier spontaneously evolved numerical models, suggesting that the wedge consists of crushed materials derived from upper and mid-crustal rocks, therefore buoyant compared to the lithospheric plates (Burov and Yamato, 2008;Warren et al., 2008;Li and Gerya, 2009). (C) A comparative plot of normal and elevated geotherm in the Himalaya-Tibet system. Due to thermal equilibrium and radioactive heat production, the geotherm is considered to be elevated by 20 Ma (Beaumont et al., 2001). (D) The viscosity profiles for our model Tibet and Himalayan wedge have been plotted against previous studies for the same system. (E) Variation of India-Asia convergence velocity with time.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 632806 5 increase in μ by two times (i.e., 3 × 10 20 to 6 ×10 20 Pa s) raises P max also by two times (2.7-5.4 kb) (Figure 2Dii). Under the set of chosen model parameters, the variation of P Dy in the wedge ranges from 2 to 6 kb. Recent field evidence also suggests that a subduction-collision zone can produce non-lithostatic pressure of large magnitudes, reaching even 15 kb (Luisier et al., 2019;Cutts et al., 2020).
The India-Asia collision history tracks a temporally nonuniform convergence motion, beginning with a high rate (∼15 cm/yr) at ∼ 50 Ma, which remarkably slowed down to ≤5 cm/yr at ≤ 22 Ma (Molnar and Stock, 2009;Copley et al., 2010). The analytical solution of lubrication theory (Eq. 1) indicates that a reducing India-Asia collision rate might have lowered the dynamic pressure in the Himalayan wedge, and resulted in a collapse of the Tibetan overriding block, forcing the deep-crustal rocks to extrude in the form of Greater Himalayan sequence, as shown earlier by Maiti and Mandal, 2020 from a simple mechanical consideration. However, to obtain a complete picture of the peak pressure-temperature condition of exhumed high-pressure rocks, its P-T-t path, shape and active phase of the orogenic channel, spatial variations of peak P-T data in the Himalayan wedge, it is necessary to develop a timedependent thermo-mechanical modeling approach, incorporating the natural complexity of the Himalaya-Tibet system and its rheological variations. In this study, we build two-dimensional thermo-mechanical numerical models within the framework of lubrication dynamics presented above.

Basic Principles
Assuming a steady-state incompressible flow, our thermomechanical modeling involves three conservation equations for mass, momentum, and energy in the following forms: where u, μ, ρ, p, F, denote velocity vector, viscosity, density, pressure, and body force term, respectively. In Eq. 4, the thermodynamic parameters T, Cp, k, Q sh represent temperature, specific heat, thermal conductivity, and heat source respectively. The mechanical part of the computation (Eqs 2, 3) calculates the velocity field, stress, strain-rate, and deformation subjected to specified mechanical boundary conditions. On the other hand, the thermal computation solves the heat balance equation (Eq. 4) on an Eulerian grid subjected to mechanical velocity field and thermal boundary conditions. The mechanical and thermal properties are coupled through the temperature-dependent viscous powerlaw rheology used in model materials. Shear heating in the model is considered to make a small contribution to the total thermal budget of a subduction system and therefore it has been neglected. It has been shown that shear heating in the lowstrength materials of the subduction system contributes to an increase in temperature of no more than 50°C (Warren et al., 2008). The values of different model parameters and material properties chosen in our models are described in Table 1. We solved the three governing equations using a finite element code (COMSOL Multiphysics, 2015), which has been benchmarked by earlier workers for studying temperature-dependent Newtonian as well as non-Newtonian rheology in fluid mechanics (van Keken et al., 2008). Moreover, the code has been used to simulate a wide range of large-scale geodynamic processes such as the mid-oceanic ridge system (Montési and Behn, 2007), subduction zone modeling (Rodríguez-González et al., 2012;Dasgupta and Mandal, 2018), magma upwelling  [Flesch et al., 2001;Liu and Yang, 2003;Copley and McKenzie, 2007;Warren et al., 2008;Li and Gerya, 2009] Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 632806 (Shahraki and Schmeling, 2012), mantle convection (He, 2014), Rayleigh-Taylor instabilities (Ghosh et al., 2020), continental deformations (Bischoff and Flesch, 2019). We have employed Arbitrary Lagrangian and Eulerian method (ALE) (Fullsack, 1995;Dasgupta and Mandal, 2018) in our timedependent model runs to advect the model upper surface with the flow and to measure the development of model topography. The most remarkable advantage of ALE is that it can explicitly track free surfaces, internal boundaries, and any moving material interfaces applicable to multiphase flow simulations. Moreover, ALE can efficiently measure large strains in simulating continuum scale deformations. The ALE method has been implemented in simulating a wide range of geodynamic processes including coupled crust-mantle dynamics and intra-cratonic deformations (Pysklywec and Cruden, 2004), slab detachment processes (Pysklywec et al., 2003), surface topography in bi-vergent subduction (Dasgupta and Mandal, 2018) and coupled crust-mantle systems (Crameri et al., 2012). Dasgupta and Mandal (2018) benchmarked the ALE scheme for computational fluid dynamics (CFD) simulations using COMSOL Multiphysics.

Model Geometry
Our initial/reference model is approximated to the Himalayan wedge geometry (discussed above) and adjoining Tibetan plateau topography with an average elevation of 5 km, which is supposed to have attained by 26 Ma in the collision history (DeCelles et al., 2007;Wang et al., 2014). The 2D model geometry, chosen in Cartesian space, has a horizontal and a vertical dimension of 2500 km and 410 km respectively, considered significantly larger than the region of interest to avoid the influence of model boundaries in the simulation results ( Figure 3A). The model replicates the Himalayan wedge in between the hinge line and the Indus-Tsangpo suture zone, with a width of approximately 300-350 km on the surface, as estimated from trans-orogen topographic profiles in the Himalayan Mountain belt , extending down to the Asian lithospheric base. We use available lithospheric cross-sections of the Himalaya-Tibet system to constrain the tapered wedge geometry (Hauck et al., 1998;Jiménez-Munt et al., 2008;Godin and Harris, 2014), and varied the subduction angle of the Indian plate between 6°to 15°w ith depth (Kumar et al., 2006;Jiménez-Munt et al., 2008;Gloaguen and Ratschbacher, 2011). The thicknesses of the Indian and overriding Asian lithosphere are chosen 110 and 130 km, respectively (Jiménez-Munt et al., 2008;Gloaguen and Ratschbacher, 2011). The Tibetan lithosphere is treated as a three-layer structure: upper to middle crust (45 km), lower crust (25 km), and mantle lithosphere (70-130 km). To reproduce the Himalaya-Tibet collision tectonics, we allowed the Indian lithosphere to subduct beneath the stationary Tibetan lithosphere (Beaumont et al., 2001), which had a stable northern boundary since 50-34 Ma (Meyer et al., 1998). The plate convergence is accommodated by the upper crust of the Indian plate, accreting into the Himalayan wedge and indenting into the Asia plate (Yakovlev and Clark, 2014).

Initial Conditions and Material Properties
We impose a set of mechanical and thermal boundary conditions in the initial/reference model ( Figure 3A). The top model boundary is subjected to an atmospheric pressure condition that allows the top surface to deform and crustal materials to extrude upward. The left end of the Indian plate is subjected to an inflow velocity u (cm/yr), whereas the right end of the Tibetan lithosphere is assigned a fixed boundary (i.e., u 0, stationary condition). The lower boundary of the mantle is permeable, which allows the material to flow out of the model domain and maintain the mass balance. The left and right vertical boundaries of the mantle are assigned to a free-slip condition.
Following the approach used in previous Himalaya-Tibet models (Jamieson et al., 2004;Rey et al., 2010), we consider two configurations of the thermal structure for our modeling ( Figure 3C). In the case of the initial model (u s 8 cm/yr, before 30 Ma), the thermal structure is defined by the normal geotherm ( Figure 3C) for continental lithosphere with T 25°C at the surface, T 400°C at 40 km depth, T 600°C at 70 km depth, and T 1,330°C at the lithosphere-asthenosphere boundary, 130 km depth (Turcotte and Schubert, 2002), which is assumed to have prevailed in the Himalaya-Tibet system earlier than 30 Ma. This inference gets support from the triaxial experiments of Rey et al. (2010) that predict the Moho temperature beneath the Tibetan plateau (at 70 km depth) must be around 500-600°C to develop a topographic elevation of 5 km in the Tibetan plateau by 30 Ma. However, the thermal gradient steepened subsequently due to crustal thickening and radiogenic heat production, and the Moho temperature beneath the plateau (at 70 km depth) possibly attained a higher value, ∼850°C at ∼25 Ma (Rey et al., 2010). In the second case (u s 5 cm/yr at 22 Ma), we thus choose an elevated geotherm ( Figure 3C) with T 25°C at surface, T 550°C at 40 km depth, T 850°C at 70 km depth, and T 1,330°C at the lithosphere-asthenosphere boundary (130 km depth). In both cases, the bottom model boundary (410 km depth) is set at T 1,500°C (Turcotte and Schubert, 2002;Li et al., 2011). The right and the left model boundaries are insulated, allowing no lateral heat outflux or influx into the model. Based on earlier numerical model data and geophysical estimates, the three principal layers of the Tibetan lithosphere: upper/middle crust, lower crust, and mantle lithosphere are assigned varying densities and viscosities (Table 1). We consider a power-law decay function (Sarkar et al., 2014) to represent decreasing viscosities (μ) of model materials with temperature (T): where μ 0 and T 0 are the reference viscosity and temperature, respectively; n is the power-law exponent, and µ is the viscosity at T. For different crustal layers, we fix μ 0 , and use n ( Table 1) to obtain the average integrated viscosities, in agreement with the bounds ( Figure 3D) estimated by several previous workers for modeling the Himalaya-Tibet system (Flesch et al., 2001;Liu and Yang, 2003;Copley and Mckenzie, 2007;Bischoff and Flesch, 2019). The temperature-dependent viscosity for each crustal layer is provided in Supplementary Figure S1A. It is to be noted that Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 632806 many previous numerical models have used exponential decay functions to account for reducing viscosity with an increase in temperature (Burov and Yamato, 2008;Warren et al., 2008;Li and Gerya, 2009). However, to keep numerical computation simple and less extensive, we have used a power-law function to obtain the viscosity-temperature regression profile in our numerical experiments. Figure 3D presents the model integrated viscosity profile in the perspective of the estimated viscosity profiles for the Himalaya Tibet system available in the literature (Liu and Yang, 2003;Copley and Mckenzie, 2007;Bischoff and Flesch, 2019). The plot suggests that our model integrated viscosity profile reasonably follows the bounds shown in the previous studies. We also consider a depth-dependent density for the Tibetan lithosphere ( Figure 3B; Supplementary Figure S1B). This is constrained from the geophysical estimates, as summarized in Jiménez-Munt et al., 2008. The values of all parameters used in the model are listed in Table 1. The 2D model domain is discretized into free triangular meshes. This leads to fine meshing (0.05-0.2 km) in the wedge and model boundaries, and a normal size (0.1-5 km) in the rest of the model areas. We have implemented an adaptive mesh refinement scheme in the model to obtain very accurate results in deforming regions during the model run time. The ALE scheme tracks the moving surfaces and measures model topographic development.
We ran a series of two-dimensional thermo-mechanical numerical model experiments ( Figure 3A), to find how reducing convergence rates, from 15 cm/yr to ≤5 cm/yr at ≤22 Ma controlled the deep-crustal extrusion processes and associated thermal structures in the Himalayan wedge. To investigate this, our initial/reference model is first subjected to a high-convergence velocity (8 cm/yr) condition, simulating the collision kinematics before 30 Ma (Copley et al., 2010). We then reduced the convergence velocity to 5 cm/yr, as applicable to the kinematic condition for 22 Ma (Molnar and Stock, 2009).

MODEL RESULTS
Model With High Convergence Velocity (u s ≥ 8 cm/yr) The reference model run with a high convergence velocity (8 cm/ yr) showed no sign of crustal extrusion (v ∼ 0) in the Himalayan wedge, but dominated by sub-horizontal flows (u 2-4 cm/ yr) in the convergence direction ( Figure 4A). The Himalayan wedge mediated the compressive stresses from the subducting plate to the overriding plate, resulting in wide-spread horizontal crustal shortening in the latter with high contractional strain rates (−6 × 10 −14 s −1 ) ( Figure 4B). The model suggests that the India-Asia collision at u s 8 cm/yr gave rise to a global crustal shortening without any appreciable focused uplift in the Himalayan wedge. At this high convergence rate, the subducting slab drags the crustal materials deep into the wedge, resulting in downward advection of isotherms, as typically reported from subduction zones ( Figure 4C). We consider a set of locations at varying depths (15-55 km) in the model wedge for quantitative analysis of the spatio-temporal changes of temperature (T) and burial/extrusion rates. T at selected locations decreases in agreement with the downward advection of isotherms ( Figure 4D). For example, at a depth of 30 km, the initial temperature (300°C) decreases to 100°C on a simulation run time of 6 Ma. We also plot spatio-temporal variations of the burial/extrusion rates in the model wedge for this high-velocity condition ( Figure 4E). The subducting slab drags crustal materials downward for burial (negative vertical velocity components) at rates of −0.2 to −0.0.3 cm/yr measured in locations of 15-55 km depth.

Model With Lower Convergence Velocity (u s = 5 cm/yr)
There is a spectacular change in the crustal flow pattern in the Himalaya-Tibet system, marked by the appearance of a welldeveloped extrusion zone in the Himalayan wedge as the collision rate is reduced to a threshold value of 5 cm/yr ( Figure 5A). At this stage, the overriding plate undergoes gravitational collapse to set in crustal flows counter to the convergence direction (∼1-1.5 cm/ yr), leading to a transition from contraction to extensional deformation at a strain rate of 2 × 10 −14 s −1 in adjacent to the wedge. On the other hand, the horizontal shortening mostly concentrates in the foreland region of the Himalayan wedge at rates of −6 × 10 −14 s −1 ( Figure 5B). We use this u s 5 cm/yr model to study the extrusion flow kinematics (discussed separately later). The crustal flow in the extrusion zone is found to occur at high velocities (3.5-4 cm/yr), which scale nearly to the subduction velocity, as also suggested by previous workers (e.g., Grujic et al., 2011). The wedge extrusion originates from a deep-crustal region (model estimate: ∼55-60 km) and meets the surface along an inclined (∼15°) trajectory ( Figure 5A).
The thermo-mechanical simulations with u s 5 cm/yr show visually contrasting isotherm structures compared to u s 8 cm/yr in the model Himalayan wedge, where the extrusion driven heat advection heavily folds the isotherms to convex upward, leading to an inverted thermal stratification ( Figure 5C), reported as  Fort, 1975;Vannay and Grasemann, 2001). In the same set of locations, as chosen for u s 8 cm/yr in the model Himalayan wedge, a quantitative analysis of the spatiotemporal changes of temperature (T) and burial/extrusion rates shows that T reverses to follow positive trends with time, where the extrusion assisted heat advection largely contributes to heating at shallow depths ( Figure 5D). For example, at a depth of 30 km, the initial temperature (300°C) increases to more than 800°C on the run time of 6 Ma ( Figure 5D). Similarly, T at 45 km and 55 km increases from ∼600 to 910°C and ∼700 to 930°C, respectively, for the same time interval since the initiation of the wedge extrusion ( Figure 5D).
The burial process is replaced by crustal extrusion (positive vertical velocity components) in those points when u s 5 cm/yr, and the extrusion rates range from 1.7 to 0.4 cm/yr ( Figure 5E). Another important model observation is that the extrusion rates, in overall, show decreasing trends with time, e.g., the highest extrusion rate drops from an initial value of ∼1.6 cm/yr to 0.6 cm/yr on a time span of 6 Ma.

Topography Development
To study how the first-order surface topography in the Himalaya-Tibet system might respond to reducing convergence rates, we calculated the vertical velocity components at the top surface of our thermo-mechanical models ( Figure 6A). The graphical plots show that the high-convergence rate gives rise to topographic upliftment mostly in the overriding plate, where the uplift rates (0.8 cm/yr) far exceed that occurring in the Himalayan wedge (0.2-0.4 cm/yr) (Figure 6Ai). But, the scenario dramatically changes following the appearance of a fully developed extrusion zone at u s 5 cm/yr when the topographic upliftment mainly concentrates in the Himalaya at 0.8-1.6 cm/ yr, leaving the adjoining part in the overriding plate to subside at rates varying from −0.7 to 0.0 cm/yr (Figure 6Aii).

Dynamic Analysis of the Wedge Extrusion
From the analysis of lubrication dynamics, presented in Conceptualization of the Himalaya-Tibet Dynamic Coupling, we showed how a decreasing convergence velocity can significantly reduce the dynamic pressure in the wedge. In our thermo-mechanical numerical models, we observe that the initiation of wedge extrusion in the Himalaya does not occur . Note, at u s 6 cm/yr, P Dy drops below the overburden pressure (P Ob ) within a small part of the wedge. But, for u s 5 cm/yr, P Dy becomes significantly less than P Ob over a large part of the wedge, resulting in a fully developed large-scale extrusion zone within the wedge.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 632806 unless the convergence rate (u s ) decreases to a threshold value ( Figure 6B). Employing the concept of flow-induced dynamic pressure, as predicted from the lubrication mechanics, we explain here the underlying dynamics of it. The analytical solution (Eq. 1) of lubrication theory, for u s 8 cm/yr the dynamic pressure in the wedge is found to be high enough (P Dy 1-5.5 kb) to overcome the overburden pressure (P Ob 1-3.5 kb) at the wedge-overriding plate interface ( Figure 6B), and thereby support the gravity load of the overriding block. The P Ob estimate takes into account both the density difference between the wedge and overriding plate, and the topographic elevation of Tibet. This dynamic condition (P Dy > P Ob ) allows the Himalayan wedge to propagate the horizontal motion into the Tibetan crust, resulting in spatially distributed contractional deformations, but without forming any extrusion zone in the Himalayan wedge ( Figure 4A). With u s reduced to 5 cm/yr, there is a significant drop in P Dy (1-2.6 kb) in the upper to middle region of the wedge (0-80 km; Figure 6B), which fails to support the overlying load, P Ob (1-2.9 kb). The unstable state (P Dy < P Ob ) eventually causes the overriding block to collapse till the stable condition, P Dy P Ob is attained in the system ( Figure 2C). The collapse forces the crustal materials in the Himalaya to extrude toward the low-pressure region along a narrow zone (discussed separately in Exhumation of Deep-Crustal Rocks), and at the same time, extensional deformations initiate in the elevated Tibetan plateau, as revealed from the horizontal flow velocity field and extensional strain rates in the models ( Figures 5A,B).

Exhumation of Deep-Crustal Rocks
In our Himalaya-Tibet collision model, the wedge extrusion originates from a depth of around 55-65 km in the wedge ( Figure 5A), and widens upward to cover a horizontal distance of ∼80 km on the surface. The extrusion zone is bounded by two crustal-scale ductile shear zones, comparable to the STDS and the HHD in the GHS ( Figure 1B) Wang et al., 2016;Carosi et al., 2018). The deeper (∼20-60 km) section of the extrusion zone runs nearly parallel to the wedge-overriding plate interface but attains an inclination of 8 to 10°while emerging to the surface. We have constructed three velocity profiles (A1, A2, and A3) across the extrusion zone in its lower (55 km), middle (35 km), and upper (20 km) segments ( Figure 7B), primarily to study the extrusion flow kinematics. The velocity profiles show nearly a parabolic (Poiseuille type) flow velocities (maximum extrusion velocity: 3-3.5 cm/yr at the axis of extrusion zone), implying that the extrusion process is driven by a pressure gradient. The bottom portions of the velocity profiles describe a Couette type of flow occurring below the wedge extrusion zone, where the down-going plate drag takes control in setting the simple shear flow kinematics. Using a Lagrangian approach, we track the material particles, located at varying depths along the extrusion zone to reconstruct the P-T-t paths of deep-crustal materials exhuming to the surface ( Figure 7C). A particle located close to the base of the extrusion zone (depth: 55 km) attains a peak metamorphic condition with a pressure and temperature of ∼17 kb and ∼820°C, respectively, and its exhumation to the surface takes a time of 5.2 Ma (Supplementary Figure S2). On the other end, material particles located at relatively shallower depths (30-40 km) show a peak P and T in the range of 8-11 kb and 500-650°C, respectively, and they are exhumed to the surface in a lesser time (2-3.5 Ma). The time scale for the exhumation of deep-crustal materials thus ranges from 2 to 5.5 Ma (Supplementary Figure  S2). From the P-T paths, it is evident that the peak P-T data of exhumed particles mostly fall in the kyanite stability field and some of them also pass through the sillimanite stability field ( Figure 7B). Our model also allows us to constrain the most active phase of the exhumation process, which occurred at 22 Ma with convergence velocity 5 cm/yr as a short-lived 4 to 6 Ma event ( Figure 5A). The wedge extrusion is expected to slow down afterward because the overriding plate collapse would progressively increase the dynamic pressure to attain a stable (P Dy P OB ) dynamic coupling between the Himalayan wedge and the overburden ( Figure 2C). For the extrusion process to remain active, the unstable condition (P Dy < P OB ) has to prevail in the Himalayan wedge, which can happen due to either a continued decrease in the Himalayan wedge viscosity or reduction of the subduction velocity further below 5 cm/yr ( Figure 2D).
Our thermomechanical model as well as lubrication theory suggest that the extrusion event in the Himalaya was accompanied by gravitational collapse of the overriding block, resulting in vertical thinning in the Tibetan part. We have estimated a temporal variation of the vertical thinning rate in the model, which indicates that the thinning rate was maximum (∼0.8 cm/yr) during the initial stage of extrusion, which gradually decreased to <0.2 cm/yr over a time span of 6 Ma ( Figure 7D). Another remarkable model finding is the transition from fast to slow exhumation process, where the rate of exhumation drops from 1.8 to 0.6 cm/yr ( Figure 7E). The model estimates are validated with the available geophysical and petrological data, which we discuss in the next section.

DISCUSSION
By combining lubrication theory and thermo-mechanical numerical simulations, the alternative tectonic model we have presented in this paper explains the extrusion mechanics of deepcrustal rocks in the Himalaya and the spatio-temporal variations of their pressure-temperature conditions. The lubrication dynamics indicates that a slow-down of the India-Asia collision rate (15 cm/yr to <5 cm/yr) appears to be the causal factor in triggering the extrusion process in the form of a southvergent extrusion zone in the Himalaya (Figures 5A, 7A), as proposed by earlier workers (Grujic et al., 1996;Vannay and Grasemann, 2001). However, unlike the previous models that account for focused surface erosion at the mountain front to show the GHS extrusion in the form of a channel (Beaumont et al., 2001), the present model suggests that the extrusion process has been entirely controlled by the coupled tectonics of the Himalayan wedge and the Tibetan plateau, where the collapse of the plateau forced the deep-crustal materials in the Himalaya to extrude along a narrow zone, without any necessary involvement of localized rheological weakening or focused surface erosion. The surface erosion probably influenced the size and shape of the surface expression of the extrusion zone ( Figure 1A), and also facilitated the removal process to exhume the extruded materials to the surface. Recently, Piccolo et al., 2018 have shown that the channel flow of hot partially melted mid-crustal rocks in the Himalaya is not sensitive to enhanced surface erosion in the mountain front, and the erosion is not a necessary condition for GHS exhumation. Moreover, several workers have reported exhumation of deep crustal materials from other types of subduction settings without any role of focused surface erosion (Brun and Faccenna, 2008;Li et al., 2011;Malusa et al., 2015;Piccolo et al., 2018). These models suggest that the extrusion of high-pressure rocks is indeed possible without any active participation of localized surface erosion in convergent boundaries. To sum up this discussion, we propose along the same line that the wedge extrusion of HP rocks in the GHS is a consequence of the collapse of the Tibetan plateau in response to the reducing India-Asia convergence velocity. The extrusion process became most active when the convergence velocity dropped down to ∼5 cm/yr at ∼22 Ma, which falls within the activity period of HHD (27-17 Ma) ( Figure 8A; Montomoli et al., 2013;Cottle et al., 2015;Khanal et al., 2015;Wang et al., 2016;Agustsson et al., 2016;Walters and Kohn, 2017;Carosi et al., 2018).

Wedge Extrusion Geometry and Kinematics
In the proposed Himalaya-Tibet model the wedge extrusion originates at a depth of around 50-60 km, and appears as a downward tapering narrow zone, bounded by an extensional shear zone (STDS) in north and a pair of thrust type shear zones (HHD and MCT) on south ( Figure 8B). Available seismic sections of the Himalaya-Tibet system also show a similar tapered-wedge geometry for the GHS, traceable down to a depth of at least 55 km; it possibly goes beyond (Hauck et al., 1998;Jiménez-Munt et al., 2008;Godin and Harris, 2014;Clift, 2017 crustal materials on the surface over a stretch of ∼80 km ( Figure 7A). Composite geological maps of the Himalaya show a large along-strike variation of the exposed width (25-120 km) of the GHS belt (Kellett et al., 2010;Spencer et al., 2012;Cooper et al., 2013). However, several lithospheric scale cross-sections of the Himalayan wedge provide a persistent thickness of the GHS to be in average 40 km (e.g., Godin and Harris, 2014) ( Figure 1B), in good agreement with our model estimate ( Figure 7A). But, our 2D model does not account for the lateral variations of the GHS width, mentioned above ( Figure 1A; Carosi et al., 2018), which we infer to be the influence of some higher-order geological factors, such as local scale faulting, folding, and up-doming of gneissic rocks (Beaumont et al., 2001;Kellett et al., 2010). In addition, according to Piccolo et al. (2018) the varying upper crustal strength can significantly influence the exhumation patterns of metamorphic rocks in orogen. All these complexities in the tectonic processes in the brittle upper crust (top 5 km), coupled with surface erosion can result in strong spatial heterogeneity in the surface manifestation of the extrusion zone.
Our thermo-mechanical model shows Poiseuille-type l flows within the extrusion zone, where the velocity magnitude is maximum (3.4 cm/yr) along the axis of the extrusion zone, decreasing to its upper and bottom boundaries. There is a narrow zone of Couette flow below the extrusion zone. By combining the two types of crustal flow, we obtain a hybrid flow in the model ( Figure 7B). Interestingly, several workers have proposed from the kinematic analyses of GHS rocks a similar hybrid flow model, dominated by the Poiseuille type (pipe-flow effect) (Grujic et al., 1996;Vannay and Grasemann, 2001;Grujic et al., 2002). The hybrid flow kinematics causes the isotherms to overturn, and ultimately gives rise to inverted metamorphic isograds ( Figure 5C), a phenomenon extensively studied over many decades (Le Fort, 1975;Kohn, 2014 and references therein).

Initiation Time of the Wedge Extrusion
The timing of rapid GHS exhumation in the Himalaya has been extensively worked out; most of the data, available in the literature, point to early Miocene (23-21 Ma) (Hodges et al., 1992;Grujic et al., 1996;Catlos et al., 2002;White et al., 2002;Daniel et al., 2003;Najman et al., 2009). For example, using monazite geochronology of the GHS rocks from the Annapurna region of Central Nepal, Corrie and Kohn (2011) have shown the attainment of their peak metamorphic conditions at 25 to 23 Ma. The Bhutan Himalaya also records the peak temperatures dating back to 23 Ma, which continued at least to 16 Ma (Kellett et al., 2010). A similar metamorphic history (26-23 Ma) is established from the Sikkim Himalaya (Sorcar et al., 2014). However, many recent P-T-t and geochronological data (U-Th-Pb on monazite) demonstrated that exhumation in the GHS started as late as 28-27 Ma  or even older in Central Nepal (41 Ma in the uppermost GHS; Carosi et al., 2016). All these data presented above converge to the fact that the main phase of exhumation in the GHS occurred in between 22 to 28-27 Ma ( Figure 8A). In earlier tectonic models the two bounding shear zones, the MCT and the STDS were considered to be synchronously active to show the exhumation of GHS rocks at 22 Ma (Burchfiel and Royden, 1985;Hodges et al., 1992, Beaumont et al., 2001. But, many recent studies seriously questioned their contemporaneous relation (Kohn, 2008;Carosi et al., 2013;Finch et al.. 2014;Weinberg, 2016;Iaccarino et al., 2017Montemagni et al., 2019;Montemagni et al., 2020). It has been shown that the STDS initiated before 24-25 Ma ( Figure 8A), as reported from western Nepal Soucy La Roche et al. 2018;Kellett et al., 2019), whereas the activity of MCT took place in a much later stage (18-17 Ma) Wang et al., 2016;Montemagni et al., 2019). It seems more probable that the HHD, which is reported to have been active from 28 to 27 Ma to 18-17 Ma Wang et al., 2016;Carosi et al., 2018) (Figure 8A), acted as a southern shear boundary of the GHS exhumation. From our numerical simulations we have shown that a fully developed wedge-extrusion zone appeared at 22 Ma when the convergence velocity reduced to 5 cm/yr (Figure 5Ai). We ran simulations with higher convergence velocities (discussed in Supplementary) to constrain the timing of extrusion initiation. Models at 6 to 5.5 cm/yr show localization of an incipient extrusion zone (Supplementary Figure S3; Figure 8B), which clearly represents a juvenile stage of the fully developed extrusion zone discussed above. It is to be noted that the Indian-Asia convergence took place at the rates of 6 to 5.5 cm/yr during 28-27 Ma ( Figure 3E). Our model thus predicts the GHS exhumation process (onset of retrograde metamorphism) to have initiated at around 27-22 Ma, and the associated P-T-t path, as reported in many recent works (Leloup et al., 2015;Zeiger et al., 2015;Wang et al., 2016;Iaccarino et al. 2017) ( Figure 8A). The wedge extrusion evolved to attain a fully developed state at 22 Ma, resulting in the most active phase of rapid exhumation process, which subsequently slowed down to a large extent at 18 Ma ( Figures 5A, 7A). The present model results also allow us to interpret that the extrusion zone in this time interval (28-18 Ma) was bounded by two crustal scale shear zones: STDS (extensional shear) in the north and HHD (thrust shear) in the south (Wang et al., 2016;Carosi et al., 2018). After 18 Ma, the southern shear boundary the wedge extrusion zone migrated southward, eventually resulting in initiation of the MCT after the most active phase of the GHS exhumation event ( Figure 8B). This late initiation of MCT in the GHS evolution has been extensively reported in the recent publications Carosi et al., 2016Wang et al., 2016;Montemagni et al., 2018) ( Figure 8A).
The present lubrication model also establishes the extrusion zone formation in the Himalaya as a crustal process coupled with the gravity collapse and extensional tectonics in the Tibetan plateau ( Figures 7A, D). Geological evidence suggests that most of the active thrusts in the central Tibet became dormant, implying a complete cessation of contractional deformations in Tibet by ∼23 Ma (Wang et al., 2014;Li et al., 2015). After 23 Ma, there was a dramatic tectonic inversion from contraction to extension, as documented from the widespread occurrence of extensional faults in the Tibetan upper crust (Wang et al., 2014;Ge et al., 2015). On the other hand, the activity of fully developed extrusion zone in the Himalayan wedge is found to prevail at ∼22 Ma, as discussed above. These two lines of observations imply a synchronous relation between the two tectonic processes, as seen in our thermo-mechanical model simulations ( Figure 7A). Previous studies also inferred the collapse of Tibet as a driving factor for the southward migration of contractional deformation from the plateau interior toward the Himalayan wedge (Wang et al., 2014;Li et al., 2015). Based on the present-day data in Southern Tibet, Ge et al. (2015) estimated the thinning strain rate for the entire thick Tibetan crust (70 km) to be around 8 nanostrain/yr, assuming that the crustal thinning began at ∼15 Ma. This estimate corresponds to a vertical thinning rate of 0.56 mm/yr since 15 Ma. According to our model estimates, the overriding Tibetan plateau underwent vertical thinning at ∼8 mm/yr rates at the initial stage (∼22-Ma) of its collapse, but weakened gradually to 0.5 mm/yr after 6 Ma ( Figure 7D), in good agreement with the estimation of Ge et al. 2015.
We also provide an estimate for the exhumation velocity (vertical extrusion rate) within the extrusion zone. The calculated exhumation velocity is found to be ∼18 mm/yr at ∼22 Ma and slowed down to ∼6 mm/yr within a time span of 6 Ma ( Figure 7E). Petrological calculations of GHS rocks also suggest for the fast exhumation (>10 mm/yr) event during the period from 22 to 17 Ma (Ganguly et al. 2000;White et al., 2002;Najman et al., 2009;Thiede et al., 2009;Grujic et al., 2011;Carrapa et al., 2014;Clift, 2017). Although the actual magnitude of calculated exhumation rates varies to some extent ( Figure 7E), they broadly agree upon a point that the exhumation in the Himalaya was rapid at ∼22 Ma, and then slowed down significantly after 17 Ma (i.e., the time of MCT initiation). Interestingly, this transition coincides with the timing of vorticity decrease during the southward migration of the deformation front, as estimated by Montemagni et al. (2020). The Himalayan wedge subsequently underwent a major reorganization of its tectonics, leading to a shift of the deformation front to the south of the MCT (White et al., 2002;Wang et al.2014).

P-T Estimates for Exhumed Metamorphic Rocks
Looking at the spatial dispositions of metamorphic grades in the Himalaya, one can find the GHS as a high-grade metamorphic belt, sandwiched between two low-grade metamorphic terrains: the Lower Himalayan Sequence (LHS) on the south and the Tethyan Himalayan Sequence (THS) on the north. In the LHS-GHS-THS package, the grade of metamorphism shows an asymmetrical spatial variation, with the highest grade located closer to the MCT (i.e. GHS-LHS boundary) (Metcalfe, 1993;Corrie and Kohn, 2011;Kohn, 2014). We validate our model results using this asymmetrical pattern. The model-calculated pressure (P) and temperature (T) are plotted as a function of horizontal distance across the channel and compared with available petrological data from the exhumed high-pressure rocks ( Figures 9A,B) in the GHS. The plots are strikingly asymmetrical; the temperature steeply increases across the MCT, (from T 550-600°C in LHS to T 700-820°C in the GHS) to attain a peak temperature (T 820°C) at a location (20 km from the MCT) in the GHS, and it then decreases at a gentle gradient toward the north (T 650-700°C at the STDS) ( Figure 9A). The pressure plot for exhumed materials shows a similar asymmetrical trend: P 6-8 kb in the LHS, which reaches a peak value of 10-15 kb within the GHS, then dropping to around 7-8 kb near STDS ( Figure 9B). It is interesting to note that the location of peak P coincides with that of peak T. We have compared our model P-T profiles with the petrological data reported for the LHS-GHS-THS litho-packet from several sections in the Himalaya, such as eastern Garhwal, central Nepal, Langtang region in Bhutan and Sikkim ( Figures 9A,B) (Kohn, 2008;Kellett et al., 2010;Spencer et al., 2012;Sorcar et al., 2014). Our model estimates tally well with the petrological calculations, where both pressure and temperature of exhumed HP rocks near the MCT region (P 8-15 kb and T 700-800°C) grossly higher than that adjacent to STDS (P 5-10 kb and T 600-750°C). This P-T variation across the GHS can be explained from the kinematics of wedge extrusion in the present model. Above the MCT, the extrusion zone exhumes crustal materials from relatively deeper depths than below the STDS ( Figure 7A).
In our model, the peak P-T data of extruded HP rocks mostly fall in the kyanite stability field, and few in the sillimanite stability field ( Figure 9C). This finding is supported by the petrological data reported from different regions of GHS in the Himalaya. Most of the petrological peak P-T data, compiled from the northwestern Himalaya, the central Himalaya (Kumaun to Langtang), and Bhutan track the kyanite stability field (Catlos et al., 2001;Daniel et al., 2003;Kohn, 2008;Corrie and Kohn, 2011;Spencer et al., 2012;Kohn, 2014). But some rocks in the Nyalum to Sikkim Himalaya, including structurally highest rocks in the Langtang region plot well within the sillimanite stability field (Kohn, 2014;Sorcar et al., 2014). The corresponding model P-T paths, on average, fall in the line of the petrological calculations (Catlos et al., 2001;Kohn, 2014;Iaccarino et al., 2015;Wang et al., 2016;Iaccarino et al., 2017;Waters, 2019) (Figure 9C). Our model peak P-T data provide an interpretation for the occurrence of eclogites (P > 1.5 GPa) in the central and eastern Himalaya (Lombardo and Rolfo, 2000;Corrie et al., 2010). However, the present study does not account for the western Himalayan eclogites (coesite bearing) of the Kaghan valley and Tso Morari nappe, which records a different exhumation history with the peak pressures, going up to When u s ≥ 6 cm/yr (before 27 Ma) there was major contractional deformation in the Tibetan plateau, which resulted in an increase of Tibetan plateau elevation. At u s ≤ 5 cm/yr (after 22 Ma) the plateau started to collapse, while Himalaya began to rise at its southern margin.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 632806 3 GPa (O'Brien et al., 2001;O'Brien, 2019). The reason is that our modeling is based on a low-dip subduction setting for the Indian lithosphere, as applicable to the India-Asia collision system since 40-45 Ma (Chemenda et al., 2000;Jiménez-Munt et al., 2008;Gloaguen and Ratschbacher, 2011), and it cannot justify the exhumation of crustal rocks from depths ≥100 km. Before 40-45 Ma, the oceanic Tethyan lithosphere most probably had a steeper subduction angle ( Figure 10A). Such steep subduction of the Tethyan lithosphere at high velocity (>15 cm/yr) ( Figures  4A,E) could drag crustal materials to greater depths (>100 km). A subduction velocity drop to 10 cm/yr during the oceanic to continental collisional tectonic transition ( Figure 10B) might favor an exhumation of deep crustal rocks, as predicted from lubrication theory ( Figure 2D) . We extrapolate our model results to suggest that a sharp velocity drop during 45-50 Ma ( Figure 9D) could facilitate the exhumation of coesite bearing eclogites (P ≈ 3 GPa) to shallower levels (40 km) in the western Himalaya (Kohn and Parkinson, 2002). During the second event (around 22 Ma) of significant subduction velocity drop (to 5 cm/yr; Figure 9D) and consequent thermal heating of the Himalayan wedge, these rocks were exhumed to the surface ( Figure 10C).

Active Phase of the Wedge Extrusion
The lubrication dynamics allows us to decipher that the collapse of the overriding plate would progressively weaken on a time span of 4-6 Ma, as discussed earlier, eventually causing the extrusion process to slow down for the same convergence velocity of 5 cm/ yr ( Figure 5A). This time-dependent weakening of the wedge extrusion, predicted from our theoretical calculations as well as numerical model, is indeed reflected in the geological records that indicate the rapid phase of GHS exhumation to phase out at around 17 Ma (Ganguly et al., 2000;White et al., 2002;Najman et al., 2009;Clift, 2017). For example, using Sr-Nd isotope analysis of sediments and 40 Ar→ 39 Ar ages of detrital micas, White et al. 2002 suggested that the rocks in lower Dharamsala groups (21-17 Ma) were derived from a rapidly exhumed high-grade metamorphic slab of GHS, whereas the Upper Dharamsala groups (17-13 Ma) of rocks were from low-grade metamorphic rocks. Based on these geological evidence, they postulated southward shifting of the active deformation zone, following weakening of the rapid exhumation process in the GHS after 17 Ma Similarly, Najman et al. (2009) showed exceptionally fast exhumation of high-grade rocks in the GHS during 21 to 17 Ma, witnessing nearly zero lag time between the cooling of exhumed rocks in GHS and their subsequent deposition in the mountain front. On the contrary, the lag time increased substantially, by around 6 Ma during the period of 17 to12 Ma, pointing to a remarkable slowdown of the exhumation process in the Himalaya after 17 Ma. The lag time then stretched more, implying a further slowdown in the exhumation process after 12 Ma (Clift, 2017). Our model suggests that the wedge extrusion, after its fully developed form at ∼22 Ma, remained active to sustain deep-crustal exhumation at high rates (>1 cm/yr) for a time duration of 4-6 Ma. After this active phase, it started to weaken as the system tended to restore its stable state. The wedge extrusion could continue to operate steadily, had the convergence velocity dropped at higher rates. It is to note here, one can see in the India-Asia collision history that the convergence velocity has dropped from ∼5 cm/yr at ∼22 Ma to ∼3.5 cm/yr at around 10 Ma (Molnar and stock, 2009;Copley et al., 2010); since then the change in the collision rate occurred at a much slow pace ( Figure 9D). Based on the lubrication dynamics, we thus suggest that the most active phase of exhumation lasted till 17 Ma, and then it started to slow down with southward migration of the deformation front and initiation of the MCT. Finally, we combine the geological history with our thermo-mechanical model results to show the evolution of the Himalaya-Tibetan tectonics in three stages in Figure 10.

CONCLUSION
From the thermo-mechanical model results, supported largely with petrological, seismic, and field data from the Himalaya-Tibet system, we conclude the following. 1) High-pressure rocks in the Greater Himalaya Sequence underwent extrusion from a deeper region in the rheologically coherent ductile Himalayan wedge, rather than from a melt weakened mid-crustal layer beneath Tibet, proposed in earlier models.
2) The wedge extrusion was driven primarily by the gravity collapse of the overriding Tibet and it became most active when the Indian-Asia convergence velocity to ∼5 cm/yr at ∼22 Ma. According to our model, the erosion does not seem to be the driving factor for deep-crustal extrusion, although it might have controlled the ultimate surface manifestation of the extruded GHS rocks.
3) The wedge-extrusion zone in its most active phase (22-17 Ma) had a downward tapered geometry, bounded by two ductile-shear zones: an extensional (STDS) zone in the north and a thrust (HHD) in the south, which subsequently migrated southward to form MCT at 17 Ma. 4) Our thermo-mechanical model successfully explains the inverted metamorphic isograds observed in the GHS. The inversion is a consequence of the hybrid flow kinematics, dominated by a Poiseuille flow within the extrusion zone and a Couette flow at its base. 5) The estimated peak pressures (10-17 kb) and temperatures (700-820°C) of exhumed GHS rocks and their asymmetric pattern, characterized by P-T peaks close to the southern boundary (above MCT) are consistent with petrologic data. 6) The lubrication dynamics presented here accounts for the fast exhumation (>1 cm/yr) event  in the GHS, which falls within the activity period of HHD (28-17 Ma) before the initiation of MCT, as reported in the literature. Finally, our model suggests that the extrusion process in the GHS slowed down further after 10 Ma as the India-Asia convergence velocity did not significantly drop anymore.

DATA AVAILABILITY STATEMENT
All the necessary information to reproduce the results are presented in the article/Supplementary Material and listed in Table 1. Petrological data used for comparison of model results are from Spencer et al. (2012), Khon (2008), Khon (2011), Sorcar et al. (2014), Kellett et al. (2010), and Daniel et al. (2003). Further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

FUNDING
This work has been supported by the DST-SERB through JC. Bose fellowship (SR/S2/JCB-36/2012) to NM and CSIR Senior Research Fellowship to GM.

ACKNOWLEDGMENTS
We thank two reviewers and Editor Guillermo Booth-Rea for their excellent comments and constructive suggestions to improve the interpretation and presentation of our model results. Our study has greatly benefitted from their valuable discussions. Commercial finite element code COMSOL Multiphysics (2015) has been used to conduct the numerical simulations presented in this article.

SUPPLEMENTARY MATERIAL
The