# Signature of transition to supershear rupture speed in the coseismic off-fault damage zone

Jorge Jara, Lucile Bruhat, Marion Y. Thomas, Solène Antoine, Kurama

Okubo, Esteban Rougier, Ares Rosakis, Charles Sammis, Y Klinger, Romain

Jolivet, et al.

► To cite this version:

Jorge Jara, Lucile Bruhat, Marion Y. Thomas, Solène Antoine, Kurama Okubo, et al. Signature of transition to supershear rupture speed in the coseismic off-fault damage zone. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2021, 477 (2255), 10.1098/rspa.2021.0364. hal-03443523

HAL Id: hal-03443523

https://hal.science/hal-03443523

Submitted on 9 Dec 2021

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
Research

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

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

Check for
updates

**Cite this article:** Jara *J et al*. 2021 Signature of transition to supershear rupture speed in the coseismic off-fault damage zone. *Proc. R. Soc. A 477:* 20210364.

https://doi.org/10.1098/rspa.2021.0364

Received: 3 May 2021

Accepted: 21 October 2021

**Subject Areas:**

geophysics, geology, mechanics

**Keywords:**

supershear earthquakes, off-fault coseismic
damage, earthquake mechanics and
observations

**Author for correspondence:**

Jorge Jara

e-mail: jara@geologie.ens.fr

Signature of transition to
supershear rupture speed in
the coseismic off-fault damage
zone

Jorge Jara<fnref n="1" />, Lucile Bruhat<fnref n="1" />, Marion Y. Thomas<fnref n="2" />,

Solène L. Antoine<fnref n="3" />, Kurama Okubo<fnref n="4" />, Esteban

Rougier<fnref n="5" />, Ares J. Rosakis<fnref n="6" />, Charles G. Sammis<fnref n="7" />,

Yann Klinger<fnref n="3" />, Romain Jolivet<fnref n="1" />,<fnref n="8" /> and Harsha S. Bhat<fnref n="1" />

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

JJ, 0000-0003-3176-0689; LB, 0000-0002-3916-4147

Most earthquake ruptures propagate at speeds below the shear wave velocity within the crust, but in some rare cases, ruptures reach supershear speeds. The physics underlying the transition of natural subshear earthquakes to supershear ones is currently not fully understood. Most observational studies of supershear earthquakes have focused on determining which fault segments sustain fully grown supershear ruptures. Experimentally cross-validated numerical models have identified some of the key ingredients required to trigger a transition to supershear speed.

© 2021 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

Electronic supplementary material is available online at https://doi.org/10.6084/m9.fishshare. c.5705327.

<footnote n="1">¹Laboratoire de Géologie, Département de Géosciences, École Normale Supérieure, CNRS, UMR 8538, PSL Université, Paris, France</footnote>

<footnote n="2">²Institut des Sciences de la Terre de Paris, Sorbonne Université, CNRS, UMR 7193, Paris, France</footnote>

<footnote n="3">³Université de Paris, Institut de Physique du Globe de Paris, CNRS, Paris 75005, France</footnote>

<footnote n="4">⁴National Research Institute for Earth Science and Disaster Resilience, 3-1 Tennnodai, Tsukuba, Ibaraki 305-0006, Japan</footnote>

<footnote n="5">⁵EES-17-Earth and Environmental Sciences Division, Los Alamos National Laboratory, Los Alamos, NM, USA</footnote>

<footnote n="6">⁶Graduate Aerospace Laboratories, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="7">⁷Department of Earth Sciences, University of Southern California, Los Angeles, CA 90089, USA</footnote>

<footnote n="8">⁸Institut Universitaire de France, 1 rue Descartes, Paris 75005, France</footnote>
However, the conditions for such a transition in nature are still unclear, including the precise location of this transition. In this work, we provide theoretical and numerical insights to identify the precise location of such a transition in nature. We use fracture mechanics arguments with multiple numerical models to identify the signature of supershear transition in coseismic off-fault damage. We then cross-validate this signature with high-resolution observations of fault zone width and early aftershock distributions. We confirm that the location of the transition from subshear to supershear speed is characterized by a decrease in the width of the coseismic off-fault damage zone. We thus help refine the precise location of such a transition for natural supershear earthquakes.

# 1. Introduction

While most earthquake ruptures propagate at speeds below the shear wave speed (sub-Rayleigh regime), some earthquakes can occasionally accelerate above the shear wave speed. Such events are known as supershear earthquakes. The earthquake rupture speed, and specifically its abrupt changes, control the high-frequency radiation [1,2]. As an earthquake goes to supershear speeds, it manifests Mach fronts that produce unusually large ground motion at distances far from the fault [3–5]. Rupture speed also governs the spatial extent of the off-fault coseismic damage zone, i.e. the volume within the crust directly surrounding the rupture that experienced mechanical damage [6–11]. Coseismic off-fault damage refers to fractures created or reactivated in the off-fault volume due to the dynamic rupture. As variations in rupture speed affect the seismic radiation and the off-fault coseismic damage, understanding the conditions for a rupture to transition from the sub-Rayleigh to the supershear regime would greatly help constrain fault properties (geometry, friction, lithology, etc.) and traction conditions that promote supershear ruptures. In addition, knowing how, why and where earthquakes attain supershear speeds would help in the reliable estimation of earthquake hazard assessment [1,12].

Whether supershear ruptures occur in nature was a matter of debate for a long time. While theoretical and numerical models demonstrated in the early 1970s that earthquakes could propagate at supershear speeds [13–19], the absence of field observation and their extreme rarity in laboratory experiments [20] first suggested that supershear earthquakes could not exist in nature. It was not until the $M_{w}$ 6.5 Imperial Valley earthquake (California, 1979) that a supershear rupture was inferred for the first time [21–23]. Pioneering laboratory experiments [24–26] together with observations from the 1999 $M_{w}$ 7.4 Izmit and the 1999 $M_{w}$ 7.2 Düzce earthquakes in Turkey [27,28], then conclusively confirmed that supershear ruptures are more common than previously expected. Supershear ruptures have now been inferred for several, albeit rare, events: the 2001 $M_{w}$ 7.8 Kunlun (China) earthquake [29–31], the 2002 $M_{w}$ 7.8 Denali (Alaska) earthquake [32,33], the 2010 $M_{w}$ 6.9 Qinghai (China) earthquake [34], the 2012 $M_{w}$ 8.6 off-Sumatra (Indonesia) earthquake [35], the 2013 $M_{w}$ 7.5 Craig (Alaska) earthquake [36], the 2013 $M_{w}$ 6.7 Okhotsk (Kamtchatka) earthquake [37,38] and most recently the 2018 $M_{w}$ 7.5 Palu (Indonesia) earthquake [39,40].

Several efforts have been made to identify and characterize these supershear earthquakes using various methods [28,31,39,41,42]. The methods are, in most cases, designed to reveal along which segment the rupture propagated at supershear velocities. The location of the transition from sub to supershear speeds is usually inferred as follows. The fault is subdivided into segments for which an average rupture speed is determined, using kinematic inversion methods. If the average rupture speed of a segment is deduced to be supershear, supershear transition is presumed to have occurred in between subshear and supershear segments. However, due to the lack of dense near-fault records, rupture speed estimates from seismological inversions and back-projection techniques rely on teleseismic records, leading to rough estimates of the transition’s location (i.e. tens of kilometres precision). Even in a situation where near-field data exists, such as during the 2002 Denali earthquake [32,33], records are still too sparse to precisely locate the transition.
Numerous theoretical and numerical models have been developed to characterize the mechanics of supershear transition. On planar faults with homogeneous stress and strength, the occurrence of supershear ruptures depends on the S ratio: S = (τp - τ0)/(τ0 - τr), where τ0 is the background shear stress and τp and τr are the peak and residual frictional strengths [15–17]. Supershear ruptures occur if S < 1.77 in two dimensions [15,17], or S < 1.19 in three dimensions [43], and the fault length is long enough for transition to occur [15]. The time and location of the transition are determined by how far S is from the Burridge–Andrews threshold of 1.77. However, natural faults do not have homogeneous stress and strength conditions before an earthquake. Numerical studies of supershear transition have also focused on spatial heterogeneity in stress and/or strength distribution, on a planar fault [43–45] and on a non-planar fault [46–48], while Kaneko & Lapusta [49] and Hu et al. [50] looked at the effect of the free surface on supershear transition.

Despite these efforts, the conditions for supershear transition in nature are still poorly understood because the imprecise location of the transition does not allow us to properly characterize the conditions on, or off, the fault that led to it. In addition, it also appears difficult to precisely capture this transition in both laboratory earthquakes (e.g. Mello [51]) and numerical models. This problem is mainly because the transition happens over very small length scales (of millimetres in the lab and over a few elements in numerical methods). A new approach that focuses specifically on the characterization of the region of transition is thus required to better understand why earthquakes accelerate to supershear speeds.

Although the modelling efforts help better understand the physics of supershear ruptures, it is not straightforward to apply that knowledge in the field, especially since we do not directly measure the stress and strength along the fault. Therefore, we propose a new approach based on the physical characterization of parameters that may promote supershear ruptures. This approach focuses on the detailed analysis of coseismic off-fault damage as a proxy to characterize the transition to supershear speeds. In this work, we cross-validate this idea by developing a theoretical model and comparing it with several numerical approaches. Then, we look at the observations available for natural supershear ruptures to corroborate this approach.

## 2. Coseismic off-fault damage around the rupture tip while approaching supershear transition: theoretical analysis

