<figure><img src="image_1.png" /></figure>

# Laboratory earthquakes decipher control and stability of rupture speeds

Received: 1 November 2022

Accepted: 12 April 2023

Published online: 27 April 2023

<figure><img src="image_2.png" /></figure>

Check for updates

Peng Dong<sup>1</sup>, Kaiwen Xia<sup>1,2,3</sup>, Ying Xu<sup>3</sup>, Derek Elsworth<sup>4</sup> & Jean-Paul Ampuero<sup>5</sup>

Earthquakes are destructive natural hazards with damage capacity dictated by rupture speeds. Traditional dynamic rupture models predict that earthquake ruptures gradually accelerate to the Rayleigh wave speed with some of them further jumping to stable supershear speeds above the Eshelby speed (~√2 times S wave speed). However, the 2018 M<sub>w</sub> 7.5 Palu earthquake, among several others, significantly challenges such a viewpoint. Here we generate spontaneous shear ruptures on laboratory faults to confirm that ruptures can indeed attain steady subRayleigh or supershear propagation speeds immediately following nucleation. A self-similar analysis of dynamic rupture confirms our observation, leading to a simple model where the rupture speed is uniquely dependent on a driving load. Our results reproduce and explain a number of enigmatic field observations on earthquake speeds, including the existence of stable subEshelby supershear ruptures, early onset of supershear ruptures, and the correlation between the rupture speed and the driving load.

The severity of seismic damage is significantly affected by rupture speed ($V_r$)<fnref n="1" />, which spans from a few hectometers per second to the P wave speed ($C_p$) of the host rock²,³, except for the forbidden zone between the Rayleigh wave speed ($C_R$) and the S wave speed ($C_S$)⁴. Particularly, supershear earthquakes with rupture speeds exceeding $C_S$⁵⁻⁹ can induce increased damage due to the passage of the Mach wave<fnref n="1" />,¹⁰,¹¹. The control on rupture speed is thus crucial for understanding earthquakes and the attendant assessment of seismic hazards.

Although seismological observations and inverse analyses have provided increasingly well-resolved source parameters for natural earthquakes, key questions remain in linking observations to fundamental theories of earthquake ruptures. A particularly mystifying aspect of earthquake rupture is the mechanism for slow-supershear ruptures propagating at theoretically unstable speeds between $C_S$ and the Eshelby speed $C_E$ (~√2$C_S$)¹²,¹³, as in the 2018 M<sub>w</sub> 7.5 Palu earthquake¹⁴, the 1999 M<sub>w</sub> 7.1 Düzce earthquake¹⁵, the 1979 M<sub>w</sub> 6.5 Imperial Valley earthquake¹⁶, and the 2017 M<sub>w</sub> 7.7 Komandorski

earthquake¹⁷. To reconcile this apparent contradiction, several models have been proposed, such as the low-velocity fault zone model¹⁶,¹⁸, faults with geometric complexity¹⁹, and mixed-mode ruptures involving oblique slip²⁰. However, these models all require special conditions and are incapable of interpreting stable subRayleigh speeds.

Additionally, how a rupture along one-dimensional homogeneous faults attains a supershear speed remains an open question. Theoretical studies²¹,²² and laboratory experiments²³ have shown that such a rupture often gradually accelerates to Rayleigh speed and that some subRayleigh ruptures can jump over the forbidden zone to a supershear speed through the Burridge-Andrews mechanism - where the jump to supershear speeds is achieved through the growth of a separate rupture at the S wave front ahead of the main rupture. Alternatively, rupture can also directly accelerate through this forbidden zone into the supershear region²⁴. However, the 2018 M<sub>w</sub> 7.5 Palu earthquake attains supershear speed at the onset (or from the first moments that a velocity can be discerned)¹⁴,²⁵, challenging such conventional understanding on supershear nucleation. The controls for

<footnote n="1">¹Institute of Geosafety, School of Engineering and Technology, China University of Geosciences, Beijing 100083, China. ²Department of Civil and Mineral Engineering, University of Toronto, Toronto, ON M5S 1A4, Canada. ³State Key Laboratory of Hydraulic Engineering Simulation and Safety, School of Civil Engineering, Tianjin University, Tianjin 300072, China. ⁴Energy and Mineral Engineering & Geosciences, G3 Center and EMS Energy Institute, Pennsylvania State University, University Park, PA 16802, USA. ⁵Géoazur, Université Côte d'Azur, IRD, CNRS, Observatoire de la Côte d'Azur; 250 rue Albert Einstein, 14 Sophia Antipolis, 06560 Valbonne, France. ✉e-mail: kaiwen.xia@utoronto.ca</footnote>
rupture speeds represent another important question. Some experi-
mental studies propose that rupture velocity depends on the initial
stress state<sup>26-28</sup> and static stress drop<sup>29,30</sup>, whereas a latest research
shows that rupture velocity is controlled by fault roughness<sup>31</sup>.

Here we use tightly-constrained laboratory experiments of spon-
taneous ruptures to address these apparent contradictions between
theoretical understanding and seismological observations of earth-
quake rupture characteristics. Our experiments show that earthquake
ruptures attain stable propagation speeds spanning widely from sub-
Rayleigh to supershear regions, congruent with field observation of
natural earthquakes. The rupture speed is dictated by the driving load
that depends on the dynamic stress drop and the normal stress. This
driving load controls the strain energy available to drive the rupture
propagation. A self-similar rupture model covering speeds from sub-
Rayleigh to supershear is developed. The model prediction of rupture
speed stability ranges is consistent with both the experimental and the
field observations.

Results

Stick-slip events on laboratory faults

Experiments are conducted on a square (500 mm × 500 mm) but thin (20 mm) polymethyl methacrylate (PMMA) plate transected by an inclined laboratory fault (Fig. 1a). The roughness of the seismogenic fault is independently controlled (Supplementary Fig. 1; see “Methods” section). The PMMA plates are held together by friction against the biaxial loading (see also Methods and Supplementary Fig. 2). With this loading configuration, the shear stress applied on the fault grows with the differential stress ($\sigma_Y - \sigma_X$). With the continuous increase in the differential stress and thus the shear stress, a series of audible stick-slip

events occur on the fault and produce a measurable static stress drop. These unstable events quasi-periodically release the stress accumulated during the loading process (Fig. 1b). As the horizontal stress ($\sigma_X$) is increased in different experiments, the quasi-periodic nature of the stick-slip events is preserved (Fig. 1c, d) while the instability becomes more violent. This observation is consistent with previous studies on saw-cut laboratory faults on rock<sup>30</sup>. For the same horizontal stress ($\sigma_X$), the peak differential stress at the onset of instability increases with the roughness of the fault. Therefore, we can conclude that the intensity of the stick-slip events is jointly controlled by the stress condition and the fault topography.

For each event, we calculate the global peak shear stress ($\tau_p$) and normal stress ($\sigma_p$) at the onset of instability by resolving $\sigma_X$ and $\sigma_Y$ on the fault plane. Experimental results show that $\tau_p$ and the ratio of global peak shear stress to normal stress ($\tau_p/\sigma_p$) increase with the normal stress and the fault roughness (Fig. 1e). A rougher fault requires higher shear stress to initiate instability under the same normal stress. For all events, the static stress drop, $\Delta\tau_s$, is nearly proportional to $\tau_p$, exhibiting a consistent relative shear stress drop ($\Delta\tau_s/\tau_p$) of ~0.2 (Fig. 1f). These results are qualitatively consistent with existing experimental results<sup>32</sup>.

Spontaneous ruptures in laboratory earthquakes

