# Seismic source dynamics, heterogeneity and friction

R. Madariaga, A. Cochard

► To cite this version:

R. Madariaga, A. Cochard. Seismic source dynamics, heterogeneity and friction. Annals of Geophysics, 1994, 37 (6), 10.4401/ag-4138. insu-03448581

HAL Id: insu-03448581

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

Submitted on 25 Nov 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
# Seismic source dynamics, heterogeneity and friction

Raúl Madariaga and Alain Cochard

Département de Sismologie, U.R.A. au C.N.R.S. 195,
Institut de Physique du Globe de Paris et Université Paris 7, France

## Abstract

It has recently been proposed by several authors that stress distribution around active faults may become critical spontaneously. Other authors believe that heterogeneity is a permanent feature of fault planes. We test these ideas on a simple but realistic fault model in the presence of non-linear rate-dependent friction. We find that if friction increases with decreasing slip rate, slip becomes unstable at low slip rates generating supersonic healing phases that lock slip prematurely. Locking of slip in turn produces stress heterogeneity. For a fault model containing a single localized asperity, Das and Okubo found that rupture starts at the asperity and propagates until it either encounters a strong barrier or the stress intensity reduces below a minimum level. Rupture in these models is completely controlled by the physics of the rupture front. Rate dependent friction changes the behavior of the fault in a fundamental way: friction can lock the fault prematurely generating supersonic healing phases. Unlike stopping phases produced by barriers, a healing phase is not a wave phenomenon; it is a direct consequence of the non-linearity of friction. In this case rise time for slip on the fault is no longer controlled by the overall size of the fault as in conventional constant-friction crack models. We find that in our asperity models it is rise time – or the healing mechanism – that controls the final size of the fault. Studying models involving several isolated asperities we find that stress heterogeneity is preserved provided that friction is strongly rate dependent. Depending on the details of the rupture process, the final state of stress on the fault can be quite complex. This behavior is not universal, it depends on the degree of rate-dependence of friction.

**Key words** earthquake – fault dynamics – seismic complexity

## 1. Introduction

A model for the accumulation of stress in the lithosphere and its sudden release by earthquakes was formulated at the beginning of this century by Reid (1910). It took almost 50 years to develop a correct elastodynamic formulation of the problem proposed in the early sixties by Haskell (1964), Maruyama (1963) and Burridge and Knopoff (1964) who laid the basis for a kinematic model of earthquakes. In this model an earthquake is described as a dislocation of the fault, although the most appropriate

term would be a displacement discontinuity or slip across the walls of the fault. Kinematic models were used to calculate the elastic waves generated by slip on the fault and are, still today, at the heart of the techniques devised to invert fault slips from seismic waveforms in the near and far-field, and from geodetic data.

Keilis-Borok (1959) studied a static model of faulting that was the first to attempt to relate the traction change across the fault, or stress drop, to slip. He established the well-known formula that relates a constant stress drop on a circular fault to the average slip on the fault. Between 1964 and 1975, Kostrov (1964, 1966, 1975) laid the foundations of the mathematical methods to study dynamic crack models, i.e. to
study the growth of rupture on a fault due to the loss of traction across an active fault. Kostrov posed the problem in the same way as Reid: earthquakes are due to slow accumulation of stress around a preexisting fault; and they are produced by the sudden triggering of a slip instability across the walls of the fault. Inside the fault, stress drops from the initial prestress to a lower value determined by the kinematic friction between the walls of the fault. Kostrov assumed that after the passage of the rupture front the traction across the fault suddenly relaxed inside the crack to a *known* value of the traction that was constant and independent of slip, slip velocity or thermodynamic state of the fault. Actually, in all these models the final traction was assumed to be zero as in mode I fractures where the final traction is exactly zero. Kostrov's models were extremely successful in providing a proper conceptual basis for the study of earthquake source dynamics but they were concerned with the continuous unilateral or bilateral growth of a rupture front in a uniformly loaded medium. A major defect of these models was that once started, rupture would never stop. The reason is that in a uniformly pre-stressed medium the stress concentration that appears near the rupture front increases indefinitely with the size of the fault. If rupture resistance is also uniform nothing can prevent the fault from running away.

Kostrov *et al.* (1969) studied the details of the energetics of faulting and showed that in order for the rupture to grow at subsonic velocities, stress drop cannot occur abruptly at the rupture front. A certain slip weakening zone has to exist where stress changes continuously as a function of slip. This model had already been proposed by Barenblatt (1964) in his study of hydraulic fractures, and by several other authors in fracture mechanics. Ida and Aki (1972) studied the consequences of simple slip weakening models for the growth of rupture.

Brune (1970) proposed a simple circular crack model in order to model the general properties of the radiation spectrum of far-field seismic waves proposed by Aki (1967) on more empirical grounds. Improving Brune's model, a numerical calculation of a circular fault with uniform stress drop was obtained by

Madariaga (1976) who assumed a simple circular rupture front running at constant velocity until it stops abruptly at a certain radius *a*. The problem of rupture arrest was not considered in this model and it was simply assumed that rupture stopped abruptly at some unbreakable barrier. An important result of these simulations was that slip continued uninterrupted until the arrival of stopping phases from the edge of the circular crack. Thus the duration of slip at any point on the fault was controlled by stopping phases, not by friction or other local properties of the fault. Although the geometry of this model is admittedly simple, it has been useful as a guide for the understanding of seismic radiation from earthquakes.

Soon after that, Das and Aki (1977a,b) introduced a model of rupture containing barriers, i.e. patches of higher rupture resistance on the fault that may slow down rupture or completely stop it. These barriers were supposed to model the intrinsic heterogeneity and discontinuity of faults as observed in the field and in the seismic waves radiated by earthquakes. An alternative model of source heterogeneity was proposed almost simultaneously by Kanamori and his colleagues (Kanamori and Stewart, 1978) who proposed that observed heterogeneity was due to asperities, patches on the fault that have not ruptured in previous seismic events and have therefore accumulated large stress compared to other places on the fault. Madariaga (1979, 1983) studied the radiation from barriers and asperities and concluded that it was not possible to distinguish between these two types of heterogeneities from the observation of seismic waves alone. In fact as has been found by numerous authors, asperities and barriers are complementary concepts that express the fact that faults evolve in a complex way, seismic ruptures propagate along faults that have already experienced previous earthquakes, encounter barriers, slow down, and continue rupturing, etc. Thus we have to consider the complete process of loading, fracturing, and fault interaction in order to understand earthquake complexity. This is a very difficult problem because of some intrinsic problems with modeling seismic ruptures and because a proper understanding of rupture requires study-
ing faulting at several length and time scales: at the shortest scale a proper understanding of the physics of rupture fronts is still lacking because laboratory data like that reported by Dieterich (1972, 1978) were obtained at very slow values of slip rate. At the next scale of complexity a good understanding of the geometry and scale of stress heterogeneity on the fault is needed in order to model observed seismic radiation. Some attempts at solving these more complex problems were made by Day (1982), Mikumo and Miyatake (1983), Das and Kostrov (1988) and many others. Their models are however intrinsically limited by the fact that most numerical methods cannot simultaneously resolve the small scale processes near the rupture front, and the large scale heterogeneities of faulting.

Because of these difficulties, rupture growth has usually been simulated assuming a numerical rupture criterion based on the maximum numerical stress ahead of the rupture front. Virieux and Madariaga (1982) and Koller *et al.* (1992) have illustrated the difficulties of this approach.

A new interest in properly modeling seismic rupture has been triggered by some results obtained by Carlson and Langer (1989) using the block and spring model proposed by Burridge and Knopoff (1967). Using a very simple rate-dependent friction model they found that a very complex state of stress develops in their model starting from slightly perturbed homogeneous conditions. In this model stress heterogeneity and a wide distribution of magnitudes is obtained in a simple homogeneous fault model by the effect of dynamics. These results have been criticized by Rice (1993), who very appropriately remarked that the Burridge-Knopoff model belongs to a category of intrinsically discrete models that lacks a continuum limit. Unfortunately doing a proper analysis of a complex fault in a continuum requires much more accurate numerical methods than those developed in the literature. The first inroads in this direction were made by Andrews (1985) who modeled the propagation of a plane fault under the control of a slip weakening model of friction. Okubo (1989) used this numerical model to study the propagation of a fault with

rate and state dependent friction of the kind proposed by Dieterich (1972). Inspired by this work, we have developed improved Boundary Integral Equation (BIE) methods that would allow us to gain sufficient control of the stress and slip velocity field on the fault in order to study slip and rate dependent friction models. As we will show in this paper, the question is still open whether a properly modeled fault can sustain the kind of heterogeneity observed by Carlson and Langer and evolve towards a regime of self-organized criticality as proposed by Bak and Tang (1989) or Sornette and Sorrette (1989).

## 2. Models of source heterogeneity

