Rupture dynamics along bimaterial interfaces: a parametric study of the shear-normal traction coupling

perturbations

is, the direction of slip in the more compliant medium (Harris & Day 2005).Asymmetric distribution of near surface damage, as observed both at the San Andreas Fault (Dor et al. 2006) and at the Anatolian fault (Dor et al. 2008), were shown to result from the different propagation along the two sides of bimaterial rupture.Asymmetric distribution of the aftershocks after a main event in the vicinity of a segment of the San Andreas Fault were also evidenced and interpreted as the signature of bimaterial rupture (Rubin 2002;Rubin & Gillard 2000;Rubin & Ampuero 2007).
Laboratory experiments also evidenced an asymmetry in the rupture propagation along bimaterial interfaces.The rupture propagation along the favoured direction was shown to tend towards an asymptotic speed intermediate between the two Rayleigh wave speeds characterizing the materials on both sides of the interface, when the contrast in shear wave velocity is not too large (Xia et al. 2005).This asymptotic speed is the generalized Rayleigh (GR) wave speed C gr characteristic of an interfacial wave that propagates along a frictionless bimaterial interface constrained to not open.For some pre-stress conditions, supershear acceleration may additionally occur along the non-favoured direction (Xia et al. 2005).This behaviour was ascribed to the slip-induced normal traction perturbations occurring during the rupture propagation along bimaterial interfaces.Sliding under variable normal pressure reveals a delayed shear traction response to sudden normal traction variations that can be described as a function of slip (Prakash & Clifton 1993;Prakash 1998).
When considering sliding along interfaces separating dissimilar linearly elastic materials, there is a coupling between interfacial slip and normal traction variations due to broken symmetry.The crack tip asymptotic fields for frictionless bimaterial interfaces were analytically studied using the classical formalism of Mushkelishvili (Muskhelishvili 1953), by for example, Williams (1959), Erdogan (1965) and Rice (1988), and of Stroh (Stroh 1962), by, for example, Suo (1990) and Yang et al. (1991).The traction ahead the crack tip, the displacement behind the crack tip and the displacement discontinuity were explicitly established, showing near-crack tip oscillatory behaviour (Williams 1959) associated with a complex stress intensity factor.An arbitrary length scale needs therefore to be introduced to define classical, real stress intensity factors (Rice 1988).For frictional bimaterial crack, using classical Coulomb friction, the oscillatory part disappears but material contrast on both sides of the interface modifies the decay of the traction in the vicinity of the crack tip (Deng 1994).A singularity weaker or stronger than the classical square-root behaviour was argued depending on the slip at crack tip (Deng 1994).A singularity stronger than the square-root singularity was discarded as it would imply a backward propagation of the rupture (Audoly 2000).
Steady propagation of slip pulse rupture mode along bimaterial frictional interfaces was also analytically investigated.For constant friction coefficient, finite-size pulses can propagate only along the favoured direction due to the slip-induced normal traction perturbations, even when the remotely imposed shear traction is smaller than the shear frictional strength of the interface (Weertman 1980).The speed of the propagating pulse was shown to be C gr , as long as the generalized Rayleigh wave exists (Weertman 1980).
When the velocity contrast is such that the generalized Rayleigh wave exists, sliding along bimaterial interface under classical Coulomb friction law is not only unstable against perturbations at all wavelengths, irrespective of the value of the friction (Martins & Simões 1995;Martins et al. 1995;Simões & Martins 1998) but it is also mathematically ill-posed due to the lack of a length (or time) scale associated with the coupling between sliding and normal traction perturbations.Steady-state slip waves are dynamically destabilized along the bimaterial interface, due to self-excited motion for a wide range of elastic contrasts and frictional conditions (Adams 1995(Adams , 1998)).The stability and the ill-posed nature of the physical problem were also investigated analysing the slip-rate response to a single-mode perturbation of the shear traction (Ranjith & Rice 2001).Unstable growth of the interfacial modes was shown to affect all wavelengths, with a growing rate inversely proportional to the wavelength, which depends on the material contrast and on the friction coefficient.
Linear stability analyses of propagating slip modes (Ranjith & Rice 2001;Rice et al. 2001;Adda-Bedia & Ben Amar 2003;Ranjith 2009Ranjith , 2014;;Brener et al. 2015) were also performed considering classical Coulomb and rate-and-state-dependent friction laws, after a regularization based on the laboratory experiments of Prakash & Clifton (1993, 1998), which introduces a monotonic memory dependence (delayed response) of the shear strength on sliding-induced normal traction perturbations.The main results of these analyses are the destabilization of different interfacial slip modes (Love, Stoneley and slip waves) and their stability criteria as a function of the frictional parameters, materials dissimilarity, wavelength and velocity of propagating modes.It was shown that for low shear wave speed ratios an admissible mode might propagate at C gr albeit it is always unstable.For larger contrasts, admissible modes may steadily propagate at a speed slightly higher than the slowest S-wave speed and stable solutions exist for a limited range of small friction coefficients (Ranjith & Rice 2001).Moreover, a finite-time response to normal traction perturbations tends to promote stabilization of larger wavenumber modes (Ranjith & Rice 2001;Adda-Bedia & Ben Amar 2003;Brener et al. 2015).
As a result, numerical simulations of bimaterial rupture propagation were shown to be mesh dependent (Cochard & Rice 2000), even in the case of slip-weakening friction (Harris & Day 1997), due to the lack of a physical length (or time) scale associated with the coupling between sliding and normal traction perturbations.Moreover unphysical numerical solutions, such as the generation of a split pulse (Andrews & Ben-Zion 1997), were observed in relation to the ill-posed nature and the instability of the physical problem.
This regularization smoothens the response of the shear strength to sudden sliding-induced variations of the normal traction, introducing a monotonic fading-memory dependence on the normal traction history.This can be viewed as a relaxation mechanism where the relaxation time varies inversely with the actual slip rate, normalized by a characteristic sliding distance.Since this regularization may limit the spontaneous development of slow nucleation such as a pulse growing from a pore pressure increase with time (Andrews & Ben-Zion 1997), a constant relaxation time mechanism was also added (Cochard & Rice 2000).
This regularization can be interpreted as low pass filter acting on the normal traction variations (Kammer et al. 2014), which should ideally damp the high frequencies responsible for the fast instability growth without attenuating the low frequency content of the normal traction perturbations.
In the case of a linear slip-weakening friction, the characteristic sliding distance associated with the relaxation is expected to be smaller than the sliding distance D c related to the frictional weakening.A constant timescale relaxation was shown to reduce the emergence of numerical noise during rupture propagation more than the slip-rate dependent relaxation (Rubin & Ampuero 2007).More recently Kammer et al. (2014) identified a critical relaxation sliding distance below which the numerical solutions become independent Downloaded from https://academic.oup.com/gji/article/209/1/48/2918650 by CNRS -ISTO user on 10 August 2022 of the regularization parameters, using a mixed regularization involving both relaxation mechanisms.
Numerical simulations confirmed the possible occurrence of sliding pulse mode in the favoured direction (Cochard & Rice 2000;Ben-Zion & Huang 2002) as well as the asymmetry in bilateral rupture-in terms of rupture acceleration, slip rate in the vicinity of the crack front, and cumulative slip in the favoured direction-in particular in the case of slip weakening (Harris & Day 1997;Shi & Ben-Zion 2006;Rubin & Ampuero 2007) and velocity weakening frictions (Ampuero & Ben-Zion 2008).In the former case, a stationary slip pulse propagating at C gr (when this asymptotic speed exists) may separate from a growing rupture, due to the inversion of the slip gradient associated with normal traction perturbations (Rubin & Ampuero 2007).In the latter case, the emergence of wrinkle-like pulses and supershear transitions were retrieved, although regularization may affect the time and mode of their occurrence (Ampuero & Ben-Zion 2008;Langer et al. 2012).Finally the supershear transition along the not favoured direction, originally found by Shi & Ben-Zion (2006), was ascribed to shear waves ahead of the rupture front (Langer et al. 2012).These waves decrease the frictional shear strength along the non-favoured direction allowing the rupture to accelerate faster than the lowest shear wave speed.
Numerical simulations also showed that the asymptotic propagation speed of expanding bimaterial rupture tends towards the generalized Rayleigh wave speed C gr , when it exists, according to the Weertman solution (Rubin & Ampuero 2007).When generalized Rayleigh wave does not exist, the asymptotic rupture speed was numerically found to be either the slower S-wave speed (Rubin & Ampuero 2007), or slightly higher than this value (Harris & Day 1997).
Slip pulses propagating at C gr , originally found in numerical solutions of bimaterial rupture under classical Coulomb friction law (Andrews & Ben-Zion 1997;Ben-Zion & Andrews 1998), were also observed after regularization of the physical problem (Cochard & Rice 2000) even though they exhibit self-sharpening and divergent features (Ben-Zion & Huang 2002).This in turns indicates that a Prakash-Clifton type of regularization does not suppress the degeneracy of the solutions nor allow the selection of a physical pulse size (Adda-Bedia & Ben Amar 2003).Numerical simulations of the regularized problem were shown to become unstable when rupture propagates long enough at the near asymptotic speed (Ben-Zion & Huang 2002).
In this study a parametric study of the Prakash-Clifton type of regularization is performed, for both a constant and a slip-rate dependent relaxation time, for the acceleration phase of the in-plane rupture propagation along bimaterial interfaces under slip weakening friction law.For both relaxation mechanisms, the regularized numerical solutions are analysed in term of their mesh-dependence and regularization-dependence properties.The near-asymptotic phase of the propagating rupture is also specifically analysed from a numerical point of view.A modified dynamic regularization is finally discussed, where the relaxation time scales as the size of the dynamic dissipation zone normalized by a constant rupture speed, leading to new insights into the coupling between the relaxation and the frictional dissipation.

