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







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


Jorge Jara^1,*, Lucile Bruhat^1, Marion Y. Thomas^2, Solène L. Antoine^3, Kurama Okubo^4, Esteban Rougier^5, Ares J. Rosakis^6, Charles G. Sammis^7, Yann Klinger^3, Romain Jolivet^1,8 and Harsha S. Bhat^1 


To appear in Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences






0em

  * Laboratoire de Géologie, Département de Géosciences, École Normale Supérieure, CNRS, UMR 8538, PSL Université, Paris, France

  * Institut des Sciences de la Terre de Paris, Sorbonne Université, CNRS, UMR 7193, Paris, France

  * Université de Paris, Institut de Physique du Globe de Paris, CNRS, 75005 Paris, France

  * National Research Institute for Earth Science and Disaster Resilience, 3-1 Tennnodai, Tsukuba 305-0006, Ibaraki, Japan

  * EES-17 – Earth and Environmental Sciences Division, Los Alamos National Laboratory, Los Alamos, NM, USA

  * Graduate Aerospace Laboratories, California Institute of Technology, Pasadena, California, 91125, USA

  * Department of Earth Sciences, University of Southern California, Los Angeles, CA 90089, USA

  * Institut Universitaire de France, 1 rue Descartes, 75005 Paris, France

  * Corresponding author: 








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 <cit.>. As an earthquake goes to supershear speeds, it manifests Mach fronts that produce unusually large ground motion at distances far from the fault <cit.>. 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 <cit.>. 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 <cit.>. 

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 <cit.>, the absence of field observation and their extreme rarity in laboratory experiments <cit.> 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 <cit.>. Pioneering laboratory experiments <cit.> together with observations from the 1999 M_w 7.4 Izmit and the 1999 M_w 7.2 Düzce earthquakes in Turkey <cit.>, 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 <cit.>, the 2002 M_w 7.8 Denali (Alaska) earthquake <cit.>, the 2010 M_w 6.9 Qinghai (China) earthquake <cit.>, the 2012 M_w 8.6 off-Sumatra (Indonesia) earthquake <cit.>, the 2013 M_w 7.5 Craig (Alaska) earthquake <cit.>, the 2013 M_w 6.7 Okhotsk (Kamtchatka) earthquake <cit.>, and most recently the 2018 M_w 7.5 Palu (Indonesia) earthquake <cit.>.

Several efforts have been made to identify and characterize these supershear earthquakes using various methods <cit.>. 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 kilometers precision). Even in a situation where near-field data exists, such as during the 2002 Denali earthquake <cit.>, 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 <cit.>. Supershear ruptures occur if S < 1.77 in 2D <cit.>, or S < 1.19 in 3D <cit.>, and the fault length is long enough for transition to occur <cit.>. 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 <cit.> and on a non-planar fault <cit.>, while Kaneko & Lapusta <cit.> and Hu et al. <cit.> 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 <cit.>) and numerical models. This problem is mainly because the transition happens over very small length scales (of millimeters 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 modeling efforts helped 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 writing a theoretical analysis 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 v ≤ c_R, where c_R is the Rayleigh wave speed. Let the origin of the polar coordinate system (r,θ) coincide with the crack tip. The near-tip stress, which also depends on the rupture speed v, is given as <cit.>,


    σ_αβ(r,θ,v)    =   K^dyn_II(v)/√(2π r)f^II_αβ(θ,v) + σ^0_αβ    α,β = 1,2
    σ_33   =   ν(σ_11 +σ_22)


where σ^0_αβ is the initial stress state and K^dyn_II is the dynamic stress intensity factor which can be approximated as


    K^dyn_II≈1 - v/c_R/√(1 - v/c_p) K_II^sta(L̂)


where c_p, the P-wave speed, is the limiting speed for a mode II crack and L̂ = ∫_0^t v(t)dt is the current crack length and c_R is the Rayleigh wave speed <cit.>. Here K_II^sta(L) = Δτ√(π L) is the static stress intensity factor of the crack and Δτ is the stress drop defined as σ^0_yx - τ_r. Here σ^0_yx is the initial shear stress on the fault and τ_r is the residual frictional strength. This solution by Freund <cit.> 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 <cit.>. This contraction has already been observed and verified experimentally for ruptures approaching the Rayleigh wave speed <cit.>. 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^dyn_II 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_DP≥ 0, where


    F_DP = √(J_2)+ptanϕ.


Here p≡ I_1/3=σ_kk/3 is the hydrostatic stress derived from the first invariant of the stress tensor, I_1, J_2=s_ijs_ji/2 corresponds to the second invariant of the deviatoric stress tensor (s_ij = σ_ij - pδ_ij) and tanϕ = 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 <cit.>, 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_DP. The closed-form expressions are quite long, but the key variables that affect the solution can be written as


    F_DP = ℱ[1√(r),
        √(L̂)(1 - v/c_R√(1 - v/c_p)),Δτ, θ,ν,tanϕ,σ^0_αβ]


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 <ref>. 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 <ref>a). The instantaneous static stress intensity factor increases with the increasing length of the crack (Figure <ref>b). However, Freund's universal function, F_II(v), monotonically decreases to 0 (Figure <ref>c). Consequently, the dynamic stress intensity factor first increases with increasing crack length but soon starts decreasing as it is dominated by F_II(v) (Figure <ref>d). 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 <ref>d). First, we compute, in Figure <ref>e, 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 <ref>f). 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 behavior 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 <cit.>, the importance, during such event, of the complex feedback between the dynamic rupture and the mechanically evolving medium was not fully considered. 

