Channel ﬂow, tectonic overpressure, and exhumation of high-pressure rocks in the Greater Himalayas

. The Himalayas are the archetype of continental collision, where a number of long-standing fundamental problems persist in the Greater Himalayan Sequence (GHS): (1) contemporaneous reverse and normal faulting, (2) inversion of metamorphic grade, (3) origin of high- (HP) and ultrahigh-pressure (UHP) rocks, (4) mode of ductile extrusion and exhumation of HP and UHP rocks close to the GHS hanging wall, (5) ﬂow kinematics in the subduction channel, and (6) tectonic overpressure, here deﬁned as TOP = P/P L where P is total (dynamic) pressure and P L is lithostatic pressure. In this study we couple Himalayan geodynamics to numerical simulations to show how one single model, upward-tapering channel (UTC) ﬂow, can be used to ﬁnd a uniﬁed explanation for the evidence. The UTC simulates a ﬂat-ramp geometry of the main underthrust faults, as proposed for many sections across the Himalayan continental subduction. Based on the current knowledge of the Himalayan subduction channel geometry and geologi-cal/geophysical data, the simulations predict that a UTC can be responsible for high TOP ( > 2). TOP increases expo-nentially with a decrease in UTC mouth width, and with an increase in underthrusting velocity and channel viscosity. The highest overpressure occurs at depths < − 60 km, which, combined with the ﬂow conﬁguration in the UTC, forces HP and UHP rocks to exhume along the channel’s hanging wall, as in the Himalayas. By matching the com-puted velocities and pressures with geological data, we constrain the GHS viscosity to be ≤ 10 21 Pa s, and the effective convergence (transpression) to a value ≤ 10 %. Variations in channel dip over time may promote or inhibit exhumation ( > or < 15 ◦ , respectively). Viscous deformable walls do not affect overpressure signiﬁcantly enough for a viscosity contrast (viscosity walls to viscosity channel) of the order of 1000 or 100. TOP in a UTC, however, is only possible if the condition at the bottom boundary is no-outlet pressure; otherwise it behaves as a leaking boundary that cannot retain dynamic pressure. However, the cold, thick, and strong lithospheres forming the Indian and Eurasian plates are a good argument against a leaking bottom boundary in a ﬂat-ramp geometry, and therefore it is possible for overpressure to reach high values in the GHS.

Abstract. The Himalayas are the archetype of continental collision, where a number of long-standing fundamental problems persist in the Greater Himalayan Sequence (GHS): (1) contemporaneous reverse and normal faulting, (2) inversion of metamorphic grade, (3) origin of high-(HP) and ultrahigh-pressure (UHP) rocks, (4) mode of ductile extrusion and exhumation of HP and UHP rocks close to the GHS hanging wall, (5) flow kinematics in the subduction channel, and (6) tectonic overpressure, here defined as TOP = P /P L where P is total (dynamic) pressure and P L is lithostatic pressure. In this study we couple Himalayan geodynamics to numerical simulations to show how one single model, upward-tapering channel (UTC) flow, can be used to find a unified explanation for the evidence. The UTC simulates a flat-ramp geometry of the main underthrust faults, as proposed for many sections across the Himalayan continental subduction. Based on the current knowledge of the Himalayan subduction channel geometry and geological/geophysical data, the simulations predict that a UTC can be responsible for high TOP (> 2). TOP increases exponentially with a decrease in UTC mouth width, and with an increase in underthrusting velocity and channel viscosity. The highest overpressure occurs at depths < −60 km, which, combined with the flow configuration in the UTC, forces HP and UHP rocks to exhume along the channel's hanging wall, as in the Himalayas. By matching the computed velocities and pressures with geological data, we constrain the GHS viscosity to be ≤ 10 21 Pa s, and the effective convergence (transpression) to a value ≤ 10 %. Variations in channel dip over time may promote or inhibit exhumation (> or < 15 • , respectively). Viscous deformable walls do not affect overpressure significantly enough for a viscosity contrast (viscosity walls to viscosity channel) of the order of 1000 or 100. TOP in a UTC, however, is only possible if the condition at the bottom boundary is no-outlet pressure; otherwise it behaves as a leaking boundary that cannot retain dynamic pressure. However, the cold, thick, and strong lithospheres forming the Indian and Eurasian plates are a good argument against a leaking bottom boundary in a flat-ramp geometry, and therefore it is possible for overpressure to reach high values in the GHS.  Grujic et al., 2011 andUnsworth et al., 2005). White line along 90 • E marks the cross section shown in (b). (b) Schematic section across the Himalayas (adapted from Grujic et al., 2011), in which the UTC stands out (GHS in red). The GHS is bounded at the top by the South Tibet Detachment (STD) and at the bottom by the Main Central Thrust (MCT). MHT is the Main Himalayan Thrust. (c) Model set-up of the UTC, with shape and dimensions similar to the natural prototype in (b). The "footwall" (moving wall) and the "hanging wall" (no-slip wall) correspond to the MCT and the STD, respectively. Apart from the later folding of both MCT and STD, the similarity between nature and model set-up is apparent.

