Thermal Weakening Friction During Seismic Slip: Experiments and Models With Heat Sources and Sinks

Experiments that systematically explore rock friction under crustal earthquake conditions reveal that faults undergo abrupt dynamic weakening. Processes related to heating and weakening of fault surfaces have been invoked to explain pronounced velocity weakening. Both contact asperity temperature Ta and background temperature T of the slip zone evolve significantly during high‐velocity slip due to heat sources (frictional work), heat sinks (e.g., latent heat of decomposition processes), and diffusion. Using carefully calibrated High‐Velocity Rotary Friction experiments, we test the compatibility of thermal weakening models: (1) a model of friction based only on T in an extremely simplified, Arrhenius‐like thermal dependence; (2) a flash heating model which accounts for the evolution of both V and T; (3) same but including heat sinks in the thermal balance; and (4) same but including the thermal dependence of diffusivity and heat capacity. All models reflect the experimental results but model (1) results in unrealistically low temperatures and model (2) reproduces the restrengthening phase only by modifying the parameters for each experimental condition. The presence of dissipative heat sinks in stage (3) significantly affects T and reflects on the friction, allowing a better joint fit of the initial weakening and final strength recovery across a range of experiments. Temperature is significantly altered by thermal dependence of (4). However, similar results can be obtained by (3) and (4) by adjusting the energy sinks. To compute temperature in this type of problem, we compare the efficiency of three different numerical approximations (finite difference, wavenumber summation, and discrete integral).