However, new modeling 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 modeling 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.




 §.§ Modeling strategies




The two numerical models presented in this study account for fracture damage-dominated brittle rheology. Earthquake ruptures induce large dynamic strain rates ( ≥ 1s^-1) around the rupture tip that cannot be fully accounted for by classical Mohr-Coulomb or Drucker-Prager plasticity, typically used low strain rates ( ≤ 10^-5s^-1)  and the two models presented below allow to account for such large dynamic strain.

Okubo et al. <cit.> explicitly model the sub-kilometric off-fault fractures using a continuum-discontinuum numerical framework. The micromechanical model <cit.> accounts for coseismically activated “microfractures” (∼10's of meters) by computing the dynamic changes in the constitutive response due to these fractures. The difference between the two models is essentially in the modeled off-fault fracture networks scale, leading to two distinct modeling strategies. In essence, the two brittle failure models are complementary, and they should, ideally, be combined. Nevertheless, they are good proxies for modeling 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., 50m to several km), Okubo et al. <cit.> use the software suite HOSSedu, developed by the Los Alamos National Laboratory (LANL) <cit.>. The numerical algorithms behind this tool are based on the combined Finite-Discrete Element Method (FDEM) proposed by Munjiza <cit.>, 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 50m, which is defined by the minimum size of the discretized mesh. For a further detailed description of the method, see Okubo et al. <cit.>.

The second modeling strategy relies on laboratory experiments <cit.> and field observations <cit.> 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 meters normal to the fault and few kilometers in depth. Following the damage constitutive laws proposed by Ashby & Sammis <cit.>, Deshpande and Evans <cit.>, Bhat et al. <cit.>, Thomas et al. <cit.> implemented an energy-based micromechanical brittle rheology, as developed by Rice <cit.>, 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 (1) the free energy Ψ^e of a solid, without flaws, deforming purely elastically and (2) 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 <cit.> for an overview). A complete description of the model can be found in Bhat et al. <cit.> and Thomas et al. <cit.>. 




 §.§ 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 2-D plane strain medium, with a 1-D right lateral fault embedded in it and loaded by uniform background stresses. The maximum compressive stress σ_1,  and the minimum compressive stress σ_3  are in the x-y plane, whereas the  intermediate principal stress σ_2 coincides with σ_zz. The fault plane makes an angle of 60^o with σ_1 with normal uniform traction (σ_yy^0) and shear traction (σ_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 σ_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 <cit.>. Static friction is set at 0.6, which corresponds to a value measured in laboratory experiments for a large range of rocks <cit.>, and dynamic friction at 0.1 as observed in high slip-rate experiments <cit.>. 

The fault length is 64km for the micromechanical model and 115km 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 6km with absorbing boundary conditions to avoid reflections interfering with the propagating dynamic rupture. For the FDEM model, we simply set the domain large enough (86km) 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 to natural or laboratory earthquakes. R_0 is defined as


    R_0 = 9π/16(1-ν)μ G/(τ_p-τ_r)^2.


Here ν is Poisson's ratio, μ is shear modulus, G is the fracture energy associated with the slip-weakening law, τ_p and τ_r are the static and dynamic frictional strengths, respectively.

Reference values for the common parameters between the two models are summarized in Table <ref>, whereas parameters specific to different modeling strategies are provided in the Supplementary Information (Tables S1 & S2).  Schematics and parameters for the simulations can also be found in 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 <ref> we display the Drucker-Prager criterion to emphasize the regions where shear sliding is likely to occur, i.e., F_DP>0 (equation <ref>), 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 <ref>a), the damage is essentially generated behind the rupture front, in the tensional quadrant. This feature has been seen by Okubo et al. <cit.> as well, using the FDEM model. Damage also occurs ahead of the rupture front, where the S-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 F_DP is positive.


