# Garnet microstructures suggest ultra-fast decompression of ultrahigh-pressure rocks

Cindy Luisier, Lucie Tajčmanová, Philippe Yamato, Thibault Duretz

► To cite this version:

Cindy Luisier, Lucie Tajčmanová, Philippe Yamato, Thibault Duretz. Garnet microstructures suggest ultra-fast decompression of ultrahigh-pressure rocks. Nature Communications, 2023, 14 (1), pp.6012. 10.1038/s41467-023-41310-w. insu-04220716

HAL Id: insu-04220716

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

Submitted on 28 Sep 2023

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>

# Garnet microstructures suggest ultra-fast decompression of ultrahigh-pressure rocks

Received: 29 November 2022

Accepted: 30 August 2023

Published online: 27 September 2023

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

Check for updates

Cindy Luisier Ⓔ<fnref n="1" />, Lucie Tajčmanová Ⓔ², Philippe Yamato Ⓔ<fnref n="1" />,<fnref n="3" /> & Thibault Duretz<fnref n="1" />

Plate tectonics is a key driver of many natural phenomena occurring on Earth, such as mountain building, climate evolution and natural disasters. How plate tectonics has evolved through time is still one of the fundamental questions in Earth sciences. Natural microstructures observed in exhumed ultrahigh-pressure rocks formed during continental collision provide crucial insights into tectonic processes in the Earth's interior. Here, we show that radial cracks around SiO₂ inclusions in ultrahigh-pressure garnets are caused by ultrafast decompression. Decompression rates of at least 8 GPa/Myr are inferred independently of current petrochronological estimates by using thermo-mechanical numerical modeling. Our results question the traditional interpretation of fast and significant vertical displacement of ultrahigh-pressure tectonic units during exhumation. Instead, we propose that such substantial decompression rates are related to abrupt changes in the stress state of the lithosphere independently of the spatial displacement.

The understanding of plate tectonics processes, such as mountain building, usually relies on petrochronological data inferred from microscale observations in exhumed rock samples<fnref n="1" />⁻<fnref n="3" />. The microscale observations reveal relics of mineral transformations related to changes in pressure (P) and temperature (T) during the tectonic evolution experienced by the rock samples. In that sense, ultrahigh-pressure (UHP) rocks—i.e., rocks that underwent extreme P conditions in the deep Earth, documented by the formation of high-pressure silica (SiO₂) or carbon polymorphs (coesite or diamond, respectively)—are the key witnesses of deep processes during burial and exhumation of the Earth's lithosphere. However, the mechanisms of exhumation of UHP units remain one of the most challenging yet unsolved questions that keep on feeding controversies over the decades⁴⁻⁸. One of the first evidence of UHP metamorphism on Earth was provided by Chopin⁹, based on the observation of coesite relics in pyrope garnets from the Dora-Maira whiteschists (Brossasco-Isasca unit, Western European Alps). Thanks to its emblematic status, the UHP Dora Maira whiteschist has been the most studied UHP rock worldwide and can serve as a proxy for refining our knowledge of the burial and exhumation cycle⁴,⁹⁻²² (Fig. 1a). The P-T and time (t) petrochronological data provide robust constrains for building up geodynamic scenarios. The P-T-t dataset (Fig. 1a) indicates that the coesite-bearing Dora-Maira

whiteschists reached peak P-T conditions of 4.3 GPa and 730 °C¹⁶ at around 35.1 ± 0.9 Ma⁴ and then a first retrograde stage of 1.1 GPa and 600 °C¹⁶ at 32.9 ± 0.9 Ma⁴. Traditionally, the geochronological and petrological data were also used to estimate an average exhumation rate that corresponds to the velocity at which the rocks are exhumed toward the surface. The time is then given by geochronological data, and the variation of depth is computed assuming that pressure is directly proportional to depth (i.e., P is considered lithostatic). Hence, for Dora Maira, an average exhumation rate of ~5 cm/year can be estimated⁴. Several exhumation mechanisms of UHP rocks have been proposed in the past based on available petrological, structural and geochronological data²³⁻²⁸. However, exhumation from mantle depth involves displacement larger than 100 km, for which geological evidence remains scarce⁷,²⁹. Alternatively, it has been suggested that the major part of the decompression of (U)HP rocks does not translate into a vertical displacement and rather corresponds to a rapid switch in the tectonic stress state within collision zones³⁰.

The iconic petrographic observations in the UHP rocks are the partial preservation of UHP coesite inclusions in garnet that were partially to completely turned into α-quartz - the low-pressure SiO₂ polymorph—during retrogression and the associated radial fractures in the garnet host. The relationship between the coesite to α-quartz phase

<footnote n="1">¹Institute of Geosciences, Goethe University Frankfurt, Frankfurt am Main, Germany. ²Institute of Earth Sciences, Heidelberg University, Heidelberg, Germany.</footnote>

<footnote n="3">³Univ Rennes, CNRS, Géosciences Rennes, Rennes, France. ✉ e-mail: cindy.luisier@gmail.com</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1 | Compilation of available petrochronological and mechanical data.**</figcaption></figure>

a Summary of the pressure and temperature (P-T) data for the Dora Maira Massif. Age data for the peak metamorphic event (4.3 GPa, 730 °C)⁴ and retrograde events¹,¹⁴ are also reported. The black line corresponds to the quartz (q)-coesite (coe) transition, and the red dashed line shows the retrograde path calculated based on the P-T data¹⁶ (red symbols). Error bars indicate the minimum and maximum values of *P* and *T* estimates. b Sketch of the elastic model explaining the two stages of pressure evolution of an inclusion (black dashed line) entrapped in a

