Signature of transition to supershear rupture speed in coseismic off-fault damage zone
Jorge Jara, Lucile Bruhat, Marion Y. Thomas, Solène Antoine, Kurama Okubo, Esteban Rougier, Ares J. Rosakis, Charles G. Sammis, Yann Klinger, Romain Jolivet, Harsha S. Bhat
SSignature of supershear transition seen in damage and aftershock pattern
Jorge Jara , Lucile Bruhat , Sol`ene Antoine , Kurama Okubo , Marion Y. Thomas , Esteban Rougier , Ares J.Rosakis , Charles G. Sammis , Yann Klinger , Romain Jolivet and Harsha S. Bhat .
1. Laboratoire de G´eologie, D´epartement de G´eosciences, ´Ecole Normale Sup´erieure, CNRS, UMR 8538, PSL ResearchUniversity, Paris, France2. Institut de Physique du Globe de Paris, Sorbonne Paris Cit´e, Universit´e Paris Diderot, CNRS, UMR 7154, Paris, France3. Department of Earth and Planetary Sciences, Harvard University, Cambridge, MA 02138, USA4. Institut des Sciences de la Terre de Paris, Sorbonne Universit´e, CNRS, UMR 7193, Paris, France5. EES-17 Earth and Environmental Sciences Division, Los Alamos National Laboratory, Los Alamos, NM, USA6. Graduate Aerospace Laboratories, California Institute of Technology, Pasadena, California, 91125, USA7. Departament of Earth Sciences, University of Southern California, Los Angeles, CA 90089, USA
Supershear earthquakes are rare but powerful ruptures with devastating consequences. How quickly an earth-quake rupture attains this speed, or for that matter decelerates from it, strongly affects high-frequency groundmotion and the spatial extent of coseismic off-fault damage . Traditionally, studies of supershear earthquakeshave focused on determining which fault segments sustained fully-grown supershear ruptures . Knowing thatthe rupture first propagated at subshear rupture speeds, these studies usually guessed an approximate locationfor the transition from subshear to supershear regimes. The rarity of confirmed supershear ruptures ,combined with the fact that conditions for supershear transition are still debated , complicates the inves-tigation of supershear transition in real earthquakes. Here, we find a unique signature of the location of asupershear transition: we show that, when a rupture accelerates towards supershear speed, the stress concen-tration abruptly shrinks, limiting the off-fault damage and aftershock productivity. First, we use theoreticalfracture mechanics to demonstrate that, before transitioning to supershear, the stress concentration aroundthe rupture tip shrinks, confining the region where damage & aftershocks are expected. Then, employing twodifferent dynamic rupture modeling approaches, we confirm such reduction in stress concentration, furthervalidating the expected signature in the transition region. We contrast these numerical and theoretical resultswith high-resolution aftershock catalogs for three natural supershear earthquakes, where we identify a smallregion with lower aftershock density near the supershear transition. Finally, using satellite optical image cor-relation techniques, we show that, for a fourth event, the transition zone is characterized by a diminution in thewidth of the damage zone. Our results demonstrate that the transition from subshear to supershear rupturecan be clearly identified by a localized absence of aftershocks, and a decrease in off-fault damage, due to atransient reduction of the stress intensity at the rupture tip. a r X i v : . [ phy s i c s . g e o - ph ] S e p hile the rupture speed of most earthquakes is limited by the speed of the Rayleigh and shear waves, earthquakescan occassionaly rupture at higher speeds. These are known as supershear earthquakes. Because the on-fault rupturespeed controls the characteristics of the resulting off-fault damage, understanding the conditions for a rupture totransition from the sub-Rayleigh to the supershear regime is critical for earthquake hazard assessment. Abrupt changesin the rupture speed indeed increases high-frequency radiation , and governs the spatial extent of the coseismic off-fault damage zone . Whether supershear ruptures are real or not has been a matter of debate for a long time. Eventhough theoretical models and laboratory experiments provided evidence that supershear ruptures exist sincethe early 1970s, it was not until the M w . Pioneering laboratory experiments together with observations from the M w M w , then conclusively confirmed that supershearruptures are in fact much more common than previously expected. Supershear ruptures have now been inferred forseveral, albeit rare, events: The M w , the M w , the M w , the M w and most recently the M w . Efforts have been made in order to identify andcharacterize these supershear earthquakes using various methodological approaches , where the methods used werein most cases, designed to reveal a posteriori which segment of the rupture propagated at supershear velocities. Thesemethods often ignore the details of the transition to supershear rupture speed. They simply limit themselves to theanalysis of full-grown supershear ruptures. When the rupture was known to propagate earlier at subshear speeds, thesupershear transition is presumed to have occurred in-between, leading to an imprecise location of this transition. Theconditions for supershear transitions in nature are still poorly understood, despite numerical efforts to characterizethe mechanics of the transition process . Here, we present an original signature of a supershear transition bycombining theoretical and numerical modeling with field observations (high-resolution aftershock catalogs and surfacedisplacement fields from optical satellite image correlation analyses).The theoretical analysis is conducted using Linear Elastic Fracture Mechanics which provides closed form so-lutions to describe the state of stress around a rupture tip. Considering a semi-infinite plain-strain crack in a 2Dhomogeneous isotropic linear medium, moving at a speed v < c R (where c R is the Rayleigh wave speed), the near-tipstress solution is controlled by the dynamic stress intensity factor, K dynII , which evolves with the rupture speed v .A Lorentz-like contraction of the stress field occurs around the rupture tip as the rupture speed approaches its limit-ing speed, c R . This contraction has also been observed and verified experimentally . For an earthquake rupturetransitioning to a supershear regime, the rupture has to first accelerate to the Rayleigh wave speed. As the rupture ap-proaches the supershear transition, K dynII monotonically decreases to zero, strongly reducing the stress concentrationat the rupture tip. Thus the off-fault domain affected by this stress concentration will also shrink (see Methods section2or details). We illustrate this effect by calculating the extent of the region where the stress state goes beyond the limitsdefined by a Drucker-Prager failure criterion (Figs. 1a and 1b). This domain, which describes the theoretical extent ofcoseismic off-fault deformation, such as fracture damage and aftershocks, will be directly affected by stress changes.As the rupture approaches the Rayleigh wave speed, with uniform (Fig. 1a) or non-uniform (Fig. 1b) rupture velocity,the stress concentration at the rupture tip eventually collapses, limiting the spatial extent of possible off-fault damageand aftershocks during the transition from the sub-Rayleigh to the supershear regime.This theoretical demonstration is then validated through two different numerical models that account for dynamicevolution of coseismic damage . In these models, unlike the theoretical development presented above, the ruptureis spontaneous and there is a feedback between off-fault damage and on-fault rupture. Both models produce in-planedynamic simulations of an earthquake rupture on a 1D right-lateral planar fault embedded in a 2D medium. A slip-weakening friction law is used to model the earthquake rupture and damage only occurs on the tensional part of the fault(bottom side of the fault on Fig. 1c and d) . The first model employs enhanced numerical algorithms for earthquakerupture allowing for spontaneous activation of off-fault fracture networks . This has been observed in experiments aswell . Once the rupture is over, we evaluate the damage pattern and the spatial variation of closeness to failure ( ∆ CF )resulting from the rupture. Regions with positive values of ∆ CF are more likely to host fracture damage and triggerfuture aftershocks. We see that the regions experiencing sub-Rayleigh and full blown supershear regimes manifesthigh local stress fluctuations (Fig. 1c). However, when the rupture is transitioning to supershear speed, ∆ CF ∼ ,and the spatial extent of the off-fault damage zone drops dramatically. The second numerical model accounts foroff-fault damage using a micro-mechanics based effective constitutive law . The resulting off-fault damage density,which is directly affected by the near-tip stress state, as illustrated by the ∆ CF , is also characterized by a suddenshrinkage of the damage zone (Fig. 1d) during the supershear transition (see Methods section for details).Assuming that the nucleation of early near-fault aftershocks (occurring around a week after the main earthquake)is mainly governed by the stress state left in the wake of the earthquake, we analyze the spatial and temporal distribu-tion of relocated aftershocks for 3 well-known supershear ruptures: The M w (Fig. 2a), the M w (Fig. 2d), and the M w (Fig. 2g). For these examples, we compute the seismic moment released by the aftershocks at a distance of less than5 km from the main fault, and, for different periods of time after the main shock: 1-3 days, 1, 2, 3 weeks and 1 month(Figs. 2 c, f, and i). In each case, we observe a region of limited extent ( ∼ . It suggests that the gap in the early aftershocks productivity (less than3 weeks after the mainshock) is mainly related to the supershear transition (see Methods section for details).As high-resolution aftershocks catalogs are not available for all the supershear earthquake ruptures, we explore anew method to investigate the spatial evolution of the damage zone width. Recent developments in satellite opticalimage analysis and sub-pixel correlation methods now allow for the detection of displacement variations due to anearthquake down to a resolution of one meter. Noting that, for earthquakes of large magnitudes ( M w ≥ ), faultzones (fault core and damage zone) are usually of metric- to kilometric-scale , we apply this method to the M w . Displacement profiles normal to the fault are computed using pre- and post-earthquake images and are used to infer the width of the fault zone (fault core and damage zone Fig. 3). While theoverall mean of the fault zone width is around 238 meters, we observe a clear localized region, ∼
11 km-long, wheredamage largely reduced down to a mean of 127 meters, located around what was previously inferred as the supersheartransition zone . As expected from numerical modeling, the supershear transition is characterized by a significantreduction in damage zone width (see Methods section for details). We acknowledge here that the aftershock catalogof Robinson and colleagues also alludes to the same conclusion. However due to the lack of high spatio-temporaldensity of aftershocks in their catalog we instead chose the above technique.Theoretical and numerical modeling of off-fault damage both show that the transition from the sub-Rayleigh tothe supershear regime is characterized by a reduction in the width of the damage zone and a paucity of aftershocks.This results from the transient shrinkage in stress concentration around the rupture tip as its velocity approaches theRayleigh wave speed. We cross-validate this predicted feature, with seismological evidence of aftershock quiescencefor the Izmit, Denali and Craig supershear earthquakes, and with geodetic observations of the damage zone for theKunlun earthquake. These results are valid for a well-developed sub-Rayleigh rupture that transitions to supershearspeeds. Recent observations from the Palu earthquake hint that the rupture might have either nucleated directly at asupershear speed, or transitioned very early on . For this particular case, further exploration is required to see ifsimilar features can be observed in the field. In conclusion, identifying an absence of aftershocks and decrease in off-fault damage allows us to pinpoint the transition from sub-Rayleigh to supershear speeds. This work provides a newframework to better characterize supershear transition zones in the field, and further explore the mechanical conditionsfor such transitions. 4 eferences [1] Das, S. The Need to Study Speed. Science , 905–906 (2007).[2] Bhat, H. S., Dmowska, R., King, G. C., Klinger, Y. & Rice, J. R. Off-fault damage patterns due to supershearruptures with application to the 2001 Mw 8.1 Kokoxili (Kunlun) Tibet earthquake.
J. Geophys. Res. , 1–19(2007).[3] Thomas, M. Y. & Bhat, H. S. Dynamic evolution of off-fault medium during an earthquake: A micromechanicsbased model.
Geophys. J. Int. , 1267–1280 (2018).[4] Okubo, K. et al.
Dynamics, radiation and overall energy budget of earthquake rupture with coseismic off-faultdamage.
J. Geophys. Res. , 1–41 (2019).[5] Vall´ee, M., Land`es, M., Shapiro, N. M. & Klinger, Y. The 14 November 2001 Kokoxili (Tibet) earthquake:High-frequency seismic radiation originating from the transitions between sub-Rayleigh and supershear rupturevelocity regimes.
J. Geophys. Res. , 1–14 (2008).[6] Ellsworth, W. L. et al.
Near-Field Ground Motion of the 2002 Denali Fault, Alaska, Earthquake Recorded atPump Station 10.
Earthq. Spectra , 597–615 (2004).[7] Bao, H. et al. Early and persistent supershear rupture of the 2018 magnitude 7.5 Palu earthquake.
NatureGeoscience , 200–205 (2019).[8] Socquet, A., Hollingsworth, J., Pathier, E. & Bouchon, M. Evidence of supershear during the 2018 magnitude7.5 Palu earthquake from space geodesy. Nature Geoscience , 192–199 (2019).[9] Bouchon, M. & Karabulut, H. The aftershock signature of supershear earthquakes. Science , 1323–1325(2008).[10] Archuleta, R. J. A faulting model for the 1979 Imperial Valley earthquake.
J. Geophys. Res. , 4559–4585(1984).[11] Bouchon, M. et al. How fast is rupture during an earthquake? New insights from the 1999 Turkey earthquakes.
Geophys. Res. Lett. , 2723–2726 (2001).[12] Robinson, D. P., Brough, C. & Das, S. The Mw 7.8, 2001 Kunlunshan earthquake: Extreme rupture speedvariability and effect of fault geometry. J. Geophys. Res. , 1–10 (2006).[13] Yue, H. et al.
Supershear rupture of the 5 January 2013 Craig, Alaska (Mw7.5) earthquake.
J. Geophys. Res. , 5903–5919 (2013). 514] Zhan, Z., Helmberger, D. V., Kanamori, H. & Shearer, P. M. Supershear rupture in a Mw 6.7 aftershock of the2013 Sea of Okhotsk earthquake.
Science , 204–207 (2014).[15] Bruhat, L., Fang, Z. & Dunham, E. M. Rupture complexity and the supershear transition on rough faults.
J.Geophys. Res. , 210–224 (2016).[16] Burridge, R. Admissible Speeds for Plane-Strain Self-Similar Shear Cracks with Friction but Lacking Cohesion.
Geophys. J. Int. , 439–455 (1973).[17] Andrews, D. J. Rupture velocity of plane strain shear cracks. J. Geophys. Res. , 5679–5687 (1976).[18] Wu, F. T., Kuenzler, H. & Thomson, K. C. Stick-slip propagation velocity and seismic source mechanism. Bull.Seism. Soc. Am. , 1621–1628 (1972).[19] Rosakis, A. J., Samudrala, O. & Coker, D. Cracks Faster than the Shear Wave Speed. Science , 1337–1340(1999).[20] Xia, K. W., Rosakis, A. J. & Kanamori, H. Laboratory Earthquakes: The Sub-Rayleigh-to-Supershear RuptureTransition.
Science , 1859–1861 (2004).[21] Passelegue, F. X., Schubnel, A., Nielsen, S., Bhat, H. S. & Madariaga, R. From Sub-Rayleigh to SupershearRuptures During Stick-Slip Experiments on Crustal Rocks.
Science , 1208–1211 (2013).[22] Liu, Y. & Lapusta, N. Transition of mode II cracks from sub-Rayleigh to intersonic speeds in the presence offavorable heterogeneity.
J. Mech. Phys. Solids , 25–50 (2008).[23] Liu, C., Bizzarri, A. & Das, S. Progression of spontaneous in-plane shear faults from sub-rayleigh to compres-sional wave rupture speeds. J. Geophys. Res. , 8331–8345 (2014).[24] Freund, L. B. Crack propagation in an elastic solid subjected to general loading-ii. non-uniform rate of extension.
J. Mech. Phys. Solids , 141–152 (1972).[25] Svetlizky, I. Fineberg, J. Classical shear cracks drive the onset of dry frictional motion. Nature , 205–208(2014).[26] Rosakis, A. J. Intersonic shear cracks and fault ruptures.
Adv. Phys. , 1189–1257 (2002).[27] Freed, A. M., B¨urgmann, R., Calais, E., Freymueller, J. & Hreinsd´ottir, S. Implications of deformation followingthe 2002 Denali, Alaska, earthquake for postseismic relaxation processes and lithospheric rheology. J. Geophys.Res. , B01401 (2006). 628] Hearn, E. H., McClusky, S., Ergintav, S. & Reilinger, R. E. Izmit earthquake postseismic deformation anddynamics of the North Anatolian Fault Zone.
J. Geophys. Res. , B08405 (2009).[29] Ding, K., Freymueller, J. T., Wang, Q. & Zou, R. Coseismic and Early Postseismic Deformation of the 5 January2013 Mw 7.5 Craig Earthquake from Static and Kinematic GPS Solutions.
Bull. Seism. Soc. Am. , 1153–1164(2015).[30] Gold, R. D. et al.
On- and off-fault deformation associated with the September 2013 Mw 7.7 Balochistanearthquake: Implications for geologic slip rate measurements.
Tectonophysics , 65–78 (2015).7 upplementary Information is available in the online version of the paper.
Acknowledgements
J.J. and R.J. acknowledge the funding of the European Research Council (ERC) under the Euro-pean Union’s Horizon 2020 research and innovation program (grant agreement 758210, project Geo4D). L.B. thanksthe funding from the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Pro-gramme (FP7/2007-2013) under REA grant agreement PCOFUND-GA-GA-2013-609102, through the PRESTIGEprogramme coordinated by Campus France. S.A. and Y.K. are partly supported by the ANR project DISRUPT (ANR-18-CE31-0012).
Author Contributions
H.S.B. conceived, designed and supervised the project. J.J. performed the aftershock analysis.L.B. and H.S.B. conducted the theoretical study. S.A. did the image correlation. K.O., E.R., M.T. and H.S.B. con-ducted the numerical analysis. J.J., L.B. and H.S.B. wrote the manuscript. All the authors contributed to the analysis,interpretation and manuscript preparation.
Author Information igure 1: Theoretical and Numerical Evidence of Supershear Transition Signature . a. & b. Domain that exceedsDrucker-Prager failure criterion for a crack of uniform rupture velocity (a) and for a crack of non-uniform rupturevelocity (b) . Both cases show that the region of potential damage shrinks as the rupture velocity v approaches theRayleigh wave speed c R . c. Changes in closeness to failure associated with Drucker-Prager failure criterion, computedusing 2D FDEM dynamic simulation including coseismic off-fault damage generated by the rupture propagation . d .Same as (c) but using an effective medium theory . The grey regions map the spatial distribution of damage densitythat occurs during the entire event while the black field record the damage that takes place at the time at which theDrucker-Prager failure criterion was computed. In both cases ( c & d ), regions of positive change are more likely tohost fracture damage and trigger future aftershocks. The transition from the sub-Rayleigh to the supershear regime isclearly denoted by a more localized and weaker stress perturbation and a consequent gap in damage.10 igure 2 igure 2: High-Resolution Aftershock Catalog Analysis . 1-month aftershock distribution for Izmit ( a ), Denali ( d ),and Craig ( g ) earthquakes, color-coded by time and indicating the respective event epicenter (color-coded star) andfocal mechanism. The black continuous line denotes the surface rupture for each event. c,f,i Cumulative aftershockseismic moment release projected on the main fault (in log scale) at different temporal scales (1-3 days, 1-2-3 weeksand 1 month), for Izmit ( c ), Denali ( f ), and Craig ( i ) earthquakes. All the aftershocks within a distance of 5 kmfrom the fault are considered in the calculation (area denoted by the black discontinuous lines in a , d and g ). Color-coded arrows indicate the different speed regimes reported for each event (green for sub-Rayleigh and orange forsupershear) , while the pink boxes indicate this work’s proposed Transition Zones. b,e,h Zoom plot of the regionproposed as Transition zone in this study for each earthquake.12 igure 3: Optical Correlation Images Analysis . a. Map of the strike slip section of the M w . b. Along-Srike fault zone width (black) and its associated uncertainty (grey), obtained from the analysis of 40 km-longprofiles, sampling the fault zone every 500 m, on the surface displacement maps. The latter is derived from correlatingpre- and post-earthquake SPOT 1-4 images. The 11 km-long red area highlighted the specific region with a mean faultwidth (red dashed line) of only 127 m compared to 238 m recorded for the rest of the rupture (red line). The latterexcludes the area where two parallel fault strands are activated and for which the fault zone is exceptionally large ( > c. Zoom of the Fig. 3 b . 13 ethods Linear elastic fracture mechanics solution
The variability of earthquake rupture speed affects the high-frequency content generation and the activation of thecoseismic off-fault damage . We develop an approximation to understand how the velocity grows during a seismicrupture using the Linear Elastic Fracture Mechanics approach. In a 2D homogeneous isotropic linear-elastic body, weconsider a semi-infinite plain-strain crack of length L , and ( r , θ ) polar coordinates centered at the crack tip. The staticnear-tip stress σ αβ field is given by σ αβ ( r, θ ) = K II √ πr f IIαβ ( θ ) (1)where K II is the static stress intensity factor and f IIαβ are universal angular functions . The static stress intensityfactor varies with the crack length such that K II = φ ( σ, τ ) √ πL (2)where φ ( σ, τ ) depends on the applied normal σ and shear τ stress. The stress at the rupture tip increases then with thecrack length.Now consider that the crack tip moves at a speed v ≤ c R , where c R is the Rayleigh wave speed. The near-tipstresses now depend on a the rupture speed v as follows σ αβ ( r, θ, v ) = K dynII ( v ) √ πr f IIαβ ( θ, v ) (3)where K dynII is the dynamic stress intensity factor. This solution is entirely analogous to the static problem. However,due to the moving coordinate system, all the fields undergo a Lorentz-like contraction, affecting both the stress intensityfactor and the angular distribution. The dynamic stress intensity factor can be approximated as K dynII ≈ − v/c R (cid:112) − v/c P K II (4)where c P is the limiting speed for a mode II crack, being c P the P-wave speed. Also K II = φ ( σ, τ ) (cid:112) πL ( t ) wherethe crack length is given by L ( t ) = (cid:82) t v ( t ) dt .Using the expressions above, we will compute the extent of the ‘yield’ region in the off-fault medium. To generalizeand illustrate the state of stress, we represent normal and shear stress by their invariants I = σ kk and J = (cid:114) s ij s ji with s ij = σ ij − I δ ij (5)Solutions from linear elastic fracture mechanics give the state of stress at the crack tip for a mode II crack: σ + σ = K II √ πr sin θ . (6)14ince σ = ν ( σ + σ ) , where ν Poisson’s ratio, we can compute the stress invariant I : I = 1 + ν K II √ πr sin θ . (7)Likewise, we compute s : s = σ (8) = K II √ πr (cid:20) cos θ (cid:18) − sin θ θ (cid:19)(cid:21) (9)The bracketed term is designated as A ( θ ) . The resulting second invariant of the stress deviator J is J = (cid:114) s , (10) = K II √ πr | A ( θ ) | . (11)Using the corresponding Drucker-Prager failure criterion F with a friction coefficient f , we define the yield region,i.e. the extent of the region where the stress state goes beyond the limits defined by the Drucker-Prager failure criterion. F = J − f I , (12) = K II √ πr (cid:20) | A ( θ ) | − f ν θ (cid:21) . (13)We label the bracketed term as (cid:112) B ( θ, ν, f ) : F ≡ K II √ πr (cid:112) B ( θ, ν, f ) . (14)We can then derive the maximum distance r DP at which the Drucker-Prager failure criterion F = 0 is reached: r DP ( θ, ν, f ) = K II B ( θ, ν, f )2 πF (15)Now consider that the crack is moving at a speed v . The stress intensity factor K II becomes K dynII , defined inequation (4). Also, one can assume that the static intensity tensor is given by K II = ∆ τ √ πL , so that we can updatethe distance r DP for a moving crack: r DP ( θ, ν, f, v ) = (1 − v/c R ) (1 − v/c P ) L ( t ) ∆ τ B ( θ, ν, f )2 F (16)For an accelerating crack, the extent of the yield region r DP is then a combined effect of the increase in L due tocrack growth and the decrease of the term (1 − v/c R ) / (1 − v/c P ) with increasing speed. Numerical dynamic simulations including off-fault coseismic damage
Several efforts has been made to understand the supershear ruptures from a dynamic perspective . However,models accounting for off-fault coseismic change damage has only been recently developed, and should now be used15o explore complex feedback that exist between this particular type of rupture dynamic and the surrounding medium.Here we employ two approaches to simulate dynamic earthquake ruptures: one based on the finite discrete elementmethod (FDEM) and another one based on a constitutive model that accounts for the dynamic evolution of elasticmoduli at high-strain rates . FDEM modeling
To reproduce coseismic damage, we use the software suite HOSSedu, developed by the Los Alamos National Lab-oratory (LANL) . The numerical algorithms behind this tool are based on the combined Finite-Discrete ElementMethod (FDEM) to produce dynamically activated off-fault fracture networks. One of the key FDEM is to alloweach individual interface between the finite elements describing the off-fault medium to have their own tensional andshear cohesion. Furthermore, these interfaces can break under appropiate stress conditions. Each broken interface isthen assimilated to a damage fracture. When the rupture propagates, this allows for a live build-up of the damagepatterns .Closeness to failure is derived using the invariant form of the Mohr-Coulomb yield criterion: F MC ≡ J R MC − I tan φ − C pII (17)where tan φ = f s is initial static friction coefficient and C pII is initial peak shear cohesion. Failure occurs when F MC ≥ . Note that when R MC = 1 , the above become the Drucker-Prager yield criterion.Here, R MC = 1 √ φ sin (cid:16) Θ + π (cid:17) + 13 cos (cid:16) Θ + π (cid:17) tan φ ; Θ = 13 cos − (cid:18) rJ (cid:19) (18)where r = (cid:18) S ij s jk s ki (cid:19) .The closeness to failure yields ∆ CF ( x , t ) = F MC ( x , t ) F MC ( x , t ) − . (19)Here F MC is the Mohr-Coulomb yield function for the initial, uniform, state of stress. Failure is more likely tooccur when ∆ CF > . (See Table S1 for all the parameters used in the simulations).During the supershear transition, the intensity and spatial extent of off-fault damage drops (Fig 1). In fact, whileboth the sub-Rayleigh and the supershear segments present high fluctuations of the stress invariants I and J , thetransitional region is clearly depleted in stress changes (Fig. S1). Micromechanics approach
The second numerical method, reflects the micro-physics of damage evolution by relating damage density to the near-tip stress state and by computing the corresponding dynamic changes of elastic propierties in the medium due to the16resence of newly formed cracks.Here, closeness to failure ∆ CF is simply derived using the Drucker-Prager yield criterion: ∆ CF ( x , t ) = F MC ( x , t ) F MC ( x , t ) − with F MC ≡ J − I tan φ, (20)where tan φ = f s = 0 . is initial static friction coefficient. Failure is more likely to occur when ∆ CF > . Aftershock Catalog Analysis
We have studied three well known supershear ruptures for which high-resolution aftershock catalogs are available, atleast 1-month after the main event. The M w had a mean rupture velocity of ∼ (Fig. 2a for 1-month aftershock locations and Figs. S3-S6 forspatiotemporal evolution of the aftershocks). In the case of the M w ∼ , but in this case, 1-year aftershocks catalog has been published (Fig. 2d for1-month aftershock locations and Figs. S7-S10 for spatiotemporal evolution of the aftershocks). Finally, the M w (Fig. 2g for 1-monthaftershock locations and Figs. S11-S14 for spatiotemporal evolution of the aftershocks). Subpixel correlation methods for damage fault width
The M w is a very well known supershear rupture, but unfortunately theseismic network in the region was poor developed, being studied only by teleseismic data. Because of it, we are notable to apply the high-resolution aftershock analysis developed in the main text. On the other hand, this earthquakewas well recorded by satellite images, that allows us to explore the size of the off-fault coseismic damage zone usingoptical correlation images techniques.We used SPOT 1 to 4 images, covering the Kunlun fault area from 1988 to 2004 with a ground resolution of10 meters, allowing for a change detection larger than a meter . The correlation of pre- and post-earthquakepaired images is done using MicMac , measuring the horizontal surface displacement between the two opticalacquisitions.We focused on the central part of the rupture, dominated by fault parallel motion, excluding the two eastern andwestern ends where fault normal and thrust components are respectively acting. Using the tool Stacking profilesfrom ENVI version 5.5.1 (Exelis Visual Information Solutions, Boulder, Colorado), 460 stacked profiles were placedperpendicular to the fault trace every 500 meters along the analysed section (stacked profiles are rectangles of 40 km-long and 1 kilometer wide in which individual slip profiles are stacked allowing to reduce noise and artefacts). On each17rofile, it is possible to observe the coseismic offset produced by the earthquake (See Supplementary Information, Fig.S15, for a profile example) where the width of this offset corresponds to the fault’s zone width (including fault coreand damage zone), that is computed over 396 clear-enough profiles along the fault trace, allowing to draw a fault zonewidth curve (Fig. 3).Two different means for the fault zone width are calculated (red and dashed red lines in Fig. 3). One is at 238 metersand corresponds to the mean of the fault zone width for the whole area of study, however removing the partitioningzone with two fault strands. The second one, at 127 meters, corresponds to the mean of a fault zone width where sliplocalizes for 11 km-long, located next to what was defined as the supershear transition using seismological data .18 eferences [31] Madariaga, R. High frequency radiation from dynamic earthquake. Ann. Geophys. , 17 (1983).[32] Freund, L. B. The mechanics of dynamic shear crack propagation. J. Geophys. Res. , 2199 (1979).[33] Williams, M. L. On the stress distribution at the base of a stationary crack. J. Appl. Mech. , 109–114 (1975).[34] Dunham, E. M., Favreau, P. & Carlson, J. M. A Supershear Transition Mechanism for Cracks. Science ,1557–1559 (2003).[35] Dunham, E. M. Conditions governing the occurrence of supershear ruptures under slip-weakening friction.
J.Geophys. Res. , 1–24 (2007).[36] Rougier, E., Knight, E. E., Lei, Z. & Munjiza, A. HOSS.edu2.x (Hybrid Optimization Software Suite - Educa-tional Version, Second Generation) (2016).[37] Klinger, Y. et al.
Earthquake damage patterns resolve complex rupture processes.
Geophys. Res. Lett. (2018).[38] Bouchon, M. et al.
Seismic imaging of the 1999 Izmit (Turkey) rupture inferred from the near-fault recordings.
Geophys. Res. Lett. , 3013–3016 (2000).[39] Dunham, E. M. & Archuleta, R. J. Evidence for a supershear transient during the 2002 Denali fault earthquake. Bull. Seism. Soc. Am. , S256—-S268 (2004).[40] Ratchkovski, N. A. et al. Aftershock Sequence of the Mw 7.9 Denali Fault, Alaska, Earthquake of 3 November2002 from Regional Seismic Network Data.
Seismol. Res. Lett. , 743–752 (2003).[41] Holtkamp, S. & Ruppert, N. A high resolution aftershock catalog of the magnitude 7.5 Craig, Alaska, earthquakeon 5 January 2013. Bull. Seism. Soc. Am. , 1143–1152 (2015).[42] Bouchon, M. & Vall´ee, M. Observation of Long Supershear Rupture During the Magnitude 8.1 KunlunshanEarthquake.
Science , 824–826 (2003).[43] Vall´ee, M. & Dunham, E. M. Observation of far-field Mach waves generated by the 2001 Kokoxili supershearearthquake.
Geophys. Res. Lett. , n/a–n/a (2012).[44] Rosu, A.-M., Pierrot-Deseilligny, M., Delorme, A., Binet, R. & Klinger, Y. Measurement of ground displacementfrom optical satellite image correlation using the free open-source software MicMac. ISPRS J. Photogramm.Remote Sens. , 48–59 (2015). 1945] Rupnik, E., Daakir, M. & Pierrot Deseilligny, M. MicMac – a free, open-source solution for photogrammetry.
Open Geospatial Data, Software and Standards , 14 (2017).[46] Tchalenko, J. S. & Ambraseys, N. N. Structural analysis of the Dasht-e Bayaz (Iran) earthquake fractures. Bull.Geol. Soc. Am. (1970).[47] Vallage, A., Klinger, Y., Grandin, R., Bhat, H. & Pierrot-Deseilligny, M. Inelastic surface deformation duringthe 2013 Mw 7.7 Balochistan, Pakistan, earthquake.
Geology , G37290.1 (2015).[48] Mitchell, T. & Faulkner, D. The nature and origin of off-fault damage surrounding strike-slip fault zones with awide range of displacements: A field study from the Atacama fault system, northern Chile. Journal of StructuralGeology31