of Geophysical Research : Solid Earth Strain localization driven by thermal decomposition during seismic shear

Field and laboratory observations show that shear deformation is often extremely localized at seismic slip rates, with a typical deforming zone width on the order of a few tens of microns. This extreme localization can be understood in terms of thermally driven weakening mechanisms. A zone of initially high strain rate will experience more shear heating and thus weaken faster, making it more likely to accommodate subsequent deformation. Fault zones often contain thermally unstable minerals such as clays or carbonates, which devolatilize at the high temperatures attained during seismic slip. In this paper, we investigate how these thermal decomposition reactions drive strain localization when coupled to a model for thermal pressurization of in situ groundwater. Building on Rice et al. (2014), we use a linear stability analysis to predict a localized zone thickness that depends on a combination of hydraulic, frictional, and thermochemical properties of the deforming fault rock. Numerical simulations show that the onset of thermal decomposition drives additional strain localization when compared with thermal pressurization alone and predict localized zone thicknesses of ∼7 and ∼13 μm for lizardite and calcite, respectively. Finally we show how thermal diffusion and the endothermic reaction combine to limit the peak temperature of the fault and that the pore fluid released by the reaction provides additional weakening of ∼20–40% of the initial strength.


Introduction
Field studies of fault zones show a hierarchical structure, with a fault core composed of ultracataclasite and fault gouge sitting within a broader damage zone [e.g., Faulkner et al., 2010]. Further investigation reveals a zone of highly localized shear on the order of 10-300 μm wide nested within the fault core [Heermance et al., 2003;Chester et al., 2003;De Paola et al., 2008;Collettini et al., 2013;Bullock et al., 2014]. These field observations are consistent with laboratory observations from high-velocity rotary shear experiments, which reveal micron-scale strain localization at slip rates of order 1 m/s. In experimental deformation tests performed at a slip rate of 1 m/s on a dry, natural clay-bearing fault gouge, Brantut et al. [2008] identified a zone of darker material ∼1-10 μm wide that, due to the lack of other indicators of deformation elsewhere in the sample, was interpreted as the main slipping zone in the experiment. In similar deformation experiments performed under wet conditions on similar natural fault zone materials, Kitajima et al. [2010] showed that a 100 μm thick zone of extremely fine grained material with a strong foliation forms at seismic slip rates. This zone is thought to have accommodated the majority of deformation in the experiment, and the foliation may indicate that the width of a single localized shear zone is much smaller than 100 μm. A more detailed discussion of these observations and further examples of micron-scale strain localization can be found in the introduction to Rice et al. [2014].
In general, strain localization should be expected in gouge undergoing thermally driven dynamic weakening. If a region is straining faster than the surrounding material, then it will experience more shear heating; more shear heating leads to faster weakening; weaker regions of the gouge layer will be more likely to accommodate subsequent deformation. Two distinct thermally driven dynamic weakening mechanisms can be considered in fluid-saturated fault rocks: thermal pressurization and thermal decomposition. Both mechanisms rely on rapid increases in pore fluid pressure leading to an overall strength decrease. Thermal pressurization is due to thermal expansion of the pore fluid and pore volume as the fluid-saturated gouge material is heated. If the heating occurs faster than the pore fluid can drain from the gouge, then the pore pressure will increase, Journal of Geophysical Research: Solid Earth 10.1002 leading to dynamic weakening [Lachenbruch, 1980;Smith, 1985, 1987]. Thermal decomposition corresponds to the chemical breakdown and devolatilization of hydrated or carbonated minerals, such as clays or calcite, which are often present in faults. Such chemical transformations provide an independent source of fluid pressure that is important at high temperatures when the reaction kinetics are fast compared to the timescale for seismic slip. High-velocity friction experiments have revealed several devolatilization reactions that can occur on timescales of a few seconds. Evidence for thermal decomposition was shown for siderite [Han et al., 2007a], calcite [Han et al., 2007b], serpentinites [Hirose and Bystricky, 2007;Proctor et al., 2014], kaolinite [Brantut et al., 2008], dolomite [De Paola et al., 2011], and gypsum [Brantut et al., 2011]. Evidence of thermal decomposition during seismic slip has also been inferred from field observations of faults [Collettini et al., 2013;Bullock et al., 2014]. In the crustal seismogenic zone these decomposition reactions are typically endothermic, and at a fixed pressure the reaction products occupy a larger volume than the reactants for undrained conditions. The combination of these two effects implies that the onset of rapid thermal decomposition leads to an increase in the pore pressure and a plateau in the maximum temperature, as shown theoretically in , Brantut et al. [2010], and experimentally in Brantut et al. [2011]. Throughout this manuscript we will refer to dynamic weakening exclusively due to thermal expansion of in situ pore fluid as thermal pressurization, and dynamic weakening due to the release of additional pore fluid during a devolatilization reaction as thermal decomposition, though what we call thermal decomposition has also been called thermochemical pressurization [Brantut et al., 2010].
The width of the deforming zone during seismic shear, which this paper attempts to constrain, is of crucial importance in theoretical models of thermally driven dynamic weakening. Lachenbruch [1980] showed that for undrained and adiabatic conditions, dynamic weakening by thermal pressurization is controlled by a critical weakening strain, so the slip-weakening distance for thermal pressurization is proportional to the deforming zone thickness. This may explain why the gouge layer thickness plays a role in determining if a rupture propagates as a crack-like rupture or slip pulse in the results of Noda et al. [2009]. Another example can be found in Garagash [2012], which showed that for steadily propagating slip pulses, thinner deforming zones lead to smaller slips and faster rupture velocities.
For thermal pressurization alone, Rice et al. [2014] used a linear stability analysis to predict how the localized zone thickness depends on the gouge properties. This analysis was complemented by the numerical simulations presented in Platt et al. [2014] that went beyond the linear regime. For strain rate localization stabilized by frictional rate strengthening alone the localized zone thickness is set by a balance between thermal pressurization, hydrothermal diffusion, and frictional strengthening. Using hydraulic and thermal parameters from Rempel and Rice [2006], which model a depth of 7 km as a typical centroidal depth for a crustal seismogenic zone, and friction data from Blanpied et al. [1998], they predicted that the localized zone is between 4 and 44 μm wide, with the smaller number assuming parameters based on experiments on undamaged gouge and the larger number representing an estimate of the effect of damage at the onset of rapid shear (e.g., microcracking). Platt et al. [2014] also showed that strain localization has a dramatic effect on the temperature and strength evolution of the gouge. As straining localizes, the frictional heating is focused into a narrower zone, leading to an acceleration in dynamic weakening and a temperature rise much larger than that predicted when strain rate localization is not accounted for. In this paper we extend the work in Rice et al. [2014] and Platt et al. [2014] to account for thermal decomposition. A linear stability analysis leads to a prediction for the localized zone thickness as a function of the gouge properties and current fault temperature, and these predictions are tested using numerical simulations. Next we show how thermal decomposition combines with thermal diffusion to limit the maximum temperature rise, and how we can estimate the temperature at which thermal decomposition operates. Finally we study the strength evolution during localization, showing that the onset of thermal decomposition leads to a sudden strength drop of ∼20-40% of the initial strength.