Coseismic off-fault damage depends on the evolution of stress during the propagation of a rupture in the crust. We use Linear Elastic Fracture Mechanics (LEFM) to provide closed-form solutions to quantify the state of stress around the rupture tip.

Consider a semi-infinite plain-strain crack that moves at speed $v \le c_R$, where $c_R$ is the Rayleigh wave speed. Let the origin of the polar coordinate system $(r, \theta)$ coincide with the crack tip. The near-tip stress, which also depends on the rupture speed $v$, is given as [52],

$$ \sigma_{\alpha\beta}(r, \theta, v) = \frac{K_{\text{II}}^{\text{dyn}}(v)}{\sqrt{2\pi r}} f_{\alpha\beta}^{\text{II}}(\theta, v) + \sigma_{\alpha\beta}^{0} \quad \alpha, \beta = 1, 2 \quad (2.1) $$

and

$$ \sigma_{33} = v(\sigma_{11} + \sigma_{22}), \quad (2.2) $$

where $\sigma_{\alpha\beta}^0$ is the initial stress state and $K_{\text{II}}^{\text{dyn}}$ is the dynamic stress intensity factor, which can be approximated as

$$ K_{\text{II}}^{\text{dyn}} \approx \frac{1 - v/c_R}{\sqrt{1 - v/c_p}} K_{\text{II}}^{\text{sta}}(\hat{L}), \quad (2.3) $$

where $c_p$, the P-wave speed, is the limiting speed for a mode II crack and $\hat{L} = \int_0^t v(t)dt$ is the current crack length and $c_R$ is the Rayleigh wave speed [52]. Here $K_{\text{II}}^{\text{sta}}(L) = \Delta\tau\sqrt{\pi\hat{L}}$ is the static stress intensity factor of the crack and $\Delta\tau$ is the stress drop defined as $\sigma_{yx}^0 - \tau_r$. Here $\sigma_{yx}^0$ is the initial shear stress on the fault and $\tau_r$ is the residual frictional strength. This solution by Freund [53] allows us to transform a static solution into a dynamic one. The rupture velocity dependence
of the stress field makes it undergo a Lorentz-like contraction affecting both the stress field and its angular distribution as the rupture speed approaches the Rayleigh wave speed, $c_R$ [52]. This contraction has already been observed and verified experimentally for ruptures approaching the Rayleigh wave speed [54]. For a crack-like rupture, there is a competing effect of the increase in stress intensity factor with increasing crack length. If it is a stable slip pulse of fixed length, then the effect of the crack length is invariant with rupture speed.

When an earthquake rupture transitions to the supershear regime, the rupture first accelerates to the Rayleigh wave speed. Within the framework of LEFM, as the rupture approaches the Rayleigh wave speed, $K_{\text{II}}^{\text{dyn}}$ monotonically decreases to zero, strongly reducing the stress concentration at the rupture tip. Thus the off-fault domain affected by this stress concentration also shrinks. We illustrate this effect by calculating the extent of the region where the stress state exceeds the Drucker–Prager failure criterion, i.e. $F_{\text{DP}} \ge 0$, where

$$F_{\text{DP}} = \sqrt{J_2} + p \tan \phi. \qquad (2.4)$$

Here $p \equiv I_1/3 = \sigma_{kk}/3$ is the hydrostatic stress derived from the first invariant of the stress tensor, $I_1, J_2 = s_{ij}S_{ji}/2$ corresponds to the second invariant of the deviatoric stress tensor ($s_{ij} = \sigma_{ij} - p\delta_{ij}$) and $\tan \phi = f$ is the static coefficient of friction. This yield criterion, effectively a smooth approximation to the Mohr–Coulomb Law using invariants of the stress tensor [55], takes into account all possible planes of slip, as opposed to Mohr–Coulomb where the potential slip planes need to be defined a priori or are optimally oriented.

Using the LEFM solutions for a mode II crack, we can compute $F_{\text{DP}}$. The closed-form expressions are quite long, but the key variables that affect the solution can be written as

$$F_{\text{DP}} = \mathcal{F} \left[ \frac{1}{\sqrt{r}}, \sqrt{\hat{L}} \left( \frac{1 - v/c_R}{\sqrt{1 - v/c_p}} \right), \Delta\tau, \theta, v, \tan \phi, \sigma_{\alpha\beta}^0 \right]. \qquad (2.5)$$

Thus, given a rupture velocity history, $v(t)$, one can then calculate the region around the fault where the Drucker–Prager yield criterion is violated.

We illustrate the solution in figure 1. Consider a fault in an elastic medium subject to an initial stress state (see table 1 for all the parameters used). Assume that the rupture is accelerating to the Rayleigh wave speed (figure 1a). The instantaneous static stress intensity factor increases with the increasing length of the crack (figure 1b). However, Freund’s universal function, $F_{\text{II}}(v)$, monotonically decreases to 0 (figure 1c). Consequently, the dynamic stress intensity factor first increases with increasing crack length but soon starts decreasing as it is dominated by $F_{\text{II}}(v)$ (figure 1d). Thus, we should expect to see an initial increase in the spatial extent of the damage followed by a reduction as the rupture approaches the limiting speed. We illustrate this effect at key moments in the rupture’s history (marked by points numbered 1 through 9 in figure 1d). First, we compute, in figure 1e, the ‘instantaneous’ damage, i.e. not accounting for the damage accumulated from the past history of the rupture around the rupture tip. We observe that the spatial extent of the domain where Drucker–Prager failure criteria is violated first increases (time points 1 through 3) and then decreases (time points 4 through 9) as the rupture accelerates. We can now superimpose all the snapshots of instantaneous damage to visualize the ‘cumulative’ damage accumulated by the rupture throughout its history (figure 1f). Note that this last exercise is purely for illustrative purposes as the LEFM solution does not really account for the history of the rupture. These limitations will be overcome in the following sections. Nevertheless, to the first order, this approach illustrates that, as the rupture approaches the Rayleigh wave speed, the width of the coseismic off-fault damage zone should decrease significantly.

It is also worth noting that in the calculations above, the off-fault medium remains elastic. Thus, the rupture is insensitive to changes in the constitutive response of the medium due to coseismic off-fault damage. To resolve these issues, we need to take advantage of numerical simulations that allow for proper constitutive descriptions of the off-fault medium, including feedback from changes in elastic properties on the rupture behaviour and that also allow for transition from subshear to supershear speeds.
<figure><img src="image_5.png" /><figcaption>**Figure 1.** LEFM solution illustrating the temporal evolution of damage around a rupture tip. All plots are normalized by their maximum values. (a) Prescribed rupture velocity history normalized by the Rayleigh wave speed. (b) Evolution of the static stress intensity factor, which is proportional to the length of the crack, $\sqrt{L(t)}$. (c) Freund's universal function. Here $c_R$ and $c_p$ are the Rayleigh and P-wave speeds, respectively. (d) Evolution of the dynamic stress intensity factor. (e) Instantaneous snapshots of the contours of the domain where the Drucker–Prager failure criterion is violated, $F_{DP} \ge 0$, at key points labelled 1 through 9 in (d). (f) Superposition of the contours where $F_{DP} \ge 0$ throughout the history of the rupture. (Online version in colour.)</figcaption></figure>

### 3. Coseismic off-fault damage around the rupture tip during supershear transition: numerical analysis

While many efforts have been made to reveal the dynamics of supershear ruptures [6,43–46,49,56], the importance, during such an event, of the complex feedback between the dynamic rupture and the mechanically evolving medium was not fully considered. However, new modelling approaches have recently emerged, allowing for dynamic activation of coseismic off-fault damage around faults and its feedback on rupture dynamics. In this section, we discuss two modelling
strategies that allow for off-fault damage associated with brittle fracture. In the following, we show that the observed damage pattern does not depend on any specific constitutive law used to account for damage. Instead, crustal damage is a universal feature that can aid in identifying the region of supershear transition in the field.

## (a) Modelling strategies

The two numerical models presented in this study account for fracture damage-dominated brittle rheology. Earthquake ruptures induce large dynamic strain rates (≥1 s⁻¹) around the rupture tip that cannot be fully accounted for by classical Mohr–Coulomb or Drucker–Prager plasticity, typically used at low strain rates (≤10⁻⁵ s⁻¹). The two models presented below allow us to account for such large dynamic strain rates. Okubo et al. [10,11] explicitly model the sub-kilometric off-fault fractures using a continuum–discontinuum numerical framework. The micromechanical model [8,9] accounts for coseismically activated ‘microfractures’ (approx. tens of metres) by computing the dynamic changes in the constitutive response due to these fractures. The difference between the two models is essentially in the scale of modelled off-fault fracture networks, leading to two distinct modelling strategies. In essence, the two brittle failure models are complementary, and they should, ideally, be combined. Nevertheless, they are good proxies for modelling off-fault damage at various scales, and we further show that they produce similar results.

To reproduce coseismic damage at relatively large scales (i.e. approx. 50 m to several km), Okubo et al. [10,11] use the software suite HOSSedu, developed by the Los Alamos National Laboratory [57,58]. The numerical algorithms behind this tool are based on the combined Finite-Discrete Element Method (FDEM) proposed by Munjiza [59], to produce dynamically activated off-fault fracture networks. The key feature here is that FDEM allows each interface between the finite elements, describing the off-fault medium, to have its own tensile and shear failure criterion. Thus, these interfaces can rupture under appropriate traction conditions. When the earthquake rupture propagates, the dynamic stress field around its tip will increase to violate the tensile or shear failure criteria leading to off-fault fracture damage. The lower limit of fracture resolution is around 50 m, which is defined by the minimum size of the discretized mesh. For a further detailed description of the method, see Okubo et al. [10,11].

