Nicolas Brantut
email: n.brantut@ucl.ac.uk

Ioannis Stefanou

Jean Sulem

Dehydration-induced instabilities at intermediate 1 depths in subduction zones 2

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 localisation, over characteristic length scales determined by the hydraulic diffusivity, the elastoplastic parameters of the rock, and the reaction rate. Typical lower bounds for the localisation length are of the order of 10 to 100 m for antigorite dehydration and deformation at 3 GPa. The fluid pressure and deformation instability is expected to induce stress build-up in the surrounding rocks forming the subducted slab, which provides a mechanism for the nucleation and propagation of intermediate-depth earthquakes.

Introduction

During prograde metamorphism in subduction zones, hydrous phases such as serpentines progressively dehydrate, forming free fluid phases at depth. Such dehydration reactions are systematically associated with a net decrease in solid volume (the reactions forming solid products denser than reactants), and with a variable change in fluid volume, the sign of which being controlled by the pressure and temperature conditions at which the reaction occurs. From initially nonporous metamorphic rocks (such as antigorite serpentinite), dehydration reactions therefore produce, at least transiently, a porous rock saturated with fluids, the rheology of which is markedly different from the original rock [e.g. [START_REF] Rutter | Dehydration and deformation of intact cylinders of serpentinite[END_REF]. The occurrence of metamorphic dehydration reactions has therefore a great impact on the stress/strain state in subduction zones.

One key specific impact of dehydration reactions is their potential to trigger unstable faulting and earthquakes, a phenomenon generally termed "dehydration embrittlement". This phenomenon corresponds to the transition from ductile to brittle deformation due to a dehydration-induced increase in pore fluid pressure. It is often thought that dehydrationembrittlement is one of the main causes of, or is at least linked to, intermediate-depth earthquakes in subduction zones [e.g., Hacker et al., 2003b]. Dehydration embrittlement has been observed experimentally [e.g., [START_REF] Raleigh | Experimental deformation of serpentinite and its tectonic implications[END_REF][START_REF] Murrell | The effect of decomposition of hydrous minerals on the mechanical properties of rocks[END_REF] and is well explained theoretically when the reaction produces an excess fluid volume, i.e., typically at relatively low pressure conditions (e.g., less than around 2.5 GPa in antigorite). Under those conditions, the excess fluid volume generated by the reaction tends to increase the pore fluid pressure, reducing the effective stress, and therefore D R A F T July 10, 2017, 12:11pm D R A F T bringing back the material into the brittle field. At higher pressure, where the total volume change of the dehydration reaction is negative, laboratory experiments indicate that dehydration embrittlement and earthquake instability might still occur [e.g. [START_REF] Jung | Intermediate depth earthquake faulting by dehydration embrittlement with negative volume change[END_REF], but the exact mechanism remains unclear. One possibility is that the reaction products are plastically weak and facilitate brittle deformation in their surroundings [e.g. [START_REF] Rutter | Dehydration and deformation of intact cylinders of serpentinite[END_REF][START_REF] Brantut | Strain localisation and slip instability in a strain-rate hardening, chemically weakening material[END_REF]. Another potential mechanism is that the porosity generated by the reaction is rapidly compacted, hence producing a pore pressure rise which could bring back the material into the brittle regime.

Coupling between mechanical compaction and dehydration reactions has been investigated in detail in the case of viscous rock rheology: in his seminal study, Connolly

[1997] (followed by [START_REF] Connolly | Compaction-driven fluid flow in viscoelastic rock[END_REF] and [START_REF] Connolly | Fluid flow in compressive tectonic settings: Implications for midcrustal seismic reflectors and downward fluid migration[END_REF], summarised recently in [START_REF] Connolly | The mechanics of metamorphic fluid expulsion[END_REF]) used a viscous compaction rheology coupled to devolatilisation reactions together with a power-law relation between porosity and permeability, and determined that compaction would drive pore fluid pressure up to near-lithostatic values, while producing intermittent upward motions of fluid (so-called porosity-waves). Despite the great success of this model for the prediction of fluid extraction from the lower crust (or along subduction zones, see [START_REF] Skarbek | Dehydration-induced porosity waves and episodic tremor and slip[END_REF]), one key assumption is that the compaction behaviour is essentially driven by a viscous creep process, and does not include the instantaneous response of the material.

Serpentinites, and more specifically the high pressure form, antigorite, are known to behave in a semi-brittle manner even at high pressure and temperature [e.g. [START_REF] Chernak | Deformation of antigorite serpentinite at high temperature and pressure[END_REF][START_REF] Proctor | Ductile to brittle" transition in thermally stable antigorite gouge at mantle pressures[END_REF], i.e., antigorite deformation systematically involves a significant degree of microcracking, and its behaviour at high pressure is similar D R A F T July 10, 2017, 12:11pm D R A F T X -5

to cataclastic flow. Furthermore, laboratory experiments [START_REF] Rutter | Dehydration and deformation of intact cylinders of serpentinite[END_REF] show that dehydrated and partially dehydrated serpentinite behave essentially like classical porous rocks, and that concepts of porous rock mechanics can be used to describe their mechanical behaviour. Therefore, a purely viscous constitutive law may not capture all the features and potential instabilities associated with coupled compaction and dehydration in antigorite.

Here, we model coupled dehydration and deformation in antigorite using as a first approximation a time-independent inelastic flow law which includes strain hardening, strain-dependent dilatancy/compaction and a porosity-dependent yield envelope. Our approach is based on the concepts typically used to model the behaviour of porous rocks [e.g Rudnicki and Rice, 1975;[START_REF] Issen | Conditions for compaction bands in porous rock[END_REF][START_REF] Wong | The brittle-ductile transition in porous rock: a review[END_REF][START_REF] Stefanou | Chemically induced compaction bands: Triggering conditions and band thickness[END_REF], and dehydration has here an indirect effect by contributing to the overall change in porosity and fluid pressure. In this framework, two types of instabilities can arise: a rate-independent bifurcation related to the constitutive behaviour of the rock, and a reaction-driven, rate-dependent instability due to the growth to small pore pressure perturbations. We specifically focus on the behaviour of serpentinite under conditions such that the dehydration reaction produces a negative total volume change, with the aim of determining whether pore pressure instabilities can occur.

Model and governing equations

2.1. Fluid mass balance and pore pressure change

During dehydration, antigorite becomes a porous aggregate, with a porosity n filled with water at a pressure denoted p f . We consider that the porosity of the aggregate is connected (at least at the scale of interest, here of the order of 10 to 100 metres), and that D R A F T July 10, 2017, 12:11pm D R A F T the fluid flows through the rock according to Darcy's law with a permeability k. The fluid pressure is modified by two independent contributions: the elasto-plastic compaction of the rock (bulk volumetric strain , taken negative in compression), and the generation of fluids from the dehydration reaction. These assumptions lead to the following governing equation for pore pressure (see full derivation in Appendix A1):

∂p f ∂t = k ηc b ∇ 2 p f - 1 c b ∂ ∂t + m 0 d (1 + ρ f ∆ r V s ) ρ f c b ∂ξ ∂t , (1)
where m 0 d is the total fluid mass that can be released by the reaction per unit rock volume, ρ f is the fluid density, ∆ r V s is the solid volume change of the reaction, ∂ξ/∂t is the reaction rate, and

c b = (1 -n)c s + nc f (2)
is an effective compressibility combining the compressibility of the solid c s , the compressibility of the fluid c f , and the porosity n. Note here that c b is not the usual storage capacity, because we did not split the volumetric strain rate into an elastic and plastic one.

Rheology

It is well established experimentally that antigorite aggregates undergo a brittle to ductile transition at confining pressures of the order of 300 to 400 MPa [START_REF] Escartín | Nondilatant brittle deformation of serpentinites: Implications for Mohr-Coulomb theory and the strength of faults[END_REF],

and that this transition depends weakly on temperature (within antigorite's stability field).

Near the dehydration temperature of antigorite, the ductile behaviour remains dominated by cataclastic flow even at mantle pressures [START_REF] Chernak | Deformation of antigorite serpentinite at high temperature and pressure[END_REF][START_REF] Proctor | Ductile to brittle" transition in thermally stable antigorite gouge at mantle pressures[END_REF], which is most likely due to the strong [001] cleavage plane and the insufficient number of independent slip systems in antigorite single crystals. Across the stability boundary

D R A F T July 10, 2017, 12:11pm D R A F T
of antigorite, a very significant volume change occurs: the solid volume tends to decrease by up to around 25%, generating a porosity occupied by pressurised water. Therefore, at pressure and temperature conditions near its stability boundary, antigorite is expected to behave very much like a ductile (cataclastic) porous rock. Such a behaviour has been well documented in lizardite by [START_REF] Rutter | Dehydration and deformation of intact cylinders of serpentinite[END_REF], and we assume here that the same behaviour applies to antigorite.

In order to describe the elasto-plastic behaviour of antigorite, we introduce the yield function f (σ, ζ), where σ is the stress tensor (here we follow the sign convention of continuum mechanics and take compressive stresses as negative) and ζ is an internal variable on which the yield cap depends. The yield function limits the elastic domain in the stress space (f < 0, see Figure 1). It is assumed that inelastic strain increments are generated when the stress state lies on the yield surface ( f = 0) and if loading is taking place. ζ can be identified to either the porosity of the material, or more directly to the reaction progress in the case of a pure chemical control over the material's strength [START_REF] Stefanou | Chemically induced compaction bands: Triggering conditions and band thickness[END_REF][START_REF] Sulem | Thermal and chemical effects in shear and compaction bands[END_REF]. While either option could be deemed acceptable in the light of the available experimental data from [START_REF] Rutter | Dehydration and deformation of intact cylinders of serpentinite[END_REF], we will develop our model assuming that the primary control on the rock's strength is given by its porosity (a robust observation in porous rocks, see [START_REF] Wong | The brittle-ductile transition in porous rock: a review[END_REF]). We therefore equate incremental changes in the internal variable ζ to irreversible (inelastic) porosity changes. We also assume, in accordance with experimental observations, that the material undergoes strain hardening. The incremental constitutive behaviour resulting from our D R A F T July 10, 2017, 12:11pm D R A F T assumptions is written as follows (see full derivation in Appendix A2):

dp = GK [(1 + (h -f β)/G)d -βdγ] h -f β + G + βµK + f βKm 0 d ∆ r V s h -f β + G + βµK dξ, (3)
dτ = GK -µd + (h -f β)/K + βµ dγ h -f β + G + βµK + f Gm 0 d ∆ r V s h -f β + G + βµK dξ, (4)
where p is the Terzaghi effective mean stress (p = p + p f , where p is the mean stress), τ is the shear stress (taken equal to the square root of the second invariant of the deviatoric stress tensor), K and G are the bulk and shear elastic moduli of the rock, respectively, γ is the shear strain (taken equal to the square root of the second invariant of the deviatoric strain tensor), h is the strain hardening coefficient, β is the dilatancy factor and f = ∂f /∂ζ is the dependency of the yield cap on porosity. We observe in equations 3 and 4 that the effect of the variation of f with porosity (terms in f ) on the mechanical behaviour is entirely captured by the modified hardening modulus h -f β. The chemical coupling appears through the product f ∆ r V s : the reaction has only an indirect effect, which is to modify the porosity. The solid volume change is always negative (porosity creation), and f is positive (the yield surface shrinks with increasing porosity), so that the overall effect of the reaction is to weaken the material.

Reaction rate

A very general formulation of mineral reaction rates is given by [START_REF] Lasaga | Fluid flow and chemical reaction kinetics in metamorphic systems[END_REF]]

∂ξ ∂t = κA rlm s|∆G| nr , ( 5
)
where κ is the temperature-dependent kinetic constant (typically following an Arrhenius law), A rlm is the specific surface area of the rate-limiting mineral, ∆G is the Gibbs energy change of the reaction, s is the opposite of the sign of ∆G, and n r is the order of the reaction. Under isothermal conditions, and for small departures from equilibrium, we can D R A F T July 10, 2017, 12:11pm D R A F T expand ∆G in terms of pore pressure; only retaining the leading order term, we have [START_REF] Wang | Effects of reaction kinetics and fluid drainage on the development of pore pressure excess in a dehydrating system[END_REF] ∆G ≈ c (p f -p eq ),

where c = ∂∆G/∂p f and p eq is the pore pressure at equilibrium. Following [START_REF] Wang | Effects of reaction kinetics and fluid drainage on the development of pore pressure excess in a dehydrating system[END_REF], we can rewrite 5 as

∂ξ ∂t = sr 0 A rlm A 0 rlm |1 -p f /p eq | nr , (7)
where r 0 is a reference reaction rate, and A 0 rlm is the surface area of the rate limiting mineral at the reference rate. The change in surface area is not very well constrained by experimental data. Hence, for simplicity, in the following we will consider the ratio

A rlm /A 0 rlm ≈ 1.
This simplification is valid only when the reaction progress is small, i.e., when depletion of the reaction is negligible (ξ 1).

Geometry and stress equilibrium

We consider a simple system made of a uniform horizontal layer of antigorite, sufficiently extended so that lateral strains can be neglected (i.e., the system is invariant in the plane of the layer). In this geometry, the vertical stress is given by p -2τ / √ 3, and stress equilibrium requires that:

∂ ∂y p - 2 √ 3 τ = 0. ( 8
)
where y denotes the vertical coordinate. The boundary condition is a constant applied vertical stress, which implies that

∂p ∂t - 2 √ 3 ∂τ ∂t = 0. ( 9
)
In such a geometry, only the vertical strain component is nonzero. Because no lateral deformation is allowed, the vertical strain is equal to the volumetric strain. Therefore, D R A F T July 10, 2017, 12:11pm D R A F T the shear and volumetric strains are related to one another by:

γ + 2 √ 3 = 0. ( 10
)
Using relations ( 9) and (10) in the incremental constitutive formulation (Equations 3, 4), and combining with the fluid mass conservation (Equation 1) and expression (7) for the reaction rate, we arrive at a single, nonlinear diffusion equation that governs the fluid pressure (see details in Appendix A3):

∂p f ∂t = M k/η 1 + c b M ∂ 2 p f ∂y 2 + M m 0 d (1/ρ f + ∆ r V s ) -f X 1 + c b M sr 0 1 - p f p eq nr , (11)
where

M = GK h -f β + G + βµK 1 + 2 √ 3 µ 1 + 2 √ 3 β + (h -f β) 1 G + 4 3K , (12)
and

X = -m 0 d ∆ r V s βK -2G/ √ 3 h -f β + G + βµK . ( 13
)
3. Parameters

Yield function

The model described above contains a number of parameters that ought to be constrained from experimental data. The yield cap can be constrained from the extensive dataset of [START_REF] Rutter | Dehydration and deformation of intact cylinders of serpentinite[END_REF] on intact and dehydrated blocks of lizardite. The data used are shown in Figure 2. For simplicity, we use a modified Cam-clay yield surface, given by

f (p , τ ) = τ -C (b + p )(p * -p ) = 0, ( 14
)
where C is the critical state line ratio, b is the tensile strength, and p * is the compaction yield pressure (following the notation of [START_REF] Wong | The brittle-ductile transition in porous rock: a review[END_REF]). Specific values of C, b and p * for intact, partially and fully dehydrated serpentinites are reported in Figure 2.

D R A F T July 10, 2017, 12:11pm D R A F T

X -11

In accordance with observations on porous sandstones, the critical compaction pressure p * decreases with increasing porosity. Zhang et al. [1990] proposed a grain crushing model in which p * scales with the porosity n as p * ∝ n -3/2 . For the serpentinite samples dehydrated at 35 MPa effective pressure and deformed at room temperature (squares in Figure 2), we find a reasonable fit with p * = 130 MPa, while the inferred porosity of the sample was around 19%. The partially dehydrated samples (diamonds in Figure 2) had a porosity of around 4%, and we determine a p * of 380 MPa. The relationship between p * and n for these two sample types does not seem compatible with the scaling proposed by Zhang et al. [1990]. Here we will use an empirical scaling p * ∝ 1/n [START_REF] Rutter | The deformation of porous sandstones; are Byerlee friction and the critical state line equivalent?[END_REF], with a constant of proportionality equal to 19 MPa. The corresponding yield caps are reported as dotted lines in Figure 2.

Obviously, the yield cap described by Equation 14 is not appropriate as the porosity approaches zero, since in that case p * diverges. However, the focus of this work is the description of rocks that are already undergoing dehydration, i.e., in which the porosity is never exactly zero. Furthermore, in absence of a more complete dataset on dehydrated and partially dehydrated serpentinite, the Cam-clay yield surface is one of the simplest yield criterion which is closed at high pressure (i.e., the material can fail by pure hydrostatic compaction). Hence, our choice for the yield function should be viewed as a first order approximation which incorporates the essential qualitative elements of the behaviour of dehydrating serpentinite: a yield cap that is closed at high pressure and that shrinks with increasing porosity.

D R A F T July 10, 2017, 12:11pm D R A F T

Mechanical and hydraulic parameters

The elastic properties of serpentinite can be obtained from Voigt-Reuss-Hill averages of the single crystal properties, and are given by [START_REF] Bezacier | Elasticity of antigorite, seismic detection of serpentinites, and anisotropy in subduction zones[END_REF]: K = 67.9 GPa, G = 38.5 GPa. The average Poisson's ratio is hence ν = 0.26. Note that we are modelling dehydrating serpentinite, and hence the average elastic properties of the rock should be made dependent upon the evolving rock mineralogical composition and porosity. However, we focus here on the initiation of the dehydration reaction, and hence expect that the pure antigorite end-member is a good approximation to the overall properties of the rock at the beginning of dehydration.

The friction coefficient µ is given by the local slope of the yield envelope (equation A24):

µ = ∂f ∂p = C 2p + b -p * 2 (p + b)(p * -p ) . (15)
In the framework of associated plasticity, we could assume that the dilatancy factor is merely equal to the friction coefficient. However, it is well known that rocks do not follow associated flow rules, and hence we shall leave the dilatancy factor β as a free parameter, and explore how it influences the stability of compaction in our model. Likewise, we will leave the hardening coefficient as a free parameter, in order to encompass the widest possible range of behaviours.

The permeability of the rock is expected to vary as a function of porosity, and hence be impacted by the compaction of the rock. However, these second order controls on permeability should only influence the behaviour of the material at large times, and not the initiation of the instability. For instance, using power-law permeability-porosity relationships, Connolly [1997] has shown that a dehydrating rock can generate travelling porosity pulses, a well known feature of nonlinear parabolic equations. Since we want

D R A F T July 10, 2017, 12:11pm D R A F T X -13
to focus on the initiation of the instability, we shall assume a constant value for the permeability, keeping in mind that this assumption should be relaxed when modelling the long-term evolution of the system. The permeability of dehydrating serpentinite was measured by [START_REF] Tenthorey | Reaction-enhanced permeability during serpentinite dehydration[END_REF], who report k = 10 -22 m 2 for the intact material and k = 10 -20 to 10 -18 m 2 during dehydration at 600 • C and 700 • C, respectively. Here we choose k = 10 -20 m 2 as a representative value for serpentinites which are dehydrating not far from equilibrium.

The viscosity of the fluid η can be precisely determined from interpolation of experimental measurements. We will use here the formulation of the International Association for the Properties of Water and Steam, 2008 (http://iapws.org/relguide/visc.pdf), and compute the appropriate viscosity at the target pressure and temperature conditions. The resulting viscosity of water ranges from η = 6.9 × 10 -5 to 4.9 × 10 -7 Pa s at p f = 2 GPa and 600 • C and p f = 5 GPa and 700 • C, respectively.

The effective compressibility of the rock, c b , does not play any role in the stability of the system; it is only required for the computation of the full numerical solutions, and acts as a scaling factor for the pore pressure rate. Here we compute c b from Equation (2), assuming that c s = 1/K.

Chemical parameters

The parameters associated with the chemical reaction can be obtained from the ther- The properties of each mineral can be extracted from thermodynamic tables [e.g. [START_REF] Holland | An internally consistent thermodynamic data set for phases of petrological interest[END_REF]Hacker et al., 2003a], and are reported in Table 1. The thermodynamic properties of water are determined as a function of pressure and temperature from the IAPWS formulation 95.

modynamic
The computation of the solid volume change ∆ r V s requires the knowledge of the molar volumes of each solid phase at the pressure and temperature conditions of each reaction.

The molar volumes are computed following the approach explained in Hacker et al. [2003a],

which is recalled in Appendix B for completeness.

The phase boundaries and net volume change of each reaction (16,17,18) are shown as a function of pressure and temperature in Figure 3. Average volume changes along each phase boundary are presented in Table 2.

In addition, the knowledge of the molar volume of antigorite allows to compute precisely 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.

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

The reaction kinetics of antigorite as a function of pressure is not well constrained by existing experimental data, which typically focus on the effect of temperature. However, the formulation 5 is general, and hence kinetic parameters obtained from experiments in which ∆G is imposed from a temperature over-or under-step should also be valid in the

D R A F T July 10, 2017, 12:11pm D R A F T X -15
case when ∆G changes due to pressure fluctuations. Here we extract kinetic parameters from the study of [START_REF] Eggler | Rate of antigorite dehydration at 2 GPa applied to subduction zones[END_REF], in which the dehydration kinetics of antigorite was determined at 2 GPa (i.e., for reaction 17) as a function of temperature. At a temperature T , the rate of antigorite dehydration is given as

rate = κ A rlm s ∆G RT nr in mol atg /cm 3 rock /s, ( 19
)
where R is the gas constant. In terms of reaction progress ξ, we rewrite equation 19 as

∂ξ ∂t = sr 0 A rlm A 0 rlm |1 -p f /p eq | nr , (20)
where

r 0 = V m A 0 rlm κ p eq c RT nr . ( 21
)
The molar volume of antigorite V m can be taken as the average along the phase boundary between antigorite and forsterite, enstatite and water, using equation B6. This yields V m ≈ 1.73 × 10 -3 m 3 /mol atg . The initial specific surface area of antigorite, A 0 rlm , depends on the grain size and shape; for square prisms with width W and length L, A 0 rlm = (2W 2 + 4LW )/LW 2 . Using W = 50 µm and L = 10 µm, we find A 0 rlm = 2.8 × 10 5 m 2 /m 3 .

The rate constant κ given by Eggler and Ehmann [2010] is 9.2 × 10 -11 mol atg /m 2 /s. The coefficient c , as defined in equation 6, is the net volume change of the reaction per unit mole of antigorite:

c = ν f M f (1/ρ f + ∆ r V s ). ( 22
)
Using the parameter value reported in Table 2 for the reaction of antigorite into forsterite and enstatite, we obtain c ≈ 8.1 × 10 -5 m 3 /mol atg . The equilibrium pressure p eq and the appropriate temperature T can be found from the phase boundary (see Figure 3). As a representative value, we choose p eq = 3 GPa and [2010] report that the exponent n r is equal to one, i.e., the kinetic is linear. Combining all the above parameters into our lumped kinetic parameter r 0 , we find a representative value of r 0 ≈ 1.47 × 10 -6 /s. (23)

Stability analysis

The system governed by Equation ( 11) involves two coupled phenomena: mechanical deformation (compaction), through the parameter M , and metamorphic effects with fluid production and porosity creation through the source term in (11). In this Section, we detail the possible sources of mechanical and chemical instabilities, and derive the key stability conditions in terms of the model parameters.

Mechanical instability

If M < 0, we immediately observe that Equation 11 is a diffusion equation with a negative diffusivity. This corresponds to an unstable system as non uniformities become more localised rather than more diffuse with increasing time. In fact, the condition M < 0 is strictly equivalent to the compaction localisation condition derived by Issen and Rudnicki [2000, their equation 19] for axisymmetric compression:

h -f β G < h crit G = - 3K 3K + 4G 1 + 2 √ 3 µ 1 + 2 √ 3 β , (24)
with the additional term -f β occurring here due to the explicit dependence of the yield cap on porosity. For a yield cap shrinking with increasing porosity, as expected and observed in dehydrating serpentinites (Figure 2), f is positive, and β is negative (due to shear-enhanced compaction at high pressure). Hence, the critical hardening modulus is reduced when the yield cap is assumed to be directly dependent on porosity; this merely D R A F T July 10, 2017, 12:11pm D R A F T reflects that the yield cap expansion due to shear hardening is in fact offset by a further yield cap expansion (resp. shrinkage) due to shear-induced compaction (resp. dilatancy).

Note, in passing, that we have derived here a slightly more general case for the compaction instability including pore fluid pressurisation effects; the stability condition, without chemical effects, is in fact M/(1 + c b M ) < 0. This has no practical consequence on the criterion because c b is usually very small.

Chemical instability and overpressure development

The case of interest here is when M > 0 (Equation ( 12)), i.e., a mechanically stable case. Equation 11 is a nonlinear diffusion equation. We examine here the case when the total volume change of the reaction is negative, for which the reaction kinetics increases when p f increases beyond p eq . Before delving into the full analysis including the effect of pore pressure diffusion, it is instructive to first analyse the undrained case; the governing equation for pore pressure reduces to

∂p f ∂t = M m 0 d (1/ρ f + ∆ r V s ) -f X 1 + c b M r 0 p f p eq -1 nr . ( 25
)
Assuming constant parameters, this equation has an analytical solution, which is

p f (t) = p eq + (p 0 f -p eq ) exp M m 0 d (1/ρ f + ∆ r V s ) -f X p eq (1 + c b M ) r 0 t (26)
if n r = 1, and

p f (t) = p eq + p -nr eq (1 -n r )r 0 M m 0 d (1/ρ f + ∆ r V s ) -f X 1 + c b M t + p 0 f -p eq 1-nr 1/(1-nr) (27)
if n r > 1, where p 0 f is the initial pore pressure in the system (a small perturbation above the equilibrium pressure).

D R A F T July 10, 2017, 12:11pm D R A F T

In both cases (linear and nonlinear kinetics), the evolution of pore pressure is an unbounded growth if

M m 0 d (1/ρ f + ∆ r V s ) -f X > 0. ( 28
)
For linear kinetics, the growth is exponential, while for nonlinear kinetics the growth corresponds to a finite time blow-up. In practice, this distinction is unimportant since pore pressure diffusion, as well as other nonlinearities not accounted for in our simplified system (such as the depletion of the reactant or the change in mechanical properties with the evolving deformation and mineralogy of the rock), are expected to strongly change the evolution of pore pressure at large times. Despite these subtleties, the solution of the undrained problem yields a key condition (Inequality 28) to observe a potential pore pressure runaway. This condition for instability can be expressed as a pair of conditions in terms of the dilatancy factor β and the dependence of the yield cap on porosity f (see Appendix C1 for details):

β < β crit = -( √ 3/2) 3 2 1 -ν 1 -2ν 1 ρ f ∆ r V s + 1 -1 , ( 29
)
f > f crit = √ 3 1 -ν 1 -2ν h -h crit 1 -β/β crit 1 + 1 ρ f ∆ r V s . ( 30
)
Using the numerical values detailed in the previous section, we remark that β crit is always positive. Hence, for the cases of interest where β < 0 (shear enhanced compaction), we always have β < β crit , and thus the condition for instability is simply f > f crit . stability between these two boundaries. For h = 0 (which seems relatively plausible for antigorite, see [START_REF] Escartín | Nondilatant brittle deformation of serpentinites: Implications for Mohr-Coulomb theory and the strength of faults[END_REF]), the stability boundaries do not overlap; the area of stability expands with decreasing (negative) values of β. Because the friction coefficient µ only appears in the expression for h crit , changing µ independently from β (which is a reasonable choice for most rocks) has a moderate impact on the stability diagram for the pore pressure instability (condition (30)), but introduces large changes for the bifurcation criterion (condition ( 24)).

The stability boundaries (h =

Using our assumed form for the yield function f , we note that f diverges as porosity decreases towards zero: in practice, f is potentially very large at the initiation of the dehydration reaction (when porosity is very low). If the applied stress is such that β < 0, as in the case at high effective pressure, then the reaction-driven instability is very probable at the onset of dehydration. This implies that the pore fluid pressure will quickly rise to decrease the effective stress; the pressure (solid and fluid) will hence equilibrate very suddenly, generating a transient pore pressure pulse.

Strain localisation and pore pressure build-up

The considerations above are restrained to the undrained system. The diffusion of pore fluids will tend to stabilise the pore pressure runaway (if it occurs) by draining the rock over a certain length scale. There is no general analytical solution available for equations of the form (11); here we restrict our analysis to a study of the stability of the system to small departures from equilibrium (which corresponds to p f = p eq ).

As stated in Section 3.3, the experimental data of [START_REF] Eggler | Rate of antigorite dehydration at 2 GPa applied to subduction zones[END_REF] show that the near-equilibrium kinetics is linear, n r = 1. In that case, the stability analysis detailed in Appendix C2 shows that pore pressure runaways are possible if f > f crit (same as D R A F T July 10, 2017, 12:11pm D R A F T condition 30) and if the spatial wavelength of the perturbation λ is such that

λ > λ crit = 2π p eq (k/η)[1 -f β/(h -h crit )] r 0 m 0 d (1/ρ f + ∆ r V s )(1 -f /f crit ) . ( 31
)
We observe that the critical wavelength depends on a number of parameters, some of them well constrained (equilibrium pressure p eq , permeability k, fluid viscosity η, and the net volume change of the reaction), and others much more poorly known (essentially, all the parameters associated with the mechanical behaviour, including f and β). However, as discussed above, the value of f at the onset of the reaction is expected to be very large since small increments in porosity have large effects on the yield cap when the rock is initially non porous. We can use this fact to our advantage by noticing that the critical wavelength λ crit tends to a constant, nonzero value for f f crit :

λ crit ∼ 2π p eq (k/η) r 0 m 0 d ∆ r V s √ 3 ν -1 2ν -1 1 β -1 -β -1 crit . (32)
Equation ( 32) provides a simple lower bound for the critical wavelength, which is, quite remarkably, independent from the hardening modulus h.

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

Based on the parameter values outlined above for antigorite, we estimate typical values for λ crit of the order of 10 to 100 m. In our model, this length scale corresponds to the characteristic width over which pore pressure builds up, and compaction (negative volumetric strain) localises.

D R A F T July 10, 2017, 12:11pm D R A F T

Numerical tests

We performed a series of numerical computations in order to explore further the effects of potential nonlinearities associated with the variations of mechanical, hydraulic and chemical parameters during deformation and reaction. We included a power-law dependence of the permeability on porosity, k ∝ ζ 3 (where we recall that ζ is the porosity of the rock minus its poroelastic variations), and accounted for depletion of antigorite by using a simple first-order approximation for the reactant surface area,

A rlm /A 0 rlm ≈ (1 -ξ)
, in Equation ( 7). The numerical method is described fully in Appendix D, and the Matlab R implementation and source code is available online at http://github.com/nbrantut/Compaction_Dehydration. The solution for pore fluid pressure, strain, reaction progress and porosity is computed within a layer of width L, with periodic boundary conditions, and an initial sinusoidal infinitesimal pore pressure perturbation is added to the homogeneous initial conditions. We chose a representative example by using an initial pore pressure of 3 GPa, an initial total mean stress of 3.61 GPa (i.e., an initial shear stress of 0.87 GPa), and an initial porosity of 3 %. Using the initial, reference parameters, we find that the critical wavelength for instability is λ crit ≈ 0.08L, so that we expect some compaction localisation (at least transiently).

The volumetric strain profile within the layer is shown as a function of time in Figure 6 on the edge (y = ±L/2). The peak compaction at the centre develops quite rapidly, initially accelerating, and then develops over a timescale of the order of t × r 0 ≈ 10, and then stabilises. At the edges of the layer, the compaction develops more slowly and a strong strain gradient develops in the initial phase of the instability; once the compaction stabilises in the centre (at around t × r 0 ≈ 20), the strain becomes more homogeneous as the edges also compact further until the whole process eventually stabilises, leaving only a slight strain heterogeneity near the centre.

The evolution of all other key variables as a function of time is shown in Figure 7.

The onset of the compaction instability is marked by a rapid acceleration of the reaction progress, porosity and pore pressure changes localised in the centre of the layer (solid lines, t×r 0 between 10 and 15). As the reaction progresses further the differences between the reaction progress, porosity, strain and pore pressure at the centre and at the edges progressively decrease. As the reaction approaches completion, the porosity stabilises at around 27% while the pore pressure remains very high, corresponding to a mean effective stress of around 0.07 GPa. At this point, no further compaction is possible because the periodic boundary conditions effectively ensure that the system is undrained.

The full stress path of each material element in the layer is drawn in Figure 8 in the effective mean stress, shear stress (p , τ ) space. As the pore pressure and porosity increase in the layer, the yield surface shrinks and the stress state evolves to maintain mechanical equilibrium. The stress paths of different elements in the layer are not exactly the same, since heterogeneities in pore pressure and stress develop between the centre (solid line)

and the edges (dashed line). As the reaction approaches completion, the effective stress becomes very small (p /σ n ≈ 0.019, i.e., p ≈ 0.07 GPa) and the material cannot compact

D R A F T July 10, 2017, 12:11pm D R A F T X -23
further because of the residual (small) nonzero strength, illustrated by the final yield surface.

Overall, the numerical simulation confirms that a compaction instability is possible, but that it is only transient as the system tends again to a homogeneous state while the reaction approaches completion. The typical timescale for instability is of the same order of magnitude as 1/r 0 , i.e., around 7 to 8 days. In our simulations, we impose periodic boundary conditions so that the layer is effectively undrained: hence, the porosity and pore pressure remain high after the reaction is completed. We performed complementary simulations using drained boundaries as another end-member case scenario, and observed that there is a long-term compaction occurring after the initial instability, and the pore pressure returns progressively towards its equilibrium value over timescales determined by the drained length across the layer.

Discussion

Model assumptions and limitations

The key assumption of the model presented here is that ductile deformation of antigorite is essentially time-independent. This approximation is justified by experimental observations of antigorite deformation at elevated pressures and temperature, showing that cataclastic mechanisms typically dominate [START_REF] Chernak | Deformation of antigorite serpentinite at high temperature and pressure[END_REF][START_REF] Amiguet | Deformation mechanisms and rheology of serpentines in experiments and in nature[END_REF][START_REF] Auzende | Deformation mechanics of antigorite serpentinite at subduction zone conditions from experimentally and naturally deformed rocks[END_REF] and that fully plastic flow is unlikely to occur due to the large crystal anisotropy and lack of available slip systems in antigorite. Furthermore, partially dehydrated serpentinite has also been shown to deform very similarly to porous sandstones [START_REF] Arkwright | Role of porosity and dehydration reactions on the deformation of hot-pressed serpentinite aggregates[END_REF][START_REF] Rutter | Dehydration and deformation of intact cylinders of serpentinite[END_REF], which motivates the use of a closed yield envelope.

D R A F T July 10, 2017, 12:11pm D R A F T However, there are also clear experimental indications that ductile flow of antigorite depends on strain rate [e.g. [START_REF] Hilairet | High-pressure creep of serpentine, interseismic deformation, and initiation of subduction[END_REF][START_REF] Amiguet | Creep of phyllosilicates at the onset of plate tectonics[END_REF], and time-dependency may not be negligible if deformation occurs over very long timescales. We can test whether time-dependent plastic flow contributes significantly to deformation by estimating the strain rates developing during the instability in our model. As mentioned in the previous Section, the characteristic timescale over which the instability develops is determined by the reaction kinetics, 1/r 0 . The typical strain achieved during the instability is of the order of a few percents, so that the strain rate is of the order of r 0 /100, which is around 10 -8 s -1 . Using that strain rate, the typical shear flow stress extrapolated from the plastic flow laws given in Amiguet et al. [2012, their figure 6] is of the order of 0.1 GPa. Therefore, our model based on time-independent ductile flow is broadly consistent with the rheology of serpentinites as long as relatively low shear stresses are considered (around 0.1 GPa).

However, under near-isostatic conditions, the shear stress is not expected to be large enough to produce significant viscous flow within the timescale of instability. The contribution of viscous flow to isostatic compaction typically scales with the inverse of the porosity [e.g. [START_REF] Wilkinson | Pressure sintering by power law creep[END_REF][START_REF] Mckenzie | The generation and compaction of partially molten rock[END_REF], so that the driving effective mean stress should be of the order of 1 GPa to achieve strain rates of around 10 -8 s -1 for a porosity of around 10%. It is therefore likely that viscous compaction in the absence of shear stresses slows down the development of the pore pressure instability.

A number of coupled deformation-dehydration physical models have been developed

based on time-dependent rheology of rocks [e.g. Connolly, 1997;[START_REF] Connolly | Compaction-driven fluid flow in viscoelastic rock[END_REF][START_REF] Connolly | The mechanics of metamorphic fluid expulsion[END_REF][START_REF] Skarbek | Dehydration-induced porosity waves and episodic tremor and slip[END_REF]. One key parameter exhibited by these models is the compaction length scale, which is related to the hydraulic diffu-D R A F T July 10, 2017, 12:11pm D R A F T sivity, compressibility and bulk viscosity of the rock. In this context, the timescale for compaction is entirely determined by the rheology of the rock. In our approach using time-independent ductile deformation, a characteristic length for pore pressure development and compaction (λ c , Equation ( 31)) also arises, associated with the hydraulic diffusivity, rheology and reaction rate. This length λ c is clearly analogue to the one obtained from viscous compaction, and corresponds to the limiting case where the timescale for compaction is determined by the reaction rate. For typical crustal metamorphic reactions, Connolly [1997] determines that the minimum compaction length is of the order of 10 to 100 m. Quite interestingly, we find that the minimum compaction length in our model is also of the order of 10 to 100 m for antigorite dehydration at intermediate depth (see Figure 5). This similarity in compaction length scale between the two types of model essentially arises from the similarity in timescales between viscous flow and reaction rate.

Here we made the assumption that the material deforms under isothermal conditions.

Although a full analysis including thermal effects is beyond the scope of the present work, we discuss here qualitatively how changes in temperature can arise and modify the behaviour of the material. Firstly, the dehydration reaction of antigorite is endothermic, which constitutes a significant heat sink. Secondly, the irreversible work done by inelastic deformation corresponds to energy dissipation and is a heat source. At 4 GPa and 625 • C, the enthalpy change of the reaction is 442 kJ per mole of antigorite (computed from the database of [START_REF] Holland | An internally consistent thermodynamic data set for phases of petrological interest[END_REF]), i.e., around 2.5 × 10 2 MJ m -3 . For a reaction rate of the order of 10 -8 s -1 , the rate of heat absorbed by the reaction is of the order of 2.5 J m -3 s -1 . The rate of work of the volumetric strain is given by ˙ p , and for an effective stress p of around 100 MPa and a characteristic strain rate of the order of 10 -8 s -1 , the D R A F T July 10, 2017, 12:11pm D R A F T rate of heat generated by inelastic deformation is of the order of 1 J m -3 s -1 . Although the exact quantities are only approximate, the rate of heat absorbed by the reaction and released by deformation are of the same order of magnitude, and might compensate each other. If significant changes in temperature occur, for instance a cooling due to the endothermic character of the reaction overcoming the heat generated by deformation, the reaction kinetics will also be modified according to the corresponding change in ∆G (Equation ( 5)). In that case, heat flow across the deformed/reacted zone might be the rate limiting process. Assuming a heat diffusivity of 10 -6 m 2 s -1 , the characteristic diffusion time across a layer of 100 m in thickness is 10 10 s, and it drops to 10 8 s for a 10 m width layer. Therefore, we expect heat flow to be limiting only when the deformation (and reaction) rate becomes faster than 10 -10 s -1 (respectively, 10 -8 s -1 ) in a thick (respectively, thin) layer. In our simulations, such strain rates are achieved transiently during the instability, so that thermal effects might affect the behaviour of the system only after the instability has initiated.

Despite the limitations outlined above, the model formulation is quite general and could be applicable to most devolatilisation reactions in subduction zones but also at shallower depths in the crust. This is especially relevant for relatively cold parts of the crust (or subduction zones), where the viscous creep rates of rocks are slow compared to reaction kinetics, so that the hypothesis of time-independent rheology would be justified. Not all reactions are expected to generate a pore pressure instability, and the surprising result here is that instabilities only arise when the total volume change of the reaction is negative (i.e., at high pressure for antigorite). Under shallow crustal conditions, most devolatilisation reactions are expected to produce a positive volume change, so instaiblities should be D R A F T July 10, 2017, 12:11pm D R A F T investigated on a case-by-case basis. In any case, the governing equation for pore pressure (Equation ( 11)) remains valid and can be used to make predictions for fluid flow in active prograde metamorphic settings, such as deep sedimentary basins (where gypsum and clay mineral dehydrate), around rising plutons and magma chambers, and of course along subduction zones.

Effective stress in subduction zones and implications for intermediatedepth earthquakes

One of the key outcome of our model is that antigorite dehydration at intermediate depth results in a very rapid build-up of pore pressure, reducing the effective mean stress towards near-zero values while maintaining a significant open porosity (at least as long as the fluids are trapped inside the dehydrating layer). This pore pressure build-up occurs despite the net negative volume change of the reaction, and is primarily driven by the collapse of porosity. Hence, we expect the effective mean stress to remain near zero throughout the regions where dehydration proceeds, independently from the net volume change associated with the reaction. While this is also expected in models involving viscous rock deformation, we find here that the path towards low effective stress states is unstable, but only in the case when the net volume change is negative 1/ρ f + ∆ r V s < 0.

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

One interesting outcome of our model is that we show that the dehydration and compaction process produces significant shear deformation in the rock, and not just pure D R A F T July 10, 2017, 12:11pm D R A F T volumetric compaction. In our approach, we used uniaxial strain boundary conditions in order to simulate a simple, tractable problem. In nature, the stress state and boundary conditions are necessarily more complex. The existence of shear stresses, even very small (0.01 to 0.1 GPa), is expected to produce significant shear strains during the pore pressure build-up. The dehydrating body of antigorite then acts as a very deformable layer or inclusion, amplifying the stresses around it. In addition, pore fluid diffusion outside the dehydrating body also contributes to decreasing the effective stress in the surrounding rocks.

The combination of relatively rapid stress amplification and effective pressure reduction provides a reasonable mechanism for the inception of brittle deformation in the surrounding ultramafic rocks (peridotites and metagabbros). Indeed, at the temperatures and pressures relevant to antigorite dehydration and at strain rates as high as 10 -8 s -1 during the pore pressure instability, peridotites have a high strength and are unlikely to accomodate deformation in a purely viscous manner, thus raising elastic stresses and favoring brittle deformation. Such a mechanism is essentially a kind of dehydration embrittlement, but in the surrounding rocks and not in the dehydrating serpentinite itself. While this process had been suggested in earlier works by [START_REF] Kirby | Localized polymorphic phase transformations in high-pressure faults and applications to the physical mechanism of deep earthquakes[END_REF] or [START_REF] Rutter | Dehydration and deformation of intact cylinders of serpentinite[END_REF],

our model provides first-order quantitative constrains on its likelyhood by establishing a closed-form stability criterion that depends on the rock rheology and reaction kinetics (Equation ( 28)).

As a final note of caution regarding our interpretations in terms of fluid pressure at intermediate depths, we recall here that the concept of effective stress makes sense only if a uniform pore pressure can be defined for a representative volume element of the rock.

D R A F T July 10, 2017, 12:11pm D R A F T X -29
This is most likely the case when porosity is larger than several percents, above the percolation threshold [START_REF] Guéguen | Introduction to the Physics of Rocks[END_REF]. However, over long timescales, the progressive drainage of the pore fluid outside of the dehydrating zones tends to allow compaction to reduce porosity; this reduction in porosity occurs concomitantly with surface diffusion and dissolution-precipitation processes that heal and seal the pore space, leaving a disconnected pore network and free fluids present only as fluid inclusions [e.g. [START_REF] Smith | Diffusional crack healing in quartz[END_REF]. In the long term, pore pressure is not a well defined concept and neither is effective stress. There are currently insufficient constrains on healing mechanisms in silicates to draw definitive quantitative conclusions on the persistence and connectivity of the pore space at intermediate depths in subduction zones, but we expect that the short timescale of the instability (typically of the order of 10 days in our simulations) and the large porosity generated by the reaction (around 20%) ensure a reasonable pore network connectivity and validates the use of the concept of effective stress.

Conclusions

We developed a model to simulate coupled deformation and dehydration of antigorite at intermediate depths in subduction zones. Our model shows that dehydration can lead to unstable pore pressure rise and deformation when the net volume change of the reaction is negative, due to a positive feedback between pore fluid pressure, compaction and dehydration rate. The stability criterion (Equation ( 30)) is controlled by the dependence of the yield envelope on porosity (parameter f = ∂f /∂ζ); using estimates for f consistent with well established mechanics of porous rocks, we find that antigorite dehydration leads to unstable deformation under typical intermediate-depth conditions in subduction zones.

D R A F T July 10, 2017, 12:11pm D R A F T Furthermore, we also show that the instability is associated with localised deformation and fluid pressure over a characteristic length scale controlled by hydraulic diffusivity, rheological parameters and reaction rate. A lower bound for this characteristic length is of the order of 10 to 100 m, commensurate with the viscous compaction length scale obtained in models using time-dependent rheologies.

Our model predicts that the typical strain rates during the instability are of the order of 10 -8 s -1 . At such rates and at the relatively cool dehydration temperature of antigorite (around 600 • C), the surrounding, chemically stable peridotites and metagabbroic rocks forming the subducted slab have a high strength, and are expected to build up elastic stresses. The elevated pore pressures associated with the dehydration reaction can therefore transiently bring the surrounding rocks back into the brittle field, thus allowing the nucleation and propagation of earthquakes. This mechanism is clearly a type of dehydration embrittlement, but is crucially based on a rapid stress transfer between weak dehydrating rocks and strong surrounding ultramafic rocks, in a manner conceptually similar to the early model of [START_REF] Kirby | Localized polymorphic phase transformations in high-pressure faults and applications to the physical mechanism of deep earthquakes[END_REF] for transformation-induced instabilities.

Appendix A: Derivation of governing equations A1. Fluid pressure For a reacting porous medium, the continuity equations for the solid skeleton and the fluid are [Coussy, 2004]

∂ ρ s (1 -n) ∂t + div ρ s (1 -n)v s = -r, (A1) ∂(nρ f ) ∂t + div(nρ f v f ) = +r, (A2) D R A F T July 10, 2017, 12:11pm D R A F T X -31
where ρ s is the density of the solid, ρ f is the density of the fluid, n is the Eulerian porosity, v s is the velocity of the solid and v f is the velocity of the fluid. In the above equations, r denotes the rate at which fluid mass is generated from solid mass. This will be our definition for the reaction rate. Neglecting gradients in ρ s and ρ f , the combination of equations A1 and A2 leads to

1 -n ρ s ∂ρ s ∂t + n ρ f ∂ρ f ∂t + div(v s ) + div n(v f -v s ) = r(1/ρ s -1/ρ f ). (A3)
We can relate the divergence of the relative fluid velocity with respect to the solid to the gradient in fluid pressure by using Darcy's law [Coussy, 2004]:

n(v f -v s ) = -(k/η)grad(p f ), ( A4
)
where p f is the fluid pressure, k is the permeability of the material, and η is the viscosity of the fluid. The divergence of the velocity of the solid is the bulk volumetric (Eulerian) strain rate:

div(v s ) = ∂ ∂t , (A5)
where is the bulk volumetric strain. The combination of relations A4 and A5 with equation A3 yields:

1 -n ρ s ∂ρ s ∂t + n ρ f ∂ρ f ∂t + ∂ ∂t - k η ∇ 2 p f = r(1/ρ s -1/ρ f ). (A6)
The variation of the fluid density can expressed as

1 ρ f ∂ρ f ∂t = c f ∂p f ∂t , ( A7
)
where c f is the compressibility of the fluid. The variation of the solid density is decomposed into two contributions: (

1 ρ s ∂ρ s ∂t = c s ∂p s ∂t + 1 ρ s ∂ρ s ∂m r, (
1 -n)c s + nc f ∂p f ∂t = k η ∇ 2 p f - ∂ ∂t + r 1 ρ f - 1 ρ s - 1 -n ρ s ∂ρ s ∂m . ( A9
)
We now need to express the evolution of the average density of the solid as a function of the reaction progress. We are interested in the following type of chemical reaction:

mineral 0 -→ i ν i mineral i + ν f fluid, ( A10
)
where ν i,f are stoichiometric coefficients. Denoting ξ the reaction progress and m 0 d the total mass of fluid that can be released by the reaction (per unit volume of rock), the solid volume is expressed as

V s = m 0 d ν f M f M 0 (1 -ξ)/ρ 0 + i (ν i M i /ρ i )ξ , (A11)
where M i is the molar mass of constituent i. The average solid mass is denoted m s . The density of the solid is ρ s = m s /V s ; hence we have:

1 ρ s ∂ρ s ∂m = 1 m s ∂m s ∂m - 1 V s ∂V s ∂m . (A12)
The conservation of mass imposes that ∂m s /∂m = -1, so that

1 ρ s ∂ρ s ∂m = 1 V s - 1 ρ s - ∂V s ∂m . ( A13

A3. Uniaxial compaction

The combination of relations 10 and 3 yields

∂p ∂t + ∂p f ∂t = GK (1 + (h -f β)/G) + 2β/ √ 3 h -f β + G + βµK ∂ ∂t - f βK h -f β + G + βµK ∂ξ ∂t . (A28)
Assuming a constant σ n over time, equation 9 yields

∂p ∂t = 2 √ 3 ∂τ ∂t , (A29)
which we combine to equation A28, making use of relation 4, to obtain

∂p f ∂t = M ∂ ∂t -f X ∂ξ ∂t , (A30)
where

M = GK h -f β + G + βµK 1 + 2 √ 3 µ 1 + 2 √ 3 β + (h -f β) 1 G + 4 3K , (A31)
and

X = -m 0 d ∆ r V s βK -2G/ √ 3 h -f β + G + βµK . ( A32
)
Now we can use equation 1 to express the volumetric strain rate:

∂ ∂t = k η ∂ 2 p f ∂y 2 + m 0 d (1/ρ f + ∆ r V s ) ∂ξ ∂t -c b ∂p f ∂t . (A33)
We finally use the expression of the volumetric strain rate given by A33 into equation A30 to obtain

∂p f ∂t = M k/η 1 + c b M ∂ 2 p f ∂y 2 + M m 0 d (1/ρ f + ∆ r V s ) -f X 1 + c b M ∂ξ ∂t . ( A34
)
Using the chemical kinetics established in Equation 7, we finally arrive at Equation ( 11) of the main text.

Assuming that the material is nominally stable, i.e, h > h crit , and considering that the reaction has a total negative volume change (1/ρ f + ∆ r V s < 0), the rhs of Inequality (C1) is a positive quantity. If the term in brackets on the lhs is negative, f would have to be also negative in order to satisfy the inequality. This is a contradiction since the material is porosity-softening and f > 0. So a first requirement for the instability to be possible is that the bracketed term is positive, which implies that

β 1 ρ f + 4G 3K + 4G ∆ r V s < 2 √ 3 -3G 3K + 4G ∆ r V s . (C2)
Considering that the material is compactant (β < 0) and that the solid volume change of the reaction is negative (∆ r V s < 0), Inequality (C2) is satisfied when either

1 ρ f + 4G 3K + 4G ∆ r V s > 0, (C3) or 1 ρ f + 4G 3K + 4G ∆ r V s < 0 and β > β crit , (C4)
where β crit is defined in Equation ( 29) of the main text. In that case, the instability criterion in terms of f is given by (see Equation (C1))

f > -(K + 4G/3)(h -h crit )(1/ρ f + ∆ r V s ) ∆ r V s (βK -2G/ √ 3) -β(K + 4G/3)(1/ρ f + ∆ r V s ) , (C5)
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 for pore pressure (11) as follows:

∂p f ∂t = M k/η 1 + c b M ∂ 2 p f ∂y 2 + M m 0 d (1/ρ f + ∆ r V s ) -f X 1 + c b M r 0 p f p eq nr . (C6) D R A F T July 10, 2017, 12:11pm D R A F T X -39
The qualitative behaviour of the system depends on the value of n r .

For n r = 1 (i.e., the dehydration reaction kinetics is approximated to be first-order), the equation is linear and we can perform a linear stability analysis to explore how the system evolves. We assume no flux conditions at the boundaries of a dehydrating layer of thickness W (which simulates a serpentinite layer embedded in an impermeable host rock). Then we look for solutions of (C6) (with n r = 1) in the form

p f = A cos(2πy/λ) exp(St), ( C7
)
where A is the amplitude of the perturbation, S is its growth factor, and λ is its wavelength. Our suggested solution must be consistent with the prescribed boundary conditions, hence we require that λ = W/k, (k = 1, 2, . . .). The perturbation is unstable

(S > 0) if (1) f > f crit and
(2) the wavelength is greater than a critical wavelength λ crit :

λ > λ crit = 2π p eq k/η r 0 [m 0 d (1/ρ f + ∆ r V s ) -f X/M ] . (C8)
After some algebra, the critical wavelength can be rewritten as in (31) in the main text.

For n r > 1, the situation is mathematically more complicated. Indeed, one can immediately observe in equation (C6) that the reaction term does not appear in a first order stability analysis (it is elevated to a power greater than 1). However, one could hope to make useful analytical predictions without resorting to a full numerical treatment. Let us assume that the perturbation has a characteristic amplitude A and a characteristic length scale L, i.e., p f (y, t) = A(t)g(y/L(t)) where g is a non dimensional function of the order of 1. The governing equation for p f is then y/L = 1. The full system of ODEs is then:

∂p f ∂t = A L 2 M k/η 1 + c b M g (y/L) + M m 0 d (1/ρ f + ∆ r V s ) -f X 1 + c b M r 0 A nr g nr (y/L). ( C9
∂ξ i ∂t = s(1 -ξ i )|p fi /p eq -1| nr , (D1) ∂p fi ∂t = M 1 + c b M 1 η∆y 2 k i+1/2 p fi+1 -p fi -k i-1/2 p fi -p fi-1 + M m 0 d (1/ρ f + ∆ r V s ) -f X 1 + c b M ∂ξ i ∂t , (D2)
∂ i ∂t = 1 M ∂p fi ∂t + f X ∂ξ i ∂t , ( D3
)

∂p i ∂t = 1 √ 3/2 + β -f /K (f /K -β) ∂p fi ∂t -f ∂ i ∂t + f m 0 d ∆ r V s ∂ξ i ∂t , ( D4
)

∂ζ i ∂t = -m 0 d ∆ r V s ∂ξ i ∂t + ∂ i ∂t - 1 K ∂p i ∂t + ∂p fi ∂t , ( D5
)

∂τ i ∂t = -β ∂p i ∂t + ∂p fi ∂t -f ∂ξ i ∂t , (D6)
where subscripts i indicate variables at point x i , and

k i±1/2 = (k i + k i±1 )/2. (D7)
In addition, a consistency check is performed by computing the total normal stress:

∂σ ni ∂t = ∂p i ∂t - 2 √ 3 ∂τ i ∂t , (

  properties of the mineral (and fluid) involved. At elevated pressure and D R A F T July 10, 2017, 12:11pm D R A F T temperature, there are three distinct dehydration reactions involving antigorite: antigorite → 4 talc + 18 forsterite + 27 water,

  h crit + f β and f = f crit ) are shown as function of f /G and β in Figure 4, where we have assumed β = µ. For h > 0, the mechanical compaction bifurcation arises only for negative values of f , while the reaction-driven compaction and pore pressure instability occurs only for positive values of f , leaving an area of D R A F T July 10, 2017, 12:11pm D R A F T

  (a). The initial ( (y, 0) = 0) and final (at t × r 0 = 100) profiles are highlighted in black, and intermediate stages are shown light grey. Over time, a net volumetric strain localisation develops around the centre of the layer, and remains there permanently at large times. The time evolution of the compaction localisation instability is better observed in Figure 6(b), which shows the volumetric strain in the centre of the layer (y = 0) and D R A F T July 10, 2017, 12:11pm D R A F T

  A8) D R A F T July 10, 2017, 12:11pm D R A F T where c s is the compressibility of the constituents of the solid skeleton, and m denotes the mass of fluid released by the chemical reaction. The last term in A8 corresponds to the evolution of the average density of the solid skeleton as the reaction proceeds (this term would be zero if the solid was transforming entirely into a fluid, without generating solid products). Combining relations A7 and A8 into equation A6 yields:

)

  The last term in parenthesis of the previous equation corresponds to the solid volume 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 D R A F T

Figure 1 .Figure 2 .

 12 Figure 1. Schematic of the yield surface in the stress space (p , τ ). The yield surface is

Figure 3 .

 3 Figure 3. Net volume change 1/ρ f + ∆ r V s , in m 3 per unit mass of water, as a function

Figure 4 .Figure 5 .-Figure 6 .

 456 Figure 4. Stability boundaries for a dehydrating material with a porosity-dependent

Figure 7 .Figure 8 .

 78 Figure 7. Volumetric strain (black), reaction progress ξ (orange), porosity ζ (green)

Table 1 .

 1 Thermodynamic properties of phases involved in the dehydration reactions of antigorite. Data from[START_REF] Holland | An internally consistent thermodynamic data set for phases of petrological interest[END_REF].

		Molar weight Molar volume Thermal expansion Bulk modulus
		M	V • m	α •	K •	∂K/∂P
	Phase	g/mol	cm 3 /mol	×10 -5 K -1	GPa
	Atg, antigorite	4536	1754.7	4.7	67.9	2.77
	Talc	379.7	136.4	3.7	41.6	6.5
	Fo, forsterite	140.7	43.7	6.1	127	5.37
	Ens, enstatite	200.8	62.6	5.1	106	8.5
	PhA, phase A	456.3	154.4	8.3	97.4	6.0
	D R A F T		July 10, 2017, 12:11pm		D R A F T

Table 2 .

 2 Average fluid mass and volume change associated with the dehydration reactions of antigorite.

	Mass of releasable fluid	Volume change
	m 0 d	1/ρ f + ∆ r V s

X -33

use relation A11 to express ∆ r V s as

Finally, we have to keep in mind that V s = 1 -n by definition, so that equation A13 becomes

The combination of the relation above with the mass balance equation A9, eventually 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

, where σ is the stress tensor and ζ is an internal variable on which the yield cap depends. ζ can be identified to either the finite porosity of the material, or more directly to the reaction progress in the case of a pure chemical control over the material's strength. These two options will be discussed later on. We also assume, in accordance with experimental observations, that the material undergoes strain hardening. For the sake of simplicity, f and g were assumed linear in terms of τ . In that case, the consistency condition for plastic loading is therefore

where h is the hardening modulus and dλ is a positive infinitesimal scalar (so-called plastic increment). If we now introduce the plastic potential g(σ), the elasto-plastic stress increment is then given by where M el is the elastic tensor and is the total strain tensor. The combination of equations A16 and A17 allows the determination of the plastic increment dλ and yields the full incremental constitutive relation

where

and

We further assume that the material is isotropic. The scalars p and τ are used here and they represent, respectively, the Terzaghi effective mean stress (i.e. the difference between the total mean stress and the pore pressure, p = tr(σ)/3 + p f ) and the shearing stress intensity. The shearing stress intensity τ is defined as the square root of the second invariant of the deviatoric part, s, of the stress tensor: τ = 1 2 s ij s ij . The Einstein summation convention is adopted and the indices i, j take values 1, 2, 3. In the (p , τ ) 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:

where K is the bulk modulus of the porous material, and G its shear modulus. For a general plastic behaviour, the derivatives of the yield function and plastic potential are expressed as follows:

where µ is the friction coefficient, and β the dilatancy factor. For the sake of simplicity, f and g were assumed linear in terms of τ . The full expression for the incremental stressstrain relation in the (p , τ ) space becomes: et al., 1990] that P * is scaled by the product of the grain size and the porosity of the rock, to the power 3/2. This dependency of f on P * , and of P * on porosity implies that the yield cap can be considered as a function of the porosity of the rock. Hence, a natural choice for the parameter ζ is the total nominal porosity (i.e., the current finite porosity resulting from deformation and reaction, minus any variations produced by elastic deformations).

In such a framework, we can write The density of a phase as a function of temperature is given by

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 pressure K = ∂K/∂P :

The total change in density as a function of pressure and temperature is finally obtained from

which implies that the molar volume is

Appendix C: Stability analysis

C1. Derivation of criterion

The inequality (28) can be rewritten as

The reaction term will be dominant, i.e., the system will be unstable, if

The condition given in (C10) is very similar to the one obtained in (C8) for the linear case, but here we see that the amplitude of the perturbation A appears in the definition of the critical length scale. Hence, according to (C10), the system is unstable only if the wavelength and the amplitude of the perturbations are large enough to overcome diffusion.

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

Here we are only interested in the behaviour at early times, because we have assumed that the parameters of the equations are constant. For further evolution of the system, the full nonlinearities should be included and it is not worth going too far mathematically with our simplified system.

Appendix D: Numerical methods

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


Tian Ma
email: matian56@sina.com

Shouhong Wang
email: showang@indiana.edu

DYNAMICAL THEORY OF THERMODYNAMICAL PHASE TRANSITIONS

Keywords: July 14, 2017 dynamical equation, thermodynamic systems, principle of dynamical transitions, thermodynamical potential, phase diagram, transition diagram, dynamical diagram, Ehrenfest classification

The objective of this paper is establish the basic dynamical theory for thermodynamic phase transitions. The main components of the theory include 1) dynamical equations of thermodynamical systems, 2) general principle of dynamical transitions, and 3) three basic theorems of thermodynamic phase transitions.

Introduction

Phase transition is a universal phenomena in most, if not all, natural systems, and refers to the transformation of the system from one state to another, as the control parameter crosses certain critical threshold. In physics, a state often refers to a "stable" solution of the mathematical model.

Most problems in the natural sciences are described by either dissipative or conservative systems. A systematic dynamic transition theory for dissipative systems has been developed by the authors through a sequence of paper starting around 2004, and the theory and its wide range of applications are synthesized in the authors' book [START_REF]Phase Transition Dynamics[END_REF]. The theory demonstrates a general principles that phase transitions of all dissipative systems can be classified into three categories: continuous, catastrophic, and random.

The goals of phase transition theories are [START_REF] Ehrenfest | Phasenumwandlungen im ueblichen und erweiterten Sinn, classifiziert nach den entsprechenden Singularitaeten des thermodynamischen Potentiales[END_REF] to determine the definition and types of phase transitions, [START_REF] Ma | Bifurcation theory and applications[END_REF] to derive the critical parameters, (3) to obtain the critical exponents, and (4) to understand the mechanism and properties of phase transitions, such as supercooled and superheated states and the Andrews critical points, etc.

There are two routes for studying phase transitions: one is Landau's approach using thermodynamic potentials, and the other is a microscopic approach using statistical theory. These routes are complementing to each other. Our dynamical transition theory follows the Landau's route, but is not a mean field theory approach. It is based, however, on the thermodynamical potential.

The objective of this paper is to provide a fundamental theory of thermodynamic phase transitions based on the dynamic transition theory and its various applications. The main ingredients of the paper are as follows.

First, for a thermodynamic system, there are three different levels of physical quantities: control parameters λ, the order parameters (state functions) u = (u 1 , • • • , u N ), and the thermodynamic potential F . These are well-defined physical quantities, fully describing the system. The potential is a functional of the order parameters, and is used to represent the thermodynamic state of the system.

In a recent paper [START_REF]Dynamic law physical motion and potential-descending principle[END_REF], we postulated a potential-descending principle (PDP) for statistical physics, which gives rise to a dynamic equation for thermodynamical systems; see (4) and (5) in Section 2 below. Then for a given basic equilibrium state ū, the deviation from the basic state ū satisfies a dissipative dynamical system given by [START_REF] Stanley | Introduction to Phase Transitions and Critical Phenomena[END_REF], which dictates the phase transition dynamical behavior of the system.

Second, with the potential functional F , the basic state ū, and the dynamical equation ( 6) at our disposal, the dynamical transition theory can then be applied, and we obtain in this paper three theorems, providing a full theoretical characterization of thermodynamic phase transitions.

Third, the first theorem, Theorem 1, states that as soon as the linear instability occurs, the system always undergoes a dynamical transition, to one of the three types: continuous, catastrophic and random. This theorem offers the detailed information for the phase diagram and the critical threshold λ c in the control parameter λ

= (λ 1 , • • • , λ N ) ∈ R N .
They are precisely determined by the following equation:

(1)

β 1 (λ 1 , • • • , λ N ) = 0.
Here β 1 is the first eigenvalue of the linear operator L λ given by ( 2)

L λ w ≡ - ∂ 2 F (ū; λ) ∂u i , ∂u j w = β 1 (λ)w,
where F is the potential functional.

Fourth, the second theorem, Theorem 2, provides the corresponding relationship between the Ehrenfest classification and the dynamical classification. The Ehrenfest classification identifies experimental quantities to determine the transition types, while the dynamical classification provides a full theoretical characterization of the transition.

This theorem establishes a natural bridge between theory and experiments. With the corresponding relationship given in this theorem, we obtain a precise and easy theoretical approach to determine the transition order, the critical exponents, and the transition diagram, which are encoded in the following transition equation:

(3)

du λ dt = β 1 (λ)u λ + a(λ)u 2 λ + b(λ)u 3 λ + h.o.t.
We remark here that the detailed information on critical exponents is entirely encoded in the potential functionals, and again, our approach is not a mean-field theory.

Also, the theorem shows that there are only first, second and thirdorder thermodynamic transitions. The third-order transition can hardly be determined by thermodynamic parameters experimentally. The dynamical transition theory, however, offers an easy and clear approach to identify the third-order transition.

Fifth, the last theorem, Theorem 3, states that both catastrophic and random transitions lead to saddle-node bifurcations, and both latent heat, superheated and supercooled states always accompany the saddle-node bifurcation associated with the first order transitions.

Theorem 2 with the transition equation 3 provides the needed information on the dynamic phase diagrams; this is important for understanding the dynamical behavior of transition states.

We further emphasize that the three theorems lead to three important diagrams: the phase diagram, the transition diagram and the dynamical diagram. These diagrams appear only to be derivable by the dynamical transition theory presented in this paper. In addition, our theory achieves the fourth goal of phase transition theories as stated in the beginning of this Introduction, which is hardly achievable by other existing theories.

The paper is organized as follows. Section 2 establishes the dynamical equation for thermodynamic equations, and Section 3 recalls the general principle for dynamical transitions. Section 4 describes the Ehrenfest classification. Section 5 states and proves the the three basic theorems for thermodynamic phase transitions.

Dynamic Equations of Statistical Physics

For a given thermodynamic system, the order parameters (state functions) u = (u 1 , • • • , u N ), the control parameters λ, and the thermodynamic potential F are well-defined quantities, fully describing the system. The potential is a functional of the order parameters, and is used to represent the thermodynamic state of the system. There are four commonly used thermodynamic potentials: the internal energy, the Helmholtz free energy, the Gibbs free energy, and the enthalpy.

In a recent paper [START_REF]Dynamic law physical motion and potential-descending principle[END_REF], we postulated the following potential-descending principle (PDP) for statistical physics:

Principle 1 (Potential-Descending Principle [START_REF]Dynamic law physical motion and potential-descending principle[END_REF]). For each thermodynamic system, there are order parameters u = (u 1 , • • • , u N ), control parameters λ, and the thermodynamic potential functional F (u; λ). For a non-equilibrium state u(t; u 0 ) of the system with initial state u(0, u 0 ) = u 0 , we have the following properties:

1) the potential F (u(t; u 0 ); λ) is decreasing:

d dt F (u(t; u 0 ); λ) < 0 ∀t > 0;
2) the order parameters u(t; u 0 ) have a limit

lim t→∞ u(t; u 0 ) = ū;
3) there is an open and dense set O of initial data in the space of state functions, such that for any u 0 ∈ O, the corresponding ū is a minimum of F , which is called an equilibrium of the thermodynamic system:

δF (ū; λ) = 0.
We have shown that PDP is a more fundamental principle than the first and second laws, and provides the first principle for describing irreversibility, leads all three distributions: the Maxwell-Boltzmann distribution, the Fermi-Dirac distribution and the Bose-Einstein distribution in statistical physics. Consequently, the potential-descending principle is the first principle of statistical physics.

Also importantly, based on PDP, the dynamic equation of a thermodynamic system in a non-equilibrium state takes the form du dt = -AδF (u, λ) for isolated systems, (4)

       du dt = -AδF (u, λ) + B(u, λ), AδF (u, λ) • B(u, λ) = 0 for coupled systems, ( 5
)
where δ is the derivative operator, B represents coupling operators, and A is a symmetric and positive definite matrix of coefficients. We refer interested readers to [START_REF]Dynamic law physical motion and potential-descending principle[END_REF] for details.

Phase transition is a universal phenomena in most, if not all, natural systems. In fact, the central problem in statistical physics is on phase transitions. A phase transition refers to the transition of the system from one state to another, as the control parameter crosses certain critical threshold. In physics, a state often refers to a "stable" solution of the mathematical model, and in statistical physics, a state refers to an equilibrium state. Consequently, in thermodynamic phase transitions, we can consider a basic equilibrium state ū, and study the dynamic transitions of this state ū. If we write u = ū + u , then the dynamic law (4) or (5) leads to an equation for the deviation order parameter u . After suppressing the primes, we arrive at the following dynamic equation for the deviation order parameter u:

(6) du dt = L λ u + G(u, λ).
Again u is the order parameter, λ is the control parameter of the system, L λ is a linear operator, and G(u, λ) is the nonlinear operator. For all thermodynamic phase transitions, ( 6) is a dissipative dynamical system, and u = 0 represents the basic equilibrium state ū.

Principle of Phase Transition Dynamics

A comprehensive and unified dynamic transition theory for dissipative systems have been developed recently by the authors with applications to a wide range of problems in nonlinear sciences; see refer [START_REF]Phase Transition Dynamics[END_REF] and the references therein. In connection to applications, the dynamic transition theory provides a systematic approach for classifying and determining the detailed information of the transition.

A thorough examination of the dynamic transition theory and its applications enable us to postulate the following general principle of phase transition dynamics.

Principle 2 (Principle of Phase Transition Dynamics [START_REF]Phase Transition Dynamics[END_REF]). Phase transitions of all dissipative systems can be classified into three categories:

(1) continuous transition if lim λ→λ 0 u λ = ū;

(2) catastrophic transition if lim λ→λ 0 u λ = ū; and

(3) random transition if both lim λ→λ 0 u λ = ū and lim λ→λ 0 u λ = ū occur.

Here ū is the basic state, and u λ are the transition states (physically, transition states correspond to local attractors).

A few remarks are now in order. First, the validity of the above principle is ensured by the dynamic transition theorem [4, Theorem 2.1.3]. In addition, all applications of the dynamic transition theory so far are in agreement with this principle. In fact, such principle is also valid for transition associated with pattern formations.

Second, this is a universal principle, applicable to phase transitions of all dissipative systems in Nature. It cannot be derived from classical bifurcation theory.

Third, it offers a guiding principle for studying phase transitions of natural systems.

Ehrenfest Classification

Paul Ehrenfest was the first who gave a definition and a classification of thermodynamic phase transitions in terms of singularities, at the critical threshold, of such thermodynamic observable parameters as heat capacity, magnetic susceptibility, etc., which are observable; see [START_REF] Ehrenfest | Phasenumwandlungen im ueblichen und erweiterten Sinn, classifiziert nach den entsprechenden Singularitaeten des thermodynamischen Potentiales[END_REF]. Phase transitions are then classified based on the behavior of the thermodynamic potentials, and were labeled by the n-th order derivative of the potential that is discontinuous at the transition.

More precisely, consider a thermodynamic system with potential functional F (u, λ). Definition 1 (Ehrenfest Classification Scheme). Let λ 0 be a critical point, and u λ are equilibrium states near λ = λ 0 . A phase transition of the system at λ = λ 0 is of the n-th order if

∂ k F (u λ , λ)
∂λ k are continuous with respect to λ at λ = λ 0 for any 0 ≤ k ≤ n -1, and

∂ n F (u λ , λ) ∂λ n is discontinuous at λ = λ 0 .
Under this scheme, phase transitions are labeled by the lowest-order derivative of the potential that is discontinuous at the transition. Firstorder phase transitions exhibit a discontinuity in the first derivative of the potential with a thermodynamic variable. In principle, there could be third, fourth, and higher-order phase transitions. In fact, we shall prove that there are only first, second and third order transitions. Now we recall the physical meaning of the derivatives of F (u, λ) on λ in statistical physics. For a thermodynamic system, the first order derivatives include the entropy: S = -∂F /∂T , and the (observable) phase volume: V = ∂F /∂p or density. The second order derivatives are the heat capacity in the constant pressure:

C p = -T ∂ 2 F ∂T 2 , the compression coefficient κ = -1 V ∂ 2 F ∂p 2 , the thermal expansion coefficient α = 1 V ∂ 2 F
∂T ∂p , and magnetic susceptibility: ∂ 2 F /∂M ∂H, where M is the magnetization of the material and H is magnetic field strength. Thus, for the first-order phase transition, the discontinuity of ∂F/∂p at critical point λ 0 = (T 0 , p 0 ) implies the discontinuity of phase volume:

(7) ∆V = V 2 -V 1 = ∂F + ∂p - ∂F - ∂p ,
and the discontinuity of ∂F/∂T implies that there is a gap between both phase entropies:

(8) ∆S = S 2 -S 1 = - ∂F + ∂T + ∂F - ∂T .
Physically, ∆S is a non-measurable quantity, hence, we always use latent heat ∆Q = T ∆S to determine the first order phase transition, and ∆Q stand for the absorbing heat for transition from phase ū to phase u λ . In the two phase coexistence situation, ∆Q and ∆V are related by the Clausius-Clapeyron equation dp dT = ∆Q T ∆V .

For second-order phase transitions, the variance of heat capacity (or specific heat), compression coefficient, the thermal expansion coefficient and magnetic susceptibility at the critical value, are important measurable physical quantities.

Three Basic Theorems of Thermodynamic Phase Transitions

5.1. First Theorem. We have proved a dynamical transition theorem [START_REF]Phase Transition Dynamics[END_REF]Theorem 2.1.3]. It shows that for [START_REF] Stanley | Introduction to Phase Transitions and Critical Phenomena[END_REF], as soon as the linear instability occurs, the system always undergoes a dynamical transition. As a direct application of this general transition theorem, we arrive immediately at the following first theorem of thermodynamic phase transitions.

Theorem 1. Consider a thermodynamic system [START_REF] Stanley | Introduction to Phase Transitions and Critical Phenomena[END_REF].

Let β 1 (λ), β 2 (λ), • • • ∈ R be eigenvalues 1 of the linear operator L λ . If (9) β i (λ)    < 0 if λ < λ 0 , = 0 if λ = λ 0 , > 0 if λ > λ 0 1 ≤ i ≤ m, β j (λ 0 ) < 0 m + 1 ≤ j,
then the system (6) always undergoes a dynamic transition to one of the three types of dynamic transitions: continuous, catastrophic and random, as λ crosses the critical threshold λ 0 .

In the physics literature, (9) is called the principle of exchange of stabilities. It is important to note that this theorem provides the detailed information for the phase diagram in the control parameter space:

(10) β 1 (λ) = 0.
For example, in the Cahn-Hilliard equation describing the binary system [START_REF]Phase separation of binary systems[END_REF][START_REF]Phase Transition Dynamics[END_REF], the control parameters are λ = (T, u 0 , L), where T is the temperature, u 0 is molar density, and L is the length scale of the system. The phase diagram is then obtained by solving (10). For fixed L, the phase diagram is given by Figure 1: For a fixed molar density u 0 , the T L phase diagram was derived and given by Figure 2: 5.2. Second Theorem. We have at our disposal the Ehrenfest classification and our dynamic classification based on Principle 1; the latter is applicable for both equilibrium and non-equilibrium phase transitions.

The second theorem of thermodynamic phase transitions shows that there are only first, second and third order transitions, and provides the relationship between the Ehrenfest and dynamic classifications. Theorem 2. For the phase transition of a thermodynamic system, there exist only first-order, second-order and third-order phase transitions. Moreover the following relations between the Ehrenfest classification and the dynamical classification hold true: second-order ⇐⇒ continuous first-order ←-catastrophic either first or third-order ←-random first-order -→ either catastrophic or random third-order -→ random with asymmetric fluctuations.

Proof. By Theorem 1, we only have to examine the three types of dynamical transitions. Before we proceed, we observe that the first eigenvalue β 1 (λ) satisfies

(11) β 1 (λ) = α(λ)(λ -λ c ), α(λ c ) = 0.
First, it is clear from the dynamical transition theory [START_REF]Phase Transition Dynamics[END_REF] that a catastrophic dynamic transition must be of first-order.

Second, consider the case where the dynamical transition at λ c be continuous. By the dynamical transition theory [START_REF]Phase Transition Dynamics[END_REF], we have shown that the transition states near λ c are given by ( 12)

u λ = A(λ -λ c ) 1/m + o((λ -λ c ) 1/m ) for λ > λ c , 0 for λ < λ c ,
where m ≥ 2 and A = 0. In fact, physically, 2 ≤ m ≤ 3; see among others [START_REF] Stanley | Introduction to Phase Transitions and Critical Phenomena[END_REF]. Also, in view of the derivation of the equation ( 6) from ( 4) and ( 5), since u is the deviation order function, we let (13)

F (u, λ) = F (u + ū, λ) -F (ū, λ),
where F is the potential functional of the thermodynamic system. Then near λ c , F takes the following form:

(14) F (u λ , λ) = 1 2 β 1 (λ)u 2 λ + h.o.t.
Therefore, in view of (12), we have

(15) F (u λ , λ) = A(λ -λ c ) 1+ 2 m + h.o.t for λ > λ c , 0 for λ < λ c .
Hence a continuous dynamical transition leads to a second-order transition in the Ehrenfest sense. Third, we consider the case where the dynamical transition at λ c is random. If the fluctuation is symmetric, then the transition must be of first-order. Otherwise, it can either be a first-order transition or the transition with the following transition states [START_REF]Phase Transition Dynamics[END_REF]:

(16) u λ = A(λ -λ c ) + h.o.t for λ > λ c , 0 for λ < λ c ,
The potential functional is then written as ( 17)

F (u λ , λ) = A(λ -λ c ) 3 + h.o.t for λ > λ c , 0 for λ < λ c .
This gives rise to a third-order transition. Theorem 1 shows that there are only three types of dynamical transitions, which lead to only first, second and third-order transitions. Therefore, there are no higher than third-order transitions in thermodynamics and the relations between the Ehrenfest classification and the dynamical classification given in the theorem hold true as demonstrated above. The proof is complete.

A few remarks are now in order. First, the Ehrenfest transition types are determined by experiments. The dynamical transition theory provides a systematic theory to easily determine the types and transition states. The above theorem establishes a natural bridge between theory and experiments.

Second, the third-order transition can hardly be determined by thermodynamic parameters experimentally. The dynamical transition theory, however, offers a simple approach to identify the third-order transition.

Third, the above theorem leads to transition diagrams, providing an important characteristic of phase transitions. In fact, the dynamical transition theory includes a systematic technical approach in determining the transition types and states, based on the central manifold reduction idea; see [START_REF]Phase Transition Dynamics[END_REF] for details. In particular, for a thermodynamic phase transition, the approach gives rise to the following transition equation, which encodes the all the needed information on transition types and states:

(18) du λ dt = β 1 (λ)u λ + a(λ)u 2 λ + b(λ)u 3 λ + h.o.t.
Then the transition types are determined by the coefficients a(λ) and b(λ) as follows:

(1) If a(λ c ) = 0, then the transition is random;

( 5.3. Third Theorem. For a thermodynamic system, there is no bifurcation point for given pressure control parameter with the temperature higher than the critical threshold. Also, given pressure, at very high temperature, the only stable state is the basic equilibrium state ū. Consequently, based on the topological degree theory [START_REF] Ma | Bifurcation theory and applications[END_REF], both catastrophic and random transitions will give rise to saddle-node bifurcations at a temperature T 0 c higher than the critical threshold T 1 c . We refer readers to Figure 4 for a better understanding of superheated and supercooled states. Also, if the transition is of the first-order, the jump of the thermodynamic potential ∆F at T 1 c ≤ T ≤ T 0 1 corresponds to the latent heat, and the two metastable states are the superheated and supercooled states of the system.

Namely, we obtain the following third theorem for thermodynamic phase transitions. Theorem 3. For a thermodynamic system (6), the following statements hold true:

(1) both catastrophic and random transitions lead to saddle-node bifurcations; and (2) both latent heat, superheated and supercooled states always accompany the saddle-node bifurcations associated with the first order transitions.

This theorem, together with (18), offers precise information on the dynamic phase diagrams. For example,Figure 4, Figure 5

Figure 1 .

 1 Figure 1. T u 0 -phase diagram for a fixed length scale L; see [3, Fig. 4] and [4, Fig. 3.23].

Figure 2 .

 2 Figure 2. TL phase diagram for fixed u 0 ; see [3, Fig. 6] and [4, Fig. 3.25].

)

  If a(λ c ) = 0, then b(λ c ) > 0 =⇒ catastrophic transition b(λ c ) < 0 =⇒ continuous transition.

Figure 3 .

 3 Figure 3. Lu 0 -transition diagram: region II is the first order transition region with latent heat, and region I is the second order transition region.

1 Figure 4 .

 14 Figure 4. Catastrophic transition: the basic state ψ = 0 in T 1 c < T < T 0 c is supercooled, the transition states Σ 2 T are superheated; see [4, Fig. 3.31 on page 202].

  and Figure 6 are such dynamical phase diagrams, derived from the dynamical transition theory.

Figure 5 .

 5 Figure 5. Dynamic phase diagram for continuous transition; see [4, Fig. 3.3 on p.138].

Figure 6 .

 6 Figure 6. Dynamical phase diagram of random transition for a 2 > 0; see [4, Fig. 3.4 on p. 139].

The linear operator L λ = -∂

F (ū;λ) ∂ui,∂uj is always self-adjoint, and all eigenvalues of L λ are real.

The work was supported in part by the US National Science Foundation (NSF), the Office of Naval Research (ONR) and by the Chinese National Science Foundation.