S E I S M I C RU P T U R E S A L O N G B I M AT E R I A L I N T E R FA C E S
Seismic waves travelling in elastic media are governed by the following equations: where ρ is the material density, v the particle velocity, σ the Cauchy stress tensor, and c is the linearized fourth-order symmetric elastic tensor.To solve the problem (1) appropriate initial and boundary conditions need to be specified.Free surface boundary conditions (zero traction) at the surface are assumed, while absorbing boundary layers, such as PML (Festa & Vilotte 2005), are introduced for unbounded physical domains.Modelling dynamic rupture requires also specific contact and friction conditions that couple kinematic and dynamic fields on the fault interface.
The fault is modelled as a zero-thickness frictional interface across which kinematic discontinuity (slip and slip rate) may occur.The interface separates two media (referred to as medium 1 and medium 2) that may have similar or dissimilar elastic properties.The interface is assumed smooth enough to define a normal vector n, here assumed as the outward normal with respect to medium 1.The traction on the fault interface is defined as,T = σ • n, where negative tractions are compressive.The traction is continuous across the fault, while the displacement u and the velocity v may have a discontinuity, denoted as the slip δu = u 2 − u 1 and the slip rate δv = v 2 − v 1, across the sliding portion of the interface.
The traction, the slip and the slip rate are classically decomposed along the normal and tangential directions.If ξ is one of these quantities, ξ = ξ n n + ξ t where the superscripts n and t refer to the normal and tangential contributions, with Frictional sliding along the fault is governed by non-smooth contact and friction conditions.In the normal direction, the Signorini contact condition is written as (2) This condition ensures the impenetrability of the two bodies when being in contact along the interface under compressive regime (T n < 0, δu n = 0).When the two bodies are no longer in contact, opening occurs (δu n ≥ 0), and the interface behaves as a free surface (T = 0).In this case, a mode I rupture is allowed.
When in contact, a frictional resistance is assumed when sliding that is governed by where μ is the actual friction coefficient and μT n the frictional strength.First condition in eq. ( 3) indicates that the non-sliding to sliding transition occurs when the modulus of the tangential traction T t lies on the surface of the actual frictional cone defined by the second condition in eq.(3).During sliding, T t remains on the surface of the friction cone and the third condition in eq. ( 3) implies that the tangential traction is collinear to the tangential slip rate.In 2-D, the last condition is simply a scalar condition.While tangential traction T t remains inside the frictional cone, no sliding occurs, that is, δv = 0.
In this study, a slip-weakening friction (Ida 1972)  Rupture dynamics along bimaterial interfaces 51 Bimaterial interface separates elastic domains with dissimilar properties.The stiffer material is the domain having the largest impedance, defined as the product of the material density by the S-wave velocity.Without loss of generality the stiffer material is assumed here to be medium 1, and the more compliant material is medium 2.
It has been established that sliding along a bimaterial interface, under classical or slip-weakening friction, is mathematically illposed and unstable against perturbations at all wavelengths and irrespective to the value of the friction coefficient, when the bimaterial contrast is such that a generalized Rayleigh wave exists.When a generalized Rayleigh wave does not exist, the problem is well-posed and stable only for a very limited range of low friction values with an upper bound that can be analytically characterized.The problem can be regularized introducing a delay between the slip-induced normal traction perturbation and the frictional shear response (Cochard & Rice 2000;Ranjith & Rice 2001).
Once regularized, the physical problem is no longer exactly the same, and becomes mathematically and numerically well posed, that is, numerical solutions do not exhibit mesh dependency.The delay of the frictional strength response to slip-induced normal traction perturbations can be accounted for by introducing an effective normal traction T n eff that is related to the normal traction by an exponential relaxation law This equation means that the relaxation is effective only where the fault is actually sliding.The effective normal traction replaces the normal traction in eq. ( 3) for the definition of the frictional strength: In the eq.( 5) t * is a characteristic relaxation time that may vary along the interface and depend on the rupture dynamics.Cochard & Rice (2000) suggested a slip-rate dependent relaxation time function t * = t * (|δv t |) defined as where δl is a characteristic sliding distance and t c is a constant timescale and the equation is valid only along the sliding portion of the fault.We also refer to δl as the relaxation slip parameter that can be scaled by the frictional slip-weakening distance D c .The relaxation time function ( 7) is the sum of two contributions.The first term in the right hand-side is the inverse of a slip-rate dependent relaxation time that can be interpreted as characterizing a loss of frictional memory of prior strength (Adda-Bedia & Ben Amar 2003; Rubin & Ampuero 2007).The second term is the inverse of a constant relaxation time that was originally introduced to avoid a singular behaviour at low slip rates, in particular during slow nucleation phase such as the one analysed by Andrews & Ben-Zion (1997), when rupture is originated by an external normal load.Nevertheless in many numerical simulations (e.g.Rubin & Ampuero 2007;Ampuero & Ben-Zion 2008;Langer et al. 2012) only the constant relaxation timescale is actually used to study bimaterial ruptures.
In this work, we numerically analyse these two mechanisms separately.We define t d = δl/|δv t | as the dynamic relaxation timescale, and t c as the constant relaxation timescale.A parametric study is performed to define numerically well-posed, as the solutions show-ing convergence for mesh refinement (Cochard & Rice 2000) and physically well-posed solutions, defined as numerical solutions that are independent of a particular choice of the regularization parameters.It is worth noting that the normal traction perturbations on the interface are due to sliding-induced perturbations and elastodynamic flux perturbations associated with bulk wave propagation.The latter variations are not expected to make the problem ill posed ahead of the rupture front, and wave transmission across non-sliding portion of the interface has to remain unchanged.Therefore the regularization has to be considered as a regularization of the sliding interfacial effects only, that is, of the sliding-induced normal traction perturbations, and it is defined only along the sliding portion of the interface, that is for |δv t | not equal to zero.
Finally, we limit the investigation to 2-D in-plane ruptures, where the tangential interface quantities are scalar.For sake of simplicity we indicate the normal traction as σ n , the effective normal traction as σ eff , the tangential traction (also referred to as shear traction) as τ , the tangential slip and slip rate as δu and δv respectively.The numerical simulations presented in this work are performed with the extension of the spectral element method (Komatitsch & Vilotte 1998) taking into account non-smooth contact and friction conditions (Festa & Vilotte 2006) whose main features are described in Appendix A.

