Dynamic non-planar crack rupture by a finite volume method

SUMMARY Modelling dynamic rupture for complex geometrical fault structures is performed through a ﬁnite volume method. After transformations for building up the partial differential system following explicit conservative law, we design an unstructured bi-dimensional time-domain numerical formulation of the crack problem. As a result, arbitrary non-planar faults can be explicitly represented without extra computational cost. On these complex surfaces, boundary conditions are set on stress ﬂuxes and not on stress values. Prescribed rupture velocity gives accurate solutions with respect to analytical ones depending on the mesh reﬁnement, while solutions for spontaneous propagation are analysed through numerical means. An example of non-planar spontaneous fault growth in heterogeneous media demonstrates the good behaviour of the proposed algorithm as well as speciﬁc difﬁculties of such numerical modelling.


I N T RO D U C T I O N
As the number of unsaturated seismograms recorded nearby the earthquake rupture zone increases dramatically, understanding the physics of the rupture process requires more and more sophisticated tools where the geometry of the ruptured surface is taken into account as well as realistic friction laws on this surface. New formulations have recently been proposed for modelling the dynamic shear crack rupture when considering the complexity of earthquake mechanisms embedded in heterogeneous crustal structure (Kame & Yamashita 1999;Ando et al. 2004;Cruz-Atienza & Virieux 2004;Huang & Costanzo 2004). Boundary integral equation formulations (Das & Aki 1977;Tada & Yamashita 1977;Andrews 1985;Bécache & Duong 1994;Tada & Madariaga 2001) provide highly accurate field estimations nearby the crack at the expense of a rather simple medium description. Finite element formulations (Day 1977;Cohen & Fauqueux 2001), especially spectral formulations (Komatitsch & Vilotte 1998;Capdeville et al. 2003;Chaljub et al. 2003), have proved to be quite accurate for handling spontaneous propagation while considering complex crack structure (Festa & Vilotte 2005). Finally, finite difference approach (Madariaga 1976;Virieux & Madariaga 1982;Olsen et al. 1997;Madariaga et al. 1998) turns out to be quite efficient for 3-D configurations at the expense of less accurate field estimations nearby the crack (Madariaga 2005). Differences between these methods depend essentially on how boundary conditions are imposed on the discrete numerical formulation.
Other numerical methods, as the finite volume (FV) approach (LeVeque 2002), have been used for wave propagation, applied to the elastodynamic equations with mitigated results (Dormy & Tarantola 1995), although new interest has been raised recently by Käser & Iske (2005). Zhang (2005) proposed an hybrid scheme by writing the integral forms of the elastic-momentum equations together with a triangular finite difference operator with good results when considering scattering by cracks in both homogeneous and heterogeneous media. Other approaches using various formulations as finite element and discontinuous Galerkin methods for dynamic rupture in elastic media could be found in Moës & Belytschko (2002) and Huang & Costanzo (2004). They rely essentially on weak formulations of boundary conditions. With that respect, we believe that our presentation based on a discrete energy conservation is new.
We propose in this paper a complete reanalysis of the FV approach with a great attention to boundary conditions. We shall apply it to the dynamic crack rupture problem for an arbitrary geometry of the crack surface. We first introduce the elastodynamic equations as well as boundary conditions to be applied on the crack surface. We build up the FV scheme as a piecewise constant description of both velocity and stress on triangular mesh in a 2-D geometry which is equivalent to a P0 discontinuous Galerkin approach (Remaki 2000;Piperno et al. 2002). Through a careful analysis of total discrete energy, we specify boundary conditions on the crack in order to insure the correct discrete energy time variation and, therefore, the system stability. Accuracy of results will be checked against selected analytical solutions when a rupture velocity is specified. Finally, we discuss the spontaneous dynamic crack rupture by considering a simple slip-weakening law. Influence of meshing structures will be analysed before we end up by different illustrations of non-planar rupture evolution in a heterogeneous medium.