garnet during the retrograde P-T path (red dashed line). The first stage (along the coesite-quartz transition) corresponds to the pressure vessel effect. The second stage of pressure decrease occurs after garnet fracturing, when the threshold of $P_{inc} = 3 \times P_{host}$ is reached (see text for explanation). c Colormap of SiO₂ density in kg·m⁻³. The red dashed line and symbols are as in (a). Black curves show the brittle-ductile transition for different deviatoric strain rates (brittle and ductile on the left and right side of the curve, respectively) computed using the garnet flow law³⁰.

transition and the development of radial fractures in garnet may provide precious constraints on the exhumation of UHP units and the associated geodynamic evolution³¹. Previous studies focusing on the preservation of coesite have traditionally explained it by the pressure vessel effect, which relies on the differential volumetric expansion of the inclusion and host minerals during decompression³¹,³² (Fig. 1b). Garnet is a strong mineral that cannot accommodate the volume increase induced by the coesite to α-quartz transformation. Therefore, high pressure is maintained on the expanding inclusion, and the reaction is prevented from running to completion. However, these models only consider the elastic behavior of garnet (i.e., no permanent deformation and hence no fracturing) and therefore, radial fracturing was traditionally constrained to occur as soon as the pressure difference between the inclusions ($P_{inc}$) and the garnet host ($P_{host}$) overcomes a factor 3, i.e., $P_{inc} \ge 3P_{host}$³³,³⁴ (Fig. 1b). Alternatively, the growth of pre-defined cracks has been modeled by combining elasticity and fracture mechanics³⁵. Previous models usually used simple geometries involving either a single spherical inclusion in a spherical host³⁴ or faceted inclusion³⁵ in a host. The role of viscous relaxation and plastic deformation has been shown to affect the pressure distribution in decompressing inclusion-hosts systems³⁶. However, such analysis has never been quantitatively applied to investigate the fracturing of coesite-bearing UHP garnets.

Recent improvements in computing tools enable us to account for complex geometries (e.g., multiple mineral phases, digitalized thin sections). In contrast with past studies, elastic, viscous and plastic properties of the phases are now also incorporated, which allows for consistently simulating the coupling between reaction-induced volume changes and deformation³⁷. Furthermore, recent studies have shown that UHP rocks, such as eclogites and, in particular, the mineral garnet, can deform in a brittle manner at high *P* conditions,

where they are expected to be viscous³⁷,³⁸. In fact, the brittle-ductile transition (BDT) of garnet for a whole range of strain rates intersects the decompression path of the UHP rocks such as the Dora Maira Massif (Fig. 1c). Interestingly, the brittle behavior of garnet can be used to constrain the range of decompression rates required to explain the petrographic observations. It is important to highlight that the term decompression path used here is purely related to the pressure decrease, i.e., independent of any spatial displacement.

Here, we combine mineral-scale thermo-mechanical numerical models with microstructural observations to quantify the geological conditions controlling radial fracture development in garnet at UHP metamorphic conditions in an inclusion-host-matrix system undergoing a large density change. Our results capture the typical patterns of radial microcracks in garnet surrounding retrogressed coesite inclusions. They challenge the traditional interpretation of a late-stage fracturing of the Dora Maria garnets based on the elastic model and show that, provided significant decompression rates, fracturing is induced by the large shear stress generated during volume change in the SiO₂ inclusions at the early stages of the phase transition.

## Results

### The Dora Maira microstructures

We selected a sample from the famous UHP Dora Maira whiteschists, consisting of garnet, phengite, kyanite, talc, quartz/coesite and rutile (Fig. 2a). The SiO₂ in the matrix is quartz and shows a polygonal texture, suggesting that it formed after retrogression from coesite⁹, whereas the SiO₂ inclusions in garnet are either fully turned into quartz or consist of coesite rimmed by quartz. Two types of microcracks are observed in garnet: (1) a set of parallel straight fractures that have the same orientation in all garnets (Fig. 2a) and (2) sets of radial fractures
<figure><img src="image_5.png" /><figcaption>**Fig. 2 |** Microscopic observations, model configuration and reference result.</figcaption></figure>

a Microphotograph of the sample in plane-polarized light showing one garnet with several SiO₂ (quartz) inclusions surrounded by radial cracks. The radial cracks are slightly bent. The matrix around garnet consists of phengite, quartz and rutile.

b Model setup consisting of a garnet with several SiO₂ inclusions in a matrix of SiO₂. A constant dilation rate is applied to the box to simulate the decompression from

the initial conditions of 4.3 GPa and 730 °C¹⁶. Pressure is computed, and temperature depends on the pressure following: $T = (1/(1.02 \times 10^{-2})) \cdot \ln(P/(2.523 \times 10^{6})) + 273.15$. c Colormap showing the logarithm of the plastic strain rate (s⁻¹) after 346 kyr of decompression to 0.5 GPa and 500 °C (dilation rate: 10⁻¹⁴ s⁻¹, mean decompression rate: 11 GPa/Myr). Dilatant shear bands concentrate into radial and curved zones, shown in orange, mimicking the radial fractures around SiO₂ inclusions.

