Fault Friction During Simulated Seismic Slip Pulses

Theoretical studies predict that during earthquake rupture faults slide at non‐constant slip velocity, however it is not clear which source time functions are compatible with the high velocity rheology of earthquake faults. Here we present results from high velocity friction experiments with non‐constant velocity history, employing a well‐known seismic source solution compatible with earthquake source kinematics. The evolution of friction in experiments shows a strong dependence on the applied slip history, and parameters relevant to the energetics of faulting scale with the impulsiveness of the applied slip function. When comparing constitutive models of strength against our experimental results we demonstrate that the evolution of fault strength is directly controlled by the temperature evolution on and off the fault. Flash heating predicts weakening behavior at short timescales, but at larger timescales strength is better predicted by a viscous creep rheology. We use a steady‐state slip pulse to test the compatibility of our strength measurements at imposed slip rate history with the stress predicted from elastodynamic equilibrium. Whilst some compatibility is observed, the strength evolution indicates that slip acceleration and deceleration might be more rapid than that imposed in our experiments.

been made to verify this compatibility by comparing high velocity fault rheology obtained in laboratory experiments with earthquake source parameters estimated from strong motion seismograms, with limited success (Chang et al., 2012;Fukuyama & Mizoguchi, 2010).
In nature, there is a feedback between the slip rate on the fault (associated with a given stress state through elastodynamic equilibrium) and the shear strength of the fault zone materials. A substantial body of work has shown that fault strength is dramatically lower at rapid deformation rates in comparison to that during slow interseismic deformation (Brantut et al., 2008;Di Toro et al., 2011;Faulkner et al., 2011;Goldsby & Tullis, 2011;Han et al., 2007;Hirose & Shimamoto, 2005). The low strength observed during rapid deformation in the laboratory can explain the propagation of ruptures in numerical models (Noda et al., 2009;Noda & Lapusta, 2013), and is quantitatively consistent with the low heat flow and shear heating estimated from borehole measurements after the 2011 Tohoku-oki earthquake (Fulton et al., 2013;Ujjie et al., 2013). High velocity friction experiments have shown that sliding velocity exerts a first order control in governing the strength of faults (Di Toro et al., 2004Goldsby & Tullis, 2011;Han et al., 2007;Tsutsumi & Shimamoto, 1997). Several physical models have been proposed to explain this weakening. Nielsen et al. (2008Nielsen et al. ( , 2010 developed a model of frictional melting in basic igneous rocks which showed good agreement with the experimental data of Hirose and Shimamoto (2005). The flash heating model (FH), based on frictional heating at microscale contacting asperities (Beeler et al., 2008;Rice, 2006), was successful in explaining experimental observations in crystalline felsic rocks in the absence of bulk frictional melting (Goldsby & Tullis, 2011). More recently in carbonate rocks, grain size sensitive creep has been proposed to explain the weakening behavior due to the presence of nanometric grains which facilitate rapid diffusive mass transfer (Demurtas et al., 2019;De Paola et al., 2015;Pozzi et al., 2018).
Although velocity exerts a direct control on fault strength, a number of studies have also demonstrated the important role of fault slip history in determining strength evolution; In particular, hysteresis in the frictional strength has been systematically observed between acceleration and deceleration phases of experiments, suggesting a change in the physical state of the fault during experiments (e.g., Goldsby & Tullis, 2011;Proctor et al., 2014;Sone & Shimamoto, 2009). A number of authors have shown that this difference can be accounted for by considering the temperature rise of the fault as a state variable in FH models (Proctor et al., 2014). By controlling the thermal evolution of gouge material during high velocity friction tests Yao et al. (2016) demonstrated the key role of temperature in high velocity weakening mechanisms. In parallel, experimental data and modeling by Passelègue et al. (2014) showed that subtle effects of temperature variations on fault rock properties, such as thermal conductivity, could lead to significant effects on fault weakening.
Given that temperature evolution is directly coupled to slip history and strength through dissipation of mechanical energy, the strength evolution and the resulting dynamics of earthquakes are expected to be controlled by thermo-mechanical feedbacks.
Here, we aim to investigate the role of slip history on the the frictional response of rocks, and test the compatibility of laboratory-derived strength evolution with elastodynamics. We use a slip history that is representative of earthquake source-time functions in the form of a so-called "modified Yoffe function" as derived by Tinti, Fukuyama, et al. (2005) (see Methods section), and test a range of initial accelerations for a fixed total slip (Figure 1). We first explore the physical mechanisms that give rise to the observed complex frictional response, and show that thermally activated mechanisms like flash heating (near the onset of the slip) and viscous creep (at late stages) are broadly consistent with the observed frictional response. This agreement confirms the key role of temperature and temperature history in the high velocity frictional behavior of rocks. We then analyze the compatibility of our experimental results with the traction evolution expected from a simple elastodynamic slip pulse model. The measured frictional response is not totally consistent with the model in that it shows more abrupt weakening at the onset of slip and too large re-strengthening at the termination of slip. These differences indicate that elastodynamics would likely produce either shorter slip pulses or sharper drops in slip rate at the tail end of pulses (self-healing).
HARBORD ET AL.

Methodology
Our experiments are analogous to slip on a single point on a fault, and we therefore need to select an appropriate slip function representative of a rupture propagating through a single point in space. In practise it is not possible to define a unique solution since fault slip history depends on the interactions between fault strength and elastodynamics, however there are several candidates we may choose to represent the seismic source (Kostrov, 1964;Nielsen & Madariaga, 2003;Tinti, Fukuyama, et al., 2005;Yoffe, 1951). The Yoffe function represents an attractive solution to model slip history due to its direct compatibility with elastodynamic rupture propagation (Nielsen & Madariaga, 2003;Tinti, Fukuyama, et al., 2005), and the fact that it shares similarity to slip histories observed in experimental studies (Berman et al., 2020;Rubino et al., 2017). It is characterized by a singular acceleration at the moving crack tip, corresponding to the crack tip stress concentration, followed by an inverse square-root decay in velocity with respect to time (Figure 1b). This results in slip that is approximately square root in time at a fixed observation locality (Figure 1a). Given that singular acceleration is unrealistic in nature, and also not possible to simulate in the laboratory, we used a regularized form of the Yoffe function presented in Tinti, Fukuyama, et al. (2005). The solution is equivalent to convolving a true Yoffe function (with singular acceleration) with a triangular function of time duration, 2t s , where s t is defined as the smoothing time. Small values of the smoothing time generate more impulsive, shorter duration events with faster initial accelerations, and increasing s t generates longer duration, less impulsive events (Figure 1). In experiments, the deconvolved timespan R t = 2 s, and maximum displacement, U tot = 1.65 m, were kept constant, to simulate seismic slip equivalent to a typical M w = 7 earthquake (Wells & Coppersmith, 1994). We varied s t from 0.05 to 0.8 s, with the rise time,   r R s t t 2t , varying between 2.1-3.6 s, which may be considered analogous to varying the rupture velocity (Cochard & Madariaga, 1994).
We utilized a slow to high velocity rotary shear apparatus (SHIVA, Di Toro et al., 2010) installed in the HPHT laboratory at INGV in Rome. The apparatus is capable of applying up to 70 kN of axial load using an electromechanical piston (Bosch-Rexroth EMC105HD), which is servo controlled at a frequency of 8 kHz. A 300 kW motor servo-controlled at 16 kHz drives the rotary motion up to 3000 RPM, with which we achieved an instantaneous acceleration of <80 m/s 2 and a deceleration of <60 m/s 2 ( Figure 1c); outside of this range machine vibrations were too strong to gather reliable data. Displacement was measured using a high resolution encoder (6,297,600 divs) for low velocity (<0.15 m/s) and a low resolution encoder (4,000 divs) for high velocity (≥0.15 m/s), the encoder-derived velocity (Figure 2 gray curves) and the imposed velocity function (Figure 2 black curves) show good agreement. Annular cohesive samples of Etna basalt and Carrara marble of 50 mm external and 30 mm internal diameter were prepared for testing, and were squared using a lathe before being ground with #80 grit prior to experimentation. All tests were performed at a normal stress of 10 MPa. Torque was measured using an S-type load cell on the stationary side of the apparatus and all data was logged at 12.5 kHz. A total of over 60 simulated slip events are presented in this study (see Text S2). In the majority of experiments slip pulses were repeated using the same sample, with the normal load kept constant during a minimum time period of at least 20 min between individual pulses. Measured frictional HARBORD ET AL.   Tinti, Fukuyama, et al. (2005). Here s t is the smoothing time (see main text).
strength was found to be highly reproducible after the first pulse (see Figure S1), which we interpret to indicate a consistent microstructure between individual pulses.

Experimental Results
All experiments show three stages of behavior typical of high velocity friction tests: (a) An elastic loading and slip-strengthening phase, which corresponds to an increase in friction coefficient from an initial value at zero velocity,  0 = 0.5-0.6 (τ 0 = 5-6 MPa), to a peak as slip rate increases,  p = 0.6-0.8 (  p= 6-8 MPa, Figures 2b and 2c); (b) A breakdown phase past the peak in frictional strength, where friction drops from  p to a minimum weakened value,  min , which is generally coeval with the peak in velocity. Values of  min are typical of high velocity friction, with values around  min = 0.05-0.2 (  (Figures 2d and 2f). At the largest values of t s = 0.8 s, during the restrengthening phase, the frictional strength in basalt reaches an apparent value of μ = 1-1.2 (10-12 MPa), before reducing to μ = 0.6-0.9 (6-9 MPa). In tests where this behavior was observed the sample often failed in a brittle manner with audible cracking coeval with the peak in friction.
When comparing between experiments we observe a clear dependence between the overall frictional behavior and the imposed smoothing time. Inspection of tests with t s = 0.1 s (with an initial acceleration rate A V /t  is of similar value to t s . In order to quantify how the smoothing time t s influences the overall mechanical behavior of the simulated faults, we now estimate key quantities relevant to the energetics of faulting.

Weakening Distance
The coseismic weakening distance is an important parameter governing the propagation of earthquake rupture, providing an indication of rupture efficiency (Ida, 1972). It may also provide insight into the weakening mechanisms active during experiments (Hirose & Shimamoto, 2005;Niemeijer et al., 2011). To estimate the weakening distance, w D , we consider the distance at which strength decreases by 95% and use the same formulation as Hirose and Shimamoto (2005) and fit data using a least squares regression. Where the reduction in strength is not monotonic (e.g., Figure 6f), we excluded results from the analysis. The fitted values range from 0.05-1 m for the presented experiments (see Table S1), and are strongly dependent on the slip function applied, but are within the range of values presented in previous studies under similar experimental conditions (Chang et al., 2012;Niemeijer et al., 2011). A clear trend is observed between t s and w D (Figure 3), given that smaller values of s t correspond to a larger acceleration demonstrating an inverse dependence on the initial acceleration. For example, for Carrara marble (Figure 3

Energy Dissipation
Following previous literature, we define the breakdown work b W [MJ m −2 ], according to the general definition of (Tinti, Fukuyama, et al., 2005), where min D is the displacement when τ = τ min . We also define restrengthening work, r W , in a similar manner, accordingly: Both of these parameters were calculated by numerical integration of the experimental shear stress record with respect to slip (see Nielsen et al., 2016). This provides a quantitative estimate of energy partitioning during experiments. We find that both b W and r W depend strongly on the impulsiveness of the Yoffe function applied.

Driving Processes of Frictional Evolution in the Presence of Complex Slip Velocity Histories
Our results show that strong variations in slip rate induce correspondingly strong variations in frictional strength, with a rapid weakening at high slip rate and significant restrengthening as slip rate decreases. In order to identify the key driving mechanisms responsible for these variations, here we test whether such variations are captured and predicted by existing physics-based high velocity friction laws.
One key experimental observation is that the minimum strength is almost systematically occurring at the peak velocity achieved during the tests (Figure 2), which corresponds to a velocity-weakening behavior of the rocks. Such a behavior is typically associated with some state evolution, whereby instantaneous changes in slip rate should induce strength increase, followed by adjustments towards a lower strength state. Here, the observation of direct correlation between peak slip rate and minimum strength indicates that this state evolution occurs over time (and slip) scales much smaller than that of the change in slip rate imposed in the experiments. This is consistent with models where state was assumed to evolve over slip distances of the order of tens of microns, indeed much shorter than the slip scales measured here (e.g., Noda et al., 2009). Therefore, we primarily focus on a description of strength that excludes considerations of short-slip state evolution.
Firstly, we focus on predicting the strength of experiments using Carrara marble. We explore two commonly proposed descriptions of strength, flash heating (Goldsby & Tullis, 2011;Proctor et al., 2014;Yao et al., 2016)   Etna basalt hosted faults, and provide a simple comparison to previous models of high velocity friction for melt accommodated weakening.

Flash Heating
Weakening by flash heating (FH) is based on the idea that contacting asperities at the sliding interface dramatically weaken at some threshold temperature (Beeler et al., 2008;Rice, 2006). High velocity experimental data obtained using simple slip rate histories has been shown to be in general agreement with this model (Goldsby & Tullis, 2011;Proctor et al., 2014;Yao et al., 2016). The shear strength is assumed to be given by f is the low velocity coefficient of friction, w f the weakened coefficient of friction and w V a critical weakening velocity that depends on the difference between the ambient fault temperature T and a critical weakening temperature w T . The critical velocity defines a threshold at which a contacting asperity spends a portion of its lifetime above the prescribed temperature w T corresponding to some weakening process, for example, mineral decomposition (see Text S2 for further details). Here we also explore the impact of temperature dependent asperity strength and size, which vary according to an asperity stress exponent, n.
Using the velocity histories imposed in the experiments, we first modeled fault strength with a fixed ambient temperature ( Figure 5, blue curves), that is, purely velocity-dependent strength. Comparison of this model directly with our experimental data shows that for all cases, purely velocity dependant strength is initially consistent with weakening behavior but diverges with increasing time. When fully accounting for the rise in background temperature, modeled by introducing the bulk heat dissipation and diffusion in the rock (e.g., Proctor et al., 2014; see Text S2), the model predictions significantly improve ( Figure 5 red curves), and strength predictions generally match initial weakening behavior during acceleration to peak velocity (t < 2s). However, the models still tend to diverge from the data at larger timescales, overestimating HARBORD ET AL.  To place the modeling results in the context of previous results we also compiled data from experiments performed at similar conditions in SHIVA, where "box-car" slip histories that is, constant acceleration to constant velocity followed by a constant deceleration to zero velocity, were employed (Violay et al., 2013(Violay et al., , 2019. Experiments run at a range of velocities are reproduced and compared to models of flash heating (Figures 5d-5f), and are shown to highlight that all model predictions tend to overlap at high constant velocity. This overlap tends to obfuscate the determination of realistic model parameters, particularly at the highest velocities, at least for the range considered here. At relatively low velocity conditions differences are observed between individual models (V = 0.3 m/s, Figure 5f). In agreement with experiments where slip rate was given by Yoffe functions, we observe that a reasonable prediction of data can be obtained when n = 1.
A consistent observation in all flash heating models is that they significantly over-predict the restrengthening behavior, and demonstrate that addition of temperature dependant asperity properties does not significantly improve the prediction of strength. An improvement may be yielded by accounting for a log normal distribution of asperity sizes which is a smooth function of velocity (Beeler et al., 2008), which may be more representative of the geometry of experimental surfaces (Candela & Brodsky, 2016).  Here we check consistency of our data with direct predictions from a GSS creep law derived from deformation of fine grain calcite aggregates at high pressure temperature conditions (Schmid et al., 1977):

Grain Size Sensitive Creep
where A is a pre-exponential factor,  is the shear rate, d is the grain size, m is the grain size exponent, E a is the activation energy, R is the gas constant and k the stress exponent (see Supporting Information S1, for detailed parameter values and modeling assumptions). Similarly to our computations using flash heating, we fully account for the background temperature evolution in the rock with temperature dependant thermal diffusivity.
Results from GSS models systematically overpredict the strength of faults at short timescales, and do not predict the initial weakening for all values of s t (Figures 6a-6c). However, from the late stages of weakening, up to the later stages of restrengthening we observe a good prediction of strength evolution. When s t = 0.1 s the restrengthening is well matched, however for larger values of s t = 0.4 and 0.8 s, GSS models systematically predict a faster restrengthening rate during the final deceleration period than experiments. At cessation of slip, as velocity decreases below ∼1 mm/s, the model predicts a complete loss of strength at all conditions consistently with the GSS flow law (Figures 6a-6c). The prediction of zero strength when compared to the experimental data suggests that GSS may no longer accommodate deformation during the final stages of slip. A grain size of 1 nm is systematically required to predict the strength of faults with Yoffe velocity history (Figures 6a-6c). Again, for comparison purposes, results are shown from Violay et al. (2013 and. With GSS models a reduced degree of convergence is observed at constant velocity, although the same grain sizes show a generally similar behavior at large timescales where sliding velocity is high (V>=1 m/s, Figures 6d and 6e). A reasonable agreement between models and data is observed for a box-car slip history at velocity of 0.3 m/s (Figure 6f), identifying the wider applicability of the creep model across the range of sliding velocities. Again, 1 nm grain sizes are required to accurately predict the strength of experiments conducted with a box car velocity history at low velocity, low normal stress conditions (Figure 6f), however use of a 100 nm grain size shows good agreement with the higher velocity experiments (Figures 6d and 6e).
The nanometric grain size required to match fault strength is probably unrealistic (De Paola et al., 2015;Pozzi et al., 2018Pozzi et al., , 2021Violay et al., 2013). However, this could be remedied by using a modified, much larger value for the preexponential factor in Equation 4; here, we decided to use an empirical estimate from an exisiting data set obtained at low strain rate, but several physical phenomena might dramatically change that value. Pre-exponential factors include contributions from grain boundary geometry (thickness and roughness) and grain boundary self-diffusion (Poirier, 1985). It is possible that the fault microstructure during initial weakening, which has been demonstrated to result from dislocation avalanches (Spagnuolo et al., 2015), may result in anhedral nanograins with larger grain boundary aspect ratio when compared to the final microstructure which is likely to have annealed during cooling of the fault. Raj and Ashby (1971) demonstrated that increases in the aspect ratio of contacting grain boundaries increases the self-diffusion coefficient, resulting in reduced yield stress, which may preclude the need for unrealistically small grain sizes.
It is also important to consider that if flash processes occur during initial fault weakening, temperature may be locally higher than predicted from GSS models. An initial flux of heat resulting from asperity-scale processes may be sustained throughout the test duration (Aretusini et al., 2021), which would allow larger grain sizes to give quantitative agreement with the experimental data. We also note that in Violay et al. (2019) the authors were able to match the final fault restrengthening of the data presented in Figure 6d by accounting for heat loss in two dimensions.

The Importance of Accurate Temperature History
Both the flash heating and grain size sensitive creep models demonstrate the importance of incorporating temperature history into the models, and shows that it is important to consider appropriate thermal HARBORD ET AL.

10.1029/2021JB022149
9 of 17 properties in model boundary conditions. This point was first highlighted by Yao et al. (2016) where experiments were conducted using sample holders of varying thermal diffusivity, demonstrating that varying temperature histories give differing strength evolutions. Here we further test this hypothesis by comparing the output of both FH and GSS models by using previously published calcite gouge experimental data obtained with a range of sample holders of varying diffusivity (Cverna, 2002). We consider, in order of increasing thermal diffusivity, grade 4 Titanium alloy (Ti90Al6V4, Pozzi et al., 2018), AISI 316 stainless steel (De Paola et al., 2015) and tungsten carbide (Smith et al., 2013). We approximate each experimental geometry as closely as possible in 1D, with the principal slip zone localized asymmetrically on the boundary between the gouge layer and the sample holder with appropriate thermal diffusivity (De Paola et al., 2015;Pozzi et al., 2018;Smith et al., 2013, see Supporting Information S1).
For both rheological models we observe that for increasing thermal diffusivity, the fault strength predictions increase, consistently with experimental observations (Figure 7). Generally FH does not predict the initial weakening behavior, predicting faster weakening in Titanium alloy ( Figure 7a) and Stainless steel (Figure 7b), and less abrupt weakening for tungsten carbide (Figure 7c). The differences in model predictions and experimental data may result from strain localisation and grain crushing that occurs during early stages of slip in gouges (Logan et al., 1992). During steady state sliding conditions the FH models are able to predict strength evolution with reasonable success, and in particular predict a slow progressive weakening with slip resulting from a progressive temperature rise, sharing similarities to the experimental data (Figures 7a-7c). The restrengthening is systematically over predicted by the flash heating models, similarly to the results shown in the previous section.
Initial weakening is predicted comparatively better for GSS than it is for FH. Similarly to FH, the GSS models also predict progressive weakening observed during constant velocity conditions (e.g., Figure 7e). The restrengthening rate predicted by models is slower than observed in experiments, and consistently with previous discussion of GSS models, strength falls to zero as slip arrests. Generally we observe that the best predictions of fault strength for the gouge experiments are obtained for the previously used parameter HARBORD ET AL.

10.1029/2021JB022149
10 of 17 Figure 7. Effects of varying thermal diffusivity in full thickness models (blue curves) with realistic sample boundary conditions compared to previously published experimental data (gray curves). Panels (a-c) are modeled using the flash heating model described in previous discussion, with fixed thermal diffusivity with curves labeled according to the asperity stress exponent used. Panels (d-f) show the same experimental data, however this time using the grain size sensitive creep model defined in the previous discussion with fixed thermal diffusivity, curves are labeled according to the grain size used in the model prediction. Thermal diffusivity increases from the left to right of the figure.
set, except at the highest normal stress conditions. Given that we may use the same parameters in the constitutive friction law (FH or GSS), it suggests that the key variable is the bulk temperature evolution. In short, reconciling these individual experimental observations is difficult without carefully considering model boundary conditions and demonstrates that it is of fundamental importance to accurately capture on and off-fault thermal boundary conditions accurately, confirming the conclusions of Yao et al. (2016).

Interplay of Weakening Mechanisms
Flash heating predictions are better at shorter timescales, whereas longer timescales demonstrate a better prediction by GSS models. In fact all FH models significantly over predict the final strengthening behavior. Taken together these observations suggest that multiple weakening mechanisms may operate simultaneously during experiments with a potential transition in dominance. At early stages when the bulk fault temperature is low, and GSS is not efficient, behavior will be dominated by asperity scale flash heating processes leading to bulk heating of the principle slip zone. However as slip and fault temperature increases, GSS deformation becomes increasingly favorable. This transition has been previously proposed by Pozzi et al. (2018) and De Paola et al. (2015), however they did not explicitly consider FH at early stages of slip, suggesting instead that the transition is simply from cataclastic processes to GSS. If FH was active during early stages of slip it is possible that the high contact temperatures during weakening may be sustained during later stages of the experiment and deformation could be accommodated by larger grain sizes. Effectively the two constitutive equations define a threshold temperature at which fault strength approaches a residual strength. In the case of FH this is given by the temperature at which a generic weakening process occurs (which could be GSS), whereas for GSS it defines the temperature at which efficient diffusive mass transfer occurs, in both the governing state variable is fault temperature.

Weakening and Restrengthening in Basalt
Weakening of basaltic experimental faults is facilitated by frictional melting, which leads to the formation of a hot low viscosity melt layer (Hirose & Shimamoto, 2005;Niemeijer et al., 2011;Rempel & Weaver, 2008;Violay et al., 2019, see Movies S3 and S4). Modeling of weakening accommodated by melting has previously addressed by Nielsen et al. (2008) and Rempel and Weaver (2008) who explicitly considered the effects of the effects of progressive melt formation during high velocity sliding. Melting of rock during frictional sliding at high velocity can be shown to result in a complex 2-stage weakening behavior, reflecting the degree of melt layer formation, with the presence of initially patchy melt leading to strengthening, followed by secondary weakening due to the formation of a continuous meltlayer (Del Gaudio et al., 2009;Hirose & Shimamoto, 2005;Rempel & Weaver, 2008). This is evident in our experiments with slower initial acceleration rates (e.g., s t = 0.8 s, Figure 2f). When acceleration is sufficiently high, then weakening is monotonic (Figure 2d), consistently with Del Gaudio et al. (2009).
Turning attention now to the restrengthening phase of basalt experiments, we observe a clear relationship between the final deceleration rate and restrengthening behavior (Figure 8). Where final deceleration is sufficiently rapid, t s = 0.1 s, then no strength overshoot is observed, and friction monotonically increases up to the end of the experiment, with  1 = 0.9. As the deceleration rate is decreased as a result of increasing t s , we observe increasing amounts of strength overshoot, and faster restrengthening rates. For the largest value of smoothing time ( s t = 0.8 s), the strength overshoot is considerable, with a coefficient of friction close to 1.4 (Figure 8a), almost twice the initial value of μ = 0.7. Such large increases in strength suggest melt solidification and cohesion of the fault, and where large overshoot was observed cracking was heard, identifying that the melt solidified and failed in a brittle manner (see Movie S4). In the limit of adiabatic instantaneous deceleration, the fault stress would instantaneously drop as a result of the Arrhenius dependance of melt viscosity (Giordano et al., 2008). However where deceleration is slow, heat diffusion dominates and significant strengthening occurs due to melt solidification. According to Nielsen et al. (2008)

Are Laboratory Friction Data Compatible With Elastodynamics?
In the previous section we analyzed the potential driving processes that produce the observed evolution of friction in response to an imposed slip history. In nature, during an earthquake, the evolution of frictional strength feeds back into the slip history due to elastodynamic stress redistribution and the requirement of stress equilibrium. To illustrate this, let us consider the elastic stress field associated with anti-plane slip along a 1d linear fault trace: where G is the shear modulus, s c is the shear wave speed, V the on fault particle velocity, x is the position along the fault, and K the dynamic load associated to points on the fault that are still slipping.
The dynamic load term in Equation 5 is composed of the integrated slip history across the entire span of the rupturing fault. Waves radiated from other points on the fault results in dynamic loading which modifies slipstress history, typically resulting in a heterogeneous slip history on the fault plane. In addition, the transfer of stress, wave propagation and rupture velocity depends also on the geometry of the rupturing fault, so that the typical non-planar geometry of natural faults will also influence the slip-stress history (e.g., Romanet et al., 2020). Therefore, slip history at a point on a fault is highly non-unique and depends on the entire integrated rupture history, and in practise there is no unique test of elastodynamic compatibility.
In order to test the compatibility of our experiments with elastodynamics, we must make several simplifying assumptions. To do this we consider a steady-state slip pulse model, where both the rupture velocity and source duration are constant. In this case the elastodynamic equilibrium in anti-plane geometry can be simplified to where  b is the ambient shear stress, L the pulse length (equivalent to the product of rupture velocity and total rise time),  G is the shear modulus multiplied by a function of the ratio of rupture speed r V and shear wave speed. Equation 6 gives the dynamic elastic stress produced by the slip rate distribution along the pulse. In our experiments, the slip rate is imposed as a function of time. Here we consider that this slip rate evolution represents the relative motion of two opposing points along a steadily propagating pulse. Choosing a constant rupture speed, we first compute the elastic stress by direct integration of Equation 6, and compare it to the measured experimental strength (for details of methodology see Viesca and Garagash (2018) and Text S3). Since strength should be equal to stress during slip to satisfy mechanical equilibrium, any deviations between predicted stress and measured strength would indicate inconsistency between the rheological behavior of the fault and our choice of imposed slip rate.
HARBORD ET AL. The stress predicted by imposing the velocity history is only qualitatively compatible with the overall evolution of strength during tests: there is an initial weakening phase, with strength decreasing until sliding occurs at constant stress after which the stress increases during final slip deceleration, although the precise magnitudes and timings do not agree. In particular, the predicted final stress increase occurs later than in experimental observations, with a comparatively smaller magnitude.
We can also use our strength measurements to predict what would be the slip rate evolution along a hypothetical pulse, that is, to determine V(x) based on    x in Equation 6, assuming this time that strength is equal to elastic stress, and compare this slip rate to the originally imposed experimental slip rate. By imposing zero slip velocity before and after the rupture interval, we also constrain the background stress  b for our hypothetical pulse (see Text S3 for details). While there are encouraging similarities between model and observation, the predicted slip rate is quantitatively inconsistent with the imposed one. In particular, the final increase in stress measured in experiments results in back slip where velocity is negative (e.g., Figure 9k). The prediction of back slip is not realistic and would not occur during spontaneous rupture.
Overall, the experimental data show limited compatibility with our simple slip pulse model. Considering that the strength is mostly controlled by slip rate (with short state evolution distances) and temperature, we expect that slip rate and strength evolution that are compatible with elastodynamics would involve abrupt changes in slip rate together with rapid strength changes, both at the rupture tip and at the cessation of slip. For instance, in Carrara marble (e.g., Figures 9a and 9d, sample S1766 g), imposing a relatively constant slip rate after initial acceleration will lead to slowly decreasing strength (due to temperature rise), which is likely to eliminate the possibility of back slip. Then, an abrupt velocity drop might be consistent with an increase in strength above the elastic stress, producing spontaneous slip arrest. Our observations of peak weakening coeval with peak velocity is partially at odds with elastodynamic models (Tinti, Fukuyama, et al., 2005;Tinti, Spudich, et al., 2005), where peak weakening occurs after peak velocity during slip rate deceleration. In contrast, slip functions and associated elastic stress in Mikumo et al. (2003) show peak weakening coeval with peak slip velocity, after which slip rate drops to a relatively constant value, which is qualitatively consistent with our previous discussion.
The results on Etna basalt further support the requirement for rapid final slip deceleration as the strength increases quickly during melt solidification, resulting in a highly unrealistic minimum compatible slip rate of ≈ −2 m/s (Figure 9l), consistent with the notion of melt "fusion" during high velocity sliding (Fialko & Khazan, 2005).

Conclusions
In this work we document results from high velocity friction experiments imposing a realistic source time history, in order to investigate how fault strength evolves during earthquakes. Simple first order observations show that the weakening distance and breakdown work are inversely dependent on the initial acceleration rate. Experimental results combined with modeling demonstrate that the high velocity strength of faults during variable velocity strongly depends on prior sliding history and temperature evolution. Carbonate built fault strength can be accurately predicted by flash heating at small time scales and grain size sensitive creep at larger timescales, provided that model boundary conditions are faithful to experimental conditions. Where flash heating is utilized to model the fault strength of carbonate built faults, then final restrengthening is always over predicted. In the case that a creep constitutive law is used there are some significant differences between requisite grain sizes for accurate strength predictions and observed grain sizes from microstructural observations (De Paola et al., 2015;Pozzi et al., 2018Pozzi et al., , 2021. This discrepancy remains unresolved, and might be due to incorrect assumptions about our choice of deformation mechanism or the estimated temperature. However, the remarkable agreement between model predictions and observations indicates that thermally activated viscous flow laws are good candidates for the rheology of faults at high velocity.
These results provide an important validation of constitutive laws of frictional strength under non-constant velocity histories, justifying their use in coupled elastodynamic models, when the temperature rise of the fault is considered (e.g., Brantut  In our experiments, we imposed a slip rate history and measured the resulting strength. In nature, there is a feedback between strength and slip rate evolution due to elastodynamic stress redistribution. We tested the consistency of our experimental data with a simple elastodynamic model, and found discrepancies, that is, the measured strength does not match the predicted elastic stress associated with the imposed slip.
HARBORD ET AL.

10.1029/2021JB022149
14 of 17 Figure 9. Experimental data compared to elastodynamic solution using steady state pulse model of rupture propagation. Dashed shear stress curves indicate solutions compatible with the imposed velocity history during an experiment, and overlay smoothed experimental data (solid curves, labeled as observed in (a)). The solid velocity history is that which is predicted from the measured evolution of shear stress (traction) during an experiment, lines are colored gray where V < 0 m/s. It is likely that the rheology of the fault gives rise to velocity changes (acceleration and deceleration) more abrupt than in our imposed source-time functions.

Data Availability Statement
All raw experimental data are available at https://zenodo.org/record/4644245.