# Continuum of earthquake rupture speeds enabled by oblique slip

Huihui Weng, Jean-Paul Ampuero

► To cite this version:

Huihui Weng, Jean-Paul Ampuero. Continuum of earthquake rupture speeds enabled by oblique slip. Nature Geoscience, 2020, 13 (12), pp.817-821. 10.1038/s41561-020-00654-4. hal-03279193

HAL Id: hal-03279193

https://hal.science/hal-03279193v1

Submitted on 23 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Continuum of earthquake rupture speeds enabled by oblique slip

Huihui Weng¹* and Jean-Paul Ampuero¹

¹Université Côte d'Azur, IRD, CNRS, Observatoire de la Côte d'Azur, Géoazur, 250 rue Albert Einstein, Sophia Antipolis, 06560 Valbonne, France

Correspondence to Huihui Weng (email: weng@geoazur.unice.fr)

Earthquake rupture speed can affect ground shaking and thus seismic hazard. Seismological observations show that large earthquakes span a continuum of rupture speeds, from slower than Rayleigh waves up to P wave speed, and include speeds that are predicted to be unstable by 2D theory. This discrepancy between observations and theory has not yet been reconciled by a quantitative model. Here we present numerical simulations that show that long ruptures with oblique slip (both strike-slip and dip-slip components) can propagate steadily at various speeds, including those previously suggested to be unstable. The obliqueness of slip and the ratio of fracture energy to static energy release rate primarily control the propagation speed of long ruptures. We find that the effects of these controls on rupture speed can be predicted by extending the 3D theory of fracture mechanics to long ruptures with oblique slip. We propose that this model provides a quantitative framework to interpret supershear earthquakes, to constrain the energy ratio of faults based on observed rupture speed and rake angle, and to relate the potential rupture speed and size of future earthquakes to the observed slip deficit along faults.
Earthquake rupture speed controls seismic energy radiation and thus affects ground shaking and seismic hazard, yet the quantitative factors controlling the rupture speed of large earthquakes are still not completely understood and the speeds of some earthquakes remain to be reconciled with basic models. In general, faster ruptures can generate stronger ground shaking, near to and far from the fault¹⁻⁴, unless competing source complexity and path effects are dominant⁵,⁶. A compilation of earthquake rupture speeds estimated from seismological observations⁷⁻¹⁴ (Fig 3a) illustrates that most earthquakes propagate at speeds slower than the shear wave speed, $v_S$, and some¹¹,¹² at speeds faster than the Eshelby speed, $v_E = \sqrt{2}v_S$ (hereafter called “fast supershear” earthquakes). Recent evidence⁷,⁸ shows that supershear earthquakes can also propagate steadily at sub-Eshelby speed (hereafter called “slow supershear”), which is unexpected from the 2D theory of fracture mechanics¹⁵. Such unexpected speeds have been reported in large earthquakes, whose ruptures are much longer than wide¹⁶,¹⁷. The propagation of ruptures with large aspect ratio has been studied theoretically in 3D in mode III, corresponding to pure-dip-slip faulting¹⁸. By extending that theory to 3D mode II ruptures, we found that slow supershear speeds are also inadmissible for long, steady, pure-strike-slip earthquakes (Methods A3), as in the 2D theory. However, natural earthquakes generally have oblique slip, with both strike-slip and dip-slip components¹⁶. A 2D computational study¹⁹ suggested that such mixed-mode ruptures can propagate at speeds between the Rayleigh wave speed $v_R$ and $v_S$, which is a “forbidden zone” for pure mode II rupture. Such speeds have been observed in 3D numerical simulations only during very short transients that would be difficult to observe in nature²⁰. Here, we show that large mixed-mode earthquakes can propagate steadily at speeds spanning the continuum of speeds observed in nature, including
"forbidden" and slow supershear speeds.

Rupture speed controlled by rake angle and energy ratio

The propagation of long mixed-mode ruptures (Fig 1) is controlled primarily by two dimensionless quantities, as deduced by dimensional analysis and confirmed by numerical simulations (Methods A1): the rake angle $\theta$ between the initial fault traction and the horizontal direction, and the energy ratio $G_c/G_0$ between dissipated and potential energies. Here, $G_c$ is the fracture energy dissipated near the rupture front and $G_0$ is the static energy release rate of mode III subshear ruptures. The latter depends on stress drop and rupture width $W$ but not on rupture length, provided rupture length is long compared to $W^{18}$. A third non-dimensional parameter, the ratio $L_c/W$ between the size of the weakening process zone at the rupture tip and the rupture width, has a secondary effect on the asymptotic rupture behavior$^{21}$. Five different rupture behaviors emerge in 3D numerical simulations as the two primary control parameters are systematically varied (Fig 2a). We first identify two large classes: self-arresting ruptures decelerate and eventually stop spontaneously, while runaway ruptures propagate unabated through the entire fault and eventually approach a steady rupture speed (Fig 2b). We further classify runaway ruptures according to their steady speed: subshear, "forbidden", slow-supershear and fast-supershear ruptures.

Remarkably, long ruptures can propagate steadily at a variety of speeds faster than the Rayleigh wave speed, even at slow supershear speeds and in the "forbidden" zone (Fig 2b). The steady speed of subshear ruptures is $v_R$ for mode II (strike-slip), $v_S$ for mode III (dip-slip), and
lies between the two for mixed mode (oblique slip)¹⁹. The steady speeds of supershear ruptures lie between $v_S$ and $v_P$ and decrease as the rake angle and energy ratio increase (Fig 2c). The same rupture behaviours are identified on the basis of the apparent horizontal speed (Methods A2), a quantity more usually constrained by seismological analyses, except that apparent horizontal speeds in the “forbidden” zone are not found (Fig S2).

The conditions separating the different rupture behaviors can be understood and quantitatively predicted by extending the theory of fracture mechanics to 3D mixed-mode long ruptures. A basic element of the theory is that the energy release rate for mixed-mode rupture is the sum of the mode II and mode III contributions²²,²³, and is of the form $G^{mix} = G_0f(\theta, v_r)$ (Methods A4). For subshear ruptures, $f(\theta, v_r)$ is independent of $v_r$. Ruptures are runaway if their energy release rate exceeds the fracture energy, $G^{mix} > G_c$, otherwise they are self-arresting. Thus the boundary between self-arresting and runaway ruptures satisfies $G_c = G^{mix}$:

$$G_c/G_0 = f(\theta, v_r) = (1 - \nu)^{-1} \cos^2 \theta + \sin^2 \theta \quad (1)$$