D Y N A M I C C R A C K P RO B L E M
We consider an isotropic, linearly elastic infinite medium containing a surface across which the displacement vector may have an unknown discontinuity while stress conditions are specified on this surface. This so-called crack problem is very different from the kinematic formulation where the displacement discontinuity is specified as introduced in seismology by Haskell (1964). Although we shall restrict our problem to 2-D geometry for illustration, the FV formulation we propose could be straightforwardly extended to 3-D geometry.
Away from the fracture surface , the medium is governed by the following linearized elastodynamic equations, where the identity matrix is denoted by I 2 , the displacement vector by u, the symmetric stress tensor by σ . The spatially varying density is denoted by ρ and Lamé coefficients by λ and µ. The superscript t denotes the transposition operation. We define the velocity vector v as the time derivative of the displacement vector u. The following velocity-stress first-order hyperbolic system, describes the propagation of elastic waves in the heterogeneous medium (Madariaga 1976;Aki & Richards 1980;Virieux 1986). An initial heterogeneous stress σ ( x, 0) = σ 0 could be defined inside the medium from previous loading histories (Virieux & Madariaga 1982). However, in this paper, we shall only consider uniform prestress conditions which can be set to zero. The crack surface ( x, t) which may have a complex geometry and which may depend also on time, will be piecewise discretized and a normal vector n is defined at each segment of the crack surface. We shall consider a frictional resistance on the crack surface. More specifically, we deal with an in-plane fracture mode, and we suppose that the contact between the two sides of the fracture is perfect. This means that no opening mechanism happens during the process, thanks to the confining pressure in the Earth crust. Inside the crack surface , the tangential stress, also called the shear stress, is assumed to drop down to its dynamic frictional level using a specific constitutive law we shall discuss later. The shear stress verifies the relationship where g is a function depending on time and a local set of a constitutive law parameters . The tangent to the surface is denoted by t. We assume that this function does not depend on the normal stress. We shall assume as well that the function g(t) is spatially invariant on the crack surface, although the numerical method we develop might handle more complex friction behaviours for other applications. Because, we allow a displacement discontinuity across the surface , we define limiting values of both the displacement and the velocity vector, respectively as : The slip U and the slip-rate V are, respectively, jumps of the tangential displacement and the tangential velocity vectors across the surface . These quantities are numerically determined. With previous notations, we have the following expressions for the slip and the slip-rate magnitudes, where the scalar product is denoted by a dot. The crack surface ( x, t) expands during the rupture process from its initial configuration 0 = ( x, 0) to the final one f = ( x, T f ) at the time T f when the rupture process stops, while waves are still propagating inside the medium. Whatever the dynamic fracture mechanism is, it depends critically on the accuracy of the elastic field estimation nearby the crack surface. It will be our main concern when considering the numerical implementation of boundary conditions.

C O N S E RVAT I V E F L U X F O R M U L AT I O N
In order to solve the elastodynamic equations by a FV method, we transform the system (3)-(4) into a conservative formulation on which we apply FV discretization. We identify the discrete total energy inside the elastic medium and its time variation related to the energy release when the crack rupture occurs. The study of the energy variation allows us to define the appropriate boundary conditions on the crack surface. This is the new feature we would like to stress in our work.