The second modelling strategy relies on laboratory experiments [60–62] and field observations [63,64] that show significant changes in elastic properties related to fracturing. Observations in the field show up to 40% coseismic reduction in P- and S-wave velocities on spatial scales of hundreds of metres normal to the fault and few kilometres in depth. Following the damage constitutive laws proposed by Ashby & Sammis [65], Deshpande & Evans [66] and Bhat et al. [67], Thomas et al. [8] implemented an energy-based micromechanical brittle rheology, as developed by Rice [68], to account for such dynamic change of bulk rheological properties during earthquakes. In short, at each equilibrium state, the Gibbs free energy density ψ of the damaged solid is defined as the sum of (i) the free energy ψ^e of a solid, without flaws, deforming purely elastically and (ii) the free energy ψ^i corresponding to the contribution of the current set of microcracks. Using thermodynamic arguments, ψ = ψ^e + ψ^i is used to derive the new stress–strain constitutive law and the changes of elastic properties in the medium at the equilibrium stage. In the model, following laboratory experiments, the evolution of ψ^i is determined by taking into account the effect of loading rate and crack-tip velocities on crack growth (see Zhang & Zhao [69] for an overview). A complete description of the model can be found in Bhat et al. [67] and Thomas et al. [8].

## (b) Model set-up

For comparison purposes, we set up the two brittle rheology models to be as close to each other as possible. In both cases, we consider a two-dimensional plane strain medium, with a one-dimensional right lateral fault embedded in it and loaded by uniform background stresses. The
maximum compressive stress $\sigma_1$, and the minimum compressive stress $\sigma_3$ are in the x-y plane, whereas the intermediate principal stress $\sigma_2$ coincides with $\sigma_{zz}$. The fault plane makes an angle of 60° with $\sigma_1$ with uniform normal traction ($\sigma_{yy}^0$) and shear traction ($\sigma_{xy}^0$) everywhere except in the nucleation zone of the micromechanical model where the shear traction is slightly above the nominal static strength. In the FDEM model, the dynamic rupture is initiated by locally decreasing the static friction, in the nucleation zone, instead of the local change of $\sigma_{xy}^0$. These different nucleation strategies do not affect the results since this study is focused on damage occurring when the rupture is dynamic (sub and supershear). For each model, the initial shear stress and the S-ratio were chosen so that the rupture transitions to supershear speed reasonably early on.

Rupture propagation along the main fault plane is governed by a slip-weakening friction law [70]. Static friction is set at 0.6, which corresponds to a value measured in laboratory experiments for a large range of rocks [71], and dynamic friction at 0.1, as observed in high slip-rate experiments [72].

The fault length is 64 km for the micromechanical model and 115 km for the FDEM model. The domain width is determined based on the numerical method and on the scale of the fracture networks accounted for in the calculations. For the micromechanical model, the width around the fault is 6 km with absorbing boundary conditions to avoid reflections interfering with the propagating dynamic rupture. For the FDEM model, we simply set the domain large enough (86 km) so that the reflections do not arrive on the fault over the computation duration. We non-dimensionalize all length scales by the static size of the process zone ($R_0$) so that the models could be compared with natural or laboratory earthquakes. $R_0$ is defined as

$$R_0 = \frac{9\pi}{16(1-\nu)} \frac{\mu G}{(\tau_p - \tau_r)^2} \quad (3.1)$$

Here $\nu$ is Poisson's ratio, $\mu$ is shear modulus, G is the fracture energy associated with the slip-weakening law, $\tau_p$ and $\tau_r$ are the static and dynamic frictional strengths, respectively.

Reference values for the common parameters between the two models are summarized in table 1, whereas parameters specific to different modelling strategies are provided in the electronic supplementary material, tables S1 and S2. Schematics and parameters for the simulations can also be found in the electronic supplementary material, figure S1.

### (c) Influence of rupture dynamics on off-fault damage

To understand the complex interaction between the main fault and the surrounding medium undergoing coseismic damage, it is necessary to simultaneously look at the rupture dynamics, the associated stress field around the fault and the triggered damage in bulk. Here we use the micromechanical model to underline the key features that arise from numerical studies. We note that the same conclusions can be drawn from the FDEM model.

In the micromechanical model, under a compressive regime, damage occurs by the growth of pre-existing 'flaws'. They represent the faults, joints, cracks, mineral twins, defects in the crystal structure, grain boundaries, etc., we observe in nature. Frictional sliding occurs on these pre-existing fractures when the shear stress overcomes the frictional resistance acting on the fracture interface. As the faces slide in opposing directions, it creates a tensile wedging force that opens wing cracks at the tips of the shear fracture. The wing cracks grow when their stress intensity factors overcome the fracture toughness. In figure 2, we display the Drucker-Prager criterion to emphasize the regions where shear sliding is likely to occur, i.e. $F_{DP} > 0$ (equation (2.4)), and therefore where we expect wing cracks to grow. Fault slip rate (white curves) and cumulative damage (greyscale) are superimposed on each snapshot of $F_{DP}$. New coseismic damage induced by the displayed stress field is highlighted in red. As the rupture propagates below the Rayleigh wave speed (figure 2a), the damage is essentially generated behind the rupture front, in the tensional quadrant. This feature has been seen by Okubo et al. [10] as well, using the FDEM model. Damage also occurs ahead of the rupture front, where the S-wave field concentrates,
<table><caption>**Table 1.** Parameters for the LEFM, micromechanical and FDEM models. *Micromechanical model. †FDEM model.</caption><thead><tr><th>parameters</th><th>description</th><th>values</th></tr></thead><tbody><tr><td>state of stress</td><td></td><td></td></tr><tr><td>-σ<sup>0</sup><sub>xx</sub></td><td>fault-parallel stress (MPa)</td><td>37.7</td></tr><tr><td>-σ<sup>0</sup><sub>yy</sub></td><td>normal stress on the fault (MPa)</td><td>60.7</td></tr><tr><td>σ<sup>0</sup><sub>yx</sub></td><td>shear stress on the fault (MPa)</td><td>19.9</td></tr><tr><td></td><td>nucleation zone* (MPa)</td><td>36.4</td></tr><tr><td>ψ</td><td>angle between σ<sub>1</sub> and fault (°)</td><td>60</td></tr><tr><td>bulk properties</td><td></td><td></td></tr><tr><td>ρ</td><td>density (kg m<sup>-3</sup>)</td><td>2700</td></tr><tr><td>c<sub>s</sub></td><td>S-wave speed (km s<sup>-1</sup>)</td><td>3.12</td></tr><tr><td>c<sub>p</sub></td><td>P-wave speed (km s<sup>-1</sup>)</td><td>5.62</td></tr><tr><td>ν</td><td>Poisson's ratio</td><td>0.276</td></tr><tr><td>fault parameters</td><td></td><td></td></tr><tr><td>f<sub>s</sub></td><td>static friction coefficient</td><td>0.6</td></tr><tr><td></td><td>nucleation zone†</td><td>0.3</td></tr><tr><td>f<sub>d</sub></td><td>dynamic friction coefficient</td><td>0.1</td></tr><tr><td>D<sub>c</sub></td><td>characteristic slip distance (m)</td><td>1.0</td></tr><tr><td>R<sub>0</sub></td><td>process zone size (m)</td><td>1057</td></tr><tr><td>S</td><td>S-ratio</td><td>1.2</td></tr></tbody></table>

although the resulting damage density is much smaller than behind the rupture front. The location of newly formed cracks corresponds to the region where $F_{DP}$ is positive. The transition to supershear velocity (figure 2b) directly impacts the generation of coseismic damage. Around $t = 6.5$ s a new pulse is generated ahead of the rupture (figure 3b). The transition to supershear velocity coincides with a decrease in the width of the damage zone (figures 2b and 3b-d). As we argued in §2, this is likely related to the decrease in the stress intensity factor as the rupture speed increases. When approaching the Rayleigh wave speed, $K_{II}^{dyn}$ becomes small, even if the total length of the fault that has ruptured ($\hat{L}$ in equation (2.3)), has increased significantly. Moreover, a snapshot of $F_{DP}$ at $t = 7.2$ s (figure 2b) shows that the first pulse is not strong enough to change the state of stress optimally, thus reducing even further $K_{II}^{dyn}$. These two factors likely explain why damage mostly occurs behind the front of the pulse propagating at sub-Rayleigh rupture speed. After $t = 7.2$ s, the second pulse ahead of the rupture front becomes stronger and the first pulse weaker (figure 3d,e). Consequently, damage now essentially occurs behind the rupture front propagating at supershear speed and in the area behind the Mach cone (figure 2c). Therefore, even if the sub-Rayleigh pulse propagates now inside the 'transition zone', it induces little to no further damage. Therefore, this local reduction in the damage zone width is expected to remain unchanged for the earthquake duration and could, in principle, be observed in the field.

