Earthquake nucleation on rough faults

Earthquake nucleation is currently explained using rate and state stability analysis, which successfully models the behavior of laboratory simulated faults with constant thickness gouge layers. However, roughness is widely observed on natural faults and its influence on earthquake nucleation is little explored. Here we conduct frictional sliding experiments with different roughness on granite samples at upper crustal conditions (30–200 MPa). We observe a wide range of behaviors, from stable sliding to stick slip, depending on the combination of roughness parameters and normal stress. Stick slip is repeatedly observed in velocity-strengthening regimes, and increases in normal stress stabilize slip; these features are not fully predicted by current stability analysis. We derive a new instability criterion that matches our observations, based on fracture energy considerations and the size of weak patches created by fault roughness.


INTRODUCTION
A central question regarding tectonic faults concerns the onset of earthquake-generating stick slip as opposed to aseismic stable sliding. This problem has been addressed in observational (Dodge et al., 1996), theoretical (Rice and Ruina, 1983), and experimental studies (Leeman et al., 2016;Scuderi et al., 2016) using the predictions of rate and state friction laws and stability analysis, where instability develops under velocity-weakening friction and low mechanical stiffness (Leeman et al., 2016;Marone, 1998;Scuderi et al., 2016). In general, most experiments have been conducted on homogeneous materials, either generating slip in a constant thickness gouge layer (e.g., Leeman et al., 2016;Scuderi et al., 2016) or on roughened cohesive rock surfaces (e.g., Passelegue et al., 2013). However, natural faults are highly heterogeneous features with variable composition, physical properties, and complex slip surface geometries (Bistacchi et al., 2011;Candela et al., 2012;Sagy et al., 2007).
In this study, we investigate the effects of heterogeneity due to the roughness of fault surfaces, and its influence on the onset of unstable sliding. Roughness is observed on faults at all scales (Bistacchi et al., 2011;Candela et al., 2012;Sagy et al., 2007), and plays a key role in faulting by controlling the size and distribution of asperities (Dieterich and Kilgore, 1994;Scholz, 1988) and the stress distribution on the fault surface (Persson, 2013;Selvadurai and Glaser, 2017). Therefore, roughness should have significant implications for both the static  and dynamic frictional strength of fault zones (Fang and Dunham, 2013), critical slip distances Ohnaka and Shen, 1999;Okubo and Dieterich, 1984), and nucleation size (Ohnaka and Shen, 1999;Okubo and Dieterich, 1984). Currently only a narrow range of conditions has been investigated experimentally (Marone and Cox, 1994;Ohnaka and Shen, 1999;Okubo and Dieterich, 1984).
Here we report results of the first systematic experimental study investigating the occurrence of frictional instability under a range of roughness and normal stress conditions. We show that the combination of these parameters controls the onset of frictional instability of faults. We then provide a novel explanation based on the interaction between maximum weak patch scaling, roughness, and normal stress.

METHODS
We performed 23 frictional sliding experiments with Westerly Granite (northeastern United States) subjected to stress conditions representative of Earth's upper crust (30 MPa ≤ σ n ≤ 200 MPa; σ n is normal stress) and loaded in a direct shear configuration to force sliding. Samples were subjected to an initial period of 1-1.5 mm run-in before carrying out velocity steps between 0.1 and 10 µm s -1 . To create different distributions of heterogeneity on the simulated faults, the sliding surfaces were axially precut and carefully polished to obtain variable degrees of roughness, characterized in terms of root mean square roughness (Z rms ) using stylus profilometry measurements (see Sections DR1 and DR2 in the GSA Data Repository 1 ). 1 GSA Data Repository item 2017311, summary of experimental conditions (sections DR1 and DR4), surface topography measurements (section DR2), mechanical dataset (section DR3), microstructural images (section DR5), derivation of Equation 4 (section DR6), and details of elastic modelling (section DR7), is available online at http://www.geosociety .org /datarepository /2017/ or on request from editing@geosociety.org.

