Sergio Chibbaro

Walid Hamouche

Corrado Maurini

Stefano Vidoli

Angela Vincenti

Chaotic and regular dynamics of a morphing shell with a vanishing-stiffness mode

Thin elastic shells are almost inextensible but easy to bend. In the presence of prestresses, geometric frustrations can produce complex elastic energetic landscapes, which have been tailored for the design of morphing structures with multiple stable equilibria or neutrally stable manifolds. We show that the co-existence of stiff and floppy modes leads to unexploited dynamical features. We build a neutrally stable cylindrical shell that under dynamical excitation alternates a chaotic behaviour with a surprisingly regular regimes with a continuous precession of the curvature axis at a constant speed. We explain the experimental findings with a minimal model, showing how the intriguing dynamics is due to the subtle coupling between the prestress, geometrical nonlinearity, material anisotropy and inertial effects. Our results shed a new light on morphing structures dynamics and can be exploited in engineering applications such as energy harvesting.

INTRODUCTION

Slender elastic structures can experience reversible, controllable large shape-changes. The design of morphing structures, i.e. structures capable to adapt their shape to different functions, have received a large interest in the last two decades, based on the observation of natural phenomena [START_REF] Dervaux | Morphogenesis of growing soft tissues[END_REF][START_REF] Forterre | How the venus flytrap snaps[END_REF] or aiming at production of technological devices [START_REF] Cao | Bistable structures for advanced functional systems[END_REF][START_REF] Kebadze | Bistable prestressed shell structures[END_REF]. Nowadays a rich gallery of shapeable multistable objects is available, spanning from origami [START_REF] Waitukaitis | Origami multistability: From single vertices to metasheets[END_REF] to large scale deployable structures [1].

Because of the interplay between elasticity and geometry, prestressed elastic shells provide a fascinating playground for morphing structure design [START_REF] Seffen | Growth and shape control of disks by bending and extension[END_REF][START_REF] Sharon | The mechanics of non-euclidean plates[END_REF]. Their shape can be controlled by prescribing through active materials or plastic prestress the target metric or curvature. They can exhibit multiple stable equilibrium configurations and a complex nonlinear behavior, which can be controlled by plastic deformation [START_REF] Seffen | Prestressed Morphing Bistable and Neutrally Stable Shells[END_REF], growth [START_REF] Chen | Nonlinear geometric effects in mechanical bistable morphing structures[END_REF][START_REF] Goriely | Differential growth and instability in elastic shells[END_REF][START_REF] Pezzulla | Geometry and mechanics of thin growing bilayers[END_REF], thermal effects [START_REF] Brunetti | A class of morphing shell structures satisfying clamped boundary conditions[END_REF] or other multiphysical couplings [START_REF] Hamouche | Multi-parameter actuation of a neutrally stable shell: a flexible gear-less motor[END_REF][START_REF] Na | Grayscale gel lithography for programmed buckling of non-euclidean hydrogel plates[END_REF]. Yet, only few works have investigated their nonlinear dynamic behavior [START_REF] Arrieta | A piezoelectric bistable plate for nonlinear broadband energy harvesting[END_REF][START_REF] Gomez | Critical slowing down in purely elastic 'snap-through' instabilities[END_REF]. Bistable structures have interesting applications for the design of energy harvesters in realistic vibration environments, due to their ability to possible self-synchronisations over a broad frequency range [START_REF] Harne | A review of the recent research on vibration energy harvesting via bistable systems[END_REF][START_REF] Stanton | Nonlinear dynamics for broadband energy harvesting: Investigation of a bistable piezoelectric inertial generator[END_REF]. It is has been shown experimentally that when the energy barrier between the energy wells is low, the potential gain for power harvesting ranges between 400% and 600% with respect to standard linear oscillators [START_REF] Cottone | Nonlinear energy harvesting[END_REF].

In this Letter, we show how the presence of soft bending modes in almost inextensible shells leads to peculiar dynamical properties. We consider thin multistable shells that are obtained by applying a large isotropic plastic curvature to an initially flat disk [START_REF] Guest | A zerostiffness elastic shell structure[END_REF]. Their equilibrium shape is almost cylindrical, and because of the thinness, the shell is almost inextensible. This kind of structures have interesting properties for shape-control [START_REF] Hamouche | Multi-parameter actuation of a neutrally stable shell: a flexible gear-less motor[END_REF] and fluid-structure [START_REF] Corsi | A neutrally stable shell in a stokes flow: a rotational taylor's sheet[END_REF] applications. Here, we demon-FIG. 1. Top: the initially flat copper disk and the two stable equilibrium shapes after the plastic deformations. Center: the current cylindrical shape is described by the curvature amplitude c and direction ϕ; the plot shows the elastic energy landscape in polar coordinates c and 2ϕ, with the two energy minima m1,2 (c = 1, ϕ = 0, π/2) and the two saddle points s1,2 (c = (1 + ε) -1 , ϕ = ±π/4). Bottom: the energy barrier between the two minima along the circle c 1, proportional to ε, a small dimensionless parameter related to material and geometrical imperfections.

strate that the interplay of the pre-stress, geometrical nonlinearities, small material anisotropies and inertial effects results into a complex dynamics alternating chaotic and regular responses. In the regular regime, the shell converts a translational excitation in a continuous precession mode characterized by a constant-speed rotation of the curvature axis, see the video in the Supplemental Material [2]. We develop a minimal two-degree-of-freedom model relating the presence of the soft precession mode to geometrical nonlinearities and the prestress, and its small stiffness to the material anisotropy. The nontrivial dynamics results from the nonlinear inertial coupling between the transverse oscillations and the soft bistable mode associated to the precessions of the curvature axis. That may open up the way to new applications of morphing structures in dynamics, as for the design of energy harvesting devices [START_REF] Arrieta | A piezoelectric bistable plate for nonlinear broadband energy harvesting[END_REF][START_REF] Cottone | Nonlinear energy harvesting[END_REF].

EXPERIMENTAL PROTOTYPE

The experimental prototype of vanishing-stiffness shell is fabricated by winding an initially flat disk of radius R = 122.5 mm and thickness h = 0.3 mm over PVC cylinders of progressively smaller diameters, up to a minimal radius of curvature 110 mm. This process induces an isotropic plastic curvature [START_REF] Hamouche | Multi-parameter actuation of a neutrally stable shell: a flexible gear-less motor[END_REF]. The disk is made of nominally isotropic copper (Cua1H14) with Young modulus Y = 124GPa, Poisson ratio ν = 0.33, mass density ρ = 9690Kg m -3 . After several winding cycles along many radial directions, the shell presents two stable cylindrical configurations of almost uniform curvature, marked by the points m 1 and m 2 in Fig. 1. The orientations of their curvature axes are mutually orthogonal, whilst their radii of curvature are approximately equal to R 0 = 178 ± 2mm. The direct manipulation of the shell reveals a stiff bending mode and a soft deformation mode associated with the rotation of the curvature axis. In the following, we will refer to the soft mode as zero-stiffness or neutrally stable modes, see also [11, 17, 18, 26? ? ] To study the system dynamics, we have clamped the center of the shell on a shaker (Bruel&Kjaer 4809) and imposed an harmonic seismic oscillation of the mounting point in the form w(τ ) = W sin 2πf τ . To measure the time-histories of curvature amplitude c(τ ) and direction ϕ(τ ), we have recorded the shell motion with a camera giving a top view of the shell. For every frame, we used an image recognition tool to fit the shell boundary to an ellipse (dashed green in the inset of Fig. 2). The direction of the minor axis of the fitted ellipse gives the value of ϕ at that instant. The curvature is inferred from the ratio of the axis lengths. For small values of the amplitude W , the shell experienced small oscillations near the stable equilibria. For both m 1 or m 2 we have measured a natural frequency f 0 = 18 ± 0.3Hz and a damping factor ζ 0.03. When increasing the forcing amplitude W the rotation of the curvature axis starts to be visible. One progressively observes transitions between the two stable equilibria, eventually associated with the continuous precession of the curvature direction ϕ or with an intermittent dynamics, see the experimental videos in [2]. The diagrams in Fig. 2 show two typical time-histories for the curvature direction ϕ. In the orange curve (f = 15.5 Hz) the transition between the stable equilibria m 1 (ϕ = 0 ≡ π) and m 2 (ϕ = π/2) is chaotic. The dark red curve (f = 16.35 Hz) corresponds to a regular precession of the curvature direction with at an almost constant speed dϕ(τ )/dτ of about π/(0.122 s) 25.7 rad/s.

UNIFORM CURVATURE SHELL MODEL

To develop an analytical model, we consider the object as an inextensible Koiter shell with curvature uniform in space [START_REF] Ciarlet | An introduction to differential geometry with applications to elasticity[END_REF]. Let (O, e 1 , e 2 , e 3 ) be an orthonormal fixed reference frame (see Figure 1). Hence, the curvature is described by a rank-one 2 × 2 symmetric tensor k(c, ϕ) = (c/R 0 ) e S (ϕ) ⊗ e S (ϕ) with e S (ϕ) = cos ϕ e 1 + sin ϕ e 2 , where c is the dimensionless curvature amplitude, ϕ the orientation of the curvature axis, and ⊗ denotes the usual tensor product. After the initial plastic forming process, the material behavior is linearly elastic, but slightly anisotropic. Denoting by 1 the identity tensor, the target plastic curvature is k p = 1/R p . The elastic misfit is k el = kk p and the elastic energy density per unit of surface is

D 2 (1 -ν) k el 2 + ν Tr(k el ) 2 + ε (k el e 1 • e 2 ) 2
where k 2 = k ij k ji and D := Y h 3 /(12(1 -ν 2 )) is the bending modulus and ε introduces the effect of a small material anisotropy, that we attribute to the anisotropic plastic hardening. Integrating over the shell surface, the elastic potential energy reads where

V(c, ϕ) = V 0 c 2 2 1 + ε 2 (1 -cos 4ϕ) -c + const, (1)
V 0 = πh 3 Y R 2 12R 2 0
. In Fig. 1 the shell elastic energy is plotted in the plane where c is the radial coordinate and 2ϕ the angular coordinate. The energy barrier between the saddle points and the minima is

