An empirically based steady state friction law and implications for fault stability.

Empirically based rate-and-state friction laws (RSFLs) have been proposed to model the dependence of friction forces with slip and time. The relevance of the RSFL for earthquake mechanics is that few constitutive parameters define critical conditions for fault stability (i.e., critical stiffness and frictional fault behavior). However, the RSFLs were determined from experiments conducted at subseismic slip rates (V < 1 cm/s), and their extrapolation to earthquake deformation conditions (V > 0.1 m/s) remains questionable on the basis of the experimental evidence of (1) large dynamic weakening and (2) activation of particular fault lubrication processes at seismic slip rates. Here we propose a modified RSFL (MFL) based on the review of a large published and unpublished data set of rock friction experiments performed with different testing machines. The MFL, valid at steady state conditions from subseismic to seismic slip rates (0.1 µm/s < V < 3 m/s), describes the initiation of a substantial velocity weakening in the 1-20 cm/s range resulting in a critical stiffness increase that creates a peak of potential instability in that velocity regime. The MFL leads to a new definition of fault frictional stability with implications for slip event styles and relevance for models of seismic rupture nucleation, propagation, and arrest.

Introduction

Friction laws based on one or more state variables are extensively used to model earthquake nucleation, seismic slip, and after slip [Ruina, 1980, 1983; Dieterich, 1981; Rice and Ruina, 1983; Gu et al., 1984; Okubo and Dieterich, 1984]. Most friction laws incorporate a dependence on slip history and have the following properties: (1) steady state dependence on velocity, (2) instantaneous dependence on velocity, and (3) evolutionary dependence on characteristic slip distance.

A typical example is the rate‐and‐state friction law (RSFL) [Dieterich, 1978; Ruina, 1983] which, in its compact form, can be stated as follows:

The friction coefficient μ is a function of sliding velocity V and of a state variable θ. The constitutive parameters a, b, and D (the characteristic slip distance) are determined by imposing velocity steps (from a fraction of µm/s up to a few cm/s) with respect to a reference slip rate V 0 and friction coefficient μ 0. The θ is often interpreted as an effective average time for the interchange of a population of asperity contacts or a time delay between a change in velocity and the corresponding change in friction coefficient. Its evolution is described by (slowness law [Dieterich, 1979]), for example,

Once the transient perturbation due to a velocity step is at completion, the steady state friction μ ss is achieved, dθ/dt = 0, and θ ss = D / V ss.. From equation (1),

According to equation (3), depending on the sign of a  −  b, either μ ss increases with increasing slip rate (rate‐hardening behavior, resulting in an unconditionally stable behavior) or μ ss decreases with slip rate (rate‐weakening behavior, potentially resulting in an unstable behavior leading to dynamic runaway).

The dependency of the steady state friction coefficient on slip rate (equation (3)) also has implications for the stability of a fault in terms of fault stiffness [e.g., Tse and Rice, 1986; Gu and Wong, 1991]. The steady state condition allows defining a critical stiffness: where σ is the effective stress normal to the fault and F(V, θ) is a generic term for inertia at a given state θ [Rice and Ruina, 1983; Ruina, 1983]. Under the assumption of negligible inertia (i.e., quasi‐static approximation) equations (3) and (4) yield k − a)σ. The ratio of the stiffness of the system k (either the natural fault plus the wall rocks or the experimental fault plus the apparatus) to k determines whether sliding is stable or unstable. For k  >  k, the system is stable but can become conditionally stable if perturbed sufficiently [Gu et al., 1984]; if k < k, sliding is always unstable. Variations in the condition of stability (k > k) may be induced if any of the constitutive parameters in equation (4) changes, i.e., the dependency of friction coefficient upon slip rate, a situation which is observed within the sliding velocity range discussed here. Experiments performed over a wide range of slip rates have shown that values of μ ss have a pronounced drop for V > 0.1 m/s [Di Toro et al., 2011; Goldsby and Tullis, 2011]. This drop corresponds to the onset of a dramatic velocity‐weakening dependence, with the consequence of substantially increasing k in equation (4) and potentially favoring the onset of unstable sliding (see section 4.2). Thus, equations (3) and (4) hold at subseismic velocities, but their validity is questionable at velocities approaching or exceeding a critical velocity V ~ 0.1 m/s.

