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

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

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 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 7.4 Izmit and the 1999 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 7.8 Kunlun (China) earthquake [29-31], the 2002 7.8 Denali (Alaska) earthquake [32,33], the 2010 6.9 Qinghai (China) earthquake [34], the 2012 8.6 off-Sumatra (Indonesia) earthquake [35], the 2013 7.5 Craig (Alaska) earthquake [36], the 2013 6.7 Okhotsk (Kamtchatka) earthquake [37,38] and most recently the 2018 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 ratio: , where is the background shear stress and and are the peak and residual frictional strengths [15-17]. Supershear ruptures occur if in two dimensions [15,17], or 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 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.

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 , where is the Rayleigh wave speed. Let the origin of the polar coordinate system coincide with the crack tip. The near-tip stress, which also depends on the rupture speed , is given as [52], and where is the initial stress state and is the dynamic stress intensity factor, which can be approximated as where , the -wave speed, is the limiting speed for a mode II crack and is the current crack length and is the Rayleigh wave speed [52]. Here is the static stress intensity factor of the crack and is the stress drop defined as . Here is the initial shear stress on the fault and 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, [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, 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. , where Here is the hydrostatic stress derived from the first invariant of the stress tensor, , corresponds to the second invariant of the deviatoric stress tensor () and 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 . The closed-form expressions are quite long, but the key variables that affect the solution can be written as Thus, given a rupture velocity history, , 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, , 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 (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.

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, . (c) Freund’s universal function. Here and are the Rayleigh and -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, , at key points labelled 1 through 9 in (d). (f) Superposition of the contours where throughout the history of the rupture. (Online version in colour.)

Table 1

Parameters for the LEFM, micromechanical and FDEM models. *Micromechanical model. †FDEM model.

parameters	description	values
state of stress
−σxx0	fault-parallel stress (MPa)	37.7
−σyy0	normal stress on the fault (MPa)	60.7
σyx0	shear stress on the fault (MPa)	19.9
	nucleation zone* (MPa)	36.4
ψ	angle between σ_1 and fault (°)	60
bulk properties
ρ	density (kg m^−3)	2700
c_s	S-wave speed (km s^−1)	3.12
c_p	P-wave speed (km s^−1)	5.62
ν	Poisson’s ratio	0.276
fault parameters
f_s	static friction coefficient	0.6
	nucleation zone†	0.3
f_d	dynamic friction coefficient	0.1
D_c	characteristic slip distance (m)	1.0
R_0	process zone size (m)	1057
S	S-ratio	1.2

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.

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.

Modelling strategies

The two numerical models presented in this study account for fracture damage-dominated brittle rheology. Earthquake ruptures induce large dynamic strain rates () around the rupture tip that cannot be fully accounted for by classical Mohr–Coulomb or Drucker–Prager plasticity, typically used at low strain rates (). 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 - and -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 of a solid, without flaws, deforming purely elastically and (ii) the free energy corresponding to the contribution of the current set of microcracks. Using thermodynamic arguments, 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 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].

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 , and the minimum compressive stress are in the – plane, whereas the intermediate principal stress coincides with . The fault plane makes an angle of with with uniform normal traction () and shear traction () 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 . 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 -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 () so that the models could be compared with natural or laboratory earthquakes. is defined as Here is Poisson’s ratio, is shear modulus, is the fracture energy associated with the slip-weakening law, and 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.

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. (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 . 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 -wave field concentrates, although the resulting damage density is much smaller than behind the rupture front. The location of newly formed cracks corresponds to the region where is positive. The transition to supershear velocity (figure 2b) directly impacts the generation of coseismic damage. Around 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, becomes small, even if the total length of the fault that has ruptured ( in equation (2.3)), has increased significantly. Moreover, a snapshot of at (figure 2b) shows that the first pulse is not strong enough to change the state of stress optimally, thus reducing even further . These two factors likely explain why damage mostly occurs behind the front of the pulse propagating at sub-Rayleigh rupture speed. After , 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.

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 ), using the micromechanical model. The rupture is nucleated around 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, (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 and the compressional and the tensional quadrants, respectively. (Online version in colour.)

Figure 3

Slip and slip rate evolution with the micromechanical model. Cumulative slip (a), and slip rate, (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.)

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

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, (see eqn (4.142) of Chen & Han [73]): where is the cohesion and is given by where where is the third invariant of the deviatoric stress tensor. Failure occurs when . Note that when and , becomes the Drucker–Prager criterion as described in equation (2.4).

Figure 4 displays the 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 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.

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

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.

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.

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 7.8 Kunlun (China) earthquake is a strike-slip event with approximately 400 km long surface rupture and a mean rupture speed between and , 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 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.

Figure 5

Optical image correlation analysis. (a) Map of the strike–slip section of the 2001 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.)

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 , 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 over 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.

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 7.4 Izmit earthquake (Turkey) had a mean rupture velocity of about [27,28], and a well-recorded three-month aftershock catalogue [41] (figure 6a for one-month aftershock locations). In the case of the 2002 7.9 Denali earthquake (Alaska), rupture velocities of about were reported [32], and 1-year long catalogue of aftershocks is available [88] (see figure 6b for one-month aftershock locations). The 2013 7.5 Craig earthquake (Alaska) was also inferred as a supershear event, with rupture velocities between 5.5 and [36] and a five-month aftershock catalogue is available [89] (see figure 6c for one-month aftershock locations).

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

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 of the fault trace (see the electronic supplementary material, for the same analysis considering ; electronic supplementary material, figure S6). We assume that each aftershock, having a seismic moment , is a circular crack and compute the slip distribution using [90]: where is the crack radius, is the shear modulus and 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 , Denali and Craig ). 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 th bin containing aftershocks is given by Note that is part of the slip distribution of the th aftershock projected on the th 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 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,d,f, 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.

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

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 , assuming that the location uncertainty of aftershocks is normally distributed. Here, corresponds to the aftershock epicentre, while 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 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).

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

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

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