In this study, ruptures are spontaneously nucleated as stress builds, similar to some existing experiments. Concurrent measurements of the spatiotemporal evolution of fault-parallel slip and the shear strain (Methods and Supplementary Fig. 3) enable the evolution of rupture propagation to be tracked (Fig. 2). These events are all nucleated at random locations away from the fault ends, indicating relatively

<figure><img src="image_3.png" /></figure>

**Fig. 1** | Stick-slip events mimicking natural earthquakes. a Schematic of laboratory earthquake model. The laboratory fault, with total length of 650 mm and observation window of 550 mm, is subjected to biaxial loading. The fault slip directions are indicated by red arrows. **b–d** Differential stress ($\sigma_Y - \sigma_X$) as a function of time at increasing horizontal stress ($\sigma_X = 1, 2, \text{and } 3 \text{ MPa}$) on extremely smooth (black), smooth (blue), and rough (red) faults. **e** Stress state, global peak shear stress $\tau_p$ as a function of peak normal stress $\sigma_p$, at the onset of instability for

events on extremely smooth (cycles), smooth (squares), and rough (triangles)
faults. The color of the symbol is coded by the horizontal stress $\sigma_X$. Dashed lines
correspond to different ratios of $\tau_p/\sigma_p$. For Static shear stress drop $\Delta\tau_s$ as a function
of peak shear stress $\tau_p$ for each event. Dashed lines indicate different values of
relative shear stress drop $\Delta\tau_s/\tau_p$. The inset illustrates the definition of static shear
stress drop $\Delta\tau_s$ and peak normal stress $\tau_p$. The symbol follows that in (e).
<figure><img src="image_4.png" /><figcaption>**Fig. 2** | Nucleation and rupture process of spontaneous events. Each subgraph shows the space-time distribution of fault parallel slip. Colors indicate the distributed measurements of fault slip obtained from optical images. The magenta lines correspond to the shear stress measured using the shear strain gauges. The</figcaption></figure>

ratio of global peak shear stress to normal stress ($\tau_p/\sigma_p$) and the measured rupture speed $V_r$ normalized by S wave speed $C_S$ (see "Methods") of each event are indicated on each subgraph.

uniform distribution of the loading as further demonstrated in Supplementary Fig. 4, and the fault strength. The distribution of shear stress along the fault may be more uniform than that in previous direct shear experiments²⁶. For events on the extremely smooth faults under low horizontal stress $\sigma_X$, the rupture initiates as a preslip zone expanding quasi-statically. This preslip zone is defined as the “nucleation zone” with its expansion representing the rupture nucleation phase. Dynamic rupture follows once the nucleation zone reaches a critical size. This observation is broadly consistent with the classic model of rupture nucleation, comprising three stages (prior to the unstable rupture), viz.: formation of the failure nucleus, quasi-static growth and acceleration³³,³⁴. With an increase in both the horizontal stress $\sigma_X$ and the fault roughness, the nucleation process becomes more complex, featuring more than one preslip zone, similar to numerical models with inhomogeneous static friction³⁵. As preslip tends to develop where the local shear stress exceeds the frictional strength, it is plausible that fault roughness similarly may result in a slightly heterogeneous distribution of the fault strength, which is initially uniform. In addition, the critical size of the nucleation zone decreases with increasing $\sigma_X$. Thus, in these experiments, the nucleation zone on the rough fault with $\sigma_X$=3 MPa is too small to be visible as limited by the spatial resolution of the measurement method.

Figure 2 also shows that the measured dynamic rupture speed (see “Methods” section) depends directly on both $\sigma_X$ and the fault roughness. Rupture on the extremely smooth faults at $\sigma_X$ = 1 MPa propagates upward at a nearly constant velocity of ~300 m/s (0.21C$_S$). In contrast, ruptures on rough faults with higher $\sigma_X$ always attain supershear speed (> C$_S$). Regardless of the speed, these ruptures attain their final propagation speed immediately following spontaneous nucleation. No intermediate phase of acceleration preceding the dynamic rupture propagation is detected in these results. Furthermore, the supershear transition by the Burridge-Andrews mechanism, predicted by numerical simulations²¹ and

confirmed by laboratory experiments with similar artificial abrupt nucleation²³, is absent. Our observations are also different from numerical results featuring direct supershear transition on sufficiently stressed faults²⁴, where the rupture speed increases continuously through the forbidden zone. Importantly, our experiments reproduce the early supershear rupture of the 2018 Palu earthquake. This spontaneous nucleation may be responsible for the absence of the speed transition proposed in the literature featuring numerical or experimental abrupt triggering. Such immediate transition to supershear has also been observed from some experiments³⁶⁻³⁸, but with limited resolution. In this study we focus on one-dimensional homogenous faults subjected to uniform loading. For one-dimensional heterogeneous faults or two-dimensional homogeneous faults, other nucleation mechanisms for supershear ruptures are possible³⁹⁻⁴¹. Furthermore, nearly steady subEshelby supershear speeds are observed (Supplementary Figs. 5 and 6), which confirms the other puzzle in the 2018 Palu earthquake. Although rupture speeds between C$_S$ and C$_E$ have been reported³⁸,⁴², the observed rupture speeds are always transient. These observations of non-steady rupture speed may be associated with the stress gradient along the fault in direct shear experiments²⁶ as first pointed out by Scholz⁴³,⁴⁴. If the rupture initiates from the edge of the direct shear specimen that suffers the most from the stress nonuniformity, then such an issue may persist even if the stress is uniform in the central part of the fault⁴⁵.

## Factors governing rupture speed

The ratio of shear stress to normal stress²⁶ and shear stress drop²⁹,⁴⁶ have been proposed as the key factors that control the rupture speed. To explore the factors governing the rupture speed, a total of 81 stick-slip events (Supplementary Data 1) are summarized in Fig. 3, enabling verifications of existing theoretical conjectures. Following Ben-David et al.²⁶, we first examine the normalized rupture velocity relative to the
<figure><img src="image_5.png" /><figcaption>**Fig. 3 | Factors controlling rupture speed.** a Dependence of rupture speed on fault roughness and the ratio of global peak shear stress to normal stress ($\tau_p/\sigma_p$), b Dependence of rupture speed on fault roughness and static stress drop $\Delta\tau_s$. The</figcaption></figure>

<figure><img src="image_6.png" /></figure>

rupture speed $V_r$ is normalized by S wave speed $C_S$. In (a) and (b), symbols represent rupture events on extremely smooth (circles), smooth (squares), and rough (triangles) faults, with color coded by horizontal stress $\sigma_X$.

ratio of $\tau_p/\sigma_p$. As shown in Fig. 3a, the rupture velocity increases with $\tau_p/\sigma_p$, which is in general consistent with their prior results. However, we observe that the specific relationship between $\tau_p/\sigma_p$ and $V_r$ also varies with $\sigma_X$. This phenomenon was either not evident or obscured within their results²⁶ where data were broadly scattered due to stress inhomogeneity. In our observations, dynamic rupture propagation speeds are steady (Fig. 2) and are absent from the theoretically predicted forbidden zone ($C_R < V_r < C_S$) (Fig. 3)—a feature not honored in existing results²⁶,²⁹,³⁰, possibly due to uncertainties in rupture speed measurements. Moreover, our results show that fault roughness plays an important role in determining dynamic rupture speed. Despite significant scatter, Passelègue et al.²⁹ proposed that the rupture speed may be related to the static stress drop following a similar argument based on field data⁴⁷. Our results also show a consistent increase in the rupture speed with the static stress drop, but the relation between these two quantities is apparently modulated by the fault roughness (Fig. 3b).

