# Extrusion of subducted crust explains the emplacement of far-travelled ophiolites

Kristóf Porkoláb, Thibault Duretz, Philippe Yamato, Antoine Auzemery, Ernst Willingshofer

► To cite this version:

Kristóf Porkoláb, Thibault Duretz, Philippe Yamato, Antoine Auzemery, Ernst Willingshofer. Extrusion of subducted crust explains the emplacement of far-travelled ophiolites. Nature Communications, 2021, 12 (1), pp.1499. 10.1038/s41467-021-21866-1. insu-03163067

HAL Id: insu-03163067

https://insu.hal.science/insu-03163067

Submitted on 9 Mar 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ARTICLE

https://doi.org/10.1038/s41467-021-21866-1

OPEN

<figure><img src="image_3.png" /></figure>

Check for updates

# Extrusion of subducted crust explains the emplacement of far-travelled ophiolites

Kristóf Porkoláb<fnref n="1" />✉, Thibault Duretz², Philippe Yamato²,³, Antoine Auzemery<fnref n="1" /> & Ernst Willingshofer<fnref n="1" />

Continental subduction below oceanic plates and associated emplacement of ophiolite sheets remain enigmatic chapters in global plate tectonics. Numerous ophiolite belts on Earth exhibit a far-travelled ophiolite sheet that is separated from its oceanic root by tectonic windows exposing continental crust, which experienced subduction-related high pressure-low temperature metamorphism during obduction. However, the link between continental subduction-exhumation dynamics and far-travelled ophiolite emplacement remains poorly understood. Here we combine data collected from ophiolite belts worldwide with thermo-mechanical simulations of continental subduction dynamics to show the causal link between the extrusion of subducted continental crust and the emplacement of far-travelled ophiolites. Our results reveal that buoyancy-driven extrusion of subducted crust triggers necking and breaking of the overriding oceanic upper plate. The broken-off piece of oceanic lithosphere is then transported on top of the continent along a flat thrust segment and becomes a far-travelled ophiolite sheet separated from its root by the extruded continental crust. Our results indicate that the extrusion of the subducted continental crust and the emplacement of far-travelled ophiolite sheets are inseparable processes.

<footnote n="1">¹Faculty of Geoscience, Utrecht University, Utrecht, Netherlands. ²Univ Rennes, CNRS, Géosciences Rennes, Rennes, France. ³Institut Universitaire de France (IUF), Paris, France. ✉email: kristof.porkolab@gmail.com</footnote>
Finding physical mechanisms that explain how dense oceanic lithosphere, referred to as ophiolite, is emplaced on top of lighter continental plates (i.e. obduction) has prompted a long-standing scientific discussion¹⁻³. Ophiolites may be accreted to continents by being scraped off from the subducting oceanic lower plate¹,⁴, or emplaced on top of the continent in a continental lower plate – oceanic upper plate subduction setting⁵,⁶. In the latter case, subduction typically initiates within the ocean resulting in the formation of high temperature-low/medium pressure (HT-LP/MP) metamorphic soles at the base of the upper plate⁷,⁸. Intra-oceanic subduction is followed by continental subduction below the oceanic upper plate when the continental margin enters the subduction zone. Such ophiolite belts exhibit high pressure-low temperature (HP-LT) metamorphic units structurally underlying the ophiolite (Fig. 1a–c). These HP-LT units predominantly consist of

continental upper crust that represents the former passive margin of the subducting continent⁹⁻¹¹. The subduction and exhumation of these units appear to be a relatively short-lived process, as evidenced by the characteristic duration of 10–30 Myr for the subduction-exhumation cycle recorded in the metamorphic rocks (Fig. 1e). The exhumed continental units are typically found between an ophiolite sheet that was transported for tens or even hundreds of kilometers over the continent (far-travelled ophiolite sheets) and the rest of the oceanic domain (open ocean or suture zone of the former ocean; Fig. 1a, b). Far-travelled ophiolite sheets are up to 10 km thick, and their average width is 50–55 km (Fig. 1e), suggesting that their dimensions are mechanically limited during emplacement.

Models of ophiolite emplacement hence have to account for (1) short-lived continental subduction below the oceanic upper plate followed by exhumation of the subducted upper crust, and (2) the

<figure><img src="image_4.png" /></figure>

**Fig. 1** Datasets collected from natural ophiolite belts. For detailed explanation and references regarding the datasets (c–e) see Supplementary Note 1 and Supplementary Tables 1 and 2. Simplified cross sections showing the structure of the Oman¹³ (a) and New Caledonia⁵⁶,⁵⁷ (b) ophiolites and their structural relationship to the underlying continental rocks. The far-travelled ophiolite sheets are separated from their oceanic root by exhumed continental units, which show a transition from non-metamorphic to high metamorphic grade. c Peak pressure-temperature (P-T) data collected from the subducted and exhumed geological units in natural ophiolite belts. Multiple data points from the same location represent the peak conditions of different structural levels. Error bars represent the uncertainty of thermodynamic calculations. d Duration of the continental subduction-exhumation cycle below the oceanic upper plate as recorded by the buried and exhumed continental formations. The duration is calculated from the arrival of the continental margin in the subduction zone to the near surface/surface exhumation of the subducted formations and is indicated by colored circles. Error bars show uncertainties associated with the duration of the subduction-exhumation cycle, where solid lines reflect well-constrained and dotted lines poorly constrained cases, respectively (see Supplementary Note 1). e Average width of far-travelled ophiolite sheets (colored circles) and the deviation from the average along strike (error bars) in natural ophiolite belts. The width was measured along 3–8 cross sections along each ophiolite belt.
observed dimensions of ophiolite sheets. Recent numerical modelling studies achieved these conditions by imposing con-vergence to reach the state of continental subduction, and sub-sequently imposing divergence to exhume the continental margin and thin the oceanic upper plate⁶,¹². However, these models do not lead to nappe formation in the subducted continental crust, which is a key feature of natural systems and seems to play an important role in the exhumation of the HP-LT continental rocks⁹,¹³,¹⁴. Moreover, stacking of continental nappes below the oceanic upper plate is critical for inducing uplift on top of the rising thrust sheet(s), which in turn may lead to gravity-driven extension in the upper plate¹⁵⁻¹⁷. Extension in the upper plate might significantly contribute to the unroofing of the subducted continental rocks as well as in separating the far-travelled ophiolite sheet from its root¹⁵,¹⁷. Hence, exploring decoupling and nappe formation mechanisms in the continental lower plate is crucial for understanding the deformation of the oceanic upper plate and far-travelled ophiolite emplacement.

Here, we combine numerical thermo-mechanical simulations of oceanic upper plate-continental lower plate subduction systems and data acquired in ophiolite belts worldwide to unravel the physical processes that explain the structure of ophiolite belts. We highlight the genetic link between the extrusion of the subducted continental crust and the emplacement of far-travelled ophiolite sheets, and further identify the key parameters controlling this process.