V 1 1 + ε , π 4 -V(1, 0) = ε 2(1 + ε) V 0 . (2)
In the ideal case of a perfectly isotropic shell (ε = 0), the elastic energy is independent of the orientation ϕ of the axis of curvature, see the inset in Fig. 1. This justifies the appellation "zero-stiffness" given in [START_REF] Guest | A zerostiffness elastic shell structure[END_REF] as the shell equilibrium is indifferent to variations of ϕ. In the physical object, the imperfections lead to positive values of the parameter ε and to a bistable behavior. Relating the value of the imperfection to the margin of stability of the two stable equilibria we have measured ε 0.032 ± 0.005 [START_REF] Hamouche | Multi-parameter actuation of a neutrally stable shell: a flexible gear-less motor[END_REF]. Under the same assumptions used to derive the elastic energy, the leading term of the kinetic energy of the shell is found to be T (c, ċ, φ) = ( ċ2 2 + 2 3 c 2 φ2 )V 0 where the superimposed dot denotes the derivative with respect to a dimensionless time t := Ω 0 τ , with Ω 0 = h 20Y /ρ/(3R 2 ). The dimensionless Lagrangian of our two degree-of-freedom model reads as

L = ċ2 2 + 2c 2 φ2 3 - c 2 2 1 + ε(1 -cos 4ϕ) 2 + c. (3)
The structural and acoustic damping is modeled by a Rayleigh dissipation function R = δ( ċ2 /2 + 2c 2 φ2 /3), where δ is a dimensionless damping coefficient. Hence, modeling the effect of the shaker as an harmonic seismic forcing of dimensionless amplitude A = 40R0W 9R 2 and frequency ω = Ω/Ω 0 , the equations of motions read as:

c + δ ċ + c 1 - ε(1 -cos 4ϕ) 2 - 4c φ2 3 = 1 + Aω 2 sin ωt, d dt 4c 2 φ 3 + δ 4c 2 φ 3 + ε c 2 sin 4ϕ = 0.
(4) We refer the reader to the Supplemental Material, section I, for a detailed derivation of the model equations.

FREE UNDAMPED DYNAMICS

We start looking at the undamped free response (δ = 0, A = 0) of the associated Hamiltonian system, for which the total energy H = T + V is conserved and constitutes a first integral. In the perfect case (ε = 0), the angular momentum p ϕ = 4c 2 φ/3 is also conserved, and the system reduces to a one-dimensional oscillator in the c variable under the action of the effective potential

V eff (c) = c 2 /2 -c + 3p 2 ϕ /8c 2 .
The system is integrable, its dynamics is diffeomorphic to a 2-dimensional torus and chaos is not possible [START_REF] Arnold | Mathematical methods of classical mechanics[END_REF][START_REF] Cencini | Chaos: from simple models to complex systems[END_REF]. Motions with precession of the curvature axis at constant speed v and curvature amplitude c v are possible for initial conditions φ(0

) = v, ċ(0) = 0, c(0) = c v = 1 1-4v 2 /3
, where c v is the minimum of the effective potential V eff . Small oscillations around this equilibrium are in the form

c(t) = c v + a sin ω v t, ϕ(t) = vt + 2av ω v c v (1 -cos ω v t) (5)
with a 1 and

ω v = V eff (c v ) = √ 1 + 4v 2 . Since ϕ is a π 2 -periodic angular variable, see (3)
, a countable family of periodic solutions (c, ϕ) n/k , with integer n and k, are possible for precession speeds

v n/k = k 4n ω v = ±1/(2 4n 2 /k 2 -1), see the Suppemental Material [2],
Section 2. This kind of periodic solutions are key to understand the free and forced response of the imperfect system.

To investigate the imperfect system, we numerically integrate the governing equations ( 4) with ε = 0.032, varying the initial condition φ(0) = v. For each v, we calculate the fractal dimension [START_REF] Cencini | Chaos: from simple models to complex systems[END_REF] of the Poincaré sections in the hyperplane ϕ = π/2, see Fig. 3. The details on the numerical methods used are reported in the supplemental material [2]. The dimension of the attractor tells that a chaotic dynamics is possible. The onset of global chaos is analyzed through the maximal Lyapunov exponent, λ max , which proves a possible chaotic behavior for sufficiently large imperfections (ε 0.01), and notably for the experimental value used. Yet, some stable quasi periodic motions are experienced, as shown by the points with vanishing fractal dimension. As shown in the insets in Fig. 3, consistently with the KAM theory [START_REF] Arnold | Mathematical methods of classical mechanics[END_REF], the lowest-order periodic orbits, (c, ϕ) n/k in (5) obtained for ε = 0, persist almost unchanged in the imperfect system [2, Fig. 4].

FORCED DYNAMICS

The free undamped response shows that the almost zero-stiffness deformation mode implies the presence of periodic solutions with continuous precessions of the curvature axis. In the damped physical system, the dissipation is not negligible and an external excitation A > 0 is necessary to observe them experimentally. Figure 4 reports the numerical results for the forced damped evolution, using the experimentally identified value of the damping factor δ 3%. For each value of the forcing amplitude A and frequency ω, we integrate numerically (4) for a sufficiently long time interval t max ∼ 10 3 and plot the ratio between the average precession speed φ and the forcing frequency ω, neglecting the initial transient response [2]. Around the basic bending frequency f 0 = √ 5h/(3πR 2 ) Y /ρ, the contour plot shows the existence of a large yellow region, where φ = ω/4. This regime corresponds to the periodic motions with a regular precession of the curvature axis. The numerical results are consistent with the experimental points E 1-6 in the plane (ω,A) for which such continuous precession regimes were observed. These periodic solutions are as the free periodic motion (5), but with the oscillating frequency ω v replaced by the forcing frequency ω. They have a well-defined precession speed v 2/2 = ω/4, corresponding to a locking condition with n = 2 c-oscillations for a full precession of the curvature axis, k = 2. Higher harmonics are in principle possible but not visible experimentally. The central region with the regular preces- sion regime is surrounded by chaotic dynamics where the curvature axis intermittently oscillates between the two energy wells as displayed by numerical and experimental trajectories, which are in a reasonable agreement in all cases. For small forcing amplitudes (case N 1 ), the result of the forced motion are small oscillations in c without rotation of the curvature axis (ϕ 0). The cubic nonlinearity renders these solutions unstable for sufficiently large amplitudes. We have performed the stability analysis of the periodic solution with ϕ = 0 , which provides the marginal curve A s (ω) in the contour plot of Fig. 4 [see 2]. For A > A s (ω), finite amplitude oscillations in ϕ are found. Depending on the available energy, they can trigger a transition from intra-well oscillations around ϕ = 0, as shown by the trajectories N 2 , to cross-well motions, including the solutions with continuous precession of the curvature axis, as E 5 ≡ N 5 and intermittent chaotic oscillations, as N 4 . A particularly intriguing feature of the nonlinear dynamics is that the transversal seismic forcing in the c variable results in regular precessions of the angular variable ϕ. The bottom of Fig. 4 shows that the precession speed, in the yellow region, is locked to the value ω/4 independent of the forcing amplitude, the imperfection value and the damping factor [2, Fig. 6]. The insets of Figure 4 shows that the experimental trajectories in red are asymetric, because of the imperfections present in the tested shell. This symmetry breaking determines the direction of the precession in the experiments [see the video in 2].

CONCLUSIONS

We have shown that prestressed cylindrical shells with floppy deformation modes exhibit an intriguing dynamics including regular the precession or the intermittent oscillations of their curvature axis. Our minimal two degrees-of-freedom model traces back the observed behavior to the nonlinear inertial coupling between the soft bistable mode associated to the rotation of the curvature axis and the stiff mode associated to the variations of the curvature amplitude. The presence of the soft mode is due to the competition between the prestress, that promotes a spherical curvature, and the geometrical nonlinear effects, that constrains the thin elastic shell to remain almost cylindrical, without assigning any preferential orientation of the curvature axis. The non-vanishing stiffness of the soft mode is related to a small material anisotropy. Breaking the rotational symmetry, this nonvanishing stiffness produces the bistable behavior and it is essential to explain the chaotic dynamics. We have characterized the regular regimes as persistent periodic solutions of the perfect system, showing that their average precession speed is precisely determined by a locking condition with the transverse oscillation forcing frequency f , being constant for a large span of amplitudes:

dϕ dτ = 2πf 4 = π 2 f. (6)
The reported behavior is independent of the length-scale and can easily be replicated. That opens up the way to the design of energy harvesting devices at microscopic and macroscopic scales exploiting almost neutrally stable manifolds. The presence of floppy modes that allows for controllable shape-change is the characteristic feature of morphing structures, as pre-stressed annular rings [START_REF] Moulton | Morphoelastic rods. part i: A single growing elastic rod[END_REF]25] or degenerate meta-materials [START_REF] Gary | Programming shape using kirigami tessellations[END_REF][START_REF] Noselli | Smart helical structures inspired by the pellicle of euglenids[END_REF]. This gives them specific properties, as the ability to support low-energy fractional excitations [START_REF] Sun | Fractional excitations in noneuclidean elastic plates[END_REF]. Our results constitute a first step in unveiling the influence of soft modes in the structural dynamics.

Our minimal model can constitute also a novel case-study for the illustration of surprising nonlinear effects in Hamiltonian mechanics.

vanishes for vanishing imperfections. In the ideal case ε = 0 of a perfectly isotropic shell, the elastic energy is flat along the circle c = 1 + ν. This justifies the appellation "zerostiffness" given in [? ] as the shell equilibrium is indifferent to variations of ϕ. In reality imperfections, in the material or in the imposed plastic deformations, leads to positive values of the parameter ε: in our shells we have measured ε 0.03 ± 0.005. Despite the fact that, when ε > 0, the curvature in the saddle is slightly different from the curvature in the minima, one can easily check that the maximal energy gap along the circle c = 1 + ν still scales as O(ε). Hence, in any configuration the stiffness to variations of the curvature direction is

∂ ϕϕ V = O(ε), but
w = h -1 p (1 -cos(cs))/c, v = v e S (ϕ), v = h -1 p (sin(cs)/c -s). (6)
We derive these expressions recalling that c but also s by means of ϕ are function of time, to get velocity fields ẇ and v. Hence we compute the shell kinetic energy by

T = 1 τ 2 c Ω ρ h ( ẇ2 + v • v) dΩ ( 7
)
where ρ is the density per unit volume and the dots indicate derivatives with respect to the dimensionless time scale t = τ /τ c . To obtain a closed form result we use power series of [START_REF] Cao | Bistable structures for advanced functional systems[END_REF] in term of c to get:

T = 3π R 6 h ρ 80τ 2 c R 2 0 1 + 50R 2 189R 2 0 c 2 + O(c 3 ) ċ2 2 + 1 + 5R 2 36R 2 0 c 2 + O(c 3 ) 2 3 c 2 φ2 . ( 8
)
Since for our shell R 0.1225 m and R 0 0.175 m, we note that the two dimensionless terms of order c 2 are both small with respect to 1 even if c = O(1). Therefore, we neglect the higher order terms in c (this assumption corresponding to the adoption of a shallowshell kinematics) and suitably rescale the kinetic energy to obtain the expression used in the paper T (c; ċ, φ) = ċ2 /2 + (2/3) c 2 φ2 . When the shell is forced by the shaker the kinetic energy is computed similarly by adding the shaker imposed motion, w M (τ ) = W sin(2πf τ ), to the transverse displacement (6) 1 . The elastic and kinetic energy can be rescaled by the same term choosing a suitable time scale, namely τ c = 3R 2 /h ρ/(20Y ). With the Young modulus and density of copper and a thickness of about h 0.3 mm we have τ c 0.01 s.

The Lagrangian in dimensionless form is therefore:

L(c, ϕ; ċ, φ) = T -V = ċ2 2 + 2 3 c 2 φ2 - c 2 2 1 + ε 1 -cos 4ϕ 2 + c (1 + ν) + A α cos(αt) ċ. ( 9
)
where

A := 40W R 0 9τ 2 c R 2
is the dimensionless amplitude and α := 2πf τ c the dimensionless frequency. The equations of motion are easily computed to get

c + c 1 - ε 2 + ε 2 cos(4ϕ) -1 -ν - 4c φ2 3 = Aα 2 sin(αt), d dt
4c 2 φ 3 = -εc 2 sin(4ϕ). ( 10
)
Using Hamilton formalism, we define the generalized coordinates q = {c, ϕ} and the conjugate momenta p = {p c , p ϕ } = { ċ, 4c 2 φ/3}. The free dynamics of our shell is governed by the Hamiltonian

H(c, ϕ; p c , p ϕ ) = p 2 c 2 + c 2 2 -c (1 + ν) + 3p 2 ϕ 8c 2 + ε c 2 4 (1 -cos 4ϕ), (11)
or by the first-order ODEs

ċ = p c , φ = 3p ϕ 4c 2 , ṗc = 1 + ν -c 1 - ε 2 + ε 2 cos(4ϕ) + 3p 2 ϕ 4c 3 , ṗϕ = -εc 2 sin(4ϕ). ( 12
)
A. The integrable case

In absence of forcing A = 0 and for ε = 0, the shell has no preferred direction of curvature and the system is integrable by quadratures. Indeed, the Hamiltionan does not depend on ϕ and, therefore, the momentum p ϕ = 4c 2 φ/3 is conserved. Moreover, the energy

E = ċ2 /2 + c 2 /2 -c (1 + ν) + 3p 2 ϕ /(8c 2 ) is conserved.
We have two first integrals and, by the Arnold theorem [? ], the motion is diffeomorphic to a 2-dimensional torus.

To solve the problem, we eliminate φ = 3p ϕ /(4c 2 ) to obtain

ċ2 = 2(E -U (p ϕ , c)), with U (p ϕ , c) := 3p 2 ϕ 8c 2 + c 2 2 -c (1 + ν). (13)
or equivalently the following equation for the curvature: Given the constant value of p ϕ this is a nonlinear ODE in c(t). It is interesting to remark that, if initially ċ(0) = 0, one possible solution is a uniform precession at constant curvature, namely c(t) = (1 + ν)/(1 -4/3v 2 0 ) and ϕ(t) = v 0 t. Here the effective potential well has been shifted from c 1 + ν due to the angular momentum (a sort of centrifugal effect).

c(t) + c(t) -(1 + ν) - 3 p 2 ϕ (t) 4 c(t) 3 = 0. (14)
We also observe that, if p ϕ = 0, then c must be strictly positive for the energy to be finite. The relation [START_REF] Dervaux | Morphogenesis of growing soft tissues[END_REF] between c and ċ is represented in Fig. 1 on the left. It allows to compute the relation between t and c by a simple quadrature

t = c c 0 1 2(E -U (p ϕ , x)) dx, (15)
The maximal and minimal values of c during these nonlinear oscillations can be computed solving for ċ = 0 into [START_REF] Dervaux | Morphogenesis of growing soft tissues[END_REF]. Figure 1 on the right displays such values. In particular the limit values of the curvature can be well approximated solving the problem with p ϕ = 0 obtaining the bounds 0 ≤ c ≤ (1 + ν) 2E + (ν 2 -1), see Fig. 1 on the right.

II. NUMERICAL METHODS

For ε > 0, the shell is bistable with preferred direction of curvature at ϕ = 0 and ϕ = π/2;

the system is no longer integrable and numerical experiments suggest an intricate dynamics between regular and irregular patterns. For ε > 0, the momentum p ϕ is not conserved but, as the system is autonomous, the total energy is still conserved, the perturbed system being quasi-integrable with two degrees of freedom, the KAM theorem assures the persistence of the integrable trajectories and, the KAM tori divides the space-phase, so that for a small enough perturbation the system should remain in a deformed integrable state.

Numerical evidence.

We have comprehensively analysed the system numerically solving the Hamilton equations corresponding to [START_REF] Corsi | A neutrally stable shell in a stokes flow: a rotational taylor's sheet[END_REF].

We have tested the free system response when the shell is initially placed at rest in the minimum m 1 and is then perturbed by a precession velocity φ(0) = v 0 . In this case, the initial total energy is E 0 = 2(1 + ν) 2 v 2 0 /3; equating this value to the elastic energy in the saddle we get v * 0 = √ 3ε/2 + O(ε 3/2 ) as the minimal initial precession speed for the system to explore the whole state space.

Starting from vanishing initial angular speed v 0 , we numerically observe small oscillations around m 1 . These motions are superpositions of oscillations with dimensionless period T = 2π in the c coordinate and oscillations with dimensionless period √ 3 ε T in the ϕ coordinate. The actual physical period for our shell is T = 2π/f 0 0.35 seconds, with f 0 = 18 Hertz the measured frequency of the small oscillations in c.

As the initial angular speed v 0 increases, these motions show an increasingly nonlinear character but the angular variable remains bounded -π/2 < 2ϕ(t) < π/2; refer to Fig. ??

to check that the other stable equilibrium position m 2 is never reached. The minimal initial energy to break this bound is related to the elastic energy in the saddle, see Eq. ( 5), and, therefore scales as O(ε).

Poincarè sections evidence

For every initial condition v 0 , we have run a numerical simulation in the interval t ∈ Already when the initial precession speed v 0 0.158, slightly above the critival value v * 0 0.15, we observe a region of chaotic motion followed by a small interval, around v 0 0.169, where the motion is quasi-periodic. From v 0 ∼ 0.2 to v 0 ∼ 0.35 the character of the motions is also quasi-periodic. A large region of chaotic motion is observed from v 0 ∼ 0.35 to v 0 ∼ 0.43; then again at v 0 0.435 the motion is almost periodic.

The three quasi-periodic motions, displayed in the insets of Fig.

?? at v (1) 0 = 0.435, v (2)
0 = 0.232, v (3
) 0 = 0.169, corresponds to local minima of the Poincarè sections volume. These are characterized by precessions at almost constant speeds φ(t) ∼ v(i) and by functional relations, between curvature and curvature direction, in the form c(ϕ) 1 + ν -h (i) cos(2 n (i) ϕ) with h (i) three real values and n (i) three integers for i = 1, 2, 3.

For instance, when v 0 = v

(t) + 1 + ε 2 - 4v 2 3 - ε 2 cos(4tv) c(t) = (1 + ν), (16)
or, changing the time scale t → τ /(4v), the Mathieu equation

c (τ ) + 2 (1 -η cos τ ) c(τ ) = c, (17)
where c := (1 + ν)/(16v 2 ), η := 3ε/(3ε -8v 2 + 6), 2 := (3ε -8v 2 + 6)/(96v 2 ) and a prime indicates differentation with respect to τ . The resonance condition of ( 17), namely = k/2 for k = 1, 2, 3, ..., leads to the resonant precession speeds 2% viscous dissipation. Results are displayed in Fig. 4: for small amplitudes, the system is trapped in a region around the inital position and cannot move freely. For an energy sufficient to overcome potential barrier, the resonance is activated and the average speed turns out to be around the Mathieu k = 2 resonance over a range of amplitudes. For large amplitudes, such that the energy landscape is basically irrelevant, the dynamics is found to be regular with average velocity given by Mathieu k = 1 resonance. Between the two resonances, there is a region of chaotic dynamics which depends on the imperfections.

v * k=1,2,3... = √ 6 + 3ε 2 √ 2 √ 1 + 3k 2 = √ 3 2 √ 1 + 3k 2 + O(ε). (18)
Remarkably, at experimental conditions it is found that the plate operate almost always with no-stifness dynamics.

FIG. 2 .

 2 FIG. 2. Two recorded histories in the forced response in terms of curvature direction ϕ(τ ) obtained with amplitude W = 2.3 mm but different forcing frequencies f . The orange curve displays an irregular pattern with quasi-periodic jumps for f 0.86f0. The dark red curve indicates a regular periodic motion with almost constant precession speed for f 0.92f0. The inset displays a shell image as recorded by the high-frequency camera lying on the vertical of the shaker: in every frame the shell boundary is fitted to an ellipse (green) to estimate the current values of c and ϕ. Videos are available in [2].

FIG. 3 .

 3 FIG. 3. Free undamped response of Eq.(4) with ε = 0.032: fractal dimension of the Poincaré sections as function of the initial precession speed (blue points), for initial conditions ϕ(0) = ċ(0) = 0, c(0) = 1, φ(0) = v. For the seven representative values of v (thin vertical lines), the insets display the trajectories in polar coordinates (c, 2ϕ). Periodic motions are marked with their ratio n/k denoting the number n of oscillations in c for passing from ϕ = 0 to ϕ = kπ/2. The gray region on the left distinguishes the minimal initial speed v * √ 3ε/2 0.15 to reach the saddle point s and to escape from the initial energy well. The Lyapunov exponent in the top inset indicates a chaotic behavior for ε > 0.008.