Within the framework of linear elastic fracture mechanics, rupture propagation is governed by the energy balance⁴⁸. According to this rule, a theoretical model has been derived for subRayleigh ruptures, where the rupture speed is correlated with $\Delta\tau_d^2/\sigma_N$, with $\Delta\tau_d$ and $\sigma_N$ being the dynamic stress drop and the normal stress (i.e. $\sigma_p$), respectively (see “Methods”)²⁸. This useful conclusion remains conjectural since they adopted a two-dimensional theory to model their three-dimensional triaxial tests. Here we examine our experimental data against their conjecture, with $\Delta\tau_d$ being quantified at high fidelity by the shear strain gauges near the fault (Fig. 1a, Supplementary Fig. 3 and inset of Fig. 4a). A consistent trend is obtained for all levels of fault roughness (Fig. 4a). Our experimental results not only verify this conjecture, but also extend it from the subRayleigh to both sub-Rayleigh and supershear regions.

## Self-similar analysis

Steady-state fracture models have been widely used to analyze dynamic rupture¹²,⁴⁹. In this model, for fractures with a given speed, the mechanical fields are assumed to remain invariant with respect to the reference-frame of the moving crack-tip. Because the mechanical fields (including the stress field) are invariant, the stress intensity factor remains the same. Constant stress intensity factor and constant fracture velocity result in constant fracture energy in this model. Cracks propagating in speed regimes $V_r < C_R$ and $C_S < V_r < \sqrt{2C_S}$ are shown to be unstable, while the speed regime $\sqrt{2C_S} < V_r < C_P$ is stable¹²,¹³. More generally, stable rupture speed is not accessible for simple

models with constant fracture energy²⁰,⁵⁰. Therefore, observations of ruptures growing bilaterally at roughly constant speeds (Supplementary Figs. 5 and 6) over a wider range of rupture speeds imply that the steady-state model cannot interpret our experiments. Furthermore, as the nucleation zone size is negligible in comparison to the fault length, and our one-dimensional homogenous laboratory faults are subjected to uniform loading, we can confidently treat our experiment as a mode II self-similar crack problem⁵¹,⁵², which has been partially adopted to analyze supershear ruptures⁴².

As detailed in “Methods”, the theoretical analysis is valid for both self-similar subRayleigh and supershear ruptures. The specific fracture energy $\Gamma_{II}$ can include not only the fracture energy at the rupture tip, but also other tip-localized dissipation⁵³. Experiments have revealed that specific fracture energy may be a function of the rupture speed⁵⁴ and rupture length⁵⁵, at least for mode I ruptures. We assume that the specific fracture energy is equivalent to the fracture energy of a rupture that obeys the classical slip-weakening law²¹. By further using the relation between the critical slip distance $D_c$ and the rupture length by Passelègue, et al.²⁸, we arrive at a unique relation between $\Delta\tau_d^2/\sigma_N$ and $\beta$, with the former being the generalized driving load and the latter being the rupture speed normalized by the P wave speed. This theoretical result nicely confirms the experimental dependence of the rupture speed on the driving load (Fig. 4a). We further fit the experimental data using the theoretical model (Fig. 4a), where the consistency between the model prediction and experimental data suggests that the laboratory earthquake ruptures are indeed self-similar.

These same data (Fig. 4a) show stable rupture speeds from ~0.2$C_S$ to near $C_R$ and from $C_U$ (~1.08$C_S$) up to $C_P$, consistent with field observations²⁰. However, such stable speed regions are not compatible with those predicted by the steady-state rupture model with constant fracture energy¹², where only $C_R$ in the subRayleigh region and supershear speeds above $C_E$ are deemed stable. A rupture is stable when an increase in rupture speed requires an increase in the driving load ($\Delta\tau_d^2/\sigma_N$ in the current case)¹². Following this rule, we can divide the rupture speed spectrum into several ranges (Fig. 4a). Ruptures with speeds falling in the ranges of $V_r < C_R$ and $C_U < V_r < C_P$ are all stable, while ruptures with speeds in the range $C_S < V_r < C_U$ are unstable. Our results confirm that even a simple homogeneous fault can allow a spontaneous rupture to propagate stably in the conventionally deemed unstable regions of $C_S < V_r < C_E$ and $V_r < C_R$. Unlike previous attempts to reconcile the contradiction between conventional theory and field observation by invoking geometrically sophisticated fault models¹⁶,¹⁸,²⁰, we provide an elegant and simple resolution here.
earthquakes (Methods). A compilation of the average rupture velocities and estimated stress drop of strike-slip earthquakes is available from the SRCMOD database⁶⁰ (Supplementary Data 2), together with fast ruptures on fault segments of some well-instrumented earthquakes (Supplementary Data 3). As seen in Fig. 4b, most earthquakes assume average rupture speeds slower than the Rayleigh wave speed. However, ruptures on certain fault segments can propagate at speeds faster than the S wave speed of the crust (-3.5 km/s). Supershear ruptures with speeds less than the Eshelby speed are indeed accessible for natural earthquakes. The speeds of natural earthquakes fall nicely into the stable speed ranges delineated by the self-similar earthquake model. While there are considerable uncertainties in the resolved source parameters, the rupture speed of natural earthquakes is positively correlated with the stress drop, which is generally consistent with our experimental results. We also note that some observations from natural earthquakes indicate that rupture speed may decrease with an increase in the stress drop, especially for earthquakes with subRayleigh rupture speeds⁶¹,⁶². However, the validity of this finding may suffer from tradeoffs and uncertainties in estimating source properties for natural earthquakes.

We use highly instrumented laboratory experiments to constrain magnitudes and mechanisms governing rupture speeds on one-dimensional laboratory faults with homogeneous strength and uniform tectonic loading. Our experimental results reveal an abrupt jump to supershear rupture immediately following a spontaneous nucleation stage, countering expectations of the Burridge-Andrews transition by a speed jump over the forbidden speed zone, or the direct transition through that zone by continuous acceleration. Secondly, contrary to prior characterizations, we observe that ruptures can propagate stably over nearly the full range of subEshelby speeds and that the rupture speed is directly conditioned by the magnitude of the driving load (Δτ<sub>d</sub><sup>2</sup>/σ<sub>N</sub>). These dynamic observations conform to a self-similar model of dynamic rupture that spans the full range of rupture speeds, which neatly explains the stable slow-supershear rupture of the 2018 M<sub>w</sub> 7.5 Palu earthquake and stable subRayleigh ruptures for other natural earthquakes. Thus, this model is of fundamental importance in the interpretation of the spectrum of earthquake rupture speeds, in defining critical controls of rupture speeds, in revealing nucleation mechanisms, and therefore in the assessment of seismic hazards. Because natural faults are more complex than laboratory faults, alternative models in the framework of linear elastic fracture mechanics coupled with velocity-dependent fracture energy⁶³ or non-constant fracture energy⁵³ may also work.

## Methods

### Laboratory faults

Polymethyl methacrylate (PMMA) has been commonly used as an analog to rock material in laboratory earthquake investigations²⁶,⁶⁴,⁶⁵. In this work, square PMMA plates with a dimension of 500 mm × 500 mm × 20 mm is used. PMMA has the physical and mechanical properties as follows: density ρ = 1190 kg/m³, the Young's modulus E = 6.24 GPa, the shear modulus μ = 2.40 GPa, and the Poisson's ratio ν = 0.3. The seismic wave velocities are C<sub>S</sub>=1.43 km/s and C<sub>P</sub>=2.40 km/s for plane stress conditions.