In a recent paper Carlson and Langer (1989) reconsidered the block-and-spring model of Burridge and Knopoff (1967) (B-K in the following). In most previous studies of the B-K model – with a few exceptions like Cao and Aki (1986) – an extremely simple frictional law was used. In this law, sometimes called Coulomb friction, it was assumed that friction is linearly proportional to the normal pressure across the fault, but independent of slip or slip velocity on the fault. Instead of this simple friction model, Carlson and Langer (1989) used a rate-dependent friction in their simulations. Starting from a very small initial heterogeneity of stress they observed that heterogeneity developed in a natural way until it reached a self-similar distribution. Slip events (earthquakes) of all sizes were observed on the same stretch of the fault and they found that these events obeyed a Gutenberg-Richter law.

In the Burridge and Knopoff (1967) model each block obeys well posed mechanical equations and reproduces, at least locally, the frictional instability that causes earthquakes. Unfortunately, this model has two severe limitations: it does not radiate seismic waves and it does not include long range elastic forces. Rice (1993) has raised an even stronger objection, namely, that this model may be intrinsically discrete so that the observed heterogeneity may be due to inadequate sampling of the friction law on the fault.
In order to test the applicability of Carlson and Langer's (1989) results to earthquakes it is clearly necessary to study more realistic earthquake models, based on the elastodynamics of rupture propagation. As a first step in that direction, we present here a study of an antiplane fracture model using a general rate-dependent model for friction. To our knowledge, the only previous work on the dynamics of faulting with rate-dependent friction is that of Okubo (1989) who used the friction law proposed by Dieterich (1972). Okubo's simulations did not show any evidence of heterogeneity of the final stress on the fault. Our results - Cochard and Madariaga (1994) - are different: using Carlson and Langer's rate dependent friction we found that the final stress on the fault becomes heterogeneous at least for certain values of the rate-dependence.

A related, and very important, aspect of rupture heterogeneity is the study of possible mechanisms for the generation of short duration pulses observed by Heaton (1990) and confirmed by recent work on the Lander's earthquake by Wald and Heaton (1994), Cohee and Beroza (1994) and Cotton and Campillo (1994). With the help of a very simple model that we can solve exactly under certain circumstances we will show that these short pulses

may arise due to the interplay of rate-dependent friction and the presence of asperities on the fault.

### 3. Mathematical statement of the dynamic faulting problem

We consider a homogeneous elastic body V with a single flat fault that lies along the x-axis. The elastic body is infinite or at least its boundaries are considered to have no influence on slip on the fault. The displacement field $u(x, t)$ in the elastic medium satisfies the elastodynamic wave equation

$$ \rho \frac{\partial^2 u}{\partial t^2} = \mu \nabla^2 u + (\lambda + \mu) \nabla \cdot \nabla u \quad (3.1) $$

without body forces and with zero initial conditions. $\rho$ is mass density, $\lambda$ and $\mu$ are the elastic constants. We also assume that radiation conditions are satisfied at infinity. The only source of dynamic motion in this medium is due to slip across the fault, so that we need to examine carefully the boundary conditions on the fault. For this purpose we define a positive and a negative side of the fault surface (see fig. 1).

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Antiplane fault model studied in this paper. Slip is measured as the displacement discontinuity between the positive and negative sides of the fault. Friction is assumed to be a function of slip velocity in the fault.</figcaption></figure>
The traction on the positive side of the fault is $T = -\sigma_{yz}$ where $\sigma$ is the stress tensor. The minus sign appears here because of the conventional definition of stress.

We assume that before the earthquake starts the fault is subject to a certain pre-stress field $\sigma_0(x, z)$ due to the superposition of loading due to plate motion and of stresses due either to slow internal deformation of the lithosphere, or produced by previous events on the fault or on neighboring faults. Contrary to the assumptions implicit in many fault models, we assume that pre-stress is very heterogeneous due to the previous seismic history of the fault. The initial traction across the fault is thus $T_0(x) = -\sigma_{0, yz}$. This pre-stress field is variable in space because the fault has already suffered previous earthquakes, patches where $T_0$ is large are called asperities by most seismologists. At time $t=0$ rupture starts and a certain zone of the x-axis slips as a consequence of the loss of traction across the fault. Let $\Gamma(t)$ designate the areas of the fault that are slipping at a given time. We call $\Delta u = u^+ - u^-$ the slip or displacement discontinuity across $\Gamma(t)$. On the slipping areas of the fault we assume that traction drops to

$$T(x, t) = T_0(x) + \Delta T(x, t) \quad \text{for } x \in \Gamma(t) \qquad (3.2)$$

where $\Delta T(x, t)$ is the traction change on the fault. Traction on the fault decreases during slip so that seismologists call $\Delta T$ the stress drop. In most previous work on dynamic faulting $\Delta T$ has been assumed to be known. This will not be the case here because we are interested in modeling the effect of non-linear rate and slip dependent friction. Thus for points $x \in \Gamma$, we assume that traction $T$ is a function of both slip $\Delta u$ and slip rate $\Delta \dot{u}$. A discussion of possible friction laws and how to solve the boundary value problem will be presented later in this paper.

On other points on the fault $x \notin \Gamma$, the boundary condition is

$$\Delta u(x, t) = 0 \quad \text{for } x \notin \Gamma(t). \qquad (3.3)$$

The solution of the elastodynamic eq. (3.1) subject to these boundary conditions is a mixed boundary value problem whose solution is very difficult because the rupture zone $\Gamma(t)$ has to be calculated as part of the solution. Except for some simple problems with known slip area $\Gamma(t)$ as a function of time, the only way to obtain solutions to this problem is to use numerical methods. Finite difference, finite elements and Boundary Integral Equations (BIE) have been used in the past. From our previous experience with finite differences (Virieux and Madariaga, 1982) we are convinced that BIE methods are the most accurate method for solving dynamic faulting problems.

### 4. Integral equations for an antiplane fault model

We look for solutions to the problem (3.1), (3.2) and (3.3) by the Boundary Integral Equation (BIE) method. This is relatively easy for the problem at hand because we know the Green function $G(x, z, t) = (2\pi\mu)^{-1}(t^2 - r^2/\beta^2)^{-1/2} H(t-r/\beta)$ in closed form (see, e.g., Aki and Richards, 1980). Here $r = \sqrt{x^2 + z^2}$ is the distance to the source point, and $\beta$ is the shear wave velocity. We start from the classical Betti representation theorem. Displacement inside the elastic body is given by

$$u(x, z, t) = \int_{\Gamma} \int_{0}^{t} \Delta u(\xi, \tau) \Sigma(x, z, \xi, t-\tau) d\tau d\xi \qquad (4.1)$$

where $\Sigma = \mu \partial G/\partial z$ is the $yz$ element of the stress tensor associated with the 2-D Green function $G$. When $z \to 0$, this equation reduces to an identity, so that this version of the representation theorem does not lead to a useful BIE.

In order to get a well-posed integral equation we calculate first the stress change $\sigma_{yz}$ due to slip $\Delta u$ in the slipping parts of the fault, and then we let $z \to 0$. In this limit $\sigma_{yz} \to -\Delta T$, the traction across the fault plane. Calculating
$\Delta T = - \mu \frac{\partial u}{\partial z}$ at $z = 0$ from (4.1) we get the following integral equation:

$$ \Delta T(x, t) = - \frac{\mu}{2\pi\beta^2} \int_{\Gamma} \int_{0}^{\tau_m} \frac{\Delta u(\xi, \tau)}{[(t-\tau)^2 - (x-\xi)^2 / \beta^2]^{3/2}} d\tau d\xi $$

where $\tau_m = \max(0, t - \|x - \xi\| / \beta)$. This looks as a very simple integral equation, unfortunately, it can not be used as it is because it is hyper-singular near the source point, when $\xi \to x$ and $\tau \to t$.

We eliminate these very strong singularities by the method proposed by Koller et al. (1992). Assuming that $\Delta u$ and its derivative with respect to $x$ (the dislocation density) are continuous, Cochard and Madariaga (1994) transformed (4.1) into the following regularized integral equation (i.e., an integral equation with an integrable kernel):

$$ \Delta T(x, t) = - \frac{\mu}{2\beta} \Delta \dot{u}(x, t) - $$

$$ -\frac{\mu}{2\pi} \int_{\Gamma} \int_{0}^{\tau_m} \frac{\sqrt{(t-\tau)^2 - (x-\xi)^2/\beta^2}}{(t-\tau)(x-\xi)} \frac{\partial}{\partial \xi} \Delta \dot{u}(\xi, \tau) d\tau d\xi \quad (4.2) $$

where, $\Delta \dot{u}$ is the slip velocity. The domain of integration of eq. (4.2) is shown in fig. 2.

The two terms in (4.2) have a very simple physical interpretation. The first one represents an instantaneous traction change produced by a corresponding change in slip velocity. If slip velocity $\Delta \dot{u}(x, t)$ is independent of position $x$ the second term cancels out and traction is directly proportional to slip velocity. The ratio $\mu/2\beta$ is the radiation impedance of the fault.

