Illuminating the physics of dynamic friction through laboratory earthquakes on thrust faults

Seismological Laboratory, Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125; Department of Earth and Environmental Sciences, Ben-Gurion University of the Negev, Beer-Sheva 84105, Israel; Graduate Aerospace Laboratories, California Institute of Technology, Pasadena, CA 91125; and Division of Engineering and Applied Science, California Institute of Technology, Pasadena, CA 91125

Large, destructive earthquakes often propagate along thrust faults including megathrusts. The asymmetric interaction of thrust earthquake ruptures with the free surface leads to sudden variations in fault-normal stress, which affect fault friction. Here, we present full-field experimental measurements of displacements, particle velocities, and stresses that characterize the rupture interaction with the free surface, including the large normal stress reductions. We take advantage of these measurements to investigate the dependence of dynamic friction on transient changes in normal stress, demonstrate that the shear frictional resistance exhibits a significant lag in response to such normal stress variations, and identify a predictive frictional formulation that captures this effect. Properly accounting for this delay is important for simulations of fault slip, ground motion, and associated tsunami excitation.
dynamic friction | laboratory earthquakes | thrust faults M any large, destructive, tsunamigenic earthquakes occur along thrust faults, such as the 2011 Mw 9.0 Tohoku earthquake in Japan (1) and the 1999 Mw 7.7 Chi-Chi earthquake in Taiwan (2). An important feature of these two earthquakes is that the maximum slip occurred close to the free surface (3)(4)(5) (Fig. 1), despite the notion that shallow portions of the faults are supposed to be frictionally stable (e.g., ref. 6). Numerical and theoretical studies have showed that the interaction of the ruptures with the free surface results in rapid variations in the fault-normal stress, σ (7-11), which feed back into the rupture process. The fault-normal stress affects the fault frictional shear resistance, which controls the fault resistance to sliding and, as a consequence, how far and how fast a fault slides. For example, the unexpectedly large slip of 60 to 80 m in the shallow portion of the 2011 Mw 9.0 Tohoku earthquake may have been caused by this phenomenon (12,13). While friction traditionally is defined as being proportional to normal stress, experimental studies (14)(15)(16)(17)(18) have shown that frictional shear resistance does not change instantaneously in response to a step changes in σ but rather evolves gradually with slip. This effect, and its proper representation in friction formulations, is critically important to understanding the dynamics of ruptures on thrust and normal faults, as it affects near-fault shaking and tsunami generation (7-10, 12, 13) Recent advances in imaging spontaneous dynamic ruptures in the laboratory (19)(20)(21)(22)(23) enable us to image and quantify the normal stress reduction near the free surface and its effect on friction for realistic fault-slip histories similar to natural thrust earthquakes. We report and interpret measurements from our laboratory earthquake experiments conducted in a setup, initially designed in refs. 13 and 24 to mimic thrust earthquakes. The original version of this setup used dynamic photoelasticity to qualitatively visualize the rupture process while particle velocimeters were employed to record the resulting ground motion and to measure fault opening and slip at only a few discrete locations on the fault and the free surface (13,24,25). The current version, of this setup, described here, uses a newly developed, high-resolution, dynamic imaging technique based on a combination of high-speed photography and digital image correlation (DIC) (19)(20)(21)(22)(23). This method allows us to monitor the entire rupture process within the field of view (FOV), to follow individual ruptures as they approach and break the free surface ( Fig. 1A and SI Appendix, Fig. S1), and to analyze this entire process through, previously unavailable, full-field measurements of displacements, velocities, strains, and stresses at time intervals of 1 μs (Materials and Methods and SI Appendix). The analog material used in the experiments results in critical length scales for instability that are one to two orders of magnitude smaller than those of rocks, allowing us to reproduce spontaneously propagating dynamic ruptures in laboratory samples of manageable size.