S I M U L AT I O N S E T U P
Rubin & Ampuero (2007) widely analysed the problem of a bimaterial growing interfacial rupture under linear slip weakening friction law, and a similar setup was chosen in this study.
The 2-D model is described in Fig. 1: a straight interface, for example, a fault, separates two blocks of dissimilar properties.The densities ρ i and the S and P-wave velocities C s i and C p i , i = 1, 2, are assigned for each block.In this configuration the expected favoured direction is towards the right, being the direction of the slip in the more compliant medium.The dynamics of the rupture is driven by the four dimensionless parameters C s 1 /C s 2 , ρ 1 /ρ 2 , ν 1 and ν 2 , with ν i the Poisson's coefficients.For all simulations ν 1 = ν 2 = 0.25 is assumed.In agreement with the analytical results of Weertman (1980), dynamic features (asymmetry of slip rate, traction evolution, etc.) are shown numerically to be mainly sensitive to the ratio γ = C s 1 /C s 2 , and less influenced by the density ratio (Ben-Zion & Andrews 1998;Rubin & Ampuero 2007).Weertman (1980) showed that a steady state slip pulse can propagate along the favoured direction of a bimaterial interface inducing shear and normal perturbations according to where the moduli ζ and ζ * having the dimension of a traction, depend on both elastic properties and rupture velocity and x is the direction of propagation of the pulse.In particular ζ decreases with the increasing rupture speed and for small contrasts of impedance a real rupture speed exists for which τ (x) = 0.This speed is the generalized Rayleigh speed (C gr ) and it is the expected asymptotic speed for the growing rupture.The explicit expression for ζ can be found in the Appendix A (eq. A2) of Rubin & Ampuero (2007).Keeping the density uniform across the two layers and assuming ν 1 = ν 2 = 0.25, C gr is shown to be real for γ < 1.359 (Harris & Day 1997).In any case, the eq.(A2) of Rubin & Ampuero (2007) allows to compute the generalized Rayleigh speed for different contrasts of density.The simulation setup for the numerical models: the fault line is the thick black line in the middle whereas its prolongations on both sides form the fictitious boundary f (see Appendix A).We also represent the position of receivers along favoured direction (to the right of nucleation) where kinematic and dynamic quantities are collected as a function of time.Within the insets the expected linear slip weakening for the strength in homogeneous case (dotted lines) is compared to the typical shapes of the strength weakening for a bimaterial simulation.
Sliding along an interface, under linear slip weakening friction and a uniform remote loading, that separates two semi-infinite dissimilar homogeneous elastic media, leads to new complexities as compared to an interface separating similar media.In the latter case, due to symmetry, there is no normal-shear traction coupling and the rupture occurs under constant normal traction, for example, the pre-stress value σ 0 n .The shear strength exhibits therefore at all points of the interface the same linear weakening as the friction.As a consequence, the fracture energy G c , defined as the area under the strength-weakening curve in a slip-traction plot, can be simply evaluated as and it is a property of the interface.In the case of a bimaterial interface, the effective normal traction σ eff dynamically changes as a response to the sliding-induced normal traction perturbations.During the acceleration phase, at receivers along the favoured direction, the increasing compressive perturbation of the normal traction ahead of the crack tip increases the strength, whereas the extensive perturbations induced by the crack front arrival make the weakening sharper and the dynamic level of shear strength will be lower than in the case of a homogeneous medium.The fracture energy G c dynamically changes according to the perturbations of the normal traction and to the relaxation parameters of the formulation, as shown in the insets of Fig. 1.
In most of the simulations and unless otherwise stated, we fixed a uniform density at a representative crustal value (ρ 1 = ρ 2 = 2700 kg m −3 ), while we investigated several values of γ , corresponding to the existence or not of a generalized Rayleigh wave speed C gr .
On the fault, the initial normal traction is set to the uniform value σ 0 n = −100 MPa, whereas the initial shear stress is τ 0 = 62.5 MPa.The linear slip weakening parameters are μ s = 0.7, μ d = 0.6 and D c = 6 mm.Accordingly the strength parameter The rupture is initiated increasing the initial tangential traction slightly above the initial static frictional strength (0.5 per cent) over a fixed length scale.The size of this patch is selected to be slightly larger than the nucleation length L c for a linear slip-weakening rupture (Uenishi & Rice 2003;Rubin & Ampuero 2007): where the effective elastic modulus for quasi-static plane strain deformation ζ depends on the shear moduli ρ i C 2 s i and the Poisson coefficients ν i of the two media and it is defined as ζ of eq. ( 8) when the rupture speed is equal to zero (Rubin & Ampuero 2007).
W is the initial slope of the slip weakening frictional strength at the nucleation: where sliding-induced perturbations on the normal tractions are neglected.With this condition, the rupture dynamically grows, mimicking the propagation of a seismic rupture without considering the initial quasi-static nucleation phase.
For numerical simulations, a regular spectral element method mesh is considered, with square elements with interpolation polynomials of degree 8 defined by 9 × 9 GLL collocation points.The maximum element size h = 12 m for all the simulations guarantees at least five points per wavelength during the rupture propagation.The Courant number for all simulations is always smaller than 0.35, ensuring stability of the second-order explicit time stepping (Komatitsch & Vilotte 1998).The fault is embedded in an infinite medium numerically modelled using the Perfectly Matching Layers (Festa & Nielsen 2003;Festa & Vilotte 2005) along the edges of the numerical domain.

