published in: Journal of Geophysical Research, vol. 105, pp. 23,765-23,789, 2000
LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
To appear in the Journal of Geophysical Research, 2000. 1

Elastodynamic analysis for slow tectonic loading with spontaneous rupture episodes on faults with rate- and state-dependent friction
Nadia Lapusta and James R. Rice
Division of Engineering and Applied Sciences and Department of Earth and Planetary Sciences, Harvard University, Cambridge, Massachusetts

Yehuda Ben-Zion
Department of Earth Sciences, University of Southern California, Los Angeles, California

Gutuan Zheng
Technology Group, IBM, Somers, New York

Abstract. We present an efficient and rigorous numerical procedure for calculating the elastodynamic response of a fault subjected to slow tectonic loading processes of long duration within which there are episodes of rapid earthquake failure. This is done for a general class of rate- and state-dependent friction laws with positive direct velocity effect. The algorithm allows us to treat accurately, within a single computational procedure, loading intervals of thousands of years and to calculate, for each earthquake episode, initially aseismic accelerating slip prior to dynamic rupture, the rupture propagation itself, rapid post seismic deformation which follows, and also ongoing creep slippage throughout the loading period in velocity-strengthening fault regions. The methodology is presented using the two-dimensional (2-D) antiplane spectral formulation and can be readily extended to the 2-D in-plane and 3-D spectral formulations and, with certain modifications, to the space-time boundary integral formulations as well as to their discretized development using finite difference or finite element methods. The methodology can be used to address a number of important issues, such as fault operation under low overall stress, interaction of dynamic rupture propagation with pore pressure development, patterns of rupture propagation in events nucleated naturally as a part of a sequence, the earthquake nucleation process, earthquake sequences on faults with heterogeneous frictional properties and/or normal stress, and others. The procedure is illustrated for a 2-D crustal strike-slip fault model with depth-variable properties. For lower values of the state-evolution distance of the friction law, small events appear. The nucleation phases of the small and large events are very similar, suggesting that the size of an event is determined by the conditions on the fault segments the event is propagating into rather than by the nucleation process itself. We demonstrate the importance of incorporating slow tectonic loading with elastodynamics by evaluating two simplified approaches, one with the slow tectonic loading but no wave effects and the other with all dynamic effects included but much higher loading rate.

1. Introduction
The purpose of this paper is to establish an efficient algorithm for elastodynamic shear rupture analysis of a fault governed by a general class of rate- and state-dependent friction laws in situations for which the total time of loading is vastly longer than the time for waves to traverse the domain of interest. Such an algorithm is needed to study slow tectonic loading processes during which there are episodes of spontaneous rapid failure in earthquakes. Investigating these processes requires a special approach, since quasi-static methods (used for calculating slow deformational processes of long duration) fail as instabilities develop, while standard elastodynamic algorithms not only use relatively small time steps but also require an increasing amount of memory and computational time at each time step to take into account all the prior deformation history, and hence they are excluded from direct implementation for investigating long-duration processes because of limitations on computing resources. Various solutions have been proposed. One of them [e.g., Okubo, 1989; Shibazaki and Matsu’ura, 1992] is to employ a quasi-static method during slow deformation and then to switch to a dynamic method once an instability starts. However, the abrupt switching from one scheme to another may disrupt the natural development of the instability, and the effects of this disruption on the further model response cannot be easily determined within this approach. Other approaches [e.g., Cochard and Madariaga, 1996; Myers et al., 1996] neglect all aseismic fault slippage, so that stressing between earthquakes is trivially modeled, and give the fault a ”kick” in the form of an abrupt small strength drop, once a critical stress has been reached somewhere. At that stage an elastodynamic algorithm calculates rupture until arrest occurs. This inevitably generates a population of small ruptures, and it requires careful study of dependence on the abrupt strength drop magnitude to separate which may be physical and which are artifacts of the abrupt drop [Cochard and Madariaga, 1996]. Still another alternative is to use a plate loading rate which is only a few orders of magnitude less than representative seismic slip rates, rather than the roughly 10 orders as for natural faults, and to use standard elastodynamic numerical methodology throughout (like in the work by Shaw and Rice [2000]). This is straightforward to implement, at least if some provision is made for dissipating wave energy, but makes it difficult to suitably model aseismic slip processes and can blur the distinction between aseismic slip before instability and small earthquakes. The developments of the present work provide an integrated numerical scheme allowing resolution of both slow and fast deformational phases, as well as the transition between them, within a single mathematical framework for 2

elastodynamics. The method enables us to perform calculations over thousands of years of slow tectonic loading, punctuated by earthquakes and the processes which lead to and follow them. Thus we can resolve aseismic slip on velocitystrengthening fault regions, advance of slip into more firmly locked zones, and slowly accelerating aseismic slippage that grows in spatial extent and will ultimately break out into an earthquake but has duration that is vastly longer than the seismic event itself. We also resolve all details of the break out of rupture, its propagation and arrest, and the transient post seismic slippage that develops. Our methodology for studying slow loading processes has two main ingredients. The first is based on the form of elastodynamic relations that we use, in which the dependence of the inertial response on prior deformation history can be truncated so that only a (fixed) part of the deformation history back from current time needs to be considered. That translates into fixed memory requirements and fixed amount of computation per each time step. It also makes the computation at each time step independent of how much time has already been simulated. The methodology is illustrated in this paper for the two-dimensional (2-D) antiplane case and uses a spectral representation of elastodynamic relations developed by Perrin et al. [1995] in which the slip distribution is represented as a Fourier series in the spatial coordinate, truncated at large order, and fast Fourier transform (FFT) methods are used. The corresponding methodology for the 2-D in-plane and 3-D cases is conceptually very similar and can be easily adopted from the one presented here using 3-D spectral elastodynamic relations developed by Geubelle and Rice [1995] and Cochard and Rice [1997]. Our algorithm can also be generalized to the (closely related) space-time boundary integral formulation. Furthermore, for situations such as elastic property heterogeneity that are not congenial to spectral or boundary integral approaches, finite difference or finite element procedures could be used, not in their conventional application to directly calculate the rupture propagation itself, but rather to calculate and numerically tabulate the convolution kernels for use in our methodology. The second ingredient is variable time stepping. The size of the time step to be made is dictated by the current values of slip velocities and parameters of the constitutive law. The smaller the slip velocities, the larger the time step, and vice versa. While the truncation of the convolutions over prior slip velocity history reduces the amount of computation required to complete one time step, the variable time stepping reduces enormously the number of time steps needed to simulate processes during the essentially aseismic phases of deformation which constitute almost all of the fault history. Throughout the computation, time steps can change by many orders of magnitude in value, allowing us to go in rel-

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
atively few steps through periods of essentially quasi-static loading, to consider more carefully the nucleation phase, and to resolve in great detail the features of the dynamic propagation during an instability. The coefficients of proportionality between the time steps and slip velocities depend on the parameters of the constitutive law as well as on numerical stability considerations that we derive here. We present the formulation for a general class of rate- and state-dependent friction laws with a positive direct velocity effect. The prototype of such laws is the experimentally derived logarithmic law of Dieterich [1979, 1981] and Ruina [1983]. The presence and size of the positive direct velocity effect for the quasi-static range of slip velocities, amply documented in such experiments, are shown to be crucial in allowing long time steps during slow deformation phases without losing stability (during such phases, velocity-strengthening parts of the fault zone are continuously slipping, producing an aseismic viscoplastic type response to which rate- and statedependent friction then reduces). The main goal of the present paper is to give the detailed description of the method in its current, much improved form. Earlier versions of the methodology were briefly outlined and/or implemented by Zheng et al. [1995], Rice and Ben-Zion [1996], and Ben-Zion and Rice [1997]. We describe the algorithm ingredients in sections 2-6. The new developments include understanding constraints on the time step during slow deformation phases and corresponding limitations on the procedure applicability (section 4), much more efficient truncation and evaluation of the convolution integrals involved (section 6), and a new procedure for updating the system in a time step (section 5). These developments allow consideration of a much wider range of the constitutive parameters, better numerical convergence of the results, and enhanced resolution in time and space with the same computational resources. The proposed methodology can be used to address a number of important issues, such as fault operation under low overall stress, interaction of dynamic rupture propagation with pore pressure development, patterns of rupture propagation in events nucleated naturally as a part of a sequence, the earthquake nucleation process, earthquake sequences on faults with heterogeneous frictional properties and/or normal stress, and others. Section 7 demonstrates the implementation of the algorithm by considering the elastodynamic response of a 2-D crustal strike-slip model, with depthvariable properties, descended from the model of Tse and Rice [1986] and studied by Rice and Ben-Zion [1996] and Ben-Zion and Rice [1997]. Considering a wider range of constitutive parameters than the range tractable for previous studies, we observe that small events appear for lower values of the state-evolution distance. The nucleation phases of the

3

small and large events are very similar, suggesting that the size of an event is determined by the conditions on the fault segments that the event is propagating into rather than by the nucleation process itself. We show how insufficient resolution in time can produce more complex slip accumulation that looks ”smooth” and plausible yet is just a numerical artifact. We also evaluate two simplified approaches, one with the slow tectonic loading but no wave effects (quasi-dynamic approach, as in the work by Rice [1993]), and the other with all dynamic effects included but much higher loading rate (like in the work by Shaw and Rice [2000]). The comparison shows that incorporating slow tectonic loading is very important for determining the true model response.

2. Elastodynamic Relation and Truncation of Convolution Integrals
As an illustration of the elastodynamic relations, let us consider a 2-D antiplane framework, in which the fault plane coincides with the -¡ plane of a Cartesian coordinate system and all particles move parallel to the direction. , and we define The only nonzero displacement is slip on the fault plane as the displacement discontinuity . The relevant shear stress on the fault plane is denoted by . It is possible to express the stress on the fault plane in terms of the slip history on the fault plane only [e.g., Cochard and Madariaga, 1994, Perrin et al., 1995] as

is the shear wave speed, is the slip rate, is the ”loading” stress (i.e., the stress that would act if the fault plane were constrained against any slip), and , incorporating stress transfers, is a linear functional of prior slip over the causality cone (i.e., all and satisfying ). The last term of (1) represents radiative damping [Rice, 1993], and the explicit extraction of the damping term from the functional allows for evaluation of without concern for singularities. In (1), most of the elastodynamic response is contained in the stress transfer functional . Cochard and Madariaga [1994] have expressed it as a double convolution integral in space and time. Perrin et al. [1995] have derived a spectral representation of as a single convolution integral in time for each Fourier mode, when representing slip and the functional as Fourier series in space. Generalizations to general slip and/or opening states in 3-D problems are given by Geubelle and Rice [1995] and Cochard and Rice [1997]. We use the spectral representation of Per-

 ¡ © § 7

 ¡ © !§ 5 '

U

 ¡ © § 7

U #¡

 ¡ © !§ 7

 ¡ © § 7

 I Q   ¡ § I    ¡ § G  9   ¡ SR© ! PF© !HF© !§ D

where