Model Derivation
In this section we derive a model for a fluid-saturated gouge material sheared between two undeforming thermoporoelastic half-spaces that allow diffusion of heat and pore fluid, the same geometry used in Platt et al. [2014]. In this one-dimensional model the only nonzero velocity component, u(y, t), is parallel to the fault zone and depends only on the coordinate perpendicular to the direction of slip y, and the time since shear commenced t. The gouge layer has a finite thickness h and the half-spaces are moved relative to each other at a kinematically imposed slip rate V, which leads to a nominal strain rate in the gouge layer oḟo = V∕h. A sketch of this geometry is shown in Figure 1. A gouge layer with a finite thickness h is sheared between two undeforming thermoporoelastic half-spaces moving relative to each other at a slip rate V, leading to a nominal strain rate oḟo = V∕h within the gouge layer. In this one-dimensional model we only account for variations in the across-fault direction y. The straining is allowed to localize within the gouge layer, as shown by the Gaussian strain rate profile sketched within the gouge layer. The width W of the zone of localized straining is then estimated as twice the root-mean-square width of the Gaussian.
Our derivation extends the model of Rice et al. [2014] to account for thermal decomposition, which is modeled using the ideas in ], and Brantut et al. [2010. For clarity we model a single reaction, but the modeling framework used is general and results are given for other decomposition reactions.  hypothesized that the short distances associated with hydrothermal diffusion make inertial effects within the gouge layer unimportant. This hypothesis was tested in Platt et al. [2014] and found to be true for typical seismogenic conditions. Based on this we use the equations for mechanical equilibrium to model the stresses within the gouge layer,

Mechanical Equilibrium
where is the shear stress in the gouge material, and n is the normal stress on the gouge layer. As in Rice et al. [2014] and Platt et al. [2014] we assume that the normal stress on the gouge layer is constant throughout shear. The assumed quasi-static behavior forces the shear stress to be constant throughout the layer, and thus, is at most a function of t.

Gouge Friction
The shear stress is linked to the normal stress using a friction coefficient f and the Terzaghi effective stress where p = p(y, t) is the local pore pressure. For a constant or rate-weakening friction coefficient, and neglecting dilatancy, only two forms of deformation satisfy mechanical equilibrium: uniform shear of the gouge layer or slip on the plane of maximum pore pressure . Small perturbations away from uniform shearing will be unstable, and the deformation will collapse to a plane. However, when the friction coefficient is rate strengthening, a finite thickness shear zone can exist.
Current high-velocity friction experiments are unable to separate out the complicated temperature and pore fluid effects to provide a friction law as a function of strain rate alone at seismic strain rates. Lacking such a friction law we assume the steady state friction law It is important to note that equations (2)-(4) link the pore pressure and strain rate profiles within the gouge layer. Locations with high pore pressures will have smaller effective stresses, corresponding to a higher strain rate for the rate-strengthening friction law assumed in this paper. This makes it crucial to understand how spatial variations in pore pressure across the gouge layer develop due to the positive feedback between frictional heating and the two thermally driven weakening mechanisms.
As discussed in Rice et al. [2014], the friction law in equation (4) neglects important effects of temperature, mineralogy, and state evolution and is unlikely to accurately describe the frictional response of gouge at the seismic slip rates considered here. However, it is important to note that the results presented in this paper will be qualitatively the same for any rate-strengthening friction law. For a guide on how to reinterpret our results for other friction laws we refer the reader to Rice et al. [2014], which showed how effective values of f o and (a − b) could be extracted from other friction laws of the form f (̇).

Conservation of Pore Fluid Mass
Defining m to be the mass of pore fluid per unit reference volume of porous material, we can write the conservation of pore fluid mass as, where q f is pore fluid flux and m d is the mass of pore fluid released by the thermal decomposition reaction per unit reference volume. For a saturated gouge m = n f , where f is the pore fluid density and n is the pore volume fraction. It follows that where we have split the porosity change into an elastic component n el and an inelastic component n in .
The new derivatives for f and the elastic porosity n el can be linked to changes in pore pressure and temperature using where T = T(y, t) is the temperature, n and f are the pore volume and pore fluid compressibilities, and n and f are the thermal expansion coefficients for pore volume and pore fluid, respectively. Platt et al. [2014] showed that dilatant effects that depend on strain rate alone are expected to have minimal impact on strain localization at seismic depths, although they may play an important role at the lower effective stresses used in high-velocity friction experiments. Motivated by this we neglect dilatancy and assume that all inelastic porosity change is due to the thermal decomposition reaction.
Denoting the mass of a chemical species x per unit reference volume of fluid-saturated gouge by m x , and the density of that chemical species by x , we can express the rate of inelastic porosity change for the decarbonation reaction in equation (1) using the rate of volume change for each of the solid phases, Next, using the molar masses M x for a chemical species x in equation (1), we can tie the volume changes to the mass of pore fluid released, Finally, we relate the pore fluid flux q f to the pore pressure gradient across the fault using Darcy's law, where k is the intrinsic permeability and f is the pore fluid viscosity.
Combining equations (5)-(8), (12), and (13) and neglecting the dependence of the hydraulic properties on pore pressure, temperature, and porosity, we arrive at Here is the storage coefficient and Λ is the ratio of pore pressure change to temperature change for thermal pressurization under undrained and adiabatic conditions [Lachenbruch, 1980]. We define the hydraulic diffusivity and the inelastic porosity created per unit mass of fluid released All three terms on the right-hand side of equation (14) have a clear physical interpretation. The first represents thermal pressurization of the pore fluid, the second term models hydraulic diffusion, and the final term models the pore pressure generated by thermal decomposition.
Reactant depletion may become important at large slips. To model this, we consider the total pore fluid mass that can be released by a decomposition reaction per unit volume of fluid-saturated gouge, m tot d . Using this, we define the reaction progress as the mass of pore fluid released divided by the total mass of pore fluid that could be released in a fully completed reaction, For this definition = 0 represents virgin material and = 1 indicates full reactant depletion. Using this definition, we can write the final term in equation (14) as Note that the total pore fluid mass m tot d that can be released during decomposition will depend on the specific reaction activated as well as the initial reactant mass fraction of the gouge. To separate these two effects, we write 1 where we have defined, Here m 100% d is the pore fluid mass per reference volume released by a completed reaction in a pure material and thus, P r is the pore pressure generated by a completed reaction of a pure reactant under undrained and isothermal conditions.
The final equation modeling the conservation of pore fluid mass is,

Conservation of Energy
Assuming that energy is generated by frictional heating in the gouge layer and absorbed by the endothermic reaction, we can write the conservation of energy as where c is the effective heat capacity per unit reference volume and ΔH is the enthalpy change associated with the generation of a unit mass of pore fluid through thermal decomposition. We will study endothermic reactions, and thus ΔH > 0. To model the heat flux, we use Fourier's law, where K is the thermal conductivity, which is assumed to be constant. Equations (24) and (25) neglect small additional terms modeling the work done by the normal stress and pore pressure, and heat transfer due to fluid flow. These are common assumptions justified in Smith [1985, 1987] for representative fault gouge permeabilities. Combining equations (24) and (25), we find where the thermal diffusivity is defined as As in the previous subsection we recast the pore fluid mass released per unit reference volume m d in terms of the reaction progress by normalizing the total mass of pore fluid released by the total amount that would be released in a completed reaction. Equation (26) becomes where The parameter E r is the net temperature change for a completed reaction in a pure material under adiabatic and isobaric conditions.

Reaction Kinetics
Finally we model the reaction kinetics, which control how fast thermal decomposition progresses. We assume a first-order reaction with an Arrhenius temperature dependence, where A is the rate constant for the reaction, Q is the activation energy for the reaction, and R is the gas constant. To recast this in terms of the reaction progress , we divide through by m tot d to find, The reaction kinetic has a sensitive dependence on temperature, with higher temperatures leading to a more vigorous reaction. For a fixed temperature a lower value of leads to a larger reaction rate, and when = 1, the reaction is complete and thus the reaction rate is zero.
The strong temperature dependence of the reaction kinetic allows us to predict when each of the dynamic weakening mechanisms will dominate. At low temperatures the reaction rate for thermal decomposition will be slow, and we expect thermal pressurization to dominate. As the temperature rises the reaction rate increases and may reach a temperature where thermal decomposition dominates. We do not expect to exceed this temperature because any increase in temperature will be absorbed by the enthalpy change of the endothermic reaction, as can be seen clearly in the numerical simulations of , Brantut et al. [2010].

Parameter Values
The model presented above is rich in parameters. In this section we will choose typical values for these parameters and discuss how well constrained each parameter is. In Appendix A we nondimensionalize the model from the previous section, showing that there are eight dimensionless parameters, each with a clear physical meaning.
The hydraulic parameters are highly variable and depend on pore pressure, temperature, and the amount of damage the surrounding material has sustained. We use the path-averaged parameters modeling a damaged material from Rempel and Rice [2006], which are based on Tables 1-3 in  and the procedures in  to account for variations in the hydraulic properties due to damage as well as pore pressure and temperature changes. This parameter set models a depth of 7 km, which is a typical centroidal depth for rupture zones of crustal earthquakes. The hydraulic diffusivity is chosen to be 6.71 mm 2 /s, the storage capacity to be = 2.97 × 10 −10 /Pa, and Λ = 0.3 MPa/K. A detailed discussion of the assumptions and laboratory measurements used to develop these parameters can be found in  and Rempel and Rice [2006].
Compared to the hydraulic parameters, the thermal parameters th and c are relatively well constrained. Following our choice of the path-averaged parameter set modeling a damaged material taken from Rempel and , we choose the effective heat capacity per unit reference volume to be c = 2.7 MPa/K, and the thermal diffusivity to be th = 0.54 mm 2 /s. Both of these fall in the typical range of values quoted in .
The frictional parameters are as variable as the hydraulic parameters. The friction law assumed here-given in equation (4)-is motivated by steady state friction values from low strain rate experiments [Dietrich, 1979], and the applicability to the rapid shear considered here is unclear. However, the analysis provided below is qualitatively similar for any rate-strengthening friction law, and Rice et al. [2014] shows how effective values of f o and (a − b) could be inferred from a general friction law f = f (̇). Understanding these limitations, we choose f o = 0.6 and a − b = 0.025, both in the observed range for low strain rate experiments on granite under hydrothermal conditions [Blanpied et al., 1998], though a wide range of other choices for f o and (a − b) could be justified.
The numerical calculations in this paper are performed for calcite decarbonation and lizardite dehydration reactions, and our results are discussed for two other reactions in section 6. We will first discuss the parameters associated with the decarbonation reaction given in equation (1) closely following . Dollimore et al. [1996] reported values of Q = 319 kJ/mol, and A = 2.95 × 10 15 s −1 for the decarbonation of calcite mixed with silica. These kinetic parameters neglect any dependence of reaction rate on the partial pressure of carbon dioxide, but more accurate models could be constructed to account for this. The sign of this effect can be understood using Le Chatelier's principle, and for a fixed temperature and reactant mass, as the partial pressure of carbon dioxide increases, the reaction rate will decrease. For the isobaric mode the enthalpy change of the reaction is equal to the activation energy [L'vov, 2002]. Thus, using the molar mass of carbon dioxide, M CO 2 = 44 g/mol, we find ΔH = 7.25 MJ/kg. The value of can be calculated using the parameter values from , leading to = 0.46 × 10 −3 m 3 /kg. Using the molecular weights and density from  and the path-averaged porosity n = 0.043 from Rempel and Rice [2006], we find m 100% d = 1140 kg/m 3 .
Choosing the fluid density is hard for decarbonation reactions in a water-saturated gouge since the in situ pore fluid is different from the fluid released by the decomposition reaction. We assume that the appropriate PLATT ET AL. LOCALIZATION AND THERMAL DECOMPOSITION    to account for damage to the gouge material at the onset of shearing and parameter changes due to changes in pore pressure and temperature. Frictional parameters are based on Blanpied et al. [1998]. A fuller discussion on the origin of the parameters can be found in Rice et al. [2014]. density is that of supercritical carbon dioxide and calculate this using the equation of state in Saxena and Fei [1987]. To determine the conditions at which to evaluate this equation of state, we must estimate the conditions at which thermal decomposition operates. We assume that thermal decomposition begins at a pore pressure of p = p a + 0.5( n − p a ), where p a is the ambient pore pressure. This is intended to crudely model a gouge that has already experienced significant dynamic weakening due to thermal pressurization before the reaction is triggered. To estimate the temperature T r at which thermal decomposition operates, we assume that all of the frictional heating is absorbed by the endothermic reaction and reactant depletion is negligible. These assumptions are consistent with the results in  and lead to To evaluate T r , we usem = 0.5 and a heating ratė= 378 MPa/ms, which corresponds to the shear stress = f o ( n − p a )∕2 and the strain rate implied by a slip rate of 1 m/s accommodated across a zone 100 μm wide. These choices lead to T r = 960 ∘ C, f = 418 kg/m 3 , E r = 3.06 × 10 3 ∘ C, and P r = 7.42 GPa for calcite decarbonation. Note that the value of E r is used to predict T r , which is then used to determine our value of P r .
Next we discuss the dehydration of lizardite into talc, olivine, and water: Llana-Fúnez et al. [2007] provide a range of kinetic parameters associated with the dehydration of intact blocks or powders of lizardite. Here we use a rate constant A = 6.40 × 10 17 s −1 and an activation energy Q = 328 kJ/mol, which correspond to the dehydration kinetics of a mixture of lizardite and brucite (originally reported in Wegner and Ernst [1983]). The reaction enthalpy is calculated using the thermodynamic software Geotab from Berman [1991], which yields ΔH = 2.56 MJ/kg. From the stoichiometry of the reaction and the densities of the reactants and products we calculate the solid volume change = 0.88 × 10 −3 m 3 /kg and the total mass of water released by the reaction m 100% d = 240 kg/m 3 . Finally, we use a procedure similar to that outlined above to determine the density of water of 267 kg/m 3 at the reaction temperature. For the dehydration of lizardite we find E r = 275 ∘ C and P r = 2.80 GPa.
Aside from the decarbonation of calcite and the dehydration of lizardite, a wide variety of other thermal decomposition reactions can be triggered during earthquake slip. Potential candidates include carbonates such as dolomite, magnesite, and siderite, as well as hydrous minerals such as gypsum and phyllosilicates (e.g., clays, serpentines, and talc). Our model requires a number of reaction parameters that are rarely available in a consistent set in the published literature. The full set of reaction parameters could be obtained for the dehydration reactions of illite-muscovite mixtures and talc. The dehydration of illite-muscovite was studied experimentally by Hirono and Tanikawa [2011], who provide all the relevant parameters needed for our model. In the case of talc dehydration, we used the kinetics reported by Bose and Ganguly [1994] and determined the enthalpy change using Geotab [Berman, 1991].
The hydraulic, frictional, and thermal parameter values are summarized in Table 1, and the parameters for the four thermal decomposition reactions are summarized in Table 2.   Bose and Ganguly [1994], reaction enthalpy from Geotab [Berman, 1991]. e Note that the values reported are per unit fluid mass released.