analyzing the synthetic seismograms, resulted in estimates of several meters. However, as noted by Ulrich et al. (2019), the second estimate is expected to yield larger values because it implicitly includes the effects of off-fault damage in the rupture energy budget . In spite of such rare highlights, dynamic weakening remains difficult to quantify based on seismological earthquake data. Hence, a number of laws with enhanced velocity-weakening have been implemented in models of seismic fault rupture (Nielsen & Madariaga, 2003;Noda et al., 2009;Zheng & Rice, 1998), relying mostly on theoretical arguments (Archard, 1959;Beeler et al., 2008;Nielsen et al., 2008;Noda et al., 2009;Rempel & Rice, 2006;Rempel & Weaver, 2008; J. R. Rice, 2006, and references therein).
On the other hand, an increasing number of well-constrained observations are being obtained in laboratory experiments performed under close to co-seismic conditions. Prakash (2008, 2012) used an impact bar to load impulsively a frictional slip surface under extreme conditions of slip rate (tens of meters per second) and normal stress (hundreds of Mega Pascals) while measuring the shear resistance to slip; they found an abrupt weakening occurring over extremely short slip distances (<1 μm) and times (<1 μs). Intermediate, more seismic-like conditions were studied (0.5-6.5 m/s, 1-50 MPa) using rotary shear machines (Di Toro et al., 2004, 2006Fondriest et al., 2013;Han et al., 2007;Hirose & Shimamoto, 2005;Mizoguchi et al., 2007;Sone & Shimamoto, 2009;Tsutsumi & Shimamoto, 1997;Violay, Nielsen, Gibert et al., 2013;Violay et al., 2015) also resulting in pronounced weakening; however, in the latter experiments, the measured weakening distances were much longer (of the order of tens of centimeters to several meters).
In the case of frictional melting, the role of temperature and frictional power in the weakening was directly modeled, and it was shown that the weakening was accelerated under larger normal stress and slip velocity (Nielsen, Mosca, et al., 2010). Theoretical arguments (Nielsen et al., 2008) predicted that the final, steadystate friction level depended on normal stress to a power 1/4, which was later confirmed by accurate experiments (Violay et al., 2014). The center of the molten layer can reach peak temperatures well above those of melting temperatures of the rock-constituent minerals (overheating), creating an ultra-thin, ultra-low viscosity slip layer whose lubricant effect is all the more efficient when slip rate and normal stress are elevated. However previous modeling of frictional melt during the transient weakening (Nielsen, Mosca, et al., 2010) relies on a complex, explicit numerical model; in addition, the model considers cases where melting is observed (in silica-built rocks) but not the cases where no melting occurs-for example, carbonate rocks (Han et al., 2007;-although, there too, considerable weakening takes place. Thermal pressurization of fluids confined to a narrow fault zone has also been invoked as a cause for profound weakening in natural faults (Rempel & Rice, 2006; J. R. Rice, 2006). Such mechanism may take place in fluid-saturated, relatively low permeability fault zones, and modeling shows that it can be compatible with the estimates of fracture energy from natural earthquakes (J. R. Rice, 2006;Viesca & Garagash, 2015). In a few cases, high-velocity experiments were conducted with fluids under drained and undrained conditions, on bare rock samples with no gouge, showing the onset of thermal pressurization only in the later phases of the slip where friction was already low (Violay et al., 2011(Violay et al., , 2015Violay, Nielsen, Gibert et al., 2013;. However, it is observed in most high-velocity friction experiments that extremely fast, efficient weakening takes place on natural rocks even in the absence of fluids. Drawing on early studies of flash heating of asperities in metal friction (Archard, 1959), Rice (J. Rice, 1999; J. R. Rice, 2006) introduced a model for rock friction at high slip rates where temperature rise is implicit and the asperities weakening is essentially related to slip velocity; such flash-weakening model was subsequently discussed by Beeler et al. (2008) and Rempel and Weaver (2008). These works were followed by a rapidly developing body of studies on thermal weakening during fast slip, from either experimental tests on specific lithologies-for example, serpentine (Proctor et al., 2014), illite and quartz gouge (Yao et al., 2016)-or theoretical modeling perspectives-for example, thermal pressurization (Brantut & Platt, 2017) and flash weakening (Brantut & Viesca, 2017).
It has been claimed that slip acceleration (Chang et al., 2012) plays a fundamental role in dynamic friction reduction; however, combined high slip velocity and high normal stress, producing high frictional power, appear to be the key requirements to induce pronounced weakening  independently of the imposed acceleration. In fact one direct consequence of elevated frictional power is to induce an elevated and localized temperature growth on the slip surface and its immediate vicinity; accordingly high NIELSEN ET AL. 10.1029/2020JB020652 2 of 31 temperature has been indicated as a likely cause of dynamic frictional weakening. On the other hand, the direct effect of temperature on the weakening may be questioned. Experiments conducted on preheated samples of dolerite (Noda et al., 2011) with the use of a furnace, and on Westerly granite, India gabbro, and quartzite (Noda et al., 2011;Passelègue et al., 2014), reveal a correlation between weakening and temperature, independently of slip velocity, although the weakening is relatively modest. Frictional weakening of pre-heated olivine samples at slow slip rates is even more modest (King & Marone, 2012). We shall discuss these results here in terms of localized versus background temperature changes.
Here, we discuss aspects of frictional contact and the effects of heating under high-velocity sliding. We define a model based on the flash-weakening formalism discussed above (Archard, 1959;Beeler et al., 2008;Rempel & Weaver, 2008;J. Rice, 1999; J. R. Rice, 2006), and extend it to include the effect of frictional heating on the background temperature, thermal diffusion, and the presence of heat sinks due to decomposition or melting. We propose an extension to account for the onset of frictional melting and shortening. We indicate a parsimonious numerical scheme for the temperature update. Assumptions and approximations are used to obtain a sufficiently elementary and uncomplicated model for practical use as a friction law in earthquake slip models, while retaining the essential behavior observed during rock test experiments under coseismic conditions. Temperature evolution is important in the weakening behavior, but its accurate numerical evaluation can be costly and inefficient over an extended number of time steps. Such cost may become prohibitive when the temperature needs to be evaluated at many different points in an extended fault model. Here we test three different temperature computation schemes and compare the efficiency and the flexibility of each (see Appendix A). First, the temperature resulting from an imposed boundary flow can be written as an analytical integral, which can be directly discretized and solved numerically. This is easy to implement, but by far the least efficient approach, and it does not allow to include variations of the parameters in time or in space. Second, a spatial Fourier transform of the diffusion equation can be written, and a discrete wavenumber equation (Noda & Lapusta, 2010) can be solved numerically. Accuracy within a few percent can be achieved with a small number of memory variables (16 or less) which are updated at each iteration and can be augmented to an arbitrary level by increasing the number of memory variables. The scheme is easily implemented and can be integrated with existing numerical codes. It is far more efficient than the direct summation, but in this method, too, it is not easy to include parameter variability in time and space. Finally, we implement the classic finite difference method with a Crank-Nicolson scheme, where variables are defined on a grid of nodes in space and time. In conditions of frictional heating, this method to solve for temperature diffusion is the most flexible, and-from our tests in either Python 3.7 or Fortran 90 on a desktop computer-it can achieve equal accuracy in less time than the two other methods, if correctly configured. However, it is slightly more costly in terms of memory than the wavenumber method. In addition, it allows to include a moving boundary to solve the Stefan problem (Carslaw & Jaeger, 1959;Nielsen, Mosca, et al., 2010) in the presence of melting, and allowing to introduce spatial and time variations in the parameters, as illustrated in Section 6.
The models are calibrated and tested against a number of selected rotary shear, high-velocity experiments performed on the SHIVA machine hosted at Istituto Nazionale di Geofisica e Vulcanologia, Roma. We used hollow cylindrical samples of 30/50 mm inner/outer diameter, respectively, machined from gabbro (as a representative of silicate-built rocks) and Carrara marble (as a representative of carbonate-built rocks). The experimental conditions cover the range from 10 to 30 MPa in normal stress and from 1 to 6.5 m/s in slip velocity.

Heat Sources and Sinks
The frictional work released per unit time (power) per unit area on the sliding interface is where V is the slip rate and τ is a macroscopic average of the shear stress. Frictional work is responsible for temperature rise, but it is in part dissipated by thermal diffusion and in part by other endothermal processes (latent heat for melting, decomposition, heat removal by fluid mass escape from the interface, surface NIELSEN ET AL.

10.1029/2020JB020652
3 of 31 energy involved in comminution, etc.). The latter processes are known to act as a buffer which inhibits the continuous rise of temperature; we argue that they have a significant effect on the background temperature and friction.
Assuming that the frictional heat rate from Equation 1 takes place on the fault surface (or within a principal slip zone of negligible thickness) at z = 0 and propagates away from the fault surface, we may write the one-dimensional thermal diffusion equation We delay to Section 5 the discussion of particular forms of Equation 4 including thermal dependence, and their compatibility with observed experiments of high velocity friction, and discuss here the nature of possible heat sinks.
The temperature-limiting effect of thermal decomposition has been modeled explicitly in the case of carbonate rocks (Sulem & Famin, 2009), gypsum (Brantut et al., 2010), and dolomite with application to the emplacement of a giant landslide (Mitchell et al., 2015). Generally, the kinetics of a reaction is accelerated exponentially with temperature (Arrhenius kinetics), as a consequence the rate of latent heat loss should increase likewise. In the case of experiments performed in the open air or in a water-filled vessel, there is an additional heat loss (Newton's law of cooling, generally considered as proportional to the temperature difference between a body and the surrounding fluid), which is enhanced by the presence of fluid convection (Acosta et al., 2018;Violay, Nielsen, Gibert, et al., 2013).
On the other hand, if the thermal dependence of diffusivity and heat capacity is considered, an increase in temperature can result in a decrease in thermal conductivity (Merriman et al., 2018). This will increase the insulating property of the rock (at least locally in the rock layer where temperature rise is substantive) inducing a feedback that enhances the temperature rise. This effect can then counteract the action of the heat sinks. This competition between the two mechanisms is further investigated in Section 6.
Heat loss due to a constant flow rate of cooling fluid may be approximated as q s = ψ (T − T i ) per unit time, assuming that fluid enters the interface at T i (ambient rock temperature) and exits at temperature T (background interface temperature); ψ is the heat capacity of the fluid times its flow rate (per unit area).
Heat loss due to decomposition processes will be represented by an exponential Arrhenius law of the form (Sulem & Famin, 2009): Here h is the thickness of the zone affected by the decomposition, (1 − n) is the remaining proportion of (unreacted) material (in first approximation 1 − n ≈ 1), and R = 8.31JK −1 mol −1 . Note that in Equation 5 T is the absolute temperature, but in the following sections T is the excess above initial temperature. In the example of decarbonation of pure calcite, indicative literature values (Sulem & Famin, 2009) are E a = 319 10 3 J/mol, A = 2.95 10 15 s −1 , L = 3190 10 3 J/kg, for activation energy, pre-exponential factor and latent heat, NIELSEN ET AL.
We note that the term q s (T(t)) in Equation 3 implicitly describes heat sinks that are distributed over a finite thickness h. Distributed heat sinks can be explicitly included using the values of temperature and temperature gradient away from the sliding surface which are derived in Equations A14-A15. However here for simplicity, we will assume that the temperature over h can be equated to that of the sliding surface T = T(z = 0) and that h is the constant; a similar approximation was also used in Mitchell et al. (2015).
In the case of frictional melt, the latent heat and melt extrusion have been taken into account explicitly to model temperature evolution and to describe frictional behavior in both steady-state (Nielsen et al., 2008) and transient conditions (Nielsen, Mosca, et al., 2010). In the case of pervasive melting, it is necessary to solve the Stefan problem with an added term of mass transport in Equation 2, a case also discussed further in Section 6.1 and in Appendix B.

Real and Nominal Contact Area
Shear and normal stress across the sliding interface are supported by local asperities whose real contact area A r represents only a small fraction α of the nominal area A n . During slip, asperity contacts coalesce, deform and disappear forming a distribution at various stages of evolution and under continuous renewal. Because rock constitutive minerals yield under a few percent of strain, within each asperity the shear stress reaches the yield point early during its contact lifetime.
After yielding each asperity deforms under either brittle failure or creep, possibly at a very high strain rate; however, it cannot support a stress value much in excess of the yield value, lest yielding would propagate from the asperity into the supporting substrate thus keeping the stress value bounded. Hence, the majority of asperities are close to the yield shear stress τ y which can be considered as an average asperity value. The bulk frictional force resisting slip can be written as F = A r τ y which, normalized by A n , yields the bulk frictional stress: α results from the ratio of applied normal stress σ n to indentation hardness σ c (or penetration hardness according to Persson, 2000, Chapter 5.1) such that but in any case α ≤ 1 so that we may write α = Min(1, σ n /σ c ).
We remark that both indentation hardness (Atkins & Tabor, 1965;Hirth & Kohlstedt, 2004;King & Marone, 2012) and yield shear stress (Raterron et al., 2004;Weidner et al., 1994) are observed to have a strong negative temperature dependence. Accordingly, an increasing temperature induces a decrease in τ y but an increase in α-see also discussion in (Hirth & Beeler, 2015, and references therein). An increase in the area of asperity contact has been documented, for example, in olivine under slow slip velocity (King & Marone, 2012). Since the thermal weakening and contact area increase have antagonistic effects, the temperature dependence of bulk friction τ is not trivial to predict. In fact experiments performed under slow slip velocity do not show a systematic or pronounced frictional drop with temperature (King & Marone, 2012;Noda et al., 2009).
However, under high slip velocity, a pronounced weakening is observed in correspondence of the temperature rise at the interface. First, we shall propose how to reconcile these two conflicting observations based on the role of slip velocity and temperature localization. Then, we shall proceed to the computation of the bulk frictional resistance of the interface based on local, temperature-dependent rheology.
NIELSEN ET AL.
10.1029/2020JB020652 5 of 31 In case that a pervasive lubricant layer develops and fills continuously the space between the asperity contacts (e.g., a pervasive melt layer  and supports the bulk of shear and normal stress, the temperature effect on α is buffered as we may consider α ≈ 1. At this point, τ y = τ is the viscous shear stress supported by the lubricant layer within a principal slip zone (PSZ), and resistance to sliding is due to the viscous shear of a thin melt layer. Though the heating is not localized at the asperity contacts, it is still localized and concentrated within a thin shear layer provided that slip is brief enough (earthquake-like duration, typically seconds) that heat diffusion away from the PSZ is reduced (close to adiabatic conditions). In a different context (no melting), Cornelio et al. (2019Cornelio et al. ( , 2020 have shown how viscous fluids permeating natural rock samples affect frictional weakening at high slip velocity by activating elastohydrodynamic lubrication.

Shear Thickness and Shear Rate
Observations on paleoseismic faults which were active at epicentral depths (≈10 km Di Toro et al., 2005) show that slip often localizes within a PSZ of a limited thickness (of the order of 100 μm or less). Active or fossil faults at moderate depths also often exhibit localized principal slip zones within a wider fault core (De Paola et al., 2008;Collettini et al., 2011;Otsuki et al., 2003;J. R. Rice, 2006;Sibson, 2003 and references therein).
Finally, laboratory experiments conducted under high stress and velocity also report the development of extremely thin PSZs either between two consolidated rock samples or within simulated or natural fault gouge. In the latter case, localization is achieved only after a critical slip of several centimeters (Pozzi et al., 2019(Pozzi et al., , 2018Smith et al., 2015). The strain rate can be equated either to the ratio of slip velocity to the thickness of the PSZ, in the presence of a pervasive lubricant layer, or in the absence thereof, to the ratio of the slip velocity to asperity height (typically ≈ 10-100 μm).
Since average seismic fault slip velocity estimated during earthquakes is typically V ≈ 1 m/s, the resulting shear strain rate within the PSZ or within the asperity contacts is extremely high (     4 4 / 10 10 V s −1 ) and is associated with a number of thermally triggered decomposition, alteration or amorphization processes (dehydration, melting, decarbonation, stress corrosion, comminution, etc.) which may directly or indirectly affect friction through the action of pressurization and/or the formation of a lubricant layer (Hirose & Shimamoto, 2005).
The consequence of such extreme localization of strain rate is that a very low equivalent viscosity should be achieved within the PSZ to result in the observed low value of shear stress during fast sliding. In addition, if the shear heating is generated within a very small thickness   4 h t (where t is the duration of the slip, indicatively 1 s), the temperature evolution is better approximated by a zero-thickness model than by the adiabatic solution within the PSZ (J. R. Rice, 2006).

Arrhenius Thermal Dependance and Flow Stress
While little is known about rock rheology at large strain rates, stress relaxation occurs through any of i different crystal plastic mechanisms (dislocation diffusion, grain boundary migration, etc.) which generally obey standard Arrhenius thermal dependence with an activation energy Q j and a power dependence on stress such that where the constants A i may include a grain size dependence for particular deformation mechanisms (e.g., diffusion plasticity). Within a given range of stress and temperature, we assume that a single mechanism will dominate and invert it to obtain NIELSEN ET AL.

10.1029/2020JB020652
6 of 31 where σ,  are the stress and the shear strain rates, respectively and T is the temperature. The exponent may be as low as n = 1 for some purely diffusive processes (Nabarro-Herring Poirier, 1985) but in most cases n > 1; for example, n = 2 for grain boundary sliding (Karato, 2008) and typical values 1.5 < n < 3 are observed at    3 4 10 10  in experiments on ceramics in brittle conditions (Lankford, 1996). The term C (Pa s 1/n ), though considered constant here for simplicity, may be strongly dependent on grain size (e.g., in the case of grain boundary sliding) among other parameters (Violay et al., 2012). A consequence of n > 1 is that the  1/ n term does not vary greatly under extremely elevated values of strain rate (e.g., the term  1/3 varies of about 25% upon a twofold increase of slip velocity from 1 to 2 m/s, assuming a shear zone of 100 μm). On the other hand, expected temperature changes of a few hundred degrees may induce a huge variation in the exponential dependence. As a consequence, under high slip velocity we can expect that the variation of τ y is primarily due to temperature changes and, for simplicity, we may neglect the variability of  1/ n to write where τ a is a reference stress T c = Q/n R and α is the real contact area and τ 0 = ατ. Here, T c is an absolute temperature (in °K). The term T i has been added to the denominator, because in the following section T is the background temperature rise with respect to the initial temperature (in the case of the experiments this will be the room temperature T = 293°K).

Local and Background Temperatures
The growth of α with temperature (due to thermal weakening of σ c ) may reduce or even surpass the weakening due to the right-hand term in Equation 11, which may explain the results that no significant weakening is observed under low slip rate even at high temperatures (King & Marone, 2012;Noda et al., 2009). However, the growth of the contact area is controlled mainly by the formation of new asperity contacts (as opposed to the growth of pre-existing ones Persson, 2000), and involves a sensibly deeper-rooted strain in a larger volume than the immediate layer below asperity contact. Consequently, the increase of α should be sensitive to the average increase of the interface temperature, or background temperature T, while the weakening of τ y will be sensitive to the local, transient temperature peak T + ΔT reached at the asperity contact, where ΔT is the transient temperature increase in the asperity ( Figure 1b). Provided that the slip rate is high, frictional heat has little time to diffuse away from the asperity during its limited contact lifetime; therefore, the excess temperature ΔT will be significant (Archard, 1959;J. R. Rice, 2006). However with continued slip, the local overheating ΔT starts to diffuse away from the asperities and gradually contributes to the rise of the background temperature T.
If the background temperature continues to rise ( Figure 1c) at some point the formation of a pervasive layer of amorphous material, melt, wear product, or viscous nanocrystalline material may fill the interstitial space between asperities, as a saturation value of α ≈ 1 is reached. In case of a continuous lubricant layer, under high velocity the thermal gradient in the vicinity of the slip zone is very steep and an extremely thin (<100 μm), overheated, and low effective viscosity layer develops ( Figure 1c). This situation has been documented in the case of frictional melt (Nielsen et al., 2008) and in the case of coseismic viscous flow in coseismic ultramylonites (Pozzi et al., 2019).
Consequently, we may consider that under high slip velocity, the growth of α is initially negligible (i.e., models of flash weakening acting in the very early stages of slip), but gradually increases with the accompanying growth of a pervasive lubricant layer, in which case the lubricant effect will compensate the increase in the contact area and the friction will not increase.
However, the transition from an initial flash heating to a fully developed lubricant layer can be complex and non-monotonic, especially at low normal stress. A relative restrengthening can occur due to the increase of the contact area ratio α with temperature rise, and the straining and elongation of the contact asperities, while voids are filled by-products of comminution, decomposition or cool melt. Microstructures corresponding to such stages were described to some extent for experiments on Gabbro under increasing slip amounts, see Hirose and Shimamoto (2005  (a) Slow velocity: the heat diffuses beyond the asperity size during the lifetime of a contact, the temperature rise is homogeneous within (h) The temperature T b is evenly distributed, so that weakening is compensated by the growth of the contact area. High velocity: (b) The asperity undergoes a local, transient temperature rise ΔT which diffuses within a limited thickness h during the short contact lifetime; on timescales of multiple contacts, the heat diffuses throughout thickness H and background interface temperature rises to T b . Temperature is unevenly distributed; the weakening due to highly localized temperature T b + ΔT surpasses the friction increase due to the growth of contact area under an average interface temperature T b . (c) A pervasive layer of overheated, lubricant material has formed with peak temperature T b + ΔT within the layer, and weakening is efficient.

Temperature Computation
We provide a number of indications that background temperature plays an important role in the weakening (Section 4). Therefore, its computation is of paramount importance in problems of thermal weakening.
The problem of heat diffusion is well known and a number of analytical and numerical solutions have been proposed. However, not all numerical solutions are equally effective, and inefficiency can be limiting in a problem where temperature needs to be computed repeatedly at many different points (as in the case of dynamic rupture models). For the type of frictional heating problem at hand, we compare the performance and advantages of each of three different numerical solutions of temperature diffusion: (1) discrete summation resolution of the integral solution, (2) a wavenumber formulation, which solution consists in the update and summation of a small number of memory variables, and (3) a finite difference, Crank-Nicolson scheme. Details about the methods and compared performance are found in Appendix A.

Weakening and Temperature Change
The particular scaling of weakening with friction, slip velocity, and frictional power systematically observed in high velocity friction experiments is compatible with a thermal signature. If temperature is the culprit, weakening should be achieved after a slip distance scaling as u c ∝ 1/(τ 2 V) and after slip duration scaling as t c ∝ 1/(τ V) 2 , as argued in Nielsen, Di Toro, and Griffith (2010). Indeed for an indicative constant value of shear stress and slip rate, an indicative solution of A1 with constant τV yields a temperature rise    c c c T V t after a time interval t c , and solving for time yields To illustrate this, we may compare two similar experiments conducted on marble in Figures'2a and 2b, where average frictional power  V differs by about a factor of two, and average  2 V by a factor of six. A similar drop to 1/3 of the initial shear stress is achieved in the two experiments at slip distances which differ by about a factor of six. Similarly, the experiments indicate that weakening time t c scales roughly as the inverse power squared (τ V) 2 as predicted.
An upper bound temperature can be obtained using Equation A1 with no heat sinks, and ρ = 2,700.0 kg/ m 3 , c = 833 J/K and κ = 0.821 10 −6 m 2 /s, for mass density, heat capacity, and thermal diffusivity of marble, respectively (Merriman et al., 2018). The estimated T curves are represented in Figure 2, showing that a comparable temperature rise (T ≈ 190°C) is achieved in both experiments at an equivalent weakening stage. These simple scaling relations seem to reinforce the idea that background temperature T exerts a strong control on the weakening.
However, at a time where weakening is already pronounced, the background temperature is still much too low (T ≈ 210°C) to trigger melting or decomposition processes (a lower bound ≈ 570°C is indicative of calcite decomposition). A similar observation can be made for the weakening of gabbro (Figure 3), with the example of experiment s555 where pervasive frictional melt formed at advanced stages (t > 1 s). Substantial weakening is observed much earlier (t ≈ 0.15 s). Using ρ = 3,000 kg/m 3 , c = 715 J/K and κ = 1.1 10 −6 m 2 /s for gabbro (Miao et al., 2014), the background temperature estimate is still only ≈ 200°C after weakening to 1/3 of the peak. The expected bulk melting temperature, about 1200°C, is achieved only later in the experiment.
Efficient weakening occurs at consistently lower background estimated temperatures than those expected to induce weakening of the material through decomposition or melting. However, local intensification of heating ΔT well above that of the background temperature T can be achieved if stress is concentrated on a fraction of the contact area only. This concentration mechanism forms the basis of the asperity flash heating (Archard, 1959; J. R. Rice, 2006), one of the frictional models discussed and tested below. In the coming sections, we will define T and ΔT and demonstrate that both are of fundamental importance in the weakening process.
NIELSEN ET AL.

Weakening of Nitrogen-Cooled Marble
In Figure 4a, we compare two experiments s876 and s880, performed on solid calcite (Carrara marble) under identical conditions, except that in s880 the rock sample was immersed in liquid Nitrogen for several seconds immediately before the experiment. This is a simple experimental test to verify that the background temperature difference has an effect in line with the prediction of basic dimensional analysis of thermal weakening. Alternatively, the temperature of the sample may have been raised before the experiment, but cooling poses lesser technical difficulties.
Our indicative estimate is that the temperature is −140 < T i < −40°C within 0.5 mm of the sample surface, at the beginning of the experimental sliding. Our estimate assumes that Newton's law of heating applies (with a poorly constrained heat transfer coefficient for air 2.5-25 W/m −2 K −1 ) within the 30 s elapsed between extraction of the sample from the nitrogen and the beginning of the experimental sliding. As a consequence, the initial temperature is lower for s880 than for s876, since the latter was initially at room temperature. If background temperature plays a role in the weakening, we expect to see some delay in the weakening for s880, which we may estimate as follows. Let T c be the temperature rise achieved in s876 after sliding for about t = 0. 70 Pa m/s V . A factor of six reduction in the weakening distance corresponds roughly to a factor of six increase in  2 V , which is expected if weakening is related to temperature increase. Similarly, the factor of two reduction in average power  V results in a factor of four increase in time (0.21 and 0.05 s, respectively) to reach equivalent weakening. A similar background temperature increase (T ≈ 190°C) above initial ambient temperature (T = 20°C) is estimated in both experiments at 1/3 weakening; however, a background temperature T ≈ 210°C is much lower than that expected to trigger weakening by decomposition reactions in calcite (about 600°C). See text for further details.
parameters for marble as in the previous section, and equating T c for both experiments, we obtain , and a delay t−t′ ≈ 0.06 s, which is roughly the delay observed in the experiment. A more accurate computation of the weakening, including the full evolution of T is shown in Figure 4b, yielding similar delay times (the full model is developed in further sections).

Test of Thermal Weakening Models
We test the fit of experimental data with two thermal weakening models: a direct Arrhenius thermal dependence, and a model flash weakening (J. R. Rice 2006) which includes both the evolving background temperature and heat sources and sinks. We discuss the differences between flash weakening and frictional melting, and test both thermal weakening models for either.
Two end-member lithologies were tested here, a calcite rock (Carrara marble) and a silicate build rock (microcrystalline gabbro). These show quite different behavior under frictional heating, as gabbro will undergo profuse melting past the initial stages of slip, but not the marble. Therefore, the micromechanics of friction are similar in the initial part of slip, but differ more widely in the later phases, in particular the recovery during the deceleration phase. The flash weakening law (Archard, 1959; J. R. Rice, 2006) was initially proposed as a kinematic treatment and did not consider the frictional recovery. Even so, the inclusion of background temperature allows to reproduces reasonably well the recovery in marble, but it over-predicts the recovery in gabbro. Recovery of marble has been considered in the context of thermal-dependent, diffusion creep plasticity . Recovery of gabbro has been analyzed in a full model of frictional melting (Nielsen, Mosca, et al., 2010), but we propose here a simplified alternative model which follows an Arrhenius thermal dependence on background temperature.
All numerical replications of the experiments are performed by imposing the experimental slip velocity history and the peak stress (i.e., static friction coefficient times normal stress). The temperature is revised at each time step based on the shear heating power as a heat source, and both thermal diffusion and dissipative heat sinks. (Considering that the friction model should be predictive, the heat source is based on the computed shear stress, not the experimentally measured shear stress.) Friction laws are based on either two or three parameters, as indicated in the text and the figures. For the solution of temperature diffusion, we include heat sinks due to endothermal phase transitions which follow (6). Rock properties (κ, c, ρ) are fixed for a given lithology (see Table 1), save for the case where thermal dependence is included for κ, c (Section 6).
We find combinations of frictional parameters and heat sink parameters T s , C s that provide a reasonable fit to the mechanical data (friction) by trial and error for each lithology. In the case of experiments s308 and s324, we conducted a systematic search on a grid of T w , B values in the flash weakening model. For each gridpoint, we computed the model friction curve and its misfit (sum of squared differences for all time steps) with the experimentally measured curve, selecting the smaller misfit value and verifying that it did provide a reasonable fit. We repeated the grid search for different combinations of T s , q s , including the case where q s = 0 (no heat sink). A trade-off is observed for the combined values T w , B, on one hand, and the combined values T s , C s , on the other hand. (See grid search example in Figure S3). This results in a reasonable fit for models over a range of parameters, only a subset of which is presented here. NIELSEN ET AL.

Background Temperature Only: Arrhenius Dependence
From the discussion and the examples of Section 4.1, it is clear that (1) weakening precedes any substantial rise of background temperature; however, (2) the background temperature still plays an important role in the weakening. We first ask the question of how a direct temperature dependence alone is capable of fitting the data where T is the background temperature.
Using Equation 11, a straightforward Arrhenius dependence for the shear stress requires adjustment of two parameters T c and τ 0 . Two additional parameters T s , C s are introduced to account for heat sinks (Equation 6) due to endothermal phase transitions. The shear value τ used in the model is the smaller of either that obtained from Equation 11 with the current temperature value, or that of the peak stress τ p = μ f σ n (where μ f is static friction coefficient and σ n is normal stress).
(1) Case of melting ( Figure 5). Parameters used for the Arrhenius thermal weakening law (Equation 11) are T c = 2700°K and τ o = 7 10 4 MPa. The parameters used for the heat sink (Equation 6) are C s = 40,10 6 Wm −2 and T s = 2000°C. We note that a strong trade-off exists between τ o and T c in the Arrhenius dependence, whereby increasing τ c can be compensated by lowering T c to achieve a very similar result; therefore, such values are purely indicative. The same remark applies to the trade-off between C s and T s . Equation 11 with a single mechanism captures some essential features of fast-slip weakening. However to accurately capture both the very initial weakening and the recovery phase, a combination of strong thermal dependence and large heat sink need to be included, to the point that the background temperature rise remains unrealistically low. This suggests that a simple model of viscous shear heating is not a realistic description of the microscale process. Alternatively, in the case of melting, the heat sink q s would implicitly incorporate the effect of heat removal by extrusion, which is not accounted for explicitly in this model-as further discussed in Section 6.1. The effect is that the computed temperature T is biased toward lower values.
(2) Case of no melting. We now use the same Arrhenius weakening model but in an attempt to reproduce the case of Carrara marble. A similar approach was adopted in Violay et al. (2019) where plastic deformation of calcite within a thin layer was assumed to follow an Arrhenius-like thermal dependence. In Pozzi et al. (2019), the steady-state friction in calcite was also interpreted in similar terms, and microstructural evidence of plastic flow was provided in support of this high-velocity deformation process. Parameters used for the Arrhenius thermal weakening law (Equation 11) are T c = 2000°K and τ o = 3 10 4 MPa. The parameters used for the heat sink (Equation 6) are C s = 50,10 6 Wm −2 and T s = 2000°C. As seen in Figure S1 the recovery of friction during the deceleration is severely underestimated, in addition, the temperature again is unrealistically low.
While the Arrhenius dependence is capable of reproducing the main features of frictional melting, it is more problematic to use it in the case of marble where no melting occurs. This is not altogether surprising, as the Arrhenius model assumes the shearing of a layer with temperature-dependent viscosity, a situation well adapted to frictional melting. However, NIELSEN ET AL.  Note. κ, ρ, c are thermal diffusivity, mass density and heat capacity, respectively.

Table 1
Rock Parameters Used Throughout the Paper Unless Otherwise Indicated in the case of marble, a model of flash weakening is likely to occur in the initial part of the slip, before the development of a continuous, high-temperature layer of low-viscosity material which can be either melted, as observed in silicate rocks (Nielsen et al., 2008), or not, as observed in carbonate rocks (Pozzi et al., 2019).

Flash Weakening With Background Temperature Evolution, Heat Sources and Sinks
We explore here a model of FWSS of heat included in the temperature estimation. Flash weakening and heating of contact asperities have been proposed as a model for high-velocity friction evolution (Archard, 1959;Beeler et al., 2008;Rempel & Weaver, 2008; J. R. Rice, 2006). There are strong experimental indications (Acosta et al., 2018;Chen & Rempel, 2014;Goldsby & Tullis, 2011;Tisato et al., 2012;Violay et al., 2014Violay et al., , 2015Violay et al., , 2011Violay, Nielsen, Gibert et al., 2013;) that this model is relevant for high-velocity experiments, in both silicate-and carbonate-built rocks, at least in the first millimeters of a slip or until melting or decomposition of the rock minerals creates an almost continuous, amorphous interstitial layer. One motivation to explore flash heating is that weakening precedes the substantial rise of the background temperature of the sliding interface (as discussed above in connection to Figures 2 and 3). Initial thermal weakening may be achieved only if local temperatures T + ΔT at asperity contacts are much higher than the background temperature T.
The FW model considers that the lifetime of asperity of linear dimension D is indicatively t c = D/V. For an asperity sheared under incipient yield stress τ c , the heating results from frictional power τ c V. Assuming that heat diffusion is mostly perpendicular to the fault, during the asperity lifetime, solution of Equation A1 with q ≈ τ c V = const. yields the local temperature rise The average strength of an asperity contact during its lifetime will be τ a = (τ r (t c − t w ) + τ c t w )/t c , where τ r is the residual shear stress supported by the weakened asperity. Assuming an asperity population with dominant dimension D, using τ p = α τ c , τ w = α τ a and noting that found (Beeler et al., 2008;Rempel & Weaver, 2008;J. R. Rice, 2006) that the effective sliding shear stress is for V > V w . The flash weakening friction is adjusted with the three parameters B (°C −2 ms −1 ), T w (°C) and τ w (Pa). The shear value τ used in the model is the smaller of either that obtained from Equations 12 and 13 with the current temperature value, or that of the peak stress τ p = μ f σ n .
In previous models, the variation of the background temperature T is often neglected in Equation 12, with the consequence that V w remains constant (Goldsby & Tullis, 2011;Noda et al., 2009). Indeed the direct numerical computation of T with classical methods can be rather costly and inefficient. However, T increases substantially after the first millimeters of slip as shown in Figure 2, and unless evolution of T is included, flash weakening fails to reproduce accurately the friction recovery observed in the experiments. One immediate evidence that friction is not purely velocity dependent is the lack of symmetry in the acceleration and deceleration phase, whereby an hysteresis loop is observed-see, for example, the τ versus V representation in Figure 6d, and also experiments reported in previous studies (Goldsby & Tullis, 2011;Proctor et al., 2014). Thus, inclusion of background temperature, which is substantially higher in the recovery phase than in the weakening phase, allows to moderate the velocity effect by acting as a state variable. This allows to obtain a hysteresis cycle where the initial weakening and the final recovery are not symmetrical, and are not purely NIELSEN ET AL.

10.1029/2020JB020652
14 of 31 velocity dependent. The higher temperature at the end of the experiment allows the friction recovery to be relatively slower, as observed.
In addition, we note that an accurate evolution of T should include both heat sources (frictional power τ V) and any significant heat sink (other than diffusion). Heat sinks are not triggered at the beginning of the slip, when the background temperature is still low; thus, they do not affect the initial temperature rise and weakening. However, as slip advances and more heat is produced, the heat sinks allow to buffer the temperature rise (Brantut et al., 2010;Sulem & Famin, 2009). We show that the temperature modulation induced by heat sinks has a sensible effect on the recovery of friction during slip deceleration.
The FWSS law is based on the full set of Equations 3, 6, 12, and 13; the only input variable is the slip velocity V(t). Input parameters are (1) the group   2 2 c D (with dimensions [L T −1 ], allowing the definition of a characteristic velocity V w ), (2) the peak stress τ p (which may be predicted using τ p = μ σ n ) and residual stress τ w (for computation of τ), and (3) C and T s (for a single dominant heat sink).
One of the peculiarities of the model described by Equation 13 is the absence of explicit dependence on normal stress. However, taking into account the evolution of background temperature T, with τ V as a heat source implicitly includes normal stress. Indeed during the initial part of the slip τ = τ p = μ s σ n where μ s is the initial friction coefficient (of the order of 0.6 before the onset of weakening), so the heat production rate is higher if the normal stress is higher. On the other hand, if the initial (peak) stress is higher under higher normal stress, temperature rise and weakening will be accelerated by a similar proportion. As a consequence, the weakening slip distance and the fracture energy may not be significantly altered by a change in normal stress. This behavior was indeed observed in a synthesis of different high-velocity friction laboratory experiments .
In Figures 7-9    in four different experiments in Figures 7-9, and in most cases the weakening, steady-state, and recovery are all reasonably well reproduced with a single set of parameters. The parameters for flash weakening (Equations 12 and 13) are T w = 800°C, τ w = 0.5 10 6 Pa, B = 0.85 10 − 6°C −2 ms −1 . The parameters used for the heat sink (Equation 6) are C s = 3 10 6 Wm −2 and T s = 1700°C. One exception, though, is experiment s324 which was performed under the most extreme frictional work rate (highest normal stress and slip velocity combination). In this case the heat sink parameter C s had to be increased by 50% (C s = 4.5 10 6 ) to reproduce correctly the recovery. One possible interpretation is that loss by combined radiation and Newton's cooling is substantial in this experiment, due to the large rate of heating, introducing additional sinks.
If the heat sinks are excluded, the final temperature is much larger (Figure 11c). This reflects on the friction, in particular, on the recovery phase which is underestimated if the same frictional parameters are maintained (Figures 11a and 11b). Note that an equally reasonable fit can be found without heat sink, but by adjusting frictional parameters for each individual experiment (Figure 11c), with a loss of generality in the model. Indeed, we could not find a single combination of parameters that would fit different experiments for a given lithology, unless we do include heat sinks in the model. Therefore, we posit that the accurate reproduction of a range of experiments, including the recovery phase, can only be achieved with background temperature evolution due to both sinks and sources. A similar observation applies to models including the thermal dependence of κ, c (see Section 6), where a reasonable fit can be obtained both with and without heat sinks provided that the frictional parameters T w , B are modified for each individual experiment.

Effect of Thermal Dependence of Diffusivity and Capacity (κ, c) on Temperature and Friction
We use experimental data from Merriman et al. (2018) on Carrara Marble to derive an empirical dependence on temperature of diffusivity and heat capacity (see supporting information). Although diffusivity and heat capacity are affected in opposite ways by the temperature rise, their combined effect still results in a significant net decrease in the conductivity (k = ρ c κ) at relatively high temperatures. However, during high-velocity frictional sliding, the high temperatures are usually reached within a small boundary layer, with a strong negative temperature gradient, because the duration of the slip is short and the diffusion distance is limited. Therefore, it is difficult to predict how important the effect of thermal dependence can be in this context, if not performing a full computation ( Figure 10).
To do so requires to take into account the inhomogeneous temperature as a function of depth z from the sample surface, and, therefore, the inhomogeneous spatial distribution of the conductivity/diffusivity parameters, in addition to their change in time. Sadly, in the wavenumber solution, the product of two spatial dependent variables (κ and  2 z T ) results in a cumbersome convolution which would nullify the advantage of its efficiency. However, the inhomogeneous solution is manageable by using a numerical method where both time and spatial steps are explicitly defined. To this end, we use a finite difference, the Crank-Nicolson diffusion scheme used previously for frictional heating problems (see Nielsen, Mosca, et al., 2010, and references therein).
Using the conditions of experiment s308, we compare the temperature and shear stress evolution between (1) a model including thermal dependence and (2) one with no thermal dependence, but with equal heat sinks in both. Finally, we repeat the simulation with thermal dependence but increase the amount of heat sink in parameter C s .
The result in terms of temperature evolution shows that under comparable frictional power curves, the temperature difference introduced by the thermal dependence is substantial (up to about 350°C at the temperature peak). If other parameters (frictional parameters, heat sinks) are kept equal, the difference between models with variable or constant thermal parameters differs largely toward the end of the frictional curve, in particular during the recovery phase (Figure 11b). The net decrease of conductivity, resulting in a higher background temperature, is sufficient to prevent the strength recovery in the flash weakening friction.
However, both models with or without κ, c thermal dependence can reproduce reasonably the frictional curve by adjusting the frictional or the heat sink parameters. For example, by increasing the heat sinks in NIELSEN ET AL.

10.1029/2020JB020652
17 of 31 the model with κ, c thermal dependence to C s = 2 10 6 MW m −2 , the temperature rise is buffered and is closer to that obtained in the model with no κ, c thermal dependence (where C s = 1.5 MW m −2 ). As a result, both models produce a similar frictional curve, compatible in both cases with the experimental observation (Figure 11a).
Finally, we show that a model with no heat sinks (C s = 0) and variable κ, c also fits the data, provided that the frictional parameters are altered to (T w , B) = (3800°C, 0.02 10 −6 ms −1 C −2 ) to compensate for the much higher temperature (Figure 11c). (Such a high value of T w largely above the weakening temperature of the rock-constitutive minerals, reinforcing the point that if heat sinks are ignored a realistic solution is not achieved.) These results were obtained by conducting simulations with 22 different combinations of (T w , B) and selecting the outcome with the best fit (also see discussion at the beginning of Section 6).
In conclusion, by adjusting the frictional parameters (T w , B) for each individual experiment, it seems possible to reach a reasonable fit for either variable or constant κ, c and either including heat sinks or not. However, a robust frictional model should reproduce a combination of different experiments with the same parameter set. This is better achieved by adding the presence of heat sinks and using a single combination of (T w , B, C s , T s ) for all experiments conducted on the same rock type. In addition, the effect of variable κ, c can be emulated in models with constant κ, c by lowering the power of heat sinks (parameter C s ).
Note that variations of mass density ρ will also occur due to thermal expansion, but are of second order in this analysis (e.g., less than 1% for granite under a temperature change of 1000°C, Richter and Simmons (1974)).

Flash Weakening Followed by the Formation of a Viscous Shear Layer
A model for frictional melting has been proposed for both steady-state (Nielsen et al., 2008) and transient (Nielsen, Mosca, et al., 2010) behavior. It accounts for the advancement of a melting boundary (solution of a Stefan problem) and the possibility that melt extrusion occurs through lateral injection veins (natural faults, Di Toro et al., 2005) or at the edges of the sample (experimental simulated faults: Niemeijer et al., 2009;Violay et al., 2014). We will only revisit some features of such a model here, and test to what extent a simpler flash weakening model differs from the melting dynamics.
In the case of frictional melt with extrusion, the thermal balance is quite different from that resulting in (2) from simple diffusive temperature. Apart from the heat loss due to phase transitions already discussed above, there is a radical change in the thermal diffusion equation with an additional convective term as a consequence of the advancement of the melting front into the rock. If we choose to attach the coordinate frame to the moving boundary, we can write: where ν is the current velocity at which the melt boundary is advancing into the solid. (We recall here that while the melt front is advancing, most of the excess melt produced is extruded, so that the shear heating remains confined in a thin melt layer.) Heat sinks due to phase transitions in this case are essentially due to melting latent heat L such that q s = ν ρ L. Thermal diffusion solutions with a moving boundary are known as the Stefan problem. It can be assumed that the boundary between melt and solid rock is at the melting temperature T m . As a consequence, ν can be computed by applying the boundary condition that T = T m and ∂ t T = 0 at the melting boundary (x = 0).
One key process in the presence of melting and extrusion is the advancement of the melting front which counteracts the advancement of thermal diffusion. In the absence of melting and extrusion, as in the case of flash weakening, the background temperature T is representative of recent frictional power dissipated NIELSEN ET AL.

10.1029/2020JB020652
18 of 31 Figure 10. Effect of thermal-dependent diffusivity and heat capacity κ, c. Temperature evolution was obtained with a similar frictional heat input, but for three different thermal models.
(3) Thermal dependence as in (1), but with the heat sink parameter increased to C s = 2.0 MW m −2 (dashed purple curve). At high temperature, the heat conductivity is lower, but the energy sinks are more effective, inducing a competing effect on the temperature rise. This allows to obtain similar results in (2) and (3)  on the fault, which induces heating within a finite thickness around the slip zone. Therefore temperature may be considered as a state variable storing the memory of the frictional heating history. However, the advancement of a melting front combined with extrusion will constantly reset the heat stored around the slip zone. The boundary will remain at the melting temperature T ≈ T m . Super-heating above T m may occur within the melt reaching a maximum at the center of the melt layer (Nielsen et al., 2008), but melt is rapidly extruded. Heat diffusion penetrates to an indicative depth   2 r z t within a given time interval t r . Given a shortening rate ν (melt front advancement) the typical time of residence of the heat in the rock adjacent to the slip zone will be t r = z/ν resulting in t r = 4 κ/ν 2 upon substitution of z. Therefore, the thermal memory of the system is reset over times of the order of t r , that is, a few seconds, assuming shortening rates of 1 mm/s and standard diffusivity values.
NIELSEN ET AL.

10.1029/2020JB020652
19 of 31  Figure 10, and the red dotted curve is temperature. In (c), we show results for a model with no heat sinks (C s = 0) and variable κ, c, but here, the substantially higher temperature requires a different combination of friction parameters (T w = 3,800, B = 0.02 10 −6 ) to fit the shear stress evolution reasonably well.
As discussed in Section 4.1 and in Figure 3, it appears that weakening precedes the bulk melting temperature. Therefore, we assume that the initial part of the weakening, before pervasive melt starts, is due to flash weakening behavior and can be modeled as such.
At later times, when pervasive melt and extrusion occur, a full steady-state condition may be reached, as predicted by the model of Nielsen et al. (2008) and experimentally confirmed by Violay et al. (2014). Such a steady state is not reached as rapidly in the absence of a migrating boundary. Only an apparent steady state is reached, later, in the presence of heat sinks q s (Equation 2), which inhibit the temperature rise is temporarily but only as long as the decomposition products are not depleted in the host rock. Therefore, the boundary migration needs to be included explicitly for an accurate temperature evolution; however, temperature diffusion in the presence of heat sinks does allow reproduce the appearance of a steady-state solution (the background temperature reaches a plateau, although the temperature diffusion continues to progress).
Finally, we note that the frictional recovery during the deceleration phase is expected to differ between flash weakening and frictional melt, although in both cases the sliding surface will undergo irreversible transformation at high slip velocity (phase transitions due to heat, roughness change through abrasion), and in both cases, the background temperature is higher at the end of sliding than at the beginning.
Ideally, in the occurrence of melting, a mixed model should be used with a transition from flash heating to frictional melting after the background temperature reaches T m . Such a transition can be rather complex, and several experiments show that there is a partial frictional recovery at that stage (Hirose & Shimamoto, 2005), although it is less pronounced in experiments performed at higher normal stress (Hung et al., 2019).
We do not implement here the mixed flash heating/frictional melting model which is a complex task and requires a separate study. However, we show to what extent frictional melt can be approximated either by a simplified Arrhenius dependence ( Figure 5 and discussion in Section 5.1), or by the flash weakening law of Section 5.2 (an example is shown in Figure S2). We find that with Arrhenius dependence, a reasonable fit of the stress curve is obtained only with models resulting in unrealistically low temperatures at the interface. On the other hand, flash weakening alone does not predict the recovery accurately in the presence of melting.
A mixed model can be developed also in the case of calcite, where an initial flash-weakening process is followed by the formation of a continuous viscous layer which undergoes high-velocity plastic shear, as posited in Violay et al. (2019) for bare surfaces in frictional contact. Where the slip initiates within a gouge, the eventual formation of a viscous layer after several mm of slip is also observed by Pozzi et al. (2019); however, the initial weakening includes a slip hardening phase. A microstructural evolution with complex processes is observed before the viscous shear is mature. This process is arguably difficult to represent with a flash weakening model.
Finally, a mixed model would also be indicated in the case where fluids are initially permeating the fault surface allow for elastohydrodynamic weakening: as slip accelerates, a transition between three lubrication regimes (boundary, mixed, and fully lubricated regimes) will occur as discussed in Cornelio et al. (2019).

Conclusions
We investigated different versions of thermal dependent, high-velocity friction models and tested how well they could be adjusted to replicate a number of experimental observations. Each friction law and each set of parameters was tested against several experiments, conducted under different normal stress and slip velocity, to verify whether the fit was robust under different conditions. We considered aspects of computational efficiency, the role of energy sinks, and the effect of thermal dependence on diffusivity and heat capacity.
The computation of temperature diffusion can be numerically costly, in particular foreseeing its use in models of an extended fault surface where T needs to be tracked at a great number of points. Therefore, we test and compare different methods: the simple discrete summation of the analytical solution, the discrete wavenumber transform, and a Crank-Nicolson finite difference scheme. We find that the two latter methods NIELSEN ET AL.
10.1029/2020JB020652 20 of 31 are comparable in terms of speed, for the problem to be solved in the case of frictional heating, and that both are much more efficient than the simple discrete summation of the analytical solution. In addition, the explicit spatial grid of the finite differences scheme allows to consider local variations of parameters. We took advantage of such a feature to investigate the thermal dependence of diffusivity and heat capacity.
We include additional heat sinks with an Arrhenius reaction rate, due to endothermal phase transitions are triggered under frictional heating. The temperature evolution is affected, in particular in the final phase of frictional recovery.
Diffusivity and heat capacity, often considered as constant, can undergo important changes when temperatures reach a substantial fraction of the melting temperature. We include such thermal dependence and show that it can have a substantial effect on the temperature and frictional evolution in problems of frictional heating. However, by altering the frictional parameters or the amount of the heat sinks, an effect similar to such thermal dependence can be simulated even in models which exclude it. Interestingly, J. R. Rice (2006) argued that the use of constant thermal properties as constant may be justified if such properties were considered as an equivalent average over a given thermal profile. However, here the thermal properties evolve with temperature, so that they will not only vary spatially but also change in time. Because the heat sinks will be activated only when temperatures are high, concomitantly with the expected drop in thermal conductivity, the time variation is best emulated by decreasing the effectiveness of sinks (C s ) as illustrated in Figures 10 and 11.
For models of flash weakening, we illustrate the existing trade-offs between frictional parameters, heat sinks, and the presence of thermal dependence of the parameters.
In addition to the flash weakening model, we investigate a simplified model including direct Arrhenius dependence of friction on background temperature. Such a model captures some of the main features of the weakening, however fails to account for the initial rapid weakening. A flash weakening model, instead, captures well the initial weakening; however, it requires to include the evolution if the background temperature to reasonably reproduce friction from start to end.
We discuss the differences between flash weakening and profuse frictional melting. It is known that both can take place within a single slip episode, with the flash heating occurring at the start. Given the differences between flash weakening and frictional melting processes, an accurate representation should include both, and then model the subtle transition from one to the other, an endeavor that we leave for future work. However, as an approximation, we show here that a flash weakening model including thermal dependence with heat sources and sinks is able to reproduce cases of frictional melting reasonably well over the limited interval of parameters of the experiments shown here in support.
Finally, we note the experiments presented here are limited to precut, cohesive rocks of two end-members (Carrara marble, gabbro) under dry conditions. One needs to consider that rupture on natural faults will develop in different lithologies, including clays. In addition, deformation of natural faults develops in more complex ways, including diffuse strain, gradual localization of slip within an initially thick fault gouge, multiple branching of rupture, and other dissipative processes which take place off-fault. Fault deformation, therefore, develops in part within a volume rather than only by slip across a surface or by shear of a thin gouge layer. The off-fault deformation and dissipation contribute to the energy balance and to the stress-strain relation, an aspect which needs to be integrated somehow with the frictional slip on the fault for realistic modeling, in ways which are not yet defined. Until then, such limitations need to be kept in mind when using laboratory-derived friction laws in earthquake models.

A3. In a Nutshell: The Temperature Update Scheme
Finally, we may sum-up the iterative scheme as follows. Upon discretization in time steps of size δt and M wavenumber steps of size δs, the memory variables θ m are updated at each time iteration n, and summed in order to obtain the current temperature T according to: where   m is the m th memory variable from the previous time step and s m is defined in Equation A10. The above temperature T can then be used to update the temperature-dependent stress according to Equation 13. If the gradient of the temperature is required it can simply be obtained by For most practical purposes, M = 16 yields a satisfying approximation (with a precision of a few percent after t = 20s). M = 32 yields a result which differs less than 1/10,000 at t = 20 s when compared to the classic solution obtained by direct discretization (2). Timing the computation example above with M = 16 (on a desktop computer with python) yields 8.44 ms for the wavenumber solution versus 132 ms for the classic time summation, and storage of 16 floating point values versus 500, that is, an estimated gain of 97% in CPU time and 94% in memory storage.

A4. Finite Differences
The third type of temperature diffusion was tested, a finite differences Crank-Nicolson algorithm. The latter is the most flexible type of computation as it allows to consider parameters which are inhomogeneous along z and also variable in time (see Section 6 on thermal-dependent diffusivity and heat capacity).
Additional tests were performed using Fortran for the three types of the computation above, and using sufficient grid points in the finite difference scheme to avoid reflections from the end boundary of the grid. A sufficient precision was used for all 3 methods to diverge by no more than 2%. The resulting timing was: (1) wave number summation CPU time: 8.66E-2 s. Classic time summation CPU time: 0.30 s. Finite differences Crank-Nicolson CPU time: 1.87E−2 s.
Time summation is by far the less effective method in terms of time and memory usage, although it is easy to implement. The performance in time is of the same order for the wavenumber summation and the finite differences (in this example, the finite difference code is even more performant than the wavenumber summation).
Therefore, it is worth conducting the temperature solution using finite differences, for both best performance and increased flexibility. For the finite difference to work correctly, a sufficiently small spatial step should be selected (in our example, we use 2E−4 m) and in time (in our example, we use 10E−2 s) to NIELSEN ET AL.

10.1029/2020JB020652
25 of 31 reflect the fastest changes expected in the heat flow regime. Once the duration of the simulation and the spatial step are determined, a sufficient number of nodes should be introduced to avoid reflections from the end boundary of the model (we use 50 nodes in this examples). To design the optimal parametrization, tests can be performed under the typical conditions required, decreasing the step size until reasonable convergence is achieved and increasing the number of steps until spurious boundary effects are negligible.
In case that the thermal computation is coupled with a rupture model, the time steps in both methods need to be harmonized.
The above tests were performed on 16 CPUs Intel(R) Xeon(R) E5-2640 v3 clocked at 2.6 GHz. Example codes for temperature computation are provided in Supplementary Materials.