AG_56.05.13_BIZZARRI_finalonline:Layout 6 ANNALS OF GEOPHYSICS, 56, 5, 2013, S0560; doi:10.4401/ag-6279 S0560 Calculation of the local rupture speed of dynamically propagating earthquakes Andrea Bizzarri Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Bologna, Bologna, Italy ABSTRACT The velocity at which a propagating earthquake advances on the fault surface is of pivotal importance in the contest of the source dynamics and in the modeling of the ground motions generation. The rupture speed (vr ) is one of the results provided by spontaneous dynamic models of ruptures, in that it is a part of the solution and it is not imposed a priori, like in non spontaneous models or in kinematic models. Since vr is numerically retrieved from the spatial distribution of the rupture times (tr ), a well– constrained value of vr in a given fault node is important. In this paper we focus on the numerical determination of vr. By comparing different numerical schemes to compute vr from tr we show that, in general, central finite differences schemes are more accurate than forward or backward schemes, regardless the order of accuracy. Overall, the most efficient and accurate algorithm is the five–points stencil method at the second–order of accuracy. These conclusions hold for homogeneous and heterogeneous configurations and for different constitutive models, such as the slip– weakening law and the rate– and state–friction governing equations. It is also shown how the determination of tr can affect vr; numerical results in- dicate that if the fault slip velocity threshold (vl ) used to define tr is too high (vl ≥ 0.1 m/s) the details of the rupture are missed, for instance the jump of the rupture front eventually occurring for 2–D supershear rup- tures. On the other hand, for vl ≤ 0.01 m/s the results appear to be sta- ble and independent on the choice of vl. Finally, it is confirmed that in the special case of the linear slip–weakening friction law the computa- tion of vr from the threshold criterion on the fault slip velocity and from the achievement of the maximum yield stress are equivalent. 1. Introduction 1.1. Scientific rationale One of the most important physical observables of an earthquake rupture is represented by the rupture speed, which quantifies how fast a rupture front is prop- agating on a fault surface. It is well known that the vari- ations in the rupture speed (and thus the earthquake acceleration) control the radiation of the seismic waves in the surrounding medium, which in turn are respon- sible of most damage from an engineering point of view [Madariaga 1983]. Many studies [Bizzarri and Spu- dich 2008, Dunham and Bhat 2008, Bizzarri et al. 2010] clarified that there are systematic differences, in terms of amplitude of the motion and of the frequency content, between ruptures which propagate at a speed lower than the S wave speed, vS (i.e., subshear events) and ruptures which accelerate above vS (i.e., supershear earthquakes). Numerical models of earthquakes are powerful tools to investigate the behavior of faulting in realistic ambient conditions and provide the opportunity to sim- ulate different shaking scenarios, in the framework of a deterministic approach. Spontaneous dynamic models simulate, as a part of the solution, the expansion through time of the broken region of a fault structure, because the rupture speed is not prior imposed — as in non spontaneous models — but it is completely con- trolled by the adopted fault governing model and by the initial conditions. Indeed, spontaneous dynamic models provide the spatial distribution of the rupture times (tr), which namely define the location of the rup- ture front during the dynamic propagation of the rup- tures. The rupture speed vr has then to be retrieved numerically from tr. This poses a relevant numerical problem, which is the main focus of the present study. In seismology, determination of vr is very impor- tant. From a purely theoretical point of view pure mode II, 2–D singular cracks are known to possess only limited values of allowable rupture speeds, namely from 0 to the Rayleigh wave speed (vR) and from vS up to the limiting speed of mode II (which is the P wave speed, vP). In this case they then predict [Freund 1979, Broberg 1989, 1999] the existence of a so–called forbid- den zone (between vR and vS). Recently, Bizzarri and Das [2012] show that in 3–D, non singular ruptures this forbidden zone disappears and that the rupture can ad- vance with rupture speeds continuously distributed be- Article history Received January 9, 2013; accepted September 6, 2013. Subject classification: Earthquake dynamics, Rupture speed, Supershear earthquakes, Rheology and friction of fault zones, Computational seismology. tween 0 and vP. This result of course requires a very ac- curate estimate of vr. Moreover, a well–constrained value of vr is impor- tant not only in spontaneous dynamic rupture models (as such considered by Bizzarri and Das [2012] and in the present study), but also in the case of kinematic models of earthquakes. Indeed, in these models the an- alytical evolution of the fault slip velocity is a priori specified and this requires the exact determination of the rupture times (see for instance Ide and Takeo [1997], among many others). Finally, it is worth mentioning that the knowledge of the rupture speed from numerical models is very rel- evant because recently it has been developed an analy- sis on possible spatial correlations existing between the various dynamic variables, such as rupture speed, total developed slip, fracture energy density, peak fault slip velocity and stress drop [Guatteri et al. 2004, Song et al. 2009, Schmedes et al. 2010, Song and Sommerville 2010, Bizzarri 2010b, 2012a]. These correlations are ex- tremely important, because they could be inserted as constraints in kinematic modeling of faults, on which current practice in seismic engineering relies. 1.2. Methodological approach In this paper we consider the simplified problem of a single, isolated, vertical, planar, strike–slip fault, embedded in a perfectly elastic medium. The fault geometry is the same as that used Bizzarri and Spu- dich ([2008]; see their Figure 3). The elastodynamic problem is solved numerically, as described in details in Bizzarri and Cocco [2005], and the parameters adopted in the present study are tabulated in Table 1. We hypothesize that the fault obeys the linear slip– weakening constitutive law [Ida 1972], which postu- lates that the shear traction on the fault drops from its maximum, yield level xu down to the residual value xf, over a characteristic slip distance d0 (see Equation 25 in Bizzarri [2011]). There is no need to overempha- size that this kind of friction model can be merely re- garded as one possibility to describe the non singular stress decay during a seismic failure, and that more elaborated governing equations can be employed (see the discussion in Bizzarri [2009, 2011]). Indeed, in Sec- tion 4 we will consider the Ruina–Dieterich law (Ruina [1983]; see also Equation 35 in Bizzarri [2011]) in order to generalize the conclusions obtained in the slip–weakening case. Moreover, we disregard here the effects on the dynamic rupture propagation of the variations of the effective normal stress, due for in- stance to the thermal pressurization of pore fluids, which has been already discussed elsewhere [Bizzarri and Cocco 2006a, 2006b]. The idealized geometry and the simplified friction law adopted here do not affect the results presented and discussed in the present study, in that they can be exactly extended to other, ANDREA BIZZARRI 2 Parameter Value Medium and discretization parameters Lamé’ s constants, m = G 35.9 GPa S wave velocity, vS 3.464 km/s Rayleigh velocity, vR 3.184 km/s P wave velocity, vP 6 km/s Eshelby velocity, vE = vS 4.899 km/s Fault length, Lf 16 km Fault width, Wf 11.6 km Spatial grid size, Dx 5 m (a) Time step, Dt 1.2 × 10−4 s (b) Coordinates of the hypocenter, H ≡ (x1H, x3H) (8,7) km Fault constitutive parameters Magnitude of the effective normal stress, vn eff 120 MPa Magnitude of the initial shear stress (pre–stress), x0 73.8 MPa Static friction coefficient, nu 0.677 ( ↔ xu = 81.24 MPa ) Dynamic friction coefficient, nf 0.46 ( ↔ xf = 55.20 MPa ) Characteristic slip–weakening distance, d0 0.4 m Table 1. Parameters adopted in the present paper. (a) The spatial discretization ensures a good resolution of the breakdown zone length (Xb). On average we have: = /Dx = 100. (b) For the adopted parameters the Courant–Friedrichs–Lewy ratio, ~CFL vSDt/Dx, equals 0.083 and a conservative estimate [e.g., Archuleta and Frazier 1978] of the critical frequency for spatial grid dispersion, facc (s) = vS/(6Dx), equals 115 Hz. As required by explicit time stepping schemes we satisfy the condition Dt ≤ Dx/(2 vP). 2 df = 3 more complicated configurations. In Sections 2 and 3 we will consider faults with homogeneous properties, while in Section 4 we will discuss the case of hetero- geneous rheology. Once the fault starts to propagate from the im- posed hypocenter H (the nucleation strategy is thor- oughly described in Bizzarri [2010a]), it spreads bilaterally over the fault releasing energy and emitting seismic radiation. In each node (x1, x3) of the fault plane we define the rupture time tr(x1, x3) as the first instant when the L2–norm of the fault slip velocity in that point, v(x1, x3), exceeds a threshold value vl. We will dis- cuss the choice of the value of vl and its effects of the computation of the rupture speed in Section 3. Once the rupture times are known, the definition of the rupture tip C at time t is straightforward: (1) which has the following numerical counterpart: (2) where the symbol denotes the discrete equivalent of a generic quantity q, the doublet (i,k) identifies a fault node ( = iDx and = (k−1)Dx, being Dx the spatial grid size) and m defines the time = m Dt, being Dt the time step. From tr one can also compute, in each node of the fault, the rupture speed vr simply as the inverse of the slowness [e.g., Bizzarri and Spudich 2008]: (3) In order to explicitly know at which rate the rup- ture tip C is enlarging in a given direction, we introduce the rupture speed vector vr, defined as it follows [see also Pulido and Dalguer 2009, Hok and Fukuyama 2011]: (4) where (5) where are the unit vectors along the strike and depth directions, respectively. The angle U mathemati- cally defines the direction of enlargement of the rup- ture front C. In particular, when U = 0, we are in the pure mode II (i.e., inplane) direction; U = r/2 corre- sponds to the pure mode III (i.e., antiplane) direction. (Note that in Equations 4 and 5 we have omitted the explicit dependence of variables on the spatial coordi- nates x1 and x3 for brevity of notations.). In this study we will explore the different possible computations of vr from Equation (3) in order to pro- vide a clear indication of the most accurate and stable computational method of the rupture speed. We will also discuss the determination of the tr, which affects, by definition, the calculation of vr. 2. Computation of the rupture speed 2.1. Numerical schemes In Figure 1 it is imaged the fault slip at t = 1.92 s re- sulting from a numerical experiments; the white por- tion of the fault represents the unbroken area and the rupture front C is also indicated as the black line marked with label “0.00”. This is the curve which sep- arates the fault nodes that are slipping and releasing en- ergy from those that are still at rest at this time. The calculation of the rupture speed vr requires the numerical computation of the spatial derivatives ap- pearing in Equation (3). In full of generality, we can write CALCULATION OF EARTHQUAKE RUPTURE SPEED Figure 1. Distribution of the fault slip developed by a 3–D synthetic earthquake at time t = 1.92 s. The adopted parameters are listed in Table 1. Due to the symmetry with respect to the hypocenter H, we report only one half of the total fault length Lf. , , for the first time,t v x x t vx x l1 31 3 ; $C =^ ^ ^h h h" , ,, i ki k t m( )m ;C = =r^ ^h h" ,L I , , .v x x x xt 1 ( , r x x r 1 3 1 3)1 3 d< < =^ ^ h h ,v cos sinvr r U U= ^ h arctg arctg x x t t x t x t r r r r 3 1 $ $ d d 2 2 2 2 /U = 1 3 p r q q q q e t v u u u u oV V qJ x1K x3K t( )mI andx x1 3V V (6) and (7) ( Just for an example, for the central finite difference scheme at second order of accuracy we have: lstart = − 1, lend = 1, c−1 = − 1/2, c0 = 0 and c1 = 1/2, so that In order to see how much the values of vr is af- fected by the choice of the numerical algorithms we test several finite differences schemes, at different or- ders of accuracy and of different type, either forward or central. We also consider a modification of the central scheme at second order of accuracy, which comes from the so–called five–points stencil [e.g., Lapidus and Pin- der 1999]: (8) and (9) A compendious summary of the considered ap- proximations is reported in Table 2. These numerical schemes, as well as the coefficient reported in Table 2, are well–known standard approaches to estimate deriv- atives and they have been introduced in a quite vast body of literature. However, a systematic comparison of these schemes in the computation of the rupture speed in the source dynamic community is presently missed. 2.2. Numerical results: example of a 3–D homogeneous rupture To compare the results obtained by adopting the different numerical approaches described in Section 2.1, instead of comparing the spatial distribution of the re- sulting rupture speed over the whole fault plane, we se- lect four profiles, which are at the hypocentral depth and aligned along the strike direction (x1), at the strike coordinate of the imposed hypocenter and aligned along the dept (x3) and along two mixed–mode direc- tions, having azimuth angles of 45° and 30° with re- spect to x1. The two former profiles represent the situations where the propagation is basically mode II and mode III, respectively. The two latter profiles are important, in that both the two non null components of slip, slip velocity and traction coexist (recall that the geometry of the problem is truly 3–D). The first and fourth profiles are indicated by black and red dashed lines, respectively, in Figure 2a of Bizzarri and Das [2012]. The results are presented in Figures 2 to 5, re- spectively. We exclude the values of vr computed within the initialization patch, where the nucleation strategy is imposed (see Section 1.2 and Bizzarri [2010a] for further details). First of all, we observe that it is not true that an increase of the order of accuracy always produces more stable results; indeed, the level of fluctuations ANDREA BIZZARRI 4 i Type of the scheme Order of accuracy Coefficients of different terms, cl – 4 – 3 – 2 – 1 0 1 2 3 4 1 forward 1 0 0 0 0 – 1 1 0 0 0 2 forward 2 0 0 0 0 – 3/2 2 – 1/2 0 0 3 central 2 0 0 0 – 1/2 0 1/2 0 0 0 4 forward 3 0 0 0 0 – 11/6 3 – 3/2 1/3 0 5 forward 4 0 0 0 0 – 25/12 4 – 3 4/3 – 1/4 6 central 4 0 0 1/12 – 2/3 0 2/3 – 1/12 0 0 7 central 8 1/280 – 4/105 1/5 – 4/5 0 4/5 – 1/5 4/105 – 1/280 8 central, from 5– points stencil 2 n/a; see Equations (8) and (9) Table 2. Different numerical schemes used to discretize Equation (1). The various multiplier coefficients {cl } of Equations (6) and (7) are also reported. For each type of approximation the first and the last non null coefficients determine lstart and lend of Equation (6), respectively. , ,x t x x x c t i l k 1 r l l l l l r 1 1 3 start end 2 2 , D + = = ^ ^h hK/ , ,x t x x x c t i k l 1 r l r l l l l 1 3 3 start end 2 2 , D + = = ^ ^h hK/ , , , .)x t x x x t i k t i k1 2 2 11 r r r 1 1 32 2 , D - - + + ^ ^ ^ h h hc mK K , , , , , x t x x x t i k t i k t i k t i k 1 4 1 1 4 1 1 1 1 1 1 r r r r r 1 1 32 2 , D + + - - - + + + - - - - ^ ^ ^ ^ ^ h h h h h K K K K , , , , , x t x x x t i k t i k t i k t i k 1 4 1 1 4 1 1 1 1 1 1 r r r r r 3 1 32 2 , D + + - - + - + - + - - - ^ ^ ^ ^ ^ h h h h h K K K K 5 CALCULATION OF EARTHQUAKE RUPTURE SPEED 0 5 10 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo c it y ( k m /s ) mode II itype = 1 - acc. = 1 0 5 10 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo c it y ( k m /s ) mode II itype = 2 - acc. = 2 itype = 3 - acc. = 2 itype = 8 - acc. = 2 0 5 10 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo c it y ( k m /s ) mode II itype = 4 - acc. = 3 0 5 10 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo c it y ( k m /s ) mode II itype = 5 - acc. = 4 itype = 6 - acc. = 4 0 5 10 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo c it y ( k m /s ) mode II itype = 7 - acc. = 8 0 5 10 3 3.5 4 4.5 5 5.5 6 Distance from the hypocenter (km) R u p tu re v e lo c it y ( k m /s ) v P v S v R vS vR forward central 0 2 4 6 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo c it y ( k m /s ) mode III itype = 1 - acc. = 1 0 2 4 6 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo c it y ( k m /s ) mode III itype = 2 - acc. = 2 itype = 3 - acc. = 2 itype = 8 - acc. = 2 0 2 4 6 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo c it y ( k m /s ) mode III itype = 4 - acc. = 3 0 2 4 6 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo c it y ( k m /s ) mode III itype = 5 - acc. = 4 itype = 6 - acc. = 4 0 2 4 6 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo c it y ( k m /s ) mode III itype = 7 - acc. = 8 0 2 4 6 3 3.5 4 4.5 5 5.5 6 Distance from the hypocenter (km) R u p tu re v e lo c it y ( k m /s ) v P v S v R vS vR forward central Figure 2. Comparison between different numerical computations of the rupture speed, as described in the main text (see Section 2.1). The profile is taken at the hypocentral depth and parallel to the strike direction (therefore the propagation is basically mode II). Blue and red curves refer to forward and central schemes, respectively. The legends report the type of the schemes, as listed in Table 2. For a better comparison, each plot reports numerical schemes having the same level of accuracy. Values of vR and vS are also reported as purple and green lines, re- spectively. The parameters are those listed in Table 1. Figure 3. The same as in Figure 2, but now for a profile along the x3–direction (i.e., the propagation is basically mode III). strongly depend on the type of the adopted numerical approach. In other words, high–order methods tend to have larger oscillations. Interestingly, the central schemes produce less fluctuations with respect to the forward schemes (compare red and blue curves in Fig- ures 2 to 5). We can also conclude that the five–points stencil (Equations 8 and 9; black curves in Figures 2 to 5) gives the best results in absolute terms; it is even bet- ter than the central schemes at fourth–order. Com- pared to the central scheme having the same order of accuracy (the second), we can also note that the five– points stencil is more stable, especially in the cases of mode III propagation (Figure 3) and mixed–modes (Fig- ure 4 and 5). Indeed, it is able to fit the overall trend, but oscillations are smaller than in the other cases. By looking at Figure 2 and 5 it emerges that, con- trarily to 2–D singular cracks, for 3–D non singular rup- tures (for which the stress release is not abrupt, but it is accomplished over a process zone of finite width), the so–called forbidden zone of rupture speeds disappears. Indeed, as already discussed in more details by Bizzarri and Das [2012], in 3–D all the values of rupture speed in between 0 and to the limiting speed of P wave speed are admissible. The results discussed above do not depend on the specific spatio–temporal discretization, but only on the ratio between Dx and Dt. This ratio basically describe the numbers of the discrete levels of vr (clearly visible as horizontal ensembles of dots in Figure 6). Remark- ably, the previous results hold also for 2–D geometries; in Figure 6 we report the results pertaining to a 2–D, pure mode II rupture, which develops on a fault still obeying the linear slip–weakening friction law. The so- lution of the elastodynamic problem is performed by using the finite difference scheme described in Bizzarri et al. [2001] and references cited therein; the parame- ters are the same as for the 3–D case (see Table 1); the only exception is the spatio–temporal discretization (now Dx = 1 m and Dt = 2.5 × 10−5 s). Note that in this case the values of Dx and Dt are different with respect to the 3–D simulation. Nevertheless, also in this case the forward schemes give more oscillating estimates of the rupture speed, compared to the predictions of the cen- tral schemes. Overall, the second–order accurate cen- tral difference approach produce the smaller oscillations (note that in the case of a 2–D rupture we are able to compute only the schemes 1 to 7 of Table 2). In conclusion, the five–points stencil method at the second–order of accuracy appears to be to most effi- cient (for a purely computational point of view) and ac- curate (in terms of oscillations) method to compute the rupture speed from the rupture times. ANDREA BIZZARRI 6 0 5 10 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) mixed - mode 45° itype = 1 - acc. = 1 0 5 10 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) mixed - mode 45° itype = 2 - acc. = 2 itype = 3 - acc. = 2 itype = 8 - acc. = 2 0 5 10 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) mixed - mode 45° itype = 4 - acc. = 3 0 5 10 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) mixed - mode 45° itype = 5 - acc. = 4 itype = 6 - acc. = 4 0 5 10 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) mixed - mode 45° itype = 7 - acc. = 8 0 5 10 3 3.5 4 4.5 5 5.5 6 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) v P v S v R vS vR forward central Figure 4. The same as in Figure 2, but now for a mixed–mode profile, having an azimuth angle of 45° with respect to the x1–direction. 7 CALCULATION OF EARTHQUAKE RUPTURE SPEED 0 5 10 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) mixed - mode 30° itype = 1 - acc. = 1 0 5 10 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) mixed - mode 30° itype = 2 - acc. = 2 itype = 3 - acc. = 2 itype = 8 - acc. = 2 0 5 10 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) mixed - mode 30° itype = 4 - acc. = 3 0 5 10 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) mixed - mode 30° itype = 5 - acc. = 4 itype = 6 - acc. = 4 0 5 10 0 1 2 3 4 5 6 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) mixed - mode 30° itype = 7 - acc. = 8 0 5 10 3 3.5 4 4.5 5 5.5 6 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) v P v S v R vS vR forward central 0 1 2 3 4 5 6 0 5 10 Distance from the hypocenter ( km ) R u p tu re v e lo c it y ( k m /s ) itype = 1 - acc. = 1 0 1 2 3 4 5 6 0 5 10 Distance from the hypocenter ( km ) R u p tu re v e lo c it y ( k m /s ) itype = 2 - acc. = 2 itype = 3 - acc. = 2 0 1 2 3 4 5 6 0 5 10 Distance from the hypocenter ( km ) R u p tu re v e lo c it y ( k m /s ) itype = 4 - acc. = 3 0 1 2 3 4 5 6 0 5 10 Distance from the hypocenter ( km ) R u p tu re v e lo c it y ( k m /s ) itype = 5 - acc. = 4 itype = 6 - acc. = 4 0 1 2 3 4 5 6 0 5 10 Distance from the hypocenter ( km ) R u p tu re v e lo c it y ( k m /s ) itype = 7 - acc. 8 vS vR forward central Figure 5. The same as in Figure 2, but now for a mixed–mode profile, having an azimuth angle of 30° with respect to the x1–direction. Figure 6. The same as in Figure 2, but now for a 2–D, purely mode II rupture. Due to the 2–D nature of the rupture we are able to compute only the schemes 1 to 7 of Table 2. The nearly vertical clouds of points appearing at about 2 km from the imposed hypocenter are numeri- cal noise, due to the jump of the rupture front occurring in 2–D. 3. Effect of minimum detectable slip velocity on rup- ture time estimates 3.1. The importance of vl In addition to the need to accurately compute the rupture speed from the rupture time array tr (discussed in Section 2.2), an accurate estimate of vr strongly relies on the accuracy of the determination of tr. As stated in Section 2.1 and in agreement with other studies [e.g., Day et al. 2005, Rubin and Ampuero 2005, Bizzarri and Spudich 2008], tr is defined in each node of the fault sur- face as the time level at which the fault slip velocity at that node first exceeds a threshold value (vl). To quan- titatively evaluate the effects of the choice of vl in the determination of tr we perform several numerical ex- periments by varying the value of the threshold veloc- ity and keeping the other parameters unchanged. We perform the simulations in 2–D, with the same param- eters as in Table 1. In Figure 7a we report the resulting values of the ruptures times, as depending on the choice of vl. In gen- eral, we observe a dependence of tr on vl especially in the correspondence of the bifurcation of the rupture tip, i.e., when the rupture jumps to the supershear regime. Indeed, in 2–D the rupture speed does not cross the so–called forbidden zone (i.e., the interval be- tween the Rayleigh wave speed and the S wave speed) and correspondingly the rupture front exhibits a jump [e.g., Bizzarri et al. 2001 and references cited therein]. The “mother” front remains sub–Rayleigh, while the “daughter” front accelerates from vS up to the com- pressional wave speed, which represents the limiting speed for this rupture mode [Burridge 1973, Das and Aki 1977, Rosakis et al. 1999, among others]. From Figure 7a it is clear that for vl = 0.1 m/s (green circles) the separation of the two rupture fronts are quite different with respect to the prediction at lower vl. In par- ticular, for vl = 0.5 m/s (purple symbols) we can not see the bifurcation, as the rupture front appears to be a con- tinuous line. Correspondingly, if one looks at the result- ing rupture speed (see Figure 7b), it is possible to see that the forbidden region disappears for vl ≥ 0.1 m/s. This is physically reasonable, because it indicates that when the threshold velocity is too high some details of the daugh- ter rupture are not captured by the criterion used to de- fine tr. We are aware of the fact that 0.5 m/s is a relatively high value to detect the rupture front; we have consid- ered it just to complete the exploration of the parameter space. Moreover, that value is not improper; in Figure 8 we plot the time evolution of the fault slip velocity for the 2–D numerical experiment in different locations of the fault. We can clearly see from that figure that even in the first stages of the rupture, when it is still subshear, the peaks of v significantly exceed 0.5 m/s. This becomes even more true as long as the rupture propagates; after birth of the supershear rupture front the peaks are even larger, as already known since Bizzarri et al. [2001]. In spite of the relative infrequency of large values of fault slip velocity (namely greater than 1 m/s) obtained in kinematic modeling of faults, it should be noted that large v is routinely found in dynamic models of fault, not only when the fault is very unstable (due for instance to the incorporation of the thermal pressurization of pore ANDREA BIZZARRI 8 Figure 7. Sensitivity of the determination of the rupture times tr from the choice of the threshold velocity vl. The simulations pertain to 2–D ruptures having the parameters of Table 1, but different val- ues of vl, as listed in the legends. (a) Distribution of tr as a function of the distance from the imposed hypocenter, with an inset report- ing a zoom in the correspondence of the rupture tip bifurcation. At a given rupture time, the fault nodes at distances larger than the lo- cation of the symbols are not rupturing (unbroken parts). (b) Re- sulting rupture speeds, compute by adopting the central scheme at second order of accuracy (which has been proved to be the best computational approach). Note that for vl ≥ 0.1 m/s the bifurcation ( jump) of the rupture front is not well resolved and correspondingly the forbidden region is penetrated (values of vR and vS are also re- ported dashed horizontal lines). The results are independent on the choice of the threshold velocity for vl ≤ 0.01 m/s. 9 fluids; e.g., Bizzarri and Cocco [2006b] or of the me- chanical lubrication; e.g., Bizzarri [2012b]), but also in presence of spatial heterogeneities (e.g., Bizzarri et al. [2010]; their Animation S3). On the other hand, for small values of vl the results are quite similar. In particular, for vl ≤ 0.01 m/s the rup- ture times are indistinguishable and independent on the choice of the threshold velocity (Figure 7a). Remarkably, the same holds also for the rupture speeds (Figure 7b). This is not surprising, because the slip velocity histories in the case of the slip–weakening friction law are character- ized by an abrupt onset (i.e., a fast accelerating phase). As an overall conclusion, we have found that the value of vl = 0.01 m/s to define the occurrence of a rupture in a given node of the fault is optimal. The re- sults presented here pertain to 2–D ruptures, but the same conclusion holds also in the case of 3–D ruptures. 3.2. Comparison between two criteria to identify the rupture front As well–known [Bizzarri et al. 2001, among others], the slip–weakening friction gives a more abrupt fault slip acceleration in the early stages of the rupture com- pared to the more elaborated rate– and state–depen- dent friction laws [e.g., Ruina 1983]. In particular, within the framework of the linear slip–weakening con- stitutive model the fault slip velocity quickly increases from zero to non zero values at the rupture onset. This is intimately due to the sliding logic; when the traction first reaches the upper yield stress xu the slip com- mences. The fact that the fault slip velocity v becomes non zero simultaneously when the traction reaches xu is peculiar of the slip–weakening constitutive models; with the rate– and state–dependent friction laws v is non null also before the rupture onset. This is impor- tant, because for the slip–weakening law we have an exact estimate of the occurrence of the rupture (and thus an exact determination of tr) against which we can compare the estimates arising from the threshold crite- rion in terms of v. To quantitatively compare the rup- ture times defined from the vl criterion (with vl = 0.01 m/s) and those defined when x first reaches xu we compute the misfit m: (10) The spatial distribution of m is reported in Figure 9 from which it is possible to see that the difference be- tween the rupture times defined in the two different ways are in general very small (at maximum we have a difference of 0.20 %). The differences are mainly con- centrated near the boundary discriminating between the sub– and the supershear regime (black line in Figure 6) and in the neighboring of the imposed hypocenter (where the nucleation is imposed and the linear slip– weakening law is not yet operating; see Bizzarri [2010a]). From the inspection of Figure 9 one also has that tr (from xu) ≤ tr (from v l ) (because of m ≥ 0); this is phys- ically reasonable, because in general first x reaches xu and subsequently v exceeds vl. 4. Results for heterogeneous configurations and dif- ferent friction laws In the previous sections (2.2 and 3) we have con- sidered the simplest case of a fault having homoge- neous properties. Indeed, we know that real earthquake CALCULATION OF EARTHQUAKE RUPTURE SPEED 1.50 1.25 1.00 0.75 0.50 0.25 0.0 2.5 5.0 7.5 10.0 12.5 15.0 Time ( s ) F a u lt sl ip v e lo ci ty ( m /s ) x 1 = 7 k m x 1 = 6 k m x 1 = 4 k m x 1 = 7 .5 5 k m x 1 = 2 km v l = 0 .5 m /s Figure 8. Time evolution of the fault slip velocity v at difference spatial coordinates (indicated near each slip velocity peak). The results per- tain to the 2–D simulation of Figure 7. The horizontal gray dashed line indicates the greater threshold value of v employed in Figure 7 to detect the location of the rupture front (vl = 0.5 m/s). , , , , .m x x t x x x x t x xt 100 (from ) (from )(from ) r v rr v 1 3 1 3 1 3 1 3 l ul = - x ^ ^ ^ ^ h h h h are heterogeneous in terms of dynamic parameters and frictional properties. A relevant body of literature deals with the numerical simulation of the dynamic propagation of ruptures on heterogeneous fault sur- faces [far of being exhaustive, Tinti et al. 2005, Laval- lée et al. 2006, Bizzarri et al. 2010, Andrews and Barall 2011], so that it is very important to quantify the relia- bility of the rupture speed computed in these situa- tions. To illustrate how the different numerical schemes described in Section 2.1 perform, we select a hetero- geneous configuration in which the initial shear stress is still aligned in the strike direction, but it is now het- erogeneous. In particular, we assume that its magni- tude x0 decays a k –1 at high wavenumbers k, namely it has a power spectral density described by Equation (21) in Bizzarri [2010b]. The spatial distribution of the x0 is reported in Figure 10a of Bizzarri [2010b] and the cor- responding vr is in Figure 10d of the same paper. In Figure 10 we compare the computations of vr with the different numerical schemes considered in this work (see Section 2.1) along the mode II profile. These results basically confirm and reinforce what we have previously found in homogeneous conditions (see Sec- tion 2.2); in general higher orders of accuracy do not correspond to small fluctuations in the resulting values of the rupture speeds and the five–points stencil aver- age (black line in Figure 10) gives more stable estimates (small oscillations) of vr. The parameters adopted in this heterogeneous model are listed in Table 1 of Bizzarri [2010b]; it is apparent that it is characterized by different values of Dx and Dt (greater than those used in homo- geneous conditions; see Table 1 of the present paper), different elastic constant (Lamè constants, P and S wave speeds), a different Courant–Friedrichs–Lewy ratio, dif- ferent values of the frictional levels and different char- acteristic slip–weakening distances. We have deliberately chosen a rather different configuration in order show that the conclusions reached in the present study do not depend on the specific configuration we use. As thoroughly discussed by Bizzarri [2011], there is no a general consensus on the most appropriate con- stitutive model to be applied to describe the traction evolution during an earthquake failure. In addition to the linear slip–weakening law considered so far, we also consider the Ruina–Dieterich friction law [e.g., Ruina 1983], in which the dependence of the fault traction is on the slip velocity and on one state variable, account- ing for the memory of the whole history of the seis- mogenic fault. This governing model (see Equation 35 in Bizzarri [2011]) is inherently different with respect to the linear slip–weakening law, but has received a lot of attention in the modeling of the seismic source processes (again, far of being exhaustive, Bizzarri and Cocco [2005], Lapusta and Liu [2009]). We have then consid- ered this model to test the different numerical schemes for the computation of vr. The considered model has the same heterogeneous initial shear stress as the pre- vious model; the parameters are listed in Table 1 of Biz- zarri [2010b] and the resulting spatial distribution of the rupture speed is imaged in Figure 10c of that paper. It is clear that the heterogeneities cause a very compli- cated evolution of the rupture, which is now charac- terized by some accelerations and further decelerations of the rupture front, which in turn has a very complex shape (compared to the homogeneous configuration). Figure 11 reports the results of the resulting rupture speed computed with the different numerical schemes. It is apparent that now the absolute level of fluctuations is even larger than in the previous heterogeneous case (Figure 10), where the rupture is governed by the linear slip–weakening law. This difference can be easily ex- plained by considering the different global rupture be- haviors in the two cases; in the simulation with the Ruina–Dieterich law the rupture heals in many patches of the fault surface, while in the slip–weakening case, it ANDREA BIZZARRI 10 Figure 9. Misfit between the rupture times defined from the vl cri- terion and those defined when x first reaches xu. The quantity m re- ported in the figure is defined in Equation (10). The simulation pertains to a 3–D rupture with the parameters are tabulated in Table 1. The transition between sub– and supershear regime is indicated with the black line. 11 CALCULATION OF EARTHQUAKE RUPTURE SPEED Figure 10. The same as in Figure 2, but now for a heterogeneous configuration, which is described in details in Section 4. Figure 11. The same as in Figure 10, but now for a heterogeneous configuration obeying to the Ruina–Dieterich friction law (see Section 4 of the main text for further details). 0 2 4 6 0 1 2 3 4 5 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) mode II itype = 1 - acc. = 1 0 2 4 6 0 1 2 3 4 5 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) mode II itype = 2 - acc. = 2 itype = 3 - acc. = 2 itype = 8 - acc. = 2 0 2 4 6 0 1 2 3 4 5 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) mode II itype = 4 - acc. = 3 0 2 4 6 0 1 2 3 4 5 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) mode II itype = 5 - acc. = 4 itype = 6 - acc. = 4 0 2 4 6 0 1 2 3 4 5 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) mode II itype = 7 - acc. = 8 0 2 4 6 2.5 3 3.5 4 4.5 5 5.5 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) v P v S v R vS vR forward central 0 2 4 6 0 1 2 3 4 5 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) mode II itype = 1 - acc. = 1 0 2 4 6 0 1 2 3 4 5 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) mode II itype = 2 - acc. = 2 itype = 3 - acc. = 2 itype = 8 - acc. = 2 0 2 4 6 0 1 2 3 4 5 Distance from the hypocenter (km) R u p tu re v e lo ci ty ( km /s ) mode II itype = 4 - acc. = 3 0 2 4 6 0 1 2 3 4 5 Distance from the hypocenter (km) R u p tu re v el o ci ty ( km /s ) mode II itype = 5 - acc. = 4 itype = 6 - acc. = 4 0 2 4 6 0 1 2 3 4 5 Distance from the hypocenter (km) R u p tu re v el o ci ty ( km /s ) mode II itype = 7 - acc. = 8 0 2 4 6 2.5 3 3.5 4 4.5 5 5.5 Distance from the hypocenter (km) R u p tu re v el o ci ty ( km /s ) v P v S v R vS vR forward central does not heal. A direct comparison of these behaviors comes from Figures 10c and 10d of Bizzarri [2010b], where it is plotted the spatial distribution of the rup- ture speed over the whole fault plane. Finally, we men- tion that in extreme cases of heterogeneities in frictional properties, to better constrain the resulting rupture speed we should improve the spatial and temporal dis- cretization of the model. The results presented in Figure 11 confirm that central schemes perform better than forward schemes (at a given order of accuracy). We note that at high or- ders of accuracy (4 and 8) the oscillations of the cen- tral schemes are greater than in the previous case (compare Figures 10 and 11), but they are generally lower than those emerging in forward schemes. More- over, we also confirm here that an increase of the order of the accuracy does not improve the results. Also in this case the five–points stencil method can be regarded as the best candidate to compute the rupture speed. 5. Conclusions In this paper we have considered the problem of the determination of the rupture speed vr characteriz- ing a fully dynamic, spontaneous earthquake. In spon- taneous models the position, as well as the shape, of the rupture front, separating the broken from the un- broken portions of a discontinuity interface, is a part of the solution and it is totally controlled by possible geometric irregularities, rheological properties of the fault structure, heterogeneities, etc. The computation of the rupture speed is performed quite routinely by the modellers of the seismic source. However, a sys- tematic, quantitative exploration of the different nu- merical method to retrieve it from the results of a numerical experiments is currently missed and this is the primary goal that the present study aims to fill. We analyze the simplistic case of a single fault, governed by the linear slip–weakening friction and with homogeneous and heterogeneous rheological proper- ties (see Section 2.2 and 4, respectively), in order to pro- vide some hints for the determination of the rupture times tr and for the computation vr. We also consider a rate– and a sate–dependent friction law (the Ruina–Di- eterich model; see Section 4) to generalize the conclu- sions of the present study. We consider a test case of a supershear event with both 2–D (pure mode II) and truly 3–D (with mixed mode II and mode III, with rake variation allowed) geometries. The results clearly indicate that the way to define tr is very important. In particular, if the fault slip veloc- ity threshold (vl) used to define tr is too high (vl ≥ 0.1 m/s) some details of the rupture process are missed. For an example, in the case of a 2–D, pure mode II rup- ture, for which it is expected to have a rupture jump from sub–Raleigh to supershear speeds [Andrews 1976, Bizzarri et al. 2001], the rupture front bifurcation is lost and the jump of the front is not properly captured (see Figure 7a). This has relevant consequences when one computed the resulting rupture speeds. On the other hand, for vl ≤ 0.01 m/s the results appear to be stable and, more importantly, independent on the choice of vl. We recall that in the framework of the linear slip– weakening governing law we have an exact estimate of the rupture time, which is formally defined as the in- stant when the traction reaches the upper yield stress and this exact estimate has been checked against those resulting from the slip velocity criteria with different threshold values. We can therefore regard the value of vl = 0.01 m/s as the optimal candidate to the used also with other constitutive models, such as the rate– and state–dependent friction laws, for which the rupture time can not be intrinsically determined by the gov- erning model itself. The choice of vl = 0.01 m/s, indicated as the pre- ferred one from the present study, is in agreement with previous papers on dynamic modeling of faults [Day et al. 2005, Bizzarri and Spudich 2008, Bizzarri et al. 2010, Bizzarri 2011]. We also emphasize that in heteroge- neous conditions some patches of the fault can poten- tially produce small stress drop, due to the frictional heterogeneities, and therefore exhibit there relatively small values of fault slip velocity. At the same time we can also expect rapid changes of rupture times on small scales, when the rupture accelerates or it slows down due to the local stress inhomogeneities. As a conse- quence high values of vl can in principle produce some bias in the determination of tr and thus of vr. (Inciden- tally, it is worth to mention that the adoption of vl = 0.1 m/s has been used in other previous papers [e.g., Belardinelli et al. 2003, Antonioli et al. 2006, Bizzarri and Belardinelli 2008] to discriminate between aseismic and seismic regimes of a fault subjected to an external stress perturbation and not to define the rupture front.) Moreover, we show that, in general, central finite differences schemes are more accurate than forward or backward schemes, regardless the order of accuracy. Overall, the most efficient and accurate algorithm is the five–points stencil method at the second–order of accuracy, because it produces small oscillations. Re- markably, these conclusions hold both for simplistic homogeneous and for more realistic heterogeneous configurations, and do not depend on the particular choice of the parameters of the model (expressed in terms of properties of the elastic medium in which the fault plane is embedded, spatial and temporal dis- cretization, constitutive parameters, ratio between the ANDREA BIZZARRI 12 13 spatial grid size and the temporal step size). More im- portantly, these results are general and hold not only for the linear slip–weakening law, but also for different constitutive model, such as the rate– and state–depen- dent friction laws. The time required by the different algorithms used to compute the rupture time (post–processing time) is in all cases negligible with respect to the total CPU time required to solve the elastodynamic problem. We re- mark that these conclusions hold for the finite differ- ence, conventional grid method employed here [Bizzarri and Cocco 2005], and a similar test would be possible also for other numerical algorithms, such as finite elements, spectral elements, discontinuous Galerkin, etc. and this comparison in those cases would be also beneficial. Acknowledgements. We would thank S. Das for intriguing discussions and two anonymous reviewer for their comments which contribute to improve the paper. The editor responsible of this paper was F. Bianco. References Andrews, D.J. (1976). Rupture velocity of plane strain shear cracks, J. Geophys. Res., 81 (32), 5679-5687. Andrews, D.J., and M. Barall (2011). Specifying initial stress for dynamic heterogeneous earthquake source models, Bull. Seismol. Soc. Am., 101 (5), 2408-2417; doi:10.1785/0120110012. Antonioli, A., M.E. Belardinelli, A. Bizzarri and K.S. Vogfjord (2006). Evidence of instantaneous dynamic triggering during the seismic sequence of year 2000 in South Iceland, J. Geophys. Res., 111, B03302; doi: 10.1029/2005JB003935. Archuleta, R.J., and G.A. Frazier (1978). 3–D numerical simulations of dynamic faulting in a half space, Bull. Seismol. Soc. Am., 68, 573-598. Belardinelli, M.E., A. Bizzarri and M. Cocco (2003). Earthquake triggering by static and dynamic stress changes, J. Geophys. Res., 108 (B3), 2135; doi: 10.10 29/2002JB001779. Bizzarri, A., M. Cocco, D.J. Andrews and E. Boschi (2001). Solving the dynamic rupture problem with different numerical approaches and constitutive laws, Geophys. J. Int., 144, 656-678; doi:10.1046/j.13 65-246x.2001.01363.x. Bizzarri, A., and M. Cocco (2005). 3D dynamic simula- tions of spontaneous rupture propagation governed by different constitutive laws with rake rotation al- lowed, Annals of Geophysics, 48 (2), 279-299. Bizzarri, A., and M. Cocco (2006a). A thermal pressur- ization model for the spontaneous dynamic rupture propagation on a three–dimensional fault: 1. Method- ological approach, J. Geophys. Res., 111, B05303; doi:10.1029/2005JB003862. Bizzarri, A., and M. Cocco (2006b). A thermal pressur- ization model for the spontaneous dynamic rupture propagation on a three–dimensional fault: 2. Trac- tion evolution and dynamic parameters, J. Geophys. Res., 111, B05304; doi:10.1029/2005JB003864. Bizzarri, A., and M.E. Belardinelli (2008). Modelling in- stantaneous dynamic triggering in a 3–D fault sys- tem: application to the 2000 June South Iceland seismic sequence, Geophys. J. Int., 173, 906-921; doi: 10.1111/j.1365-246x.2008.03765.x. Bizzarri, A., and P. Spudich (2008). Effects of supershear rupture speed on the high−frequency content of S waves investigated using spontaneous dynamic rup- ture models and isochrone theory, J. Geophys. Res., 113, B05304; doi:10.1029/2007JB005146. Bizzarri, A. (2009). What does control earthquake rup- tures and dynamic faulting? A review of different competing mechanisms, Pure Appl. Geophys., 166 (5-7), 741-776; doi:10.1007/s00024-009-0494-1. Bizzarri, A. (2010a). How to promote earthquake rup- tures: Different nucleation strategies in a dynamic model with slip–weakening friction, Bull. Seismol. Soc. Am., 100 (3), 923-940; doi:10.1785/0120090179. Bizzarri, A. (2010b). On the relations between fracture energy and physical observables in dynamic earth- quake models, J. Geophys. Res., 115, B10307;doi:10. 1029/2009JB007027. Bizzarri, A., E.M. Dunham and P. Spudich (2010). Co- herence of Mach fronts during heterogeneous super- shear earthquake rupture propagation: Simulations and comparison with observations, J. Geophys. Res., 115, B08301; doi:10.1029/2009JB006819. Bizzarri, A. (2011). On the deterministic description of earthquakes, Rev. Geophys., 49, RG3002; doi: 10.10 29/2011RG000356. Bizzarri, A. (2012a). Rupture speed and slip velocity: What can we learn from simulated earthquakes?, Earth Plan. Sci. Lett., 317-318, 196-203; doi:10.1016/ j.epsl.2011.11.023. Bizzarri, A. (2012b). The mechanics of lubricated faults: Insights from 3-D numerical models, J. Geophys. Res., 117, B05304; doi:10.1029/2011JB008929. Bizzarri, A., and S. Das (2012). Mechanics of 3–D shear cracks between Rayleigh and shear wave rupture speeds, Earth Plan. Sci. Lett., 357-358, 397-404; doi:10.1016/j.epsl.2012.09.053. Broberg, K.B. (1989). The near–tip field at high crack velocities, Int. J. Fracture, 39, 1-13. Broberg, K.B. (1999). Cracks and Fracture, Academic Press, New York, 752 pp. Burridge, R. (1973). Admissible speeds for plane–strain self–similar shear cracks with friction but lacking co- CALCULATION OF EARTHQUAKE RUPTURE SPEED hesion, Geophys. J. R. Astron. Soc., 35, 439-455. Das, S., and K. Aki (1977). A numerical study of two–di- mensional spontaneous rupture propagation, Geo- phys. J. Roy. Astr. Soc., 50, 643-668. Day, S.M., L.A. Dalguer, N. Lapusta and Y. Liu (2005). Comparison of finite difference and boundary inte- gral solutions to three–dimensional spontaneous rupture, J. Geophys. Res., 110, B12307; doi:10.1029/ 2005JB003813. Dunham, E.M., and H.S. Bhat (2008). Attenuation of radiated ground motion and stresses from three–di- mensional supershear ruptures, J. Geophys. Res., 113, B08319; doi:10.1029/2007JB005182. Freund, L.B. (1979). Mechanics of dynamic shear crack–propagation, J. Geophys. Res., 84 (B5), 2199- 2209. Guatteri, M., P.M. Mai and G.C. Beroza (2004). A pseudo-dynamic approximation to dynamic rupture models for strong ground motion prediction, Bull. Seismol. Soc. Am., 94 (6), 2051-2063; doi:10.1785/01 20040037. Hok, S., and E. Fukuyama (2011). A new BIEM for rup- ture dynamics in half−space and its application to the 2008 Iwate−Miyagi Nairiku earthquake, Geo- phys. J. Int., 184, 301-324; doi:10.1111/j.1365-246X. 2010.04835.x. Ida, Y. (1972). Cohesive force across the tip of a longi- tudinal−shear crack and Griffith’s specific surface energy, J. Geophys. Res., 77 (20), 3796-3805. Ide, S., and M. Takeo (1997). Determination of consti- tutive relations of fault slip based on seismic wave analysis, J. Geophys. Res., 102, 27379-27391; doi: 10.1029/97JB02675. Lapidus, L., and G.F. Pinder (1999). Numerical solution of partial differential equations in science and engi- neering, John Wiley & Sons, 677 pp. Lapusta, N., and Y. Liu (2009). Three–dimensional boundary integral modeling of spontaneous earth- quake sequences and aseismic slip, J. Geophys. Res., 114, B09303; doi:10.1029/2008JB005934. Lavallée, D., P. Liu and R.J. Archuleta (2006). Stochastic model of heterogeneity in earthquake slip spatial distributions, Geophys. J. Int., 165, 622-640; doi: 10.1111/j.1365-246X.2006.02943.x. Madariaga, R. (1983). High frequency radiation from dynamic earthquake fault models, Annales Geo- physicae, 1, 17-23. Pulido, N., and L.A. Dalguer (2009). Estimation of the high–frequency radiation of the 2000 Tottori ( Japan) earthquake based on a dynamic model of fault rupture: Application to the strong ground mo- tion simulation, Bull. Seismol. Soc. Am., 99 (4), 2305-2322; doi:10.1785/0120080165. Rosakis, A.J., O. Samudrala and D. Coker (1999). Cracks faster than the shear–wave speed, Science, 284, 1337-1340. Rubin, A.M., and J.-P. Ampuero (2005). Earthquake nu- cleation on (aging) rate and state faults, J. Geophys. Res., 110, B11312; doi:10.1029/2005JB003686. Ruina, A.L. (1983), Slip instability and state variable fric- tion laws, J. Geophys. Res., 88, 10359-10370; doi:10. 1029/JB088iB12p10359. Schmedes, J., R.J. Archuleta and D. Lavallée (2010). Cor- relation of earthquake source parameters inferred from dynamic rupture simulations, J. Geophys. Res., 115, B03304; doi:10.1029/2009JB006689. Song, S.G., A. Pitarka and P. Sommerville (2009). Ex- ploring spatial coherence between earthquake source parameters, Bull. Seismol. Soc. Am., 99, 2564-2571; doi:10.1785/0120080197. Song, S.G., and P. Sommerville (2010). Physics-based earthquake source characterization and modeling with geostatistics, Bull. Seismol. Soc. Am., 100 (2), 482-496; doi:10.1785/0120090134. Tinti, E., A. Bizzarri and M. Cocco (2005). Modeling the dynamic rupture propagation on heterogeneous faults with rate– and state–dependent friction, An- nals of Geophysics, 48 (2), 327-345. *Corresponding author: Andrea Bizzarri, Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Bologna, Bologna, Italy; email: bizzarri@bo.ingv.it. © 2013 by the Istituto Nazionale di Geofisica e Vulcanologia. All rights reserved. ANDREA BIZZARRI 14 << /ASCII85EncodePages false /AllowTransparency false /AutoPositionEPSFiles false /AutoRotatePages /None /Binding /Left /CalGrayProfile (Dot Gain 20%) /CalRGBProfile (sRGB IEC61966-2.1) /CalCMYKProfile (U.S. Web Coated \050SWOP\051 v2) /sRGBProfile (sRGB IEC61966-2.1) /CannotEmbedFontPolicy /Warning /CompatibilityLevel 1.3 /CompressObjects /Tags /CompressPages true /ConvertImagesToIndexed true /PassThroughJPEGImages true /CreateJobTicket false /DefaultRenderingIntent /Default /DetectBlends true /DetectCurves 0.1000 /ColorConversionStrategy /LeaveColorUnchanged /DoThumbnails false /EmbedAllFonts true /EmbedOpenType false /ParseICCProfilesInComments true /EmbedJobOptions true /DSCReportingLevel 0 /EmitDSCWarnings false /EndPage -1 /ImageMemory 1048576 /LockDistillerParams true /MaxSubsetPct 100 /Optimize false /OPM 1 /ParseDSCComments true /ParseDSCCommentsForDocInfo true /PreserveCopyPage true /PreserveDICMYKValues true /PreserveEPSInfo true /PreserveFlatness true /PreserveHalftoneInfo false /PreserveOPIComments false /PreserveOverprintSettings true /StartPage 1 /SubsetFonts true /TransferFunctionInfo /Apply /UCRandBGInfo /Preserve /UsePrologue false /ColorSettingsFile (None) /AlwaysEmbed [ true /AndaleMono /Apple-Chancery /Arial-Black /Arial-BoldItalicMT /Arial-BoldMT /Arial-ItalicMT /ArialMT /CapitalsRegular /Charcoal /Chicago /ComicSansMS /ComicSansMS-Bold /Courier /Courier-Bold /CourierNewPS-BoldItalicMT /CourierNewPS-BoldMT /CourierNewPS-ItalicMT /CourierNewPSMT /GadgetRegular /Geneva /Georgia /Georgia-Bold /Georgia-BoldItalic /Georgia-Italic /Helvetica /Helvetica-Bold /HelveticaInserat-Roman /HoeflerText-Black /HoeflerText-BlackItalic /HoeflerText-Italic /HoeflerText-Ornaments /HoeflerText-Regular /Impact /Monaco /NewYork /Palatino-Bold /Palatino-BoldItalic /Palatino-Italic /Palatino-Roman /SandRegular /Skia-Regular /Symbol /TechnoRegular /TextileRegular /Times-Bold /Times-BoldItalic /Times-Italic /Times-Roman /TimesNewRomanPS-BoldItalicMT /TimesNewRomanPS-BoldMT /TimesNewRomanPS-ItalicMT /TimesNewRomanPSMT /Trebuchet-BoldItalic /TrebuchetMS /TrebuchetMS-Bold /TrebuchetMS-Italic /Verdana /Verdana-Bold /Verdana-BoldItalic /Verdana-Italic /Webdings ] /NeverEmbed [ true ] /AntiAliasColorImages false /CropColorImages true /ColorImageMinResolution 150 /ColorImageMinResolutionPolicy /OK /DownsampleColorImages true /ColorImageDownsampleType /Bicubic /ColorImageResolution 300 /ColorImageDepth -1 /ColorImageMinDownsampleDepth 1 /ColorImageDownsampleThreshold 1.10000 /EncodeColorImages true /ColorImageFilter /DCTEncode /AutoFilterColorImages true /ColorImageAutoFilterStrategy /JPEG /ColorACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /ColorImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000ColorACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000ColorImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasGrayImages false /CropGrayImages true /GrayImageMinResolution 150 /GrayImageMinResolutionPolicy /OK /DownsampleGrayImages true /GrayImageDownsampleType /Bicubic /GrayImageResolution 300 /GrayImageDepth -1 /GrayImageMinDownsampleDepth 2 /GrayImageDownsampleThreshold 1.10000 /EncodeGrayImages true /GrayImageFilter /DCTEncode /AutoFilterGrayImages true /GrayImageAutoFilterStrategy /JPEG /GrayACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /GrayImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000GrayACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000GrayImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasMonoImages false /CropMonoImages true /MonoImageMinResolution 1200 /MonoImageMinResolutionPolicy /OK /DownsampleMonoImages true /MonoImageDownsampleType /Bicubic /MonoImageResolution 1200 /MonoImageDepth -1 /MonoImageDownsampleThreshold 1.08250 /EncodeMonoImages true /MonoImageFilter /CCITTFaxEncode /MonoImageDict << /K -1 >> /AllowPSXObjects false /CheckCompliance [ /None ] /PDFX1aCheck false /PDFX3Check false /PDFXCompliantPDFOnly false /PDFXNoTrimBoxError true /PDFXTrimBoxToMediaBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXSetBleedBoxToMediaBox true /PDFXBleedBoxToTrimBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXOutputIntentProfile (None) /PDFXOutputConditionIdentifier () /PDFXOutputCondition () /PDFXRegistryName (http://www.color.org) /PDFXTrapped /Unknown /CreateJDFFile false /SyntheticBoldness 1.000000 /Description << /ENU (Use these settings to create PDF documents with higher image resolution for high quality pre-press printing. The PDF documents can be opened with Acrobat and Reader 5.0 and later. These settings require font embedding.) /JPN /FRA /DEU /PTB /DAN /NLD /ESP /SUO /NOR /SVE /KOR /CHS /CHT /ITA >> >> setdistillerparams << /HWResolution [2400 2400] /PageSize [595.000 842.000] >> setpagedevice