where $\nu = 0.25$ is Poisson's ratio, which is in good agreement with our 3D dynamic simulations results (Fig 2a). The values of $f(\theta, v_r)$ at $v_r = v_E$ and $v_r = v_S$ are determined theoretically (Methods A4) and, in combination with the steady energy balance $G_c = G^{mix}$, the boundaries between slow and fast supershear ruptures and between slow supershear and “forbidden” ruptures are well predicted (Fig 2a).

Ruptures with oblique slip can propagate steadily at slow supershear and “forbidden” speeds because their rupture fronts are not vertical but tilted (Fig 2d & S1b). A kinematic model that cap-
tures purely geometrical effects, considering an expanding elliptical front with obliquely oriented
major axis (Fig S4), qualitatively explains the occurrence of unexpected speeds on long faults but
also shows substantial discrepancies with the dynamic model (Methods A6). Fracture dynamics
theory provides a mechanical explanation for the existence of steady rupture speeds in the “forbid-
den” zone. While the mode III contribution to the energy release rate is negative in the “forbidden”
zone, in a tilted mixed-mode rupture front it is compensated by the positive mode II contribution
(Methods A4), thus enabling a steady energy balance $G^{mix} = G_c$.

## Model consistency with observed supershear rupture speeds

The theory developed here provides a new interpretive framework for supershear earthquakes that suggests a method to constrain the energy ratio $G_c/G_0$ of faults based on observations of earthquake rupture speed and rake angle. Model and observations can be compared in terms of rupture speed, rake angle and energy ratio (Fig 3a). All the supershear earthquakes observed so far have rake angles lower than 60° and a continuum of rupture speeds up to $v_P$. The basic model explains these earthquake observations and constrains the energy ratios of faults to lie between 0.5 and 0.89. For energy ratios smaller than 0.5, supershear speeds are, in theory, allowed over a wider range of rake angles (dashed curves in Fig 3a) but have not been observed in nature. A recent example of slow supershear rupture is the 2018 Mw7.5 Palu earthquake, which was inferred to propagate steadily at a sub-Eshelby speed ~4.1 km/s⁸,²⁴. Considering the rakes constrained by different studies of the Palu earthquake (~25°²⁴, ~6-15°¹⁷, and ~15-17° from USGS and gCMT), such slow supershear rupture requires an energy ratio between 0.75 and 0.85. An alternative interpretation
of the unusual speed of this earthquake assumes the presence of a low velocity fault zone²⁵, which
remains to be confirmed by local fault studies. The 2013 Mw6.7 Okhotsk deep earthquake⁹ and
the 1999 Mw7.5 Turkey Izmit earthquake¹⁰ were estimated to propagate at Eshelby speed. This
requires values of rake angle and energy ratio near the boundary between slow and fast supershear
ruptures. The rake angle of these two events are very close to mode II ruptures⁹,¹⁶. Thus, if these
ruptures have a steady Eshelby speed, their energy ratio should be around 0.89 (Fig 3a); if their
speed is not steady and the rupture comprises both super-Eshelby (fast supershear) and subshear
segments, this value is an upper bound on the energy ratio for the fault segments with super-
Ehshelby speed. An example of fast supershear rupture is the 2001 Mw8.1 Kunlun earthquake¹¹.
An intermediate portion of the rupture had super-Eshelby speed ~5 km/s and rake ~10°¹⁶, which
requires 0.7 < G<sub>c</sub>/G<sub>0</sub> < 0.8.

The model presented here also explains the continuum of earthquake rupture speeds, ranging
from slower than Rayleigh wave speeds to P wave speed (Fig 3b). For subshear runaway rup-
tures, steady propagation at speeds arbitrarily lower than the shear wave speed requires the fracture
energy to increase with rupture speed, which can result from velocity-dependent friction¹⁸. Oth-
erwise, subshear runaway ruptures accelerate to a rake-dependent steady speed between vR and
vS and, for a given rupture length, their average rupture speed increases from 0 to vS as the en-
ergy ratio decreases. In the “forbidden”, slow-supershear and fast-supershear regimes, ruptures can
propagate steadily at speeds between vR and vP, even in the absence of velocity-dependent friction:
they are stable because the velocity-dependence of energy release rate can stabilize perturbations
of rupture speed (Methods A3).
Implications for physics-based seismic hazard assessment

The fracture mechanics theory of long ruptures developed here provides a physics-based frame-
work to relate the time-dependent seismic hazard along large faults to quantities that can be con-
strained and monitored by geophysical observations, such as seismic coupling (Fig 4). A rupture
potential Φ was introduced by Weng and Ampuero¹⁸ to infer the arrest distance of long dip-slip
(modue III) ruptures with a given spatial distribution of G_c/G_0 along strike. We heuristically extend
their definition to mixed-mode long faults as:

$$
\Phi(L_1, L_2) = \int_{L_1}^{L_2} (1 - G_c / G^{mix}) dL / W \quad (2)
$$

where G^mix = G_0f(θ, v_r) is the energy release rate for subshear ruptures and W is the rupture width. The rupture potential serves to anticipate the final size of a rupture: a rupture can propagate over the entire fault segment [L_1, L_2] only if Φ(L_1, L_2) > 0, i.e., if the average of the mixed-mode energy ratio G_c/G^mix along the segment is < 1. In addition, if G_c/G^mix is much smaller than 1, such as in the slow-supershear and fast-supershear regimes in Fig. 3b, the rupture of the entire fault segment can be supershear. Therefore, two source properties that can strongly affect seismic hazard, namely rupture length and speed, can be assessed from estimates of the rake angle θ and the energy ratio G_c/G_0 along the fault. The rake angle can be estimated from geodetic data. We propose below an approach to estimate the energy ratio at each along-strike location on long faults.

On the one hand, G₀ on long faults is approximately related to final slip *D* by G₀ = CμD²/W, where C is a geometrical factor of order 1 (Methods A5). On the other hand, fracture energy G_c can be estimated from scaling relations as a function of final slip *D*. Such rela-
tions have been derived over a wide range of earthquake sizes by different approaches: dynamic earthquake modeling<sup>26-28</sup>, laboratory experiments<sup>29</sup>, and seismological methods such as kinematic source inversion<sup>30,31</sup> (Fig 4a). As a crude first-order approximation, we seek a scaling relation of the form $G_c \approx BD^n$. Scale-independent $G_c$ corresponds to $n = 0$ (e.g. linear slip-weakening friction; Fig. S5a), while scale-dependent $G_c$ corresponds to $n > 0$ (e.g., power-law slip-weakening friction; Fig. S5b). Theoretical models with off-fault inelastic dissipation<sup>2,32</sup> lead to $n = 1$ and for thermal pressurization<sup>30</sup> $n = 2/3$. As we focus here on large earthquakes, we only consider the data with $D > 0.1$ m and estimate $n$ simply based on the observations. We ignore the data of kinematic source inversions which are likely to overestimate the fracture energy due to their over-smoothing of the slip rate function<sup>27</sup>. Least squares regression gives $n = 0.7$ and $B = 3$ (the units of $G_c$ and $D$ are $MJm^{-2}$ and $m$, respectively).