The second term contains all the elastodynamic interactions between different points on the fault. This term contains both long range elastic interactions – inversely proportional to $x - \xi$ – and the local wave field of the wave-front singularities that propagate with the shear wave speed ($\beta$) along the fault. It is interesting to note that just as in the static problems (Aki and Richards, 1980) traction change is related to the gradient of slip not to slip itself. In the theory of dislocations $b = \partial_x \Delta u$ is usually

<figure><img src="image_3.png" /><figcaption>Fig. 2. Domain of integration for the solution to the Boundary Integral Equation (BIE) for the antiplane crack (4.2).</figcaption></figure>
called Burger's vector density. In integral eq. (4.2) the elastodynamic interactions are proportional to the rate of Burger's vector rather than Burger's vector. This is understandable since elastic waves are emitted only by temporal changes in the dislocation distribution. Integral eq. (4.2) is regular everywhere except near $\xi \to x$ where it has to be interpreted in the usual sense of a Cauchy integral. This Cauchy-like singularity is associated with the static stress field of a dislocation as can be shown letting $\Delta u(x, t) \to \Delta u(x) H(t)$ in (4.2).

Let us remark that boundary integral equations of the type (4.2) were introduced in the study of faulting by Burridge (1969) who did not use the form (4.2), but a numerical approximation derived from the more common displacement formulation of the boundary integral equations. The displacement formulation was developed independently by Hamano in an unpublished work, and was extensively studied by Das and Aki (1977a), Andrews (1985) and Kostrov and Das (1989). In the displacement formulation eq. (4.2) is inversed, so that $\Delta u$ is given as an integral function of the change in traction across the fault:

$$ \Delta u(x, t) = -\frac{1}{2\pi\mu} $$

$$ \int_{-\infty}^{\infty} \int_{0}^{\tau_m} \frac{1}{\sqrt{(t-\tau)^2 - (x-\xi)^2/\beta^2}} \Delta T(\xi, \tau) d\tau d\xi. \quad (4.3) $$

In this formulation the unknowns appear both on the left hand side of the equation ($\Delta u$) and under the integral, because $\Delta T$ is not known outside the fault. To solve this equation Das and Aki (1977a) used (4.3) as an integral equation in order to compute $\Delta T$ outside the fault; and then used this value of $\Delta T$ to compute $\Delta u$ inside the fault. Although both formulations are equivalent, we consider that the BIE (4.2) is preferable for the study of non-linear stress-slip rate friction laws, because it does not require computing stress outside the slipping patches $\Gamma$ of the fault.

We can now pose the boundary value problem

for the points $x \in \Gamma(t)$ (inside the crack). Since the integral in (4.2) extends over points inside the crack, the second boundary condition (3.3) is automatically satisfied by our integral equation. The stress boundary condition (3.2) can be easily written using (4.2), so that for points $x \in \Gamma$ we get

$$ T(x, t) = T_0(x) - \frac{\mu}{2\beta} \Delta \dot{u}(x, t) - \frac{\mu}{2\pi} K(x, t) * \frac{\partial}{\partial x} \Delta \dot{u}(x, t) \quad (4.4) $$

where $K(x, t) = \sqrt{t^2 - r^2/\beta^2}/(\alpha t)$ designates the kernel of eq. (4.2), the star $*$ indicates time and space convolution, and $T_0(x)$ is the stress field before rupture starts. Although it is not essential for the solution to the boundary value problem, we can compute the stress field in the locked areas of the fault, $x \notin \Gamma(t)$, as the sum of the preexisting stress level $T_0$ and the additional stress due to slip on the portion of the fault that is already sliding:

$$ T(x, t) = T_0(x) - \frac{\mu}{2\pi} K(x, t) * \frac{\partial}{\partial x} \Delta \dot{u}(x, t) \quad (4.5) $$

that is we remove the instantaneous response. Since the current point $x$ is not included in the convolution, the last term is regular and easy to compute numerically. Given the traction $T$, (4.4) and (4.5) provide a complete set of equations for the computation of stress and slip on the fault plane. Unfortunately, as we discuss in the next section, the solution to the problem may be very difficult if the friction law that relates $T$ to slip and slip-rate is non-trivial.

## 5. Friction models

As remarked in the introduction almost all the numerical models used to study dynamic faulting in a continuum used an extremely simple model of friction that assumed that slip de-
velops on the fault plane once the traction be-
comes larger than a certain threshold $T_u$, usu-
ally called static friction. Once this threshold is
overcome, traction decreases instantaneously to
a constant kinematic friction that is indepen-
dent of slip or slip rate. Inside the fault then
the integral eq. (4.4) reduces to a much simpler
one because $T(x, t)$ is constant equal to the
kinematic friction $T_k$. The solution of eq. (4.4)
is still quite difficult because of the nature of
mixed boundary value problems, but it is nev-
ertheless a linear integral equation relating
stress drop and to slip gradient.

This simple friction law contains a number
of physical inconsistencies that have been dis-
cussed by numerous authors in the literature,
see e.g. Burridge (1973). The most important
problem is that if friction drops instantaneously
to the kinematic friction once slip starts, then
there is no energy dissipation at the rupture
front, and as a consequence the rupture front
either does not grow at all, or it moves at the
terminal velocity – shear wave velocity for an-
tiplane cracks. For all other values of the rup-
ture velocity, stress near the rupture front be-
comes singular violating the assumption that
stress can not be larger than the threshold $T_u$.
This of course is not what is observed in real
earthquakes: careful studies of rupture velocity
from broad band data indicate that average
rupture velocities are about 70-75% of the
shear wave velocity for most earthquakes
which show clear evidence of directivity (see
e.g. Campos *et al.*, 1994).

The solution to this paradox is well known:
a finite slip weakening zone must exist near
the rupture front where a finite amount of en-
ergy is dissipated in order to make the fault
grow. The energy released at the rupture front
is sometimes considered to be an intrinsic ma-
terial property although this has never been
proven for cracks or faults. Some authors – fol-
lowing Aki (1987) – propose that the size of
the process zone is of the order of a few hun-
dred meters or less and it is responsible for
$f_{max}$, the high frequency cut-off observed in
near field acceleration spectra. Other authors
propose that the weakening zone is smaller and
that $f_{max}$ is due to attenuation. In either case the
slip weakening zone is so small that its effects

on seismic radiation are very difficult to ob-
serve.

One may think that since slip weakening is
such a localized phenomenon it could be ne-
glected in a first approximation. This is what is
usually done in fracture dynamics. For the pur-
pose of modeling large scale fracture phenom-
ena, the rupture front reduces to a point and the
maximum stress criterion $T_u$ is replaced by an
energy release rupture criterion. In numerical
modeling of faulting, however, it is not possi-
ble to resolve the stress field near the rupture
front. For this reason, Das and Aki (1977a,b)
proposed a simpler numerical criterion that
they called the Irwin rupture criterion. In this
model the crack grows if the numerically com-
puted stress immediately ahead of the rupture
front is larger than a certain number. As we
showed by careful numerical experiments in
Virieux and Madariaga (1982) and Koller *et al.*
(1992), the Irwin criterion is grid size depen-
dent. The reason is simple: any reasonable
model of the rupture front has to include a
length-scale, that measures the size of the slip-
weakening zone. This length-scale, indepen-
dent of the overall size of the fault, appears ex-
plicitly in the energy release rate rupture crite-
ria, or in the more detailed models of the slip
weakening zone discussed by Ida and Aki
(1972). In numerical computations the length
scale is by default the grid size.

Thus the peak numerical stress criterion be-
longs to the class of intrinsically discrete mod-
els justly criticized by Rice (1993). More accu-
rate models of the slip-weakening zone ahead
of the rupture front were modeled by Andrews
(1985) who showed that the size of the slip
weakening zone reduces in size as the rupture
front accelerates. This effect makes the correct
modeling of the rupture front dynamics even
more difficult.

More realistic friction laws were proposed
by Dieterich (1972), Rice and Ruina (1983)
based on experimental studies of friction of
rocks under low slip rates. They proposed that
friction is a non-linear function of slip velocity
and a number of hidden thermodynamic vari-
ables that describe the *state* of the fault at the
time of the earthquake. To our knowledge the
only numerical model where a rate and state
dependent friction law of the Dieterich type has been used is that of Okubo (1989) who studied several rupture problems. Okubo applied Dieterich's friction model only to low slip velocities; for high slip velocities he assumed a constant friction level as in all previous work on dynamic faulting. Implementation of Dieterich-Ruina models is not easy because of the numerical difficulties for properly modeling processes at widely different time scales. The Dieterich-Ruina friction is important during slip initiation at very low slip rates. Actually, Okubo (1989) found that at high rupture speeds, the Dieterich-Ruina dynamic friction behaves like the more common slip weakening models. If this were to be confirmed by further numerical modeling it would mean that most of the phenomenology of friction with slip weakening carries over to the case of Dieterich-Ruina law.

In the previous models of friction the emphasis

was put on the initial loss of traction of rock surfaces for very low values of slip velocity. In fact, frictional behavior at low speeds has a fundamental role in controlling the initiation of slip and the propagation of a rupture front. But this is not the only aspect of friction that is important for the study of dynamic faulting: the behavior of friction for large values of slip rate is equally important and more determinant of the large scale features of faulting as we will show in the rest of this paper.

