Research ArticleGEOPHYSICS

The equation of motion for supershear frictional rupture fronts

See allHide authors and affiliations

Science Advances  18 Jul 2018:
Vol. 4, no. 7, eaat5622
DOI: 10.1126/sciadv.aat5622


The rupture fronts that mediate the onset of frictional sliding may propagate at speeds below the Rayleigh wave speed or may surpass the shear wave speed and approach the longitudinal wave speed. While the conditions for the transition from sub-Rayleigh to supershear propagation have been studied extensively, little is known about what dictates supershear rupture speeds and how the interplay between the stresses that drive propagation and interface properties that resist motion affects them. By combining laboratory experiments and numerical simulations that reflect natural earthquakes, we find that supershear rupture propagation speeds can be predicted and described by a fracture mechanics–based equation of motion. This equation of motion quantitatively predicts rupture speeds, with the velocity selection dictated by the interface properties and stress. Our results reveal a critical rupture length, analogous to Griffith’s length for sub-Rayleigh cracks, below which supershear propagation is impossible. Above this critical length, supershear ruptures can exist, once excited, even for extremely low preexisting stress levels. These results significantly improve our fundamental understanding of what governs the speed of supershear earthquakes, with direct and important implications for interpreting their unique supershear seismic radiation patterns.


The onset of frictional sliding occurs through the dynamic propagation of rupture fronts (1). Local slip occurs only once these fronts rupture the microcontacts that define the rough contacting surfaces that separate sliding bodies. Propagation speeds of idealized rupture fronts that are driven by singular stresses at their tips are limited by the Rayleigh wave velocity (2, 3). However, when regularization of these singularities is taken into account, theoretical (4) and numerical (5) work has predicted the existence of supershear cracks, a class of cracks that propagate beyond the shear wave speed and may approach the longitudinal wave speed. While supershear ruptures have since been observed in laboratory friction experiments (610) and inferred in natural earthquakes (1114), it is not completely understood what determines their speed. Understanding what determines the speed of supershear earthquakes has important implications, as their propagation speeds strongly influence the structure of their radiated waves. The supershear radiation structure, which is fundamentally different from sub-Rayleigh earthquake radiation (15), has important consequences for resulting seismic hazards (16).

While the conditions for the transition of ruptures from sub-Rayleigh to supershear regime have been studied extensively (5, 1722), the question of what determines the evolution of supershear rupture speed is still entirely open. In the sub-Rayleigh regime, using the framework of linear elastic fracture mechanics (LEFM), the equation of motion for accelerating ruptures has been constructed via the superposition of solutions for constant-velocity (singular) shear cracks (23). This solution has recently been validated experimentally for frictional ruptures (2, 3). In the supershear regime, however, this approach fails (24, 25), as the elastic fields in the crack tip vicinity are coupled to the crack velocity history.

Here, we construct an approximate equation of motion for supershear cracks by using a particular solution given by Broberg (26). This equation of motion builds on previous work on energy flux for supershear cracks (2628) but provides a direct link between crack length and rupture speed. It therefore enables us to predict supershear crack speeds for various nonuniform systems. We will show that this equation of motion provides a good description of both experiments and numerical observations of frictional supershear ruptures that have been proposed to model natural earthquakes. Our results demonstrate the broad applicability and the scope of these approximate solutions. These results demonstrate how predictions based on fracture mechanics can be successfully extended to the supershear regime.


Sustained supershear propagation

Our experimental setup consists of two poly(methyl methacrylate) (PMMA) plates of same thickness that are pressed together by an applied normal load (Fig. 1A). We use both dry and boundary-lubricated frictional interfaces (3). Shear forces FS are applied, and slip events are either nucleated spontaneously (7) or induced (3) by a slight out-of-plane shear perturbation at x ≈ 0. A high-speed camera at 580,000 frames/s records the dynamic changes of the real contact area A(x, t). Figure 1B shows an example of a frictional rupture event that nucleates at x ≈ 0, accelerates, and transitions at x ≈ 50 mm to supershear speeds (Cf > CS), which it maintains until reaching the leading edge.

Fig. 1 Experimentally measured supershear rupture velocities.

(A) Two PMMA blocks are used in a stick-slip friction experiment. We consider the elastic medium to be 2D with a quasi 1D frictional interface and measured the full 2D tensorial strains along and ≈3.5 mm above the frictional interface (blue squares). (B) Real contact area A(x, t) measurements along the interface, normalized before the event, show rupture nucleation at x ≈ 0, acceleration, and transition to supershear at x ≈ 50 mm. (C) Shear stress variations, τ = σxyEmbedded Image, relative to the rupture tip arrival time, ttip, measured at x = 105 mm. Arrows denote the dynamic stress drop τ0 and peak shear strength τp. We obtain τp from the measured cohesive zone size. The gray curve schematically depicts the interfacial shear stress evolution. (D) Measured rupture velocities for two examples of rupture events in experiments with a dry interface. (B), (C), and the blue example in (D) correspond to the same rupture event.