The resulting relation between energy ratio and slip is: $G_c/G_0 = BWD^{n-2}/C\mu$. The spatial distribution of slip deficit rate along a fault can be inferred from geodetic observations<sup>33-35</sup>. Given an estimate of slip deficit at a future time, a worst-case scenario (largest possible magnitude) is obtained by assuming all the slip deficit is released by a single large earthquake, i.e., $D$ is set equal to the slip deficit. Because $G_0$ depends more strongly than $G_c$ on $D$ ($n < 2$), the energy ratio $G_c/G_0$ decreases with increasing slip deficit $D$. Thus the condition for runaway ruptures (equation (1)) predicts that fault segments need to accumulate a certain critical slip deficit $D^{run}(\theta)$ to become capable of hosting long runaway ruptures, otherwise they can only host self-arresting ruptures. Combining the scaling relation of energy ratio versus slip with equation (2) allows to infer an estimate of the largest possible rupture size from a slip deficit distribution. As an illus-
tration, the time-dependent evolution of the segmentation of the central Andes subduction zone in Chile predicted by the model is shown in Fig 4, and yields a reasonable estimate of return time of a 1960-like mega-earthquake of ~360 yrs (250 – 500 yrs, accounting for model uncertainties). This physics-based estimate of runaway recurrence time is consistent with the kinematic estimate of ~300 years based on slip budget³³. Similarly, a minimum slip deficit value $D^{sup}(\theta)$ is required for steady supershear ruptures (Methods A4). The model also implies that, on a given fault, supershear earthquakes should have larger slip than subshear ones, which is supported by a compilation of observations of previous supershear earthquakes showing that rupture portions that developed supershear speeds tended to host the largest slip area in each earthquake¹¹,³⁶⁻³⁸.

Previous results³⁹⁻⁴¹ suggest that the theory developed here may be extended to account for complexities with length scales smaller than the seismogenic width. The effects of low-velocity damaged zones on 3D shear ruptures can be explained by a reduced effective shear modulus³⁹, heterogeneous fracture energy can be represented by an equivalent homogeneous $G_c^{40}$, and fault roughness can be accounted for in $G_0^{41}$. These topics warrant further investigation in the context of 3D shear ruptures.

While some model inputs currently have considerable uncertainties, continued efforts to establish robust scaling relations between fracture energy and slip and to better constrain slip deficit accumulation, from synergistic developments of frictional theories, laboratory experiments and geophysical observations, should allow to integrate the concepts presented here into earthquake hazard assessment. Capitalising on advances in geodesy, seismology and rock mechanics, the pro-
posed analysis would allow to partition a fault into segments with different potential behaviors in
future earthquakes: self-arresting or runaway, subshear or supershear. By accounting for the fi-
nite width of seismogenic zones and the obliqueness of earthquake slip, our findings quantitatively
reconcile the observations of earthquake rupture speeds with the basic theory of rupture dynamics
while opening new avenues for physics-based seismic hazard assessment.
## Methods

**A1. Dynamic rupture simulations.** We set 3D dynamic rupture simulations with oblique slip on a long fault with finite seismogenic width W embedded in an unbounded, linear elastic, homogeneous medium. We use a computational domain large enough to avoid the effects of the reflected waves from the domain boundaries within the simulation time. We assume a Poisson's ratio $\nu$ of 0.25. The shear modulus and S wave speed of the medium are denoted $\mu$ and $v_S$, respectively. The P wave speed, the Eshelby speed, and the Rayleigh wave speed are $v_P = \sqrt{3}v_S$, $v_E = \sqrt{2}v_S$, and $v_R = 0.92v_S$, respectively.

We use the linear slip-weakening friction law with slip-weakening distance $d_c$, static strength $\tau_s$, and dynamic strength $\tau_d$. This is the most simple friction law adopted in computational earthquake dynamics, and allows to prescribe a constant fracture energy $G_c = 0.5d_c(\tau_s - \tau_d)$ (Fig S5a). The strength values are also fixed because the fault normal stress is constant due to the symmetries of the problem. For a pure-dip-slip fault (rake angle of 90°), Weng and Ampuero¹⁸ demonstrated that the key parameter that controls the evolution of rupture speed is the energy ratio $G_c/G_0^{III}$, where the energy release rate is $G_0^{III} = \lambda_{III}\Delta\tau^2W/\mu$ and $\Delta\tau = \tau_0 - \tau_d$ is the nominal stress drop and $\lambda_{III}$ a geometric factor of order 1. The definition of the mode II energy ratio $G_c/G_0^{II}$ is the same⁴² except for the value of the geometric factor $\lambda_{II}$. The energy ratio for purely mode II or purely mode III (assuming the same stress drop $\Delta\tau$) can be written as:

$$ 
\begin{aligned}
\frac{G_c}{G_0^{II}} &= \frac{1}{2\lambda_{II}} \frac{L_c}{W} \left[ \frac{\Delta\tau}{\tau_s - \tau_d} \right]^{-2} \\
\frac{G_c}{G_0^{III}} &= \frac{1}{2\lambda_{III}} \frac{L_c}{W} \left[ \frac{\Delta\tau}{\tau_s - \tau_d} \right]^{-2},
\end{aligned}
\quad (3) $$
where

$$L_c = \frac{\mu d_c}{\tau_s - \tau_d} \qquad (4)$$

is a characteristic frictional length proportional to the static cohesive zone size$^{43}$. The value of $\lambda_{III}$ was determined analytically and validated numerically$^{18}$: $0.96/\pi$ for a deep buried fault (infinite space, like considered here), $1.92/\pi$ for a surface-breaking fault in a half-space, and between $0.96/\pi$ and $1.92/\pi$ for a buried fault in a half-space. Here, we found numerically for mode II ruptures on a deep buried fault that $\lambda_{II} \approx 0.96/\pi/(1 - \nu)$, which is similar to the value 0.43 obtained by Weng and Yang$^{42}$. Then we have $\lambda_{II}/\lambda_{III} = (1 - \nu)^{-1}$. In the main text, we denote $G_0 = G_0^{III}$, and thus $G_0^{II} = (1 - \nu)^{-1}G_0$. To prescribe the energy ratio $G_c/G_0$, we fix the value of the cohesive ratio $L_c/W = 0.25$ and vary the stress ratio $\Delta\tau/(\tau_s - \tau_d)$. Note that here we denote $\Delta\tau$ the absolute amplitude of stress drop. The minimum value of the energy ratio is proportional to the cohesive ratio, $G_c/G_0 \propto L_c/W$, and is obtained when the stress drop $\Delta\tau$ equals the strength drop $\tau_s - \tau_d$ (in such extreme case, the P wave from the hypocenter can trigger the rupture of the entire fault, enabling rupture at the P wave speed for all mixed-mode ruptures). Since we consider oblique slip with rake angle $\theta$ (the direction between the initial traction vector and the horizontal direction), the initial shear stress, whose amplitude is $\tau_0$, has an along-strike component $\tau_0 \cos \theta$ and along-dip component $\tau_0 \sin \theta$. Exploiting the symmetries of the problem, we only need to simulate rake angles between 0° and 90°. Other values $\theta'$ between -180° and 180° can be mapped to the 0-90° range as $\theta = \min(|\theta'|, 180 - |\theta'|)$. If the absolute initial stress $\tau_0$ is too small compared to the stress drop $\Delta\tau$, the slip direction may be time-dependent inside the cohesive zone$^{19}$ and thus the actual fracture energy may be larger than $G_c$. To have full control on the actual value of the
fracture energy, we set up a relatively large initial stress, $\tau_0/\Delta\tau \approx 10$.