spreading from the SiO₂ inclusions (Fig. 2a). The latter form networks of fractures that either linkup to other SiO₂ inclusions or spread toward the sides of the garnet and eventually connect with the matrix. The fracture density is higher at the inclusion's corners as well as in the linkup regions than on the curved sides of inclusions, as documented and discussed in detail by Whitney et al. ³⁵.

## Thermo-mechanical numerical modeling

A full description of the code used in this study is provided in the “Methods” section, together with a detailed presentation of the model setup (also summarized in Fig. 2b). The use of visco-elasto-plastic rheological models allows for the emergence of dilatant frictional-plastic shear bands in the garnet. These structures share similarities with radial microcracks observed in natural samples and are further qualified as radial fractures. Figure 2c shows the colormap of the plastic strain rate at the conditions where the matrix reaches 1 GPa (i.e., after 3.3 GPa of decompression). The high strain rate focuses on thin localized bands radiating from the SiO₂ inclusions and propagating between inclusions as well as toward the garnet rim. The patterns produced by the shear bands are in good agreement with the observed radial fractures in natural samples: (1) the shear bands are either straight or bent, and (2) the density of the shear bands is higher at the inclusion's corners, as expected based on previous studies³⁵.

The results of the thermo-mechanical numerical models are presented in Fig. 3. During the first stage of decompression, starting in the coesite stability field at 4.3 GPa and 730 °C, the matrix, garnet and inclusions undergo a progressive decrease in *P* and *T* until reaching the conditions of the phase transition (i.e., 2.95 GPa and ~700 °C; Fig. 3a). First, the matrix SiO₂ undergoes the coesite/α-quartz phase transition. At a constant dilation rate, the transition is not immediate and acts as a pressure buffer in the matrix, which remains at the pressure of the phase transition (i.e., 2.95 GPa, Fig. 4a, b). The SiO₂ inclusions remain in the coesite field as the volume increase in the matrix pressurizes the garnet and inclusions. This effect maintains the pressure of inclusions slightly above the level of the phase transition (Fig. 3b). When the entire matrix SiO₂ is turned into quartz, decompression of the garnet and inclusions resumes. At this stage, the decompression paths of the garnet, matrix and inclusions split. The decompression in the matrix progresses normally with a continuous *P-T* decrease along the decompression path, and the coesite in the inclusions can start the phase transition (Fig. 3c). At this point, the geometry of the inclusions as well as the dilation rate control the stress and pressure evolution in

the garnet. The dilation rate is the first-order parameter that controls the decompression rate (Fig. 4c) and that most affects the deformation style and the pressure field in the model. The relative amount of matrix and garnet exerts a minor, second-order influence on the resulting decompression rate (see Supplementary Fig. 1). Hereafter, we present the results based on different dilation rates.

For dilation rates higher than and equal to 10⁻¹⁴ s⁻¹, we achieve a decompression rate of at least 8.0 GPa/Myr (Fig. 4). The garnet behaves in a brittle manner so that the stress increase in the garnet surrounding the transforming inclusions reaches the yield stress as soon as quartz starts forming, which triggers the nucleation of radial fractures (Figs. 2c, 3b). Therefore, the first fractures form early on the decompression path when the matrix reaches a *P* of ca. 2.6 GPa, i.e., soon after the phase transition in the matrix is accomplished (Fig. 3b). The first fractures form linkups between neighboring SiO₂ inclusions by propagating the fracture tips toward the closest inclusions. Then, newly formed radial fractures propagate from the inclusion-garnet boundary toward the garnet sides (Fig. 3c).

If the phase transition of SiO₂ inclusions is disabled (Fig. 4a), i.e., coesite density follows a linear equation of state during decompression, the onset of plastic yielding occurs at lower *P-T* conditions (ca. 1.6 GPa). This confirms the role of the SiO₂ phase transition in building up stresses in the host garnet. The phase transition fosters plastic yielding at high pressure, while differential expansion alone would trigger yielding only at later stages.

For a dilation rate of 10⁻¹⁵ s⁻¹, the decompression rate reaches 1.0 GPa/Myr. Diffuse plastic deformation occurs at a late stage (~1.5 GPa) once all SiO₂ turned to quartz. If the dilation rate is lower than 10⁻¹⁵ s⁻¹, garnet behaves in a ductile manner from the onset of the decompression path until ca. 2.0 GPa. The activation of viscous creep prevents the development of large differential stress around the SiO₂ inclusions so that the formation of quartz is enhanced without fracturing and goes to completion. Overall, localized plastic deformation is observed only for dilation rates greater or equal to 10⁻¹⁴ s⁻¹ and decompression rate of at least 8.0 GPa/Myr (Fig. 4b), and the onset of fracturing occurs consistently between 2.5 and 2.8 GPa.

The model results reproduce the first-order characteristics (i.e., shape, dimensions, location) of the radial fractures observed in natural samples. However, in our simulations, the timing of fracturing and the conditions in which it occurs differ from previous studies. In our models, plastic yielding occurs in response to a stress increase at the interface between the SiO₂ inclusions and the garnet host due to the
<figure><img src="image_6.png" /><figcaption>**Fig. 3 | Time evolution of the numerical model of inclusion-host decompression.** Model results for successive time steps for a decompression rate of 11 GPa/Ma. The first column shows P-T conditions of the matrix (black star) and each SiO₂ inclusion (white diamonds) along the retrograde path (gray dashed line). The black dashed line corresponds to the brittle-ductile transition for garnet at the strain rate of 10⁻¹⁴ s⁻¹, as in Fig. 1c. The second column shows the logarithm of the plastic strain</figcaption></figure>