Finally, it is worth noticing that while the surrounding medium records the transition to supershear rupture (via the damage zone width), the final displacement along the fault, that geodesy could provide, carries no such information. As one can see in figure 3a, the slip continues to accumulate in the 'transition zone' after the passage of the rupture, while the damage zone width remains unchanged (figure 2c).
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Damage and Mohr–Coulomb Failure with the micromechanical model. Temporal evolution of a dynamic rupture occurring on a right-lateral fault embedded in granite (solid grey line at $y = 0$), using the micromechanical model. The rupture is nucleated around (0, 0) denoted by the yellow star. Sub-figures display snapshots of the Drucker–Prager yield criterion (equation (2.4)) when the rupture propagates at subshear (a) and supershear velocities (c), and during the transition (b). Fault slip rate (white curves) and the cumulative damage density, *D* (grey scale) are superimposed. Instant damage happening in relation to the displayed stress field is underlined in red. The yellow star denotes the nucleation patch and *C* and *T* the compressional and the tensional quadrants, respectively. (Online version in colour.)</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Slip and slip rate evolution with the micromechanical model. Cumulative slip (a), and slip rate, $\dot{\delta}$ (b-e) on the fault are displayed at various time steps. Coloured curves correspond to the dynamic simulation with the damage evolution law, dotted grey curves represent a simulation with the same parameters but within a pure elastic medium. The grey box gives the location of the area where less damage is observed. (Online version in colour.)</figcaption></figure>

### (d) Comparison between the different modelling strategies

We now perform a similar exercise using the FDEM model and choose parameters such that the micromechanical and FDEM models are as close to each other as possible. We adopt a slightly different yield criterion to reflect on the way damage is accounted for in this numerical method (see §3a and Okubo et al. [10]). We evaluate the potential regions of failure, in shear, by calculating the invariant form of the Mohr-Coulomb yield function, $F_{MC}$ (see eqn (4.142) of Chen & Han [73]):

$$F_{MC} = R_{MC} \sqrt{f_2} + p \tan \phi - c, \quad (3.2)$$

where c is the cohesion and $R_{MC}$ is given by

$$R_{MC} = \sin\left(\theta + \frac{\pi}{3}\right) \sec\phi + \frac{1}{\sqrt{3}} \cos\left(\theta + \frac{\pi}{3}\right) \tan\phi, \quad (3.3)$$

where

$$\cos 3\theta = \frac{3\sqrt{3}}{2} \frac{J_3}{J_2^{3/2}} \quad \text{and} \quad J_3 = \frac{1}{3} s_{ij} s_{jk} s_{ki}, \quad (3.4)$$
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Damage and Mohr–Coulomb Failure with the FDEM model. Dynamic rupture occurring on a right-lateral fault embedded in granite, using the FDEM model. The figure displays the invariant form of the Mohr–Coulomb yield criterion (equation (3.2)) when the rupture propagates at supershear velocities. Black lines give the spatial distribution of off-fault fractures that occurs during the entire event. Signature of the transition from the sub-Rayleigh to the supershear regime is highlighted by the gap in off-fault damage. (Online version in colour.)</figcaption></figure>

where $J_3$ is the third invariant of the deviatoric stress tensor. Failure occurs when $F_{MC} \ge 0$. Note that when $R_{MC} = 1$ and $c = 0$, $F_{MC}$ becomes the Drucker–Prager criterion as described in equation (2.4).

Figure 4 displays the $F_{MC}$ criteria when the rupture propagates at supershear velocities and the cumulative damage pattern resulting from the entire coseismic rupture. On top of the area undergoing $F_{MC} > 0$ in relation to the rupture propagation along the main fault plane, we observe positive values at the tip of the discretized off-fault fractures, far away from the rupture front. This feature is not so evident within the micromechanical model because the fractures are homogenized in the constitutive law (§3a). We also observe a shrinkage in the spatial extent of damage when the rupture transitions to supershear, as seen earlier with the micromechanical model.

In addition to the above-described numerical experiments using brittle rheology, we acknowledge that Templeton & Rice [7] first explored the extent and distribution of off-fault plasticity during supershear rupture (see fig. 13 of Templeton & Rice [7]). They conducted a two-dimensional in-plane dynamic rupture modelling with Drucker–Prager elasto-plastic constitutive law. In their model, the degree of damage is inferred from the accumulation of plastic strain. Although plasticity, as a proxy for damage, essentially holds only at low-strain rates [69,74], they nevertheless pointed out a 'remarkable contraction' in the damage zone width associated with supershear transition.

We can thus conclude, quite confidently, that this reduction in damage zone width (a gap in off-fault fractures) is a universal characteristic of supershear transition and is insensitive to the constitutive law used to model damage.

## 4. Natural observations of coseismic off-fault damage

Theoretical and numerical models, with three different rheological descriptions of off-fault damage [7,8,10], all suggest that the region affected by the stress field around the rupture tip shrinks during the supershear transition. This shrinkage of the stress field results in a narrow off-fault damage zone. We now look for a direct (or indirect) signature of this reduced damage zone in the case of natural supershear earthquakes, verifying their location with corresponding kinematic models that invert for rupture speeds.
Vallage et al. [75] and Klinger et al. [76] showed that there is a one-to-one relationship between the features of the displacement field, around the fault, of an earthquake and off-fault fracture damage. Using the FDEM method described earlier, Klinger et al. [76] showed that the observed, spatially diffuse, displacement field around a fault is due to displacement accommodated by off-fault fractures (see fig. 3 in [76]). This feature is in stark contrast to the sharp displacement field obtained, around a fault, when displacement is only accommodated by the main fault. Thus deviations from this sharp displacement field would allow us to characterize the width of the damage zone around a fault.

We can also make a reasonable assumption that the newly created/reactivated off-fault fractures are likely to host early off-fault aftershocks (based on the observed stress field). We should then expect a reduction in the spatial extent of the early off-fault aftershocks at the location of the supershear transition. A more thorough analysis would involve continuing the simulations described earlier for up to a week, or so, after the earthquake. However, computational limitations make this beyond the scope of the present study.

## (a) Optical image correlation to observe off-fault coseismic damage

Recent developments in satellite optical image analysis and sub-pixel correlation methods allow for detecting displacement variations due to an earthquake down to sub-metric resolutions (e.g. [75,77–79]). These methods enable characterizing the surface rupture geometry, the amount of surface displacement and the width of the zone affected by this displacement after an earthquake, also referred to as the ‘fault zone width’ [80]. The fault zone width reflects the lateral extent of rock damage, and secondary deformation features around the fault [81,82]. In the field, the fault zone width corresponds to the last deformation features observed when moving away from the fault core [81]. The same definition is employed in geodetic studies providing a more detailed characterization of this region due to the density and resolution of the observations [76].

The 2001 $M_{w}$ 7.8 Kunlun (China) earthquake is a strike-slip event with approximately 400 km long surface rupture and a mean rupture speed between 3.3 and 3.9 km s$^{-1}$, which is larger than the shear wave speed and hence reported as a supershear rupture [29,30]. This earthquake was well recorded by satellite images, allowing us to investigate the spatial distribution of coseismic off-fault damage by focusing on the fault zone width. We use SPOT-1 to SPOT-4 images, covering the 2001 Kunlun fault area from 1988 to 2004 with a ground resolution of 10 m, allowing for change detection of less than 1 m [75,83]. We focus on the central part of the rupture dominated by fault parallel motion [84,85], obtaining through optical image correlation the horizontal displacement field for the earthquake east-west and north-south components (electronic supplementary material, figure S2). The maximum strike-slip displacement for this area of the rupture is around 8 m [84,85].

Prior to correlation, SPOT images are corrected from viewing angle and topography using a unique 90 m SRTM digital surface model. The pre-and post-earthquake corrected paired images are correlated using the MicMac package [83,86], which measures the horizontal surface displacement between the acquisitions (see electronic supplementary material, figure S2 to see the data employed in this study). The MicMac method enables preserving the resolution of the input images, that is 10 m, since it measures the displacement at every pixel location in the pre-earthquake image. The size of the pixel pattern considered in the pre-earthquake images, and that is searched for in the post-earthquake ones, is 5 pixels, so 50 m. A regularization of 0.3 is used to force MicMac to consider realistic displacement values relative to neighbouring pixels, even if this procedure assumes pixels with a lower correlation score (regularization 0 in Micmac means taking the best correlation score without regarding on the coherence of surrounding pixels). Output displacement maps show displacement in the east-west and north-south directions. We measure 460 stacked profiles perpendicular to the fault trace in the displacement map (electronic supplementary material, figure S2) every 500 m to analyse variations in the surface displacements and fault zone width. Individual slip profiles are stacked on rectangles of 40 km-long and 1 km-wide, in order to average out noise and artefacts (Stacking profiles from ENVI v. 5.5.1; Exelis
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Optical image correlation analysis. (a) Map of the strike-slip section of the 2001 *M*<sub>w</sub> 7.8 Kunlun earthquake (China), where P1 denotes the transition zone reported for the event from seismological far-field data [31]. (b) Along-Strike fault zone width (black) and its associated uncertainty (grey), obtained from the analysis of 40 km long profiles, sampling the fault zone every 500 m, on the surface displacement maps. The latter is derived from correlating pre- and post-earthquake SPOT-1 to SPOT-4 images. The 11 km long red area highlights a region with a mean fault width (red dashed line) of only 127 m compared with 238 m recorded for the rest of the rupture (red line). The latter excludes the area where two parallel fault strands are activated and for which the fault zone is exceptionally large (greater than 1000 m). (c) Zoom of figure 5b. (Online version in colour.)</figcaption></figure>

Visual Information Solutions, Boulder, Colorado). On each profile, displacement is estimated in the fault-parallel direction, observing the coseismic offset produced by the earthquake (see electronic supplementary material, figure S3 for a profile example). This procedure allows us to evaluate the fault zone width in the strike-slip direction, over 396 profiles (figure 5b). Field measurements of fault zone widths are at least a few tens to a few hundred metres in the case of this earthquake [84,85]. Thus, we infer that the fault zone width can be estimated using displacement maps, employing satellite image correlation.

