# Dynamic non-planar crack rupture by a finite volume method

Mondher Ben Jemaa, N. Glinsky-Olivier, V. M. Cruz-Atienza, J. Virieux, Serge Piperno

► To cite this version:

Mondher Ben Jemaa, N. Glinsky-Olivier, V. M. Cruz-Atienza, J. Virieux, Serge Piperno. Dynamic non-planar crack rupture by a finite volume method. Geophysical Journal International, 2007, 171 (1), pp.271-285. 10.1111/j.1365-246X.2006.03500.x. insu-00355237

HAL Id: insu-00355237

https://insu.hal.science/insu-00355237

Submitted on 20 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamic non-planar crack rupture by a finite volume method

M. Benjemaa,<sup>1,2</sup> N. Glinsky-Olivier,<sup>1</sup> V. M. Cruz-Atienza,<sup>2,3</sup> J. Virieux<sup>2</sup> and S. Piperno<sup>1</sup>

Accepted 2007 May 15. Received 2007 April 28; in original form 2006 January 9

## SUMMARY

Modelling dynamic rupture for complex geometrical fault structures is performed through a finite 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 fluxes and not on stress values. Prescribed rupture velocity gives accurate solutions with respect to analytical ones depending on the mesh refinement, 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 specific difficulties of such numerical modelling.

**Key words:** boundary conditions, dynamic rupture, finite volume approach, friction law, numerical methods, seismic source.

## 1 INTRODUCTION

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.

<footnote n="1"><sup>1</sup>INRIA-ENPC, CAIMAN project, Sophia-Antipolis, France. E-mail: Mondher.Ben.Jemaa@sophia.inria.fr</footnote>

<footnote n="2"><sup>2</sup>CNRS, Géosciences Azur, Sophia-Antipolis, France</footnote>

<footnote n="3"><sup>3</sup>Department of Geological Sciences, San Diego State University, USA</footnote>
## 2 DYNAMIC CRACK PROBLEM

We consider an isotropic, linearly elastic infinite medium containing a surface $\Gamma$ 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 $\Gamma$, the medium is governed by the following linearized elastodynamic equations,

$$ \rho \frac{\partial^2 \vec{u}}{\partial t^2} = \overline{\mathrm{div}} \underline{\sigma} \quad (1) $$