Regularization with slip-rate dependent relaxation time
We refer to a dynamic relaxation timescale when t * = t d = δl/|δv| in eq. ( 7), where the relaxation sliding distance δl is the only parameter that controls the regularization.It is chosen in the range (2%D c −100%D c ) to ensure that the regularization of normal perturbations induced by the propagating rupture mainly occurs in the vicinity of the crack front.
The regularization depends on the tangential slip rate δv, and therefore the relaxation time is local and point-wise varying along the interface.Variations of the relaxation time are relevant around the crack tip and within the dissipation zone, that is, where the energy is dissipated according to the shear strength weakening and where the slip rate sharply increases to its maximum.Beyond this maximum, the slip rate decreases towards an almost constant value outside the cohesive zone almost leading to a constant time regularization with a larger relaxation time.
Convergence of the numerical solutions is first considered.Once regularized, the physical problem is no longer exactly the same as it originally was.Therefore, convergence needs to be studied as a function of the mesh size h and of the relaxation sliding distance Rupture dynamics along bimaterial interfaces 53 δl respectively.Convergence for mesh refinement is a condition for numerically well-posed solutions in the sense of Cochard & Rice (2000).This is analysed by numerically investigating the maximum value of the mesh size h, at a fixed δl, below which the solutions do not depend anymore on the discretization length within the dispersion error.Then the numerically stable solutions, that is, mesh size independent solutions, are investigated as a function of the relaxation sliding distance δl.An upper limit for the sliding distance parameter δl max is expected below which numerical solutions do not depend on δl, within the dispersion error (Kammer et al. 2014).This convergence is here referred to as physical convergence of the solutions.Both convergences are measured using the L 2 -norm of the difference between the numerical solution and a reference one obtained using the finest mesh and the smallest regularization parameter.Since the regularization introduces a time delay in the kinematic and dynamic fields, the traces are first aligned by crosscorrelation before computing the L 2 -norm difference.The relative error is finally obtained by normalizing this result by the L 2 -norm of the reference solution.
The comparison between the different models is presented both in space and time domains.In the space domain the comparison is based on the slip rate that allows to identify the position of crack tip, to characterize the rupture speed, and to analyse the asymmetry between the two directions of propagation of the crack.In the time domain, the comparison is based on the effective normal traction σ eff recorded at receivers located along the interface, at increasing distances from the initiation zone, as indicated in Fig. 1.Both quantities are representative of the rupture dynamics during the whole acceleration phase.
First a model for which the generalized Rayleigh wave speed exists (γ = 1.18,C s 2 = 2.620 km s −1 and C s 1 = 3.092 km s −1 ) is investigated.
Fig. 2(a) shows the slip-rate profile at time step t 0 = 0.12 s for δl = 2%D c which is the smallest value used in this study.Fig. 2(b) is a zoom of Fig. 2(a) around the rupture front.We analysed the convergence for different mesh sizes, h = 2, 3, 4, 6, 12 m.
For coarser meshes (h > 4 m) the rupture is faster in both directions as compared to finer grids (Figs 2a and d).Additionally for δl = 2%D c and h > 4 m strong oscillations occur, also producing pathological effects in the not favoured direction, such as multiple pulses due to continuous arresting and restarting of the rupture (Fig. 2a).These results hold for all the acceleration phase of the rupture.
Figs 2(c) and (d) are the same representation of Figs 2(a) and (b) at the same time step, for δl = 10%D c .In this case, slip-rate oscillations for the coarser meshes are considerably damped (Fig. 2c).Nevertheless when zooming around the rupture front (Fig. 2d), the rupture in the coarser meshes is still in advance as compared to the  slip-rate evolution observed in finer grids.For the two values of δl showed in Fig. 2, mesh convergence is achieved when h ≤ 4 m.
The same convergence condition is obtained when analysing σ eff as a function of time.In Fig. 3, the effective normal traction is represented as a function of time at receiver 5.For δl = 2%D c (Fig. 3a), the coarser meshes (h = 6 m and h = 12 m) clearly show oscillations, whose features are similar to those retrieved by Kammer et al. (2014) for the slip rate, while the maximum value of the effective traction occurs earlier in time (Fig. 3b).For finer meshes (h ≤ 4 m) the curves overlap (Figs 3a and b) as a result of the mesh convergence.
These features are preserved also for larger δl, for which the oscillations of σ eff are more and more damped (Figs 3c and d).
The maximum value of the element size, below which stable solutions are observed, is finally found to be almost independent of the specific value of δl, in the explored range of variation of this parameter.In this analysis numerically well-posed solutions are found for h ≤ 4 m.This is slightly different from the results of Kammer et al. (2014), who consider an arresting slip pulse, for which smaller δl values require finer meshes.
Having characterized for different δl the maximum element size, below which solutions are mesh independent, the dependence of the solution on the sliding distance parameter δl of the dynamic regularization is investigated.
Supporting Information Fig. S1 shows the slip-rate profile at time step t 0 = 0.12 s, while the inset shows a zoom of the same profile around the rupture front in the favoured direction.The figure highlights the expected asymmetry of the bimaterial rupture propagation with larger differences around the rupture front in the favoured direction.Here, convergence of the maximum amplitude of slip rate is achieved for δl ≤ 20%D c .Above this bound, the maximum amplitude decreases as δl increases.The position of the rupture front is more sensitive to the regularization and becomes independent of δl when δl ≤ 10%D c .Above this bound, the rupture speed decreases as δl increases.In the favoured direction, the normal traction perturbation changes sign moving from a compressive regime ahead of the rupture front to an extensive regime behind the rupture front.Increasing δl corresponds to increasing the dynamic relaxation time which is no more able to properly resolve the sharp variation of the normal traction at the rupture front.The regularization damps the high frequency energy of the propagating rupture within the dissipation zone, decreasing the maximum amplitude of the slip rate and preventing fast acceleration of the rupture.In the not favoured direction, an opposite behaviour is observed as a function of δl, although the effect is less pronounced as compared to the favoured direction.The normal traction perturbation is now extensional ahead of the rupture front and compressive behind the rupture front.
The so-called physical convergence of the solutions can be also shown in time domain looking at the variations of σ eff .Fig. 4(a) shows the evolution of σ eff at receiver 5.The curves superimpose before the arrival of the rupture, as expected since ahead of the rupture front the normal perturbation is not regularized along the Downloaded from https://academic.oup.com/gji/article/209/1/48/2918650 by CNRS -ISTO user on 10 August 2022 non-sliding interface.Behind the rupture front, the curves are different and depend on δl.These differences are enhanced in the zoom around the maximum value of the effective normal traction (Fig. 4b) where physical convergence of the solutions can be inferred from the superposition of the curves.Supporting Information Fig. S2( In the case of the dynamic timescale relaxation, numerical solutions are shown to become independent of the regularization pa-rameter during the whole acceleration phase, that is, of the sliding distance δl, for all δl below an upper bound that can expressed as a fraction of D c .This results from the fact that the relaxation mechanism is an adaptive low-pass filter of the normal traction, with a cut-off varying as a function of the slip rate (Kammer et al. 2014).The dynamic relaxation timescale hence provides adaptive cut-off frequencies along the rupture.As the rupture accelerates towards the asymptotic speed, the slip rate at the rupture front sharply increases, also increasing the cut-off frequency of the filter.To clarify this interpretation the amplitude spectrum of σ eff is shown in Fig. 5 at the same receiver analysed in Fig. 4(a).In Supporting Information Figs S3(c) and (d), the amplitude spectra of σ eff at receivers 2 and 8 are plotted.Different physical and numerical frequencies associated to the propagating rupture can be identified in the spectra, and are marked with dotted lines in Fig. 5.The lowest frequency is related to the largest physical timescale associated with the normal traction perturbations and it usually corresponds to a first slope break in the spectrum (black dotted line).At the end of the nucleation phase, this timescale is associated to the waves generated during the initiation phase.As the rupture accelerates, this timescale becomes shorter and shorter.Since the energy balance is mostly controlled by the size of the weakening region of the frictional strength, the first corner frequency in the spectra is associated with the duration of the weakening process at a given point of the interface.At shorter timescale within the dissipation zone, the frequency linked to the coupling between tangential and normal traction perturbations can be identified (red dotted line).This coupling timescale is estimated as the delay between the maximum of slip rate and the rupture front.The timescales associated with the weakening process and the normal/shear traction coupling are shown in Fig. 6.Since instability of bimaterial rupture comes from the high-frequency coupling between normal and shear tractions, the cut-off frequency associated with the dynamic regularization has to lie between the characteristic frequency of the normal traction variation and the coupling frequency to preserve the specific timescales of the propagating rupture, while damping the unstable frequencies generated by the coupling.In Fig. 5 this cut-off frequency, for δl ≤ 10%D c is marked by the magenta dotted line.Finally, the frequency associated with the numerical resolution can be identified as the maximum one well resolved by the discretization size, which for spectral element method depends on the smallest velocity and the minimum number of collocation points (∼5) per wavelength required to accurately resolve the dominant wavelength of propagating seismic waves (Komatitsch & Vilotte 1998).It is marked by a green dotted line.
During the acceleration phase, the cut-off frequencies associated with the weakening and the coupling processes increase while the dissipation zone shrinks.As a result the amplitude of the slip rate increases and the normal traction perturbation follows the same evolution as the slip rate (Fig. 6).Indeed we argue that a dynamic relaxation filter, which adapts the cut-off frequency to slip rate variations, is able to properly filter the normal traction close to the crack front all along the acceleration phase.To better understand what drives the physical convergence of the numerical solutions in the case of a dynamic regularization of the slip-weakening friction law, the evolution of the effective normal traction σ eff versus that of the actual normal traction σ n is analysed along the sliding portion of the interface.In Fig. 7 the difference σ n − σ eff is plotted as a function of slip, for all the receivers indicated in Fig. 1.For a fixed sliding distance parameter δl, the slip value δu * at which σ n − σ eff is zero, behind the crack front, does not seem to depend on the position of receiver and thus appears to be independent of the slip rate and the rupture speed during the acceleration phase, although the maximum of the difference between σ n and σ eff increases as the crack grows up.When δl ≤ 10%D c , it is found that δu * < D c , as shown in Fig. 7 The regularization provides physically convergent numerical solutions when the sliding distance over which it is effective, remains smaller than the sliding distance over which the dissipation takes place.For δl > 10%D c , the filter operates over a sliding distance scale larger than the dissipation and as a result it filters the physical scale that we need to be resolved during the rupture propagation.For a linear slip weakening friction law, these results are naturally scaled by D c , which is the characteristic finite sliding distance over which the weakening of the friction takes place.
The same convergence analysis was also performed for a higher impedance contrast (γ = 1.80) at which the generalized Rayleigh wave speed does not exist.This is done varying the shear wave speeds while keeping constant the density and the effective shear modulus ζ .This allows to start the rupture using the same initiation length as before.
For a larger γ the rupture is faster along the favoured direction and slower along the opposite direction enhancing the rupture asymmetry.The mesh dependence and the δl analysis show similar results.Supporting Information Fig. S4 shows the evolution of σ eff in time (Supporting Information Fig. S4a) and frequency (Supporting Information Fig. S4b).Physical convergence is again achieved when δl ≤ 10%D c .Since for γ = 1.80 the asymptotic propagation phase is approached more rapidly, the results are plotted for sake of clarity at a closer receiver (receiver 3) where the rupture is still far from the asymptotic phase.The characteristic slip δu * associated with the regularization appears now to depend on the receiver position during the acceleration phase (Supporting Information Fig. S4c).However it is worth noting that the normal traction perturbations are getting stronger than in the previous case (σ n − σ eff is about 8 times larger at receiver 5).Resolving the slip at which σ n − σ eff is zero becomes harder due to the strong normal traction perturbation at crack front.Nevertheless it is still observed that as long as δu * < D c the physical convergence is achieved both in time and frequency.This condition holds until the rupture approaches the asymptotic speed (see Section 5).

Regularization with constant relaxation time
When a constant relaxation timescale t c is used, the regularization of sliding-induced normal tractions perturbations do not depend anymore on the slip rate and is therefore independent of the rupture dynamics.
The parametric study is performed varying t c in the range 1.2 × 10 −4 ≤ t c ≤ 6.0 × 10 −3 .Defining t c = δl δv * , as proposed by Cochard & Rice (2000), corresponds to a range of variation of δl between 2%D c and 100%D c for δv = 1 ms −1 , allowing to directly compare the constant timescale regularization to the dynamic one, as analysed in the previous subsection.
As for the dynamic timescale regularization, two parametric analyses are performed to study the mesh dependence and the physical convergence.All the results presented here refer to the case γ = 1.18 with the same elastic properties as used in the Section 4.1.
The mesh dependence analysis provides very similar results as compared to the dynamic timescale analysis.The analysis can be summarized looking at the slip rate profiles at time step t 0 = 0.12 s in the space domain for the smallest value of t c (t c = 1.2 × 10 −4 s, Supporting Information Figs S5a and b) and for t c = 6.0 × 10 −4 s (Supporting Information Figs S5c and d).As for the dynamic timescale regularization, the coarser meshes (h = 6 m and h = 12 m) provide mesh-dependent solutions.In the first case, spurious oscillations are observed (zoom in Supporting Information Fig. S5b) while in the second case these oscillations are damped within the rupture (Supporting Information Fig. S5d).Although the simulations qualitatively yield the same numerical results as for the dynamic timescale regularization, the oscillations are more pronounced in the case of a constant timescale regularization.Solutions still become mesh independent for h ≤ 4 m during all the acceleration phase, and therefore numerically well-posed.
The influence of the constant relaxation timescale is investigated by varying t c .In contrast with the dynamic relaxation timescale, numerical solutions never become independent of the regularization parameter and therefore they correspond all to different physical problems.In this respect a constant relaxation timescale regularization never achieves a physical convergence.This can be observed both for the kinematic and dynamic fields, in space and time domains.
In Fig. 8 σ eff is plotted as a function of time at two receivers.Although solutions appear to be physically convergent below some value of t c at the beginning of the acceleration phase (receiver 2, Fig. 8a), this no more the case farther away from the initiation phase (receiver 8, Fig. 8c).
This lack of physical convergence can be explained in comparison with the results obtained for the dynamic regularization.An equivalent δl max eq can be estimated, all along the rupture propagation, in the case of a constant relaxation time as δl max eq = δv max • t c , where δv max is the maximum slip rate observed at each point of the fault, this value being reached soon after the rupture front passes through the point.For the receivers of Fig. 8, considering t c = 6 × 10 −4 s the equivalent sliding distance at receiver 2 is δl max eq#2 ∼ 6 × 10 −4 m = 10%D c while for receiver 8 it is δl max eq#2 ∼ 2.4 × 10 −3 m = 40%D c .The equivalent sliding distance increases as the rupture moves away from the initiation zone, owing to the sharpening of the slip rate.For this specific t c , δl max eq#8 is well outside the range of sliding distances that were shown to guarantee a physical convergence.
Analysing the solutions in the frequency domain (Figs 8b and  d), a fixed relaxation timescale implies a fixed cut-off frequency for the regularization all along the acceleration phase (yellow dotted lines in Figs 8b and d).Conversely the physical and the coupling frequencies increase with the rupture growth (respectively black and red dotted lines in Figs 8b and d).Since the size of the dissipation zone goes to zero as the rupture tends towards the asymptotic speed, there will be always a position on the fault beyond which the timescale associated to the shear strength weakening process will become shorter than t c .In this case, the regularization will overfilter the energy balance during the weakening process, providing no longer physical convergent solutions.
An effect deriving from the lack of physical convergence is the slower acceleration of the rupture towards the asymptotic speed as compared to dynamically regularized solutions.In Fig. 9  physical convergence no longer holds, the rupture speeds differ more and more due to the excess of filtering in the case of a constant relaxation time.
In light of this parametric study, the classical Prakash-Clifton regularization (eq.( 5)) where t * is given by eq. ( 7) can be rewritten as where f d = |δv|/δl is a dynamic frequency and f c is a constant cutoff frequency.The former provides a self-adaptive regularization, which leads to mesh independent and physically convergent solutions when the element sizes and the sliding distance parameters are within a finite range.The latter still provides numerically wellposed solutions, that is, mesh independent solutions, but does not lead to physically convergent solutions.When f c is small enough and the absolute value of the slip rate is large enough to guarantee f d f c , the numerical solutions are similar to those provided by the dynamic timescale regularization.This occurs in most of the simulations using the complete Prakash-Clifton regularization.Nevertheless, when f c ≥ f d the numerical solutions, although numerically well-posed, depend now on the parameters chosen for the regularization and therefore they no longer correspond to the same physical problem.As proposed by Cochard & Rice (2000) adding a constant timescale relaxation help in modelling very slow nucle-ation processes, even though the choice of f c becomes nucleation dependent.

