Cracks, pulses and macroscopic asymmetry of dynamic rupture on a bimaterial interface with velocity-weakening friction

We study in-plane ruptures on a bimaterial fault governed by a velocity-weakening friction with a regularized normal stress response. Numerical simulations and analytical estimates provide characterization of the ranges of velocity-weakening scales, nucleation lengths and background stresses for which ruptures behave as cracks or pulses, decaying or sustained, bilateral or unilateral. With strongly velocity-weakening friction, ruptures occur under a wide range of conditions as large-scale pulses with a preferred propagation direction, that of slip of the more compliant material. Such ruptures have macroscopic asymmetry manifested by signiﬁcantly larger seismic potency and propagation distance in the preferred direction, and clearly quantiﬁed by the directivity ratio derived from the second order moments of the spatiotemporal distribution of slip rate. The macroscopic rupture asymmetry of the large-scale pulses stems from the difference in the criticality conditions for self-sustained propagation in each rupture direction, induced by the asymmetric normal stress changes operating in bimaterial interfaces. In contrast, crack-like ruptures show macroscopic asymmetry under restrictive conditions. The discussed mechanism is robust with respect to regularization parameters, ranges of stress heterogeneities and a proxy for off-fault yielding and should operate similarly for crustal-scale rupture pulses even in the absence of velocity-weakening. Small-scale pulses, driven by the bimaterial normal stress reduction at the scale of the process zone, can detach from the rupture front of the large-scale pulses that propagate in the preferred direction. However, their occurrence depends on the relaxation scale in the regularization of the normal stress response and their development can be hindered by off-fault yielding.