## 5.1. Rate-dependent friction and slip instability

In their study of the B-K model, Carlson and Langer (1989) used a slip velocity dependent friction law that was initially proposed by Burridge and Knopoff (1967). This friction model, shown in fig. 3, contains no internal

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Rate-dependent friction law used by Carlson and Langer (1989) and the trajectory of a point of the fault in the stress-slip velocity plane. The dynamic unloading lines labeled by *S* have a slope $-\mu/(2\beta)$. When *S* = *S*<sub>crit</sub> slip becomes unstable and healing occurs instantaneously.</figcaption></figure>
state variable nor scale length:

$$T(\Delta\dot{u}) = T_u \frac{V_0}{V_0 + \Delta\dot{u}} \quad (5.1)$$

where T is the absolute traction, $\Delta\dot{u}$ is the slip velocity, $V_0$ is a reference slip rate that determines the rate of slip velocity weakening in the model and $T_u$ is the traction threshold or static friction.

Friction law (5.1) was originally proposed by Burridge and Knopoff (1967) to simulate friction at high slip rates. As shown by Cochard and Madariaga (1994) for a numerical model, (5.1) is very unstable at low values of slip rate, both during the passage of the rupture front and during slip arrest. In order to understand the nature of the instability of (5.1), let us rewrite our BIE (4.4) in the form

$$T[\Delta\dot{u}(x, t)] = -\frac{\mu}{2\beta}\Delta\dot{u}(x, t) + S(x, t) \quad (5.2)$$

where we introduce

$$S(x, t) = T_0(x) - \frac{\mu}{2\pi} K(x, t) * \frac{\partial}{\partial x} \Delta\dot{u}(x, t)$$

a term that contains the effect of pre-stress and the elastodynamic interactions of the current stress field with past values of slip velocity.

For slipping points, eq. (5.2) is an explicit non-linear algebraic equation for the slip velocity $\Delta\dot{u}$ of that point. This equation may have one, none or several solutions and requires careful analysis. The solution of (5.2) at different instants of the slip cycle may better be discussed by referring to the graphical solution shown in fig. 3. We notice first that the right hand side of eq. (5.2) represents a family of straight lines of slope $-\mu/2\beta$ parameterized by S. The left hand side, on the other hand, is simply the friction law (5.1). From a geometrical point of view, the solution of the integral eq. (5.2) consists in determining the intersection of the straight lines parameterized by S with the friction law.

Let us follow the trajectory of a point of the fault on the stress-slip velocity plane. When the rupture front approaches this point, the sum of interactions S increases rapidly so that the total stress T computed from (4.5) approaches the threshold. Slip starts when T reaches the friction threshold $T_u$. At that point $\Delta\dot{u} = 0$. This is indicated by the diamond in fig. 3 and corresponds to $S = T_u$. At that instant of time, the fault is in unstable equilibrium and traction drops instantaneously to a value determined by the intersection of the line $S = T_u$ with the friction law. This is indicated by the black dot in the same figure. There is no energy dissipation in this process and as a consequence rupture propagates at the shear wave velocity. If rupture velocity were lower than $\beta$, a stress concentration would appear immediately ahead of the rupture front, violating the maximum stress condition $T < T_u$. As already discussed above the solution to this problem, already identified by Rice (1993), is to introduce the small characteristic length of the region where stress drops continuously from its peak value $T_u$ to the kinetic friction at high slip rate. We are aware that the friction model (5.1) is incorrect but as will be shown in the next section using exact solutions, we do not believe that this deficiency of the friction law (5.1) affects our main results qualitatively.

For later times, the fault element under consideration receives information from the other places of the fault so that S changes resulting in a translation of the straight line parallel to itself and consequently to a variation of $\Delta\dot{u}$ along the friction curve (e.g. the circle in fig. 3). There are actually two intersections between the straight line and the friction curve but the only relevant one is that with lower stress.

## 5.2. Slip arrest

As slip velocity decreases, S will eventually reach a critical value when the straight line labeled $S_{crit}$ becomes tangent to the friction law. At this time, this point becomes unstable because for $S < S_{crit}$ there are no real solutions of eq. (5.2). $\Delta\dot{u}$ drops instantaneously to 0
<figure><img src="image_5.png" /><figcaption>Fig. 4. Modified rate dependent friction law used in order to avoid the low slip rate instability of the friction law in fig. 3.</figcaption></figure>

and the fault healing instantaneously. At the same time, the absolute traction T increases to $S_{crit}$ as is indicated by the rectangle in fig. 3.

Thus there are two instabilities in this problem, the first produces a sudden stress drop after rupture initiation, the second during slip healing due to the inability of the fault to continue slipping when friction is increasing with reduced slip-velocity. The latter mechanism was discussed by Heaton (1990). In order to avoid the numerical problems due to the healing instability in fig. 3, we have modified the friction law (5.1) in order to reduce the numerical effects of the instability. We take as shown in fig. 4

slope $-2\mu/\beta$. The advantage of this friction model is that it avoids the instability of healing without changing the qualitative properties of the solution. With this modification, in particular, we avoid the annoying numerical noise produced by the sudden jump of stress at the moment of healing. The noise appears because in the numerical computations the healing of different elements of the fault occurs at discrete instants of time so that two neighboring fault elements that heal at two different time intervals end up having a different state of stress. With the law (5.3), stress increases continuously as slip rate decreases to zero.

$$T(\Delta\dot{u}) = S_{sp} \frac{V_x}{V_x + \Delta\dot{u}} \quad \text{for } \Delta\dot{u} > 0 \quad (5.3)$$

where $V_x$ is chosen so that the derivative $dT/d\Delta\dot{u}$ at the origin is greater than the critical

## 6. The origin of short slip events (Heaton pulses)

In the previous discussion we tried to emphasize that faults under rate-dependent friction are unstable: traction on the fault jumps
abruptly at the beginning of slip. For large slip rates friction reaches a finite limit and the fault slips at almost constant stress. As slip velocity decreases slip becomes unstable again because stress tends to increase with reduced slip rate. A major effect of this stress increase during slip deceleration is that it can produce premature healing of the fault. Recently Heaton (1990) observed that the source time history of several earthquakes, determined from the inversion of near field accelerograms, was very short. He proposed that rate-dependent friction was the cause of this short rise time. In the following we will demonstrate that this is indeed the case but that rate-dependent friction by itself is not enough: narrow pulses develop only if the stress field itself is heterogeneous, in particular if it contains localized asperities before the earthquake starts.

We will demonstrate this using a simple but exact solution. Solving crack problems exactly is very difficult but thanks to earlier work by Kostrov (1966) and Burridge (1969),

Madariaga (1983) found an exact solution to the near field of some antiplane fault models under general time-independent stress drop.

### 6.1. Faulting of an asperity under constant friction

We consider in a first step faulting under constant rate-independent friction. In order to obtain the exact solution we assume that a finite crack appears suddenly at time $t = 0$ along a stretch of length $\ell$ of the fault. After the fault is nucleated we let it grow at constant prescribed rupture speed. Traction inside the rupture zone drops instantaneously to a constant kinematic friction level $T_k$. Although it is possible to solve the problem of an arbitrary rupture velocity we prefer to use a prescribed rupture velocity in order to obtain very simple exact solutions. Let us assume, as shown in fig. 5 that the fault is loaded by an initial stress field such that traction is exactly at the friction

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Model of a simple fault with a localized pre-stress concentration or asperity of width $2a$. Rupture starts abruptly from an initial instantaneous rupture of width $2\ell$ and then grows at constant rupture velocity $v$. In the applications, the kinematic friction was $T_k = 0$.</figcaption></figure>
threshold $T_u$ inside an asperity of width $2a$ near the center of the fault. Outside the asperity the initial stress $T_0 \le T_u$. This simple model which is a slight modification of another one studied by Okubo (1989) represents a crude model of an asperity: a region on the fault where the stress is higher than its surroundings where it has been reduced by previous earthquakes.

For simplicity's sake we assume also that the initial rupture is centered inside the asperity. As shown in fig. 5, at time $t=0$ we then let rupture start abruptly running at constant rupture velocity $v$. As rupture grows bilaterally it breaks through the asperity and keeps growing into the area of lower pre-stress. If the pre-stress $T_0$ outside the asperity is only slightly larger than the final kinematic friction $T_k$, it is well known that the fault will propagate along the fault plane forever and that points near the center of the fault will keep slipping also forever. Unless, of course, the rupture front hits an unbreakable barrier.

We can solve for the stress and velocity fields for this fault as long as the waves emitted by one tip do not reach the other one. For longer times waves diffracted by the rupture front must be taken into account. The computation that we will describe is exact until the arrival of these diffractions, for later times our solution is still a reasonable approximation because diffracted waves are quite small at high rupture velocities. If one wants to make an exact computation, it would be better to resort to numerical solutions.

