The fracture energy of ruptures driven by flash heating

We present a model for dynamic weakening of faults based on local flash heating at microscopic asperity contacts coupled to bulk heating at macroscopic scale. We estimate the fracture energy G associated with that rheology and find that for constant slip rate histories G scales with slip δ as G∝δ2 at small slip, while G∝δ1/2 at large slip. This prediction is quantitatively consistent with data from laboratory experiments conducted on dry rocks at constant slip rate. We also estimate G for crack‐like ruptures propagating at constant speed and find that G∝δ2/3 in the large slip limit. Quantitative estimates of G in that regime tend to be several orders of magnitude lower than seismologically inferred values of G. We conclude that while flash heating provides a consistent explanation for the observed dynamic weakening in laboratory experiments with kinematically imposed slip, its contribution to the energy dissipation during earthquakes becomes negligible for large events when considering the elastodynamic coupling between strength and slip evolution.


Introduction
The dynamics of earthquakes is primarily controlled by the balance between the available elastic strain energy (i.e., the prerupture stress level along the fault), the energy radiated away from the fault, and the fracture energy G consumed to advance the rupture front. The fracture energy, far from being a material constant, depends on how the fault weakens during slip, and hence is ultimately controlled by the physical processes responsible for fault weakening.
In the context of seismology and shear rupture propagation with complex friction laws, G is generally defined by the integral of the shear strength change over the local slip [e.g. Kanamori and Heaton, 2000;Abercrombie and Rice, 2005]. Therefore, G integrates potentially complex strength evolution with slip, slip rate, and other evolving physical variables. Far-field seismological observations provide constraints on the magnitude of G, which is typically derived from estimates of moment magnitude, stress drop and radiated energy [e.g. Abercrombie and Rice, 2005;Viesca and Garagash, 2015]. However, disentangling the details of stress, slip or slip rate evolution from G, is generally not possible (or at least not unequivocally) with seismological data alone [e.g. Guatteri and Spudich, 2000]. Hence, the physical mechanisms giving rise to earthquake propagation remain only accessible through a single integrated quantity, the fracture energy.
One approach to circumvent this issue and identify the underlying physics of dynamic weakening is to make predictions of G based either on empirical laboratory data  or theoretical analysis [Rice, 2006;Viesca and Garagash, 2015], and examine if and how the resulting scaling of G with other source parameters (typically, fault slip) matches with independent seismological estimates. In other words, the key question is: what does the observed scaling of G tell us about the physics of rupture ? Such an approach has been remarkably successful in identifying thermal pressurisation as a potentially ubiquitous weakening mechanism, compatible with earthquake data over a very wide range of magnitudes [Rice, 2006;Viesca and Garagash, 2015]. Thermal pressurisation is a mechanism by which faults weaken due to an increase in pore fluid pressure on the fault plane driven by frictional heating. In a purely empirical approach, Nielsen et al. [2016] have shown that the fracture energy derived from laboratory friction experiments, almost regardless of the experimental conditions, is in fact consistent with that of earthquakes. These experimental results highlight the potential nonuniqueness of the weakening mechanisms responsible for the scaling of G with slip: indeed, most of the experimental data used by Nielsen et al. [2016] were obtained on dry rocks, in a setup that essentially precludes the efficiency of thermal pressurisation.
Overall, a key question is to determine what features of weakening mechanisms are essential to reproduce the scaling of G derived from seismological or experimental data, and whether weakening mechanisms other than thermal pressurisation could also be viable candidates to explain the fracture energy of earthquakes.
Here, we tackle this issue by exploring in detail the fracture energy associated with weakening by flash heating, which is a theoretically and experimentally documented weakening mechanism occurring at the onset of seismic slip [e.g. Rice, 1999Rice, , 2006Beeler et al., 2008;Goldsby and Tullis, 2011;Passelègue et al., 2016;Brantut et al., 2016]. We first present an updated flash heating model, which includes progressive weakening due to bulk frictional heating, and then compute the associated fracture energy under either imposed slip rate or within an elastodynamic crack model. We then discuss the resulting scaling of G with slip, and compare it to experimental and earthquake data. Finally, we extract several general conclusions about how fracture energy should scale with slip for ruptures driven by thermal weakening processes.