In addition, the stresses σij are calculated from strains measured every 1 μs at multiple locations along (and slightly above) the frictional interface (2). Figure 1C shows the variation of the shear stress Embedded Image (Embedded Image is the residual frictional resistance) and the dynamic stress drop τ0, as the rupture passes by the strain measurement location. At the interface, we expect the shear stress to increase to the peak shear strength of the frictional interface (see definition of τp in Fig. 1C). Although this highly localized increase of stress is not directly measured by our (slightly) off-fault measurements, τp as well as the fracture energy Г (the dissipated energy per unit area), were obtained independently (2, 3) by analyzing the sub-Rayleigh crack propagation regime for this system. Briefly, Г was extracted from the 1/r LEFM singular form of the stresses, and τp was obtained from measurements of the cohesive zone size (details are provided in Materials and Methods). While inhomogeneities along the frictional interface are reflected in spatial variations of both and Г and τp, we find that Г and τp are independent of the rupture velocity. We therefore use these inferred values in analyzing the behavior of the system in the supershear regime.

The rupture front velocity Cf(l), inferred from the measurements of A(x, t) for the event presented in Fig. 1B, is shown in Fig. 1D (blue curve), along with an additional rupture event (red curve) whose high propagation speed approaches the longitudinal wave speed CL. Here, l denotes the rupture length (as ruptures are nucleated at x ≈ 0), and Embedded Image. What dictates Cf(l) and how it relates to τ0 and the frictional properties of the interface is the central question investigated in this article.

We now briefly describe how LEFM can be used to provide quantitative predictions for the speeds of supershear frictional rupture fronts. Within the wake of frictional ruptures, the bodies are always in partial contact [σyy(x) ≠ 0], and the frictional resistance, Embedded Image, opposes sliding. The problem of frictional rupture fronts, propagating within an interface separating identical materials, can be mapped to the stress-free conditions that define the mode II crack problem (29) by using the linearity of the governing equations and defining the shear stress variations, Embedded Image (see definition of τ, τ0, and τp in Fig. 1C). The stress fields near the tip of an idealized supershear crack with no cohesive zone are singular (26, 30, 31), τ ~ K/(xl)g, where K is the stress intensity factor and where, in contrast to sub-Rayleigh ruptures, the singular exponent g depends on Cf [with g(Cf/CL) ≤ 1/2]. K has been calculated explicitly only for a handful of particular examples, such as a semi-infinite crack subject to a pair of concentrated shear forces (27) or for a symmetric bilateral crack subject to uniform remote shear stress (26). In both cases, cracks were assumed to expand at constant velocity. A consequence of the singular description of supershear cracks is the vanishing energy flux into the crack tip for Cf ≠2CS. It was shown (31), however, that when a cohesive zone is introduced, a finite region where these singularities are regularized to a finite τp value, the requirement for a positive energy flux to the crack tip is fulfilled for any Cf > CS. By combining the cohesive zone model with the constant-velocity singular crack solution, the energy flow to the crack tip per unit area (the energy release rate) G can be calculated (2628). Here, we follow Broberg’s bilateral crack solution (26), where G was expressed in terms of Cf, τ0, l, and the cohesive zone size xc. Here, motivated by recent experimental results (2), which showed that the interface shear strength τp is roughly independent of Cf (in contrast to the Cf dependence of xc), we rederive G in terms of and Cf/CL, l, τ0, and τp. Balancing G and Г yields a crack propagation criterionEmbedded Image(1)where μ is the shear modulus and Embedded Image(Cf/CL) and Embedded Image are known functions. The shape of the stress distribution within the cohesive zone is carried within Embedded Image [see Materials and Methods for the derivation and definition of Embedded Image(Cf/CL) and Embedded Image].

Studies of supershear crack arrest (24) showed that steady-state singular fields are not established instantaneously. Therefore, a solution of an accelerating crack cannot be constructed by the superposition of constant-velocity crack solutions, as had been performed previously for sub-Rayleigh propagation (23). In this sense, while Eq. 1 was strictly derived to describe cracks propagating at a constant velocity, we hypothesize that its applicability can be extended for “slowly” accelerating cracks, where velocity history can be neglected. As Eq. 1 implicitly links Cf to l, it will provide an equation of motion for supershear cracks.

The salient property of the supershear crack growth criterion (Eq. 1) is its dependence on the cohesive zone characteristics, in striking contrast to sub-Rayleigh brittle cracks (9, 30). In addition, note the emergence of a characteristic length scale Embedded Image, which is proportional to the static cohesive zone size. Finally, we note that the construction of Eq. 1 describes velocity response to spatial changes in interface properties [Г(x) and τp(x)], as for sub-Rayleigh cracks. The rupture response to nonuniform spatial profiles, τ0(x), is not captured by Eq. 1. We will consider this question in a later section of this article.