Here to shed light onto the change of frictional regime which occurs at the transition from subseismic to seismic slip conditions, we investigate a large data set of published and unpublished data of friction experiments performed with rotary shear, biaxial and triaxial testing machines on silicate‐bearing, carbonate‐bearing, and quartz‐rich rocks (Table 1 and supporting information Table S1). Our aim is to define a μ ss curve whose range of validity extends from subseismic (V ≲ 0.1 m/s) to seismic slip rates (V ≳ 0.1 m/s) and discuss its implications for fault stability.

Table 1

Best Fit Parameters (a − b, V, and p) Resulting From Equation (5) (MFL, Figure 2) for a Given Range of σ, V ss, and Lithologya

Family	Rock Type	Referenceb	V ss (m s−1)	σn (MPa)	a − b	Vc (m s−1)	p
Quartz rich	Sandstone Novaculite	Di Toro et al. [2004, 2006 also unpublished data]	1E−10 ÷ 6.5	5 ÷ 18.7	−0.013 ÷ 0.002	0.016 ÷ 0.04	0.42 ÷ 0.86
Carbonate bearing	Calcite	This study, Han et al. [2007, 2010], Shimamoto and Logan [1981], Morrow et al. [2000] and also Di Toro et al. unpublished data	1E−7 ÷ 6.5	5 ÷ 200	−0.0007 ÷ 0.01	0.06 ÷ 0.12	0.78 ÷ 1.07
	Dolomite	This study, Shimamoto and Logan [1981], and Weeks and Tullis [1985]
Silicate bearing	Gabbro	Niemeijer et al. [2011], Nielsen et al. [2008], and Marone and Cox [1994]	3E−8 ÷ 6.5	0.5 ÷ 75	0.009 ÷ 0.035	0.02 ÷  0.093	0.44 ÷ 0.64
	Tonalite	Di Toro et al. [2004, 2006, also unpublished data]
	Monzodiorite	Mizoguchi and Fukuyama, [2010]
	Peridotite	Del Gaudio et al. [2009]
	Granite	Dieterich [1978] and Di Toro et al. [2004, 2006]

V = 1 µm/s; μ = μ (V).

When not original, data reported in Figure 2 are from the following sources (see also Table S1).

Figure 2

The MFL fit of experimental data under room humidity conditions for the three main rock categories: quartz‐rich, silicate‐bearing, and carbonate‐bearing rocks. Within the silicate‐bearing rocks cathegory (red dots) granites are in chess pattern and monzodiorites in stripes. Within the silicate‐bearing rocks cathegory (red dots) granites are in chess pattern and monzodiorites in stripes. Each rock category is pertinent to a specific fault lubrication style (see section 4). Experimental data and best fit parameters used in the MFL for each individual lubrication style are listed in Table 1. Given the data scattering for silicate‐bearing rocks in the intermediate velocity range (Figure 1b), the 95% interval of confidence is reported for completeness (dashed red lines).

Methods

We conducted experiments with the rotary shear apparatus SHIVA (Slow to High Velocity Apparatus), which imposes slip rates up to 6.5 m/s and normal stresses up to 40 MPa on solid rock samples with 30/50 mm internal/external diameter [Di Toro et al., 2010; Niemeijer et al., 2011]. We performed experiments with single (Figure 1a) and multiple velocity steps (Figure 1b) on silicate‐bearing (microgabbro) and calcite‐bearing (Carrara marble, 99% calcite) rocks. Samples were mounted inside an aluminum jacket, embedded within resin, and their surface rectified with a lathe [Nielsen et al., 2012]. The maximum height of the initial surface roughness was less than 10 µm for both rocks. The sample preparation procedure and the stiffness of SHIVA (k = 0.074 MPa/µm) resulted in a thorough reproducibility of the experimental data (Figure 1a).

Figure 1

Evolution of friction coefficient with slip and slip rate in experiments performed with SHIVA: (a) velocity control at seismic slip rates (Carrara marble, s308 and s307) and (b) velocity stepping in the intermediate velocity regime (microgabbro, s672) with typical pronounced oscillations of the friction coefficient at V = 0.1 m/s. Here the average value of the steady state friction coefficient was determined in the last 50% of cumulated slip resulting in a standard deviation Δμ ss = 0.2.