The square PMMA plate is cut diagonally into two identical triangular plates, using a computer-numerical control (CNC) engraving and milling machine. Then the cut faces are polished to remove the machining lines. Subsequently, the polished surfaces are roughened by sandblasting using white fused alumina particles for a specified duration, pressure, and blasting distance. Three kinds of grain size are utilized in blasting, resulting in planes with three kinds of roughness as quantified by a stylus surface profilometer. The measured root mean square (RMS) heights are -0.1 µm, -1 µm, and -5 µm, respectively (Supplementary Fig. 1a). We obtain the power spectral density (PSD) of the roughness by performing fast Fourier transform (Supplementary

Fig. 1b). For a self-affine surface, the PSD value p follows a power law with the wavelength λ

$$p = C\lambda^{\alpha} \quad (1)$$

where C and α are constants. The value of α is 1.15, 1.87, 2.62 for the extremely smooth, smooth, and rough fault, respectively. It can be seen that a rougher fault has a higher α value, consistent with existing studies⁵⁶.

### Experimental setup and testing procedure

Four independently controlled hydraulic rams are used to apply simulated tectonic loads to the fault model (Supplementary Fig. 2). First, the horizontal stress (σ<sub>χ</sub>) and the vertical stress (σ<sub>γ</sub>) are increased to a specified stress level at the same rate of 0.01 MPa/s. After a hold of 15 minutes, σ<sub>χ</sub> is held constant at the specified stress level while the vertical rams are switched to displacement-controlled mode at a rate of 1 µm/s. With the movement of vertical rams, σ<sub>γ</sub> increases quasi-statically and so does the differential stress (σ<sub>γ</sub>-σ<sub>χ</sub>). This biaxial loading configuration can achieve a relatively uniform normal stress and shear stress along the fault, reducing the influence of stress gradients on the nucleation process and subsequent rupture behaviors⁴³. During experiments, both σ<sub>χ</sub> and σ<sub>γ</sub> are controlled and measured by load cells to a resolution of 0.005 MPa. The displacement of the rams is controlled and measured with a resolution of 1 µm. Both stress and displacement are recorded at a sampling rate of 100 Hz.

### Diagnostic methods

We use a high-speed camera (Photron FASTCAM SAL1) to capture the image of the speckle pattern on the front side of the specimen during rupture nucleation and propagation (Supplementary Figs. 2 and 3a). Two LED lamps provide stable high intensity illumination for imaging. The camera is triggered by the acoustic emission signal by a piezoelectric acoustic sensor mounted on the specimen. The framing speed of the camera is set as 50,000 frames/s and the exposure time as 2 µs. The digital images recorded have a spatial resolution of 1024 × 112 pixels. The field of view is 600 mm in length and thus one pixel represents a physical dimension of about 600 µm. The digital image correlation (DIC) method is used to reveal the spatiotemporal evolution of the fault-parallel slip during the rupture (Supplementary Fig. 3b). This technique has been successfully applied in recent laboratory earthquake experiments⁶⁶,⁶⁷. Following an error analysis of DIC⁶⁸, we determine that the full-field displacement measurement has a standard deviation of 2 µm.

In addition to the full-field displacement measurement using the DIC method, seven shear strain gauges are mounted at 8 cm separation on the backside of the specimen along the fault to record the shear strain (Supplementary Fig. 3c). The strain gauge is 1 mm away from the linear fault. The strain signals are amplified by an amplifier (gain -200, bandwidth -1 MHz) and sampled at a rate of 1 MHz by a digital acquisition system. The digital acquisition is also triggered by the acoustic emission signal. Therefore, the strain signals are synchronized with the digital images. A schematic of the shear strain gauge is shown in Supplementary Fig. 3d. The shear strain gauge has two components, one component is in the vertical direction and the other is in the horizontal direction. Denoting the measured strain components as ε<sub>1</sub> and ε<sub>2</sub>, the shear strain (ε<sub>XY</sub>) along the fault is

$$\varepsilon_{XY} = 1/2(\varepsilon_{1} - \varepsilon_{2}) \quad (2)$$

### Determination of rupture velocity

Taking advantage of the simultaneous measurements of the full-field displacement and local shear strains, we can determine the rupture speeds by two independent methods. First, the full-field DIC
displacement measurement captures the slip distribution along the fault at different instants in time. The rupture front is identified as the point where the fault-parallel particle velocity $\dot{u}$ exceeds a threshold $\dot{u}_{th}$. We first calculate the full-field fault-parallel particle velocity $\dot{u}$ by temporal differentiation at points along the fault, at an interval of -0.6 mm (=1 pixels) (Supplementary Fig. 5). Since the error in displacement is -2 µm, the error in particle velocity is -0.1 m/s (= 2 µm / 20 µs). We set $\dot{u}_{th} = 0.2$ m/s as the threshold for the passage of rupture front, allowing the arrival time of the rupture to be determined. Because the slip velocity increases sharply as the rupture passes, the arrival time is insensitive to the small change in $\dot{u}_{th}$. Finally, we can obtain the average rupture velocity within 20 µs by temporal differentiation. Thus, the error in rupture velocity is less than 30 m/s.

Second, from the local shear strain measurement at seven locations along the fault, the arrival time of the rupture front can be picked from the shear strain record with an error of only 1 µs. The rupture velocity between neighboring strain gauge arrays is computed as the interval (80 mm) divided by the difference in arrival time. For super-shear ruptures propagating at velocities of -2 km/s, the difference in arrival time is 40 µs. Thus, the relative error in rupture velocity is 2.5% (-50 m/s); for subRayleigh ruptures propagating at a velocity of -1 km/s, the error is smaller (-25 m/s).

The along-fault rupture velocity calculated from these two independent methods is shown in Supplementary Figs. 5 and 6. These independently evaluated velocities match closely and with nearly constant velocity - as there is no sharp accelerations or decelerations. The variation in rupture velocity is small and usually within -0.1 km/s— likely resulting from small variations in applied fault stress (Supplementary Fig. 4). The average rupture velocity is computed by averaging the local rupture velocity along strike.

## Static and dynamic stress drop

In this work, static stress drop $\Delta\tau_s$ is defined as the global drop in the macroscopic shear stress resolved on the fault (inset in Fig. 1f) - similar to previous studies$^{69}$. However, static stress drop may not accurately quantify the change in shear stress associated with the rupture. This is because static stress drop includes the contributions from both the specimen and the load frame "machine stiffness" during unloading$^{37,59}$. Furthermore, static stress drop is an integration of stress change resulting from multiple ruptures during a stick-slip event. The free boundary at the end of the fault is subject to surface breakout when a rupture arrives$^{31}$, potentially allowing backward-propagating re-ruptures, as shown in Supplementary Fig. 5. Thus, the globally measured static stress drop $\Delta\tau_s$ cannot represent the stress drop accompanying the first rupture during a stick-slip event. To isolate the stress change of the first rupture from the global value, we measure the local dynamic stress drop $\Delta\tau_d^{local}$ of the first rupture from the evolution of the shear stress, as shown in the inset in Fig. 4a. This treatment is in accordance with previous studies$^{30,37}$. The dynamic stress drop $\Delta\tau_d$ is achieved by averaging the local dynamic stress drop along strike.

For natural earthquakes, because the faults are embedded in a nearly infinite medium (the crust), boundary-induced repeated ruptures are likely absent. Thus, the dynamic stress drop $\Delta\tau_d$ of laboratory earthquakes is equivalent to the "seismological" stress drop $\Delta\tau$ in the investigation of natural earthquakes, especially for ruptures on simple faults.

## Self-similar model for subRayleigh and supershear ruptures