§ © ¡ © 3¥ 2£01)© §('  ¡ & §¥ " §¥  © ¡ © !£¤ %© ¡ © #!£¤ ©  $    ¡§  ¡§ ©  ©   ¡ D B 9A $   ¡ 7 6   ¡ § 5    ¡ § E© § ¨C@%© !§ 8(© !#3' 4© !('
is the shear modulus,

 

 © § ¥ © ¡ ¢ ¨¦¤

B

U Y a¡ $ ¡ `X U V4 § B YW $  U © #  U ¡§

 ¡ © !§ 7

 

 T¢

¡ £  ¢

(1)

4

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
Figure 1
1

rin et al. [1995] for the illustration here and write

W(p)

Figure 1. Convolution kernel lation in the 2-D antiplane case.

(4)

where and with . As any other boundary integral formulation, the spectral representation (2), (3) or (2), (4) gives as a functional of , because can be expressed in terms of the , are related to , and can be expressed in terms of by the inverse Fourier transform. We emphasize that the spectral representation, in comparison with space-time boundary integral formulations, is very advantageous from the computational point of view. The matrix of convolution integrals, implied by a spacetime formulation after discretization in space, is replaced in the spectral approach by a diagonal matrix, once the FFT

 ¡ © § 7

„

where is the Bessel function of the first kind of order one. Equations (1)-(3) are referred to as the ”displacement” representation of the elastodynamic relations. An analogous ”velocity” representation can be obtained by combining (1) and (2) with the result of integrating (3) by parts, giving

is used to transform from to and then the inverse FFT from to . Even though the spectral approach uses a larger number of degrees of freedom than needed for the domain of interest itself, the drastic reduction in the number of time convolutions significantly shortens the computation of the stress transfer functional , which is the most time-consuming stage of the analysis. We show this in Appendix B, where we further discuss the relation between spectral and space-time formulations. Note that Cochard and Rice [1997] showed how to reformulate the spectral method to rigorously eliminate the replications, but that requires far more complex calculations of the convolution kernels and still twice more degrees of freedom than needed for the domain of interest. If the convolution integrals in (3) or (4) had to be computed in full, the algorithm would be impractical for investigation of long deformational processes. Evaluation of the convolution integrals is the most computationally demanding part of the elastodynamic analysis and may take more than 99% of the total computational time [Perrin et al., 1995]. Fortunately, truncation of the convolutions is possible, which significantly reduces the overall computational time. If the duration of the physical problem is much longer than the time required for elastic waves to traverse the spatial domain of the system, it is not necessary to keep examining the influence of displacements of points on the failure surface at all prior times. This is reflected in rapid decay of

  § q

t

 ¡ © !§ 7  ¡§ © 

(3)

pm  q§

where is the length of the fault domain under consideration, replicated periodically. The replication distance has to be chosen several times larger than the domain over which rapid faulting takes place, to assure that there is negligible influence of waves arriving from the periodic replicates of the rupture process. For adaptation to our numerical procedure, (even) will be some large number of FFT sample points used to discretize this domain. Also, coefficients and are complex in general, with respective conjugates and , but take real values for and , so that the representations involve degrees of freedom. To satisfy the elastodynamic wave equation, and are related by

  § q

~

  § q

ge y k wy QR f§ xv tr q§ m z y w )s  p u kQ #SR § t m Vd © U s U ji § q G  U  B Y q Y† § ˜ — k  $ Y

q

  § q

t

“ •“ ’ ” AQ

‘

„

t

©‘ ‰ˆ SA

© U s U ji § q k  $

 ¡§ © !

 ¡ © !§ 7

  § q

‘ & r „3Rg c fc b sq ie  q †‡© ‚s€…v u § q   ¡ ƒ© !§ 7 x p ig ec hc fdb t ‰   § q “ •“ ’ ” „   § q  ¡§ © !   § q „ t A l1 § q 9 $  Y q Y† A9 q Y† –6 q q &„   § q

© ‚ € xv u  s#yw § q
(2)

ge § t „ U ge  hfd — ˜ ™ A 9 –) § q $  U  B Y q †Y § Y q Y†

& r t3Rg c fc b sq ie   ¡§ )©  p ig ec hc fdb  ¡ © !§ 7 m |  !§ t }{ t k  on § q G “ •F’ ”“ q „ t „q &t

0.8

0.6

0.4

0.2

0

-0.2 10 20 p 30 40 50

for the velocity formu-

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
the kernels of the convolution integrals for both the velocity and displacement formulations. At large values of , the kernels (shown in Figure 1) and have amplitude for an oscillation (at circular frequency decay like ) that averages in time to zero. To truncate the convolutions, we define the elastodynamic time window as the time interval , where, in the computationally most efficient version of our methodology, may be different for different Fourier modes. is chosen such that the contribution to the functional of the deformation history occurring prior to time has negligible effect on the simulation results. Since only the effects need be infrom the current time backward to cluded in the dynamic response, the convolution integrals are truncated by computing them only within the elastodynamic time window defined. This transforms (3) and (4) into

5

the framework of the space-time representation for , as we briefly discuss in Appendix B. Note that the combination of (1), (2), and (6) with (no convolution) would amount to static calculation of stress transfers, then corresponding to the ”quasi-dynamic” procedure of Rice [1993], also discussed by Ben-Zion and Rice [1995] and Rice and Ben-Zion [1996]. Because of the retention of the radiation term of inertial elastodynamics, as in (1), the quasi-dynamic procedure allows solutions to exist during instabilities; the solutions would not exist in a formulation with no damping term, which we usually call quasi-static.

3. Constitutive Laws and Space Discretization
Constitutive laws used here are rate- and state-dependent friction laws developed to incorporate experimental observations [Dieterich, 1979, 1981; Ruina, 1983]. These laws include dependence of strength on slip velocity and on an evolving state variable (or variables) which characterizes asperity contacts, thus allowing for loss of strength in rapid slip and for subsequent rehealing so that repetitive failures can occur. The laws have been successfully used to explain various aspects of stable and unstable sliding between elastic solids [Ruina, 1983; Rice and Ruina, 1983; Gu et al., 1984; Tullis and Weeks, 1986] as observed in the laboratory. Also, they have been used to model earthquake phenomena, including nucleation, ductile and brittle crustal slip regions, spatio-temporal slip complexities, and earthquake aftershocks [e.g., Tse and Rice, 1986; Stuart, 1988; Okubo, 1989; Horowitz and Ruina, 1989; Rice, 1993; Dieterich, 1992, 1994; Perrin et al., 1995; Ben-Zion and Rice, 1995, 1997; Rice and Ben-Zion, 1996; Stuart and Tullis, 1995; Tullis, 1996; Boatwright and Cocco, 1996]. A formulation of such laws which assumes constant normal stress and one state variable, to record dependence on slip history, is of the general form

(5)

and

(6)

where is the state variable and , , and depend on space variables and time. The rate- and state-dependent constitutive laws as usually formulated, based on laboratory observations, have the following properties. If slip velocity is held constant, the state variable and hence the stress evolve toward constant values, called steady-state values and denoted and , respectively, where satisfies and is given by . All laws of the class (7) reduce, for near , to , where has dimensions of slip and can be interpreted as a characteristic



™

 f '  D § f ' 



z  D § f šr w 4CQ D § ˜—#Shsk  $ ™ $ – kQ  D § !    f © D § • f ' ‹   ” f © D § Ž   !   D § fw’‘!    

 “f '

D

© #© D § …#Sa#k Ž  kQ © ¦a© D § Œ ' ‹

'



D

respectively. The discussion of the truncation implementation is given in section 6. In view of the truncation procedure, the velocity formulation has an important advantage over the displacement formulation. As pointed out by Perrin et al. [1995], the first (algebraic) term in (4) and (6) corresponds to the final static elastic stress, and the remaining integral term corresponds to wave-mediated stress transfer carrying the elastodynamic effects. The isolation of the static term is important in our computational procedure where we truncate the remaining convolution integral. During slow deformational periods is small, the static term conwhere tributes most to . During all deformation phases, the velocity formulation with truncation (1), (2), and (6), unlike the displacement formulation with truncation (1), (2), and (5), ensures that regardless of the way the convolution integral is truncated, the long-term static stress field (that is, the stress field after passage of all waves) due to slip up to the time is always exactly represented. Thus, the velocity formulation should be used for long deformational histories, although the displacement formulation can also be useful, for example, to study individual events. The same separation into static and dynamic parts, with truncation of the convolution on time within the dynamic part, may be carried out in

 ¡ © § 7

‚ ‰

9  CA § Q D B

t m ˆ‡d— A †6 © U s U ji § q G  U  B Y q † Y § k  $ 9 Y q Y†

Q § A R q

U #h U ‰i § q k  $

‚

p

 ‚ ƒ£ § $

7

 ‚ …„ § $

t

9 Y q Y† $

z €V€ w © ‚  $

p  q§ Qp

t „ U  g e wV— ˆ‡d A n†) § q 9$  U  B Y q †Y § Y q Y†

g he

  § q



t

A Š… § q 9 $  Y q Y†

  § q

i Rg „

Q £p w|

  § q G

pm  q§ t  „

B Y q Y†

‚

(7a) (7b)

6

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
at steady state given by , . The origin of the well-known result (9a) is important for our consideration of variable time stepping in section 4. We review it in Appendix A and use it to restrict the size of time steps. Such a stability dependence on the system stiffness in the case of steady-state velocity weakening has important implications for the proper space discretization, imposing an upper bound on a spatial element size in numerical modeling. To demonstrate this, let us continue with the antiplane example. Selecting equally spaced sample points along the domain of length , we discretize the domain into space , , elements (also called ”cells”) , . The discretized formulation deals with slips and shear stresses at the sample points, taken at the cell centers . Each of the cells has the effective stiffness (defined as reduction in due to elastic interactions with the surroundings for unit slip at the same sample point), given by . Here is a modeldependent constant, of order unity. For example, when using the cellular basis set for slip (i.e., calculating the as if the slip were locally uniform in each cell) like in the work by Rice [1993] and for the spectral basis set of (2). The values of cited apply for cells whose distance from any free surface is many times and, in the spectral case, for . If we start with steady quasi-static sliding of the whole (taking this rate to be sufficiently domain at slip rate small so that the dynamic effects are negligible) and slightly perturb the motion of one cell while maintaining steady sliding of the other cells, we get that the linearized response to the perturbation is governed by the same system of equations as for the spring-slider model, with the spring stiffness replaced by the effective stiffness of the cell. Hence the perturbation will grow if or, equivalently, , and the perturbation will decay if or, equivalently, . This property defines the critical cell size . As emphasized by Rice [1993], the growth of the perturbation on one cell, while the others continue the steady sliding, would imply that the cell is capable of failing independently of the surrounding cells, which would make the results dependent on the numerical discretization. Hence, to insure that the perturbation on a single cell decays, so that each cell can fail only as a part of larger space segment, the mesh should be refined enough for space element size to be much smaller than the critical cell size . In other words, the condition

(8)

where the precise definition of

is

(9b)

z D SR D § f ' k D w Crà kQ ÅÄ  9 ¢† · $ † •¡ CQ ¹ Ÿ µ 9

Here and in the following, notation means that the expression in the brackets has to be evaluated

™

·

©#|œW—WƒÂ1ÁCQ Ÿ µ Àµ Ÿµ

(9a)

(10)

ˆQA

 ¹·

9 ¢¡ † •sCQ 9 • Ÿ µ · ¢¡ † •sCQ ·  µ

¢ •¡ †

v

·

 CQ 9 @l† µ ·

µ

µ

µ CQ

º wQ ˆ j· 

9

· ¸¶†

Ÿ  D

·

9 ¢ ¡ †9 W µ ·  •s¿¾sQ ½f† ¢¡ † · W •sCQ }¼µ

|

»” €“ •“ ’

v'

This property is sometimes called ”direct velocity dependence” or ”direct effect” and is well established experimentally. As discussed in section 4, the presence and size of this direct effect are essential for our numerical procedure to be efficient in simulating processes of long duration. Hence our method is not immediately applicable to other constitutive laws, such as slip-weakening laws (which emerge as the limit case here for rapid slip if and are, in the limit, independent of ), which do not have that property. Stability of steady frictional sliding, governed by the constitutive laws of type (7) with the properties discussed above, has been extensively investigated [Ruina, 1983; Rice and Ruina, 1983; Dieterich, 1992; Gu et al., 1984; Ranjith and Rice, 1999], particularly for single degree of freedom elastic systems. Such systems are generically represented by a spring-slider model, in which a rigid block is attached to a spring of stiffness and slides on a frictional surface, with the other end of the spring moving at the imposed rate . Linear stability analysis of such a system, perturbed about steady-state sliding at the rate , as in the work by Ruina [1983], shows that the sliding is always stable for friction with steady-state velocity strengthening, while for friction with steady- state velocity weakening, there exists a critical value of the spring stiffness such that perturbations from steady-state sliding grow in time for systems with and decay in time for systems with . For rates sufficiently small so that the inertia effects can be ignored, the critical stiffness is given by

 ¶³

g & µ³ h´ v ¡ z v ¡ © v w¡

v'

† ¡  Q¦µ8$ v !§ v' A v “ •“ ’ Q ƒµ “ •“ ’ © ­yž © s| ” ” žž © A ‘

D

“ •“ ‘ ’ ”

slip distance required for evolution to the steady state ( is also sometimes denoted by or ). The state variable is usually chosen in such a way that , particularly if is to be interpreted as a characteristic lifetime of the asperity population on the contact surfaces; is then interpreted as a measure of the sliding distance required to establish a new population of asperity contacts, and is assumed to be independent of . The law (7) is said to exhibit steady-state velocity weakening if and steady-state velocity strengthening if . If the slip velocity is suddenly increased or decreased, the stress simultaneously increases or decreases; that is, instantaneous positive viscosity is incorporated in (7) through the requirement

Ÿ    Ÿ  XD § fw²‘ XD 

° ¨±§ ® ¬­¬ « r ©« ¨i§ r § z w ž R¨ § ¯­« r « ª¨ § r § z #C#© D § ¦I § Q D $ w V™   IQ  Ž  ° ®¬¬ ©

Ÿ  D ¢¡ •C†  †

Ÿ  D



™

D

¢¡ † W •s£8†

©

D Q ™   €œfw t › ' D )Ž Q ™ ž ¨W

™ ¨%§ r §3¦ ™ $¤  ¢¡ D C D § f ' k D •¥ƒ•s† kQ 

Ÿ XD

W D CQ f ' k k  oD SQ ! ' k k

‹

DI  #© D § £I ‹

¢¡ •s†

› hk

†

D

D



LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
should hold, where the maximum is sought over all allowable slip rates . The critical cell size is directly related to the nucleation size of model earthquakes, that is, to the size of the patch that initiates unstable, dynamic slip. In simulations with the constitutive laws of the type discussed here, rapid, dynamic break-out of an instability is always preceded by quasi-static slipping of a small zone. The size of this zone just before the dynamic instability is what we call the nucleation size. Changing , and consequently , changes the nucleation size in an essentially linearly proportional manis not only a ner in all simulations we have done. Thus very important numerical parameter, it is also a crucial physical parameter. Note that expression (10) is derived neglecting inertial effects. Rice and Ruina [1983] have shown (by considering the spring-slider model with mass) that inclusion of the inertial in (9a) and hence decreases in (10), effects increases which requires making cell size even smaller. They have also analyzed uniform slip (at constant slip velocity) between elastic continua, with spatial perturbation of the type and found that (1) there exists a critical wavelength such that for smaller wavelengths the perturbation is stable and for larger ones it is unstable and (2) the value of the critical wavelength decreases appreciably with increase in the slipping velocity. As confirmed by simulations, proper resolution of high slip velocities during dynamic instabilities requires the statically estimated critical cell size to be discretized by tens and sometimes (e.g., in the case of strong velocity weakening) even hundreds of cells . Proper discretization is further discussed in sections 4 and 7 and by Zheng and Rice [1998]. The need for such fine discretization stems from the rateand state-dependent friction laws that we discuss here. As already mentioned, these laws are supported by experimental evidence at low , and their feature of state evolution over a slip distance is supported by the concept of a characteristic slip required for renewal of the asperity contact population. The laws produce high slip velocities near the rupture tips and incorporate small characteristic slip distances to be resolved there, and hence require fine discretization in space and time. The discretization constraints may be possible to relax by using modified forms of the friction laws, depends on for example, in which in the law for and increases significantly for the seismic range of . Taking proportional to at high slip rates would be equivalent to having state evolve over a characteristic time (rather than over a characteristic slip distance), as for the velocityweakening range of the ad hoc type of friction law used by Shaw and Rice [2000]. Such modifications and their influence on the qualitative features of the simulation results still have to be explored.

7

are the characteristic slip distance, where , , and the current slip velocity, and a prescribed parameter for the ith cell of the discretized domain (introduced earlier by , , , ), respectively, and the minimum is sought over all the cells. The choice of parameters depends on the constitutive law and stability considerations as explained below. Criterion (11) allows us to adjust the evolution time stepping during a simulation based on current slip velocities of the cells, so that slip in a time step does not exceed a fraction of the characteristic slip distance of the friction law, the fraction being prescribed by for cell . The adaptive time step can be enormously longer than the time for waves to propagate over the space domain during periods of slow, essentially quasi-static, loading, before unstable rupture begins, and can be very small during spontaneous failure, spanning up to 10 orders of magnitude in value in some of the simulations that we have done. In between occurrences of dynamic ruptures, when slip rates are very small, we would like evolution time stepping to be as large as possible without compromising the algorithm accuracy and stability. A rather insightful constraint on the time steps at low slip rates can be derived by considering, as in our motivation for existence of the critical cell size, quasi-static stability of perturbed motion of a single cell with continuing steady sliding of the other cells at velocity . If the grid is properly refined, then the perturbation on a single cell dies away, as we have already considered. Demanding that our time discretization preserves this property, we get a condition for the size of the time step allowed, which provides a constraint for , from (11). We derive this constraint in Appendix A by analyzing, as a simple

Ê Ë“  É

ŸD

Ê Ë“  É

v y g “ •“ ’ Q ̵ “ •“ ’ © ­yž © #`„³ h² v ¡ z v ¡ © & v w¡ ” ” žž © |  µ³ A ‘

Ê © z v D Q v ™ v wy ­Ç 1 ˓  É ÈÃ

“ •Í’ © ­yž © s†%³ v y ” “ žž © |  A

³

vy

vy

vD v™

™

D

ٵ

D

ٵ

µ

ٵ

ٵ

kQ #Shsk

µ

ٵ

™

™

D

¢ •¡ †

D

™

D

4. Variable Evolution Time Step
Simulating truly slow loading while capturing details of occasional rapid failures requires varying evolution time steps. Our time step selection criterion is based on two observations. First, we recognize that the slower the particle velocities in a rupture process are, the longer the time steps should become, and vice versa. Second, to assure proper integration of the constitutive law during the calculation, we would like the relative displacement in each time step to be small compared to the characteristic slip evolution distance . To fulfill both of the above requirements, the time step from one updating of field variables (slip velocity, stress, etc.) to another, which we call the evolution time step , is chosen as (11)

¢ •¡ ‘  ‘   SA § ­Æ Q ˆ ÈÇ ™

8

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
in purely quasi-static phases of the analysis. The condition also reveals that if the direct effect is decreased, then the time steps should be chosen smaller as well. That is why the efficiency of our algorithm, which relies on using long adaptive time steps during quasi-static loading periods, depends on the size of the positive direct effect (quantified by ) for the quasi-static range of slipping velocities. Note that (12) or (14) are not applicable to the case , as their derivation (Appendix A) assumes nonzero . Moreover, for a certain range of (very small) values of , the inertial effects become comparable to the direct effect even for small sliding velocities and can no longer be ignored. The Rice and Ruina [1983] inertial analysis indicates that linearized perturbation to steady-state sliding of all wavelengths are unstable in the case of . In practice, it is possible to simulate a single dynamic event in the case with , apparently without numerical instability (A. Cochard, private communication, 1999), possibly due to the low rates of growth of the instability for the highest modes. However, the time stepping has to be so small that long deformation histories are excluded from consideration. Since our algorithm relies on using long adaptive time steps during quasi-static loading periods, it can be efficiently used only for the constitutive laws that exhibit the experimentally verified positive direct effect as in (8). On the basis of conditions (12), we choose parameters , (used in selection of the evolution time steps (11)) as

(12a)

if

and

(13)

if

(14b)

µ CQ

9

·

v™

Ú ƒ• É ÙØ

Conditions (12) or (14) give an estimate of the required time stepping for low slip velocities; if these conditions are not met, cell-by-cell instabilities arise which either make the simulations impossible or corrupt the results. Deriving these restrictions is an important development in the methodology, as they explained and eliminated many of the numerical difficulties that we had. Note that when the condition (12a) or (14a) is applicable (which is often the case since large is required for proper space discretization), it implies that if we refine the grid (taking smaller cell size and hence larger in (14a) or larger in (12a)), then we have to decrease the time stepping as well, even

and subscript denotes the value of the corresponding quantities for the cell and is the single-cell stiffness, . The term 1/2 enters (15) to enforce the condition that for each cell the slip in every time step is not larger than half of the characteristic slip distance . The time step selection criterion (11) and (15) captures the essence of our variable time-stepping scheme, but in actual simulations we modify the criterion slightly to reconcile the variability in time steps with the necessity to compute convolution integrals, uniformly discretized in time. To store deformation histories in an efficient way, we introduce a time parameter, , which is the minimum value of the

ž

Ó

(14c)

Î Ÿv v ¨† $ ™

i

Ÿv

† ³ ³ Î Î Î Ÿv Ÿv Ñ º $v vÐ v $ŸÏ ™ † |

vÐ

†

žµ $ Ÿµ

i

Ó

Î µ Ñ ÎŸ º | ¹$ Ÿ µ Ÿ $ŸÏ|

µ CQ Ÿ µ

µ CQ Ÿ µ

Ð

µ

Ð

if

, with

if

, where (15c)

× A©v Ÿ

|

Îv ™¨† $ ÖÈ Ã „| ryÇ … v y v $ŸÏ

WvÐ

if

and

and

Î

v $ Ÿ Ï!Î § $ v ™† ÖÈ Ã r­Ç … v y

(14a)

Ÿv × €A ©  |

Ÿv

“ ”•“ ’ © ­yž © #Õ4³ v y žž © |  A

As the derivation in Appendix A shows, constraints (12) are applicable to both steady-state velocity strengthening and steady-state velocity weakening with a sufficiently dense grid. For the latter case we have and , and hence (12) can be rewritten in a more insightful form

Ÿ

Î

Î

¢¡ †  ™ Q •CŠ8C Ÿ

© ° ¨i§ ® ¬y¬ « r ©« ¨%§ r § z D ICQ#© D §(' I D w  Ÿ  ŸÏ

Î

ž z D k Q D §  ' wD ¨ § r § C f k

Ï $ Ÿ !§

| $ µ Ô8nCQ Ÿ µ §  Î Ÿ Ÿ

ӟµ

µ

$ „|

Î

Ñ ÒŸ D

Ï $ Ÿ § Ÿ D  É  ™

™

| W µ µ  ¢¡ † Q l‘sQ Ÿ œƒ•Cs3†

É

ŸÏ $Ÿ

Î

Î

Ÿ

and

and

are given by

Ÿ

Ÿ

™ ¨†

Î

$

i

Ó

Ÿ

ΠΠΟ Ÿ Ѻ $ $ŸÏ ™ † |

Ð

Ð

if

, where (12c)

Î ŸÎ

ŸÎ

(12b)

(15b)

Ÿ

Î

Ÿ

Ÿ

model case, explicit integration of the governing equations with a constant time step . From that we obtain:

Î

Ÿ Ó Ÿ Î

Ñ ÒŸ D Î ™ † É $ „| $ŸÏ ™

Ï $ Ÿ Χ j† Ÿ D  É $™  Ÿ ™

É

WÐ WÐ Î

(15a)

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
evolution time step allowed (and also the discretization interval for computing the convolution integrals, as explained in section 6). It is selected as a fraction of the time needed for elastic waves to traverse a spatial element, in the form

9

to the desired degree, whereas sufficient increase in does, provided other parameters, such as elastodynamic time windows, replication period, etc., are chosen appropriately.

(16)

We insist that every time step we take be an integer multiple and not smaller than . That is, we first comof pute the (tentative) time step using criterion (11) and (15) and then convert it into a multiple of : int

5. Updating Scheme: Advancing One Evolution Time Step
Let us consider how the values of field variables are updated over one evolution time step. We will use the velocity formulation without truncation in this section, for generality, and consider truncation of the convolution integrals in section 6. We continue the antiplane case with the domain discretized into cells , , , . Suppose that the discretized values of slip , slip velocity , state variable , stress , and state rate at cell centers are known at time for all and that the slip velocity history is known for all prior time , , where is the beginning of the deformation process considered. When using the spectral formulation, we also assume that the Fourier coefficients of the slip distribution are known at time and note that the velocity history need only be available in the for . To Fourier domain, as the values of advance the field values by one evolution time step and to determine all the quantities just mentioned at the end of that step, we proceed in the spirit of a second-order Runge-Kutta procedure as follows: 1. Determine the evolution time step to be made using criterion (11) and (15) - (17). 2. Make first predictions of the values of slip and state , based on known values at , as variable at time

(18)

É

  § v D

| ê —$ @X³ &u “ •F’ Q ”“ åq è e c b å v é – å c f±Ëg q i µæ hç ª¡

gåè q& u c fc±“ Ëg Íq ’ v ­Ryižéž ž ˜ e b ”•“ å © © A #fTæ ©|  v   É ƒ § Ÿ  6

  § v D

If a better resolution in time is desired, it is often advantageous, for better stability and faster convergence of the results, to keep or and to increase (and hence ), rather than to keep and to decrease , even though increasing is more costly in terms of computational time and memory. If there are numerical oscillations or other features in the simulation that point to an inadequate resolution, possibly in time, decreasing without increasing does not always solve the numerical problems

3. Make a corresponding first prediction of the functional, using slip prediction (19) and treating the slip rates as if they were constant through the time step and equal to . To implement this in the spectral formulation, we first compute the Fourier coefficients of and . To represent FFT operations, we shift the coordinate origin so that the cell centers are at , , and define and (unless used as sub, where

v   É ) § Ÿ 7 6



  U 

Ê Ë“ Ü  É



 r É

v © ä § v D  É d § v  ) É d § Ÿ  6  6

ž  § v G   É ± § v 1) É d § Ÿ v   6   6