RU P T U R E B E H AV I O U R T O WA R D S A S Y M P T O T I C S P E E D
The numerical analysis of bimaterial rupture propagation under linear slip-weakening friction law was so far focused on the acceleration phase.Results of Section 4 show that dynamic regularization, with a slip-rate dependent relaxation time, leads to meshindependent solutions that are also independent of the regularization parameter under certain conditions.
In this section, the asymptotic regime of the rupture propagation is numerically investigated using the previous results.
For in-plane rupture propagation along a planar interface, under linear-weakening friction, separating isotropic linear elastic bodies made of identical materials, the asymptotic rupture speed is either the Rayleigh wave speed in the sub-shear regime or the P-wave speed in the supershear regime (Burridge 1973;Andrews 1976).
For an interface separating different materials, the predicted asymptotic rupture speed depends on the velocity contrast between the two materials.When the generalized Rayleigh wave exists along the interface, that is, S-wave velocity ratio C s 1 /C s 2 < 1.359 when the materials have the same density, the predicted asymptotic rupture speed is the generalized Rayleigh wave speed C gr .When Downloaded from https://academic.oup.com/gji/article/209/1/48/2918650 by CNRS -ISTO user on 10 August 2022 for dynamic and constant timescales: while a convergence among different timescales is still detectable, the acceleration is equivalent to that deriving from dynamic timescale models.Conversely when the timescale is too large to properly regularize the problem the rupture speeds differ more and more and the acceleration is slower for the constant timescale.this latter wave does not exist, Ranjith & Rice (2001) found that an unstable steady-state mode might propagate at a speed slightly higher than the shear wave speed C s 2 of the more compliant material.Numerical solutions produced contrasting results either in agreement with these results (Harris & Day 1997) or suggesting an asymptotic rupture speed close to C s 2 (Rubin & Ampuero 2007).

Interface with no existing generalized Rayleigh wave
Numerical solutions of the dynamically regularized physical problem are investigated here for large γ values, with similar (ρ 1 = ρ 2 ) or dissimilar (ρ 1 = ρ 2 ) densities.For all the considered γ parameters, generalized Rayleigh wave does not exist, that is, eq.(A2) in Rubin & Ampuero (2007) does not have real roots.
The mesh discretization and the regularization parameters are chosen to be such that mesh independency and physical convergence are satisfied following results of Section 4.1.
Seven different contrasts of impedance γ , between 1.5 and 2.1, have been explored first varying the shear wave speed of the more compliant material for ρ 1 = ρ 2 = 2700 kg m −3 and C s 1 = 4058 ms −1 .Finally, the P-wave velocity of the two materials was selected such that ν 1 = ν 2 = 0.25.
In Fig. 10(a) the numerical rupture speed, normalized to C s 2 , along the favoured direction of propagation is shown for different values of γ as a function of the propagation distance.Asymptotic rupture speed is shown to be a function of γ and in all cases it is shown to tend towards values larger than C s 2 .As such C s 2 is not the asymptotic speed for these simulations, and the rupture can accelerate to larger values as suggested by Ranjith & Rice (2001).The asymptotic speed increases as γ increases, as shown in Fig. 10(b) where error bars represent the standard deviation for the average rupture speed determination.
The capability of the rupture to accelerate towards a speed larger than the shear wave velocity of the more compliant medium generates peculiar effects.As the rupture proceeds at a sub-shear speed, the emitted radiation, although asymmetric, exhibits the classical pattern with P and S waves ahead of the rupture.In Fig. 11(a) the kinetic energy field around the fault is shown at a fixed time step when the rupture speed is lower than C s 2 , and most of the energy is confined in the more compliant material.When the rupture accelerates towards a speed larger than C s 2 , classical Mach cone can be recognized in the half space occupied by the more compliant material (Fig. 11b).This acceleration also generates a change in the normal traction perturbation pattern along the favoured direction.Compressive normal traction perturbation ahead of the rupture front and extensional normal traction perturbation at the rupture front still hold, but S-waves emitted behind the rupture generate now a small extensional normal traction perturbation behind the dissipation zone, see Fig. 11(c), where the perturbation is evidenced by a black circle.Henceforth slip rate and normal traction perturbations at the crack front remain almost constant during all this phase and the dissipation zone can be resolved with enough grid points, that is, at least three points.
For each of the investigated values of γ , variations of the density contrast have been explored such that C gr does not exist.Some of the results are summarized in Supporting Information Figs S6(a) and (b).Supporting Information Fig. S6(a) shows the evolution of the rupture speed for three different ratios γ = C s 1 /C s 2 (γ = 1.7, 1.9, 2.1) with a fixed density contrast of 1.2 (ρ 1 = 3240 kg m −3 , ρ 2 = 2700 kg m −3 ).As discussed before, the asymptotic speed increases as γ increases and is  always larger than the shear wave speed C s 2 .Supporting Information Fig. S6(b) shows now the evolution of the rupture speed for four different contrasts of density, ρ 1 /ρ 2 = 1.0, 1.2, 1.4, 1.6, with ρ 2 = 2700 kg m −3 in all cases, at a fixed value of γ = 1.9.The asymptotic speed is a function of the density contrast, that is, it increases as the density contrast increases, and in all cases the rupture accelerates again towards an asymptotic speed larger than the S-wave speed of the more compliant material.
It is worth to note here that dynamic timescale regularization of the physical problem lead to mesh-independent numerical solutions that are physically convergent during the whole rupture propagation including the asymptotic phase.The asymptotic rupture speed is found to be always larger than the shear wave speed C s 2 , independently of the regularization parameter in the parameter domain of physical convergence estimated in Section 4.1.

