3 D dynamic simulations of spontaneous rupture propagation governed by different constitutive laws with rake rotation allowed

In this work we present a 3D Finite Difference numerical method to model the dynamic spontaneous propagation of an earthquake rupture on planar faults in an elastic half-space. We implement the Traction-at-Split-Nodes fault boundary condition for a system of faults, either vertical or oblique, using different constitutive laws. We can adopt both a slip-weakening law to prescribe the traction evolution within the breakdown zone or rateand state-dependent friction laws, which involve the choice of an evolution relation for the state variable. Our numerical procedure allows the use of oblique and heterogeneous distribution of initial stress and allows the rake rotation. This implies that the two components of slip velocity and total dynamic traction are coupled together to satisfy, in norm, the adopted constitutive law. The simulations presented in this study show that the rupture acceleration to super-shear crack speeds occurs along the direction of the imposed initial stress; the rupture front velocity along the perpendicular direction is slower than that along the pre-stress direction. Depending on the position on the fault plane the orientation of instantaneous total dynamic traction can change with time with respect to the imposed initial stress direction. These temporal rake rotations depend on the amplitude of initial stress and on its distribution on the fault plane. They also depend on the curvature and direction of the rupture front with respect to the imposed initial stress direction: this explains why rake rotations are mostly located near the rupture front and within the cohesive zone. 279 ANNALS OF GEOPHYSICS, VOL. 48, N. 2, April 2005 and arrest as well as the comprehension of the mechanisms controlling slip duration (i.e. healing of slip) and the earthquake energy balance. This progress was possible because of the results achieved in laboratory experiments on fault friction and rock mechanics, in direct observations of earthquake and faulting episodes and in numerical simulations of earthquake dynamics. Stimulated by the collection of high-quality original observations on these phenomena, further investigations are needed to shed light on the numerous remaining open questions. In this context, improving the capability to perform robust numerical simulations of earthquake rupture dynamics represents a fundamental task that has to be pursued to achieve a complete understanding of the earthquake processes. The numerical algorithms widely adopted to simulate earthquake dynamics can be divided in

and arrest as well as the comprehension of the mechanisms controlling slip duration (i.e.healing of slip) and the earthquake energy balance.This progress was possible because of the results achieved in laboratory experiments on fault friction and rock mechanics, in direct observations of earthquake and faulting episodes and in numerical simulations of earthquake dynamics.Stimulated by the collection of high-quality original observations on these phenomena, further investigations are needed to shed light on the numerous remaining open questions.In this context, improving the capability to perform robust numerical simulations of earthquake rupture dynamics represents a fundamental task that has to be pursued to achieve a complete understanding of the earthquake processes.
The numerical algorithms widely adopted to simulate earthquake dynamics can be divided in-

Introduction
The understanding of earthquake rupture propagation and the seismic wave generation process is an important task that has focused scientific research in recent years.Substantial progress has been achieved in this field, allowing the modeling of rupture initiation, propagation Mailing address: Dr. Andrea Bizzarri, Istituto Nazionale di Geofisica e Vulcanologia, Sede di Bologna, Via D. Creti 12, 40128 Bologna, Italy; e-mail: bizzarri@bo.ingv.it to three main groups: Boundary, Domain and Hybrid methods.There are distinct advantages and limitations on each of these numerical approaches and different opinions exist within the scientific community on the choice of the best numerical strategy to perform robust and accurate simulations.Boundary Elements are certainly more accurate than Domain methods; however, they are much less flexible to simulate medium and source heterogeneities.In other words, while in principle Boundary methods are applicable to complex media, they are far from being sufficiently computationally efficient; Domain methods are not as accurate as the Boundary methods in the case of relatively simple structural models, but they are appropriate to simulate complex and more realistic situations.The Domain methods include the Finite Element (FE) and Finite Difference (FD) approaches.The FE can account for complex boundary conditions at the free surface or internal material discontinuity, but the implementation of the complex boundary conditions is not a trivial task in this context.On the other hand, a clear advantage of the FD method is its relative computational efficiency.
Different numerical algorithms have been proposed to solve the elasto-dynamic equation in a 3D continuum medium and to simulate earthquake ruptures (Mikumo and Miyatake, 1978;Day, 1982a,b;Virieux and Madariaga, 1982;Das and Kostrov, 1983;Quin, 1990;Olsen et al., 1997;Fukuyama and Madariaga, 1998, among many others); they both belong to Boundary or Domain methodologies.Some of them have the property, like Finite Element and Finite Difference schemes, that all components of displacement and displacement rate (or (particle) velocity) are defined at the same position in space.With the adoption of such a conventional-grid scheme a fault surface can be represented as by a class of split nodes, represented as disjoint grid points in contact.This approach is called the Traction-at-Split-Nodes (TSN) Fault Boundary Conditions (FBC) and it has been implemented in 2D by Andrews (1973) and in 3D by Day (1977Day ( , 1982a,b),b), Archuleta and Day (1980), Andrews (1999).It will briefly be described in the next section.An alternative is to represent a fault as a yielding distributed through a region of finite thickness; if the independent variables are displacement rate and stress, rather than displacement and displacement rate, yielding can be treated with the Stress-Glut (SG) FBC.Andrews (1999) tested these two methods for a 3D fault without constitutive laws and for a circular crack propagating with an imposed and constant rupture velocity (equal to the Rayleigh wave speed).The conclusion of that paper is that for non-spontaneous models and for adequately refined grid size, both TSN and SG methods converge to the same results.
A main advantage of using the TSN FBC is that it can easily allow the introduction of a constitutive law relating the total dynamic traction to fault friction.The possibility to compute fault slip, fault slip velocity and traction at the same grid point allows the specification of either slipor rate-and state-dependent constitutive laws, as we will discuss in the following sections.
Another important assumption commonly required in numerical procedures is the choice of the amplitude and the direction of the initial stress on the fault plane.Many of the existing numerical procedures assume, for simplicity, that the initial stress is assigned only in one direction (Olsen et al., 1997;Aochi et al., 2000); in this way slip becomes automatically a scalar and only one component is non-null.Other papers (Fukuyama and Madariaga, 1998;Fukuyama et al., 2003) consider the friction law in a vectorial formulation, in the sense that the constitutive law is written independently for each components of the solution, but also in these models the initial stress is non-null only in one direction.
In this study we present and discuss the main features and several numerical simulations performed with a fully (namely truly) 3D spontaneous dynamic rupture model in which the crack propagation and the dynamic traction evolution are governed by assigned constitutive laws.We use the TSN FBC, which assumes an explicit displacement discontinuity between the two-halves of the fault, and the initial stress is a vector on the fault plane.We therefore allow both components of slip to be non zero; they are mutually coupled by the assumed constitutive relation.

