Supershear Rupture, Daughter Cracks, and the Definition of Rupture Velocity

Rupture velocity is a fundamental feature of earthquake behavior. How rupture velocity is defined can therefore affect our understanding of earthquake physics. Based on dynamic rupture simulations, we compare the tangent and average rupture velocities calculated via one‐ and two‐dimensional gradients. Two strike‐slip scenarios with free‐surface‐induced supershear ruptures are presented within a homogeneous and a depth‐dependent stress regime, respectively. Although both scenarios produce a daughter crack that propagate over the seismogenic depth, in the depth‐dependent case the 1D and 2D rupture velocities capture different features: A 1D horizontal gradient measurement implies a supershear rupture, while a 2D gradient measurement reveals sub‐Rayleigh rupture propagation everywhere except very close to the free surface. A large area on the fault with the 1D horizontal tangent rupture velocity does not necessarily produce an observable Mach front, which arises the importance of rupture velocity definitions in supershear analysis.

A perhaps subtle issue embedded in the above discussion is the formal definition of rupture velocity. In his 3D dynamic rupture modeling of earthquakes with heterogeneous pre-stress (Day, 1982), notes two ways to calculate rupture velocity. The first is the local or "tangent" rupture velocity, which can be calculated at every point on the fault as the reciprocal of the gradient of the rupture time at that point. The second is the average or "secant" rupture velocity, which is simply the distance of a point from the hypocenter divided by the rupture time at that point (Day, 1982). shows that these two calculations can give quite different results for the same earthquake, including supershear speed in the former but not in the latter. Using 3D dynamic rupture models of faults with both homogeneous and depth-dependent initial stress, we explore these and other methods of calculating rupture velocity in the presence of free-surface-induced supershear rupture. We find that in some cases, different methods of calculating rupture velocity can give very different results as to whether the rupture is propagating at supershear speed. These different methods of calculation may be linked to different physical observables in real earthquakes, with implications for the detection and interpretation of supershear rupture.