We measured the steady state friction coefficient at a constant velocity as depicted in Figure 1. Experiments are summarized in Table S1 together with data derived from a review of either published or unpublished data sets. The entire data set is synthetized in Table 1 with reference list.

Results

We studied the velocity dependence of μ ss and grouped the investigated lithologies under three rock categories (Tables 1 and S1): quartz‐rich rocks (novaculite and sandstone), silicate‐bearing rocks (microgabbro, monzodiorite, granite, and peridotite), and carbonate‐bearing rocks (calcite and dolomite). The frictional response to a velocity step varies slightly with rock category but preserves similar features (Figure 2). Within the entire velocity spectrum it is possible to distinguish three regimes for frictional sliding:

Low velocity (V < 10 − 4m/s): μ ss has a log linear dependence with V ss compatible with the RSFL (equation (3)) with a − b > 0 (velocity hardening) for silicate‐bearing rocks, a − b < 0 (velocity weakening) for quartz‐rich rocks, and a − b ~ 0 (velocity‐neutral) for carbonate‐bearing rocks.

Intermediate velocity (10−4  <  V ≲ 10− 2 m/s): the dependence of μ ss with V is negligible in case of silicate‐ and carbonate‐bearing rocks. However, in the experiments performed with SHIVA, this velocity range is associated to the onset of pronounced oscillations around a steady state value (Figure 1b). The amplitude of the oscillations Δμ (~0.2, e.g., Figure 1b) dampens down after either an increase toward seismic slip rates (V > 10−2 m/s) or a decrease to subseismic ones (V < 10−4 m/s). The Δμ is smaller than the one measured in experiments performed in the same velocity range on monzodiorite (e.g., Δμ ~ 0.5 [Mizoguchi and Fukuyama, 2010]), probably because of the higher stiffness of SHIVA.

High velocity (V ≳ 10 − 2 m/s): strong and negative velocity dependence of μ ss in all rock categories.

The velocity dependence of μ ss within the three regimes can be described by a nonlinear least squares empirical fit which results in a modified (rate‐and‐state) friction law (modified RSFL (MFL), solid curves in Figure 2): where V 0  =  1 µm /s and μ ≈  0.7 are reference values often used in the RSFL literature. Parameters V and p were determined with a best fit procedure resulting in the range of variability reported in Table 1. The a − b value was estimated from the evolution of μ ss in the μ ss‐log V plane, i.e., a  −  b  =  Δμ ss / Δlog(V ss). The best fit parameters resulting from equation (5) will be discussed in the next section.

The data scattering for silicate‐bearing rocks was large. The standard deviation is shown in Figure 1 for the experiment performed with SHIVA in the intermediate velocity regime but was not available for a large number of the published experiments. We account for the large data scattering of the silicate‐bearing rocks category by relaxing the best fit constrains to the 95% interval of confidence (dashed red lines in Figure 2).

Discussion

Rate‐Dependent Steady State Curve