$$ \underline{\sigma} = \lambda \mathrm{div} \vec{u} I_2 + \mu [\vec{\nabla} \vec{u} + t'(\vec{\nabla} \vec{u})], \quad (2) $$

where the identity matrix is denoted by $I_2$, the displacement vector by $\vec{u}$, the symmetric stress tensor by $\underline{\sigma}$. The spatially varying density is denoted by $\rho$ and Lamé coefficients by $\lambda$ and $\mu$. The superscript $t$ denotes the transposition operation. We define the velocity vector $\vec{v}$ as the time derivative of the displacement vector $\vec{u}$. The following velocity-stress first-order hyperbolic system,

$$ \rho \frac{\partial \vec{v}}{\partial t} = \overline{\mathrm{div}} \underline{\sigma} \quad (3) $$

$$ \frac{\partial \underline{\sigma}}{\partial t} = \lambda \mathrm{div} \vec{v} I_2 + \mu [\vec{\nabla} \vec{v} + t'(\vec{\nabla} \vec{v})], \quad (4) $$

describes the propagation of elastic waves in the heterogeneous medium (Madariaga 1976; Aki & Richards 1980; Virieux 1986). An initial heterogeneous stress $\underline{\sigma}(\vec{x}, 0) = \underline{\sigma}_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 $\Gamma(\vec{x}, t)$ which may have a complex geometry and which may depend also on time, will be piecewise discretized and a normal vector $\vec{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 $\Gamma$, 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

$$ t' \underline{\sigma}(t, x) \vec{n} = g(t, \Psi) \quad \forall x \in \Gamma, \quad (5) $$

where $g$ is a function depending on time and a local set of a constitutive law parameters $\Psi$. The tangent to the surface is denoted by $\vec{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 $\Gamma$, we define limiting values of both the displacement and the velocity vector, respectively as :

$$ \vec{u}^{\pm}(\vec{x}, t) = \lim_{\epsilon \to 0} \vec{u}[\vec{x} \pm \epsilon \vec{n}(\vec{x}), t] \quad (6) $$

$$ \vec{v}^{\pm}(\vec{x}, t) = \lim_{\epsilon \to 0} \vec{v}[\vec{x} \pm \epsilon \vec{n}(\vec{x}), t]. \quad (7) $$

The slip $\mathcal{U}$ and the slip-rate $\mathcal{V}$ are, respectively, jumps of the tangential displacement and the tangential velocity vectors across the surface $\Gamma$. These quantities are numerically determined. With previous notations, we have the following expressions for the slip and the slip-rate magnitudes,

$$ \mathcal{U} = (\vec{u}^{+} - \vec{u}^{-}) \cdot \vec{t} \quad (8) $$

$$ \mathcal{V} = \frac{\partial \mathcal{U}}{\partial t} = (\vec{v}^{+} - \vec{v}^{-}) \cdot \vec{t}, \quad (9) $$

where the scalar product is denoted by a dot.

The crack surface $\Gamma(\vec{x}, t)$ expands during the rupture process from its initial configuration $\Gamma_0 = \Gamma(\vec{x}, 0)$ to the final one $\Gamma_f = \Gamma(\vec{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.

## 3 CONSERVATIVE FLUX FORMULATION

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.

### 3.1 Finite volume equations

Over the elastic domain we consider in this paragraph, the stress tensor $\underline{\sigma}$ can be split into a trace tensor $\underline{s} = 1/2 \mathrm{tr} \underline{\sigma} I_2$ and a deviatoric tensor $\underline{d} = \underline{\sigma} - \underline{s}$ (i.e. $\mathrm{tr} \underline{d} = 0$). For a 2-D geometry, we may consider equivalently the two numbers $T = (\sigma_{xx} + \sigma_{zz})/2$ and $T' = (\sigma_{xx} - \sigma_{zz})/2$ and the shear-stress $\sigma_{xz}$. We write the system (3)–(4) into the form

$$ \rho \frac{\partial v_x}{\partial t} = \frac{\partial (T + T')}{\partial x} + \frac{\partial \sigma_{xz}}{\partial z} \quad (10) $$

$$ \rho \frac{\partial v_z}{\partial t} = \frac{\partial \sigma_{xz}}{\partial x} + \frac{\partial (T - T')}{\partial z} \quad (11) $$

$$ \frac{1}{\lambda + \mu} \frac{\partial T}{\partial t} = \frac{\partial v_x}{\partial x} + \frac{\partial v_z}{\partial z} \quad (12) $$

$$ \frac{1}{\mu} \frac{\partial T'}{\partial t} = \frac{\partial v_x}{\partial x} - \frac{\partial v_z}{\partial z} \quad (13) $$

$$ \frac{1}{\mu} \frac{\partial \sigma_{xz}}{\partial t} = \frac{\partial v_z}{\partial x} + \frac{\partial v_x}{\partial z} \quad (14) $$

which can be written into a pseudo-conservative form

$$ \Lambda \vec{W}_t = \mathrm{div} \vec{\mathcal{F}}(\vec{W}), \quad (15) $$

where the subscript $t$ means the temporal derivative. The vector $\vec{W} = t'(v_x, v_z, T, T', \sigma_{xz})$ has to be estimated. The diagonal matrix
$\Lambda = \text{diag}(\rho, \rho, \frac{1}{\lambda+\mu}, \frac{1}{\mu}, \frac{1}{\mu})$ contains the material description, while the transformation $\mathcal{F}$ is a linear function of $\vec{W}$ given by

$$ \vec{\mathcal{F}} = (\mathcal{F}_1, \mathcal{F}_2) $$

where

$$ \mathcal{F}_1(\vec{W}) = {}^t(T + T', \sigma_{xz}, v_x, v_x, v_z) $$

and

$$ \mathcal{F}_2(\vec{W}) = {}^t(\sigma_{xz}, T - T', v_z, -v_z, v_x). $$

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 pre-specified 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 $\mathcal{T}_i$. Assuming both the solution and the medium's characteristics constant in each cell, we obtain by the Green theorem,

$$ \int_{\mathcal{T}_i} \Lambda_i \vec{W}_i = \int_{\partial \mathcal{T}_i} \vec{\mathcal{F}}(\vec{W}) \vec{n} \, dS, \qquad (16) $$

where $\partial \mathcal{T}_i$ designs the boundary of the cell $\mathcal{T}_i$, $\Lambda_i$ contains the values of elastic parameters inside the cell $\mathcal{T}_i$ and $\vec{n}$ is the unitary outwards normal vector. The eq. (16) is approximated by,

$$ A_i \Lambda_i (\vec{W}_i)_{\mathcal{T}_i} = \sum_{\mathcal{T}_j \in V(\mathcal{T}_i)} \Phi_{ij}, \qquad (17) $$

where $A_i$ is the surface of the cell $\mathcal{T}_i$. The expression $(\vec{W}_i)_{\mathcal{T}_i}$ is an approximation of $\vec{W}_i$ inside the cell $\mathcal{T}_i$. For each cell $\mathcal{T}_i$, the set of neighbouring cells is denoted by $V(\mathcal{T}_i)$. The numerical flux integral across the interface $\mathcal{T}_{ij} = \mathcal{T}_i \cap \mathcal{T}_j$ between $\mathcal{T}_i$ and $\mathcal{T}_j$ is denoted by $\Phi_{ij}$. The scheme is conservative as the following equality

$$ \Phi_{ij} + \Phi_{ji} = 0 \qquad (18) $$

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 $\Gamma(\vec{x}, t)$. We can write the following expression,

$$ \int_{\mathcal{T}_{ij}} \vec{\mathcal{F}}(\vec{W}) \vec{n} \, dS \approx \Phi_{ij} \equiv \vec{\mathcal{F}}\left(\frac{W_i + W_j}{2}\right) \int_{\mathcal{T}_{ij}} \vec{n} \, dS \\ = \frac{\vec{\mathcal{F}}(W_i) + \vec{\mathcal{F}}(W_j)}{2} \vec{n}_{ij}, $$

where arithmetic means of fields are used in this flux evaluation along the edge between two cells. The edge normal $\vec{n}_{ij} = \int_{\mathcal{T}_{ij}} \vec{n} \, dS$

<figure><img src="image_2.png" /><figcaption>Figure 1. Two contiguous cells within the mesh. $\rho_i$ and $\tilde{\Lambda}_i$ are the local elastic properties of the medium in the cell $\mathcal{T}_i$. The unknown $\vec{W}$ is assumed constant in each cell, and $\Phi_{ij}$ design the fluxes between $\mathcal{T}_i$ and $\mathcal{T}_j$.</figcaption></figure>

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 $\mathcal{T}_i$ to $\mathcal{T}_j$ when time integration is performed for the $\mathcal{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 $\vec{\gamma} = {}^t(T, T', \sigma_{xz})$ the stress part of the vector $\vec{W}$. Using this notation, the previous flux integral is split into

$$ \Phi_{ij} = \Phi(W_i, W_j) = {}^t[\Phi_{\vec{v}}(\gamma_i, \gamma_j), \Phi_{\vec{\gamma}}(v_i, v_j)] \qquad (19) $$

with

$$ \Phi_{\vec{v}} = {}^t(\Phi_{v_x}, \Phi_{v_z}) \quad \text{and} \quad \Phi_{\vec{\gamma}} = {}^t(\Phi_T, \Phi_{T'}, \Phi_{\sigma_{xz}}) $$

given by the following expressions,

$$ \Phi_{v_{xij}} = \frac{T_i + T_j + T'_i + T'_j}{2} n_{x_{ij}} + \frac{\sigma_{xz_i} + \sigma_{xz_j}}{2} n_{z_{ij}} \qquad (20) $$

$$ \Phi_{v_{zij}} = \frac{\sigma_{xz_i} + \sigma_{xz_j}}{2} n_{x_{ij}} + \frac{T_i + T_j - T'_i - T'_j}{2} n_{z_{ij}} \qquad (21) $$

$$ \Phi_{T_{ij}} = \frac{v_{x_i} + v_{x_j}}{2} n_{x_{ij}} + \frac{v_{z_i} + v_{z_j}}{2} n_{z_{ij}} \qquad (22) $$

$$ \Phi_{T'_{ij}} = \frac{v_{x_i} + v_{x_j}}{2} n_{x_{ij}} - \frac{v_{z_i} + v_{z_j}}{2} n_{z_{ij}} \qquad (23) $$

$$ \Phi_{\sigma_{z_{ij}}} = \frac{v_{z_i} + v_{z_j}}{2} n_{x_{ij}} + \frac{v_{x_i} + v_{x_j}}{2} n_{z_{ij}}. \qquad (24) $$

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

$$ \rho_i v_i^{n+\frac{1}{2}} = \rho_i v_i^{n-\frac{1}{2}} + \frac{\Delta t}{\mathcal{A}_i} \sum_{j \in V(i)} \Phi_{\vec{v}}(\gamma_i^n, \gamma_j^n) \qquad (25) $$

$$ \tilde{\Lambda}_i \gamma_i^{n+1} = \tilde{\Lambda}_i \gamma_i^n + \frac{\Delta t}{\mathcal{A}_i} \sum_{j \in V(i)} \Phi_{\vec{\gamma}} \left( v_i^{n+\frac{1}{2}}, v_j^{n+\frac{1}{2}} \right) \qquad (26) $$
with the discrete matrix $\tilde{\Lambda}_i = \text{diag}(\frac{1}{\lambda_i+\mu_i}, \frac{1}{\mu_i}, \frac{1}{\mu_i})$. Replacing the fluxes by their respective expressions, we obtain the following 2-D discrete system of five equations,

$$ \rho_i v_i^{n+\frac{1}{2}} = \rho_i v_i^{n-\frac{1}{2}} + \frac{\Delta t}{\mathcal{A}_i} \sum_{j \in V(i)} N_{ij} \frac{\gamma_i^n + \gamma_j^n}{2} \quad (27) $$

$$ \tilde{\Lambda}_i \gamma_i^{n+1} = \tilde{\Lambda}_i \gamma_i^n + \frac{\Delta t}{\mathcal{A}_i} \sum_{j \in V(i)} t N_{ij} \frac{v_i^{n+\frac{1}{2}} + v_j^{n+\frac{1}{2}}}{2}, \quad (28) $$

with the following discrete geometric matrix $N_{ij} = (n_{xij} \ n_{xij} \ n_{zij})$ depending on cell edges. This discrete time-evolution 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 $\Delta 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.

## 3.2 Energy consideration

Let us consider the energy $E$ of the system defined by

$$ E = \underbrace{\int_{\Omega} \frac{1}{2} \rho \| \vec{v} \|^2 + \int_{\Omega} \frac{1}{2} t \vec{\sigma} \cdot \vec{\epsilon}}_{E_c} - \underbrace{\int_{\Omega} \frac{1}{2} t \vec{\sigma} \cdot \vec{\epsilon}}_{E_m}, \quad (29) $$

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:

$$ E_m = \int_{\Omega} \frac{1}{2} t \vec{\sigma} \cdot C \vec{\sigma} = \int_{\Omega} \frac{1}{2} t \vec{\gamma} \tilde{\Lambda} \vec{\gamma}. \quad (30) $$

Thus, the total energy in the continuum,

$$ E = \int_{\Omega} \left( \frac{1}{2} \rho \| \vec{v} \|^2 + \frac{1}{2} t \vec{\gamma} \tilde{\Lambda} \vec{\gamma} \right), \quad (31) $$

is transformed into the discrete total energy,

$$ \mathcal{E}^n = \frac{1}{2} \sum_i \mathcal{A}_i \left[ \rho_i t (v_i^{n+\frac{1}{2}}) v_i^{n-\frac{1}{2}} + t (\gamma_i^n) \tilde{\Lambda}_i \gamma_i^n \right], \quad (32) $$

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.

## 4 NUMERICAL BOUNDARY CONDITIONS

Let us now consider the crack surface $\Gamma$ 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

$$ \Delta \mathcal{E}^n = \sum_{\substack{i,j,l \\ \partial \mathcal{T}_i \cap \partial \mathcal{T}_j \subset \Gamma}} \frac{\Delta t}{2} \left( \vec{F}_{ij}^{n+\frac{1}{2}} + \vec{F}_{ji}^{n+\frac{1}{2}} \right) \cdot \vec{n}_{ij}, \quad (33) $$

where the summation is over the crack edges (see Appendix B). Each energy variation $\Delta \mathcal{E}_{ij}^n$ is related to the cell $\mathcal{T}_i$ towards the cell $\mathcal{T}_j$. The local discrete energy rate vector $\vec{F}$ is given by

$$ \vec{F}_{ij}^{n+\frac{1}{2}} = \left\{ \begin{aligned} & [T_i^{(n+\frac{1}{2})} + T_i^{(n+\frac{1}{2})}] v_{x_j}^{n+\frac{1}{2}} + \sigma_{xz_i}^{(n+\frac{1}{2})} v_{z_j}^{n+\frac{1}{2}} \\ & - [\sigma_{xz_i}^{(n+\frac{1}{2})} v_{x_j}^{n+\frac{1}{2}} + [T_i^{(n+\frac{1}{2})} - T_i^{(n+\frac{1}{2})}] v_{z_j}^{n+\frac{1}{2}} \end{aligned} \right\} $$

while the normal vector $\vec{n}_{ij} = (n_{xij} \ n_{zij})$ is oriented as usual from the cell $\mathcal{T}_i$ towards the cell $\mathcal{T}_j$. The stress estimation at half-integer time steps is obtained through the following time averaging expression

$$ \gamma_i^{(n+\frac{1}{2})} = \frac{\gamma_i^n + \gamma_i^{n+1}}{2}. $$

If we assume continuity of velocity and stress fields through these crack edges, the discrete total energy time variation will be zero.

### 4.1 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 $\vec{n}_{ij} = (n_{xij} \ n_{zij})$, with $n_{xij} = \pm 1$ depending on which side of the crack line we are. The energy variation is then reduced to the following expression,

$$ \Delta \mathcal{E}^n = \frac{\Delta t}{2} \sum_{\substack{i,j \subset \Gamma \\ \partial \mathcal{T}_i}} \left\{ \sigma_{xz_i}^{(n+\frac{1}{2})} v_{x_j}^{n+\frac{1}{2}} + [T_i^{(n+\frac{1}{2})} - T_i^{(n+\frac{1}{2})}] v_{z_j}^{n+\frac{1}{2}} \right. \\ \left. + \sigma_{xz_j}^{(n+\frac{1}{2})} v_{x_i}^{n+\frac{1}{2}} + [T_j^{(n+\frac{1}{2})} - T_j^{(n+\frac{1}{2})}] v_{z_i}^{n+\frac{1}{2}} \right\} n_{z_{ij}}. \quad (34) $$

Without any stress boundary condition, this discrete energy should be kept constant. The variation $\Delta \mathcal{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

$$ \sum_{\substack{i,j \subset \Gamma \\ \partial \mathcal{T}_i}} \left[ \sigma_{xz_i} v_{x_j} + (T_i - T_i') v_{z_j} + \sigma_{xz_j} v_{x_i} + (T_j - T_j') v_{z_i} \right] n_{z_{ij}} = 0 \quad (35) $$

for every couple $(i,j)$ such that $\partial \mathcal{T}_i \cap \partial \mathcal{T}_j \subset \Gamma$.

Consider now two adjacent cells $\mathcal{T}_i$ and $\mathcal{T}_j$ sharing a common crack edge. Since, we allow discontinuities because of the local rupture, the fluxes integral $\Phi_{ij}$ and $\Phi_{ji}$ through the segment $\mathcal{T}_{ij} = \partial \mathcal{T}_i \cap \partial \mathcal{T}_j \subset \Gamma$ 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 $\Phi_{ij}^*$ instead of $\Phi_{ij}$ for the cell $\mathcal{T}_i$ and $\Phi_{ji}^*$ instead of $\Phi_{ji}$ for the cell $\mathcal{T}_j$. We should define rules how to estimate these new flux integrals $\Phi_{ij}^*$ and $\Phi_{ji}^*$ from local variables from cells $\mathcal{T}_i$ and $\mathcal{T}_j$.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Two contiguous cells above and below the fault Γ. [AB] design the same local horizontal crack edge between cells $\mathcal{T}_i$ and $\mathcal{T}_j$. Dashed lines specify fictitious cells.</figcaption></figure>

Due to the discontinuities we will introduce, the flux integrals $\Phi_{ij*}$ and $\Phi_{ji*}$ across the crack does not verify $\Phi_{ij*} + \Phi_{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

$$ \Delta\mathcal{E}_{ji*}^n + \Delta\mathcal{E}_{ij*}^n = 0. \quad (36) $$

In other words, the following equality

$$ \sigma_{xz_i} v_{x_j} + (T_i^* - T_{j^*}) v_{z_j} + \sigma_{xz_j} v_{x_{j^*}} + (T_j - T_{j^*}) v_{z_{j^*}} - [\sigma_{xz_i} v_{x_{j^*}} + (T_i - T_{j^*}) v_{z_{j^*}} + \sigma_{xz_{j^*}} v_{x_i} + (T_{j^*} - T_{j^*'}) v_{z_{j^*}}] = 0 \quad (37) $$

must be verified at each time $t$ and for all $i$ and $j$ such that $\partial \mathcal{T}_i \cap \partial \mathcal{T}_j \subset \Gamma$ considering fictitious cells $\mathcal{T}_{j^*}$ and $\mathcal{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 $\mathcal{T}_{j^*}$ cell by assigning specific values of both $\mathcal{T}_i$ and $\mathcal{T}_j$ cells through

$$ v_{x_{j^*}} = v_{x_i} \quad (38) $$

$$ v_{z_{j^*}} = v_{z_j} \quad (39) $$

$$ T_{j^*} = T_j \quad (40) $$

$$ T'_{j^*} = T'_{j} \quad (41) $$

$$ \sigma_{xz_{j^*}} = -\sigma_{xz_i}. \quad (42) $$

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,

$$ \sigma_{xz_{j^*}} = -\sigma_{xz_i} + 2g \quad (43) $$

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 $\Phi_{ij*}$ across $\mathcal{T}_{ij}$ by,

$$ \Phi_{v_{x_{j^*}}} = \frac{T_i + T_j + T'_i + T'_j}{2} n_{x_{j^*}} + g n_{z_{j^*}} \quad (44) $$

$$ \Phi_{v_{z_{j^*}}} = g n_{x_{j^*}} + \frac{T_i + T_j - T'_i - T'_j}{2} n_{z_{j^*}} \quad (45) $$

$$ \Phi_{T_{j^*}} = v_{x_i} n_{x_{j^*}} + \frac{v_{z_i} + v_{z_j}}{2} n_{z_{j^*}} \quad (46) $$

$$ \Phi_{T'_{j^*}} = v_{x_i} n_{x_{j^*}} - \frac{v_{z_i} + v_{z_j}}{2} n_{z_{j^*}} \quad (47) $$

$$ \Phi_{\sigma_{xz_{j^*}}} = \frac{v_{z_i} + v_{z_j}}{2} n_{x_{j^*}} + v_{x_i} n_{z_{j^*}}. \quad (48) $$

We may proceed similarly for the other flux integral $\Phi_{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 $\Phi_{ij*}$ and $\Phi_{ji*}$, as shown above. Unambiguous rules have been elaborated for boundary conditions across an horizontal crack segment.

## 4.2 Local arbitrary oriented crack edge

Let us consider a local crack edge making an angle $\theta$ 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 $\bar{\sigma}$ and $\bar{v}$, using the transformation matrix $P_{\theta}$ between coordinate systems. The following standard relationships could be deduced:

$$ \bar{\sigma} = P_{\theta}^{-1} \sigma P_{\theta} \quad \text{and} \quad \bar{v} = P_{\theta}^{-1} v. $$

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 $\bar{\sigma}_{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 $\tau$ 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.

# 5 SELF-SIMILAR CRACK WITH CONSTANT RUPTURE VELOCITY

For a self-similar planar crack with a bilateral propagation at a constant velocity, Kostrov (1964) has obtained an analytical solution
for a sub-Rayleigh rupture velocity. In the 2-D geometry, the slip velocity follows the analytical expression,

$$v_x(t, x, 0) = C \text{Re} \left( \frac{t/x}{\sqrt{t^2/x^2 - v_r^{-2}}} \right),$$

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 \alpha$ and we consider a Poissonian medium where the ratio between P and S velocities is $\sqrt{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 $\tau_f$, when the point ruptures with the prescribed rupture velocity. Stresses are normalized by the stress drop $\tau_0 - \tau_f$, where $\tau_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:

$$\rho \frac{\partial \vec{v}}{\partial t} = \overrightarrow{\text{div}} \sigma \qquad (49)$$

$$\frac{\partial \sigma}{\partial t} = \lambda \operatorname{div} \vec{v} I_2 + \mu [\vec{\nabla} \vec{v} + l'(\vec{\nabla} \vec{v})] + \eta \{\lambda \operatorname{div} \vec{v} I_2 + \mu [\vec{\nabla} \vec{v} + l'(\vec{\nabla} \vec{v})]\}, \quad (50)$$

where $\vec{v}$ is the time derivative of the velocity vector and $\eta$ 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 $\eta$ is determined numerically and turns out to be $0.5\Delta 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.

## 5.1 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.

## 6 DYNAMIC RUPTURE USING STRESS THRESHOLD CRITERION

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,b, 1985; Mikumo & 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;
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Comparison of the numerical (circles) and analytical (solid lines) solutions for the self-similar dynamic crack growth problem for (a) fault slip and (b) shear stress flux integral. Each panel correspond to a given crack inclination θ with respect the horizontal axis.</figcaption></figure>

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:

$$ \tau_c = \begin{cases} \tau_u - (\tau_u - \tau_f) \frac{U}{\delta_0} & 0 \le U \le \delta_0 \\ \tau_f & U \ge \delta_0. \end{cases} \qquad (51) $$

The shear traction fluxes on the fault are bounded above by τ_c. We then verify the following jump conditions on the rupture surface

$$ \tau \le \tau_c \qquad (52) $$

$$ (\tau - \tau_c) V = 0, \qquad (53) $$

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.
<table><caption>**Table 2.** Different meshes and associated numerical quantities.</caption><thead><tr><th>Mesh</th><th><i>h</i> (m)</th><th>Δ <i>t</i> (10<sup>-3</sup> s)</th><th>Vertex number</th><th>Triangles number</th></tr></thead><tbody><tr><td>M1</td><td>100</td><td>7</td><td>4347</td><td>8452</td></tr><tr><td>M2</td><td>50</td><td>3</td><td>6037</td><td>11820</td></tr><tr><td>M3</td><td>25</td><td>1.5</td><td>14425</td><td>28536</td></tr><tr><td>M4</td><td>20</td><td>1.2</td><td>24681</td><td>49012</td></tr><tr><td>M5</td><td>10</td><td>0.6</td><td>54309</td><td>108180</td></tr><tr><td>M6</td><td>5</td><td>0.4</td><td>96573</td><td>192692</td></tr></tbody></table>

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, 5 s 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.

## 6.1 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

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Comparison between analytical (solid line) and numerical (circle) solution for the self-similar constant velocity rupture of (a) the shear stress and (b) the slip rate, at four points located on the crack. Left-hand panels are for $η = 0$ (i.e. without damping coefficient), while right-hand panels are for $η = 0.5Δt$. Time Envelope Misfit (TEM) and Frequency Envelope Misfit (FEM) are two misfit criteria for quantitative comparison of seismograms (see Kristekova et al. 2006, for details).</figcaption></figure>

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
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison of the numerical (circles) and analytical (solid lines) shear stress solutions for the self-similar dynamic crack growth problem for six different meshes.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Slip-weakening friction law. The curve represents the total shear stress as a function of cumulative slip.</figcaption></figure>

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 $\tau_u = 1.4$ MPa; $\tau_f = -2$ MPa and $\delta_0 = 0.25$ m with an initial shear stress $\tau_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 $\tau_f$. Due to

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Numerical solutions at a fault point located at 6 km from the end of the nucleation zone.</figcaption></figure>

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.

Fig. 10 shows the superposition of the slip rate on a point located at the middle of the spontaneous fault (left-hand panel) and the fault length (right-hand panel) as function of the time for five different meshes. The solutions are clearly dependent on mesh refinement. One can see that the mesh refinement induces the convergence of the different computed solutions to the finest one. The rupture times
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Velocity seismograms computed at seven points located around the fault, for six different orientations (θ = 0°, 9°, 18°, 27°, 36° and 45°).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Superposition of the slip (a) and the slip rate (b) computed in the middle point of the spontaneous rupture for six different fault inclination angles (θ = 0°, 9°, 18°, 27°, 36° and 45°).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Superposition of the slip rate computed in the middle point of the spontaneous rupture (right-hand side) and the fault length (left-hand side) for five different meshes.</figcaption></figure>

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*<sub>c</sub> 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 τ<sub>u</sub>. 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*<sub>c</sub>, which should be greater than eight for making the solution rather independent of the mesh definition.

## 6.2 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
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Difference in time of rupture as a function of the mesh size (left-hand side) and the number of mesh segments inside the cohesive zone (right-hand side), relative to a reference solution.</figcaption></figure>

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 $\tau_u$, we tested two cases for which only $\tau_u$ is changed. For the first case, we choose $\tau_u = 1.4$ MPa: the rupture propagates in a subshear regime. For the second case $\tau_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: $\tau_f = -2$ MPa and $\delta_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 $\tau_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.

## 7 NON-PLANAR FAULT GEOMETRY

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.

### 7.1 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: $\tau_u = 1.5$ MPa; $\tau_f = -3.3$ MPa and $\delta_0 = 0.05$ m with an initial shear stress $\tau_0 = 0$ MPa. The nucleation zone is governed by the following parameters: $\tau_u = \tau_0 = 0$ MPa; $\tau_f = -10$ MPa and $\delta_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.

## 8 CONCLUSION

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 self-similar 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
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Comparison of the shear stress and the slip rate solutions obtained by our finite volume scheme and the finite difference method (Cruz-Atienza & Virieux 2004), in three points located on the crack surface. The rupture propagates at (a) subshear regime and (b) supershear regime.</figcaption></figure>

by performing a simulation for a non-planar fault geometry embed-
ded 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.

ACKNOWLEDGMENTS

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
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Rupture length function of the time for two cases. The left-hand panel shows that the crack propagates at a subshear regime, while the right-hand panel shows that the crack propagates at a supershear regime.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**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).</figcaption></figure>

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.

## REFERENCES

Aki, K. & Richards, P., 1980. *Quantitative Seismology: Theory and Methods*, W.H. Freeman & Co, San Francisco.

Ando, R., Tada, T. & Yamashita, T., 2004. Dynamic evolution of a fault system through interactions between fault segments, *J. geophys. Res.*, **109**, doi:10.1029/2003JB002665.

Andrews, D., 1976a. Rupture propagation with finite stress in antiplane strain, *J. geophys. Res.*, **81**, 3575–3582.

Andrews, D., 1976b. Rupture velocity of plain strain shear cracks, *J. geophys. Res.*, **81**, 5679–5687.

Andrews, D., 1985. Dynamic plane-strain shear rupture with a slip weakening friction law calculated by a boundary integral method, *Bull. seism. Soc. Am.*, **75**, 1–12.

Aochi, H. & Fukuyama, E., 2002. Three-dimensional nonplanar simulation of the 1992 landers earthquake, *J. geophys. Res.*, **107**, 4.1–4.12.

Bécache, E. & Duong, T.H., 1994. A space-time variational formulation for the boundary integral equation in a 2D elastic crack problem, *J. geophys. Res.*, **28**, 141–176.

Benjemaa, M., Piperno, S. & Glinsky-Olivier, N., 2006. *Étude de stabilité d'un schéma volumes finis pour les équations de l'élastodynamique en maillages non structurés*, INRIA Sophia Antipolis. RR-5817, France.

Capdeville, Y., Chaljub, E., Vilotte, J.-P. & Montagner, J.-P., 2003. Coupling the spectral element method with a modal solution for elastic wave propagation in realistic 3D global Earth models, *Geophys. J. Int.*, **152**, 34–68.

Chaljub, E., Capdeville, Y. & Vilotte, J.-P., 2003. Solving elastodynamics in a fluid-solid heterogeneous sphere: a parallel spectral element approximation on non-conforming grids, *J. Comput. Phys.*, **187**, 457–491.

Cohen, G. & Fauqueux, S., 2001. 2D elastic modelling with efficient mixed finite elements, in *Extended Abstracts*, Eur. Ass. Expl. Geophys.

Cruz-Atienza, V.M., 2006. Rupture dynamique des faille non-planaires en différences finies, *PhD thesis*, University of Nice — Sophia Antipolis, France.

Cruz-Atienza, V.M. & Virieux, J., 2004. Dynamic rupture simulation of non-planar faults with a finite-difference approach, *Geophys. J. Int.*, **158**, 939–954.

Cruz-Atienza, V.M., Virieux, J. & Aochi, H., 2007. 3d finite-difference dynamic-rupture modelling along non-planar faults, *Geophysics*, in press.

Das, S. & Aki, K., 1977. A numerical study of two-dimensional spontaneous rupture propagation, *Geophys. J. R. astr. Soc.*, **50**, 643–668.

Das, S. & Kostrov, B., 1987. On the numerical boundary integral equation method for three-dimension dynamic shear crack problems, *J. Appl. Mech.*, **54**, 99–104.

Day, S., 1977. Finite element analysis of seismic scattering problems, *PhD Dissertation*, **54**, 99–104.

Day, S., 1982. Three-dimensional simulation of spontaneous rupture: the effect of nonuniform prestress, *Bull. seism. Soc. Am.*, **72**, 1881–1902.

Day, S. & Ely, G.P., 2002. Effect of a shallow weak zone on fault rupture: Numerical simulation of scale-model experiments, *Bull. seism. Soc. Am.*, **92**, 3022–3041.

Day, S., Dalguer, L.A., Lapusta, N. & Liu, Y., 2005. Comparison of finite difference and boundary integral solutions to three-dimensional spontaneous rupture, *J. geophys. Res.*, **110**, B12307.

Dormy, E. & Tarantola, A., 1995. Numerical simulation of elastic wave propagation using a finite volume method, *J. geophys. Res.*, **100**, 2123–2133.

Festa, A. & Vilotte, J.-P., 2005. Spectral element simulations of dynamic rupture along kinked faults, *EGU Meeting*, pp. SRef-ID: 1607–7962/gra/EGU05-A-05122.

Fezoui, L., Lantéri, S., Lohrengel, S. & Piperno, S., 2005. Convergence and stability of a discontinuous galerkin time-domain method for the 3D heterogeneous maxwell equations on unstructured meshes, *ESAIM*, **39**, doi: 10.1051/m2an:2005049.

Harris, R., Archuleta, R. & Day, S., 1991. Fault steps and the dynamic rupture process: 2-D numerical simulations of a spontaneously propagating shear fracture, *Geophys. Res. Lett.*, **18**, 893–896.

Haskell, N., 1964. Total energy and energy spectral density of elastic wave radiation from propagating faults, *Bull. seism. Soc. Am.*, **54**, 1811–1842.

Huang, H. & Costanzo, F., 2004. On the use of space-time finite elements in the solution of elasto-dynamic fracture problems, *Int. J. Fract.*, **127**, 119–146.

Ida, Y., 1972. Cohesive force across the tip of a longitudinal-shear crack and griffith's specific surface energy, *J. geophys. Res.*, **77**, 3796–3805.

Kame, N. & Yamashita, T., 1999. Simulation of the spontaneous growth of a dynamic crack without constraints on the crack tip path, *Geophys. J. Int.*, **139**, 345–358.

Käser, M. & Iske, A., 2005. Ader schemes on adaptive triangular meshes for scalar conservation laws, *J. Comput. Phys.*, **205**, 486–508.

Knopoff, L. & Ni, X., 2001. Numerical instability at the edge of a dynamic fracture, *Geophys. J. Int.*, **147**, F1–F6.

Komatitsch, D. & Vilotte, J.-P., 1998. The spectral element method: an efficient tool to simulate the seismic response of 2-D and 3-D geological structures, *Bull. seism. Soc. Am.*, **88**, 368–392.
Kostrov, B., 1964. Selfsimilar problems of propagation of shear cracks, *PMM*, **28**, 889–898.

Kristekova, M., Kristek, J., Moczo, P. & Day, S., 2006. Misfit criteria for quantitative comparison of seismograms, *Bull. seism. Soc. Am.*, **96**(5), 1836–1850.

LeVeque, R., 2002. *Finite Volume Methods for Hyperbolic Problems*, Cambridge University Press, UK.

Madariaga, R., 1976. Dynamics of an expanding circular fault, *Bull. seism. Soc. Am.*, **66**, 639–666.

Madariaga, R., 2005. Seismic energy radiation from dynamic faulting, *EGU Meeting*, pp. SRef–ID: 1607–7962/gra/EGU05–A–05775.

Madariaga, R., Olsen, K. & Archuleta, R., 1998. Modeling dynamic rupture in a 3d earthquake fault model, *Bull. seism. Soc. Am.*, **88**, 1182–1197.

Mikumo, T. & Miyatake, T., 1978. Dynamical rupture process on a three-dimensional fault with non-uniform frictions and near-field seismic waves, *Geophys. J. R. astr. Soc.*, **54**, 417–438.

Moës, N. & Belytschko, T., 2002. Extended finite element method for cohesive crack growth, *Eng. Fract. Mech.*, **69**, 813–833.

Olsen, K.B., Madariaga, R. & Archuleta, R., 1997. Three-dimensional dynamic simulation of the 1992 landers earthquake, *Science*, **278**, 834–838.

Palmer, A. & Rice, J., 1973. The growth of slip surfaces in the progressive failure of overconsolidated clay slopes, *Proc. R. Soc. Lond.*, **A332**, 537.

Piperno, S., Remaki, M. & Fezoui, L., 2002. A non-diffusive finite volume scheme for the 3d maxwell equations on unstructured meshes, *SIAM J. Numer. Anal.*, **39**, 2089–2108.

Remaki, M., 2000. A new finite volume scheme for solving maxwell system, *COMPEL*, **19**(3), 913–932.

Tada, T. & Madariaga, R., 2001. Dynamic modelling of the flat 2-D crack by a semi-analytical BIEM scheme, *Int. J. Numer. Methods Eng.*, **50**, 227–251.

Tada, T. & Yamashita, T., 1977. Non-hypersingular boundary integral equations for two-dimensional non-planar crack analysis, *Geophys. J. Int.*, **130**, 269–282.

Virieux, J., 1986. P-SV wave propagation in heterogeneous media, velocity-stress method, *Geophysics*, **51**, 889–901.

Virieux, J. & Madariaga, R., 1982. Dynamic faulting studied by a finite difference method, *Bull. seism. Soc. Am.*, **72**, 345–369.

Zhang, J., 2005. Elastic wave modeling in fractured media with an explicit approach, *Geophysics*, **70**, 75–85.

## APPENDIX A:

We shall express the mechanical energy with respect to the stress components $\vec{\gamma}$. Because the stress tensor is symmetric, it is somehow easier to write it down in a vectorial form as

$$ \vec{\sigma} = \begin{pmatrix} \sigma_{xx} \\ \sigma_{zz} \\ \sigma_{xz} \end{pmatrix} \qquad \text{and it is straightforward to check the following equality} \\ \qquad \vec{\sigma} = \underline{M} \vec{\gamma}, \tag{A1} $$

where

$$ \underline{M} = \begin{pmatrix} 1 & 1 & 0 \\ 1 & -1 & 0 \\ 0 & 0 & 1 \end{pmatrix}. $$

We need to introduce as well the deformation tensor

$$ \underline{\epsilon} = \begin{pmatrix} \epsilon_{xx} & \epsilon_{xz} \\ \epsilon_{zx} & \epsilon_{zz} \end{pmatrix}, $$

where

$$ \epsilon_{xx} = \frac{\partial u_x}{\partial x}, \quad \epsilon_{xz} = \epsilon_{zx} = \frac{1}{2} \left( \frac{\partial u_x}{\partial z} + \frac{\partial u_z}{\partial x} \right) $$

and

$$ \epsilon_{zz} = \frac{\partial u_z}{\partial z}. $$

Because the $\underline{\epsilon}$ is also symmetric, we can write it down in vectorial form as

$$ \vec{\underline{\epsilon}} = \begin{pmatrix} \frac{\partial u_x}{\partial x} \\ \frac{\partial u_z}{\partial z} \\ \frac{\partial u_x}{\partial z} + \frac{\partial u_z}{\partial x} \end{pmatrix}. $$

For an elastic medium, the generalized Hooke's law links deformation and stress through the linear relationship

$$ \vec{\underline{\epsilon}} = \underline{C} \vec{\sigma} \tag{A2} $$

with

$$ \underline{C} = \frac{1}{\mu} \begin{pmatrix} \frac{\lambda + 2 \mu}{4(\lambda + \mu)} & -\frac{\lambda}{4(\lambda + \mu)} & 0 \\ -\frac{\lambda}{4(\lambda + \mu)} & \frac{\lambda + 2 \mu}{4(\lambda + \mu)} & 0 \\ 0 & 0 & 1 \end{pmatrix}. $$

Therefore, the mechanical energy defined as $E_m = \int_{\Omega} \frac{1}{2} \underline{\epsilon} \cdot \vec{\underline{\epsilon}}$ should be expressed with deformation or stress only. Because we consider a crack problem, we express the mechanical energy with stress components. By eqs (A1) and (A2), we find

$$ \begin{aligned} E_m &= \int_{\Omega} \frac{1}{2} \underline{\epsilon} \cdot \underline{C} \underline{\sigma} \\ &= \int_{\Omega} \frac{1}{2} \underline{\gamma} \cdot \underline{\Lambda} \underline{\gamma} \\ &= \int_{\Omega} \frac{1}{\lambda + \mu} T^2 + \frac{1}{\mu} T'^2 + \frac{1}{\mu} \sigma_{xz}^2, \end{aligned} $$

where the matrix $\underline{\Lambda}$ stands for $\underline{M} \underline{C} \underline{M}$ (as $\underline{M}$ is symmetric).

## APPENDIX B:

The discrete total energy inside the cell $T_i$ at time $n$, expressed as

$$ \mathcal{E}_i^n = \frac{1}{2} A_i (\rho_i {}^t (v_i^{n+\frac{1}{2}}) v_i^{n-\frac{1}{2}} + {}^t (\gamma_i^n) \tilde{\Lambda}_i \gamma_i^n), \tag{B1} $$

could vary in time when we apply boundary conditions. The variation of this total energy between times $n$ and $n + 1$ is

$$ \begin{aligned} 2 \Delta \mathcal{E}_i^n &= 2 (\mathcal{E}_i^{n+1} - \mathcal{E}_i^n) \\ &= A_i \rho_i {}^t (v_i^{n+\frac{1}{2}}) (v_i^{n+\frac{3}{2}} - v_i^{n-\frac{1}{2}}) \\ &\quad + A_i {}^t (\gamma_i^{n+1} + \gamma_i^n) \tilde{\Lambda}_i (\gamma_i^{n+1} - \gamma_i^n) \end{aligned} $$

thanks to the symmetry of the matrix $\tilde{\Lambda}_i$. From eqs (27) to (28), we may deduce

$$ \frac{2 \Delta \mathcal{E}_i^n}{\Delta t} = \sum_{j \in V(i)} {}^t(v_i^{n+\frac{1}{2}}) N_{ij} [\gamma_i^{(n+\frac{1}{2})} + \gamma_j^{(n+\frac{1}{2})}] \\ + {}^t[\gamma_i^{(n+\frac{1}{2})}] {}^t N_{ij} (v_i^{n+\frac{1}{2}} + v_j^{n+\frac{1}{2}}). \tag{B2} $$

Normal vectors inside a cell should verify the following consistent relationship,

$$ \int_{\partial T_i} 1.N \, dS = \int_{T_i} \vec{\nabla} 1 \, dxdz = 0 \quad \Rightarrow \quad \sum_{j \in V(i)} N_{ij} = 0. \tag{B3} $$
Thus the local discrete energy variation is simplified into the expression,

$$ \frac{2 \Delta \mathcal{E}_i^n}{\Delta t} = \sum_{j \in V(i)} \left[ \begin{aligned} & v_i^{n+\frac{1}{2}} N_{ij} \gamma_j^{[n+\frac{1}{2}]} + v_i^{[n+\frac{1}{2}]} N_{ij} v_j^{n+\frac{1}{2}} \end{aligned} \right] \quad (B4) $$

and allows us to estimate the total discrete energy variation over cells by the following expression,

$$ \begin{aligned} 2 \Delta \mathcal{E}^n &= 2 (\mathcal{E}^{n+1} - \mathcal{E}^n) \\ &= \Delta t \sum_i \sum_{j \in V(i)} \left[ \begin{aligned} & v_i^{n+\frac{1}{2}} N_{ij} \gamma_j^{(n+\frac{1}{2})} + v_i^{[n+\frac{1}{2}]} N_{ij} v_j^{n+\frac{1}{2}} \end{aligned} \right] \\ &= \Delta t \sum_{\text{interfaces}} \left[ \begin{aligned} & v_i^{n+\frac{1}{2}} N_{ij} \gamma_j^{(n+\frac{1}{2})} + v_i^{[n+\frac{1}{2}]} N_{ij} v_j^{n+\frac{1}{2}} \end{aligned} \right] \end{aligned} $$

$$ \begin{aligned} & + v_i^{n+\frac{1}{2}} N_{ji} \gamma_i^{(n+\frac{1}{2})} + v_i^{[n+\frac{1}{2}]} N_{ji} v_i^{n+\frac{1}{2}} \end{aligned} \\ & + \Delta t \sum_{\text{external interfaces}} \left[ \begin{aligned} & v_i^{n+\frac{1}{2}} N_{ij} \gamma_j^{(n+\frac{1}{2})} + v_i^{[n+\frac{1}{2}]} N_{ij} v_j^{n+\frac{1}{2}} \end{aligned} \right] \end{aligned} $$

Since, we have this conventional identity $N_{ji} = -N_{ij}$ and since the velocity vector is discontinuous across crack surface $\Gamma$, we may deduce the time variation of the total discrete energy as

$$ \mathcal{E}^{n+1} = \mathcal{E}^n + \frac{\Delta t}{2} \sum_{i,j \in \Gamma} \left[ \begin{aligned} & v_i^{n+\frac{1}{2}} N_{ij} \gamma_j^{(n+\frac{1}{2})} + v_i^{[n+\frac{1}{2}]} N_{ij} v_j^{n+\frac{1}{2}} \end{aligned} \right] \quad (B5) $$

relating to the elastic energy release into the medium along the crack edge.