The numerical methodology
The numerical algorithm presented in this paper is based on a 3D Finite Difference (FD) scheme written by J. Andrews.A grid of nodes is introduced in the 3D space and the solution of the fundamental elasto-dynamic equation is obtained by stepping the main physical quantities explicitely through time by calculating the net force (formally the restoring force) acting at all nodes.The discretization is made by using the quadrilateral isoparametric elements (Hughes, 1987) with all edges parallel to the axes of the Cartesian coordinate system; stress is not uniform inside an element and the fundamental physical variables are displacement and force at nodes.All components of displacement, displacement rate, acceleration and nodal force are defined in the same grid point; components of displacement, acceleration and force are known at integer time levels, while those of displacement rate are known at half-integer time levels.From the displacement components at each node of the medium we compute the strain tensor components at integration points and then stress tensor (using the Hooke's Law for an isotropic medium with the small displacement approximation adopted) and net force components at nodes.Local forces are calculated using the 8-points Lobatto integration.Displacement, displacement rate and the net force are iterated explicitly through time.Particle acceleration is derived from the net force acting on each node, simply applying the second law of mechanics.Finite differences in space are formulated to be equivalent to finite elements and therefore the numerical algorithm can be considered either as a Finite Element or as a Finite Difference scheme.This formulation is mathematically equivalent to the local stiffness matrix, but it is more efficient.
In fact the number of finite difference equations is 8 times greater than in a simple staggeredgrid scheme, in which the different components of the fundamental variables (displacement rate and stress) are defined in different points in space.In the latter scheme it is very easy to increase the order of accuracy to fourth-order in space, while keeping it second-order in time.
There is no doubt that the fourth-order staggeredgrid method is the most efficient and accurate scheme to propagate seismic waves through the medium outside the fault plane.However, accurate implementation of FBC is not simple with a fourth-order staggered-grid scheme.In fact, in or- (stresses are assumed to be negative for compression: σ n and pfluid are positive numbers).A particle located on the negative side is therefore affected by a total normal traction Σ Σ n + Σ Σ p fluid .Following the Terzaghi effective stress law (Terzaghi et al., 1996;Wang, 2000), the fault friction τ is proportional to the modulus σ n eff of the effective normal traction Σ Σ = Σ Σ n − Σ Σ p fluid = − σ n eff n ˆ= − (σn − pfluid) n ˆ.In this paper we neglect pore fluid pressure changes and therefore we consider constant normal stress.In the corner of each panel are represented the fundamental building elements: a equilateral triangle in 2D case and a parallelepiped in 3D one.a b der to simulate earthquake ruptures the FBC is more important than propagation of seismic waves to large distances.Our numerical simulations show that this method is about 12 times more expensive in time machine than the staggered-grid scheme of Madariaga (1976).
The numerical algorithm presented here has second-order accuracy in space and second-order accuracy in time, except at the starting and at the healing of slip, where it is first-order.The geometric configuration is reported in fig.1b for a single vertical fault plane.Our numerical algorithm allows the simulation of several faults, including dip-slip faulting mechanisms.Solutions depend on two spatial coordinates (along strike and dip directions) and time [u i(x1, x3, t), u being the fault slip (i.e. the displacement discontinuity)] and it has two non-null components (u 1 and u3); in this way we have a fully (namely truly) 3D fault model.The medium is parameterized with rectangular box elements in the x1-x3 plane (suitable for vertical faults, see fig.1b) and with right prisms having rectangular cross section on the x 1-x3 plane (suitable for dip-slip faults).The x3 = 0 plane is the free surface.
All outputs of the code should be regarded as average over an area centered on the point that extends a half-grid interval in each space dimension and extends in time a half time step before and after the nominal time of the output.The convergence and stability conditions of the proposed numerical procedure are discussed in detail in the Appendix.

The fault constitutive relations
In the case of constant sliding velocity on a surface, we can define the coefficient of friction µ as the ratio between the lateral (shear) friction force and the externally applied normal force (e.g., Scholz, 1990).From the first Amonton's Law it is also the ratio of the modulus τ of the shear traction T (T = τ T ˆ, where T ˆis the unit vector expressing the local direction of the shear traction on the surface) and the modulus σ n eff of normal traction Σ Σ (Σ Σ = − σ n eff n ˆ, n ˆ being the unit vector normal to the surface).In our geometry the fault is a plane perpendicular to the x2 axis of the Cartesian coordinate system Ox 1 x 2 x 3 (see fig. 1a,b) and therefore n ˆis collinear to x ˆ2 ≡ ≡ (0,1,0), i.e. the second vector of the canonical orto-normal basis.Indicating with symbol σ22 the normal component of stress tensor, follows that Σ Σ = σ22n ˆ= − σ n eff n ˆ.In this paper stresses are assumed to be negative for compression.Moreover, if we denote with σ21 and σ23 the shear components of the stress tensor applied on the positive side of the fault surface, we can write the fault shear traction modulus as In full of generality a governing law is the explicitation of the analytical dependencies of the fault friction τ (also called fault strength) on some physical observables , , , , , , , , , where u and v are the fault slip (i.e. the displacement discontinuity) and the fault slip velocity modula, respectively; Ψ Ψ = (Ψ 1, …, ΨN) is the state variable vector (Ruina, 1983); T is the temperature, accounting for ductility, plastic flow, rock melting and vaporization; H is the humidity (Dieterich and Conrad, 1984); λc is the characteristic length of the fault surface, accounting for roughness and topography of asperity contacts (Ohnaka and Shen, 1999;Ohnaka, 2003) and eventually responsible of mechanical lubrication (Brodsky and Kanamori, 2001); h is the material hardness; g is the gouge, accounting for surface consumption and gouge formation during sliding episodes (e.g.Marone et al., 1990;Marone and Kilgore, 1993;Mair and Marone, 1999;Mair et al., 2002); Ce is the chemical environment; σ n eff is the value of the effective normal stress, that in general may varies in time (Linker and Dieterich, 1992;Prakash, 1998) and is expressed as the difference existing between the value of the reference normal stress σn (due to the regional tectonic loading) and value of the pore fluid pressure pfluid (Andrews, 2002;Bizzarri andCocco, 2004, 2005a,b).
In this study we implement either a Slip-Weakening (SW) constitutive law (Ida, 1972;Palmer and Rice, 1973;Andrews, 1976a,b) or a Rate-and State-dependent (RS) friction formulation (Dieterich, 1979a,b;Beeler et al., 1994;Bizzarri and Cocco, 2003 and references therein) under the assumption of constant normal stress (i.e.σ n eff = σn).The fault friction is updated through time in any position during rupture propagation in a way that depends on the adopted constitutive relation.In all cases we assume that total shear traction is collinear with the fault slip velocity at each time step at any point on the fault during dynamic rupture propagation (i.e.T ˆ= v/v).

The slip-weakening law
We use in this study the linear SW law in the linear, or classical, formulation proposed by Ida (1972) and Andrews (1976a,b), in which the dynamic traction (τ) degrades linearly for increasing slip (u) from the maximum yield stress (τu) to the kinetic frictional level (τf) over a characteristic slip distance (d0) The decrease of traction for increasing slip (the slip-weakening behaviour) has been experimentally observed by Ohnaka and Yamashita (1989) and eq.( 3.3) has been often adopted in numerical modeling of earthquake ruptures (Andrews, 1994;Olsen et al., 1997;Aochi et al., 2000, among many others).Recently, several papers have attempted to infer the critical slip-weakening distance from real earthquakes (Ide and Takeo, 1997;Guatteri and Spudich, 2000;Mikumo et al., 2003): these studies propose values of d0 of the order of 0.1-1 m.According to this constitutive law fault strength depends only on slip.
As stated before, we recall that τ in eq. ( 3.3) is intended as the modulus of the dynamic shear traction vector T. Because we assume that dynamic traction is collinear with slip velocity at each time step, we derive the direction cosines from slip velocity (see eq. ( 5) in Andrews, 1999) and use them to compute the two components of the dynamic traction.Slip velocity and dynamic traction have been updated through time following the numerical strategy described in Andrews (1999).In other words, in our approach we associate the fault constitutive relation to the modulus τ of the total dynamic traction T and from that we compute the two traction components from the direction cosines obtained from fault slip velocity.

The rate-and state-dependent friction laws
Rate-and state-dependent friction laws assume that the fault friction depends on the fault slip velocity v and on one state variable Ψ (Dieterich, 1979a,b).In this study we use the slowness friction law (see Beeler et al., 1994), which is defined as follows: eff and L are the constitutive parameters; the former accounts for the direct effect of friction and the latter is the characteristic length over which the state variable evolves.v * and µ * are reference values, v is the modulus of slip velocity and σ n eff is the value of the effective normal stress that in this paper is assumed to be constant.In Bizzarri and Cocco (2005a,b) we have generalized the method in order to solve the problem with temporally varying effective normal stress.The second equation is currently named the evolution law.Cocco and Bizzarri (2002) have shown that RS constitutive laws yield the slip-weakening behavior; Bizzarri and Cocco (2003) have associated SW and RS frictional parameters and extensively discuss the relationship existing between the characteristic lengths of the two different constitutive laws.
Also in this case the constitutive law defined in eq.(3.4) is associated to the modulus τ of total dynamic traction T; it defines the frictional behavior as a function of the state variable and the modulus of slip velocity.We use a 3D generalization of the Rosembrok Stiff Integration (Press et al., 1992), previously adopted in our 2D fault models (Bizzarri et al., 2001).This in-tegration method performs better than the Runge-Kutta algorithm also in 3D case.We solve the following set of differential equations to get the state variable and the two independent components of fault slip velocity (v 1 and v3): where ui are the components of the slip, α is the differential acceleration for unit stress drop and Li are the component of the dynamic load that will give no change of differential acceleration.We therefore calculate the director cosines from slip velocity and, assuming that traction is collinear to slip rate, we obtain the two components of total dynamic traction.
As already stated above, our numerical procedure allows the calculation of the two components of slip, slip velocity and total dynamic traction at the same positions on the fault plane for each time step.The two independent components of the relevant physical quantities are coupled through the adopted constitutive law.In this study, the adopted constitutive law is considered an analytical relation between the modulus of traction and that of slip or slip velocity (and the scalar state variable).This is different from the assumption made by other authors (Fukuyama and Madariaga, 1998, for instance), who specify distinct constitutive relations independently of each component of slip, and then consider for simplicity only one non-null components of initial stress.In this study we use an initial stress defined as a vector with two non-null components, and our numerical strategy guarantees the possibility to perform general 3D simulations which allows the rake to rotate in time.

Results with the slip-weakening constitutive equation
In this section we compare different numerical simulations performed using our 2D and 3D FD codes, both based on the TSN FBC.The geometry of these configurations is represented in fig.1a,b (fig.1a for the 2D case and fig.1b for the 3D one).Domain boundary conditions are the same and the initial stress in the 3D model is oriented along the x 1 axis (rake 0°); in this case the initial component of the shear stress is non zero only along the x1 axis.The nucleation strategy is in the same for both the configurations: we impose an initial time-weakening friction, until the rupture begin to propagate spontaneously (for further details see Bizzarri et al., 2001).Medium and constitutive parameters are listed in table I.In fig.2a we plot the slip obtained from the 2D numerical simulation as a function of time and the spatial coordinate along the fault.Figure 2b,c displays the total slip pattern obtained from the 3D numerical simulation: fig.2b shows the total slip as a function of the x 3 coordinate (taking x1 equal to the coordinate of the nucleation point), while fig.2c shows the slip as a function of x1 (fixing x3 as the coordinate of the nucleation point).Figure 2b describes the anti-plane mode of propagation, while fig.2c identifies the in-plane mode.

Parameter
Value As the strength parameter S (Das and Aki, 1977a,b) is lower than the critical value of 1.77, the crack accelerates to super-shear velocity and it bifurcates.This behavior was first analytically demonstrated by Burridge (1973) and subsequently numerically modeled by Andrews (1976b) for a 2D in-plane spontaneous dynamic rupture.The crack tip bifurcation might not be representative of most real-world observations, that show the prevalence of the sub-shear ruptures in earthquakes.However, there are some laboratory observations (e.g., Rosakis et al., 1999) in which a super-shear crack tip velocity was measured.Xia et al. (2004) enumerates a series of earthquakes that exhibit supershear rupture speeds (1979 Imperial Valley EQ; 1992 Landers EQ; 2002 Denali EQ; Kunlunshan EQ, 14 November 2001).Additionally, the interested reader can find a discussion in a recent paper by Bouchon et al. (2001).Although the crack tip bifurcation has already been modeled in 2D numerical simulations of spontaneous dynamic rupture governed by different governing equations (see for instance Bizzarri et al., 2001), there are no truly 3D numerical experiments showing this behavior with rateand state-dependent constitutive equations.Moreover, it is important to emphasize that we have chosen this set of parameters and modeled super-shear rupture models because rake rota-tions are much more pronounced in these cases (as we will discuss in Section 5).In fig.2a we clearly distinguish two rupture fronts: the first traveling at Rayleigh-wave velocity and the second one accelerating from S-to P-wave speed.Also in the 3D simulation we observe the crack tip bifurcation, but it occurs only in one component (in-plane) and the distinction between the two rupture fronts is less evident.The transition to super shear rupture velocity is more gradual in the 3D simulations, as shown in fig.2c.Moreover, the 2D simulations show that the cohesive zone thickness at the crack bifurcation increases; this behavior is not confirmed in our 3D simulations and the thickness of the cohesive zone is different between the two components plotted in fig.2a-c.
In fig.3a,b we plot the solutions for the 2D and for the 3D models in a fault point located on the x1 axis 18 spatial steps from the nucleation point.Figure 3a shows the traction versus slip behavior, representing the imposed constitutive law, while in fig.3b we plot the phase diagram (i.e. the traction as a function of the slip velocity).It is important to remark here that in 3D simulations we plot the modulus of the vectors.The phase diagram shown in fig.3b for the 3D fault is more similar to that obtained with rate-and state-dependent friction laws (see next section and Bizzarri et al., 2001).The temporal evolution of traction, fault slip and slip velocity are shown in fig.4a,b, for both the 2D model (fig.4a) and for the 3D one (fig.4b).We emphasize that, for this particular value of the strength parameter S, the peak of slip velocity coincides with the minimum of total traction, in agreement with Bizzarri (2003) and Tinti et al. (2004).For this particular configuration the value of the characteristic slip-weakening distance (d0) can be inferred form the slip at the time of slip velocity peak, as proposed by Mikumo et al. (2003).
We also observe that for 3D simulations the slip acceleration and the weakening phase in traction evolution are shorter than that in 2D; this is due to the difference in the dynamic load caused by the fault points that are already slipping.Fig. 5a-c.The same as fig.2a-c, but adopting the Dieterich-Ruina governing equation.All parameters are listed in table II.Also with the rate-and state-dependent friction law we observe that the crack tip accelerate asymptotically up to the P-wave velocity.(3.4)).Medium and constitutive parameters are listed in table II.For particular configurations, the dynamic propagation of a rupture obeying to RS friction laws is also characterized by a crack tip bifurcation, with a secondary front that accelerate up to the P-wave velocity (Okubo, 1989;Bizzarri et al., 2001).This is represented in fig.5a for the 2D fault model and in fig.5b,c for the 3D simulation.In the figures the black arrow identify the point at which such a bifurcation occurs.We emphasize that in the case of Dieterich-Ruina Law the transition from sub-to super-shear velocity is gradual in both fault model.These calculations confirm the results of the comparison between 2D and 3D simulations with the SW law, pointing out that the rupture propagation along the direction of imposed initial traction (in-plane in the simulations here discussed) is more similar to the 2D solution than the rupture propagation in the anti-plane direction (shown in figs.3b and 5b).