The analysis of the profiles highlights three main domains in fault zone width. The first domain, which represents the majority of the rupture, exhibits a mean value of 238 ± 80 m, and it corresponds to the strike-slip section of the rupture where only one fault zone is observed at the surface (figure 5; electronic supplementary material, figures S4A and S5). The second domain presents a width larger than 1 km lasting for about 15 km. This section shows a complex rupture pattern, where the rupture splits into two fault zones evolving into slip-partitioning between pure strike-slip along the southern fault strand and pure normal faulting along the northern strand. The slip-partitioning lasts for about 70 km eastward [84,87] (figure 5; electronic supplementary material, figures S4B and S5). The third domain is located right before the onset of the slip-partitioning region, where the average fault zone width drops to 127 ± 39 m over 11 km long segment. This value is about half of the mean fault zone width of 238 m calculated for the whole fault zone (ignoring the slip-partitioning area). This area co-locates with the region where the rupture is inferred, from teleseismic data, to have transitioned to the supershear regime [29–31]. The onset of slip-partitioning cannot cause a reduction in the damage zone width because the normal faulting motion was quite likely triggered by the dynamic stress field of the main
supershear segment [6]. This means that the overall dynamics in that region was still dominated by the rupture on the strike-slip segment. With these observations, and considering the theoretical analysis and numerical modelling conducted in the previous sections, we interpret this localized fault zone width reduction as a shrinkage of the off-fault damage zone associated with the rupture transitioning to supershear speeds (e.g. figure 2). Hence, this along-strike feature seen by the optical correlation analysis represents a natural observation of the supershear transition zone.

## (b) Distribution of early aftershocks as a proxy for off-fault coseismic damage

Assuming that the nucleation of early aftershocks is mainly governed by the stress state left in the wake of the earthquake, we test whether the extent of coseismic off-fault damage in nature manifests in the spatial distribution of early off-fault aftershocks. The logic is that considering the state of stress in the wake of the earthquake, weakened regions will preferentially host early aftershocks, and thus we should observe a region with less events where the rupture transitioned to supershear. When re-analysing the regions where established supershear ruptures transitioned from sub-Rayleigh to supershear speeds, we might find local minima in the spatial extent of early aftershocks. If this region coincides with kinematic observations of supershear transition, we would have an independent and more precise location of supershear transition.

We analyse three reported supershear ruptures for which high-resolution aftershock catalogues are available, at least up to one month after the main event. The 1999 $M_{w}$ 7.4 Izmit earthquake (Turkey) had a mean rupture velocity of about 4.8 km s⁻¹ [27,28], and a well-recorded three-month aftershock catalogue [41] (figure 6a for one-month aftershock locations). In the case of the 2002 $M_{w}$ 7.9 Denali earthquake (Alaska), rupture velocities of about 5.5 km s⁻¹ were reported [32], and 1-year long catalogue of aftershocks is available [88] (see figure 6b for one-month aftershock locations). The 2013 $M_{w}$ 7.5 Craig earthquake (Alaska) was also inferred as a supershear event, with rupture velocities between 5.5 and 6.0 km s⁻¹ [36] and a five-month aftershock catalogue is available [89] (see figure 6c for one-month aftershock locations).

We observe a paucity in aftershocks where the ruptures have transitioned to supershear, as highlighted by the pink boxes in figure 7b,d,f. We quantify such a lack of aftershocks by computing the cumulative moment density released by aftershocks, considering those located within $l_{a} = 5$ km of the fault trace (see the electronic supplementary material, for the same analysis considering $l_{a} = 2.5$ km; electronic supplementary material, figure S6). We assume that each aftershock, having a seismic moment $M_{0}$, is a circular crack and compute the slip distribution $\delta(r)$ using [90]:

$$ \delta(r) = \frac{24\Delta\tau}{7\pi\mu} \sqrt{\left(\frac{7M_0}{16\Delta\tau}\right)^{2/3} - r^2}, \quad (4.1) $$

where r is the crack radius, $\mu = 30$ GPa is the shear modulus and $\Delta\tau$ is the stress drop, assumed to be 3 MPa. We bin the aftershock crack radius along-strike, based on the minimum aftershock magnitude for each case (Izmit $M_{w}$ = 2.4, Denali $M_{w}$ = 1.5 and Craig $M_{w}$ = 1.7). Thus, the binning sizes employed to discretize the main fault are 90 m for Izmit, 32 m for Denali and 40 m for Craig earthquake. With the slip distribution, the cumulative seismic moment density of the ith bin containing N aftershocks is given by

$$ m_c^{(i)} \approx \mu \sum_{j=1}^{N} \delta_j^{(i)}. \quad (4.2) $$

Note that $\delta_{j}^{(i)}$ is part of the slip distribution of the jth aftershock projected on the ith fault segment. We also note that we do not compute the full moment density tensor as we do not have focal mechanisms for all the aftershocks. As we are interested in relative spatial variation (along the strike of the main fault) of this quantity, the above approximation is adequate. We compute the seismic moment density over two periods following the mainshock: one and three weeks (figure 7a,c,e). Focusing on the region where the rupture is expected to transition to
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Aftershock catalogues. One-month aftershock distribution for Izmit (a), Denali (b) and Craig (c) earthquakes, colour-coded by time. Pink stars indicate the mainshock epicentres with focal mechanism. The continuous black lines denote the surface rupture for each event, while the dashed ones indicate a distance of 5 km from the fault for each event. Brown dashed boxes are the target regions to explore the supershear transition based on published kinematic models (Bouchon *et al.* [27] for a, Ellsworth *et al.* [32] for b and Yue *et al.* [36] for c). The pink boxes indicate this work's proposed transition zones. (Online version in colour.)</figcaption></figure>

the supershear regime (based on kinematic models), we systematically observe a small area characterized by a reduced seismic moment density and lack of aftershocks (figure 7b, df, pink boxes). The extent of these regions is different for each earthquake cf., 3.6 km for Izmit, 6 km for Denali and 3.2 km for Craig. The difference in the smoothness of the cumulative moment curves for different earthquakes is due to the number of aftershocks available for the calculation. In the Craig earthquake case, we have 68 (83) aftershocks at a distance of 5 km from the main fault in one week (three weeks), 2396 (4265) for Denali and 641 (1118) for the Izmit earthquake.