Stress drop inside the asperity is simply $\Delta T_a = T_u - T_k$ and outside the asperity it is $\Delta T_0 = T_0 - T_k$. The final traction everywhere on the rupture zone is $T_k$.

Assuming that the rupture front moves at constant rupture velocity $v$, the exact velocity field inside the fault is given by

$$ \Delta \dot{u}_i (x, t) = \frac{2\Delta T_a \beta}{\mu} V_{fi} (x - \ell, t, v) + \frac{2\Delta T_a v}{\pi \mu} [V_0 (x - \ell, t, v) - V_d (x - \ell, t, v, a - \ell) + $$

$$ +V_0 (-\ell - x, t, v) - V_d (-\ell - x, t, v, \ell - a)] + \\ + \frac{2\Delta T_0 v}{\pi \mu} [V_d (x - \ell, t, v, 0) + V_d (-\ell - x, t, v, 0)] \quad (6.1) $$

where

$$ V_{fi} = 1 \quad \text{for} \quad -\ell + \beta t \le x \le \ell - \beta t $$

and

$$ V_{fi} = V_i (x - \ell, t, v) - V_i (-\ell - x, t, v) \quad \text{otherwise} $$

with

$$ V_i (x, t, v) = 1 \quad \text{for} \quad x < -\beta t $$

$$ = \left[ 1 - \sin^{-1} \left( \frac{(x + \beta t) + (x - vt)}{(1 + v/\beta) \beta t} \right) \right] \\ \text{for} \quad -\beta t \le x < vt \quad (6.2) $$

and $V_i(x, t, v) = 0$ for $x \ge vt$. We defined the functions

$$ V_0 (x, t, v) = \frac{1}{(1 + v/\beta) \sqrt{vt - x}} \sqrt{x + \beta t} \\ \text{for} \quad -\beta t \le x < vt $$

with $V_0(x, v, t) = 0$ outside this range; and

$$ V_d (x, t, v, a) = \frac{1}{(1 + v/\beta)\sqrt{vt - x}} \sqrt{\frac{v}{\beta}} \sqrt{(x - a) + \beta(t - a/v)} \\ \text{for} \quad a(1 + v/\beta) - \beta t \le x < vt $$

with $V_d(x, t, v, a) = 0$ outside this range. Let us remark that $V_d$ is the radiation from the edge of
the loaded zone located at $x = a$. The expressions for $V_0$ and $V_d$ are the same except for the factor $\sqrt{v/\beta}$ that reduces the radiation from the border of the jump in stress drop. This factor appears whenever there is a jump in stress drop on the fault.

Similarly we can compute the stress field outside the crack

$$ T(x, t) = \frac{2\Delta T_a}{\pi} [T_0(x-\ell, t, v) - T_d(x-\ell, t, v, a-\ell) + \\ +T_0(-\ell-x, t, v) - T_d(-\ell-x, t, v, a-\ell)] + \\ +\frac{2\Delta T_0}{\pi} [T_d(x-\ell, t, v, 0) + T_d(-\ell-x, t, v, 0)] \quad (6.3) $$

where

$$ T_0(x, t, v) = \sqrt{\frac{\beta t - x}{x - vt}} - \tan^{-1} \sqrt{\frac{\beta t - x}{x - vt}} $$ 
for $vt < x \le \beta t$

and $T_0(x, t, v) = 0$ for $x > \beta t$, is the stress concentration ahead of the rupture front produced by the initial instantaneous rupture at time $t = 0$. The other function

$$ T_d(x, t, v, a) = \sqrt{\frac{v}{\beta}} \sqrt{\frac{\beta(t - a/v) - (x - a)}{x - vt}} - \\ - \tan^{-1} \sqrt{\frac{v}{\beta}} \sqrt{\frac{\beta(t - a/v) - (x - a)}{x - vt}} $$

for $vt = x \le \beta t + a(1 - v/\beta)$. $T_d(x, t, v, a) = 0$ for $x > \beta t + a(1 - v/\beta)$, represents the radiation from the stress drop jump at $x = a$. We observe as shown by Madariaga (1983) that ruptures radiate only when there are changes in rupture velocity (the term $T_0$) or when there are abrupt changes in stress drop (the term $T_d$).

## 6.2. Exact solution for a fault with uniform load and constant friction

In order to compare exact solutions to slip velocity on a fault with asperities to the well-known crack models, we compute first the slip velocity field on a fault loaded by a uniform pre-stress field. In order to use the analytical solutions of the previous section we assume that rupture starts from a finite patch of size $2\ell = 2a$ that fails instantaneously at time $t = 0$, and then grows at constant rupture velocity $v$. The only difference between this model and the usual self-similar crack model of Kostrov (1964) is that here rupture starts from a finite patch instead of a point. Since initial stress and friction are uniform, referring to fig. 5 we choose $T_u = T_0 = 1$, $T_k = 0$. Rupture velocity is taken as $0.75\beta$, and the results are normalized by $\mu = 1$, $\beta = 1$ and the size of the initial rupture $a = \ell = 1$.

As shown in fig. 6, for a uniformly loaded fault, slip velocity is strongly concentrated near the rupture front and remains very large at all times. Slip velocity is infinite at the rupture front, in order to represent it in fig. 6 we smoothed the singularity near the rupture front. This is just a local correction that does not affect results. Once the fault has grown large enough, slip velocity near the center of the crack has lost the memory of the finite initial rupture and looks almost like that of Kostrov's self-similar crack. In this model slip velocity near the crack center is of the order of $2\beta\Delta T_a/\mu = 2$ which is the slip velocity for an instantaneous rupture of the fault over the entire x-axis. Slip in this model never locks and the rupture extends indefinitely along the fault plane. This is the well known defect of classical fracture models: once a crack starts in a uniformly loaded medium, nothing can stop it except for an unbreakable barrier.

## 6.3. Exact solution for an isolated asperity with constant friction

We consider now a simple model of a localized asperity. We model the asperity as shown in fig. 5. A highly stressed region or asperity of length $2\ell = 2a = 4$ is ready to break at time
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Slip velocity as a function of position and time for a uniformly loaded fault. Rupture occurs instantaneously on a central patch and then grows at constant rupture velocity $v = 0.75\beta$. Pre-stress is uniform and kinematic friction is independent of slip rate. In order to prepare this figure we smoothed the slip velocity singularity near the rupture front.</figcaption></figure>

$t=0$. Stress on the asperity is taken as $T_u = 1$. The asperity is surrounded by a region where pre-stress $T_0 = 0.2T_u$ and without loss of generality we assume that $T_k = 0$, so that $\Delta T_a = T_u$ and $\Delta T_0 = 0.2T_u$. Rupture velocities are normalized by $\beta = 1$ and stresses by $T_u/\mu = 1$. Slip velocities are scaled by $T_u/\mu\beta = 1$.

Figure 7 represents the evolution of slip velocity for the single asperity model when friction is not rate-dependent. As in the previous figure slip velocity has been locally smoothed near the rupture front in order to reduce the singularity. Rupture of the asperity occurs instantaneously. After the initial rupture of the asperity, we assume that rupture propagates at a constant velocity $v = 0.75\beta$. As observed in

fig. 7 at the beginning slip velocity in the asperity is high and uniform. Once the rupture extends beyond the asperity, slip velocity near the center of the fault decreases continuously. For later times, two slip velocity events propagate trailing the rupture fronts and velocity near the center of the fault decreases steadily towards a minimum slip velocity of the order of $2\beta\Delta T_0/\mu$, i.e. the slip velocity for an instantaneous rupture of a fault initially loaded with a pre-stress field of intensity $\Delta T_0$. As long as $\Delta T_0$ is significantly less than $\Delta T_a$, slip velocity near the center of the fault decreases faster than near the borders. The slip velocity field in fig. 7 is completely different from that of fig. 6; for the asperity model slip velocity consists
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Slip velocity as a function of position and time for a fault model containing a single pre-loaded asperity. Pre-stress in the asperity equals $T_u$, the peak frictional stress. Outside the asperity $T_0 = .20T_u$. Rupture occurs instantaneously on a central patch and then grows at constant rupture velocity $v = 0.75\beta$. Kinematic friction was assumed to be independent of slip rate. In order to prepare this figure we smoothed the slip velocity singularity near the rupture front.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Stress as a function of position and time for a fault model containing a single pre-loaded asperity. Pre-stress in the asperity equals $T_u$, outside the asperity $T_0 = .20T_u$. Rupture occurs instantaneously on a central patch and then grows at constant rupture velocity $v = 0.75\beta$. Kinematic friction was assumed to be independent of slip rate. Stress was smoothed in the vicinity of the rupture front in order to avoid singularities in the figure.</figcaption></figure>
in two propagating velocity pulses of width comparable to the size of the initial asperity (2a). The amplitude of these rupture pulses decreases very quickly after rupturing through the asperity, but at longer distances their amplitude increases slowly. These velocity pulses are still rather large because the friction law used here was the classical rate-independent friction model.