rate. Blue areas correspond to the parts of the garnet that reached the plastic yield. The third column shows the relative pressure (in GPa) calculated as the pressure difference relative to the matrix. The time steps correspond to: a the SiO₂ inclusions, garnet and matrix SiO₂ reach the coesite-quartz transition; b the onset of garnet fracturing after the complete retrogression of the coesite in the matrix; c the matrix reaches a decompression stage of 0.4 GPa and 500 °C.

phase transition, whereas the elastic model assumes that fracturing is triggered by a pressure difference of a factor 3 between the inclusion and the garnet³³,³⁴. Moreover, in elastic models, fracturing occurs at lower temperatures (<400 °C) so that the kinetics of the coesite-quartz transformation is slow enough to prevent its total completion. In our model, the stress generated by the phase transition produces a pressure variation between the inclusions and garnet host that reaches ca. 2 GPa at the time when fracturing occurs. This happens much earlier and at higher P-T conditions (ca. 2.6 GPa, 680 °C in the matrix) than for elastic models. The complex geometry of inclusions also favors the development of shear stresses in the host in comparison with the purely spherical shapes considered in the previous elastic models. In our case, the preservation of coesite cannot be attributed to sluggish kinetics because the coesite-quartz transformation is almost instantaneous at high T³⁹,⁴⁰, and another mechanism must be invoked to explain the preservation of coesite in the natural samples. To our knowledge, there is no available flow law for coesite. However, based on the observation of highly fractured coesite relicts in the SiO₂ inclusion cores, one could speculate that coesite is more resistant than quartz and hence behaves in a brittle manner, as garnet does, during the decompression. If this is true, the high mechanical strength of

coesite could allow it to remain at high pressure and prevent a total transformation into quartz. Furthermore, elastic models suggest that the pressure vessel effect is released soon after fracturing occurs. In our models, the pressure vessel effect remains after fracturing so that a pressure gradient of up to a factor 2 between the matrix and the inclusions persists toward the late stage of the decompression. The effect of a 3D geometry on the stress distribution around the inclusions remains to be investigated.

## Discussion

The decompression rates based on geochronological data of Dora Maira⁴,⁴¹ between the peak conditions of 4.3 GPa and the retrograde conditions of 1.1 GPa translate into decompression rates varying between 0.8 GPa/Myr (maximum duration of ~4 Ma) and 8 GPa/Myr (minimum duration of ~0.4 Ma). Our results are thus consistent with the fastest decompression rates estimated from petrochronology and independently report the occurrence of ultra-fast decompression rates. Furthermore, the results suggest that radial fractures are formed early after the onset of decompression, i.e., after the phase transition has occurred in the matrix and while the phase transition is still occurring within the inclusions.
<figure><img src="image_7.png" /><figcaption>**Fig. 4 | Systematic investigation of host fracturing during decompression.**</figcaption></figure>

a Pressure (P) evolution of the matrix versus time (t) for a dilation rate of 10⁻¹⁴ s⁻¹. Simulations with phase transition activated either in the inclusions and matrix or only in the inclusions reach the yield stress at higher P conditions than simulations where phase transition is not activated or when phase transition occurs only in the matrix. The lines denote viscoelastic-dominated deformation, and crosses indicate the occurrence of localized frictional deformation. b Pressure evolution of the

matrix versus time for different dilation rates. In contrast with runs with dilation rates higher than 10⁻¹⁴ s⁻¹ (red, light blue and purple lines with crosses), a dilation rate of 10⁻¹⁵ s⁻¹ induces diffuse plastic deformation (green line with circles), and dilation rates lower than 10⁻¹⁵ s⁻¹ caused no frictional plastic deformation (yellow, dark blue lines). c Relationship between dilation rate and the mean decompression rate, the blue dot indicates conditions for diffuse plastic deformation, and the yellow crosses indicate conditions for prominent shear banding.

Based on our results, a decrease in the decompression rate leads to the displacement of the BDT of garnet toward lower P-T conditions on the decompression path (Fig. 1c). As a result, it postpones the onset of fracturing by delaying the moment at which the brittle field is reached. Consequently, a decrease in decompression rate during the decompression would prevent the formation of radial fractures if they are not already formed earlier along the decompression path. Similarly, a slow decompression at such high T would also lead to an efficient viscoelastic relaxation of the P in the inclusions, hence preventing any fracturing from occurring. Therefore, it is essential that the first stage of decompression occurs under fast decompression rates ≥8.0 GPa/Myr in order to predict the petrographic observations. Traditionally, the decompression path is related to the spatial displacement of geological units and exhumation rates (e.g., cm/yr). However, as suggested by Yamato and Brun³⁰, a fast decompression rate can be explained by a sudden change in the state of stress from compression to extension within the lithosphere. The substantial decompression rates of 8 GPa/Myr obtained from the thermo-mechanical modeling support this hypothesis because of the lack of solid evidence for such an enormous spatial (vertical or oblique) displacement⁷.