In order to evaluate the robustness of our findings, we perform an uncertainty analysis of the seismic moment density estimation by accounting for errors in the location of aftershocks in each catalogue. We employ a multivariate normal distribution {$\mathbf{x} = \mathcal{N}_2(\mu_a, \sigma_a)$, assuming that the location uncertainty of aftershocks is normally distributed. Here, $\mu_a$ corresponds to the
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Spatio-temporal seismic moment density evolution. Cumulative seismic moment density projected on the main fault at different temporal scales (one to three weeks), for Izmit (a), Denali (c) and Craig (e) earthquakes. All the aftershocks within a distance of 5 km from the fault are considered in the calculation (area denoted by the black discontinuous lines in figure 6a–c). Colour-coded arrows (at the top in a, b and c) indicate the different speed regimes reported for each event (green for sub-Rayleigh and orange for supershear) [32,36,41], while the stars denote the epicentre of each earthquake and the arrows indicate the direction of the rupture. The pink boxes highlight our proposed transition zone, also observed in a map view in (b) for Izmit, (d) for Denali and (f) for Craig earthquakes. (Online version in colour.)</figcaption></figure>

aftershock epicentre, while $σ_a$ is the covariance matrix of the location. This procedure allows us to generate, randomly, 10 000 synthetic catalogues and repeat the along-strike moment density evaluation with time for each of the 10 000 catalogues, deriving the mean and the standard deviation (1 − σ and 10 − σ ) of the spatio-temporal evolution of seismic moment density (red
<figure><img src="image_12.png" /></figure>

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

<figure><img src="image_14.png" /><figcaption>**Figure 8.** High-resolution aftershock catalogue statistical analysis. Aftershock seismic moment density projected on the main fault for one week, considering all the earthquakes at a distance of 5 km from the fault for Izmit (a), Denali (b) and Craig (c) earthquakes. The black curve is the mean seismic moment density on the main fault ± the standard deviation (1 - σ) computed from 10 000 synthetic random catalogues. Pink boxes highlight our proposed transition zone. (Online version in colour.)</figcaption></figure>

areas in figure 8 for one week, denoting the mean of the seismic moment density on the main fault ± the standard deviation). Such analysis suggests that despite potential changes in the spatio-temporal distribution of seismic moment density, our conclusions are robust (see the electronic supplementary material for a comparison between one and three weeks using 1 − σ and 10 − σ ; electronic supplementary material, figures S7–S9).

We consider three weeks as the maximum reasonable period to investigate the evolution of
the seismic moment for these earthquakes to avoid potential effects of post-seismic deformation.
Afterslip is proposed as a mechanism driving the aftershock triggering (e.g. [91,92]). For each studied event, the inferred transitional region co-locates with areas where some authors report afterslip occurrence [93–95]. Therefore, the observed gap in the early aftershocks productivity (less than three weeks after the mainshock) is mainly related to the mainshock rupture. It is worth noting that results for the $M_w$ 7.8 Kunlun earthquake reported by Robinson *et al.* [30], employing Harvard CMT solutions, also allude to the same conclusion. However, due to the lack of high spatio-temporal density of aftershocks in their catalogue, we are more confident in our results using optical correlation techniques to characterize the transition zone.

We note that in the distribution of moment plotted for the entire fault profile (figure 7), local gaps in aftershock density exist and do not necessarily indicate a supershear transition. We hypothesize that such gaps could also be related to abrupt local changes in rupture speeds, off-fault medium strength and fault geometrical complexities, which warrants further investigation.

To conclude, we show that the reduction in the spatial extent of early aftershocks (less than or equal to three weeks after the mainshock), and the cumulative seismic moment density around the fault (less than or equal to 5 km on either side of the fault) is associated with the transition to supershear speed. This type of analysis also helps us identify, more precisely than kinematic models, the region where supershear transition occurs.

## 5. Conclusion and discussion

Using theoretical arguments and numerical models that account for coseismic off-fault damage, we have shown that supershear transition is characterized by a significant reduction in the width of the off-fault damage zone. This feature is due to a Lorentz-like contraction of the spatial extent of the stress field around a rupture tip. We then cross-validated this phenomenon with natural observations of coseismic off-fault damage zone width using image correlation and analysis of aftershock catalogues. We confirm that supershear transition is indeed characterized by a significant reduction in the width of the off-fault damage zone.

Our results are in general agreement with the published kinematic models for the Izmit [28,96], Denali [32,33], Craig [36] and Kunlun [29–31] earthquakes, relative to the location where the rupture accelerates to supershear speeds. Our inferred location of the transition zone can sometimes be different from those found in kinematic inversions, but only by a few kilometres.

While the theoretical and numerical models are still idealized (modelled on a one-dimensional planar fault with uniform traction and frictional strength), the reduction in the width of the coseismic off-fault damage zone, related to supershear transition, is observed for natural earthquakes (figures 5 and 7). This approach could thus be used to refine and precisely define the zone of supershear transition inferred from coarser kinematic models. Moreover, as the zone of supershear transition is now narrowed down to a few kilometres, it could also be used to guide future fieldwork to study the segment of the fault where such a transition occurred. A better understanding of the physical conditions for supershear transition might help us foresee the location of future supershear transition, although, as of now, this remains a difficult task.

The lack of aftershocks on the segment experiencing supershear rupture has been debated in the past. Bouchon & Karabulut [41], for instance, already concluded that the entire supershear segment collocated with a region of aftershock quiescence. To explain this, they claimed that the stress and strength conditions were more homogeneous on the supershear segments [41], impeding large ground accelerations near the fault [28,32]. Our work, by contrast, focuses on characterizing the region where the rupture transitions to supershear speed, before the said supershear segment.

The results of this study are valid for a well-developed sub-Rayleigh rupture that transitions to supershear speeds. However, the 2018 $M_w$ 7.5 Palu (Indonesia) earthquake might have either nucleated directly at a supershear speed, or transitioned very early [39,40,97,98]. The model we develop here might provide insights to understand whether an earthquake can nucleate and propagate directly at supershear speeds, and if so, why that would be the case.
**Data accessibility.** All the catalogues and numerical simulation results employed in this work have been obtained from published works, cited in the main text and electronic supplementary material.

**Authors' contributions.** H.S.B. conceived of and designed the study. H.S.B. developed the LEFM solution. M.Y.T. conducted damage mechanics-based numerical modelling and K.O. performed FDEM-based numerical modelling. S.A. conducted the image correlation analysis. J.J. analysed the aftershock catalogue with inputs from M.Y.T. and H.S.B. J.J. and L.B. drafted the first version of this manuscript, which was subsequently modified with inputs from H.S.B., M.Y.T., K.O., S.A., Y.K., A.J.R. and R.J. All authors read and approved the manuscript.

**Competing interests.** We declare we have no competing interests.

**Funding.** H.S.B. acknowledges the European Research Council grant PERSISMO (grant no. 865411) for partial support of this work. J.J. and R.J. acknowledge the funding of the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 758210, project Geo4D). L.B. thanks the funding from the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme (FP7/2007-2013) under REA grant agreement no. PCOFUND-GA-GA-2013-609102, through the PRESTIGE programme coordinated by Campus France. S.A. and Y.K. are partly supported by the ANR project DISRUPT (ANR-18-CE31-0012). R.J. acknowledges funding from the Institut Universitaire de France. Portions of this research were obtained using resources provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under contract no. 89233218CNA000001. This publication was approved for unlimited release under LA-UR-21-24016. SPOT images are from the ISIS programme from CNES.

## References

1. Madariaga R. 1983 High frequency radiation from dynamic earthquake. *Ann. Geophys.* **1**, 17.

2. Spudich P, Frazer LN. 1984 Use of ray theory to calculate high-frequency radiation from earthquake sources having spatially variable rupture velocity and stress drop. *Bull. Seismol. Soc. Am.* **74**, 2061–2082. (doi:10.1785/BSSA0740062061)

3. Bernard P, Baumont D. 2005 Shear Mach wave characterization for kinematic fault rupture models with constant supershear rupture velocity. *Geophys. J. Int.* **162**, 431–447. (doi:10.1111/gji.2005.162 issue-2)

4. Dunham EM, Archuleta RJ. 2005 Near-source ground motion from steady state dynamic rupture pulses. *Geophys. Res. Lett.* **32**, 1–4. (doi:10.1029/2004GL021793)

5. Dunham EM, Bhat HS. 2008 Attenuation of radiated ground motion and stresses from three-dimensional supershear ruptures. *J. Geophys. Res.: Solid Earth* **113**, 1–17. (doi:10.1029/2007JB005182)

6. Bhat HS, Dmowska R, King GCP, Klinger Y, Rice JR. 2007 Off-fault damage patterns due to supershear ruptures with application to the 2001 $M_{w}$ 8.1 Kokoxili (Kunlun) Tibet earthquake. *J. Geophys. Res.* **112**, B06301. (doi:10.1029/2006JB004425)

7. Templeton EL, Rice JR. 2008 Off-fault plasticity and earthquake rupture dynamics: 1. Dry materials or neglect of fluid pressure changes. *J. Geophys. Res.* **113**, B09306. (doi:10.1029/2007JB005529)

8. Thomas MY, Bhat HS, Klinger Y. 2017 Effect of brittle off-fault damage on earthquake rupture dynamics. In *Fault zone dynamic processes: evolution of fault properties during seismic rupture* (eds MY Thomas, TM Mitchell, HS Bhat), pp. 255–280. American Geophysical Union.

9. Thomas MY, Bhat HS. 2018 Dynamic evolution of off-fault medium during an earthquake: a micromechanics based model. *Geophys. J. Int.* **214**, 1267–1280. (doi:10.1093/gji/ggy129)

10. Okubo K, Bhat HS, Rougier E, Marty S, Schubnel A, Lei Z, Knight EE, Klinger Y. 2019 Dynamics, radiation, and overall energy budget of earthquake rupture with coseismic off-fault damage. *J. Geophys. Res.: Solid Earth* **124**, 11771–11801. (doi:10.1029/2019JB017304)

11. Okubo K, Rougier E, Lei Z, Bhat HS. 2020 Modeling earthquakes with off-fault damage using the combined finite-discrete element method. *Comput. Part. Mech.* **7**, 1057–1072. (doi:10.1007/s40571-020-00335-4)

12. Das S. 2007 Geophysics: the need to study speed. *Science* **317**, 905–906. (doi:10.1126/science.1142143)

13. Burridge R. 1973 Admissible speeds for plane-strain self-similar shear cracks with friction but lacking cohesion. *Geophys. J. Int.* **35**, 439–455. (doi:10.1111/j.1365-246X.1973.tb00608.x)
14. Hamano Y. 1974 Dependence of rupture time history on the heterogeneous distribution of stress and strength on the fault plane. *Eos Trans. AGU* **55**, 352.

15. Andrews DJ. 1976 Rupture velocity of plane strain shear cracks. *J. Geophys. Res.* **81**, 5679–5687. (doi:10.1029/JB081i032p05679)

16. Das S. 1976 A numerical study of rupture propagation and earthquake source mechanism. PhD thesis, Massachusetts Institute of Technology.

17. Das S, Aki K. 1977 A numerical study of two-dimensional spontaneous rupture propagation. *Geophys. J. R. Astron. Soc.* **50**, 643–668. (doi:10.1111/j.1365-246X.1977.tb01339.x)

18. Burridge R, Conn G, Freund LB. 1979 The stability of a rapid mode II shear crack with finite cohesive traction. *J. Geophys. Res.* **84**, 2210–2222. (doi:10.1029/JB084iB05p02210)

19. Freund LB. 1979 The mechanics of dynamic shear crack propagation. *J. Geophys. Res.* **84**, 2199–2209. (doi:10.1029/JB084iB05p02199)

20. Wu FT, Kuenzler H, Thomson KC. 1972 Stick-slip propagation velocity and seismic source mechanism. *Bull. Seismol. Soc. Am.* **62**, 1621–1628. (doi:10.1785/BSSA0620061621)

21. Olson AH, Apsel RJ. 1982 Finite faults and inverse theory with applications to the 1979 Imperial Valley earthquake. *Bull. Seismol. Soc. Am.* **72**, 1969–2001. (doi:10.1785/BSSA07206A1969)

22. Archuleta RJ. 1984 A faulting model for the 1979 Imperial Valley earthquake. *J. Geophys. Res.: Solid Earth* **89**, 4559–4585. (doi:10.1029/JB089iB06p04559)