Linear Stability Analysis
In this section we predict the localized zone thickness using a linear stability analysis. To make progress analytically, we linearize the reaction kinetic about = 0 and a current fault temperature T = T f , leading to Given that the Arrhenius factor has a strong dependence on temperature, such a linearization will have a very limited range of validity. However, performing the linear stability analysis with the linearized reaction, kinetic above is equivalent to performing the linear stability analysis with the Arrhenius reaction kinetic and then freezing the coefficients in the resulting time-dependent linear system. This means that the linearized reaction kinetic is valid provided that perturbations in temperature are small, which is expected to be true at the onset localization. Thus, despite the rather crude approximation made when linearizing a highly nonlinear function, we will find that the linearized analysis does convey some key qualitative features observed in the more precise nonlinear solutions presented later in this paper.
Inserting the linearized reaction kinetic into equations (23) and (28), we arrive at As in Rice et al. [2014] we now perturb about the solution for uniform shearing, where the uniform shear solution is denoted by a subscript 0. This is done by setting, 10.1002/2014JB011493 wherē0(t) = n − p 0 (t) is the effective stress for uniform shear and we have assumed thaṫ0 is equal to the nominal strain ratėo. Somewhat surprisingly we do not need to solve for the uniform solution since it does not enter into the final linearized system for perturbations in p and T.
Substituting (38) into the model and linearizing, we find that Next we assume that the perturbation is proportional to a Fourier mode with a wavelength , This simplifies equations (39) tō0 Eliminating the only time-dependent term in the system,̄0(t), we arrive at a linear system with constant coefficients, Equations (42) can be solved by assuming pore pressure and temperature perturbations of the form A nontrivial solution to the linear system exists only when  Tables 1 and 2, a reactant mass fractionm = 0.5, and a nominal strain ratėo = 10, 000 s −1 . The horizontal dotted lines show LT pT and HT pT for both materials. The vertical lines show the location of the temperature T r predicted by equation (32) assuminġ= 378 MPa/ms. As expected we see that at low temperatures, the critical half wavelength is equal to LT pT and for high temperatures the critical half wavelength is equal to HT pT , with a smooth transition between the two regimes occurring at intermediate temperatures. Our prediction for the temperature at which thermal decomposition operates at lies in this intermediate temperature regime, so it is unlikely that the high-temperature limit of the linear stability analysis will provide a good quantitative prediction for the localized zone thickness.
Equation (44) determines the growth rate s of a perturbation with a given wavelength , allowing us to determine the stability of the uniform shear. Whenever the real part of s is positive, the perturbations will grow unstably, and whenever the real part of s is negative, the perturbation will decay. The critical wavelength that separates growing and decaying perturbations in p and T, which we call pT following the notation in Rice et al. [2014], occurs when the real part of s is zero. This critical wavelength will be used to predict a localized zone thickness.
We can identify two physically instructive limits from equation (44), one for low temperatures where thermal decomposition is negligible, and the other for high temperatures where the thermal decomposition dominates thermal pressurization. To study the low-temperature (LT) limit, we set 1 = 0, corresponding to a reaction rate so slow that thermal decomposition can be neglected. We recover the system of equations analyzed in Rice et al. [2014], and the critical wavelength for perturbations in p and T is given by This critical wavelength is set by a balance between frictional rate strengthening, thermal pressurization, and hydrothermal diffusion.
Next we study the high-temperature (HT) limit, where thermal decomposition dominates thermal pressurization. Numerical solutions of (44) show that when the real part of s is zero, the imaginary component of s is also zero. This allows us to find a closed form solution for pT by setting s = 0 and neglecting the thermal diffusion term, which is equivalent to assuming that at high temperatures, the endothermic reaction eliminates temperature gradients much faster than thermal diffusion. Equation (44) then becomes which can be solved to find Interestingly the critical wavelength is independent of any reaction kinetic parameters (i.e., A and Q) and the reactant mass fraction. The reaction controls the localized zone width through the parameters E r and P r . We see that the endothermic nature of the reaction acts to widen the localized zone, while the pore pressure generated by the reaction acts to thin the localized zone.
Next we test the above predictions by finding the critical wavelength pT numerically for a wide range of values of T a . Figure 2 shows how the critical wavelength varies for calcite and lizardite using the parameters in Tables 1  and 2, a reactant mass fractionm = 0.5, and a strain ratėo = 10, 000 s −1 , which is equivalent to a slip rate of 1 m/s accommodated across a zone 100 μm wide. For comparison we show the low-and high-temperature