The ability to model the visco-elasto-plastic rheology of garnet reveals important insights into the link between phase transition and deformation in inclusion-host systems. Here, the density change associated with the SiO₂ phase transition in inclusions of complex geometries has the potential to develop considerable shear stress, leading to fracturing of the garnet upon fast decompression. The results confirm that under the published range of decompression rates for the first stage of the decompression of the UHP Dora Maira Massif, garnet is fracturing at still high P-T conditions of >2 GPa and >600 °C, unlike previously thought based on elastic models. The major implication is that the onset of the decompression must be associated with a rapid switch in the stress state of the buried continental lithosphere, favoring the orogenic wedge model, in which nappe stacking triggers transient stress variations²⁶.

## Methods

### Code description

We use the 2D numerical code MDOODZ³⁸,⁴² based on the finite difference method/marker-in-cell approach, which is suitable for the study deformation of heterogenous visco-elasto-plastic materials. The

latest developments incorporate the effect of compressibility and progressive mineral transformations on rheology and strain localization³⁷. The equation of momentum is formulated as:

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

where $τ_{ij}$ denotes the deviatoric strain rate tensor components, P is the pressure and $x_i$ corresponds to the spatial coordinates. Deviatoric stress is calculated using a Maxwell, visco-elasto-viscoplastic rheological model (e.g., ref. 43):

$$ \dot{\epsilon}_{ij} = \dot{\epsilon}_{ij}^{e} + \dot{\epsilon}_{ij}^{v} + \dot{\epsilon}_{ij}^{vp} \qquad (2) $$

where $\dot{\epsilon}_{ij}$ corresponds to the deviatoric strain rate tensor components, and the e, v and vp superscripts respectively correspond to elastic, viscous and viscoplastic. The deviatoric elastic strain rate tensor is assumed to relate linearly with respect to the objective deviatoric stress rate:

$$ \dot{\epsilon}_{ij}^{e} = \frac{d\tau_{ij}}{dt} \qquad (3) $$

Viscous creep strain rate follows a power law relation:

$$ \dot{\epsilon}_{ij}^{v} = 2^{-n} A e^{-\frac{E+PV}{RT}} \tau_{\parallel}^{n-1} \tau_{ij} \qquad (4) $$

where A is a material parameter, E is the activation energy, V is the activation volume, $τ_{\parallel}$ is the second deviatoric stress invariant and n is the power law index.

We assume that garnet yielding can be modeled using a pressure-dependent and dilatant frictional plastic rheological model. The emergence of shear failure patterns is hence modeled using a visco-plastic Drucker-Prager model. The yield (F) and potential (Q) functions are expressed as:

$$ F = \tau_{\parallel} - C \cos \varphi - P \sin \varphi - \dot{\lambda}\eta^{vp} \qquad (5) $$

$$ Q = \tau_{\parallel} - P \sin \psi \qquad (6) $$
where $\varphi$ is the friction angle, $\psi$ is the dilation angle, C is the cohesion and $\eta^{\text{vp}}$ is the Kelvin element viscosity$^{44}$. This model includes temporal regularization, which allows for successful satisfaction of force equilibrium at each time step. The viscoplastic plastic strain rate is expressed as:

$$ \dot{\varepsilon}_{ij}^{\text{vp}} = \dot{\lambda} \frac{\partial Q}{\partial \tau_{ij}} \qquad (7) $$

where $\dot{\lambda}$ is the rate of the plastic multiplier, which is determined using an analytical return mapping procedure based on trial yield function and material parameters (see details in ref. 42). Frictional shear banding only occurs if the pressure-dependent yield stress is attained. At such temperatures (>600 °C), this limit can be reached in strong hosts like garnet but not in weaker (e.g., omphacite, Supplementary Fig. 2).

The model is compressible and the continuity equation is expressed as:

$$ \frac{d \ln(\rho)}{dt} = - \frac{\partial v_i}{\partial x_i} \qquad (8) $$

The simulations thus account for volumetric strains induced by elasto-plastic deformations and by phase transitions. The density ($\rho$) is function of pressure and temperature. For phases that involve phase transitions (i.e., SiO$_2$), the equilibrium density ($\rho$) is pre-computed using Gibbs free energy minimization and stored in look-up tables. For other materials, the following equation of state is applied:

$$ \rho = \rho_0 e^{\beta P - \alpha T} \qquad (9) $$

where $\rho_0$ is the reference density, $\beta$ is the incompressibility modulus and $\alpha$ is the thermal expansivity. Here, the use of a constant compressibility for garnet is a simplifying assumption$^{45}$. We have tested several values of compressibility ($K = \beta^{-1}$). These variations did not have any significant effect (Supplementary Fig. 3), which is also in line with a recent study by Moulas et al.$^{46}$.

In case of phase transitions, the density evolves as:

$$ \frac{d\rho}{dt} = \frac{\rho^{\text{eq}} - \rho}{t^k} \qquad (10) $$

where $t^k$ is the transformation kinetic time, which is expressed as:

$$ t^k = \ln\left(\frac{1}{3}\right) (-2S\dot{g})^{-1} \qquad (11) $$

where S is the grain boundary area. The growth rate is:

$$ \dot{g} = k_0 T e^{-\frac{Q_g}{RT}} \left(1 - e^{-\frac{\Delta G}{RT}}\right) \qquad (12) $$

where $\Delta G$ is the free energy change of reaction computed using Gibbs free energy minimization, $Q_g$ is the activation energy for growth and $k_0$ is a constant$^{39}$. Note that at the considered temperature, growth kinetic rates are much larger than considered decompression rates. Kinetics has thus a negligible effect as transformation rates are quasi-instantaneous.

## Model configuration