In fig. 8 we present the stress field for the asperity model. As shown in this figure, inside the asperity stress drops instantaneously from $T_u = 1$ to $T_k = 0$ and remains constant forever because of the constant kinematic friction assumed in this model. Outside the asperity stress drop is lower, changing from an initial value $T_0 = 0.20$ down to zero after the passage of the rupture front. Since we assumed a constant subsonic rupture velocity we observe a stress concentration ahead of the rupture front.

## 6.4. Exact solution for rate-dependent friction

We will compare the solution of the previous section for constant kinematic friction $T_k = 0$ with solutions for a particular rate dependent friction. Ideally we would like to solve for Carlson and Langer's friction (5.1). But we cannot solve these non-linear problems exactly. For this reason we adopt a much simpler friction law, in which friction is constant for all slip velocities higher than a threshold $V_0$. When the slip velocities $\Delta\dot{u}$ decrease below $V_0$, slip on the fault suddenly locks or heals. Thus we write the non-linear friction law

$$T = T_k \quad \text{for} \quad \Delta\dot{u} \ge V_0 \quad (6.4)$$

where $T_k$ is a constant kinematic friction for slip rates larger than $V_0$. When slip rates are lower than $V_0$, the fault is assumed to heal instantaneously. This sudden healing simulates the healing instability produced by the friction law (5.1) when slip velocity is close to the critical value shown by the diamond in fig. 3. As discussed previously, because of the first term in (4.2), when the fault heals stress increases instantaneously by the amount $\Delta T = \mu/(2\beta)V_0$.

Although this friction law is non-linear, it is possible to obtain an exact solution because we can find the instant of time when slip locks for simple fault models. In particular we can obtain an exact solution for the simple asperity model discussed in the previous section. As shown in fig. 7 for constant kinematic friction, slip velocity decreases rapidly behind the rupture pulse. Slip velocity first reaches the critical slip rate $V_0$ at the center of the fault. The $V_0$ slip velocity isoline then spreads from the center of the fault at a supersonic apparent velocity. It is this property that allows us to solve for the slip velocity for the friction law (6.4): because healing propagates supersonically different points on the fault do not interact, i.e. waves from the initial healing point arrive at other fault points after they have locked in turn. This shows that slip arrest is strictly a local phenomenon due to the nonlinear rate dependence of the friction. In fig. 9 we show the effect of rate-dependent friction (6.4) on the slip velocity field for the single asperity model. This figure can be compared to fig. 6 where friction was considered to be constant for all values of slip rate. The main difference between these two figures is that slip velocity has locked prematurely due to the healing mechanism included in (6.4). The slip velocity field consists now in two narrow slip velocity pulses propagating bilaterally, and leaving behind a central region that has slipped and healed in a relatively short time controlled by the limiting slip rate $V_0$. The corresponding stress field is presented in fig. 10. Stress is computed exactly outside the fault and is frozen at a constant value inside the fault after healing. Unfortunately stress inside the fault cannot be computed exactly after healing.

The most striking result of fig. 9 is the narrow rupture event running behind the rupture front. This narrow pulse is similar to those observed by Heaton (1990) in the modeling of several recent Californian earthquakes. Because the healing event is supersonic, the slip pulses which propagate at the subsonic rupture velocity become increasingly narrower with time. Eventually the healing pulse catches up with the rupture front producing the complete arrest of rupture without intervention of any
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Exact solution for slip velocity as a function of position and time for a fault model containing a single asperity under rate-dependent friction. Pre-stress in the asperity is equal to $T_u$, the peak frictional stress. Outside the asperity $T_0 = .20T_u$. Rupture occurs instantaneously on a central patch and then grows at constant rupture velocity $v = 0.75\beta$. Slip velocity was smoothed in the vicinity of the rupture front in order to avoid singularities in the figure.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Exact solution for stress as a function of position and time for a fault model containing a single asperity under rate-dependent friction. Pre-stress in the asperity equals $T_u$, outside the asperity $T_0 = .20T_u$. Rupture occurs instantaneously on a central patch and then grows at constant rupture velocity $v = 0.75\beta$. Stress was smoothed in the vicinity of the rupture front in order to avoid singularities in the figure.</figcaption></figure>
unbreakable barrier. Thus we have found a mechanism for stopping seismic ruptures that does not depend on the presence of barriers on the fault. Although law (6.4) is a particular case of slip-rate weakening friction, the results we have obtained apply to most slip-rate weakening friction laws as will be shown in next section using numerical models. Our results are exact, therefore they do not suffer from the limitations of other computations done with so-called intrinsically discrete models. The rupture front in our example has the correct inverse-square root singularity of all crack models of faulting. Ideally one would like to model a slip-weakening or slip-rate and state dependent rupture front but there is no way at present to obtain exact solutions to such a problem. In reality, of course, the rupture front will not run at constant speed, but since rupture will be subsonic anyway, the healing phase will always catch up with the rupture front, producing a narrow Heaton pulse and eventually stopping the rupture front.

## 7. Numerical solution to the BIE on the fault

In order to discretize the integral equation (4.2) we proceed as in Cochard and Madariaga (1994). First, we introduce the following discretization of the slip velocity field:

$$ \Delta \dot{u}(x, t) = \sum_{j, m} V_{j, m} d(x, t; x_j, t_m) $$

where $d(x, t; x_j, t_m)$ is the simple box-car function:

$$ d(x, t; x_j, t_m) = \begin{cases} 1 & \text{if } x_j \le x < x_{j+1} \\ 0 & \text{and } t_m \le t < t_{m+1} \end{cases} $$

$$ d(x, t; x_j, t_m) = 0 \quad \text{otherwise.} $$

We obtain a matrix equation by collocation of the integral equation at a series of knots located inside each boundary element. These collocation points are defined by the coordinates

$x_i + 1/2 \Delta x$ and times $t_n + \epsilon \Delta t$ with $0 < \epsilon < 1$. In most applications we use $\epsilon = 0.999 < 1$. Writing $T_{i,n} = T(x_i + 1/2 \Delta x, t_n + \epsilon \Delta t)$, we get the following discrete boundary integral equation:

$$ T_{i, n}(V_{i, n}) = -\frac{\mu}{2\beta} V_{i, n} + S_{i, n} \quad (7.1) $$

where we have written

$$ S_{i, n} = T_{0i} - \frac{\mu}{2\beta} \sum_{m=0}^{n-1} \sum_{j} V_{j, m} K_{i-j, n-m} \quad (7.2) $$

Since $T_{0i}$ is the initial traction of element $i$, the term $S_{i, n}$ contains the pre-stress and the elasto-dynamic interactions of the current stress field with past values of slip velocity.

In (7.1) $T_{i, n}(V_{i, n})$ is the friction law that relates the total stress at the collocation point $i$ to the slip rate $V$ of the element. In order to solve the integral equation, we distinguish as in eqs. (4.4) and (4.5) between boundary elements that are in the process of slipping, from boundary elements that are locked, for which either the rupture front has not yet arrived or that have already slipped and healed. For non-slipping elements we set $V_{i, n} = 0$ in eq. (7.1), and we use it to compute stress in the element.

## 8. Numerical results

In order to verify that the numerical results for other friction laws are qualitatively the same as those obtained with the exact solution for the peculiar friction law (6.4), we will first show a numerical solution for the rupture of a single isolated asperity. In Cochard and Madariaga (1994) we already presented numerical solutions for this problem using the friction law (5.1). Unfortunately, those simulations presented some numerical noise produced by the healing instability already discussed with reference to fig. 3. In order to avoid numerical noise due to healing, we use the modified friction law (5.3) in which we choose the low slip
<figure><img src="image_12.png" /><figcaption>Fig. 11. Two models of pre-stress on an antiplane fault containing a single or a twin pre-loaded asperity. The rest of the fault is supposed to be free from pre-stress.</figcaption></figure>

rate behavior so as to avoid the healing insta-
bility. In the first numerical model we study
the rupture of the single isolated asperity
shown in fig. 11.

8.1. *Single asperity model*

For the simulation we used a single asperity
containing 100 elements. For clarity, the slip
velocity history presented in fig. 12 shows only
one line every ten grid elements. In all the
computations, the non-dimensional threshold
T<sub>u</sub>/μ was fixed at 5.0. We used the friction law
defined in (5.3) with S<sub>sp</sub>/(μ) = 1.1 and V<sub>x</sub> =
0.75βS<sub>sp</sub>/μ. In the numerical results slip veloc-
ity was normalized by βT<sub>u</sub>/μ, stress by μ and
lengths are given in units of grid size.