Journal of Geophysical Research: Solid Earth
10.1002/2014JB011493 limits LT pT and HT pT for both materials using horizontal dotted lines. We see that the numerically calculated critical wavelength agrees with the appropriate limit for extreme values of T f , and in the intermediate region we see a smooth transition between one critical wavelength and the other.
Finally, to determine where we expect typical temperatures during thermal decomposition to lie with respect to the high-and low-temperature limits, we plot the reaction temperature T r estimated in equation (32) for both materials using vertical dashed lines. We see that T r lies in the intermediate temperature regime, and thus, the simple formula in equation (47) may not be a good prediction for the localized zone thickness when thermal decomposition is active.

Predicting a Localized Zone Thickness
It is important to note that the critical wavelengths LT pT and HT pT depend on the strain ratėo. Following the procedure in Rice et al. [2014] we now eliminatėo from the two critical wavelengths to find the linear stability analysis (LSA) prediction for the localized zone thickness W LSA as a function of the gouge properties and the slip rate V. We set For the high-temperature limit this leads to the formula and in the low-temperature limit we find As shown in Rice et al. [2014], the linear stability analysis presented in this section can be specialized for a gouge layer of thickness h sheared between rigid, impermeable, and thermally insulating blocks moving relative to each other with a slip rate V. In this case the width W LSA corresponds to the widest possible gouge layer that can be sheared uniformly. These boundary conditions are different from the geometry used in the numerical simulations, but we will show that the linear stability analysis is still able to predict important features seen in the numerical simulations. It should also be noted that to predict the localized zone thickness we have used the critical half wavelength separating growing and decaying perturbations in pore pressure and temperature, not the critical half wavelength that controls perturbations in strain rate. However, Rice et al. [2014] showed that for (a − b) ≪ f o , the two wavelengths are almost equivalent, so the use of pT to predict the localized zone thickness is justified.
As shown in Figure 2, the reaction temperature T r predicted in equation (32) does not fall in the high-temperature regime. Motivated by this we now develop a more complicated prediction for the localized zone thickness in the intermediate temperature range between the high-temperature and low-temperature limits. As before we set s = 0 in equation (44), leading to a quadratic equation for 2 Using the definitions in equation (48), we turn this quadratic into an equation for the localized zone thickness in the intermediate regime W int , As expected, in the high-temperature limit (i.e., 1 2 → ∞) the final term in equation (52)

Journal of Geophysical Research: Solid Earth
This formula is more cumbersome than that given in equation (49) but in the next section we will show that it provides predictions that agree more closely with the results of numerical simulations. However, the more accurate prediction comes at a price and we now must know the kinetic parameters A and Q as well as an estimate of the current fault temperature T f . Equations (49) and (53) are the key results of this study and provide a framework to understand the different physical balances that control the localized zone thickness when thermal decomposition is active.

Shear of a Finite Width Layer
In this section we solve numerically for a gouge layer with a finite width h sheared between two undeforming thermoporoelastic half-spaces that conduct heat and pore fluid moving relative to each other with a slip rate V, the same geometry assumed in Platt et al. [2014]. A sketch of this geometry is shown in Figure 1. At each time step the pore pressure and temperature are updated using equations (23), (28), and (31). To update the shear stress, we require one additional condition. As in Platt et al. [2014] we use which forces the total straining within the gouge layer to equal the total slip rate V accommodated across the gouge layer.
The initial conditions are set to the ambient conditions p = p a and T = T a , and a uniform strain ratė=̇o throughout the gouge layer. To be consistent with the parameters in Rempel and , which are intended to model a depth of 7 km, we choose p a = 70 MPa and T a = 210 ∘ C. This is equivalent to an assumed geotherm of 30 ∘ C/km and a hydrostatic pore pressure gradient of 18 MPa/km.
Note that the geometry used in the numerical simulations is different from the impermeable and thermally insulating boundary conditions assumed in the linear stability analysis. However, as shown in Platt et al. [2014], this is not expected to matter when deformation localizes to a zone much narrower than the gouge layer thickness because the physical balances that control strain rate localization in our simulations will be exactly the same as in the linear stability analysis. Furthermore, hydrothermal diffusion from the gouge layer into the adjacent half-spaces introduces small variations away from the initially uniform pore pressure and temperature profiles, with the largest pore pressures and temperatures near the center of the gouge layer. Strain rate localization naturally develops from this initial perturbation, which has a wavelength comparable to the gouge layer thickness, and thus, we do not need to seed our calculations with a small initial perturbation away from uniform straining.
During the initial stages of deformation the reaction rate is slow, making thermal decomposition negligible. For certain gouge properties the maximum temperature within the gouge layer may eventually become large enough to trigger thermal decomposition. Throughout this section we will focus on this transition from thermal pressurization to thermal decomposition and the behavior of the system after thermal decomposition is triggered. The behavior before thermal decomposition is triggered, where dynamic weakening occurs due to thermal pressurization alone, was analyzed in Platt et al. [2014].
A simple test to determine if thermal decomposition will be triggered in our simulations is to compare the maximum temperature rise for a gouge layer undergoing thermal pressurization alone with the temperature predicted by equation (32). If the two temperatures are comparable or the prediction from equation (56) is larger than the value from equation (32), then it is likely that thermal decomposition will be triggered. All simulations reported here were designed to trigger thermal decomposition, though , indicating that thermal decomposition can be neglected during the initial stages of deformation. Eventually thermal decomposition becomes important anḋm ax increases to a new peak valuėT D peak . Following the peak̇m ax decays, but the values are always above those for thermal pressurization alone. The minimum and maximum strain rates used to calculate Δṫare shown by the black plus and black cross.
we performed other simulations with a larger value of Λ and found that thermal decomposition was rarely triggered.
We will begin by discussing how thermal decomposition drives strain localization during seismic shear, move on to show how thermal diffusion and the endothermic reaction limit the peak temperature, and end by illustrating how the onset of thermal decomposition leads to a sudden strength drop.