Results with Dieterich-Ruina Law
As done in fig.3a,b, we plot in fig.6a,b the behavior of dynamic traction as a function of fault slip (fig.6a) and slip velocity (fig.6b).Both the slip-weakening curve and the phase diagram obtained with the rate-and state-dependent law defined in eq.(3.4) does not change substantially between 2D and 3D calculations.The equivalent SW distance d 0 eq (Cocco and Bizzarri, 2002) is the same for 2D and 3D simulations.These results corroborate and generalize the conclusion of Bizzarri and Cocco (2003), who used a 2D fault model obeying to RS constitutive laws.
In fig.7a,b we show the temporal evolution of dynamic traction, slip and slip velocity, similarly to what illustrated in fig.4a,b for the SW model.These calculations corroborate and generalize the results of Tinti et al. (2004), obtained with 2D fault model, according to which for a crack obeying to rate-and state-dependent friction laws the maximum peak slip velocity does not correspond in time to the minimum value of the traction; in this case if we measure the equivalent slip-weakening distance d 0 eq from the slip at the time of slip velocity we can underestimate it up to the 50%.

Temporal rake rotations during dynamic rupture propagation
In this section we present the results of several simulations performed by assuming an oblique initial stress direction.The sketch depicted in fig.8 illustrates a general configuration of a 3D simulation in which the crack is propagating on the fault plane and, at a selected time step, the rupture propagation direction does not coincide with the slip velocity orientation at the same time step.The direction of crack growth at any point on the fault plane (solid blue vector in the figure) determines the local in-plane and anti-plane components of motion, which in general do not coincide with the strike-and dip-slip components (#1 and 3 in fig.8 and shown with solid green vectors) of sliding.Slip and slip velocity vectors are collinear with dynamic traction, because this is imposed in the solution of the elasto-dynamic equation.
We have performed several 3D dynamic simulations using the parameters listed in table III and assuming an initial stress oriented at 45°on the fault plane.We solve the elasto-dynamic equation using both a slip-weakening and a rateand state-dependent friction law.In fig.9a-f we present the results of the simulation performed by using the SW law: we show the slip velocity and the total dynamic traction evolution on the fault plane at three different time steps.Figure 10a-f shows a similar calculation obtained using RS and a slowness evolution equation.As expected, these simulations show that the rupture propa- gates faster in the direction of the initial stress (45°for these calculations) and slip velocity is largest in this direction.The two constitutive laws yield similar results and the main difference concerns the nucleation duration and therefore the dynamic propagation is not synchronous.In both cases we observe that the crack tip bifurcation (with a primary front traveling at the Rayleigh wave velocity and a secondary one traveling at    the P-wave velocity) and the acceleration to super-shear rupture velocity occurs only in the local in-plane direction, that corresponds to the direction of the imposed initial traction (see also Bizzarri, 2003).We have plotted the rake variations on the fault plane in figs.11a and 12a for the two constitutive behaviors here considered; the rake variation is measured with respect to the initial orientation of the initial stress, that for these simulations is taken equal to 45°.These figures represent the changes in slip velocity directions (which correspond to the orientation of the instantaneous dynamic traction) in the positive quadrant of the fault plane x 1-x3 at a fixed time step.Figure 11a-d shows the results obtained with the slip-weakening law, while fig.12a,b shows similar calculations for the Dieterich-Ruina (ageing) Law.These simulations show that within the broken region (the area inside the crack tip) there is no spatial variation of the rake in both cases, which means that the dynamic traction is always collinear with the initial stress.The rake rotation occurs mostly within the cohesive zone and in correspondence of the rupture front bifurcation.We emphasize that this behavior is general and independent of the particular value of the initial stress orientation (Bizzarri, 2003).We plot in fig.11c,d the corresponding cases shown in fig.11a,b, but assuming an initial pre-stress directed along the x 1 axis (i.e. initial rake of zero degrees).We can observe that the behavior is quite similar.We would emphasize, however, that simulations performed with an initial traction having an arbitrary orientation with respect to the x1 axis are not simply a Cartesian rotation on the fault plane of the solutions obtained for a case in which the initial stress is directed along the x 1 axis, because of the presence of the free surface.In figs.11b and 12b we plot the time history of the rake (formally the dynamic traction orientation at different time steps) for three positions indicated in panels (a) of each figure by the stars.These fault positions are located at a distance of 2 and 3 km from the nucleation point in figs.11a-d and 12a,b, respectively.For both the figures, point #1 is located along the x1 axis, #2 is located in the direction of the initial stress and #3 lies on the x3 axis.The vertical grey lines in panels (b) indicate the duration of the breakdown process.We emphasize that we do not have any temporal variation on the rake in position #2: this confirms the previous findings that along the direction of the initial stress the slip velocity does not change its orientation.This position corresponds to a case in which the rupture front is perpendicular to the initial stress direction (or equivalently the local rupture propagation direction is collinear with the initial stress, see fig. 8).We can clearly observe that the temporal variation of the rake is concentrated in the cohesive zone, but also before a variation occurs (see locations #1 and #3).In particular, we observe for location #1 that the compo-nent 3 decreases before the breakdown process, then it begins to increase inside the cohesive zone and reaches its maximum in correspondence of the end of breakdown process; it therefore decreases again and reaches a final value that is equal to the initial one.This is a general behavior, numerically modeled both with the slip-weakening law (fig.10b) and with the Dieterich-Ruina constitutive equation (fig.12b).We emphasize, however, that the variation of the rake is more evident in the case of the slipweakening law.

Dependence on the absolute level of stresses
Andrews (1994) simulated the mixed-mode dynamic propagation of an earthquake rupture using a Boundary Integral Equation method in which the crack tip is an infinite line (not closed) propagating in the x1 direction; the solutions of the dynamic problem depend only on the x1 coordinate but have two non-null components (1 and 3).Simulations were done using an initial stress oriented at 45°.Although in these calculations the cohesive zone is not well resolved (only three points within the breakdown duration), he found that the two components of slip velocity (in-plane and anti-plane) interact with each other and the rake changes with time; he sug-Fig.12a,b.The same as fig.11a-d, but assuming a Dieterich-Ruina Law.Vertical lines in panel b) have to be regarded as the time step at wich are reached the equivalent yield and kinetic levels (for a detailed discussion about the correspondence of frictional parameters of different governing laws see Cocco andBizzarri, 2002 andBizzarri andCocco, 2003).gested that the variation of slip direction decreases with increasing the absolute initial stress level.Our simulations, performed with a better resolution of the cohesive zone, confirm the results of Andrews (1994) and generalize them because we consider a more complex coupling between the two local components of slip velocity, having a fully 3D dynamic rupture propagation and we consider also what happens before the rutpure onset.We performed several simulations using the initial and constitutive parameters listed in table IV, an initial stress direction oriented at 45°(as in previous calculations and in Andrews, 1994) and a strength parameter S always equal to 0.8.The strength excess and the dynamic stress drop is the same among the four simulations performed with the parameters listed in table IV.We measured the rake variation as a function of time in a fault position located on the x 1 axis and at a distance of 18 steps form the nucleation point.We show in fig.13 the rake variations as a function of time.This figure clearly shows that the amount of rake variation depends on the absolute values of stress parameters: the change in direction of slip velocity is larger for smaller values of initial stress.In particular, we observe that the rake variation decreases as the ratio of the breakdown stress drop (defined as the difference between the upper yield stress and the kinetic frictional level) and the initial stress decreases.The same direct proportionality exists between the ratio of the dynamic stress drop (defined as the difference between the initial stress and the kinetic frictional level) and the initial stress.This is in agreement with the results of Spudich (1992) and with Guatteri and Spudich (1998).We point out, however, Fig. 13.Temporal rake evolution (in a fault point located at 18 units from the nucleation point and along the x1 axis) during the propagation of different cases where the slip-weakening law has been adopted.For each color, representing a different numerical simulation, the vertical lines marks the start and the end of the breakdown process.The parameters are listed in table IV. that in our simulations the variation of these stress ratios are determined by variations of the initial stress amplitude.Increasing the initial stress amplitude by a factor 10 reduces the rake rotation from 44% to 11% of the initial direction.For all these simulations the largest rake change occurs within the cohesive zone, in agreement with previous results.

