A minimal model for the onset of slip pulses in frictional rupture
Kjetil Thøgersen, Einat Aharonov, Fabian Barras, François Renard
AA minimal model for the onset of slip pulses in frictional rupture
Kjetil Thøgersen a , ∗ Einat Aharonov b , Fabian Barras a , and Fran¸cois Renard a,c The Njord Centre, Departments of Physics and Geosciences, University of Oslo, 0316 Oslo, Norway a Institute of Earth Sciences, The Hebrew University, Jerusalem, 91904, Israel b Universit´e Grenoble Alpes, Universit´e Savoie Mont Blanc,CNRS, IRD, IFSTTAR, ISTerre, 38000 Grenoble, France c (Dated: December 22, 2020)We present a minimal one-dimensional model for the transition from crack-like to pulse-like prop-agation of frictional rupture. In its non-dimensional form, the model depends on only two freeparameters: the non-dimensional pre-stress and an elasticity ratio that accounts for the finite heightof the system. The model contains stable slip pulse solutions for slip boundary conditions, andunstable slip pulse solutions for stress boundary conditions. Results demonstrate that the existenceof pulse-like ruptures requires only elastic relaxation and redistribution of initial pre-stress. Thenovelty of our findings is that pulse-like propagation along frictional interfaces is a generic elasticfeature, whose existence is not sensitive to particular rate- or slip-dependencies of dynamic friction. Frictional rupture, the process by which relative slidingmotion starts along the interface between two contactingsolid surfaces, is of prime importance in the description ofvarious systems in physical sciences and engineering. Ex-amples range from the squealing of brake pads and bushbearings [1, 2], corrugation and wear of mechanical com-pounds [3], to the earthquake cycle along crustal faults [4]and the surge of glaciers [5]. A frictional rupture is char-acterized by its rupture speed and its rupture mode. Therupture speed can vary from slow, through sub-Rayleighup to super-shear [6–8]. The rupture mode can be de-scribed by analogy to the dynamics of shear cracks, a be-havior referred to as crack-like dynamics [9]. Frictionalrupture can also be pulse-like rupture, where the slippingportion of the interface is spatially and temporally con-strained to the vicinity of the propagating rupture tip.Since the seminal work of Heaton [10] in 1990, thepulse-like rupture mode has been successfully used to ra-tionalize the short slip duration observed in seismic inver-sions and describe the propagation of both fast [11, 12]and slow [13] earthquakes. Pulse-like ruptures have alsobeen reproduced in laboratory experiments [14–18]. Theubiquity of pulse-like rupture observed in frictional sys-tems has motivated the development of theoretical andnumerical models to investigate the conditions support-ing the emergence of slip pulses. Slip-pulses have beenfound in systems with a large variety of boundary con-ditions and domain approximations [19–26] and frictionlaws including from rate-and-state friction [22, 23, 26],velocity-weakening friction [19, 20], slip-dependent fric-tion [27] and Coulomb friction [25].Several studies have demonstrated the important con-trol of the initial stress distribution along the inter-face, also called pre-stress, on slip pulses. Pulses widenwith increasing pre-stress, approaching crack-like rupture[17, 28]. Heterogeneous stress distributions can cause slippulses [29]. Stress barriers can cause cracks to arrestand create a backward propagating arresting front thatcauses a pulse in the opposite direction [30]. Also, fault geometry may play a role. Large earthquakes have beenshown to favour a transition from crack-like to pulse-likepropagation [28]. Large aspect ratios may favour the for-mation of slip pulses [31], and the seismogenic depth canlimit the width of earthquake slip pulses [32, 33]. Re-cent studies have also demonstrated the inherent insta-bility of slip-pulses for different friction constitutive laws[22, 24, 26, 34].Despite these extensive studies, the fundamental set ofconditions required to develop slip pulse in frictional sys-tems remains debated. Several studies have highlightedthe need of velocity-weakening friction to obtain slip-pulses, and there is a growing consensus that geometri-cal effects and initial pre-stress may play a role, but thatthe existence of slip-pulses requires either strong velocity-weakening friction [11, 16, 17, 35] or a coupling betweennormal stress and slip due to, for example, bimaterial ef-fects [16] or thermal pressurization of pore fluids alongthe frictional interface [36].In light of the many different models and interpre-tations of slip pulses, our goal is to identify the mini-mum ingredients necessary for their appearance. Here,we introduce a one-dimensional model which, in non-dimensional form, contains only two free parameters; thenon-dimensional initial shear pre-stress, ¯ τ , and a ratio ofelastic moduli, ¯ γ , that governs the elastic relaxation ofshear stress due to a finite system size. We demonstratethat the onset of slip pulses can be explained in full by theelastic relaxation or redistribution of initial pre-stress. Inour simplified model, the existence of slip pulses does notrely on a coupling between slip and normal stress, nor ve-locity weakening friction. The stability of the pulse-likepropagation is set by the boundaries. For stress bound-ary conditions, pulses are sensitive to small perturba-tions in the pre-stress along the interface, which controltheir shape. For displacement boundary conditions, themodel predicts a transition from crack-like propagationto steady-state pulses for constant pre-stress. a r X i v : . [ phy s i c s . g e o - ph ] D ec Loading Nucleation and propagation Arrest H
We solve one-dimensional elasticity by integrating overa small system thickness H , assuming an infinitely rigidsubstrate. In the case of earthquakes, H would corre-spond to the distance perpendicular to the fault beyondwhich elastic parameters do not vary, i.e. the thicknessof the damage zone. This approximation is useful in thesense that the system is primed for frictional rupture atlarge aspect ratios, for which observations have shown tofavour slip pulses. Note that a similar approach has beenused previously by e.g. Bouchbinder et al. [37]. Here, weuse two types of boundary conditions at the top surface,i) displacement boundary condition and ii) stress bound-ary condition. In non-dimensional form (SupplementalMaterial), linear elasticity with Amontons-Coulomb fric-tion for the part of the interface that is moving reducesto ¨¯ u = ∂ ¯ u∂ ¯ x − ¯ γ ¯ u + ¯ τ ± − ¯ β ∂ ˙¯ u∂ ¯ x , (1) while the static friction threshold in non-dimensionalform can be written as | ∂ ¯ u∂ ¯ x − ¯ γ ¯ u + ¯ τ | ≥
1. Here,¯ τ ± = τ /σ n ∓ µ k µ s − µ k (2)where σ n is the normal stress, µ s and µ k are the staticand dynamic friction coefficients, and τ is the initial shearstress along the interface.¯ γ = 2 Gλ + 2 G (3)where G is the shear modulus, λ is the first Lam´e pa-rameter. For ¯ γ > γ = 0 corresponds to a stress boundary condition.¯ x = x/X , ¯ u = ( u − u ) /U , ¯ t = t/T are the dimen-sionless position, slip, and time respectively. The char-acteristic length and time scales are chosen so that anon-dimensional front speed ¯ v c = 1 corresponds to thewave speed of the material; X = H , U = 2 H µ s p − µ k pλ +2 G (or U = H µ s p − µ k pλ +2 G ), T = H q ρλ +2 G . To limit oscillationson the scale of the discretization, we also introduce asmall damping term − ¯ β ∂ ˙¯ u∂ ¯ x where we use ¯ β = 10 − . Thetransition from dynamic to static friction occurs whenthe local slip velocity reaches zero. We solve the systemfor varying ¯ τ (¯ x ) and ¯ γ . It is worth noting that the fi-nite difference discretization of equation 1 is identical tothe dimensionless form of the classical Burridge-Knopoffspring block model under the additional constraint thatthe spatial discretization d x equals the thickness of thesystem d x = H (d¯ x = 1). THE CRACK-PULSE TRANSITION
We initialize the shear pre-stress with a (small) Gaus-sian function with maximum ¯ τ = 1 and standard devia-tion 1, so that rupture onset occurs in the middle of thedomain (FIG. 1). For a symmetric pre-stress distributionwith ¯ τ >
0, we observe bi-directional propagation of ei-ther a crack (¯ γ = 0) or a crack that transitions into twoslip pulses (¯ γ > τ > τ < γ = 0 and for ¯ γ > τ = 0 in one directionand ¯ τ < Steady-state slip pulse
For ¯ γ > τ >
0, the steady-state solution ofdynamic frictional rupture is pulse-like. In the limit of¯ β = 0, this steady-state solution for a slip pulse can be -1 -0.5 0 0.5 100.20.40.60.811.2 -100 -50 0 00.511.522.533.5-20 0 20 -0.500.51 FIG. 2.
Left : steady-state solution of the pulse velocity com-pared to the analytical solution for the shear pre-stress pro-file given in the inset (the Gaussian function with standarddeviation 1 in the center of the domain is used to nucleaterupture). The figure also shows the solution for non-zero ¯ β . Right : steady-state solution for pulses with ¯ γ = 0 that arekept in steady-state because they propagate in a region of¯ τ = 0. The shear pre-stress distribution at the onset of slipdetermines the shape of the pulse. calculated in closed form and the derivation is provided inthe Supplemental Material. The slip velocity is a cosine˙¯ u = 1 √ ¯ γ cos p ¯ γ (1 − ¯ τ )¯ τ ¯ x ! , (4)where the maximum pulse velocity is ˙¯ u max = √ ¯ γ , and thepulse width is ¯ W = π ¯ τ √ ¯ γ (1 − ¯ τ ) . To obtain this solutionwe also used the steady-state front velocity ¯ v c = √ − ¯ τ which has been obtained previously by Amundsen et al. [38]. FIG. 2 shows the steady-state pulse velocity for ¯ γ > γ = 0, which correspond to displacement andstress boundary conditions applied on the top surface,respectively. For displacement boundary condition at thetop surface (¯ γ > γ = 0, a unique solution of a steady-state pulse does notexist because under uniform stress ¯ L trans and ¯ T trans tendto infinity when ¯ γ →
0. However, under non-uniformpre-stress conditions, pulses can still nucleate at ¯ γ = 0.In this case, the pulse size will be entirely determined bythe rupture history which is controlled by the initial pre-stress distribution, and is only stable when propagatingin a region of ¯ τ = 0, i.e. the pulse grows when ¯ τ > τ < Crack-pulse transition for ¯ γ > For displacement boundary condition applied at thetop surface, we can estimate the transition from crack -101-10-50510 0 50 100-100-80-60-40-20020406080100 -101
FIG. 3.
Left : Spatio-temporal slip velocity using ¯ γ >
Right : Spatio-temporal slip velocity for ¯ γ = 0 with anon-symmetric pre-stress distribution (left graph). The figurecontains the definitions of the transition time ¯ T trans and thetransition length ¯ L trans . to pulse using the analytical steady-state solution for theslip pulse. The transition time ¯ T trans is defined by thetime it takes to reach a slip of ¯ u trans = τ ¯ γ at the pointof nucleation. If we assume that ∂ ¯ u∂ ¯ x is small, the equa-tion of motion for the nucleation point can be approx-imated as ¨¯ u (¯ t ) ≈ − ¯ γ ¯ u (¯ t ) + ¯ τ , which under the initialconditions ¯ u (0) = 0 and ˙¯ u (0) = 0 has the solution ¯ u (¯ t ) = ¯ τ ¯ γ (1 − cos( √ ¯ γ ¯ t )). ¯ T trans is found from ¯ u ( ¯ T trans ) = 2 ¯ τ ¯ γ where we use the first solution ¯ T trans = π √ ¯ γ . Using ¯ T trans ,we find ¯ L trans using the front speed ¯ v c :¯ L trans = 2 ¯ T trans ¯ v c = 2 π p ¯ γ (1 − ¯ τ ) . (5)FIG. 3 (left) shows an example of a crack-pulse transitionat ¯ γ >
0, in agreement with the predictions derived inthis section.
Crack-pulse transition for ¯ γ = 0 Pulses at ¯ γ = 0 are inherently unstable, and the pulseshape is not unique. Pulses only keep a steady shapeif propagating in a region where the available elasticstrain energy exactly corresponds to the dissipation inthe pulse, which for Amontons-Coulomb friction corre-sponds to ¯ τ = 0. This can be seen from equation 1 whichreduces to the one-dimensional wave equation when ¯ τ = 0and ¯ γ = 0. Any small perturbation in ¯ τ can cause thepulse to change its shape, vanish or transition to a crack.For more complicated friction laws or when solving thesystem in more than one dimension, it will in generalbe very unlikely to nucleate a pulse where dissipationexactly matches the available elastic strain energy. For¯ γ = 0, the onset of a steady-state pulse requires that thefollowing two conditions are met together. i) The pre-stress is such that an initial crack arrests in one directionand not the other. This requires an asymmetry in thepre-stress. ii) The pre-stress in the propagating directionreaches zero at some point. For ¯ γ = 0, L trans is deter-mined by the crack size when the crack arrests in onedirection. The position of front arrest can be approxi-mated by noting that stress balance at zero accelerationrequires that ¯ τ = 0 in the ruptured region. Assuming nu-cleation of rupture is located at ¯ x = 0 and the arrestingoccurs in the negative direction, we can write Z ∞ x - ¯ τ ( x ) = 0 (6)which can be used to determine ¯ x - for any ¯ τ (¯ x ). Thetransition time is given by the time it takes to propagatefrom ¯ x = 0 to ¯ x - , which can be found if the rupture speed¯ v c is known. ¯ T trans = | ¯ x - |h ¯ v c i - (7)The total crack size is¯ L trans = | ¯ x - | + | ¯ x + | (8)where ¯ x + = Z ¯ T trans ¯ v c, + d¯ t (9)In the arresting direction, ¯ τ has to be negative. Then, wedo not have an analytical expression for the front speed,but we instead use h ¯ v c i ± ≈ T trans ∼ | ¯ x − | .and ¯ L trans ∼ | ¯ x − | , where ¯ x − is found from equation 6.Note that this will lead to a slight underestimation of¯ T trans . FIG. 3 (right) shows an example of a crack-pulsetransition at ¯ γ = 0, in agreement with the predictionsderived in this section.We have calculated the scaling relations for ¯ γ > γ = 0 over a wide range of ¯ γ and ¯ τ . While the non-dimensional elasticity ratio ¯ γ is expected to be found ina limited range close to 2 / γ = 0, we systematicallyvary the initial width ¯ W ¯ τ of the pre-stress perturbationto obtain different transition lengths and times. FIG. 4shows the data collapse using the scaling relations de-rived above. The model and analytical results are ingood agreement for the prediction of the transition fromcrack-like to pulse-like propagation. DISCUSSION
The present study introduces a one-dimensional two-parameter model for the onset of slip-pulses which can -0.8 -0.45 101520-0.8 -0.45 10152010 -3 -1 -2 -1 -3 -1 -2 -1 FIG. 4. Scaling relations of ¯ L trans and ¯ T trans for the crack-pulse transition for both ¯ γ > γ = 0 (rightpanels). For ¯ γ = 0, we use an asymmetric stress distribution(top right inset), while for ¯ γ = 0 we use a symmetric stressdistribution (top left inset). The face color shows the stress¯ τ or ¯ τ ,l , while the color of the rim shows ¯ γ or ¯ W ¯ τ . be understood from elastic relaxation and redistributionof pre-stress only. Therefore, in this minimal model, slippulse requires neither a coupling between slip and normalstress nor velocity-weakening friction.In our model, the onset of slip-pulse frictional ruptureoccurs through a transition from crack-like to pulse-likerupture at a crack size that depends on pre-stress andboundary conditions. For stress boundary conditions,the transition length is controlled in full by the initialpre-stress distribution, and the transition occurs becausethe crack arrests in one direction and not the other. Thepre-stress in the propagating direction is set so that dissi-pation in the pulse is equal to the available elastic strainenergy. For slip boundary conditions, the transition fromcrack to a self-similar pulse is given by the system width H , the elasticity ratio ¯ γ as well as the ratio of rupturespeed v c and speed of sound v s ; L trans = 2 πH q λ +2 G G v c v s .The pulse width is W = π ¯ τ H q λ +2 G G v c v s with a maximumvelocity ˙ u max = ∆ τ √ Gρ , where ∆ τ = ( µ s − µ k ) σ n is thestress drop. The model predicts that the onset of slipoccurs at shorter propagation lengths in systems withlarge aspect ratios, and low initial pre-stress. Applied toearthquake mechanics, these results are consistent withi) the observation that large earthquakes favor slip-pulses[28], ii) the observation that large aspect ratios seem tofavor the formation of slip pulses [31], iii) the observationthat increasing pre-stress widens slip-pulses so that theyapproach crack-like solutions [17], and iv) the predictionthat the seismogenic depth may limit the width of slippulses [32, 33].Recent studies have demonstrated the inherent insta-bility of slip-pulses for different kinds of friction consti-tutive laws [22, 24, 26, 34]. These models reported theexistence of a steady-state pulse solutions located at thesharp transition between growing pulses, whose spatialextent increases during propagation, and decaying pulses,whose spatial extent progressively shrinks and eventu-ally arrests the pulse. The slip pulses in our simplifiedmodel have the same intrinsic instability only when stressboundary conditions are used (¯ γ = 0). Then, slip pulsesonly keep their shape when the available elastic strainenergy in the form of pre-stress is exactly equal to thedissipation in the pulse. For our choice of Amontons-Coulomb friction law, this corresponds to ¯ τ = 0, whenthe equation of motion reduces to the one-dimensionalwave equation. In this case, the net stress change is zeroif the pulse velocity profile is constant, and the net stresschange is non-zero only if the pulse accelerates or de-celerates. This is consistent with the findings of; pulsescan propagate rapidly with negligible net stress change[40]. Any small perturbation in ¯ τ will make the slip pulsedecrease or grow. The velocity profile of these unstablepulses at ¯ γ = 0 are not unique in our model. This non-uniqueness is likely a feature of our specific choice of fric-tion law. For different friction laws, unique solutions maybe likely found through couplings between sliding veloc-ity, slip and pulse dissipation. Interestingly, a slip bound-ary condition (¯ γ >
0) stabilizes the slip pulses throughelastic relaxation. Steady-state slip pulses are then foundfor any pre-stress ¯ τ ∈ (0 , τ willsimply transition pulses from one steady-state solutionto another.In view of the number of features of slip pulses cap-tured by the minimal model presented here, we proposethat pulse-like propagation along frictional interfaces isa likely generic feature. While different friction lawswill yield different velocity profiles, crack-pulse transi-tion lengths, propagation speeds and pulse widths, theexistence of slip pulses does not rely on details in thefriction constitutive law in our model. Instead, the pre-stress and the boundary conditions alone can produce slippulse solutions along frictional interfaces. This interpre-tation has important consequences. In particular, in lab-oratory experiments and natural earthquakes, one shouldbe careful with interpreting observations of slip pulses asproof of velocity weakening friction. Preventing the ap-pearance of slip pulses is in many situations desirableas they are the origin of squealing noise in train breaksand several industrial processes [1, 2]. Our findings alsodemonstrate the prospect of controlling the onset of slippulses as well as their stability through i) manipulation of the pre-stress at the frictional interface and ii) controlof the system boundary conditions. ∗ [email protected][1] N. Kinkaid, O. O’Reilly, and P. Papadopoulos, Journalof Sound and Vibration , 105 (2003).[2] G. Chen, J. Lv, Q. Zhu, Y. He, and X. Xiao, Wear , 1788 (2019), 22nd International Conferenceon Wear of Materials.[3] K. Kato, Wear , 151 (2000).[4] C. H. Scholz, Nature , 37 (1998).[5] K. Thøgersen, A. Gilbert, T. V. Schuler, and A. Malthe-Sørenssen, Nature communications , 1 (2019).[6] K. Thøgersen, H. A. Sveinsson, D. S. Amundsen,J. Scheibert, F. Renard, and A. Malthe-Sørenssen, Phys-ical Review E , 043004 (2019).[7] O. Ben-David, G. Cohen, and J. Fineberg, Science ,211 (2010).[8] J. K. Trømborg, H. A. Sveinsson, J. Scheibert,K. Thøgersen, D. S. Amundsen, and A. Malthe-Sørenssen, Proceedings of the National Academy of Sci-ences , 8764 (2014).[9] I. Svetlizky and J. Fineberg, Nature , 205 (2014).[10] T. H. Heaton, Physics of the Earth and Planetary Inte-riors , 1 (1990).[11] K. Chen, J.-P. Avouac, S. Aati, C. Milliner, F. Zheng,and C. Shi, Nature Communications , 1 (2020).[12] J. Galetzka, D. Melgar, J. F. Genrich, J. Geng, S. Owen,E. O. Lindsey, X. Xu, Y. Bock, J.-P. Avouac, L. B. Ad-hikari, B. N. Upreti, B. Pratt-Sitaula, T. N. Bhattarai,B. P. Sitaula, A. Moore, K. W. Hudnut, W. Szeliga,J. Normandeau, M. Fend, M. Flouzat, L. Bollinger,P. Shrestha, B. Koirala, U. Gautam, M. Bhatterai,R. Gupta, T. Kandel, C. Timsina, S. N. Sapkota, S. Ra-jaure, and N. Maharjan, Science , 1091 (2015).[13] S. Michel, A. Gualandi, and J.-P. Avouac, Nature ,522 (2019).[14] T. Baumberger, C. Caroli, and O. Ronsin, Physical re-view letters , 075509 (2002).[15] H. Shlomai and J. Fineberg, Nature communications ,1 (2016).[16] H. Shlomai, D. S. Kammer, M. Adda-Bedia,and J. Fineberg, Proceedings of the Na-tional Academy of Sciences , 27 (2010).[18] G. Lykotrafitis, A. J. Rosakis, and G. Ravichandran,Science , 1765 (2006).[19] J. Carlson and J. Langer, Physical Review A , 6470(1989).[20] J. M. Carlson, J. S. Langer, and B. E. Shaw, Reviews ofModern Physics , 657 (1994).[21] A. E. Elbanna and T. H. Heaton, Geophysical JournalInternational , 1797 (2012).[22] E. A. Brener, M. Aldam, F. Barras, J.-F. Molinari,and E. Bouchbinder, Physical review letters , 234302(2018).[23] G. Zheng and J. R. Rice, Bulletin of the SeismologicalSociety of America , 1466 (1998). [24] N. Brantut, D. I. Garagash, and H. Noda, Journal ofGeophysical Research: Solid Earth , 8998 (2019).[25] S. Nielsen and R. Madariaga, Bulletin of the Seismolog-ical Society of America , 2375 (2003).[26] A.-A. Gabriel, J.-P. Ampuero, L. A. Dalguer, and P. M.Mai, Journal of Geophysical Research: Solid Earth (2012).[27] Z. Wu and Y. Chen, (1998).[28] D. Melgar and G. P. Hayes, Geophysical Research Letters , 9691 (2017).[29] S. M. Day, G. Yu, and D. J. Wald, Bulletin of the Seis-mological Society of America , 512 (1998).[30] E. Johnson, Geophysical Journal International , 125(1990).[31] S. B. Nielsen, J. Carlson, and K. B. Olsen, Journal ofGeophysical Research: Solid Earth , 6069 (2000).[32] J. P. Ampuero and X. Mao, Fault zone dynamic pro-cesses: Evolution of fault properties during seismic rup-ture , 243 (2017). [33] K. Bai and J.-P. Ampuero, Journal of Geophysical Re-search: Solid Earth , 10 (2017).[34] H. Noda, E. M. Dunham, and J. R. Rice, Journal ofGeophysical Research: Solid Earth (2009).[35] N. Beeler and T. Tullis, Bulletin of the Seismological So-ciety of America , 1130 (1996).[36] D. Garagash, Journal of Geophysical Research: SolidEarth (2012).[37] E. Bouchbinder, E. A. Brener, I. Barel, and M. Urbakh,Phys. Rev. Lett. , 235501 (2011).[38] D. S. Amundsen, J. K. Trømborg, K. Thøgersen,E. Katzav, A. Malthe-Sørenssen, and J. Scheibert, Phys-ical Review E , 032406 (2015).[39] S. Ji, S. Sun, Q. Wang, and D. Marcotte, Journal ofGeophysical Research: Solid Earth (2010).[40] G. C. McLaskey, B. D. Kilgore, and N. M. Beeler, Geo-physical Research Letters , 7039 (2015). upplemental material: A minimal model for the onset of slip pulses in frictionalrupture Kjetil Thøgersen a , ∗ Einat Aharonov b , Fabian Barras a , and Fran¸cois Renard a,c The Njord Centre, Departments of Physics and Geosciences, University of Oslo, 0316 Oslo, Norway a Institute of Earth Sciences, The Hebrew University, Jerusalem, 91904, Israel b Universit´e Grenoble Alpes, Universit´e Savoie Mont Blanc,CNRS, IRD, IFSTTAR, ISTerre, 38000 Grenoble, France d (Dated: December 22, 2020) a r X i v : . [ phy s i c s . g e o - ph ] D ec ˆ x (0)
EQUATION OF MOTION
To build the minimal 1D model, we start with the 2D system presented in the FIG. S1. We start with the two-dimensional problem, following the approach of ? ]. We consider a long linear elastic strip with height H in contactwith an infinitely rigid substrate. H is assumed to be much smaller than any length scale over which all otherquantities vary in the x-direction. The equation of motion is ρ ∂ ~u∂t = ∇ · σ , (1)where the stress tensor σ is given by σ = λ tr( ∇ ~u ) I + G [ ∇ ~u + ( ∇ ~u ) T ] . (2) G is the shear modulus, λ is the second Lam´e coefficient, and ~u is the displacement. In this 2D system, the boundaryconditions are described in FIG. S1. At the bottom interface, the shear stress is set by the frictional force f f . On thetop surface, two kinds of boundary conditions are considered, either imposed shear stress or imposed displacement.To obtain a one-dimensional approximation, we integrate the equation of motion over a finite thickness H , corre-sponding to the red rectangle in FIG. S1: Hρ h ∂ ~u∂t i y = Z H ∇ · σ d y = ∂∂x Z d x Z H ∇ · σ d x d y = ∂∂x I S σ ~n d S (3)= ∂∂x (d x σ ~n | y = H + d x σ ~n | y =0 + H h σ ~n | x =d x i y + H h σ ~n | x =0 i y ) (4)= σ ~n | y = H + σ ~n | y =0 + ∂∂x H [ h σ xx i y , h σ xy i y ] T (5)Here, S denotes the surface of the volume over which the integration is calculated, and ~n is the normal vector. Theboundary conditions can now be set through the surface tractions [ˆ σ x ( H ) , ˆ σ y ( H )] T ≡ σ ~n | y = H and [ˆ σ x (0) , ˆ σ y (0)] T ≡ σ ~n | y =0 . We further assume u y = 0 which results in the last term being reduced to ∂∂x H [ σ xx , σ x ( x, , t ) is set by the friction law.At the top interface, we either set a constant traction, or a constant displacement (i.e. no slip boundary).Now we insert for u x ( x, t ) and assume u y ≈ λ and G , and we obtain Hρ ∂ h u x i y ∂t = H ∂ h σ xx i y ∂x + ˆ σ x ( x, H, t ) − ˆ σ x ( x, , t ) , (6)where we have introduced a minus sign on ˆ σ x ( x, , t ) in order to consider the f f as positive in equation 11. Insertingfor σ xx , we obtain Hρ ∂ h u x i y ∂t = H ( λ + 2 G ) ∂ h u x i y ∂x + ˆ σ x ( x, H, t ) − ˆ σ x ( x, , t ) . (7)The friction law at the interface can be applied through ˆ σ x (0 , t ), and ˆ σ x ( H, t ) is the driving stress. y = 0 defines thefrictional interface location. We use two variations of the boundary condition ˆ σ x ( x, H, t ).First, we follow what was used by ? ] and apply a shear stress boundary condition. ˆ σ x ( x, H, t ) is then set as aconstant function τ ( x ) that does not vary in time. The equation of motion can then be written as Hρ ∂ h u x i y ∂t = H ( λ + 2 G ) ∂ h u x i y ∂x + τ − ˆ σ x ( x, , t ) . (8)To be able to apply a friction law at y = 0, we need to make an additional approximation that relates the averageslip with the slip at the interface. When we later will apply Amontons-Coulomb friction at y = 0, it is sufficient toassume sign( ˙ u ( y = 0)) = sign( h ˙ u i ). Alternatively, for velocity-dependent friction laws an alternative approach wouldbe to approximate ˙ u ( y = 0) ≈ h ˙ u i .Second, it is also possible to introduce a constant displacement (i.e. no slip) boundary condition at y = H ; ∂u x ( H ) ∂t = 0. This can be rewritten in terms of a driving stress which is set by a relative displacement at the topboundary and the frictional interface u x ( x, H, t ) − u x ( x, , t ). Then, we need to introduce an approximate relationbetween the displacement at y = H and the displacement at y = 0, which we do by linearizing the displacementprofile. The shear stress ˆ σ x ( x, H, t ) can then be rewritten in terms of the imposed displacement at the top boundary,matching the traction and the shear stress in the volume close to the interface G ∂u x ∂y | y = H . Hρ ∂ h u x i y ∂t = H ( λ + 2 G ) ∂ h u x i y ∂x + u x ( x, H, − u x ( x, , t ) H G − ˆ σ x ( x, , t ) . (9)To evaluate the relative slip of the frictional interface and the top layer we can no longer use the averaged value h u x i y ,but instead change the slip of interest to u x ( x, , t ).12 Hρ ∂ u x ( x, , t ) ∂t = H λ + 2 G ∂ u x ( x, , t ) ∂x − H λ + 2 G ∂ u x ( x, H, ∂x + u x ( x, H, − u x ( x, , t ) H G − ˆ σ x ( x, , t ) . (10)where ∂ u x ( x,H, ∂x is a constant function that arises if the the displacement at the top interface is not set as constant.At y = 0, we apply a frictional boundary condition. In principle, this means that we adopt the linearization in anapproximate sense. i.e. small deviations from linearity can cause a difference in shear stress between y = 0 and y = H . The mathematical formulation is also equivalent to treating the y-coordinate as quasi-static, where therelevant difference in shear stress reduces to the difference between the shear stress at y → y = H ), and the kinetic shear stress at the boundary ( y = 0). As our friction law we choose an Amontons-Coulombfriction law f f = ˆ σ x ( x, , t ) at the interface, which can be written as f f ( ≤ µ s σ n if ˙ u = 0= µ k σ n if | ˙ u | > σ n is the normal stress (which equals σ yy ), with the additional criteria that a transition from static to dynamicfriction occurs when the local shear stress is greater than µ s σ n , and a transition from dynamic to static friction occurswhen the local velocity reaches zero. In the following, we will write u x ( x, , t ) as u . Dimensionless formulation
For the rest of the document, and in the main text, we drop the subscript x . We start with the equation ofmotion using a displacement boundary condition at y = H . Following the approach from ? ] We assume aninitial displacement field u ( x, , u ( x, H, u = u ( x, , t ) − u ( x, , τ :¨ u = λ + 2 Gρ ∂ u ∂x − GρH u − ρH ˆ σ x ( x, , t ) + 2 ρH τ (12)We then define the dimensionless variables ¯ u = u U , ¯ t = tT and ¯ x = xX so that¨¯ uU/T = λ + 2 Gρ UX ∂ ¯ u∂ ¯ x − GUρH ¯ u − ρH ˆ σ x ( x, , t ) + 2 ρH τ, (13)where the derivative is now taken with respect to ¯ t and ¯ x . We can further manipulate the expression to arrive at¨¯ u = λ + 2 Gρ T X ∂ ¯ u∂ ¯ x − GT ρH ¯ u − T U ρH ˆ σ x ( x, , t ) + T U ρH τ, (14)Next, we reduce the number of parameters by selecting λ +2 Gρ T X = 1, which leads to T = X q ρλ +2 G . This automaticallyobeys dimensionless wave speed equal to 1 in the limit of small GT ρH ;¯ v s = s λ + 2 Gρ TX = s λ + 2 Gρ r ρλ + 2 G = 1 . (15)Choosing U = 2 H µ s p − µ k pλ +2 G leads to T ρU H ( τ − f f ) = X H τ − ˆ σ x µ s σ n − µ k σ n . (16)If we now select the length scale X = H , we obtain the familiar parameter [ ? ? ] T ρU ( τ − ˆ σ x ) = τ /σ n ∓ µ k µ s − µ k ≡ ¯ τ ± (17)which applies when the interface is sliding, and where ± corresponds to the sign of the velocity. In this paper onlypositive velocities appear, which means we can simplify to ¯ τ ≡ ¯ τ + . We can then define the non-dimensional equationof motion ¨¯ u = ∂ ¯ u∂ ¯ x − ¯ γ ¯ u + ¯ τ . (18)where ¯ γ = 2 GT ρH = 2 GH ρρH ( λ + 2 G ) = 2 G ( λ + 2 G ) . (19)To remove oscillations at the grid scale we introduce a damping term ¯ β ¨¯ u = ∂ ¯ u∂ ¯ x − ¯ γ ¯ u + ¯ τ − ¯ β ∂ ˙¯ u∂ ¯ x . (20)Equation 20 applies for the sliding part of the interface, and dynamic friction is included through ¯ τ . The staticfriction threshold in dimensionless units reduces to | ∂ ¯ u∂ ¯ x − ¯ γ ¯ u + ¯ τ | ≥ u = ∂ ¯ u∂ ¯ x + ¯ τ . (22)where the dimensionless scaling parameters are the same as before apart from U which changes by a factor two dueto the different way of integrating the thin layer for the two different boundary conditions: U = H µ s p − µ k pλ +2 G .In the manuscript, these two variations as ¯ γ > γ = 0, refer to the two types of boundary condition applied atthe top surface of the domain. It is also worth noting that the finite difference scheme of the dimensionless equationof motion derived here is identical to the classical Burridge-Knopoff spring block model with leaf springs (see e.g. ? ]) under the assumption that the discretization ∆ x equals the thickness of the system ∆ x = H . ANALYTICAL SOLUTION FOR STEADY STATE SLIP PULSE ( ¯ β = 0 ) To obtain the steady state solution for a slip pulse we solve¨¯ u = ∂ ¯ u∂ ¯ x + ¯ τ − ¯ γ ¯ u (23)under the assumptions ˙¯ u ( − ¯ W /
2) = ˙¯ u ( ¯ W /
2) = 0, ¯ u ( − ¯ W /
2) = ∆¯ τ / ¯ γ , and ¯ u ( ¯ W /
2) = 0, where ¯ W = WH is thedimensionless pulse width. For a constant front speed ¯ v c we can relate ¯ u to ˙¯ u and ¨¯ u :˙¯ u = − ¯ v c ∂ ¯ u∂ ¯ x (24)where the minus sign arises from the assumption of propagation in the positive direction. Then,¨¯ u = ¯ v c ∂ ¯ u∂ ¯ x , (25)which leads to ¯ v c ∂ ¯ u∂ ¯ x = ∂ ¯ u∂ ¯ x + ¯ τ − ¯ γ ¯ u, (26)which can be rewritten as ∂ ¯ u∂ ¯ x (cid:0) ¯ v c − (cid:1) + ¯ γ ¯ u − ¯ τ = 0 . (27)The general solution is ¯ u (¯ x ) = c sin √ ¯ γ ¯ x p ¯ v c − ! + c cos √ ¯ γ ¯ x p ¯ v c − ! + ¯ τ ¯ γ (28)To find the constants c , c , as well as the pulse width ¯ W and the stress drop ∆¯ τ , we use the following boundaryconditions: ¯ u ( ¯ W /
2) = 0, ¯ u ( − ¯ W /
2) = ∆¯ τ ¯ γ , where ∆¯ τ is the stress drop. ˙¯ u ( ¯ W /
2) = 0, ˙¯ u ( − ¯ W /
2) = 0. In addition,we set ¯ v c = 1 / √ − ¯ τ which is exact in the limit ¯ γ → ? ], but is also a good approximation for ¯ γ ’ u (¯ x ) = ¯ τ ¯ γ − sin( p ¯ γ (1 − ¯ τ )¯ τ ¯ x ) ! (29)for ¯ x ∈ [ − ¯ W / , ¯ W / W is the pulse width¯ W = π s ¯ v c − γ = π ¯ τ p ¯ γ (1 − ¯ τ ) (30)The pulse velocity profile is ˙¯ u (¯ x ) = 1 √ ¯ γ cos p ¯ γ (1 − ¯ τ )¯ τ ¯ x ! (31)where the maximum pulse velocity is ˙¯ u max = 1 √ γ . (32)The acceleration profile is ¨¯ u (¯ x ) = 1¯ τ sin p ¯ γ (1 − ¯ τ )¯ τ ¯ x ! (33)In addition, the stress drop is given by ∆ τ = 2¯ τ (34)which can also be understood from symmetry arguments. FIG. S2 summarizes the pulse solution in steady state.FIG. S3 shows pulse solutions for ¯ γ = 0 expanding on the data shown in the main text. -1 -0.5 0 0.5 100.10.20.30.40.50.60.70.80.91 -1 -0.5 0 0.5 100.10.20.30.40.50.60.70.80.91 -1 -0.5 0 0.5 1-1-0.8-0.6-0.4-0.200.20.40.60.81 FIG. S2. Steady state solution of slip, slip velocity, and acceleration for ¯ γ > Left : Non-dimensional slip in the pulse insteady state.
Middle : Non-dimensional slip velocity in the pulse in steady state.
Right : Non-dimensional acceleration insidethe propagating pulse in steady state. -100 -50 0050100150 -100 -50 000.511.522.533.5 -100 -50 0-0.100.10.20.30.40.5 -20 0 20 -0.500.51
FIG. S3. Steady state slip, slip velocity and acceleration/stress for a selection of pulses at ¯ γ = 0 with prestress shown in theinset. Left : Non-dimensional slip in the pulse in steady state.
Middle : Non-dimensional slip velocity in the pulse in steadystate.
Right : Non-dimensional acceleration inside the propagating pulse in steady state.
SCALING RELATIONSHIPS FOR SYSTEMATIC VARIATION OF THE DAMPING PARAMETER ¯ β In the main text, we present a set of scaling relationships for the transition from crack-like to pulse like rupture.FIG. S4 shows how the data collapse for ¯ γ > β . FIG. S5 shows how the datacollapse for ¯ γ = 0 is affected by the damping parameter ¯ β . -3 -1 -2 -1 -3 -1 -2 -1 -3 -1 -2 -1 -3 -1 -2 -1 -3 -1 -2 -1 -3 -1 -2 -1 FIG. S4. Scaling relationships for the transition from crack-like to pulse-like rupture for a variety of damping parameters ¯ β for¯ γ > -0.8 -0.45 10 15 20-0.8 -0.45 10 15 20 -0.8 -0.45 10 15 20-0.8 -0.45 10 15 20 -0.8 -0.45 10 15 20-0.8 -0.45 10 15 20 FIG. S5. Scaling relationships for the transition from crack-like to pulse-like rupture for a variety of damping parameters ¯ β for¯ γγ