Rupture dynamics is governed by energy balance. For mode II shear rupture with constant specific fracture energy (i.e., total energy dissipation), the rupture will accelerate towards the limiting speed $C_R$, and some of them could suddenly jump into supershear regime with the final speed approaching the other limiting speed $C_P^{20}$. However, a constant rupture speed has been observed in seismological observations$^{17}$. This constant rupture speed requires an increase in

dissipated energy with growing rupture length$^{53}$. Several mechanisms have been proposed to provide this increase, including off-fault yielding$^{53}$, provision of specific friction laws$^{70,71}$ and slip-rate-dependent fracture energy$^{63}$. Experiments have revealed that for constant rupture speed, the dissipated energy increases with rupture length, at least for mode I ruptures$^{55}$.

In a self-similar rupture model, an inherent theoretical consequence is that the specific fracture energy $\Gamma_{ll}$ must increase with the rupture length$^{72}$. $\Gamma_{ll}$ is a generalized concept that may contain energy consumed at the rupture tip (fracture energy) and off-fault inelastic dissipation$^{53}$. The proportionality between $\Gamma_{ll}$ and rupture length is consistent with dynamic rupture models with off-fault dissipation$^{50,53,72}$. Also, this proportionality is consistent with field observations that "seismological" fracture energy (breakdown work) scales with the rupture size$^{73}$.

According to the theory of dynamic fracture mechanics$^{52}$, the energy release rate for a mode II self-similar rupture is

$$G_{ll}(a, \Delta\tau_d, \beta) = G_{ll}^{s}(a, \Delta\tau_d, 0)g_{ll}(\beta) \quad (3)$$

where $\beta = V_r/C_p$ is the normalized rupture speed, $g_{ll}(\beta)$ is a universal function that describes the effect of rupture speed on the energy release, $G_{ll}^s(a, \Delta\tau_d, 0) = \frac{n(\Delta\tau_d)^2a}{4(1-k^2)\mu}$ is the energy release rate for a static mode II rupture with $\Delta\tau_d$ being the dynamic stress drop, $k$ the ratio between $C_S$ and $C_p$ for plane stress conditions, $a$ the crack half length, and $\mu$ the shear modulus.

For self-similar rupture, the specific fracture energy $\Gamma_{ll}$ increases with the rupture length$^{72}$. Another essential feature of specific fracture energy is its dependency on rupture velocity, which has been verified in experiments, at least for mode I ruptures$^{54,74}$. Such dependency is also implied in the cell model$^{72}$. Taking the above features into account, we obtain a generic expression for the specific fracture energy

$$\Gamma_{ll} \propto af(\beta) \quad (4)$$

where $f(\beta)$ is an unknown function that describes the effect of rupture speed on the fracture energy.

We further assume that specific fracture energy is equivalent to energy dissipation for a rupture that obeys the classical slip-weakening law, enabling $\Gamma_{ll}$ to be written as

$$\Gamma_{ll} = \frac{1}{2} \Delta\tau_b D_c \quad (5)$$

where $D_c$ is the critical slip distance, $\Delta\tau_b (=\tau_p - \tau_r)$ the breakdown stress drop$^{75}$, $\tau_p$ denotes the peak shear stress and $\tau_r$ the residual shear stress in the slip-weakening model. Further, $\Delta\tau_b = \tau_p - \tau_r = \sigma_N(f_p - f_r)$, where $\sigma_N$ is the normal stress(i.e. $\sigma_p$ in our experiments), and $f_p$ and $f_r$ the peak and the residual friction coefficients, respectively.

A slip-weakening model with scale-dependent $D_c = \eta a$ has been utilized for analyzing earthquake ruptures$^{53,70}$; such assumption can be considered as a lumped representation of the fracture energy, off-fault dissipation, and non-edge-localized dissipation. Thus, Eq. (5) can be expressed as

$$\Gamma_{ll} = \frac{1}{2} \sigma_N (f_p - f_r) \eta a \quad (6)$$

where $\eta$ is a dimensionless parameter that may be dependent on the rupture speed (section 4.5 in reference 53). So far, the methodology is similar to that used to model laboratory earthquake experiments in triaxial tests$^{28}$.
Comparing Eqs. (4) and (6), $\Gamma_{II}$ can be rewritten as

$$ \Gamma_{II} = \sigma_N f(\beta)a \quad (7) $$

where $f(\beta) = \frac{1}{2}(f_p - f_r)\eta$ is introduced to jointly consider the effect of rupture speed on the breakdown friction drop and $\eta$. Combining Eqs. (3) and (7), the energy balance $G_{II} = \Gamma_{II}$ leads to

$$ f(\beta) = \frac{\pi}{4(1-k^2)\mu} g_{II}(\beta) \frac{\Delta\tau_d^2}{\sigma_N} \quad (8) $$

For a subRayleigh rupture, the analytical expression for $g_{II}(\beta)$ is (equation 6.9.92 in reference 52)

$$ g_{II}(\beta) = \frac{2(1-k^2)\sqrt{(k^2 - \beta^2)^3} R(\beta)}{k\beta^2[g_2(\beta)]^2} \quad (9) $$

where $k = C_S/C_P$, $R(\beta) = 4k^3\sqrt{1-\beta^2}\sqrt{k^2-\beta^2} - (2k^2-\beta^2)^2$ is the Rayleigh function, and $g_2(\beta) = \beta^{-2}[8k^2(k^2-\beta^2)]\mathbf{E}(\sqrt{1-\beta^2}) - 4k^2(k^2-\beta^2)\mathbf{K}(\sqrt{1-\beta^2}) - \beta^{-2}[8k^4 - 8k^2\beta^2 + \beta^4]\mathbf{E}(\sqrt{1-\beta^2/k^2}) + (4k^2 - 3\beta^2)\mathbf{K}(\sqrt{1-\beta^2/k^2})$, in which $\mathbf{K}$ and $\mathbf{E}$ are complete elliptic integrals of the first and second kind.

The values of $f(\beta)$ as a function of $\beta$, evaluated from our experimental results, are shown in Supplementary Fig. 7. In particular, $f(\beta)$ approaches zero as $\beta$ is close to 0.56 (corresponding to $C_R$). We fit the experimental results by the following function in the subRayleigh region:

$$ f_1(\beta) = 3.925 \left(1 - \left(\frac{\beta}{C_R/C_P}\right)^{8.129}\right) \times 10^{-6} \quad (10) $$

This equation automatically ensures that $f_1(C_R/C_P)$ is zero. Combining Eqs. (5) and (6), we can evaluate the relationship between $\Delta\tau_d^2/\sigma_N$ and $V_r$. This function matches well our experimental data at subRayleigh speeds (Fig. 4a).

For a supershear rupture, the analytical expression for the energy release rate is more involved. In the framework of linear elastic fracture mechanics, there would be no energy flux into the crack tip for mode II cracks propagating at an intersonic (supershear) velocity, except when $V_r = C_E$, the Eshelby speed. By introducing a Barenblatt-type cohesive zone, the energy flux into the crack tip becomes positive and all supershear rupture speeds become physically acceptable. According to Broberg$^{52,76}$, $g_{II}(\beta)$ for a supershear self-similar rupture with a cohesive zone at the crack tip (equation 6.9.124 in reference 52 and equation 68 in reference 76) is

$$ g_{II}(\beta) \approx 4(1 - k^2)B(\beta)\Gamma_D(\gamma)\left(\frac{d}{a}\right)^{1-2\gamma} \quad (11) $$

