TY - JOUR
AU - Vilotte,, Jean-Pierre
AB - Summary Earthquake ruptures often develop along faults separating materials with dissimilar elastic properties. Due to the broken symmetry, the propagation of the rupture along the bimaterial interface is driven by the coupling between interfacial sliding and normal traction perturbations. We numerically investigate in-plane rupture growth along a planar interface, under slip weakening friction, separating two dissimilar isotropic linearly elastic half-spaces, and we perform a parametric study of the classical Prakash–Clifton regularization, for different material contrasts. In particular the mesh-dependence and the regularization-dependence of the numerical solutions are analysed in this parameter space. When the regularization involves a slip-rate dependent relaxation time, a characteristic sliding distance is identified below which numerical solutions no longer depend on the regularization parameter, that is, they are physically well-posed solutions. Such regularization provides an adaptive high-frequency filter of the slip-induced normal traction perturbations, following the dynamic shrinking of the dissipation zone during the acceleration phase. In contrast, a regularization involving a constant relaxation time leads to numerical solutions that always depend on the regularization parameter since it fails in adapting to the shrinking of the process zone. Dynamic regularization is further investigated using a non-local regularization based on a relaxation time that depends on the dynamic length of the dissipation zone. Such reformulation is shown to provide similar results as the dynamic timescale regularization proposed by Prakash–Clifton when the slip rate is replaced by the maximum slip rate along the sliding interface. This leads to the identification of a dissipative length scale associated with the coupling between interfacial sliding and normal traction perturbations, together with a scaling law between the maximum slip rate and the dynamic size of the process zone during the rupture propagation. Dynamic timescale regularization provides mesh-independent and physically well-posed numerical solutions during the acceleration phase towards an asymptotic speed. When generalized Rayleigh wave does not exist, numerical solutions are shown to tend towards an asymptotic velocity higher than the slowest shear wave speed. When the generalized Rayleigh wave speed exists, numerical solutions tend towards this velocity becoming noisier and noisier as the rupture progresses. In this regime regularization dependent, unstable finite-size pulses may be generated. Numerical solutions, Earthquake dynamics, Earthquake source observations, Computational seismology 1 INTRODUCTION Seismic ruptures often propagate along interfaces between dissimilar materials. Yet a complete theoretical and numerical understanding of the dynamics and stability of rapidly propagating rupture along frictional bimaterial interfaces is lacking. Field observations, laboratory experiments, analytical studies, and numerical simulations evidenced a number of specific features of rupture propagation along interfaces separating dissimilar materials, as compared to interfaces separating similar materials. Frictional rupture propagating along a planar bimaterial interface leads to slip-induced normal traction perturbations, as a result of symmetry breaking, which depends on the material property contrast and on the frictional behaviour of the interface. Geological observations at the San Andreas Fault highlighted favoured direction of propagation during large rupture events, that 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 Cgr 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 Cgr, 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, 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 2009, 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 Cgr 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. Mesh dependence can be removed (Cochard & Rice 2000) when the friction is modified with the Prakash–Clifton regularization (Prakash & Clifton 1993; Prakash 1998). 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 Dc 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 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 Cgr (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 Cgr, 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 Cgr, 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. 2 SEISMIC RUPTURES ALONG BIMATERIAL INTERFACES Seismic waves travelling in elastic media are governed by the following equations: \begin{equation}\!\!\left\{ {\begin{array}{@{}*{1}{l}@{}} {\rho {{\bf \dot{v}}} = \vec{\nabla } \cdot {{\boldsymbol \sigma }}}\\ {{{\boldsymbol{\dot{\sigma }}}} = {{\bf c}}:\vec{\nabla }{{\bf v}}} \end{array}} \right.,\end{equation} (1) where ρ is the material density, v the particle velocity, |${{\boldsymbol \sigma }}$| 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,|${{\bf T}} = {{\boldsymbol \sigma }} \cdot {{\bf 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 = u2 − u1 and the slip rate δv = v2 − v1, 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, ξ = ξnn + ξt where the superscripts n and t refer to the normal and tangential contributions, with ξn = ξ · n and ξt = ξ −(ξ · n)n. 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 \begin{equation}\!\!\left\{ {\begin{array} {c} {\delta {u^n} \ge 0}\\ {{T^n} \le 0}\\ {\delta {u^n}{T^n} = 0.} \end{array}} \right.\end{equation} (2) This condition ensures the impenetrability of the two bodies when being in contact along the interface under compressive regime (Tn < 0,  δun = 0). When the two bodies are no longer in contact, opening occurs (δun ≥ 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 \begin{equation}\!\!\left\{ {\begin{array} {c} {\left( {\left| {{{{\bf T}}^t}} \right| + \mu {T^n}} \right)\left| {\delta {{{\bf v}}^t}} \right| = 0}\\ {\left| {{{{\bf T}}^t}} \right| + \mu {T^n} \le 0}\\ {{{\rm T} ^t} \cdot \delta {{{\bf v}}^t} = \left| {{{{\bf T}}^t}} \right|\left| {\delta {{{\bf v}}^t}} \right|} \end{array}} \right.\end{equation} (3) where μ is the actual friction coefficient and μTn the frictional strength. First condition in eq. (3) indicates that the non-sliding to sliding transition occurs when the modulus of the tangential traction Tt lies on the surface of the actual frictional cone defined by the second condition in eq. (3). During sliding, Tt 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 Tt remains inside the frictional cone, no sliding occurs, that is, δv = 0. In this study, a slip-weakening friction (Ida 1972) is assumed. During sliding, the friction drops linearly from a static value μs to a dynamic value μd over a characteristic sliding distance Dc according to \begin{equation}\mu = \max \left\{ {\left( {{\mu _s} - \frac{{{\mu _s} - {\mu _d}}}{{{D_c}}}\left| {\delta {u^t}} \right|} \right);{\mu _d}} \right\}.\end{equation} (4) 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 ill-posed 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_{\rm{eff}}^n$| that is related to the normal traction by an exponential relaxation law \begin{equation} \left\{ \begin{array} {c@\quad l} \frac{{\partial T_{{\rm{eff}}}^n}}{{\partial t}} = \frac{1}{{{t^*}}}\left( {{T^n} - T_{{\rm{eff}}}^n} \right) & {\rm if}\,\,\,\left| {\delta {{\bf v}}} \right| \ne 0\\ T_{\rm{eff}}^n = {T^n} & {\rm if}\,\,\left| {\delta {{\bf v}}} \right| = 0 \end{array} \right.\end{equation} (5) 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: \begin{equation}\!\left( {\left| {{{{\bf T}}^t}} \right| + \mu T_{\rm{eff}}^n} \right)\left| {\delta {{{\bf v}}^t}} \right| = 0;\quad \left| {{{{\bf T}}^t}} \right| + \mu {T_{{\rm{eff}}}} \le 0.\end{equation} (6) 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*(|δvt|) defined as \begin{equation}\frac{1}{{{t^*}}} = \frac{{\left| {\delta {{{\bf v}}^t}} \right|}}{{\delta l}} + \frac{1}{{{t_c}}}\quad {\rm{if}}\,\,\left| {\delta {{\bf v}}} \right| \ne 0.\end{equation} (7) where δl is a characteristic sliding distance and tc 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 Dc. 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 td = δl/|δvt| as the dynamic relaxation timescale, and tc as the constant relaxation timescale. A parametric study is performed to define numerically well-posed, as the solutions showing 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 |δvt| 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. 3 SIMULATION SETUP 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 |$\gamma = {C_{{s_1}}}/{C_{{s_2}}}$|⁠, and less influenced by the density ratio (Ben-Zion & Andrews 1998; Rubin & Ampuero 2007). Figure 1. Open in new tabDownload slide 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 1. Open in new tabDownload slide 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. 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 \begin{equation}\!\left\{ {\begin{array}{@{}*{1}{c}@{}} {\Delta \tau \left( x \right) = \frac{{\bar{\zeta }}}{{2\pi }}\int_{ - \infty }^\infty {\frac{{{\rm{d}}\delta u/{\rm{d}}x}}{{x - s}}{\rm{d}}s} }\\ {\Delta {\sigma ^n}\left( x \right) = {\zeta ^*}\frac{{{\rm{d}}\delta u}}{{{\rm{d}}x}}} \end{array}} \right.\end{equation} (8) where the moduli |$\bar{\zeta }$| 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 |$\bar{\zeta }$| 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 (Cgr) and it is the expected asymptotic speed for the growing rupture. The explicit expression for |$\bar{\zeta }$| 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, Cgr 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. 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 σ0n. The shear strength exhibits therefore at all points of the interface the same linear weakening as the friction. As a consequence, the fracture energy Gc, defined as the area under the strength-weakening curve in a slip-traction plot, can be simply evaluated as Gc = 0.5Dc(μs − μd)σn0, 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 Gc 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 Cgr. On the fault, the initial normal traction is set to the uniform value σ0n = −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 Dc = 6 mm. Accordingly the strength parameter |$s = \frac{{{\mu _s}{\sigma _0}^n - \tau }}{{\tau - {\mu _d}{\sigma ^n}_o}} = 3$|⁠. 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 Lc for a linear slip-weakening rupture (Uenishi & Rice 2003; Rubin & Ampuero 2007): \begin{equation}{L_c} = 1.118\frac{{\zeta '}}{W}\end{equation} (9) where the effective elastic modulus for quasi-static plane strain deformation ζ' depends on the shear moduli |${\rho _i}{C^2}_{{s_i}}$| and the Poisson coefficients νi of the two media and it is defined as |$\bar{\zeta }$| 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: \begin{equation}W = \frac{{{\sigma _0}^n\left( {{\mu _s} - {\mu _d}} \right)}}{{{D_c}}}\end{equation} (10) 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. 4 REGULARIZATION: A PARAMETRIC STUDY 4.1 Regularization with slip-rate dependent relaxation time We refer to a dynamic relaxation timescale when t* = td = δ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%Dc − 100%Dc) 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 δ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 δlmax  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 L2-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 cross-correlation before computing the L2-norm difference. The relative error is finally obtained by normalizing this result by the L2-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\,{\rm{km}}\,{{\rm{s}}^{ - 1}}$| and |${C_{{s_1}}} = 3.092\,{\rm{km }}\,{{\rm{s}}^{ - 1}}$|⁠) is investigated. Fig. 2(a) shows the slip-rate profile at time step t0 = 0.12 s for δl = 2%Dc 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. Figure 2. Open in new tabDownload slide Analysis of the slip rate for grid refinement. At the same time step for δl = 2%Dc: (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 2. Open in new tabDownload slide Analysis of the slip rate for grid refinement. At the same time step for δl = 2%Dc: (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). 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%Dc 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%Dc. 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%Dc (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. Figure 3. Open in new tabDownload slide σeff as a function of time for all used grid sizes and for two different relaxation slip values. (a) δl is equal to 2%Dc: 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%Dc: 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 3. Open in new tabDownload slide σeff as a function of time for all used grid sizes and for two different relaxation slip values. (a) δl is equal to 2%Dc: 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%Dc: 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). 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 t0 = 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%Dc. 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%Dc. 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 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(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%Dc. Supporting Information Fig. S2(b) shows how the normalized time-shift decreases when decreasing δl. Curves for δl ≤ 10%Dc, 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 Cgr, 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%Dc. The figure also indicates that the rupture is accelerating towards the expected asymptotic speed Cgr. 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%Dc. Figure 4. Open in new tabDownload slide 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%Dc. (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 Cgr. The non-convergent solutions are clearly slower than the convergent ones. Figure 4. Open in new tabDownload slide 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%Dc. (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 Cgr. The non-convergent solutions are clearly slower than the convergent ones. In the case of the dynamic timescale relaxation, numerical solutions are shown to become independent of the regularization parameter 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 Dc. 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. Figure 5. Open in new tabDownload slide 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%Dc). The numerical limit (green dashed line) is related to the mesh with h = 4 m. Figure 5. Open in new tabDownload slide 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%Dc). The numerical limit (green dashed line) is related to the mesh with h = 4 m. 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%Dc 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. Figure 6. Open in new tabDownload slide 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%Dc. The coupling time and the physical time interval from which the respective frequencies are inferred are explicitly shown (cf. Fig. 5). Figure 6. Open in new tabDownload slide 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%Dc. The coupling time and the physical time interval from which the respective frequencies are inferred are explicitly shown (cf. Fig. 5). 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%Dc, it is found that δu* < Dc, as shown in Fig. 7(a) where δu* ∼ 95%Dc for δl = 10%Dc. When δl = 15%Dc, δu* ∼ Dc, whereas for larger values of δl, δu* > Dc, for example, see Fig. 7(b) where δu* ∼ 110%Dc for δl = 50%Dc. 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%Dc, 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 Dc, which is the characteristic finite sliding distance over which the weakening of the friction takes place. Figure 7. Open in new tabDownload slide σn − σeff as a function of slip. Warm colour full lines in panel (a) represent simulations with δl = 10%Dc (convergent solutions) at two receivers, whereas cold colour full lines in panel (b) are relative to δl = 50%Dc (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 7. Open in new tabDownload slide σn − σeff as a function of slip. Warm colour full lines in panel (a) represent simulations with δl = 10%Dc (convergent solutions) at two receivers, whereas cold colour full lines in panel (b) are relative to δl = 50%Dc (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). 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%Dc. 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* < Dc the physical convergence is achieved both in time and frequency. This condition holds until the rupture approaches the asymptotic speed (see Section 5). 4.2 Regularization with constant relaxation time When a constant relaxation timescale tc 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 tc in the range 1.2 × 10− 4 ≤ tc ≤ 6.0 × 10− 3. Defining |${t_c} = \frac{{\delta l}}{{\delta {v^*}}}$|⁠, as proposed by Cochard & Rice (2000), corresponds to a range of variation of δl between 2%Dc and 100%Dc 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 t0 = 0.12 s in the space domain for the smallest value of tc (tc = 1.2 × 10− 4 s, Supporting Information Figs S5a and b) and for tc = 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 tc. 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 tc 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). Figure 8. Open in new tabDownload slide 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 tc. 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 tc = 6.0 × 10− 4s) 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 8. Open in new tabDownload slide 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 tc. 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 tc = 6.0 × 10− 4s) is fixed (yellow dotted lines). The physical (black dotted lines) and the coupling frequencies (red dotted lines) increase as expected with the crack growth. This lack of physical convergence can be explained in comparison with the results obtained for the dynamic regularization. An equivalent |$\delta l_{\rm{eq}}^{\max }$| can be estimated, all along the rupture propagation, in the case of a constant relaxation time as |$\delta l_{\rm{eq}}^{\max } = \delta {v_{\max }} \cdot {t_c}$|⁠, where δvmax  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 tc = 6 × 10− 4 s the equivalent sliding distance at receiver 2 is |$\delta l_{{\rm{eq}}\# 2}^{\max }\,\, \sim 6 \times {10^{ - 4}}\,{\rm{m}} = 10\% {D_c}$| while for receiver 8 it is |$\delta l_{{\rm{eq}}\# 2}^{\max }\,\, \sim 2.4 \times {10^{ - 3}}\,{\rm{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 tc, |$\delta l_{{\rm{eq}}\# 8}^{\max }$| 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 tc. In this case, the regularization will over-filter 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 the instantaneous rupture speed (normalized to Cgr) is shown as the rupture propagates along the fault for the constant (blue line) and dynamic (red line) relaxation time mechanisms. For fixed tc, the solutions are initially superimposed. During the crack growth, when 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. Figure 9. Open in new tabDownload slide Acceleration of the rupture towards the asymptotic speed (Cgr) 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 9. Open in new tabDownload slide Acceleration of the rupture towards the asymptotic speed (Cgr) 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. In light of this parametric study, the classical Prakash–Clifton regularization (eq. (5)) where t* is given by eq. (7) can be rewritten as \begin{equation}\frac{{\partial {\sigma _{{\rm{eff}}}}}}{{\partial t}} = \left( {{f_d} + {f_c}} \right)\left( {{\sigma _n} - {\sigma _{{\rm{eff}}}}} \right),\end{equation} (11) where fd = |δv|/δl is a dynamic frequency and fc is a constant cut-off 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 well-posed solutions, that is, mesh independent solutions, but does not lead to physically convergent solutions. When fc is small enough and the absolute value of the slip rate is large enough to guarantee fd ≫ fc, 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 fc ≥ fd 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 nucleation processes, even though the choice of fc becomes nucleation dependent. 5 RUPTURE BEHAVIOUR TOWARDS ASYMPTOTIC SPEED 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 mesh-independent 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 Cgr. When 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). 5.1 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\,{\rm{m }}{{\rm{s}}^{ - 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. Figure 10. Open in new tabDownload slide 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 10. Open in new tabDownload slide 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 γ. 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. Figure 11. Open in new tabDownload slide 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 11. Open in new tabDownload slide 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). For each of the investigated values of γ, variations of the density contrast have been explored such that Cgr 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 |$\gamma = {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. 5.2 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 Cgr. In particular the acceleration phase for γ = 1.18 and Cgr = 2570 ms− 1 was analysed in Section 4.1. The rupture speed is shown to tend towards the asymptotic Cgr 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 Cgr, 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 Cgr 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. Figure 12. Open in new tabDownload slide Slip rate profile at time t0 = 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 12. Open in new tabDownload slide Slip rate profile at time t0 = 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). During the initial phase of the instability, a slip pulse emerges (red curve in Fig. 12), which might be triggered by the occurrence 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 Cgr, 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 tc have been explored spanning two orders of magnitude. While instability can be delayed when increasing tc, that is, the larger is tc 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 Cgr 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 mesh-independent but depends on the constant relaxation time tc and it occurs earlier for smaller tc (Fig. 13a). Figure 13. Open in new tabDownload slide Pulse onset for different constant timescale: the smallest tc 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 tc. In all cases the pulse size is pretty constant during its propagation and when a larger tc is used the average size of the pulse is quite larger. Figure 13. Open in new tabDownload slide Pulse onset for different constant timescale: the smallest tc 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 tc. In all cases the pulse size is pretty constant during its propagation and when a larger tc is used the average size of the pulse is quite larger. 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 tc leads to a smaller pulse size as shown in Fig. 13(b). In particular defining the non-dimensional parameter |${\tilde{L}_p}$| as the pulse size normalized by the element size h, we found |${\tilde{L}_p} \approx 0.75$| for tc = 3.0 × 10− 4 s, |${\tilde{L}_p} = 1.50$| for tc = 6.0 × 10− 4 s, and |${\tilde{L}_p} = 2.50$| for tc = 1.2 × 10− 3s. 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 generalized Rayleigh wave exists, and mathematically ill-posed (Adams 1995; Martins et al. 1995; Simões & Martins 1998; Ranjith & Rice 2001). 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. 6 MODIFIED REGULARIZATION WITH DISSIPATION LENGTH DEPENDENT RELAXATION TIME 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 Dc. This aspect suggests that the instability of the solutions for a bimaterial interface under slip-weakening 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: \begin{equation}\frac{{\partial {\sigma _{{\rm{eff}}}}}}{{\partial t}} = \frac{1}{{{t_{{L_d}}}}}\left( {{\sigma ^n} - {\sigma _{{\rm{eff}}}}} \right)\quad {t_{{L_d}}} = \frac{{\beta {L_d}}}{{{V^a}}}\end{equation} (12) where now Ld is the size of dissipation zone, Va 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 well-posed 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 t0 = 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 slip-rate 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* < Dc for β ≤ 0.10 (Fig. 14c), while δu* > Dc 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. Figure 14. Open in new tabDownload slide 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 14. Open in new tabDownload slide 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). 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 δvmax . This value always occurs in the vicinity of the rupture front and it defines a characteristic velocity scale: \begin{equation}\frac{{\partial {\sigma _{{\rm{eff}}}}}}{{\partial t}} = \frac{1}{{{t_d}}}\left( {{\sigma ^n} - {\sigma _{{\rm{eff}}}}} \right)\quad {t_d} = \frac{{\delta l}}{{\delta {v_{\max }}}}.\end{equation} (13) The numerical solutions for the two regularizations superimpose as it can be observed in the slip rate profile at time step t0 = 0.13s (Fig. 15). Since the relaxation slip can be also expressed as δl = βDc the equivalence between the two regularizations also implies that the dissipation length scales as the inverse of the maximum slip rate: \begin{equation}{L_d} \propto \frac{1}{{\delta {v_{\max }}}}.\end{equation} (14) Figure 15. Open in new tabDownload slide Slip rate profiles at a fixed time step for dissipation length scale (Ld in the legend) and maximum slip rate scale (δvmax  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. Figure 15. Open in new tabDownload slide Slip rate profiles at a fixed time step for dissipation length scale (Ld in the legend) and maximum slip rate scale (δvmax  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. 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 well-posed 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 Ld smaller than the size of the dissipation zone. 7 CONCLUSIONS 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%Dc. 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 now replaced by the maximum slip-rate δvmax  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 Cgr. However numerical solutions exhibit increasing spurious oscillations and become unstable irrespective of the value of the regularization parameter. 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 Cgr, 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 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 1996, 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). Acknowledgments Part of the simulations in this work were performed on the S-CAPAD resources of IPG Paris with the support of Geneviíve Moguilny. The authors acknowledge Harsha Bhat Suresh (ENS Paris) for the fruitful scientific discussions about the paper topics and the two anonymous reviewers for their constructive review and discussions. The authors also acknowledge the Université Franco-Italienne and the project VINCI, as well as the project ReLUIS 2016 for the funding support. REFERENCES Adams G.G. , 1995 . Self-excited oscillations of two elastic half-spaces sliding with a constant coefficient of friction , J. Appl. Mech. , 62 , 867 – 872 . Google Scholar Crossref Search ADS WorldCat Adams G.G. , 1998 . Steady sliding of two elastic half-spaces with friction reduction due to interface stick-slip , J. Appl. Mech. , 65 , 470 – 475 . Google Scholar Crossref Search ADS WorldCat Adda-Bedia M. , Ben Amar M., 2003 . Self-sustained slip pulses of finite size between dissimilar materials , J. Mech. Phys. Solids , 51 , 1849 – 1861 . Google Scholar Crossref Search ADS WorldCat Ampuero J.-P. , Ben-Zion Y., 2008 . Cracks, pulses and macroscopic asymmetry of dynamic rupture on a bimaterial interface with velocity-weakening friction , Geophys. J. Int. , 173 ( 2 ), 674 – 692 . Google Scholar Crossref Search ADS WorldCat Andrews D.J. , 1976 . Rupture velocity of plane strain shear cracks , J. geophys. Res. , 81 , 5679 – 5687 . Google Scholar Crossref Search ADS WorldCat Andrews D.J. , Ben-Zion Y., 1997 . Wrinkle-like slip pulse on a fault between different materials , J. geophys. Res. , 102 , 553 – 571 . Google Scholar Crossref Search ADS WorldCat Atkinson C. , 1977 . On stress singularities and interfaces in linear elastic fracture mechanics , Int. J. Fract. , 13 , 807 – 820 . Google Scholar Crossref Search ADS WorldCat Audoly B. , 2000 . Asymptotic study of the interfacial crack with friction , J. Mech. Phys. Solids , 48 , 1851 – 1864 . Google Scholar Crossref Search ADS WorldCat Bar-Sinai Y. , Spatschek R., Brener E.A., Bouchbinder E., 2014 . On the velocity-strengthening behavior of dry friction , J. geophys. Res. , 119 ( B3 ), 1738 – 1748 . Google Scholar Crossref Search ADS WorldCat Bar-Sinai Y. , Spatschek R., Brener E.A., Bouchbinder E., 2015 . Velocity-strengthening friction significantly affects interfacial dynamics, strength and dissipation , Sci. Rep. , 5 , 7841 , doi:10.1038/srep07841 . Google Scholar Crossref Search ADS PubMed WorldCat Benveniste Y. , Miloh T., 2001 . Imperfect soft and stiff interfaces in two-dimensional elasticity , Mech. Mater. , 33 , 309 – 323 . Google Scholar Crossref Search ADS WorldCat Ben Zion Y. , Andrews D.J., 1998 . Properties and implications of dynamic rupture along a material interface , Bull. seism. Soc. Am. , 88 , 1085 – 1094 . OpenURL Placeholder Text WorldCat Ben-Zion Y. , Huang Y., 2002 . Dynamic rupture on an interface between a compliant fault zone layer and a stiffer surrounding solid , J. geophys. Res. , 107 ( B2 ), 1 – 13 . Google Scholar Crossref Search ADS WorldCat Brener E.A. , Weikam M., Spatschek R., Bar-Sinai Y., Bouchbinder E., 2015 . Dynamic instabilities of frictional sliding at a bimaterial interface , J. Mech. Phys. Solids , 89 , 149 – 173 . Google Scholar Crossref Search ADS WorldCat Burridge R. , 1973 . Admissible speeds for plane-strain self-similar shear cracks with friction but lacking cohesion , Geophys. J. R. astr. Soc. , 35 , 439 – 455 . Google Scholar Crossref Search ADS WorldCat Carcione J.M. , 1996 . Elastodynamics of a non-ideal interface: application to crack and fracture scattering , J. geophys. Res. , 101 ( B12 ), 28 177–28 188 . Google Scholar Crossref Search ADS WorldCat Carcione J.M. , 2001 . Wave fields in real media: wave propagation, in anisotropic, anelastic and porous media , in Handbook of Geophysics Exploration, Seismic Explorations , Vol. 31, eds Helbig K., Treitel S., Pergamon . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Chaljub E. , Komatitsch D., Vilotte J.-P., Capdeville Y., Valette J.-P., Festa G., 2007 . Spectral-element analysis in seismology , Adv. Geophys. , 48 , 365 – 419 . OpenURL Placeholder Text WorldCat Cochard A. , Rice J.R., 2000 . Fault rupture ill-posedness between dissimilar materials: ill-posedness, regularization, and slip-pulse response , J. geophys. Res. , 105 ( B11 ), 25 891–25 907 . Google Scholar Crossref Search ADS WorldCat Dalguer L.A. , Day S.M., 2007 . Staggered-grid split-node method for spontaneous rupture simulation , J. geophys. Res. , 112 ( B2 ), 1 – 15 . Google Scholar Crossref Search ADS WorldCat Del Piero G. , Raous M., 2010 . A unified model for adhesive interfaces with damage, viscosity and friction , Eur. J. Mech. A , 24 ( 4 ), 496 – 507 . Google Scholar Crossref Search ADS WorldCat Deng X. , 1994 . An asymptotic analysis of stationary and moving cracks with frictional contact along bimaterial interfaces and in homogeneous solids , Int. J. Solids Struct. , 31 ( 17 ), 2407 – 2429 . Google Scholar Crossref Search ADS WorldCat Dor O. , Rockwell T.K., Ben-Zion Y., 2006 . Geologic observations of damage asymmetry in the structure of the San Jacinto, San Andreas and Punchbowl faults in southern California: a possible indicator for preferred rupture propagation direction , Pure appl. Geophys. , 163 , 301 – 349 . Google Scholar Crossref Search ADS WorldCat Dor O. , Yildirim C., Rockwell T.K., Ben-Zion Y., Emre O., Sisk M., Duman T.Y., 2008 . Geologic and geomorphologic asymmetry across the rupture zones of the 1943 and 1944 earthquakes on the North Anatolian Fault: possible signals for preferred earthquake propagation direction , Geophys. J. Int. , 173 , 483 – 504 . Google Scholar Crossref Search ADS WorldCat Erdogan F. , 1965 . Stress distribution in bonded dissimilar materials with cracks , J. Appl. Mech. , 32 , 403 – 410 . Google Scholar Crossref Search ADS WorldCat Festa G. , 2004 . Slip imaging by isochrones back projection and source dynamics with spectral element methods , PhD Thesis , Università degli Studi di Bologna “Alma Mater Studiorum”, Bologna (Italy) & Università degli Studi di Napoli “Federico II” , Napoli, Italy . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Festa G. , Nielsen S.B., 2003 . PML absorbing boundaries , Bull. seism. Soc. Am. , 93 ( 2 ), 891 – 903 . Google Scholar Crossref Search ADS WorldCat Festa G. , Vilotte J.-P., 2005 . The Newmark scheme as velocity-stress time-staggering: an efficient implementation for spectral element simulations of elastodynamics , Geophys. J. Int. , 161 ( 3 ), 789 – 812 . Google Scholar Crossref Search ADS WorldCat Festa G. , Vilotte J.-P., 2006 . Influence of rupture initiation on the intersonic transition: crack-like versus pulse like modes , Geophys. Res. Lett. , 33 ( 15 ), L15320 , doi:10.1029/2006GL026378 . Google Scholar Crossref Search ADS WorldCat Gurtin M.E. , Murdoch A.I., 1975 . A continuum theory of elastic material surfaces , Arch. Ration. Mech. Anal. , 57 , 291 – 323 . Google Scholar Crossref Search ADS WorldCat Harris R.A. , Day S.M., 1997 . Effects of a low-velocity zone on a dynamic rupture , Bull. seism. Soc. Am. , 87 , 1267 – 1280 . OpenURL Placeholder Text WorldCat Harris R.A. , Day S.M., 2005 . Material contrast does not predict earthquake rupture propagation direction , Geophys. Res. Lett. , 32 , L23301 , doi:10.1029/2005GL023941 . Google Scholar Crossref Search ADS WorldCat Ida Y. , 1972 . Cohesive force across the tip of a longitudinal-shear crack and Griffith's specific surface energy , J. geophys. Res. , 77 , 3796 – 3805 . Google Scholar Crossref Search ADS WorldCat Kammer D.S. , Yartebov V.A., Anciaux G., Molinari J.-F., 2014 . The existence of a critical length scale in friction , J. Mech. Phys. Solids , 63 , 40 – 50 . Google Scholar Crossref Search ADS WorldCat Kim C.I. , Schiavone P., Riu C.-Q., 2011 . The effect of surface elasticity on an interface crack in plane deformations , Proc. R. Soc. A , 467 , 3530 – 3549 . Google Scholar Crossref Search ADS WorldCat Komatitsch D. , Vilotte J.-P., 1998 . The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures , Bull. seism. Soc. Am. , 88 , 368 – 392 . OpenURL Placeholder Text WorldCat Kozdon J.E. , Dunham E.M., 2013 . Rupture to the trench: dynamic rupture simulations of the 11 March 2011 Tohoku earthquake , Bull. seism. Soc. Am. , 103 , 1275 – 1289 . Google Scholar Crossref Search ADS WorldCat Langer S. , Olsen-Kettle L., Weatherley D., 2012 . Identification of supershear transition mechanisms due to material contrast at bimaterial faults , Geophys. J. Int. , 190 , 1169 – 1180 . Google Scholar Crossref Search ADS WorldCat Lapusta N. , Rice J.R., Ben-Zion Y., Zheng G., 2000 . Elastodynamic analysis for slow tectonic loading with spontaneous rupture episodes on faults with rate- and state-dependent friction , J. geophys. Res. , 105 , 23 765–23 789 . Google Scholar Crossref Search ADS WorldCat Marone C.J. , Scholz C.H., Bilham R., 1991 . On the mechanics of earthquake afterslip , J. geophys. Res. , 96 ( B5 ), 8441 – 8452 . Google Scholar Crossref Search ADS WorldCat Martins J.A.C. , Simões F.M.F., 1995 . On some sources of instability/ill-posedness in elasticity problems with Coulomb friction , in Contact Mechanics , pp. 95 – 106 , eds Raous M., Jean M., Moreau J.J., Plenum Press . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Martins J.A.C. , Guimaraes J., Faria L.O., 1995 . Dynamic surface solutions in linear elasticity and visco-elasticity with frictional boundary conditions , Trans. ASME. J. Vib. Acoust., 117 , 445 – 451 . OpenURL Placeholder Text WorldCat Mishuris G. , 2004 . Imperfect transmission conditions for a thin weakly compressible interface. 2D problems , Arch. Mech. , 56 ( 2 ), 103 – 115 . OpenURL Placeholder Text WorldCat Muskhelishvili N.I. , 1953 . Some Basic Problems of the Mathematical Theory of Elasticity , P. Noordhoff . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Prakash V. , 1998 . Frictional response of sliding interfaces subjected to time varying normal pressures , Trans. ASME , J. Tribol. , 120 , 97 – 102 . Google Scholar Crossref Search ADS WorldCat Prakash V. , Clifton R.J., 1993 . Time resolved dynamic friction measurements in pressure-shear , in Experimental Techniques in the Dynamics of Deformable Solids , AMD-Vol. 165 , pp. 33 – 48 , ed. Ramesh K.T., ASME . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Ranjith K. , 2009 . Destabilization of long-wavelength Love and Stoneley waves in slow sliding , Int. J. Solids Struct. , 46 , 3086 – 3092 . Google Scholar Crossref Search ADS WorldCat Ranjith K. , 2014 . Instabilities in dynamic anti-plane sliding of an elastic layer on a dissimilar elastic half-space , J. Elast. , 115 , 47 – 59 . Google Scholar Crossref Search ADS WorldCat Ranjith K. , Rice J.K., 2001 . Slip dynamics at an interface between dissimilar materials , J. Mech. Phys. Solids , 49 , 341 – 361 . Google Scholar Crossref Search ADS WorldCat Rice J.R. , 1988 . Elastic fracture mechanics concepts for interfacial cracks , J. Appl. Mech. , 55 , 98 – 103 . Google Scholar Crossref Search ADS WorldCat Rice J.R. , Lapusta N., Ranjith K., 2001 . Rate and state dependent friction and the stability of sliding between elastically deformable solids , J. Mech. Phys. Solids , 49 , 1865 – 1898 . Google Scholar Crossref Search ADS WorldCat Rojas O. , Dunham E.M., Day S.M., Dalguer L.A., Castillo J.E., 2009 . Finite difference modelling of rupture propagation with strong velocity-weakening friction , Geophys. J. Int. , 179 ( 3 ), 1831 – 1858 . Google Scholar Crossref Search ADS WorldCat Ru C.-Q. , 2010 . Simple geometrical explanation of Gurtin-Murdoch model of surface elasticity with clarification of its related versions , Sci. China , 53 , 536 – 544 . OpenURL Placeholder Text WorldCat Rubin A.M. , 2002 . Aftershocks of microearthquakes as probes of the mechanics of rupture , J. geophys. Res. , 107 , 2142 , doi:10.1029/2001JB000496 . Google Scholar Crossref Search ADS WorldCat Rubin A.M. , Ampuero J.-P., 2007 . Aftershock asymmetry on a bimaterial interface , J. geophys. Res. , 112 , B05307 , doi:10.1029/2006JB004337 . OpenURL Placeholder Text WorldCat Rubin A.M. , Gillard D., 2000 . Aftershock asymmetry/rupture directivity among central San Andreas fault microearthquakes , J. geophys. Res. , 105 , 19 095–19 109 . Google Scholar Crossref Search ADS WorldCat Shi Z. , Ben-Zion Y., 2006 . Dynamic rupture on a bimaterial interface governed by slip-weakening friction , Geophys. J. Int. , 165 , 469 – 484 . Google Scholar Crossref Search ADS WorldCat Simões F.M.F. , Martins J.A.C., 1998 . Instability and ill-posedness in some friction problems , Int. J. Eng. Sci. , 36 , 1265 – 1293 . Google Scholar Crossref Search ADS WorldCat Stroh A.N. , 1962 . Steady state problems in anisotropic elasticity , J. Math. Phys. , 41 , 77 – 103 . Google Scholar Crossref Search ADS WorldCat Suo Z. , 1990 . Singularities, interfaces and cracks in dissimilar anisotropic media , Proc. R. Soc. A , 427 , 331 – 358 . Google Scholar Crossref Search ADS WorldCat Uenishi K. , Rice J.R., 2003 . Universal nucleation length for slip-weakening rupture instability under nonuniform fault loading , J. geophys. Res. , 108 ( B1 ), 2042 , doi:10.1029/2001JB001681 . Google Scholar Crossref Search ADS WorldCat Weeks J.D. , 1993 . Constitutive laws for high-velocity frictional sliding and their influence on stress drop during unstable slip , J. geophys. Res. , 98 ( B10 ), 17637 , doi:10.1029/93JB00356 . Google Scholar Crossref Search ADS WorldCat Weertman J. , 1980 . Unstable slippage across a fault that separates elastic media of different elastic constants , J. geophys. Res. , 85 , 1445 – 1461 . OpenURL Placeholder Text WorldCat Williams M.L. , 1959 . The stresses around a fault or crack in dissimilar media , Bull. seism. Soc. Am. , 49 , 199 – 204 . OpenURL Placeholder Text WorldCat Xia K. , Rosakis A.J., Kanamori H., Rice J.R., 2005 . Laboratory earthquakes along inhomogeneous faults: directionality and supershear , Science , 308 , 681 – 684 . Google Scholar PubMed OpenURL Placeholder Text WorldCat Yang W. , Suo Z., Shih C.F., 1991 . Mechanics of dynamic debonding , Proc. R. Soc. A , 433 , 679 – 697 . Google Scholar Crossref Search ADS WorldCat SUPPORTING INFORMATION Supplementary data are available at GJIRAS online. Figure S1. Slip rate profiles at a fixed time step for decreasing δl: the figure shows the expected typical bimaterial asymmetry. The inset shows the convergence for small δl in terms of maximum amplitude and position of the crack front. Figure S2. Physical convergence criterion in time domain based on the L2-norm using cross-correlation is shown. We defined the distance between a signal s(t − τ) and a reference signal r(t) as: |$D = \int_{{t_1}}^{{t_2}} {{{| {s( {t - \tau } ) - r( t )} |}^2}{\rm{d}}t} $|⁠, where τ is the time delay obtained maximizing the cross-correlation function: |$c( t ) = \int_{{t_1}}^{{t_2}} {s( {t + \xi } )r( \xi )} \,{\rm{d}}\xi $|⁠. Normalizing D by the L2-norm of the reference signal yields the relative error e such that |${e^2} = \scriptstyle\frac{D}{{\int_{{t_1}}^{{t_2}} {{r^2}( t ){\rm{d}}t} }}$|⁠. For δl ≤ 10%Dc this error is lower than 0.05 as shown in panel (a). Panel (b) shows how the normalized time-shift decreases to zero when decreasing δl. Figure S3. Physical convergence for decreasing relaxation slip in time domain. σeff is shown at receivers 2 (a) and 8 (b). The induced perturbations are huger and sharper moving away from the nucleation and they are smoothed for increasing δl. Panels (c) and (d) show the amplitude spectra of the quantity in panels (a) and (b). Figure S4. (a) σeff as a function of time at receiver 3 when Cgr does not exist: the inset shows that convergence is still achieved when δl ≤ 10%Dc. The convergence is still driven by the introduced filtering and the panel (b) shows the amplitude spectra for decreasing δl. (c) Even for higher γ, δu* is fixed by dynamic timescale and it is lower than Dc when δl ≤ 10%Dc. Figure S5. Grid refinement in space domain for constant timescale, slip rate at the same time step for tc = 1.2 × 10− 4s (a) and (b) and tc = 6.0 × 10− 4s (c) and (d). (a) For smaller tc, when solutions are not convergent, strong oscillations of slip rate can emerge up to pathological effects (e.g. stop and go of rupture). Those effects can boost the rupture producing a spurious acceleration of the rupture front. The black square indicates the zoom around the crack front (b). (c) Even for higher tc, for which the oscillatory effects are damped, the solutions for coarser meshes do not converge with those obtained from finer grids. When solutions converge position of crack front and amplitude of the maximum coincide. The range of mesh convergence is the same found for dynamic timescale. The black square indicates the zoom around the crack front (d). Figure S6. Acceleration of rupture when Cgr does not exist. In panel (a) three different ratios γ are used, whereas the ratio ρ1/ρ2 is fixed at 1.2. In panel (b) four different ratios ρ2/ρ1 are used, whereas the ratio γis fixed at 1.9. Figure S7. Grid refinement in space domain for dissipation length scale, slip rate for two different parameters β when a dissipation length scale is used: the features for all models are pretty the same obtained for dynamic and constant timescales (panels a and c). The black square indicates the zoom around the crack front (panels b and d). Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. APPENDIX A: SPECTRAL ELEMENT METHOD FOR RUPTURE DYNAMICS The spectral element method solves the general problem of the elastodynamics (eq. 1) using a weak formulation. This is obtained by multiplication of the eq. (1) by an admissible displacement w which is zero at the boundary of the physical domain Ω and by integration over Ω: \begin{equation}\int_\Omega {\rho {{\bf w\dot{v}}}\,{\rm{d}}\Omega } = - \int_\Omega {\nabla {{\bf w}}:{{\boldsymbol \sigma }}} + \int_{{\Gamma _N}} {{{\bf wT}}\,{\rm{d}}\Gamma } .\end{equation} (A1) The latter term comes from integration by parts of the divergence of the stress tensor and ΓN is the portion of the boundary over which Neumann conditions are applied. In this work we assume ΓN = Γe∪Γ, where Γe is the Earth surface and Γ represents the fault. The free surface condition is naturally accounted for, considering zero the contribution of the last integral in eq. (A1) coming from Γe. To account for the faulting process, we decompose the physical domain Ω into two regions Ω1 and Ω2, such that Ω1∪Ω2 = Ω, Ω1∩Ω2 = Γ∪Γf, the measure of Ω1∩Ω2 is zero, and Γf is an arbitrary prolongation of Γ, needed to reach the external boundary of the physical domain Ω (Fig. 1). Across Γf continuity of kinematic and dynamic quantities is ensured. We separately solve eq. (A1) in the two domains, where the interface ΓN = Γ∪Γf is common to Ω1 and Ω2. The domains are discretized in quadrangular or hexahedral elements for 2-D or 3-D regions, the displacement and velocity are then approximated within each element by Lagrange polynomials defined on Gauss–Lobatto–Legendre (GLL) collocation points and numerical integration is based on the GLL quadrature formula. Finally the problem is reduced to the following algebraic system (Komatitsch & Vilotte 1998): \begin{equation}{{{\bf M}}_i}{{{\bf \dot{v}}}_i} = {{\bf F}}_i^{{\mathop{\rm int}} } + {{\bf B}}_i^T{{{\bf T}}_i}.\end{equation} (A2) Here M is the mass matrix, which inherits diagonality from numerical orthogonality of the Lagrange polynomials, |${{{\bf F}}^{{\mathop{\rm int}} }}$| is the vector of internal forces, which depends on the local displacement, and BT is interface matrix related to boundary Γ∪Γf. The details of the approximations introduced by the spectral element method and yielding eq. (A2) can be found in Komatitsch & Vilotte (1998) and Chaljub et al. (2007). We also assemble the contributions of the mass matrix and the internal forces on the fictitious boundary Γf, reducing the contribution of the interface matrix to the fault Γ. The eq. (A2) can be discretized in time using a displacement-velocity Newmark scheme. We adopt the formulation suggested in Festa & Vilotte (2005), where a dissipative scheme is used to damp the high frequency spurious oscillations generated by the propagating rupture. Time discretization at time (p + 1)Δt, where Δt is the time step, yields: \begin{equation}{{{\bf M}}_i}{{\bf v}}_i^{(p + 1)} = {{{\bf M}}_i}{{\bf v}}_i^{(p)} + \Delta t\left( {{{\bf F}}_i^{{\mathop{\rm int}} \,(p + 1)} + {{\bf B}}_i^T{{\bf T}}_i^{(p + 1)}} \right).\end{equation} (A3) After inversion of the mass matrix, subtracting eq. (A3) for the medium 2 by the respective equation for the medium 1 leads to \begin{equation}\delta {{\bf v}}_{}^{(p + 1)}\, =\, \delta {V^{(p + 1)}} - {{{\bf C}}_\Gamma }{{{\bf T}}^{(p + 1)}},\end{equation} (A4) where we used the position T = T1 = −T2, δv(p) is the slip rate at any GLL point on the fault at time pΔt and |$\delta {V^{(p + 1)}} = \delta {{{\bf v}}^{(p)}} + \Delta t\,( {{{\bf \tilde{M}}}_2^{ - 1}{{\bf \tilde{F}}}_2^{{\mathop{\rm int}} \,(p + 1)} - {{\bf \tilde{M}}}_1^{ - 1}{{\bf \tilde{F}}}_1^{{\mathop{\rm int}} \,(p + 1)}} )$|⁠, where the tilde operator represents the restriction of mass matrices and internal forces to the fault. The operator |${{{\bf C}}_\Gamma }$| is defined as \begin{equation}{{{\bf C}}_\Gamma } = \Delta t\left( {{{{{\bf \tilde{M}}}}^{ - 1}}_2\, + {{{{\bf \tilde{M}}}}_1}^{ - 1}} \right){{{\bf B}}^T}.\end{equation} (A5) An analogous relationship can be obtained for the displacement discontinuity (Festa 2004) \begin{equation}\delta {{\bf u}}_{}^{(p + 1)}\,=\, \delta {U^{(p + 1)}} - {{{\bf Q}}_\Gamma }{{{\bf T}}^{(p + 1)}},\end{equation} (A6) where |$\delta {U^{(p + 1)}} = \delta {{{\bf u}}^{( p )}} + \Delta t\delta {{{\bf v}}^{( p )}} + \frac{1}{2}\Delta {t^2}( {{\bf \tilde{M}}}_2^{ - 1}{{\bf \tilde{F}}}_2^{{\mathop{\rm int}} \,(p + 1)} -{{\bf \tilde{M}}}_1^{ - 1}$||${{\bf \tilde{F}}}_1^{{\mathop{\rm int}} \,(p + 1)})$| and |${{{\bf Q}}_\Gamma } = \frac{1}{2}\Delta t{{{\bf C}}_\Gamma }$|⁠. The diagonality of the matrices |${{{\bf Q}}_\Gamma }$| and |${{{\bf C}}_\Gamma }$| 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 |$\sigma _T^{n\,\,(p)} = \sigma _0^{n\,\,} + {T^{n\,\,(p)}}$| and |$\tau _T^{\,\,(p)} = \tau _0^{\,\,} + {T^{t\,\,(p)}}$|⁠, where |$\sigma _0^n$| and τ0 represent the pre-stress field components, and Tn and Tt 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: \begin{equation}\bar{\sigma }_T^{n,K(p + 1)} = {\left( {Q_\Gamma ^{KK}} \right)^{ - 1}}\,\delta {U_{}}^{n,K\,\,(p + 1)} + \sigma _0^{n,K}.\end{equation} (A7) If |$\bar{\sigma }_T^{n,K\,\,(p + 1)} \le 0$|⁠, the Signorini condition is automatically satisfied and |$( {0,\bar{\sigma }_T^{n,K\,\,(p + 1)}} )$| is the solution for the couple |$( {\delta u_{}^{n,K\,\,(p + 1)},\sigma _T^{n,K\,\,(p + 1)}} )$|⁠. In the case for which eq. (A7) yields |$\bar{\sigma }_T^{n,K\,\,(p + 1)} > 0$|⁠, the normal traction is forced to be zero (⁠|$\sigma _T^{n,K\,\,(p + 1)} = 0$|⁠) and |$\delta u_{}^{n,K\,\,(p + 1)} = Q_\Gamma ^{KK}\,\,\bar{\sigma }_T^{n,K\,\,(p + 1)}$|⁠. From the Signorini condition, when opening occurs (⁠|$\sigma _T^{n,K\,\,(p + 1)} = 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): \begin{equation}\sigma _{{\rm{eff}}}^{(p + 1)} = {\left( {1 + \frac{{\Delta t}}{{{t^*}}}} \right)^{ - 1}}\left[ {\sigma _{{\rm{eff}}}^{(p)} + \frac{{\Delta t}}{{{t^*}}}\sigma _{}^{n\,\,(p + 1)}\,} \right].\end{equation} (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 δvt  (p + 1) = 0. Accordingly, the corresponding tangential traction from (A4) is \begin{equation}\bar{\tau }_T^{K\,\,(p + 1)} = \tau _0^K + {\left( {C_\Gamma ^{KK}} \right)^{ - 1}}\,\delta {V_{}}^{t,K\,\,(p + 1)}.\end{equation} (A9) If this quantity is below the threshold |$| {\bar{\tau }_T^{K\,\,(p + 1)}} | \le - \mu ( {\delta u_{}^{t,K\,\,(p + 1)}} )\,\sigma _{eff,T}^{K\,\,(p + 1)}$|⁠, the couple |$( {0,\bar{\tau }_T^{K\,\,(p + 1)}} )$| is the solution for |$( {\delta {v^{t,K\,\,(p + 1)}},\tau _T^{K\,\,(p + 1)}} )$| and the point K is at rest. Otherwise, the traction is at the threshold |$( {| {\bar{\tau }_T^{K\,\,(p + 1)}} | = - \mu ( {\delta u_{}^{t,K\,\,(p + 1)}} )\sigma _{{\rm{eff}},T}^{K\,\,(p + 1)}} )$| and the corresponding solution for tangential slip rate is: \begin{equation}\delta {v^{t,K\,\,(p + 1)}} = \delta {V_{}}^{t,K\,\,(p + 1)} + C_\Gamma ^{KK}\mu \left( {\delta u_{}^{t,K\,\,(p + 1)}} \right)\sigma _{{\rm{eff}},T}^{K\,\,(p + 1)}.\end{equation} (A10) Eqs (A8) and (A10) have to be simultaneously solved. When using a regularization based on the constant timescale, eq. (A8) can be solved separately from the Coulomb problem. In this case, we sequentially solve the normal traction, the relaxation scheme and the Coulomb projection. When t* depends on the actual value of the tangential slip rate, as for the dynamic regularization, the two equations are coupled and they are iteratively solved. Author notes * " Now at: Istituto Nazionale di Geofisica e Vulcanologia, Rome, Italy. © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
TI - Rupture dynamics along bimaterial interfaces: a parametric study of the shear-normal traction coupling
JF - Geophysical Journal International
DO - 10.1093/gji/ggw489
DA - 2017-04-01
UR - https://www.deepdyve.com/lp/oxford-university-press/rupture-dynamics-along-bimaterial-interfaces-a-parametric-study-of-the-bP9fYheyzG
SP - 48
VL - 209
IS - 1
DP - DeepDyve
ER -