To investigate the applicability of the LEFM equation of motion, we supplement the experiments with two-dimensional (2D) numerical simulations of dynamic bilateral shear cracks using a spectral boundary integral method (32). The modeled system consists of two linear elastic semi-infinite half spaces with an interface governed by a linear slip weakening friction law (see Materials and Methods) and τ0 applied along the interface. Poisson ratio v = 0.35 and plane strain are assumed (CS/CL = 0.48). We use a weak patch (τ0p = 0.8) around the nucleation center to promote a direct transition (20) to supershear propagation. To test the predictions of Eq. 1 for various τ0p values, we fix the value of τ0p to 0.2 to 0.7 after a finite propagation distance (see Fig. 2A, top) by increasing τp and Г (l0 remains constant). This leads to an instantaneous response of the crack.

Fig. 2 Comparing theoretical predictions of supershear crack velocities with numerical simulations and experimental measurements.

(A) Top: Spatially uniform τ0 and nonuniform τp profiles are considered in simulations. The imposed τ0p profiles are shown. Bottom: Colors represent the crack velocities Cf(l) corresponding to the stress profiles in the top panel. τp is low, dashed lines, near the point of nucleation (l/l0 = 0) and increases at distances 50 for τ0p > 0.2 and at l/l0 = 130 for τ0p = 0.2 (red). Black solid lines denote predictions of Eq. 1. For sufficiently low τ0p (red example), two solution branches to Eq. 1 exist (dashed line shows the unstable solution). (B) Top: Profiles of the measured shear stress τ0 and estimated interfacial strength τp in experiments. Bottom: Measured rupture velocities for rupture events with stresses shown in the top panel. Theoretical predictions according to Eq. 1 are shown (solid lines); average τ0 values are used, l0 = μГ/Embedded Image = 2 mm. Dashed lines correspond to the error estimates of l0 and τp. (C) Cf, measured at l = 120 mm (averaged over ±10 mm), for multiple experiments is plotted with respect to the measured τ0 profiles averaged over the same interval. Solid dots and open circles indicate experiments performed with dry and boundary-lubricated interfaces, respectively. Black lines are the theoretical predictions with the estimated errors.

Figure 2 demonstrates that Eq. 1 provides an excellent description of both the numerical simulations and the experiments with no adjustable parameters. Lower τ0 are associated with slower ruptures, as observed in early numerical work (33) and in sub-Rayleigh propagation (3). In our simulations, we observe transient discrepancies for the τ0p = 0.2 case due to the strong localized increase of the interfacial strength (l/l0 = 140 to 170). We expect these discrepancies when strong crack accelerations are present, as was previously discussed. We note that other types of transition to the supershear regime, for example, supershear seed crack or transition through a secondary crack, have been applied and result in equivalent observations (see Materials and Methods). In our experiments, rupture nucleation, acceleration, and the supershear transition are sensitive to the high stress gradients in the block edge vicinity and are beyond the scope of the current study. Instead, we concentrate on the sustained supershear propagation within a central region where these edge effects are negligible. Note the spatial variation of the τp obtained profile and the accompanying response of the equation of motion (Fig. 2B).

Figure 2C demonstrates the generality of the experimental results, where we consider both dry and boundary-lubricated interfaces, each with significantly different values of Г and τp (3, 34). Here, we compare the predictions of Eq. 1 at a fixed spatial position for multiple rupture events driven at different levels of τ0p. Slight discrepancies for the boundary lubrication case might result from nonuniformity of τ0 along the interface, which is not taken into account here.

In systems with uniform τ0, τp, and Г, the sub-Rayleigh–to–supershear transition [the Burridge-Andrews mechanism (46, 9, 18)] occurs only along weak interfaces above a critical value of τ0p, and no supershear crack can exist for lower values of τ0p. These critical values depend on CS/CL; τ0p = 0.31 and 0.37 for CS/CL = 0.48 (simulations) and 0.57 (experiments), respectively. Both our simulations and experiments reveal strong hysteretic behavior; once a supershear crack has been excited [for example, the transition is facilitated through spatial nonuniformity (17, 19, 35)], supershear cracks may propagate in regions of strong interfaces, well below the critical τ0p values given by the homogeneous Burridge-Andrews transition mechanism. For example, the green and red simulated examples (Fig. 2A) correspond to τ0p = 0.3 and 0.2, respectively. Similarly, the particularly slow laboratory supershear rupture, CfEmbedded Image (blue example in Fig. 2B), corresponds to τ0p < 0.37, suggesting that, in this case, the transition has not been induced by the Burridge-Andrews mechanism but rather by favorable heterogeneities (19) at the block edge.