Dynamic Rupture Interaction with the Free Surface
The experiments illustrate in vivid detail how the thrust ruptures interact with the free surface while dynamically capturing all of the quantitative information needed to constrain fault-normal stress variations and friction formulations. In a representative experiment (experiment 1, Fig. 2), a rupture approaches the free surface with a peak particle velocity magnitude of 4 m/s at a time t = 55 μs after nucleation. Note that the velocity map so far is symmetric with respect to the fault line. The velocity field of the propagating rupture shows a transition from a dominant faultnormal motion ahead of the rupture front to a fault-parallel motion behind it. As the rupture arrives at the free surface (t = 62 μs), the symmetry in velocity field is broken, and the peak particle velocities increase to 8 and 4.5 m/s in the hanging and foot walls, respectively. At a later time (t = 75 μs), sliding Significance Our study explores a major challenge in earthquake sciencethe dynamics of thrust earthquake ruptures as they interact with the Earth's surface, which is relevant to some of the most destructive earthquakes that have ever occurred. The work illustrates key features of the complex dynamic behavior associated with this interaction using an ultrahigh-speed imaging technique. The interaction leads to large and rapid reductions in normal stress. However, the frictional shear resistance does not decrease instantaneously with normal stress, as typically assumed, but experiences a significant delay. Such delay is important for a range of earthquake source problems that involve rapid normal stress variations.
continues with smaller levels of particle velocity. Note that the motion of the hanging wall is mostly vertical, while the footwall has a large horizontal component. The supershear rupture just described is followed by another rupture, traveling at a sub-Rayleigh speed. This sub-Rayleigh rupture originated the supershear one (26), and subsequently trails behind it (22,24,27). As the trailing-Rayleigh signature arrives to the free surface (t = 101 μs), particle velocities increase, with larger values at the hanging wall. The particle velocities of the footwall are almost completely horizontal near the free surface. Full-field images of the fault-parallel velocity, shear stress, τ, and normal stress, σ (Fig. 3A), provide further insight into the rupture process. As the rupture propagates through the FOV, τ shows a symmetric pattern with respect to the fault. The shear stress increases above the initial stress level (τ 0 = 6.4 MPa) at the rupture tip and decreases to smaller values behind it. The normal stress and fault-parallel velocity show asymmetric pattern with respect to the fault.

Temporal Evolution of Stresses and Slip Rate near the Free Surface
The pronounced reduction in the fault-normal stress σ and the associated response of shear resistance τ are captured by plots of time histories of τ, σ, and the slip rate V at a point on the fault close to the free surface (Fig. 3B). As the rupture arrives, σ initially increases from σ = 11.6 to σ ∼ 13.5 MPa (t = 59 μs), and then decreases to σ ∼ 8 MPa over about 11 μs. At the arrival of the trailing Rayleigh rupture (t = 102 μs), there is additional reduction to a minimum of σ ∼ 5 MPa. The shear stress increases to τ ∼ 7 MPa at the arrival of the supershear rupture, then decreases rapidly to τ ∼ 4 MPa over about 3 μs. At later stages (t > 70 μs), there are only small variations in τ, with a total decrease of 1 MPa over 70 μs. The slip rate shows two peaks of V = 12 and 8 m/s, corresponding to the arrival of the supershear rupture and the trailing Rayleigh signature, respectively.
Interestingly, supershear transition later in the rupture process leads to a more intense near-surface rupture. An experiment (experiment 2, Fig. 4B) under similar loading conditions to those in experiment 1, but with a weaker nucleation process, have resulted in a later supershear transition and a later arrival of the rupture to the free surface (t ∼ 82 μs). Because of the later transition, most of the rupture energy in experiment 2 is still carried by the trailing Rayleigh signature rather than the (newly formed) supershear rupture front. Thus, the peak in V and reduction of σ associated with the first arrival of the rupture to the free surface are lower than in experiment 1. However, as the supershear rupture is reflected, it interacts with the arriving trailing Rayleigh, resulting in a significantly more intense trailing Rayleigh in experiment 2 than both the supershear rupture front and the trailing Rayleigh in experiment 1, with peak V of 16 m/s and reduction of σ to a minimum of 2 MPa.
To create different types of ruptures and gather a suite of data suitable for constraining friction formulations, we have conducted experiments at different initial levels of normal stress σ 0 (SI Appendix, Table S1). In the already presented experiment 1, σ 0 = 11.6 MPa. The temporal evolution of τ, σ, and V near the free surface in experiments under a prestress level of σ 0 = 7.6 and 7. In the experiments under lower initial normal stress, σ is completely released, for a short period of time, as the trailing Rayleigh signature arrives to the free surface ( Fig. 3C and SI Appendix, Fig. S2). This indicates a potential local opening of the fault. However, the shear resistance never goes to zero and the normal stress returns to compressive values shortly after, indicating that the actual opening may not have occurred at least for the set of loading parameters and the rupture scenario of supershear rupture approaching the free surface as reported here. Note that fault opening near the point where rupture meets the surface was suggested in previous experiments based on laser velocimetry measurements of the fault opening velocity (13) but under much lower initial compression resulting in a very different rupture scenario than the one investigated here (i.e., purely sub-Rayleigh rupture in ref. 13 vs. supershear rupture with sub-Rayleigh signature here).

