S Duprat-Oualid
email: sylvia.duprat-oualid@univ-rennes1.fr

P Yamato

P Pitra

Major role of shear heating in intracontinental inverted metamorphism: Inference from a thermo-kinematic parametric study

Keywords: Inverted metamorphism, Shear heating, Numerical modelling, Thermal properties of the lithosphere

Inverted metamorphism corresponds to the stacking of high-temperature metamorphic units structurally on top of lower-temperature units and is commonly observed along main thrusts in major orogens. Yet, in spite of many existing models, the origin and preservation of the metamorphic inversion in intracontinental collision belts is still debated. In this study, we use a crustal-scale 2D thermo-kinematic model in order to investigate the key parameters controlling the inversion of the geothermal gradient at crustal scale. Our results confirm that the kinematic framework strongly impacts the thermal evolution around the thrust. Erosion velocity and thermal conductivity of rocks are two parameters that control the spatial location of the thermal perturbation and the intensity of inversion, respectively. However, even in extreme kinematic configurations, i.e., convergence velocities > 3cm.yr -1 and relatively high thrust dip angles ~ 30˚, the thermal inversion is fleeting and thrust temperatures cannot reach the high temperature peak values characteristic of natural occurences (> 600˚C) if shear heating is not taken into account. Conversion of mechanical energy into heat represents a main contribution to the thermal budget along main crustal shear zones. It leads to high temperature conditions in the thrust zone and our results attest that it is the only process that allows the preservation through time of an intense thermal inversion. Our quantification shows that shear heating is much more efficient than other processes such as accretion and surface denudation and is compatible with the observations of inverted metamorphism in the Himalayan or Variscan belts, for example. This comparison with natural occurrences suggests that the formation and preservation of intracontinental inverted metamorphism require shear zone viscosity values of the order of 10 20 -10 21 Pa.s for convergence velocities between 1 and 3 cm.yr -1 .

Introduction and state of the art