Minimal length for supershear propagation

The LEFM-based equation of motion describes two distinct shapes. For high prestress levels, for example, τ0p = 0.8 (blue example in Fig. 2A), Cf monotonically increases with increasing l. As τ0p is decreased, however, the solution bifurcates into two branches, for example, τ0p = 0.2 (red example in Fig. 2A), one of which appears to be unstable (dashed line in Fig. 2A). The accelerating branch (the physical solution) only appears at a finite crack length lm, predicting that no supershear cracks can exist for l < lm. Thus, lm represents a critical length for supershear propagation, similar to Griffith’s length for the sub-Rayleigh regime (31). We have confirmed this prediction for several values of τ0p by varying the length of the weak patch in our simulation (Fig. 3A). If the weak patch exceeds lm, then the crack obeys the supershear equation of motion in the following stronger region. However, in cases where the weak patch ends at l < lm, the crack transitions instantaneously to the sub-Rayleigh regime. From this point, the crack remains at the sub-Rayleigh regime even for l > lm, where the supershear solutions exist. This demonstrates the coexistence of the sub-Rayleigh and supershear propagation solutions and the hysteretic transition between the two.

Fig. 3 Minimal lengths exist for supershear crack propagation under uniform loading.

(A) Crack simulations for nonuniform τ0p profiles (see Fig. 2A, top), where the size of the weak patch is varied. Dashed and solid curves indicate propagation within the weak and strong regions, respectively. τ0p values in the figure correspond to the strong region (in the weak region, τ0p = 0.8). Supershear propagation in the strong region cannot be sustained for l < lm (lm is indicated by stars), where the supershear–to–sub-Rayleigh transition occurs. (B) lm obtained numerically (green line) from Eq. 1 shows that no minimal length exists for τ0p ≳ 0.45 (for fixed CS/CL). Normalization by l0 and τp is applied for convenient comparison to (A). Stars correspond to the values denoted in (A). The dashed line is the analytic approximation based on Eq. 2.

Figure 3B describes the dependence of lm on τ0. No critical length exists for sufficiently high stress levels. For τ0p ≲ 0.45, a finite value of lm emerges and lm increases rapidly with decreasing τ0. We find empirically that Cf(lm) ≈ Embedded Image; hence, by using Eq. 1, lm can be approximated byEmbedded Image(2)Figure 3B shows a comparison of Eq. 2 with the exact solution of Eq. 1. This approximation reveals that when it exists, the scale of lm is determined by the Griffith length for crack nucleation lc ~ μГ/Embedded Image (lm ≈ 3.7lc for CS/CL = 0.48). We emphasize that the minimal length presented here defines whether the supershear solution does or does not exist, which is fundamentally different from the question of what determines the sub-Rayleigh–to–supershear transition and whether it actually occurs (for example, Andrew’s transition length, which is, in a uniform setup, always larger than lm).

Supershear propagation along interfaces with nonuniform loading

Both changes in frictional strength (barriers) and nonuniform prestress (asperities) affect fault ruptures. These heterogeneities influence not only rupture propagation (17) but also produce different ground motions (36). In Fig. 4A, we show simulations that model nonuniform loading to illustrate the differences in the effect of barriers versus asperities for supershear propagation; changes in τ0p are now induced by varying τ0, in contrast to the simulations shown in Fig. 2A, where τ0p profiles were varied by spatial changes in τp and Г.

Fig. 4 Comparing theoretical predictions of supershear rupture velocities to numerical simulations for nonuniform prestress.

(A) Top: Spatially nonuniform τ0 profiles with uniform τp profiles are considered. The imposed τ0p profiles are shown. Bottom: Rupture velocities Cf(l) corresponding to the stress profiles in the top panel. Top and bottom: A high-stress nucleation patch (τ0p = 0.8 region) facilitates a direct supershear transition. Black solid lines show theoretical predictions for the nonuniform loading described by Eq. 3. Note the difference of Cf(l) profiles with those in Fig. 2A for similar τ0p profiles. Predicted locations for supershear–to–sub-Rayleigh transition are marked by gray triangles. (B) Comparison of supershear–to–sub-Rayleigh transition positions measured in 50 simulations, with the predicted transition position based on Eq. 3. In the considered τ0(x) profiles, the size of the high-stress nucleation patch is modified [see yellow curve in (A)]. Dashed line with slope 1 and is given for reference. Colors indicate the minimal needed increase in τ0 within the low-stress region that would negate the predicted supershear–to–sub-Rayleigh transitions. Gray triangles denote the simulations shown in (A).