Finite volume equations
Over the elastic domain we consider in this paragraph, the stress tensor σ can be split into a trace tensor s = 1/2 trσ I 2 and a deviatoric tensor d = σ − s (i.e. tr d = 0). For a 2-D geometry, we may consider equivalently the two numbers T = (σ xx + σ zz )/2 and T = (σ xx − σ zz )/2 and the shear-stress σ xz . We write the system (3)-(4) into the form which can be written into a pseudo-conservative form = diag (ρ, ρ, 1 λ+µ , 1 µ , 1 µ ) contains the material description, while the transformation F is a linear function of W given by This function definition is only for compact notation of the eq. (15). Let us underline that the right-hand side (RHS) of this equation does not include medium properties description. In other words, specific parameters describing heterogeneities are grouped on the left-hand side (LHS) of the eq. (15) which allows a non-ambiguous use of centred space scheme as we will see in following paragraphs.
The elastic medium is divided in triangular cells in such a way that the crack surface coincides with edges of specific cells at any time. By anticipation, we may consider that the crack ruptures on a prespecified surface related to a mechanically weak zone of the Earth crust. Therefore, the initial meshing of the entire medium could be such that any evolution of the crack surface will match numerical edges of cells, an easy problem compared to new fractures on a fresh material.
The eq. (15) is integrated over each finite control surface (volume in 3-D geometry) or cell T i . Assuming both the solution and the medium's characteristics constant in each cell, we obtain by the Green theorem, where ∂T i designs the boundary of the cell T i , i contains the values of elastic parameters inside the cell T i and ñ is the unitary outwards normal vector. The eq. (16) is approximated by, where A i is the surface of the cell T i . The expression ( W t ) T i is an approximation of W t inside the cell T i . For each cell T i , the set of neighbouring cells is denoted by V (T i ). The numerical flux integral across the interface T i j = T i ∩ T j between T i and T j is denoted by ij . The scheme is conservative as the following equality is verified, thanks to the convention for normal vector orientation for edges of each cell. For cells having connection with the external edges of the domain, one must consider flux integral which may require specific attention as Absorbing conditions as studied by Benjemaa et al. (2006), which will be applied without discussion in this paper focused on the numerical crack boundary implementation. We use a centred scheme for a numerical approximation of the flux integral between contiguous cells for elements without edges on the crack surface ( x, t). We can write the following expression, where arithmetic means of fields are used in this flux evaluation along the edge between two cells. The edge normal is now isolated from field quantities W i and W j which may vary spatially and in time. The index ij specifies the direction from T i to T j when time integration is performed for the T i cell (Fig. 1). These centred numerical fluxes fulfil the conservative property we want to verify when the medium is continuous.
Let us denote γ = t (T, T , σ xz ) the stress part of the vector W . Using this notation, the previous flux integral is split into with v = t ( vx , vz ) and γ = t ( T , T , σxz ) given by the following expressions, For temporal integration, we use a leap-frog scheme where velocity is discretized at half-integer time steps and stress at integer time steps, which gives us the following discrete scheme Replacing the fluxes by their respective expressions, we obtain the following 2-D discrete system of five equations, with the following discrete geometric matrix N i j = ( n x i j n x i j n z i j n z i j −n z i j n x i j ) depending on cell edges. This discrete timeevolution system has a rather low number of arithmetic operations. A CFL criterion, for which the stability of this explicit time scheme can be proved in the general case of an unstructured mesh, will depend on the smallest triangle of the mesh. In other words, the discrete time increment t must be bounded by a value depending on the highest wave velocity and the smallest space path, which is taken for our unstructured mesh as the shortest height among all heights of triangles inside the mesh. A quite attractive feature would be a local time step. For more details on how to estimate this CFL value, we refer to the work of Benjemaa et al. (2006) as notations are quite tedious when considering unstructured meshes.

Energy consideration
Let us consider the energy E of the system defined by where E c is the kinetic energy of the system and E m its mechanical energy (see Appendix A for more details). By considering Hooke's law, we may express the mechanical energy with respect to stresses which are quantities specified at crack boundaries: Thus, the total energy in the continuum, is transformed into the discrete total energy, expressed at integer time steps n where the stress is estimated. Once the rupture process has stopped, we have verified that this discrete total energy inside the medium is kept constant, which is a consequence of the joint use of the leap-frog time-scheme and centred numerical fluxes (Fezoui et al. 2005), for the time and space discretizations of the conservative flux formulation (15). How numerically the total energy varies in time during the rupture process is described now.

N U M E R I C A L B O U N D A RY C O N D I T I O N S
Let us now consider the crack surface which coincides with edges of specific cells at any time. Let us remind that solutions could be discontinuous through these edges. Therefore, for these cells, a specific flux estimation should be performed on those edges belonging to the crack surface, called from now on crack edges. Using the scheme (27)-(28), the discrete total energy time variation is given by where the summation is over the crack edges (see Appendix B). Each energy variation E n i j is related to the cell T i towards the cell T j . The local discrete energy rate vector F is given by is oriented as usual from the cell T i towards the cell T j . The stress estimation at half-integer time steps is obtained through the following time averaging expression If we assume continuity of velocity and stress fields through these crack edges, the discrete total energy time variation will be zero.