Geological setting
Based on metamorphic grade and structural style, four units and the major faults separating them were distinguished by Gansser (1964). These are from bottom to top in Fig. 1: Sub-Himalayan Sequence (SHS -unmetamorphosed Tertiary rocks), Main Boundary Thrust (MBT), Lesser Himalayan Sequence (LHS -low-grade metamorphic rocks), Main Central Thrust (MCT), Greater Himalayan Sequence (GHShigh-grade metamorphic rocks), South Tibetan Detachment (STD), and Tethyan Sedimentary Sequence (TSS -unmetamorphosed to weakly metamorphosed rocks). All the main faults are north-dipping thrusts, except the STD that also dips to the north but is a normal fault.
The GHS shows patchy occurrences of eclogites close to the STD Ganguly et al., 2000;O'Brien et al., 2001;Groppo et al., 2007;Corrie et al., 2010;Kellett et al., 2013;Sorcar et al., 2014;Zhang et al., 2015;Fig. 1a). Recent petrologic studies provide estimates for spatiotemporal variations in pressure (P ) and temperature (T ) in the GHS. The peak metamorphic conditions are T ∼ 760 • C and P ≥ 1.5 GPa for eclogitization in the Bhutan Himalayas . Peak conditions with T = 670 • C and P ≥ 1.5 GPa were reported for the Nepal Himalayas (Corrie et al., 2010). On the other hand, an estimate of the metamorphic peak at P = 2.7-2.9 GPa and T = 690-750 • C from coesite-bearing eclogites in the western Himalayas was provided by O'Brien et al. (2001). The eclogites have been, in part, overprinted by regionally more extensive granulite facies conditions of 800 • C at ∼ 1 GPa Ganguly et al., 2000;Groppo et al., 2007;Zhang et al., 2015). PT-time paths suggest exhumation of these high-grade rocks under nearly isothermal decompression after peak metamorphic conditions (Ganguly et al., 2000;Groppo et al., 2007;Sorcar et al., 2014). Using cooling rates, the exhumation history of the high-grade rocks was interpreted as a two-stage event by Ganguly et al. (2000), marked by exhumation at a rate of 15 mm yr −1 to a depth of 15 km, followed by slow exhumation at a rate of 2 mm yr −1 to a depth of at least 5 km, which occurred broadly in Miocene times Corrie et al., 2010;Kellett et al., 2013;Sorcar et al., 2014;Warren et al., 2011;Rubatto et al., 2013).
The exhumation mechanics of GHS rocks is one of the most debated issues in the Himalayas (and elsewhere where HP and UHP rocks outcrop), having led to a variety of tectonic models that postulate channel flow by topographic forcing (Wobus et al., 2005;Beaumont et al., 2001) or transpression (Grujic et al., 1996). Grujic et al. (1996) first proposed the GHS in the Bhutan Himalayas as deep crustal ductile rocks extruded between the MCT and the STD. Numerical models have integrated geological, tectonic, geophysical, metamorphic, and rheological data to provide possible explanations for the exhumation process. The models postulate a channel flow of low-viscosity rocks in the middle to lower crust, driven by a topographic pressure gradient, to account for the extrusion dynamics of high-grade metamorphic rocks in the GHS (Wobus et al., 2005;Beaumont et al., 2001). The channel flow model can also explain the coeval reverse and normal kinematics along the MCT and STD, respectively (Fig. 1b). However, as Grujic et al. (2011) pointed out, these models cannot "predict the exhumation of lower orogenic (> 50 km, i.e. > 1.4 GPa) crustal material" in their basic form. To overcome this limitation, an alternative exhumation mechanism was proposed by Grujic et al. (2011), with additional tectonic forcing (transpression) by the impingement of strong Indian crust into the already weak lower crustal granulitized eclogites below southern Tibet. However, previous models do not comprehensively address the mechanics of overpressure leading to the formation of eclogites (Schulte-Pelkum et al., 2005), and their focused exhumation close to the STD.
Given that the current models do not fully explain the observations in the GHS, in this study we couple eastern Himalayan geodynamics with numerical simulations to show how one single model, upward-tapering channel (UTC) flow, as in the current eastern Himalayas (Fig. 1b), can be used to find a unified explanation for the following persisting problems: (1) contemporaneous reverse and normal faulting, (2) inversion of metamorphic grade, (3) origin of HP and UHP rocks, (4) mode of ductile extrusion and exhumation of HP and UHP rocks close to the GHS hanging wall (STD), (5) flow kinematics in the subduction channel, and (6) tectonic overpressure.

Premises
We model channel flow with a linear viscous fluid by the Navier-Stokes equation with body force (gravity); therefore, pressure in the channel depends on the viscosity and velocity configuration. Most critically, the velocity field depends on channel geometry and conditions applied at the boundaries (e.g. Marques et al., 2018). Ultimately, tectonic overpressure (TOP) can only exist if the channel walls are strong enough. Therefore, when investigating pressure in a viscous channel, one has to take into account four fundamental issues: -Viscosity. The viscosity term in the Navier-Stokes equation depends on a number of parameters, all of which are incorporated in the Arrhenius term in a constitutive equation. Therefore, the modeller has two options when investigating the effects of viscosity on pressure: either use a full constitutive equation and test all the parameters in the Arrhenius term or simply and directly vary the magnitude of the viscosity. We chose the second option in our numerical simulations, since our focus is the assessment of parameter variations on the development of overpressure and flow configuration.
-Geometry of the channel. Given that flow configuration inside the channel plays a critical role in the pressure distribution, we tested three main shapes of the channel: parallel-sided (parallelepiped), and upward-(similarly to Marques et al., 2018) or downward-tapering channels.
-Boundary conditions. The conditions at the boundaries can either promote or inhibit TOP, because they control the flow pattern and the pressure retention inside the channel. Therefore, we tested different velocity configurations applied at the underthrusting (foot) wall (simple or simple + pure shears), and different conditions at the boundaries like slip, no-slip, or outlet pressure.
-How the walls of the pressure vessel react to internal pressure. Under particular applied boundary conditions, the Navier-Stokes equation produces TOP in an upward-tapering channel that can reach values orders of magnitude greater than observed in nature; therefore we will discuss the theoretical values in view of the current knowledge on natural HP and UHP rocks. The discussion of channel flow is similar to discussing a pressure vessel with an overpressured fluid inside: one has to investigate the conditions to produce overpressure inside the vessel (the channel in the prototype), and simultaneously the strength of the vessel walls (the lithosphere in the prototype) to support the internal pressure without failure (by brittle or viscous yield). We will therefore discuss the strength of the channel walls in view of the current knowledge about the Indian (footwall) and Eurasian (hanging wall) lithospheres, especially in terms of thickness and strength.
This study builds on the conceptual work by Marques et al. (2018) on tectonic overpressure, in which the main conclusions are that TOP depends critically on boundary conditions (e.g. an upward-tapering channel can produce large TOP, whereas an outlet condition at the bottom prevents TOP from developing) and on critical parameters like strain rate and viscosity.
Given the above premises, we investigated the conditions under which overpressured rocks can form and be exhumed in a prototype like the Himalayas: geometry of the channel, conditions at the boundaries, applied velocities, and viscosity. Based on the numerical simulations and the current knowledge of the Himalayas, we discuss the theoretical values of overpressure, the obtained exhumation velocities, the most likely viscosity of the subducted rocks, and finally the effects of the strength of the channel walls on overpressure.