Governing equations
The constitutive law governing frictional weakening by flash heating has been derived in detail by Rice [2006] and Beeler et al. [2008]. Here we develop a model for flash heating that is modified from the original formulation: following the steps initially outlined by Rempel [2006], and further developed by Brantut and Platt [2017] [see also Proctor et al., 2014;Yao et al., 2015], we include here the dependence of flash heating on the fault bulk temperature, and extend the flash heating model to gouge. We first recall the general form of the shear strength evolution governed by flash heating [Rice, 1999[Rice, , 2006Beeler et al., 2008]: where τ is the strength, τ 0 is the initial frictional strength of the fault, V w is a critical weakening slip rate (temperature-dependent), and V is the slip rate. Equation (1) is valid only for V ≥ V w , and τ is assumed to remain constant and equal to τ 0 at lower slip rates [Rice, 2006;Beeler et al., 2008]. Here we do not include any "residual" strength level, and assume that the strength at high slip rate approaches zero. The weakening slip rate is given by [Rice, 2006]: where α is the thermal diffusivity of the fault rock, ρc is its heat capacity , D is the asperity contact diameter, τ c is the asperity contact shear strength, and T w is the critical weakening temperature. Equation (2) states that the weakening temperature is a function of the ambient temperature T , i.e., the background temperature of asperities before they start sliding. The evolution of the shear strength of the fault is therefore given by the evolution of both the slip rate V and fault zone temperature T . The latter is a bulk average over many particles and asperities, and is governed by the heat equation: where y is the spatial coordinate normal to the fault surface, andγ is the distributed shear strain rate in the fault gouge. For a Gaussian strain rate distribution across the fault, Equation 3 has the following solution for temperature evolution at y = 0 [Carslaw and Jaeger, 1959]: where w is a measure of the shear zone thickness. The assumption of a Gaussian strain rate profile with constant w is a conservative one, since strain could further localise inside the gouge due to thermal weakening. Therefore, the temperature rise and resulting weakening are likely lower bound estimates. Using the expression (1) for shear stress, we obtain an expression for temperature that does not include any direct dependence on slip rate: In our assumption of distributed strain rate over a finite thickness (and not bare surface contact), we implicitly extend the flash heating model to an ensemble of frictional contacts distributed over the fault thickness. This generalisation has been developed by Rempel [2006] and Brantut and Platt [2017], who showed that the model would hold provided that V w is modified by a factor proportional to the number of contacts within the fault thickness.

Solutions
In order to estimate the evolution of shear stress, and therefore of fracture energy, as a function of cumulated slip, one needs to solve Equation 5. Three informative end-member solutions can be found analytically.
At early times, while the fault effective thickness w √ 2π remains large compared to the thermal boundary layer width √ αt, heating is mostly adiabatic and Equation (5) simplifies to: Combining with expression (2) for the weakening velocity, we obtain the following solution for temperature: where . The time t A w corresponds to the time required to heat a layer of thickness √ 2πw from T = T 0 up to the weakening temperature T w . Equation (7) depends only on time and is not affected by the slip rate history on the fault.
At large times, when the thermal boundary layer becomes much wider than the shear zone thickness, shear heating is essentially concentrated on an infinitely narrow width, which acts as a line source. Under those conditions, Equation (5) simplifies to: A natural characteristic time in this regime is the following: with time during flash heating for a range of timescale ratios t A w /t SP w . The evolution of temperature (a) is independent from the slip rate history. In the computation of strength we assumed V = V w0 . The black curves correspond to the full numerical solution for each ratio of timescales. The dashed blue lines correspond to the asymptotic solutions in the adiabatic regime, the dashed green lines are the asymptotic solutions for the slip-ona-plane regime in the large time limit, and the dashed orange lines are the asymptotic solution for the slip-on-a-plane regime in the small time limit. asymptotic solution, valid for t t SP w , is given by (see Supplementary Materials, Section 1) If we want to insist that the whole flash heating process occurs in the slip-on-a-plane limit, which is relevant for instance for bare rock frictional surfaces or when the two time scales t A w and t SP w are very different, we can also determine a simple asymptotic form for T (t) in the small time limit. For t A w t t SP w , we find (see Supplementary Materials, Section 1) that the temperature is well approximated by Overall, for shear over a finite thickness, we observe that the temperature, and hence the strength evolution, is controlled by only two characteristic timescales, and therefore by only one nondimensional parameter, namely t A w /t SP w . This ratio of timescales controls the dominant thermal regime of the fault zone.
A set of numerical solutions of the general problem, computed using the spectral in space, finite-difference in time method given by Noda and Lapusta [2010], are shown in Figure 1(a) for a range of ratios t A w /t SP w , along with the asymptotic solution derived above. The corresponding evolution of strength, computed using (1) and V = V w0 , is given in Figure 1(b). We observe a gradual decrease in strength over time, due to the reduction in V w (T ) induced by the macroscopic heating of the fault.

Fracture energy
Based on our strength computations, we can now make predictions for the fracture energy associated with flash heating. Here we use the generalised definition of G given by Abercrombie and Rice [2005]: where δ is the slip. Since the strength depends directly on the slip rate history, we also expect the fracture energy to do so. In the following, we analyse how G scales with slip using two models for slip rate evolution, one with constant slip rate, and one derived from elastodynamics.

Analysis using constant slip rate
In a first approximation, we use a simple assumption of constant slip rate to compute G. In this case, analytical formulae can be derived for G(δ ) in the three asymptotic cases outlined in the previous Section. In the adiabatic regime, a direct computation of (13) using δ = V t yields: In the slip-on-a-plane approximation, in the small slip limit (i.e., V t A w δ V t SP w ), an approximate form for the fracture energy is (see Supplementary Materials, Section 2) where δ SP w = V w0 t SP w . Finally, a slip-on-a-plane, large slip approximation is given by (see Supplementary Materials, Section 2): (slip on a plane, large slip), which is valid for δ V t SP w . The results for G(δ ) in the general case (computed numerically) and the approximate analytical solutions are shown in Figure 2. Essentially, we find that there is a switch from G ∝ δ 2 at small slip to G ∝ √ δ at large slip; if the two timescales t A w and t SP w are separated, an intermediate regime arises where G ∝ δ 3/2 . This behaviour is completely analogous to the scaling given by Rice [2006] for weakening by thermal pressurisation.
The reason for the similarity between flash heating and thermal pressurisation in this context is the fact that, at small slip, both mechanisms are essentially similar to linear slip-weakening; whereas at large slip, both mechanisms are dominated by a thermal weakening mechanisms controlled by a thermal (and/or hydraulic) diffusive boundary layer.

Analysis using a dynamic crack model
The kinematic approach outlined above gives initial insight into the scaling of fracture energy with slip. However, using a constant slip rate is a simplification, inconsistent with the mechanics of rupture propagation in which slip rate evolves in concert with strength behind the rupture tip.
For a semi-inifinite shear crack propagating at constant speed, the elastodynamic equilibrium requires that where x is the position from the rupture tip, V r is the rupture speed and µ * is an elastic shear modulus which depends on the mode of rupture and on the rupture speed [Rice, 1980]. The elastodynamic stress (17) has to be consistent with the strength on the fault given by the flash heating process (Equation 1). Therefore, the slip rate, stress and temperature histories are coupled and have to be determined simultaneously. Far from the crack tip, for large slip, an asymptotic analysis of the coupled system (17), (1) and (9)    energy: Notably, G scaling with a 2/3-power law in slip was also found by Viesca and Garagash [2015] for dynamic ruptures driven by thermal pressurisation.

Comparison with Laboratory and Earthquake Data
The theoretical results outlined in the previous Section can be compared to laboratory data obtained at high constant (imposed) slip rate. Figure 3 shows the fracture energy compilation of Nielsen et al. [2016] as a function slip, plotted together with the theoretical predictions for flash heating using a realistic set of parameter values [Brantut and Platt, 2017]: constant slip rate V = 1 m/s, critical time t SP w = 0.1 s, critical slip rate V w0 = 0.2 m/s, and nominal stress τ 0 = 20 MPa. The flash heating model reproduces the shift in trend, modelled here as a transition between δ 3/2 at small slip and δ 1/2 at large slip.
The theoretical results for the dynamic crack-tip problem can be used to see whether we can also explain earthquake fracture energy data with flash heating only. This is attempted in Figure 4, which shows the G vs. δ data compiled by Viesca and Garagash [2015] together with the large slip asymptote obtained from the dynamic steady-state crack analysis. We used a rupture speed V r equal to around 90% of the shear wave speed, so that the elastodynamic shear modulus µ * is reduced by approximately a factor 2 compared to its static value.
Even though the set of parameters used in the simulation are similar to that used to fit the laboratory data, the fracture energy predicted by the model remains much smaller than for earthquakes over a significant range of slip (δ > 10 −2 m). Beyond 10 mm slip, the fracture energy from flash heating is much smaller than for earthquakes, implying that other dissipation processes dominate (e.g., thermal pressurisation). This is consistent with previous estimates for the relative contribution of flash heating vs. thermal pressurisation [Brantut and Rice, 2011]. In Figure 4 we used a value for V w0 consistent with bare rock surfaces; in gouge, V w0 should be divided by a factor commensurate with the number of contacts within the gouge thickness. This modification would increase the value of fracture energy, but has little quantitative impact due to the weak sensitivity of G to V w0 (power −1/3). Furthermore, the correction of V w0 for gouge at large slip is likely small because of the potential strain localisation occurring during the early stages of slip.
At small slip distances, the earthquake data are consistent with a scaling of G with δ 2 , as for instance produced by a linear slipweakening friction law (as long as the constant residual friction is not reached) or by thermal pressurisation of pore fluids. Note that slip-weakening is not necessarily incompatible with the physics of flash heating, since there must be a critical slip distance beyond which asperities start to weaken at high speed [e.g. Noda et al., 2009;Brantut and Rice, 2011;Viesca and Garagash, 2015]. Here, we did not include any slip-weakening process in our flash heating model in order to explore the properties of the thermal weakening process alone. It appears that the involvement of an element of slip weakening at small slip distances within the flash heating framework is necessary to produce a more realistic scaling of fracture energy with slip.
Despite the similarity in the scaling of G with slip between thermal pressurisation and flash heating at large slip (G ∝ δ 2/3 ), any quantitative estimate using realistic parameter values reveals that flash heating has a negligible contribution in G.

Discussion and Conclusion
Laboratory data are explained quantitatively well by the flash heating model, but not earthquake data, especially for slip distances larger than a few millimetres. By contrast, Viesca and Garagash [2015] have shown that thermal pressurisation is able explain earthquake data over 9 orders of magnitude in slip. This difference is due mostly to the low heat diffusivity of rocks, which make critical weakening times and distances for flash heating very short compared to those linked with thermal pressurisation (see Supplementary Materials, Section 3). Indeed, the characteristic slip distance associated with thermal pressurisation is governed by the hydraulic diffusivity of fault rocks, which is widely variable and typically orders of magnitude larger than their thermal diffusivity. In addition, the large slip rates arising in the elastodynamic crack model tend to  induce a faster strength reduction than in constant slip rate cases, therefore producing overall lower fracture energies. Despite the quantitative discrepancies between flash heating and thermal pressurisation, for both processes the fracture energy at large slip scales with δ 1/2 at constant slip rate and with δ 2/3 for propagating ruptures. This similarity is not coincidental; in fact, any thermal weakening process for which temperature remains bounded at large times would produce similar scalings of G with slip. Indeed, the large slip asymptote is obtained by observing that (1) V (x) and τ(x) both decrease with same power x λ far from the crack tip (Equation 17, see Viesca and Garagash [2015] and Supplementary Materials), and (2) the integral in (4) approaches a constant, finite temperature for sufficiently large times. These requirements imply that λ = −1/4, from which G ∝ δ 2/3 is deduced.
By contrast, the apparent stronger scaling of G with δ 2 at small slip merely reflects a linear slip-weakening process. Thermal pressurisation under adiabatic, undrained conditions is a likely possibility, but may not be the only one. For instance, a regularised flash heating process including an intrinsic critical slip distance for asperities to weaken [Noda et al., 2009] would also be a possibility.
From a phenomenological point of view, brittle fracture of intact rocks is also characterised by a slip-weakening process [e.g. Ohnaka, 2003]; so does rate-and-state friction at moderate slip rates (i.e., in the absence of healing or state recovery). Any of these phenomena is compatible with seismological estimates of G(δ ) for small slip (typically δ 1 cm).
In summary, the most general conclusion that can be drawn from the comparison of friction models and seismological constraints of fracture energy is that seismic slip occurs by a succession or combination of physical processes which (1) initially resemble linear slip-weakening, and (2) progressively become dominated by diffusion across the fault. In other words, the progressive change in scaling of G vs. slip with increasing slip imply that shear work dissipation occurs more and more outside the fault, either due to thermal or hydraulic diffusion (as in flash heating or thermal pressurisation), or alternatively by off-fault damage (as explained by Andrews [1976Andrews [ , 2005; Nielsen et al. [2016]).
The theoretical developments presented here show that great care is required when comparing friction models (empirical or physics-based) to earthquake data: except for the purely slipdependent friction laws, boundary conditions in terms of slip rate history generally have an impact on the strength evolution and on the resulting fracture energy. Dynamic steady-state rupture models [e.g. Garagash, 2012;Viesca and Garagash, 2015] provide a useful tool to circumvent the shortcomings of assuming a priori slip rate histories, without having to resort to computationnally intensive numerical elastodynamic simulations.