U§ q G

t

 É ( 6

  § q

t

The parameter determines how fine our resolution in time is. To understand how to choose (or from (16)), let us consider how the evolution time step during dynamic instability relates to . We recognize from (1) that a characteristic slip velocity of order is induced by an abrupt dynamic stress drop . Hence, to resolve the characteristic slip distance of the friction law, has to be of the order . At the same time, using our constraints (10) on the grid spacing , we can express , where . From the above formulae, is comparable to . The conhas to be tens or hundreds stant is of order unity, in order to properly discretize the critical cell size , and the ratio can be considerably smaller than unity. This suggests that the smallest required time step is comparable to . We have found that if other parameters, most notably , are chosen appropriately, the standard ”sampling” choice of , giving , produces stable and satisfactory results for the cases we considered; can also be successfully used in most cases. Note that (16) can be rewritten as

(19)

‘

(17)

U  —   U  “ •Í’ © ­yž © sâã³ ” “ žž © |  A    § v G   § v '   § v    § v D   § v  “ •“ ’ Q ” g  ‘ âµ “ •“ ’ © ­yž © s@ͳ h` v ¡ z v ¡ © & v w¡ ” žž © |  µ³ A

Â

” •“ ¡  É ” Ê Ë“  É Ú ƒØ Û Ù ' XÉ µ

Ú ƒØ Û Ù µ µ CQ Ÿ †‘Â

” •“ ¡  É ” Ê Ë“ É ' XÉ Q ÔXÉ á' Ÿ µ µ µ sQ Ÿ œn · ' ÍÉ ” •“ • É Q ˓  É ”¡ Ê · á' Q ÔXÉ Â zD á' kQ SR D § f ' k D w Cr“` ÍÉ ÅÄ Ã $  9 á' ”¡  ÔXÉ B  § €™ –à B  § Q Ÿ ¶ B ¦¶ ” •“ • É Q · µ Qµ

” •“ ¡ É ” A Q •

™

9 '  XÉ B § Q ™

Ú ƒ£Û ÙØ

” •“ É ” Qµ B ¦† ¡ 

Ú ƒØ  É Ù

ž   § Q Ÿ sÚ ƒ3… ¦sÚ ƒ£ŒÍÚ ƒ• É µ ÙØÛ B BQµ Ù ØÛ  Ù Ø

Ú ƒ• É ÙØ

ž hÚ ƒ• #˓ Ú ƒ• É ÞCr1 ˓ ÒÉ Ê Ü ß Ù Ø ÉÊ ‰© Ù Ø Ý ÅÄ Ã

ž ¦sÚ ƒ3… ” •“ • É Ú ƒ£ŒÍÚ ƒ• É ”¡ ÙØÛ  ÙØ BQµ Ù ØÛ  ÙØ ‰Ú ƒ• É Ê Ù © z Ú ƒØ  É Q ˓  É w  ºQ wS| Ê Ë“  É Ú ƒ• É ÙØ Q A w|  “ •“ ’ ” Ù Q| A SÕ Ú ƒØ Û Ú ƒ• É ÙØ ºQ|  Ù Ø Sw¾ÒÚ ƒ3Û Ê  ˓ ‰ Ú ƒ3Û ÙØ Ê Ë“  É Â µµ CQ Ÿ HP 9 XÉ Q' B Ú ƒ• É ÙØ

10 script), to get

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
(The superscript double asterisks could be dispensed with at this point, but is useful for comparison to other methods.) 6. Make a corresponding prediction of the functional, using the and treating the slip rates as if they were constant throughout the time step at , consistently with updating slip in stage 5. The specific steps are analogous to (20) in stage 3. Note that the second term in the brackets of (20b) will be the same in this stage, and hence it can be computed just once, in stage 3, and stored for use here. This is computationally very advantageous, since this term incorporates most of the convolution evaluation. 7. Make final predictions and of the slip rate and state rate, similar to stage 4, using the new predictions for the state variable from (22) and of the last stage for the functional. the result That is, solve the equation like (21), but with superscripts double asterisks, to find , and then compute from (7b). 8. Declare the values of field quantities , , , to be equal to the predictions with the superscript double asterisks. Compute the corresponding , if needed, from (7a). Store slip values of stress velocity history for the time interval , to be used in future convolution evaluations, as for in . In the spectral formulation, store instead the history of the Fourier coefficients as for in , since they are actually used in the convolutions, and set . Finally, return to stage 1 to advance through the next time step. This scheme is second-order accurate in for the slip and state variable, assuming that the predictions of the functional are computed accurately enough (N. Lapusta, Ph.D. thesis in preparation, 2000). We use the midpoint integration scheme to compute the convolution integrals. Note that while it is important to know the order of accuracy of an updating scheme, the actual performance also depends on other things, such as stability characteristics, balance of terms that achieves most error cancellation, etc. Ultimately, the most important thing is the ability of a scheme to produce numerically stable simulations with results convergent through space grid reduction and better time resolution, which can often be checked only by actually doing the simulation. Earlier studies used different updating schemes that provided a foundation for the development of the scheme described here. The scheme used by Rice and Ben-Zion [1996] and Ben-Zion and Rice [1997] proceeds through stages 1 to 4 and declares the values of the field variables at time to be equal to the predictions , , , and , all of which are first-order ac-

Then, using (4), we get the Fourier coefficients of the prediction of the functional

(20b)

(21)

and then solving (21) for . We find using Newton-Rhapson search with as the first guess. Once are obtained, the corresponding state rates can be readily found from (7b) as . 5. Calculate the final prediction of slip and state variable at time by

(22)

v v   É C § Ÿ   É S § Ÿ  6  6

É

v   É ( § Ÿ G  6

v   É d § XD 6Ÿ   É S 6

t

4. Find predicted slip rates corresponding to the predicted state from (19) and functional from the last stage. This is done by equating stress (1) to the strength (7a) to get

Q   É „ § Ÿ q G 6 A 

z  É ¨E w 6 ©

 É i § Ÿ q 6Ÿ t 6 „ § q G §

U

(20c)

6 # § v   É # § v   6

t ) É  § q  6 U t z  É €E w 6 © n U  § q G v 6 d § ŸXD 6d § v wD  U  § v D   z  É aE w 6 ©

  É % § v ' 6

  É E § v G   É  § v D  É 6  6 

Within the brackets of (20b), the second term can be computed since the slip velocity history is known. The third term is an approximation of the convolution on the time interval corresponding to the current time step. We then obtain the prediction of the functional through an inverse FFT as

v   É C § Ÿ G  6Ÿ

v   É C § XD 6 ŸŸ

v   É „ § ÍD 6 ŸŸ

v Q   É i §  D i § v D § 6Ÿ 6 A R v   É 4 § vRŸ 7 6Ÿ   É d § Ÿ  6Ÿ t

(20a)

  É

v 6Ÿ F § Ÿ 7

v   É 4 § RŸ  6Ÿ

v   É ( § Ÿ G  6Ÿ

z É A Q 

v v   É d § Ÿ E É d § XD § Ž 6 © 6Ÿ v  ¶ É Í § Ÿ G  6 v   É F §  D 6Ÿ v   § Ÿ D v v  É 6d § ŸÍD   É ± § XD 6Ÿ v v ò  É d § Ÿ E É d § Ÿ D  Œ 6 © 6 ñ‹

ž ¦ U s U  Y q Y† § m ˜ d — kB ë t § m˜ U #h U j$± É ( § q G  U  k 6 B Y q Y† t   É d § Ÿ q ¤ 6

v © z  É d § Ÿ D d § v wD  6 6

ž Ëz  É d § Ÿ v G  ± § v G  w A ± § v … É d § Ÿ v  6 6 6   6Ÿ É

vB v v   É ( § Ÿ D C9 A i É d § Ÿ 8d É d § 3' 6 $ 6 7 6 65

& r Cq „ ð ï ìí ì g v ˯î¯ib v ž   É d § Ÿ q q & è 6  ) É d § Ÿ 7 6 p ig ec Rhc f±b

t t t ž  § q G  É 6± § q  4 É ( § Ÿ q 6 g rv t §q p 4 G c f±b ec v  É  § Ÿ  6 v    É € § XD 6Ÿ vv © ä § D q è t  § ˜ q G $ ë ˜ š— $ ë ("

A d § v   É d § Ÿ v  6  6Ÿ É

A 9 1 É ( § Ÿ q $  6 Y q Y†  É  6 „

v   É 4 § Ÿ 7 6

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS

11

If is the same for all , then for the highest mode is times longer than for the lowest mode. Since the convolution kernel decays (Figure 1) and is usually a large number (spanning values from 512 to 65536 in the simulations we have done), much of the computation for the highest modes has negligible contribution. To save computational time and memory, we examined use of time windows which are mode-dependent and significantly shorter for the higher modes. In the current implementation, two parameters determine the window sizes. One of them is , the time window for the lowest mode , and the other is , the ratio of the for the highest and the lowest modes. Once these parameters are selected, the for all the modes are determined as

‘

ú 

‚ | § 

Let us consider the evaluation of the truncated convolution integrals in the velocity formulation (1), (2), and (6). The elastodynamic window (introduced in section 2) can be selected the same for all Fourier modes, or it can be mode-dependent. We have studied both approaches. Let

‚ù

 w| § ‚ 

6. Evaluation of the Truncated Convolution Integrals

That is, the kernel window for the highest mode is times longer than the window for the lowest mode, and the kernel windows for the modes in between vary linearly with . The corresponding time windows can be found from (23). Note that for , this approach is equivalent to the one with the constant . The advantage arises from the fact that far smaller can be used; acceptable values can be as low as 4 for some problems, including our implementation examples in section 7. We determine the parameters and by trial and error, starting with an educated guess and then comparing the results with the ones for smaller and larger values of both parameters, until convergence is reached. A useful parameter for making the initial guess is the time for elastic waves to propagate through the domain of size treated

‚

Q Bˆ ‘ a Y ‰ Y (SA

“ •F’ ”“ | ž Ó ¹wY Y § ¹$ A Q | $ ‰ |$ ¹‰‚ ù | ¹$ | $ ¹wY ‰ Y § § ÷ A  w| ‚ $ %

”“ ’ © ä| § ‚ ÷ ‚ ù ƒ A Q “ •F§ ‚ ÷ 

‚÷

‚ù ‚  “ •F’ ™‚ ù ”“  AQ

“ •F’ ”“

 ` B Y q Y†

”“ ’  A Q “ •Í!§ ‚÷  ‰ § ¦ ‚ ‰ ž  § ‚• Q §Y § ÷ ‰  ‘ (SA Y ‰ ²) ‰ ‚ Bˆ  ‰§ ‚ ÷ ”“ Q “ •F’ sY  wA zQ “ •F!§  “ •Í’ ”“ ’ ‚ ”“ | £Y B ˆz © ‘  A Q Q § ‚ B ˆ ‘ w| ¦ (SA ©

‚ù  | § ‚ ÷ ”’  “ •“ § ‚ ÷ AQ Ñ 6 4| | § ‚ h ‘ Q (SA §    Bˆ ” Q “ •“ ’ 6  § ÷  ‰ § ‚ ÷  “ •F§ ‚ ÷ (| ‚ ”“ ’ AQ

‰

‚ù

©E| § ‚¦h ‘ Q B(SA § )w| § ‚ ÷    ˆ ‚÷

‰

‚

Y ‰Y

‚ w| § ¦

| øàY ‰ Y

‰Y ‰ Yw

curate (N. Lapusta, Ph.D. thesis in preparation, 2000). The only departure in this ”incomplete” scheme is in stage 2, where the value of the state was computed using not (19), but through exact integration of (7b) assuming that the slip velocity was constant in time throughout the step. Another updating scheme was originally developed by Morrissey and Geubelle [1997] for the case of constant evolution time steps and constitutive laws without state variables and described by them as a ”semi-implicit velocity formulation”, ”with delay”, ”discretized kernel”, and ”convolutions by trapezoid rule”. It incorporates steps similar to stages 1-5, and then uses predictions , , , and as the values of the field variables at time . It also approaches differently the convolution evaluation, using a trapezoidal rule and delay in kernel. The delay in kernel, discussed in detail by Morrissey and Geubelle [1997], was introduced as an empirical step by Cochard and Madariaga [1994] to smooth numerical oscillations in slip velocity right behind the rupture front, but at the cost of reducing the slip velocities at the tips of the propagating disturbances. The present scheme performs better than both of the above mentioned schemes in the cases that we considered. It does not use the delay, captures more accurately the (high) slip velocities at the rupture tips, and has essentially the same stability performance. However, it is more costly in terms of the computational time (but not memory, which is often the primary limitation), mostly because it uses another pair of FFT transforms at the stage 6. Each of the FFT transforms requires floating point operations. The other time-consuming computation is the evaluation of convolution integrals, which, as considered in section 6, requires to operations, depending from on the truncation procedure used. Clearly, if the truncation scheme used requires operations, then the extra FFTs do not make much difference in terms of the cpu time. However, if the more efficient truncation procedure can be used in the problem at hand, then the number of operations for the FFTs and for the convolutions can have comparable orders of magnitude, in which case the advantage of the full scheme 1-8 has to be weighted against the increase in the computational time.

denote the length of the elastodynamic time window for mode . Keeping the window the same for all Fourier modes simplifies the procedure, but it is not a very efficient choice, for the following reason. The ar, gument of the convolution kernel, depends on the mode number and varies in the ranges for the lowest (spatially nonuniform) mode and for the highest mode . We call the length of these ranges kernel windows , so that (23)

Y ‰ Y § ‚ ’Ò ‰ § ‚  

  v É s § 6  É  § Ÿ G  6 v v v    É à §  D  É Þ § Ÿ   É Þ § Ÿ  6Ÿ  6Ÿ  6Ÿ “ •“ ’ § ”  i ró “ •“ ’ § ” ”“ ’ ö õô ”“  i !ró ò  “ •Í!§ i sy¨“ •F’ ñ ó ‚  ” “ ’ ö õô ” “ ò  “ •Í!§ i #­¨“ •Í’ ñ ó

(24)

12

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
is a (25) overall cpu time due to the mode-dependent is very significant. The analogous reduction arises in memory requirements, as much less deformation history has to be stored for higher modes.

in the spectral formulation. In general, good initial guess. We select in the form

y

z

with depth only, δ = δ (z,t)

24 km

-24 km < z < 0: Fault zone with depth-variable properties; rate- and state-dependent friction law applies

-96 km < z < -24 km: Moving substrate; slip imposed at uniform rate of 35 mm/year

Figure 2. A vertical strike-slip fault in an elastic half-space (like in the work by Rice [1993]).

¦

¥¥¥ ¥¥¥ ¥¥¥ ¤¥ ¤¥ ¤¥ ¤¤¤ ¤¤¤

¥¥¥ ¥¥¥ ¥¥¥ ¤¥ ¤¥ ¤¥ ¤¤¤ ¤¤¤

Slip constrained to vary

x

 

 ¡ © § 5 '

ÿþ ý ü º A  䕓 ƒû

 “¢

”  ªý D $ D §

 ¡§ © 

¥¥¥ ¥¥¥ ¥¥¥ ¤¥ ¤¥ ¤¥ ¤¤¤ ¤¤¤

¡  !§ 5 '

§





¥¥¥ ¥¥¥ ¥¥¥ ¤¥ ¤¥ ¤¥ ¤¤¤ © ¤¤¤

”  Ôý D €© !!§ $  ¡§ 

¨

¢  £¡

¡

” Ôý D

where, from our experience, of 1 to 4 are sufficient values for most problems. In order to use standard procedures for computing the using the (conconvolution integrals, we discretize stant) time interval introduced in (16). Computing the field values using evolution time steps from (17), , we store the field values which are multiples of needed for convolution evaluation on a uniform time grid of spacing , repeating each value times. The lengths of elastodynamic time windows determine how many values have to be stored from the current value back in time for each Fourier mode. This array of the stored deformation history contains the values needed for the discretized convolution in the useful format (spaced ), which simplifies and speeds up the computation. by At each evolution time step, the array of stored history is updated by writing newly computed values over the values that have moved outside the elastodynamic time windows . Hence only a fraction of the array is typically updated at each time step, and the assignment for each value in the array is done only once. In earlier implementations [Zheng et al., 1995; Rice and Ben-Zion, 1996; Ben-Zion and Rice, 1997], the field values were stored at (variable) evolution time steps, and then time-consuming search and mapping routines were employed in every time step to transfer values from the set of the evolution time steps to the uniformly spaced array required for the calculation of the convolution integrals. Obtaining this array was the most expensive part of the computation, being up to 10 times more time-consuming than the multiplications needed to evaluate convolutions. The current procedure reduces the cpu time for updating the array of the stored history by orders of magnitude, so that the cpu time for computing convolution integrals is essentially determined by the cpu time required to perform multiplications. Let us estimate the order of magnitude of this latter cpu time. When the time windows are the same for each Fourier mode, they are taken to be of the order of the time for elastic waves to propagate through the domain of interest. In this case, the convolution evaluation, at each time step, requires floating point operations for each Fourier mode and operations altogether (for all modes). In the case of dependent on the Fourier modes according to (23)-(24), an upper bound on the number of operations for the mode scales as , and the number of operations for all modes is at most . Since the typical values for are thousands to tens of thousands, and can be as small as 4, the reduction in the

7. Implementation Example
7.1. Formulation of 2-D Model and its Response To demonstrate how the ideas outlined in the previous sections are combined to produce long-duration simulations, let us consider elastodynamic response of a 2-D depthvariable fault model (Figure 2), to which earlier implementations of related procedures have already been applied [Rice and Ben-Zion, 1996; Ben-Zion and Rice, 1997]. In this model [Rice, 1993], a vertical strike-slip fault with depthvariable properties is embedded in an elastic half-space. The fault is driven below depth km with a plate rate of mm/yr. In the shallower zone, governed by a constitutive law, the slip is calculated as a function of depth and time (variations with along-strike distance are not included in this 2-D model). This model can be mathematically described as the antiplane problem discussed in previous sections. It proves convenient to express the formulae in terms of variables and , in which case , the stress which would act if the plane were constrained against any slip, becomes independent of time and equal to the initial stress . Hence, following relations (1) and

‚

ú (

 ‰§ ‚ 

“ •Í’ ”“ ò ”“ ’ ô ”“  “ •Í!§ È 3“ •F’ ‚ ù ñ ó ù§ Y ‰ ¨Q “ •F’ ‚ ró Y ”“

ú ¼Òw| § ‚    Ê˓ Ü  É

 ‰ §  ‚ Ê Ú ƒØ  É Q ˓ Ü  É Ù

 ‰ §  ‚

© B Q ‘ ‚ … ú  ‚ 1)| § ‚     w| § ‚  ‚  Ú ƒØ  É Ù ‚Â Ú ƒØ  É Ù

‚ “ •“ ’ § ”  i !ró ” ’§  “ •“ !ró

‚ù

Ú ƒØ  É Ù ‰

Ú ƒ• É ÙØ

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
(2), we write

13

(26)

(27)

(28a)

except that we regularize the first of these near as discussed below and, to allow of either sign, replace with in (28b). In (28), is the effective normal stress, is the value of friction coefficient at the reference velocity m/s, and are frictional parameters, and , as before, is the characteristic slip distance. As (28) indicates, , , , and vary with depth but not with time. Examples of dynamic modeling with addition of



D

(28b)

A drawback of the logarithmic form (28a) is that the stress is not defined for . The logarithmic form was derived from purely empirical considerations to match experimental observations [Dieterich, 1979, 1981; Ruina, 1983]. However, it has a theoretical basis, in that such a form would result if the direct velocity effect is due to stress biasing of the activation energy in an Arrhenius rate process at contact junctions, at least in the range for which forward microscopic jumps, in the direction of shear stress, are overwhelmingly more frequent than backward jumps. Such interpretation seems implicit in the work by Chester

·

º SQ ˆ j· 

“ •Í’ Q ”“ A

00 ( 

§ ž h $€9 0 º

™ $ SQRh1” § ( •C† 0  ¢¡ Ÿ Ï )( Ÿ 0

µ CQ

 6

9

™

·

µ CQ

“ •Í’ Q XSQ ”“ A 6º ˆ

2

ˆ

ž º 3 j·

9·

· ¶‰†

ٵ

Î

4ž 532¨

“ •F’ ”“

 Š·

·

–º ¹wQ ˆ

º SQ ˆ

ٵ

 –·

with the relation between and given by (6). We use GPa and km/s. We select the replication distance in the direction in the following way. In addition to the region where we wish to simulate slip, we include in the substrate (mantle) region , where a plate velocity of mm/yr is imposed. That means there, so FFT sample points in that part of the domain create zero padding. Finally, to model the free surface at , we map the slip and slip velocity from the region to the region as even functions. Hence, the spatial domain to consider is in the direction, and or 192 km for km. As noted before, because of the Fourier series representation (27), we actually solve a problem where this domain is periodically repeated along , but the zones of potential rapid slip accumulation, which coincide with the steady-state velocity-weakening regions of , are separated enough (by at least ) to prevent significant influence of spatial replications on each other. For our examples here, we take the constitutive law in one of the standard laboratory-derived forms of rate- and statedependent friction (7) with the Dieterich-Ruina version of state variable evolution

power law creep at depth and other features are given by Rice and Ben-Zion [1996] and Ben-Zion and Rice [1997]. The assumed variation of and with depth, like in the work by Rice [1993], is shown in Figure 3a. This is consistent with the experimentally determined temperature dependence of by Blanpied et al. [1991, 1995] for granite under hydrothermal conditions, as mapped by them into a depth variation based on a San Andreas fault geotherm. Variation of with depth is discussed below. The effective normal stress is assumed in this example to vary with depth in a way that incorporates high fluid over pressurization at depth, according to MPa. In this distribution (shown in Figure 3b), is equal to overburden minus hydrostatic pore pressure at shallow depth (up to about 2.6 km), with transition to lithostatic pore pressure gradient with 50 MPa offset at depth. Figure 3b also shows the initial stress, which is the same for all the cases considered here. Since the friction law (28) is a particular (and widely used) version of (7), all our conclusions from sections 3 and 4 hold with , . The critical stiffness (9a) becomes . To find the critical cell size corresponding to (10), we need to determine the coefficient for our present model in the expression for the single-cell stiffness . As mentioned in section 3, . That was determined by performing the static elastic spectral analysis of unit slip at a single FFT sample point (and all its replications) and equating the resulting stress reduction there to . When we ignore the presence of a free surface, the calculation is straightforward and gives (where , which arises from the periodic replication, can be ignored for the large values of we consider). A full calculation, which we did numerically, also includes unit slip in the mirror cell and in all its periodic replicates. It resulted in values of very close to for all the cells except for the one adjacent to the free surface, where . Using , we estimate the critical cell size by

0 z © † C ä“ CQ ¡ | §  ž w yÇ …Ò0 6 AÈ Ã

'¢ & $  %





$ ä ‘§ 0 ™

 

z ÿþ ý ü S © !䕓 iû ‘ $ w

©

ÿý þ䕓 ü û ‘ zª!þªý•“ ü û © w ÿ º

¡

¡ ¦  § ™ 5 ô § È ¡  ¡ © §  D



ÿþ ý ü û !䕓 ih

 

D

& r ti e 3Rg c fc b Cq © ú g s…v  § q ” ) ªý D %©  $  ¡§ ‚ €x u p ig ec R•c fdb B 9A $   ¡ 7 6 ¡   ¡§ ©” ¦ ªý D %© § D § s@±© § Vd !§ 5 ' )© (' $  ¡



© ú g C…v  § q ‚€x u ¡ © 0

¡  !§ ™  I $| „†  ¡  ¡ © !§ h© !§ D  ¡  © !§ I

ÿþ ý ü $ ÿþ ý ü º z !䕓 iû Ô© !ª•“ ƒû „$ w



  § q

 #6

‘

„

5 D ¡ È ô  !§  ¡ © § D



D

t B  § q & r „3Rg c fc b Cq ie   ¡ ƒ© § 7 p ig ec R•c fdb ™ z S ÿþ ý © !ª•“ ü û „$ w º ¡ 0 ™ 5 9 |– D ž3  5 7 Y ¡ Ô© !§ D Y º ¡A

¢  £ Ԕý D  

  

ÿþ ý  !ª•“ ü û z ÿþ ý !䕓 ü û © !䕓 ü û „$ w º ÿþ ý º

§  "6 5 f7¤  !¡§ !© ¡(' 0

” ” ªý D 4© ! $  ¡§ D

  


 ¡ © § D

$ ‘

9

(29)

14

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS

a

0 (a - b) a -5

b

0 normal stress initial shear -5

Depth (km)

Depth (km)

-10

-10

-15

-20

-0.01

0

0.01 0.02 0.03 Value of (a - b) and a

0.04

and (like in the work by Rice [1993]), consistent Figure 3. (a) Depth-variable distribution of frictional parameters with the measured temperature and inferred depth variation of of Blanpied et al. [1991, 1995] for granite under hydrothermal conditions. (b) Depth-variable distribution of the effective normal stress (solid line) and initial shear stress (dashed line). and Higgs [1992] and Chester [1994] and is more explicitly proposed by Brechet and Estrin [1994] and Baumberger [1997]. To account in a simple way for backward jumps, which could not be neglected near , we solve (28a) for , identifying the factor , which then appears as the stress biasing of forward jumps, and replace it by to account for backward jumps too. This procedure, used by Rice and Ben-Zion [1996] and Ben-Zion and Rice [1997], replaces (28a) with the regularized form tional grid and wish to keep uniformly high for good numerical resolution. As is well known, it is not presently feasible to do computations with values of chosen in consistency with laboratory values of (which are typically in the few micron range and give in the range of 1 m). We thus take larger and but do keep small compared to other feature sizes in the model, such as the seismogenic depth, which in our model corresponds to the steady-state velocity weakening region, extending over km. We prescribe in the simulations, using values approximately equal to 0.94 km and km in the examples shown. The distribution of the characteristic slip distance is assigned from (29) based on the desired for cells with steady-state velocity weakening. For the above values of , this results in of order of millimeters to tens of millimeters, much larger than laboratory values but, as explained, needed to keep the size of possible to resolve. Such selection of has to be adjusted in the regions close to transition from the velocity-weakening to the velocity-strengthening friction, where, owing to near-zero values of , it results in very small values of . Resolution of these values would require extreme refinement of time stepping during dynamic rupture, when slip velocities are large, as follows from the time selection criterion (11). Thus we increase in those regions, effectively changing (increasing) the nucle-

(30)

™

$ h à§

 Q

™

ٵ

™

During forward sliding at rates of order the modification to is of order or less, where , and so this is a negligible change from (28a). As discussed before, the critical cell size from (29) is a very important parameter, both for the physics and numerics of the problem. In principle, can vary with depth, since it is determined by the depth-variable frictional properties and stress. In all our examples, however, we attempt to make uniform throughout the velocity-weakening depth range. The motivation is that we have a uniform computa-

ٵ

™

| A x–

ٵ

ٵ

ٵ

¢ž £  A ¨

™

µ sQ Ÿ µ

ٵ

 ºQ ž ÒShº 3

ٵ

4

0

™

™



ٵ

$ ä T§ $ ª X§

   9

8
¦ ž „Ó 5  ™Q 4Sh D § È ô Ÿµ  '  d0 Q !§ D Ÿµ

5 §A @  0 ' A d)CQ f§ &Å @ 5 D 5Q 7 A $ PÅ 5 Ñ  ADHE  #6 5 7 &Å @ D A ¤ I­ÈÇÆ GFCÄ D

B5 &Å @ A

z (0 Q ' $ § 

B5 &Å @ ¦d)CQ f'§ &Å w@  A $0 A

 ¡ © § 0

B

7
  ¡§ © (' 5 º– 1X Q 7 Ÿµ ¡  !§ 5 '

-15

-20

0

20

40 60 Stress (MPa)

80

100

 5

D

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS

15

a

0 h* = 0.94 km h* = 0.23 km -5

b

0

-5

Depth (km)

-10

Depth (km)

-10

-15

-20

0

10 20 30 40 Characteristic slip distance L (mm)

50

Figure 4. (a) Depth-variable distribution of the characteristic slip distance of the friction law, for the two cases considered. (b) Modification in the critical cell size caused by the adjustments (increases) to in the regions next to the transition from the velocity-weakening to the velocity-strengthening friction; is defined only for the velocity-weakening part of the fault zone. ation size there. As for the velocity-strengthening regions, keeping a particular value of is not of a concern, since nucleation cannot happen there. However, it is the resolution of in the velocity-strengthening regions that usually controls the size of time steps during essentially quasistatic phases of deformation. This is not surprising since during the quasi-static phases the velocity-weakening regions are stuck with near-zero velocities, while the velocitystrengthening regions are creeping with (much larger) slip velocities close to the plate velocity. Hence, while assigning in the velocity-strengthening regions, it is practical to take into consideration time step constraints (12). Keeping in mind the restrictions discussed and aiming for a simple and continuous distribution of , we use, for the examples in this paper, the distribution of shown in Figure 4a. Such an assignment of keeps constant and unmodified from 13.5 to 4 km depth. Figure 4b demonstrates the modification of caused by the adjustments in . Before we consider choices of numerical parameters, let us have a look at the response this model produces and our simulations are able to capture. Figures 5 and 6 show parts of the slip accumulation for the cases and 0.235 km, respectively. Figure 7 shows the maximum slip velocity histories for these two cases. From the data, as well as from the assorted slip velocity output, we notice that the velocitystrengthening region at the bottom of the fault is creeping, with roughly the plate velocity of 35 mm/yr, whereas the velocity-weakening region accumulates slip through dynamic failure events. The velocity-strengthening region near the free surface is thin and gets broken by strong ruptures coming from the bottom of the fault segment, but it also exhibits some creeping. Even on these plots, nucleation zones km), can be distinguished (especially for the case corresponding to several dashed lines plotted on top of each other. The slip accumulates quasi-statically at those regions for a long time, then the slip accelerates and dynamically expands from there. Notice that the nucleation size scales with , and in both cases the nucleation size is times larger than . The model earthquakes generally nucleate at the bottom transition between the velocity-weakening and velocity-strengthening regions, because it is there that the creeping region transfers stresses to the seismogenic depth, loading it up. By looking at the distance between the tips of the dashed lines, which are 1 s apart in time, we can estimate the rupture propagation velocities. These range from slower ones right after the nucleation up to almost the shear wave speed of 3 km/s further in the rupture development. Many of the model earthquakes are large and reach the free surface, sending a wave of slip back to depth. The provided examples show that our simulation algorithm deals very well with slow loading of the fault with the equivalent of the millimeters per year plate rate, slow

ž º¨ Ÿµ

¢ ž T–  

4

™

ٵ

™

ٵ

ٵ

S
ž º¨ Ÿµ

4

™

ٵ

¡  !§ ™

ٵ

ٵ

™



™

ٵ

R
™ ™

-15 h* = 0.94 km h* = 0.23 km -20

0.1

1

10

100

1000

Modified value of h* (km)

16

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS

0

-5

Depth (km)

-10

-15

-20

15

18

Figure 5. Accumulation of slip versus depth for the case km, . The solid lines are plotted every 5 years. The dashed lines are plotted above 18 km depth every second if the maximum velocity anywhere on the fault exceeds 0.001 m/s. The model response consists of large, essentially periodic events rupturing the whole fault. The nucleation size is km.

h* = 0.23 km, properly resolved run

0

-5

Depth (km)

-10

-15

-20

15

18

21 Slip (m)

24

27

30

Figure 6. Accumulation of slip versus depth for the case km, . The solid lines are plotted every 5 years. The dashed lines are plotted above 18 km depth every second if the maximum velocity anywhere on the fault exceeds 0.001 m/s. The model response consists of a repeated pair of larger and smaller events. The nucleation size is km. In comparison with Figure 5, note the difference in the system behavior and the change in the nucleation size as is changed.

A

ž

  T–

ž ¨

–

ٵ

ºµ fÌsQ Ÿ µ

XXYYXYXYXY XXYYXYXYXY XXYYXYXYXY XXYYXYXYXY XXYYXYXYXY XXYYXYXYXY XXYYXYXYXY XXYYXYXYXY XXYYXYXYXY XXYYXYXYXY XXYYXYXYXY XXYYXYXYXY XXYYXYXYXY XXYYXYXYXY XXYYXYXYXY XXYYXYXYXY XXYYXYXYXY XXYYXYXYXY XXXXX

ºµ …ÁsQ Ÿ µ

VW VW VW VW VW VW VW VW VW VW VW VW VW VW VW VW VW VW VW VWVWVVWW VWVWVVWW VWVWVVWW VWVWVVWW VWVWVVWW VWVWVVWW VWVWVVWW VWVWVVWW VWVWVVWW VWVWVVWW VWVWVVWW VWVWVVWW VWVWVVWW VWVWVVWW VWVWVVWW VWVWVVWW VWVWVVWW VWVWVVWW VWVWVVWW
h* = 0.94 km, properly resolved run 21 Slip (m) 24 27 30

¢ž 0  A ¨

ž º¨ Ÿµ

4

 ٵ

U U

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
h* = 0.94 km, properly resolved run h* = 0.23 km, properly resolved run

17

a
10 0.1

b

10

Maximum velocity (m/s)

1e-05

1e-07

1e-09 500 600 700 800 900

Simulated time (years)

Figure 7. Maximum slip velocity on the fault as a function of time for (a) km and (b) km. Individual rupture events collapse onto a straight line on this timescale of hundreds of years. The maximum slip velocities reached during dynamic ruptures are of the order of 10 m/s for larger events and 1 m/s for smaller events in the Figure 7b. The maximum slip velocity in between the dynamic events is m/s, which corresponds to the plate velocity of 35 mm/yr. Read 1e-05 as . event nucleation, and dynamic slip with kilometers per second rupture velocities and meters per second slip velocities. The results shown are well-resolved numerically. We now discuss the choices of numerical parameters that produce these simulations. 7.2. Parameter Selection for Well-Resolved Simulations As explained in sections 3 and 4, the choice of the parameter is crucial for the simulation stability, accuracy, and tractability. Together with the selection of and , this parameter determines the number of the discretization points required (and hence the problem size) and the cell size through and . If from the expression for the minimum time step (16) is kept unchanged, the parameter can also be considered as controlling the smallest time discretization allowed. For the examples here, we find that produces stable results that are closely matched by the results obtained with (hence convergence through grid reduction is achieved). provides less satisfactory resolution. As an illustration, consider Figure 8, where the slip rate history during a part of the second event of the sequence with km is shown for the point at 3 km depth. Although in a different problem the value of required may

be different, the above consideration illustrates how to approach selection of . We use of 1/2, giving . For the convolution truncation, we choose the elastodynamic time windows by specifying or with , and , confirming these values by varying the parameters to make sure the results do not depend on their choice. Note that since the replication period in our model is much larger than the region failing dynamically, is quite a large window for this model, and almost identical results can be obtained by using . Larger values of the elastodynamic time windows in this model tend to increase the maximum velocities achieved, by a small amount, but do not noticeably influence other features, such as the amount and distribution of slip or the rupture velocities. The fact that can be chosen as small as 4 in this problem provides a huge computational gain. Let us illustrate that for our example with km. To achieve , we need elements along the replication period. This is also the number of Fourier modes for which we have to compute the convolution integrals at each time step (and hence to store slip velocity history sufficiently back in time). For mode , with

ú  A  | § ‚  

ú  ‚ (Œw| § ¦

ú (Õ | §  ‚

” •“ É ” A Q ¡

¢ž 0  A 3

”“ | ¨  “ •Í’ žA º¨  Ÿµ

4 4!

Ù  Ú ƒØ  É

| £ ‰

ú  A  | §  ‚

ٵ

µ µ sQ Ÿ 1ÁÂ

º o²‚ ù

ž º¨ Ÿµ

‚ù

4

Ú ƒØ Û Ù

ºµ µ ²‰sQ Ÿ Õ Â

ú ( A )w| § ¦  ‚

BQ‘

ú š

`

0.001

Maximum velocity (m/s)

0.1

0.001

ٵ

ʵ SQ ٠lĵ

a b&

| †–

“ •Í’ ”“

Â

Ÿ CQ …ƒsQ  “ •F’ ”“ µ ‘Â  µ ‘

ºµ µ fšCQ Ÿ âÌÂ

Â

µ µ sQ Ÿ †“Â

AÂ

µ

`

1e-05

1e-07

1e-09 500 600 700 800 900

Simulated time (years)

 dVÂ

ž º¨

4

c b&

Ú ƒ£Û ÙØ

‘

 ٵ

|

18
6

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
h*/h = 80 h*/h = 40 h*/h = 20 5

4

Slip velocity (m/s)

Figure 8. Slip velocity history at 3 km depth as a function of time for the second event in the sequence with km. Zero time is chosen arbitrarily for plotting convenience. The resolution gives essentially the same results as . The case is less sufficiently resolved. and , the number of values stored at spacing is . If of a constant window were used for all modes, this would be the number of required values for each Fourier mode, and we would need two arrays of the size , one for the Fourier coefficients and the other for the values of the discretized kernel. At each time step we would have to use these arrays to compute the convolution integrals. With and mode-dependent time windows, the number of values needed to be stored and used at each time step decreases significantly for higher Fourier modes. We pack more than one mode into most columns of the arrays, getting, for this particular example, the size of the arrays down to , which results in a very substantial reduction, by more than a factor of 300, of memory and cpu time for doing convolutions. This is consistent with our order-ofmagnitude considerations in section 6. Our evolution time step selection follows criterion (11) and (15)-(17). With the choice of other parameters dis, cussed, the minimum evolution time step allowed, is s for the case km and s for the case km. We can understand why such small values are needed by recognizing, as we did establishing criterion (11), that slip in one time step must be comparable (and preferably smaller) than the characteristic slip distance of the friction law, to resolve the state variable evolution. Since the simulated is of order millimeters (Figure 4a) and

h* = 0.94 km, properly resolved run

1e+06

Evolution time step (s)

10000

100

1

0.01 500 600 700 Simulated time (years) 800

Figure 9. Values of evolution time steps (in seconds) plotted as a function of the simulated time in years for the case km. Variable time stepping works well, making the time step span more than 8 orders of magnitude in this simulation, consistently with the changes in the slip velocity.

ž º3

4

g ghghgh

hgh

 ٵ

g ghghgh

hgh

g ghghgh

hgh

g ghghgh

hgh

g ghghgh

hgh

g ghghgh

hgh

g ghghgh

hgh

g ghghgh

hgh

g ghghgh

hgh

g ghghgh

hgh

g ghghgh

hgh

g ghghgh

hgh

g ghghgh

hgh

g ghghgh

hgh

g ghghgh

hgh

g ghghgh

hgh

ºµ l‘CQ Ÿ µ

g ghghgh

hgh

g ghghgh

hgh

g ghghgh

hgh

ž º¨  Ÿµ

4

ž | s¨ Ú ƒ• É ÙØ

   52 f

z

| A¨

4  D  C2 © A   w 

©A

–

”“  “ ”•“F’ º@TsÚ ÙƒØ3Ûsú( § R| § ‚ “ •Í’ Q Ú ƒ• É ÙØ Q  ÙØ A S| lÚ ƒ£Û

ž º¨

4

µŸ A ÁCQ µ

 ٵ

e

3

2

1

8

8.5

9 Simulated time (s)

9.5

10

™

¢ž 0  A ¨

z S

 A D  52 © A   w 

Ÿµ ž º #3 –

µ ÁsQ Ÿ µ
º ½

‚ù ™

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
the largest slip velocities are of order meters per second, it is clear that the smallest time steps should be of the order of one thousandth of a second, as we have here. Note that once the dynamic propagation begins, large slip velocities at the rupture fronts determine the time step size, and it is , so almost always equal to the minimum time step that the rapid event propagation is essentially modeled using constant time steps. Preventing the time step from becoming smaller than may also mean that the state variable evolution is occasionally not ideally resolved right at the rupture peak, but in all the cases we have checked, not resolving only occasionally at the very tip of the rupture does not change the results in any significant way. For the slow deformation periods in between dynamic rupture events, the time steps taken are quite large. Figure 9 shows the values of evolution time steps for the case km. We see that the time steps span more than eight orders of magnitude in this simulation. Conditions (12) or (14) do a very good job restricting the time steps during slow deformation. If they are violated, the computation becomes corrupted. We can show this by relaxing the conditions several times, that is, by selecting the time step using (11) and (15) - (17) as before, but with coefficients increased by a certain factor, although still insisting that is not larger that 1/2. Figure 10a shows the maximum velocity for the case km and a factor of 2 increase, and Figure 10b for a factor of 5 increase, of the time steps in the sense discussed above. Comparing to Figure 7a, which shows the well-resolved response, we see that numerical instabilities start to appear at slow sliding velocities for modestly increased time steps (Figure 10a), while for the factor of 5 increase the response looks very complex, with numerous events of different maximum velocities (Figure 10b), all of which are artifacts of the improper time discretization. Figure 11 shows slip accumulation for a factor of 3 increase, m with some small ”events” appearing (like the one at slip) and aperiodic large events. Those features are caused by the improper resolution in time; the true response in this case is the periodic sequence of large events shown in Figure 5. Improper space discretization also produces artificial complexification of the model response, as discussed by Rice [1993] and Ben-Zion and Rice [1995, 1997].

19

8. Discussion
As the considered examples show, the algorithm presented here is capable of rigorous treatment of long-duration deformation histories with continuing aseismic creep slippage in velocity-strengthening fault regions throughout the loading period, with gradual nucleation of model earthquakes followed by dynamic propagation of ruptures, and

with rapid post seismic deformation after such events. The algorithm is formulated for general rate- and state-dependent friction laws, and the positive direct effect observed experimentally and represented by such laws is decisive for its success during long intervals with essentially quasi-static response and aseismic slip. The algorithm employs a number of important ideas. Separation of the stress transfer functional into static and dynamic parts localizes the effects of the prior deformation history in convolution integrals on slip velocity with rapidly decaying kernels. Truncation of these convolutions is justified by rapid decay in time of the convolution kernels and allows us to simulate long processes without the necessity to deal with all prior deformation history at each time step. Variable time stepping makes the number of time steps during slow deformation periods numerically manageable while still capturing the details of both the nucleation and dynamic propagation phases. Proper space and time discretization ensures reliability of the results which can be verified through space and time grid refinement. The methodology has been presented using the 2-D antiplane spectral formulation. The described procedures can be readily extended to the 2-D inplane and 3-D spectral formulations. They can also be applied at least in part to discretized models based on spacetime boundary integral formulations, as we briefly discuss in Appendix B, where we further suggest that such formulations could be founded on kinematic modeling input from more versatile methods like finite difference, possibly being practical in cases for which the spectral diagonalization of convolutions does not apply. The numerically most challenging parts of the algorithm are calculations of the convolution integrals and, in the spectral formulations, fast Fourier transforms (FFTs). If the elastodynamic time windows used to truncate the convolutions are long and truncation according to the mode does not shorten much the time windows for higher modes, then the convolution evaluations take most of the computational time and even minor optimizations of convolution evaluations can be very beneficial. If, on the other hand, the time windows are much shorter for higher modes, as we have been able to use for our depth-variable example here, then the FFTs start to use a comparable fraction of the computational time, and an efficient FFT routine can make a significant difference. It is important to ensure that the results of the simulation do not depend on the discretization and other numerical parameters. For example, we could conclude that the model response is complex (Figure 10b) or that the events are aperiodic (Figure 11), whereas better resolution in time leads, in this case, to a periodic sequence of large events (Figures 5 and 7a). The verification of the independence of the results on numerics can be done through establishing convergence

vy

4 †– |

vy

Ú ƒØ  É Ù

Ú ƒ• É ÙØ ž º¨

4

 ٵ

™

ž º¨

4

 ٵ

20

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
h* = 0.94 km, twice larger steps h* = 0.94 km, 5 times larger steps

a
10 0.1

b

10

Maximum velocity (m/s)

0.001

Maximum velocity (m/s)

0.1

0.001

1e-05

1e-05

1e-07

1e-07

1e-09 500 600 700 800 900

1e-09 400 420 440 460 480

Simulated time (years)

Simulated time (years)

Figure 10. Maximum velocity on the fault for km and insufficient resolution in time. (a) increased by a factor of 2, (b) increased by a factor of 5 (in both cases, is not allowed to exceed 1/2). In Figure 10a, numerical instabilities start to show; Figure 10b shows seemingly very ”complex” behavior, which is actually just a numerical artifact.
h* = 0.94 km, 3 times larger time step 0

-5

Depth (km)

-10

-15

-20

15

18

21 Slip (m)

24

27

30

Figure 11. Accumulation of slip versus depth for the case km with insufficient resolution in time (y increased by a factor of 3). The solid lines are plotted every 5 years. The dashed lines are plotted above 18 km depth every second if the maximum velocity anywhere on the fault exceeds 0.001 m/s. The response differs from that of the properly resolved run in Figure 5. Nonperiodic large events and even some smaller events appear, all of which are, in this case, artifacts of the inadequate time discretization.

pq pq pqpqpq pqpqpq
v

pq pq pqpqpq pqpqpq

pq pq pqpqpq pqpqpq

vy

pq pq pqpqpq pqpqpq

pq pq pqpqpq pqpqpq

pq pq pqpqpq pqpqpq

pq pq pqpqpq pqpqpq

pq pq pqpqpq pqpqpq

pq pq pqpqpq pqpqpq

pq pq pqpqpq pqpqpq

pq pq pqpqpq pqpqpq 4

ž º¨  Ÿµ

pq pq pqpqpq pqpqpq i

pq pq pqpqpq pqpqpq

pq pq pqpqpq pqpqpq 4

vy ž º3  Ÿµ

pq pq pqpqpq pqpqpq

pq pq pqpqpq pqpqpq

pq pq pqpqpq pqpqpq

pq pq pqpqpq pqpqpq

pq pq pqpqpq pqpqpq

pq pq pqpqpq pqpqpq

pq pq pqpqpq pqpqpq r

pq pq pqpqpq pqpqpq

vy

i

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
of the results (or at least of their qualitative features) as the parameters of the simulation are refined. This is often necessary even when the output looks smooth and plausible, as it can still be qualitatively different from the true response of the model (e.g., Figure 11 versus Figure 5). Usually, the plots of slip velocity or stress reveal much more about the numerical stability and convergence than their slip counterparts. The examples presented are based on a rather simple model, with the fault properties uniform throughout the seismogenic (steady-state velocity-weakening) zone. The response consists of periodic sequences of events (Figures 5 and 6), and such a regular response has allowed us to concentrate on developing a rigorous and efficient numerical procedure. We verify, under the conditions of much better resolution and wider parameter range, the result of Rice and BenZion [1996] and Ben-Zion and Rice [1997] that the dynamic effects alone are not sufficient to produce event complexity. Moreover, we find periodic response in some cases where the earlier studies have found, evidently due to insufficient numerical resolution, chaotic sequences of large events. In our consideration, the characteristic slip distance of the rate- and state-dependent friction is taken to be much larger than the laboratory values to achieve numerical tractability. To predict the model behavior with the laboratory-derived values of , it is important to observe trends as we decrease . Comparing Figures 5 and 6, with in the case of Figure 6 being 4 times smaller than that of Figure 5, we see that the reduction in the characteristic slip distance introduced small events at the base of the seismogenic zone. Further twofold reduction in produces a sequence of one large and one small event much like Figure 6 (although with smaller slip per event), but it is possible that still much smaller values of would introduce more elaborate sequences, with more small events. To produce realistically complex behavior, additional features have to be introduced. We would expect that adding strong fault heterogeneities in the form of (highly) nonuniform normal stress and/or frictional properties would naturally complexify the model response. Accounting for shear heating is also very important. It would introduce pore pressure development, which would add a second weakening mechanism due to the evolution of the effective normal stress. Interaction of two weakening mechanisms has complexified events sequences for a model studied by Shaw and Rice [2000], in a certain parameter range for the ad hoc type of friction law with two slip-weakening distances used. Another consequence of the shear heating would be temperature-induced time variations in the frictional properties, which may also contribute to event complexity. These problems, as well as other important problems such as the

21

earthquake nucleation process or patterns of rupture propagation in events nucleated naturally as a part of a sequence, can be studied within the methodology presented in this paper. Comparison of the larger and the smaller events for the case km supports the view that large events are just small events that run away and shows how the runaway can be prevented by prior stress release. Let us consider the 3-D plots of slip and slip velocity for individual events shown in Figures 12 and 13. The distribution of slip before each of the shown model earthquakes (Figure 12) reflects the slip in previous events, the creeping velocity-strengthening regions on both ends of the fault segment, and a clear nucleation zone, which actually extends long back in time. The slip velocity plots (Figure 13) show how the dynamic events develop. Once an event nucleates, two rupture fronts propagate in the opposite directions. One of them is arrested in the velocity-strengthening region at the bottom of the fault. In the case of the larger event (Figure 13a), the other rupture front reflects off the free surface and runs down, rerupturing the seismogenic depth, with dynamic waves of slip propagating on the surface of the rupture. The spikes on the rupture front are an artifact of the outputting and plotting procedure; for a given space location as a function of time and for a given moment in time as a function of space, the slip velocity profiles are smooth. The plotting procedure also reduces the maximum slip velocities achieved, owing to insufficient resolution of the image surface. In the case of the smaller event (Figure 13b) the rupture gets arrested long before it reaches the free surface. From the slip distribution in Figure 12b we notice that the smaller event fails to advance into the region of larger slip (and hence higher stress release) left by the previous (larger) event. Comparing Figures 12a and 12b, as well as Figures 13a and 13b, we notice that the slips and slip velocities during and right after the nucleation of the smaller event look just like the ones for the beginning of the larger event. This means that observing signals from the nucleation and beginning of such an event, we would not be able to tell whether the final size of the event will be large or small. We can use the developed methodology, which incorporates both truly slow tectonic loading and all dynamic effects, to evaluate simplified approaches. Let us consider two such approaches: (I) a procedure with truly slow tectonic loading but with a part of the dynamic effects (namely, dynamic stress transfers) ignored, and (II) a procedure with all dynamic effects incorporated but much faster loading. , To get a procedure of type I, we take which coincides with the quasi-dynamic approximation used by Rice [1993] and Ben-Zion and Rice [1995]. Figure 14 shows the results for the case km. Comparing

 ‚ …¶ | § ¦ ‚

ž º¨

4

 ٵ

¢ž £  A ¨

 ٵ

™

™

™

™

™

™

22
28.8

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS

a

Slip, m

24.0

25.2

26.4

27.6

64 60
-4

56

Ti

-12 -16 -20

De

, pth

-8

km

m e, s
52 48

b
Sep 20 18:53

Slip, m

25.2

26.4

27.6

28.8

30.0

GMT

Figure 12. Slip in individual20 17:48 from the sequence in Figure 6, km, for (a) a larger event and (b) the following GMT Sep events smaller event. The time axis spans 20 seconds, with zero time chosen arbitrarily for plotting convenience. The slip axis spans 6 m in both cases. Notice the clear nucleation zone that extends much further back in time. The smaller event in Figure 12b looks just like the beginning of the larger event in Figure 12a; it stops by not being able to advance into the higher slip/lower stress region in the middle of the fault. This supports the idea that large events are small events that run away due to favorable stress/strength conditions on the fault.

‚ xw u t £C€yIvbs

64 60
-4

m Ti

-12 -16 -20

De

, pth

-8

km

56

e, s
52 48

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS

23

3.0

a

Slip velocity, m/s

0.0

0.6
64 60 5 Tim 6 e,

1.2

1.8

2.4

s

52 48

3.0

b

GMT

Sep 20 18:59

Slip velocity, m/s

0.0

0.6

1.2

1.8

2.4

-2

0

-1

6

D -12 ep th
-4 -12 -16 -20

De

p

k th,

Figure 13. Slip velocity inSep 20 20:07 events as in Figure 8, km, for (a) a larger event and (b) the following smaller GMT the same event. The spikes on the rupture front are artifacts of the outputting and plotting procedures; the slip velocity is actually smooth, as confirmed by plots of slip velocity at particular depths as a function of time and at particular times as a function of depth. Notice the dynamic effects in the form of waves, including the one reflected from the free surface. As in Figure 12, the smaller event in Figure 13b looks just like the beginning of the larger event in Figure 13a.

, k -8 m
-8

m

-4

‚ xw u t 0†y …„ƒs

64 60 56

Ti

m e, s
52 48

24
0

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
h* = 0.94 km, quasidynamic run (Tw = 0)

-5

Depth (km)

-10

Figure 14. Accumulation of slip versus depth for the case km with (i.e., the ”quasi-dynamic” approximation). The solid lines are plotted every 5 years. The dashed lines are plotted above 18 km depth every second if the maximum velocity anywhere on the fault exceeds 0.001 m/s. There is no reflected front of slip from the free surface which is present when the wave-mediated stress transfers are included (Figure 5). The slip per event and the slip and rupture velocities are smaller.
h* = 0.94 km, much larger loading velocity 0

-5

Depth (km)

-10

-15

-20

5

8

11 Slip (m)

14

Figure 15. Accumulation of slip versus depth for the case km with the much larger loading velocity mm/year. The solid lines are plotted every years. The dashed lines are plotted above 18 km depth every second if the maximum velocity anywhere on the fault exceeds 0.001 m/s. (The initial conditions could be tuned so that the slip in the middle of the lower velocity strengthening region is roughly the same as at its bottom. The qualitative features of the simulation are independent of the initial conditions.) Notice that, compared to the case of tectonic loading (• m/s, Figure 5), the nucleation phase is very different and the slip per event is more than twice smaller.

d ˜‚ u – 0w | ™£b$—I•

‚ £u $—–

wu‘ T“’

‰ xw u t Cˆ IT„ƒs

‰ xw u t Cˆ IT„bs

df ˜ ƒ&w | e‚

‡ ”

-15

-20

15

18

21 Slip (m)

24

27

30

17

20

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
with Figure 5, we see that the wave effects disappear as they should (there is no slip wave reflecting off the free surface), the rupture velocities are slower (the tips of the dashed lines are closer to one another), and the slip velocities are smaller too, as is the accumulated slip per event. We have examined the possibility that the quasi-dynamic calculations might approximately reproduce the slip per event as in Figure 5 if we could make the slip velocities faster during dynamic events. To that end, we did a series of studies in which the numerical value of in the radiation damping term of elastodynamic relation (26) was varied in size. Reduction to approximately half the proper value did give rupture propagation and slip velocities during events that tend to be of comparable order to those of the proper dynamic simulation (Figure 5), but there was very little effect on the slip per cycle, which remained much as in Figure 14. This suggests that the greater slip in Figure 5 is significantly due to the discussed reflected wave effect at the free surface, a feature which could not appear in the quasi-dynamic simulation. However, there is no qualitative difference in the system behavior, probably due to the relative simplicity of this model’s response. If the dynamic effects were more complex, their elimination may have made a more significant difference. Finally, we consider a procedure in which the plate loading rate is only a few orders of magnitude less than representative seismic slip rates, which allows to use standard elastodynamic numerical methodology throughout (like in the work by Shaw and Rice [2000]). To this end, we change the loading velocity from mm/yr in our implementation example above to mm/yr, for the case km, keeping other parameters the same. The resultant slip accumulation is shown in Figure 15, which is very different from the model response with the tectonic-like loading velocity (Figure 5). The increase in loading velocity has totally altered the nucleation process and location and m resulted in more than twice smaller slip per event (r slip per event, compared to m slip per event in the case of mm/yr). The loading velocities mm/yr and mm/yr also give quite altered (although, naturally, less so) nucleation, and slips per event of and 2.3 m, respectively. The evaluation of the simplified approaches shows that while accounting for the proper dynamic response can be very important, simulating slow tectonic loading is also crucial for uncovering the true model response. The results also demonstrate that even though simplified approaches may be unavoidable in some cases, they have limitations that can be uncovered and remedied only within a more general approach like the one presented here.

25

Appendix A: Derivation of Critical Stiffness and Time Discretization Constraints Dictated by Constitutive Law
in To derive expression (9) for the critical stiffness quasi-static slip, let us consider a spring-slider system moving steadily at rate , with state and frictional resistance . If we impose a perturbation on the sliding velocity and use , , and to denote the corresponding perturbations of other quantities, then from the friction law (7) we can write the linearized constitutive response to the perturbation as

where , , , and are used to denote the partial derivatives of the functions in (7): , , , with the derivatives evaluated at steady state, and given by (9b). The definitions of and are equivalent to the more familiar , . At steady state we have , and thus the sign of determines whether the system exhibits steady-state velocity weakening or strengthening. and can, in general, depend on the steadystate velocity , but in the most commonly used DieterichRuina logarithmic forms of the rate- and state-dependent friction laws, they are chosen as constants times effective normal stress (as we do in our implementation example) and are often denoted by and . The (quasi-static) elastic response to the perturbation is given by

where is the spring stiffness. Eliminating by combining (A1a) and (A1b), substituting , and defining, for simplification, nondimensional quantities , , and , we obtain the following simple system of two linear differential equations:

–  • t t – t †Š – i†—‰t g tŠ ‘ t k ‘t g u  •Cu• it vž‚Šº“ iw u 0¹x ‡„ ‡  © g ª žP5œkƒ•tg‚€¸€† • § –‡ ¨ •· i „· u ‘ « ® « © ª ¨ • ™ i gl„l™ © § u t ¶t vt†Š ¢µ´³©›²±¯°¯lœª® t­¬ll‚©ž’– 5'k0¢¡«Ž ƒ•)bŒ• ¦qžŠ t t †Š – t  $kt “ ‘¦w’5'k0¢¡ƒ•¥ƒ™ t  £b5'k0¢¡ƒ•†…™ u •™ i Ž g¤ u ™ i Ž g š t CŸkžŠu’“• ™Cœk0œp›•gvb™ •t iŽ š  –t “ t t †Š Ž– 9˜Œ‡ • tt •—‹ ‡  t  •“ • ’ƒ‡ t  ”t ‘u ‹ QŒt • Š u ‡Ž… t ‡ • t ‰‡ „

(A1b)

“ k ‡ˆ u ½

(A2b)

~}†{ |

– t Š‘ x Éʞ5žtitŠ 'k Š5žkt ȖǑt g ††5CkŠŠtt k “ »{{ ‘ ÆÄà Åu t “ ·k Ž Á ¾½ u Á†0·C5e¾0·· à k½ À À– t £‡k  t  u ¾

‡„

‡ˆ ‡  ‡„ $iuƒ•g ‚€€ T†…„ t „ u• t  ti ƒ•tg‚€€fu t  t u ‡ ˆ ¼‡ • t• u “ œk¸¿f¦ ‘u ‡Ž ˆ »–{ ’‰‡ „ Š

{

‡•

iw 0u ‘ hg ‰x |

v ˜ ‚ †w | u£u $—– •

d ˜ ‚ u – †w | q£!p$—I• ‚ u – 0on$—I•

zx | r y ˜ ‚ Cw | x£wu $—– • ‚ 0Tu t—– • s ITr x

myk C†lj

(A1a)

‰ xw u t Cˆ IQƒs

(A2a)

26

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
with defined by (A3). Considering eigenvalues (A6), we conclude the following. If , which corresponds to the situation with steady-state velocity weakening and improperly sparse grid (giving ), the eigenvalues have absolute values larger than unity, and the perturbation grows, consistently with the continuous case, regardless of the chosen time step . For , which corresponds to either steady-state velocity strengthening, or steady-state velocity weakening with a sufficiently dense grid, the continuous case predicts decay of the perturbation. Indeed, the discretized case follows this stable behavior if the time step is sufficiently small. That is, in the case , the absolute values of the eigenvalues is less than unity, and hence the perturbation decays, if the (dimensional) time step satisfies the conditions written in the text as equations (12).

(A3)

For velocity strengthening friction, , we have and the perturbation always decays. For velocity weakening friction, , the behavior depends on the sign of the expression in brackets of (A3), which defines the critical stiffness (A4)

(A5a)

or

(A5b)

with

(A5c)

(A6)

ð

éå ê

i¢h¿g é å  ¢i‚¿g å „

We would like this discrete dynamic system to have the same stability behavior as the continuous (in time) system (A2) does. The behavior of such discrete dynamic systems as (A5) can be understood by finding the eigenvalues of the governing matrix and comparing their absolute values to unity. The eigenvalues of matrix are given by

where subscripts and denote values discretized in space. Here is the static elastic stiffness matrix, and is a matrix of convolution kernels. These could emerge as the result of discretizing a boundary integral formulation, or could be determined by standard elastodynamic finite difference or finite element calculations. The calculations would be purely kinematic, imposing a step in slip at each node singly (say at node ), over a single elastodynamic time step, and calculating the stress histories at all other nodes . By analyzing the using (1) and (B1), the could be extracted numerically as the long-term limit of the stress changes and the determined from the transient response and tabulated. The relation (B1) is analogous to (2) and (4), and most ideas of the presented algorithm may be directly applied, including the truncation of convolutions, the schemes to advance through a time step, and the time step selection pro-

Ú¨¢i‚¿g  § ¨ òê § ñð Ô ÚŽï¨î°†·¢³i°î¿‚g ¿ ‘ ì í ë ‘ ¿ æª é ‘ u å é  ˆ ±i°îÈ´h¿g é å  "9™¢ih¿g é ˆ é å ê § ȗp¢ih¿g 5ä ±ç¯ è

h–Œã €

and the perturbation grows in time for and decays in time for , as claimed in (9a). To derive the constraint (12) on the evolution time stepping, we analyze, as a simple model case, explicit integration ) of the governing (with a constant time step system for the perturbation (A2). In the system (A2), inertia effects are ignored, but they are negligible at low slip rates. The discretized system is given by

Appendix B: Remarks on Extension of Methodology to Cases With No Translational Invariance
Our methodology can be extended to general cases which lack translational invariance, to accommodate problems such as a fault oblique to the free surface or heterogeneity of bulk material properties (like a layered Earth structure). If a boundary integral equation is discretized with sample points (nodes) over the fault domain, the stress transfer functional will still be related to the slip by an expression which can be put in the velocity form (B1)

wÍ Ö1Ë Ñ  – wÐÏ}~†{ ¬u “ 'k¸ižŒIt ά{ Ì | ‘ { t Š– ‘ g ‘ w Í tŠ‘ ÐÏ$i†…t g ¢i‚¿g å „

¿ Ñ wÍ —6Ë ¿ t• u  “ œk $ѐ¹Ñ

ñ

wÌ âáË

Qualitative behavior of such systems depends on the eigenvalues of the matrix . If the real parts of all the eigenvalues are negative, the solution vanishes with increasing time, and if at least one eigenvalue has a positive real part, the solution grows without limit. When the largest real part is zero, long-time periodic motion results. In the system (A2), the perturbation grows (exponentially) for and decays (exponentially) for , where

Ë

– Š  t 6‘ ‹ Ù t Š » Š k x É tž5k„¢Ñ5$iŠ v†$Ñt¥t – ‘ g †5†vCÑ Ñ 6{ ‘ ÚÙØ’× t †k  t Š k  “ {“ Å u Ô Ž Á ÒÒ × Á ÔlÓCÒÓ Ò ¶¾½ À Tu lC¶¾½ À Á Ñ i ‘ÔÓ Ž ÁÖÒÒ Ã ÕÑCœkiIÒÒ ½ nlCÒÒ g ¶¾½ À u C'kIu¾6nlC¶¾ƒ½g À ‘ÔÓ

wÌ —“Ë

¿ t•u “ œktÑuTÑ ~}| e{ {Í ~}| !#{ {Ì ©ª Ž « ¬)© É P5'k ›•“ g h€€ 0Œ• ~‘Æu t 6‘“ t – p~}C{ u| •· i „· Å Š – w Í tŠ‘ – 1ÏtižŒIt g wÍ (Ë wÐÌÏ$it†elt g Š‘ x É †Š6Ït – u Å{ t Š uË t ‘“ ‘ “ u wÍ 1ÎË Ž à „Ñ …t {Š ‘ Ü Ë Ù‰ †Ë Ùy ‘ º‹ Ù u Ü Ô­ Û ßÞ Ý “ × ×

Ã

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS

27

relation (B1) to

(B2)

where

Acknowledgments. These studies were supported by the Southern California Earthquake Center (SCEC), by USGS grant 99-HQ-GR-0025 to Harvard University, and by a Blaise Pascal International Research Chair of the Foundation of Ecole Normale Superieure, Paris, to JRR. SCEC is supported by NSF Cooperative

h–Œã €



summations on and extend over and on over . One may hope that truncation times for can be chosen shorter for the higher wave numbers represented by, say, the maximum of or . In implementing such an approach, one can determine and directly, without determining and first. The response to the slip history , which is a step function in time and Fourier mode (or mode of the wavelet expansion used) in space, would be determined, and hence values corresponding to would be known. Then

i h–€ 52g Ü 1)('h–€ 5 ã 0& ã

ÿ i gÒ £ ‘ u ùÒ Q’¬!CÎ¥ ê ~¨¢ih¿g  § ä

h–Ü € ã Ü  h–# € ã h–Ü € ã Ü 

h–€ ã5"“ — !ã u Û

h–Ü € ã

Ú¨¢ih¿g ¥  § ÿ i gÒ£ ‘ gÒ£ u g ù QBln¢i‚¿)lž¢i‚¿†Ò  ¢i‚¿†Ò ¥  CÎ¥¥ ê gù ùÒ ¢i‚¿)l£ gÒ h–$’¶‚–Œò˜ h–Ü € ã #€ ã u € ã
 é f å  u é å  Ú¨¢i‚¿g  § 

cedure. One exception is the notion of shorter truncation windows for higher mode numbers, which naturally arises in the spectral formulation and greatly reduces storage requirements and execution time compared to a uniform truncation window for all modes, as discussed in section 6. Unfortunately, there is no directly similar notion for the general case (B1). The following concept may, however, provide the extension of the spectral formulation to the general case (B1), enabling the possibility of shorter truncation windows for some time convolutions in (B1). We expect that even in the absence of translational invariance, if we read into (B1) a step-in-time slip distribution which varies spatially like , or like in 3-D modeling of a fault which spans two space dimensions and , then for high , or , there would generally be rapid decay of the transient stress in approach to the static limit. In such a case, the convolution could be truncated earlier than for lower or . Hence, the following procedure can be explored, illustrated further for the 2-D case. The slips can be represented as a linear combination of independent basis functions, each with progressively higher wave number features. While some wavelet expansion may turn out to be a preferred way of accomplishing that, a simple illustration is provided by using the Fourier sum like in (3) for the and . Using the notation of section 5 for this representation and its inversion, with and , we can transform the

fast transforms would be done on these computed values, to convert them to histories . From (B2), components of and would be given as and , for the value of used in the slip history selected. The procedure would have to be repeated for all values of . For cases without translational invariance, there is no simple transformation to diagonalize the matrix from (B1), and the matrix from (B2) will generally not be diagonal. For the problems that have translational invariance (or can be mapped to such by addition of a mirror imin the space-time formulation (B1) age), the matrix is still nondiagonal, even though it acquires a special form with . However, in this case the spectral formulation considered in section 2 diagonalizes the matrix of the convolution kernels. This translates into substantial reduction in the number of time convolutions and makes the spectral formulation of a translationally invariant problem (such as (2) and (4)) computationally much more efficient than the space-time formulation of the same problem, or than any formulations of problems that lack translational invariance. This is true despite the additional cost of the spectral formulation, which needs more degrees of freedom. Suppose that the replication distance of the spectral formulation is times larger than the domain which we would like to simulate (so that ), and that the time for computation of one truncated time convolution integral scales linearly with the number of space elements involved (as is true when all truncation windows are comparable to the time for the elastic waves to propagate through the simulated region). In a formulation with a nondiagonal matrix of convolution kernels, we need time convolutions to compute the values of the (discretized) functional for all nodes. This requires the number of floating point operations proportional to . In the spectral formulation, such as (3) and (7), due to the diagonalization, we need only time convolutions, one per each Fourier mode. This translates into the number of operations proportional to . The FFTs used in the spectral formulation add the number of operations proportional to , so that the overall number of operations in the spectral formulation is still proportional to . Since is usually small number, such as 4 to 8, and typical values for are 512 to 65536, comparing and immediately shows that the spectral representation of a translationally invariant problem is computationally much more efficient.

h–Ü € ã Ü  h–€ 5%Ï — !ã ã u

¢ih¿g å ÿ ÿ Cä ¢ih¿g CÒ  ù ùÒ †Îê é ˆ ¢ih¿uC  CÎ gù ùÒ ¥Ò ¥ÿ ê ýý ý ý ¢i‚¿g †Ò  i ÿ ™n’ Ž g‘ ù Ÿ¨†0h–Œ'㥠&C0‚–Œ„ §‘ yk € Žyk € ã ÿ Ÿ¨h–Œ'㎠٠§ € ñð Ž Ô ù f é © ~¨¢ih¿g é å  Ž é å ê § å †  é è å è )Ú¨¢i‚¿CÒ ¥   C‰¥ ê § Ò u g ù Ž ùÒ Ž ¨§î †¢· îi ‚¿g ù  ½ îi È´h¿g CÒ ¥  "9‘ ¿ ¿ ‘ ù ìíë

 

ü

¢i‚¿gù¶½†CÒÎ¦ê¥ ¤ ù è ‘u¡¢i‚¿’l£ ù gÒ ¡ ù ½ Ô ù f é   ù ¢u é ˆ

û

¡ å5ƒäå Ò   å ¢u Ò £ ¢i‚¿g é ä ¢ih¿g é ˆ

i µ ¥)Ó ² ‚õg Çó ú ù ø÷ö å ô

‚–€ ã ¢ih¿g ‚–Œã € ÿ"‹ ˆå Ü Ü { þ ý {ý ÿ"‹ Ü Ü { þ ý {ý i ‚õg Çó ¸÷±öå ô

28

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
1995. Gu, J.-C., J. R. Rice, A. L. Ruina, and S. T. Tse, Slip motion and stability of a single degree of freedom elastic system with rate and state dependent friction, J. Mech. Phys. Solids, 32, 167196, 1984. Horowitz, F. G., and A. Ruina, Slip patterns in a spatially homogeneous fault model, J. Geophys. Res., 94, 10,279-10,298, 1989. Morrissey, J. W., and P. H. Geubelle, A numerical scheme for mode III dynamic fracture problems, Int. J. Numer. Meth. Eng., 40, 1181-1196, 1997. Myers, C. R., B. E. Shaw, and J. S. Langer, Slip complexity in a crustal-plane model of an earthquake fault, Phys. Rev. Lett., 77, 972-975, 1996. Okubo, P. G., Dynamic rupture modeling with laboratory-derived constitutive relations, J. Geophys. Res., 94, 12,321-12,335, 1989. Perrin, G., J. R. Rice, and G. Zheng, Self-healing slip pulse on a frictional surface, J. Mech. Phys. Solids, 43, 1461-1495, 1995. Ranjith, K., and J. R. Rice, Stability of quasi-static slip in a single degree of freedom elastic system with rate and state dependent friction, J. Mech. Phys. Solid, 47, 1207-1218, 1999. Rice, J. R., Spatio-temporal complexity of slip on a fault, J. Geophys. Res., 98, 9885-9907, 1993. Rice, J. R., and Y. Ben-Zion, Slip complexity in earthquake fault models, Proc. Natl. Acad. Sci. U.S.A., 93, 3811-3818, 1996. Rice, J. R., and A. L. Ruina, Stability of steady frictional slipping, J. Appl. Mech., 50, 343-349, 1983. Ruina, A. L., Slip instability and state variable friction laws, J. Geophys. Res, 88, 10,359-10,370, 1983. Shaw, B. E., and J. R. Rice, Existence of continuum complexity in the elastodynamics of repeated fault ruptures, J. Geophys. Res., in press, 2000. Shibazaki, B., and M. Matsu’ura, Spontaneous processes for nucleation, dynamic propagation, and stop of earthquake rupture, Geophys. Res. Lett., 19, 1189-1192, 1992. Stuart, W. D., Forecast model for the great earthquakes at the Nankai Though subduction zone, Pure Appl. Geophys., 126, 619-641, 1988. Stuart, W. D., and T. E. Tullis, Fault model for preseismic deformation at Parkfield, California, J. Geophys. Res., 100, 24,07924,099, 1995. Tse, S. T., and J. R. Rice, Crustal earthquake instability in relation to the depth variation of frictional slip properties, J. Geophys. Res., 91, 9452-9472, 1986. Tullis T. E., Rock friction and its implications for earthquake prediction examined via models of Parkfield earthquakes, Proc. Natl. Acad. Sci. U.S.A., 93, 3803-3810, 1996. Tullis, T. E., and J. D. Weeks, Constitutive behavior and stability of frictional sliding of granite, Pure Appl. Geophys., 124, 383-414, 1986. Zheng, G., and J. R. Rice, Conditions under which velocityweakening friction allows a self-healing versus a crack-like mode of rupture, Bull. Seismol. Soc. Am., 88, 1466-1483, 1998. Zheng, G., Y. Ben-Zion, J. R. Rice, and J. Morrissey, A new method for calculating elastodynamic response over long stressing histories containing rapid rupture episodes (abstract), Eos Trans.

Agreement EAR-8920136 and USGS Cooperative Agreements 1408-0001-A0899 and 1434-HQ-97AG01718. This is SCEC contribution 511. We are grateful to John Morrissey for participation in earlier phases of this research and to Alain Cochard for extensive discussion and comments on the manuscript.

References
Baumberger, T., Contact dynamics and friction at a solid-solid interface: material versus statistical aspects, Solid State Commun., 102, 175-185, 1997. Ben-Zion, Y., and J. R. Rice, Slip patterns and earthquake populations along different classes of faults in elastic solids, J. Geophys. Res., 100, 12,959-12,983, 1995. Ben-Zion, Y., and J. R. Rice, Dynamic simulations of slip on a smooth fault in an elastic solid, J. Geophys. Res., 102, 17,77117,784, 1997. Blanpied, M. L., D. A. Lockner, and J. D. Byerlee, Fault stability inferred from granite sliding experiments at hydrothermal conditions, Geophys. Res. Lett., 18(4), 609-612, 1991. Blanpied, M. L., D. A. Lockner, and J. D. Byerlee, Frictional slip of granite at hydrothermal conditions, J. Geophys. Res., 100, 13,045-13,064, 1995. Boatwright, J., and M. Cocco, Frictional constraints on crustal faulting, J. Geophys. Res., 101, 13,895-13,909, 1996. Brechet, Y., and Y. Estrin, The effect of strain rate sensitivity on dynamic friction of metals, Scri. Metall. Mater., 30, 1449-1454, 1994. Chester, F. M., Effects of temperature on friction: Constitutive equations and experiments with fault gouge, J. Geophys. Res., 99, 7247-7261, 1994. Chester, F. M., and N. Higgs, Multimechanism friction constitutive model for ultrafine quartz gouge at hypocentral conditions, J. Geophys. Res., 97, 1859-1870, 1992. Cochard, A., and R. Madariaga, Dynamic faulting under ratedependent friction, Pure Appl. Geophys., 142, 419-445, 1994. Cochard, A., and R. Madariaga, Compexity of seismicity due to highly rate dependent friction, J. Geophys. Res., 101, 25,32125,336, 1996. Cochard, A., and J. R. Rice, A spectral method for numerical elastodynamic fracture analysis without spatial replication of the rupture event, J. Mech. Phys. Solid, 45, 1393-1418, 1997. Dieterich, J. H., Modeling of rock friction, 1. Experimental results and constitutive equations, J. Geophys. Res., 84, 2161-2168, 1979. Dieterich, J. H., Constitutive properties of faults with simulated gouge, in Mech. Behavior Crustal Rocks, edited by N. L. Carter et al, Geophys. Monogr. Ser., vol. 24, pp. 103-120, AGU, Washington, D.C., 1981. Dieterich, J. H., Earthquake nucleation on faults with rate- and state-dependent strength, Tectonophysics, 211, 115-134, 1992. Dieterich, J. H., A constitutive law for rate of earthquake production and its application to earthquake clustering, J. Geophys. Res., 99, 2061-2618, 1994. Geubelle, P. H., and J. R. Rice, A spectral method for 3D elastodynamic fracture problems, J. Mech. Phys. Solids, 43, 1791-1824,

LAPUSTA ET AL.: DYNAMIC ANALYSIS OF SLOWLY LOADED FAULTS
AGU, 76(46), Fall Mtg. Suppl., F405, 1995.

29

3 4

N. Lapusta and J. R. Rice, Division of Engineering and Applied Sciences and Department of Earth and Planetary Sciences, Harvard University, Cambridge, MA 02138. (lapusta@esag.harvard.edu; rice@esag.harvard.edu) Y. Ben-Zion, Department of Earth Sciences, University of Southern California, Los Angeles, CA 90089-0740. (benzion@topaz.usc.edu)

G. Zheng, Technology Group, IBM, Somers, NY 10589. (gutuan@us.ibm.com)
Received October 18, 1999; July 13, 2000. revised May 31, 2000; accepted

A This preprint was prepared with AGU’s L TEX macros v5.01, with the extension package ‘AGU3 ’ by P. W. Daly, version 1.5e from 1997/11/18.



Add Comment

To add a comment, you need to login or signup (you can signup without leaving this page)