The transition to supershear velocity (Figure <ref>b) directly impacts the generation of coseismic damage. Around t=6.5s a new pulse is generated ahead of the rupture (Figure <ref>b). The transition to supershear velocity coincides with a decrease in the width of the damage zone (Figures <ref>b and <ref>b-d). As we argued in section <ref>, 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 (L̂ in equation <ref>), has increased significantly. 
Moreover, a snapshot of F_DP at t=7.2s (Figure <ref>b) 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.2s, the second pulse ahead of the rupture front becomes stronger and the first pulse weaker (Figure <ref>d-e). Consequently, damage now essentially occurs behind the rupture front propagating at supershear speed and in the area behind the Mach cone  (Figure <ref>c). 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 <ref>a, the slip continues to accumulate in the “transition zone" after the passage of the rupture, while the damage zone width remains unchanged (Figure <ref>c).




 §.§ Comparison between the different modeling 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 section <ref> and Okubo et al. <cit.>). 
We evaluate the potential regions of failure, in shear, by calculating the invariant form of the Mohr-Coulomb yield function, F_MC (see equation 4.142 of Chen & Han <cit.>):  


    F_MC = R_MC√(J_2) + ptanϕ - c


where c is the cohesion and R_MC is given by:  


    R_MC =sin( Θ + π3)ϕ+ 1√(3)cos( Θ + π3)tanϕ


where


    cos3Θ = 3√(3)2J_3J_2^3/2  ;  
    J_3 = 13s_ijs_jks_ki


where J_3 is the third invariant of the deviatoric stress tensor.

Failure occurs when F_MC≥ 0. Note that when R_MC = 1 and c=0, F_MC becomes the Drucker-Prager criterion as described in equation <ref>.






Figure <ref> 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 (section <ref>). 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 et al. <cit.> first explored the extent and distribution of off-fault plasticity during supershear rupture (see Figure 13 of Templeton & Rice et al. <cit.>). They conducted a 2D in-plane dynamic rupture modeling 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 <cit.>, 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 <cit.>, 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. <cit.> and Klinger et al. <cit.> 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. <cit.> showed that the observed, spatially diffuse, displacement field around a fault is due to displacement accommodated by off-fault fractures (see Figure 3 in <cit.>). This feature is in stark contrast with 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., Vallage et al. <cit.>, Barnhart et al. <cit.>, Delorme et al. <cit.> and Ajorlou et al. <cit.>). 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” <cit.>. The fault zone width  reflects the lateral extent of rock damage, and secondary deformation features around the fault <cit.>. In the field, the fault zone width corresponds to the last deformation features observed when moving away from the fault core <cit.>. 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  <cit.>. 

The 2001 M_w 7.8 Kunlun (China) earthquake is a strike-slip event with ∼ 400km long surface rupture and a mean rupture speed between 3.3-3.9 km/s, which is larger than the shear wave speed and hence reported as a supershear rupture <cit.>. 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 10m, allowing for change detection of less than 1m <cit.>. We focus on the central part of the rupture dominated by fault parallel motion <cit.>, obtaining through optical image correlation the horizontal displacement field for the earthquake east-west and north-south components (Figure S2). The maximum strike-slip displacement for this area of the rupture is around 8 m <cit.>.