An inverted metamorphic sequence is characterised by an upward intensification of metamorphism, typically an increase of the metamorphic peak temperature. Inverted metamorphism is mostly related to major thrusts, where convergence velocities are high. The major occurrences can be ascribed to one of t hree distinct geological contexts: (i) metamorphic soles beneath ophiolites (e.g., [START_REF] Williams | Metamorphic aureoles beneath ophiolite suites and alpine peridotites; tectonic implications with west Newfoundland examples[END_REF][START_REF] Jamieson | P-T paths from high temperature shear zones beneath ophiolites[END_REF][START_REF] Hacker | The role of deformation in the formation of metamorphic gradients: Ridge subduction beneath the Oman Ophiolite[END_REF][START_REF] El-Naby | Evolution of the Pan-African Wadi Haimur metamorphic sole, Eastern Desert, Egypt[END_REF]; (ii) oceanic subduction zones (e.g., [START_REF] Toksöz | Temperature Field and Geophysical Effects of a Downgoing Slab[END_REF][START_REF] Platt | Metamorphic and deformational processes in the Franciscan Complex, California: Some insights from the Catalina Schist terrane[END_REF]; (iii) major thrusts in continental collision (e.g., [START_REF] Lefort | Himalayas : The collided range -Present knowledge of the continental arc[END_REF][START_REF] Kohn | P-T-t data from central Nepal support critical taper and repudiate largescale channel flow of the Greater Himalayan sequence[END_REF].

The present study is focused on the third geological setting. The best documented example is the inverted metamorphic sequence associated with the Main Central Thrust (MCT) zone in the southern Himalayas (e.g., [START_REF] Frank | Relations between metamorphism and orogeny in a typical section of the Indian Himalayas[END_REF][START_REF] Lefort | Himalayas : The collided range -Present knowledge of the continental arc[END_REF][START_REF] Sinha-Roy | Himalayan main central thrust and its implications for Himalayan inverted metamorphism[END_REF][START_REF] Arita | Origin of the inverted metamorphism of the lower Himalayas, Central Nepal[END_REF][START_REF] Hubbard | Thermobarometric constraints on the thermal history of the Main Central Thrust Zone and Tibetan Slab, eastern Nepal Himalaya[END_REF][START_REF] Jain | Inverted metamorphism in the intracontinental ductile shear zone during Himalayan collision tectonics[END_REF][START_REF] Henry | Kinematic, thermal and petrological model of the Himalayas: constraints related to metamorphism within the underthrust indian crust and topographic elevation[END_REF][START_REF] Harrison | A late Miocene-Pliocene origin for the central Himalayan inverted metamorphism[END_REF][START_REF] Harrison | A model for the origin of Himalayan anatexis and inverted metamorphism[END_REF][START_REF] Harrison | The origin of Himalayan anatexis and inverted metamorphism : Models and constraints[END_REF][START_REF] Caddick | Burial and exhumation history of a Lesser Himalayan schist: Recording the formation of an inverted metamorphic sequence in NW India[END_REF][START_REF] Kohn | P-T-t data from central Nepal support critical taper and repudiate largescale channel flow of the Greater Himalayan sequence[END_REF], but other examples were described in the Variscan belt [START_REF] Burg | Inverted metamorphic zonation and large-scale thrusting in the Variscan Belt: an example in the French Massif Central[END_REF][START_REF] Arenas | High-pressure micro-inclusions and development of an inverted metamorphic gradient in the Santiago Schists (Ordenes Complex, NW Iberian Massif, Spain): evidence of subduction and syncollisional decompression[END_REF][START_REF] Štípská | Inverted metamorphic zonation in a basement-derived nappe sequence ; eastern margin of the Bohemian Massif[END_REF][START_REF] Ballèvre | Palaeozoic history of the Armorican Massif : Models for the tectonic evolution of the suture zones[END_REF][START_REF] Pitra | Inverted metamorphic field gradient towards a Variscan suture zone (Champtoceaux Complex, Armorican Massif, France)[END_REF], in the Caledonian belt [START_REF] Andreasson | Occurrence and significance of inverted metamorphic gradients in the western Scandinavian Caledonides[END_REF][START_REF] Mason | Inverted isograds at Sulitjelma, Norway: the result of shear-zone deformation[END_REF][START_REF] Watkins | Geothermometry and geobarometry of inverted metamorphic zones in the W central Scottish Dalradian[END_REF][START_REF] Johnson | A discussion of possible heat sources during nappe stacking: the origin of Barrovian metamorphism within the Caledonian thrust sheets of NW Scotland[END_REF] or in the Appalachian belt [START_REF] Camiré | Development of inverted metamorphic gradient in the internal domain of the Taconian belt, Gaspé Peninsula[END_REF]. The pressure and temperature conditions within such thrust zones attest to a medium pressure, medium-high temperature metamorphism. As an example, inverted metamorphic sequences in both the Himalayas and the Variscan belt have recorded peak temperatures between 500 and 700˚C and peak pressures between 8 and 11 kbar (e.g., [START_REF] Burg | Inverted metamorphic zonation and large-scale thrusting in the Variscan Belt: an example in the French Massif Central[END_REF][START_REF] Macfarlane | An evaluation of the inverted metamorphic gradient at Langtang National Park, Central Nepal Himalaya[END_REF][START_REF] Guillot | An overview of the metamorphic evolution of central Nepal. In : "Geology of Nepal[END_REF][START_REF] Kohn | P-T-t data from central Nepal support critical taper and repudiate largescale channel flow of the Greater Himalayan sequence[END_REF][START_REF] Pitra | Inverted metamorphic field gradient towards a Variscan suture zone (Champtoceaux Complex, Armorican Massif, France)[END_REF][START_REF] Corrie | Metamorphic history of the central Himalaya, Annapurna region, Nepal, and implications for tectonic models[END_REF]. Natural intracontinental metamorphic inversions are thus characterised by inverted thermal gradients between 10 and 50˚C.km -1 (e.g., [START_REF] Kohn | P-T-t data from central Nepal support critical taper and repudiate largescale channel flow of the Greater Himalayan sequence[END_REF][START_REF] Pitra | Inverted metamorphic field gradient towards a Variscan suture zone (Champtoceaux Complex, Armorican Massif, France)[END_REF].

Such intense thermal perturbations over a thickness of several kilometres necessarily imply geodynamic processes at crustal or even lithospheric scale. Several numerical models have been developed, most of them constrained by or aiming to fit data from the Himalayan belt (e.g., [START_REF] Jamieson | Tectonic assembly of inverted metamorphic sequences[END_REF][START_REF] Jamieson | Crustal channel flows: 2. Numerical models with implications for metamorphism in the Himalayan-Tibetan orogen[END_REF][START_REF] Henry | Kinematic, thermal and petrological model of the Himalayas: constraints related to metamorphism within the underthrust indian crust and topographic elevation[END_REF][START_REF] Bollinger | Mountain building in the Nepal Himalaya: Thermal and Kinematic model[END_REF][START_REF] Kohn | P-T-t data from central Nepal support critical taper and repudiate largescale channel flow of the Greater Himalayan sequence[END_REF].

They suggest two alternative hypotheses to explain the metamorphic inversion.

On the one hand, the inverted zonation is the result of a post-metamorphic deformation of a preexisting "normal" metamorphic sequence. This may result from late thrusts cutting through the initial metamorphic sequence (e.g., [START_REF] Brunel | Etude pétro-structurale des chevauchements ductiles himalayens sur la transversale de l'Everest-Makalu (Népal oriental)[END_REF], passive deformation of metamorphic isograds within a ductile shear zone (e.g., [START_REF] Jain | Inverted metamorphism in the intracontinental ductile shear zone during Himalayan collision tectonics[END_REF][START_REF] Grujic | Ductile extrusion of the Higher Himalayan Crystalline in Bhutan : evidence from quartz microfabrics[END_REF][START_REF] Gibson | Deformation-induced inverted metamorphic field gradients: an example from the southeastern Canadian Cordillera[END_REF] or crustal-scale folding (e.g., [START_REF] Searle | Thermal model for the Zankstar Himalaya[END_REF][START_REF] Stephenson | Inverted metamorphism and the Main Central Thrust : field relations and thermobarometric constraints from the Kishtwar Window, NW Indian Himalaya[END_REF].

On the other hand, the inverted zonation may result from a temporary spatial perturbation of crustal isotherms [START_REF] Lefort | Himalayas : The collided range -Present knowledge of the continental arc[END_REF]Peacock, 1987a;[START_REF] Jamieson | Tectonic assembly of inverted metamorphic sequences[END_REF][START_REF] Jamieson | Crustal channel flows: 2. Numerical models with implications for metamorphism in the Himalayan-Tibetan orogen[END_REF][START_REF] Burg | Viscous heating allows thrusting to overcome crustalscale buckling: numerical investigation with application to the Himalayan syntaxes[END_REF]. In this hypothesis, metamorphism is synchronous with the deformation event. Some authors propose that inverted metamorphism occurs due to a diffusive heat transfer across the major thrust from the hotter upper unit to the colder subjacent unit (e.g., [START_REF] Lefort | Himalayas : The collided range -Present knowledge of the continental arc[END_REF][START_REF] Shi | Two-dimensional modeling of the P-T-t paths of regional metamorphism in simple overthrust terranes[END_REF][START_REF] England | The interpretation of inverted metamorphic isograds using simple physical calculations[END_REF].

Particularly, the thermal inversion may be due to a "channel flow", i.e., the extrusion of deep hot crustal rocks above a colder plate (e.g., [START_REF] Grujic | Ductile extrusion of the Higher Himalayan Crystalline in Bhutan : evidence from quartz microfabrics[END_REF][START_REF] Jamieson | Tectonic assembly of inverted metamorphic sequences[END_REF][START_REF] Jamieson | Crustal channel flows: 2. Numerical models with implications for metamorphism in the Himalayan-Tibetan orogen[END_REF]. This model presents the advantage of reproducing the ranges of metamorphic pressures and temperatures characterising the outcroping sections of the MCT [START_REF] Beaumont | Himalayan tectonics explained by extrusion of a low-viscosity crustal channel coupled to focused surface denudation[END_REF].

However, it requires strong erosion localised on the thrust front and specific internal physical properties in terms of viscosity and angle of friction. Furthermore, it fails to reproduce some important metamorphic and geochronological records on both sides of the MCT [START_REF] Kohn | P-T-t data from central Nepal support critical taper and repudiate largescale channel flow of the Greater Himalayan sequence[END_REF].

Beyond the first-order role of the kinematic framework and thermal diffusion, several studies have specifically focused on the contribution of radiogenic heat, including the accretion of radioactive material across the active thrust, and on the role of erosion [START_REF] Royden | The steady state thermal structure of eroding orogenic belts and accretionary prisms[END_REF][START_REF] Ruppel | Pressure-temperature-time paths from two-dimensional thermal models: Prograde, retrograde, and inverted metamorphism[END_REF][START_REF] Huerta | The interdependence of deformational and thermal processes in mountain belts[END_REF][START_REF] Huerta | The thermal structure of collisional orogens as a response to accretion, erosion, and radiogenic heating[END_REF][START_REF] Huerta | The effects of accretion, erosion and radiogenic heat on the metamorphic evolution of collisional orogens[END_REF][START_REF] Bollinger | Mountain building in the Nepal Himalaya: Thermal and Kinematic model[END_REF].

Accretion models allowing the accumulation of highly radioactive material over a duration longer than 30 Myr can lead to significant increase of temperature within the hanging wall (e.g., [START_REF] Huerta | The interdependence of deformational and thermal processes in mountain belts[END_REF][START_REF] Huerta | The thermal structure of collisional orogens as a response to accretion, erosion, and radiogenic heating[END_REF]. However, such studies involve disputable initial assumptions.

As an example, the orogenic accretionary wedges of these models (e.g., [START_REF] Huerta | The thermal structure of collisional orogens as a response to accretion, erosion, and radiogenic heating[END_REF] are deeply rooted at 30 to 60 km depths. Reasonable accretion velocities associated with lower, realistic erosion velocities lead to wedges several hundreds of kilometers wide at the surface. Outcrouping rocks across such wide zones should consequently be characterised by a continuous intense deformation and none of the known frontal thrusts displays a such configuration. Moreover, in most cases, strong accretion alone cannot reproduce the high temperatures observed in natural inverted metamorphic sequences. On the contrary, accretion leads to progressive cooling of the active thrust plane to steady temperatures lower than 400˚C at 30 km depth with no significant thermal inversion [START_REF] Royden | The steady state thermal structure of eroding orogenic belts and accretionary prisms[END_REF]. Concomitant surface denudation characterised by an erosion velocity greater than the accretion velocity is necessary to raise the temperatures in the thrust zones (e.g., [START_REF] Royden | The steady state thermal structure of eroding orogenic belts and accretionary prisms[END_REF][START_REF] Bollinger | Mountain building in the Nepal Himalaya: Thermal and Kinematic model[END_REF].

In the absence of accretion processes, moderate erosion velocities can lead to thrust temperatures of the order of 600˚C and to the development of inverted isotherms comparable to those observed in natural inverted metamorphic sequences [START_REF] Royden | The steady state thermal structure of eroding orogenic belts and accretionary prisms[END_REF]. However, these conclusions require a thrust activity lasting for more than 70 Myr and exhumation of ultrahigh pressure metamorphic rocks from more than 100 km depth [START_REF] Royden | The steady state thermal structure of eroding orogenic belts and accretionary prisms[END_REF]. Such implications are not compatible with natural cases where rocks recorded metamorphic peak pressures of 8 to 11 kbar. From a mechanical approach, erosion catalyses the development of main orogenic structures (e.g., [START_REF] Burg | Viscous heating allows thrusting to overcome crustalscale buckling: numerical investigation with application to the Himalayan syntaxes[END_REF], and some models highlight the important influence of both the erosion and the exhumation on the thermal perturbation (e.g., [START_REF] Jamieson | Tectonic assembly of inverted metamorphic sequences[END_REF][START_REF] Beaumont | Crustal channel flows: 1.Numerical models with applications to the tectonics of the Himalayan-Tibetan orogen[END_REF][START_REF] Bollinger | Mountain building in the Nepal Himalaya: Thermal and Kinematic model[END_REF][START_REF] Kohn | P-T-t data from central Nepal support critical taper and repudiate largescale channel flow of the Greater Himalayan sequence[END_REF].

The role of heat production by conversion of mechanical energy into thermal energy in a major shear zone (shear heating) was also considered to complete the thermal budget and the possible development of thermal inversion around a thrust [START_REF] Minear | Thermal regime of a downgoing slab[END_REF][START_REF] Toksöz | Temperature Field and Geophysical Effects of a Downgoing Slab[END_REF][START_REF] Lefort | Himalayas : The collided range -Present knowledge of the continental arc[END_REF][START_REF] Graham | Thermal regimes and regional metamorphism in the vicinity of overthrust faults: an example of shear heating and inverted metamorphic zonation from southern California[END_REF][START_REF] Scholz | Shear heating and the state of stress on faults[END_REF][START_REF] England | Pressure-temperature-time paths of regional metamorphism. I. Heat transfer during the evolution of regions of thickened continental crust[END_REF][START_REF] Pavlis | The Role of Strain Heating in the Evolution of Megathrusts[END_REF][START_REF] Molnar | Temperatures, Heat Flux, and Frictional Stress Near Major Thrust Faults[END_REF][START_REF] England | The interpretation of inverted metamorphic isograds using simple physical calculations[END_REF][START_REF] Burg | The role of viscous heating in Barrovian metamorphism of colli-sional orogens: thermomechanical models and application to the Lepontine Dome in the Central Alps[END_REF][START_REF] Burg | Viscous heating allows thrusting to overcome crustalscale buckling: numerical investigation with application to the Himalayan syntaxes[END_REF], but without real quantification. [START_REF] Kidder | Shear hesating not a cause of inverted metamorphism[END_REF] specifically refute the importance of shear heating under high convergence velocity (~10 cm.yr -1 ) and attribute the thermal inversion to accretion process.

The demonstration provided by [START_REF] Kidder | Shear hesating not a cause of inverted metamorphism[END_REF] that shear heating is not a cause of inverted metamorphism in accretion domains is based on an oceanic subduction beneath an active arc [START_REF] Kidder | High temperatures and inverted metamorphism in the schist of Sierra de Salinas, California[END_REF]. The very high convergence velocity (~10 cm.yr -1 ) and the high initial geothermal gradients characterising the overriding forearc [START_REF] Kidder | High temperatures and inverted metamorphism in the schist of Sierra de Salinas, California[END_REF] and the subducting oceanic plate easily and rapidly lead to intense thermal inversion.

With low shear heating (~10 -5 W.m -3 implying very low viscosities of about 10 17 Pa.s), heat advection and conduction dominate the thermal budget: both the overriding plate and the slab are rapidly cooled and the subducting material is progressively and slowly heated.

Consequently, the rocks accumulated in the accretion zone are representative of the thermal state of the thrust plane. Understandingly, the thermal peaks are thus colder from the top to the base of the accretion area and such thermal inversion i s compatible with the thermochronologic data from the Pelona Schist of southern California [START_REF] Kidder | Shear hesating not a cause of inverted metamorphism[END_REF].

Considering a stronger shear heating along the active thrust, which is superposed on the overall cooling pattern, decreases the thermal inversion intensity recorded in the accreted sediments. Nevertheless, such shear heating assumptions are not generally applicable because very low shear stress and viscosities are involved to balance the high strain rate (10 -12 to 10 -11 s -1 ).

Previous studies of syn-deformational intracontinental inverted metamorphism used various models with very different degrees of complexity (e.g., [START_REF] Shi | Two-dimensional modeling of the P-T-t paths of regional metamorphism in simple overthrust terranes[END_REF][START_REF] England | The interpretation of inverted metamorphic isograds using simple physical calculations[END_REF][START_REF] Jamieson | Crustal channel flows: 2. Numerical models with implications for metamorphism in the Himalayan-Tibetan orogen[END_REF]. However, none of the proposed models is commonly accepted to give a general explanation to the development of syn-deformation inverted metamorphism. Whereas analytical studies provide good mathematical solutions for understanding the respective influence of kinematics and heat diffusion on the thermal evolution, the other factors involved (e.g., erosion, various heat sources, variable kinematic configuration, thermo-dependence laws) are difficult to address analytically but can be solved and quantified numerically.

In the present paper, our goal is not to fit a particular natural case (e.g., the Himalayas). Rather, using a systematic approach, we aim to explore and provide a detailed synthesis of the different conditions required for the formation and preservation of crustalscale inverted metamorphism along one intracontinental thrust. We define and compare the relative importance of each one of the parameters and processes involved: kinematic setting, thrust geometry, erosion velocity, rock properties (in particular the thermal diffusivity), accretion velocity, radiogenic heat production and shear heating. The goal is to infer (i) the conditions required to locally invert the thermal field; (ii) the exact impact of each one of these parameters on the intensity, the location and the duration of the inversion and (iii) the conditions required to preserve the thermal inversion through time. Finally, we discuss how and why shear heating is crucial for the formation and preservation of inverted thermal gradients under realistic kinematic settings.

Numerical model

Initial setup

In order to test independently the influence of the various parameters and processes on the thermal evolution of a crustal thrust, a simple numerical model is defined in which each component can be controlled. The 2D-model setup (Fig. 1a) is focused on the first 80 km of a continental lithosphere, which includes a 30 km thick crust. In order to study the thermal evolution of a major continental thrust zone, a thrust cutting through the whole continental crust with a dip angle θ is simulated by imposing a velocity contrast between the downgoing plate and the upper plate, moving with a velocity V lp and V up , respectively (Fig. 1). The convergence velocity V th is consequently given by V th = V lp + V up . The thrust is mainly represented by a straight segment, the dip angle of which, however, progressively decreases towards both the base of the crust and the surface in order to avoid any advection errors due to local strong velocity direction contrasts. Since the present study focuses on the slanted part of the thrust and because the temperature evolution around the thrust and in the deepest levels of the fault is not noticeably affected by the configuration of the lower block velocity field beneath the first 30 km (see Appendix A,Fig. A.1), the downgoing plate intentionally lies flat under the upper crust to simplify the boundary conditions and the isostasy computation.

The initial geotherm is defined in the whole model according to the following equation:

T(x,z,t 0 ) = T surface + z⋅ q base k + T source (x,z) (1)
where T surface , q base , k and T source are the surface temperature, the mantle heat flux, the thermal conductivity and the possible local heat sources, respectively. Initially, the only implemented heat source is the radiogenic heat production. It is considered as a permanent property of the material and is defined according to the following decreasing exponential law from [START_REF] Turcotte | Geodynamics[END_REF]:

Q(z,t 0 ) = Q 0 + exp -z z r       (2)
where Q 0 is the surface radiogenic heat production, and z r corresponds to the length scale for the decrease in Q with depth [START_REF] Turcotte | Geodynamics[END_REF] controlling the curvature of the vertical profile of Q. Eq. ( 2) allows to easily test a wide range of possible distributions of radiogenic production and both parameters Q 0 and z r , were independently tested. In the presence of radiogenic heat sources, Eq. ( 1) is then only z-dependent and becomes:

T(z,t 0 ) = T surface + z⋅ q base k + Q 0 ⋅ z r 2 k ⋅ 1 -exp -z z r             (3)

Computation methods

Both kinematic and thermal computations are based on a regular Eulerian grid with a space resolution of 1 km. During the thrust activity, at each time step, advection and thermal diffusion are realised separately in order to minimise the numerical diffusion of temperature inside a non-static environment by using the marker-in-cells method [START_REF] Gerya | Introduction to numerical geodynamic modelling[END_REF]. The advection of rock properties and temperature is ensured by a field of initially randomly distributed particles (see Appendix A for more details). Three kinematic parameters can be distinguished here (Fig. 1b). The first one is the burial velocity Vb of the downgoing plate corresponding to the vertical component of the velocity field V lp and defined by:

Vb = V lp • sin(θ) (4)
The second parameter defines the vertical accretion velocity a. As in previous studies (e.g., [START_REF] Royden | The steady state thermal structure of eroding orogenic belts and accretionary prisms[END_REF][START_REF] Huerta | The interdependence of deformational and thermal processes in mountain belts[END_REF][START_REF] Huerta | The thermal structure of collisional orogens as a response to accretion, erosion, and radiogenic heating[END_REF][START_REF] Huerta | The effects of accretion, erosion and radiogenic heat on the metamorphic evolution of collisional orogens[END_REF], accretion is realised at each time step by adding a horizontal velocity component a' (Fig. 1b) in the whole model so that both the active thrust position and the convergence velocity remain fixed. Both parameters are related according to:

a' = a • tan(θ) (5)
The third kinematic parameter corresponds to the maximum erosion velocity Ve linked to the upper plate velocity V up by the following equation:

Ve = V up • sin(θ) (6)
The erosion velocity Ve controls the upper plate velocity V up in the models (Fig. 1b, see

Appendix A for more details).

Compared to the tectonic advection times, the isostatic balance is reached instantaneously [START_REF] Crough | Isostatic rebound and power-law flow in the asthenosphere[END_REF]. Thus, at each time step, the isostatic response is calculated in the entire model thereby producing a realistic topography associated with the thrust.

Then, for the thermal part, the following heat equation is solved by using the implicit finite difference method on the Eulerian grid. In order to limit numerical diffusion, the Eulerian computation is coupled with intermediate thermal interpolations with the Lagrangian particles (see Appendix A for more details):

ρ⋅ Cp⋅ ∂T ∂t - ∂ ∂x k ∂T ∂x       - ∂ ∂z k ∂T ∂z       = Q + Hs (7)
The heat equation (Eq. ( 7)) is composed by a heat diffusion part where ρ, Cp and k are respectively the density, the heat capacity and the thermal conductivity (left side of Eq. ( 7))

and by a heat production part (right side of Eq. ( 7)) where Q and Hs represent the radiogenic heat production and the shear heating, respectively.

The parameters ρ, Cp, and k define the capacity of the material to diffuse heat. This property can be reduced to one parameter, the thermal diffusivity κ, defined as κ = k/(ρ.Cp).

Each one of these three parameters depends on the lithology, but also on the physical environment of the rock (e.g., pressure, temperature). Various empirical laws describe the temperature and pressure dependence of these parameters in rocks (e.g., [START_REF] Zoth | Handbook of Terrestrial Heat Flows Density Determination[END_REF][START_REF] Seipold | The variation of thermal transport properties in the Earth's crust[END_REF][START_REF] Seipold | Temperature dependance of thermal transport properties of crystalline rocks -a general law[END_REF][START_REF] Abdulagatov | Effect of Pressure and Temperature on the Thermal Conductivity of Rocks[END_REF][START_REF] Hofmeister | High-Pressure Geoscience Special Feature: Pressure dependence of thermal transport properties[END_REF][START_REF] Whittington | Temperature-dependent thermal diffusivity of the Earth's crust and implications for magmatism[END_REF]. In the present study, we tested independently the influence of a range of realistic values for each parameter on the thermal evolution around the thrust (see Table 1). First, these values were considered constant in the crust and in the mantle. In a second time, their temperature and pressure dependence is considered and their influence on the thermal field investigated. These results are discussed in section 3.3.2.

Shear heating Hs results from the transformation of mechanical energy into heat during deformation and is recognised as a major component of the total energy budget (e.g., [START_REF] Brun | Strain heating and thermal softening in continental shear zones: a review[END_REF][START_REF] Lachenbruch | Heat Flow and Energetics of the San Andreas Fault Zone[END_REF][START_REF] Chester | Fracture surface energy of the Punchbowl fault, San Andreas system[END_REF]. Shear heating is a function of the strain rate and stress tensors ( ij and σ ij respectively) and is computed according to the following equation:

Hs = σ xx ⋅ ˙ ε xx + σ zz ⋅ ˙ ε zz + 2⋅ σ xz ⋅ ˙ ε xz (8)
The strain rate tensor is defined as ij = 0.5 [(∂V i /∂j) + (∂V j /∂i)] where ∂V i /∂j corresponds to the variation of the velocity of the component i in the direction j and the same convention applies for ∂V j /∂i. Because shear zones related to inverted metamorphism can be several kilometres thick (e.g., Grugic et al., 1996), ductile rheology is generally attributed to the whole thrust (e.g., [START_REF] Searle | Thermal model for the Zankstar Himalaya[END_REF][START_REF] Jamieson | Tectonic assembly of inverted metamorphic sequences[END_REF][START_REF] Kidder | High temperatures and inverted metamorphism in the schist of Sierra de Salinas, California[END_REF].

In the present study, the shear zone thickness is space resolution dependent and is thus approximately of about 1 km. Shear heating is computed in a ductile manner, where stress and strain rate are linked by the effective viscosity η of the material constituting the thrust:

σ ij = 2⋅ η⋅ ˙ ε ij (9)
Viscosity depends on temperature, strain rate and creep parameters (e.g., [START_REF] Goetze | Stress and temperature in the bending lithosphere as constrained by experimental rock mechanics[END_REF][START_REF] Carter | Flow properties of continental lithosphere[END_REF][START_REF] Ranalli | Rheological stratification of the lithosphere[END_REF]. It tends to decrease with increasing strain rate and temperature [START_REF] Meissner | Weakness of the lower continental crust: a condition for delamination, uplift, and escape[END_REF]. Inversely, both the rock deformation and the associated shear heating depend on the viscosity (e.g., [START_REF] Burg | Viscous heating allows thrusting to overcome crustalscale buckling: numerical investigation with application to the Himalayan syntaxes[END_REF][START_REF] Hartz | Toasting the jelly sandwich : The effects of shear heating on lithospheric geotherms and strength[END_REF]. However, in order to better define the first order role of viscosity, we simplified the problem by using a constant viscosity η all within the thrust (see discussion in section 3.4).

Model results

Burial velocity effects

In order to test the validity of our numerical model, we first realised a comparative study with analytical solutions from [START_REF] England | The interpretation of inverted metamorphic isograds using simple physical calculations[END_REF]. The results are described in Appendix B.

In extreme kinematic configurations (e.g., V th = 5 cm.yr -1 , θ = 30˚), which favour the development of thermal inversion, the thermal evolution around the thrust and in the whole model can be best characterised. In this case, because the initial geotherm is linear, the thermal inversion starts on a large part of the thrust with a relatively homogeneous intensity.

Then, the affected area progressively deepens and widens and the maximum intensity is localised towards the deepest levels of the thrust. On each point of the thrust, the vertical temperature gradient ∂T/∂z first decreases, reaches a minimum value (negative when thermal inversion occurs) and then increases again until a steady state when the temperature advection and the heat diffusion are completely balanced. Nevertheless, even under such extreme kinematic settings, none of the (V th , θ) combinations allows the preservation of the thermal

Erosion effects

In our model, erosion can be activated and is only controlled by the erosion velocity Ve (Eq. ( 6), Fig. 1 and Table 1). In this case, contrary to the previous analysis (section 3.1, Appendix B), the thrust velocity V th is given by a combination of the lower and the upper plates velocities (V lp and V up , see Fig. 1b). Consequently, in order to test the influence of the erosion velocity independently of the kinematic setting, either the downgoing plate velocity

V lp or the thrust velocity V th have to be fixed. Here, we chose to keep constant the thrust velocity V th for two reasons. (i) Kinematic simulations having the same vertical component Vz of the thrust velocity (Fig. 1b) are characterised by very similar thermal evolutions (see Appendix B). Each given Vz value corresponds to a group of (V th , θ) combinations. This allows us to easily make connections and distinctions between the thermal evolution around the thrust and, on the one hand, the erosion velocity, and on the other hand, the convergence velocity and the thrust dip angle. In comparison with the MCT in the Himalayas which is characterised by a high downgoing velocity of the Indian plate of ~ 2 cm.yr -1 for the last 10 to 20 Myr (e.g., [START_REF] Bilham | GPS measurements of present-day convergence across the Nepal Himalaya[END_REF][START_REF] Lavé | Active folding of fluvial terraces across the Siwaliks Hills, Himalayas of central Nepal[END_REF] and a dip angle of ~ 15 ̊ (e.g., [START_REF] Kaneko | Timing of Himalayan ultrahigh-pressure metamorphism: sinking rate and subduction angle of the Indian continental crust beneath Asia[END_REF], Vz values of the order of 0.5 cm.yr -1 are realistic and greater values are excessive. (ii) A further parametric analysis concerning the shear heating (section 3.4) takes into account the velocity contrast on the thrust as a major parameter. Controlling the thrust velocity instead of the downgoing plate velocity is thus more practical to compare the combined effects of the kinematic setting, the erosion velocity and the shear heating. The distribution of the velocity field on both sides of the thrust can be described as the ratio between the erosion velocity Ve and the vertical component Vz of the thrust velocity (see Fig. 1).

Erosion and motion of the upper plate, for a constant V th , does not lead to significant variations of the reached intensity (Fig. 2a) and of the thickness of the affected area (Fig. 2b) of the thermal inversion during the thrust activity. However, it strongly influences its space and time locations. Increasing erosion and hence displacement of the upper plate tends to localise the thermal inversion at shallower depths along the thrust (Fig. 2). Taking out crustal material by erosion leads to progressive cooling of the remaining shallowest material controlled by surface thermal conditions. The maximum thermal inversion intensity is not restricted to the base of the slanted part of the thrust but moves to the surface with increasing erosion velocities. However, because erosion velocities are not sufficiently high to efficiently balance the thermal diffusion, even for extreme cases, thermal inversion cannot occur at very superficial levels along the thrust.

The time of maximum inversion is not significantly changed, even if a slight delay can be noticed. The duration of the thermal inversion is widely extended with erosion. However, if for extreme kinematic and erosion values thermal inversion can be preserved until the steady state, the most realistic kinematic setting tested (corresponding to Vz = 0.5 cm.yr -1 ) never leads to a thermal inversion, whatever the applied erosion velocity (Fig. 2).

Hence, when the displacement along the thrust is entirely ensured by the downgoing plate (Ve = 0), temperature on each point along the thrust is continually decreasing until the steady state. When the upper plate velocity V up approximates the lower plate velocity V lp , (i.e., Ve ≈ 0.5Vz), the thermal evolution along the thrust tends to be reversed. Thus, for Ve > 0.5Vz, temperature at a given point of the thrust starts to increase in the beginning of the thrust activity (Fig. 2a, case Ve = 2Vz/3) due to the fast exhumation of the deep hot portions of the upper plate. This trend increases with the vertical component Vz of the thrust velocity. After reaching a local temperature peak, maximum thermal inversion is quickly reached and the temperature along the whole thrust decreases until the steady state. Actually, the thrust profile flattening in depth limits the temperature increase because the deepest rocks involved in the thrust come from the base of the upper continental block and cannot thus be hotter than the initial Moho temperature.

To summarise, erosion tends to localise the thermal inversion towards shallower depths of the thrust and significantly increases the thermal inversion duration. Nevertheless, it does not affect the intensity and the thickness of the thermal inversion. Furthermore, even in extreme cases, the initial temperature increase on the thrust is limited (maximum 25˚C) and is far from sufficient to account for the high temperatures known in major thrust zones associated with inverted metamorphism.

Rock properties effects

In the lithosphere, conduction is generally the predominant heat transport mechanism ahead of convection and radiation [START_REF] Lachenbruch | Heat Flow and Energetics of the San Andreas Fault Zone[END_REF][START_REF] Abdulagatov | Effect of Pressure and Temperature on the Thermal Conductivity of Rocks[END_REF]). As shown previously in section 2.2, ρ, Cp and k constitute the three rocks properties involved in the heat diffusion equation. Moreover, as defined by Eq. ( 3), the initial geotherm in the lithosphere is not only controlled by k and the radiogenic heat sources, but also by the basal mantle heat flux q base . These thermal parameters are here considered independently and the impact of the radioactive character of the crust on the thermal evolution around the thrust is also discussed.

Heat diffusion properties and the mantle heat flux

The ranges of realistic values (Table 1) for each thermal parameter depend strongly on the lithology. On the one hand, the k values are generally comprised between 2 W.m -1 .K -1 and 4 W.m -1 .K -1 [START_REF] Clauser | Thermal conductivity of rocks and minerals[END_REF][START_REF] Turcotte | Geodynamics[END_REF]. This high degree of variability of more than 100 % proportionally controls the thermal diffusivity κ with the same order of magnitude. On the other hand, the heat capacity Cp varies up to ~ 50 %, from 800 to 1200 J.kg -1 .K -1 [START_REF] Waples | A review and evaluation of specific heat capacities of rocks, minerals, and subsurface fluids. Part 1: Minerals and nonporous rocks[END_REF], involving maximum differences for κ of ~ 30 %.

Depending on the lithology, the crustal rock density ρ varies from ~ 2400 kg.m -3 for sedimentary rocks to ~3300 kg.m -3 for basic eclogites as an extreme case, but is generally comprised between 2500 and 3000 kg.m -3 with an associated variability of ~ 20%. Heat diffusion is therefore mainly controlled by the variation of thermal conductivity, and the influence of the variations of density and heat capacity are less significant (see Appendix C,

Fig. C.1). This has two reasons:

First, the hotter is the geotherm (obtained for a low thermal conductivity in the crust or a high mantle heat flux), the higher is the thermal inversion (Fig, 3a,c). The basal heat flux, q base , only involved in the initial and boundary conditions, influences the intensity of the thermal inversion proportionally to the magnitude of its variability. For example, between q base = 20 mW.m -2 and q base = 40 mW.m -2 , the initial thermal gradient is two times higher and the thermal inversion intensity reached on each point of the thrust is doubled (Fig. 3c).

Proportionally to a q base increase, the intensification of the inversion with decreasing k is more significant. From k = 3 to 2 W.m -1 .K -1 , the thermal inversion is approximately two times higher (Fig. 3a). Thus, the initial thermal structure of the lithosphere, defined by the distribution of both q base and k, constitutes a key parameter controlling the intensity reached by the thermal inversion. Nevertheless, the location of the inversion, both in space and time, is mostly insensitive to either q base or k (Fig. 3).

Second, the strong correlation between κ and k leads the latter to strongly influence the heat diffusion through the thrust. A low thermal conductivity decreases the heat diffusion and thus favours thermal inversion during thrusting. For this reason, the influence of both k and q base on the thermal evolution around the thrust is different (Fig. 3). A decrease of k induces a more significant intensification of the inversion than an increase of q base (Fig. 3a,c).

Furthermore, contrary to the intensification of the basal heat flux, a decrease of k induces a widening of the thrust part affected by the inversion, a thickening of the area where the thermal inversion occurs (Fig. 3b,d) and a continuation in time of the thermal inversion (Fig.

3).

Temperature and pressure dependences

Thermal properties also vary as a function of temperature and pressure. Various experiments have been realised and some empirical laws have been proposed in order to define the thermal dependence of the thermal conductivity k [START_REF] Zoth | Handbook of Terrestrial Heat Flows Density Determination[END_REF][START_REF] Horai | The effect of pressure on the thermal conductivity of silicate rocks up to 12 kbar[END_REF][START_REF] Clauser | Thermal conductivity of rocks and minerals[END_REF][START_REF] Seipold | Temperature dependance of thermal transport properties of crystalline rocks -a general law[END_REF][START_REF] Voosten | Influence of temperature on thermal conductivity, thermal capacity and thermal diffusivity for different types of rock[END_REF][START_REF] Abdulagatov | Effect of Pressure and Temperature on the Thermal Conductivity of Rocks[END_REF][START_REF] Whittington | Temperature-dependent thermal diffusivity of the Earth's crust and implications for magmatism[END_REF] and heat capacity Cp [START_REF] Waples | A review and evaluation of specific heat capacities of rocks, minerals, and subsurface fluids. Part 1: Minerals and nonporous rocks[END_REF][START_REF] Whittington | Temperature-dependent thermal diffusivity of the Earth's crust and implications for magmatism[END_REF]. Generally, with increasing temperature, k decreases, whereas Cp increases. Conversely, thermal conductivity tends to increase with increasing pressure [START_REF] Horai | The effect of pressure on the thermal conductivity of silicate rocks up to 12 kbar[END_REF][START_REF] Seipold | Depth dependance of thermal transport properties for typical crustal rocks[END_REF][START_REF] Kukkonen | Temperature and pressure dependencies of thermal transport properties of rocks: implications for uncertainties in thermal lithosphere models and new laboratory measurements of high-grade rocks in the central Fennoscandian shield[END_REF][START_REF] Osako | Simultaneous measurements of thermal conductivity and thermal diffusivity for garnet and olivine under high pressure[END_REF][START_REF] Abdulagatov | Effect of Pressure and Temperature on the Thermal Conductivity of Rocks[END_REF][START_REF] Hofmeister | High-Pressure Geoscience Special Feature: Pressure dependence of thermal transport properties[END_REF]. However, the pressure effect is limited compared to the influence of temperature [START_REF] Abdulagatov | Effect of Pressure and Temperature on the Thermal Conductivity of Rocks[END_REF].

The reported variation of each thermal property as a function of temperature and pressure lies within the range of constant values tested in our numerical simulations (Table 1).

Only the variability of k has a significant impact on the thermal evolution around the thrust.

However, k depends on both the lithology and the physical environment, and this impact is difficult to quantify. Indeed, in spite of the numerous laboratory experiments (e.g., [START_REF] Čermák | Thermal conductivity and specific heat of minerals and rocks, in Landolt-Börnstein: Numerical Data and Functional Relationships in Science and Technology, New[END_REF][START_REF] Osako | Simultaneous measurements of thermal conductivity and thermal diffusivity for garnet and olivine under high pressure[END_REF][START_REF] Abdulagatov | Effect of Pressure and Temperature on the Thermal Conductivity of Rocks[END_REF], none of the proposed empiric laws [START_REF] Zoth | Handbook of Terrestrial Heat Flows Density Determination[END_REF][START_REF] Clauser | Thermal conductivity of rocks and minerals[END_REF][START_REF] Seipold | Temperature dependance of thermal transport properties of crystalline rocks -a general law[END_REF][START_REF] Voosten | Influence of temperature on thermal conductivity, thermal capacity and thermal diffusivity for different types of rock[END_REF][START_REF] Abdulagatov | Effect of Pressure and Temperature on the Thermal Conductivity of Rocks[END_REF][START_REF] Whittington | Temperature-dependent thermal diffusivity of the Earth's crust and implications for magmatism[END_REF] appears to be undisputable (Fig. 4). Diverse equations relative to the same lithology appear inconsistent with one another (Fig. 4c,d). Furthermore, they are associated with large uncertainties (e.g., Seipold, 1998, Fig. 4a) and rarely integrate both temperature and pressure dependences in the same time [START_REF] Kukkonen | Temperature and pressure dependencies of thermal transport properties of rocks: implications for uncertainties in thermal lithosphere models and new laboratory measurements of high-grade rocks in the central Fennoscandian shield[END_REF][START_REF] Abdulagatov | Effect of Pressure and Temperature on the Thermal Conductivity of Rocks[END_REF].

Regarding only the vertical distribution of the theoretical lithologies (after [START_REF] Turcotte | Geodynamics[END_REF], the thermal conductivity increases from the superficial levels (sediments) to the mantle rocks. On the other hand, within each lithology, k decreases significantly with temperature (and hence depth; Fig. 4). This progressive drop may reach 30-50 % from 0 to 600 ˚C. Concerning the pressure influence, k increases with depth in the shallowest crustal levels due to the reduction of the pore volume [START_REF] Abdulagatov | Effect of Pressure and Temperature on the Thermal Conductivity of Rocks[END_REF].

Beyond ~ 100 MPa, k tends to reach a linear behaviour [START_REF] Horai | The effect of pressure on the thermal conductivity of silicate rocks up to 12 kbar[END_REF][START_REF] Clauser | Thermal conductivity of rocks and minerals[END_REF][START_REF] Abdulagatov | Effect of Pressure and Temperature on the Thermal Conductivity of Rocks[END_REF], and is supposed to evolve following the thermodependence law appropriate for the considered lithology.

To sum up, within the range of realistic values, both Cp and ρ have negligible effects on the thermal evolution of the lithosphere, and the diffusive transport of heat is controlled by the variability of k. The temperature and pressure dependences do not change the order of variability, but they can lead to different thermal behaviours inside the lithosphere [START_REF] Maierová | The effect of variable thermal diffusivity on kinematic models of subduction[END_REF]. In particular, the dependence of the k value on temperature for a given lithology and pressures > 100 MPa favours the thermal inversion in depths. These physical and lithologic dependences are generally not taken into account in the numerical models although the values allocated to the thermal conductivity have to be considered with caution [START_REF] Nabelek | Strain heating as a mechanism for partial melting and ultrahigh temperature metamorphism in convergent orogens: Implications of temperature-dependent thermal diffusivity and rheology[END_REF].

Radiogenic heat production

As for the heat diffusion properties (section 3.3.1, Appendix C), numerical simulations were realised in an extreme kinematic context (V th = V lb = 5 cm.yr -1 and θ = 30˚) in order to quantify the effect of the radiogenic heat production. They show that although the radiogenic heat production contributes to the initial geothermal profile, it does not influence the thermal inversion around the thrust in terms of intensity and localisation (Appendix C, Fig. C.2).

The previous analysis revealed that the burial velocity, the erosion velocity and the thermal conductivity of rocks are the main parameters controlling the thermal evolution around a crustal scale thrust. Nevertheless, considering realistic kinematic settings (V th < 3 cm.yr -1 ), none of these parameters allows the preservation of an established thermal inversion in time.

Shear heating effects

Shear heating constitutes a part of the conservation of energy and a study of the thermal evolution around a thrust during its activity requires to take into account the relative displacement between the two crustal blocks as a source of heat [START_REF] Toksöz | Temperature Field and Geophysical Effects of a Downgoing Slab[END_REF][START_REF] Brun | Strain heating and thermal softening in continental shear zones: a review[END_REF][START_REF] Pavlis | The Role of Strain Heating in the Evolution of Megathrusts[END_REF]. By considering shear heating in the viscous way, two parameters have to be examined: the strain rate and the effective viscosity η within the thrust (section 2.2, Eq. ( 8) and ( 9)). Since the velocity field is not mechanically computed in our numerical models but imposed, both and η can be considered independently. The initial kinematic conditions used here are based on realistic cases for continental collision: the maximum dip angle θ of the thrust set to 15˚ and the velocity thrust values between 1 and 3 cm.yr -1 , implying strain rates from 3.10 -13 to 10 -12 s -1 .

In rocks, viscosity is temperature-and strain rate-dependent. Considering an undeformed crust, the viscosity field is characterised by a decrease from ~ 10 25 Pa.s in the cold upper crust to ~ 10 21 Pa.s in the hot lower crust (e.g., [START_REF] Burg | The role of viscous heating in Barrovian metamorphism of colli-sional orogens: thermomechanical models and application to the Lepontine Dome in the Central Alps[END_REF]. The interdependence of strain rate and viscosity tends to a steady state and, as an example, with a convergence velocity of 2 cm.yr -1 , the viscosity drops to values of 10 19 to 10 21 Pa.s within major shear zones [START_REF] Burg | The role of viscous heating in Barrovian metamorphism of colli-sional orogens: thermomechanical models and application to the Lepontine Dome in the Central Alps[END_REF]. Consequently, we chose to apply this range of viscosity values, constant in time and space along the thrust. In our kinematic models, the strain field is only localised along the thrust (Fig. 1) thus limiting the frictional heat production to the shear zone. In order to highlight the effects of shear heating, no additional heat source is taken into account, and the erosion is not considered so that the downgoing plate only ensures the displacement along the thrust.

The theoretical amount of heat produced (Eq. ( 8)) along the thrust is very sensitive to even small variations of viscosity (Fig. 5c). Considering the lowest viscosity (η = 10 19 Pa.s)

and strain rate ( = 3.17 . 10 -13 s -1 ), the resulting frictionnal heat produced (Hs ≈ 4 µW.m -3 ) is of the same order of magnitude as the maximum realistic radiogenic heat production (Table 1). A viscosity increase induces a proportional increase of shear heating along the thrust. For example, an increase of 10 0.25 Pa.s will typically result in an elevation of ~ 10 0.25 = 178 % of the produced heat (Fig. 5c). This way, (i) shear heating can efficiently balance the diffusive cooling of the overriding plate. For high viscosity values, the thermal gradient above the thrust can be elevated beyond the initial crustal geotherm. Consequently, if the heating produced by shearing becomes higher than the cooling driven by the underthrusting of cold material, the surface heat flow at the top of the upper block increases during thrusting [START_REF] Hartz | Toasting the jelly sandwich : The effects of shear heating on lithospheric geotherms and strength[END_REF]. (ii) The thermal budget at the top of the downgoing crust has two sources: the transfer by diffusion of the heat inherited in the upper block and the shear heating along the thrust. Consequently, in presence of shear heating, the maximum thermal perturbation appears and develops under the thrust, whereas it straddles the thrust when shear heating is not considered (Fig. 5a,b).

As shown in Fig. 6, the amount of heat produced by shearing leads to an intense thermal inversion characterised by a significant thickening of the area affected and to a long duration. Unlike the other situations considered, shear heating allows the preservation of the thermal inversion in time until the steady state under realistic kinematic settings (V th = 1-3 cm.yr -1 ; Fig. 5b and6).

Thermal inversion intensities, even at shallow depths, can exceed 100˚C.km -1 and temperatures extend up to more than 1000˚C (Fig. 7). Both ranges of temperatures between 500 and 700˚C and of thermal inversion intensities from 10 to 50˚C.km -1 are reached between 20 and 30 km depth for different combinations of viscosity and convergence velocity (Fig. 8).

Consequently, shear heating alone can reproduce and preserve in time thermal features comparable to natural occurences of inverted metamorphic sequences. Furthermore, our model suggests that, in these cases, viscosity values within the thrust zone have to range between 10 20 to 10 21 Pa.s (Fig. 8).

In short, for small variations of viscosity or strain rate across a thrust (constant width of the sheared area), the thermal inversion is strongly influenced by shear heating in terms of intensity, space and time location. Although shear heating does not necessary imply an inversion, it plays a major role in the preservation of the perturbation in time.

Shear heating, accretion and radiogenic heat production competition

Accretion processes were implemented while keeping shear heating in action along the active thrust. As previously, the thrust dip angle θ was set to 15˚. The convergence velocity V th was fixed at 2 cm.yr -1 and no erosion process was considered. The vertical accretion velocities a tested were 1 and 2 mm.yr -1 and are comparable to the values used in previous works [START_REF] Royden | The steady state thermal structure of eroding orogenic belts and accretionary prisms[END_REF][START_REF] Huerta | The interdependence of deformational and thermal processes in mountain belts[END_REF][START_REF] Huerta | The thermal structure of collisional orogens as a response to accretion, erosion, and radiogenic heating[END_REF][START_REF] Huerta | The effects of accretion, erosion and radiogenic heat on the metamorphic evolution of collisional orogens[END_REF][START_REF] Bollinger | Mountain building in the Nepal Himalaya: Thermal and Kinematic model[END_REF][START_REF] Kidder | Shear hesating not a cause of inverted metamorphism[END_REF]. In order to evaluate the thermal influence of accumulation of radioactive material in the accretion prism, an elevated accretion velocity (a = 2 mm.yr -1 ) was considered and two radiogenic profiles (Eq. ( 2) and ( 3)) characterised by a common surface radiogenic heat production Q 0 = 2 µW.m -3 but differing in the specific depth z r were tested : (i) a realistic configuration with z r = 15 km and (ii) an extreme case with z r = 30 km such that the whole continental crust is highly enriched in radioactive elements. Mantle flux was adapted to each thermal configuration in order to ensure a realistic initial temperature at the base of the crust between 420 and 500˚C.

Results are synthesised in Fig. 9 in terms of maximum thermal inversion intensity and maximum temperatures reached in the thrust zone after the thermal steady state is attained.

Whether we take into account accretion and radiogenic heat production or not, no significant thermal inversion occurs if the viscosity η is lower than 10 20.25 Pa.s. Low viscosity implying low amounts of shear heating (Eq. ( 8) and ( 9)), it follows that, in the framework of realisitc kinematic and thermal settings, it is impossible to reproduce metamorphic inversions comparable to natural observations by considering accretion processes only. Conversely, if η exceeds 10 21 Pa.s, temperatures and geothermal gradients within the thrust zone are excessively elevated whatever the depth considered. Thus, under a given kinematic context, the range of viscosities leading to realistic thermal inversions is relatively narrow, of the order of 10 20 to 10 21 Pa.s (Fig. 8 and9).

The concomitance of accretion and shear heating implies that material previously deformed and warmed within the thrust zone is accumulated in the hanging wall while subjacent volumes are in turn sheared and heated. For low amounts of shear heating, upward heat transfert by accretion leads to a smaller pertubation on the geothermal gradient around the active thrust (Fig. 9a) and no change in the maximum temperatures (Fig. 9b). Conversely, when shear heating is significant, accretion of hot material efficiently reduces the heat diffusion from the thrust zone to the overriding plate. Thermal inversion in the active thrust zone is consequently more intense and this trend is more pronounced in depth. However, because the material heated in the shear zone is continuously moved upwards out of the active shear zone, the maximum temperatures are significantly lower (e.g., about 50˚C lower with a 2 mm.yr -1 accretion velocity and η = 10 20.5 Pa.s. and almost 100˚C lower with a 1 mm.yr -1 accretion velocity and η = 10 20.75 Pa.s.) (Fig. 9b).

Accretion of radioactive material associated with shear heating has no significant effect on the maximum intensity of thermal inversion, even if extremely radioactive crustal material is considered (Fig. 9a). Because the top of the downgoing plate is highly radioactive, the maximum temperatures reached along the thrust are obviously hotter than without radiogenic heating.

In short, whatever the accretion velocity or the radiogenic heat production, even in extreme cases, none of these processes can lead to thermal inversion without shear heating in the framework of realistic kinematic and thermal settings. Furthermore, maximum thermal features in the thrust zone are more influenced by small viscosity variations than by high changes of the accretion velocity and of the amount of radiogenic heating.

Shear heating and erosion

On the one hand, surface denudation enhances the localisation of the thermal inversion at shallower depths and can lead to increased temperatures along the thrust (see section 3.2).

On the other hand, shear heating strongly influences the crustal thermal budget and highly intensifies the thermal inversion (see sections 3.4 and 3.5). Consequently, numerical simulations taking into account both processes are necessary to complete our parametric study (Fig. 9). Here again, tests were performed with a thrust dipping 15˚ and a convergence velocity V th of 2 cm.yr -1 . Accretion and radiogenic heat production were switched off.

Viscosities ranged between 10 19 and 10 21 Pa.s and erosion velocity was set to Ve = 2 mm.yr -1 .

This results in 38% of the thrust velocity V th being ensured by the upper plate velocity V up .

When the hanging wall velocity partially ensures the convergence (i.e., erosion is taken into account), maximum temperatures attained i n and along the thrust zone are significantly rised (Fig. 9b). As an example, at 20 km depth, a 2 mm.yr -1 erosion velocity induces an increase of the maximum temperature of 80˚ for η = 10 20.25 Pa.s, of 110˚ for η = 10 20.5 Pa.s and of 140˚ for η = 10 20.75 Pa.s. This thermal jump increases with the viscosity.

More generally, the maximum temperatures recorded along the thrust for a such erosion velocity (Ve = 2 mm.yr -1 ) show an increase of 10 to 18% at 30 km depth, of 20 to 25% at 20 km depth and 19 to 32% at 10 km depth.Thermal influence of erosion is thus significantly intensified when shear heating is important, whereas it is almost negligible at low frictionnal heat production. However, if exhumation of the overthrusting rocks favours the thermal raise, its influence is inconsiderable with regard to the thermal perturbation intensity (Fig. 9a).

Erosion velocities greater than 2 mm.yr -1 yield similar results and follow the same trend.

To sum up, concomitance of erosion processes with shear heating along a major thrust in realistic configurations (i.e., Ve < Vb) does not lead to significant variations in the maximum thermal inversion intensity. But the maximum temperatures reached along the thrust are highly increased by several tens of degrees.

Discussion and conclusions

Through the present parametric study, we tested and compared the relative influences of the different processes and parameters involved in the thermal evolution in space and time of the continental crust around a major thrust zone:

(i) the kinematic configuration of the thrust, i.e., the thrust dip angle, the convergence velocity, and the distribution of the thrust velocity on both sides of the shear zone (expressed here as the erosion velocity);

(ii) the thermal rock properties controlling the heat diffusion across the rocks, i.e., the thermal conductivity, the heat capacity and the density;

(iii) the heat sources, i.e., the basal mantle flux, the internal radiogenic heat production, including the accretion of radioactive material, and the mechanical energy converted into heat along the thrust during its activity (shear heating).

The initial kinematic configuration plays a major role. W h ereas the thermal perturbation around the thrust increases in depth with increasing burial velocity of the downgoing plate, erosion localises the perturbation towards shallower depths of the crust and significantly increases its duration. However, for identical convergence velocities, any changes in the erosion velocity do not impact the intensity and the thickness of the thermal inversion contrary to the thermal conductivity. Nevertheless, whereas a limited and transient thermal inversion can be achieved in extreme kinematic configurations, inversion is never observed in the range of reasonable (reported from nature) values of these parameters. Among all the processes and parameters tested, only shear heating allows the preservation of an intense thermal inversion through time under realistic kinematic schemes.

The thermal conductivity problem

The thermal conductivity, in the range of realistic values (i.e., between 2 and 4 W.m - 1 .K -1 ), has a major influence on the heat transfer by diffusion and on the initial geothermal gradient. Consequently, it strongly influences the local thermal gradient during the thrust activity, and when thermal inversion occurs, the intensity of this perturbation (Appendix C).

As a rule, the initial thermal structure of the lithosphere is an important parameter to take into account in the models. However, tests comparing the thermal evolution of crusts with an identical initial temperature field but different thermal conductivities (the slope of the initial geotherm is thus balanced by the amount of mantle heat flux) attest to the significant effect of k on the intensity and duration of the thermal inversion (Appendix C).

Whereas the thermal conductivity (or thermal diffusivity) is generally considered to be constant in numerical models, it strongly depends on both the lithology and the physical environment, more particularly the temperature. The choice of the k values [START_REF] Nabelek | Strain heating as a mechanism for partial melting and ultrahigh temperature metamorphism in convergent orogens: Implications of temperature-dependent thermal diffusivity and rheology[END_REF][START_REF] Maierová | The effect of variable thermal diffusivity on kinematic models of subduction[END_REF] and corresponding empirical laws chosen in the numerical thermal models are capital and can lead to significant differences of the thermal behaviour of the lithosphere. For pressures higher than 100 MPa, the pressure effect is negligible and the thermal influence dominates with a strong decrease of k of ~ 1 W.m -1 .K -1 for temperatures ranging from the surface to the base of the crust. Empirical thermo-dependence laws for various rock types were developed for atmospheric pressures. Empirical laws describing the evolution of k as a function of temperature under pressures higher than 100 MPa should give a better estimation of the thermal conductivity field in the crust for a given lithology.

Furthermore, metamorphic reactions occurring along and across the thrust, the varying amount of fluids in the rocks and the possible associated chemical changes linked to their circulation and/or chemical diffusion have also an impact on the evolution of the thermal conductivity of a given rock. These parameters are difficult to integrate in the models at present, but may have significant effects.

Shear heating, the key parameter

Shear heating, as defined in this paper and in the majority of other numerical models, considers that the whole energy required to overcome the fault friction is converted into thermal energy. However, this is a simplification and other energy transfers should be taken into account and their influence on the formation and preservation of the inverted thermal gradient tested. (i) The energy dissipated during a continuous deformation is not entirely converted into heat but can also be used to change the physical properties of the rock inside the sheared area (e.g., the rock microstructure; [START_REF] Brun | Strain heating and thermal softening in continental shear zones: a review[END_REF]. (ii) Metamorphic reactions occurring within and astride the shear zone can significantly consume or release heat [START_REF] Oxburgh | The physico-chemical behaviour of the descending lithosphere[END_REF][START_REF] Anderson | Thermal model of subduction with dehydration in the downgoing slab[END_REF][START_REF] Devore | Relations between subduction, slab heating, slab dehydration and continental growth[END_REF]Peacock, 1987aPeacock, , b, 1990)). (iii) Potential strain and gravitational energy [START_REF] Kanamori | Energy Partitioning During an Earthquake[END_REF]) can be transformed into pressurisation heat without displacement along the fault. Although these aspects should be taken into account in order to refine the general picture, they are not expected to modify the first-order message delivered by the modelling presented here.

Indeed, our results show that without shear heating, the occurence of thermal inversion requires extreme kinematic configurations (e.g., [START_REF] Kidder | Shear hesating not a cause of inverted metamorphism[END_REF], i s systematically transient (a few million years) and the inversion paroxysm is reached astride the thrust. In such a case, in order to preserve the trace of the thermal inversion in metamorphic parageneses (i.e., to avoid that any other metamorphic event significantly overprints the inverted record), exhumation of the whole sequence needs to be initiated during the short period of time concerned by the inversion (some Myr) and to be extremely fast with a very high associated erosion velocity. On the other hand, under realistic kinematic configurations for intracontinental thrusts (up to 3 cm.yr -1 ), shear heating is the only parameter which allows the preservation of the thermal inversion in time until the thermal steady state. It also significantly influences the spatial location and extent of the thermal inversion: (i) frictionnal heating along the thrust causes the localisation of the inversion area under rather than astride the thrust, within the uppermost levels of the lower plate, and (ii) the more frictionnal heat is produced, the thicker is the area affected by the thermal inversion leading to inversion zones up to tens of kilometers thick. The amount of heat produced depends on the local strain rate and the shear zone viscosity (Eq. ( 8)). The values of both parameters, and particularly the effective viscosity, strongly influence the crustal scale thermal perturbation, especially in terms of thermal inversion intensity.

The present study attests that shear heating is the key-parameter and allows to reproduce thermal features coming from natural inverted metamorphic sequences (Fig. 8).

Both the high metamorphic peak temperatures (500-700˚C at 8-11 kbar; e.g., [START_REF] Guillot | An overview of the metamorphic evolution of central Nepal. In : "Geology of Nepal[END_REF][START_REF] Pitra | Inverted metamorphic field gradient towards a Variscan suture zone (Champtoceaux Complex, Armorican Massif, France)[END_REF][START_REF] Corrie | Metamorphic history of the central Himalaya, Annapurna region, Nepal, and implications for tectonic models[END_REF] and the thermal inversion intensities (generally comprised between 10 and 50˚C.km -1 ; e.g., [START_REF] Kohn | P-T-t data from central Nepal support critical taper and repudiate largescale channel flow of the Greater Himalayan sequence[END_REF][START_REF] Pitra | Inverted metamorphic field gradient towards a Variscan suture zone (Champtoceaux Complex, Armorican Massif, France)[END_REF]) can be reproduced.

In parallel, our results suggest that the development and preservation of such inverted metamorphic zoning in realistic kinematic and thermal setting requires that the viscosity of the material involved has to be between 10 20 and 10 21 Pa.s.

Applicability to the Himalayan inverted metamorphism

The MCT in the Himalayas is characterised by a dip angle of 15˚ (e.g., Kaneko et al., and a downgoing velocity of the Indian plate of 2 cm.yr -1 . (e.g., [START_REF] Bilham | GPS measurements of present-day convergence across the Nepal Himalaya[END_REF][START_REF] Lavé | Active folding of fluvial terraces across the Siwaliks Hills, Himalayas of central Nepal[END_REF]. High metamorphic temperature peaks within the MCT reach 600 to 650˚C and pressure peaks are of 8 to 11 kbar (e.g., [START_REF] Macfarlane | An evaluation of the inverted metamorphic gradient at Langtang National Park, Central Nepal Himalaya[END_REF][START_REF] Guillot | An overview of the metamorphic evolution of central Nepal. In : "Geology of Nepal[END_REF][START_REF] Corrie | Metamorphic history of the central Himalaya, Annapurna region, Nepal, and implications for tectonic models[END_REF], corresponding to depths of about 30-40 km (assuming that these pressures correspond to lithostatic pressures). In agreement with the kinematical constraints, Fig. 8 shows that shear zone viscosity of 10 20.5 to 10 20.6 Pa.s is required to reproduce such maximum temperatures. The corresponding steady state thermal inversion intensity is 20 to 30˚C.km -1 , hence matching the estimates extracted from [START_REF] Kohn | P-T-t data from central Nepal support critical taper and repudiate largescale channel flow of the Greater Himalayan sequence[END_REF].

By taking into account erosion and accretion of radioactive rocks, as proposed by previous studies (e.g., [START_REF] Huerta | The thermal structure of collisional orogens as a response to accretion, erosion, and radiogenic heating[END_REF][START_REF] Bollinger | Mountain building in the Nepal Himalaya: Thermal and Kinematic model[END_REF] the range of consistent viscosities is slightly extended to 10 20.35 to 10 20.65 Pa.s and the corresponding maximum thermal inversion intensity spreads between 20 and 40˚C.km -1 (Fig. 9). Whether or not accretion and surface denudation are considered, the thermal features of the inverted metamorphic zoning across the MCT can simply be reproduced by using shear heating and involving viscosities in the shear zone close to those commonly expected (e.g., [START_REF] Burg | The role of viscous heating in Barrovian metamorphism of colli-sional orogens: thermomechanical models and application to the Lepontine Dome in the Central Alps[END_REF].

Hence, our model, without pretending to explain the inverted metamorphism related to the MCT, shows that the downgoing of the Indian plate under the Himalayan belt with kinematical constraints from previous works associated with shear heating can reconcile: (i) the simultaneous increase of pressure and temperature conditions within the lower unit; (ii) the metamorphic temperature peaks recorded within the thrust zone as maximum temperatures reached and representative of the thermal steady state; (iii) the order of thermal inversion intensity and (iv) the location of such a thermal inversion mainly below and across the thrust zone. Furthermore it does not contradict the possible concomittance of accretion processes and suface denudation that, in a very lesser extent compared to shear heating, modify the characteristics of the thermal inversion.

Thermal computations: heat diffusion equation ( 7) is solved by using the implicit finite difference method on the Eulerian grid. Then, the computed temperature difference on nodes ∆T between the current and the previous time step is interpolated from the grid to the particles. Finally, the new temperature is computed by adding the interpolated ∆T to the previous particle temperature. By ensuring the advection by particles, this method has the advantage of defining precisely the new temperature field directly on the particles.

Appendix B: Numerical model validation [START_REF] Molnar | Temperatures, Heat Flux, and Frictional Stress Near Major Thrust Faults[END_REF] and [START_REF] England | The interpretation of inverted metamorphic isograds using simple physical calculations[END_REF], by using a series of simplifications, proposed mathematical derivations of the fundamental equations controlling the heat transfer in the lithosphere. This way, without additional heat sources, they defined the tectonic conditions for the occurrence of a thermal inversion on each point of the thrust such as:

Z f ⋅ V th ⋅ sin(θ ) π ⋅ κ ≥1 (B.1)
where Z f is the depth of the considered point on the fault. This can be used to provide the minimum depth Z min where thermal inversion appears:

Z min = π ⋅ κ V th ⋅ sin(θ ) = π ⋅ κ Vz (B.2)
Considering a point in the downgoing plate located at the thrust front, it will reach the depth Z min after a time tz min : In spite of an overestimation of the minimum depth Z min of about 5 km, the slope of the linear regression law between Z min and Vz -1 is very close to the coefficient π•κ from Eq. (B.2) (Fig.

tz min = Z min Vz = π ⋅ κ Vz 2 (B.

B.1b

). In the same way, the numerical times t min show also a good agreement with the trend defined by the analytical solutions for tz min (Eq. (B.3) and Fig. B.1c,d). Shifts between the analytical laws and linear regressions are essentially due to the smoothed thrust profile we used (Fig. 1).

Advection effects on the thermal evolution within a thrust zone can be thus described as a function of the vertical component of the thrust velocity (Fig. 1). In Maximum temperatures are drawn for values of 500, 600 and 700˚C (continuous lines).

Maximum inverted geothermal gradients are represented from 10 to 50˚C.km -1 (dashed lines).

Dark green areas represent the combinations ( V th , η) for which both thermal characteristics are in the range of values considered as "realistic". Hatched areas indicate the combinations ( V th , η) for which no thermal inversion occurs. (continuous lines). The numerical simulations were realised at a convergence velocity V th set to 2 cm.yr -1 , with a thrust dip angle θ = 15˚ and a mantle heat flux q base = 40mW.m -3 . The grey areas correspond to the ranges of inverted thermal gradients and metamorphic peak temperatures characterising natural intracontinental inverted metamorphic sequences (see sections 3.5 and 3.6 for details about the radiogenic, accretion and erosion conditions).

Hatched area in (a) delimits the range of geothermal gradients not inverted. off. In order to highlight the thermal characteristics differences between both models, computations were realised under extreme kinematic conditions: V th = 5 cm.yr -1 and θ = 30˚.

[top] Snapshots and corresponding isotherms (in ˚C) of both models after a 1 Myr convergence (amount of shortening S = 50 km), i.e., when the thermal inversion reached its paroxysm.

[bottom] Differences of characteristic values of the thermal evolution along the thrust between the "Collision mode" and the "Subduction mode": temperature, vertical thermal gradient and thickness of the thermal inversion, respectively from left to right. For each characteristic value, the difference between both models is minor compared to the absolute values. imposed velocity field is circular in order to make the disc rotate around its center. Time step is defined according to the Courant criterion applied on the maximum local velocity, i.e., on the disc border. For example, a Courant criterion of 3 thus maximises the advection to a thirdcell by time step. A reference marker is localised on this border. Velocity field on particles was computed according to four different modes: on the one hand, direct computation is only characterised by linearly interpolating the velocity field from the nodes to the particles, while on the other hand, advection is realised according to the Runge-Kutta approximations of order 2, 3 or 4. Tests were also realised for a range of particles amount, from 1 to 100 per cell.

Number of particles does not impact on the effective advection. b) Offset after a complete revolution. Every test leads to a certain numerical error, i.e., the reference marker does not keep constant its distance to the disc center. While direct computation is characterised by a terminal offset of about 10 to 40 %, the Runge-Kutta methods allow a significant improvement of the particles advection leading to total offset of about 10 -2 to 10 -3 % for the order 2 and 3, and for less than 10 -4 % for the order 4. In the present test, the offset is positive for the direct computation and the 2 nd order of Runge-Kutta methods, i.e., the disc tends to dilate, whereas it is negative for the other two Runge-Kutta orders. This explains why the 2 nd and 3 rd orders are very similar in absolute. The Courant criterion also plays a main role in numerical advection. In the present test, increasing it with a factor 2 assures the reduction of the space error (offset) of a factor 2 for the simple computation, of a factor 8 for the Runge-Kutta methods of order 2 and 3, and of a factor 32 for the 4 th order. c) Relative time of normalised to the maximum time (i.e., for the 4 th order of the Runge-Kutta approximation, with a Courant criterion of 4 and 100 particles/cell). Results are given in percents. Time of computation linearly increases with the Courant criterion, with the amount of particles and also with the order of the used Runge-Kutta method. Even if there is an important gap between the calculation times of the direct computation and the 2 nd order of the Runge-Kutta approximation, the benefit of space advection precision significantly prevails.

The present test reveals an advection error (offset) of only 0.0098 % and an associated computation time more than 2 times lower than for the 3rd order and approximately 4 times better than for the 4 th order. Consequently, these observations lead us to use a Courant criterion of 2 associated with the 2 nd order of the Runge-Kutta method for a good compromise between the advection approximation and the computation time. (background). h) Numerical t min as a function of κ.Vz -2 (coloured circles). The linear regression is represented by the continuous line. The dashed straight line corresponds to the analytical linear law after [START_REF] England | The interpretation of inverted metamorphic isograds using simple physical calculations[END_REF]. Steady state is more rapidly reached in the shallowest part of the thrust, and thermal inversion duration in depths does not exceed 6 Myrs.

2 3 4 k [W.m -1 .K -1 ] k [W.m -1 .K -1 ] T [˚C] a b 1 2 3 4 1 2 3 4 k [W.m -1 .K -1 ] [b ] [a]

Figure 7

no inversion no inversion no inv.

no inversion no inversion no inv.

  inversion in time beyond the first 6 Myr (see Appendix B, Fig. B.2).

  3) Whereas England and Molnar considered θ constant along the entire slanted part of the fault, θ represents here the maximum dip angle of the thrust and an exact fit with the analytical solutions (B.2) and (B.3) cannot be expected and the thermal inversion is slightly underestimated.The numerical minimum depth min and start time t min of the thermal inversion as a function of V th and θ show a good agreement with the general trend of the analytical values given by or deduced fromEngland and Molnar (1993) (Fig. B.1a, c). This consistency is confirmed by the linear regressions (Fig.B.1b, d) resulting from the set of numerical values.

  Fig. B.1e, g, the results for a range of numerical simulations with different (Vz, k) combinations are treated in the same way as previously. Even if the slopes a relative to the linear regression laws corresponding to the minimum depth Z min and time t min (Fig. B.1f, h) of the beginning of the thermal inversion are rather far from the coefficient π predicted by England and Molnar (1993) (Eq. (B.2) and (B.3)), both the high similarity between our two regression slopes a and the very good correlation coefficients r confirm the tendencies of the analytical equations.

Figure Captions : Fig. 1 .

 :1 Figure Captions:

Fig. 2 .

 2 Fig. 2. Thermal evolution along the thrust as a function of the vertical component Vz of the

Fig. 3 .

 3 Fig. 3. Thermal evolution along the thrust as a function of the vertical component Vz of the

Fig. 4 .

 4 Fig. 4. Empirical thermo-dependence laws for k for various lithologies. a) Seipold (1998).

Fig. 5 .

 5 Fig. 5. a) Influence of the shear heating Hs on the spatial evolution of the geothermal gradient

Fig. 6 .

 6 Fig. 6. Impact of shear heating on the thermal evolution along the thrust during the

Fig. 7 .Fig. 8 .

 78 Fig. 7. a) Maximum thermal inversion intensity reached beneath the thrust zone in presence of

Fig. 9 .

 9 Fig. 9. a) Maximum thermal inversion intensity and b) maximum temperature from the thrust

Fig. A. 1 .

 1 Fig. A.1. Comparison between two velocity field configurations: "Collision mode", in which

Fig. A. 2 .

 2 Fig. A.2. Influence of the computation mode of the velocity field on particles mode and of the

Fig. B. 1 .

 1 Fig. B.1. Role of the kinematic conditions (V th and θ) and the thermal conductivity k:

Fig. B. 2 .

 2 Fig. B.2. Thermal evolution for an extreme kinematic setting: V th = 5 cm.yr -1 , θ = 30˚. a)

Fig. C. 1 .

 1 Fig. C.1. Rock thermal properties and mantle heat flux effects on the thermal inversion under extreme kinematic conditions (V th = 5 cm.yr -1 , θ = 30˚). a) and e) Heat capacity Cp effects. b)

Fig. C. 2 .

 2 Fig. C.2. Radiogenic heat production effects on the thermal inversion under extreme kinematic conditions: V th = 5 cm.yr -1 and θ = 30˚. a) and d) Surface radiogenic heat

  Figure 1

  inversion: 6.01 ˚C.km-1 Maximum inversion: 9.98 ˚C.km-1 Maximum inversion: 4.65 ˚C.km-1 Maximum inversion: 8.17 ˚C.km-1 Maximum inversion: 4.59 ˚C.km-1 Maximum inversion: 8.62 ˚C.km -

  Figure 4

  Figure 5

  Figure 8

Figure 13 :Combined

 13 Figure 13: Appendix Figure B.2

Table Captions :Table 1 .

 Captions1 Parameters tested in the model.

Acknowledgments

We thank L. Le Pourhiet, T. Duretz and D. May for fruitful discussions and their numerical advices. Stimulating discussions with M. Ballèvre and P. Gautier at the onset of this work were greatly appreciated. We also thank P. Tackley and D. May for the material support (ETH Brutus cluster). This paper benefited from thoughtful and careful revisions from S. Guillot, S. Schmalholz and an anonymous reviewer.

Figure 3 Maximum thickness: 15.96 km Maximum thickness: 12.63 km Maximum thickness: 9.92 km Maximum thickness: 15.82 km Maximum thickness: 17.64 km Maximum thickness: 13.95 km Maximum thickness: 3.77 km Maximum inversion: 2.02 ˚C.km -1 Maximum inversion: 14.71 ˚C.km -1 Maximum inversion: 6.01 ˚C.km -1 Maximum inversion: 2.97 ˚C.km -1 Maximum inversion: 5.82 ˚C.km -1 Maximum inversion: 9.97 ˚C.km -1 Maximum inversion: 19.89 ˚C.km -1 Maximum thickness: 15.96 km Maximum thickness: 17.6448 km Maximum thickness: 13.9501 km Maximum inversion: 9.97 ˚C.km -1 Maximum inversion: 6.67 ˚C.km -1 Maximum thickness: 15.80 km Maximum thickness: 15.80 km Maximum inversion: 13.34 ˚C.km -1 Figure 6 Maximum thickness: 12.70 km Maximum thickness: 8.67 km Maximum thickness: 11.75 km Maximum thickness: 3.57 km Maximum thickness: 7.59 km Maximum thickness: 16.74 km Maximum thickness: 3,7684 km Maximum thickness: 18.72 km Maximum thickness: 12.65 km Maximum thickness: 21.76 km Maximum inversion: 10.49 ˚C.km -1 Maximum inversion: -5.77 ˚C.km -1 Maximum inversion: 3.11 ˚C.km -1 Maximum inversion: 4.81 ˚C.km -1 Maximum inversion: 2.03 ˚C.km -1 Maximum inversion: 0.66 ˚C.km -1 Maximum inversion: 16.87 ˚C.km -1 Maximum inversion: 9.81 ˚C.km -1 Maximum inversion: 5.85 ˚C.km -1

Maximum inversion: 0.59 ˚C.km -1 Maximum thickness: 2.60 km Maximum thickness: 7.98 km Maximum inversion: 2.77 ˚C.km -1

Appendix A: Computation methods complements

Continental collision vs. continental subduction: The velocity field configuration of the downgoing plate beneath the first 30 km has a negligible influence on the thermal evolution around the thrust (Fig. A.1).

Interpolation computations between nodes end particles: They are realised by using a weight interpolation scheme [START_REF] Gerya | Introduction to numerical geodynamic modelling[END_REF] between the Lagrangian particles and the Eulerian nodes. The interpolation from particles to nodes is realised on 1 cell by using at least 25 particles per cell. This constitutes a good compromise between the time needed for interpolation calculations and the generated numerical error [START_REF] Yamato | Numerical modelling of magma transport in dykes[END_REF]. Advection computations: At each time step, the velocity field corresponding to the tectonic advection (both lp and V up ) is computed on each particle. The velocity direction of each particle is defined as parallel to the nearest segment of the thrust. In order to avoid a sharp variation of the velocity field on particles, the thrust profile was intentionally smoothed at the border of its slanted part (Fig. 1). The magnitude and the sign of each particle velocity are defined by its location relatively to the thrust. Particles are then advected by using the Runge-Kutta method. A Courant criterion maximising the advection to a half-cell by time associated the second-order Runge-Kutta method leads to an accurate advection while limiting the computation time (Fig. A.2). The differences between the fourth, third and second orders of the Runge-Kutta method are negligible.

Erosion computation: To simulate erosion, each particle advected above the topography is then numerically considered as part of the atmospheric domain by changing each one of the properties it carries. This way, erosion does not perturb the topography and consequently the crustal thickness that controls the vertical displacements linked to the isostatic balancing. The mass balance is assured by adding the corresponding amount of material laterally to the upper plate.

Thermal conductivity k thermo-dependence laws as a function of the lithology, from [START_REF] Zoth | Handbook of Terrestrial Heat Flows Density Determination[END_REF] and Clauser and Huenges (1995) (T in ˚C):