Numerical modelling
We modelled the subduction channel, as illustrated in Fig. 1c, with an incompressible linear viscous fluid. The assumptions of incompressibility and linearity considerably simplify the model, and constitute standard procedure in many geophysical and geodynamic problems (cf., e.g. Ranalli, 1995;Turcotte and Schubert, 2014). The set-up simulates a flat-ramp geometry of the main underthrust faults, as shown in many cross sections of the Himalayas, in particular the one shown in Fig. 1b. For steady-state flow of a viscous incompressible Newtonian fluid at very low Reynolds number, the dynamic Navier-Stokes equations reduce to the Stokes approximation, which is the basis of the COMSOL code for computational fluid dynamics used here (COMSOL, 2018).

Boundary conditions and model set-up
The boundary conditions were as follows (see Fig. 1c, and Methods in the Appendix for further details): (1) slab-parallel velocity (U ) applied on the underthrusting (foot)wall (2 to 20 cm yr −1 ; Feldl and Bilham, 2006;DeMets et al., 2010), and fixed hanging wall; (2) viscosity (η) between 10 19 and 10 22 Pa s (Beaumont et al., 2001;England and Houseman, 1989;Copley et al., 2011); (3) channel dip α (15-30 • ); (4) channel mouth width W m = 25 to 100 km, and width at the channel's base W b = 150 or 200 km, from which we define W * m = W m /W b ; (5) constant density of the material in the channel (2800 kg m −3 ). Given the viscosity contrast between the footwall and hanging wall of the GHS and channel material, the channel walls were assumed undeformable in the first simulations, except when testing the effects of nonrigid walls on overpressure.
The metamorphic processes occur in response to the total isotropic stress, called dynamic pressure, which is a sum of the tectonic (Stokes) and lithostatic pressures (ρgz, where ρ is density, g is gravitational acceleration, and z is depth). We evaluate the dynamic pressure to explain the occurrence of high-pressure rocks in the GHS, and we define an overpressure factor (TOP) as the non-dimensional ratio between dynamic and lithostatic pressures (Figs. A1 and A2). For a better understanding of overpressure in a UTC, we carried out a parametric study of TOP as a function of η, W m , α, U , and effective convergence velocity (transpression; see Methods in the Appendix for details). The prime focus of our investigation concerned the simulations with U = 5 cm yr −1 , α = 20 • , W m = 100 km, and W b = 150 km, which represent the most common and conservative values. We then use the numerical results to constrain the viscosity, pressure, and velocity in the channel, consistent with current geological data and estimates.