We prescribe a time-dependent weakening over the nucleation zone of size $L/W = 2$ to nucleate unilateral ruptures at prescribed speeds. Rupture propagation becomes spontaneous outside the nucleation zone. To study steady supershear ruptures, without focusing on the supershear transition, we set the nucleation speed as $1.1v_S$ or $1.414v_S$. Tests show that the value of the nucleation speed does not affect the steady-state supershear speed (Fig S3). To study self-arresting and runaway ruptures, we use a sub-Rayleigh nucleation speed of $0.5v_S$.

We use the spectral element software SPECFEM3D$^{44-47}$ for the dynamic simulation. All the simulations are conducted on a medium-scale computing cluster with 64 cores and 384 GB memory. We set the time step based on the Courant-Friedrichs-Lewy stability condition. To guarantee sufficient numerical resolution, we set a grid size much smaller than the characteristic frictional length, i.e., $L_c/\Delta x = 10$. We also test a few models with refined grid, $L_c/\Delta x = 20$, and find their results are the same.

**A2. Calculations of rupture speed** We compute two types of rupture speed: depth-averaged real speed $v_r$ and apparent horizontal speed $v_r^{hor}$ (Fig S1a). The real speed is computed at each point on the fault from the gradient of rupture time $t(x_1, x_3)$

$$v_r^{real}(x_1, x_3) = \frac{1}{\sqrt{(\partial t/\partial x_1)^2 + (\partial t/\partial x_3)^2}} \quad (5)$$
where $x_1$ and $x_3$ are the along-strike and along-dip distances, respectively. The horizontal apparent speed is based on the horizontal gradient of rupture time:

$$v_r^{hor}(x_1, x_3) = \frac{1}{\partial t/\partial x_1} \qquad (6)$$

We averaged the real speed and apparent horizontal speed along depth at each along-strike position.

**A3. Energy release rate for mode II supershear rupture** The 2D theory predicts the energy release rate of supershear ruptures has the following form$^{48}$:

$$G = g(v_r) \frac{\Delta\tau^2 L}{\mu} \left(\frac{\Lambda}{L}\right)^{p(v_r)} \qquad (7)$$

where $g(v_r)$ and $p(v_r)$ are known functions of rupture speed, $L$ is the rupture propagation distance, and $\Lambda$ is the size of the dynamic cohesive zone, $\Lambda \propto L_c$. In general, $g(v_r)$ depends on the shape of the slip-weakening curve$^{48}$, but in this study the friction law is fixed. In 2D, $G$ increases from 0 at $v_r = v_S$ to its peak value at $v_r = v_E$, and then decreases to 0 at $v_r = v_P$. As $p(v_r) < 1$ for all speeds between $v_S$ and $v_P$, $G \propto L^{1-p(v_r)}$ is a monotonously increasing function of $L$. Hence, for a constant fracture energy $G_c$, the rupture speed $v_r$ approaches the P wave speed as $L$ grows. Only if the fracture energy is scale-dependent in the form $G_c \propto L^{1-p(v_r)}$ can steady supershear ruptures exist. Otherwise the only admissible steady speed is the P wave speed. For elongated ruptures in 3D, the theory by Weng and Ampuero$^{18}$ predicts that $G$ saturates when the rupture reaches a finite width $W$; it becomes a function of $W$ instead of $L$. Here, we make heuristic modifications to equation (7) by replacing $L$ with $W$:

$$G = g(v_r) \left(\frac{\Lambda}{W}\right)^{p(v_r)} G_0. \qquad (8)$$
Here $g(v_r)$ differs from the one in the 2D theory by a geometrical factor of order 1. The energy balance $G = G_c$ gives

$$ \frac{G_c}{G_0} = g(v_r) \left(\frac{\Lambda}{W}\right)^{p(v_r)} \qquad (9) $$

We suppose that, like in the 2D case, the right side of the above equation also increases from 0 at $v_r = v_S$ to its peak at $v_r = v_E$, and then decreases to 0 at $v_r = v_P$. This equation of motion of mode II long ruptures predicts that supershear propagation is stable if the energy ratio is below the maximum of the right side of equation 9, which is numerically estimated as $g(v_E) \approx 0.9$ (note that$^{48}$ $p(v_E) = 0$). If $G_c/G_0 < 0.9$, there are two mathematical solutions of this equation of motion, one with speed between $v_S$ and $v_E$ and the other between $v_E$ and $v_P$. Only the latter is stable, because the velocity-decreasing energy release rate provides a negative feedback to any perturbation of rupture speed, which stabilizes steady ruptures. In our 3D purely mode II dynamic simulations, we only observe steady supershear ruptures at speeds between $v_E$ and $v_P$ (Fig 2a), which is well explained by the heuristic equation of motion.

**A4. Energy release rate for mixed-mode rupture** For mixed-mode ruptures in 3D faults with finite width $W$, we use a reduced-dimensionality (2.5D) model to derive the energy release rate. The 2.5D model has been previously shown to be a very good approximation of the 3D elongated rupture model for pure dip-slip faulting$^{18}$. It assumes that the rupture front is nearly vertical. This assumption is reasonable because in the 3D dynamic simulations, the angles of mixed-mode rupture fronts are quite small (<10°) for fast supershear, sub-shear runaway and self-arresting ruptures (Fig 2d). For slow supershear and “forbidden” speeds, the rupture front tilt is substantial and its effects cannot be ignored, as explained below.
The energy release rate is the rate of mechanical energy flow into the rupture tip per unit rupture advance. The stress drop vector (fault-parallel traction change) is approximately parallel to the slip vector, because we focus on situations with little rake rotation. The total energy release rate for a mixed-mode rupture is the sum of the mode II and III contributions<sup>22,23</sup>, which are associated to the along-strike $\Delta\tau_{str} = \Delta\tau \cos\theta$ and along-dip $\Delta\tau_{dip} = \Delta\tau \sin\theta$ components of stress drop, respectively:

$$G^{mix} = G^{II} \cos^2 \theta + G^{III} \sin^2 \theta \qquad (10)$$

where $G^{II}$ and $G^{III}$ denote the energy release rates of purely mode II and III ruptures, respectively, that would prevail if both modes had the same stress drop $\Delta\tau$. Equation (10) stems from the energy release rate of each mode being proportional to the square of the associated stress intensity factor<sup>22,23</sup> and from the principle of linear superposition of stress intensity factors which, for a straight crack, results in $K_{II}$ depending only on $\Delta\tau_{str}$ and $K_{III}$ only on $\Delta\tau_{dip}$. The equation also applies to some curved front situations. The stress intensity factors at any point along a static circular rupture front of radius $a$ are<sup>49</sup>

$$K_{II} \propto \Delta\tau \sqrt{a} \cos \omega; \quad K_{III} \propto \Delta\tau \sqrt{a} \sin \omega \qquad (11)$$

where $\omega$ is the angle between the slip direction and the local rupture propagation direction. The expressions have a similar form at the major axis tip of an elliptical rupture, which can be set horizontal and $a = W$ to draw an analogy to the 2.5D model.
Based on 2.5D models (Methods A1 and A3):

$$G^{II} = \begin{cases} (1 - \nu)^{-1} G_0, & \text{if } v_r < v_R \\ G_{fb}^{II}, & \text{if } v_R < v_r < v_S \\ g(v_r)(\frac{\Delta}{W})^{p(v_r)} G_0, & \text{if } v_S < v_r < v_P \end{cases} \quad (12)$$

and

$$G^{III} = \begin{cases} G_0, & \text{if } v_r < v_S \\ G_{Sup}^{III}, & \text{if } v_r > v_S \end{cases} \quad (13)$$

where $G_0 = \lambda_{III} \Delta \tau^2 W / \mu$ and $v_r$ is the depth-averaged real speed. The 2D analytical solutions of $G_{fb}^{II}$ and $G_{Sup}^{III}$ depend on different mathematical assumptions: inverse-square-root singularity or no singularity of stress at the rupture tip$^{15,19,50}$. In the 2.5D model, under the same mathematical assumptions, the singular solution yields the same expressions as for sub-Rayleigh ruptures but with negative values, $G_{fb}^{II} = -(1 - \nu)^{-1} G_0$ and $G_{Sup}^{III} = -G_0$, and the non-singular solution $G_{fb}^{II} = 0$ and $G_{Sup}^{III} = 0$. Within the ranges bracketed by these two end-member analytical solutions, the values of $G_{fb}^{II}$ and $G_{Sup}^{III}$ appropriate to describe our simulations are determined empirically as explained below.

The mixed-mode energy release rate is of the form $G^{mix} = G_0 f(\theta, v_r)$. Self-arresting ruptures occur if the energy release rate of mixed-mode steady ruptures is too small to match the fracture energy, $G_c > G^{mix}$. Runaway ruptures near the boundary with self-arresting ruptures have sub-Rayleigh speeds and almost vertical fronts (< 5°; Fig 2d). Thus the theoretical boundary between self-arresting and runaway ruptures corresponds to the condition $G_c = G^{mix}$ evaluated at
sub-Rayleigh speeds ($v_r < v_R$):

$$f(\theta, v_r) = G_c/G_0 = (1 - \nu)^{-1} \cos^2 \theta + \sin^2 \theta. \quad (14)$$

Note that at sub-Rayleigh speeds $f(\theta, v_r)$ does not depend on $v_r$. The theoretical boundary between fast supershear and slow supershear ruptures is obtained by evaluating the energy balance $G_c = G^{mix}$ at $v_r = v_E$:

$$G_c = g(v_E)G_0 \cos^2 \theta + G_{Sup}^{III} \sin^2 \theta \quad (15)$$

where $(\frac{\Lambda}{W})^{p(v_E)} = 1$ because $p(v_E) = 0$ and we know that $g(v_E) = 0.9$ (Methods A3). We find that if we set an empirical value of $G_{Sup}^{III}/G_0 = -0.4$ (amid the two end-member analytical solutions) in the theoretical equation (15) the resulting equation fits well the Eshelby boundary from our 3D dynamic simulations (Fig 2a):

$$f(\theta, v_E) = G_c/G_0 = 0.9 \cos^2 \theta - 0.4 \sin^2 \theta. \quad (16)$$

For the boundary between supershear and subshear regimes, the theoretical relation is

$$G_c = g(v_{S+}) (\frac{\Lambda}{W})^{p(v_{S+})} G_0 \cos^2 \theta + G_{Sup}^{III} \sin^2 \theta \quad (17)$$

where $v_{S+}$ is a rupture speed slightly larger than the S wave speed. Near this boundary, the rupture front tilts severely, ~30° (Fig 2d), thus the effects of tilted rupture front needs to be considered. The first term on the right side of equation (17) is positive and the second term is non positive. The first term on the right side of equation (17) need to be sufficient to support the dissipated fracture energy. As the term $g(v_{S+})$ is quite small according to the theory (Methods A3), we suppose that the geometrical effect of tilted front enlarges the size of the “apparent cohesive zone” along the
strike direction to make the first term sufficiently large. We find that if we set empirical values of $g(v_{S+})(\frac{\Lambda}{W})^{p(v_{S+})} = 0.9$ (same as equation (16)) and $G_{Sup}^{III}/G_0 = 0.0$ (one end-member analytical solution) in the theoretical equation (17) the resulting equation fits well the supershear boundary for small rake angle (Fig 2a)

$$f(\theta, v_S) = G_c/G_0 = 0.9 \cos^2 \theta. \qquad (18)$$

For larger rake angle, equation (18) underestimates the energy release rate due to the even larger tilt of the rupture front (>30°, Fig 2d).

**A5. Scaling relation of energy release rate** On long faults, the depth-averaged static energy release $G_0(x)$ is related to final slip $D(x, z)^{51}$ by

$$G_0(x) = \frac{1}{2} \int_0^W \Delta\tau(x, z) D(x, z) dz/W \qquad (19)$$

where $x$ and $z$ are along-strike and along-dip distances, respectively. To first order, $G_0(x) = \frac{1}{2}\Delta\tau(x)D(x)$ (Fig. S5) and $\Delta\tau(x) = 2C\mu D(x)/W$, thus this equation is approximated as

$$G_0(x) = \frac{C\mu D(x)^2}{W} = \frac{1}{4C} \frac{\Delta\tau(x)^2 W}{\mu} \qquad (20)$$