Localized Zone Thickness
In this subsection we will study how the localized zone thickness evolves when thermal decomposition is triggered. Following Platt et al. [2014] we define the maximum strain rate within the gouge layer to bė Because the total straining in the layer is fixed by the slip rate V (see equation (55)),̇m ax can be used as a proxy for the localized zone thickness, with a larger value oḟm ax indicating a thinner localized zone. Figure 3 shows hoẇm ax evolves for the thermal decomposition of calcite and lizardite. This plot was generated using the parameters in Tables 1 and 2, a gouge layer thickness h = 1 mm, a slip rate V = 1 m/s, and a reactant mass fractionm = 0.5. For comparison the solution from Platt et al. [2014] that neglects thermal decomposition and models thermal pressurization alone (i.e., E r = P r = 0) is shown by the black dashed curve. As expected our results initially match the calculation for thermal pressurization alone, corresponding to the initial stages of deformation when the reaction progresses so slowly it can be neglected. When thermal decomposition is triggered, we see thaṫm ax rises to a new peak before decaying. We find that throughout the simulation the shape of the strain rate profile is well described by a Gaussian function, in agreement with the results of Platt et al. [2014] for thermal pressurization alone.
We use the Gaussian shape oḟand the peak strain rate after thermal decomposition is triggereḋT D peak to estimate the localized zone thickness W in the numerical simulations, assuming that W is equal to twice the root-mean-square width of the Gaussian. Integrating condition (55) assuming the Gaussian-shaped strain rate profilėg and that the localized zone thickness is much less than h, we find that If the localized zone thickness is comparable to the gouge layer thickness, then equation (59) is not valid, though a more complicated formula can be found that depends on h, V, anḋT D peak . Figure 4 shows a plot of the strain rate profile at peak localization for the simulation modeling the decarbonation of calcite shown in Figure 3 alongside the Gaussian function given in equation (58). The solid black line indicates where the localized zone thickness is measured when we assume that W is equal to twice the root-mean square width of the Gaussian. We see that twice the root-mean-square width may not be the  Table 1 and the calcite parameters in Table 2, a reactant mass fractionm = 0.5, a slip rate V = 1 m/s, and a gouge layer thickness h = 1 mm. Straining localizes to a zone a few tens of microns wide, and we see great agreement between the numerical simulation and the Gaussian fit. The horizontal lines show the two ways to infer a width from the Gaussian function. The upper black line shows where the width is measured assuming that W is equal to twice the root mean square of the Gaussian, and the lower black line shows where the localized zone thickness is measured when we assume the localized zone thickness is equal to 2W. best measure of the localized zone thickness, and if we integrate equation (58), we find that only ∼ 68% of the deformation occurs between y = −W∕2 and y = +W∕2. A better estimate of the deforming zone thickness may be 2W, and this region of the Gaussian accommodates ∼ 95% of the total straining.
Next we investigate how the localized zone thickness depends on the gouge layer thickness and ambient fault temperature. Figure 5 shows W as a function of the gouge layer thickness h for the parameters in Tables 1 and 2, a reactant mass fractionm = 0.5, and a slip rate V = 1 m/s. We see that the localized zone thickness does not change much as h changes from 100 μm to 1750 μm, replicating the behavior observed in Platt et al. [2014] for thermal pressurization alone. This weak dependence of W on the gouge layer thickness suggests that the localized zone thickness is controlled by the gouge properties and not the initial width of the deforming zone. The small increase in W observed for the smallest values of h is thought to be due to the localized zone thickness becoming comparable to the gouge layer thickness. Figure 5 also shows the dependence of W on the ambient temperature T a . We observe that the localized zone thickness does not vary dramatically as the ambient temperature varies from 150 ∘ C to 420 ∘ C, which is to be expected because this range of ambient temperatures is much lower than the temperature at which thermal decomposition operates.
Having shown that the localized zone thickness when thermal decomposition is active depends weakly on the initial conditions, we now study how W varies with the material properties of the gouge. This parameter Figure 5. A plot showing how the localized zone thickness W depends on the gouge layer thickness h and ambient fault temperature T a for calcite and lizardite. These simulations were performed using the parameters in Tables 1 and 2, a reactant mass fractionm = 0.5, a slip rate V = 1 m/s. The simulations varying T a use a gouge layer thickness h = 0.5 mm. We see that the localized zone thickness is almost independent of the gouge layer thickness. From this we can conclude that the localized zone thickness is controlled by the gouge properties and not the initial thickness of the deforming zone, in agreement with the conclusions from Platt et al. [2014] for strain localization driven by thermal pressurization alone. Furthermore, we see that W is almost independent of T a , which is to be expected since the temperature at which thermal decomposition operates does not depend on the ambient fault temperature.  Tables 1 and 2, a reactant mass fraction m = 0.5, a slip rate V = 1 m/s, and a gouge layer thickness h = 0.5 mm. For comparison we also show the linear stability prediction from equation (49) with the dotted curves, the prediction from equation (60) evaluated using the peak temperature from the numerical simulations with the dashed curves, and the prediction from equation (60) evaluated using the temperature from equation (32) assuminġ= 252 MPa/ms with the dash-dotted curves. The predictions from equation (60) give the best agreement with the numerical simulations, especially when the peak temperature from the numerical simulations is used to evaluate (60).

10.1002/2014JB011493
sweep, shown by the solid curves in Figure 6, covers all the dimensionless parameters in the model except for T I (see Appendix A), which was studied in Figure 5. In each plot one parameter is varied while the remaining parameters are fixed to the values in Tables 1 and 2, a reactant mass fractionm = 0.5, a slip rate V = 1 m/s, and a gouge layer thickness h = 0.5 mm.
We compare the localized zone thicknesses observed in numerical simulations with the linear stability predictions from section 4. First, we use the high-temperature limit from the linear stability analysis, given in equation (49) and shown by the finely dashed curves in Figure 6. We see that the predictions from the high-temperature limit of the linear stability analysis are in qualitative agreement with the localized zone thickness predicted by the numerical simulations, with curves representing the analytic prediction and numerical simulations having roughly similar shapes. However, the quantitative agreement between the two is often quite poor, with equation (49) consistently predicting localized zone thicknesses that are a factor of ∼2-3 smaller than those observed in the numerical simulations. This can be understood by looking at Figure 2, which shows that the endothermic reaction caps the maximum temperature at a value that is less than the lower bound of the high-temperature regime, and thus, W HT is not a good approximation.
Next we fit our simulations using the formula and 1 and 2 are given in equation (35). This is based on the linear stability prediction for the intermediate temperature regime (given in equation (53)) with the pore pressure generated P r replaced by (P r − ΛE r ). This change is made because setting s = 0 in the linear stability analysis removes the effects of thermal pressurization, but inserting equation (28) into equation (23), we see that when the thermal pressurization is accounted for the total pore pressure rise in a completed reaction is P r −ΛE r . For all parameters used in this paper P r > ΛE r and the reaction acts as a pore pressure source.
To evaluate the formula in equation (60), we must assume a current fault temperature T f . In Figure 6 this is done in two ways. First, we use the peak temperature from the numerical simulations, shown by the coarsely dashed curves. In addition, we use the prediction T r from equation (32), shown by the lines with alternating short and long dashes, assuminġ= 252 MPa/ms. This power density is equivalent to an effective stress equal to half of the ambient effective stress, a friction coefficient of 0.6, and a slip rate of 1 m/s accommodated across a deforming zone 150 μm wide. This value oḟhighlights the extreme frictional heating rates produced during seismic slip that make thermal pressurization and thermal decomposition such effective dynamic weakening mechanisms.
We see that the more general formula given in equation (60) provides a much better quantitative fit to the numerical simulations than the simple high-temperature asymptote W HT . Using a single-fitting parameter (the numerical factor of 0.55 in equation (60)), we get good agreement with a parameter sweep over seven dimensionless parameters for both calcite and lizardite. The best fit is obtained when we set T f to be the peak temperature from the simulations, though using the temperature predicted by equation (32) often still gives reasonable agreement.
As shown in Figure 3,̇T D peak is not achieved instantly when thermal decomposition is triggered. Instead,̇m ax increases smoothly from the value predicted for thermal pressurization alone to the new peak value over a finite time. To quantify the time taken for localization to occur after decomposition is triggered, we define Δtṫ o be the time between the local minimum iṅm ax and the second maximuṁT D peak . These points are shown by a black plus and a black cross in Figure 3. Studying how Δṫvaries in the parameter sweeps that led to Figure 6, we find that Δṫincreases as the localized zone thickness decreases. This means that more intense localization develops faster than less intense localization.
Finally we study the decay from the peak strain rate shown in Figure 3. The simulations leading to Figure 6 show that larger values oḟT D peak , and thus smaller values of W, correspond to more rapid decay after the peak strain rate, where we have used the peak value of −̈to measure the speed of decay. This can be seen in Figure 7. A plot showing the evolution of the maximum temperature T max for calcite and lizardite. These simulations were performed using the parameters in Tables 1 and 2, a reactant mass fractionm = 0.5, a slip rate V = 1 m/s, and a gouge layer thickness h = 1 mm. For comparison the solution from Platt et al. [2014] for thermal pressurization alone (i.e., E r = P r = 0) is shown by the dashed black line. Initially our simulations agree with the simulations from Platt et al. [2014], indicating that thermal decomposition can be neglected during the initial stages of deformation. Eventually thermal decomposition becomes important, and T max rises to a new peak before settling onto a slowly decaying plateau. As in Brantut et al. [2010], thermal decomposition leads to a capping of the maximum temperature rise below a typical melting temperature. Figure 3, which shows thaṫm ax decays more rapidly for lizardite than calcite. Decay from the peak strain rate indicates that the localized zone thickens with increasing shear. Thickening of the localized zone makes it hard to describe the localized zone throughout a seismic event using a single width and also means that materials that have different localized zone thicknesses immediately after decomposition is triggered could have very similar thicknesses during the later stages of shear. This can be seen near the end of the simulations in Figure 3 where calcite and lizardite have similar values oḟm ax .