Flow patterns
The model UTC shows two main layers, one flowing downward due to applied underthrusting motion in the footwall, and another flowing upward and so inducing relative normal faulting on the hanging wall ( acts as an internal large-scale shear zone with curved geometry and thrust motion. The upward flowing layer shows, at shallow depth, a maximum velocity ≈ 0.5 × 10 −9 m s −1 , i.e. ∼ 16 mm yr −1 . The line of flow convergence separates crustal materials of contrasting pressures: one towards the footwall with P < 1.5 GPa the other towards the hanging wall with P > 1.5 GPa (red curve in left panel in Fig. 2d), which is the pressure at which eclogite formation is possible at −30 km. Overall, the flow pattern shows that significantly overpressured rocks (TOP > 2.) can be exhumed rapidly through a narrow region close to the hanging wall of the channel, which corresponds to the STD in the Himalayas and where HP and UHP rocks have been found.

Dynamic pressure and overpressure
Model results are presented as colour maps (Fig. 2) and graphs (Fig. 3), the latter showing the effects of several parameters on overpressure in the subduction channel.
Varying W m with other parameters remaining constant and η = 10 21 Pa s shows that the UTC develops overpressure in the entire range of W m /W b = W * m = 25/150 to 100/150 km (Fig. 3a). TOP is inversely proportional to W * m , and can be as high as 10 for W * m = 0.17 at depths between 20 and 60 km, with the highest TOP at 20 km depth.
TOP is sensitive to α in a UTC under a given set of values for W m , η and U (Fig. 3b). The results plotted in Fig. 3b show TOP > 1 for 15 • < α ≤ 30 • . TOP is maximal at α = 20-25 • , reaching 1.7 at depths between 40 and 60 km.
The plot in Fig. 3c shows increase in TOP with increase in U , from TOP ≈ 1.5 at U = 2-5 cm yr −1 (current Indian velocity), to TOP ≈ 11 when U = 20 cm yr −1 (Indian velocity at 60-70 Ma).
The simulations show a near-exponential variation in TOP with η ( Fig. 3d), which we use to constrain the viscosity in the Himalayan collision zone. Given that the code is based on the Stokes equation, viscosity and velocity play a fundamental role on the development of TOP. However, the flow configuration is also critical, because velocity depends on the divergence of the velocity gradient in Stokes' equation. Furthermore, the flow configuration also critically depends on the boundary conditions, therefore some conditions favour the development of TOP (e.g. narrow channel mouth in an upward-tapering channel) and others prevent it (e.g. an outlet condition at the bottom boundary).
Above we presented numerical simulations for η = 10 21 Pa s, typically applicable to the Himalayan tectonic setting. However, we ran additional simulations with different viscosities, and a set of results is presented for a viscosity of 10 22 Pa s (Fig. 4); η = 10 22 Pa s induces much higher overpressure, especially when the mouth width decreases, and when the underthrusting velocity increases to velocities that have been estimated to exist at 60-70 Ma.
Taken together, the results shown in Figs. 3 and 4 place constraints on the factors affecting overpressure. Extremely high values of TOP are obtained for η > 10 21 Pa s, U > 5 cm yr −1 , and W * m < 0.50. Varying channel dip (α) involves significant changes in the flow pattern, as shown in Fig. 5. For α = 15 • , the channel is dominated by downward flow, setting in a large-scale vortex in the deeper level, and does not show conspicuous zones of ductile extrusion, which only occurs when α > 15 • .
Besides the results obtained for a channel base width of 150 km, and variable mouth width, we also evaluated the effects of the channel base width on flow patterns and pressure distribution, by running a set of numerical simulations with a base width of 200 km. The channel flow shows similar patterns in the two cases, and small variations in pressure.

Effects of transpression on overpressure and flow
We ran a set of simulations to investigate how much a transpressional movement across the viscous channel might influ- ence the magnitude of tectonic overpressure and, especially, velocity at the channel mouth (extrusion velocity). Transpression in the numerical models was introduced by setting the magnitude of horizontal velocity in excess of that corresponding to the underthrusting movement, i.e. transpression was set by adding an extra horizontal velocity component that made the velocity vector less steep than the moving subduction footwall. Figure 6 shows a plot of TOP as a function of transpression, represented as the ratio between horizontal velocity and non-transpressional horizontal component (ca. 1.49 × 10 −9 m s −1 ). The numerical results indicate that (1) transpression has appreciable effects on overpressure, especially if transpression is large (> 20 %) and (2) trans- pression has great effects on extrusion velocity, as shown in Figs. 6 and 7.

Viscous deformable walls
We used a similar modelling approach to evaluate the magnitude of overpressure in subduction channels confined by deformable walls, a model condition that closely replicates the actual mechanical setting in the Himalayas. This model allows for both channel walls to deform viscously, thus raising the question of how much overpressure they can retain inside the channel. We developed the deformable wall models with a channel geometry similar to that in rigid wall models, as shown in Fig. 8a. The footwall and the hanging wall of the channel were rheologically modelled as a viscous material, which provides a good approximation for simulation of long-term (millions of years) rheology of the lithosphere. Several earlier workers have used viscous rheology to model continental-scale deformation for an India-Tibet collision. The assumed viscosity values of the cold Indian craton range from 10 23 to 10 25 Pa s (e.g. Jiménez-Munt and Platt, 2006;Yang and Liu, 2013), whereas that of Himalayan subducted material ranges between 10 20 and 10 21 Pa s (e.g. Liu and Yang, 2003;Copley and Mckenzie, 2007). The viscosity  ratio (viscosity walls to viscosity channel) is therefore of the order of 10 2 to 10 5 . In our modelling we chose a conservative value of the viscosity ratio equal to 10 3 , where the walls and channel viscosities are 10 23 and 10 20 Pa s, respectively. We constrained the model boundaries with kinematic conditions as in the reference model with rigid walls. The lateral and the top boundaries of the footwall were subjected to a velocity of 4 cm yr −1 sub-parallel to the channel, whereas the lateral vertical boundaries of the hanging wall were fixed with zero horizontal velocity components, leaving the vertical component unconstrained. Its top boundary was also left unconstrained, allowing the material to extrude upward freely. The wallchannel interfaces had a no-slip condition.
Model results show channel flow patterns quite similar to those observed in rigid wall models. The extrusion occurs along a region close to the hanging wall in the form of a Poiseuille flow (Fig. 8a). It is noteworthy that the footwall undergoes little or no deformation, although being deformable. The entire footwall underthrusts by translational motion parallel to the channel. We calculated both the dynamic and the static pressures along the channel axis, and plotted them as a function of depth (Fig. 8b). Similarly to rigid wall models, the dynamic pressure here exceeds the static pressure by nearly 1.5 GPa. For example, the static pressure at a depth of 60 km is about 1.5 GPa, whereas the corresponding dynamic pressure reaches 3 GPa. The pressure plots clearly suggest that subduction channels with deformable walls can also give rise to large tectonic overpressures. The viscosity ratio (viscosity walls to viscosity channel) is therefore of the order of 10 2 to 10 5 . In our modelling we chose a conservative value of the viscosity ratio equal to 10 2 , where the walls and channel viscosities are 10 23 and 10 21 Pa s, respectively. For a viscosity ratio of 10 3 , the deformable wall models are found to be mechanically identical to rigid wall models. We also used a lower viscosity contrast of 10 2 , and found that even at this relatively low contrast there is significant overpressure in the subduction channel.

Condition at the bottom boundary
This is a critical boundary condition because it is directly related to the retention of overpressure. When we assign an outlet pressure (calculated lithostatic pressure at the depth of the bottom wall) to the bottom wall, TOP does not develop in the whole channel (Fig. 9).

Comparison with previous work
The occurrence of TOP has received much attention in the geological literature (e.g. Rutland, 1965;Mancktelow, 1993Mancktelow, , 1995Mancktelow, , 2008Petrini and Podladchikov, 2000;Schmalholz and Podladchikov, 2013;Schmalholz et al., 2014b). TOP has been argued to exist in both hard   (Mancktelow, 1993) and soft (Mancktelow, 1995) layers, and its occurrence has been predicted by force balance considerations independent of rheology (Schmalholz and Podladchikov, 2013;Schmalholz et al., 2014a, b). We have previously explored (Marques et al., 2018) the occurrence of TOP in higher viscosity layers intercalated in lower viscosity layers (layer-parallel shortening of a rheologically stratified lithosphere), and in a lower viscosity layer between higher viscosity walls (subduction zone).
Previous work has investigated the occurrence of TOP at all scales: (1) local variations in pressure (e.g. Mancktelow, 1993;Tenczer et al., 2001;Taborda et al., 2004;Marques et al., 2005aMarques et al., , b, c, 2014Podladchikov, 2003, 2004;Ji and Wang, 2011;Schmalholz and Podladchikov, 2014;Tajčmanová et al., 2014Tajčmanová et al., , 2015Angel et al., 2015), which in many cases is the natural consequence of the use of Stokes flow in the model, similarly to the numerical model used in the present study; and (2) TOP in subduction zones (e.g. Li et al., 2010;Reuber et al., 2016). Given the great dependence of pressure on geometry, boundary and ambient conditions, and flow pattern, we cannot compare our Stokes flow models directly with the cited self-consistent geodynamic models, because in these the controlling parameters are combined with many other variables and parameters that act simultaneously and change with time. Therefore, we analysed, separately, the effects of the various parameters and boundary conditions on pressure in order to gain a better understanding of the effects of each of them.
TOP has been investigated as a function of the tectonic environment (e.g. Stüwe and Sandiford, 1994;Petrini and Podladchikov, 2000;Vrijmoed et al., 2009;Pleuger and Podladchikov, 2014;Schmalholz et al., 2014a), and geometrical effects on TOP have also been addressed (e.g. Schmalholz and Podladchikov, 1999;Moulas et al., 2014), e.g. in downward-tapering (e.g. Mancktelow, 1995Mancktelow, , 2008 and references therein) and parallel-sided subduction channels, which have been argued to be the most appropriate configurations to model natural subduction zones. However, given the complexity and unsteady nature of subduction zones, the subduction channel can adopt all possible configurations, and the strictly parallel-sided configuration should be considered an exception rather than a rule, especially if we consider the 3-D, non-cylindrical, nature of subduction zones. Previous models have used two of the three main possible configurations of a subduction channel: parallel-sided and downward tapering, which have been shown to produce TOP < 3 (e.g. Li et al., 2010;Reuber et al., 2016). Here we investigated a different channel geometry, the upward-tapering channel. In fact, the parallel-sided geometry corresponds to W * m = 1, which can thus be considered an end-member of the UTC. Therefore, we can compare numerical results of overpressure obtained for parallel-sided and UTC channels, by looking at the graph where we vary W * m (Fig. 3a). Our best explanation for this effect is that the narrower the mouth the higher the flow confinement, which results in increased velocity gradient in the channel flow, and therefore the dynamic pressure.
The formation and exhumation of high-(HP) and ultrahigh-pressure (UHP) rocks is a persisting fundamental problem, especially regarding UHP rocks. The problem is even greater if one assumes that pressure estimated from paleopiezometry can be converted directly to depth, because then the UHP rocks must be exhumed from great depths. Several models have been proposed for the exhumation of HP and UHP rocks in several orogens (e.g. Hacker and Gerya, 2013;Warren, 2013;Burov et al., 2014a, b): channel flow (e.g. England and Holland, 1979;Mancktelow, 1995;Grujic et al., 1996;Beaumont et al., 2001Beaumont et al., , 2009Burov et al., 2001;Raimbourg et al., 2007;Gerya et al., 2008;Warren et al., 2008;Li and Gerya, 2009); eduction (e.g. Andersen et al., 1991;Kylander-Clark et al., 2012); buoyancy-driven crustal delamination and stacking (e.g. Chemenda et al., 1995Chemenda et al., , 1996Sizova et al., 2012); microplate rotation (e.g. Hacker et al., 2000;Webb et al., 2008); trans-mantle diapirism (e.g. Stöckhert and Little et al., 2011;Gordon et al., 2012); and slab rollback (e.g. Brun and Faccenna, 2008;Faccenda et al., 2009;Vogt and Gerya, 2014;Malusà et al., 2015). No model has so far provided a complete and unique explanation. The UTC model presented here is a new potential model to explain the exhumation of HP and UHP rocks, because it shows that it is possible to form rocks recording HP or UHP at depths < 60 km and to exhume them to the surface as a consequence of the flow configuration in the UTC.
Regarding the discrepancy between previous estimates of possible values of overpressure and ours, we call attention to two factors: (1) we use a subduction channel geometry, the UTC, not investigated previously; and (2) the values reported here are very large only for small W * m , or U > 5 cm yr −1 , or η > 10 21 Pa s. In other words, for relatively small tapering (W * m ), average plate tectonics velocities, and reasonable viscosities, the numerical results reported here for overpressure are not excessive, but nevertheless still very important as a factor for depth overestimation. The values used for the controlling parameters, W * m , α, and η are conservative; in fact, the model channel in Fig. 1c shows rather small tapering as compared with the cross section in Fig. 1b, but, nevertheless, the model overpressure is still quite high, especially at low depth.

Meaning and applicability of the numerical results
The numerical simulations reported here clearly discriminate the conditions favourable or unfavourable to the development of high TOP. The conditions that favour high TOP shallow in the subduction channel are the following: upwardtapering geometry, high viscosity (> 10 20 , which also means relatively low temperature), strong channel walls, general shear (i.e. simple + pure shears), low subduction angle, noslip condition at the bounding lateral walls, and a no-outlet condition at the bottom wall. All these conditions do not need to act simultaneously to generate TOP. We conclude that if during the unsteady evolution of a subduction zone, the boundary conditions, geometrical configuration, and ambient conditions meet the favourable model setting reported here, then high TOP can develop. Otherwise, only small TOP can be expected. In great contrast, the single action of low viscosity, downward-tapering geometry, weak channel walls, or outlet pressure at the bottom wall can prevent the development of TOP, or even promote underpressure.
We analysed the consistency between the numerical results and geological and geophysical data to constrain the most probable viscosity and pressure, at the same time satisfying a reasonable velocity at the channel mouth (i.e. exhumation rates; Fig. 7). On the one hand, the viscosity of rocks comprising the lithosphere can vary between 10 19 and 10 23 Pa s. On the other hand, overpressure is sensitive to the viscosity within the UTC, increasing rapidly with an increase in viscosity. Additionally, from the values shown in Figs. 3 and 4, the formation of HP rocks can occur at very shallow levels if η = 10 21 Pa s. However, despite the relatively wide range of possible viscosity values, η > 10 21 Pa s, combined with other favourable conditions in a Himalayan UTC yields overpressures > 8. This means that, for η = 10 22 Pa s, a rock metamorphosed at 50 km depth would record a total pressure equivalent to the lithostatic pressure at a depth of 400 km, which is not acceptable on the basis of our current knowledge of subduction zone dynamics. Therefore, we propose that the viscosity in the subduction channel is probably in the range 10 20 ≤ η ≤ 10 21 Pa s, in agreement with the estimates for Himalayan subducted material (between 10 20 and 10 21 Pa s) by Liu and Yang (2003) and Copley and Mckenzie (2007).
The UTC simulations show that there is no need for gravitational collapse, buoyancy-controlled crustal exhumation, or orogen-perpendicular pressure gradient induced by a topographic gradient to explain simultaneous reverse and normal fault kinematics in the MCT and STD, or inverse metamorphic grade or exhumation of HP rocks. We conclude that flow in a UTC, without the need for topography or density contrasts, can be responsible for these three simultaneous and seemingly paradoxical processes in the Himalayas.
An important question regarding TOP in nature still persists: why do we not see TOP in all subduction zones around the globe? On the one hand, our simulations indicate that roll-back subduction (transtension, in opposition to the favourable transpression) is unfavourable for the development of TOP. In contrast, collision-type subduction zones, like the Himalayas, with intervening old, cold, and strong lithospheres are favourable for TOP. On the other hand, the recognition of TOP depends on methods and analytical technology, as shown by the most recent literature on petrology. There is growing evidence that TOP is recorded by minerals, as shown by Tajčmanová et al. (2014Tajčmanová et al. ( , 2015, Moulas et al. (2013, and Angel et al. (2015). Constraints from host-inclusion elasticity show that TOP can greatly depart from lithostatic pressure; Angel et al. (2015) showed that deviations from lithostatic pressure in excess of 1 GPa can be readily produced in quartz inclusions within garnet in metamorphic rocks.

Comparison between model and nature
Inspired by the cross section of the natural upward-tapering channel shown in Fig. 1b, we investigated the effects of this geometry on TOP, and use it to find new explanations to the problems raised by the Himalayan geodynamics.
Given our incomplete knowledge of natural prototypes and the limitations of modelling very complex systems, we must distinguish between the theoretically and naturally possible values of overpressure. The study here reported for a UTC shows that relevant parameters like channel mouth width (W * m ), subduction dip (α), underthrusting velocity (U ), and viscosity (η) can produce very high overpressure; however, these theoretically possible values must be constrained by the current knowledge of the Himalayas, in particular exhumation velocities and spatial distribution, occurrence of HP and UHP rocks, and strength of the lithosphere bounding the subduction channel. Despite the natural constraints imposed by our knowledge of the current Himalayas, one cannot ignore that, under specific boundary conditions, geometrical configurations and parameter sets that could have existed in the past (e.g. much higher subduction velocity), high values of overpressure are theoretically possible, which should guide us in the search of new evidence in the natural prototype.
Previous models can explain channel flow, but neither account for the exhumation of HP rocks (Rubatto et al., 2013) nor the exhumation velocities  reported from the Himalayas. Our UTC model provides an alternative explanation for the pressure required for eclogite metamorphism (Hetényi et al., 2007;Zhang et al., 2014), and the process of rapid exhumation. For exhumation by extrusion to occur in the subduction channel, the flow pattern inside the channel must have a specific configuration, as in the UTC. In such a velocity configuration, underthrusting and exhumation on the channel footwall add to produce enhanced overthrusting on the MCT, and above the MCT along the line of flow reversal. Conversely, exhumation (upward flow) on the hanging wall is greater than underthrusting and produces relative normal fault displacement on the STD, not because the block to the north of the STD (hanging wall) moves down but because the rocks south of the STD (footwall) move up due to exhumation by extrusion.
Previous channel flow models can explain the exhumation mechanism; however, they leave a number of problems unaddressed. Here we raise some of these issues, pointing to our UTC model as a unifying model to explain the GHS evolution: 1. The classical channel flow model assumes that the entire GHS crustal mass thrusts up along the MCT, with concomitant normal motion on the STD (Poiseuille flow). However, recent studies have shown large-scale thrusts within the GHS Larson et al., 2015), suggesting a more complex kinematics of the extrusion process. The UTC model we propose here shows flow partitioning in the channel, leading to thrust-type shear localization within the model GHS.

2.
A typical channel flow model fails to explain the occurrence of HP rocks (> 1.5 GPa) close to the STD. Our UTC model yields an asymmetrical flow pattern in which HP or UHP materials extrude along a narrow zone located close to the STD.
3. The assumption of lithostatic pressure raises two main problems: (i) a conceptual problem, because the subduction channel is dynamic, therefore the lithostatic and dynamic pressures are not identical (e.g. Yamato and Brun, 2017); and (ii) a practical problem, because the exhumation velocities are calculated on the basis of depth estimated from ρgz (where z is depth), and not normalized by the overpressure. For instance, con-version of 2 GPa to depth using a static assumption (ρgz) yields a depth of ca. 70 km for a rock density of 2900 kg m −3 . However, the UTC flow develops an overpressure of the order of 2 at much smaller depths, and thereby yields lower exhumation rates, as compared to those calculated from petrologic modelling. Estimated metamorphic paths should reflect the shape of the isotherms in the subduction channel, which must have a relationship with velocity in order to carry cold rocks to depth, and preserve the HP and UHP mineral parageneses during exhumation.
4. Model velocities in the channel and at the channel mouth must be consistent with the values reported in the literature. Assuming lithostatic pressure, an exhumation rate of ∼ 15 mm yr −1 to a depth of at least 15 km was estimated by Ganguly et al. (2000). An estimate of 22-44 mm yr −1 , and increasing linearly with depth, was provided by Grujic et al. (2011). According to the UTC dynamic model, the assumption of lithostatic pressure where TOP = 2 yields an overestimation of the exhumation velocity by a factor of 2. If this is the case, then the velocity estimates have to be divided by 2 (15/2 = 7.5 mm yr −1 , and 33/2 = 16.5 mm yr −1 ). Our UTC model shows a high velocity layer with the materials flowing upward at a rate of 16 mm yr −1 at a depth of ca. 40 km, which is thus in agreement with the estimated average exhumation. The velocity map in Fig. 6 reveals variations in exhumation rates with depth, as predicted for the GHS in the Sikkim Himalayas by Ganguly et al. (2000), who showed that the exhumation was rapid (15 mm yr −1 ) to a depth of 15 km, and then decreased to ca. 2 mm yr −1 until a depth of 5 km. These values estimated for exhumation in the GHS constrain the theoretical values of overpressure numerically obtained by varying the amount of transpression. Transpression values > 10 % imply velocities at the mouth (exhumation) much higher than estimated for the GHS, therefore we conclude that transpression must be very limited (< 10 %).

A critical issue regarding overpressure in a subduction
channel is the strength of the channel walls to support high overpressure values. One of the most debatable boundary conditions in our modelling is the use of rigid walls. For this discussion, we can compare the subduction channel to a pressure vessel, in which the resistance of the vessel to internal pressure depends on two main parameters: the strength of the vessel (the lithosphere hosting the subduction zone), and the thickness of the pressure vessel walls (hoop stress). In nature, if the walls of the pressure vessel (subducting and overlying lithospheres) are old and cold, which is the case in the Himalayan collision, then their mechanical strength can be very high. If, additionally, the cold and strong lithosphere is thick, then the walls of the subduction chan-nel can support high overpressure, as indicated by the numerical results with viscous deformable walls. Given that the Indian plate and the TSS above the STD are almost undeformed (attesting to the rigidity contrast between footwall and hanging wall of the GHS) and thick, the channel walls were assumed undeformable in the reference simulations. In order to investigate the effects of viscous deformable walls on tectonic overpressure, we used viscosity contrasts (viscosity of channel walls to viscosity of subduction channel) down to 100, which are well within the accepted values of lithosphere viscosity (up to 10 23 Pa s) and subducted material (down to 10 19 Pa s). These simulations indicate that viscosity contrasts of 1000 or 100 do not change significantly the overpressure obtained with rigid walls. Another critical issue in overpressure build-up is the condition at the bottom boundary: if an outlet pressure is assigned to the bottom wall, then this boundary behaves as a leaking boundary that cannot retain dynamic pressure. However, the cold, thick, and strong lithospheres that comprise the Indian and Eurasian plates are a good argument against a leaking bottom boundary in a flat-ramp geometry such as the Himalayan collision zone. If, for some reason, the channel walls become weaker, in the brittle or viscous regimes, then the walls will yield and not be able to support large TOP.
6. In order to explain the non-linear variation in overpressure with channel dip (α) we need to analyse the variations in channel flow patterns with increasing α (Fig. 5). For low α values (15 • ), the underthrusting motion drags materials to a larger extent into the downward flow, and produces a large vortex in the deeper channel, where the curl dominates the flow field. Consequently, the dynamic pressure remains low. Note that flow divergence increases the dynamic pressure. With increasing α (20 • ) the flow pattern is characterized by the development of an extrusion channel on the hanging wall side, along which the material extrudes upward with flow convergence at the mouth. Such a negative divergence in the flow builds overpressure on the hanging wall side (Fig. 2d). With further increase in α the extrusion channel widens, and causes the overpressure to drop, as it happens in a pipe flow. This is the reason why the overpressure has a maximum at α around 20-25 • .
7. Inverted metamorphic grade has not been explained by previous models, but the UTC can provide an explanation if one considers the flow pattern shown in Fig. 2b. HP and UHP rocks can be exhumed by two flow cells, both inverting metamorphism because low-grade rocks go down close to the footwall, and high-grade rocks are exhumed close to the hanging wall.

Conclusions
The UTC model integrates and provides a robust physical explanation for a number of landmark features in the Greater Himalayan geodynamics, such as simultaneous reverse and normal faulting (channel flow), inversion of the metamorphic grade in the GHS, and exhumation of HP and UHP rocks along a narrow conduit close to the STD. Viscous flow in a UTC involves dynamic pressures in excess of lithostatic pressure, resulting in significant overpressure by a factor of more than 1.5, even at depths as shallow as 40 km. The UTC model predicts high-pressure (> 1.5 GPa) metamorphism of underthrusted rocks, e.g. eclogitization, to occur above 60 km depth. The UTC model shows that the GHS is segmented broadly into two sub-terrains with contrasting pressures: wide southern and narrow northern terranes, with pressures less and greater than 1.5 GPa, respectively. It further shows that temporal variations in channel dip may promote (α > 15 • ) or inhibit (α < 15 • ) exhumation.
Overpressure increases with increase in U , from TOP ≈ 1.5 for U = 2-5 cm yr −1 (current Indian velocity) to TOP ≈ 11 when U = 20 cm yr −1 (Indian velocity at 60-70 Ma), which means that in the past all the dynamic processes discussed here may have been enhanced. We tested different model setups (e.g. parallel walls) and boundary conditions (e.g. slip or no-slip condition at bounding walls) but they do not reproduce the prototype. The UTC model shows that tectonic pressure alone can drive the extrusion of HP rocks by channel flow. Viscous deformable walls do not affect overpressure significantly for viscosity contrasts (viscosity walls to viscosity channel) of the order of 1000 or 100. If, during the subduction process, the mouth width, dip, velocity, viscosity, or the conditions at the boundaries change in space and time, then TOP will change accordingly, and the exhumation mechanism (flow in the channel) and exhumation depth will also change.
TOP in a UTC is only possible if the condition at the bottom boundary is not outlet pressure; otherwise it behaves as a leaking boundary that cannot retain dynamic pressure. However, the cold, thick, and strong lithospheres that comprise the Indian and Eurasian plates are a good argument against a leaking bottom boundary in a flat-ramp geometry, which means that overpressure can build up to high values in the GHS. The argument does not apply if the channel is "open" at the bottom, because overpressure cannot be retained. This could be the case in subduction zones where there is no evidence for return flow and exhumation concomitant with subduction.
The numerical results reported here show that, under specific boundary conditions, geometrical configurations, and parameter sets, high values of overpressure are theoretically possible, which should guide us in the search of new evidence in the natural prototype to prove or disprove the natural existence of high overpressure.
Data availability. Data used in this article can be found in the references cited in Sects. 1 and 4. Figure A1. Evolution of dynamic and lithostatic pressures in a UTC with η = 10 21 Pa s and ρ = 2800 kg m −3 . The ratio dynamic pressure/lithostatic pressure corresponds to the overpressure factor (TOP).

A1 Boundary conditions and model set-up
The boundary conditions needed to complete the mathematical formulation for numerical simulations were as follows: (1) slab-parallel velocity applied on the underthrusting wall, consistent with the horizontal velocity of the Indian plate (5 cm yr −1 ; DeMets et al., 2010); (2) slip condition on (parallel to) the bottom boundary (Nábělek et al., 2009); (3) noslip condition on the hanging wall; (4) outlet condition with 1 atm pressure at the channel mouth; (5) gravity applied to the whole channel (∼ 9.8 m s −2 ); and (6) constant density of the material in the channel = 2800 kg m −3 (no phase changes in the models), representing the association felsic (mostly) and mafic granulites carrying the eclogite pods. Given that the Indian plate and the TSS above the STD are almost undeformed, attesting to the rigidity contrast between footwalls and hanging walls of the GHS, the channel walls were assumed undeformable in the simulations, except those testing the effects of viscous walls. In order to investigate flow kinematics and dynamic pressure in the channel, we varied the following parameters: (1) channel viscosity (η), (2) underthrusting velocity (U ), (3) channel dip (α), (4) channel mouth width (W m ), and (5) viscosity of channel walls. The viscosity in the channel was varied between 10 19 and 10 22 Pa s to cover a broad spectrum of crustal viscosities, as reported in the literature (Beaumont et al., 2001;England and Houseman, 1989;Copley et al., 2011). The current convergence rate between India and Eurasia has been estimated to be of the order of 5 cm yr −1 ; however, given the wide range of estimated velocities (Feldl and Bilham, 2006;DeMets et al., 2010), we ran numerical simulations varying U between 2 and 20 cm yr −1 (6.34 × 10 −10 to 6.34× −9 m s −1 in the model). Channel dip was varied between 15 and 30 • , which broadly covers the geometry of the GHS shown in different geological sections. We assumed W m = 25 to 100 km, and W b (width at the channel's base) = 150 or 200 km, from which we define W * m = W m /W b . We tested a viscosity con-trast (viscosity of channel walls to viscosity in the channel) of 1000 to investigate the effects of viscous deformable walls on overpressure. Despite varying all these parameters, the prime focus of our investigation concerned the simulations with U = 5 cm yr −1 , α = 20 • , W m = 100 km, and W b = 150 km, as they represent the most common and conservative values regarding published data. We then use the numerical results to constrain the viscosity, pressure, and velocity in the channel, consistent with current geological data and estimates.
The metamorphic processes occur in response to the total isotropic stress, called dynamic pressure, which is a sum of the tectonic (Stokes) and lithostatic pressures (ρgz, where ρ is density, g is gravitational acceleration, and z is depth; Figs. A1 and A2). The dynamic pressure results from the viscous flow driven by tectonic stresses in the gravity field. Using the present mechanical model, we evaluate the dynamic pressure to explain the occurrence of high-pressure rocks in the GHS, as a consequence of dynamic pressure in excess of lithostatic pressure at a given crustal depth. We define an overpressure factor (TOP) as the non-dimensional ratio between dynamic and lithostatic pressures. For a better understanding of overpressure in a UTC, we carried out a parametric study of TOP as a function of η, W m , α, U , and effective convergence velocity (horizontal velocity component > U ).

A2 Mathematical formulation
The mathematical model used in the present work is based on the Navier-Stokes equations for two-dimensional steadystate incompressible viscous flows: ρ ∂u ∂t + u + ∇u = −∇p + η∇ 2 u + F, where u is the velocity vector, p the pressure, ρ the density, η the dynamic viscosity, and F the external body force (gravity); ρ and η are constant. Then, defining the scaled variables x = x/L, u = u/U , p = p/P , and t = t/T , in terms of the characteristic length L, velocity U , pressure P , and time T = L/U , Eqs. (A1) and (A2) become where Re = ρU L/η and Eu = P /ρU 2 are, respectively, the Reynolds and Euler numbers. For flows at low characteristic velocity U and high viscosity η, inertial terms Eu and Re in Eq.
(3) become negligible. We thus obtain the Stokes approximation of the momentum equation for quasi-static (creeping) flows, which in dimensional form and under a gravity field reads as −∇p + η∇ 2 u + F = 0. (A5) Figure A2. Overpressure in the UTC under the velocity field shown in Fig. 3.
The Stokes equations were solved on the 2-D domain illustrated in Fig. 1c, which was filled with an incompressible viscous linear fluid. The flow equations, with the boundary conditions specified, were solved in the primitive variables u ≡ (u, v) and p over a finite element mesh, using the algorithm for incompressible Stokes flows implemented in COM-SOL.
Author contributions. FOM had the idea and designed the project, ran part of the numerical simulations, analysed and discussed the results, and wrote the text; NM and GR contributed to the design of the project, to the analysis and discussion of the results, and to the writing; SG ran part of the numerical simulations and contributed to the discussion of the results and to the writing; SB contributed to the discussion of the results and to the writing.
Competing interests. The authors declare that they have no conflict of interest.