Interface with an existing generalized Rayleigh wave
The asymptotic phase of in-plane bimaterial rupture is now investigated for values of γ small enough to ensure the existence of generalized Rayleigh wave with wave speed C gr .In particular the acceleration phase for γ = 1.18 and C gr = 2570 ms −1 was analysed in Section 4.1.The rupture speed is shown to tend towards the asymptotic C gr wave speed independently of the contrasts of velocity and density when the physical problem is dynamically regularized, using slip rate dependent relaxation time.
However when the rupture speed approaches the generalized Rayleigh wave speed, the rupture propagation becomes numerically unstable and numerical solutions are rapidly polluted by spurious oscillations.Despite the dynamical regularization, this instability occurs for all the regularization parameters that lead to physical convergent solutions.
The instability is generated by the continuous shrinking of the dissipation zone, driven by the increasing singular behaviour of the normal traction perturbations at the rupture front.
As the rupture speed tends towards C gr , the absolute amplitude of the compressive normal traction perturbation just ahead of the rupture front increases, and as such the shear strength of the interface is increasing.Immediately behind the rupture front, the normal traction, within the dissipation zone, perturbation changes sign becoming more and more extensional.This variation of the normal traction perturbation occurs over a shrinking dissipation zone leading to higher and higher frequency variations as the rupture propagates close to the C gr wave speed.As a result the energy release rate increases with a continuous transfer of energy to shorter wavelengths.This leads to a singular behaviour of the slip rate, and the problem cannot be regularized using slip-rate dependent relaxation time.
As the dissipation zone shrinks, the numerical mesh can no more properly resolve it.The numerical solution is no longer able to propagate the high frequencies generated at the rupture front and increasing spurious oscillations rapidly pollute the whole simulation.This is shown in Fig. 12, where the number of grid points actually resolving the dissipation zone is shown as a function of position of the rupture front.
During the initial phase of the instability, a slip pulse emerges (red curve in Fig. 12), which might be triggered by the occurrence Downloaded from https://academic.oup.com/gji/article/209/1/48/2918650 by CNRS -ISTO user on 10 August 2022 of a slip gradient minimum, behind the dissipation zone, associated with the change of sign of the normal traction perturbation (Rubin & Ampuero 2007) but more probably by the spurious oscillations that can unload the interface below the frictional strength level.
Results of Section 4 show that dynamical regularization, with a slip-rate dependent relaxation time, leads to physical convergent solutions as long as it adaptively filters slip-induced normal traction perturbations without modifying the energy balance and the energy scale transfer driving the rupture propagation.The latter is associated with the size of the dissipation zone.However, when the adaptive cut-off frequency of the filter moves beyond the mesh cut-off frequency, aliasing effects occur together with spurious numerical oscillations (Festa & Vilotte 2005).
As expected, a larger sliding distance parameter δl allows the rupture to propagate over a longer distance as a result of increased smoothening of the slip-induced normal traction perturbations and stronger attenuation of the spurious oscillations.Nevertheless as the rupture speed gets close to C gr , emergence of spurious oscillations cannot be avoided leading to unstable solutions independently of the value of the regularization parameter δl.
It is worth to note here that shrinking of the dissipation zone is also observed in the case of an interface separating identical materials or of a bimaterial interface with no generalized Rayleigh wave, when the rupture approaches the predicted asymptotic speed.However in both cases numerical solutions lead to an asymptotic rupture speed that tends from below to the predicted one, as the mesh is refined, and always end up with a dissipation zone that remains resolved by at least three grid points.In contrast when generalized Rayleigh wave exists, numerical solutions behave quite differently as a result of the singular behaviour of the normal traction perturbation, which changes sign at the rupture front.
In the case of isotropic linear elastic similar bodies sliding along interfaces separating them there is no coupling between interfacial slip and normal traction variations, and the slip gradient decreases monotonically away behind the rupture tip when initial conditions along the interface are uniform.In contrast in the case of a bimaterial interfaces due to broken symmetry such coupling exists but there is a qualitative difference depending on the existence of generalized wave along the interface.In the case where generalized Rayleigh wave does not exist, the rupture speed becomes larger than the S-wave speed of the more compliant medium and therefore elastodynamic coupling between interfacial slip and normal traction perturbations mainly affect the region behind the rupture front, as previously discussed, and this coupling does not lead to a singular behaviour at the rupture front nor produces a local minimum of the slip gradient behind the dissipation zone.
In an attempt to mitigate the singular behaviour observed for bimaterial interface when generalized Rayleigh wave exists, slip-rate dependent time regularization might be switched to a constant time regularization as the rupture approaches the asymptotic speed.This switch can be numerically triggered based on the actual resolution of the dissipation zone, that is, when the number of grid points in the dissipation zone becomes less than 5 corresponding to half of the element size for these simulations, or by introducing an upper cut-off limit for the slip rate dependence.
Different values of the constant relaxation time t c have been explored spanning two orders of magnitude.While instability can be delayed when increasing t c , that is, the larger is t c the longer the rupture can propagate, it can never be suppressed.In those simulations a metastable slip pulse is also generated, which propagates at almost the C gr wave speed during a certain time before the simulation finally blows up.
The behaviour of the rupture when approaching the asymptotic speed has been further investigated regularizing the physical problem with a constant relaxation time.Results of Section 4 show that constant relaxation time leads to a fixed cut-off frequency and that there is no physical convergence of the numerical solutions during the acceleration phase, that is, solutions always depend on the regularization parameter.Nevertheless, rupture speed still tends towards the asymptotic generalized Rayleigh wave speed.When approaching the asymptotic speed, a propagating finite-size pulse emerges and it propagates during a certain time along the interface before simulation finally blows up.The onset of the pulse is meshindependent but depends on the constant relaxation time t c and it occurs earlier for smaller t c (Fig. 13a).
Once the pulse is generated, the size of the pulse remains constant during propagation.The size selection depends on the filter cut-off frequency introduced by the constant relaxation time regularization.A smaller t c leads to a smaller pulse size as shown in Fig. 13(b).In particular defining the non-dimensional parameter L p as the pulse size normalized by the element size h, we found L p ≈ 0.75 for t c = 3.0 × 10 −4 s, L p = 1.50 for t c = 6.0 × 10 −4 s, and L p = 2.50 for t c = 1.2 × 10 −3 s.It is worth to recall that we have 9 grid points for each element: this ensures that the pulse size is always well resolved.The pulse size appears also to be mesh dependent and slightly decreases, as the mesh is refined.
The pulse size selection seems to be associated to competing effects between constant relaxation time regularization and mesh discretization as the size of the dissipation zone continuously decreases.The same behaviour occurs for a dynamical scale regularization, but in this case as adaptive cut-off frequency of the regularization filter moves beyond the mesh frequency cut-off, spurious oscillations blur the numerical solutions.
The instability of the pulse, which emerges for all the adopted relaxations, owes to the singular behaviour of the slip-induced normal traction perturbations and of the slip rate that a Prakash-Clifton law type cannot regularize.
It has long been recognized that sliding along a bimaterial interface under classical Coulomb friction law is unstable against perturbations at all wavelengths and irrespective of the value of the friction coefficient, when the bimaterial contrast is such that Pulse size for different models obtained for constant timescale regularizations.The sizes, at different time steps, are compared for different t c .In all cases the pulse size is pretty constant during its propagation and when a larger t c is used the average size of the pulse is quite larger.
When the physical problem is well-posed with a Prakash-Clifton type of law, numerical solutions (Ben-Zion & Huang 2002) still exhibit self-sharpening and divergent behaviour of wrinkle-like pulse propagating along bimaterial interface after long enough propagation distance, in agreement with analytical stability analysis (Ranjith & Rice 2001;Adda-Bedia & Ben Amar 2003).According to the results of this study, this feature also emerges when slip pulse detaches at the end of the acceleration phase of a growing rupture propagating along a bimaterial under linear slip-weakening law.

M O D I F I E D R E G U L A R I Z AT I O N W I T H D I S S I PAT I O N L E N G T H D E P E N D E N T R E L A X AT I O N T I M E
The dynamic timescale was shown to provide physical convergent solutions as long as the relaxation of the normal traction perturbations occurs over a sliding distance slip length smaller than the frictional slip weakening distance D c .This aspect suggests that the instability of the solutions for a bimaterial interface under slipweakening friction law arises from sliding-induced normal traction perturbations at the scale of the frictional weakening zone, where most of the dissipation and wave emission take place.
For this reason an alternative regularization is proposed in order to link the relaxation timescale to the size of the dissipation length.We can describe this new regularization by the following equation: where now L d is the size of dissipation zone, V a is a reference rupture speed, which is selected to be the expected asymptotic rupture velocity, and β is a non-dimensional parameter used to perform a parametric analysis, by analogy with the dynamic timescale.The eq. ( 12) still provides a dynamic, adaptive relaxation timescale, since the size of the dissipation zone is shrinking during the acceleration of the rupture leading to smaller and smaller relaxation times.In contrast to the slip-rate dependent relaxation timescale, this relaxation timescale is non-local and is related to a characteristic length scale of the rupture.Even though such a formulation is not numerically convenient since it requires the determination of the size of the dissipation zone at each time step, it is worth to explore it here for a better understanding of the connection between the normal traction regularization and the frictional dissipation.
This regularization provides mesh independent numerically wellposed solutions as mesh is refined (h ≤ 4m, where h is the element size); conversely the coarsest meshes show spurious oscillations.This is shown here looking at the slip rate profile at time step t 0 = 0.12s, both for the smallest used value of β(β = 0.01, see Supporting Information Figs S7a and b) and for β = 0.10 (Supporting Information Figs S7c and d).
Then, we investigated the physical convergence, that is, numerical solutions independent of the regularization parameter.The numerical solutions were again compared in time and frequency domains to check the influence of the parameter β parameterization on the rupture dynamics during the acceleration phase.Similarly to the dynamic timescale, the convergence is achieved for small β values, for which σ eff becomes independent of the parameter β at all receivers.In Fig. 14(a) the time evolution of the effective normal traction is plotted for different β values at receiver 5.The inset in this figure shows that the convergence is achieved for β < 0.10.As for sliprate dependent relaxation timescale, within β < 0.10 the cut-off frequency of the filter associated with the regularization is located between the characteristic scale of the dissipation process and the coupling scale, during all the rupture acceleration.Fig. 14(b) shows the amplitude spectra for σ eff at receiver 5.For β < 0.10 the spectra superimpose in the low frequency range.Figs 14(c) and (d) show that this regularization mechanism also leads to a characteristic slip δu * , with δu * < D c for β ≤ 0.10 (Fig. 14c), while δu * > D c for β > 0.10 (Fig. 14d).The similarity between the results obtained when regularizing the solutions with the dissipation length and the slip rate, suggests performing a direct comparison between the two kinds of regularization.
The slip-rate dependent relaxation time formulation is modified such that the relaxation time is now dependent on the maximum slip rate along the interface at a given time δv max .This value always occurs in the vicinity of the rupture front and it defines a characteristic velocity scale: The numerical solutions for the two regularizations superimpose as it can be observed in the slip rate profile at time step t 0 = 0.13s (Fig. 15).Since the relaxation slip can be also expressed as δl = β D c Downloaded from https://academic.oup.com/gji/article/209/1/48/2918650 by CNRS -ISTO user on 10 August 2022 the equivalence between the two regularizations also implies that the dissipation length scales as the inverse of the maximum slip rate: This scaling law characterizes the rupture dynamics and it suggests that the dynamics of a growing in-plane rupture along a planar bimaterial interface under slip-weakening friction is mainly driven by the bimaterial coupling between interfacial slip and normal traction perturbations at the scale of the actual size of the dissipation zone.Thus provides further evidence that regularization based on slip-rate dependent relaxation time only leads to physically wellposed solutions irrespective of the value of δl when the sliding distance over which relaxation occurs is such that the relaxation occurs over a finite dissipation length scale L d smaller than the size of the dissipation zone.