As shown in the slip velocity plot of fig. 12,
healing starts spontaneously from the center of
the fault, just as in the exact calculations of
fig. 9 using the rate limited friction law (6.4).
Contrary to the exact solution, however, slip
velocity in the numerical simulation decreases
continuously to zero, eliminating the jump in
slip rate across the healing phase. Healing
spreads supersonically from the center of the
fault until it catches up with the rupture front.
Once the healing phase reaches the rupture
front, rupture stops spontaneously. No barrier
or local variation in stress field is needed, rup-
ture stops simply as a consequence of the en-
ergy dissipation produced by local rate-depen-
dent healing. The overall result is qualitatively
the same as for the exact solutions. This con-
firms the existence of the self-healing pulses
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Numerical solution for the slip velocity as a function of position and time produced during the rupture of a fault containing a single localized asperity under rate dependent friction (5.3).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Numerical solution for the stress as a function of position and time produced during the rupture of a fault containing a single localized asperity under rate dependent friction (5.3).</figcaption></figure>
proposed by Heaton (1990), at least for the particular case of faults loaded by isolated asperities and for rate-dependent friction.

The stress field inside fault shown in fig. 13 is somewhat different from that obtained in the exact solution of fig. 10. Inside the fault stress changes continuously, increasing as slip rate decreases. This smooth variation continues after the slip arrest. The net effect is an apparent supersonic healing phase which is identified in fig. 13 by a change in slope of the time variation of the stress field. As expected the stress field is quite smooth after healing, although some low amplitude oscillations are still apparent. The oscillations appear because healing jumps from cell to cell in the numerical simulations producing small but perceptible differences in the final stress field. This numerical noise is much weaker than in our previous simulations using unstable friction (5.1) (Cochard and Madariaga, 1994). When the healing phase approaches the rupture front, stress after healing becomes noisy because of the interference of numerical noise in the slip velocity field near the rupture front with the healing phase. Qualitatively, however, the exact results and the numerical ones confirm the generation of short duration pulses by rate dependent friction in faults containing localized stress heterogeneities.

## 8.2. Twin-asperity model

A much more interesting model is a finite fault pre-loaded by two asperities as shown at the bottom of fig. 11. The first asperity is assumed to be at the rupture threshold $T_u$, and the second one is almost at the rupture threshold so that it is ready to break when the rupture arises. These asperities are 100 and 50 elements large, respectively, and are separated by a low pre-stress segment containing 100 elements. Since in this model there is no way to stop rupture, we have assumed that two unbreakable barriers are located at the two ends of the fault at a distance of 250 grid elements from the center. The non-dimensional threshold $T_u/\mu$ was fixed at 5.0 and, in the first simulation, we assumed constant kinematic friction $T_k = 0$.

Stress outside the asperities was exactly 0. Slip velocity normalized by $\beta$ is shown in fig. 14 and the corresponding stress field normalized by $\mu$ in fig. 15. As seen from these two figures the main effect of rupture of the two asperities is to transfer the pre-stress accumulated in the asperities towards the borders of the fault. The stress concentrations in fig. 15 are due to this stress transfer. The stress transfer from the asperities to the stress concentration around the borders of the fault is amplified by stress overshoot near the center of the asperities. This overshoot is clearly seen in fig. 15 in the form of two narrow V-shaped «channels» in stress field. Overall duration of slip near the center of the fault is clearly controlled by the stopping phases that are clearly seen crossing the fault after rupture has stopped at the ends of the fault.

An interesting observation is that the rupture of the second asperity produces a slip pulse that crosses the entire fault and re-triggers slip near the abscissa $x = -200$. Thus even for simple friction pre-stress heterogeneity produces new phenomena that were not included in classical crack models.

In the following simulation we studied again the twin asperity model, but this time we used the friction law defined in (5.3) with $S_{sp}/\mu = 1.8$ and $V_x = 1.5\beta S_{sp}/\mu$. Slip velocity is normalized by $\beta T_u/\mu$, stress by $\mu$. $T_u/\mu = 5$ in these simulations. The evolution of slip velocity and traction in the presence of rate-dependent friction are shown in figs. 16 and 17, respectively.

The most important differences between this simulation and the previous one with constant friction is that the global arrest of movement is spontaneous and due only to the friction law. Rupture is completely contained inside the boundaries of the fault. The healing phases apparent in fig. 16 catch up with the rupture front and stop the rupture before it reaches the boundaries of the fault at $x = \pm 250$. Another important difference with fig. 14 is that slip duration is no longer controlled by the overall size of the fault but by the friction law itself. This can occur in a very complex way, as seen in fig. 16 rupture of the second asperity produces a new rupture event that propagates for a
<figure><img src="image_15.png" /><figcaption>Fig. 14. Numerical solution for the slip velocity as a function of position and time produced during the rupture of a fault containing two localized asperities under constant rate-independent kinematic friction.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15. Numerical solution for the stress as a function of position and time produced during the rupture of a fault containing two localized asperities under constant rate-independent kinematic friction.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Numerical solution for the slip velocity as a function of position and time produced during the rupture of a fault containing two localized asperities under rate dependent friction (5.3).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Numerical solution for the stress as a function of position and time produced during the rupture of a fault containing two localized asperities under rate dependent friction (5.3).</figcaption></figure>
while into the already broken parts of the fault, regenerating slip. In fig. 14 this event crosses the whole fault, here it is arrested by friction before it has the time to reach the center of the fault.

Due to the interference between the slip histories of two asperities, the final stress field at the end of rupture presents a very strong heterogeneity. The initial stress localized in the asperities has been transferred to several places along the fault. In fact the stress field inside the fault has decreased much less than in the case of constant kinematic friction. Two strong stress concentrations now appear near $x = 50$ and $x = -100$, both are due to the arrest of slip events emanating from each of the asperities. The stress field is weaker in the direction towards the second asperity. It looks as if the result of rupture was to transfer stress into the less loaded areas of the fault. The details of this stress transfer are strongly friction-law dependent. They demonstrate, however, that depending on the details of pre-stress distribution and the friction law, the final state of stress may be very complex. This is the main point we wanted to stress in this work.

## 9. Discussion and conclusion

From exact and numerical solution of a model of faulting with a single initial stress asperity, we conclude that faulting under rate-dependent friction presents a number of fundamental differences with respect to the well known models of propagation of a fault under constant friction. The main new result is that rate-dependent friction makes fault healing unstable, so that slip velocity drops abruptly to zero on the arrival of any stopping phase issued from a barrier or the edges of the fault. Fault healing spreads from the center of the fault forming a healing phase, the hyperbolic curve where the slip velocity jumps to zero in fig. 9. This phase propagates bilaterally at supersonic speed and finally catches up with the rupture front. Thus the whole rupture process may stop spontaneously before the rupture has had the time to reach the edges of the fault. These results were first obtained by Cochard

and Madariaga (1994) using the very unstable friction model (5.1). We have found that the supersonic healing phases appear also for other friction laws and we showed this exactly for the friction model (6.4).

Similar results were found numerically for the friction model (5.3), a variation of (5.1) that has the advantage that slip velocity and traction vary continuously at healing. The unwanted numerical oscillations observed by Cochard and Madariaga (1994) have practically disappeared. The only unsatisfactory element that remains is the lack of length scale of the friction criterion. The main effect of this deficiency is that rupture velocity depends on the grid size. Although we are convinced that a proper treatment of the rupture front would not qualitatively change the results presented in this paper, an improved friction law will be introduced in future work.

From the solution of the slightly more complex case of two interacting stress asperities, we find that stress heterogeneity may be spontaneously maintained on the fault plane. This is due to the instability of healing already mentioned above. This result is different from that of Okubo (1989) who found in his simulations that, for a single asperity, the state of stress inside the fault at the end of rupture was very smooth. This result is probably due to the use of a constant friction beyond a limiting slip velocity. It is not clear from his paper what friction law was used during the healing process. This point should also be the object of future work.

We must emphasize that our fault model contains no material heterogeneities at all. The rupture threshold is homogeneous along the fault axis (except at the ends of the fault which in fact do not play any part in the rupture history when the friction is rate-dependent) and so are the characteristics of the friction law. What we call an asperity is not a fault patch with a higher rupture threshold as is sometimes assumed, it is just an area where the *initial stress* is high compared with the surroundings. Presumably, this stress heterogeneity has been left from a previous rupture event. It is clear that, had we included heterogeneities of the rupture threshold, we could have also observed short
rise times and heterogeneous distributions of
traction after the rupture as Das and Kostrov
(1988) did.

We cannot quantitatively compare our two-
dimensional results with real data. However,
our numerical simulations are in qualitative
agreement with Heaton's (1990) observation of
short slip velocity pulses. He attributed this to
the non-linear rate-dependence of friction, just
as in our examples. The rate dependence of
friction is not sufficient to produce short slip
pulses, an additional mechanism is needed in
order to reduce slip velocity fast enough. No
level of rate dependence will produce short
pulses in the homogeneous stress model pre-
sented in fig. 6. In our models, slip velocity de-
creases because of the localized loading of the
fault at the asperities. Indeed, in the single as-
perity model of figs. 7 and 9, slip velocity nat-
urally decreases even with constant friction.
The decrease in slip rate depends on the value
of the pre-stress outside the asperity: the lower
the pre-stress outside the asperities, the more
rapidly slip velocity decreases. Thus, rise time
depends not only on the amount of rate-depen-
dence of friction, but also on the difference be-
tween the pre-stress inside and outside the as-
perities. Ideally one would like to quantify this,
but friction being non linear, there is no simple
way to relate the pulse width to asperity size
and rate-dependence of friction. For example,
one can have a relatively high rate-dependence
in the friction law, and yet not be able induce
the early arrest of slip if the pre-stress outside
the asperity is not too low.