Heterogeneous distribution of pre-stress on the fault plane
All the simulations presented and discussed in previous sections were performed in a homogeneous configuration: all medium and constitutive parameters are uniform.In this section we will present a result obtained by imposing a SW law and a distribution of initial stress highly heterogeneous (fig.14), characterized by a high initial stress patch in the central portion of the fault plane that corresponds to the nucleation patch.Such a distribution is obtained by interpolating a slip distribution of the 1999 ChiChi, Taiwan, earthquake proposed by Yagi (http://www.eic.eri.u-tokyo.ac.jp/yuji/taiwan/ aiwan.html).In this case the initial stress is directed only in the x1 direction.The other constitutive and model parameters are listed in table V. Figure 15a-f shows the 3D snapshots of slip and slip velocity calculated at three distinct time intervals.This figure points out that both the slip and the slip velocity distributions on the fault plane are heterogeneous.Rupture is arrested by areas of very low initial stress and accelerates in areas of high pre-stress.We plot in fig.16 the rake variation, with respect to the initial value of 0°, on the fault plane at a fixed time.We can observe that the main variations of the slip direction are concentrated in the rupture front, and near the barrier around 6875 m, where the rupture propagation is inhibited due to a very low value of the initial stress.We conclude that both the heterogeneity of initial stress and its absolute amplitudes control the complexity of the crack propagation causing local rotations or deviations of the slip direction from the initial stress orientation.

Discussions and conclusive remarks
In this study we have presented a fully (namely truly) 3D Finite Difference method to solve the elasto-dynamic equation in order to simulate the spontaneous rupture propagation of an earthquake rupture obeying to different con- stitutive laws.We use the Traction-at-Split-Nodes technique which allows the evaluation of slip, slip velocity and total dynamic traction at the same position on the fault plane.We use both a slip-weakening and a rate-and state-dependent constitutive laws to model the dynamic rupture propagation on an assumed fault plane with an oblique orientation or a heterogeneous distribution of the initial stress.Our numerical procedure allows the adoption of heterogeneous distribution of constitutive parameters on the fault plane.We will not present here the results of simulations with heterogeneous distribution of frictional; we will discuss in Tinti et al. (2005) for a 2D fault model.
The 3D simulations of dynamic rupture propagation performed in this study by using either a slip-weakening or a rate-and state-dependent constitutive law confirm the results presented in previous investigations which used 2D dynamic solutions of the elasto-dynamic equation: in particular, the scaling relation between the two length-scale parameters of the adopted constitu-   tive formulation (d0 and L, see Bizzarri and Cocco, 2003, and references therein) as well as the temporal evolution of dynamic traction and slip velocity (Tinti et al., 2004).However, our simulations generalize these results and point out that the rupture front acceleration to super-shear speeds occurs along the imposed direction on initial stress and the crack velocity is slower in the perpendicular direction.Simulations performed with an oblique direction of initial stress confirm this conclusion.Depending on the position on the fault plane the direction of instantaneous dynamic traction can be different from the imposed direction of pre-stress.Because we assume that total dynamic traction is collinear with slip velocity, this variation of the instantaneous orientation is associated to a rake change.These rake rotations are localized in the proximity of the rupture front and within the cohesive zone.
Our calculations show that the amount of rake rotation depends on the absolute value of initial stress: for a constant dynamic or breakdown stress drop, larger amplitudes of initial stress yield smaller rotations of instantaneous dynamic traction around the initial imposed direction.Thus, in general we expect that rake rotations might depend on the ratio between breakdown (or dynamic) stress drop and initial stress.Preliminary simulations with heterogeneous pre-stress on the fault plane suggest that Appendix.Convergence and stability conditions.
In this Appendix we summarize the convergence and stability relations that have to be satisfied in order to correctly resolve the dynamic traction and the slip velocity evolution during the dynamic processes.
The first condition that has to be satisfied is the temporal and spatial resolution of the cohesive zone where total friction decreases from the maximum yield value to the kinetic level.For a 2D case the requirement of resolution of this characteristic time duration and spatial scale consists to impose where T b (τ 1 ) and T b (τ 3 ) are the breakdown zone times calculated for the time evolution of the components 1 and 3 of traction, respectively.We refer to a vertical fault, as depicted in fig.1b.The break-in non-uniform conditions the temporal rake rotation might be more diffuse on the fault.This study aimed to show the applicability and the flexibility of our 3D Finite Difference numerical approach to simulate the nucleation and the dynamic propagation of an earthquake rupture as well as to investigate the constitutive behavior.The possibility to include the heterogeneity of constitutive parameters allows the modeling of the healing of slip and the dynamic traction evolution.The adoption of the Traction-at-Split-Nodes technique allows further implementations aimed to include the effects of variable normal stress perturbations and porefluid effects on the friction coefficient.down zone times in a fixed fault point (x 1 * , x 2 * ) can be expressed as (see eq. (A.5) and its analytical derivation in Bizzarri et al., 2001)  where φ is the rake calculated in the actual fault point.These relations are true both for the slip-weakening law (eq.(3.3)) and for the rate-and state-dependent friction laws (eq.(3.4)).In the latter case quantities Tb and X b have to be regarded as equivalent breakdown zone time and as equivalent breakdown zone extension, respectively, in the sense discussed in Cocco and Bizzarri (2002) and in Bizzarri and Cocco (2003).
The second condition that have to be considered states that the spatial and time steps are coupled by the general condition (Andrews, 1985;Fukuyama and Madariaga, 1998;Bizzarri et al., 2001, among many different others) where d is the maximum distance between two neighboring nodes [d = (∆x 1 2 + ∆x 2 2 + ∆x 3 2 ) 1/2 ].In other words eq.(A.4) requires that no coupling exists between first neighbors.This condition is common to Boundary Integral Equation methods and to Finite Difference approaches (Bizzarri et al., 2001).
In the case of rate-and state-dependent friction laws there is an additional requirement, introduced by Rice (1993), that in full of generality it can be expressed as k diag >> k cr , where k diag is the diagonal term of stiffness matrix and kcr is the critical stiffness.The condition kdiag >> kcr corresponds to impose that locally each single element of the discretized fault is conditionally stable (Scholz, 1990).This avoids that a single node may fail independently of the neighbors (numerical noise and artificial complexity) and guarantees that the discrete medium can be considered as a continuum.
In our 3D fault model, the local stiffness is expressed as kdiag = 1/C = (ρ∆x 2 )/(4∆t 2 ), where C is the local compliance (Andrews, 1985;Bizzarri et al., 2001;Bizzarri and Cocco, 2003) and represents the proportionality constant between instantaneous traction and dynamic slip.The critical stiffness can be expressed as (b − a)σ n eff /L (Ranjith and Rice, 1999), where the constitutive parameters a, b, σ n eff and L have been assigned.When k diag = k cr , we have the critical grid size because the time step and the spatial grid size are related by the Courant-Friedrichs-Levy (CFL) ratio w CFL (w CFL = v S ∆t/∆x, see Fukuyama and Madariaga, 1998;Bizzarri et al., 2001).The Rice's condition can be therefore expressed as or, alternatively, as .6)

Fig
Fig. 1a,b.Schematic representation of geometry of the 2D fault model (a) and of the 3D one (b).The stars represent the points of the nucleation, blue arrows indicate the local crack enlargement direction, while red arrows indicate the direction of propagation.In the 2D case (a) the rupture propagates only in the x1 direction and all solutions are independent on the x3 spatial coordinate.In both cases (2D and 3D) particles placed on the positive side of the fault plane, oriented through the unit vector n ˆ// x ˆ2 ≡ ≡ (0,1,0), exert on the particles on the negative side a normal traction Σ Σ n = − σ n n ˆand a pore fluid pressure Σ Σ p fluid = − pfluidn(stresses are assumed to be negative for compression: σ n and pfluid are positive numbers).A particle located on the negative side is therefore affected by a total normal traction Σ Σ n + Σ Σ p fluid .Following the Terzaghi effective stress law(Terzaghi et al., 1996;Wang, 2000), the fault friction τ is proportional to the modulus σ n eff of the effective normal traction Σ Σ = Σ Σ n − Σ Σ p fluid = − σ n eff n ˆ= − (σn − pfluid) n ˆ.In this paper we neglect pore fluid pressure changes and therefore we consider constant normal stress.In the corner of each panel are represented the fundamental building elements: a equilateral triangle in 2D case and a parallelepiped in 3D one.
Fig.2a-c.Slip obtained for simulations adopting the linear slip-weakening law for a 2D fault (a) and for a 3D fault (b and c).All medium and constitutive parameters are listed in table I. Arrows identify the point where the bifurcation occurs: a primary front propagate at the R-wave velocity and a secondary one accelerate from S-to P-wave speed.The cohesive zone (where the traction degrades from the maximum yield stress τu down to the kinetic frictional level τf) is also indicated by a contouring in all plots.

Figure
Figure 5a-c illustrates the comparison between 2D and 3D simulations performed by us-

Fig
Fig. 3a,b.Solutions at the fault point located at 18 units form the nucleation point, along the x1 axis: a) traction versus slip; b) phase portrait (or phase diagram), i.e. traction versus slip velocity.Solid square are obtained form the 2D fault calculation, while open circles are obtained from the fully 3D calculation.All parameters are the same as fig.2a-c.
1.5 m) ∆x1 = ∆x2 = ∆x3 0.01 m Initial rake 0°F ig.6a,b.The same as fig.3a,b, but for the case shown in fig.5a-c.Solutions are plotted at a fault point located at 2.0 m from the nucleation point.