Note that in the context of LEFM, spatial variations in τ0 are fundamentally different from those in τp (see Materials and Methods). Equation 1 assumed uniform prestress and cannot be strictly applied to nonuniform loading configurations. We therefore need to generalize Eq. 1 for spatially varying τ0. To accomplish this, the information about the τ0(x) profile should be incorporated within K. We assume K = κ(Cf/CL)Ks0, l, g) and hypothesize that Ks is a weighted integral functional of τ0(x) that can be calculated in a manner analogous to the static stress intensity factor used for sub-Rayleigh cracks (see Materials and Methods). This yieldsEmbedded Image(3)which coincides with Eq. 1 for uniform loading. Figure 4A shows predicted rupture speeds based on Eq. 3. Equation 3 generally compares well to simulations with minor discrepancies that appear in regions immediately following nonuniform areas. Note that even though the τ0p profiles appear to be very similar to the nonuniform strength configurations in Fig. 2A, the dynamics are significantly different, reflecting the influence of the spatial variation of τ0. For instance, while Eq. 2 predicts that at τ0p ≈ 0.15 no supershear cracks can exist for l/l0 < 200 (see Fig. 3), Eq. 3 predicts sustained supershear propagation at τ0p = 0.1 (see violet example in Fig. 4).

Only for extremely low background stress levels do we observe the transition to sub-Rayleigh propagation (see yellow and brown curves in Fig. 4A). At the crack length la of this transition, similar to the minimal supershear crack length lm (Fig. 3), the energy flux to the crack tip is insufficient to drive supershear cracks. Figure 4B shows a systematic comparison between predicted supershear to sub-Rayleigh transitions and measured results from simulations for 50 different prestress configurations. To evaluate the criticality of prestress configurations, we computed the additional background stress Δτ0 required to prevent the transition to sub-Rayleigh (shown in color in Fig. 4B) and sustain supershear propagation. Simulations with lower prediction accuracy have lower Δτ0 and hence exhibit an extreme sensitivity to values of the prestress driving the transition. The prediction works generally well, although low values of Δτ0p result in underpredicted la for la/l0 > 180. Thus, small changes in the background prestress, for instance due to wave radiation from the crack, which have been neglected in the current theoretical description, should have considerable effects on the rupture dynamics and could explain the prediction discrepancies.


Our measurements and simulations have shown that the energy balance at the rupture tip (Eqs. 1 and 3) provides quantitative predictions for the evolution of Cf(l). These results are general so long as several necessary assumptions are satisfied. For instance, a region in the rupture tip’s vicinity should exist where the stresses are singular τ ~ 1/(xl)g. While typically this is assumed to be true when lxc, this condition is hard to meet for τ0 → τp, where g → 0. Furthermore, in all rupture events considered here, shear waves that originate at the rupture nucleation (9, 18) trail behind the supershear rupture tip. However, one might consider more complex scenarios. Rough faults, for example, may lead to rapid acceleration and deceleration of sub-Rayleigh ruptures (35) that would result in stress wave radiation. If the sub-Rayleigh–to–supershear transition were to occur at later times, then the supershear rupture will eventually catch up with these waves, and our assumption of time-independent loading will be violated. Finally, in our experiments, there is no significant variation of Embedded Image after the passage of the rupture tip (37), and frictional ruptures can be mapped to simple cracks with traction-free faces. Slip pulses (15, 38), which may result because of strong velocity weakening of friction (39) and rapid healing of the interface, are fundamentally different from crack-like slip fronts and would require a different predictive theoretical model.

Natural faults are complex entities that include significant heterogeneity of fracture energy, friction laws, and stresses. Furthermore, complicated fault geometry, together with significant challenges in measuring the speeds of supershear earthquakes (11, 12, 14) and determining fault properties, often hinder comparisons with modeling. In some cases, qualitative and sometimes quantitative features of natural earthquakes can be successfully compared with idealized numerical models (13) and laboratory experiments (40). This fracture mechanics–based picture presented in this article is one of these, as it provides new and significant fundamental understanding of supershear rupture dynamics. These ideas may be potentially generalized to further account for the complexities associated with natural earthquakes.

Our results not only highlight how key parameters such as τ0, τp, and Г control supershear earthquakes but also provide a tool for quantitative predictions of supershear rupture speeds. These results have important implications, as the propagation speeds of supershear earthquakes dictate the amplitudes of the stress fields along the radiated Mach fronts (15). These amplitudes, while disappearing at Cf = Embedded Image, greatly increase as CfCL. These subtle changes of the supershear earthquake speeds, therefore, have significant impact on off-fault damage and provide a key parameter for predicting the serious seismic hazards that these earthquakes impose.