where $C$ is a geometrical factor of order 1 and $\Delta\tau(x)$ and $D(x)$ are the depth-averaged stress drop and slip, respectively. We retrieve $G_0(x) = G_0^{III}(x)$ by taking $C = 1/(4\lambda_{III})$ for mode III and $G_0(x) = G_0^{II}(x)$ by taking $C = 1/(4\lambda_{II})$ for mode II (Method A1), respectively. These two values of $C$ agree within 4% with previous analytical solutions of static energies in full space, i.e., $C = 2/\pi/(1 - \nu)$ for mode III$^{52}$ and $C = 2/\pi$ for mode II$^{53,54}$, respectively. For a mixed-mode
rupture,

$$ 
\begin{aligned}
\Delta\tau^2 &= \Delta\tau_{str}^2 + \Delta\tau_{dip}^2 \\
\Delta\tau_{str} &= \frac{(1 - \nu)\pi\mu}{4W} D_{str} \\
\Delta\tau_{dip} &= \frac{\pi\mu}{4W} D_{dip}.
\end{aligned}
\qquad (21) $$

Since we focus on situations with little rake rotation, we have

$$ 
\begin{aligned}
D_{str} &= D \cos \theta \\
D_{dip} &= D \sin \theta.
\end{aligned}
\qquad (22) $$

Therefore, the factor relating stress drop $\Delta\tau$ and final average slip $D$ for mixed-mode rupture is

$$ C = \frac{\pi}{4} \sqrt{(1 - \nu)^2 \cos^2 \theta + \sin^2 \theta} \qquad (23) $$

For a half-space with rupture breaking the surface, $C$ is divided by two¹⁸. For a buried fault, the value of $C$ lies between the values of the full-space and half-space models. Our conceptual conclusions hold qualitatively for all values of $C$. The quantitative effect of $C$ on our results is encapsulated in equation 20; in particular, it affects the value of the critical slip deficits for runaway and supershear ruptures, discussed in the main text.

**A6. Kinematic model** To investigate the possible reasons of unexpected steady speeds, we compare the dynamic models with a simple kinematic model designed to capture purely-geometric effects (Fig S4a). We assume that a supershear rupture extends as an elliptical front propagating at the P wave speed along its major axis and at the S wave speed along its minor axis, the limiting speeds for mode II and III ruptures, respectively. The elliptical fronts are truncated to lie inside the seismogenic portion of the fault. The rake angle is the angle between the major axis of the ellipse
and the strike direction. We vary the rake angle and compute the depth-averaged real speed, the
horizontal speed and the depth-averaged rupture propagation angle (angle of the real rupture speed
relative to the horizontal direction).

We find that the basic geometrical effects of tilted elliptical front represented in the kinematic
model only account for part of the dynamic simulation results. The rupture speeds decrease as the
rake angle increases (Fig S4b and S4c), as in the dynamic models, but there are also important
discrepancies between the two models. Beyond a rake of 20-30°, the speed of the dynamic models
decreases faster than that of the kinematic models. An eventual drop to sub-Rayleigh speeds is
only found in the dynamic models. Furthermore, the variability of the real speed across the depth
profile is larger in the dynamic model than in the kinematic model. The dependency of the rupture
propagation angle as a function of rake angle (Fig S4d) is totally different between the two models.
Also, the variability across depth of the rupture angle is much larger in the dynamic model than
in the kinematic model, which means the curvature of the dynamic front is larger than that of the
kinematic elliptical front. The rupture angle drops to less than 10° once the speed drops below the
Rayleigh speed, only in the dynamic models.

1. Scala, A., Festa, G. & Del Gaudio, S. Relation between near-fault ground motion impulsive signals and source parameters. *Journal of Geophysical Research: Solid Earth* **123**, 7707–7721 (2018).

2. Gabriel, A., Ampuero, J., Dalguer, L. & Mai, P. M. Source properties of dynamic rupture pulses with off-fault plasticity. *Journal of Geophysical Research* **118**, 4117–4126 (2013).
3. Dunham, E. M. & Bhat, H. S. Attenuation of radiated ground motion and stresses from three-dimensional supershear ruptures. *Journal of Geophysical Research: Solid Earth* **113** (2008).

4. Andrews, D. Ground motion hazard from supershear rupture. *Tectonophysics* **493**, 216–221 (2010).

5. Bizzarri, A., Dunham, E. M. & Spudich, P. Coherence of mach fronts during heterogeneous supershear earthquake rupture propagation: Simulations and comparison with observations. *Journal of Geophysical Research* **115** (2010).

6. Vyas, J. C., Mai, P. M., Galis, M., Dunham, E. M. & Imperatori, W. Mach wave properties in the presence of source and medium heterogeneity. *Geophysical Journal International* **214**, 2035–2052 (2018).

7. 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).

8. Bao, H. *et al*. Early and persistent supershear rupture of the 2018 magnitude 7.5 palu earthquake. *Nature Geoscience* **12**, 200–205 (2019).

9. Zhan, Z., Helmberger, D. V., Kanamori, H. & Shearer, P. M. Supershear rupture in a Mw 6.7 aftershock of the 2013 sea of okhotsk earthquake. *Science* **345**, 204–207 (2014).

10. Bouchon, M. *et al*. How fast is rupture during an earthquake? new insights from the 1999 turkey earthquakes. *Geophysical Research Letters* **28**, 2723–2726 (2001).
11. Bouchon, M. & Vallée, M. Observation of long supershear rupture during the magnitude 8.1 kunlunshan earthquake. *Science* **301**, 824–826 (2003).

12. Yue, H. et al. Supershear rupture of the 5 january 2013 craig, alaska (Mw 7.5) earthquake. *Journal of Geophysical Research* **118**, 5903–5919 (2013).

13. Hicks, S. et al. Back-propagating super-shear rupture in the 2016 M7.1 Romanche transform fault earthquake. *Accepted for publication in Nature Geoscience* with doi: 10.1038/s41561-020-0619-9 Preprint at https://eartharxiv.org/bs7na/ (2019).

14. Wang, D. & Mori, J. Short-period energy of the 25 april 2015 Mw 7.8 nepal earthquake determined from backprojection using four arrays in europe, china, japan, and australia. *Bulletin of the Seismological Society of America* **106**, 259–266 (2016).

15. Burridge, R. Admissible speeds for plane-strain self-similar shear cracks with friction but lacking cohesion. *Geophysical Journal International* **35**, 439–455 (1973).

16. Mai, P. M. & Thingbaijam, K. Srcmod: An online database of finite-fault rupture models. *Seismological Research Letters* **85**, 1348–1357 (2014).