Limiting of Peak Temperature
Next we look at the temperature evolution in the gouge layer. To do this, we define the maximum temperature to be ] .
(62) Figure 7 shows the evolution of T max for the same parameters used to generate the results shown in Figure 3. For comparison we also include the solution from Platt et al. [2014] for thermal pressurization alone, which is shown by the dashed black line in Figure 7. We see that the onset of thermal decomposition initially causes the maximum temperature rise to increase faster than for thermal pressurization alone, a surprising result for an endothermic reaction. This is due to the additional strain rate localization that accompanies the onset of the reaction, focusing frictional heating into a narrower zone. However, the reaction kinetic and thermal diffusion quickly catch up, leading to a peak in T max followed by a gradual decay. This limiting of the temperature is qualitatively similar to the results in Brantut et al. [2010] for a uniformly sheared layer with a thickness between 1 mm and 10 mm, though our peak temperature is higher because straining is more localized in our model, and thus, frictional heating is more intense.
To quantitatively study the maximum temperature rise when thermal decomposition is triggered, we define the peak temperature as Using the parameter sweeps from Figure 6, we plot the dependence of T peak on a range of parameters, as shown in Figure 8. Alongside the numerical simulations we plot the predictions from equation (32) evaluated witḣ= 252 MPa/ms. We see an overall good agreement between the numerical simulations and equation (32). The maximum difference between the two temperatures is typically around 50 ∘ C, though larger discrepancies are seen for the smallest values of E r and A.
To understand the differences between the numerical results and equation (32) we study the magnitude of the three terms on the right-hand side of equation (28) c , The first term models frictional heating, the second term models thermal diffusion, and the final term models the endothermic reaction. At the peak temperature the time derivative of T is zero so these three terms must  Tables 1 and 2, a reactant mass fraction m = 0.5, a slip rate V = 1 m/s, and a gouge layer thickness h = 0.5 mm. For comparison we include the temperature predictions from equation (32) assuminġ= 252 MPa/ms. We see good agreement between our numerical simulations and the simple formula to estimate the temperature at which thermal decomposition operates, with typical discrepancies of ∼ 50 ∘ C.
sum to zero. Physically this means that at the peak temperature the frictional heating is exactly balanced by thermal diffusion and the endothermic reaction. Figure 9 shows how these three terms vary with E r and th for the simulations modeling the decarbonation of calcite shown in Figure 8, alongside the heating rate corresponding tȯ= 252 MPa/ms that was inserted into equation (32) to fit the simulations shown in Figure 8. We observe that for all the simulations shown here, thermal diffusion is more important than the endothermic reaction. Other parameter sweeps show that in almost all simulations thermal diffusion is a factor of 2-3 larger than the reaction, and thus, we conclude Figure 9. A plot showing how the magnitude of frictional heating, thermal diffusion, and the endothermic reaction at peak temperature vary with E r and th for calcite. These plots were generated using the parameters in Tables 1 and 2, a reactant mass fractionm = 0.5, a slip rate V = 1 m/s, and a gouge layer thickness h = 0.5 mm. The black dashed line shows the heating rate corresponding to the value oḟ= 252 MPa/ms used to fit the numerical simulations in Figure 8, where we assumed that frictional heating exactly balances the endothermic reaction. However, this figure shows that thermal diffusion plays a larger role than the reaction in limiting the maximum temperature. In both parameter sweeps the magnitude of the frictional heating and thermal diffusion terms increases as the localized zone thins, and the units in this plot reinforce the extreme heating rates associated with micron-scale strain rate localization.
that thermal diffusion is more important than thermal decomposition in limiting the maximum temperature. This large contribution from thermal diffusion explains why the value oḟthat agrees with the numerical simulations is considerably smaller than the values oḟobserved in the simulations. Micron-scale localization makes thermal diffusion efficient, and the endothermic reaction only needs to offset a percentage of the frictional heating. However, we emphasize that it may not be appropriate to extrapolate this conclusion to other parameter values where the localized zone thickness is much wider than the few tens of microns we observe because the efficiency of thermal diffusion drops rapidly as the localized zone thickness increases, and the endothermic reaction may need to offset all of the frictional heating.
Studying the dependence of the three terms shown in Figure 9 on other parameters allows us to find two general trends that may explain the deviations between (32) and the numerical results. First, for all parameters we see that the magnitude of the frictional heating and thermal diffusion terms increase as the localized zone thickness decreases. These increases largely offset and we see a modest positive correlation between peak temperature and localized zone thickness, indicating that thermal diffusion is decreasing slightly faster than frictional heating as W increases. This can be seen in the subplots of Figure 8 showing the dependence on (a − b), hy , and P r . Second, any parameter change that causes thermal decomposition to be triggered earlier during shear tends to increase the peak temperature above that predicted by equation (32). This trend can be understood by noting that if thermal decomposition is activated earlier, then thermal pressurization will contribute less dynamic weakening, and thus, frictional heating will be more vigorous when the peak temperature is achieved, which equation (32) suggests should lead to a larger peak temperature. This trend can be observed in the subplots of Figure 8 showing the dependence on A and Q, where we see that equation (32) underpredicts the numerically observed value at high A and overpredicts at low A.
Following the peak temperature we see a gradual decrease in the maximum temperature, coinciding with the thickening of the localized zone described in the previous subsection. During this gradual cooling the magnitude of all three terms in equation (28) fall. This is to be expected since frictional heating and thermal diffusion are largely controlled by the width of the deforming zone, and the reaction rate is controlled by the maximum temperature. The ratio of the reaction term to thermal diffusion and the ratio of the reaction term to frictional heating both decay with increasing slip, so as expected thermal decomposition becomes less important as the maximum temperature decays. In a few simulations we observed a gradually increasing temperature after thermal decomposition was triggered, instead of the gradually decreasing temperature seen in Figure 7, with this being particularly common for lizardite.

Strength Drop Due To Thermal Decomposition
In this subsection we study how the onset of thermal decomposition alters the shear strength evolution of the gouge layer. Figure 10 shows the shear strength evolution for calcite and lizardite for the same parameters [2014] that considers dynamic weakening from thermal pressurization alone (i.e., E r = P r = 0) is shown by the dashed black line. Initially our simulations agree with the simulations from Platt et al. [2014], indicating that thermal decomposition can be neglected during the initial stages of deformation. Eventually thermal decomposition becomes important, and the shear strength drops below that predicted for thermal pressurization alone.The locations of the stresses used to calculate the strength drop associated with thermal decomposition are indicated by the black plus symbols.
as those used in Figures 3 and 7. We see that the onset of thermal decomposition leads to a rapid acceleration in dynamic weakening, followed by a return to more gradual weakening.
Platt et al. [2014] showed that for thermal pressurization alone the strength evolution after localization is in good agreement with the Mase-Smith-Rice slip on a plane solution Smith, 1985, 1987;]. The shear strength evolution after thermal decomposition is triggered obviously does not agree with the slip on a plane solution, but the weakening rate −d ∕dt is found to be in reasonable agreement with the slip on a plane solution. Figure 11 shows the weakening rate for calcite and lizardite alongside the weakening rate for the slip on a plane solution. We clearly see a large increase in the weakening rate at the onset of thermal decomposition, but at later times the weakening rate is comparable to that predicted by the slip on a plane solution. Because the weakening rate returns to a value comparable to the value for the slip on a plane solution, weakening due to thermal decomposition can be crudely described as a discrete strength drop coinciding with the onset of the reaction.
Next we quantify how this strength drop depends on the gouge properties. To do this, we first define the strength before thermal decomposition to be the stress at the local minima in the weakening rate associated with the onset of decomposition. Next we define the time at which thermal decomposition stops being important as the moment at which the separation between the weakening rate and the slip on a plane solution is the same as it was before thermal decomposition was triggered. The strength after thermal decomposition is defined as the strength at the time when thermal decomposition stops being important. These two values are used to define the strength drop associated with thermal decomposition Δ , and this strength drop is equivalent to integrating across the large peak in the weakening rate associated with thermal decomposition seen in Figure 11. For clarity we use plus signs to indicate the strength before and after thermal decomposition in the lizardite simulation shown in Figure 10 and use dashed lines to show Δ . Figure 12 shows how the strength drop associated with thermal decomposition varies with the parameters in the model. We see that typical strength drops are between 0.2 and 0.4 of the initial strength 0 , meaning that in these simulations thermal decomposition is as important as thermal pressurization in controlling the total coseismic strength drop of the gouge layer. For the parameter sweeps over E r , P r we see that Δ increases as the localized zone thickness after thermal decomposition is triggered decreases, which is unsurprising since a more vigorous reaction drives more severe localization. It is hard to extend this conclusion to the parameter sweeps over th , hy , and (a − b) because these parameters also influence the evolution of the system before the reaction is triggered, or the parameter sweeps over A and Q since these parameters control the temperature at which the reaction is triggered. This may indicate that Δ is not the perfect variable to measure impact of thermal decomposition or alternatively that the balance between thermal pressurization and thermal decomposition is largely controlled by the amount of slip that occurs before the reaction is triggered and not the properties of the reaction itself. For each individual parameter sweep we observe that larger strength drops occur over shorter times. Finally, we highlight the significant drops in Δ observed when for low values ofm, which we believe are caused by reactant depletion becoming important at low initial reactant mass  Smith, 1985, 1987;. During the initial stages of deformation the two solutions agree, and we see a first spike in weakening rate associated with the onset of localization driven by thermal pressurization. Eventually thermal decomposition is triggered, and we see a second spike in weakening rate, before the two numerical solutions return to a weakening rate comparable to the slip on a plane solution at large slips. The second spike is much larger for lizardite, corresponding to the larger strength drop. This plot shows how weakening due to thermal decomposition can be related to previous solutions for pore fluid weakening and emphasizes the extreme weakening rates associated with the onset of thermal decomposition.
fractions. This conclusion is supported by the fact that the drop in Δ at lowm is more pronounced for lizardite, which has a lower value of E r and thus will be more prone to depletion.