Fig
Fig. 7a,b.The same as fig.4a,b, but for the Dieterich-Ruina Law.All parameters are those of figs.5a-c and 6a,b.The fault point is the same as fig.6a,b.The vertical grey line marks the time step at which the friction is at its final level (the equivalent frictional level); this instant does not coincide with the time step at which the slip velocity is at its maximum value.See the text for a detailed discussion of physical implications.

Fig. 8 .
Fig. 8. Sketch representing a time snapshot of the positive quadrant (with respect to the nucleation point) of the fully 3D fault.The star indicates the nucleation point, as in fig.1a,b; the solid red line identifies the instantaneous position and shape of the crack tip, that defines the contour of the broken area (gray region).In a generic point of the tip a blue arrow indicates the local crack enlargement direction, that in general does not coincide with the local slip velocity direction (red arrow), that is collinear with the dynamic traction.Black arrows represent the local in-plane and the local antiplane modes of propagation, that in our fault model are coupled through the constitutive law.These modes are not always coincident with the components along the x1 and x3 axis.
Table III.Medium and constitutive parameters for 3D numerical simulations shown in figs.9a-f, 10a-f, 11a-d and 12a,b.Configuration for the slip-weakening law (figs.9a-f and 11a-d) refers to a real-world fault, while input parameters for the Dieterich-Ruina Law (figs.10a-f and 12a,b) are typical of a laborato-