17. Song, X. et al. Geodetic observations of the 2018 mw 7.5 sulawesi earthquake and its implications for the kinematics of the palu fault. *Geophysical Research Letters* **46**, 4212–4220 (2019).

18. Weng, H. & Ampuero, J. The dynamics of elongated earthquake ruptures. *Journal of Geophysical Research: Solid Earth* **124**, 8584–8610 (2019).
19. Andrews, D. Dynamic growth of mixed-mode shear cracks. *Bulletin of the Seismological Society of America* **84**, 1184–1198 (1994).

20. Bizzarri, A. & Das, S. Mechanics of 3-d shear cracks between rayleigh and shear wave rupture speeds. *Earth and Planetary Science Letters* **357**, 397–404 (2012).

21. Ampuero, J. & Mao, X. Upper limit on damage zone thickness controlled by seismogenic depth. *Fault Zone Dynamic Processes: Evolution of Fault Properties During Seismic Rupture* **227**, 243 (2017).

22. Freund, L. *Dynamic fracture mechanics* (Cambridge university press, 1998).

23. Irwin, G. R. Analysis of stresses and strains near the end of a crack transversing a plate. *Trans. ASME, Ser. E, J. Appl. Mech.* **24**, 361–364 (1957).

24. Fang, J. et al. The 2018 mw 7.5 palu earthquake: A supershear rupture event constrained by insar and broadband regional seismograms. *Remote Sensing* **11**, 1330 (2019).

25. Oral, E., Weng, H. & Ampuero, J. P. Does a damaged-fault zone mitigate the near-field impact of supershear earthquakes? -application to the 2018 Mw 7.5 palu, Indonesia earthquake. *Geophysical Research Letters* **47**, e2019GL085649 (2020).

26. Mai, P. et al. On scaling of fracture energy and stress drop in dynamic rupture models: Consequences for near-source ground-motions. *Earthquakes: Radiated Energy and the Physics of Faulting* 283–293 (2006).
27. Weng, H. & Yang, H. Constraining frictional properties on fault by dynamic rupture simulations and near-field observations. *Journal of Geophysical Research* **123**, 6658-6670 (2018).

28. Gallovič, F., Valentová, L., Ampuero, J. & Gabriel, A. Bayesian dynamic finite-fault inversion: 1. method and synthetic test. *Journal of Geophysical Research: Solid Earth* **124**, 6949-6969 (2019).

29. Nielsen, S. et al. G: Fracture energy, friction and dissipation in earthquakes. *Journal of Seismology* **20**, 1187–1205 (2016).

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

31. Tinti, E., Cocco, M., Fukuyama, E. & Piatanesi, A. Dependence of slip weakening distance (dc) on final slip during dynamic rupture of earthquakes. *Geophysical Journal International* **177**, 1205–1220 (2009).

32. Andrews, D. Rupture dynamics with energy loss outside the slip zone. *Journal of Geophysical Research* **110** (2005).

33. Moreno, M. et al. Heterogeneous plate locking in the South-Central Chile subduction zone: Building up the next great earthquake. *Earth and Planetary Science Letters* **305**, 413–424 (2011).

34. Jolivet, R., Simons, M., Agram, P., Duputel, Z. & Shen, Z. Aseismic slip and seismogenic coupling along the central san andreas fault. *Geophysical Research Letters* **42**, 297–306 (2015).
35. Villegas-Lanza, Juan Carlos and Chlieh, Mohamed and Cavalié, O and Tavera, Hernando and Baby, Patrice and Chire-Chira, J and Nocquet, J-M et al. Active tectonics of peru: Heterogeneous interseismic coupling along the nazca megathrust, rigid motion of the peruvian sliver, and subandean shortening accommodation. *Journal of Geophysical Research: Solid Earth* **121**, 7371–7394 (2016).

36. Pondard, N., Armijo, R., King, G. C. P., Meyer, B. & Flerit, F. Fault interactions in the sea of marmara pull-apart (north anatolian fault): earthquake clustering and propagating earthquake sequences. *Geophysical Journal International* **171**, 1185–1197 (2007).

37. Lin, A. et al. Co-seismic strike-slip and rupture length produced by the 2001 ms 8.1 central kunlun earthquake. *Science* **296**, 2015–2017 (2002).

38. Perrin, C., Manighetti, I., Ampuero, J., 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. *Journal of Geophysical Research: Solid Earth* **121**, 3666–3685 (2016).

39. Weng, H., Yang, H., Zhang, Z. & Chen, X. Earthquake rupture extents and coseismic slips promoted by damaged fault zones. *Journal of Geophysical Research* **121**, 4446–4457 (2016).

40. Roux, S., Vandembroucq, D. & Hild, F. Effective toughness of heterogeneous brittle materials. *European Journal of Mechanics-A/Solids* **22**, 743–749 (2003).

41. Vandembroucq, D. & Roux, S. Mode iii stress intensity factor ahead of a rough crack. *Journal of the Mechanics and Physics of Solids* **45**, 853–872 (1997).
42. Weng, H. & Yang, H. Seismogenic width controls aspect ratios of earthquake ruptures. *Geo-physical Research Letters* **44**, 2725–2732 (2017).

43. Day, S., Dalguer, L., Lapusta, N. & Liu, Y. Comparison of finite difference and boundary integral solutions to three-dimensional spontaneous rupture. *Journal of Geophysical Research* **110**, B12307 (2005).

44. Ampuero, J. Etude physique et numérique de la nucléation des séismes. *PhD Thesis, University of Paris VII, France* (2002).

45. Komatitsch, D. & Vilotte, J.-P. The spectral element method: an efficient tool to simulate the seismic response of 2d and 3d geological structures. *Bulletin of the seismological society of America* **88**, 368–392 (1998).

46. Galvez, P., Ampuero, J.-P., Dalguer, L. A., Somala, S. N. & Nissen-Meyer, T. Dynamic earthquake rupture modelled with an unstructured 3-d spectral element method applied to the 2011 m 9 tohoku earthquake. *Geophysical Journal International* **198**, 1222–1240 (2014).

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

48. Broberg, K. Intersonic bilateral slip. *Geophysical Journal International* **119**, 706–714 (1994).

49. Tada, H., Paris, P. C., Irwin, G. R. & Tada, H. *The stress analysis of cracks handbook*, vol. 130 (ASME press New York, 2000).
50. Freund, L. & Clifton, R. On the uniqueness of plane elastodynamic solutions for running cracks. *Journal of elasticity* **4**, 293–299 (1974).

51. Kanamori, H. & Rivera, L. Energy partitioning during an earthquake. *Earthquakes: Radiated Energy and the Physics of Faulting* 3–13 (2006).

52. Starr, A. Slip in a crystal and rupture in a solid due to shear. In *Mathematical Proceedings of the Cambridge Philosophical Society*, vol. 24, 489–500 (Cambridge University Press, 1928).