Local horizontal crack
Without loss of generality, we may consider a horizontal segment of the crack surface with respect to the Cartesian coordinate system. Any other segment orientation could be considered by performing a local coordinate rotation as we shall see. Therefore, the normal vector of this horizontal segment is n i j = ( 0 nz i j ), with n z i j = ±1 depending on which side of the crack line we are. The energy variation is then reduced to the following expression, Without any stress boundary condition, this discrete energy should be kept constant. The variation E n must vanish for every time step n. For the sake of simplicity, we omit the temporal superscript and write down the energy conservation as Consider now two adjacent cells T i and T j sharing a common crack edge. Since, we allow discontinuities because of the local rupture, the fluxes integral ij and ji through the segment T i j = ∂T i ∩ ∂T j ⊂ could not be estimated through relations (17)-(21). The centred space scheme has to be modified for handling such field discontinuities and we must check that the eq. (35) will be verified when specific homogeneous boundary conditions are applied. For these reasons, we consider new flux integrals, say i j * instead of ij for the cell T i and ji * instead of ji for the cell T j . We should define rules how to estimate these new flux integrals i j * and ji * from local variables from cells T i and T j . Due to the discontinuities we will introduce, the flux integrals i j * and ji * across the crack does not verify i j * + ji * = 0 as for a continuous medium. Nevertheless, the total energy variation must remain equal to zero when specific homogeneous boundary conditions are applied. A simple way consists in verifying the local condition In other words, the following equality must be verified at each time t and for all i and j such that ∂T i ∩ ∂T j ⊂ considering fictitious cells T j * and T i * (see Fig. 2). The minus sign comes from the orientation of the normal crack vector.
Since we are interested in the in-plane fracture mode with no opening mechanism, the normal velocity component must be continuous, while the tangential velocity component is discontinuous. This leads us to the definition of local variables of the fictitious T j * cell by assigning specific values of both T i and T j cells through Similar equalities must be verified for the fictitious cell i * by replacing indexes j * by i * and i by j. Let us also remark that the eqs (40)-(42) are nothing but the continuity of the traction vector through the crack.
This definition enables us to conserve the discrete total energy as we sum up over the crack surface when no boundary conditions are specified through a local conservation over each segment. In order to satisfy the boundary condition (5), we modify the last eq. (42) as, prescribing the shear stress average value across the crack to g. The factor 2 comes from the centred scheme. The same equality must be taken for the cell i * . Following standard centred scheme, we are now able to estimate the flux i j * across T i j by, We may proceed similarly for the other flux integral ji * by inverting indexes i and j. Let us remark that fictitious cells are not specified in the computer code keeping memory management simple. Only related rules are applied for deducing the appropriate variable values when defining the split fluxes i j * and ji * , as shown above. Unambiguous rules have been elaborated for boundary conditions across an horizontal crack segment.

Local arbitrary oriented crack edge
Let us consider a local crack edge making an angle θ with respect to the Cartesian coordinate system (x, z). We may define a local Cartesian coordinate system (x , z ) with the axis x along the crack direction and express both the stress tensor and the velocity vector in the local Cartesian system, respectively, denoted as σ and v, using the transformation matrix P θ between coordinate systems. The following standard relationships could be deduced: We may apply the boundary condition (37) on this local crack edge which is now horizontal in this new Cartesian coordinate system. On this crack edge, the tangential shear stress σ xz is assumed to drop down to its dynamic friction level. More precisely, its flux integral through the local crack edge drops down to the dynamic level, while the cell value is already representative of the elastic medium response. We shall note this crack edge shear stress τ for clear distinction with the shear stress value itself inside the cell. Therefore, any crack shape could be considered as well by performing this local transformation for each individual edge of the complex crack surface: the crack surface is discretized as a subsequent edges at any time, making necessary the knowledge of the crack geometry before rupture initiation. More sophisticated strategies could be developed with adaptative remeshing as the crack surface expands. This important numerical investigation is left to further works.