where $d$ is the size of the cohesive zone, $B(\beta)$ and $\Gamma_D(\gamma)$ are known functions: $B(\beta) = \frac{\beta^2 F_1(\beta) e^{-2\gamma(1/\beta)}}{2^{2-2\gamma} \gamma^2 N^2} \frac{\sin\pi\gamma}{(1+\beta)^{2\gamma}}$, $\Gamma_D(\gamma) = 2\left[\int_0^{\frac{1}{u^{1-\gamma}}\frac{D(-u)}{du}} du\right]^{-2} w_D(\gamma)$, with

$$ \gamma(s) = \frac{1}{\pi} \tan^{-1} \frac{4k^3 \sqrt{1 - 1/s^2} \sqrt{1/s^2 - k^2}}{(1/s^2 - 2k^2)^2} $$

$$ I_0(s) = \text{p.v.} \, 2 \int_1^{1/\beta} \frac{\xi \gamma(\xi) - 1/\beta \gamma}{\xi^2 - s^2} d\xi $$

$$ F_1 = \frac{\beta^2 \sin\pi\gamma}{4k^2\sqrt{1 - \beta^2}} $$

$$ N = \int_1^{1/\beta} \frac{M\left(\frac{1}{\beta}\right) \left(\frac{2/\beta}{1/\beta - s}\right)^{\gamma} - M(s) \left(\frac{1/\beta + s}{1/\beta - s}\right)^{\gamma/\beta s}}{1/\beta - s} ds + \frac{M\left(\frac{1}{\beta}\right) \left(\frac{2}{\beta}\right)^{\gamma}}{\gamma(1/\beta - 1)^{\gamma}} $$

$$ M(s) = \frac{\sin[\pi\gamma]e^{-2I_0(s)}}{1/\beta + s} \left(\frac{s-1}{s+1}\right)^{\gamma/\beta s} $$

$$ w_D(\gamma) = \int_0^1 \xi^{1-\gamma} D(-\xi) \left\{ \int_0^1 \frac{D(-\xi) - D(-u)}{u^{1-\gamma}(u-\xi)} du + \int_1^\infty \frac{D(-\xi) du}{u^{1-\gamma}(u-\xi)} \right\} d\xi \\ D(-\xi) = 1 - \xi $$

$\gamma = \gamma(1/\beta)$ is implied when the argument is omitted.

We estimated the ratio between $d$ and $a$ based on the shear stress change during the dynamic rupture process, as shown in Supplementary Fig. 8. Generally, $d/a$ falls in the range between 0.03 and 0.05. As shown in Supplementary Fig. 9, a small change in $d/a$ does not significantly change the normalized energy release rate$^{76}$. Therefore, we assume that the ratio between $d$ and $a$ is a constant ($\approx 0.04$).

Combining Eqs. (8) and (11), we recover a relationship between $\beta$ and $f(\beta)$ for supershear ruptures, which is illustrated in the Supplementary Fig. 7. We fit the data by the following relation:

$$ f_2(\beta) = \frac{0.0009 \left[ \left( \frac{\beta}{c_S/c_P} \right)^{0.0013} - 1 \right]}{1 - \beta^{14.42}} \quad (12) $$

Note that the boundary condition that $f_2(C_S/C_P)$ is equal to zero is enforced by this equation. After substituting $f_2(\beta)$ into Eq. (8), we evaluate the relationship between $\Delta\tau_d^2/\sigma_N$ and $V_r$ numerically. The result is shown as the magenta dashed line in the supershear region in Fig. 4a.

## Estimates of rupture speed and stress drop for natural earthquakes

The data for natural earthquakes are mainly adopted from the SRCMOD database of finite-fault rupture models$^{60}$. SRCMOD is an open-source database including more than 300 finite fault rupture models. For the models included in SRCMOD, we only consider the strike-slip events that can be approximately regarded as mode II ruptures (Supplementary Data 2). Except for the stress drops ($\Delta\tau$), key source parameters, including the geometry of the fault, the seismic moment ($M_0$), and the rupture velocity ($V_r$) are provided in the finite fault rupture models. The analytical relation between the fault rupture length ($L_r$), fault rupture width ($W_r$) and the seismic moment ($M_0$)$^{77}$ allows one to estimate the stress drop ($\Delta\tau$) and to explore the correlation between $V_r$ and $\Delta\tau$. For a circular fault with radius $r$ embedded in homogeneous material, the analytical solution of the stress drop is

$$ \Delta\tau = 7M_0/16r^3 \quad (13) $$
Whereas for a rectangular fault with large aspect ratio ($L_r \gg W_r$), the stress drop becomes

$$ \Delta\tau = 2M_0/\pi L_r W_r^2 \quad (14) $$

However, the exact application range of (14) is ambiguous. Recently, Luo et al.⁷⁸ proposed a more unified model that considers the transition of scaling relation between seismic moment ($M_0$) and rupture area (A). According to this model, we can estimate the stress drop from $M_0$ and rupture dimension ($L_r$ and $W_r$) as

$$ \Delta\tau = \begin{cases} C\left(\frac{L_r}{W_r}\right) \frac{M_0}{A^{3/2}} & L_r \le W_r \\ C\left(\frac{L_r}{W_r}\right) \frac{M_0}{A W_r} & L_r > W_r \end{cases} \quad (15) $$

where $C\left(\frac{L_r}{W_r}\right) = \frac{\pi}{8} \left(7 - \frac{7-16/\pi^2}{1+\left(2W_r/L_r\right)^2}\right)$.

Noting that the rupture dimension in SRCMOD is always larger than the real rupture area on the fault. This oversized rupture dimension leads to an underestimate of stress drop when using Eq. (15). One way to improve the estimation is to trim the finite fault model to better approximate the true source dimension. Several methods have been proposed to determine the effective rupture area⁷⁹,⁸⁰. Here, we use a simple method to obtain the effective rupture area. In this method, the sub-faults with cumulative slip less than 5% of the maximum fault slip are eliminated. In addition, we assume that the aspect ratio of the finite fault model is representative of the real one. Therefore, combining the area of the trimmed fault and the original fault aspect ratio enables the rupture length ($L_r$) and rupture width ($W_r$) to be determined. Then, following Eq. (15), we can estimate stress drop for the compiled strike slip earthquakes in SRCMOD. The results for the estimated stress drop combined with other source parameters are listed in Supplementary Data 2.

Furthermore, we compile data for several well-instrumented earthquakes (Supplementary Data 3) that featured supershear ruptures. For these events, we are concerned with the rupture speed and the stress drop on the fast-rupture segments. These data are derived directly from the rupture models of these earthquakes as reported in the literature.

## Data availability