Results for the twin asperity model show
that the final state of stress after rupture is not
simply related to the initial stress as is usually
assumed. Clear dynamic effects are observed
in figs. 16 and 17. The final state of stress
after an event plus the loading due to plate
motion should serve as the initial state of the
following event on the fault. Will the stress
distribution become spontaneously heteroge-
neous as in the box spring models of Carlson
and Langer? Preliminary results show that this
is actually the case if friction is highly rate-de-
pendent. This will also be the subject of further
work.

Acknowledgements

This work was sponsored by Groupement
Scientifique Sismique, a joint project of
CNRS, IFP, and Elf Aquitaine. All computa-
tions were carried out in the Connexion Ma-
chine CM5 of the Centre National de Calcul
Parallèle en Sciences de la Terre. We thank
Patrick Stoclet of this center for his invaluable
help in developing the parallel FFT algo-
rithms.

REFERENCES

AKI, K. (1967): Scaling law of seismic spectrum, *J. Geophys. Res.*, **73**, 5359-5376.

AKI, K. (1987): Magnitude-frequency relation for small earthquakes: a clue to the origin of $f_{max}$ of large earthquakes, *J. Geophys. Res.*, **92**, 1349-1355.

AKI, K. and P.G. RICHARDS (1980): *Quantitative Seismology* (W.H. Freeman, New York).

ANDREWS, D.J. (1985): Dynamic plane-strain shear rupture with a slip-weakening friction law calculated by a boundary integral method, *Bull. Seismol. Soc. Am.*, **75**, 1-21.

BAK, P. and C. TANG (1989): Earthquakes as a self-organized critical phenomenon, *J. Geophys. Res.*, **94**, 15,635-15,637.

BARENBLATT, G. I. (1964): On some general concepts of the mathematical theory of brittle fracture, *J. App. Math. Mech.*, **28**, 778-792.

BRUNE, J.N. (1970): Tectonic stress and the spectra of seismic shear waves from earthquakes, *J. Geophys. Res.*, **75**, 4997-5009.

BURRIDGE, R. (1969): The numerical solution of certain integral equations with non-integrable kernels arising in the theory of crack propagation and elastic wave diffraction, *Phil. Trans. R. Soc.*, **A265**, 363-381.

BURRIDGE, R. (1973): Admissible speeds for plane-strain self-similar shear cracks with friction but lacking cohesion, *Geophys. J. R. Astron. Soc.*, **35**, 439-455.

BURRIDGE, R. and L. KNOPOFF (1964): Body force equivalents for seismic dislocations, *Bull. Seismol. Soc. Am.*, **54**, 1,875-1,888.

BURRIDGE, R. and L. KNOPOFF (1967): Model and theoretical seismicity, *Bull. Seismol. Soc. Am.*, **57**, 341-371.

CAMPOS, J., R. MADARIAGA, J. NABÈLEK, B. BUKCHIN and A. DESCHAMPS (1994): Faulting process of the 20 June 1990 earthquake from broadband records, *Geophys. J. Int.*, **118**, 31-46.

CAO, T. and K. AKI (1986): Seismicity simulation with a rate- and state-dependent friction law, *Pure Appl. Geophys.*, **124**, 487-513.

CARLSON, J.M. and J.S. LANGER (1989): Mechanical model of an earthquake fault, *Phys. Rev. A*, **40**, 6470-6484.

COCHARD, A. and R. MADARIAGA (1994): Dynamic faulting under rate-dependent friction, *Pure Appl. Geophys.*, **142**, 419-445.
COHEE, B.P. and C.BEROZA (1994): Slip distribution of the 1992 Landers earthquake and its implications for earthquake source mechanics, *Bull. Seismol. Soc. Am.*, **84**, 692-712.

COTTON, F. and M. CAMPILLO (1994): Frequency domain inversion of strong motions: application to the 1992 Landers earthquake, *Pure Appl. Geophys.* (in press).

DAS, S. and K. AKI (1977a): A numerical study of two-dimensional spontaneous rupture propagation, *Geophys. J. R. Astron. Soc.*, **50**, 643-668.

DAS, S. and K. AKI (1977b): Fault plane with barriers: a versatile earthquake model, *J. Geophys. Res.*, **82**, 5658-5670.

DAS, S. and K. KOSTROV (1988): An investigation of the complexity of the earthquake source time function using dynamic faulting models, *J. Geophys. Res.*, **93**, 8035-8050.

DAY, S.M. (1982): Three-dimensional finite difference simulation of fault dynamics: rectangular faults with fixed rupture velocity, *Bull. Seismol. Soc. Am.*, **72**, 705-727.

DIETERICH, J.H. (1972): Time-dependent friction as a possible mechanism for aftershocks, *J. Geophys. Res.*, **77**, 3771-3781.

DIETERICH, J.H. (1978): Time-dependent friction and the mechanics of stick-slip, *Pure Appl. Geophys.*, **116**, 790-806.

HASKELL, N.A. (1964): Total energy and energy spectral density of elastic wave radiation from propagating faults, *Bull. Seismol. Soc. Am.*, **54**, 1,811-1,841.

HEATON, T.H. (1990): Evidence for and implications of self-healing pulses of slip in earthquake rupture, *Phys. Earth Planet. Inter.*, **64**, 1-20.

IDA, Y. and K. AKI (1972): Seismic source time function of propagating longitudinal-shear cracks, *J. Geophys. Res.*, **77**, 2034-2044.

KANAMORI, H. and G.S. STEWART (1978): Seismological aspects of the Guatemala earthquake of February 4, 1976, *J. Geophys. Res.*, **83**, 3427-3434.

KEILIS-BOROK, V.L. (1959): On the estimation of the displacement in an earthquake source and of source dimensions, *Ann. Geofis.*, **12**, 205-214.

KOLLER, M.G., M. BONNET and R. MADARIAGA (1992): Modeling of dynamical crack propagation using regularized time-domain boundary integral equation, *Wave Motion*, **16**, 339-366.

KOSTROV, B.V. (1964): Self-similar problems of propagation of shear cracks, *J. App. Math. Mech.*, **28**, 1077-1087.

KOSTROV, B.V. (1966): Unsteady propagation of longitudinal shear cracks, *J. App. Math. Mech.*, **30**, 1241-1248.

KOSTROV, B.V. (1975): On the crack propagation with variable velocity, *Int. J. Fract.*, **11**, 47-56.

KOSTROV, B.V. and S. DAS (1989): *Principles of Earthquake Source Mechanics* (Cambridge University Press, Cambridge, U.K.).

KOSTROV, B.V., L.V. NIKITIN and L.M. FLITMAN (1969): The mechanics of brittle fracture, *Mekhanika Tverdogo Tela*, **4**, 112-125.

MADARIAGA, R. (1976): Dynamics of an expanding circular fault, *Bull. Seismol. Soc. Am.*, **66**, 639-666.

MADARIAGA, R. (1979): On the relation between seismic moment and stress drop in the presence of stress and strength heterogeneity, *J. Geophys. Res.*, **84**, 2243-2250.

MADARIAGA, R. (1983): High frequency radiation from dynamic earthquake fault models, *Annales Geophysicae*, **1**, 17-23.

MARUYAMA, T (1963): On the force equivalents of dynamic elastic dislocations with reference to the earthquake mechanism, *Bull. Earthquake. Res. Inst., Tokyo University*, **41**, 467-486.

MIKUMO, T. and T. MIYATAKE (1983): Numerical modeling of space and time variations of seismic activity before major earthquakes, *Geophys. J. R. Astron. Soc.*, **74**, 559-583.

OKUBO, P.G. (1989): Dynamic rupture modeling with laboratory-derived constitutive relations, *J. Geophys. Res.*, **94**, 12321-12335.

RICE, J.R. (1993): Spatio-temporal complexity of slip on a fault, *J. Geophys. Res.*, **88**, 9885-9907.

RICE, J.R. and A.L. RUINA (1983): Stability of steady frictional slipping, *J. App. Mech.*, **50**, 343-349.

REID, H.F. (1910): The mechanics of the earthquake, in *The California Earthquake of April 18, 1906, Report of the State Investigation Commission*, edited by CARNEGIE INSTITUTE OF WASHINGTON (Washington D.C.), vol. 2.

SORNETTE, D. and A. SORNETTE (1989): Self-organized criticality and earthquakes, *Europhys. Lett.*, **9**, 197-202.

VIRIEUX, J. and R. MADARIAGA (1982): Dynamic faulting by a finite difference method, *Bull. Seismol. Soc. Am.*, **72**, 345-369.

WALD, D.J. and T.H. HEATON (1994): Spatial and temporal distribution of slip for the 1992 Landers, California earthquake, *Bull. Seismol. Soc. Am.*, **84**, 668-691.