The best fit model (MFL) of the entire data set indicates the existence of three velocity‐dependent frictional regimes with similar characteristics within each rock category. A critical velocity V exists that determines the transition from the intermediate to the high‐velocity frictional regime, which is likely dependent on material properties [Dieterich and Linker, 1992], as indicated by the variability of V with lithology (Table 1). The decay of μ ss with increasing slip rate in the high‐velocity regime (V > V) is described by the parameter p ranging from 0.4 to 1.1 depending on rock category. The two parameters, within each rock category, have a correspondence with the most accredited weakening mechanism operating at seismic velocities in the experiments (see Niemeijer et al. [2012] for a review). For cohesive silicate‐bearing rocks (microgabbro, granite, tonalite, peridotite, etc.), weakening is promoted by frictional melts which are produced easily at seismic slip rates [e.g., Tsutsumi and Shimamoto, 1997; Di Toro et al., 2004, 2006, 2011; Spray, 2005; Hirose and Shimamoto, 2005]. Theoretical interpretation of lubrication operated by friction melts predicts p ~ 0.5 [Nielsen et al., 2008, 2010], in agreement with the data presented here. Interestingly, V ranges from 1 to 9 cm/s in equation (5) to fit the experimental data for silicate‐bearing rocks. The estimated V is consistent with the critical velocities for the activation of flash heating and melting [Goldsby and Tullis, 2002, 2011; Violay et al., 2014] due to thermal softening at the asperity scale [see also Carlson and Langer, 1989; Rice, 2006; Tullis and Goldsby, 2003; Rempel, 2006; Beeler et al., 2008]. For quartz‐rich rocks (novaculite) in the absence of melting, fault lubrication has been associated to the amorphization and hydration of quartz in the presence of room humidity at the highly stressed contacts where strained Si‐O‐Si bonds are highly reactive (silica gel lubrication [Goldsby and Tullis, 2002; Di Toro et al., 2004; Nakamura et al., 2012; Kirkpatrick et al., 2013]). Hydrated amorphous silica was found on the slip surface of experimental faults of single crystals and polycrystalline quartz sheared at V > 3 mm/s [Hayashi and Tsutsumi, 2010; Nakamura et al., 2012]. This velocity is almost in the range of the V reported in Table 1. For carbonate‐bearing rocks weakening has been associated with decarbonation, carbon amorphization, and thermally activated grain size‐dependent processes [Han et al., 2007, 2010; Verberne et al., 2014; De Paola et al., 2015; Green et al., 2015; Spagnuolo et al., 2015]. Also, for this rock category, V (ranging from 8 to 13 cm/s) is compatible with the critical slip rate measured in experiments performed on calcitic marbles and limestones [e.g., Tisato et al., 2012]. For these rocks, the critical velocity has been interpreted in the framework of the fast‐moving dislocation theory [Spagnuolo et al., 2015].

The friction coefficient data of silicate‐built rocks (e.g., granite and gabbro) are highly scattered due to the following:

Mechanical properties of their forming minerals and respective melts (i.e., fracture toughness, melting temperature, and melt viscosity [Spray, 1993, 2010]).

The nonlinear dependence of shear stress with normal load leading to a nonconstant friction coefficient over the range of normal loads considered here (0.5 to 75 MPa). This is the case for frictional melt lubrication, where the shear stress has a power law dependency on normal stress with an exponent ranging from 0.25 to 0.5 [Nielsen et al., 2008; Niemeijer et al., 2011; Violay et al., 2014].

The dependence of the friction coefficient on slip. In Figure 2, we included all the friction data independently of the cumulated slip. However, for rough surfaces of gabbro, Marone and Cox [1994] reported μ = 0.45 at slip rates ranging from 1 to 10 µm s−1 after a slip < 50 mm and μ = 0.7 after cumulative slip of 53 mm.

Overall, the empirical steady state friction law equation (5), tested over a large data set with different rocks and machines, satisfies the theoretical prediction of an enhanced velocity‐weakening condition which was introduced for dynamic fracture simulations [Zheng and Rice, 1998] and in models of seismic rupture propagation [e.g., Noda et al., 2009, Dunham et al., 2011] by setting p  =  1.

Frictional Stability and Critical Stiffness Under Variable Rate Dependence