The model configuration is characterized by a 2D medium of 5 × 5 mm, consisting of a matrix SiO$_2$, containing a garnet crystal with several SiO$_2$ inclusions (Fig. 2b). The geometry is based on the digitalized thin section microphotography in Fig. 2a. The density-pressure relationship of SiO$_2$ is calculated based on the Holland and Powell 2011 database$^{47}$

and the density evolution of garnet follows the equation of state. Note that the variation of the molecular volume of $\alpha$-quartz in Holland and Powell 2011$^{47}$ is equal to the one in their 1998 database$^{48}$, which does not significantly differ from other estimates in the literature within the temperature range of interest (730–600 °C)$^{49}$. The kinetics of the coesite to quartz phase transition is taken into account in the models following Mosenfelder and Bohlen$^{39}$ ($Q_g = 269 \pm 26$ kJ/mol, $k_0 = 26.537$ m/s/K). Creep parameters are based on laboratory-derived flow laws: garnet creep follows the flow law of Ji and Martignole$^{50}$, and SiO$_2$ creep was modeled using dry quartz experimental data$^{51,52}$ (see Supplementary Table 1 for material parameters).

The initial conditions are set to the peak Alpine P-T conditions of 4.3 GP and 730 °C, according to Hermann$^{16}$. The decompression of the model from the peak conditions follows the retrograde P-T path constrained by the peak and retrograde P-T data from Hermann$^{16}$ and Rubatto and Hermann$^{4}$, as shown in Fig. 1a. The temperature (unit Kelvin) is assumed constant at the scale of the sample and is updated at each time step using the following regression of the P-T data:

$$ T = 273.15 + \frac{1}{a} \log_{10} \left( \frac{P \times 10^{-9}}{b} \right) \qquad (13) $$

where $a = 1.02 \times 10^{-2}$ and $b = 2.523 \times 10^{-3}$.

Progressive dilation and associated decompression are applied via kinematic boundary conditions. The velocity component normal to each boundary ($v_i$) is updated at each time increment as:

$$ v_i = x_i \times \nabla v^{\text{BG}} \qquad (14) $$

where $x_i$ is the current boundary coordinate and $\nabla v^{\text{BG}}$ is the applied dilation rate, which is assumed constant with time. For materials involving no phase transition, the decompression rate can be directly related to the dilation rate. However, for materials involving a phase transition, the pressure remains buffered until the phase transition is complete. For this reason, the modeled decompression rates vary with time, although the dilation rate is constant (see Fig. 4a). The provided decompression rates are thus averaged over the decompression paths. We model purely isotopic dilation of the sample, and we ensure that boundary conditions do not induce deviatoric strains. To this end, the following non-zero out-of-plane strain rate is applied:

$$ \frac{\partial v_y}{\partial y} = \frac{1}{2} \left( \frac{\partial v_x}{\partial x} + \frac{\partial v_z}{\partial z} \right) \qquad (15) $$

All internally modeled deviatoric strains (e.g., shear bands) are thus solely induced by stress concentrations and caused by internal sample heterogeneity. The development of shear bands could be further enhanced if tectonic deviatoric deformation fields were additionally considered (e.g., pure or simple shear). Our models thus provide conservative estimates on the development of shear failure. In order to further verify our 2D models, we have run 1D host-inclusion-matrix decompression models with radial symmetry (Supplementary Fig. 4). We have tested two types of boundary conditions: (1) constant dilation rate and (2) constant decompression rates. Both types of models show that fast decompression rates (>5 GPa/Myr) or fast dilation rates (>5 × 10$^{-15}$ s$^{-1}$) are needed to reach frictional yielding of garnet, thus confirming the presented 2D results.

## Data availability

The presented results are based on numerical simulations and do not involve source data.

## Code availability