I N T RO D U C T I O N
The juxtaposition of materials with dissimilar elastic properties is a general feature of tectonic faults in a variety of contexts and scales. Subduction interfaces bring into contact oceanic and continental crusts, normal and reverse faults offset vertical stratifications, large strike-slip faults displace different crustal blocks and oceanic transforms juxtapose rocks of different ages. At the scale of the inner architecture of a fault zone, bimaterial interfaces associated with rock damage are present with various degrees of sharpness. Failure along a bimaterial interface can also be relevant at the microscopic scale where grain boundaries are potential fracture surfaces and in the context of glacial quakes with slip at the interface between ice and rock (Deichmann et al. 2000;Weertman 2005;Danesi et al. 2007). Moreover, bimaterial rupture is a prevalent mode of failure in composite engineering materials: impact-induced delamination of composite laminates, decohesion of adhesive joints in bonded structures, particle debonding in reinforced elastomers, fibre/matrix debonding and pull-out in fibre-reinforced materials, delamination of thin-film/substrate interfaces, etc.
In-plane ruptures on planar bimaterial interfaces have remarkable dynamic properties that may be relevant to many issues of basic and applied science (e.g. Ben-Zion 2001). In contrast to slip between similar media, slip along a planar bimaterial interface generates dynamic changes of normal stress σ that depend on the slip gradient, material properties, rupture velocity and direction of rupture propagation (Weertman 1980). For subshear ruptures, the change of σ at the tip propagating in the direction of slip of the compliant solid, referred to as the 'preferred direction', is tensile and produces a dynamic reduction of the frictional strength. In contrast, the change of σ at the tip propagating in the opposite direction is compressive, leading to dynamic strengthening. For supershear ruptures the senses of changes of σ are reversed (Weertman 2002;Shi & Ben-Zion 2006). The amplitudes of the near-tip changes of σ increase with propagation distance along the bimaterial interface, due to a continual transfer of energy to shorter wavelengths (Adams 1995; Ben-Zion & Huang 2002). In the following we refer to dynamic change of σ due to in-plane slip on a planar interface as the 'bimaterial effect'.
The theoretical analysis by Weertman (1980) suggested the existence of a peculiar rupture mode: even in the absence of frictional weakening, a slip pulse can propagate along a bimaterial interface in the preferred direction without significant energy dissipation by friction. Andrews & Ben-Zion (1997) confirmed this prediction through numerical simulations of self-sustained pulses, triggered by transient pore pressure release, and referred to such a rupture mode as 'wrinkle-like pulse'. Subsequent studies with constant coefficient of friction attempted to clarify various aspects of the wrinklelike pulse, including the conditions for its existence (Ben-Zion & Andrews 1998), conditions leading to ill-posedness and regularization (Adams 1995;Cochard & Rice 2000;Adda-Bedia & Ben Amar 2003), effects of low velocity fault zone layer and multiple possible rupture surfaces (Ben-Zion & Huang 2002;Brietzke & Ben-Zion 2006) and effects of off-fault plastic yielding (Ben-Zion & Shi 2005). The characteristics of the wrinklelike pulse suggest that earthquakes on bimaterial faults tend to propagate in a preferred direction determined by the material contrast. Several lines of evidence were proposed to support this view, including asymmetric pattern of immediate small aftershocks (Rubin & Gillard 2000), asymmetric damaged fault zone rocks observed in the field (Dor et al. 2006a(Dor et al. ,b, 2008 and asymmetric fault zone damage inferred from seismic trapped and head waves (Lewis et al. 2005(Lewis et al. , 2007. Numerical simulations of ruptures on a bimaterial interface governed by slip-weakening friction produced results that depend on the employed nucleation procedure. Simulations with localized strong nucleation of the type used by Andrews & Ben-Zion (1997) generated ruptures that evolved to unidirectional wrinkle-like pulses similar to those generated with slip-independent friction (Shi & Ben-Zion 2006). In contrast, relatively large smooth nucleation process that lead to runaway ruptures in a homogeneous solid generated bilateral cracks with a superposed wrinkle-like pulse at the preferred direction (Harris & Day 1997, 2005Rubin & Ampuero 2007). Laboratory experiments also showed a diversity of rupture styles along bimaterial interfaces (Anooshehpoor & Brune 1999;Xia et al. 2005) that probably reflect differences in material and frictional properties, nucleation conditions and rupture size.
Some aspects of these apparently contrasting results have been the subject of recent debate (Andrews & Harris 2005;Harris & Day 2005;Ben-Zion 2006a,b). Nevertheless, it is clear that the bimaterial effect can lead to a wrinkle-like pulse, independent of the presence of frictional weakening and of the nucleation procedure employed. It is also clear that the wrinkle-like pulse can have important effects on earthquake ruptures, whether it propagates on its own or remains confined to the crack tip region of a bilateral crack. Rubin & Ampuero (2007) showed that the superposed wrinkle-like pulse spawned upon interaction with barriers from the crack tip that propagates in the preferred direction is consistent with the body of observations cited above, without requiring predominantly unidirectional rupture.
As pointed out in the discussion sections of Andrews & Ben-Zion (1997), this paper and other works, unilateral rupture is not expected to be realized stricto sensu in natural bimaterial faults due to heterogeneities and various perturbations from the basic wrinklelike pulse model. Instead, a given bimaterial fault may generate events with different degrees of observable 'macroscopic source asymmetry', which may be quantified for instance by the spatial skewness of the final slip (Shi & Ben-Zion 2006). The focus of this paper is on a physical mechanism that can plausibly produce a statistical propensity for 'macroscopically asymmetric ruptures', rather than on whether all bimaterial ruptures are strictly unilateral. An important related question is what signatures of the bimaterial effect are accessible for seismological and laboratory observations. Our results indicate that the bimaterial effect manifests itself not only through the emergence of small-scale wrinkle-like pulses, but more significantly through the asymmetric destabilization of largescale rupture pulses.
The simulations with slip-weakening friction by Harris & Day (1997, 2005 and Rubin & Ampuero (2007), although macroscopically bilateral, produced very prominent asymmetry of slip velocities at the opposite rupture fronts. Under friction that depends only on slip, this strong asymmetry of slip velocities can not manifest itself into macroscopic asymmetry of the final slip. However, laboratory experiments (Tsutsumi & Shimamoto 1997;Di Toro et al. 2004;Mizoguchi et al. 2006;Yuan & Prakash 2007) and theoretical considerations of thermal and fluid weakening mechanisms (Sibson 1973;Noda & Shimamoto 2005;Rice 2006) imply that the friction coefficient reduces dramatically with high slip rates relevant for seismic rupture. Ben-Zion (2006a) suggested that velocity-dependent friction may produce a positive feedback between frictional weakening and the bimaterial effect, leading to larger stress drop in the preferred direction and macroscopically asymmetric rupture. In this paper we examine in detail this possibility, through numerical simulations of 2-D in-plane ruptures on a bimaterial interface. The calculations employ a spectral boundary integral method and a rateand state-dependent friction law with strong velocity dependence. The parameter-space study covers effects associated with velocityweakening, nucleation zone size and heterogeneous initial stress.
The remainder of the paper is organized as follows. In Section 2, we introduce the assumptions of the model, the governing equations and the numerical method to solve them. The assumed friction law contains slip-weakening or velocity-weakening as limit cases, depending on the relaxation scale in the state evolution law. The steady-state friction coefficient is inversely proportional to sliprate, to mimic the weakening mechanisms thought to operate on natural faults at high slip velocities. In Section 3, we examine the behaviour of ruptures triggered by a slightly overstressed nucleation zone. We first study the conditions that lead to crack-like and pulselike ruptures in a homogeneous solid with uniform initial stress. We then show that pulse-like ruptures on a bimaterial interface develop strong macroscopic asymmetry, manifested by significantly larger seismic potency and propagation distance in the preferred direction. We assess the statistical robustness of this rupture asymmetry with respect to stress heterogeneities and a proxy for off-fault yielding, and quantify it by a seismologically observable measure based on second order spatio-temporal moments of slip rate. In Section 4, we discuss the implications of our results in relation to the relevant observational evidence and to additional physical processes. Several Appendices provide analytical results that are used to verify the numerical calculations and interpret the results.
to operate on natural faults at high slip rates, while keeping a simple mathematical formulation. We adopt a velocity-and state-dependent friction law equivalent to that employed by Cochard & Madariaga (1996) and Shaw & Rice (2000). The friction coefficient is given by and is defined by the following parameters: a static friction coefficient μ s , a characteristic velocity scale V c and two positive coefficients, α and β, quantifying a direct effect and an evolution effect, respectively. The second term in eq. (1) is a direct velocity strengthening mechanism with a regularizing effect. The third term is a weakening mechanism, with a state variable θ governed bẏ where τ c is a characteristic timescale over which the friction coefficient evolves to its steady-state value If α < β, the steady-state friction coefficient weakens asymptotically as 1/V at slip velocities much faster than V c , approaching its The relaxation time τ c tunes the weakening term between two limit cases: slip-weakening and velocity-weakening (Shaw & Rice 2000). If τ c is much longer than the typical timescale of fluctuation of the state variable, eq. (2) becomesθ ≈ V /τ c , implying that θ is proportional to slip and that the evolution term of the friction coefficient is effectively slip-weakening, with characteristic slipweakening distance Conversely, if τ c is short, the relaxation to steady state is fast and the friction is effectively velocity-weakening, with characteristic velocity scale V c . In our simulations we keep D c constant and vary V c , so the transition from slip-weakening to velocity-weakening is obtained when V c is increased. The remaining friction parameters are fixed as in Table 1. The frictional and elastic parameters determine a slip-velocity scale V dyn that is relevant for the nucleation of instabilities and required numerical discretization (Appendix A). With the parameters in Table 1, V dyn = 0.23 m s −1 . For the values of V c studied here, the critical wavelength (defined in Appendix A) ranges from 3.5 to 5.2 m, whereas for pure slip-weakening friction it would be 1.9 m.

Normal stress response
There is currently no consensus or definitive experimental constraints on the normal stress response of natural faults at high slip rates (Ben-Zion 2001; Rice et al. 2001;Rubin & Ampuero 2007). It is, however, recognized that a characteristic time or slip scale of normal stress response is required to regularize the possible illposedness of the bimaterial rupture problem (Cochard & Rice 2000;. Following Rubin & Ampuero (2007), the fault strength is assumed to be proportional to an effective normal stress σ * which follows a regularized evolution law: where V * and δ σ are constitutive parameters, with values given in Table 1. To emphasize the role of large-scale pulses generated by the velocity-weakening friction, we deliberately set δ σ ≈ D c . This choice minimizes the bimaterial effect-it delays the onset of small-scale wrinkle-like pulses spawned by large-scale rupture fronts (Rubin & Ampuero 2007). To facilitate proper grid resolution, in a number of simulations we assumed insteaḋ

Background stress and heterogeneities
Natural faults are believed to contain heterogeneities on a broad range of scales, although their properties and nature are poorly constrained at short scales. We assume a uniform initial normal stress σ 0 and explore the effect of heterogeneous initial shear stresses τ 0 (x). Faults with velocity-weakening friction can generate complex seismicity patterns even under uniform friction properties (Cochard & Madariaga 1996), so we do not find it essential at present to consider heterogeneities of friction parameters. The nominal static and dynamic strength are denoted τ s = σ 0 μ s and τ d = σ 0 μ d , respectively. We define the background shear stress τ 0 as the spatial average of τ 0 (x). These quantities are used to define the strength excess parameter The value of τ 0 is chosen low enough to ensure that rupture speed remains sub-Rayleigh; we set S = 1.7, which in homogeneous media leads to supershear transition at hypocentral distances well beyond the limits of our modelled fault (Andrews 1976).
Following Ampuero et al. (2006) and Ripperger et al. (2007), plausible τ 0 (x) distributions are generated as stochastic fields with the following power spectrum shape where k is the spatial wavenumber along the fault, std is the standard deviation normalized by the nominal strength drop, L 1 is the von Karman correlation length, L 2 is a Gaussian smoothing length and H is the Hurst exponent. We consider four different sets of spectral parameters, A-D, defined in Table 2. Examples are shown for each set in Fig. 1. For the purposes of comparison, set A is viewed in two different ways. Sets A-C are pink noise (H = 0, power spectrum ∝ 1/|k|) with increasing correlation length. Sets A and D are white noise (H = −1/2, flat power spectrum), but small-scale correlation is still present through L 2 . These cases are representative of the stresses (H = 0) that are Table 2. Heterogeneous initial stress parameter sets: von Karman correlation length L 1 , Gaussian smoothing length L 2 and Hurst exponent H. The symbol × means the parameter is irrelevant. The two equivalent views of set A are related to comparison groupings ABC and AD, respectively.  Table 2, with normalized standard deviation std = 0.25 and background stress τ 0 = 82.81 MPa. The higher and lower dashed lines indicate the nominal yield and dynamic strength, τ s and τ d , respectively. theoretically expected for self-similar seismicity with constant average stress drop (Andrews 1980), and of stresses (H ≈ 0 or H < 0) derived from kinematic source slip inversions (Tsai 1997;Mai & Beroza 2002). The normalized standard deviation, std, is varied from case to case.

Nucleation procedure
The nucleation conditions are known to affect the rupture style, and in particular, the selection between crack-like and pulse-like ruptures (Nielsen & Madariaga 2003;Festa & Vilotte 2006;Shi et al. 2008). To minimize the number of arbitrary assumptions, while keeping a reasonable computational time, we adopt the following nucleation procedure: we set on an initial patch of size L nuc the initial stress slightly above the nominal static strength, with a very small triggering overstress factor = 10 −5 . We explore a range of values of L nuc , somewhat larger than the theoretical quasistatic critical wavelength (see Appendix A).

Numerical method
Spontaneous dynamic ruptures are calculated by the spectral boundary integral equation method (SBIEM) developed by Cochard & Rice (2000) and parallelized by Rubin & Ampuero (2007). The Fortran/MPI code, now modified to include velocity-weakening friction, is freely available through the Orfeus Seismological Software Library (http://www.orfeus-eu.org/links/softwarelib.htm). The spatial discretization, x = 0.15 m, is taken smaller than the characteristic length with maximum growth rate in a linear stability analysis (Appendix A), and smaller than an estimate of the process zone size for bimaterial crack-like ruptures (Rubin & Ampuero 2007).

Rupture behaviour in homogeneous medium with uniform initial stress
To provide a context for the calculations on bimaterial faults, we first discuss results for ruptures in a homogeneous solid with uniform initial stresses. Dynamic ruptures on velocity-weakening faults can propagate as pulses or as cracks and can be either decaying or sustained (Cochard & Madariaga 1994Perrin et al. 1995;Beeler & Tullis 1996;Zheng & Rice 1998;Nielsen & Carlson 2000;Nielsen & Madariaga 2003). Fig. 2 shows several examples of simulated ruptures, and Fig. 3 summarizes the rupture style of 315 numerical simulations for ranges of values of the characteristic velocity V c and the size of the nucleation patch L nuc (the other parameters being fixed as in Table 1). The propagation distances generated at a given time after rupture initiation by the different sets of parameters are summarized in Fig. 4. It is known that pulse-like rupture can appear in faults governed by velocity-weakening friction and that, with the nucleation procedure adopted here, ruptures under pure slip-weakening friction propagate exclusively as cracks. Among the sustained ruptures in Fig. 3, we distinguish pulses from cracks by the presence or absence of healing in the hypocentre. As expected from the transition of frictional behaviour as a function of V c , sustained ruptures are crack-like at low V c (predominantly slip-weakening) and pulse-like at high V c (predominantly velocity-weakening). The crack/pulse transition is explored in more detail in     value of V c , the healing time becomes infinite (Fig. 5) signifying a transition to crack-like rupture. The critical V c is higher if the background stress τ 0 is raised or, to a minor extent, if L nuc is made larger (Fig. 6). The dependence on τ 0 is in quantitative agreement with the theoretical analysis of Nielsen & Carlson (2000), summarized in Appendix D.
At a fixed V c value, a transition from decaying to sustained ruptures occurs when L nuc is increased: beyond a critical value of L nuc the rupture pulses become unstoppable (Figs 3, 4 and 7). This may be understood by considering the speeds of the rupture and healing fronts. The variation among cases in the speed of the rupture front is mainly controlled by the energy available from the overstressed nucleation region, which grows with increasing size L nuc . The rupture front initially follows the behaviour expected for a dynamic crack: initial acceleration due to the nucleation impulse, followed by deceleration when leaving the nucleation area due to the lower available stress drop, and re-acceleration due to the increasing energy release rate roughly proportional to the crack length (Fig. 7). The rupture front speed reaches a peak soon after the onset of pulse-like rupture and then decelerates again because from then on, the energy release scales with the pulse width, which is a smaller length scale than the earlier crack size. The healing front evolves quite independent of L nuc at distances far from the nucleation region, slowing down to a stable speed of the order of 85 per cent of the Rayleigh speed in the examples of Fig. 7. The transition from decaying to sustained pulse is presumably controlled by the ratio between the healing front speed and peak rupture front speed. Gentle nucleation (small L nuc ) produces slower rupture than healing front speed, the pulse width shrinks, the energy release decreases and the pulse collapses. Conversely, vigorous nucleation (large L nuc ) pushes the rupture speed closer to or beyond the asymptotic healing front speed, leading to a sustained pulse.
At higher V c values, the transition from decaying to sustained rupture requires a larger energy supply from the nucleation process (larger L nuc ). As shown in Appendix C, both the effective fracture energy and the dynamic stress drop, two competing quantities that control rupture propagation, decrease with increasing V c (Fig. C1) so it is not evident why ruptures are easier to stop at larger V c (Fig. 8). Further analysis of the energy balance, as developed in Appendix C, shows that the effect of the reduction of dynamic stress drop eventually overcomes the reduction of fracture energy as a function of V c , favouring rupture arrest. Moreover, we observe that healing fronts tend to be faster at larger V c , which also favours rupture arrest.
If the background stress τ 0 is increased, ruptures become harder to stop due to the larger available energy and the transition to runaway ruptures occurs at smaller L nuc .   Figure 9. Rupture styles on a velocity-weakening bimaterial fault, with same parameters as in Fig. 2. Pulse-like ruptures (left-hand and middle panels) are strongly asymmetric. By the last snapshot, pulse-like (middle panel) and crack-like (right-hand panel) sustained ruptures spawn a small-scale wrinkle-like pulse from the rupture front that propagates in the preferred direction (positive).

Rupture behaviour in a bimaterial interface with uniform initial stress
We ran a similar set of 315 simulations on a bimaterial fault with uniform initial stress, 18 per cent contrast in S-wave speed and 40 per cent contrast in shear moduli. The assumed properties for the two elastic media are summarized in Table 3. With these settings, there is no significant density contrast, the effective shear impedance, defined through the harmonic average of the shear impedances of the two media (c s1 /μ 1 + c s2 /μ 2 ) −1 , is practically the same as in the previous section, and generalized Rayleigh waves do exist (Weertman 1980;. All the remaining parameters are fixed as in the previous section. The main rupture styles on the bimaterial fault are illustrated in Fig. 9 and summarized in Fig. 10. The styles previously described for a fault in homogeneous media are also present here. However, the bimaterial faults tend to produce highly asymmetric pulse-like ruptures, which develop larger slip in the 'preferred' rupture direction and decay in the opposite direction. Whereas in these simulations the crack/pulse transition is not significantly affected by the material contrast, the decaying/sustained rupture transition is strongly affected. The transition from cracks to pulses depends on the balance between radiation damping and the steady-state friction response (Appendix D). Because the effective shear impedance is chosen to be essentially the same as in the homogeneous case, and the normal stress changes in the hypocentre are small due to the vanishing slip gradient there, the value of V c at the transition from cracks to pulses remains similar to that of the homogeneous case.
In contrast, the transition from decaying to sustained ruptures splits into two branches, one for each rupture direction (compare Fig. 10 to Fig. 3). At a given value of V c , the triggering sizes L nuc required for pulses to run away in the 'preferred' direction are smaller than in the homogeneous case. Within these pulses, the bimaterial effect reduces the normal stress, which increases the dynamic stress drop and makes these rupture fronts harder to stop. The fracture energy also increases, proportionally to the increasing strength drop, providing a competing effect. However, the increase of static energy release rate, which is proportional to the square of the increasing dynamic stress drop, appears to dominate in the energy balance. Overall, the bimaterial effect enhances rupture acceleration in the 'preferred' direction and the transition to runaway ruptures is very sensitive to this effect. All the ruptures promoted to self-sustained (runaway) events on bimaterial faults are asymmetric pulse-like ruptures propagating in the 'preferred' direction. Conversely, pulses running in the opposite direction are instead stabilized by the increase of normal stress induced therein by the bimaterial effect. To become unstoppable, they require stronger energy supply from the nucleation (larger L nuc ) than in the homogeneous case. This has also been observed in pulse-like ruptures under slip-weakening friction triggered by strong localized nucleation conditions (Shi & Ben-Zion 2006).
A crude estimate of the dynamic stress drop is derived in Appendix E, from known relations for steady-state bimaterial pulses (Weertman 1980), and its implications on the bimaterial effect described above are discussed. In particular, the amplification of the dynamic stress drop is stronger when the rupture approaches the limiting speed, which is the generalized Rayleigh wave speed or the slowest shear wave speed if the former does not exist. This mechanism is excited sooner, at lower rupture speeds, for moderate material contrasts, c s2 /c s1 ≈ 30 per cent. This is consistent with the different strengths of the bimaterial effect, and a maximum around 35 per cent contrast, found in the early simulations of Ben-Zion & Andrews (1998). Bimaterial ruptures show varying degrees of asymmetry in rupture velocity, slip rate, propagation distance and seismic potency (Fig. 11). The most prominent asymmetry is found in pulse-like ruptures that are sustained in the preferred direction but decay in the opposite direction. In contrast, the macroscopic asymmetry in the crack-like ruptures is weak. Strong asymmetry might appear in decaying cracks that propagate over long distances, but these can only occur under restrictive conditions, requiring particularly small V c and low background stress (large ratios of final size to nucleation size predicted in Fig. B2 require large S values).
As observed by Rubin & Ampuero (2007) under slip-weakening friction and gentle nucleation conditions, the predominantly bilateral cracks found here at low V c spawn a small-scale wrinkle-like pulse from the crack front propagating in the 'preferred' direction. This feature is found here also for the pulse-like ruptures, as illustrated in Fig. 12: a 'small-scale pulse' is nucleated at the leading edge of the 'large-scale pulse' that runs in the 'preferred' direction. The mechanism is controlled by the increasing normal stress changes induced by the steepening of the slip gradient within the rupture front process zone, due to the Lorentz contraction (Rubin & Ampuero 2007). As the asymptotic behaviour near the rupture front and its dependence on rupture speed is the same for pulses and cracks, it is not surprising that this mechanism is also promoted at the edge of pulse-like ruptures.

Ruptures on a bimaterial fault with heterogeneous initial shear stress
We selected from the previous simulations, conditions under which rupture is pulse-like and asymmetric: L nuc = 17 m and V c = 0.19 m s −1 . In order to assess statistically the robustness of the asymmetry, we added a heterogeneous component to the initial stress τ 0 (x), as described in Section 2.3. For each set in Table 2, thirty realizations of the randomized initial stress profile were generated. Each realization was rescaled to explore three different values of the standard deviation, std = 0.05, 0.15 and 0.25. The initial stress was shifted in space to locate its maximum at the centre of the fault, and an overstressed nucleation patch was set as described in Section 2.4. To avoid multiple nucleation spots, the occasional stress values above the static yield strength outside the central nucleation patch were reset below yield by a small fraction (2 per cent) of the nominal strength drop. For each of these τ 0 (x) we performed a  second simulation with a spatially flipped version of the initial stress, τ 0 (− x). All the simulations include a limit slip-rate V tough = 10 m s −1 as a proxy for off-fault yielding (Andrews 2005), and stopping barriers (σ 0 = 300 MPa) located at 200 m from the hypocentre. Examples corresponding to a single realization of stress heterogeneity from set D, with different levels of std, are shown in Fig. 13 together with their spatially flipped versions. In cases (b-d) the macroscopic asymmetry in the 'preferred' direction persists, independently of the amplitude of the heterogeneities. In cases (f-h) the asymmetry attenuates as the heterogeneities increase. For case (h) a subevent is triggered in the opposite direction. The fate of each rupture, especially at high std, is strongly dependent on the particular details of the initial stress distribution. This is expected for pulselike ruptures, because they are controlled by the fluctuations at the local scale of the pulse width, whereas it should be less marked for crack-like ruptures.   Table 2). There are roughly five times more ruptures with strong asymmetry towards the 'preferred' direction than ruptures with weak or opposite asymmetry. Because we considered also symmetric versions of each initial stress realization, the corresponding histograms for a fault in homogeneous medium would be exactly symmetric. The ratio of asymmetry in our bimaterial simulations, taken as a whole, is consistent with the 3:1 asymmetry in the occurrence of immediate aftershocks observed by Rubin & Gillard (2000). However, this agreement may be fortuitous, as the potency ratio in our models depends on other model parameters. Fig. 15 provides histograms of potency ratios for different levels of stress heterogeneity (std) in the set D. For weak heterogeneity with std = 0.05, the behaviour is strongly asymmetric, as in the reference uniform case. For larger std, the troughs of the initial stress field are large enough to significantly affect the rupture propagation, and some ruptures manage to run almost unilaterally in the opposite direction. Statistically, the macroscopic asymmetry in the 'preferred' direction becomes less pronounced with increasing std, but it is still significant (ratio 2:1) at std = 0.25. Table 4 summarizes the overall asymmetry of seismic potency as a function of std for all sets A-D, quantified by the ratio of the number of events with larger seismic potency in the 'preferred' direction to the number of events with larger seismic potency in the opposite direction. Increasing the correlation length and the amplitude (std) of heterogeneities attenuates the macroscopic asymmetry. The effect of correlation length is more dramatic on white noise heterogeneities (compare sets A and D). The asymmetry persists up to larger std when the correlation length is shorter.

Seismological measures of macroscopic asymmetry
Testing the statistical predictions associated with ruptures on a bimaterial interface in the context of natural faults requires data that reflect properties of a large number of earthquake ruptures. Indirect tests may be done with catalogue analysis (Rubin & Gillard 2000) and observations of rock damage in fault zone structures (Dor et al. 2006a(Dor et al. ,b, 2008Lewis et al. 2005Lewis et al. , 2007. It would be desirable to obtain direct tests on imaged properties of earthquake ruptures. While a number of large earthquakes have been recorded by dense near-field networks, the more numerous small earthquakes are generally not accessible for detailed seismological imaging of source properties. However, the degree of source asymmetry may be measured by the first and second order moments of the spatio-temporal distribution of seismic moment rate density (McGuire et al. 2001). To avoid the ambiguity of the seismic moment on bimaterial faults (Ben-Zion 2003; Ampuero & Dahlen 2005), we work here with the seismic potency.
One observable indicator of rupture directivity is the azimuth dependence of the duration τ of the apparent source time function, which can be determined for large events from teleseismic observations and for smaller events by empirical Green's function deconvolution (McGuire et al. 2002;McGuire 2004). Here τ is defined as the square root of the second central moment of the apparent source time function (analogous to eq. F3). It is a function of the take-off angle θ, which is the angle from the positive fault axis direction to the direction of the ray leaving from the source centroid towards the station. In the far-field and under the Fraunhofer approximation (Silver 1983): whereμ (i, j) are the second order central moments (Appendix F) of the spatio-temporal distribution of slip rate V (x, t), and c is the wave speed at the source corresponding to the selected seismic phase. The second term in eq. (10) is indicative of macroscopic source asymmetry; it is related to the apparent velocity of the instantaneous centroid and thus combines the effects of rupture directivity and of asymmetric slip distribution. A 'directivity ratio ' (McGuire et al. 2002;McGuire 2004) is defined as with values in the range −1 to 1. High absolute values, for example, |DR| > 0.5, imply strong source asymmetry. Here DR = 1 indicates strong rupture asymmetry towards the 'preferred' direction, and DR = −1 asymmetry towards the opposite direction. Fig. 16 shows histograms of DR values generated in simulations with different amplitudes of initial stress heterogeneities (std) from set D. In comparison to the potency ratio (Fig. 15), the parameter DR is a more sensitive and informative measure of macroscopic source asymmetry. For low level heterogeneities std = 0.05, all events become unilateral in the preferred propagation direction. The bias and skewness of the distribution for std = 0.15 is amplified and it becomes clear that most cases are unilateral in the 'preferred' direction, but a few cases have evolved into almost unilateral rupture in the 'opposite' direction. At std = 0.25, most ruptures have strong asymmetry in either direction (|DR| > 0.5), although the skewness is still clear. A tendency for unilateral ruptures, but with no preferred  Figure 16. Histograms of the directivity ratio (DR = ± 1 indicates strong asymmetry in the 'preferred' and opposite rupture direction, respectively) in bimaterial ruptures with heterogeneous initial stress, for three values of std indicated in each plot. The ratio of the number of events with DR > 0.5 to the number of events with DR < −0.5 is also indicated below std.
direction, was also found in a similar set of simulations performed in homogeneous medium, for sufficiently high std. Fig. 16 also indicates for each std an overall measure of the tendency for strong asymmetry, defined as the ratio of the number of events with DR > 0.5 to the number of events with DR < −0.5. This ratio is summarized in Table 5 as a function of std for all the heterogeneity sets of Table 2. Again, the macroscopic source asymmetry is statistically enhanced by short correlation length and small amplitude in the stress heterogeneities, and it persists significantly at std = 0.25.
We repeated the simulations of set A assuming lower values of the velocity-weakening scale to favour crack-like rupture. Simulations with V c = 0.1 and 0.05 m s −1 resulted in a narrower range of values of directivity ratio, with |DR| < 0.5 even when std = 0.25. This suggests that, also in presence of heterogeneous stress, significant macroscopic asymmetry is not efficiently generated under slip-weakening friction. We also noted that larger std increases the range of values of DR, and that the few cases with |DR| > 0.5 are large-scale pulses where the healing phases are generated by the interaction with the heterogeneities.

D I S C U S S I O N
We summarize here the properties of the different rupture pulses found in our simulations and their implications on earthquake mechanics, and discuss some potentially important effects that are not examined in this paper and should be the subject of future work.

Small-scale bimaterial pulses and large-scale hybrid pulses
Previous studies have shown that in-plane ruptures on a bimaterial interface with constant friction coefficient propagate in a preferred direction as wrinkle-like pulses generated by dynamic changes of the normal stress (Weertman 1980 . Wrinkle-like pulses were also generated on a bimaterial interface governed by slip-weakening friction in cases initiated with strong localized nucleation (Shi & Ben-Zion 2006). However, calculations with slip-weakening friction initiated by relatively large smooth nucleation produced bilateral cracks with a small-scale superposed wrinkle pulse in the preferred direction (Harris & Day 1997, 2005Rubin & Ampuero 2007). The simulations of this work combine two mechanisms for pulse generation, velocity-weakening friction and the bimaterial effect. In bimaterial interfaces we find two different types of rupture pulses. 'Large-scale pulses' are primarily generated by the velocityweakening friction, as in homogeneous media, but are enhanced in the preferred propagation direction by the bimaterial effect. We refer to these also as 'hybrid pulses', because of the involvement of both mechanisms. Their asymmetry is generated by the bimaterial effect on the dynamic stress drop, which favours the transition to self-sustained rupture in the preferred direction and to decaying rupture in the opposite direction. 'Small-scale pulses' detach from the main rupture front that propagates in the preferred direction, progressively or at rupture arrest on abrupt barriers. Their size is comparable or smaller than the size of the process zone of the parent rupture front. These pulses are exclusively generated and sustained by the bimaterial effect, similarly to the wrinkle-like pulses predicted under constant friction, and we refer to these also as 'pure wrinkle-like pulses' or 'pure bimaterial pulses'. Rubin & Ampuero (2007) analysed the conditions that lead, in bilateral crack ruptures under slip-weakening friction, to the generation of small-scale wrinkle-like pulses. The mechanism operates by nucleation of a healing front at the tail of the process zone in the preferred direction, where the slip gradients are large. This spawned wrinkle-like pulse smears the stress concentrations at the crack arrest tip, inducing asymmetry in the final stress distribution. This mechanism was proposed by Rubin & Ampuero (2007) to explain the asymmetry of microearthquake triggering observed by Rubin & Gillard (2000) on segments of the San Andreas fault that have a material contrast. The arguments of Rubin & Ampuero (2007) are independent on whether the main rupture propagates as a crack or as a pulse, and our calculations with velocity-weakening friction show that the small-scale daughter pulse can also detach from a largescale pulse. Our stopping ruptures also produce asymmetric static Coulomb stress changes near the final rupture tips that are more favourable for aftershock triggering near the tip in the 'preferred' direction, and this asymmetry does not correlate significantly with the macroscopic source asymmetry quantified by DR. Rubin & Ampuero (2007) noted that the rupture distance at which the wrinkle-like pulse appears is sensitive to details of the regularization of the constitutive normal stress response. This aspect of fault rheology has not been constrained unequivocally by laboratory experiments and remains unknown at geological scales. The asymmetries generated by the small-scale pulse are favoured when σ * follows closely σ at the scale of the slip-weakening process zone (eq. 5 with δ σ D c ). Here we deliberately minimized this effect by setting δ σ ≈ D c , to isolate the bimaterial effect on the large-scale pulse while keeping a regularized rheology. Note that Andrews & Harris (2005) also used δ σ ≈ D c in their slip-weakening simulations, which resulted in insignificant bimaterial effect.
In contrast, the large-scale pulse is largely independent on the details of the regularization of the normal stress response, and the bimaterial effect on this pulse is significant for a wide range of regularization timescales. In our simulations, δ c /V * = 0.035 s is longer than the typical model rise times and much longer than the timescale of the process zone, and yet the hybrid large-scale pulse produces macroscopic source asymmetry that is robust with respect to regularization details, range of stress heterogeneities and a proxy for off-fault yielding. The effects generated by the large-scale asymmetric pulse are naturally stronger than those associated with the small-scale pulse.

Effect of off-fault yielding
Ben-Zion &  showed that spontaneous generation of offfault yielding regularizes the evolution of the wrinkle-like pulse and limits the slip and rupture velocities. Rubin & Ampuero (2007) noted that off-fault yielding tends to reduce slip gradients in the process zone, and counteracts the development of the high normal stress changes required to generate the small-scale pulse. To generate significant asymmetry with these pulses while minimizing off-fault stress requires increasing the velocity contrast and decreasing the friction drop, for example, more than 20 per cent velocity contrast for μ = 0.1. In contrast, the large-scale velocity-weakening pulses require smaller slip gradients and can operate without large off-fault stresses. Andrews (2005) suggested imposing a limit on slip rate in elastic simulations as a proxy for off-fault dissipation. Our simulations of Section 3.3 incorporate this 'velocity toughening' procedure. The results of Figs 14-16 confirm that the large-scale pulses and associated asymmetric ruptures are not sensitive to this feature. However, off-fault yielding produces a 'damage' length scale and other effects (Ben-Zion & Shi 2005) that can not be accounted for by this procedure and may be important for the behaviour of natural ruptures. Several recent studies found with direct geological mapping (Dor et al. 2006a(Dor et al. ,b, 2008 and inversions of seismic fault zone waves (Lewis et al. 2005(Lewis et al. , 2007 macroscopic damage asymmetry in the structure of large faults. The width of these asymmetric damage zones can reach several hundreds of metres. The robustness and large-scale asymmetry of the hybrid pulses simulated in this paper suggest that the observed damage asymmetries may reflect a statistical preference for the occurrence of asymmetric ruptures on those fault zones.

3-D effects
The transition between decaying and sustained ruptures is very sensitive to small changes in the available energy. The dynamic changes of normal stress induced by the bimaterial effect are confined to relatively small areas of high slip gradient. In pulse-like ruptures, the control on rupture properties by the spatial variability of dynamic stress drop operates very locally, at the scale of the pulse width. In contrast, crack-like ruptures are sensitive to a weighted average of stress drop over a larger scale. Due to their higher responsiveness, pulse-like ruptures are affected more efficiently by the bimaterial effect than crack-like ruptures. Velocity-weakening friction is not the only possible mechanism to generate large-scale pulses, and we expect the asymmetry mechanism studied here to generally operate on rupture pulses, regardless of their origin.
In faults with large aspect ratio, ruptures naturally evolve into pulses, even in the absence of velocity-weakening friction (Day 1982;Johnson 1992). The bimaterial effect on these 3-D pulses can lead to predominately unilateral crustal scale ruptures, as suggested by recent 3-D calculations on bimaterial faults governed by slipweakening friction and with uniform initial stress (Dalguer & Day 2007).
Pulses in slip-weakening faults can also be generated in 3-D by strong heterogeneities (Oglesby & Day 2002). The asymmetric destabilization mechanism might similarly operate on these pulses. In our 2-D simulations with low V c (slip-weakening dominant), increasing the amplitude of the heterogeneities favours the emergence of large-scale pulses and results in a larger range of macroscopic rupture directivity (DR). However, strong heterogeneities also tend to reduce the asymmetry of large-scale pulse ruptures. It is uncertain at this stage which effect prevails in highly heterogeneous slip-weakening faults. Brietzke et al. (2007) showed that growing wrinkle-like pulses can also be generated in 3-D calculations from an initial mixedmode rupture on a fault with a constant coefficient of friction. Once a 3-D rupture front has developed a large radius of curvature, the asymptotic fields near the rupture tip are essentially 2-D. Thus, under more realistic friction, the mechanism of small-scale pulse detachment from a main rupture front should also exist in 3-D whenever rupture velocity and slip gradients are not limited too soon by the saturation effect of the finite width of the seismogenic crust.

Role of nucleation
The most prominent bimaterial effect found here is the difference in the criticality conditions for self-sustained pulse propagation in each rupture direction. This is an efficient mechanism for macroscopic source asymmetry in large-scale pulses that are close to their stability limit. This near-critical state can be plausibly achieved also under nucleation conditions that are more natural than assumed in this work. For instance, when nucleation is driven by tectonic loading near a stress concentration (Uenishi & Rice 2003), the slipping patch is expected to reach a critical state very slowly. More detailed modelling is required to explore the effects of smoother nucleation conditions on properties of bimaterial ruptures. For this purpose the friction law can be modified at low slip rates to match the classical rate-and-state framework, for which there are ample laboratory constraints.
Ampuero & Rubin (2008) showed that under rate-and-state friction with the slip law for the evolution of the state variable, pulses can develop spontaneously during slow, quasi-static nucleation. They also observed, in quasi-static simulations on bimaterial faults and in results of a linearized analysis (unpublished work), a tendency of the nucleation slip patch to migrate in the preferred direction. The combination of spontaneous pulse-like nucleation with a bimaterial effect might be an efficient and natural mechanism to initiate asymmetric large-scale dynamic pulses. However, this possibility has not been fully explored yet.
The nucleation of natural earthquakes may be completely different from the smooth picture described above (Ben-Zion 2006a,b). On a fault with strong heterogeneity at all scales, earthquake initiation might operate instead by a cascade process in which small subevents dynamically trigger larger ones. In this scenario the nucleation process, when seen from the large scale of the main shock, might appear as abrupt and, if vigorous enough it could push the fault beyond criticality in both rupture directions, counteracting the asymmetry mechanism. On the other hand, strong localized nucleation processes propagating at subshear velocities are expected to produce macroscopically asymmetric ruptures (Shi & Ben-Zion 2006). At this time it is uncertain which type or types of nucleation mechanisms are prevalent in nature.

Stress and strength heterogeneities
In models of seismicity that are inherently discrete, in the sense of Rice (1993) and Rice & Ben-Zion (1996), and that exhibit intermittent criticality, the statistical properties of the heterogeneities of fault stress and strength are not stationary but evolve in time throughout the seismic cycle. Typically, their correlation length and standard deviation increase at the approach of a large earthquake in a process associated with increasing short-scale fluctuations and decreasing large-scale roughness . Microseismicity patterns provide observational constraints on the mechanics of bimaterial faults (Rubin & Gillard 2000;Rubin 2002). We have found that rupture pulses are statistically sensitive to the amplitude (std) and correlation length of stress heterogeneities. This suggests that the intensity of the macroscopic source asymmetry that could be observed on the background seismicity might not be a stationary property of a bimaterial fault. More quantitative information about the statistical properties of heterogeneities near the end of the 'seismic cycle' of natural faults is needed to fully understand the implications of our results for large earthquakes.
To eventually stop all our simulated ruptures we arbitrarily placed strong barriers at a fixed hypocentral distance. This stopping procedure has an impact on the statistics of asymmetry, as they limit the value of DR for those ruptures that otherwise would have continued their unilateral propagation. We have verified that the asymmetry in Fig. 16 is still significant when we consider the subset of ruptures that stop spontaneously before reaching the barriers. Moreover, the steepness of the stopping barriers controls how far the small-scale pulse propagates beyond the final tips of the macroscopic rupture, which is an essential part of the aftershock asymmetry mechanism studied by Rubin & Ampuero (2007). Although it would be more natural to prescribe the initial stress in such a way that all ruptures stop spontaneously, in consistency with the organization of the fault state in an earthquake cycle model, this has not been attempted here.
Ideally, a statistical characterization of the spatial distribution of the fault state prior to large earthquakes could be consistently obtained from a complete simulation that resolves the short coseismic timescales as well as the long interseismic timescales, and that spontaneously generates a broad range of event magnitudes and heterogeneity length scales (Cochard & Madariaga 1996;Shaw & Rice 2000). Doing this while maintaining proper resolution of the dynamic nucleation stages is computationally challenging (Ben-Zion & Rice 1997;Lapusta et al. 2000), and such modelling has not been attempted yet on bimaterial faults. However, if the large-scale hybrid pulses are more essential than the small-scale pure wrinkle pulses, as this work suggests, the numerical grid resolution might not need to be as dense as previously thought, and this type of simulations might be achievable in the near future. The generation of off-fault yielding is also expected (Ben-Zion & Shi 2005) to produce a largescale smoothing effect that regularizes the bimaterial ruptures and reduces the grid size requirements.

Relevance to observations in other contexts
Laboratory experiments on bimaterial interfaces have revealed a diversity of rupture styles: unilateral or bilateral, subshear or supershear, crack-like or pulse-like. These experiments have mainly imaged the motion of the rupture front. Only recently the local slip velocity at selected points on the interface has been recorded and rupture pulses properly identified (Lykotrafitis et al. 2006). Anooshehpoor & Brune (1999) found in sliding experiments with different foam-rubber blocks and gradual loading that large ruptures propagate only in the positive direction, with properties compatible with the predictions for wrinkle-like pulses, whereas small events have various propagation directions and transient properties.
In the experiments of Xia et al. (2005), ruptures were initiated by an abrupt explosive wire, similar to the strong localized procedure used in the simulations of Andrews & Ben-Zion (1997) and later related works, but with essentially infinite propagation speed. The results produced rupture fronts in the preferred direction having the generalized Rayleigh wave speed, and rupture fronts in the opposite direction with velocities that range from sub-Rayleigh to near the slower P-wave speed. Shi & Ben-Zion (2006) showed that strong localized nucleation with subshear velocities produce subshear ruptures that tend to evolve for ranges of conditions towards unilateral wrinkle pulses in the preferred direction, whereas nucleation with supershear velocities excite also a weak front of supershear rupture in the opposite direction. As discussed earlier, with velocity-weakening friction a smooth nucleation process favours large-scale asymmetric pulses, while abrupt triggering over a sufficiently large scale can reduce the tendency for asymmetry. Laboratory observations of the effects studied in this work will require relatively gradual initiation procedures, with subshear propagation velocities, and resolution of the amount of slip generated in the opposite propagation directions.
Deep seismic events in glaciers have been shown to occur near the bottom of the glaciers (Deichmann et al. 2000;Danesi et al. 2007), which is a natural bimaterial (rock/ice) or trimaterial (rock/till/ice) configuration. Most of the glacier bottom interface is creeping and the seismicity is highly clustered in space, including identifiable repeating events. These features are similar to the repeating earthquakes on asperities surrounded by fault creep, as observed in bimaterial sections of the San Andreas fault and as expected near the brittle-ductile transitions of large faults. Deep icequakes are an interesting analogue to these earthquakes, with the advantage of an easier access to remote and direct observations. Rudnicki & Rice (2006) examined the interaction between changes of normal stress generated by large-scale contrast of elastic solids and fluid effects associated with small-scale contrasts of permeability across the fault. The effects generated by these two types of contrasts can augment each other or compete, depending on the senses of contrasts. When the two mechanisms produce competing changes, the bimaterial effects associated with the large-scale elastic contrast are likely to be larger for rupture velocities near the generalized Rayleigh wave speed. Yamashita (2007) found with quasi-static calculations that a contrast of elastic and diffusivity properties across a fault lead to strongly asymmetric distribution of afterslip. It would be interesting to examine the interaction between large-scale bimaterial contrast and small-scale fluid-related effects in simulations that account also for off-fault damage, which is expected to modify considerably the permeability and diffusivity properties.
Additional important effects that should be considered in future works include ruptures on dipping bimaterial faults (Shi & Brune 2005), effects of complex 3-D velocity structures (Archuleta et al. 2006), and interaction of the bimaterial effects with thermoelastic instability (Afferrante & Ciavarella 2007). We also note that bimaterial contrasts might be important for interpretation of geodetic data across large faults (Le Pichon et al. 2005;Fialko 2006;Wdowinski et al. 2007).

C O N C L U S I O N S
We perform a systematic parameter-space study of properties of inplane ruptures on a bimaterial interface governed by a velocity-and state-dependent friction law. The study is done primarily through numerical simulations and is augmented by analytical estimates. We characterize the range of velocity-weakening scales, nucleation parameters and initial stress distributions for which ruptures behave as cracks or pulses, decaying or sustained, and quantify the symmetry properties of the ruptures with several measures. The simulations produce large-scale pulses associated with strongly velocityweakening friction and small-scale secondary wrinkle-like pulses generated by bimaterial effects. The latter are sensitive to various poorly constrained details. In contrast, the former occur under more general conditions and lead to macroscopic rupture asymmetry that is robust with respect to regularization parameters, ranges of stress heterogeneities and a proxy for off-fault yielding. The macroscopic rupture asymmetry is manifested by significantly larger seismic potency and propagation distance in the preferred direction, similar to what was found by Shi & Ben-Zion (2006) with strong nucleation and slip-weakening friction. The rupture asymmetry is clearly quantified by the directivity ratio calculated from the second moment of the spatio-temporal distribution of slip rate. Since friction is generally velocity-dependent, and rupture pulses occur under lower stress than cracks, ruptures on bimaterial interfaces are expected to be associated with strong macroscopic asymmetry. Several additional effects, not considered in this work, may enhance or impede the generation of large scale pulses and macroscopic rupture asymmetry on bimaterial interfaces. It is important to clarify such effects with additional theoretical studies, and to test the predictions associated with ruptures on bimaterial interfaces with detailed field and laboratory observations.

A C K N O W L E D G M E N T S
JPA was funded by SPICE, a Marie Curie Research and Training Network in the 6th Framework Program of the European Commission. YBZ acknowledges support from the Southern California Earthquake Center.

A P P E N D I X A : N U C L E A T I O N L E N G T H A N D N U M E R I C A L R E S O L U T I O N
A usual indicator of the shortest length-scale that must be resolved by a numerical simulation is obtained by linear stability analysis (Rice 1993). The analysis by Rice et al. (2001) is general enough to include the friction law employed here. The correspondence between the classical rate-and-state frictional parameters a, b and L and the parameters α, β and D c in our friction law is In homogeneous media, the smallest wavelength that can become unstable after a perturbation of uniform steady-state slip is where c R is the Rayleigh velocity of the medium. We have derived this approximate result for mode II by analogy to the exact results for mode III presented by Rice et al. (2001), replacing μ by  Rice et al. (2001) and from our eq. (A4), respectively. μ/(1 − ν) and c S by c R . We have checked numerically (Fig. A1) that it is accurate enough (better than 5 per cent for ν = 0.25) for our purposes.
This critical wavelength increases monotonically as a function of V , from its quasi-static limit of λ cr at low velocities (V V c ), to its limits at high velocities ( To guarantee numerical stability the grid step size x must be smaller than the quasi-static length scale of eq. (A6). As noted by Rice et al. (2001), in absence of a direct frictional effect (α → 0) the critical wavelength vanishes (because V dyn → 0) and the problem becomes ill-posed. We have set the value of α so that this length is well resolved by our numerical grid. In bimaterial faults we estimate the quasi-static critical length by replacing μ by the effective modulus μ defined by Rubin & Ampuero (2007) based on relations of Weertman (1980).

A P P E N D I X B : T R A N S I T I O N T O R U N A W A Y C R A C K S
For crack ruptures of half-size a larger than the size of the process zone we can estimate the final rupture size by an energy balance, under the framework of fracture mechanics and Griffith's rupture criterion. Rupture stops when the energy release rate at zero rupture speed, G(a), is exactly balanced by the fracture energy, G c , We are assuming here a smooth spatial distribution of G c , otherwise the condition above should be an inequality (≤). The crack arrest position a is stable if dG da (a) < 0.
In general, G is a complicated function of rupture history and does not coincide with the energy release rate of a quasi-static crack. However, for the purpose of this semi-quantitative discussion we  Figure B2. Final rupture size a as a function of asperity size L nuc . Each curve corresponds to a different value of the strength excess parameter S = 1, 1.5, . . . , 5. Lengths are normalized by the Griffith length L c associated to the background stress drop, which is proportional to fracture energy. simply assume that G is given by the static stress intensity factor K 0 as where η is a dynamic overshoot factor. We adopt η = 4/π , the exact value in mode III for a stationary crack under dynamic plane wave loading (Ing & Ma 1997) and for bilateral cracks with uniform stress drop (Leise & Walton 2001). For an asperity with initial stress ≈ τ s and half-size L nuc ≤ a, under background stress τ 0 : Solving the equations above for a gives an estimate for the final rupture size. Fig. B1 shows a near-perfect match between the results of numerical simulations and the theoretical estimate, for the special case τ 0 − τ d = 0. For higher τ 0 , the appropriateness of the simplifying assumption B3 degrades and the prediction systematically overestimates the rupture size, but the general trends discussed next still hold. The predicted final size as a function of asperity half-size L nuc for a range of values of the background strength excess parameter S defined in eq. (7) is plotted in Fig. B2. The lengths in Fig. B2 are normalized by the Griffith length associated to the background stress drop: For each S there is a maximum L nuc above which the ruptures become unstoppable. This critical point is given by the minimum of K 0 (a). In spontaneously arrested cracks, the range of ratios of final to initial rupture size, a/L nuc , becomes wider if S is increased (decreasing the background stress τ 0 ).

A P P E N D I X C : F R A C T U R E E N E R G Y A N D D Y N A M I C S T R E S S D R O P
Following Rubin & Ampuero (2005) we estimate the fracture energy spent by our adopted friction law by mimicking the slip rate evolution in the process zone by an abrupt velocity jump, from 0 to constant V . The critical fracture energy G c is defined as the amount of energy dissipated by the frictional forces minus the work of  Table 1.  Figure C2. Ratio κ of available energy to critical fracture energy as a function of the velocity-weakening parameter V c (non-dimensional parameter v) for the values of background stress (non-dimensional parameter s) indicated in the legend. All quantities are normalized as described in the text. residual steady-state friction. Assuming a constant normal stress σ gives Taking θ (0) = 0, the state variable evolves as After integration of (C2) using (C3) we get: For an estimate of the slip rate scale V inside a crack-like rupture we use V dyna , defined in Appendix D through the balance between steady-state strength and radiation damping: where s and v are normalized measures of τ 0 and V c , respectively, also defined in Appendix D.
The estimate of G c obtained by combining eqs (C4) and (C5) is plotted in Fig. C1 as a function of V c for different values of τ 0 . At a fixed background stress τ 0 , the fracture energy decreases with increasing V c , due to the decrease of V dyna .
To quantify how easily a rupture front can be stopped, the key parameter is the ratio of available strain energy to G c , the parameter κ introduced by Madariaga & Olsen (2000). For a rupture pulse of width W and minimum dynamic stress τ dyna , We assume τ dyna = σ μ ss (V dyna ). The dynamic stress drop, τ 0 − τ dyna , plotted in Fig. C1, is also a decreasing function of V c , so it is not immediately apparent which term will dominate the energy balance. The parameter κ, normalized by is plotted in Fig. C2 as a function of V c . At sufficiently high V c it decreases: in the energy balance the effect of the reduction of dynamic stress drop with increasing V c overcomes the effect of the reduction of G c . This suggests that ruptures are easier to stop at larger V c and offers an interpretation for the shape of the boundary between decaying and sustained ruptures in Fig. 3.

A P P E N D I X D : R U P T U R E S T Y L E T R A N S I T I O N ( P U L S E / C R A C K )
Here we show that the crack/pulse transition observed in our simulations is consistent with the analysis by Zheng & Rice (1998) and Nielsen & Carlson (2000), and in better quantitative agreement with  Figure D1. Rupture style as a function of background stress τ 0 and characteristic velocity-weakening parameter V c , for a nucleation patch size fixed at L nuc = 10 m. The solid and dashed curves are pulse/crack transition boundaries predicted by critical T and critical H values, respectively. Open circles represent conditions under which the pulse/crack transition is observed in our numerical simulations, diagnosed by the divergence of the hypocentral healing time. The open square corresponds to a pulse rupture but at this or higher stress levels the transition to crack is masked by the transition to supershear rupture. the latter. We discuss on that basis the effect of the background stress level.
Crack-like ruptures are generally associated with relatively high background stress. Zheng & Rice (1998) defined a characteristic stress value τ pulse as the value of background stress for which the radiation damping line in a stress/velocity diagram (V , τ ) is tangent to the steady-state strength curve, that is, the largest τ that satisfies In their 2-D antiplane rupture simulations, Zheng & Rice (1998) obtained sustained pulses only if τ 0 was larger than a threshold value τ arrest slightly smaller than τ pulse . Zheng & Rice (1998) showed that the crack rupture mode is not viable when τ 0 < τ pulse , and argued that the transition from pulse-like to crack-like rupture occurs when the velocity-weakening rate in the bulk of the slip zone is small compared to the fault impedance (the radiation damping coefficient). They quantified this condition by the following dimensionless parameter T, defined when τ 0 > τ pulse : where V dyna is an estimate of the slip rate in the bulk of the slipping zone, given by the largest slip velocity V such that In the 2-D antiplane simulations reported by Zheng & Rice (1998), crack-like rupture required small values of T and the crack/pulse transition, although not sharp, occurred close to T = 0.5. Nielsen & Carlson (2000) introduced another dimensionless parameter, H, defined as the ratio of V c to the virtual value of V c that, for a given τ 0 , would make the radiation damping line tangent to the steady-state strength curve, that is, the smallest value of V c that satisfies (The dependence on V c is implicit in μ ss .) In their 3-D acoustic simulations, crack-like rupture required high H values and a sharp transition of rupture mode occurred at H ≈ 0.3. We introduce the following normalized notations for τ pulse , τ 0 and V c , respectively, For the specific friction law adopted here the quantities introduced above are In our simulations the characteristic velocity at the crack/pulse transition is V c ≈ 0.071 m s −1 when L nuc = 10 m, and increases weakly with L nuc . This transitional V c value corresponds to T ≈ 0.198 and H ≈ 0.496. Based on either of these critical T or H values, the crack/pulse transition boundary for different levels of background stress τ 0 is predicted to be, respectively: The predictions based on T or H do not differ significantly at low background stress (s < 0.4) and are both consistent with our simulations, as shown in Fig. D1. Both imply that the transitional V c increases with increasing τ 0 . At high τ 0 the prediction based on H agrees better with our simulations. For higher τ 0 (s > 0.55) the pulse/crack transition is masked by the transition to supershear rupture: the divergence of the hypocentral healing time is frustrated by the arrival of waves from the point of nucleation of the supershear front. The pulse/crack transition at such high stress values is seen more clearly in the antiplane mode (no supershear transition), a matter that we do not pursue here.

A P P E N D I X E : B I M A T E R I A L E F F E C T O N P U L S E S T A B I L I T Y
The stability of a rupture pulse is affected by the parameter κ defined in eq. (C6): large values of κ promote runaway rupture, small values promote rupture arrest. In this appendix, we derive a rough estimate of the bimaterial effect on the value of κ, from known relations for steady-state bimaterial pulses. Despite its crudeness, this estimate allows discussion of general trends of the transition between sustained and decaying pulses, as a function of the contrast of rock properties, dynamic friction coefficient and background stress.
For given pulse width W and critical slip D c , the parameter κ is proportional to the ratio of the squared dynamic stress drop τ to the strength drop τ s . From relations given by Weertman (1980), the normal stress change and shear stress drop within a (steady-state)  Figure E1. Contours of the ratio (c) defined in eq. A7, as a function of S-wave velocity contrast and rupture speed c (normalized by the limiting speed, c lim , which is the generalized Rayleigh or the slowest S-wave speed). The dynamic stress drop in the 'preferred' rupture direction is amplified by the bimaterial effect as (c) increases towards 1/μ d . pulse of width W and characteristic slip D, running in a bimaterial fault with speed c can be roughly estimated as σ bim ≈ ∓μ * (c)D/W (E1) The sign ∓ depends on the direction of rupture propagation, and μ * and μ are bimaterial parameters, given by Weertman (1980) and Cochard & Rice (2000), that depend on the elastic properties and density of the two rocks, and on the rupture speed c. The dynamic stress drop and the strength drop can be written respectively as where τ hom and τ s hom are the analogous quantities in the homogeneous case. Combining the equations above we obtain: where (c) and S is the nominal strength excess parameter defined in eq. (7). For given pulse width W and critical slip D c , the ratio of κ bim in a bimaterial fault to κ hom in a homogeneous medium is In defining this ratio we assume that the harmonic average of the shear moduli of the bimaterial case is equal to the shear modulus of the homogeneous case, that is, (1/μ 1 + 1/μ 2 )/2 = 1/μ. Due to the crudeness of the estimates involved in its derivation, eq. (E8) has to be taken only as an indicative order of magnitude, that will need to be refined by further work. We nevertheless explore its implications next. Fig. E1 shows that for the usual range of shear wave speed contrasts and for rupture velocities slower than the appropriate limiting speed c lim (the generalized Rayleigh wave speed c GR or the slower shear wave speed) is a positive and monotonically increasing function of c, which becomes infinite at c = c GR when this speed exists. Thus, according to eq. (E8), the bimaterial effect increases the value of κ in the 'preferred' direction and decreases it in the opposite direction. The predicted effect is enhanced by faster rupture speed, by lower background stress (higher S) and by higher dynamic friction μ d .
In particular, κ is strongly amplified in the 'preferred' direction when The reduction of κ in the opposite direction is not as pronounced. Fig. E1 shows that for typical values of dynamic friction (μ d ≤ 0.6) strong amplification ( ≥ 1.67) can only be achieved at small to moderate material contrasts (c s2 /c s1 < 1.8) and at large rupture speeds (c ≈ c lim ). For instance, for μ d = 0.6 and 30 per cent velocity contrast eq. (E9) is met at c ≈ 0.97 c GR , but requires higher c/c lim if the material contrast is different. Rupture speeds more typical of natural earthquakes (c/c lim < 0.9) produce moderate values (<0.6). However, rupture pulses are also responsive to small changes of κ, and a moderate increase of , well before the extreme condition of eq. (E9) is met, can turn a decaying pulse into a runaway pulse. At moderate rupture speeds (c/c lim < 0.9) Fig. E1 shows that increasing the velocity contrast (even beyond the domain of existence of c GR ) systematically increases , enhancing the transition to runaway pulses in the 'preferred' direction.