Methods
We use the 3D dynamic finite difference method (Z. Zhang et al., 2014), which has been validated following the SCEC/USGS Dynamic Rupture Code Verification exercises (Harris et al., 2009(Harris et al., , 2018, to model earthquakes on vertical planar faults of size 60 km along strike and 15 km in depth, which intersect the surface of a homogeneous half-space. We assume linearly elastic material behavior and slip-weakening friction (Andrews, 1976;Palmer & Rice, 1973). For a given point on the fault, rupture starts when the shear stress  0 increases to the yielding strength     u s n , and linearly decreases to the residual stress     d d n as the slip increases to the characteristic slip-weakening distance 0 d , which equals 0.4 m in this study. The static and dynamic friction coefficients  s and  d are 0.677 and 0.525, respectively. The P wave speed is 6 km/s, and the S wave speed is 3.464 km/s, with density 2,670 kg/m 3 . To allow comparison of our results with prior work, we focus on and duplicate two of the models of Kaneko and Lapusta (2010) by applying a time-varying shear stress perturbation (Equation A1 of Supporting Information) while performing a new analysis. Our prior work has shown that the results are sensitive to the physical depth of nucleation, but not to the details of the artificial nucleation procedure (Hu et al., 2019). The first model is a homogeneous stress model whose initial normal stress is 120 MPa, while the normal stress for the depth-dependent stress model is min[1.0 + 16.2z, 120] MPa, with z equals to the depth in kilometers ( Figure S1). The relative fault is 1.6 for the two models, which prevents supershear transition via Burridge-Andrews mechanism (Andrews, 1976;Das & Aki, 1977;Day, 1982;Dunham, 2007). A new analysis on the rupture speed of the three homogeneous stress models (presented by Hu et al., 2020) are also shown. We focus on the free-surface-induced supershear mechanism, although the results have implications for supershear rupture in general, especially for the cases with a large angle of the daughter crack propagation with respect to the horizontal direction. We use a spatial grid of 50 m and a time step of 0.0025 s; grid size independence is verified by good comparison of the results with selected models with half the spatial grid and time step.
To calculate rupture speed, we utilize four different definitions.
The 2D tangent rupture speed is calculated following Bizzarri and Das (2012): where   , r t y z is the rupture time for a point on the fault, and and h is the grid spacing. The 1D horizontal tangent rupture speed is calculated using: The 2D secant (average) rupture speed is calculated using: Where Dist_hypo is the distance from the point to the hypocenter.
Finally, the 1D horizontal average rupture speed is calculated using: Where Dist is the horizontal distance between the observer ( , y z) and the point ( 0 , y z ). In their work on free-surface-induced supershear rupture, Kaneko and Lapusta (2010) and 1d r V to characterize their ruptures as supershear or subshear. Weng and Ampuero (2020) used an average of 2d r V and 1d r V to characterize rupture speed in their study of supershear rupture on oblique faults. Here we compare the results for our four definitions and explore the implications of their differences, which may highlight the challenge of interpreting inferences of supershear rupture.

Results
We compare slip rate snapshots, rupture time distributions, and different methods of calculating rupture speed for the homogeneous and depth-dependent stress models in Figure 1. The slip rate snapshots of the two models (Figures 1a and 1f) match the results of Figures 6 and 3 of Kaneko and Lapusta (2010), respectively. Both models produce a daughter crack that propagates out ahead of the main rupture front, as can be seen by the thin zone of high slip rate to the right of the main crack in the rupture snapshots, as well as by sharp kinks in the rupture time contours. Both daughter cracks eventually propagate down to the bottom of their respective faults. However, there is a qualitative difference between the two daughter cracks. In the homogeneous stress case, the daughter crack develops a largely vertical shape while propagating along strike. In contrast, the daughter crack in the depth-dependent stress case retains a diagonal shape with the crack leading at the free surface as it propagates to the right. Figures 1b-1e, the homogeneous stress model shows a free-surface-induced supershear transition by all measures of rupture speed, and as the daughter crack propagates, the area experiencing supershear rupture speed gradually grows to encompass the entire depth of the fault. Nonetheless, important differences exist between the different rupture speed estimates. The 2D gradient 2d r V shows that the area undergoing supershear rupture does not necessarily correspond precisely to the area of the fault for which the daughter crack arrives earlier than the parent crack (the daughter crack arrives first for points on the fault above the sharp kink in rupture contours). In other words, even if the daughter crack causes a point on the fault to rupture ahead of the main sub-Rayleigh crack, that does not necessarily imply that the propagation speed is supershear at that point. Conversely, the supershear area defined by the 1D horizontal gradient 1d r V precisely corresponds to the boundaries of the area in which the daughter crack arrives before the parent crack. Finally, both average rupture speeds producing a higher estimate of rupture speed over more of the fault.

As shown in
Differences between the various rupture speed calculation methods for the depth-dependent stress model are more significant. For this model the 2D gradient 2d r V indicates a rupture that is almost entirely subshear (Figures 1g-1j). Close examination reveals that there is a very thin (<200 m in depth) zone of horizontally propagating supershear rupture immediately below the surface (Figure 2b), but the zone is too small to be seen in the plots in Figure 1g. The daughter crack is revealed to be a sub-Rayleigh crack propagating diagonally downwards from the free surface. The orientation of the rupture contours for the daughter crack is similar to the wavefronts of a converted horizontal P to downward S wave along the fault ( Figure S2). In contrast to the almost entirely sub-Rayleigh 2d r V , the 1D horizontal gradient 1d r V indicates a large area propagating at supershear speed-the entire area for which the daughter crack arrives before the parent crack. Inspection of the rupture time contours indicates the reason for this supershear estimate: because the daughter crack is propagating at a large angle with respect to the horizontal, the rupture time contours are farther apart in the along-strike direction than in the direction normal to the contours (i.e., the direction of the 2D gradient), leading to a larger apparent velocity in that direction ( Figures S3-S4). The effect is analogous to the large apparent horizontal phase velocity seen on the surface of the Earth from a diagonally incident seismic wave. The 2D secant (average) rupture speed secant r V and the 1D horizontal average ave1d r V both show a significant amount of supershear rupture on the fault, but less than the 1D horizontal gradient 1d r V . It should be noted that in this model, the secant rupture speed displays a larger area of supershear rupture than the 2D gradient speed, which is the opposite of the effect observed by (Day, 1982), in which the 2D gradient speed displayed supershear rupture speed where the secant speed did not.
More information about the spatial patterns of the various rupture speed estimates can be gained by examining their distributions along strike, as shown in Figure 2. In the homogeneous stress case, at the free surface all rupture speeds with the exception of secant r V are supershear and above the Eshelby speed, which is the lower limit for the stable steady supershear regime ( (Burridge et al., 1979). Rupture speeds at a depth of 0.2 km are similar to those at the free surface. At the hypocentral depth of 7.5 km, 1d r V has an almost step-function increase from subshear to supershear as it passes to the part of the fault ruptured initially by the daughter crack, while the other measures of rupture velocity more gradually increase to supershear speed along strike. In particular, 2d r V eventually increases to above the Eshelby speed. In contrast to the homogeneous stress case, the depth-dependent stress model displays clear supershear rupture at the free surface for all rupture velocity definitions except for 2d r V , which appears to vary between supershear and subshear along strike. ave1d r V is above the Eshelby speed, while other definitions are below. An experiment with a smaller spatial grid of 25 m ( Figure S5) indicates that there is a very thin zone of clear supershear rupture right at the free surface by the 2d r V definition, but this speed is still below the stable Eshelby speed. At 0.2 km depth, 2d r V is clearly subshear, while the other definitions indicate supershear rupture. The thin, surficial nature of the supershear rupture is why there is no observable blue color for this model in Figure 1g. At hypocentral depth we can see greater differentiation between the different rupture velocity calculations: only 1d r V becomes clearly supershear (but sub-Eshelby) for areas ruptured first by the daughter crack.
HU ET AL.  . Dashed black lines depict the shear wave velocity and the Eshelby speed ( 2 S V ). The upper limit of the vertical axis indicates the P wave velocity.
An issue of considerable practical importance in supershear rupture observations is the existence of the shear Mach front, which is generated by superposition of S-wave fronts (Bernard & Baumont, 2005;Dunham & Archuleta, 2005;Dunham & Bhat, 2008;Mello et al., 2010) and is observed in laboratory crack experiments (Rosakis, 1999;Rubino et al., 2019;Xia et al., 2004). Numerical simulations show that besides the shear Mach front, a Rayleigh Mach front also propagates to the far field with little attenuation (Dunham & Bhat, 2008). Figure 3 displays snapshots of fault slip rate and strike-perpendicular ground velocity for the homogeneous (a and b) and depth-dependent (c and d) stress cases. The homogeneous stress case produces clear Mach fronts on the Earth's surface for both the Shear and Rayleigh waves; these fronts also are observed on both the fault itself and on a vertical cross-section perpendicular to strike. The depth-dependent stress case produces almost no perceptible Mach front in the fault-perpendicular component, although a very weak shear Mach front is perceptible in the fault-parallel component ( Figure S6). Moreover, the Mach cone half-angle to the fault plane which is an unstable supershear rupture speed (Mello et al., 2010). As shown in Figures 3b and 3d, the shear Mach cone half-angle is  40.2 , while the Rayleigh Mach cone half-angle is 35.5°; both  values are lower than  / 4, representing a stable supershear speed. However,  is 48.6° for the depth-dependent stress case (Figure 3d), which emphasizes the unstable supershear velocity.

Discussion and Conclusions
Our numerical results imply that the definition of rupture speed based on the two-dimensional rupture time gradient ( 2d r V ) is most consistent with traditional ground-motion-based inferences of supershear rupture, such as the formation of a Mach front and stronger FP than FN initial pulses (Aagaard & Heaton, 2004;Dunham & Bhat, 2008). However, the inferences of supershear rupture speeds can be interpreted to have different consequences, depending on how the rupture speed is calculated. For example, back-projection studies of long along-strike ruptures tend to trace the movement of the zone of maximum energy radiation in map view (Bao et al., 2019;Ishii et al., 2005;Walker & Shearer, 2009). Associating this zone of maximum energy radiation with the location of the rupture front, it is straightforward to take the horizontal distance HU ET AL.
10.1029/2021GL092832 6 of 10 than with 2d r V . For the homogeneous stress model, there is little distinction between the amount of the fault that experiences supershear rupture between the different rupture velocity definitions. This model is unequivocally supershear by all definitions, and would also likely be indicated to be supershear by all observables. However, the depth-dependent model is more equivocal. It indeed does have at least a small zone of supershear rupture by all definitions, but the area of the fault that is indicated to be supershear differs and depends on the definition. As shown in Figure S7, this more equivocal case does not display strongly the near-source ground motion signatures (Mach front, strong FP pulse) expected from supershear rupture. One can imagine that the famous Pump Station 10 record in the 2002 Denali Fault earthquake (Dunham & Archuleta, 2004) might not have produced evidence of supershear rupture if it had had a stress distribution similar to our hypothetical depth-dependent case.
The depth-dependent case has an extremely small percentage of the fault experiencing supershear rupture by the 2d r V definition, which explains why it does not generate a significant Mach front or strong FP pulse: Due to its small area and low stress drop, the supershear 2d r V zone radiates little energy and does not strongly affect the ground motion. Of course, it is possible that in the real world, stable friction near the free surface, or velocity-strengthening behavior at shallow depth (Kaneko et al., 2008) could either eliminate this surficial supershear rupture. Our depth-dependent stress model produces a striking diagonally downward-directed rupture from the daughter crack at the surface. Tang et al. (2020) also noticed this feature on dipping strike-slip faults with largely subshear rupture, and argued that the rupture which was not directed purely in the mode II direction suppressed supershear rupture (analogous to the mixed-models of Andrews, 1994 andAmpuero, 2020). Our interpretation, in contrast, is that the diagonal rupture fronts are a consequence of the supershear rupture right at the surface leading the subshear rupture at depth; the diagonal rupture fronts in our interpretation are a symptom, rather than a cause, of the subshear rupture at depth. Hu et al. (2019Hu et al. ( , 2020 have shown that the existence of a daughter crack in strike-slip faults cannot guarantee a sustained supershear rupture, since the daughter crack is possible to be absorbed by the following main crack, which turns into an unsustained supershear rupture. Our simulations show even if the "supershear" (in the 1D gradient view) daughter crack can propagate to the fault bottom, its near field ground motion could still present sub-Rayleigh characteristics as revealed in the 2D gradient. Thus, it is interesting to compare the two metrics of rupture velocities in characterizing unsustained supershear ruptures. As shown in Figure 4, we compare the 1D and 2D rupture velocity distributions of the three rupture scenarios HU ET AL.   Figure 1 of Hu et al. (2020). The black lines mean the rupture front contours with a 1.0 s time interval. The ruptures are nucleated with an high stress over an elliptical region (semimajor axis of 3 km and semiminor axis of 1 km) at different hypocenter depths of (a, d) 5 km, (b, e) 7.5 km, and (c, f) 10 km. Initial stresses and friction parameters are listed in Table S1. in a homogeneous stress model (Figure 1 of Hu et al., 2020). The hypocenter depth is the only different simulation parameter among the three cases which transits from a sub-Rayleigh rupture (Figure 4d) to a sustained supershear rupture (Figure 4f) as the hypocenter depth increases from 5 to 10 km. The difference between the 1D and 2D rupture velocities in the sustained supershear case in the homogeneous stress model (Figures 4c and 4f) is similar as discussed in Figure 1. However, for the unsustained supershear cases, both 1D and 2D metrics can unveil the fading away process of the supershear daughter crack, although slight differences still exist, especially when the angle of the tangent direction of the rupture front deviates from the vertical direction.
Thus, the presence of the leading daughter crack makes it tempting to identify the fault areas for which it arrives first as undergoing supershear rupture, but it is possible that this supershear speed is only obtained via one-dimensional gradient or average velocity definitions, not two-dimensional local gradients like 2d r V , especially when the daughter crack is downward propagating. More information than just the presence of a daughter crack is necessary to identify widespread supershear rupture.
In conclusion, we argue that the most physically meaningful definition of supershear rupture is based on 2d r V : earthquake ruptures with large areas that are supershear by this definition will be considered supershear using a range of observables. However, other definitions of supershear rupture may well be useful for other seismological techniques, such as back projection or surface wave polarization rotation. In reality, fault zones (Huang et al., 2016), fault roughness (Bruhat et al., 2016), off-fault yielding (Bhat et al., 2007), velocity-strengthening friction near the Earth's surface (Kaneko et al., 2008), sediment effect (Xu et al., 2021), and topography effect (Z. Zhang et al., 2016) may make rupture fronts complicated, which further complicates rupture velocity behavior and thus the differences between 1d r V and 2d r V . However, no matter which definition of rupture velocity is chosen, we argue that the existence of a Mach front produced by an earthquake is evidence of a true supershear rupture.