P Agard

P Yamato

M Soret

C Prigent

S Guillot

A Plunder

B Dubacq

A Chauvet

P Monié

Plate interface rheological switches during subduction infancy: control on slab penetration and metamorphic sole formation

Keywords: Subduction, metamorphic sole, rheology, plate interface, slab dehydration, mechanical coupling

Subduction infancy corresponds to the first few million years following subduction initiation, when slabs start their descent into the mantle. It coincides with the transient (yet systematic) transfer of material from the top of the slab to the upper plate, as witnessed by metamorphic soles welded beneath obducted ophiolites. Combining structure-lithology-pressure-temperature-time data from metamorphic soles with flow laws derived from experimental rock mechanics, this study highlights two main successive rheological switches across the subduction interface (mantle wedge vs. basalts, then mantle wedge vs. sediments; at ~800°C and ~600°C, respectively), during which interplate mechanical coupling is maximized by the existence of transiently similar rheologies across the plate contact. We propose that these rheological switches hinder slab penetration and are responsible for slicing the top of the slab and welding crustal pieces (high-then low-temperature metamorphic soles) to the base of the mantle wedge during subduction infancy. This mechanism has implications for the rheological properties of the crust and mantle (and for transient episodes of accretion/exhumation of HP-LT rocks in mature subduction systems) and highlights the role of fluids in enabling subduction to overcome the early resistance to slab penetration.

Introduction

