Dehydration‐induced instabilities at intermediate depths in subduction zones

We formulate a model for coupled deformation and dehydration of antigorite, based on a porosity‐dependent yield criterion and including shear‐enhanced compaction. A pore pressure and compaction instability can develop when the net volume change associated with the reaction is negative, i.e., at intermediate depth in subduction zones. The instability criterion is derived in terms of the dependence of the yield criterion on porosity: if that dependence is strong, instabilities are more likely to occur. We also find that the instability is associated with strain localization, over characteristic length scales determined by the hydraulic diffusivity, the elasto‐plastic parameters of the rock, and the reaction rate. Typical lower bounds for the localization length are of the order of 10 to 100 for antigorite dehydration and deformation at 3 GPa. The fluid pressure and deformation instability is expected to induce stress buildup in the surrounding rocks forming the subducted slab, which provides a mechanism for the nucleation and propagation of intermediate‐depth earthquakes.

fluids from the dehydration reaction. These assumptions lead to the following governing 86 equation for pore pressure (see full derivation in Appendix A1): where m 0 d is the total fluid mass that can be released by the reaction per unit rock volume,

Rheology
It is well established experimentally that antigorite aggregates undergo a brittle to duc-95 tile transition at confining pressures of the order of 300 to 400 MPa [Escartín et al., 1997], 96 and that this transition depends weakly on temperature (within antigorite's stability field).

97
Near the dehydration temperature of antigorite, the ductile behaviour remains dominated where p is the Terzaghi effective mean stress (p = p + p f , where p is the mean stress), τ

Reaction rate
A very general formulation of mineral reaction rates is given by [Lasaga and Rye, 1993] 136 ∂ξ ∂t = κA rlm s|∆G| nr , where κ is the temperature-dependent kinetic constant (typically following an Arrhenius 137 law), A rlm is the specific surface area of the rate-limiting mineral, ∆G is the Gibbs energy where c = ∂∆G/∂p f and p eq is the pore pressure at equilibrium. Following Wang and 143 Wong [2003], we can rewrite 5 as where r 0 is a reference reaction rate, and A 0 rlm is the surface area of the rate limiting 145 mineral at the reference rate. The change in surface area is not very well constrained 146 by experimental data. Hence, for simplicity, in the following we will consider the ratio 147 A rlm /A 0 rlm ≈ 1. This simplification is valid only when the reaction progress is small, i.e., 148 when depletion of the reaction is negligible (ξ 1). 149

Geometry and stress equilibrium
We consider a simple system made of a uniform horizontal layer of antigorite, sufficiently 150 extended so that lateral strains can be neglected (i.e., the system is invariant in the plane 151 of the layer). In this geometry, the vertical stress is given by p − 2τ / √ 3, and stress 152 equilibrium requires that: where y denotes the vertical coordinate. The boundary condition is a constant applied 154 vertical stress, which implies that In such a geometry, only the vertical strain component is nonzero. Because no lateral 156 deformation is allowed, the vertical strain is equal to the volumetric strain. Therefore, : and  Figure 2. For simplicity, we use a modified Cam-clay yield surface, given by where C is the critical state line ratio, b is the tensile strength, and p * is the compaction 165 yield pressure (following the notation of Wong and Baud  in Figure 2), we find a reasonable fit with p * = 130 MPa, while the inferred porosity of 172 the sample was around 19%. The partially dehydrated samples (diamonds in Figure 2 The friction coefficient µ is given by the local slope of the yield envelope (equation A24): In the framework of associated plasticity, we could assume that the dilatancy factor is 197 merely equal to the friction coefficient. However, it is well known that rocks do not follow 198 associated flow rules, and hence we shall leave the dilatancy factor β as a free parameter, 199 and explore how it influences the stability of compaction in our model. Likewise, we 200 will leave the hardening coefficient as a free parameter, in order to encompass the widest 201 possible range of behaviours.

202
The permeability of the rock is expected to vary as a function of porosity, and hence resulting viscosity of water ranges from η = 6.9 × 10 −5 to 4.9 × 10 −7 Pa s at p f = 2 GPa 220 and 600 • C and p f = 5 GPa and 700 • C, respectively.

221
The effective compressibility of the rock, c b , does not play any role in the stability of     Table 2.

236
In addition, the knowledge of the molar volume of antigorite allows to compute precisely 237 the total potential mass of water releasable by each reaction, m 0 d . Because antigorite is not very compressible, the main factor influencing m 0 d is the stoichiometry of the reactions.

239
The computed averages of m 0 d for each reaction are shown in Table 2.

240
The reaction kinetics of antigorite as a function of pressure is not well constrained by 241 existing experimental data, which typically focus on the effect of temperature. However, where R is the gas constant. In terms of reaction progress ξ, we rewrite equation 19 as where The molar volume of antigorite V m can be taken as the average along the phase boundary Using the parameter value reported in Table 2 for the reaction of antigorite into forsterite 259 and enstatite, we obtain c ≈ 8.1 × 10 −5 m 3 /mol atg . The equilibrium pressure p eq and the 260 appropriate temperature T can be found from the phase boundary (see Figure 3). As a with the additional term −f β occurring here due to the explicit dependence of the yield reflects that the yield cap expansion due to shear hardening is in fact offset by a further 281 yield cap expansion (resp. shrinkage) due to shear-induced compaction (resp. dilatancy).

282
Note, in passing, that we have derived here a slightly more general case for the com- Assuming constant parameters, this equation has an analytical solution, which is if n r = 1, and if n r > 1, where p 0 f is the initial pore pressure in the system (a small perturbation above 295 the equilibrium pressure).

296
D R A F T July 10, 2017, 12:11pm For linear kinetics, the growth is exponential, while for nonlinear kinetics the growth 299 corresponds to a finite time blow-up. In practice, this distinction is unimportant since 300 pore pressure diffusion, as well as other nonlinearities not accounted for in our simplified Appendix C1 for details): Using the numerical values detailed in the previous section, we remark that β crit is always 309 positive. Hence, for the cases of interest where β < 0 (shear enhanced compaction), we 310 always have β < β crit , and thus the condition for instability is simply f > f crit .  of the form (11); here we restrict our analysis to a study of the stability of the system to 333 small departures from equilibrium (which corresponds to p f = p eq ).

334
As stated in Section 3.3, the experimental data of Eggler and Ehmann [2010] show that 335 the near-equilibrium kinetics is linear, n r = 1. In that case, the stability analysis detailed 336 in Appendix C2 shows that pore pressure runaways are possible if f > f crit (same as 337 D R A F T July 10, 2017, 12:11pm . (31) We observe that the critical wavelength depends on a number of parameters, some of them 339 well constrained (equilibrium pressure p eq , permeability k, fluid viscosity η, and the net initially non porous. We can use this fact to our advantage by noticing that the critical 345 wavelength λ crit tends to a constant, nonzero value for f f crit : (32) Equation (32) provides a simple lower bound for the critical wavelength, which is, quite 347 remarkably, independent from the hardening modulus h.

348
The value of λ crit is shown in Figure 5 as a function of f and β for the parameters 349 relevant to the dehydration of antigorite into enstatite and forsterite. As expected, the 350 wavelength tends to the constant given by (32) at large values of f , and we confirm that 351 this limit value has only a mild dependence (square root) on the dilatancy factor β.

Numerical tests
We performed a series of numerical computations in order to explore further the ef-357 fects of potential nonlinearities associated with the variations of mechanical, hydraulic 358 and chemical parameters during deformation and reaction. We included a power-law 359 dependence of the permeability on porosity, k ∝ ζ 3 (where we recall that ζ is the 360 porosity of the rock minus its poroelastic variations), and accounted for depletion of 361 antigorite by using a simple first-order approximation for the reactant surface area, (7). The numerical method is described fully in with periodic boundary conditions, and an initial sinusoidal infinitesimal pore pressure 367 perturbation is added to the homogeneous initial conditions. We chose a representative 368 example by using an initial pore pressure of 3 GPa, an initial total mean stress of 3.61 GPa 369 (i.e., an initial shear stress of 0.87 GPa), and an initial porosity of 3 %. Using the initial, 370 reference parameters, we find that the critical wavelength for instability is λ crit ≈ 0.08L, 371 so that we expect some compaction localisation (at least transiently).

372
The volumetric strain profile within the layer is shown as a function of time in Figure   373 6(a). The initial ( (y, 0) = 0) and final (at t × r 0 = 100) profiles are highlighted in black,

Model assumptions and limitations
The key assumption of the model presented here is that ductile deformation of antig- (see Figure 5). This similarity in compaction length scale between the two types of model 456 essentially arises from the similarity in timescales between viscous flow and reaction rate.

457
Here we made the assumption that the material deforms under isothermal conditions. viscous rock deformation, we find here that the path towards low effective stress states is 505 unstable, but only in the case when the net volume change is negative 1/ρ f + ∆ r V s < 0.

506
This instability arises because the compaction tends to increase pore pressure, and takes 507 the system further away from equilibrium. If the net volume change from the reaction is 508 positive, any compaction and pore collapse would increase pore pressure and bring the 509 system back to equilibrium, unless unrealistic amounts of shear-induced dilatancy occur.

510
One interesting outcome of our model is that we show that the dehydration and com-511 paction process produces significant shear deformation in the rock, and not just pure This is most likely the case when porosity is larger than several percents, above the per-536 colation threshold [Guéguen and Palciauskas, 1994]. However, over long timescales, the 537 progressive drainage of the pore fluid outside of the dehydrating zones tends to allow com-538 paction to reduce porosity; this reduction in porosity occurs concomitantly with surface 539 diffusion and dissolution-precipitation processes that heal and seal the pore space, leaving D R A F T July 10, 2017, 12:11pm D R A F T where ρ s is the density of the solid, ρ f is the density of the fluid, n is the Eulerian porosity, 574 v s is the velocity of the solid and v f is the velocity of the fluid. In the above equations, 575 r denotes the rate at which fluid mass is generated from solid mass. This will be our 576 definition for the reaction rate. Neglecting gradients in ρ s and ρ f , the combination of 577 equations A1 and A2 leads to We can relate the divergence of the relative fluid velocity with respect to the solid to 579 the gradient in fluid pressure by using Darcy's law [Coussy, 2004]: where p f is the fluid pressure, k is the permeability of the material, and η is the viscosity 581 of the fluid. The divergence of the velocity of the solid is the bulk volumetric (Eulerian) 582 strain rate: where is the bulk volumetric strain. The combination of relations A4 and A5 with 584 equation A3 yields: The variation of the fluid density can expressed as where c f is the compressibility of the fluid. The variation of the solid density is decomposed We now need to express the evolution of the average density of the solid as a function 594 of the reaction progress. We are interested in the following type of chemical reaction: where ν i,f are stoichiometric coefficients. Denoting ξ the reaction progress and m 0 d the 596 total mass of fluid that can be released by the reaction (per unit volume of rock), the 597 solid volume is expressed as where M i is the molar mass of constituent i. The average solid mass is denoted m s . The 599 density of the solid is ρ s = m s /V s ; hence we have: The conservation of mass imposes that ∂m s /∂m = −1, so that The last term in parenthesis of the previous equation corresponds to the solid volume 602 change of the reaction, which we denote ∆ r V s . Keeping in mind that m = m 0 d ξ, we can D R A F T July 10, 2017, 12:11pm Finally, we have to keep in mind that V s = 1 − n by definition, so that equation A13 The combination of the relation above with the mass balance equation A9, eventually 607 leads to the governing equation (1) for pore fluid pressure.

A2. Constitutive behaviour
We assume that the dehydrating rock is elasto-plastic. We introduce the yield function 609 f (σ, ζ), where σ is the stress tensor and ζ is an internal variable on which the yield 610 cap depends. ζ can be identified to either the finite porosity of the material, or more 611 directly to the reaction progress in the case of a pure chemical control over the material's 612 strength. These two options will be discussed later on. We also assume, in accordance 613 with experimental observations, that the material undergoes strain hardening. For the 614 sake of simplicity, f and g were assumed linear in terms of τ . In that case, the consistency 615 condition for plastic loading is therefore where h is the hardening modulus and dλ is a positive infinitesimal scalar (so-called 617 plastic increment). If we now introduce the plastic potential g(σ), the elasto-plastic stress summation convention is adopted and the indices i, j take values 1, 2, 3. In the (p , τ ) 630 space, the stress vector σ is defined as: Likewise, the strain vector can be described by the volumetric strain and shear strain γ: The elastic tensor is written then: D R A F T July 10, 2017, 12:11pm D R A F T where K is the bulk modulus of the porous material, and G its shear modulus. For a 634 general plastic behaviour, the derivatives of the yield function and plastic potential are 635 expressed as follows: where µ is the friction coefficient, and β the dilatancy factor. For the sake of simplicity, 637 f and g were assumed linear in terms of τ . The full expression for the incremental stress-638 strain relation in the (p , τ ) space becomes: In expressions A25 and A26, the factor ∂f /∂ζ corresponds to the dependency of the

648
In such a framework, we can write Using this expression for ζ into Equations (A25) and (A26), we finally arrive at the incremental constitutive formulation of Equations (3) and (4) in the main text.
D R A F T July 10, 2017, 12:11pm D R A F T

A3. Uniaxial compaction
The combination of relations 10 and 3 yields Assuming a constant σ n over time, equation 9 yields which we combine to equation A28, making use of relation 4, to obtain and Now we can use equation 1 to express the volumetric strain rate: We finally use the expression of the volumetric strain rate given by A33 into equation A30

658
to obtain Using the chemical kinetics established in Equation 7, we finally arrive at Equation (11) 660 of the main text.

661
Appendix B: Volume change for antigorite dehydration D R A F T July 10, 2017, 12:11pm where ρ • is the density under standard conditions (at T = T 0 = 25 • C), and The density as a function of pressure is given by where is the linear strain calculated from the bulk modulus K and its derivative with 665 pressure K = ∂K/∂P : The total change in density as a function of pressure and temperature is finally obtained The inequality (28) can be rewritten as D R A F T July 10, 2017, 12:11pm D R A F T Assuming that the material is nominally stable, i.e, h > h crit , and considering that the 671 reaction has a total negative volume change (1/ρ f + ∆ r V s < 0), the rhs of Inequality (C1) 672 is a positive quantity. If the term in brackets on the lhs is negative, f would have to be 673 also negative in order to satisfy the inequality. This is a contradiction since the material 674 is porosity-softening and f > 0. So a first requirement for the instability to be possible 675 is that the bracketed term is positive, which implies that Considering that the material is compactant (β < 0) and that the solid volume change of 677 the reaction is negative (∆ r V s < 0), Inequality (C2) is satisfied when either or where β crit is defined in Equation (29) of the main text. In that case, the instability 680 criterion in terms of f is given by (see Equation (C1)) which is exactly the same as Equation (30) after some rearrangements.

C2. Linear analysis
Denoting p f the small perturbation of p f above p eq , we rewrite the governing equation 683 for pore pressure (11) as follows: where A is the amplitude of the perturbation, S is its growth factor, and λ is its wave-691 length. Our suggested solution must be consistent with the prescribed boundary con-692 ditions, hence we require that λ = W/k, (k = 1, 2, . . .). The perturbation is unstable 693 (S > 0) if (1) f > f crit and (2) the wavelength is greater than a critical wavelength λ crit : After some algebra, the critical wavelength can be rewritten as in (31) in the main text.

695
For n r > 1, the situation is mathematically more complicated. Indeed, one can imme-696 diately observe in equation (C6) that the reaction term does not appear in a first order 697 stability analysis (it is elevated to a power greater than 1). However, one could hope to 698 make useful analytical predictions without resorting to a full numerical treatment. Let us 699 assume that the perturbation has a characteristic amplitude A and a characteristic length 700 scale L, i.e., p f (y, t) = A(t)g(y/L(t)) where g is a non dimensional function of the order 701 of 1. The governing equation for p f is then D R A F T July 10, 2017, 12:11pm The condition given in (C10) is very similar to the one obtained in (C8) for the linear 704 case, but here we see that the amplitude of the perturbation A appears in the definition 705 of the critical length scale. Hence, according to (C10), the system is unstable only if the 706 wavelength and the amplitude of the perturbations are large enough to overcome diffusion.

707
In any case, we note that the reaction will only dominate the system at early times; one can 708 show that diffusion cannot be neglected everywhere when the system evolves with time.

709
Here we are only interested in the behaviour at early times, because we have assumed The numerical solution of the fully coupled, nonlinear system is obtained by discretising the governing equation for pore pressure (11) in space using a centered finite difference stencil, and then solving for reaction progress, pore pressure, volumetric strain, total mean stress, porosity and shear stress as a coupled system of ordinary differential equations (ODEs). In practice, we normalise the governing equations by using the magnitude of the imposed total normal stress σ n as the stress scale, the reaction rate 1/r 0 as the time scale, and the thickness of the antigorite layer (denoted L). We use a centered finite difference approximation of the second-order spatial derivatives of pore pressure, with a grid defined by points y i = i∆y, and implement periodic boundary conditions at the edges y = 0 and D R A F T July 10, 2017, 12:11pm D R A F T y/L = 1. The full system of ODEs is then: where subscripts i indicate variables at point x i , and In addition, a consistency check is performed by computing the total normal stress: and verifying a posteriori that it remains constant throughout space and time. All       Figure 4. Stability boundaries for a dehydrating material with a porosity-dependent yield function. The dilatancy factor β has been assumed equal to the friction coefficient µ (i.e., we used the so-called associated plasticity assumption) in order to limit the number of free parameters. Black curves correspond to stability boundary for the pore pressure runaway (condition 30), and blue curves correspond to stability boundary for conventional (mechanical) compaction instability (condition 24).  In this simulation, the total normal stress is σ n = 3.71 GPa, so that the initial effective mean stress is 0.61 GPa and the initial shear stress is 0.09 GPa.
D R A F T July 10, 2017, 12:11pm D R A F T