Finally, while supershear earthquakes are typically associated with highly stressed faults, our results support suggestions (19) that supershear earthquakes may exist even at extremely low stress levels, if favorable heterogeneities for supershear transition are present. The proposed equation of motion provides exact conditions for which these low stress earthquakes may or may not sustain supershear propagation and can explain supershear–to–sub-Rayleigh transitions such as those in (19). This work, furthermore, provides fundamental insights for understanding how the crucial interplay between fault roughness and stress levels (35) may govern supershear rupture.


Experimental system

Our experiments were conducted using PMMA plates (ρ ≈ 1170 kg/m3) of dimensions 200 mm × 100 mm × 5.5 mm (top block) and 240 mm × 100 mm × 5.5 mm (bottom block) in the x, y, and z directions, respectively (see Fig. 1A). The blocks were pressed together by applying ≈5 MPa of nominal pressure. The contacting surfaces were cleaned with distilled water and isopropyl alcohol and then dried for about 2 hours (termed here dry interfaces). We also conducted experiments in the boundary-lubrication regime (3, 34), where contacting surfaces were coated by a thin layer of lubricant (silicon oil with kinematic viscosity v ≈ 100 mm2/s). Material shear, CS, and longitudinal, CL, wave speeds were obtained by measuring the time of flight of ultrasonic pulses, yielding CS = 1345 ± 10 m/s and CL = 2700 ± 10 m/s, corresponding to the plane strain hypothesis (εzz = 0). Using these measured values, we calculated that for plane stress conditions (σzz = 0), CL ≈ 2333 m/s (CS/CL ≈ 0.58) and CR ≈ 1237 m/s. The corresponding dynamic shear modulus is μ = 2.1 GPa (9). The three components of the 2D strain tensor were continuously and simultaneously measured by Vishay 015RJ rosette strain gages every 1 μs at 16 to 19 spatial locations, ~3.5 mm above the frictional interface (2).

In previous work (2), it was found that the strains in the vicinity of rupture tips, propagating below CR, are well described by the square-root singular LEFM solutions, originally derived to describe brittle shear cracks. This comparison of the solutions with the measured strains provides a measure of the fracture energy Г, which was found to be roughly independent of rate. Furthermore, it was suggested that measurements of the contact area A variations provide a direct measurement of the cohesive zone size, xc (2). The theoretical interface strength τp is defined as the difference between the actual interface strength and the residual stress Embedded Image (see Fig. 1C). τp is obtained from measured values of Г and xc using the simplest cohesive zone models (29, 41). For the current set of experiments, we measured slight spatial inhomogeneities of the frictional interface. In the dry case, for 50 mm < x < 100 mm, Г ≈ 1.5 J/m2 and τp ≈ 1.3 MPa, while for 100 mm < x < 150 mm, Г ≈ 2.5 J/m2 and τp ≈ 1.6 MPa. This effective step increase of Г was used to compare our measurements to the predictions of Eq. 1 in Fig. 2B. In the boundary-lubrication regime, for 100 mm < x < 150 mm, Г ≈ 10 J/m2 and τp ≈ 2.66 MPa. Details can be found in the Supplemental Material of (3). In general, the measured values of Г and τp obtained for the sub-Rayleigh regime were used here for supershear ruptures.

Numerical simulations

Our numerical results were generated by solving elastodynamic equations with the spectral boundary integral method (32, 42). The propagation of a dynamic crack between two half spaces was modeled with a cohesive-type approach to describe the tractions along the weak interface. An explicit time integration was applied. The spectral formulation of the tractions and displacements at the interface results in a periodic setup. In our simulations, we used a replication length of 1.2 m with a discretization of 8192 nodes. The half-space material is linear elastic. To compare with the experiments, we applied a dynamic value of the elastic modulus E = 5.65 GPa, Poisson’s ratios v = 0.35, and density ρ = 1170 kg/m3 and used a plane strain assumption (CS/CL ≈ 0.48). The interface tractions were governed by a linear slip-weakening cohesive law, τ(δ) = τp(1 – δ/dc) for 0 < δ < dc, which imposes a strength that decreases linearly with slip δ from a peak value τp to zero over a characteristic slip weakening distance dc. Reference values applied in the uniform setup were τp = 1.0 MPa and dc = 2 μm, which led to a fracture energy of Г = 1.0 J/m2. Rupture nucleation was triggered via slowly propagating a seed crack of imposed velocity 0.1CR [following (9)] through the nucleation zone. In the nucleation zone, the value of τp was gradually reduced to zero (over a length of ≈6 mm). Once the seed crack reached a critical distance lc (the Griffith length), rupture acceleration initiated and the ruptures propagated dynamically.