Predictions for Other Common Fault Materials
In this section we use the results from the previous section to make predictions for the peak temperature and localized zone thickness for the four materials listed in Table 2.
First, we predict the maximum temperature during an earthquake-or other rapid slip events such as landslides where thermal decomposition might be triggered [Mitchell et al., 2015]-using equation (32). We use the parameters from Tables 1 and 2, and a heating rate oḟ= 252 MPa/ms. This leads to the predictions shown in Table 2, and we find that the dehydration of talc and the illite/muscovite mixture limits the peak temperature at much higher values than those predicted for the decarbonation of calcite and the dehydration of lizardite. Note that these predictions are the temperatures at which the endothermic reaction proceeds fast enough to offset all of the frictional heating, and it is possible that thermal decomposition may begin to alter the shear strength evolution before the temperature reaches T r and that other physical mechanisms may limit the peak temperature rise to a value lower than our predictions for T r .
Next we predict the localized zone thickness using the high-temperature limit given in equation (49). These predictions are shown in Table 2 for the parameters in Tables 1 and 2, a reactant mass fractionm = 0.5, and a slip rate V = 1 m/s. The localized zone thicknesses predicted for the other dehydration reactions are similar to the predictions for lizardite, with values of about a micron.
Finally we predict the localized zone thickness for the four thermal decomposition reactions using the formula given in equation (60), which is motivated by the linear stability analysis in the intermediate regime and gives the best fit to the numerical simulations. To evaluate this formula, we use equation (32) to estimate the current temperature of the deforming gouge. Using the parameters in Tables 1 and 2, and assuming a reactant mass fractionm = 0.5 and a slip rate V = 1 m/s leads to the predictions given in Table 2. We observe that these predictions are larger than the predictions from the high-temperature limit W HT , as was observed in the numerical simulations shown in section 5. For all four thermal decomposition reactions we predict that the localized zone thickness is approximately 10 μm wide.

Localized Zone Thickness During Seismic Shear
In this manuscript we showed how the localized zone thickness is expected to change during seismic shear. Thermal decomposition can be neglected during the initial stages of deformation and localization is driven by thermal pressurization alone. In this limit the localized zone thickness is set by a balance between thermal pressurization, hydrothermal diffusion, and frictional rate strengthening, and for a fixed slip rate the localized zone thickness can be predicted using the analysis in Rice et al. [2014] and Platt et al. [2014]. At high  Tables 1 and 2, a reactant mass fractionm = 0.5, a slip rate V = 1 m/s, and a gouge layer thickness h = 0.5 mm. We see that a typical strength drop at the onset of thermal decomposition is 0.2-0.4 0 . Comparing with Figure 6 we see that larger stress drops are associated with smaller values of W.
temperatures thermal decomposition provides more weakening than thermal pressurization, and we predict that the maximum strain rate in the gouge layer increases to a new peak value before decaying, indicating that the onset of thermal decomposition drives additional strain rate localization. Our observations agree with the results for strain localization driven by thermal pressurization and thermal decomposition in an elastoplastic Cosserat material presented in Veveakis et al. [2012], which also showed additional localization at the onset of thermal decomposition.
We used a linear stability analysis to quantitatively predict the localized zone thickness as a function of the fault temperature. As expected, at low temperatures we recover the predictions from Rice et al. [2014], which studied strain localization driven by thermal pressurization alone. At high temperatures the localized zone thickness is independent of the fault temperature, and the formula for localized zone thickness has a simple form that is independent of the reactant mass fraction and the reaction kinetics. The reaction controls the localized zone thickness only through the parameters E r and P r . For fault temperatures between the highand low-temperature limits we solved for the localized zone thickness using Cardano's formula for the roots of a cubic equation, leading to a more complicated formula than the simple solution in the high-temperature limit. This formula shows a weak dependence on the reactant mass fraction and reaction kinetics and requires a current fault temperature to be specified.
We tested our analytic predictions using numerical simulations. Performing a parameter sweep over all relevant dimensionless parameters, we found that the more general cubic formula makes more accurate predictions than the simpler formula valid in the high-temperature limit. This is because the endothermic nature of the reaction limits the peak fault temperature to a value below the region where the high-temperature limit is valid. Based on this we conclude that the best predictions for localized zone thickness when thermal decomposition is active are given by equation (60). However, this means we must know the reaction kinetics and hope that the peak fault temperature is well approximated by equation (32), which is only the case if we can estimate how to offset the power densitẏto account for losses by thermal diffusion. When the reaction kinetics are unknown, a prediction for the localized zone thickness can still be made using the simpler formula in equation (49), though this systematically underpredicts the localized zone thickness observed in the numerical simulations by up to an order of magnitude.
The ubiquity of carbonates and hydrated clays in mature faults and the large temperature rises expected during an earthquake suggest that thermal decomposition is likely triggered during most large earthquakes. This suggests that it may be more appropriate to compare the predictions from equation (60) with field and laboratory observations of micron-scale strain localization than the low-temperature limit studied in Rice et al. [2014] and Platt et al. [2014]. The localized zone thicknesses predicted in this paper are in good agreement with the majority of observations of strain localization, and a detailed discussion of these observations can be found in the introduction of Rice et al. [2014]. When comparing with field and laboratory observations it may be more appropriate to use 2W to estimate the width of the localized zone, since only 68% of the deformation occurs between y = −W∕2 and y = +W∕2.
Depending on the extent of grain size reduction or amorphization due to comminution and thermal decomposition, the thinnest localized zone thicknesses predicted in this paper may be comparable to a typical grain size in the gouge layer. This means that for the very thinnest localized shear zones the size of individual grains may be an important localization limiter. There are several ways to predict a localized zone thickness in this limit, as discussed in Rice et al. [2014] and Platt et al. [2014]. One option, which is based on a wide body of research on localization in granular systems, is to set the localized zone thickness equal to ∼ 10 − 20d 50 , where d 50 is the grain size such that 50% by weight of the particles have larger size. Another option is to extend the model presented in this paper to account for the motion of individual grains. This might be done using a higher order continua or gradient theory that models the inertia of individual grains, and examples of how these models interact with thermal and pore fluid effects can be found in Vardoulakis [2002], Sulem et al. [2011], andVeveakis et al. [2012].
Our model makes many simplifications that may alter our quantitative predictions significantly, though we expect the results to be qualitatively unchanged with the localized zone thickness set by a balance between thermal decomposition, frictional rate strengthening, and diffusion. First, we assume that the gouge properties are constant and approximate the expected changes with pore pressure and temperature using the path-averaging approach from . Rempel and Rice [2006] suggested that this is a reasonable approximation for most parameters, but that the changes in hydraulic diffusivity accompanying pore pressure changes may be important. Since thermal decomposition can elevate pore pressures close to the normal stress, it is possible that the hydraulic diffusivity at peak localization is much larger than the value we assumed, leading to a localized zone thickness that is much wider than our predictions. As noted in , the solid volume change accompanying thermal decomposition will also impact the hydraulic parameters, and we expect this porosity change to increase hy and lower P r . Both of these changes will act to widen the localized zone. Since limited depletion has occurred at the moment when peak localization is achieved, we do not expect this to alter the peak localized zone thickness, but it may lead to significant widening of the localized zone as the reactant is depleted.
Equation (60) shows that the localized zone thickness depends more sensitively on f o than any other parameter in the model. This means that other dynamic weakening mechanisms that alter the friction coefficient-such as flash heating and the low friction coefficients associated with nanoparticles-may lead to localized zones that are wider than our predictions. If we crudely approximate these dynamic weakening effects by assuming a lower friction value of f o = 0.2, then we predict that the localized zone thickness will increase by almost an order of magnitude.
Our results also show that, even though the localization is controlled by spatial variations in pore pressure generated by the positive feedback between frictional heating and the two dynamic weakening mechanisms, the localized zone thickness during thermal decomposition is expected to be insensitive to changes in the ambient pore pressure. However, the ambient pore pressure will control the temperature rise during thermal pressurization alone and thus may control the strain rate and strength evolution by determining if thermal decomposition is activated or not.
One major caveat that must be attached to this work is the assumption of a fixed kinematically applied slip rate. In a dynamically propagating rupture we expect the slip rate to vary by at least an order of magnitude along the fault, with the largest slip rates at the rupture tip. Our formulas for the localized zone thickness suggest that these variations in slip rate will lead to significant changes in the localized zone thickness during an earthquake. However, Figure 3 shows that localization develops over a finite slip of a few millimeters, and thus, it is not appropriate to just evaluate equation (60) as a function of V in a dynamic rupture simulation. Properly testing the effects of a variable slip rate requires a new study that imposes V(t).
Finally, it is important to note that micron-scale localization also occurs in rotary shear experiments performed at slip rates of ∼ 10 μm/s [Yund et al., 1990;Beeler et al., 1996], and the model presented here cannot explain these observations. If another mechanism drives strain rate localization during nucleation, then it may be more appropriate to reinterpret h as the thickness of the deforming zone at the moment thermal pressurization and thermal decomposition become important.