23. Spudich P, Cranswick E. 1984 1979 Imperial Valley earthquake using a short baseline. *Bull. Seismol. Soc. Am.* **74**, 2083–2114. (doi:10.1785/BSSA0740062083)

24. Rosakis AJ, Samudrala O, Coker D. 1999 Cracks faster than the shear wave speed. *Science* **284**, 1337–1340. (doi:10.1126/science.284.5418.1337)

25. Xia KW, Rosakis AJ, Kanamori H. 2004 Laboratory earthquakes: the sub-rayleigh-to-supershear rupture transition. *Science* **303**, 1859–1861. (doi:10.1126/science.1094022)

26. Passelègue FX, Schubnel A, Nielsen S, Bhat HS, Madariaga R. 2013 From sub-rayleigh to supershear ruptures during stick-slip experiments on crustal rocks. *Science* **340**, 1208–1211. (doi:10.1126/science.1235637)

27. Bouchon M, Toksöz N, Karabulut H, Bouin MP, Dietrich M, Aktar M, Edie M. 2000 Seismic imaging of the 1999 Izmit (Turkey) rupture inferred from the near-fault recordings. *Geophys. Res. Lett.* **27**, 3013–3016. (doi:10.1029/2000GL011761)

28. Bouchon M, Bouin MP, Karabulut H, Toksöz MN, Dietrich M, Rosakis AJ. 2001 How fast is rupture during an earthquake? New insights from the 1999 Turkey earthquakes. *Geophys. Res. Lett.* **28**, 2723–2726. (doi:10.1029/2001GL013112)

29. Bouchon M, Vallée M. 2003 Observation of long supershear rupture during the magnitude 8.1 Kunlunshan earthquake. *Science (New York, N.Y.)* **301**, 824–826. (doi:10.1126/science.1086832)

30. Robinson DP, Brough C, Das S. 2006 The $M_w$ 7.8, 2001 Kunlunshan earthquake: extreme rupture speed variability and effect of fault geometry. *J. Geophys. Res.* **111**, B08303. (doi:10.1029/2005JB004137)

31. Vallée M, Landès M, Shapiro NM, Klinger Y. 2008 The 14 November 2001 Kokoxili (Tibet) earthquake: high-frequency seismic radiation originating from the transitions between sub-Rayleigh and supershear rupture velocity regimes. *J. Geophys. Res.* **113**, B07305. (doi:10.1029/2007JB005520)

32. Ellsworth WL et al. 2004 Near-field ground motion of the 2002 Denali fault, Alaska, earthquake recorded at pump station 10. *Earthquake Spectra* **20**, 597–615. (doi:10.1193/1.1778172)

33. Dunham EM. 2004 Evidence for a supershear transient during the 2002 Denali fault earthquake. *Bull. Seismol. Soc. Am.* **94**, S256–S268. (doi:10.1785/0120040616)

34. Wang D, Mori J. 2012 The 2010 Qinghai, China, earthquake: a moderate earthquake with supershear rupture. *Bull. Seismol. Soc. Am.* **102**, 301–308. (doi:10.1785/0120110034)

35. Wang D, Mori J, Uchide T. 2012 Supershear rupture on multiple faults for the $M_w$ 8.6 off Northern Sumatra, Indonesia earthquake of April 11, 2012. *Geophys. Res. Lett.* **39**, 6–10.

36. Yue H, Lay T, Freymueller JT, Ding K, Rivera L, Ruppert NA, Koper KD. 2013 Supershear rupture of the 5 January 2013 Craig, Alaska ($M_w$ 7.5) earthquake. *J. Geophys. Res.: Solid Earth* **118**, 5903–5919. (doi:10.1002/2013JB010594)

37. Zhan Z, Helmberger DV, Kanamori H, Shearer PM. 2014 Supershear rupture in a $M_w$ 6.7 aftershock of the 2013 Sea of Okhotsk earthquake. *Science* **345**, 204–207. (doi:10.1126/science.1252717)
38. Zhan Z, Shearer PM, Kanamori H. 2015 Supershear rupture in the 24 May 2013 $M_w$ 6.7 Okhotsk deep earthquake: additional evidence from regional seismic stations. *Geophys. Res. Lett.* **42**, 7941–7948. (doi:10.1002/grl.v42.19)

39. Bao H, Ampuero JP, Meng L, Fielding EJ, Liang C, Milliner CWD, Feng T, Huang H. 2019 Early and persistent supershear rupture of the 2018 magnitude 7.5 Palu earthquake. *Nat. Geosci.* **12**, 200–205. (doi:10.1038/s41561-018-0297-z)

40. Socquet A, Hollingsworth J, Pathier E, Bouchon M. 2019 Evidence of supershear during the 2018 magnitude 7.5 Palu earthquake from space geodesy. *Nat. Geosci.* **12**, 192–199. (doi:10.1038/s41561-018-0296-0)

41. Bouchon M, Karabulut H. 2008 The aftershock signature of supershear earthquakes. *Science* **320**, 1323–1325. (doi:10.1126/science.1155030)

42. Walker KT, Shearer PM. 2009 Illuminating the near-sonic rupture velocities of the intracontinental Kokoxili $M_w$ 7.8 and Denali fault $M_w$ 7.9 strike-slip earthquakes with global P wave back projection imaging. *J. Geophys. Res.* **114**, B02304. (doi:10.1029/2008JB005738)

43. Dunham EM. 2007 Conditions governing the occurrence of supershear ruptures under slip-weakening friction. *J. Geophys. Res.* **112**, B07302. (doi:10.1029/2006JB004717)

44. Dunham EM, Favreau P, Carlson JM. 2003 A supershear transition mechanism for cracks. *Science (New York, N.Y.)* **299**, 1557–1559. (doi:10.1126/science.1080650)

45. Liu T, Lapusta N. 2008 Transition of mode II cracks from sub-Rayleigh to intersonic speeds in the presence of favorable heterogeneity. *J. Mech. Phys. Solids* **56**, 25–50. (doi:10.1016/j.jmps.2007.06.005)

46. Shi Z, Day SM. 2013 Rupture dynamics and ground motion from 3-D rough-fault simulations. *J. Geophys. Res.: Solid Earth* **118**, 1122–1141. (doi:10.1002/jgrb.50094)

47. Bruhat L, Fang Z, Dunham EM. 2016 Rupture complexity and the supershear transition on rough faults. *J. Geophys. Res.: Solid Earth* **121**, 210–224. (doi:10.1002/jgrb.v121.1)

48. Hu F, Xu J, Zhang Z, Chen X. 2016 Supershear transition mechanism induced by step over geometry. *J. Geophys. Res.: Solid Earth* **121**, 8738–8749. (doi:10.1002/2016JB013333)

49. Kaneko Y, Lapusta N. 2010 Supershear transition due to a free surface in 3-D simulations of spontaneous dynamic rupture on vertical strike-slip faults. *Tectonophysics* **493**, 272–284. (doi:10.1016/j.tecto.2010.06.015)

50. Hu F, Oglesby DD, Chen X. 2019 The sustainability of free-surface-induced supershear rupture on strike-slip faults. *Geophys. Res. Lett.* **46**, 9537–9543. (doi:10.1029/2019GL084318)

51. Mello M. 2012 Identifying the unique ground motion signatures of supershear earthquakes: theory and experiments. PhD thesis, California Institute of Technology.

52. Freund LB. 1990 *Dynamic fracture mechanics*. Cambridge Monographs on Mechanics. Cambridge, UK: Cambridge University Press.

53. Freund LB. 1972 Crack propagation in an elastic solid subjected to general loading—II. Non-uniform rate of extension. *J. Mech. Phys. Solids* **20**, 141–152. (doi:10.1016/0022-5096(72)90007-5)

54. Svetlizky I, Fineberg J. 2014 Classical shear cracks drive the onset of dry frictional motion. *Nature* **509**, 205–208. (doi:10.1038/nature13202)

55. Drucker DC, Prager W. 1952 Soil mechanics and plastic analysis or limit design. *Q. Appl. Math.* **10**, 157–165. (doi:10.1090/qam/1952-10-02)

56. Gabriel AA, Ampuero JP, Dalguer LA, Mai PM. 2013 Source properties of dynamic rupture pulses with off-fault plasticity. *J. Geophys. Res.: Solid Earth* **118**, 4117–4126. (doi:10.1002/jgrb.50213)

57. Rougier EE, Knight EE, Lei Z, Munjiza A. 2016 OSS.edu2.x (Hybrid Optimization Software Suite—257. Educational Version, Second Generation).

58. Knight EE, Rougier E, Lei Z, Euser B, Chau V, Boyce SH, Gao K, Okubo K, Froment M. 2020 HOSS: an implementation of the combined finite-discrete element method. *Comput. Part. Mech.* **7**, 765–787. (doi:10.1007/s40571-020-00349-y)

59. Munjiza A. 2004 *The combined finite-discrete element method*. Chichester, UK: John Wiley & Sons, Ltd.

60. Walsh JB. 1965 The effect of cracks in rocks on Poisson's ratio. *J. Geophys. Res.* **70**, 5249–5257. (doi:10.1029/JZ070i020p05249)

61. Walsh JB. 1965 The effect of cracks on the compressibility of rock. *J. Geophys. Res.* **70**, 381–389. (doi:10.1029/JZ070i002p00381)