Frictional Response under Rapid Normal Stress Variations
Our experimental measurements indicate that the shear resistance does not obey the traditionally assumed proportionality to the normal stress but evolves gradually, as suggested in refs.
14-16 and 18. This delay is directly observed in plots of the effective friction, τ/σ, vs. slip near the free surface (Fig. 4). For experiment 1, the effective friction initially increases to τ/σ ∼ 0.6 and then decreases with slip to τ/σ ∼ 0.35 at slip of about 25 μm. At larger levels of slip, when the impinging rupture is reflected at the free surface, σ decreases, and because of the delayed response of the frictional shear resistance τ, the ratio τ/σ increases back to a value of 0.6 at a slip of 120 μm. The friction τ/σ gradually decreases at larger slip, but as the trailing Rayleigh arrives and σ temporarily decreases, τ/σ increases again to a peak of 0.7, and later drops to 0.4.
The measurements of τ, σ, and V along the interface enable testing different formulations of frictional shear resistance, as well as constraining their parameters. We find that friction formulations without the delayed evolution of shear stress in response to normal stress changes cannot fit our experimental measurements.
Friction Model 1. In this model, we test a formulation featuring rate-and-state friction with enhanced weakening but without accounting for delayed shear stress response to variations in normal stress. The formulation assumes that the frictional shear resistance is proportional to σ as follows: where f is the friction coefficient. Previous experiments performed in the same setup (19) showed that, for constant σ, the frictional behavior is consistent with a combined formulation of rate-and-state (RS) friction (28,29) with aging evolution law, enhanced with flash-heating weakening (FH) (30)(31)(32)(33) in the form of: where f* is the friction coefficient at the reference velocity V*, a and b are RS friction parameters, L RS is the characteristic slip for the state variable evolution, V w is the weakening slip velocity, and f w is the residual friction coefficient. Note that, at steady state, L RS =θ → V. We discretize Eqs.  Fig. 3A) in experiment 1 and find that with L RS = 2 μm, the friction formulation captures the reduction of τ=σ at slip smaller than 25 μm. However, at larger slip, as σ decreases, the modeled response is significantly below the observed response because the formulation does not account for the delayed response (Fig. 4D). Prakash and Clifton (PC) constitutive law (15,16,36), the shear resistance is given by the following: The function ψ evolves exponentially with slip to the new value of σ as where L PC is another evolution distance. We discretize Eqs. 2-5 in time and search for the frictional parameters that together with the measured values V (t) and σ(t) fit the observed frictional response. Eqs. 4 and 5 require the values of ψ t for each time, which is obtained by discretizing Eq. 5 in time as follows: where the initial value of ψ t at the arrival of the rupture is ψ t = σ 0 . We use the frictional parameters as in the previous case and estimate L PC using a grid search, in which L PC ranges between 0.1 (equivalent to a negligible delay) and 10,000 μm (a large delay). The inclusion of the delayed response in the formulation, with L RS and L PC as fitting parameters, enables to obtain a much better fit between the experimental and modeled evolution of τ=σ in experiment 1 (Fig. 4D and SI Appendix, Fig. S6). A value of L PC = 1,000 μm, which is three orders of magnitude larger than the value of L RS = 2 μm, the RS evolution distance, is needed in order to capture the observed delay in shear resistance. The root mean square (RMS) of the differences between the observed and modeled responses (dashed black lines in Fig. 5) decreases significantly from value of 0.17 to a value of 0.06 between L PC = 10 and 1,000 μm, with most of the reduction between L PC = 30 and 600 μm, whereas, for L PC > 1,000 μm, there is only minor change in the RMS values. However, the modeled values of τ=σ is consistently under the observed values, suggesting that, in addition to the delayed response, the value of f w in our experiment is slightly larger than that estimated in ref. 19. Note that the values of V w = 1.1 m/s and f w = 0.27 for the weakening parameters were obtained in ref. 19 by fitting the steady-state frictional behavior of Homalite at different slip rates, using experimental data from experiments under different initial normal stresses (σ 0 ) (SI Appendix, Fig. S3).
Friction Model 3. In this model, we improve the fit with the observed response by considering a formulation of rate-and-state friction with enhanced weakening and Prakash-Clifton law, featuring weakening parameters that depend on normal stress. This formulation is consistent with high-speed friction experiments on gabbro that were performed under different normal loads and showed that V w (37) and f w (38) decrease with σ in the form of power laws: and A B C D E F Fig. 4. Experimental evidence, fitting, and subsequent prediction of pronounced delay in shear resistance response to rapid normal stress variations. (A-C) Temporal evolution of τ (red), σ (black), and V (blue) near the free surface (location marked in Fig. 3A) for experiments 1 to 3. Although experiments 1 and 2 were performed under similar loading conditions, a later supershear transition in experiment 2 leads to a more intense near-surface rupture with larger peak in V and larger reduction of σ. (D) Fitting of the measured effective friction τ/σ (experiment 1) for three models: 1) enhanced-weakening RS friction without account for a delayed response to variations in σ (blue); 2) enhanced-weakening RS friction that accounts for the delayed response by the Prakash-Clifton law (purple); and 3) enhanced-weakening rate-and-state friction with the Prakash-Clifton law and weakening parameters that depend on normal stress (red). The experimental data are best fit by friction model 3 with the Prakash-Clifton evolution distance that is two to three orders of magnitude larger than that of rate-and-state friction. (E and F) Comparison of the measured and predicted values of τ/σ for experiments 2 and 3 and friction model 3. The parameters constrained with the data in experiment 1 allow us to predict the nontrivial friction evolution in experiments 2 and 3, as well as experiments 4 to 6 (SI Appendix, Fig. S7).
where a v , b v , a f , and b f are the power laws coefficients. Each of the experiments in refs. 37 and 38 examined the steady-state frictional behavior under constant sliding velocity and normal stress. To account for rapid stress variations observed in our experiments, we modify Eqs. 7 and 8 such that V w and f w evolve gradually in response to a change in the normal stress as and We assume here that the same form of delay (ψ) affects all of the frictional properties that vary with σ; thus, we avoid using additional fitting parameters. We discretize Eqs. 2-5, 9, and 10 in time and search for a set of parameters a v , b v , a f , b f , and L PC , that together with the RS friction parameters above and the measured values V (t) and σ(t) fit the observed frictional response, similarly to the analysis carried out for the previous two models. As before, we estimate L PC using a grid search, in which L PC ranges between 0.1 (equivalent to a negligible delay) and 10,000 μm (a large delay). We do not simply perform a grid search over the parameters a v , b v , a f , and b f , but consider only combinations of the parameters that agree with the experimental data in ref. 19 (SI Appendix). This formulation does fit better the observed behavior in experiment 1, with the best fit obtained for evolution distances of L PC = 370 μm and L RS = 2 μm and power law parameters of a f = −0.34 and b f = 85 (Figs. 4D and 5). With these values, f w varies between 0.33 and 0.36 (Eq. 10), and the modeled effective friction fits the observed friction at slip of 25 to 200 μm better than model 2, although the additional improvement of the fit from model 2 to model 3 is significantly smaller than that from model 1 to model 2. Note that the differences between model 2 and model 3 should increase for the experiments performed under lower value of σ 0 (Eqs. 9 and 10). Plots of the RMS of vs. L PC for different values of a f and a v (Fig. 5 A and B, respectively) provide more insight into the differences between the models and the effect of using different model parameters. The RMS decreases from 0.18 for model 1 to 0.06 for model 2 with (L PC ≥ 1,000 μm) and further to 0.045 for model 3, with a large effect of the parameters a f , b f , and L PC (Fig. 5A) and negligible effect of the parameters a v and b v (Fig. 5B). The RMS value of 0.045 is obtained for four combinations of the parameters a v , b v , a f , and b f , with the following parameter ranges: a f = −0.346 to −0.334, b f = 72 to 95, a v = −1.55 to −0.14, and b v = 7 to 8 × 10 9 . The values of L PC that give the smallest RMS value for each of the four combination ranges between 350 and 415 μm. Now that we have identified a formulation and its parameters that capture the measurements of experiment 1, we verify its predictive value by applying it to experiments 2 to 6, without any changes in the parameters. Remarkably, we find that the formulation in model 3, together with the parameters constrained with the data in experiment 1, allows prediction of the friction evolution near the free surface in the other experiments ( Fig. 4 and SI Appendix, Fig. S7). Accounting for the normal stress dependence of the enhanced-weakening friction parameters (V w and f w through Eqs. 9 and 10, respectively) in the PC law mostly affect experiments 3 to 6 (SI Appendix, Figs. S7 and S8), which were performed under lower σ 0 than experiment 1, and consequently experienced smaller reductions in the friction coefficient. For experiments 1, 2, 3, and 5, the decreases in RMS of the differences between the observed and modeled frictional response between models 1 and 2 is significantly larger than that between models 2 and 3, while for experiments 4 and 6, the difference between the RMS obtained with models 1 and 2 is comparable to or smaller than that between models 2 and 3 (SI Appendix, Fig. S8).