While v = 0.35 occurs in serpentinized mantle material (43) as, for instance, along the Denali fault (44) that hosted a supershear earthquake (12, 13, 45), v = 0.25 is a common value for granite. Additional numerical results were obtained with v = 0.25 and plane strain boundary conditions (CS/CL = 0.577). Comparisons with the theoretical predictions are shown in fig. S1. The results are qualitatively identical to Fig. 2. The main quantitative difference is that the ruptures accelerate over shorter distances to considerably higher speeds, that is, Cf/CL > 0.9. Therefore, the differences in rupture speed for very different prestress levels are already relatively small shortly after transition.

Here, we used spatially nonuniform τ0p profiles. Direct supershear transition was triggered by using high τ0p levels for very short crack lengths. Supershear propagation was then studied as cracks entered into regions with lower values of τ0p (for example, brown curve in Fig. 5). We have also tested the quality of the theoretical prediction for supershear ruptures after transition through the well-known Burridge-Andrews mechanisms (orange curve in Fig. 5). The rupture first propagates in the sub-Rayleigh regime until a radiated shear wave ahead of the crack tip nucleates a secondary crack that then propagates at supershear speeds. A third simulated rupture was initiated by a seed crack that propagated at supershear speeds already during the nucleation procedure. All three tested supershear transition mechanisms lead to supershear ruptures with propagation speeds that are, after brief transient differences, quantitatively well described by our theoretical model. The observed discrepancies within this transient period, we believe, are related to the history dependence, discussed by Huang and Gao (24), which was neglected in the derivation of Eq. 1.

Fig. 5 The effect of sub-Rayleigh–to–supershear transition on the equation of motion of supershear cracks.

Three different transition mechanisms are considered. Top: The first setup (brown curve) has a spatially nonuniform τ0p profile with reduced local τp for l/l0 < 50 (see main text for details). Two additional examples have spatially uniform τ0p profiles (orange and red curves). Bottom: Colors represent the crack velocities Cf(l) corresponding to the stress profiles in the top panel. Brown curve indicates continuous crack acceleration to supershear speeds (direct transition) within a weakened nucleation (nucl.) patch (high τ0p level). Orange curve indicates sub-Rayleigh rupture transitions at l/l0 ≈ 65 to supershear speed through the Burridge-Andrews (BA) mechanism (4, 5). Red curve indicates an imposed supershear seed crack leads to a self-sustained supershear crack propagation. The black dashed line denotes theoretical prediction for a spatially uniform prestress level (τ0p = 0.45).

Theoretical framework

In what follows, we briefly describe how LEFM can be used to provide quantitative predictions for supershear rupture propagation following Broberg (26). The self-similar solution (26, 31) for a symmetrical bilateral singular crack expanding at a constant speed under uniform remote shear stress τ0 providesEmbedded Image(4)where l is the crack half length and κ(Cf/CL) a known function.

A nonrealistic consequence of the singular description is the vanishing energy flux into the crack tip for CfEmbedded Image. It was shown (31), however, that when a cohesive zone is introduced, a finite region where these singularities are regularized, the requirement for a positive energy flux to the crack tip is fulfilled for any Cf > CS. The energy flux is given byEmbedded Image(5)Here, μ is the shear modulus, τp (the peak shear strength) is defined in Fig. 1, xc is the cohesive zone size, and Embedded Image(Cf/CL) is a known function containing information about the shape of the stress distribution within the cohesive zone. Furthermore, as small-scale yielding is assumed (that is, lxc and τp ≫ τ0), the shear stress of these cohesive zone models recovers the singular form at large distances from the crack tip (31)Embedded Image(6)where κc(g) is another known function. The complete form of Eqs. 4, 5, and 6 are given, respectively, by Eqs. 52, 66, and 64 in (26), with slight changes in the notations: τp, xc, τ0, and g denote τD, d, Embedded Image, and γ.

G can be expressed in terms of τ0 and l, by comparing the two independent forms for τ (Eqs. 4 and 6). This intermediate asymptotic matching was also used by Huang and Gao (24) and Antipov et al. (28). In (26), G was determined as a function of Cf/CL, l, τ0, and xc [Eq. 68 in (26)] and given byEmbedded Image(7)where B(Cf/CL) and ГD(g) are known functions where the shape of the stress distribution within the cohesive zone is contained within ГD(g). Here, motivated by recent experimental results (2), which showed that τp is roughly independent of Cf (in contrast to the Cf dependence of xc), we rederived G in terms of Cf/CL, l, τ0, and τpEmbedded Image(8)Embedded Image(Cf/CL) and Embedded ImageD(g) are given byEmbedded Image(9)Embedded Image(10)where I0, N, f1, and ωD are given explicitly in (26) and can be calculated once Cf/CL and CS/CL are specified. Note that the cohesive zone properties are only contained within Embedded ImageD(g). D(ω) defines the shape of the stress distribution within the cohesive zone; for ξ = xxtip < 0, the shear stress gradually decreases according to a prescribed shear stress profile, τ(ξ) = τpD(ξ/xc). We chose Embedded Image to describe the experimental results. Our numerical simulations use a more simply implemented linear slip-weakening cohesive law τ(δ) (see “Numerical simulations” section). For simplicity, we modeled simulated ruptures by a linear spatial stress distribution within the cohesive zone through D(ξ/xc) = 1 + ξ/xc for –1 < ξ/xc < 0. Although some deviations between the models exist [for example, linear spatial profiles of D(ξ) result in nonlinear slip laws τ(δ), see Fig. 16 in (15)], we explicitly verified that these deviations have negligible effects on the resulting equation of motion.