The role of the stiffness in the dynamics of systems controlled by friction has been extensively discussed in theoretical and numerical studies [e.g., Dieterich, 1979; Rice and Ruina, 1983; Putelat et al., 2011; Belardinelli and Belardinelli, 1996; Tse and Rice, 1986; Gu and Wong, 1991]. The fault + wall rock system was described by spring‐slider models where a sliding mass is loaded through a spring by a driver moving at constant velocity V ss. The frictional instability is triggered if dμ ss / dV  <  0 and k ≤ k, where k, is the critical stiffness of the apparatus or of the wall rocks in nature [Rice and Ruina, 1983]. The k is defined in equation (4) under the quasi‐static approximation which implies that velocity is the fastest evolving variable of the system and it responds instantaneously (negligible inertia) to a variation of the boundary conditions [Gu et al., 1984; Belardinelli and Belardinelli, 1996]. Assuming that the quasi‐static approximation holds, we compute k according to equation (5) for V ranging from subseismic to seismic slip rates for σ = 10 MPa and D = 1 µm, 10 µm, and 100 µm. The k (V) is plotted in Figure 3a together with the machine stiffness of SHIVA (k = 0.074 MPa/µm) and of other machines (0.1 MPa/µm [Byerlee and Brace, 1968], 0.04 MPa/µm [Dieterich, 1979], and 0.105 MPa/µm [Marone et al., 1990]). For silicate‐ and carbonate‐bearing rocks given a typical laboratory D ~1 µm k (V) intersects the machines stiffness in the gray shaded area at V ~ 10−4–10−3 m/s, switching the experimental fault + machine system from stable (k) to unstable (k) sliding. As velocity further increases and once the intermediate velocity regime is overcome, dμ ss/dV begins to wane and k decreases until eventually k so that the system returns to a condition of stability in the high‐velocity regime. We note that at high velocity the inertial terms may cease to be negligible [Baumberger et al., 1994]. In this regime, the inertial terms are essentially stabilizing, provided that both dμ ss /dV and ∂μ/∂V are negative, which is what we observe experimentally at high velocities. If the critical stiffness is less but close to the stiffness of the apparatus (intermediate velocity regime), permanently sustained oscillations occur as observed for the example in Figure 1b and previous experiments [e.g., Tullis and Weeks, 1986; Mizoguchi and Fukuyama, 2010], in agreement with theoretical linearized stability analysis [e.g., Gu et al., 1984].

Figure 3

The critical stiffness k from equation (4) versus machine (or natural) stiffness k for the three rock categories. (a) Laboratory faults: k is computed using the MFL (equation (5)) for σ = 10 MPa. Dashed, solid, and dash‐dotted lines are for D = 1 µm, 10 µm, and 100 µm, respectively. The k = 0.074 MPa/µm for SHIVA is in blue solid line; a range of other machine stiffness corresponds to the gray shaded area. For k < k frictional sliding is stable; for k ≈ k and V ~ 0.1 m/s (intermediate velocity regime) theoretical stability analysis predicts the onset of oscillations, consistently with experimental observations (Figure 1b). For k > k (V > V) the friction coefficient drops, and inertial effect may become dominant. (b) Natural faults: k results from the extrapolation of the MFL to natural conditions (σ = 100 MPa, D = 0.5 mm). A range of k for natural faults (red shaded area) was computed assuming a fault patch radius between 5 and 50 m and shear stiffness of G = 25 GPa.

In general, from Figure 3a, given a D value, the intersection between k and k indicates the velocity regime where frictional instability can occur. Conversely, given a value for stiffness, a − b, and a velocity regime where oscillations occur, it is possible to determine an approximate value for D. However, D is not necessarily constant in the three velocity regimes; thus, the estimated D is the ultimate value prior to the onset of the frictional oscillations. In Figure 3a each k curve represents one stability status of the system given an initial population of asperities whose renewal requires a distance D If D increases with velocity (or slip) the intersection between the k and k curves shifts toward higher velocities.

The evolution of effective k throughout the slip velocity spectrum, as described above, has a number of implications for fault stability. For example, it can explain why the realization of friction experiments at either high (~1 m/s or more) or low (<1 cm/s) velocity is less challenging than at intermediate velocity. Indeed, at intermediate velocity the stability, and hence control of the system, are more difficult to achieve due to the peak in k (Figure 1b). In addition, the velocity dependence of k indicates that variability in slip patterns (stable sliding, slow slip pulses, and stick slip) can be expected on a given fault, as observed in the experiments here, and will possibly contribute to the variety of slip behavior observed on natural faults.

Implications for Fault Stability