The data generated in this study have been performed using the code MDoodz7.0 [https://zenodo.org/record/8112302]. This code is freely
accessible. All the information needed to install and use it, as well as
any updates, can be found here: https://github.com/tduretz/
MDOODZ7.0.

References

1. Stern, R. J. The evolution of plate tectonics. *Philos. Trans. R. Soc. A*
   *Math. Phys. Eng. Sci. 376*, 20170406 (2018).

2. Korenaga, J. Crustal evolution and mantle dynamics through Earth history. *Philos. Trans. R. Soc. A Math. Phys. Eng. Sci.* **376**, 20170408 (2018).

3. Holder, R. M., Viete, D. R., Brown, M. & Johnson, T. E. Metamorphism and the evolution of plate tectonics. *Nature* **572**, 378–381 (2019).

4. Rubatto, D. & Hermann, J. Exhumation as fast as subduction?
*Geology* **29**, 3–6 (2001).

5. Ague, J. J. *Current and Future Challenges in Metamorphic Petrology*, Vol. 51, No. 5 (GSA Abstracts, 2019).

6. Hacker, B. R. & Gerya, T. V. Paradigms, new and old, for ultrahigh-pressure tectonism. *Tectonophysics* **603**, 79–88 (2013).

7. Schenker, F. L. et al. Current challenges for explaining (ultra)high-pressure tectonism in the Pennine domain of the Central and Western Alps. *J. Metamorph. Geol.* **33**, 869–886 (2015).

8. Warren, C. J. Exhumation of (ultra-)high-pressure terranes: concepts and mechanisms. *Solid Earth* **4**, 75–92 (2013).

9. Chopin, C. Coesite and pure pyrope in high-grade blueschists of the Western Alps: a first record and some consequences. *Contrib. Mineral. Petrol.* **86**, 107–118 (1984).

10. Chopin, C., Henry, C. & Michard, A. Geology and petrology of the coesite-bearing terrain, Dora Maira massif, Western Alps. *Eur. J. Mineral.* **3**, 263–292 (1991).

11. Compagnoni, R., Hirajima, T. & Chopin, C. Ultra-high-pressure metamorphic rocks in the Western Alps. *Ultrahigh Pressure Metamorphism* (eds Coleman, R. G. & Wang, X.) 206–243 (Cambridge University Press, 1995).

12. Ferrando, S., Groppo, C., Frezzotti, M. L., Castelli, D. & Proyer, A. Dissolving dolomite in a stable UHP mineral assemblage: evidence from Cal-Dol marbles of the Dora-Maira Massif (Italian Western Alps). *Am. Mineral.* **102**, 42–60 (2017).

13. Gauthiez-Putallaz, L., Rubatto, D. & Hermann, J. Dating prograde fluid pulses during subduction by in situ U–Pb and oxygen isotope analysis. *Contrib. Mineral. Petrol.* **171**, 1–20 (2016).

14. Gebauer, D., Schertl, H.-P., Brix, M. & Schreyer, W. 35 Ma old ultrahigh-pressure metamorphism and evidence for very rapid exhumation in the Dora Maira Massif, Western Alps. *Lithos* **41**, 5–24 (1997).

15. Groppo, C. et al. What's in the sandwich? New P–T constraints for the (U)HP nappe stack of southern Dora-Maira Massif (Western Alps). *Eur. J. Mineral.* **31**, 665–683 (2019).

16. Hermann, J. Experimental evidence for diamond-facies meta-
morphism in the Dora-Maira massif. *Lithos* **70**, 163–182 (2003).

17. Hirajima, T. & Compagnoni, R. Petrology of a jadeite-quartz/coesite-
almandine-phengite fels with retrograde ferro-nyböite from the
Dora-Maira Massif, Western Alps. *Eur. J. Mineral.* **5**, 943–956 (1993).

18. Kienast, J. R., Lombardo, B., Biino, C. & Pinardon, J. L. Petrology of very-high-pressure eclogitic rocks from the Brossasco-Isasca Complex, Dora-Maira Massif, Italian Western Alps. *J. Metamorph. Geol.* **9**, 19–34 (1991).

19. Massonne, H.-J. & Schreyer, W. Stability field of the high-pressure assemblage talc + phengite and two new phengite barometers. *Eur. J. Mineral.* **1**, 391–410 (1989).

20. Nowlan, E. U., Schertl, H.-P. & Schreyer, W. Garnet-omphacite-phengite thermobarometry of eclogites from the coesite-bearing unit of the southern Dora-Maira Massif, Western Alps. *Lithos* **52**, 197–214 (2000).

21. Schertl, H.-P., Schreyer, W. & Chopin, C. The pyrope-coesite rocks and their country rocks at Parigi, Dora Maira Massif, Western Alps:

detailed petrography, mineral chemistry and PT-path. *Contrib. Mineral. Petrol.* **108**, 1–21 (1991).

22. Sharp, Z. D., Essene, E. J. & Hunziker, J. C. Stable isotope geochemistry and phase equilibria of coesite-bearing whiteschists, Dora Maira Massif, western Alps. *Contrib. Mineral. Petrol.* **114**, 1–12 (1993).

23. Beltrando, M., Compagnoni, R. & Lombardo, B. Ultra-high-pressure metamorphism and orogenesis: an Alpine perspective. *Gondwana Res.* **18**, 147–166 (2010).

24. Butler, J. P., Beaumont, C. & Jamieson, R. A. The Alps 2: controls on crustal subduction and (ultra)high-pressure rock exhumation in Alpine-type orogens. *J. Geophys. Res. Solid Earth* **119**, 5987–6022 (2014).

25. Malusà, M. G. et al. Contrasting styles of (U)HP rock exhumation along the Cenozoic Adria-Europe plate boundary (Western Alps, Calabria, Corsica). *Geochem. Geophys. Geosyst.* **16**, 1786–1824 (2015).

26. Schmalholz, S. M. & Duretz, T. Shear zone and nappe formation by thermal softening, related stress and temperature evolution, and application to the Alps. *J. Metamorph. Geol.* **33**, 887–908 (2015).

27. Schmalholz, S. M. & Schenker, F. L. Exhumation of the Dora Maira ultrahigh-pressure unit by buoyant uprise within a low-viscosity mantle oblique-slip shear zone. *Terra Nova* **28**, 348–355 (2016).

28. 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).

29. Pleuger, J. & Podladchikov, Y. Tectonic overpressure may reconcile the structural and petrological records of the Adula nappe (Central Alps). *EGU Gen. Assem. Conf. Abstr.* **16**, 15622 (2014).

30. Yamato, P. & Brun, J. P. Metamorphic record of catastrophic pressure drops in subduction zones. *Nat. Geosci.* **10**, 46–50 (2016).

31. Gillet, P., Ingrin, J. & Chopin, C. Coesite in subducted continental crust: P-T history deduced from an elastic model. *Earth Planet. Sci. Lett.* **70**, 426–436 (1984).

32. Rosenfeld, J. L. & Chase, A. B. Pressure and temperature of crys-
tallization from elastic effects around solid inclusions in minerals?
*Am. J. Sci.* **259**, 519–541 (1961).