## Results

**Modelling strategy.** We designed 2D thermo-mechanical numerical simulations governed by momentum, mass, and heat conservation equations and a visco-elasto-plastic rheological model (see Methods section for details of numerical modeling techniques and further model setup description). A total plate convergence velocity of 3 cm/yr is achieved by prescribing constant normal inflow velocities of $|V_{in}| = 1.5$ cm/yr along the upper 140 km of the two model sides. Mass conservation is satisfied by gradually increasing outflow below 140 km (Fig. 2a). A convergence velocity test was performed to explore the role of this parameter (Supplementary Note 2 and Supplementary Fig. 1). The top boundary of the model is a true free surface¹⁸. The initial model geometry is inspired by reconstructions of pre-obduction geodynamic settings where intra-oceanic subduction is initiated relatively close (<400 km) to the continental passive margin, which arrives to the subduction zone after ~10 Myr of oceanic subduction¹⁹⁻²¹. Subduction is designed to initiate along an inclined weak zone (equivalent to an oceanic detachment) close to the mid-ocean ridge to achieve a right-dipping subduction zone and a thermally young ophiolite front, in agreement with natural ophiolite belts²². A simplified continental passive margin geometry is implemented by linearly decreasing crustal thickness over the distance of 200 km (Fig. 2a). The continental basement is divided in two parts (upper and lower crust) constituted of different materials (Tables 1 and 2), which introduces a decoupling level at the base of the upper crust (Fig. 2a). A resolution test was performed to verify the robustness of the model results (Supplementary Fig. 2).

**General evolution of the reference model.** Intra-oceanic subduction initiation is followed by the subduction of oceanic lithosphere until 10 Myr. The passive margin of the continent then starts to subduct below the oceanic upper plate (i.e. obduction; Fig. 2b), and experiences HP-LT metamorphism up to eclogite facies P-T conditions (Fig. 2c). The 3-km-thick sedimentary cover of the passive margin partially subducts with the rest of the continental lithosphere, but is largely stacked and

accreted to the front of the oceanic upper plate (Fig. 2). After reaching eclogite facies conditions, the burial velocity of the upper crust decreases from ~3 cm/yr at 18 Myr to near zero at 23 Myr. This is the moment when the upper crust starts to decouple from the lower crust and the lithospheric mantle. Decoupling results in the localization of a major reverse-sense shear zone along which the subducted upper crust is extruded upwards and leftwards (Fig. 2d). The oceanic upper plate undergoes gravity-driven extension when being pushed up by the extruding thrust sheet. Extension leads to the breaking of the oceanic upper plate, which enables the continental upper crust to be rapidly extruded to the surface, separating a 50-km-wide and maximum 13-km-thick far-travelled ophiolite sheet from the rest of the oceanic lithosphere (Fig. 2d, e).

**Ductile nappe formation: the onset of exhumation.** Burial of continental crust below a denser oceanic plate leads to a progressive increase of the buoyancy force that resists subduction despite the continuously imposed plate convergence. From ~21 Myr (i.e. after 11 Myr of continental subduction) the velocity of the subducting continental upper crust becomes near zero, while the lower crust and the lithospheric mantle still subduct at ~1.5 cm/yr (Fig. 3a). This kinematic setting leads to an increase of deviatoric stress within the lower crust, which meets the conditions for strain localization by thermal softening²³. As a result, a major reverse-sense shear zone develops and exhibits strain rates of >$10^{-13}$ s⁻¹ (Fig. 3a). The shear zone is initially low-angle, and propagates along the base of the upper crust, which is the weakest horizon in the continental crust (Figs. 2a and 3a). The shear zone thus facilitates decoupling between the crustal layers, which is further accommodated by distributed deformation (folding) at a strain rate of ~$10^{-14}$ s⁻¹ in the entire subducted upper crust.

**Interplay between crustal extrusion and upper plate necking.** Upper crustal decoupling and nappe formation results in uplift, which triggers extension of the oceanic upper plate over a 50-km-wide zone (Fig. 3a). Extension leads to the formation of normal faults and to the normal-sense reactivation of the original plate boundary thrust (Fig. 3a).

Further strain localization results in the connection of the two flat reverse-sense shear zone segments by a steeper ramp segment. This structure separates very-low-grade to non-metamorphic upper crust from the low-grade to eclogite-facies upper crust (Fig. 3b). The localization of the main shear zone facilitates the acceleration of upper crustal extrusion (from 0.5–1 cm/yr until 24.5 Myr to 1.5–3 cm/yr from 25 Myr), which is accommodated by increasing displacement along the right-and left-dipping extensional shear zones and leads to the necking of the oceanic upper plate (Fig. 3b, c). Through this process a sheet of oceanic lithosphere (the future far-travelled ophiolite) gets disconnected from the oceanic plate along the left-dipping normal fault, which joins the upper flat segment of the main thrust in depth (Fig. 3b). Subsequently, the far-travelled ophiolite sheet is emplaced on top of the continent and transported to the left along the upper flat thrust segment (Fig. 3c). This process is driven by the progressive extrusion of the continental crust.