From Figure 3a, when D ≥ 100 µm, k curves rarely intersect the machine stiffness zone implying that experimental faults with D greater than the upper limit observed in experiments are in a condition of stable slip over the entire velocity spectrum. However, in the case of natural faults, all other variables in equations (4) and (5), and the stiffness of the system should be rescaled. Though the stiffness of natural faults varies sensibly throughout the evolution of the fault zone with repeated slip events [Griffith et al., 2009], a rough estimation for the stiffness right before the earthquake derives from [Eshelby, 1957] using the simplest approximation of a buried circular crack of radius r and shear modulus G and assuming that the earthquake instability does not propagate beyond the nucleation patch. We can use indicatively G = 25 GPa (e.g., granite) and r ranging from 5 to 50 m resulting in k = [0.7–7] 10−3 MPa/µm per unit fault area, up to 1  order of magnitude smaller than the stiffness of SHIVA. The conditions of stability for faults can be extrapolated as reported in Figure 3b, computing equation (4) with σ = 100 MPa (approximatively the lithostatic load at 4 km depth) and D = 0.5 mm. Again, oscillations are expected when k is close to k. As a result of the scaling shown in Figure 3b, it appears that a natural fault patch is not necessarily locked into either stable or unstable conditions. A given fault may evolve within three different regimes, depending on slip velocity. If slip rate gradually increases (for example, under the effect of an increased tectonic load) the fault will first encounter a stable slip episode which will last as long as k (in Figure 3b, slip rate below 10−4 m/s for the case where D = 0.5 mm, red solid curve). As the fault slip accelerates, k may increase to a point where it approaches fault stiffness k. Then the fault may respond to small perturbations by oscillations or stick‐slip motion, a sign that it is close to instability limit k ~ k. As the fault accelerates further to k, it may eventually dynamically weaken and enter a seismic transient episode. The above sequence of events may represent a type of precursory pattern, with accelerating slip and moderate seismic unrest preceding an earthquake rupture. We note that in addition to the evolution of k described above, another well‐known mechanism for the transition from stable to unstable sliding is the growth of the slip patch radius r, with the consequent decrease in k according to equation (6).

Conclusions

A modified friction law (MFL, equation (5)) describes the steady state conditions over the entire velocity spectrum where we identify three velocity regimes for frictional behavior. Two parameters describe the three regimes, the critical velocity V and the exponent p (Table 1) which are material dependent and likely related to specific physicochemical processes occurring at the asperity scale and within the slipping zone.

First, according to the MFL, critical conditions for laboratory fault stability (k and D) are not constant but evolve with slip rate, and their evolution can justify slip accelerations on materials expected to slip stably (e.g., a − b > 0 for silicate‐bearing rocks in Figure 3a). Consequently, the parameter a − b alone cannot be intended as a watershed between the unconditionally stable and potentially unstable frictional behaviors. Noda and Lapusta [2013] achieved similar conclusions for fault stability using rapid shear heating of pore fluids to combine the stable slip to the fast coseismic weakening.

Second, the MFL describes two possible steady state conditions for a given μ ss, at least in the case of silicate‐bearing rocks of the present data set (Figure 2), one in the low‐velocity and one in the high‐velocity regime. Consequently, experimental faults slipping stably at low velocities can jump into the high‐velocity regime under a “velocity kick” [Scholz, 1998] that forces the fault to overcome a critical slip or stiffness.

Third, k has a pronounced peak at intermediate slip velocities (i.e., V  ~  V) and for D = 10 µm the estimated k is larger than the stiffness of all rotary shears. The k peak explains why experiments performed in the intermediate velocity regime are extremely demanding resulting often in sample failure and do require machines with very high stiffness for control purposes. However, at higher velocities (V > 0.1 m/s) k drops again below the typical machine stiffness, and stable behavior is expected. Consequently, the high‐velocity regime represents a state of stable sliding in the sense of perturbation theory, resulting in manageable experimental control.

The scaling of stability conditions to natural faults is discussed under the opportune scaling of stiffness and D. In nature, at seismogenic depths (e.g., σ = 100 MPa), oscillations attributed to frictional instabilities are expected to occur in a range of slip rates (10−4–10−3 m/s) that is lower than the range (10−3–10−2 m/s) of experimental faults. The estimates given here are qualitative because of the strong assumptions made on natural fault stiffness.

Differently from previous theoretical studies the MFL here proposed is based on experiments performed on different rock types with several machines. The MFL predicts a variable critical stiffness that describes the behavior of faults under dynamic loading conditions in association with the occurrence of at least three frictional regimes and slip event styles. Accordingly, under the opportune scaling, the MFL is effective in explaining the plethora of events observed during the seismic cycle and might be applied to both observational seismology and earthquake modeling.

Table S1

Click here for additional data file.