Understanding subduction initiation, in both space and time, has been a challenge since the advent of plate tectonics [START_REF] Dewey | Ophiolite obduction[END_REF][START_REF] Regenauer-Lieb | The Initiation of Subduction: Criticality by Addition of Water?[END_REF][START_REF] Gurnis | Evolving force balance during incipient subduction[END_REF]. What is referred to as "subduction initiation" in the literature encompasses two different concepts and periods: (i) how and where subduction nucleates (i.e., what triggers the beginning of subduction; e.g., [START_REF] Regenauer-Lieb | The Initiation of Subduction: Criticality by Addition of Water?[END_REF][START_REF] Stern | Subduction initiation: spontaneous and induced[END_REF], and (ii) how subduction proceeds over the first few million years of its history ("subduction infancy"; [START_REF] Stern | Subduction zone infancy: examples from the Eocene Izu-Bonin-Mariana and Jurassic California arcs[END_REF].

This study focuses on subduction infancy, when a newly born slab starts its descent into the mantle and when the thermal regime of the subduction zone progressively cools down before reaching steady-state (e.g., [START_REF] Syracuse | The global range of subduction zone thermal models[END_REF]Plunder et al., 2015;Figs. 1a,b). The only rock remnants of this elusive geodynamic step are thin (~10-500 m) metamorphosed slivers of oceanic crust (metamorphic soles; [START_REF] Williams | Metamorphic aureoles beneath ophiolite suites and Alpine peridotites: tectonic implications with west Newfoundland examples[END_REF][START_REF] Wakabayashi | Spatial and temporal relationships between ophiolites and their metamorphic soles: A test of models of forearc ophiolite genesis[END_REF] found beneath pristine, 100-1000 km long, ≤10-15 km thick fragments of oceanic lithosphere emplaced on top of continents as ophiolites [START_REF] Coleman | Tectonic setting for ophiolite obduction in Oman[END_REF][START_REF] Nicolas | Structures in Ophiolites and Dynamics of Oceanic Lithosphere (Petrology and Structural Geology)[END_REF]Fig. 1c).

Metamorphic soles correspond to upper crustal material from the downgoing slab (with variable proportions of basalts and pelagic sediments; [START_REF] Spray | Possible causes and consequences of upper mantle decoupling and ophiolite displacement[END_REF][START_REF] Boudier | Shear zones, thrusts and related magmatism in the Oman ophiolite : Initiation of thrusting on an oceanic ridge[END_REF] and have long been recognized as formed during the first few My of intra-oceanic subduction (Fig. 1b; [START_REF] Dewey | Ophiolite obduction[END_REF][START_REF] Spray | Possible causes and consequences of upper mantle decoupling and ophiolite displacement[END_REF][START_REF] Dewey | The sole of an ophiolite: the Ordovician Bay of Islands Complex[END_REF]. Their formation would result from heat transfer from the upper plate mantle and/or shear heating when the slab enters the mantle and heats up [START_REF] Dewey | Ophiolite obduction[END_REF][START_REF] Hacker | Simulation of the Metamorphic and Deformational History of the Metamorphic Sole of the Oman Ophiolite[END_REF]. Explaining how such thin metamorphosed tectonic slivers of oceanic crust get welded ("underplated") to the upper plate along hundreds of km (e.g., Oman, Turkey: Hacker and Gnos, 1997;[START_REF] Çelik | Early-Middle Jurassic intra-oceanic subduction in the İzmir-Ankara-Erzincan Ocean[END_REF] is essential for understanding mechanical coupling during subduction infancy (and possibly during later subduction), but has so far remained enigmatic (Jamieson, 1981;[START_REF] Dewey | The sole of an ophiolite: the Ordovician Bay of Islands Complex[END_REF]. This problem is herein addressed by (i) compiling worldwide characteristics of metamorphic soles (i.e., lithologies, internal organization, thicknesses, thermobarometric constraints), augmented by refined estimates for their pressure-temperature (P-T) conditions of formation using thermodynamic modelling and by (ii) calculating effective viscosities of materials present along the plate interface from known rheological properties for the crust and mantle (i.e., peridotite, basalt, sediment, serpentinite).

This study reveals the existence of rheological switches across the subduction interface, and proposes that these changes in rheological properties control slab penetration into the mantle and the formation of metamorphic soles during subduction infancy. This mechanism has implications for effective rheologies of the crust and mantle and for the general understanding of accretion processes and early slab dynamics.

Metamorphic soles: the record of subduction infancy

Metamorphic sole constitution

The main characteristics (i.e., structural position, lithologies, constitution) and P-T conditions of metamorphic soles worldwide are reviewed in figure 2 and Table 1. This synthesis shows that metamorphic soles are ubiquitous beneath non-metamorphosed ophiolites (e.g., Oman, Turkey, Papua, Newfoundland) and share similar characteristics regardless of the ophiolite or the detailed geological/geodynamical setting [START_REF] Spray | Possible causes and consequences of upper mantle decoupling and ophiolite displacement[END_REF][START_REF] Wakabayashi | What constitutes «emplacement of an ophiolite ?: Mechanisms and relationship to subduction initiation and formation of metamorphic soles[END_REF]. Radiometric ages of metamorphic soles and ophiolites generally fall within 1-2 My (e.g., Fig. 1d for Oman; Hacker et al., 1996;[START_REF] Rioux | Tectonic development of the Samail ophiolite : High-precision U-Pb zircon geochronology and Sm-Nd isotopic constraints on crustal growth and emplacement[END_REF], suggesting the existence of a still warm (i.e., > 1000°C) upper plate mantle near the subduction interface.

Metamorphic soles comprise ~10 to ~500 m thick (Figs. 2b-c) highly strained and metamorphosed crustal rocks where amphibolitized metabasalt dominates, together with increasing proportions of pelagic metasediment structurally downwards (mainly metaradiolarite, with intercalations of metatuff and metapelite downwards). Vertically, metamorphic soles exhibit an inverted metamorphic sequence with isograds subparallel to the basal peridotite foliation [START_REF] Spray | Possible causes and consequences of upper mantle decoupling and ophiolite displacement[END_REF]. They grade steeply from thin high temperature (HT) granulite/amphibolite facies lithologies adjacent to the overlying peridotites (> 700-850°C;Fig. 2d;e.g., McCaig, 1983;Jamieson, 1986) to thicker amphibolite/greenschist facies low temperature soles (LT; ~550-650°C; Table 1). This temperature trend is not continuous, however, since the structurally and thermally lower LT sole is arguably formed later, and at lower pressure than the HT sole, by successive stacking of increasing amounts of metasediment (e.g., Malpas, 1979;[START_REF] Casey | Initiation of subduction zones along transform and accreting plate boundaries, triple-junction evolution, and forearc spreading centres: implications for ophiolitic geology and obduction[END_REF]Jamieson, 1986;Fig. 2d).

Whenever radiometric constraints are available, HT soles are coeval or slightly older than LT soles, yet within ~2 My (e.g., for Oman: Fig. 1d; Hacker et al., 1996;[START_REF] Roberts | Geochronological constraints on the metamorphic sole of the Semail ophiolite in the United Arab Emirates[END_REF]. HT soles and/or LT soles may be missing in places but wherever both are observed, and not disturbed by obvious later tectonics, HT soles are overlying LT soles.

High deformation in the HT and LT soles is marked by mylonites and complex recumbent folding. This deformation, however, is commonly less conspicuous in HT sole mafic amphibolites than in the LT soles (Jamieson, 1981; this study), due to the extent of recrystallization (e.g., Oman) and/or to the lack of lithological heterogeneities. Wherever (rarely) preserved, stretching lineations in the HT soles and LT soles strike differently (e.g., Newfoundland; [START_REF] Dewey | The sole of an ophiolite: the Ordovician Bay of Islands Complex[END_REF], suggesting that boundary conditions and/or accretion dynamics may have been modified during sequential underplating. Some earlier workers, on the basis of rare gabbroic occurrences (associated with dunite and locally intercalated between the mantle and the HT metamorphic amphibolite; Jamieson, 1981), suggested that the whole metamorphic sole could represent a metamorphosed, overturned limb of ophiolite crust (with sediments, basalts and gabbros from bottom to top; see also [START_REF] Wakabayashi | What constitutes «emplacement of an ophiolite ?: Mechanisms and relationship to subduction initiation and formation of metamorphic soles[END_REF]. This interpretation is however unlikely: (i) in contrast with the ophiolite crust, gabbros are extremely rare in the soles (and may represent small-scale intrusions in the mantle), (ii) HT mafic amphibolites tend to have a distinctive nature/geochemical signature (i.e., MORB transitional to OIB or E-MORB; [START_REF] Dewey | The sole of an ophiolite: the Ordovician Bay of Islands Complex[END_REF] and (iii) the overturned limb hypothesis fails to explain why pressure conditions in the metasedimentary LT sole are lower than in the HT sole located above (Gnos, 1998; section 2.2).

The mantle rocks immediately above the metamorphic sole are also highly deformed (Fig. 2d), showing m-to hm-scale deformation patterns consistent with those observed in the underlying HT sole [START_REF] Boudier | Shear zones, thrusts and related magmatism in the Oman ophiolite : Initiation of thrusting on an oceanic ridge[END_REF], pressure estimates equivalent to those of HT sole peak metamorphism (Jamieson, 1981;McCaig, 1983) and porphyroclastic to ultramylonitic textures formed in the temperature range of 1100°C down to ~700°C [START_REF] Boudier | Shear zones, thrusts and related magmatism in the Oman ophiolite : Initiation of thrusting on an oceanic ridge[END_REF][START_REF] Michibayashi | The role of pre-existing mechanical anisotropy on shear zone development within oceanic mantle lithosphere: an example from the Oman ophiolite[END_REF]Linckens et al., 2011b). This suggests that the base of the ophiolite mantle deformed and cooled during subduction infancy and that, from a mechanical point of view, the metamorphic sole should be considered as a threefold stack with, from bottom to top, the LT sole, the HT sole and the base of the ophiolitic mantle sequence (hereafter noted as: LTsole\HTsole\basal peridotites).

P-T conditions of metamorphic soles

Published P-T estimates for metamorphic sole formation (Fig. 3a; Table 1) spread along a high to medium T/P gradient, which partly arises from the diversity and variable precision of thermobarometric methods used. New P-T estimates for HT soles are provided here using thermodynamic modelling. P-T phase diagrams for fixed bulk rock composition (pseudosections) were calculated to constrain the conditions of HT sole formation using the Gibbs-free-energy minimization software THERIAK/DOMINO (de Capitani and Petrakakis, 2010; with the updated database of Holland and Powell (1998); tcdb55cc2d.bs) with the following solution models [START_REF] Diener | A new thermodynamic model for clino-and orthoamphiboles in the system Na2O-CaO-FeO-MgO-Al2O3-SiO2-H2O-O[END_REF] for amphibole, [START_REF] Green | An order-disorder model for omphacitic pyroxenes in the system jadeite-diopside-hedenbergite-acmite, with applications to eclogitic rocks[END_REF] for clinopyroxene, [START_REF] White | Progress relating to calculation of partial melting equilibria for metapelites[END_REF] for orthopyroxene, Holland et al. (1998) for chlorite, [START_REF] Baldwin | Modelling of mineral equilibria in ultrahigh-temperature metamorphic rocks from the Anápolis-Itaucu Complex, central Brazil[END_REF] for plagioclase and Holland and Powell (1998) for garnet. P-T conditions for LT soles are difficult to assess, unfortunately, owing to the high variance of the assemblages and uncertainties in thermodynamic models.

A representative pseudosection (Fig. 3b) was calculated within the chemical system Na 2 O-CaO-FeO-MnO-MgO-Al 2 O 3 -SiO 2 -H 2 O for a typical mafic amphibolite (YE1302b; Lycian ophiolite, W. Turkey) with garnet, clinopyroxene, plagioclase and amphibole. The corresponding sample comes from the region of Salda (South Western Turkey; N037°46'44" E029°56'21"), where garnet-clinopyroxene amphibolites are passing downwards, away from the contact with the overlying peridotite, to amphibolites and then to greenschist facies rocks. Spinel is replaced by plagioclase in the nearby peridotite, suggesting former equilibration of the rock at pressures > 0.8-1 GPa. At the thin section scale, garnet and clinopyroxene are intimately intergrown (with globular shape inclusions of clinopyroxene in garnet). Amphibole and plagioclase are found both in the matrix and as inclusions in garnet and clinopyroxene, and therefore formed during peak conditions. Clinopyroxene has diopside compositions with Mg# ~ 0.80 and a Ca content of 0.88 -0.91 per formula unit (p.f.u.). The pyrope content in garnet ranges between 0.20-0.28 p.f.u. and plagioclase (when preserved) has an anorthite fraction of 0.2-0.3. P-T conditions for the garnet-clinopyroxene-amphibole-plagioclase peak assemblage of sample YE1302b, using mineral isopleths, are 1.08 ± 0.1 GPa and 780 ± 40 °C. Water amounts were set to match amphibole modes (50-60 vol%). TiO 2 was neglected in the calculation as it enters mainly accessory minerals (rutile, titanite) and amphibole, in which titanium is not accounted for [START_REF] Diener | A new thermodynamic model for clino-and orthoamphiboles in the system Na2O-CaO-FeO-MgO-Al2O3-SiO2-H2O-O[END_REF]. The MnO content reproduces the observed garnet chemistry and ensures consistent garnet and clinopyroxene Fe-Mg exchange.

A similar range of estimates was obtained with phase diagrams for metamorphic soles from Turkey (in Kütahya; boxes in Fig. 3a; Table 1) and Oman (in Sumeini and Khubakhib; dashed boxes, Fig. 3a), pointing to definitely high P values for the upper HT soles when compared to the spread of published estimates (i.e., black dots in Fig. 3a). This conclusion is strengthened by the similarity of available mineral assemblages and compositions worldwide. These estimates are consistent with the presence of only subordinate amounts of melts in the HT soles (Gnos, 1998), suggesting that temperatures do not significantly exceed amphibolite dehydration melting (in agreement with predicted melt fractions < 5-10 vol% at 850°C, depending on pressure; [START_REF] Green | Activity-composition relations for the calculation of partial melting equilibria for metabasic rocks[END_REF].

Significance of metamorphic soles within thermal subduction regimes

The above compilation shows that conditions for metamorphic sole formation are remarkably similar worldwide and characterized by the accretion of HT soles that are thinner, more mafic, accreted earlier, at greater depths, at almost invariant P-T conditions (800 ± 50°C at 1.0 ± 0.2 GPa) and always on top of LT soles (equilibrated at ~600 ± 50°C at 0.5 ± 0.1 GPa).

These estimates are compared to peak conditions reached by oceanic rocks during later subduction, namely those for HT and LT eclogites (Fig. 3a; see [START_REF] Agard | Exhumation of oceanic blueschists amd eclogites in subduction zones: Timing and mechanisms[END_REF] for a review).

Their contrasting P-T conditions exemplify the change in the subduction thermal gradient through time, from warm to cold, from peak burial conditions of ~800°C and 1.0 GPa to ~550°C and 2.5 GPa. Figure 3a shows that during subduction infancy and subsequent cooling, three rock types successively form from the upper crust of the down-going slab: (1) metamorphic soles within the first 1-2 My, along the hotter gradient, (2) HT oceanic eclogites postdating initiation by ~5 My, commonly exhumed in serpentinite mélanges with counter-clockwise P-T paths (arrows in Fig. 3a; e.g., Wakabayashi, 1990;Garcia-Casco et al., 2006), (3) LT oceanic eclogites formed after 5-10 My and, whenever exhumed, mostly as continuous tectonic slices.

Contrary to metamorphic soles, the fraction of eclogites exhumed worldwide is only rarely accreted to an upper plate ophiolite mantle. Figure 3a also shows that accretion largely overlaps the domain where serpentine is not stable.

Most HT metamorphic soles do not show later HP-LT metamorphism overprint, such as would be expected from progressive cooling if these rocks had remained at their depth of formation. The implication is that the HT soles formed at ~1 GPa must have been accreted and partly exhumed (together with the deformed base of the ophiolite mantle) during subduction infancy (i.e., before significant cooling of the subduction thermal regime), and later underplated below the undeformed ophiolite mantle.

Rheology of the plate interface and mechanical coupling during subduction infancy

A generic mechanical process is needed to explain how slices of crustal material from the slab get accreted to the mantle wedge (i) only during subduction infancy and (ii) in such a uniform manner. Accretion (or "underplating") of any tectonic slice across the subduction interface requires the combination of:

(i) an increase in mechanical coupling across the plate contact beyond some threshold, in order to preferentially localize strain and relative displacement elsewhere within the slab, along some other physical discontinuity (e.g., basalt vs. sediment, or basalt vs. sheeted dykes; [START_REF] Kimura | Peeling oceanic crust in subduction zones[END_REF][START_REF] Dewey | The sole of an ophiolite: the Ordovician Bay of Islands Complex[END_REF], (ii) some deformation mechanisms allowing for effective slicing within the slab. Whether slicing takes place through slow creep, fluid-mediated slip (such as slow slip events) or repeated regular earthquakes is unknown to date and beyond the scope of the present study.

Noteworthily, mechanical coupling is maximum when rheologies on both sides of the plate interface are similar: strain is otherwise localized and the interface is decoupled.

Effective viscosities of plate interface material constrained by rock mechanics

Flow laws derived from experimental rock mechanics (Table 2) are used to estimate the effective viscosity (of plate interface materials as a function of temperature (T). Figure 4a shows this dependency for rocks expected to lie at the base of the upper plate (peridotite, serpentinite) and at the top of the slab (basalt, sediment), using a strain rate ( ˙ e ) derived from natural constraints (10 - 13 s -1 for Oman; Linckens et al., 2011b) and the following formula:

where n, A, Q and R correspond to the power-law exponent, the material constant, the activation energy and the gas constant, respectively.

The reader is referred to [START_REF] Karato | Rheology of the Earth's mantle: A historical review[END_REF] and Hirth and Kohldstedt (2015) for reviews of flow laws of relevance for the mantle. These were derived from experiments on dry or wet olivine for dislocation creep (e.g., Hirth and Kohlstedt, 2003) and for grain size dependent deformation mechanisms such as diffusion creep (Hirth and Kholstedt, 2003;Faul and Jackson, 2007) and dislocation-accommodated grain boundary sliding (disGBS; Hirth andKohldstedt, 2003: Hansen et al., 2011). Important independent constraints on mantle deformation along the plate interface come from natural data on basal peridotites and mantle shear zones rooting in the deformed mantle base of the ophiolite [START_REF] Boudier | Shear zones, thrusts and related magmatism in the Oman ophiolite : Initiation of thrusting on an oceanic ridge[END_REF]Linckens et al., 2011a,b). The evolution from 1100°C in porphyroclastic peridotites to ~700°C in ultramylonites shows that (i) olivine grain size decreases from ~2 mm to 10-50 m and (ii) deformation mechanisms evolve at ~800-850°C from dislocation creep to grain size sensitive creep, diffusion creep being the most likely dominant mechanism in localized shear bands (the boundary between diffusion creep and disGBS is not well known however; Linckens et al., 2011b).

Flow laws for metamorphic sole mafic rocks are scarce in comparison and can be approached by using wet/dry diabase and mafic granulite. Flow laws are even scarcer for metasediments (Table 2). Given the abundance of metaradiolarites in the LT soles, the flow law of quartz is considered here as representative. Felsic granulite (dashed purple curve; Fig. 4a) may represent an equivalent to a strongly metamorphosed metasediment.

Importantly, at any given temperature and (at least as a first order approximation) regardless of which flow law is used, the mechanical strength of peridotite (i.e., dry or wet olivine: thick green curves and green overlay, respectively; Fig. 4a) is greater than that of mafic oceanic crust (blue overlay in Fig. 4b). Sediments are weaker than basalt, yet stronger than serpentinite. This is emphasized in figure 4b, where averages of flow laws are shown for the mantle (i.e., dislocation creep for dry and wet olivine and diffusion creep or disGBS with a grain size of 30 m), for the mafic crust ("basalt") and for serpentinite. It is emphasized that changing the strain rate modifies the absolute value of η but does not affect the relative position of these curves (see supplementary figure S1). The strength contrasts existing between these different materials therefore seem independent of strain rate.

Rheological switches across the plate interface

Slab dehydration (see Faccenda, 2014 for a review) can be anticipated to be critical for mechanical coupling during subduction infancy. Fluids released into the nascent mantle wedge (Fig. 4c) will induce serpentinization of the mantle below ~550-650°C (depending on pressure: Fig. 3a; [START_REF] Ulmer | Serpentine stability to mantle depths and subduction-related magmatism[END_REF]. These fluids will be stored in progressively lesser amounts deeper down, first as hydrous phases such as chlorite ± amphibole ± phlogopite ± talc, then at T > ~850-900°C as fluid inclusions or point defects in nominally anhydrous minerals such as olivine or pyroxene [START_REF] Hirth | The stress dependence of olivine creep rate: Implications for extrapolation of lab data and interpretation of recrystallized grain size[END_REF].

Contrary to the general increase of material strength with cooling (Fig. 4b), cooling of hydrated mantle wedge peridotites will therefore progressively weaken the mantle wedge towards serpentinite rheology (from point i to f; Fig. 4b) by (i) absorption of OH in olivine ("wet" olivine) at 1000-900°C and changes in deformation mechanisms from dislocation creep to grain size sensitive creep at T < ~850-900°C (i.e., diffusion creep and disGBS), (ii) formation of weaker hydrated minerals and eventually (iii) serpentinization. Noteworthily, the viscosity of an even mildly serpentinised peridotite (> 15%) approaches that of pure serpentine [START_REF] Escartin | Strength of slightly serpentinised peridotites: implications for the tectonics of oceanic lithosphere[END_REF].

Although the exact rheological path from peridotite to serpentinite cannot be precisely quantified presently, an important conclusion is that mantle wedge viscosities will cross over the curves for metabasalt and then metasediment (Fig. 4b): these two "rheological switches" (i.e., when the mantle wedge first becomes weaker than basalts, then than sediments) most probably occur at T ~800°C (where deformation by diffusion creep and disGBS have similar effective viscosities; Figs. 4a,b) and ~600°C, respectively. The temperature and viscosity values for the HT rheological switch (Fig. 4d) both match the temperature (T~800-850°C; [START_REF] Boudier | Shear zones, thrusts and related magmatism in the Oman ophiolite : Initiation of thrusting on an oceanic ridge[END_REF]Linckens et al., 2011a) and viscosity (10 20-21 Pa.s; Linckens et al., 2011b;[START_REF] Tasaka | Influence of mineral fraction on the rheological properties of forsterite+enstatite during grain size sensitive creep: 3. Application of grain growth and flow laws on peridotite ultramylonite[END_REF] inferred from the high strain mylonitic to ultramylonitic deformation of adjacent banded peridotites from the base of the ophiolite.

Model for slab penetration into the mantle and metamorphic sole formation

Evolution of rheological contrasts across the plate interface during subduction infancy

Figure 4d shows the evolution of the rheological contrast across the plate interface during incipient slab penetration, as the top of the slab progressively heats up, weakens and dehydrates.

The top of the slab is first considered as essentially made of mafic crust (i.e., basalt; Fig. 4c), sediments being probably scarce at the start of intra-oceanic subduction, far away from the continent (Fig. 1a). At shallow depths and for temperatures below ~550-600°C, the subducting oceanic (basaltic) crust is juxtaposed against an incipiently serpentinised mantle on top (as schematized in Fig. 4c). A sharp viscosity contrast (exists on either side of the subduction plane (strong basalt vs. weak mantle wedge; Fig. 4d).

As depth and temperature increase, the subducting crust weakens and progressively encounters a warmer, stronger, unserpentinised and less hydrated mantle wedge (' on Fig. 4d), where grain size sensitive deformation mechanisms take over (Linckens et al., 2011a;Hirth and Kohldstedt, 2015). The viscosity contrast across the interface reverses (weak basalt vs. stronger mantle wedge) at T > ~800°C, once the subducting crust is juxtaposed against an almost dry peridotite rheology (" on Fig. 4d).

A similar evolution can be envisioned if sediments are present on top of the slab, with a viscosity reversal (sediment vs. mantle wedge) taking place at ~600°C ("second rheological switch"; Fig. 4d). Although the extent to which basalts and sediments harden as a result of prograde mineral transformations is unknown, flow laws indicate that these rheological switches exist, even considering the change from basalt to mafic granulite or from sediment/quartzite to felsic granulite (Fig. 4a). They will also take place regardless of the age of the overriding lithosphere, although at different depths (the warmer the lithosphere, the shallower the rheological switches).

Rheological switches bear major consequences on mechanical coupling during early slab penetration:

(i) at any given time the strength of the mantle wedge will increase downwards (Fig. 4c), so that the slab can be expected to face greater resistance to penetration with depth (at least until mantle melting occurs);

(ii) as the thermal regime cools with time, the domain where serpentine is stable ("serp. front" in Fig. 4c) will expand continuously downwards.

Metamorphic sole formation linked to strong interplate mechanical coupling

Effective viscosities of the lower plate crust and upper plate mantle therefore converge and switch, during subduction infancy, across restricted T windows (Fig. 4d; with P and time generally ~ 1 GPa and < 2Myr, these are restricted P-T-t windows too). We propose that the detachment and accretion of metamorphic soles is triggered by peaks of interplate mechanical coupling associated with rheological switches, that distribute deformation over large, km-scale bands across the subduction interface (e.g., [START_REF] Yuen | Shear deformation zones along major transform faults and subducting slabs[END_REF] and localize strain further into the slab where/if a sharper viscosity contrast exists (e.g., [START_REF] Kimura | Peeling oceanic crust in subduction zones[END_REF].

Slab penetration and metamorphic sole formation are tentatively reconstructed in figure 5 in three major steps (isotherms and depth-time trajectories are from thermo-kinematic models detailed in supplementary material; Figs. S2-4):

(1) Plate interface mechanical strength peaks during the first rheological switch (metamorphosed basalt vs. metasomatized upper mantle; Fig. 5a) at T~750-850°C, leading to the formation of essentially mafic HT metamorphic soles (dot 1a, Fig. 4d). Detachment on a weaker horizon within the slab probably takes place at the transition between hydrothermalized/weakened basaltic layers and drier basalts below, and/or between basalts and sheeted dykes (if present), accounting for the general lack of metagabbros in metamorphic soles.

(2) Accretion of the metasedimentary-rich LT soles takes place at lower T (and P) conditions at ~550-650°C and ~0.5 GPa, after the partial exhumation of the HT sole (Fig. 5b). This corresponds to the second rheological switch (dot 1b, Fig. 4d). Two different scenarios of mechanical coupling can however be envisioned: (i) between the sediments and the incipiently serpentinised, weakening upper mantle (option 1; Fig. 5d) or (ii) between the sediments and the base of the HT sole (option 2; Fig. 5d), whose viscosity will be shifted from basalt to sediment as the sedimentary column increases (blue arrow, Fig. 4b; a similar shift is expected if strain rate increases or if basalt from the HT sole get hydrated by water released from sediments below). Option 2 is supported by the general lack of (serpentinised) mantle between the LT and HT sole (see section 2). Option 1 implies that the mantle above the LT sole is tectonically removed by the exhumation of the (more resistant)

HT sole and basal peridotites (Fig. 5d). Option 1 would nevertheless explain the occurrence of LT soles directly beneath the mantle. Accretion of the LT sole indicates, in any case, that the detachment horizon within the slab is located within the sedimentary pile.

The fact that HT soles are thinner than LT soles could result from larger amounts of accumulated strain and later ductile thinning, from the peeling of thinner slices during the first rheological switch and/or from longer duration of accretion during the LT episode.

(3) As the thermal regime of the subduction continues to decrease (Fig. 5c), incoming sediment and basalt remain stronger than the increasingly serpentinised mantle wedge (Fig. 4d): the plate interface progressively 'unzips' downwards. LT eclogites, which start forming within the refrigerated subduction zone (Fig. 5c), are less likely to get mechanically coupled to the weakened upper plate. This can be the reason why they are rarely exhumed [START_REF] Agard | Exhumation of oceanic blueschists amd eclogites in subduction zones: Timing and mechanisms[END_REF], in agreement with their location in the serpentine stability field ("no accretion" domain in figure 3a). (i) The mantle wedge acts as a buttress which progressively softens with time: it only transiently peels off the slab crust during the first My, thereby forming metamorphic soles, then progressively loses strength with cooling/serpentinization until full decoupling;

Impact of subduction cooling on interplate mechanical coupling

(ii) Strain localizes with time in shallower decoupling horizons within the slab (Fig. 5e): accretion affects the top of the mafic crust (stage a), then only the sediments on top (stage b) and, when the plate interface becomes decoupled after a few My (stage c), accretion is restricted to shallow, neartrench infill (e.g., unmetamorphosed Hawasina units found in Oman beneath the metamorphic sole and the ophiolite; Searle and Malpas, 1980;[START_REF] Searle | Subduction zone metamorphism during formation and emplacement of the Semail ophiolite in the Oman Mountains[END_REF], as in present-day accretionary wedges (e.g., Nankai prism). The width of the plate interface shear zone may therefore decrease with time.

Discussion

The proposed mechanism for metamorphic sole formation and slab penetration (Fig. 5) explains why, through two main steps of accretion, metamorphic soles form and get accreted at remarkably similar P-T-time conditions worldwide (section 2.3) and across a transient period of subduction lifetime only (i.e., when rheological switches take place).

Uncertainties on P-T-time conditions, viscosity estimates and strain rate

P-T-time constraints

P-T estimates (800 ± 50°C at 1.0 ± 0.1 GPa for the early thin HT sole and 600 ± 50°C at ~0.5 ± 0.1 GPa for the late, thicker LT sole; Fig. 3) depend on the accuracy of thermodynamic models for pyroxene and amphibole. These models were considerably refined but Ti-and Fe 3+ -substitutions in amphiboles are complex and mafic melts are notoriously difficult to model (Diener and Powell, 2012 and references therein;[START_REF] Green | Activity-composition relations for the calculation of partial melting equilibria for metabasic rocks[END_REF]. Larger uncertainties on pressure conditions for LT soles arise from lesser constraints in the high-variance greenschist to epidote amphibolite fields.

The major uncertainty therefore lies in the exact pressure gap between the HT and LT soles, and whether there might be a continuum in between.

Age constraints for HT and LT metamorphic soles tightly cluster within 1-2 My (Table 1; Fig. 1d). The short duration of the process and/or apparent synchronicity might be exaggerated by the fact that 40 Ar/ 39 Ar ages on amphibole may represent cooling ages for the HT soles (i.e., below 550°C), whereas age constraints for the LT sole could represent crystallization ages. But the fact that the LT soles are found below the HT ones (whom underlie the mantle base of the ophiolite) leaves little doubt that they were accreted afterwards. Thermo-mechanical modelling of a cooling subduction zone also shows that cooling lasts no more than a few My [START_REF] Duretz | Thermo-mechanical modeling of the obduction process based on the Oman Ophiolite case[END_REF].

Viscosities and strain rates

One of the largest uncertainties with experimental flow laws is that they are performed at (and extrapolated from) conditions orders of magnitude faster than nature (typically 10 -4 -10 -6 s -1 versus 10 -12 -10 -15 s -1 in nature; Burov, 2011; [START_REF] Hirth | The stress dependence of olivine creep rate: Implications for extrapolation of lab data and interpretation of recrystallized grain size[END_REF]. Field evidence, however, suggests that mafic, clinopyroxene and plagioclase bearing lithologies are ~2 orders of magnitude weaker than dry peridotite at ~700-800°C (e.g., [START_REF] Homburg | Investigation of the strength contrast at the Moho: A case study from the Oman Ophiolite[END_REF], in agreement with the respective position of calculated flow laws (Fig. 4b,c). Viscosity estimates for the first, HT rheological switch (~10 20-21 Pa.s; Fig. 4d) are also in remarkable agreement with values deduced from natural observations on strained peridotites from the (ultra)mylonitic base of the ophiolite mantle (10 20-21 Pa.s; Linckens et al., 2011b;[START_REF] Tasaka | Influence of mineral fraction on the rheological properties of forsterite+enstatite during grain size sensitive creep: 3. Application of grain growth and flow laws on peridotite ultramylonite[END_REF].

Another uncertainty comes from the simplifying assumption of considering similar strain rates for estimating the viscosities of all lithologies (Fig. 4b), whereas faster strain rates can be expected where strain localization takes place. Regardless of strain rate, however, "rheological switches" and converging mechanical behaviour between the mantle and basalt/sediment will take place across the interface: as the mantle wedge evolves progressively from dry to serpentinised, its viscosity curve is bound to cross that of basalt then of sediment (Fig. 4b,d), and this order is independent of strain rate (e.g., from 10 12 to 10 14 s -1 ; Fig. S1). This conclusion is not modified by shear heating either, which only changes the temperature field and shifts the T window of metamorphic sole formation to shallower depths, but not the characteristic "S-shape" of the isotherms (Fig. S5; so does the thermal state of the upper plate: Figs. S2b,c).

Complex feedbacks may nevertheless exist, since similar viscosities on each side will increase interplate mechanical coupling, which will in turn decrease strain rates, hence probably decrease shear heating (also depending on the width onto which deformation is distributed). Further testing of this scenario using self-consistent fully coupled thermomechanical models is therefore needed, but there are major challenges:

(i) refined rheologies for sediments, variably hydrated basalts and gabbros, taking into account the influence of hydrothermal alteration, progressive metamorphic recrystallization, water loss or plagioclase content on crustal rocks are unknown (e.g., [START_REF] Getsinger | Amphibole fabric formation during diffusion creep and the rheology of shear zones[END_REF]. The same holds true for the rheology of the mantle wedge (as yet unconstrained, as is fluid migration inside; [START_REF] Faccenda | Water in the slab: A trilogy[END_REF], for polyphase lithologies that also probably exist along the plate interface, for the impact of subordinate amounts of melt in the HT soles (<~10%), or for assessing the influence of pressure or the extent to which frictional energy is converted into heat (which also depends on rheological laws).

(ii) appropriate spatial resolution (i.e., down to ~10 m) is required in self-consistent viscoelasto-plastic geodynamic models in order to localize strain and (progressively) slice and detach pieces from the slab, in addition to reaching sufficient temporal resolution (e.g., van Dinther et al., 2014).

Further modelling will help constrain the depth and duration of these processes. The simple thermo-kinematic models used to derive the isotherms of figure 5 (see supplementary material) suggest that for a set of realistic velocities, initial thermal age and slab dip, incoming crustal rocks may cross the temperature range of HT sole formation (~750-850°C) at depths of ~25-35 km and that accretion of individual slices may last on the order of 0.3-0.4 My (Fig. S3). Although inferred from simplified models, these constraints point to the possible formation of metamorphic soles across a range of depths, possibly accounting for some of the scatter observed in the P-T estimates for metamorphic soles worldwide (Fig. 3a).

Accretion and exhumation during subduction infancy

Accretion of HT soles and incorporation of HT eclogites in serpentinite mélanges

Accretion of mafic HT soles will last less than shown in the reconstruction and probably stop (or decrease) once a significant amount of sediments reaches the trench and becomes involved in the plate interface (Fig. 5b): since incoming sediments reach their rheological switch with the mantle wedge at lower T than basalts (Figs. 4c,d), their arrival along the plate interface will indeed localize strain and thus deactivate basalt accretion.

As a result of progressive cooling of the subduction zone, potential accretion of HT soles further downdip can be predicted for a few more My after subduction initiation (Fig. 5b). Such rocks, however, are not accreted below ophiolites and/or not exhumed, which could be due to rock densities exceeding mantle values (thermodynamic modelling shows that this will be the case for basalt at 800°C for P > 1.2 GPa), to dynamics associated with melting at depth [START_REF] Faccenda | Water in the slab: A trilogy[END_REF] and/or to the resumption of full coupling between the plates [START_REF] Syracuse | The global range of subduction zone thermal models[END_REF] dragging down these rocks irreversibly.

In contrast, later and further downdip, the mantle wedge will get colder, more hydrated and therefore more buoyant while transiently remaining fairly strong (i.e., still mechanically strongly coupled). HT eclogites and the heterogeneously hydrated/weakened mantle wedge may thus reach broadly equivalent viscosities towards the end of subduction infancy, thereby favouring mechanical coupling, rock mixing and fast (buoyant) joint exhumation. This could explain the anticlockwise, short-lived exhumation of HT eclogites in serpentinite mélanges at depths of 50-60 km and ~5 My after subduction nucleation (Figs. 3a, 5c; e.g., Franciscan complex or Serpentinite mélange of Cuba:

Garcia [START_REF] García-Casco | High pressure metamorphism of ophiolites in Cuba[END_REF].

Exhumation of HT soles and basal peridotites: the depth conundrum

This study highlights the contrast between the juxtaposition of HT soles and basal peridotites at ~25-35 km (assuming purely lithostatic P estimates) and the final ophiolite thickness (≤10-15 km).

Whether this can be explained by mantle thinning [START_REF] Casey | Initiation of subduction zones along transform and accreting plate boundaries, triple-junction evolution, and forearc spreading centres: implications for ophiolitic geology and obduction[END_REF][START_REF] Dewey | The sole of an ophiolite: the Ordovician Bay of Islands Complex[END_REF] or relative exhumation (or both) has been a matter of speculation (e.g., the "conundrum of Samail"; Hacker and Gnos, 1997). Mantle thinning would have to be concentrated within the basal peridotites (500-1000 m thick at present) as the rest of the mantle section is mostly undeformed [START_REF] Ceuleneer | Mantle flow patterns at an oceanic spreading centre: the Oman peridotites record[END_REF][START_REF] Nicolas | Accretion of Oman and United Arab Emirates ophiolite -Discussion of a new structural map[END_REF], and does not explain how metamorphic soles get accreted.

Based on the pressure difference between the HT and LT soles and on the similar P-T conditions retrieved from HT soles and basal peridotites (e.g., Jamieson, 1981;McCaig, 1983; this study), the reconstruction of figure 5 depicts their relative exhumation with respect to the rest of the overlying oceanic lithosphere, thanks to buoyancy/rheology contrasts (and through successive stacks: HTsole\peridotite, then LTsole\HTsole\peridotite). The lower density of the slices would account for their exhumation with respect to the overlying mantle, while the rheology (viscosity) constrasts above and below the slices would favour strain localization. Relative exhumation is also supported by the existence of rare blueschist facies overprints on HT soles underlain by blueschist facies rocks (e.g., Turkey; Plunder et al., 2015[START_REF] Plunder | Metamorphic sole formation, emplacement and blueschist-facies overprint: early subduction dynamics witnessed by western Turkey ophiolites[END_REF], indicating that, in some cases, the HT sole is not immediately exhumed and stagnates at depth (which would not be the case if the mantle was systematically thinned).

An alternative explanation to the depth conundrum could be that the ~1 GPa pressure estimate for the HT soles corresponds to overpressure arising from strongly coupled lithologies (i.e., dry mantle against basalt; McCaig, 1983). Depths attributed to the deformation of HT soles and basal peridotites could then be reduced by up to a factor of 2 [START_REF] Petrini | Lithospheric pressure-depth relationship in compressive regions of thickened crust[END_REF], matching both the depths/pressures of LT soles and final ophiolite thickness (~10-15 km). Whether such rocks, affected by strong ductile deformation, may sustain excess dynamic pressure is unclear and could be tested with numerical models (with limitations discussed in section 5.1.2). Overpressure would not, however, affect the conclusions of this study regarding mechanical coupling, accretion or the existence of rheological switches.

Are metamorphic soles formed beneath supra-subduction ophiolites?

While the proposed mechanism emphasizes the importance of rheology during subduction infancy, the nature and genesis of ophiolites is still a matter of debate (see [START_REF] Rioux | Tectonic development of the Samail ophiolite : High-precision U-Pb zircon geochronology and Sm-Nd isotopic constraints on crustal growth and emplacement[END_REF] for a recent discussion), with authors in favour of a MORB-type pre-existing lithosphere [START_REF] Nicolas | Structures in Ophiolites and Dynamics of Oceanic Lithosphere (Petrology and Structural Geology)[END_REF][START_REF] Nicolas | Accretion of Oman and United Arab Emirates ophiolite -Discussion of a new structural map[END_REF] and others supporting an entirely supra-subduction origin [START_REF] Stern | Subduction zone infancy: examples from the Eocene Izu-Bonin-Mariana and Jurassic California arcs[END_REF]. In the first hypothesis, formation of oceanic lithosphere shortly predates intra-oceanic subduction (via near-ridge, detachment or transform fault inversion; [START_REF] Boudier | Shear zones, thrusts and related magmatism in the Oman ophiolite : Initiation of thrusting on an oceanic ridge[END_REF], while in the second hypothesis intra-oceanic spontaneous subduction (often assumed to take place at transform faults) triggers the formation of supra-subduction lithosphere by mantle upwelling (e.g., [START_REF] Stern | Subduction zone infancy: examples from the Eocene Izu-Bonin-Mariana and Jurassic California arcs[END_REF]. The interplay between lithology/rheology, T (and P) and mechanical coupling of our proposed mechanism could in principle operate through either inversion of small oceanic basins or mantle upwelling following spontaneous subduction initiation. The second scenario, however, less easily explains the systematic juxtaposition of HT soles onto LT soles (and the sharp lithological divide between them), as well as the slightly younger and more dispersed ages of the metamorphic soles when compared to those of the ophiolite crust (Fig. 1d; Wakabayashi andDilek, 2000, 2003). Metamorphic sole formation may thus question the popular spontaneous subduction initiation model [START_REF] Stern | Subduction zone infancy: examples from the Eocene Izu-Bonin-Mariana and Jurassic California arcs[END_REF].

Rheological implications

An important inference from this model is that the effective rheologies of the mantle and crust are similar at ~800°C and ~1 GPa in the presence of fluids. This is an important anchor point for experimental rock mechanics, supporting the validity of extrapolations from laboratory experiments performed orders of magnitude faster than in nature. Although the viscosity-temperature path of the mantle wedge (from i to f; Fig. 4b) or location of the first rheological switch are not determined

with precision yet, their location (at ± 50°C and ~1 order of magnitude in viscosity) further supports the convergence of the experimental grain size dependent flow laws at ~800°C (i.e., at the temperature of HT sole formation; Fig. 4d).

Figure 5 also highlights the importance, for strain localization during subduction infancy, of lubrication by hydrous phases [START_REF] Regenauer-Lieb | The Initiation of Subduction: Criticality by Addition of Water?[END_REF]Dymkova and Gerya, 2013) and grain size reduction (Linckens et al., 2011b;Hirth and Kohldstedt, 2015). Viscosity estimates for sole formation (Fig. 4d) match those inferred for the nucleation of throughgoing low viscosity shear zones (~10 20 Pa.s; [START_REF] Regenauer-Lieb | The Initiation of Subduction: Criticality by Addition of Water?[END_REF]. Although mantle wedge rheology likely depends on the extent of serpentinization and higher T ductile deformation, the similarity of ophiolite soles worldwide suggests that regional variations (e.g., incoming material, convergence rates, effective fluid release) are averaged out.

The presence of a warm oceanic lithosphere (i.e., <~3 My and/or rejuvenated; [START_REF] Hacker | Simulation of the Metamorphic and Deformational History of the Metamorphic Sole of the Oman Ophiolite[END_REF][START_REF] Duretz | Thermo-mechanical modeling of the obduction process based on the Oman Ophiolite case[END_REF] supplementary material) appears to be an important requirement for metamorphic sole formation. Their formation is probably inhibited for older lithospheres, not because they are too cold (warm conditions will be met deeper down) but because the plate interface will be too decoupled. On the other hand, previous studies suggested that one-sided intraoceanic subduction initiation requires significant mantle cooling and subsequent weakening [START_REF] Crameri | A free plate surface and weak oceanic crust produce single-sided subduction on Earth[END_REF]. Too strong mechanical coupling, during subduction infancy, would favour double-sided subduction as in the Archean (Sizova et al., 2010 and references therein). This potentially explains why ophiolite soles older than the Neo-proterozoic are missing.

These findings also bear important implications for mechanisms of sediment or seamount underplating (i.e., transient episodes of strong coupling may control, in long-lived subduction systems, the potential accretion/exhumation of HP-LT rocks, including eclogites; [START_REF] Kimura | Peeling oceanic crust in subduction zones[END_REF][START_REF] Agard | Plate acceleration: The obduction trigger? Earth Planet[END_REF], variably hydrated mantle wedge rheologies [START_REF] Faccenda | Water in the slab: A trilogy[END_REF] and for geochemical budgets of volatile fluxing during (hot) subduction [START_REF] Ishikawa | Trace element characteristics of the fluid liberated from amphibolite-facies slab: Inference from the metamorphic sole beneath the Oman ophiolite and implication for boninite genesis[END_REF][START_REF] Sizova | Subduction styles in the Precambrian: Insight from numerical experiments[END_REF].

Conclusions

Combining structure-lithology-P-T-time data from metamorphic soles with flow laws derived from experimental rock mechanics, we herein (i) outline the existence of two major, systematic rheological switches across the subduction interface (mantle wedge vs. basalts, then mantle wedge vs. sediments), (ii) propose that they control slab penetration and the successive formation of HT then LT metamorphic soles and (iii) provide a tentative reconstruction of metamorphic sole accretion during subduction infancy.

This reconstruction provides a generic explanation for the ubiquitous formation of metamorphic soles and emphasizes how slab progression is hindered, during subduction infancy, by progressive changes in the mechanical properties of the cooling plate interface, until the interface becomes fully decoupled. Metamorphic sole formation and accretion would not so much result from a transient HT event (i.e., an 'ironing' effect), but from the existence of transiently similar rheologies and strong coupling during subduction infancy.

This study sheds light on early slab dynamics and on the role of transient mechanical coupling along the plate interface (i.e., such transient episodes may also control accretion/exhumation of HP-LT rocks in mature subduction systems) and provides a testable hypothesis for thermo-mechanical models. The proposed mechanism also strengthens the applicability of experimentally-derived flow laws for the mantle and mafic crust: although calibrated at much higher strain rates, they appear to be in remarkable agreement with natural data at ~800°C (and ~1 GPa). and of the mantle immediately above the slab (Fig. 4c) when the slab descends during subduction infancy. The rheological contrast across the plate interface between the basaltic crust and mantle wedge evolves from  to ' and ". Evolution along the black arrows shows the existence of two rheological switches during subduction infancy, when slab crustal rocks reach viscosities comparable to that of the mantle wedge (dots 1a and 1b; see also Fig. 3a). The temperature range for the formation of HT and LT soles is indicated. place first order constraints on the thermal regime of our tentative reconstruction of slab penetration into the mantle and metamorphic sole formation (Fig. 5). The advantage of such a 2D thermokinematic model, as described in detail in Duprat-Oualid et al. (Fig. S2;2013), is that it includes a minimum of parameters (i.e., initial thermal ages of the oceanic lithosphere, convergence velocities and slab dips; Fig. S3), in contrast with thermo-mechanical modelling, which is beyond the scope of this study. It also provides a first grip on the depth and duration of the rheological switches outlined in our study (see Fig. S4 below).

Model design and numerical code description

The morphology of the model is shown in figure S2a. Incoming lower plate material, with a characteristic initial age for the oceanic lithosphere (T age ), is continuously buried at a constant convergence velocity (V) along a subduction plane with prescribed dip angle (θ), to simulate the subduction of an initially flat oceanic lithosphere down to the mantle. The model box is 85 km high, and the width of the model is adjusted depending onθto ensure a constant grid resolution of 1 km in all models.

The initial temperature is computed following a half space oceanic geotherm equation (e.g., [START_REF] Turcotte | [END_REF], expressing the temperature T as a function of depth (z). Such a thermal profile depends on the thermal properties of the rocks (i.e., k,ρ and Cp), on the thermal age T age of the oceanic lithosphere and on the difference between the temperature of the mantle T m and the surface T s , fixed at 1350°C and 0°C, respectively. Thermal properties used both for this initial thermal profile and during the simulations are: k = 2 W.m -1 .K -1 , ρ = 3000 kg.m -3 and Cp = 1000 J.kg -1 .K -1 .

At each time step, the code solves the following heat thermal equation: , (Eq. 1)

where ρ, Cp, and k are the density, the heat capacity and the thermal conductivity, respectively. Q, which corresponds to the heat production part (due, for instance, to radiogenic heat production or shear heating), is set to 0. This point is however discussed below and in the discussion section.

This equation (Eq. 1) is solved by using the implicit finite difference method on the Eulerian grid (e.g., Gerya, 2010). Temperature at the top of the model is fixed at 0°C, while the left and right sides of the model are subject to insulating boundary conditions (i.e., ). The boundary condition at the bottom of the model is that of a constant flux (i.e., ). The temperature T computed on the nodes is then interpolated on markers and advected through the model by following the imposed velocity field (see Gerya, 2010 and[START_REF] Duprat-Oualid | Major role of shear heating in intracontinental inverted metamorphism: Inference from a thermo-kinematic parametric study[END_REF].

Parametric study

The influence of the convergence velocity V, subduction dip angle θ and initial age of the oceanic lithosphere T age (controlling the initial thermal profile, Eq. 2) has been tested through a parametric study. For the sake of clarity, we focus the following description on the range of temperatures associated with HT sole formation (for which refined constraints on their P-T conditions of formation also exist; Fig. 3b).

Results show that the kinematic parameters (V,θ) control the time at which the tip of the slab reaches a given temperature (Fig. S2c), while the initial thermal age of the oceanic lithosphere controls the depth where this occurs (Fig. S2b). Considering a velocity of 3 cm.yr -1 and a dip of 30 °for the subducting slab (V.sin(θ) is 1.5 cm.yr -1 ), a very low initial T age is required to reach 700°C in less than ~2 My (Fig. S2b). This is why our reference experiment uses a T age of 3 Ma. In such a configuration, our models predict that the temperature of ~750°C will be reached at around 20 km depths (Fig. S2d).

Results from modelling

The initial thermal age, convergence velocity and slab dip are set to 3 Ma, 3 cm.yr -1 and 30° in the reference model, so that the tip of the oceanic crust located on top of the slab reaches 750°C and 20 km depth after ~1.5 My (and 1000°C at ~25 km depth after 1.75 My; Fig. S2d). This configuration is the one adopted in Figure 5.

Thermo-kinematic models constrain the evolution of temperature through time along the plate interface for given depths, shown here at 10-20-30-40 km for the reference experiment (Fig. S3a).

T-time trajectories of evenly spaced incoming rocks (i.e., at the slab tip and then every 30 km away, for points A, B, C and D; Fig. S1d) are shown on figure S3b.

The evolution of temperature through time along the plate interface and T-time trajectories (Fig. S3) show that a rock starting to subduct 1 My after subduction initiation (point B; Fig. S4) will be in the temperature range of formation of the HT sole and mafic/ultramafic rheological switch (~750-850°C) between 2.8 and 3.2 My at 27-32.5 km depth. The depth and duration over which incoming rocks remain in this temperature interval, for a range of realistic upper plate thermal ages and slab dips (Fig. S4), is 25±7 km (across a few km on either side of the thrust) for no more than 0.3-0.4 My.

4. Influence of shear heating. Shear heating (Hs) corresponds to the transformation of mechanical energy into heat during deformation and depends on the strain rate and on the viscosity. For a shear zone with thickness h, exhibiting a homogeneously distributed shear velocity V and an effective mean viscosity η, Hs can be expressed as (e.g., [START_REF] Duprat-Oualid | A dimensional analysis to quantify the thermal budget around lithospheric-scale shear zones[END_REF]: , In our simulations, we can thus compute this heat production for a given viscosity and a given strain rate as presented in figure S5. Results show that Hs strongly affects the thermal evolution around shear zones, especially for high viscosity values. The temperature conditions for the formation of the metamorphic sole are in this case located at shallower depth than in experiments which do not include shear heating. This results are obtained for constant viscosity and do not account for the T-dependency of viscosity: in nature, however, almost all materials (see Fig. 4b) present a decrease in viscosity with increasing temperature. These simulations therefore provide an estimate of the maximum influence of shear heating, reality being probably between situations shown in figures S2d and S5.

The important result is that, even if shear heating highly modifies the thermal field evolution around the shear zone, it does not alter the pattern of the isotherms. Figure S5 indeed shows that whatever the viscosity of the shear zone (and therefore the amount of heat produced by shear heating), the global "S-shape" of the isotherms, and therefore the existence of rheological switches, is always reproduced. Consequently, although shear heating can impact the depth of metamorphic sole formation, it does not affect the processes related to plate interface rheological switches that we propose in this study. Figure S5. Influence of shear heating on the thermal evolution of the subduction zones. Plots are provided here after 2 My (see details in the supplementary material section 4 above). The top-left panel corresponds to the thermal field obtained for the reference model, with no shear heating. In the 3 others panels, heat production by shear heating is implemented assuming a shear zone thickness of 1 km and a constant viscosity (10 20 , 10 20.5 and 10 21 Pa.s, respectively). As the viscosity of rocks is strongly temperature dependent (Fig. 4b), these simulations are not fully consistent but nevertheless useful since the viscosity values considered here lie in the range inferred from natural data (10 20 Pa.s: Tasaka et al., 2013; 10 21 Pa.s: Linckens et al., 2011). Black curves correspond to the 700 and 800°C isotherms. The black dot is located at a depth of 30 km for reference.

Figure 5

 5 Figure 5 highlights how cooling dramatically impacts slab penetration during subduction

  Figure1(a) Geodynamic setting of metamorphic sole formation during subduction infancy, following intraoceanic subduction initiation. The later geodynamic evolution will lead to continental subduction and obduction s.s. (i.e., emplacement of the oceanic lithosphere onto continental lithosphere; after[START_REF] Agard | Plate acceleration: The obduction trigger? Earth Planet[END_REF]; (b) Close-up view of Fig.1a: the formation and accretion of metamorphic soles imply a shift of the subduction interface during subduction initiation (from thrust 1 to thrust 2); (c) simplified geological map highlighting the striking continuity of the metamorphic sole beneath the mantle of the Oman ophiolite (modified after[START_REF] Nicolas | Accretion of Oman and United Arab Emirates ophiolite -Discussion of a new structural map[END_REF]; (d) age constraints for metamorphic sole formation along the Oman ophiolite (Hb: hornblende; Ms: white mica; Zr: zircon; see Table1for references). Radiometric ages for the ophiolite are shown for comparison (after[START_REF] Rioux | Tectonic development of the Samail ophiolite : High-precision U-Pb zircon geochronology and Sm-Nd isotopic constraints on crustal growth and emplacement[END_REF]; see discussion in section 5.2).

  Figure3(a) Compilation of pressure-temperature data for HT (black) and LT (white) metamorphic soles (see Table1; all are enclosed in the grey shaded area). P-T conditions for representative eclogites (HT: purple; LT: white) are also given for comparison (after[START_REF] Agard | Exhumation of oceanic blueschists amd eclogites in subduction zones: Timing and mechanisms[END_REF]. Arrows outline counterclockwise P-T paths for HT eclogites. Boxes in the HT sole domain correspond to phase diagram calculations in this study. Plain boxes: results for two sample locations in Turkey (Kütahya and Yesilova); dashed boxes: results for two locations in Oman (Sumeini and Khubakhib). The diamond-shaped box with a thicker contour correponds to the estimate for the HT sole sample from Yesilova described in the text and in figure 3b. Samples have the following molar proportions (%): Kutahya: Si (46.92), Al (16.14), Fe (9.98), Mg (10.08), Ca (10.49), Na (6.39) / Yesilova: Si (44.29), Al (16.77), Fe (7.68), Mn (0.45), Mg (15.52), Ca (10.37), Na (4.14)/ Sumeini: Si(43.24) Al(15.66) Fe(10.46) Mg(11.37) Ca(16.81) Na(2.46) / Khubakhib: Si(44.85) Al (15.67) Fe(8.82) Mg(11.51) Ca (14.36) Na (4.79). (b) Phase diagram calculated in the Na 2 O-CaO-MnO-FeO-MgO-Al 2 O 3 -SiO 2 -H 2 O chemical system for a MORB-type HT metamorphic sole (YE1302b; Lycian ophiolite, W. Turkey). The studied equilibrium peak assemblage (i.e. garnet-clinopyroxene-amphiboleplagioclase) is stable within a P-T field refined using mineral isopleths (e.g., pyrope content in garnet: X Prp ; Calcic component in clinopyroxene: Ca Cpx ) to 1.08 ± 0.1 GPa and 780 ± 40 °C. Abbreviations for minerals with solid solutions: Amph: amphibole; Cpx: clinopyroxene; Ep: epidote; Grt: garnet; Opx: orthopyroxene; Plg: plagioclase; Qtz: quartz.

Figure 4

 4 Figure4(a) Calculated temperature dependency of effective viscosities () for key plate interface lithologies (peridotite, basalt, sediments and serpentine; see Table2 and section 3.1 for details) for an average strain rate of 10 -13 s -1 (after the Oman example;Linckens et al., 2011b). Note the systematic location of the dry and wet olivine dislocation creep flow laws above those for mafic rocks (i.e., basalt, blue overlay). Abbreviations: disGBS: dislocation-accommodated grain boundary sliding; diff creep: diffusion creep; FJ07:Faul and Jackson, 2007; HK03: Hirth and Kohlstedt, 2003. Diffusion creep flow laws have been calculated for two different olivine aggregate grain sizes (10 and 50 µm) corresponding to the extreme values estimated byLinckens et al. (2011b) in Oman mantle ultramylonites. (b) Plot of averages of published flow laws (after Fig.4a) showing the qualitative viscosity decrease of a cooling and hydrating mantle wedge, from i to f (see text). Note that the viscosity of incoming crust will tend to decrease (from basalt to sediment) with increasing amounts of sediments and/or with addition of water released from progressive sediment dehydration. Stars mark the inception of melting of basalt and sediments(after Kessel et al., 2005 and White et al., 2007, respectively). (c) 2D sketch of slab penetration during subduction infancy, emphasizing the importance of fluid liberation from the slab. The dashed green line bounds the stability field of serpentine ("serp. front"; blue dots suggest that only part of this domain may be effectively hydrated). Fluid storage capacity of the mantle wedge decreases downwards (see section 3.2). Background isotherms are taken after 1 My from thermo-kinematic modelling (see supplementary material). Abbreviations: s: sediments (shown here in red, as in Fig.1b), b: basaltic layer, g: gabbroic layer, m: mantle, serp.: serpentinite; (d) Black arrows outline evolving rheologies of the crust on top of the slab (basalt, then sediment) and of the mantle immediately above the slab (Fig.4c) when the slab descends during subduction infancy. The rheological contrast across the plate interface between the basaltic crust and mantle wedge evolves from  to ' and ". Evolution along the black arrows shows the existence of two rheological switches during subduction infancy, when slab crustal rocks reach viscosities comparable to that of the mantle wedge (dots 1a and 1b; see also Fig.3a). The temperature range for the formation of HT and LT soles is indicated.

Figure S2 .

 S2 Figure S2. Thermo kinematic modelling of the temperature field associated with slab penetration into the mantle during subduction infancy; (a) Set-up of the numerical model. The initial temperature field is flat and depends on the thermal age of the oceanic lithosphere (T age ). V corresponds to the convergence velocity; (b) and (c): Variations, as a function of time and initial oceanic T age , of the depth (left panel) and burial velocity (V.sinθ; right panel) required for the tip of the slab to reach the temperature of 700°C. Depths of metamorphic sole formation are chiefly controlled by the initial T age of the oceanic lithosphere (left panel), while kinematic parameters control the time at which this happens (right panel); (d) Thermal structure of the reference model after 1.75, 2.75 and 3.75 Myrs, which was used in figure 5.

Figure S3 .

 S3 Figure S3. Thermal modelling of slab penetration during subduction infancy. (a) Evolution of slab temperatures along crust-mantle boundaries at depths of 10-20-30-40 km (thin curves). Red line: Ttime trajectory of the surface temperature of point A (tip of the slab). Orange overlay: T range of HT sole formation. Incipient melting of basaltic crust is given for an average value between wet and dry MORB (Kessel et al., 2005). See text for deformation mechanisms (disGBS: dislocationaccomodated grain boundary sliding; Linckens et al., 2011a; Hirth and Kohlstedt, 2015); (b) Thick red lines give T-time trajectories for evenly spaced points along the top of the slab (A: slab tip; points A to E are spaced every 30 km; see Fig. S2d). Temperatures along the slab Moho are shown with dashed blue lines. Open black circles indicate when a given point enters the conditions of metamorphic sole formation (black arrows outline approximate duration under these conditions).

Figure S4 .

 S4 Figure S4. Duration and depth over which incoming rocks remain between 750 and 850°C, during the first rheological switch and HT sole formation (Fig. 4d), depending on upper plate thermal age and slab dip. For the reference experiment, incoming crustal rocks may cross the temperature range of HT sole formation (~750-850°C) at depths of ~25-35 km: accretion of individual slices may last on the order of 0.3-0.4 My, while ~1-3 My would be required to fully accrete the first incoming 100 km of crust on the top of the slab (A to C or D; see text, section 5.2.i for further discussion).

Table 2 843

 2 Published

	Lithology	Authors	n	A	Q	m
					(MPa -n )*(s -1 )	(kJ/mol)
	Slab crust
	Sediments
	Quartz		Ranalli and Murphy, 1987	3	6.80E-06
	Wet quartzite	Kirby and Kronenberg,1987	2.3	3.20E-04
	Biotite schist	Kronenberg et al., 1990	18	1.20E-30
	Schists		Shea and Kronenberg, 1993	31	1.30E-67
	Felsic granulite	Wilks and Carter, 1990	3.1	8.00E-03
	Mafic crust
	Basalt		Shelton and Tullis, 1981
			and Hacker and Christie, 1990	3.5	1.00E-04
	Wet diabase	Shelton and Tullis, 1981	3.4	2.00E+04
	Dry diabase	Mackwell et al., 1998	4.7	8.00E+00
	Diabase	Van Hunen and Van den Berg, 2008	3.4	2.21E-04
	Microgabbro	Wilks and Carter, 1990	3.5	4.85E+04
	Mafic granulite	Wilks and Carter, 1990	4.2	1.40E+04
	Mantle wedge
	Dry olivine	Chopra and Paterson, 1981	3.6	4.50E+00
		"	Chopra and Paterson, 1984	3	1.00E+04
		"	Chopra and Paterson, 1984	3.5	2.50E-04
		"	Karato et al., 1986	3.5	5.40E+00
		"	Karato and Wu, 1993	3.5	2.42E+05
		"	Karato and Rubie, 1997	3	2.40E+05
		"	Bussod et al., 1993	3.5	1.12E+05
		"	Hirth and Kohlstedt, 2003	3.5	1.10E+05
		"	Karato and Jung, 2003	3	1.26E+06
		"	Li et al., 2006	3	4.57E+03
		"	Kawazoe et al., 2009	3.5	1.10E+05
	Diffusion creep	Hirth and Kohlstedt, 2003	1	1.5E+09		3
	Diffusion creep	Faul and Jackson, 2007	1.4	2E+10		2
	DisGBS	Hirth and Kohlstedt, 2003	3.5	6.5E+03		2
	Wet olivine	Chopra and Paterson, 1981	4.4	2.76E+02
		"	Chopra and Paterson, 1981	4	2.00E+03
		"	Karato et al., 1986	3	1.50E+06
		"	Evans and Kohlstedt, 1995	4.5	2.60E+00
		"	Hirth and Kohlstedt, 1996	3.5	4.88E+06
		"	Mei and Kohlstedt, 2000	3	4.57E+03
		"	Mei and Kohlstedt, 2000	3	5.01E+02
		"	Karato and Jung, 2003	3	3.63E+00
		"	Karato and Jung, 2003	3	7.94E+02
		"	Hirth and Kohlstedt, 2003	3.5	9.00E+01
		"	Hirth and Kohlstedt, 2003	3.5	1.60E+03
		"	Mc Donnell et al., 1999*	2.14	9.10E+03		3
	Serpentinite	Hilairet et al., 2007	5.8	2.51E-13	20.8
	842	"	Hilairet et al., 2007	3.8	2.51E-09	12.1

flow laws and material creep parameters used to compute the curves of figure 4b. n, A, Q 844 and m correspond to the stress exponent, pre-exponential factor activation energy and grain size 845 exponent, respectively (*calculated with X H2O =0.05 wt%). DisGBS: dislocation-accommodated 846 grain boundary sliding. See supplementary material for references.

Acknowledgements

This work was funded through the Agence nationale de la recherche (ANR-10-BLAN-0615) and through Institut Universitaire de France (grant to P.A.). We are grateful to colleagues from the O:NLAP ANR project for fruitful and enthusiastic discussions. An earlier version of this manuscript benefited from the constructive remarks of our friend and colleague E. Burov, who sadly passed away, and A. Schubnel. This work also benefited from discussions with scientists from the ZIP project (REA grant agreement no. 604713 (ZIP "Zooming In between Plates") from the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme FP7/2007-2013).

APPENDIX. SUPPLEMENTARY MATERIAL A. Effective viscosities of plate interface material constrained by rock mechanics

Figure S1. Influence of strain rate on the effective viscosities of selected lithologies at two different temperatures (i.e., 500 and 650°C; same colors as in figure 4a, except for felsic granulite). Note that the respective positions of the dry mantle, basalt, quartz and serpentine remain unchanged for strain rates > ~10 14 s -1 . Note the stronger dependence on strain rate of serpentine and biotite schist, which are almost independent of temperature (as seen in Fig. 4a). Abbreviations: H07: Hilairet et al., 2007;[START_REF] Karato | Rheology of the upper mantle: a synthesis[END_REF]K90: Kronenberg et al., 1990;RM87: Ranalli and Murphy, 1987;ST81: Shelton and Tullis, 1981;WK90: Wilks and Carter, 1990.


Davide Faranda

Bérengère Dubrulle

Maria Emanuele Flavio

Pons

Statistical early-warning indicators based on Auto-Regressive Moving-Average processes

We address the problem of defining early warning indicators of critical transition. To this purpose, we fit the relevant time series through a class of linear models, known as Auto-Regressive Moving-Average (ARMA(p, q)) models. We define two indicators representing the total order and the total persistence of the process, linked, respectively, to the shape and to the characteristic decay time of the autocorrelation function of the process. We successfully test the method to detect transitions in a Langevin model and a 2D Ising model with nearest-neighbour interaction. We then apply the method to complex systems, namely for dynamo thresholds and financial crisis detection.

Many experimental or natural systems undergo critical transitions -sudden shifts from one to another dynamical regime. In some instances, e.g. global changes in climate science, species extinction in ecology, spin glasses, it is of crucial importance to build early warning indicators, i.e. estimates of the transition threshold based on finite time-series corresponding to situations where the bifurcation did not happen yet. The statistical approach to this issue traditionally involves so-called indicators of criticality [START_REF] Scheffer | [END_REF]2]. They are based on specific properties of ideal statistical systems (such as the Langevin or Ising model) near the transition: critical slowing down, modifications of the auto-correlation function or of the fluctuations [3], increase of variance and skewness [4], diverging susceptibility [5][6][7], diverging correlation length. However, it is known that, in some complex systems, these indicators fail to detect the transition: in spin glasses, no diverging correlation length has been found so far, and one has to resort to finer statistical tools (such as four point dynamical susceptibility [8,9]) to detect transitions. In addition, traditional early warning indicators may be inapplicable in datasets containing a small number of observations (see e.g. [10,11]), which is usually the case in several applications where the experiment is unique (as in financial or climate time series), difficult to repeat or to sample for a long time (as for atmospheric measurements, laboratory turbulence, etc). This suggests that indicators based on single statistical properties may not be sufficient for detecting transitions in complex systems, so that one should rather consider all the information contained in the finite-time series.

The main idea of the present letter is therefore to introduce a class of indicators of critical transitions based on a statistical model for the observed data when approaching a tipping point. To be interesting for applications, these indicators must satisfy certain properties: i) they must generalize the well-know indicators based on single statistical properties and ii) they must be built using a statistical model that is simple to implement and works for limited data set. Here, we show that these properties are satisfied for indicators based on the auto-regressive moving-average processes of order p, q ARMA(p, q), modeling a time series X t (λ), experiencing a transition at λ = λ c . In the first part of the paper, we first recall some basics on ARMA(p, q) modeling and define corresponding early-warning indicators. We then check that these indicators are able to detect the transition in simple theoretical models, such as Langevin double-well model or Ising model. We then apply our indicators to the analysis of complex systems for the detection of turbulent dynamo thresholds and financial crisis.

Theory. Let us consider a series of observations X t of an observable with unknown underlying dynamics, controlled by a parameter λ. We further assume that for λ < λ c the time series X t (λ) represents a stationary phenomenon. The critical threshold λ c is defined through the condition that for λ ≥ λ c , there is a bifurcation, in the sense that there exists no smooth transformation of the physical measure through the transition. Since X t (λ) is stationary, we may then model it by an ARMA(p(λ), q(λ)) process such that for all t:

X t (λ) = p i=1 φ i (λ)X t-i + ε t + q j=1 θ j (λ)ε t-j (1)
with ε t ∼ W N (0, σ 2 ) -where W N stands for white noise -and the polynomials φ

(z) = 1 -φ 1 z t-1 -• • • -φ p z t-p and θ(z) = 1 -θ 1 z t-1 -• • • -θ q z t-q
, with z ∈ C, have no common factors. Notice that, hereinafter, the noise term ε t will be assumed to be a white noise, which is a very general condition [START_REF] Brockwell | Time series: theory and methods[END_REF]. For a general stationary time series, this model is not unique. However there are several standard procedures for selecting the model which fits at best the data. The one we exploit in this paper is the Box-Jenkis procedure [START_REF] Box | Time Series Analysis: Forecasting and Control[END_REF]. We chose the lowest p and q such that residuals of the series filtered by the process ARMA(p, q) are delta correlated. This fixes p and q, and thus our statistical model. There are other model selection procedures based on information criteria (Bayesan or Akaike information criteria). Unfortunately, in our case none of them gives clear indications for discriminating the model to be used as it is not clear which range p and q must be tested to get reliable results. Intuitively, p and q are related to memory lag of the process, while the coefficients φ i (λ) and θ i (λ) represent the persistence: the higher their sum (in absolute value), the slower the system is forgetting its past history. In the sequel, we now present early warning indicators based on these parameters.

Early-warning indicators. Far from the transition, the time series of a generic physical observable can be described by an ARMA(p, q) model with a reasonably low number of p, q parameters and coefficients. On the other hand, for λ → λ c , the critical value corresponding to a transition, the statistical properties (such as the shape and/or the persistence of the autocorrelation function) of the system change, leading to different characteristics of the ARMA(p, q) model which can describe the data series or to an inadequacy of the model itself. Specifically, several changes in the dynamics occur near the transition, as the system is allowed to explore a larger portion of the phase space with higher excursions in the direction of the new stable state. First, the distributions of the observables become skewed towards the maxima or the minima, depending on the direction of the shift. Second, the system may experience a critical slowing down with diverging memory effects. This phenomenon is traditionally quantified by the autocorrelation function (ACF) of the time series X t defined (assuming a zero-mean process) as:

ACF(h) = E[X t+h X t ]/E[X 2 t ]. (2)
Far from the transition, the ACF tends to be 0 after a finite number of lags h. As λ → λ c , critical increase of memory of the system makes h → ∞. The ARMA(p(λ), q((λ)) model of the corresponding time series will then be characterized by two properties:

• p i=1 |φ i | and 1 j=1 |θ j | increase for λ → λ c as the direct consequence of h → ∞.
• p + q increases for λ → λ c because of additional new time scales associated to the trajectories moving towards the potential barrier between the two attracting states.

This rather simple observation allows us to define two indicators: O(λ) = p(λ) + q(λ), which diverges for λ → λ c , and the total persistence of the process:

R(λ) = p i=1 |φ i (λ)| + q i=1 |θ i (λ)|
which also show a divergent behavior at the transition, unless O(λ) = 1. In this latter case R(λ) = φ 1 → 1 for λ → λ c . These indicators present several advantages with respect to the traditional ones reported, for example, in [START_REF] Scheffer | [END_REF]. First, by computing the functional form for p(λ) and q(λ) and for the coefficients φ i (λ) and θ i (λ) one has also an effective statistical toy model for describing the phenomenon and to produce data with analogous statistical properties. This may be very useful for series or data which can hardly be reproduced (laboratory experiments) or integrated by new observations (climate datasets, stock market titles). Second, if several series at different λ's are available, one can extrapolate the characteristics of the process at not yet measured λ's. This property can be very useful for devising new experiments knowing the possible location of the transitions. Third, if the transition is marked by the fact that R(λ) → O(λ) rather than by a divergence of O(λ), one may argue that the potential landscape for the observable X does not change significantly when approaching the transition and therefore a Langevin reduction to a double well system is possible. If, on the contrary, the order changes significantly approaching the transition, such a low dimensional reduction is not appropriate and one should be very careful in pursuing such a model as shown, for a relevant climatic example, in [START_REF] Lucarini | [END_REF].

A toy model for critical transitions. We start considering a classical system featuring bistability under the effect of random noise, i.e. dX = -V (X)dt + dW (3) with potential V (X) = aX 4 -bX 2 +λX, where λ ≥ 0 and W is a Wiener process with unit variance. We consider system (3) for values of λ such that, in the deterministic limit, it features two stable fixed points ( X1 < 0 and X2 > 0) and an unstable fixed point X. The asymptotic behavior of the system can be assessed in terms of the solution of a Fokker-Plank equation [START_REF] Risken | Fokker-Planck Equation[END_REF]. Here we are rather interested in the finite-time behavior and we consider only the simulations such that the noise does not push the system across the bifurcation, i.e. the system is confined in one of the two wells. For this system we perform the following numerical experiments: for each value of λ we compute an ensemble of 500 trajectories X(λ) finding, for each of them, the best ARMA(p(λ), q(λ)) in the sense specified by the Box-Jenkins procedure [START_REF] Box | Time Series Analysis: Forecasting and Control[END_REF].

Then, O(λ) and R(λ) have been computed, here • stands for the ensemble average. In Fig. 1-a we report the results of this analysis, which clearly show that the average order is not affected in this case and O(λ) 1, whereas the transition is well highlighted by R(λ) which approaches O(λ) = 1 for λ → λ c . There is a simple way to understand this behavior by linking the orders p, q and of the coefficients φ i to the autocorrelation function ACF (see [START_REF] Brockwell | Time series: theory and methods[END_REF]-Chapter 3 for more details). The orders are linked to the functional form of the ACF whereas the values of φ i depend on the decay rate. In the case of system given by Eq. 3, the shape of the ACF is exponential both far from the transition (Fig 1-b) and when approaching it (Fig 1-c). However, in the latter case, the decay is much slower, this causing the increase R(λ) → 1.

The Ising model. As a second test of the indicator, we consider a classical 2-D Ising dynamics with a nearestneighbor interaction on a square lattice of size L. At each site j, a discrete spin is allowed to have two values σ j ∈ {+1, -1}. The energy of the configuration is given by the Hamiltonian:

H = -J neighbors σ i σ j (4)
under the interaction J. We consider only the case J = 1 and evolve the system by using Metropolis algorithm [START_REF] Haario | [END_REF].

A second order phase transition is expected at the temperature T = T c = 2/ ln(1 + √ 2) 2.269. To apply our early warning indicators, we performed 100 simulations for a square lattice of size L = 256 at different T > T c . We checked that our results do not depend sensitively on the size of the lattice, provided that L > 128. After discarding 100000 time iterations necessary to reach a clearly identifiable stationary state, for each temperature, an ensemble of 15 time series consisting of 200000 time units of H(t, T ) is analyzed with the procedure described above. Stationarity has been tested performing a Dickey-Fuller test on each time series. The results for O and R are reported in Fig. 2. It is evident that O and R increase when moving towards the critical temperature T c . In this case, not only the persistence of the correlations R, but also the number of terms O necessary to describe the process increases. This means that the transition cannot be modeled by a simple Langevin equation as other time scales become important. In other words, this transition is associated to a non-trivial unknown potential landscape.

An example of complex system. Up to now we have analyzed toy systems, extensively studied both analytically and numerically and for which the threshold are analitically predictable. However, interesting systems, such as turbulence or finance, lie on another level of complexity and one naturally wonders whether the technique provides reliable results. We focus on the data of the Von Kármán Sodium (VKS) experiment, a successful attempt to get a transition to dynamo in a laboratory turbulent liquid-metal experiment. The control parameter for the transition is the magnetic Reynolds number Rm. The interesting characteristic of this experiment is that several dynamo and no dynamo configurations have been obtained by changing the material of the impellers and of the cylinder [7,17]. Here we focus on two different configurations: (i) one producing a well-documented stationary dynamo at Rm ≈ 44, thereby providing a fair test of our method and (ii) one that failed to produce the dynamo within the accessible values of Rm. The time series is constructed using the modulus of the magnetic field |B|(t) as a function of time t, measured within six detectors in the core of the vessel. From this, we extract the quantities O and R, averaged over the six detectors. The results are plotted in Fig. 3: the main figure for the configuration (i), the inset for the configuration (ii). They depend quantitatively on the sensors chosen, but not qualitatively as the transition is always detected at the same Rm. The transition is very net in terms of divergence of O and R and can be located at Rm f = 47, when the dynamo is observed. In the non-dynamo case, no sign of transition is visible.

Financial crisis. We conclude by discussing the performances of the ARMA early warning indicators applied to the EUR/USD exchange rate hourly datasets (Figure 4-b). The chosen observable is the log-return of the time series, a quantity commonly examined in finance as the series obtained this way do not contain non-stationarities:

R(t) = log(X t ) -log(X t-1 ).
Here X t is the EUR/USD exchange hourly rate. Monthly values for O and R have been extracted from the time series of R(t) and results are shown in figure 4-a. Our analysis can be safely performed on these series since they are stationary, as it results from the Dickey-Fuller test [18]. The technique clearly points to three distinct warnings (marked by the red dotted lines). Interestingly, they are followed after a few months delay, by official warnings of the European Central Bank (ECB). The first warning corresponds to the Sub-prime American crisis, the second to the Greek crisis, and the third one to the Irish crisis. The crisis for the real market falls immediately after the ECB announcements. If we compare these results with the ones arising from physical systems, the warning seems to appear too early. We may argue that indicators which provide similar warnings are available also at the ECB. The time between the early warning discover and the ECB announcements may serve to the ECB for trying corrections and avoid an immediate financial crisis which is announced only when the crisis itself is unavoidable. Similar behaviors have been discovered for early warning indicators applied to financial datasets, as reported in [19,20] .

Discussion. In this work we have proposed a new method to detect early warnings of critical transitions via a statistical approach which allows also for incorporat- ing the information of several statistical indicators analysed in [2]. Here, we exploit the properties of linear (i.e. stationary and invertible) autoregressive moving-average processes, denoted ARMA(p(λ), q(λ)), being λ the system control parameter. More specifically, we have defined two indicators representing the total order and the total persistence of the stochastic process. An increase of the former is indicative of the impossibility to represent the data series in a parsimonious thus leading to the idea that the linearity hypothesis fails and the decay of the autocorrelation function of the process is no longer exponential; an increase of the latter is linked to a longer decay time of the correlation, and can be due to the increase of the total order or just of the persistence itself. The two phenomena are very different and, up to our knowledge, the traditional indicators exploit only the increase of the memory of the system (critical slowing down) to identify the threshold λ c . Here, the possibility of detecting substantial modifications to the shape of the autocorrelation function via the change in the orders p, q, allows for understanding whether reductions to simple low dimensional models are relevant or not for describing the dynamics. We have combined these two indicators to detect critical transitions both in models and in real systems. In all the cases considered, the behavior of the indicators has shown to be an effective way to investigate the proximity of the system to a critical transition; thus, they seem to be a useful tool to study critical transitions, since their estimation involves well-known, standard statistical techniques characterized by a low computational cost and applicable to relatively short time series. The application to finance seems promising. It would be interesting to extend this approach to other financial indicators as well as to climate data. On a theoretical level, one could use the technique to understand how transitions are modified when systems originally in equilibrium are driven out of equilibrium by forcing-dissipation mechanisms, starting from conceptual toy model of outof-equilibrium Ising dynamics [21,22].

FIG. 1 :

 1 FIG. 1: (a): O and R for the system defined in Eq. 3. Each error-bar represents the average of 15 realizations and the standard deviation of the mean. b) ACF for (λ-λc)/λc = -0.271. c) ACF for (λ -λc)/λc = -0.05 .

FIG. 2 :

 2 FIG.2:O and R for the system defined in Eq. 4, L=256. Each error-bar represents the average of 15 realizations and the standard deviation of the mean.

FIG. 3 :

 3 FIG. 3: Averaged O and R for the Von Karman -Sodium experiment. Solid lines refer to the experiments for increasing values of Rm. Inset: same as the main figure but for a configuration where the dynamo has not been observed.

FIG. 4 :

 4 FIG. 4: Upper panel: Average O and R for the log-return series R(t) of the EUR/USD hourly exchange rate. Central panel: EUR/USD hourly exchange rate. Lower panel: R(t). Red dotted lines correspond to early warning of the crisis. Black continuous lines correspond to actions taken by the ECB.