RESULTS
All experiments show initial elastic loading followed by frictional rollover, where the contacting surfaces begin to slide (see Fig. DR1 and the summary of results in Sections DR3 and DR4). Once past this initial stage, the frictional strength remains relatively constant and a steady state is reached (typically requiring a displacement of 0.75-1.5 mm). The full spectrum of frictional sliding behaviors is observed, from stable sliding to seismic stick slip, across the range of experimental conditions. In several experiments, it was possible to determine the rate and state friction parameters a and b by modeling the frictional data to load-point velocity stepping during stable sliding episodes (Section DR3). Figure 1 shows examples of typical slip dynamics observed in different experiments (full plots are provided in Section DR3).
When normal stress is increased to 100 MPa, smooth faults (Z rms ≤ 4.3 µm) are observed to become fully unstable with repetitive fast stick slip instabilities (Fig. 1C). Fast slip is confirmed by observations of frictional melting in scanning electron microscopy (SEM) imaging of slip surfaces ( Fig. 2C; Section DR5). Intermediate roughness surfaces (Z rms = 8 µm) show marginal stability with velocity-weakening to neutral friction accompanied by slow stress drops upon increases in velocity (Fig. 1D). Rougher faults (Z rms ≥ 18.4 µm) are stable throughout the course of experimentation with velocity-strengthening friction, and abundant cataclasis observed in SEM imaging ( Fig. 2D; Section DR5).
For σ n > 100 MPa, sliding shows a wider spectrum of behaviors, with some unexpected results. At 150 MPa, smooth faults (Z rms ≤ 4.3 µm) remain unstable with repetitive fast stick slip cycles. All rougher faults (Z rms > 4.3 µm) are marginally stable, with evidence of fast stress drops nucleating spontaneously (without a velocity kick) or upon stepwise velocity increases, in spite of velocity-strengthening friction measurements ( that increasing the normal stress to 200 MPa results in the smoothest faults (Z rms ≤ 4.3 µm) becoming marginally stable (Fig. 1F). Similar behavior is also observed on intermediate roughness faults (4.3 < Z rms < 28.2 µm), which are stable with velocity-neutral to velocity-strengthening friction. Unexpectedly, given the consistently velocity-strengthening friction at lower normal stresses, the roughest fault (Z rms = 28.2 µm) is unstable with repetitive dynamic stick slip (Section DR3).
The variety of slip behaviors observed is summarized in Figure 2A, where points correspond to experimental conditions (Z rms σ n ), allowing approximate definition of differing frictional domains. Two characteristic trends emerge in the data (Fig. 2). First, there is a transition from stable to unstable and marginally stable slip as normal stress is increased, in accordance with the predictions of rate and state (Marone, 1998;Rice and Ruina, 1983), with the transition at increasingly higher normal stress as faults become rougher.
Second, with further normal stress increase up to 200 MPa, instability is suppressed on all but the roughest fault, which becomes unstable ( Fig.  2A). The occurrence of spontaneous rupture nucleation in a velocity-strengthening regime for several experiments and the second trend of the stabilizing effect of normal stress are not predicted using a standard stability analysis (Marone, 1998;Rice and Ruina, 1983).

DISCUSSION
We now discuss these results in light of rupture stability criteria and develop a theoretical model based on roughness-induced weak fault patches. Studies of natural fault surfaces show that faults have a characteristic self-affine roughness, described by a power density spectrum over nine orders of magnitude (from 10 -4 to 10 5 m): where α is the amplitude scaling (in m 3 ), k is the inverse wavelength (in m -1 ), H is the Hurst exponent, and k 0 is a normalizing factor (here 1 m -1 ) with α = 10 -3 to 10 -1 m 3 and H = 0.6-0.8 (Bistacchi et al., 2011;Candela et al., 2012). At shorter length scales of <1-50 µm, this scaling diminishes and becomes isotropic as a result of plastic yielding at asperity contacts . From stylus profilometry measurements of our pre-experimental and postexperimental surfaces, we have a corner inverse wavelength, k min , identified using Fourier analysis (see Section DR2), above which surfaces obey self-affine scaling.   (Section DR6). Therefore, Z rms is selected as a single representative parameter of the surface statistics in Figure 2, and throughout the discussion.
Onset of rupture propagation can be interpreted either (1) in the context of rate and state dependent friction (Marone, 1998;Rice and Ruina, 1983), when a sliding patch reaches a critical size; or (2) as the consequence of stress concentration around a weak patch that may propagate unstably according to fracture energy considerations (Griffith, 1921), which have been adapted to the problem of shear cracks and earthquake faulting (Andrews, 1976;Ida, 1972).
According to criterion 1, stability is controlled by the ratio of the mechanical stiffness K f to the frictional stiffness K c , defined as , where a and b are rate and state friction dimensionless parameters and D c is the critical slip distance. When the stiffness criterion K f K c < 1 is satisfied, instability can develop; other wise, sliding is conditionally stable (Marone, 1998;Rice and Ruina, 1983). In the case of tectonic faults embedded in an elastic medium, K f represents the stiffness of the fault and can be expressed as the shear modulus, h is the linear fault dimension, and C is a dimensionless shape factor (Rubin and Ampuero, 2005). The stiffness criterion allows definition of the minimum dimension h*, of a slip patch required for instability to develop. Rate and state friction laws and Equation 2 provide effective tools to model the full spectrum of fault slip behaviors observed across relatively homogeneous sliding interfaces such as gouge-dominated faults (Leeman et al., 2016;Scuderi et al., 2016). However, stability criterion 2, based on fracture energy, surmises the presence of a pre-existing flaw or weak patch of finite size. Material flaws are inherent in Griffith's (1921) original crack theory, and in the case of tectonic faults we may equate them to an elastic bridge between asperities (Figs. 3A and 3B). For earthquake nucleation, instability arises when the growth of the weak patch is energetically favorable, requiring that the strain energy release balances or exceeds the work done against residual frictional strength. According to a model of dynamic slip weakening, fault friction drops over a characteristic distance δ c (Andrews, 1976;Ida, 1972); here we extend the significance of δ c to low sliding velocities.
Criterion 2 allows the definition of a critical length, L c = C G c p r 0 r ( ) 2 at which a shear crack undergoes unstable failure (Andrews, 1976), where τ 0 is the static shear stress on the fault, τ p = µ p σ n is the peak stress, µ p is the peak friction coefficient, and τ r = µ r σ n is the shear stress after weakening where µ r is the sliding or weak friction coefficient. To enhance similarity with h* (Equation 2), we derive a lower bound length estimate for L c by assuming that the stress state on the fault is close to the peak stress during experiments (i.e., τ 0 ≈ τ p ), yielding Although the critical patch length h* and L c share some scaling similarities, they may differ by orders of magnitude: a -b is usually small (typically −0.01) while we can expect µ p -µ r to be quite large. Here the frictional strength at asperities equates to µ p , and within elastic bridges or zones of reduced asperity density it equates to µ r . We adopt µ p -µ r = 0.2, as suggested by the observations of Selvadurai and Glaser (2017) on rough surfaces, which show that stress fluctuations can be as much as 40% of the peak stress.
Estimates of nucleation size for experiments showing stick-slip instability at moderate normal stress using Equation 2 for a fault of Z rms = 3.6 µm [a -b = −0.003, D c = 5 µm, G = 50 GPa (estimated from loading curves), and σ n = 100 MPa] yield h* values of ~1 m, ~2 orders of magnitude larger than the size of samples utilized. Following Rubin and Ampuero (2005), if we neglect the rate parameter a, we find that h b * reduces by a factor of 3 (a = 0.005), which is still an order of magnitude larger than the size of the sample. If we calculate the nucleation length using Equation 3, we obtain L c = 1.25-3.75 cm (with G = 50 GPa, σ n = 100 MPa, µ p -µ r = 0.2, and δ c = 5-15 µm). Values of L c obtained are in agreement with other studies that posit that the nucleation length should be smaller than the sample length (L 0 = 4 cm) for laboratory stick slip occurrence (Ohnaka and Shen, 1999;Okubo and Dieterich, 1984;Passelegue et al., 2013). The values estimated here for δ c = 0.05 k min -1 (Ohnaka and Shen, 1999) are consistent with previous estimates using high-frequency strain gauges (Okubo and Dieterich, 1984) and those predicted by numerical modeling of elastic surface closure (see Section DR7). This is in contrast to values of D c obtained during velocity steps that do not show a systematic dependence on roughness. The onset of instability observed at higher normal stress for increasing roughness is in accord with δ c ∝ k min -1 . While a stability criterion based on fracture energy (e.g., Equation 3) can explain the onset of stick slip during our experiments at low to moderate normal stress (30-150 MPa), the surprising observation that slip instability is suppressed at higher normal stress indicates the presence of some limiting process ( Fig. 2A). As we discuss in the following, this behavior could be explained by considering the microphysical properties of contact asperity distribution in relation to fault zone roughness, and the associated stress heterogeneity (Scholz, 1988).
Previous studies imaged the distribution of frictional contacts with increasing normal stress and varying surface roughness in transparent materials (Dieterich and Kilgore, 1994;Selvadurai and Glaser, 2017). With increasing normal stress, contact asperities increase in number and grow, as shown in Figure 3 (Section DR6). Theory indicates that stress and asperity sizes will follow a power law distribution for a self-affine surface under load (Scholz, 1988). The asperity bridging length, λ c , which is the maximum supportable elastic length or bridge between asperities, is also shown to decrease as λ c ∝ σ n -2 for a self-similar surface (Scholz, 1988). Here we consider a generalization of this result to any self-affine surface with 0 < H < 1 as is a scaling factor of length dimension (Section DR6). From measurements of experimental fault surfaces, the Hurst exponent is typically 0.6-0.9 above k min , yielding λ c ∝ σ n -2.5 -σ n -10 . In comparison, Equation 3 gives L c ∝ σ n -1 , demonstrating that as normal stress increases, the bridging length will decrease at a faster rate than that of the nucleation length (Fig. 4). In extreme cases, bridges of length scale λ c may represent voids, as shown in Figure 3, but more generally represent zones of reduced normal stress, or weak patches of low stiffness, filled with undercompacted gouge. These weak zones can act as stress concentrators and initiate rupture, provided that L c < λ c , as shown in Figure 4. However, with increasing normal stress, the bridges will gradually be closed and the maximum open patch will decrease until λ c < L c , and rupture nucleation is no longer possible in accordance with our experimental observations (Fig. 2). In general, instability leading to rupture nucleation vin our experiments is only observed when the nucleation length L c satisfies the condition L c < λ c and L c < L 0 . Conversely, the conditions L c > L 0 at lower normal stress inhibits stick slip; λ c < L c at higher normal stress leads to stable sliding (Fig. 4).

CONCLUSIONS
Our findings also have implications for the larger scale behavior of natural fault zones. In principle, the model shown in Figure 4 suggests that the transition from seismic to aseismic faulting may be controlled by the stabilizing influence of increasing normal stress upon asperities, as originally suggested by Brace and Kohlstedt (1980), in addition to currently accepted temperature-induced rheological changes (Scholz, 1988). In addition to this our results qualitatively support observations of subduction zone seismicity, where rough seafloor topography is observed to be related to creeping behavior, and smooth seafloor topography is related to seismicity and large earthquake nucleation (Wang and Bilek, 2014).
Our results highlight the key role of fault heterogeneity in earthquake nucleation. On larger scales such heterogeneity includes fault jogs, compositional contrasts, and fluid injection in addition to fault roughness, as suggested herein. These results complement rate and state friction stability analysis, providing a physical framework to include the complexity of roughness in earthquake nucleation models.