**Key parameters controlling crustal extrusion and far-travelled ophiolite emplacement.** Extensive tests have been performed on the reference model to determine key parameters that control the extrusion of the upper crustal thrust sheet, which is instrumental for the necking of the upper plate and the emplacement of far-travelled ophiolite sheets (Fig. 4). In particular, we evaluated the impact of the crustal rheology of the subducting plate and shear heating on the process of nappe formation.
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Model setup and evolution of the reference model. The compositional evolution and associated *P*-*T* paths are presented in **b-e** on the left and right side, respectively. Colored circles are the positions of selected particles for recording *P*-*T* conditions and correspond to the colored lines on the left-side *P*-*T* diagrams. **a** Model setup showing the initial distribution of compositional domains, the velocity boundary condition, and the strength profile (second stress invariant vs depth, assuming a transcurrent setting) for the continental domain. For thermal and rheological parameters of the model compositions see Table 1. Layering in the continental upper crust corresponds to passive strain markers. White lines are isotherms from 200 to 1200 °C. **b** Early-stage continental subduction at 13.9 Myr after the subduction of the left-side ocean between 0 and 10 Myr. **c** Late-stage continental subduction after the initiation of upper crustal decoupling. **d** Breaking of the oceanic upper plate and accelerated extrusion of the subducted upper crust. **e** Emplacement of the far-travelled ophiolite sheet further leftwards on top of the continent, dictated by crustal extrusion.</figcaption></figure>
<table><caption>**Table 1 Thermal and rheological parameters used for different compositions in the reference model.**</caption><thead><tr><th></th><th>ρ (kg m<sup>-3</sup>)</th><th>k (W m<sup>-1</sup> K<sup>-1</sup>)</th><th>Q<sub>r</sub> (W m<sup>-3</sup>)</th><th>α (K<sup>-1</sup>)</th><th>C (MPa)</th><th>φ (°)</th><th>A (Pa<sup>-1</sup> s<sup>-1</sup>)</th><th>n</th><th>Q (J mol<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Sedimentary cover (mica)</td><td>2700</td><td>2.55</td><td>2.9e-6</td><td>3.0e-5</td><td>10</td><td>15</td><td>1.0e-138</td><td>18</td><td>51.0e3</td></tr><tr><td>Continental upper crust (Westerly granite)</td><td>2750</td><td>2.8</td><td>1.65e-6</td><td>3.0e-4</td><td>10</td><td>30</td><td>3.1623e-26</td><td>3.3</td><td>186.5e3</td></tr><tr><td>Continental lower crust (mafic granulite)</td><td>2900</td><td>2.8</td><td>1.65e-6</td><td>3.0e-4</td><td>10</td><td>30</td><td>8.8334e-22</td><td>4.2</td><td>445.0e3</td></tr><tr><td>Oceanic crust (Maryland diabase)</td><td>2900</td><td>3.0</td><td>1.0e-10</td><td>3.0e-5</td><td>10</td><td>30</td><td>3.2e-20</td><td>3.0</td><td>276.0e3</td></tr><tr><td>Lithospheric mantle (dry olivine)</td><td>3300</td><td>3.0</td><td>1.0e-10</td><td>3.0e-5</td><td>10</td><td>30</td><td>1.1e-16</td><td>3.5</td><td>530.0e3</td></tr><tr><td>Asthenosphere (dry olivine)</td><td>3300</td><td>3.0</td><td>1.0e-10</td><td>3.0e-5</td><td>10</td><td>30</td><td>1.1e-16</td><td>3.5</td><td>530.0e3</td></tr><tr><td>Weak zone (serpentinite)</td><td>2900</td><td>3.0</td><td>1.0e-10</td><td>3.0e-5</td><td>0</td><td>30</td><td>4.4738e-38</td><td>3.8</td><td>8.9e3</td></tr></tbody></table>

The heat capacity ($C_p$) and the compressibility ($\beta$) were set to 1050 J kg$^{-1}$ K$^{-1}$ and $10^{-11}$ Pa$^{-1}$ for all compositions, respectively. Rheological parameters (pre-exponential factor (A), stress exponent (n), and creep activation energy (Q)) are set according to flow laws of mica$^{59}$, Westerly granite$^{24}$, mafic granulite$^{25}$, Maryland diabase$^{26}$, dry olivine$^{60}$, and serpentinite$^{61}$, from top to bottom in the table. Other material properties are $\rho$, $k$, $Q_r$, $\alpha$, $C$, and $\phi$.
$\rho$ density, $k$ thermal conductivity, $Q_r$ radiogenic heat production, $\alpha$ coefficient of thermal expansion, $C$ cohesion, $\phi$ friction angle.

<table><caption>**Table 2 Thermal and rheological parameters used for different compositions in the Mod 1 model variant.**</caption><thead><tr><th></th><th>ρ (kg m<sup>-3</sup>)</th><th>k (W m<sup>-1</sup> K<sup>-1</sup>)</th><th>Q<sub>r</sub> (W m<sup>-3</sup>)</th><th>α (K<sup>-1</sup>)</th><th>C (MPa)</th><th>φ (°)</th><th>A (Pa<sup>-1</sup> s<sup>-1</sup>)</th><th>n</th><th>Q (J mol<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Sedimentary cover (mica)</td><td>2700</td><td>2.55</td><td>2.9e-6</td><td>3.0e-5</td><td>10</td><td>15</td><td>1.0e-138</td><td>18</td><td>51.0e3</td></tr><tr><td>Continental upper crust (wet quartzite)</td><td>2750</td><td>2.8</td><td>0.5e-6</td><td>3.0e-4</td><td>10</td><td>30</td><td>5.0717e-18</td><td>2.3</td><td>154.5e3</td></tr><tr><td>Continental lower crust (felsic granulite)</td><td>2900</td><td>2.8</td><td>0.5e-6</td><td>3.0e-4</td><td>10</td><td>30</td><td>2.0095e-21</td><td>4.2</td><td>243.0e3</td></tr><tr><td>Oceanic crust (Maryland diabase)</td><td>2900</td><td>3.0</td><td>1.0e-10</td><td>3.0e-5</td><td>10</td><td>30</td><td>3.2e-20</td><td>3.0</td><td>276.0e3</td></tr><tr><td>Lithospheric mantle (dry olivine)</td><td>3300</td><td>3.0</td><td>1.0e-10</td><td>3.0e-5</td><td>10</td><td>30</td><td>1.1e-16</td><td>3.5</td><td>530.0e3</td></tr><tr><td>Asthenosphere (dry olivine)</td><td>3300</td><td>3.0</td><td>1.0e-10</td><td>3.0e-5</td><td>10</td><td>30</td><td>1.1e-16</td><td>3.5</td><td>530.0e3</td></tr><tr><td>Weak zone (serpentinite)</td><td>2900</td><td>3.0</td><td>1.0e-10</td><td>3.0e-5</td><td>0</td><td>30</td><td>4.4738e-38</td><td>3.8</td><td>8.9e3</td></tr></tbody></table>

The heat capacity ($C_p$) and the compressibility ($\beta$), were set to 1050 J kg$^{-1}$ K$^{-1}$ and $10^{-11}$ Pa$^{-1}$ for all compositions, respectively. Rheological parameters (pre-exponential factor (A), stress exponent (n), and creep activation energy (Q)) are set according to flow laws of mica$^{59}$, wet quartzite$^{25}$, felsic granulite$^{25}$, Maryland diabase$^{26}$, dry olivine$^{60}$, and serpentinite$^{61}$, from top to bottom in the table. Other material properties are $\rho$, $k$, $Q_r$, $\alpha$, $C$, and $\phi$.
$\rho$ density, $k$ thermal conductivity, $Q_r$ radiogenic heat production, $\alpha$ coefficient of thermal expansion, $C$ cohesion, $\phi$ friction angle.

In the reference model, a decoupled crustal rheology is achieved by using the Westerly granite flow law$^{24}$ for the continental upper crust, and a mafic granulite flow law$^{25}$ for the lower crust (Fig. 2a). Using a stronger Maryland diabase flow law$^{26}$ for the upper crust results in a more coupled, stronger crustal rheology (Fig. 2a). This prevents upper crustal decoupling after reaching eclogite-facies conditions and inhibits subsequent nappe formation and extrusion (Fig. 4a). Instead, the upper crust is subducted to greater depth and enhanced underplating occurs compared to the reference model.

A more coupled, but relatively weak crustal rheology can be achieved by decreasing the strength of the continental lower crust (Figs. 2a and 4b). Such a rheology results in the decoupling of the entire continental crust from the lithospheric mantle rather than the decoupling of the upper crust from the lower crust (Fig. 4b). This leads to distributed folding and thrusting in the lower plate rather than localized nappe formation and extrusion of the subducted upper crust. Our results thus indicate that both strong (Fig. 4a) and weak (Fig. 4b) coupled crustal rheologies inhibit upper crustal extrusion and associated the emplacement of far-travelled ophiolites. We further tested different types of decoupled crustal rheological models to determine the effect of smaller compositional or thermal differences. The results show that slight variations in the thermal and material properties lead to different timing and thus position of crustal decoupling, which is reflected in different amounts of underplated continental upper crust below the oceanic upper plate (Supplementary Fig. 3). Rheology and coupling of the continental crust are hence key factors that control nappe formation, upper crust extrusion, and upper plate necking.

The degree of coupling in the subducting continental crust is also affected by shear heating, a process that dissipates mechanical energy in the form of heat, facilitating shear zone formation in the lithosphere$^{23,27}$. In our reference model shear heating results in heat production that can be up to two orders of

magnitude larger than the average radiogenic heat production in the continental crust (Fig. 5b, see Table 1 for radiogenic heat production ($Q_r$) values). This effect is transient and has limited spatial extent, but it induces thermal softening of the ductile material and may trigger the formation of a shear zone. This is documented by a reduction of the effective viscosity at the base of the incipient nappe at the boundary of the lower and upper continental crust and further along the steeper ramp section (Fig. 5). This allows an upper crustal nappe to decouple from the rest of the subducting lithosphere and extrude upwards and leftwards along the shear zone (Figs. 3b and 5). The model which does not include shear heating (Fig. 4c) shows that strain localization and nappe formation fails to initiate, or is significantly delayed compared to the reference model (Fig. 3). Without shear heating, the coupling between the upper crust and the subducting lithosphere is too high to allow for nappe formation and its subsequent extrusion. Instead, buoyancy force leads to the underplating of the subducted upper crust. This demonstrates that shear heating is essential for strain localization and nappe formation in our reference model.

**Comparison to natural ophiolite belts.** Our models display many first-order features of natural ophiolite belts that are related to continental subduction below an oceanic upper plate (Fig. 6). It produces a far-travelled ophiolite sheet as the structurally highest unit, separated from its root by buried and then exhumed HP-LT continental rocks, and underlain by accreted non-metamorphic to low-grade sedimentary cover units (Figs. 1a, b and 6). The tectonic window of exhumed continental crust shows a gradual decrease in metamorphic grade from blueschist-eclogite facies to non-metamorphic towards the far-travelled ophiolite sheet (Fig. 6). The reference model has a relatively warm initial continental geotherm ($T_{MOHO} = 510$ °C), and therefore shows higher peak metamorphic grades than the Mod 1 variant, which has an initially colder continent ($T_{MOHO} = 375$ °C) (Fig. 6a, b, respectively).
<figure><img src="image_6.png" /><figcaption>**Fig. 3 Kinematic and dynamic evolution of continental upper crust extrusion and upper plate necking.** Top panels show the compositional domains with the interpreted main structures (shear bands in the brittle domain are interpreted as faults), the middle panels provide information of the strain rate with velocity arrows, while the bottom panels on the second stress invariant with the orientation of the most compressive principal stress (σ₁). The color fill of the middle and bottom panels are based on scientific colormaps⁵⁸. The orientations of σ₁ on bottom panels are shown by black lines, horizontal for compression, and vertical for extension. The boundaries of different compositional domains presented on the top panels are highlighted with white contour lines on the lower right panels. **a** Decoupling of the upper crust and initiation of upper plate extension. **b** Further localization of main thrust and extensional accommodation structures that break the upper plate. **c** Emplacement of far-travelled ophiolite sheet by further crustal extrusion.</figcaption></figure>

As such, the reference model closer resembles natural sites where the tectonic window largely consists of metamorphosed crust. Such a case is observed in New Caledonia or Western-Central Anatolia, where the majority of the exposed continental crust reached *P* > 1 GPa and *T* between 300 and 600 °C¹⁴,²⁸⁻³⁰. On the other hand, Mod 1 fits better to natural examples where large

fractions of the tectonic window exhibit no or very low-grade metamorphic overprint, with relatively minor volumes of high-*T* blueschist and eclogite facies rocks (e.g. the Oman or Lesser Caucasus ophiolite belts¹²,³¹). In both cases, the highest-grade continental rocks are separated from the oceanic lithosphere to the right by a major normal-sense shear zone that accommodated
the extrusion of the subducted continental crust (Fig. 6). The ophiolite sheet and the underlying accreted sedimentary units are resting on the upper flat thrust segment of the main shear zone system that further transported the upper plate units (far-travelled ophiolite sheet and accreted sediments) leftwards during the extrusion of the continental crust. The normal-sense shear zone that separated the far-travelled ophiolite sheet from its root

<figure><img src="image_7.png" /><figcaption>**Fig. 4** Effects of changing key parameters compared to the reference model setup. a Model compositions at 29 Myr with stronger upper crust rheology (Maryland diabase instead of Westerly granite)²⁴,²⁶. For the initial strength profile of the continental domain see Fig. 2a. b Model compositions at 29 Myr with weaker lower crust rheology (Maryland diabase instead of mafic granulite)²⁵,²⁶. For the initial strength profile of the continental domain see Fig. 2a. c Model compositional domains at 29 Myr with shear heating switched off.</figcaption></figure>

during the extrusion joins this upper flat thrust segment. The prograde P–T ratio and peak P–T conditions recorded by the deeper parts of the subducted upper crust in our reference model (4 MPa/°C and ~2 GPa at 500–600 °C, respectively) are average values compared to those of natural cases (Figs. 1c and 2d and Supplementary Table 1). The duration of the continental subduction-exhumation cycle (15–20 Myr) is in agreement with the well-constrained ophiolite belts like Oman or New Caledonia, and appears to be slightly shorter than the sites where subduction and/or exhumation of the continental formations is poorly dated or debated (e.g. Brooks Range, Hellenides-Dinarides, and Southern Ural; see Fig. 1e, Supplementary Note 1, and Supplementary Table 2). Observations such as post-subduction ductile to brittle extensional deformation¹³,³², or the coexistence of opposite shear sense directions in the former passive margin units¹¹,³³ fit well with the structural evolution of our model, which involves top-left shearing (thrusting) during burial and both top-left and top-right shearing during exhumation (normal faulting; Figs. 3 and 6). The width of the far-travelled ophiolite sheet (50 km) predicted by the model(s) agrees very well with the average width (53 km) of natural ophiolite belts (Fig. 1e and Supplementary Note 1).

Hence, the structure, dimensions, composition, and P-T-t conditions of our model closely resemble those of natural ophiolite belts. The geometry and time evolution of our model also allows comparison with the currently active continental subduction of the Australian continental margin below the oceanic Banda arc. The Australian continental margin has been subducting for ~10 Myr³⁴, which roughly equals the duration of upper crustal subduction in our reference model. Reconstruction of the stacked and accreted sedimentary cover shows that 215–230 km of continental lithosphere has been subducted, but only the upper-most 2 km of sedimentary cover was accreted to the upper plate³⁵. Hence the present-day structure and dimensions of the Australian continental subduction are very similar to those of our reference model in the 20–23-Myr snapshots (Fig. 2c).

## Discussion

Our results have important implications for the dynamics of ophiolite emplacement in oceanic upper plate-continental lower plate subduction systems, and thus may be widely applied to understand geological observations in natural ophiolite belts. While plate kinematic changes may play an important role in initiating intra-oceanic subduction³⁶, or in the cessation of contraction at continental subduction zones⁶, we show that the emplacement of far-travelled ophiolite sheets can result from syn-convergence, buoyancy-driven decoupling, and upward extrusion of the subducted continental upper crust accommodated by the necking of the oceanic upper plate. Crustal extrusion and the

<figure><img src="image_8.png" /><figcaption>**Fig. 5** The effect of shear heating. a Plot of visco-plastic dissipation of mechanical energy in the form of heat (Qs; left panel) and effective viscosity (ηeff, right panel) at the 21.6 Myr snapshot of the reference model, prior to nappe formation in the subducting continental crust. b Plots of visco-plastic dissipation of mechanical energy in the form of heat (left panel) and effective viscosity (right panel) at the 23.7 Myr snapshot of the reference model, after nappe formation at in the subducted continental crust. Color fills of the plots are based on scientific colormaps⁵⁸.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 6** Near-surface composition and structure of final model snapshots. The peak metamorphic zonation of the exhumed continental upper crust is plotted for a the reference model and b the Mod 1 model variant (for model parameters see Table 2, for initial strength profile see Supplementary Fig. 3). The panels are a combination of calculated metamorphic domains for the continental upper crust, compositional domains for the sedimentary cover and the oceanic lithosphere, and line interpretations of the main shear zones. The peak metamorphic zones are defined based on the recorded maximum pressure (P) and temperature (T) values of the upper crustal particles throughout the entire model evolution, and hence represent the theoretical maximum metamorphic grade reached by the given crustal domains. In both cases, a transition from non-metamorphic to eclogite-grade metamorphic rocks are observed within the upper crustal tectonic window. The reference model (Fig. 6a) has a higher initial continental geothermal gradient (TMOHO = 510 °C), therefore it shows higher peak metamorphic grades compared to the model (Mod 1) with colder Moho temperature (TMOHO = 375 °C; Fig. 6b).</figcaption></figure>

resulting necking of the upper plate in our model are consistent with cases studies from Oman¹⁷ and New Caledonia¹⁵. Extrusion of the subducted upper crust requires nappe formation. In agreement with previous studies²⁶,³⁰, the results show that shear heating is an important mechanism that facilitates strain localization and nappe formation. The precise reproduction of smaller scale nappes (nappe thickness of several kilometers) which is often observed in case of continental subduction⁹,¹³ would require very high-resolution numerical modeling and built-in heterogeneities inside the upper crust to localize shear zones at multiple horizons³⁷⁻³⁹. Our results also indicate that crustal decoupling and exhumation may take place with different timing and position in the subducted continent depending on the rheology of the continental crust⁴⁰,⁴¹. Variations in thermal or compositional properties thus might control the surface preservation (exhumation) or the subduction and recycling of different types of continental passive margins (e.g. magmatic or magma-poor passive margins)⁴²⁻⁴⁴.

Numerous natural ophiolites show evidence for subduction zone magmatism in the upper plate following intra-oceanic subduction initiation, which results in thermally younger, thus thinner upper plates⁸,⁴⁵⁻⁴⁷. Our model does not account for such effects and hence may overestimate upper-plate thickness. Thinner upper plates would result in flatter continental subduction, which may lead to crustal decoupling and upper plate necking further away from the ophiolite front. If so, the size of the resulting far-travelled ophiolite sheets would be comparable to

the widest far-travelled ophiolite sheets in Oman and Anatolia (80–150 km; see Supplementary Note 1).

The current subduction of the Australian plate below the oceanic Banda arc provides an exciting example for a prospective future ophiolite belt. As more than 200 km of continental crust has already subducted below the oceanic plate³⁵, it most likely has reached eclogite facies conditions. Based on our model, decoupling of the Australian upper crust has already begun, or will begin in the geological near future. If decoupling is followed by nappe formation, the extrusion of the upper crust and simultaneous necking of the oceanic upper plate may lead to the emplacement of far-travelled ophiolites.

## Methods

**Numerical modelling.** The above-presented modeling results were obtained by solving the coupled set of non-linear thermo-mechanical equations. The steady-state momentum equation, the heat transfer equation and the incompressible mass conservation equations are formulated as:

$$ \frac{\partial \tau_{ij}}{\partial x_j} - \frac{\partial P}{\partial x_i} = -\rho g_i \qquad (1) $$

$$ \rho c_p \frac{DT}{Dt} = \frac{\partial}{\partial x_i} \left( k \frac{\partial T}{\partial x_i} \right) + Q_r + Q_s \qquad (2) $$

$$ \frac{\partial v_i}{\partial x_i} = 0 \qquad (3) $$

where v is the velocity vector, T is the temperature, k is the thermal conductivity, ρ is the density, c_p is the heat capacity at constant pressure, Q_r is the radiogenic heat
production, τ is the deviatoric stress tensor, ˙e is the deviatoric strain rate tensor, P is the pressure, and g is the gravity acceleration vector. The visco-plastic dissipation or shear heating, Qs, is expressed as:

$$Q_s = \tau_{ij} (\dot{\epsilon}_{ij} - \dot{\epsilon}_{ij}^e) \quad (4)$$

where $\dot{\epsilon}_{ij}^e$ the elastic portion of the deviatoric strain rate.

The density field evolves according the following equation of state:

$$\rho = \rho_0 (1 - \alpha(T - T_0))(1 + \beta(P - P_0)) \quad (5)$$

where $\rho_0$ is the reference density, $\alpha$ is the thermal expansivity, $\beta$ is the compressibility, $T_0$ and $P_0$ are the reference temperature and pressure which were respectively set to 0 °C and 105 Pa. We use the Boussinesq approximation; hence both density changes in the mass conservation equation and adiabatic heating are neglected.

The effective viscosity ($\eta$) relates to the deviator stress and strain rate tensor:

$$\tau_{ij} = 2\eta\dot{\epsilon}_{ij} = \left(\frac{1}{\eta}\right)^{-1} \left(\frac{1}{\eta^e} + \frac{1}{\eta^p}\right)^{-1} \dot{\epsilon}_{ij} \quad (6)$$

using a visco-elasto-plastic rheological model:

$$\dot{\epsilon}_{ij} = \dot{\epsilon}_{ij}^v + \dot{\epsilon}_{ij}^e + \dot{\epsilon}_{ij}^p \quad \text{where} \quad \dot{\epsilon}_{ij}^{\text{dis}} = \dot{\epsilon}_{ij}^{\text{dis}} + \dot{\epsilon}_{ij}^{\text{Peierls}} \quad (7)$$

where the v, e, and p superscripts stand for viscous, elastic and plastic and the superscripts dis and Peierls correspond to dislocation and Peierls creep mechanisms.

The viscous strain rate account for contributions of both Peierls and dislocation creep. The contribution of dislocation creep is expressed as:

$$\dot{\epsilon}_{ij}^{\text{dis}} = \dot{\epsilon}_{\text{II}}^{\text{dis}} \frac{\tau_{ij}}{\tau_{\text{II}}} = \left(2A^{\frac{1}{n}} f e^{\frac{Q}{k_B T}}\right)^{-n} \frac{\tau_{ij}}{\tau_{\text{II}}^n} \quad (8)$$

where A is a pre-factor, Q is the activation energy, n is the stress exponent, R is the universal gas constant, and f is a correction factor$^{48}$. The subscripts II stand for the square root of the second tensor invariant. The Peierls mechanism is taken into account in the mantle using the regularized formulation$^{49}$. The Peierls strain rate is expressed as:

$$\dot{\epsilon}_{ij}^{\text{Peierls}} = \frac{\epsilon_{\text{II}}^{\text{Peierls}} \tau_{ij}}{\tau_{\text{II}}} \quad (9)$$

where the second invariant strain rate is spelled as:

$$\left\{ \begin{array}{l} \dot{\epsilon}_{\text{II}}^{\text{Peierls}} = (2A^{\text{Peierls}})^{-s} \frac{\tau_{ij}}{\tau_{\text{II}}^s} \frac{\tau_{ij}}{\tau_{\text{II}}} \\ A^{\text{Peierls}} = f^{\text{Peierls}} \gamma \sigma^{\text{Peierls}} \left( E^{\text{Peierls}} e^{\frac{(1-\gamma)^2 Q^{\text{Peierls}}}{RT}} \right)^{-1} \\ s = \frac{Q^{\text{Peierls}}}{RT} (1-\gamma)^{(q-1)q\gamma} \end{array} \right. \quad (10)$$

where s is an effective temperature-dependent stress exponent, QPeierls is the activation energy (=540 kJ/mol), $\sigma^{\text{Peierls}}$ is the Peierls stress (=8.5 × 10⁹ Pa), $E^{\text{Peierls}}$ (=5.7 × 10¹¹ s⁻¹), q (=2.0), and $\gamma$ (=0.1)$^{50}$.

The elastic strain rate is written as:

$$\dot{\epsilon}_{ij}^{e} = \frac{\tau_{ij} \widehat{\tau}_{ij}}{2G \tau_{\text{II}}} \quad (11)$$

where $\tau_{ij}$ is the corotational time derivative of the stress tensor (Jaumann rate) and G is the shear modulus (=10¹⁰ Pa).

The plastic strain rate takes the form of:

$$\dot{\epsilon}_{ij}^{p} = \dot{\epsilon}_{\text{II}}^{p} \frac{\tau_{ij}}{\tau_{\text{II}}} \quad \text{with} \quad \dot{\epsilon}_{\text{II}}^{p} = \frac{F}{2\eta^{ve}} \quad \text{and is defined only if} \quad F = \tau_{\text{II}} - C \cos \varphi - P \sin \varphi > 0 \quad (12)$$

where $\varphi$ is the friction angle and C is the cohesion (see Table 1 for reference model values). The coefficient $\eta^{ve}$ is a visco-elastic coefficient dependent on both viscous and elastic moduli and on the time discretization$^{51}$. No plastic strain softening was applied in the presented simulations.

The temperature is fixed at both the surface of the model (0 °C) and the lower boundary (1330 °C). No heat flows through the right and left boundaries. A plate convergence rate of 3 cm/year is achieved by prescribing constant normal inflow velocities of $|V_{\text{in}}| = 1.5$ cm/year along the upper 140 km of the two model sides, while mass conservation is satisfied by gradually increasing outflow below 140 km. The shear stress is set to zero along the left, right, and lower boundaries. The upper boundary is a true free surface$^{18}$, thus its position evolves in response to tectonic loading. The initial topography is set to 0 km. The initial temperature field is computed by solving the heat transfer equation assuming steady state and neglecting visco-plastic dissipation using the reference thermal parameters (Table 1). During this initialization step, quasi-adiabatic mantle conditions are reached by assuming an artificially high conductivity within the asthenosphere.

The thermo-mechanical equations are solved using the finite difference/marker-in-cell technique$^{52}$. The global linearized system of mechanical equations is solved with a direct-iterative scheme that combines both Powell-Hestenes iterations and Cholesky factorization$^{53}$. Due to the non-linear nature of the considered

rheological model, non-linear iterations are required at both local and global levels. At the local level, an exact partitioning of the strain rate tensor is obtained via successive Newton iterations$^{54,55}$. At the global level, Picard iterations are used to best-satisfy mechanical equilibrium equations (to an absolute tolerance of 10⁻⁶ and within a maximum of 20 iterations).

**Model geometry.** The computational domain is a cross section of 1330 × 410 km. The model resolution is 1 km in both directions. The initial compositional geometry is inspired by reconstructions of pre-obduction geodynamic settings. It contains an oceanic domain (660 km wide) with a spreading ridge and a tilted weak zone in the center that ensures left-dipping intra-oceanic subduction initiation. The thermal structure of the oceanic lithosphere is calculated by applying a half-space cooling age model from 1.5 Myr at the center to 50 Myr at the edges of the ocean. The oceanic crust is 6 km thick and is overlain by a layer of uppermost sedimentary cover that linearly thickens from the ridge (0 km) towards the right edge of the continental domain (3 km). The transition from the continent to the ocean is defined by a passive margin geometry where the continental upper and lower crust linearly thin from 30 km to 5 km over the distance of 200 km. The uppermost sedimentary cover layer has a constant 3 km thickness over the continental domain.

## Data availability

Due to the very large size of the files, the model output data are available upon request from the corresponding author.

## Code availability

The numerical code used in this study is available from the main code developer (Thibault Duretz, e-mail: thibault.duretz@univ-rennes1.fr) upon reasonable request.

Received: 17 July 2020; Accepted: 15 February 2021; Published online: 08 March 2021

## References

1. Coleman, R. Plate tectonic emplacement of upper mantle peridotites along continental edges. *J. Geophys. Res.* **76**, 1212–1222 (1971).

2. Dewey, J. Ophiolite obduction. *Tectonophysics* **31**, 93–120 (1976).

3. Agard, P. et al. Obduction: why, how and where. Clues from analog models. *Earth Planet. Sci. Lett.* **393**, 132–145 (2014).

4. Oxburgh, E. Flake tectonics and continental collision. *Nature* **239**, 202–204 (1972).

5. Boudier, F., Ceuleneer, G. & Nicolas, A. Shear zones, thrusts and related magmatism in the Oman ophiolite: initiation of thrusting on an oceanic ridge. *Tectonophysics* **151**, 275–296 (1988).

6. Duretz, T. et al. Thermo-mechanical modeling of the obduction process based on the Oman ophiolite case. *Gondwana Res.* **32**, 1–10 (2016).

7. Agard, P. et al. Plate interface rheological switches during subduction infancy: control on slab penetration and metamorphic sole formation. *Earth Planet. Sci. Lett.* **451**, 208–220 (2016).

8. van Hinsbergen, D. J. et al. Dynamics of intraoceanic subduction initiation: 2. Suprasubduction zone ophiolite formation and metamorphic sole exhumation in context of absolute plate motions. *Geochem. Geophys. Geosyst.* **16**, 1771–1785 (2015).

9. Kilias, A. et al. Alpine architecture and kinematics of deformation of the northern Pelagonian nappe pile in the Hellenides. *Austrian J. Earth Sci.* **103**, 4–28 (2010).

10. Yamato, P. et al. New, high-precision P–T estimates for Oman blueschists: implications for obduction, nappe stacking and exhumation processes. *J. Metamorphic Geol.* **25**, 657–682 (2007).

11. van Hinsbergen, D. J. et al. Tectonic evolution and paleogeography of the Kırşehir Block and the Central Anatolian Ophiolites, Turkey. *Tectonics* **35**, 983–1014 (2016).

12. Hässig, M., Rolland, Y., Duretz, T. & Sosson, M. Obduction triggered by regional heating during plate reorganization. *Terra Nova* **28**, 76–82 (2016).

13. Searle, M. P. Structural geometry, style and timing of deformation in the Hawasina Window, Al Jabal al Akhdar and Saïh Hatat culminations, Oman Mountains. *GeoArabia* **12**, 99–130 (2007).

14. Pourteau, A. et al. Neotethys closure history of Anatolia: insights from 40Ar–39Ar geochronology and P–T estimation in high-pressure metasedimentary rocks. *J. Metamorphic Geol.* **31**, 585–606 (2013).

15. Lagabrielle, Y., Chauvet, A., Ulrich, M. & Guillot, S. Passive obduction and gravity-driven emplacement of large ophiolitic sheets: The New Caledonia ophiolite (SW Pacific) as a case study? *Bull. Soc. Géol. Fr.* **184**, 545–556 (2013).
16. Agard, P., Searle, M. P., Alsop, G. I. & Dubacq, B. Crustal stacking and expulsion tectonics during continental subduction: P-T deformation constraints from Oman. *Tectonics* **29**, TC5018 (2010).

17. Chemenda, A. I., Mattauer, M. & Bokun, A. N. Continental subduction and a mechanism for exhumation of high-pressure metamorphic rocks: new modelling and field data from Oman. *Earth Planet. Sci. Lett.* **143**, 173–182 (1996).

18. Duretz, T., May, D. A. & Yamato, P. A free surface capturing discretization for the staggered grid finite difference scheme. *Geophys. J. Int.* **204**, 1518–1530 (2016).

19. Searle, M. & Cox, J. Tectonic setting, origin, and obduction of the Oman ophiolite. *Geol. Soc. Am. Bull.* **111**, 104–122 (1999).

20. Maffione, M. & van Hinsbergen, D. J. Reconstructing plate boundaries in the Jurassic neo-Tethys from the east and west Vardar ophiolites (Greece and Serbia). *Tectonics* **37**, 858–887 (2018).

21. Maffione, M., van Hinsbergen, D. J., de Gelder, G. I., van der Goes, F. C. & Morris, A. Kinematics of Late Cretaceous subduction initiation in the Neo-Tethys Ocean reconstructed from ophiolites of Turkey, Cyprus, and Syria. *J. Geophys. Res. Solid Earth* **122**, 3953–3976 (2017).

22. Maffione, M. et al. Dynamics of intraoceanic subduction initiation: 1. Oceanic detachment fault inversion and the formation of supra-subduction zone ophiolites. *Geochem. Geophys. Geosyst.* **16**, 1753–1770 (2015).

23. Kiss, D., Podladchikov, Y., Duretz, T. & Schmalholz, S. M. Spontaneous generation of ductile shear zones by thermal softening: Localization criterion, 1D to 3D modelling and application to the lithosphere. *Earth Planet. Sci. Lett.* **519**, 284–296 (2019).

24. Hansen, F. & Carter, N. in *The 24th US Symposium on Rock Mechanics (USRMS)* (American Rock Mechanics Association, 1983).

25. Ranalli, G. *Rheology of the Earth* (Springer Science & Business Media, 1995).

26. Carter, N. L. & Tsenn, M. C. Flow properties of continental lithosphere. *Tectonophysics* **136**, 27–63 (1987).

27. Duretz, T., Schmalholz, S. & Podladchikov, Y. Shear heating-induced strain localization across the scales. *Philos. Mag.* **95**, 3192–3207 (2015).

28. Potel, S., Mählmann, R. F., Stern, W., Mullis, J. & Frey, M. Very low-grade metamorphic evolution of pelitic rocks under high-pressure/low-temperature conditions, NW New Caledonia (SW Pacific). *J. Petrol.* **47**, 991–1015 (2006).

29. Brovarone, A. V. & Agard, P. True metamorphic isograds or tectonically sliced metamorphic sequence? New high-spatial resolution petrological data for the New Caledonia case study. *Contrib. Mineral. Petrol.* **166**, 451–469 (2013).

30. Okay, A. Jadeite–chloritoid–glaucothane–lawsonite blueschists in north-west Turkey: unusually high P/T ratios in continental crust. *J. Metamorph. Geol.* **20**, 757–768 (2002).

31. Searle, M., Warren, C., Waters, D. & Parrish, R. Structural evolution, metamorphism and restoration of the Arabian continental margin, Saïh Hatat region, Oman Mountains. *J. Struct. Geol.* **26**, 451–473 (2004).

32. Rawling, T. J. & Lister, G. S. Large-scale structure of the eclogite–blueschist belt of New Caledonia. *J. Struct. Geol.* **24**, 1239–1258 (2002).

33. Porkoláb, K. et al. Cretaceous-Paleogene tectonics of the Pelagonian zone: inferences from Skopelos island (Greece). *Tectonics* **38**, 1946–1973 (2019).

34. Keep, M. & Haig, D. W. Deformation and exhumation in Timor: distinct stages of a young orogeny. *Tectonophysics* **483**, 93–111 (2010).

35. Tate, G. W. et al. Australia going down under: quantifying continental subduction during arc-continent accretion in Timor-Leste. *Geosphere* **11**, 1860–1883 (2015).

36. Agard, P., Jolivet, L., Vrielynck, B., Burov, E. & Monie, P. Plate acceleration: the obduction trigger? *Earth Planet. Sci. Lett.* **258**, 428–441 (2007).

37. Kiss, D., Duretz, T. & Schmalholz, S. Tectonic inheritance controls nappe detachment, transport and stacking in the Helvetic Nappe System, Switzerland: insights from thermo-mechanical simulations. *Solid Earth* **11**, 287–305 (2020).

38. Duretz, T. et al. The importance of structural softening for the evolution and architecture of passive margins. *Sci. Rep.* **6**, 1–7 (2016).

39. Reuber, G., Kaus, B. J., Schmalholz, S. M. & White, R. W. Nonlithostatic pressure during subduction and collision and the formation of (ultra) high-pressure rocks. *Geology* **44**, 343–346 (2016).

40. Willingshofer, E., Sokoutis, D., Luth, S., Beekman, F. & Cloetingh, S. Subduction and deformation of the continental lithosphere in response to plate and crust-mantle coupling. *Geology* **41**, 1239–1242 (2013).

41. Vogt, K., Matenco, L. & Cloetingh, S. Crustal mechanics control the geometry of mountain belts. Insights from numerical modelling. *Earth Planet. Sci. Lett.* **460**, 12–21 (2017).

42. Tugend, J. et al. Reappraisal of the magma-rich versus magma-poor rifted margin archetypes. *Geol. Soc. Lond. Spec. Publ.* **476**, 23–47 (2020).

43. Franke, D. Rifting, lithosphere breakup and volcanism: Comparison of magma-poor and volcanic rifted margins. *Mar. Petrol. Geol.* **43**, 63–87 (2013).

44. Yamato, P., Burov, E., Agard, P., Le Pourhiet, L. & Jolivet, L. HP-UHP exhumation during slow continental subduction: Self-consistent thermodynamically and thermomechanically coupled model with application to the Western Alps. *Earth Planet. Sci. Lett.* **271**, 63–74 (2008).

45. Tamura, A. & Arai, S. Harzburgite–dunite–orthopyroxenite suite as a record of supra-subduction zone setting for the Oman ophiolite mantle. *Lithos* **90**, 43–56 (2006).

46. Meffre, S., Aitchison, J. C. & Crawford, A. J. Geochemical evolution and tectonic significance of boninites and tholeiites from the Koh ophiolite, New Caledonia. *Tectonics* **15**, 67–83 (1996).

47. Pearce, J. A., Lippard, S. & Roberts, S. Characteristics and tectonic significance of supra-subduction zone ophiolites. *Geol. Soc. Lond. Spec. Publ.* **16**, 77–94 (1984).

48. Schmalholz, S. M. & Fletcher, R. C. The exponential flow law applied to necking and folding of a ductile layer. *Geophys. J. Int.* **184**, 83–89 (2011).

49. Kameyama, M., Yuen, D. A. & Karato, S.-I. Thermal-mechanical effects of low-temperature plasticity (the Peierls mechanism) on the deformation of a viscoelastic shear zone. *Earth Planet. Sci. Lett.* **168**, 159–172 (1999).

50. Evans, B. & Goetze, C. The temperature variation of hardness of olivine and its implication for polycrystalline yield stress. *J. Geophys. Res. Solid Earth* **84**, 5505–5524 (1979).

51. Moresi, L., Dufour, F. & Mühlhaus, H.-B. A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials. *J. Comput. Phys.* **184**, 476–497 (2003).

52. Gerya, T. V. & Yuen, D. A. Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties. *Phys. Earth Planet. Inter.* **140**, 293–318 (2003).

53. Räss, L., Duretz, T., Podladchikov, Y. Y. & Schmalholz, S. M. M2Di: Concise and efficient MATLAB 2-DS Stokes solvers using the Finite Difference Method. *Geochem. Geophys. Geosyst.* **18**, 755–768 (2017).

54. Popov, A. & Sobolev, S. V. SLIM3D: a tool for three-dimensional thermomechanical modeling of lithospheric deformation with elasto-visco-plastic rheology. *Phys. Earth Planet. Inter.* **171**, 55–75 (2008).

55. Schmalholz, S. M. & Duretz, T. Impact of grain size evolution on necking in calcite layers deforming by combined diffusion and dislocation creep. *J. Struct. Geol.* **103**, 37–56 (2017).

56. Patriat, M. et al. New Caledonia obducted Peridotite Nappe: offshore extent and implications for obduction and postobduction processes. *Tectonics* **37**, 1077–1096 (2018).

57. Maurizot, P., Robineau, B., Vendé-Leclerc, M. & Cluzel, D. *Introduction to New Caledonia: Geology, Geodynamic Evolution And Mineral Resources* (Geological Society of London, 2020).

58. Crameri, F. Scientific colour-maps. Zenodo https://zenodo.org/record/3596401#.YDSeDugzbcS (2018).

59. Kronenberg, A. K., Kirby, S. H. & Pinkston, J. Basal slip and mechanical anisotropy of biotite. *J. Geophys. Res. Solid Earth* **95**, 19257–19278 (1990).

60. Hirth, G. & Kohlstedt, D. *Rheology of The Upper Mantle And The Mantle Wedge: A View From The Experimentalists*. Vol. 138, p. 83–106 (Geophysical monograph-american geophysical union, 2003)..

61. Hilairet, N. et al. High-pressure creep of serpentine, interseismic deformation, and initiation of subduction. *Science* **318**, 1910–1913 (2007).

## Acknowledgements

The research leading to these results has received funding from the European Union's MSCA-ITN-ETN Project SUBITOP 674899. Numerical simulations were performed on the Utrecht University cluster Eejit. Thanks are due to Jeroen van Hunen and Cedric Thieulot for discussions that helped to initiate the project.

## Author contributions

K.P. conceived the research idea. T.D. and P.Y. designed the thermo-mechanical numerical code. A.A, P.Y., T.D., and K.P. designed the model setup. K.P. conducted the numerical simulations and interpreted the results together with P.Y., T.D., and E.W. All authors discussed the results and interpretations, and contributed to writing the paper.

## Competing interests

The authors declare no competing interests.
**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-021-21866-1.

**Correspondence** and requests for materials should be addressed to K.P.

**Peer review information** *Nature Communications* thanks William Griffin, Thomas Lamont and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_10.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2021