Fig.
Fig.10a-f.The same as fig.9a-f, but using a Dieterich-Ruina governing law.Medium and constitutive parameters are reported in table III.
Fig.10a-f.The same as fig.9a-f, but using a Dieterich-Ruina governing law.Medium and constitutive parameters are reported in table III.
Fig.9a-f.3D snapshots of the solutions of an enlarging crack with an initially oblique pre-stress vector (45°).a), b) and (c): slip velocity snapshots; d), e) and (f): total dynamic traction.The slip-weakening law has been adopted and the constitutive parameters are listed in table III.

Fig. 14 .
Fig. 14.Initial shear stress distribution (component along x1) on the fault plane for a numerical experiment with spatial heterogeneity.This map refer to the inferred distribution for the 1999 ChiChi, Taiwan earthquake.

Fig
Fig.15a-f.Slip (panels a, b and c) and slip velocity (panels d, e and f) time snapshots obtained for a numerical simulation in which the pre-stress is that plotted in fig.14and the fault is governed by a slip-weakening law.Parameters are listed in table V.

Fig. 16 .
Fig. 16.Spatial distribution of the rake variation (with respect to the initial value of 0°) resulting at the end of simulation corresponding to the results shown in fig.15a-f.The black contour lines superimposed to the rake variation represent the traction values at the same time step at which the rake variation is calculated.

Table I .
Medium and constitutive parameters for numerical experiments in figs.2a-c, 3a,b and 4a,b where a linear slip-weakening governing law has been imposed on the fault plane.All quantities are in non-dimensional units.The strength parameter S (Das and Aki, 1977a,b) is 0.8.

Table II .
Medium and constitutive parameters for numerical experiments in figs.5a-c, 6a,b and 7a,b where a Dieterich-Ruina Law has been adopted.All quantities are typical of a laboratory scale.

Table IV .
Medium and constitutive parameters for 3D simulations plotted in fig.13where a linear slipweakening law has been assumed.All quantities are in nondimensional units.The strength parameter S is 0.8 for all numerical experiments.

Table V .
Medium and constitutive parameters for simulation results reported in figs.15a-f and 16.The fault is governed by a linear slip-weakening law.All quantities refer to a real-world fault.