62. Faulkner DR, Mitchell TM, Healy D, Heap MJ. 2006 Slip on 'weak' faults by the rotation of regional stress in the fracture damage zone. *Nature* **444**, 922–925. (doi:10.1038/nature05353)
63. Hiramatsu Y, Honma H, Saiga A, Furumoto M, Ooida T. 2005 Seismological evidence on characteristic time of crack healing in the shallow crust. *Geophys. Res. Lett.* **32**, L09304. (doi:10.1029/2005GL022657)

64. Froment B, McGuire JJ, van der Hilst RD, Gouédard P, Roland EC, Zhang H, Collins JA. 2014 Imaging along-strike variations in mechanical properties of the Gofar transform fault, East Pacific Rise. *J. Geophys. Res.: Solid Earth* **119**, 7175–7194. (doi:10.1002/2014JB011270)

65. Ashby MF, Sammis CG. 1990 The damage mechanics of brittle solids in compression. *Pure Appl. Geophys.* **133**, 489–521. (doi:10.1007/BF00878002)

66. Deshpande V, Evans A. 2008 Inelastic deformation and energy dissipation in ceramics: a mechanism-based constitutive model. *J. Mech. Phys. Solids* **56**, 3077–3100. (doi:10.1016/j.jmps.2008.05.002)

67. Bhat HS, Rosakis AJ, Sammis CG. 2012 A micromechanics based constitutive model for brittle failure at high strain rates. *J. Appl. Mech.* **79**, 031016. (doi:10.1115/1.4005897)

68. Rice JR. 1975 Continuum mechanics and thermodynamics of plasticity in relation to microscale deformation mechanisms. In *Constitutive equations in plasticity* (ed. A Argon). Cambridge, MA: MIT Press.

69. Zhang QB, Zhao J. 2013 Effect of loading rate on fracture toughness and failure micromechanisms in marble. *Eng. Fract. Mech.* **102**, 288–309. (doi:10.1016/j.engfracmech.2013.02.009)

70. Palmer AC, Rice JR. 1973 The growth of slip surfaces in the progressive failure of over-consolidated clay. *Proc. R. Soc. Lond. A* **332**, 527–548. (doi:10.1098/rspa.1973.0040)

71. Byerlee J. 1978 Friction of rocks. *Pure Appl. Geophys.* **116**, 615–626. (doi:10.1007/BF00876528)

72. Wibberley CAJ, Yielding G, Di Toro G. 2008 Recent advances in the understanding of fault zone internal structure: a review. *Geol. Soc. Lond. Spec. Publ.* **299**, 5–33. (doi:10.1144/SP299.2)

73. Chen WF, Han DJ. 2007 *Plasticity for structural engineers*. New York, NY: Springer US.

74. Davis R, Selvadurai A. 2005 *Plasticity and geomechanic*. Cambridge, UK: Cambridge University Press.

75. Vallage A, Klinger Y, Grandin R, Bhat HS, Pierrot-Deseilligny M. 2015 Inelastic surface deformation during the 2013 Mw 7.7 Balochistan, Pakistan, earthquake. *Geology* **43**, G37290.1. (doi:10.1130/G37290.1)

76. Klinger Y et al. 2018 Earthquake damage patterns resolve complex rupture processes. *Geophys. Res. Lett.* **45**, 210–279. (doi:10.1029/2018GL078842)

77. Barnhart WD, Gold RD, Hollingsworth J. 2020 Localized fault-zone dilatancy and surface inelasticity of the 2019 Ridgecrest earthquakes. *Nat. Geosci.* **13**, 699–704. (doi:10.1038/s41561-020-0628-8)

78. Delorme A, Grandin R, Klinger Y, Pierrot-Deseilligny M, Feuillet N, Jacques E, Rupnik E, Morishita Y. 2020 Complex deformation at shallow depth during the 30 October 2016 Mw 6.5 Norcia earthquake: interference between tectonic and gravity processes? *Tectonics* **39**, e2019TC005596. (doi:10.1029/2019TC005596)

79. Ajorlou N, Hollingsworth J, Mousavi Z, Ghods A, Masoumi Z. 2021 Characterizing near-field surface deformation in the 1990 Rudbar earthquake (Iran) using optical image correlation. *Geochem. Geophys. Geosyst.* **22**, e2021GC009704. (doi:10.1029/2021GC009704)

80. Milliner CWD, Dolan JF, Hollingsworth J, Leprince S, Ayoub F, Sammis CG. 2015 Quantifying near-field and off-fault deformation patterns of the 1992 M w 7.3 Landers earthquake. *Geochem. Geophys. Geosyst.* **16**, 1577–1598. (doi:10.1002/2014GC005693)

81. Rockwell TK, Lindvall S, Dawson T, Langridge R, Lettis W, Klinger Y. 2002 Lateral offsets on surveyed cultural features resulting from the 1999 Izmit and Duzce earthquakes, Turkey. *Bull. Seismol. Soc. Am.* **92**, 79–94. (doi:10.1785/0120000809)

82. Mitchell TM, Faulkner DR. 2009 The nature and origin of off-fault damage surrounding strike-slip fault zones with a wide range of displacements: a field study from the Atacama fault system, northern Chile. *J. Struct. Geol.* **31**, 802–816. (doi:10.1016/j.jsg.2009.05.002)

83. Rupnik E, Daakir M, Deseilligny MP. 2017 MicMac—a free, open-source solution for photogrammetry. *Open Geospat. Data Softw. Stand.* **2**, 14. (doi:10.1186/s40965-017-0027-2)

84. Klinger Y, Xu X, Tapponnier P, Woerd JVD, Lasserre C, King G. 2005 High-resolution satellite imagery mapping of the surface rupture and slip distribution of the Mw 7.8, 14 November 2001 Kokoxili earthquake, Kunlun fault, Northern Tibet, China. *Bull. Seismol. Soc. Am.* **95**, 1970–1987. (doi:10.1785/0120040233)
85. Klinger Y, Michel R, King G. 2006 Evidence for an earthquake barrier model from $M_w$ 7.8 Kokoxili (Tibet) earthquake slip-distribution. *Earth Planet. Sci. Lett.* **242**, 354–364. (doi:10.1016/j.epsl.2005.12.003)

86. Rosu AM, Pierrot-Deseilligny M, Delorme A, Binet R, Klinger Y. 2015 Measurement of ground displacement from optical satellite image correlation using the free open-source software MicMac. *ISPRS J. Photogramm. Remote Sens.* **100**, 48–59. (doi:10.1016/j.isprsjprs.2014.03.002)

87. King G. 2005 Slip-partitioned surface breaks for the $M_w$ 7.8 2001 Kokoxili earthquake, China. *Bull. Seismol. Soc. Am.* **95**, 731–738. (doi:10.1785/0120040101)

88. Ratchkovski NA et al. 2003 Aftershock sequence of the $M_w$ 7.9 Denali fault, Alaska, earthquake of 3 November 2002 from regional seismic network data. *Seismol. Res. Lett.* **74**, 743–752. (doi:10.1785/gssrl.74.6.743)

89. Holtkamp S, Ruppert N. 2015 A high resolution aftershock catalog of the magnitude 7.5 Craig, Alaska, earthquake on 5 January 2013. *Bull. Seismol. Soc. Am.* **105**, 1143–1152. (doi:10.1785/0120140179)

90. Udias A, Madariaga R, Buforn E. 2014 *Source mechanisms of earthquakes*. Cambridge, UK: Cambridge University Press.

91. Hsu YJ, Simons M, Avouac JP, Galeteka J, Sieh K, Chlieh M, Natawidjaja D, Prawirodirdjo L, Bock Y. 2006 Frictional afterslip following the 2005 Nias-Simeulue earthquake, Sumatra. *Science* **312**, 1921–1926. (doi:10.1126/science.1126960)

92. Ross ZE, Rollins C, Cochran ES, Hauksson E, Avouac JP, Ben-Zion Y. 2017 Aftershocks driven by afterslip and fluid pressure sweeping through a fault-fracture mesh. *Geophys. Res. Lett.* **44**, 8260–8267. (doi:10.1002/2017GL074634)

93. Freed AM, Bürgmann R, Calais E, Freymueller J, Hreinsdóttir S. 2006 Implications of deformation following the 2002 Denali, Alaska, earthquake for postseismic relaxation processes and lithospheric rheology. *J. Geophys. Res.* **111**, B01401. (doi:10.1029/2005JB003894)

94. Hearn EH, McClusky S, Ergintav S, Reilinger RE. 2009 Izmit earthquake postseismic deformation and dynamics of the North Anatolian fault zone. *J. Geophys. Res.* **114**, B08405. (doi:10.1029/2008JB006026)

95. Ding K, Freymueller JT, Wang Q, Zou R. 2015 Coseismic and early postseismic deformation of the 5 January 2013 $M_w$ 7.5 Craig earthquake from static and kinematic GPS solutions. *Bull. Seismol. Soc. Am.* **105**, 1153–1164. (doi:10.1785/0120140172)

96. Bouchon M. 2002 Space and time evolution of rupture and faulting during the 1999 Izmit (Turkey) earthquake. *Bull. Seismol. Soc. Am.* **92**, 256–266. (doi:10.1785/0120000845)

97. Zhang Y, Chen YT, Feng W. 2019 Complex multiple-segment ruptures of the 28 September 2018, Sulawesi, Indonesia, earthquake. *Sci. Bull.* **64**, 650–652. (doi:10.1016/j.scib.2019.04.018)

98. Li Q, Zhao B, Tan K, Xu W. 2020 Two main rupture stages during the 2018 magnitude 7.5 Sulawesi earthquake. *Geophys. J. Int.* **221**, 1873–1882. (doi:10.1093/gji/ggaa115)