Prior to correlation,  SPOT images are corrected from viewing angle and topography using a unique 90m SRTM digital surface model. The pre-and post-earthquake corrected paired images are correlated using the MicMac package <cit.>, which measures the horizontal surface displacement between the acquisitions (see Supplementary Information, Figure S2 to see the data employed in this study). The MicMac method enables preserving the resolution of the input images, that is 10m, 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 50m. A regularization of 0.3 is used to force MicMac to consider realistic displacement values relative to neighboring 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 (Figure S2) every 500 m to analyze variations in the surface displacements and fault zone width. Individual slip profiles are stacked on rectangles of 40km-long and 1km-wide, in order to average out noise and artifacts (Stacking profiles from ENVI version 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 Supplementary Information, Figure S3 for a profile example). This procedure allows us to evaluate the fault zone width in the strike-slip direction, over 396 profiles (Fig. <ref>b). Field measurements of fault zone widths are at least a few tens to a few hundred meters in the case of this earthquake <cit.>. 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 ± 80m, and it corresponds to the strike-slip section of the rupture where only one fault zone is observed at the surface (Figures <ref>, S4A , S5). The second domain presents a width larger than 1km lasting for about 15km. 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 <cit.> (Figures <ref>, S4B, S5). The third domain is located right before the onset of the slip-partitioning region, where the average fault zone width drops to 127 ± 39m over  11km long segment. This value is about half of the mean fault zone width of 238m 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 <cit.>. 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 <cit.>. 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 modeling 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 <ref>). 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-analyzing 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 analyze three reported supershear ruptures for which high-resolution aftershock catalogs 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.8km/s <cit.>, and a well-recorded 3-month aftershock catalog <cit.> (Figure <ref>a for 1-month aftershock locations). In the case of the 2002 M_w 7.9 Denali earthquake (Alaska), rupture velocities of about 5.5km/s were reported <cit.>, and 1-year long catalog of aftershocks is available <cit.> (Figure <ref>b for 1-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 to 6.0km/s <cit.> and a 5-month aftershock catalog is available <cit.> (Figure <ref>c for 1-month aftershock locations). 

We observe a paucity in aftershocks where the ruptures have transitioned to supershear, as highlighted by the pink boxes in Figures <ref>b, d, and f. We quantify such lack of aftershocks by computing the cumulative moment density released by aftershocks, considering those located within l_a=5km of the fault trace (see the Supplementary Information for the same analysis considering l_a=2.5km, Figure S6). We assume that each aftershock, having a seismic moment M_0, is a circular crack and compute the slip distribution δ(r) using <cit.>: 


    δ(r)    =   24Δτ7πμ√((7M_016Δτ)^2/3-r^2)


where r is the crack radius, μ = 30GPa is the shear modulus, and Δτ is the stress drop, assumed to be 3MPa. 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 90m for Izmit, 32m for Denali, and 40m for Craig earthquake. With the slip distribution, the cumulative seismic moment density of the i^th bin containing N aftershocks is given by


    m^(i)_c  ≈μ∑_j=1^Nδ^(i)_j


Note that δ^(i)_j is part of the slip distribution of the j^th aftershock projected on the i^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 density 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: 1 and 3 weeks (Figures <ref>a, c, and 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 (Figures <ref>b, d, and f, pink boxes). The extent of these regions is different for each earthquake c. f., 3.6km for Izmit, 6km for Denali, and 3.2km 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 5km from the main fault in 1 week (3-week), 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 catalog. We employ a multivariate normal distribution {𝐱} = 𝒩_2(μ_𝐚, σ_𝐚), assuming that the location uncertainty of aftershocks is normally distributed. Here, μ_𝐚 corresponds to the aftershock epicenter, while σ_𝐚 is the covariance matrix of the location. This procedure allows us to generate, randomly, 10000 synthetic catalogs and repeat the along-strike moment density evaluation with time for each of the 10000 catalogs, deriving the mean and the standard deviation (1-σ and 10-σ) of the spatio-temporal evolution of seismic moment density (red areas in Figure <ref> for 1-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 Supplementary Information for a comparison between 1 and 3-week using 1-σ and 10-σ, Figures S7, S8, and S9). 
 





We consider 3-week as the maximum reasonable period to investigate the evolution of the seismic moment for these earthquakes to avoid potential effects of postseismic deformation. Afterslip is proposed as a mechanism driving the aftershock triggering (e.g., Hsu et al. <cit.> and Ross et al. <cit.>). For each studied event, the inferred transitional region co-locates with areas where some authors report afterslip occurrence <cit.>. Therefore, the observed gap in the early aftershocks productivity (less than 3-week 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. <cit.>, employing Harvard CMT solutions, also allude to the same conclusion. However, due to the lack of high spatio-temporal density of aftershocks in their catalog, 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 <ref>), 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 (≤ 3-week after the mainshock), and the cumulative seismic moment density around the fault (≤ 5km 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 catalogs. 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 <cit.>, Denali <cit.>, Craig <cit.>, and Kunlun <cit.> 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 kilometers. 

While the theoretical and numerical models are still idealized (modeled on a 1D 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 <ref> and <ref>). 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 kilometers, 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 <cit.>, 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 <cit.>, impeding large ground accelerations near the fault <cit.>. Our work, in 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 <cit.>. 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 ACCESS
All the catalogs and numerical simulation results employed in this work have been obtained from published works, cited in the main text and Supplementary Information.



§ AUTHOR CONTRIBUTIONS
'HSB conceived of and designed the study. HSB developed the LEFM solution. MYT conducted damage mechanics based numerical modeling and KO performed FDEM based numerical modeling. SA conducted the image correlation analysis. JJ analysed the aftershock catalog with inputs from MYT and HSB. JJ and LB drafted the first version of this manuscript that was subsequently modified with inputs from HSB, MYT, KO, SA, YK, AJR and RJ. All authors read and approved the manuscript'.



§ FUNDING
HSB acknowledges the European Research Council grant PERSISMO (grant 865411) for partial support of this work. JJ and RJ acknowledge the funding of the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement 758210, project Geo4D). LB 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 PCOFUND-GA-GA-2013-609102, through the PRESTIGE program coordinated by Campus France. SA and YK 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 program from CNES.



SUPPLEMENTARY INFORMATION













§ NUMERICAL MODELING OF OFF-FAULT COSEISMIC DAMAGE



Specific parameters used in the FDEM model are provided in table <ref> and those used in the micromechanics based model are provided in table <ref>.





















































naturemag