Figure 6 shows g(Cf/CL) and the product of Embedded Image(Cf/CL) and Embedded ImageD(g) [Embedded ImageD(g) was calculated for a linear cohesive zone profile], so that the predictions of the supershear equation of motion (Eq. 1 in the main text) can be easily reproduced. At Cf = Embedded Image, the singularity exponent g = 1/2 and ГD(1/2) = Embedded ImageD(1/2) = 1. Therefore, Eqs. 7 and 8 coincide, and G(Embedded Image) is independent of the cohesive zone characteristics. Using these observations and explicitly verifying that Cf(lm) ≈ Embedded Image, the approximated form for a critical length for supershear propagation, lm, given by Eq. 2 in the main text, is easily obtained.

Fig. 6 Calculated g(β) and Embedded Image(Cf/CL)Embedded ImageD(g).

g(β) and the product of Embedded Image(Cf/CL) and Embedded ImageD(g), as appearing in Eq. 1, are provided for CS/CL = 0.48 (top) and CS/CL = 0.577 (bottom). Embedded ImageD(g) was calculated assuming a linear cohesive zone model.

A crucial difference between Eqs. 4 and 6 exists. While Eq. 6 adapts rather quickly to spatial changes of interface properties, Eq. 4 (and therefore Eq. 1) does not account for nonuniform τ0 due to the underlying assumption of the solution. In what follows, we consider supershear propagation along interfaces with nonuniform spatial profiles of τ0(x).

We assume for supershear cracks the condition of small-scale yielding, τ ~ K/(xl)g, where K is the stress intensity factor. The value of K explicitly incorporates the information about the τ0(x) profile. The fundamental solution for semi-infinite supershear cracks propagating at a constant speed and subjected to a pair of suddenly applied concentrated shear forces is given by Eq. 30 in (27)Embedded Image(11)For a general τ0(x) profile, K = Embedded Image(Cf/CL)Ks is obtained by integrating the fundamental solution following Freund’s superposition method (30), which results in Embedded Image(12)Note the analogy with semi-infinite sub-Rayleigh cracks, where Ks(l) ~ Embedded Image.

As far as we know, an explicit solution analogous to Eq. 12 (which was derived for semi-infinite cracks) for propagating supershear bilateral cracks in our simulations has not been formulated. Therefore, we suggest that the sub-Rayleigh static stress intensity factor Ks(l) ~ Embedded Image can be extended to supershear propagation byEmbedded Image(13)The normalization factor Embedded Images(g) = Embedded Image ensures that the generalized form of the shear stresses given byEmbedded Image(14)recovers Eq. 4 for spatially uniform τ0(x) profiles.

It is worth emphasizing again that Eqs. 11 and 12 are derived for cracks propagating at a constant speed. Previous work (24) has shown that a rigorous solution for accelerating cracks cannot be constructed by the superposition of these constant-velocity solutions, as had been performed previously for sub-Rayleigh propagation (23). In this sense, we hypothesized that these equations can be applied to problems with slowly varying crack speed. As the kernel in Eqs. 12 and 14 is singular, we assumed that the main contribution to the integral comes from the region sl. Therefore, assuming that Cf varies slowly, we calculated the integral with a value of g(Cf) that corresponds to the local crack speed at the crack position l. Finally, Cf(l) is determined numerically by Eq. 3: For any crack length l, Cf is evaluated so that Eq. 3 is satisfied. Although this is in no way a rigorous derivation and some of the approximations here are uncontrolled, our conjecture is in good agreement with the simulations and provides a useful tool for describing supershear cracks subjected to nonuniform loading.


Supplementary material for this article is available at

Fig. S1. Comparison of theoretical predictions of supershear crack velocities with numerical simulations for various shear strength levels for v = 0.25 (plain strain).

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank M. Adda-Bedia for fruitful discussions and his help. Funding: We also acknowledge the support of both the US-Israel Binational Science Foundation (grant no. 2016950) and Israel Science Foundation (grant no. 1523/15). Author contributions: D.S.K. performed the numerical work. I.S. performed the experimental and theoretical work. All authors contributed to the analysis and writing of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. The computer code used for simulations, as well as the additional data related to this paper, may be requested from the authors.
View Abstract

Navigate This Article