Limiting of Peak Temperature
In addition to studying how thermal decomposition drives strain localization, we also studied the evolution of the maximum temperature within the gouge layer. This builds on previous work by  and Brantut et al. [2010Brantut et al. [ , 2011 that showed how the endothermic decomposition reaction can limit the maximum temperature rise, possibly explaining the frequent lack of pseudotachylytes on mature faults. Figure 7 shows that thermal decomposition is initially unimportant and the maximum temperature rise follows the solution for thermal pressurization alone from Platt et al. [2014]. When thermal decomposition becomes important, the maximum temperature within the gouge layer begins to rise faster than for thermal pressurization alone. This is a surprising result for an endothermic reaction but can be understood by realizing that the pore pressure generated by the reaction is driving additional localization, focusing frictional heating into a narrower zone. Eventually the reaction kinetic becomes fast enough to offset the additional heating, and we see a peak temperature followed by a gradual decay. This gradual decay is due to the strength drop that accompanies the onset of decomposition gradually lowering the total frictional heating that the reaction has to offset.
While Brantut et al. [2010] showed that the endothermic reaction caps the maximum temperature rise, they did not provide a way to predict how this temperature will change with the gouge properties or reaction triggered. In this paper we estimated the peak temperature rise by assuming it occurs when the reaction progresses fast enough to offset all frictional heating. This highlights that the peak temperature is controlled by the kinetics and is not well estimated by the temperatures from equilibrium phase diagrams. Our estimates for the peak temperature were tested using numerical simulations. Performing a parameter sweep over all relevant dimensionless parameters we showed that our estimate is generally accurate to within ∼ 50 ∘ C when we assume a fixed frictional heating equal to a 50% strength drop and a localized zone that is 150 μm wide. From this we conclude that equation (32) can be used to estimate peak temperatures when thermal decomposition is active.
These simulations also allowed us to study the role of thermal diffusion in limiting the maximum temperature. We find that, in general, thermal diffusion, which occurs rapidly for micron-scale deforming zones, is more important than thermal decomposition in limiting the maximum temperature. However, this conclusion may not extrapolate to other parameter values, and it is possible that for higher values of hy or lower values of f o , both of which lead to wider localized zones, thermal diffusion would be unimportant in limiting the peak temperature. Note that the importance of thermal diffusion contradicts the assumptions that went into equation (32), and it may be more appropriate to consider the endothermic reaction offsetting a percentage of the frictional heating when evaluating equation (32). This can be seen in Figure 8, where we found the good agreement between equation (32) and the numerical simulations by using a value oḟlarger than that observed in the numerical simulations.
It is important to note that our results are based on a large extrapolation in the reaction kinetics, and any change in A or Q will alter our results. One important physical process that is neglected here is the interaction between the pore fluid pressure and the reaction kinetics. We expect any increase in pore pressure to slow the reaction rate, which may replace the gradual decay after the peak temperature with a gradual increase.
Our predictions for talc and the illite/muscovite mixture show that thermal decomposition may not always preclude melting. However, it is likely that, on the timescales associated with seismic slip, melting is partially controlled by the kinetics, as was shown to be the case for thermal decomposition. This means that it may not be sufficient to just compare the predictions from equation (32) with a typical equilibrium melting temperature, and instead, a melting temperature should be estimated by comparing the melting kinetics with a typical seismic slip duration. Quantitative predictions for a wider range of materials is made difficult due to the lack of data to constrain the reaction kinetics.

Impact on Dynamic Weakening
Previous work by Brantut et al. [2010] showed that the onset of thermal decomposition leads to a rapid pore pressure increase and thus accelerated dynamic weakening. Our final focus in this paper was to study how the magnitude of this strength drop is controlled by the gouge properties.
As with the localized zone thickness and maximum temperature, the shear strength evolution initially follows the solution for thermal pressurization alone from Platt et al. [2014]. This means that the initial weakening follows the solution for uniform shear under undrained and adiabatic conditions from Lachenbruch [1980], and after the first strain rate localization driven by thermal pressurization, the shear strength follows the Mase-Smith-Rice slip on a plane solution Smith, 1985, 1987;. The onset of thermal decomposition is accompanied by an acceleration in dynamic weakening, leading to a lower shear strength than the Mase-Smith-Rice slip on a plane solution. While the shear strength evolution no longer follows the slip on a plane solution, the weakening rate −̇does approach that predicted by the slip on a plane solution at large slips.
Comparing the weakening rate from our numerical simulations and the slip on a plane solution, we were able to quantify the strength drop associated with the onset of thermal decomposition. Typical strength drops are ∼ 20-40% of the initial fault strength, though we see significant variations in the parameter sweep shown in Figure 12. In general, larger strength drops are associated with more intense localization, and the larger stress drops also occur over shorter slips. From this we conclude that the strength drop due to thermal decomposition is comparable to the strength drop from thermal pressurization. Assuming that flash heating can be modeled by instantaneously reducing the friction coefficient from ∼ 0.6 to ∼ 0.2 at the rupture tip, we expect flash heating to account for ∼ 70% of the coseismic strength drop with thermal pressurization and decomposition each accounting for ∼ 15% of the strength drop. However, this conclusion relies on a crude model for flash heating, and it is unclear how efficient flash heating is when deformation is distributed in a gouge material.
As discussed in section 7.1, it is important to remember that our model assumes a fixed kinematically applied slip rate. To truly determine how much of the coseismic strength drop is due to thermal decomposition requires a dynamic rupture code that couples the strength evolution on the fault surface to an elastodynamic model for the material adjacent to the fault.

Conclusions
In this manuscript we used a model for deformation in a fluid-saturated gouge layer to study seismic strain localization driven by thermal decomposition. Combining a linear stability analysis with numerical simulations, we predicted the localized zone thicknesses as a function of the fault properties, showing that when thermal decomposition dominates thermal pressurization, this thickness is set by a balance between thermal decomposition, hydraulic diffusion, and frictional rate strengthening.

10.1002/2014JB011493
In addition, we studied how the endothermic reaction combines with thermal diffusion to limit the temperature rise during an earthquake, producing an estimate for how the peak temperature depends on reaction properties. For the materials studied here this peak temperature is controlled by the reaction kinetics and is typically much larger than the equilibrium phase transition temperature.
Next we studied how the onset of thermal decomposition accelerates dynamic weakening, showing that the onset of decomposition leads to a rapid strength drop of ∼20-40% of the initial fault strength. The weakening rate after the onset of decomposition is shown to be roughly approximated by the slip on a plane solution for weakening driven by thermal pressurization, though thermal decomposition always leads to shear strengths that are lower than those predicted by thermal pressurization alone. A parameter sweep shows that larger strength drops at the onset of decomposition are associated with more intense strain localization.
Our results were used to predict the peak temperature and localized zone thickness for four different thermal decomposition reactions. We predict localized zone thicknesses between ∼ 7 and ∼ 13 μm, and peak temperatures between 885 and 1733 ∘ C. Based on these predictions we conclude that thermal decomposition drives micron-scale strain localization, but not all thermal decomposition reactions will limit the peak temperature below a typical melting temperature. We are grateful for support by NSF Geophysics Program grant EAR-1315447, with supplemental funding by the Southern California Earthquake Center. SCEC is funded by NSF cooperative agreement EAR-1033462 and USGS cooperative agreement G12AC20038. The SCEC contribution number for this paper is 1965. N.B. acknowledges the support of the UK Natural Environment Research Council through grants NE/G016909/1 (to P. G. Meredith) and NE/K009656/1 (to NB). Part of this work originated during a research visit of NB to Harvard University in 2011, and financial support from a Harvard gift fund in support of Geomechanics Research is acknowledged. We are grateful to Robert C. Viesca, Dmitry I. Garagash, and Alex Schubnel for useful discussions.