C O N C L U S I O N S
In this study, a systematic numerical parametric study of in-plane growing rupture propagating along a bimaterial planar interface under linear-weakening friction was performed when the physical problem is regularized with a Prakash-Clifton type of law.
The existence of mesh-independent and physically well-posed numerical solutions, the latter defined here as numerical solutions not depending on the regularization parameters, was investigated as a function of the regularization parameters themselves.
When regularization involves a dynamic time inversely proportional to the slip rate, numerical solutions are shown to become mesh independent as the mesh is refined, and physically well-posed, irrespective of the value of the impedance contrast, when relaxation sliding distance is smaller than the frictional weakening sliding distance, that is, δl ≤ 10%D c .The regularization can be interpreted as an adaptive low-pass filter that smoothens high frequency normal traction perturbations without significantly affecting the energy scale transfer that drives the rupture propagation.Slip-rate dependence of the relaxation time allows for an adaptive regularization all along the acceleration phase following the elastodynamic Lorentz contraction of the dissipation zone.
A relaxation time proportional to the evolving size of dissipation zone, normalized by a reference rupture speed, leads to an adaptive filter that smoothens the slip induced normal traction perturbations at the scale of the dissipation zone.Numerical solutions for this regularization are mesh-independent and physically well posed when the relaxation occurs over a length scale smaller that the dissipation zone.They are shown to superimpose to those obtained using a slip-rate dependent relaxation time when the local slip-rate is Downloaded from https://academic.oup.com/gji/article/209/1/48/2918650 by CNRS -ISTO user on 10 August 2022 now replaced by the maximum slip-rate δv max along the sliding interface, which always occurs in the vicinity of the rupture front.Therefore the size of the dissipation zone is shown to be inversely proportional to the maximum amplitude of slip rate all along the acceleration phase of the rupture.
When regularization involves a constant relaxation timescale, numerical solutions are shown to become mesh-independent as mesh is refined but not physically well-posed.In this case, the cut-off frequency of the filter associated to the regularization is fixed during the whole acceleration phase, whereas both amplitude and frequency of the normal traction perturbations increase with the propagation distance together with a continuous transfer of energy to shorter wavelengths as a result of the shrinking of the dissipation zone.As a result, numerical solutions might be physically well posed at the beginning of acceleration phase, but sooner or later the regularization will over-filter the physical short wavelengths that drive the rupture propagation.
For interfaces, across which bimaterial contrasts are such that no generalized Rayleigh wave exists, the rupture is shown to tend towards an asymptotic speed higher than the shear wave speed in the more compliant material.Peculiar effects are numerically observed such as emission of a half Mach cone in the more compliant medium, variations in the normal traction perturbations pattern along the favoured direction, together with a small amplitude extensive normal traction perturbation behind the rupture tip that modifies the slip gradient behind the dissipation zone.Elastodynamic coupling between interfacial sliding and normal traction perturbations mostly occurs behind the crack front and does not significantly perturb the normal traction variations ahead of the rupture.Therefore the rupture may propagate at an almost constant speed during a long time.This feature only depends on the shear velocity contrast between the two materials while the asymptotic speed slightly increases as the density contrast increases, as previously pointed by Rubin & Ampuero (2007).The asymptotic rupture speeds are in agreement with the analytical results of Ranjith & Rice (2001).
When bimaterial contrasts are such that generalized Rayleigh wave exists, the rupture asymptotically tends to C gr .However numerical solutions exhibit increasing spurious oscillations and become unstable irrespective of the value of the regularization pa-rameter.This instability is associated with the singular behaviour of the slip-induced normal traction perturbations, and of the slip rate at the rupture front, in relation with the complete shrinking of the dissipation zone.This singular behaviour cannot be regularized by a slip-rate dependent relaxation since the upper cut-off frequency of the associated adaptive filter moves beyond the mesh cut-off frequency.When a constant relaxation time is used before the emergence of the instability a metastable finite-size slip pulse can be observed.The onset of this slip pulse is mesh-independent but depends as well as the pulse size on the regularization parameter.As such this slip pulse cannot be considered as a physically well-posed solution.
Results of this study confirm that a Prakash-Clifton type of law under linear slip-weakening friction law does not regularize the numerical problem as the rupture speed tends towards the asymptotic generalized Rayleigh wave speed C gr , and this remains a critical issue.
Other studies have argued that rate-and-state or more complex slip velocity dependent friction laws may regularize the physical and numerical problem.In a quasi-static regime analysis, it was shown (Rice et al. 2001) that instantaneous strengthening response, can overcome, at small slip velocities, the bimaterial destabilizing effect associated with the coupling between interfacial slip and normal traction perturbation, This may occur even if the interface is steadily rupturing under a slip weakening law.
Despite recent progress (Lapusta et al. 2000;Kozdon & Dunham 2013), further analytical and numerical investigations of the full elastodynamic stability analysis in the framework of rate-and-state friction models remain to be done especially in the high slip velocity regime.It has also been shown that dry frictional interfaces generally become velocity strengthening over some range of slip velocities (e.g.Marone et al. 1991;Weeks 1993;Bar-Sinai et al. 2014, 2015).Such a velocity-strengthening effect is expected to play a stabilization role for bimaterial sliding as recently argued by Brener et al. (2015) in the framework of generalized rate-and-state friction models together with finite-time response to normal traction perturbations.However this remains to be further numerically investigated for accelerating rupture along a bimaterial interface and extended propagation distance.
Other approaches have been proposed (Harris & Day 2005) introducing an additional interfacial viscous dissipation term to regularize the bimaterial problem, extending previous works on homogeneous rupture problem (Dalguer & Day 2007;Rojas et al. 2009).Although the introduction of an artificial interfacial viscous damping limits the emergence of spurious oscillations at high frequency, how the energy transfer to smaller resolved wavelengths modifies the rupture dynamics still needs further investigation.
Another direction is to introduce a new length scale in the physical problem.It is well known that the region near a fault experiences a local environment different from the bulk.A more accurate approach would incorporate a description of the separate mechanics of the material interface.One possibility is to refine the problem using version of the surface model proposed by Gurtin & Murdoch (1975) as argued by Ru (2010) and Kim et al. (2011).When the existence of a thin intermediate zone more compliant than the two dissimilar materials is taken into account (e.g.Ben-Zion & Huang 2002) another possibility would be to consider the effect of the mechanical properties of such soft imperfect, or non-ideal, interfaces, (e.g.Atkinson 1977;Benveniste & Miloh 2001).
In the latter case when the thickness of the soft linearly elastic intermediate zone is essentially much smaller than the rupture propagation distance, asymptotic analysis (e.g.Mishuris 2004) shows Downloaded from https://academic.oup.com/gji/article/209/1/48/2918650 by CNRS -ISTO user on 10 August 2022 that displacement field is no more continuous across the interface ahead of the rupture tip and that tractions become proportional to the displacement discontinuity.Soft imperfect interfaces can be then replaced by imperfect transmission conditions, removing square-root singularity at the crack tip.Viscoelastic or dissipative intermediate zone would lead to fractional energy loss across the imperfect interface (e.g.Carcione 1996Carcione , 2001)), which together with the additional length scale may play a stabilization role for bimaterial sliding in the framework of a material interface law (e.g.Del Piero & Raous 2010).) and Q = 1 2 tC .The diagonality of the matrices Q and C implies that the eqs (A4) and (A6) are local to the collocation points on the fault and they can be directly combined with the Coulomb and Signorini conditions respectively.
Detailing the formulation for a 2-D inplane rupture, the total tractions as projected onto the normal and tangential components are σ  p) , where σ n 0 and τ 0 represent the pre-stress field components, and T n and T t are the normal and tangential components of the dynamic traction T.
To solve the Signorini condition at the point K on the fault we consider the total normal traction at zero normal slip: ) is the solution for the couple (δu n,K ( p+1) , σ n,K ( p+1) T ).In the case for which eq. (A7) yields σ n,K ( p+1) T > 0, the normal traction is forced to be zero (σ n,K ( p+1) T = 0) and δu n,K ( p+1) = Q K K σ n,K ( p+1) T .
From the Signorini condition, when opening occurs (σ n,K ( p+1) T = 0) the two sides of the fault behave as a free surface.When the interface is still in contact, the Coulomb condition should be verified to find the tangential slip rate.
Nevertheless, for a bimaterial interface, the Coulomb law is coupled with the effective normal traction, which is related to the normal traction through the discrete version of eq. ( 5): p+1) . (A8) Since the static contributions of σ n and σ eff coincide, eq.(A8) can be solved for the dynamic part of the effective normal traction.The tangential projection of eq. ( A4) is finally coupled with the Coulomb condition that can be solved by analogy with the procedure adopted for the Signorini condition.We start assuming that the tangential slip velocity δv t ( p+1) = 0. Accordingly, the corresponding tangential traction from (A4) is If this quantity is below the threshold | τ K ( p+1) T | ≤ −μ(δu t,K ( p+1) ) σ