25 FIG. 4 .

 254 FIG.4. Top-left: Contour plot for the ratio between the asymptotic precession speed φ and ω as a function of the forcing amplitude A and frequency ω found with direct integration of Eqns. (4) for ε = 0.032 and δ = 3%. The yellow region with φ = ω/4 corresponds to the regular precession of the curvature axis and is consistent with the experimental observation of regular motions (points E1-6). Top-right: Trajectories in the polar coordinates (c, 2ϕ) at different frequencies, corresponding to the points N1-4: numerical (black) and experimental (dark-red); the corresponding time histories are reported in the top inset. Bottom: Cross section of the contour-plot above for ω = 0.92 and for three different imperfection values.

  the stiffness to variation of the curvature is ∂ cc V = O(1); the actual shell is bistable with the energy barriers O(ε) or O(1) depending on the path connecting the wells. The kinetic energy is computed once known the velocity field in terms of Lagrangian parameters. Let x = x 1 e 1 + x 2 e 2 the generic position in the shell reference domain Ω and let s(τ ) := x/R 0 • e S (ϕ(τ )) a dimensionless abscissa in the curvature direction at time τ . Referring to Fig. ??, it easily seen that the transverse w and inplane v displacement fields respectively are

FIG. 1 :

 1 FIG. 1: The integrable case. Left) Eqn.(13) in the plane ( ċ, c) for several values of the initial energy E and F 0.05. The vertical line labels the value (1 + ν). Right) The maximal and minimal values of c as functions of the initial energy E. The bounds obtained for p ϕ → 0 remain valid when p ϕ is increased.