S E L F -S I M I L A R C R A C K W I T H C O N S TA N T RU P T U R E V E L O C I T Y
For a self-similar planar crack with a bilateral propagation at a constant velocity, Kostrov (1964)   for a sub-Rayleigh rupture velocity. In the 2-D geometry, the slip velocity follows the analytical expression, where the factor C, called Kostrov constant, depends on the rupture velocity v r . This analytical solution may be used for the validation of the numerical solutions. We have selected a rupture velocity of 0.5 α and we consider a Poissonian medium where the ratio between P and S velocities is √ 3. Medium properties are taken from the first line of Table 1. At a given point of the crack plane, the flux integral value of the shear stress drops abruptly from the pre-stress level to the dynamic frictional stress value, say τ f , when the point ruptures with the prescribed rupture velocity. Stresses are normalized by the stress drop τ 0 − τ f , where τ 0 is the initial state of stress which gives a dimensionless stress drop equal to the unity.
Figs 3(a) and (b) show, respectively, the comparison between numerical and analytical solutions of the slip and the shear stress evaluation in four equidistant points placed along the crack plane for six different inclination angles with respect to the horizontal Cartesian axis. The numerical solution for the slip follows the analytic solution very closely whatever is the fault inclination. The slip is exactly equal to zero before the arrival of the rupture front and then increases hyperbolically as predicted by the theoretical study. The shear stress is also well modelled, especially the relaxation induced after the S wave arrives at the recorded point. For these simulations, we have selected numerical parameters of Table 2 using the discretization M4, where h denotes, hereafter, the mesh size along the crack.
Short period oscillations in the shear stress are observed due to the discrete stepwise progress of the fracture front. Dissipation terms can be introduced to control these spurious numerical oscillations (Virieux & Madariaga 1982;Knopoff & Ni 2001). We rewrite the system (3)-(4) as follows: where˙ v is the time derivative of the velocity vector and η is a damping coefficient to be determined. Due to our method which assumes that unknowns are piecewise-constant in each cell, we cannot add in a simple way a spatial second order derivative to the eq. (49) as usually used. We propose to add another term in the RHS of the eq. (50) which is exactly equivalent to the addition of a Laplacian term, but have the advantage to avoid spatial derivative computations. Fig. 4 shows the comparison between analytical and numerical solutions of the shear stress for the self-similar constant velocity rupture obtained without and with the damping term. The coefficient η is determined numerically and turns out to be 0.5 t for the M4 discretization of Table 2. Applying this artificial damping will be case-dependent because we have to avoid distortions in the build-up of the physical singularity of the shear-stress.

Mesh influence
One very interesting feature of the FV formulation is the capability of using simple unstructured triangular meshes. This allows us to describe quite accurately the geometry of the fault surface, especially when the geometry is non-planar. Realistic source geometries will modify quite significantly rupture behaviour as well as slip history over the fault surface (Aochi & Fukuyama 2002;Ando et al. 2004). Another advantage of considering unstructured meshes lies in the fact that one can refine the mesh nearby the fault zone in order to increase accuracy in field estimations at the expense of a fine time step which is until now selected globally for the entire medium. The mesh size along the crack h must be small enough to evaluate both shear stress concentration before the rupture as well as the finite peak associated with the S-wave motion. Fig. 5 shows the comparison between analytical and numerical solutions computed for different meshes, using medium parameters of Table 1. The different meshes were generated automatically by setting the segment length at external boundaries of the grid as well as on the fault surface. These meshes have the same mesh size at external boundaries and only the mesh size along the crack h varies, our objective being to find out the dependence of the numerical solutions on the mesh refinement around the rupture surface when neglecting the damping term. Various informations about the different meshes are given in Table 2. One can easily notice that the peak due to the S-wave travelling ahead the singularity may disappear if the mesh size along the fault is not enough refined. Moreover, one may notice that singular values depend critically on the mesh definition. We may see that it will not affect spontaneous rupture solutions which are now investigated.