Raw data for Figs. 1, 3 and 4 are provided in the Supplementary Data 1-3. The source data for all other figures and the original laboratory data have been deposited in Figshare under accession code (https://doi.org/10.6084/m9.fishshare.18865673.v2). The natural earthquakes data is available at SRCMOD database (http://equake-rc.info/srcmod/).

## References

1. Das, S. The need to study speed. Science 317, 905-906 (2007).

2. Gomberg, J., Wech, A., Creager, K., Obara, K. & Agnew, D. Reconsidering earthquake scaling. Geophys. Res. Lett. 43, 6243-6251 (2016).

3. Kanamori, H. & Brodsky, E. E. The physics of earthquakes. Rep. Prog. Phys. 67, 1429-1496 (2004).

4. Freund, L. B. Dynamic fracture mechanics. (Cambridge University Press, 1990).

5. Archuleta, R. J. A faulting model for the 1979 Imperial Valley earthquake. J. Geophys. Res. 89, 4559-4585 (1984).

6. Bouchon, M. et al. How fast is rupture during an earthquake? New insights from the 1999 Turkey Earthquakes. Geophys. Res. Lett. 28, 2723-2726 (2001).

7. Wang, D., Mori, J. & Uchide, T. Supershear rupture on multiple faults for the Mw8.6 Off Northern Sumatra, Indonesia earthquake of April 11, 2012. Geophys. Res. Lett. 39, L21307 (2012).

8. Socquet, A., Hollingsworth, J., Pathier, E. & Bouchon, M. Evidence of supershear during the 2018 magnitude 7.5 Palu earthquake from space geodesy. Nat. Geosci. 12, 192-199 (2019).

9. Bouchon, M. & Vallee, M. Observation of long supershear rupture during the magnitude 8.1 Kunlunshan earthquake. Science 301, 824-826 (2003).

10. Andrews, D. J. Ground motion hazard from supershear rupture. Tectonophysics 493, 216-221 (2010).

11. Bernard, P. & Baumont, D. Shear Mach wave characterization for kinematic fault rupture models with constant supershear rupture velocity. Geophys. J. Int 162, 431-447 (2005).

12. Samudrala, O., Huang, Y. & Rosakis, A. Subsonic and intersonic shear rupture of weak planes with a velocity weakening cohesive zone. J. Geophys. Res. 107, 2170 (2002).

13. Burridge, R., Conn, G. & Freund, L. B. The stability of a rapid mode II shear crack with finite cohesive traction. J. Geophys. Res. 84, 2210-2222 (1979).

14. Bao, H. et al. Early and persistent supershear rupture of the 2018 magnitude 7.5 Palu earthquake. Nat. Geosci. 12, 200-205 (2019).

15. Bouin, M. P., Bouchon, M., Karabulut, H. & Aktar, M. Rupture process of the 1999 November 12 Düzce (Turkey) earthquake deduced from strong motion and Global Positioning System measurements. Geophys. J. Int 159, 207-211 (2004).

16. Huang, Y., Ampuero, J.-P. & Helmberger, D. V. The potential for supershear earthquakes in damaged fault zones - theory and observations. Earth Planet. Sci. Lett. 433, 109-115 (2016).

17. Bao, H. et al. Global frequency of oceanic and continental supershear earthquakes. Nat. Geosci. 15, 942-949 (2022).

18. Oral, E., Weng, H. & Ampuero, J. P. Does a damaged-fault zone mitigate the near-field impact of supershear earthquakes?-Application to the 2018 7.5 Palu, Indonesia, earthquake. Geophys. Res. Lett. 47, L085649 (2020).

19. Okuwaki, R., Hirano, S., Yagi, Y. & Shimizu, K. Inchworm-like source evolution through a geometrically complex fault fueled persistent supershear rupture during the 2018 Palu Indonesia earthquake. Earth Planet. Sci. Lett. 547, 116449 (2020).

20. Weng, H. & Ampuero, J.-P. Continuum of earthquake rupture speeds enabled by oblique slip. Nat. Geosci. 13, 817-821 (2020).

21. Andrews, D. J. Rupture velocity of plane strain shear cracks. J. Geophys. Res. 81, 5679-5687 (1976).

22. Burridge, R. Admissible speeds for plane-strain self-similar shear cracks with friction but lacking cohesion. Geophys. J. R. Astron. Soc. 35, 439-455 (1973).

23. Xia, K., Rosakis, A. J. & Kanamori, H. Laboratory earthquakes: The sub-Rayleigh-to-supershear rupture transition. Science 303, 1859-1861 (2004).

24. Liu, C., Bizzarri, A. & Das, S. Progression of spontaneous in-plane shear faults from sub-Rayleigh to compressional wave rupture speeds. J. Geophys. Res. Solid Earth 119, 8331-8345 (2014).

25. Voosen, P. Indonesian earthquake broke a geologic speed limit. Science 363, 570 (2019).

26. Ben-David, O., Cohen, G. & Fineberg, J. The dynamics of the onset of frictional slip. Science 330, 211-214 (2010).

27. Xu, J., Zhang, H. & Chen, X. Rupture phase diagrams for a planar fault in 3-D full-space and half-space. Geophys. J. Int. 202, 2194-2206 (2015).

28. Passelègue, F. X. et al. Initial effective stress controls the nature of earthquakes. Nat. Commun. 11, 5132 (2020).

29. Passelègue, F. X., Schubnel, A., Nielsen, S., Bhat, H. S. & Madariaga, R. From sub-Rayleigh to supershear ruptures during stick-slip experiments on crustal rocks. Science 340, 1208-1211 (2013).

30. Passelègue, F. X. et al. Dynamic rupture processes inferred from laboratory microearthquakes. J. Geophys. Res. Solid Earth 121, 4343-4365 (2016).
31. Xu, S. et al. Fault strength and rupture process controlled by fault surface topography. *Nat. Geosci.* **16**, 94–100 (2023).

32. Thompson, B. D., Young, R. P. & Lockner, D. A. Premonitory acoustic emissions and stick-slip in natural and smooth-faulted Westerly granite. *J. Geophys. Res.* **114**, B02205 (2009).

33. Latour, S., Schubnel, A., Nielsen, S., Madariaga, R. & Vinciguerra, S. Characterization of nucleation during laboratory earthquakes. *Geophys. Res. Lett.* **40**, 5064–5069 (2013).

34. Ohnaka, M. & Kuwahara, Y. Characteristic features of local break-down near a crack-tip in the transition zone from nucleation to unstable rupture during stick-slip shear failure. *Tectonophysics* **175**, 197–220 (1990).

35. Albertini, G., Karrer, S., Grigoriu, M. D. & Kammer, D. S. Stochastic properties of static friction. *J. Mech. Phys. Solids* **147**, 104242 (2021).

36. Fukuyama, E., Xu, S., Yamashita, F. & Mizoguchi, K. Cohesive zone length of metagabbro at supershear rupture velocity. *J. Seismol.* **20**, 1207–1215 (2016).

37. Xu, S. et al. Strain rate effect on fault slip and rupture evolution: Insight from meter-scale rock friction experiments. *Tectonophysics* **733**, 209–231 (2018).

38. Svetlizky, I., Albertini, G., Cohen, G., Kammer, D. S. & Fineberg, J. Dynamic fields at the tip of sub-Rayleigh and supershear frictional rupture fronts. *J. Mech. Phys. Solids* **137**, 103826 (2020).

39. Kaneko, Y., Lapusta, N. & Ampuero, J. P. Spectral element modeling of spontaneous earthquake rupture on rate and state faults: Effect of velocity-strengthening friction at shallow depths. *J. Geophys. Res.* **113**, B09317 (2008).

40. Bruhat, L., Fang, Z. & Dunham, E. M. Rupture complexity and the supershear transition on rough faults. *J. Geophys. Res. Solid Earth* **121**, 210–224 (2016).

41. Dunham, E. M., Favreau, P. & Carlson, J. M. A supershear transition mechanism for cracks. *Science* **299**, 1557–1559 (2003).

42. Kammer, D. S., Svetlizky, I., Cohen, G. & Fineberg, J. The equation of motion for supershear frictional rupture fronts. *Sci. Adv.* **4**, eaat5622 (2018).

43. Scholz, C., Molnar, P. & Johnson, T. Detailed studies of frictional sliding of granite and implications for the earthquake mechanism. *J. Geophys. Res.* **77**, 6392–6406 (1972).

44. Scholz, C. H. *The Mechanics of Earthquakes and Faulting*. 2nd edn, (Cambridge University Press, 2002).

45. Svetlizky, I. & Fineberg, J. Classical shear cracks drive the onset of dry frictional motion. *Nature* **509**, 205–208 (2014).

46. Wu, F. T., Thomson, K. C. & Kuenzler, H. Stick-slip propagation velocity and seismic source mechanism. *Bull. Seism. Soc. Am.* **62**, 1621–1628 (1972).

47. Robinson, D. P., Brough, C. & Das, S. The Mw7.8, 2001 Kunlunshan earthquake: Extreme rupture speed variability and effect of fault geometry. *J. Geophys. Res* **111**, B08303 (2006).

48. Svetlizky, I., Kammer, D. S., Bayart, E., Cohen, G. & Fineberg, J. Brittle fracture theory predicts the equation of motion of frictional rupture fronts. *Phys. Rev. Lett.* **118**, 125501 (2017).

49. Freund, L. B. The mechanics of dynamic shear crack propagation. *J. Geophys. Res.* **84**, 2199–2209 (1979).

50. Andrews, D. J. Rupture dynamics with energy loss outside the slip zone. *J. Geophys. Res.* **110**, B01307 (2005).

51. Andrews, D. J. Dynamic plane strain shear rupture with a slip-weakening friction law calculated by a boundary integral method. *Bull. Seism. Soc. Am.* **75**, 1–21 (1985).

52. Broberg, K. B. *Cracks and fracture*. (Academic Press, 1999).

53. Templeton, E. *Effects of inelastic off-fault deformation on the dynamics of earthquake rupture and branch fault activation*, Harvard University (2009).

54. Rosakis, A. J., Duffy, J. & Freund, L. B. The determination of dynamic fracture toughness of AISI 4340 steel by the shadow spot method. *J. Mech. Phys. Solids* **32**, 443–460 (1984).

55. Ravi-Chandar, K. *An experimental investigation into the mechanics of dynamic fracture*, California Institute of Technology (1982).

56. Sagy, A., Brodsky, E. E. & Axen, G. J. Evolution of fault-surface roughness with slip. *Geology* **35**, 283–286 (2007).

57. Perrin, C., Manighetti, I., Ampuero, J.-P., Cappa, F. & Gaudemer, Y. Location of largest earthquake slip and fast rupture controlled by along-strike change in fault structural maturity due to fault growth. *J. Geophys. Res. Solid Earth* **121**, 3666–3685 (2016).

58. Ohnaka, M. & Shen, L.-F. Scaling of the shear rupture process from nucleation to dynamic propagation: Implications of geometric irregularity of the rupturing surfaces. *J. Geophys. Res.* **104**, 817–844 (1999).

59. Morad, D., Sagy, A., Tal, Y. & Hatzor, Y. H. Fault roughness controls sliding instability. *Earth Planet. Sci. Lett.* **579**, 117365 (2022).

60. Mai, P. M. & Thingbaijam, K. K. S. SRCMOD: An online database of finite-fault rupture models. *Seismol. Res. Lett.* **85**, 1348–1357 (2014).

61. Chounet, A., Vallée, M., Causse, M. & Courboulex, F. Global catalog of earthquake rupture velocities shows anticorrelation between stress drop and rupture velocity. *Tectonophysics* **733**, 148–158 (2018).

62. Tan, Y. & Helmberger, D. Rupture directivity characteristics of the 2003 Big Bear sequence. *Bull. Seism. Soc. Am.* **100**, 1089–1106 (2010).

63. Weng, H. & Ampuero, J. P. Integrated rupture mechanics for slow slip events and earthquakes. *Nat. Commun.* **13**, 7327 (2022).

64. Rubinstein, S. M., Cohen, G. & Fineberg, J. Detachment fronts and the onset of dynamic friction. *Nature* **430**, 1005–1009 (2004).

65. Selvadurai, P. A. & Glaser, S. D. Laboratory-developed contact models controlling instability on frictional faults. *J. Geophys. Res. Solid Earth* **120**, 4208–4236 (2015).

66. Rubino, V., Rosakis, A. J. & Lapusta, N. Understanding dynamic friction through spontaneously evolving laboratory earthquakes. *Nat. Commun.* **8**, 15991 (2017).

67. Rubino, V., Rosakis, A. J. & Lapusta, N. Spatiotemporal properties of sub-rayleigh and supershear ruptures inferred from full-field dynamic imaging of laboratory experiments. *J. Geophys. Res. Solid Earth* **125**, e2019JB018922 (2020).

68. Rubino, V., Rosakis, A. J. & Lapusta, N. Full-field ultrahigh-speed quantification of dynamic shear ruptures using digital image correlation. *Exp. Mech.* **59**, 551–582 (2019).

69. Leeman, J. R., Saffer, D. M., Scuderi, M. M. & Marone, C. Laboratory observations of slow earthquakes and the spectrum of tectonic fault slip modes. *Nat. Commun.* **7**, 11104 (2016).

70. Andrews, D. J. Rupture models with dynamically determined breakdown displacement. *Bull. Seism. Soc. Am.* **94**, 769–775 (2004).

71. Andrews, D. J. Rupture propagation with finite stress in antiplane strain. *J. Geophys. Res.* **81**, 3575–3582 (1976).

72. Broberg, K. B. Constant velocity crack propagation-dependence on remote land. *Int. J. Solids Struct.* **39**, 6403–6410 (2002).

73. Viesca, R. C. & Garagash, D. I. Ubiquitous weakening of faults due to thermal pressurization. *Nat. Geosci.* **8**, 875–879 (2015).

74. Paxson, T. L., Lucas, R. A., Broberg, K. B., Kobayashi, A. & Rosenfield, A. R. An experimental investigation of the velocity characteristics of a fixed boundary fracture model. 415-426 (Noordhoff International Publishing, 1973).

75. Ohnaka, M. A constitutive scaling law and a unified comprehension for frictional slip failure, shear fracture of intact rock, and earthquake rupture. *J. Geophys. Res.* **108**, 2080 (2003).

76. Broberg, K. B. Intersonic bilateral slip. *Geophys. J. Int* **119**, 704–714 (1994).

77. Shearer, P. M. *Introduction to Seismology*. (Cambridge University Press, 2009).

78. Luo, Y., Ampuero, J.-P., Miyakoshi, K. & Irikura, K. Surface rupture effects on earthquake moment-area scaling relations. *Pure Appl. Geophys.* **174**, 3331–3342 (2017).
79. Hayes, G. P. The finite, kinematic rupture properties of great-sized earthquakes since 1990. *Earth Planet. Sci. Lett.* **468**, 94–100 (2017).

80. Mai, P. M. & Beroza, G. C. Source scaling properties from finite-fault-rupture models. *Bull. Seism. Soc. Am.* **90**, 604–615 (2000).

## Acknowledgements

This study was supported by National Natural Science Foundation of China (NSFC) grants no. 42174061 and 42141010 (K.X.), and Project funded by China Postdoctoral Science Foundation grant no. 2022TQ0319 (P.D.), and the Natural Sciences and Engineering Research Council of Canada (NSERC) through the Discovery grant no. 72031326 (K.X.). We thank Prof. Hiroo Kanamori from Caltech, Prof. Shamita Das from Oxford University, Prof. Emily Brodsky from UC Santa Cruz, and Dr. Yingjun Deng from Tianjin University for discussions and comments. We thank Shiqing Xu, Allan Rubin, the other anonymous reviewer, and the Editor for their constructive comments.

## Author contributions

K.X. supervised the study. P.D. performed the experiments. P.D. & Y.X. carried out the analysis. P.D. wrote the paper. K.X. edited it. K.X., P.D., D.E. and J.-P.A. contributed to the revision of the paper. All authors contributed to the data interpretation.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at

https://doi.org/10.1038/s41467-023-38137-w.

**Correspondence** and requests for materials should be addressed to Kaiwen Xia.

**Peer review information** : *Nature Communications* thanks Shiqing Xu, Allan Rubin and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023