Figure 1 .
Figure1.The simulation setup for the numerical models: the fault line is the thick black line in the middle whereas its prolongations on both sides form the fictitious boundary f (see Appendix A).We also represent the position of receivers along favoured direction (to the right of nucleation) where kinematic and dynamic quantities are collected as a function of time.Within the insets the expected linear slip weakening for the strength in homogeneous case (dotted lines) is compared to the typical shapes of the strength weakening for a bimaterial simulation.

Figure 2 .
Figure 2. Analysis of the slip rate for grid refinement.At the same time step for δl = 2%D c : (a) when solutions are not convergent strong oscillations of slip rate emerge up to pathological effects (e.g.stop and go of rupture).Those effects can boost the rupture producing spurious acceleration of the rupture front.The black square indicates the zoom around the crack front (b).(c) Even for highest δl, for which the oscillatory effects are damped, the solutions for coarsest meshes do not converge to those obtained from finest ones.When solutions converge, position of crack front and amplitude of the maximum coincide.The black square indicates the zoom around the crack front (d).

Figure 3 .
Figure 3. σ eff as a function of time for all used grid sizes and for two different relaxation slip values.(a) δl is equal to 2%D c : the coarser meshes show strong oscillations and they are not convergent with the results coming from the finer grids.The zoom shows the same quantities around the crack front (b).(c) δl is equal to 10%D c : the coarser meshes show less strong oscillations but they are still not convergent to the results coming from the finest meshes.The zoom again shows the same quantities around the crack front (d).

Figure 4 .
Figure 4. Physical convergence for decreasing relaxation slip in time domain.σ eff is shown at receiver 5 (a).The induced perturbations are huger and sharper moving away from the nucleation and they are smoothed for increasing δl.The black square in (a) indicates the zoom around the crack front (b) the convergence of maximum amplitude for σ eff is evident for δl ≤ 10%D c .(c) Rupture speed along the favoured direction as a function of distance from nucleation, the overlapping of convergent curves is evident as the capability of the rupture to accelerate almost up to C gr .The non-convergent solutions are clearly slower than the convergent ones.
a) shows the decrease of the relative error as δl decreases.The relative error becomes smaller than a 5 per cent threshold value when δl ≤ 10%D c .Supporting Information Fig. S2(b) shows how the normalized time-shift decreases when decreasing δl.Curves for δl ≤ 10%D c , are mesh independent and therefore numerically well-posed solutions, but they depend on the regularization parameter and therefore each of them are numerical solutions of different physical problems.Fig. 4(c) represents the rupture speed, normalized by C gr , during the acceleration phase along the favoured direction as a function of the distance from the centre of the initiation patch.The physical convergence is achieved again for δl ≤ 10%D c .The figure also indicates that the rupture is accelerating towards the expected asymptotic speed C gr .The same results hold during the whole acceleration phase and Supporting Information Figs S3(a) and (b) show the evolution with time of the effective normal traction at receivers 2 and 8. Gathering the results for the whole rupture propagation, solutions can be shown to be physically convergent for δl ≤ 10%D c .

Figure 5 .
Figure5.Physical convergence for decreasing relaxation slip in frequency domain: amplitude spectra for the perturbations described in Fig.4.The characteristic frequencies of the physical and numerical problem are explicitly reported as an example (dashed lines).The magenta dotted line is the characteristic cut-off frequency deriving from the relaxation mechanism (for δl = 10%D c ).The numerical limit (green dashed line) is related to the mesh with h = 4 m.

Figure 6 .
Figure 6.Normal traction perturbations accordingly to slip rate variations at the same receiver of Figs 4 and 5.The physical content of the two variations is the same for the two quantities and it increases moving away from the nucleation.The plot refers to the same simulation performed with a mesh size h = 4 m and δl = 10%D c .The coupling time and the physical time interval from which the respective frequencies are inferred are explicitly shown (cf.Fig. 5).

Figure 7 .
Figure 7. σ n − σ eff as a function of slip.Warm colour full lines in panel (a) represent simulations with δl = 10%D c (convergent solutions) at two receivers, whereas cold colour full lines in panel (b) are relative to δl = 50%D c (non-convergent solutions) at the same receivers.The zero crossing recorded at other receivers is plotted with dashed lines and respectively with warm and cold colours in panels (a) and (b).

Figure 8 .
Figure 8. Physical non-convergence for decreasing relaxation slip in time domain.σ eff is shown at the receiver 2 (a) and 8 (c) and the lack of convergence is evidenced for constant timescale even for small t c .Even when solutions are similar at the beginning of acceleration phase (inset of panel a) the differences increase with the crack growth as shown in the inset of panel (c).The physical non-convergence is shown in the frequency domain at the same receivers (panels b and d).The non-convergence of solutions can be argued by the increasing difference among the low-frequency part of amplitude spectra.The cut-off frequency (related to t c = 6.0 × 10 −4 s) is fixed (yellow dotted lines).The physical (black dotted lines) and the coupling frequencies (red dotted lines) increase as expected with the crack growth.

Figure 9 .
Figure9.Acceleration of the rupture towards the asymptotic speed (C gr ) for dynamic and constant timescales: while a convergence among different timescales is still detectable, the acceleration is equivalent to that deriving from dynamic timescale models.Conversely when the timescale is too large to properly regularize the problem the rupture speeds differ more and more and the acceleration is slower for the constant timescale.

Figure 10 .
Figure 10.Acceleration of rupture, along favoured direction, for high contrasts of impedance obtained varying C s 2 .(a) Under these conditions the rupture can accelerate towards speeds higher than C s 2 .(b) Average speed (normalized to C s 2 ) during stationary phase as a function of increasing contrast of impedance γ .

Figure 11 .
Figure 11.Kinetic energy field before (a) and after (b) the acceleration of the rupture beyond C s 2 .The emission of S-waves behind the rupture front in more compliant medium (Mach-cone) is evident in panel (b).Perturbation of normal traction (c) at the same time step of panel (b): the extensive effect due to the S wave is evidenced (black circle).

Figure 12 .
Figure 12.Slip rate profile at time t 0 = 0.19 s just before the blowing up of numerical simulation (red curve) is plotted together with the number of points in dissipation zone as function of the position of crack front.When the dissipation zone shrinks at less of 3 points a slip pulse is generated and after few iterations the models blow up (grey curve).

Figure 13 .
Figure 13.Pulse onset for different constant timescale: the smallest t c leads to a faster generation of the pulse (a) after a shorter propagation distance.(b)Pulse size for different models obtained for constant timescale regularizations.The sizes, at different time steps, are compared for different t c .In all cases the pulse size is pretty constant during its propagation and when a larger t c is used the average size of the pulse is quite larger.

Figure 14 .
Figure 14.Convergence analysis for decreasing β in time domain.(a) Variations of σ eff as a function of time are plotted for receiver 5.The inset of panel (a) shows the overlapping of the curves and the physical convergence for β ≤ 0.10.(b) Convergence analysis for decreasing β in frequency domain for the same receiver in (a): the physical convergence of solutions for β ≤ 0.10 is still due to the overlapping of amplitude spectra in low frequency band as for the dynamic timescale.(c,d) σ n − σ eff as a function of slip: warm colour full lines in panel (c) represent simulations with β = 0.10 (convergent solutions) at two receivers, whereas cold colour full lines in panel (d) are relative to β = 0.50 (non-convergent solutions) at the same receivers.The zero crossing recorded at other receivers is plotted with dashed lines and respectively with warm and cold colours in panels (c) and (d).

Figure 15 .
Figure 15.Slip rate profiles at a fixed time step for dissipation length scale (L d in the legend) and maximum slip rate scale (δv max in the legend).The inset shows that the two regularizations are convergent in the sense of crack front position and maximum of slip rate amplitude.

Fint
the Signorini condition is automatically satisfied and (0, σ n,K ( p+1) T is assumed.During sliding, the friction drops linearly from a static value μ s to a dynamic value μ d over a characteristic sliding distance D c