[0, 10 3 T

 3 FIG. 2: Relevant Poincarè sections in the three-dimensional space (p ϕ , c, p c ). These reduce to single points when φ(0) = 0.169, 0.232, 0.435, respectively the orange, red and brown points. They reduce to two black (resp. three black) points when φ(0) = 0.3 and φ(0) = 0.349, to a cloud of green points φ(0) = 0.39 and to the green curve when φ(0) = 0.5.

Figure

  Figure ?? plots, on a logarithmic axis, the volume of the smallest right parallelepiped containing S v 0 . As Poincarè sections densely fill regions of the phase-space for chaotic motions and shrink to points when the motions are periodic, Fig. ?? allows for an overview of the system response for a wide range of initial precession speeds; therein v 0 ranges from 0.1 to 0.5 with steps of 2 * 10 -5 . The insets display the time-histories of c(t) (radial coordinate) and 2ϕ(t) (angular coordinate) depicted as black curves over the elastic energy contours of Fig. ??.

  average precession speed is about v(1) 0.326 and n (1) = 2 as the curvature is minimal twice (in m 1 and m 2 ) and maximal twice (in the saddles). Instead, at v 0 = 0.232 we have v(2) 0.179 and n (2) = 3, whilst at v 0 = 0.169 one observes v(3) 0.104 and n (3) = 4. It seems that the lower is the initial precession speed, depends on ϕ within an error o(ε) and we are left with c

Figure 3 FIG. 3 : 4 FIG. 4 :

 3344 Figure 3 displays the relation between the Mathieu pulsation and the precession speed v, comparing the critical values (18) with the numerical ones found in Fig. ??.

Sergio Chibbaro Sorbonne Université, CNRS, UMR 7190, Institut Jean Le Rond d'Alembert, F-75005 Paris, France Corrado Maurini Sorbonne Université, CNRS, UMR 7190, Institut Jean Le Rond d'Alembert, F-75005 Paris, France Stefano Vidoli Dipartimento di Ingegneria Strutturale e Geotecnica, Sapienza Università di Roma, Italy Angela Vincenti Sorbonne Université, CNRS, UMR 7190, Institut Jean Le Rond d'Alembert, F-75005 Paris, France * * Your e-mail address; Your web page

I. DERIVATION OF THE MODEL

To the derive the shell elastic energy we consider the shell as inextensible. For sufficiently thin shells this hypothesis is well justified, see [? ]. Starting from a flat natural configuration and recalling classical results in the theory of shells, see e.g. [? ], the inextensibility hypothesis implies that the Gaussian curvature is everywhere constant and vanishing. Hence the curvature field k is a field of rank-one 2 × 2 symmetric matrices, namely k = k(c, ϕ) = c e S (ϕ) ⊗ e S (ϕ) for some real c and some angle ϕ. Here e S = (cos ϕ, sin ϕ) means the curvature direction, see Fig. ??. Justified by the experimental evidence, we further assume that the curvature c and the curvature direction ϕ do depend on time but are spatially uniform; thus, at every instant τ , the shell configuration is cylindrical and described by the two Lagrangian parameters c(τ ) ∈ [0, ∞) and ϕ(τ ) ∈ [0, 2π). The shell being inextensible, the elastic energy V is reduced to the bending contribution only; assuming the shell as linearly elastic we have

with k = kp 1 the isotropic plastic deformation imposed by the winding process and D is the bending stiffness tensor

Here D := Y h 3 /(12(1 -ν 2 )) is the bending stiffness, Y the Young modulus, h the shell thickness, 1 the 2 × 2 identity tensor and Tr the trace operator. Due to ε parameter, the material is orthotropic with Poisson ratio ν: indeed, the material shear stiffness is

where we have introduced the dimensionless curvature c := c/R 0 . Once dropped the constant term and rescaled, (3) reads to

in the paper. The elastic energy V has two wells

and two saddle points s : (ϕ = ±π/4, c = (1 + ν)/(1 + ε). The energy gap between the saddle points and the minima

the more are the curvature oscillations in every complete precession of the curvature axis.

Remark also that in all the cases v(i) < v (i) 0 ; this is to be expected because the system is conservative and the initial precession speed is the maximal possible, being applied in the system at rest in m 1 where E e is minimal.

Figure ?? displays a contourplot of the average precession speed v ≡ φ with respect to the initial precession speed v 0 and the imperfection parameter ε. When ε is vanishing, the system is regular and the dependence of the average velocity with respect to the initial one is smooth and monotonically increasing; the angle variable ϕ is always cyclic. However, for values of the perturbation parameter higher than ε ∼ 10 -2 , the average precession speed could suddenly drop; see, for instance, the region with almost triangular shape in the center of the plot. Therein the precession speed changes sign often during the observation period and this results in small overall value of its average. The dashed curves indicates the isolines where v = 0.16, 0.24 and 0.43: at these speeds the motion seems to avoid chaos until very large imperfections values.

To corroborate quantitatively the onset of chaotic behavior at the critical value ε ∼ 10 -2 , we have computed the Lyapunov spectrum [? ? ]. Being the system autonomous and hamiltonian, only the largest Lyapunov exponent, say λ max , is relevant. We start with the system at rest in the saddle, give random initial conditions and use the numerical algorithm is observed for ε 0.008, this picture being only marginally modified changing the initial energy over a large spectrum of values. For smaller values, it is found consistently that the system is not chaotic and the value of λ max can be considered equal to zero, being of the same order of the numerical error.

III. DERIVATION OF THE MATHIEU EQUATION

a. Analytical estimates. The numerical evidence suggests a prominent role of the precessions within the zero-stiffness manifold. Particular values of the precession speed seem to be more stable and avoid chaotic behaviors, regardless of the value of the imperfection parameter. To understand the precise mechanism, we resort to a perturbative analysis and seek for solutions in the form ϕ(t) = vt + o(ε). For such an ansatz the Lagrangian does not


Luisier
email: philippe.yamato@gmail.com

Reaction-induced volume change triggers brittle failure at eclogite facies conditions

Keywords: volume change, High-Pressure metamorphism, frictional plasticity, numerical modelling

Highlights:

-Densification due to eclogitization reaction can generate significant shear stresses -Degree of reaction overstepping impacts the amount of generated shear stresses -Reaction-induced stresses can trigger brittle failure at high-pressure conditions -Negative volume change and associated brittle failure enhance reaction propagation

Introduction

When buried at depth in convergent zones, rocks are subjected to pressure (P) and temperature (T) changes. Depending on their chemical composition and fluid content, the mineral assemblages and proportions evolve toward more thermodynamically stable assemblages at increasing P-T conditions [START_REF] Spear | Metamorphic phase equilibria and pressure-temperature-time paths[END_REF]. Eclogitization of the oceanic and continental crust involves a large increase in solid density and fluid release, which have significant impacts on tectonic deformation and long-term lithosphere dynamics (e. g., [START_REF] Kirby | Metastable mantle phase transformations and deep earthquakes in subducting oceanic lithosphere[END_REF][START_REF] Hacker | Subduction factory 2. Are intermediate-depth earthquakes in subducting slabs linked to metamorphic dehydration reactions?[END_REF][START_REF] Jung | Intermediate-depth earthquake faulting by dehydration embrittlement with negative volume change[END_REF]Hetenyi et al., 2007;2021;[START_REF] Malvoisin | Sustainable densification of the deep crust[END_REF]. In addition, there is growing evidence supporting a genetic link between eclogitization reactions and short-term geological processes, such as the generation of seismic events at intermediate depth (Fig. 1a, b, Hacker et al., 2003;[START_REF] Jung | Intermediate-depth earthquake faulting by dehydration embrittlement with negative volume change[END_REF][START_REF] Austrheim | Pseudotachylytes generated during seismic faulting and eclogitization of the deep crust[END_REF]Hetenyi et al., 2007;[START_REF] Nakajima | Intermediatedepth earthquakes facilitated by eclogitization-related stresses[END_REF][START_REF] Shi | Lower-crustal earthquakes in southern Tibet are linked to eclogitization of dry metastable granulite[END_REF][START_REF] Incel | Reaction-induced embrittlement of the lower continental crust[END_REF]. However, it is still unclear and there is no consensus yet on the mechanism responsible for this behaviour. Indeed, even the strongest mineral phases are expected to deform in a ductile rather than brittle manner above 700°C and 1 GPa of confining pressure, and for a strain rate of 10 -14 s -1 (e.g. dry plagioclase flow laws, [START_REF] Rybacki | Dislocation and diffusion creep of synthetic anorthite aggregates[END_REF]2004). We thus face a chicken-and-egg problem: it is hard to decipher whether brittle deformation is induced by metamorphic reactions or if the reaction itself triggers embrittlement at high pressure. Several studies have addressed this issue by proposing either of these scenarios. On one hand, reaction-induced embrittlement has been proposed to be caused by (i) the change in strength of the newly formed reaction products (grainsize reduction, nucleation and/or growth of weaker phases, [START_REF] Stünitz | Weakening and strain localization produced by syn-deformational reaction of plagioclase[END_REF][START_REF] John | Generation of intermediate-depth earthquakes by self-localizing thermal runaway[END_REF][START_REF] Thielmann | Intermediate-depth earthquake generation and shear zone formation caused by grain size reduction and shear heating[END_REF][START_REF] Incel | Laboratory earthquakes triggered during eclogitization of lawsonite-bearing blueschist[END_REF], (ii) dehydration of the transforming rocks leading to hydrofracturing through a pore fluid pressure increase [START_REF] Hacker | Subduction factory 2. Are intermediate-depth earthquakes in subducting slabs linked to metamorphic dehydration reactions?[END_REF][START_REF] Jung | Intermediate-depth earthquake faulting by dehydration embrittlement with negative volume change[END_REF][START_REF] Hetényi | Density distribution of the India plate beneath the Tibetan plateau: Geophysical and petrological constraints on the kinetics of lower-crustal eclogitization[END_REF]2021;[START_REF] Alvizuri | Source mechanism of a lower crust earthquake beneath the Himalayas and its possible relation to metamorphism[END_REF], and (iii) volume change [START_REF] Mcgarr | Seismic moments and volume changes[END_REF][START_REF] Pennington | The effect of oceanic crustal structure on phase changes and subduction[END_REF][START_REF] Kirby | Metastable mantle phase transformations and deep earthquakes in subducting oceanic lithosphere[END_REF][START_REF] Malvoisin | Control of serpentinisation rate by reactioninduced cracking[END_REF]2020). On the other hand, it has been proposed that brittle fracturing could be triggered by transient high stress conditions over very short timescales, leading to external fluid infiltration and subsequent transformation of the granulite into eclogite (e.g., [START_REF] Jamtveit | Earthquake-induced transformation of the lower crust[END_REF]. In the latter case, brittle failure at intermediate depths is induced in response to remote seismic events occurring in the seismogenic zone at shallower levels [START_REF] Jamtveit | Earthquake-induced transformation of the lower crust[END_REF], while in the other aforementioned scenarios, brittle failure occurs in situ as a consequence of a reaction-induced stress perturbation. Whatever the mechanism responsible for embrittlement, it is also plausible that the deformation of the lower continental crust results from a cyclic interplay between fracturing-induced reaction and reaction-induced fracturing involving successive stages of fluid pumping and release. Another mechanism that has been proposed to explain embrittlement is volume change. Although embrittlement triggered by diminishing solid density reactions (e. g., serpentinization, melting) has already been suggested (e. g., [START_REF] Kelemen | Reaction-driven cracking during retrograde metamorphism: Olivine hydration and carbonation[END_REF][START_REF] Vrijmoed | An alternative model for ultra-high pressure in the Svartberget Fe-Ti garnet-peridotite, Western Gneiss Region, Norway[END_REF][START_REF] Schmalholz | 2D Hydro-Mechanical-Chemical Modeling of (De) hydration Reactions in Deforming Heterogeneous Rock: The Periclase-Brucite Model Reaction[END_REF][START_REF] Malvoisin | Achieving complete reaction while the solid volume increases: A numerical model applied to serpentinisation[END_REF] and studied using numerical models [START_REF] Malvoisin | Control of serpentinisation rate by reactioninduced cracking[END_REF][START_REF] Ulven | Fracture initiation during volume increasing reactions in rocks and applications for CO2 sequestration[END_REF], densification has received less attention. Moreover, due to complex numerical issues, only recently could frictional plasticity be robustly included into models that account for non-linear visco-elasto-plastic deformation of heterogeneous and compressible materials (e. g., [START_REF] Duretz | Modeling Lithospheric Deformation Using a Compressible Visco-Elasto-Viscoplastic Rheology and the Effective Viscosity Approach[END_REF].

Lower crustal rocks consist of amphibolite to granulite facies rocks, whose densification is mainly controlled by the transformation of plagioclase into omphacite (Fig. 1c-e). It is associated with an important density change of ~ 400-500 kg.m -3 (Fig. 1f) and with a transient variations of material strength [START_REF] Bras | Transient weakening during the granulite to eclogite transformation within hydrous shear zones (Holsnøy, Norway)[END_REF]. Another major characteristic of densification reactions in fluid undersaturated lower crustal rocks, such as amphibolites and granulites, is that they usually take place far from the predicted thermodynamic equilibrium (e.g., [START_REF] Austrheim | Eclogitization of lower crustal granulites by fluid migration through shear zones[END_REF][START_REF] Wayte | A TEM study of disequilibrium plagioclase breakdown at high pressure: the role of infiltrating fluid[END_REF][START_REF] Wain | Metastability of granulites and processes of eclogitisation in the UHP region of western Norway[END_REF][START_REF] Jackson | Metastability, mechanical strength, and the support of mountain belts[END_REF][START_REF] Hetényi | Density distribution of the India plate beneath the Tibetan plateau: Geophysical and petrological constraints on the kinetics of lower-crustal eclogitization[END_REF][START_REF] Hetényi | Metamorphic transformation rate over large spatial and temporal scales constrained by geophysical data and coupled modelling[END_REF]. Metastability of such rocks in subduction zones has been discussed in several studies [START_REF] Lund | Earthquakes in the deep continental crust-insights from studies on exhumed high-pressure rocks[END_REF]Young and Kylander-Clark, 2015;[START_REF] Palin | Subduction metamorphism in the Himalayan ultrahigh-pressure Tso Morari massif: an integrated geodynamic and petrological modelling approach[END_REF] and fast rock transformation, after a delay corresponding to a reaction overstep of several kbars, has been proposed as a cause of intermediate-depth earthquakes [START_REF] Dennis | Earthquakes resulting from metastable phase transitions[END_REF][START_REF] Jamtveit | Disequilibrium metamorphism of stressed lithosphere[END_REF][START_REF] Incel | Reaction-induced embrittlement of the lower continental crust[END_REF][START_REF] Hetényi | Density distribution of the India plate beneath the Tibetan plateau: Geophysical and petrological constraints on the kinetics of lower-crustal eclogitization[END_REF]2021).

Here, we investigate the mechanical response to densification of the lower crust using thermomechanical numerical modelling. The aim of this study is therefore to quantify the amount of shear stress that can be generated by the solid volume change induced by eclogitization, by taking into account reaction pressure overstepping (Fig. 2 a,b). We consider the eclogitization as a densification reaction and voluntary ignore the potential strength variations induced by this transformation (i.e. reaction softening/hardening). Our study can then be considered as an endmember case, in which adding a softening mechanism would further amplify strain localisation that is already generated by densification. We show that, in deforming rocks, the stress generated during densification can indeed overcome the yield stress, hence generating frictional plastic yielding in rocks that would otherwise be in the ductile realm.

Methods

Modelling strategy and validation

We designed numerical models to simulate the progress of eclogitization of a radial inclusion enclosed in a granulitic material (Fig. 2c,d) placed at metastable conditions (i.e., at a background pressure (P BG ) larger than the pressure of the reaction (P r ) predicted by equilibrium thermodynamics) for a fixed temperature (Fig. 2a,b). The initial degree of overstepping is set by initially placing the whole model at a background pressure of 500 MPa (for the reference model) above the pressure of reaction (i.e., P BG = 2.0 GPa and P r = 1.5 GPa).

First, a one-dimensional (1D) model (Fig. 2c) involving linear viscous rheology was used to extract the maximum shear stress obtained during the densification of a segment of the model (r inc ) placed in metastable conditions (i.e., P BG > P r ), assuming no background strain rate (see section 2.3). A constant linear viscous rheology (η = 10 22 Pa.s) is used throughout the model. Parameters are given in Table 1. This model allowed us to perform an intensive parametric study on the following parameters: degrees of overstepping (ΔP m ), density variation (Δρ), shapes of density evolution (d Pr ), reaction kinetics (τ k ), shear modulus (G), bulk modulus (K), and linear viscosities (η). Results of this parametric study are presented in Supplementary Figure 1. Two-dimensional (2D) models (Fig. 2d), built upon the 1D models, were used to incorporate the effect of background strain rate and to account for visco-elasto-viscoplastic (V-E-VP) rheologies (see section 2.4), allowing to better quantify the amount of brittle failure and its spatial evolution.

Both 1D and 2D models were tested against an analytical solution based on the approach of [START_REF] Lee | Self-induced fracture generation in zircon[END_REF] see section 2.2. and Supplementary Figure 2). Matlab scripts to solve the analytical solution and to reproduce the 1D model are provided in the Supplementary Materials. The results of the 1D and 2D models are identical and both converge with increasing grid resolution (see Supplementary Figure 3).

Analytical solution

The analytical solution was developed following the method presented in [START_REF] Lee | Self-induced fracture generation in zircon[END_REF].

To allow for a direct comparison with the 2D models, the original solution was modified to account for cylindrical coordinates. Displacement, pressure and stress fields are caused by the instantaneous density change (from ρ i to ρ f ) of a cylindrical inclusion of radius r inc embedded in an elastic matrix.

The associated volume change generates a self-induced stress, which is proportional to the density change through a constant C that can be expressed as:

(1)

The radial displacement vector (u r ), in cylindrical coordinates, can then be computed analytically for a variable radial distance r using the following equation:

(2) where A and B, that depend on C, are constants. The latter may be determined from the application of boundary and transmission conditions associated with material interfaces (i.e. inclusion for r ≤ r inc and matrix for r > r inc , see script SphereInHole_cylindrical.ipynb for more details). The inclusion and the matrix are both considered elastic with similar material properties (i.e. bulk modulus K=40 GPa and shear modulus G=30 GPa). However, the proposed analytical solution could also be used for different elastic properties between the reacting inclusion and the matrix. The Matlab script corresponding to this analytical solution is provided in supplementary file (AnalyticalSolution.m). It allows to test the validity of both the 1D and the 2D model used in our study (See Supplementary Figure 2).

Description of the 1D model

The 1D model used in this study solves the continuity and momentum equations in the [r, p] polar coordinates. In this coordinate system, we consider the volume change affecting a radial inclusion (with a radius r inc ) embedded in a radial domain (Fig. 2c). In the absence of body forces and inertia, the continuity and momentum equations are expressed as:

(3) and:

(4

)
where ρ is the density, V the velocity, r the distance from the center of the model, and and are the normal stress in r and p direction, respectively. By definition, total stresses ( and )

can be decomposed as:

(5

)
where P is the pressure and and , the deviatoric stresses. Note that here stresses are defined as negative in compression. Material constituting the inclusion and the matrix are defined with visco-elastic rheology (see parameters in Table 1) considering a linear viscosity (η) and a shear modulus (G). Deviatoric stresses can therefore be computed as:

, ( 6
)
where and are the effective deviatoric strain rate components, and η ve corresponds to an effective visco-elastic modulus expressed as follow:

, ( 7
)
where Δt corresponds to the time step used in the model for time discretisation. The effective viscoelastic deviatoric strain rate components are expressed as:

, ( 8
)
and corresponding to the deviatoric stresses obtained in the previous time step (see [START_REF] Duretz | Modeling Lithospheric Deformation Using a Compressible Visco-Elasto-Viscoplastic Rheology and the Effective Viscosity Approach[END_REF]. The driving force in the model is only controlled by the volume change associated with the density change. Motion occurs because the initial pressure (P BG ) is higher than the pressure of reaction (Pr). The reaction progress (see Fig. 2b), depends on pressure P and is expressed as a cumulative distribution function (CDF) and takes the following form:

(9)

where P r is the "pressure of reaction" (reaction line on Fig. 2a) and d Pr corresponds to a standard deviation of the CDF controlling the density change sharpness (Fig. 2b). At a given P, the value of X eq computed using ( 9) therefore corresponds to a value at equilibrium. In order to introduce a reaction kinetics, the value of X considered is computed implicitly using:

(10) t k being the kinetic parameter that defines the characteristic reaction time. The reference density is then computed at each time step as:

(11)

where and corresponds to the initial reference density (2850 kg.m -3 for the granulite) and the final reference density (3250 kg.m -3 for the eclogite), respectively. Density introduced in continuity equation includes pressure dependency and therefore corresponds to:

(12

)
where β = 1/K, is the compressibility (see Table 1). Note that thermal expansivity is not considered and temperature is kept constant throughout the simulations. At each time step, the equations (3-4) are discretized using a staggered finite difference scheme and solved using a pseudo-transient solving strategy (e.g. [START_REF] Duretz | Resolving thermomechanical coupling in two and three dimensions: spontaneous strain localization owing to shear heating[END_REF]. The mesh is advected in a Lagrangian manner following the obtained velocity field. The Matlab script corresponding to this code is provided in supplementary file (Numerics1D_polar.m). It allows to run the 1D reference model and can thus be used to reproduce Figures 3g-i.

Description of the 2D model

The 2D thermo-mechanical simulations were carried out using the code MDoodz6.0, which accounts for a V-E-VP rheological model [START_REF] Duretz | Modeling Lithospheric Deformation Using a Compressible Visco-Elasto-Viscoplastic Rheology and the Effective Viscosity Approach[END_REF]. The continuity and momentum equations are solved in 2D Cartesian [x,y] coordinates, neglecting both gravity and inertial terms, such as:

(13) (14)
The heat transfer equation is expressed as:

(15

)
where k is the conductivity and C p the heat capacity. corresponds to viscoplastic dissipation.

Note that heat transfer does not affect our results: the temperature remains constant at the considered spatio-temporal scale.

The partial differential equations are discretised using a staggered finite difference scheme and markers-in-cell method [START_REF] Gerya | Characteristics-based marker-in-cell method with conservative finitedifferences schemes for modeling geological flows with strongly variable transport properties[END_REF]. The non-linear mechanical equations are solved using a full Newton-Raphson method with analytical derivation of the Jacobian (see [START_REF] Duretz | Modeling Lithospheric Deformation Using a Compressible Visco-Elasto-Viscoplastic Rheology and the Effective Viscosity Approach[END_REF]. Successive Newton corrections are computed using a direct-iterative scheme that relies on Powell-Hestenes iterations, Cholesky factorization of the symmetrised Jacobian operator and Generalized Conjuguate Residual [START_REF] Räss | Concise and efficient MATLAB 2-DS tokes solvers using the Finite Difference Method[END_REF]. The optimal Newton step is computed via a direct line-search globalization procedure [START_REF] Duretz | Shear heating-induced strain localisation across the scales[END_REF]. Marker positions are evolved forward in time using a fourth order (in space) Runge-Kutta scheme. The Courant number is set to where n, A, and Q are the dislocation creep parameters of the material (see Table 1) and R is the universal gas constant (R = 8.314510 J•K -1 •mol -1 ). F pwl corresponds to the correction factor (e.g., [START_REF] Schmalholz | The exponential flow law applied to necking and folding of a ductile layer[END_REF] for invariant formulation relative to the type of experiments used for calibration (here axial compression), such as:

. ( 18
)
In contrast to the 1D simulations, the 2D simulations include the effect of frictional rheology. To this end a viscoplastic Drucker-Prager rheological model is employed 35 . The yield function is expressed as:

(19)
where is the deviatoric second stress invariant, P is the pressure, C is the cohesion and is the friction angle. The viscoplastic parameter, , is set such that the overstress ( ) is on the order of 0.1 MPa (i.e. η vp = 5 × 10 18 Pa.s for = 10 -14 s -1 ) in our simulations. Tests of convergence for 2D models involving plasticity can be found in Supplementary Figure 5. The effect of plastic dilation is neglected and the plastic flow potential (Q) thus simply reads as . The reaction progress law and reaction kinetics are included in a similar manner as in the 1D code described above (see eq. 9 to 12).

Choice of parameters for reference models

In this study, eclogitization is restricted to a progressive densification process (Fig. 2b). The choice of parameters for defining a reference model is based on data from granulite to eclogite transformation, that has been intensively studied in the Bergen Arcs. Based on recent chemical data and P-T estimates performed in Holsnøy [START_REF] Bhowany | Phase equilibria modelling constraints on P-T conditions during fluid catalysed conversion of granulite to eclogite in the Bergen Arcs, Norway[END_REF], we set P r for the granulite to eclogite transformation at 1.5 GPa for a temperature of 680°C. A value of 700 MPa was chosen as it seems to be in good agreement with the density evolution that encompass several reactions occurring with a P increase (Bras et al., 2021, Fig. 1f). Thermodynamic modelling providing the Fig. 1f was performed using Perple_X_6.9.1 [START_REF] Connolly | Multivariable phase diagrams; an algorithm based on generalized thermodynamics[END_REF]2009) and the internally consistent thermodynamic database hp62ver.dat [START_REF] Holland | An improved and extended internally consistent thermodynamic dataset for phases of petrological interest, involving a new equation of state for solids[END_REF]. The solution models are for feldspar (Furhman and Lindsley, 1988), garnet [START_REF] White | New mineral activity-composition relations for thermodynamic calculations in metapelitic systems[END_REF], clinopyroxene [START_REF] Green | An order-disorder model for omphacitic pyroxenes in the system jadeite-diopside-hedenbergite-acmite, with applications to eclogitic rocks[END_REF], clinoamphibole [START_REF] Green | Activity-composition relations for the calculation of partial melting equilibria in metabasic rocks[END_REF], epidote [START_REF] Holland | An improved and extended internally consistent thermodynamic dataset for phases of petrological interest, involving a new equation of state for solids[END_REF], white mica [START_REF] White | New mineral activity-composition relations for thermodynamic calculations in metapelitic systems[END_REF], and melt [START_REF] Green | Activity-composition relations for the calculation of partial melting equilibria in metabasic rocks[END_REF]. Bulk composition used corresponds to sample HOL7C_2014 [START_REF] Bhowany | Phase equilibria modelling constraints on P-T conditions during fluid catalysed conversion of granulite to eclogite in the Bergen Arcs, Norway[END_REF]. A difference in density of 400 kg.m -3 using a reference density of 2850 kg.m -3 for the granulite (X=0) and of 3250 kg.m -3 for the eclogite (X=1) appears also in good agreement with previous studies [START_REF] Austrheim | Eclogitization of lower crustal granulites by fluid migration through shear zones[END_REF]Hetenyi et al., 2007). The kinetic parameter t k has been fixed to 100 yrs as there are arguments that granulite to eclogite transformation is rather a fast process [START_REF] Malvoisin | Sustainable densification of the deep crust[END_REF]. The choice of the anorthite dry [START_REF] Rybacki | Dislocation and diffusion creep of synthetic anorthite aggregates[END_REF] was stimulated by the fact that initial material corresponds to dry granulite. K and G values [START_REF] Motra | Influence of loading and heating processes on elastic and geomechanical properties of eclogites and granulites[END_REF]. For 2D model involving plasticity, cohesion and friction angle were set to 50 MPa and 30°, as they correspond to an extrapolation of [START_REF] Byerlee | Friction of rocks[END_REF] data at high pressure. The strain rate used in the 2D reference model is 10 -14 s -1 . in a radial configuration around the inclusion (Fig. 3a-c). A fast, transient and important shear stress increase, larger than 500 MPa, occurs in the vicinity of the reacting area at a time close to the onset of reaction (Fig. 3a). After reaching a maximum value (here of ~0.6 GPa), shear stresses then decrease following a Maxwell-type relaxation [START_REF] Dabrowski | Viscous relaxation of grain-scale pressure variations[END_REF] until they become null once equilibrium density is reached.

Results

Results from the reference model

Simultaneously, the densification reaction induces a pressure decrease in both the matrix and the inclusion (Fig. 3 d-f). In the inclusion, there is first a rapid pressure drop (see Fig. 3d compared to the 2.0 GPa initial pressure in both inclusion and matrix, and green line in Fig. 3g-h), followed by a consecutive pressure increase to reach the same pressure as in the matrix (Fig. 3g,h). The pressure decrease in the matrix is smoother (black line in Fig. 3g-h). It is associated to the increase in matrix volume accommodating the shrinking of the inclusion, to satisfy mass conservation. An important point to note is that the reaction stops as soon as both the inclusion and the matrix reached the same density, since our model is only driven by density change, which is not necessarily the one predicted for the reaction to be complete (i.e., X=1).

As the value of the maximum shear stress reached (max(τ II )) constitutes a useful metric that indicates whether reaction-induced volume change can lead to rock failure, we have tracked this value through a complete parametric study (see Supplementary Figure 1). Figure 4a shows the effect of pressure overstepping and density change associated to the metamorphic reaction on the value of max(τ II ). The results show that density variations produce shear stress. They also show that when a metamorphic reaction is associated with a change in density larger than a few tens of kg.m -3 , pressure overstepping always leads to shear stress generation, and that the larger the density change is, the higher the effect of the overstepping is.

The proportion of reacting material is also a key factor (Fig. 4b). Maximum shear stress is reached when the amount of initial reacting material is small (<10%) in comparison to the total amount of material. If the initial reacting area is higher than 20-30%, the maximum shear stress reached is much less important.

Influence of densification on frictional plastic shear banding

The mechanical consequences of volume changes in an isoviscous host-inclusion system are different from that in host-inclusion systems involving viscosity contrasts and subjected to shear (e.g. [START_REF] Schmid | Analytical solutions for deformable elliptical inclusions in general shear[END_REF][START_REF] Moulas | Stress field associated with elliptical inclusions in a deforming matrix: Mathematical model and implications for tectonic overpressure in the lithosphere[END_REF][START_REF] Luisier | Metamorphic pressure variation in a coherent Alpine nappe challenges lithostatic pressure paradigm[END_REF]. Figure 5 compares pressure and shear stress fields obtained for models including reaction-induced volume change (Fig. 5a,b) with models involving a viscosity ratio of 1000 between inclusion and matrix (Fig. 5 c,d), under a background strain rate of 10 -14 s -1 .

In Figure 5a, inclusion and matrix have the same viscosity and frictional plasticity is not activated.

In this case, there is no significant change in pressure around the inclusion. The shear stress (τ II ) generated by the negative volume change of the inclusion is radial (see Fig. 3a), but because of the pure shear background strain rate (horizontal shortening and vertical stretching), the shear stress field is different in vertical and horizontal direction. τ II values are higher in the vertical direction where shear stresses due to the volume change and shear stresses due to stretching are in the opposite sense. For the same reason, τ II values are smaller in the horizontal direction where shortening accompanies the shrinking of the inclusion. When frictional plasticity (i.e. pressure dependent) is accounted for (Fig. 5b), we observe that plastic yield stress in matrix is reached in the area undergoing stretching (where τ II values are high) whereas it is not reached in the area undergoing shortening (where τ II values are low). The regions undergoing plastic deformation (contoured in white in Fig. 5b) are located close to the reacting inclusion in the stretching direction.

Figure 5c displays pressure and shear stress field obtained for a model using the same parameters as in Fig. 5a, except that (i) there is no volume change and (ii) inclusion and matrix have different viscosities (viscosity 1000 times weaker in the inclusion). In this case, the shear stress field is identical in vertical and horizontal directions, but the pressure field is not. Due to the presence of a weak inclusion, a pressure increase occurs in the matrix close to the inclusion in the vertical direction (parallel to stretching). In the horizontal direction, the pressure decreases. Pressure and shear stress fields are therefore drastically different than in the previous case involving volume changes without viscosity contrast. In the "viscosity-contrast" case, shear stress is maximum where the pressure is the lowest in the matrix, close to the inclusion in horizontal direction. These areas constitute favourable locations for frictional plastic yielding. In the vertical direction, because of locally higher pressure values, a larger shear stress is required to reach the yielding. This is the reason why the parts undergoing plastic deformation (white contours in Fig. 5d) are located close to the reacting inclusion in the shortening direction, in the case where plasticity is switched on.

In Figure 5, viscosities are considered stress dependent with a stress exponent n=3. Even if n is known to affect relaxation time and pressure evolution [START_REF] Dabrowski | Viscous relaxation of grain-scale pressure variations[END_REF], performing the same models considering a constant viscosity with n=1 does not drastically change the results concerning matrix deformation over the considered time scale (see Supplementary Figure 4).

Brittle failure in ductile realm

Our results show that depending on ΔP m , the amount of shear stress produced (max(τ II )) by considering only reaction-induced volume change is significant (Fig. 4). For a given material, at a fixed pressure-temperature-strain rate (P-T-) conditions in the ductile field, shear stress (τ v )

supported by this material is below the yield stress value (τ y ) required to fail (Fig. 6a). By definition, this difference in stress Δτ (with Δτ = τ y -τ v ) can be overcome by (1) increasing strain rate, (2) decreasing pressure and/or (3) decreasing temperature. We here show that shear stresses resulting from reaction-induced volume change can also overcome Δτ and thus trigger yielding. The occurrence of ductile and brittle deformation is therefore strongly linked to volume changes associated with metamorphic reactions.

Figure 6b displays the value of Δτ as a function of P and for dry anorthite [START_REF] Rybacki | Dislocation and diffusion creep of synthetic anorthite aggregates[END_REF]. Δτ isocontours are presented in black. The brittle-ductile transition (BDT, yellow line) corresponds to Δτ = 0. Conditions of the parametric study performed for different degrees of overstepping (right vertical axis) and different background strain rates considering a density difference Δρ of 400 kg.m -3 (i.e., roughly as in the granulite to eclogite transformation) are also presented (see symbols). White lines correspond to isovalues of the amount of anorthite matrix that undergone brittle failure in case of reaction-induced volume change. An illustration of the effect of the background strain rate is presented in Fig. 6c-g. Results show that the field where yield stress is achievable due to reaction-induced volume change covers a large part of the domain where ductile deformation is expected for dry anorthite if metamorphic transformation were not activated (Fig. 6b). By accounting for the stress induced by volume change, the position of the brittle-ductile boundary can be shifted towards lower strain rate values, reducing the required strain rate necessary to fail the rock up to two orders of magnitude. Destabilization of a ductile metastable granulite into eclogite anywhere in this domain can lead to brittle failure.

Discussion and conclusion

Our results show that eclogitization-induced densification occurring under metastable conditions can, by itself, generate sufficient shear stress to fail the rocks at high-pressure conditions. Brittle failure in rocks can therefore be explained by densification reactions, without even considering the rheological changes induced by the reaction. In light of these results, it is now possible to further discuss the applicability of our model on natural cases and the broader implications of this study on the strength of the continental lithosphere subjected to metamorphic reactions. These points will be addressed in the following, after a discussion on possible perspectives concerning our modelling approach.

Modelling perspectives

In the present study we intentionally did not consider the potential rheological changes associated with the metamorphic transformation. Modification of strength is transient and only caused by densification. This choice has been made to isolate and study in detail the effect of densification alone. One important consequence of this is the value of pressure obtained in the transforming area, which is lower than the pressure in the unreacted matrix. However, in models considering only strength difference between inclusion and matrix (as in Fig. 5c,d), the weaker the inclusion is compared to its matrix, the higher the pressure difference will be between the inclusion and matrix [START_REF] Schmid | Analytical solutions for deformable elliptical inclusions in general shear[END_REF][START_REF] Moulas | Stress field associated with elliptical inclusions in a deforming matrix: Mathematical model and implications for tectonic overpressure in the lithosphere[END_REF][START_REF] Luisier | Metamorphic pressure variation in a coherent Alpine nappe challenges lithostatic pressure paradigm[END_REF]. We can therefore speculate that by taking into account a viscosity decrease of the inclusion in addition to volume change during the transformation would lead to higher pressure inside the inclusion compared to models that only consider volume changes. Depending on the magnitude of such a viscosity decrease, the pressure inside the inclusion could even be higher than the pressure in the surrounding matrix. In such a case, this pressure increase would extend the degree of overstepping of the rock. This is something that can be expected for the eclogitization process in specific cases, such as the granulite to eclogite transformation, as it is classically proposed that the eclogite is weaker, at least transiently, than its granulitic protolith [START_REF] Labrousse | Initiation of crustal-scale thrusts triggered by metamorphic reactions at depth: Insights from a comparison between the Himalayas and Scandinavian Caledonides[END_REF][START_REF] Bras | Transient weakening during the granulite to eclogite transformation within hydrous shear zones (Holsnøy, Norway)[END_REF][START_REF] Kaatz | Widening of hydrous shear zones during incipient eclogitization of metastable dry and rigid lower crustholsnøy, western Norway[END_REF].

Our models assume that the third direction of space is infinite. Our analytical solution (see section 2.2) shows that the stresses induced by volume change ("self-induced stresses") in cylindrical coordinates corresponds to 2/3 of the self-induced stresses obtain in spherical coordinates [START_REF] Lee | Self-induced fracture generation in zircon[END_REF][START_REF] Gillet | Coesite in subducted continental crust: PT history deduced from an elastic model[END_REF]. Hence, considering the eclogitization process in 3D would therefore lead to higher self-induced stresses than in 2D. In that sense, our estimates of the ability to reach brittle failure can then be taken as minimum.

We assume conservation of mass and thus consider a closed system. A variation of density therefore only impacts volume change. However, it is well known that any replacement reaction, whether it is a volume increase or reduction, involves the formation of intra-crystalline porosity along with the product minerals, as evidenced by pseudomorph microtextures [START_REF] Putnis | Mineral replacement reactions: from macroscopic observations to microscopic mechanisms[END_REF]2009;[START_REF] Putnis | Transient porosity resulting from fluid-mineral interaction and its consequences[END_REF]. In the case of an open system subjected to input of externally derived fluids, such a porosity can be filled by product minerals, thus limiting or even inhibiting the volume change [START_REF] Centrella | Coupled mass transfer through a fluid phase and volume preservation during the hydration of granulite: An example from the Bergen Arcs, Norway[END_REF].

Similarly, mineral replacement can occur via a process of isovolumetric interface-coupled dissolution-precipitation during open system fluid-rock interaction [START_REF] Putnis | Mineral replacement reactions[END_REF][START_REF] Centrella | Coupled mass transfer through a fluid phase and volume preservation during the hydration of granulite: An example from the Bergen Arcs, Norway[END_REF]. From this point of view, our model constitutes an end-member that can however be applied in cases where mass balance is fulfilled [START_REF] Centrella | Textural and chemical evolution of pyroxene during hydration and deformation: a consequence of retrograde metamorphism[END_REF].

Applications and Implications

The degree of overstepping a rock can sustain before reacting plays a major role on its mechanical behaviour. Indeed, our results show that it extends the brittle field toward lower strain rate values (Fig. 6). Many examples have been reported where the plagioclase breakdown reaction is thought to have been inhibited because of the lack of free fluid [START_REF] Austrheim | Eclogitization of lower crustal granulites by fluid migration through shear zones[END_REF][START_REF] Wayte | A TEM study of disequilibrium plagioclase breakdown at high pressure: the role of infiltrating fluid[END_REF][START_REF] Rubie | Disequilibrium during metamorphism: the role of nucleation kinetics[END_REF][START_REF] Jackson | Metastability, mechanical strength, and the support of mountain belts[END_REF]Young and Kylander-Clark, 2015;[START_REF] Jamtveit | Disequilibrium metamorphism of stressed lithosphere[END_REF][START_REF] Palin | Subduction metamorphism in the Himalayan ultrahigh-pressure Tso Morari massif: an integrated geodynamic and petrological modelling approach[END_REF]. The preservation of metastable plagioclase at high pressure can generate pressure overstepping of the order of several 100's of MPa as already suggested by [START_REF] Wayte | A TEM study of disequilibrium plagioclase breakdown at high pressure: the role of infiltrating fluid[END_REF] for the Allalin metagabbro (Alps) and [START_REF] Hetényi | Density distribution of the India plate beneath the Tibetan plateau: Geophysical and petrological constraints on the kinetics of lower-crustal eclogitization[END_REF]2021) for the Indian lower crust (beneath the Tibetan Plateau), and sometimes up to pressure higher than 1.0 GPa (e.g., [START_REF] Wain | Metastability of granulites and processes of eclogitisation in the UHP region of western Norway[END_REF].

Densification of rocks triggered by mineral reactions that take place far from the equilibrium constitutes therefore a mechanism for explaining the nucleation of mechanical instabilities at depth during the burial of a dry continental lithosphere and it is in good agreement with kinetics laws (e.g. Cahn, 1956). In the case of a densification reaction such as the eclogitization of plagioclase-bearing rocks (e.g. Norvegian or Himalayan granulites, Alpine gabbros), it is hence appropriate to study the effects of volume change on the initiation of such instabilities.

Assessing the seismic potential of such a failure occurring at intermediate depth in continental subduction zones is challenging and remains under the scope of this study. The question regarding their potential seismogenic nature remains open and could be answered by using dynamic rupture models (Gabriel et al., 2020, Thomas and[START_REF] Thomas | Dynamic evolution of off-fault medium during an earthquake: a micromechanics based model[END_REF] and/or laboratory experiments [START_REF] Schubnel | Deep-focus earthquake analogs recorded at high pressure and temperature in the laboratory[END_REF].

Our results also show that the mechanical effect of reaction-induced volume change on deformation is important even if only a small amount of material is transformed. Reaction-induced shear stress is indeed maximum when the initial amount of material allowed to react is small (<10%, see Fig. 4b).

Maximum shear stress values then gently decrease when the ratio between the initial volume of reacting material and total volume of rocks increase (Vol inc /Vol tot in Fig. 4b), reaching a plateau of low shear stresses for values of Vol inc /Vol tot ~30-40%. This is in good agreement with the absence of seismicity observed in subducting continental crust where eclogitization progress exceed 30-40%

(see Fig. 1a,b). This could explain why pseudotachylytes are only observed in partially transformed granulites from Holsnøy (Fig. 1d-f), where the percentage of eclogitization is below ~40% [START_REF] Boundy | Structural development and petrofabrics of eclogite facies shear zones, Bergen Arcs, western Norway: implications for deep crustal deformational processes[END_REF].

In addition, our results show that even if frictional yielding is reached, the amount of material that undergoes plastic deformation is generally low (i.e., in the order of 1% or less, see Fig. 6).

However, this should not be overlooked. When yielding occurs, the pressure in the damaged area drops. In contrast, pressure increases at the tip of the area that has undergone plastic deformation (Fig. 5b). The development of such a pressure gradient could be responsible for a transfer of fluid that may favourably circulate through the damaged zone. Because transformation of granulite into eclogite is mainly driven by the presence of fluids, we can speculate that, if even a small amount of fluid is present in the granulite, it will favourably be collected toward the damage zone and then allow transformation of this area. Figure 7 displays a model where this hypothesis is tested, by allowing areas that have undergone plastic deformation to react (i.e., by allowing density change as it is the case for the initial reacting inclusion). It illustrates progressive eclogitization by a runaway mechanism, which lasts as long as the reaction is not limited by the lack of reactive product. In that case, eclogitization of the granulite self-propagates in the direction perpendicular to the shortening direction, forming finger-like structure of eclogitized material as described in the Holsnøy area (e.g., so-called "static eclogites", [START_REF] Putnis | Preservation of granulite in a partially eclogitized terrane: Metastable phenomena or local pressure variations?[END_REF]Zertani et al., 2019[START_REF] Jamtveit | Accelerated hydration of the Earth's deep crust induced by stress perturbations[END_REF]2016).

This study has therefore strong implications for our understanding of the state of stress of the continental lithosphere in convergence zones. In such geodynamic contexts, rocks buried at depth are mainly subjected to metamorphic reactions leading to an increase in density. d Pr = 700 MPa) but inclusion and matrix have the same material property (Table 1). Viscosity present a stress exponent (n=3). In contrast, in (c) and (d), volume change due to the reaction is switched-off and the only difference between the inclusion and the matrix corresponds to a viscosity ratio (reference viscosity in the inclusion is 10 3 times weaker than in the matrix, and n=3). White contours delimitate the parts of the matrix that reached the yield stress and have undergone brittle failure. 1) at T = 680°C. Symbols correspond to simulations performed using the 2D model (see section 2.4.) at T = 680°C for P r = 1.5 GPa and d Pr = 700

MPa. Results (extracted after t = 1000 yrs) allows to define the ratio between the amount of matrix area that reached the yield stress (A brittle , red colours in (c-g)) and the total amount of matrix material (A matrix ). Areas depicted in orange and red correspond to zones that reached a plastic strain of 10 -5 and 10 -3 , respectively.

White dashed lines in (b) correspond to isocontours of this ratio (in %). BDT = Brittle-ductile transition. Five different symbols are used to group the similar results: Grey dots are used for models similar to (c), i.e.

where brittle features are not significant; + signs are used for models as in (d) and (e), x signs are used when shear zones start to propagate as in (f); stars are used when well developed shear zones propagate as in (g); red stars correspond to models belonging to the brittle domain as illustrated in Supplementary Figure 5.

0. 2 .

 2 The implementation of the V-E-VP rheological model follows Yamato et al. (2019) and Duretz et al. (2021). The contribution of each deformation mechanism is computed iteratively to satisfy the following Maxwell model (e.g., Popov & Sobolev, 2008): (16)where the v, e and vp superscripts refer to the viscous, elastic, and visco-plastic (frictional) contributions to the total strain rate. In contrast to the 1D model, a dislocation creep flow law is considered in the 2D model. The flow stress ( ) is thus stress and thermal dependent and is expressed as:(17)

Figure 3

 3 Figure 3 displays the results obtained in absence of background strain rate and for a linear viscosity set to 10 22 Pa.s. It shows that reaction-induced volume change affects both shear stress and pressure

Figures

  Figures and captions:

Figure 1 Fig. 1

 11 Figure 1

Figure 2 Fig. 2

 22 Figure 2

Figure 3 Fig. 3

 33 Figure 3

  (i) as follows: (g) Evolution of pressure with reaction progress (X) in matrix (black line) and in the reacting area (green line). Circle, asterisk, plus and multiplication symbols corresponds to different times (0.5, 10, 50, 100 and 250 kyrs, respectively). Blue line is for a rock in thermodynamic equilibrium. (h) Evolution of pressure with time (t) for matrix (black line) and reacting area (green line). (i) Evolution with time of maximum of shear stress (max(τ II )) recorded in the matrix.

Figure 5 Fig. 5

 55 Figure 5

Figure 6 Fig. 6

 66 Figure 6

Figure 7 Fig. 7

 77 Figure 7

reference model (with background strain rate and visco-elasto-plastic rheologies)**

  Our results show that the nucleation of brittle instabilities due to the densification reaction taking place far from

		Yamato, P., Duretz, T., & Angiboust, S. Brittle/ductile deformation of eclogites: insights from Table 1 Parameters used in reference models
		numerical models. Geochemistry, Geophysics, Geosystems, 20 (7), 3116-3133 (2019).
	Young, D. J., & Kylander-Clark, A. R. C. Does continental crust transform during eclogite facies Symbol Meaning Unit Value
	equilibrium can generate important transient stresses, hence extending the domain where frictional 1D reference model (without background strain rate and using linear viscosity)
	P BG	metamorphism?. Journal of Metamorphic Geology, 33 (4), 331-357 (2015). Background pressure GPa	2.0
	P r	plastic yielding can occur. Pressure of reaction	GPa	1.5
	ΔP m d Pr ρ 1 ρ 2 t k K	Pressure overstepping (P BG -P r ) Zertani, S., Labrousse, L., John, T., Andersen, T. B., & Tilmann, F. The interplay of Eclogitization GPa 0.5 and deformation during deep burial of the lower continental crust-A case study from the Bergen GPa 0.7 P r standard deviation Initial density (before transformation) kg.m -3 2850 Final density (after total transformation) kg.m -3 3250 Kinetic parameter* s 3.1558×10 9 Arcs (Western Norway). Tectonics, 38 (3), 898-915 (2019). Bulk modulus GPa 80
	G	Shear modulus	GPa	40
	η	Viscosity	Pa.s	10 22
	2D BG	Background strain rate	s -1	10 -14
	T BG A	Temperature Pre-exponential factor (dislocation creep)***	°C Pa -n .s -1	680 5.0119×10 -6
	n	Exponent (dislocation creep)***	-	3
	E a	Activation energy (dislocation creep)***	kJ.mol -1	648
	C	Cohesion	MPa	50
	φ	Friction angle	°	30
	η vp C p k	Viscoplastic viscosity Heat capacity Conductivity	Pa.s J.kg -1 .K -1 W.m -1 .K -1	5 ×10 18 1050 2.3
	*see Supplementary Fig. 1 for more details **Parameters that are the same in both 1D and 2D model (i.e.
	elastic and metamorphic reaction properties) are not recalled for 2D reference model. ***Values from
	Rybacki and Dresen (2000) for anorthite dry.

Table 1

 1

). For comparison, results from the 1D model (see section 2.3. and Table 1 for parameters) are presented in (g) to

[START_REF] Bjørnerud | Pressure-related feedback processes in the generation of pseudotachylytes[END_REF]. (e) Granulite fully transformed into eclogite with garnet (red), omphacite (green) and zoisite-kyanite-quartz (yellowish assemblage). (f) Possible evolution of rock's density with pressure during the granulite to eclogite metamorphic transformation. Black stars correspond to density data extracted from Perple_X_6.9.1[START_REF] Connolly | Multivariable phase diagrams; an algorithm based on generalized thermodynamics[END_REF] 2009) at T = 680°C (see details in section 2.5.). Blue curve corresponds to a possible fit of the density evolution using a cumulative distribution function with a mean pressure (P r ) of 1.35 GPa and a standard deviation (d Pr ) of 0.7 GPa as defined in Fig.2band considering reference densities of 2850 kg.m -3 and 3290 kg.m -3 for granulite and eclogite, respectively.

Acknowledgements

P.Y. thanks the Institut Universitaire de France for financial support. We also thank Alexander Webb and an anonymous reviewer for their comments and the editorial handling which contributed to the improvement of the final version of the manuscript.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

All numerical model outputs are available from the corresponding author upon request. They can also be reproduced by using the scripts provided in supplementary materials.

Code availability

Matlab scripts for the analytical solution and the 1D model are provided in supplementary materials. The 2D code is available here: https://github.com/tduretz/MDOODZ6.0 (input file for the 2D reference model is given in supplementary materials). For further requests on the numerical code please contact the authors.

Declaration of competing interest

Vol inc /Vol tot = πr inc 2 /L 2 .