53. Knopoff, L. Energy release in earthquakes. *Geophysical Journal International* **1**, 44–52 (1958).

54. Haskell, N. Total energy and energy spectral density of elastic wave radiation from propagating faults. *Bulletin of the Seismological Society of America* **54**, 1811–1841 (1964).

55. Dziewonski, A. M. & Anderson, D. L. Preliminary reference earth model. *Physics of the earth and planetary interiors* **25**, 297–356 (1981).

**Correspondence** Correspondence and requests for materials should be addressed to Huihui Weng (email: weng@geoazur.unice.fr).

**Acknowledgements** The open-source software SPECFEM3D used in our 3-D dynamic rupture simulations is available from the Computational Infrastructure for Geodynamics at https://geodynamics.org/cig/software/specfem3d/. This work was supported by the French government through the Investments in the Future project UCAJEDI (ANR-15-IDEX-01) managed by the French National Research Agency (ANR). We thank Diego Molina for providing the central Andes coupling model in
digital form.

**Author contributions** H.W. designed and carried out the numerical experiments, and analysed the numerical results. H.W. and J.-P.A. developed the theoretical model, interpreted the results, and wrote the paper.

**Competing Interests** The authors declare that they have no competing interests.
**Figure 1 Earthquake rupture propagation on a long fault with oblique slip.** The fault has a finite width *W* and slips at an oblique angle (rake) relative to the strike direction. The inset shows the propagation of a tilted rupture front on a fault; arrows show the local direction of rupture propagation (real speed vector).

**Figure 2 Rupture propagation controlled by energy ratio and rake angle.** (a) Five different rupture behaviours (see legend) found in 3D dynamic rupture simulations depending on energy ratio *G*<sub>c</sub>/*G*<sub>0</sub> and rake angle θ. Black curves are theoretical estimates of boundaries between certain behaviors. For subshear and self-arresting ruptures, only cases near the theoretical boundary were simulated. (b) Normalized depth-averaged rupture speed *v*<sub>r</sub>/*v*<sub>S</sub> versus normalized distance *L*/*W*, color coded by rake angle, for *G*<sub>c</sub>/*G*<sub>0</sub> = 0.63. *v*<sub>R</sub>, *v*<sub>S</sub>, *v*<sub>E</sub>, and *v*<sub>P</sub> are the Rayleigh, S, Eshelby and P wave speeds, respectively. (c) Normalized steady rupture speed and (d) real speed angle as a function of *G*<sub>c</sub>/*G*<sub>0</sub> and θ. Only fast, slow supershear, and “forbidden” speeds (see legend in (a)) are shown in (c).

**Figure 3 Comparison between observed and simulated speeds and synopsis of rupture behaviors.** (a) Symbols: rupture speeds and rake angles of global earthquakes, color coded by aspect ratio, compiled from references⁷⁻¹⁴, SRCMOD¹⁶, and USGS. Error bars indicate uncertainties when reported in references. Speeds are normalized by *v*<sub>s</sub> values compiled from references⁷⁻¹⁴ and PREM⁵⁵. Solid and dashed curves: rupture speeds from simulations (*G*<sub>c</sub>/*G*<sub>0</sub> > 0.5) and conceptual sketch motivated by theory (*G*<sub>c</sub>/*G*<sub>0</sub> < 0.5,
Methods A1), respectively. (b) Bottom: summary of the five rupture behaviours. Top: rupture speeds versus $G_c/G_0$ at a fixed rake (conceptual sketch), along profile AA' in the bottom plot. Purple and black curves: steady-state and distance-averaged (at a fixed rupture length and initial rupture speed) rupture speeds, respectively. Gray box shows all possible values of averaged speeds.

**Figure 4 Conceptual implications for time-dependent seismic hazard assessment.**

(a) Empirical estimates of fracture energy $G_c$ versus final slip $D^{26–30}$ (symbols; error bar indicates uncertainty when available), power-law fit (black thick line), examples of theoretical $G_c – D$ relation³⁰ (blue curve) and $G_0 – D$ relations (black thin lines). (b) Spatial distribution (left) and depth-averaged profile (right) of slip deficit rate (seismic coupling ratio³³ times plate convergence rate ~66 mm/yr) of the southern Andes subduction zone, Chile; epicenter (red star) and rupture area (green curve) of the 1960 Valdivia earthquake.

(c) Elapsed time since the 1960 earthquake needed to accumulate the critical slip deficit for runaway rupture, $D_{run}$. The colored band accounts for uncertainties of the $G_c – D$ scaling relation. Blue lines indicate potentially-runaway segments.
<figure><img src="image_2.png" /><figcaption>Figure 1</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 2</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 3</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>b</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 4</figcaption></figure>

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

<figure><img src="image_11.png" /></figure>
# Continuum of earthquake rupture speeds enabled by oblique slip

Huihui Weng¹ and Jean-Paul Ampuero¹

¹Université Côte d'Azur, IRD, CNRS, Observatoire de la Côte d'Azur, Géoazur, 250 rue Albert

Einstein, Sophia Antipolis, 06560 Valbonne, France

## Contents

### 5 Supplementary Figures

Figure S1.

Figure S2.

Figure S3.

Figure S4.

Figure S5.
<figure><img src="image_12.png" /><figcaption>**Figure S1:** (a) The definition of real speed and apparent horizontal speed. (b) The shape (left), distributions of real speed (middle) and horizontal speed (right) of steady rupture fronts across the depth (colored symbols coded by rake angle).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure S2:** (a) Normalized depth-averaged horizontal speed $v_r^{hor}$ (colored curves coded by rake angle) as a function of normalized distance $L/W$ from the 3D dynamic rupture simulations with $G_c/G_0 = 0.63$. (b) Dependencies of normalized depth-averaged horizontal speed on energy ratio and rake angle. (c) The same as (b), but for the real speed</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure S3:** Normalized depth-averaged speeds as a function of normalized distance *L/W* from the 3D dynamic rupture simulations with different nucleation speeds.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure S4:** (a) Rupture contours of a kinematic model with oblique slip whose rupture extends as an elliptical front propagating at the P wave speed along its major axis and at the S wave speed along its minor axis. The rake angle is the angle between the major axis of the ellipse and the strike direction. (b) The comparison of depth-averaged real speed between the kinematic and dynamic models with $G_c/G_0 = 0.63$ versus rake angle. The definition of real speed and horizontal speed are the same for both the kinematic and dynamic models. (c) The comparison of horizontal speed versus rake angle. (d) Dependency of depth-averaged real speed angle on rake angle for both kinematic and dynamic models.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure S5:** Sketches of $G_c$ and $G_0$ for linear (a) and power-law (b) slip-weakening friction laws.</figcaption></figure>