Conclusions
Our study focuses on the conceptual frictional behavior of thrust faults for cases with rapid variations in normal stress. We present experimental measurements featuring rapid normal stress reduction as the ruptures interact with the free surface, with a temporary complete release for experiments under the initial compressive stresses of less than 5.7 MPa. Larger reduction in normal stress and larger slip rates for a rupture with late supershear transition compared to that with early transition suggest that the former may be more destructive in the case of thrust faults, in contrast with current assumptions that well-developed supershear ruptures are more damaging. Our findings clearly demonstrate the delay between normal stress changes and the corresponding changes in frictional resistance for the laboratory setting in an analog material, with important implications for the dynamics of thrust earthquakes near the free surface. In particular, our results indicate that the delay in shear resistance response to variations in normal stress is associated with an evolution distance that is two to three orders of magnitude larger than that of rate-and-state friction. Such delay is important in other earthquake source problems that involve rapid normal stress variations, such as seismic slip on nonplanar faults, on bimaterial faults, and in the presence of shear-heating-induced pressurization of pore fluids. Our past studies have demonstrated that similar conceptual features uncovered in our experiments are highly relevant to real faults in rocks. However, to establish the connection to real earthquakes, including the time delay, one would B A need to examine real rock materials, develop the corresponding friction relations, and then conduct a simulation of the real fault scenario. Our next steps would be directed toward conducting such experiments with real powdered rock on the interface, as found in the core of natural faults.

Materials and Methods
We study how dynamic shear ruptures on inclined frictional interface interact with the free surface through full-field measurements of displacements, velocities, strains, and stresses associated with the ruptures, obtained from six experiments (SI Appendix, Table S1). The ruptures are initiated by a local pressure release in a small region near the center of the interface, and, as they propagate upward, a target area near the free surface is monitored with an ultrahigh-speed camera to produce 128 images with a resolution of 400 × 250 pixels 2 and temporal sampling of 1 million frames per s (Fig. 1A  and SI Appendix, Fig. S1). We produce evolving in-plane displacement maps from the images using local DIC, filter the high-frequency displacement noise with a non-local-means filter (39)(40)(41), and use a postprocessing algorithm (21) to ensure continuity of tractions across the interface. We calculate the strains from the displacement fields using a finite difference approximation and compute the stress changes from the strain fields using the standard planestress linear elastic constitutive equations. The actual stresses are obtained from the stress changes by adding the initial stresses. A more detailed description of the laboratory setup and full-fields analysis is given in SI Appendix and in refs. 19 and 20.
Data Availability. The time series of slip, sliding velocity, and shear and normal stresses near the free surface are available at the Caltech Research Data Repository: https://data.caltech.edu/records/1405.