33. Van der Molen, I. The shift of the α-β transition temperature of quartz associated with the thermal expansion of granite at high pressure. *Tectonophysics* **73**, 323–342 (1981).

34. Van der Molen, I. & Van Roermund, H. L. M. The pressure path of solid inclusions in minerals: the retention of coesite inclusions during uplift. *Lithos* **19**, 317–324 (1986).

35. Whitney, D. L., Cooke, M. L. & Du Frane, S. A. Modeling of radial microcracks at corners of inclusions in garnet using fracture mechanics. *J. Geophys. Res. Solid Earth* **105**, 2843–2853 (2000).

36. Zhang, Y. Mechanical and phase equilibria in inclusion–host sys-
tems. *Earth Planet. Sci. Lett.* **157**, 209–222 (1998).

37. Yamato, P., Duretz, T., Bäisset, M. & Luisier, C. Reaction-induced volume change triggers brittle failure at eclogite facies conditions. *Earth Planet. Sci. Lett.* **584**, 117520 (2022).

38. Yamato, P., Duretz, T. & Angiboust, S. Brittle/ductile deformation of eclogites: insights from numerical models. *Geochem. Geophys. Geosyst.* **20**, 3116–3133 (2019).

39. Mosenfelder, J. L. & Bohlen, S. R. Kinetics of the coesite to quartz transformation. *Earth Planet. Sci. Lett.* **153**, 133–147 (1997).

40. Perrillat, J. P., Daniel, I., Lardeaux, J. M. & Cardon, H. Kinetics of the Coesite–Quartz transition: application to the exhumation of ultrahigh-pressure rocks. *J. Petrol.* **44**, 773–788 (2003).

41. Bonnet, G., Chopin, C., Locatelli, M., Kylander-Clark, A. R. C. & Hacker, B. R. Protracted subduction of the European hyper-extended margin revealed by rutile U-Pb geochronology across the Dora-Maria massif (W. Alps). *Tectonics* **41**, e2021TC007170 (2022).
42. Duretz, T., de Borst, R. & Yamato, P. Modeling lithospheric deformation using a compressible visco-elasto-viscoplastic rheology and the effective viscosity approach. *Geochem. Geophys. Geosyst.* **22**, e2021GC009675 (2021).

43. de Borst, R. & Duretz, T. On viscoplastic regularisation of strain-softening rocks and soils. *Int. J. Numer. Anal. Methods Geomech.* **44**, 890–903 (2020).

44. Duretz, T., de Borst, R. & Le Pourhiet, L. Finite thickness of shear bands in frictional viscoplasticity and implications for lithosphere dynamics. *Geochem. Geophys. Geosyst.* **20**, 5598–5616 (2019).

45. Angel, R. J., Gilio, M., Mazzucchelli, M. & Alvaro, M. Garnet EoS: a critical review and synthesis. *Contrib. Mineral. Petrol.* **177**, 54 (2022).

46. Moulas, E., Podladchikov, Y., Zingerman, K., Vershinin, A. & Levin, V. Large-strain elastic and elasto-plastic formulations for host-inclusion systems and their applications in thermobarometry and geodynamics. *Am. J. Sci.* **323**, 2 (2023).

47. Holland, T. J. B. & Powell, R. An improved and extended internally consistent thermodynamic dataset for phases of petrological interest, involving a new equation of state for solids. *J. Metamorph. Geol.* **29**, 333–383 (2011).

48. Holland, T. J. B. & Powell, R. An internally consistent thermodynamic data set for phases of petrological interest. *J. Metamorph. Geol.* **16**, 309–343 (1998).

49. Gerya, T. V., Maresch, W. V., Podlesskii, K. K. & Perchuk, L. L. Semi-empirical Gibbs free energy formulations for minerals and fluids for use in thermodynamic databases of petrological interest. *Phys. Chem. Miner.* **31**, 429–455 (2004).

50. Ji, S. & Martignole, J. Ductility of garnet as an indicator of extremely high temperature deformation. *J. Struct. Geol.* **16**, 985–996 (1994).

51. Ranalli, G. Rheology of the Earth (Chapman & Hall, 1995).

52. Kirby, S. H. & Kronenberg, A. K. Rheology of the lithosphere: selected topics. *Rev. Geophys.* **25**, 1219–1244 (1987).

## Acknowledgements

This work was supported by the Goethe University of Frankfurt, an Alexander von Humboldt research grant to C.L. and a Wilhelm and Else Heraeus Foundation visiting professorship to P.Y. The Heidelberg preparation lab is acknowledged for high-quality thin-section preparation. We thank Taras Gerya for his detailed and constructive reviews.

## Author contributions

All authors participated equally in the design of the study. C.L.: Formal analysis, investigations, numerical simulations & visualization,

writing—original draft & editing; L.T.: Writing original draft & editing and revision; P.Y.: Involvement in the numerical code development, writing—editing and revision; T.D.: Numerical code development, writing—editing and revision, supervision.

## Funding

Open Access funding enabled and organized by Projekt DEAL.

## 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-023-41310-w.

**Correspondence** and requests for materials should be addressed to Cindy Luisier.

**Peer review information** *Nature Communications* thanks Taras Gerya and the other anonymous reviewers for their contribution to the peer review of this work. A peer review file is available.

**Reprints and permissions 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.

**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) 2023