D Y N A M I C RU P T U R E U S I N G S T R E S S T H R E S H O L D C R I T E R I O N
An important issue in seismology is the study of the stress conditions on faults before and during earthquakes, and the inference of a constitutive law that characterizes the material response to the applied stress. The friction constitutive relationship represents the governing equation of the failure process, and relates the stress field with fault slip and slip-rate among other physical parameters. The constitutive relationship is a key element of the dynamic descriptions of the seismic source which is based on models that satisfy the elastodynamic equations (Andrews 1976a(Andrews ,b, 1985Mikumo & Miyatake 1978;Day 1982;Das & Kostrov 1987;Harris et al. 1991). In the framework of fracture mechanics, an earthquake may be considered as a dynamically propagating shear crack that radiates seismic waves. The resulting motion on the fault is strongly related to the shear stress drop. Hence, the slip evolution depends on the failure criterion, the constitutive properties and the initial stress conditions on the fault surface, apart from fault geometry and medium properties. In contrast with the physically consistent dynamic models, kinematic models are widely accepted as a good description of the seismic source (Haskell 1964) prescribing the displacement history of motion a priori, without an explicit attempt to investigate the physical causes of the rupture process.
In our model, the constitutive relationship on the fault surface is assumed to be a slip-weakening (SW) friction law (Ida 1972  Palmer & Rice 1973), which is completely characterized by the yield stress τ u , the dynamic frictional stress τ f and the slip-weakening distance δ 0 (Fig. 6). Thus the frictional strength τ c is given by such a constitutive law, which may be written as follows: The shear traction fluxes on the fault are bounded above by τ c . We then verify the following jump conditions on the rupture surface which also prevent retrograde fault motion and allow rupture healing (see Day et al. 2005, for details about these jump conditions). As a result, a positive fault dislocation always takes place whenever the shear traction τ exceeds the fault strength τ c . Otherwise, the fault remains locked. Following eq. (51), rupture begins in a given point if τ exceeds the yield stress τ u . The fault strength then drops down to the dynamic level τ f as the slip grows over the critical distance δ 0 .  Fig. 7 displays several phase diagrams for one point located at 6 km from the end of the nucleation zone. For this simulation, we have used the following parameters: τ u = 1.3 MPa; τ f = −3.3 MPa and δ 0 = 0.4 m with an initial shear stress τ 0 = 0 MPa. For the initiation of the unilateral rupture, we impose the rupture in a 2 km long nucleation zone at one extremity of a 12 km fault. In this nucleation zone, the shear stress drops abruptly to the final level τ f . Again, the medium properties are given by the first line of Table 1.
After initiation, the rupture front propagates at subshear velocity. As its length increases, it becomes supershear and finally approaches the P-wave velocity. Due to the choice of the constitutive friction parameters, the rupture front exhibits the so-called bifurcation (Andrews 1976b): the rupture front jumps from a subshear to a supershear regime. Fig. 7(d) shows clearly that the observational point lies in a region where the rupture front has reached the supershear regime. One can see that the slip-rate peak (around 3 s) arrives before the slip-rate perturbation due to the S wave (around 3.5 s) travelling behind the rupture front. Figs 7(a) and (b) show the evolution of the shear stress as a function of the slip and slip-rate, respectively. One can note that the linear constitutive law is respected. Finally, Fig. 7(c) shows the slip history according to time. We note that there is no slip before the arrival of rupture front. Slope variations in the slip function around 3.5, 5s and 5.8 s correspond to the direct S wave and two back propagating P-and S-waves arrest pulses.
In order to check that our method does not depend on the fault orientation with respect to the Cartesian reference axis when rupture propagates spontaneously, we compare seismograms computed around the same spontaneous rupture case with different source orientations. Fig. 8 shows the superposition of the velocity components in seven points located around the fault, for six different fault orientations. We see a good agreement between all signals. For comparison, velocity components are expressed in the local reference frame (x , z ). The constitutive values used for this test case are : τ u = 1.7 MPa; τ f = −2 MPa and δ 0 = 0.2 m with an initial shear stress τ 0 = 0 MPa. The nucleation zone is 1.5 km long and lies on the left extremity of a 6 km fault. It is governed by the following parameters : τ u = τ 0 = 0 MPa; τ f = −10 MPa and δ 0 = 0.02 m. Fig. 9 shows the superposition of the slip and the slip-rate in the middle point of the spontaneous rupture region, for various fault orientations. Good estimates of the latter are also seen independently of the orientation of the fault. Furthermore we clearly see the P-stopping phase that abruptly changes the slip around 2 s after initiation (Fig. 9). Fig. 9(b) shows that the rupture front travels at a supershear regime. This is due to the choice of the constitutive friction parameters (Das & Aki 1977). Only small numerical oscillations are found, suggesting that quite stable solutions have been constructed.

Convergence of the spontaneous crack solution
We study again the influence of the mesh on numerical solutions. An essential requirement (even though not sufficient) for an accurate numerical method is that numerical solutions become independent of grid size. Since no theoretical solution is available for the spontaneous crack problem, we look for mesh refinements that yield rupture history independent of numerical discretization. During the steady crack propagation, the only physical length in our problem is the size of the fault region lying just behind the rupture front where the shear stress has not reached the dynamic friction level. This zone, known as the cohesive zone, is the place in which the breakdown process happens. For this reason, its correct sampling is a fundamental requirement for accurate rupture estimates. In order to yield numerical solutions independent of grid discretization, the number of mesh segments inside the cohesive zone, N c , must be kept large enough. This quantity N c represents the main numerical parameter controlling the convergence of crack solutions. Of course, the smaller the grid size along the fault h is, the bigger N c is. So we expect to find some minimal value for N c that assures numerical convergence.
For the spontaneous crack propagation with a slip weakening friction law, the cohesive zone is variable during the rupture and there is no a priori estimation of the numerical mesh density for adequate description of this weakening law. Hence, the quantity N c we consider in the following represents the average of all quantities N c along the spontaneous fault before the crack stops. The constitutive parameters for this simulation are τ u = 1.4 MPa; τ f = −2 MPa and δ 0 = 0.25 m with an initial shear stress τ 0 = 0 MPa. For the initiation of the unilateral rupture, we impose the rupture in a 2 km long nucleation zone at one extremity of a 6 km fault. In this nucleation zone, the shear stress drops abruptly to the final level τ f . Due to this choice of the constitutive parameters, the fault propagates at a subshear regime. This case is more suitable for the determination of the cohesive zone than the supershear regime for which this process zone is extremely variable.  converge toward the same value when the mesh size become smaller or equal than 50 m. Fig. 11 shows the root mean square (rms) of the rupture time difference (in percentage) as a function of the mesh size along the crack h for the left-hand panel, and as a function of the number of mesh segments inside the cohesive zone N c for the right-hand panel. The plotted circles represent the rms difference of rupture times relative to a finest mesh. The rupture time of a point on the fault plane is defined as the time at which the shear stress exceeds the yield stress τ u . We used a numerical solution computed with mesh size h = 10 m as reference solution. The rms differences follow a power law with estimated exponent between 1.8 and 2.1. These results seem in agreement with those found in Day et al. (2005). Although this study is not yet quantitative since no independent reference solution was considered, these numerical comparisons allow us to estimate the minimal number of mesh segments inside the cohesive zone, N c , which should be greater than eight for making the solution rather independent of the mesh definition.

Comparison with a finite difference method
As we pointed out in the previous section, the convergence of solutions when the mesh is refined is not sufficient to guarantee the accuracy of the numerical result when no theoretical solution exists (Day & Ely 2002). The comparison of different numerical methods can be helpful for this kind of problems. In what follows, we proceed to a comparison of the FV method we propose with a finite difference method introduced and validated by Cruz-Atienza (2006) and Cruz-Atienza et al. (2006) using a complete different numerical implementation of boundary conditions based on a strong treatment inside elements neighbouring the crack surface. This comparison validates the boundary conditions of our FV approach since both methods give quite similar solutions. Let us consider a dynamic crack problem with the slip weakening law (51). We select again medium properties using the first line of Table 1 and numerical parameters using the line M5 of Table 2. Since the crack propagation velocity depends on the material strength τ u , we tested two cases for which only τ u is changed. For the first case, we choose τ u = 1.4 MPa: the rupture propagates in a subshear regime. For the second case τ u = 0.5 MPa: the rupture propagates in a supershear regime. The other constitutive parameters are the same for both cases and are the following : τ f = −2 MPa and δ 0 = 0.25 m. For rupture initiation, we impose a 2 km nucleation zone in which the shear stress drops abruptly to a final level τ f . Figs 12(a) and (b) show, respectively, a comparison between the numerical solutions of the shear stress and the slip rate for the two test cases. The three observation points are located at L/2, 2L/3 and 3L/4, where L = 4000 m is the final spontaneous fault length (out of the nucleation zone). Results were obtained with h = 10 m for the two methods. Solutions are fairly similar for both rupture cases. Only small differences on the peak of the slip rate can be noted mainly in the subshear case. In the supershear case, both the time history and the amplitude of the seismograms are almost identical.
The rupture velocity is a critical parameter that strongly depends on the local properties of the solution. The crack tip evolves similarly for both numerical methods. Fig. 13 provides a quite impressive agreement between the FV method and the FD method for the dynamic crack tip position as a function of the time. We are confident not only in the convergence of our numerical scheme but also in the precision of the numerical solution we have obtained.

N O N -P L A N A R FAU LT G E O M E T RY
Triangular unstructured meshes used along this study allow us to consider both planar and non-planar fault geometries. Moreover, all variables of the system are computed in the same control volume while the rupture geometry is defined by pre-selected edges which may break or not. Let us underline that this discretization of the fault will depend on the mesh we use but the fault line will be sampled by edges and not by staircases which allows far more flexibility than the one proposed by FD methods. The FV scheme is quite adapted for both heterogeneous media and complex structures of faulting.
As an illustration of this flexibility, a complex test case dealing with a spontaneous rupture crossing a heterogeneity is now presented.

Complex fault geometry in heterogeneous media
Let us consider a non-planar fault propagating in a heterogeneous medium. The rupture crosses a low velocity zone LVZ during its evolution. Table 1 shows the different elastic properties of both media. The fault is governed by the SW friction law (51). Fig. 14 shows a snapshot of the horizontal velocity v x at 4 s after initiation. The rupture is 14.3 km long with a 1 km long nucleation zone located at the left edge of the fault, while the LVZ (circular dashed line) has a diameter of 4 km. The various constitutive parameters used for this simulation are : τ u = 1.5 MPa; τ f = −3.3 MPa and δ 0 = 0.05 m with an initial shear stress τ 0 = 0 MPa. The nucleation zone is governed by the following parameters : τ u = τ 0 = 0 MPa; τ f = −10 MPa and δ 0 = 0.02 m.
The unilateral rupture propagates rightwards at supershear velocity. The crack tip velocity reaches the S-wave velocity and approaches the P wave one. Results are similar to those presented by Cruz-Atienza & Virieux (2004), showing that the LVZ has important consequences. We have found that the crack tip velocity abruptly decreases inside such a zone, due to the direct relationship between the SW friction law and the elastic properties of the medium. We have also found an important increase of the slip and slip rate functions inside this zone. Furthermore, we can clearly identify in Fig. 14 the reflected P and S waves on the interface between the two media, especially inside the LVZ where back-propagating trapped waves are generated.

C O N C L U S I O N
A new flexible FV method to simulate the spontaneous growth of an in-plane shear crack has been presented. Thanks to an appropriate change of variables, all parameters of the medium are grouped on the left hand side of the elastodynamic equations and integration can be made even if the medium contains heterogeneities. The study of a suitable discrete expression of energy allows us to define the appropriate fracture boundary conditions to be imposed on the crack surface. The shear stress flux integral is set instead of the shear stress itself when applying boundary conditions. This makes the fracture to have no thickness, so both crack blocs only interact through the fracture traction vector. Consequently, different elastic properties at both sides of the fracture can be properly considered (e.g. bi-materials cracks). Numerically, this corresponds to the weak treatment of boundary conditions (set on fluxes), instead of a strong treatment of boundary conditions (set on elastic fields values). Spurious high frequency content in elastic fields at the vicinity of the crack tip is reduced as we reduce the mesh size near the fault surface, and the solution has a smoother behaviour as appropriate dissipation terms are added. Unstructured triangular meshes allow us to do this without important supplementary memory requirement. The comparison between numerical and analytical solutions for the selfsimilar constant velocity case, as well as comparisons made with an independent finite difference method (Cruz-Atienza & Virieux 2004) for different spontaneous rupture cases, revealed a very good agreement between the solutions, validating thus our approach for any kind of rupture geometry. The study of the influence of the mesh refinement on the numerical solutions shows that solutions are accurate enough if at least eight fault segments are found to be inside the cohesive zone. Finally, we illustrate the robustness of our method by performing a simulation for a non-planar fault geometry embedded in a heterogeneous medium. Results are in agreement with our expectations about the presence of low velocity zones during the dynamic rupture propagation. The FV method we propose seems to be a good alternative to the widely used finite difference and finite element methods, due to its geometrical flexibility and low computational cost. It is proved that this method is a second order space accuracy over a structured mesh (see Remaki 2000, for instance). The use of discontinuous Galerkin methods (DG), which can be thought as a FV methods of higher order, will improve the solution accuracy and should be investigated in future works.
Our study has been restricted to 2-D space domain, although the extension to 3-D space domain of the numerical rupture model does not require supplementary theoretical considerations. Only the computer task is more intensive.

A C K N O W L E D G M E N T S
We thank S. Lanteri for useful discussions and helpful comments during the work.
The authors are grateful to J.P. Ampuero for his constructive remarks that helped to improve the content of the paper and his careful checking of the accuracy of the solution. Raul Madariaga The left-hand panel shows that the crack propagates at a subshear regime, while the righthand panel shows that the crack propagates at a supershear regime. Figure 14. Snapshot of the horizontal particle velocity v x four second after rupture initiation. The non-planar fault grows through a circular low velocity zone (LVZ). Spontaneous rupture propagating rightwards is governed by the slip weakening friction law (eq. 51).
is thanked as an associate editor for positive remarks. We acknowledge partial support of scientific program QSHA PROJET ANR-05-CATT-011. This paper is a contribution of the UMR Géosciences Azur 6526. may deduce Normal vectors inside a cell should verify the following consistent relationship, and allows us to estimate the total discrete energy variation over cells by the following expression, Since, we have this conventional identity N ji = − N ij and since the velocity vector is discontinuous across crack surface , we may deduce the time variation of the total discrete energy as