ML13317B472
| ML13317B472 | |
| Person / Time | |
|---|---|
| Site: | San Onofre  |
| Issue date: | 12/31/1979 |
| From: | Day S SYSTEMS, SCIENCES & SOFTWARE |
| To: | |
| Shared Package | |
| ML13317B473 | List: |
| References | |
| TASK-03-06, TASK-3-6, TASK-RR SSS-R-80-4295, NUDOCS 8005010572 | |
| Download: ML13317B472 (75) | |
Text
SYSTEMS, SCIENCE AND SOFTWARE SSS-R-80-4295 THREE-DIMENSIONAL FINITE DIFFERENCE SIMULATION OF FAULTDYNAMICS S. M, DAY FINAL REPORT SPONSORED BY:
NATIONAL AERONAUTICS AND SPACE ADMINISTRATION AMES RESEARCH CENTER M oFFETT FIELD, CALIFORNIA 94035 CONTRACT No, NAS2-10459 S PROJECT NO, 12035 REGULATORY goXE FIL CODY DECEMBER 1979 P. 0.
BOX 1620, LA JOLLA, CALIFORNIA 92038, TELEPHONE (714) 453-0060 8005O1O
- SYSTEMS, SCIENCE AND SOFTWARE SSS-R-80-4295 THREE-DIMENSIONAL FINITE DIFFERENCE SIMULATION OF FAULT DYNAMICS S. M. DAY FINAL REPORT SPONSORED BY:
NATIONAL AERONAUTICS AND SPACE ADMINISTRATION AMES RESEARCH CENTER MOFFETT FIELD, CALIFORNIA 94035 CONTRACT No, NAS2-10459 S PROJECT No1 12035 DECEMBER 1979 MAIL ADORESS:
P. O. BOX 1620, LA JOLLA, CALIFORNIA 92038 TELEPHONE: (714) 453-0060 STREET ADDRESS:
3398 CARMEL MOUNTAIN ROAD, SAN DIEGO, CALIFORNIA 92121
TABLE OF CONTENTS Section Page I.
INTRODUCTION.
1 II.
FIXED RUPTURE VELOCITY MODELING......
3
2.1 INTRODUCTION
3
2.2 DESCRIPTION
OF THE CALCULATIONS.
3 2.3 NUMERICAL RESULTS.....
.6 2.4 APPLICATION TO GROUND MOTION ESTIMATION..
22 2.4.1 Stress Drop Scaling.....
22 2.4.2 Scaling With Fault Width.
23 2.4.3 Model Plausibility.....
25 III.
SPONTANEOUS RUPTURE MODELING....
31
3.1 INTRODUCTION
31 3.2 FAULT CONSTITUTIVE MODEL...
31
3.3 DESCRIPTION
OF THE CALCULATIONS...
31 3.4 NUMERICAL RESULTS FOR UNIFORM STRESS 35 3.5 NUMERICAL RESULTS FOR NONUNIFORM STRESS................... 43 3.6 DISCUSSION..
58 IV.
SUMMARY
61 V.
REFERENCES.
62 APPENDIX A:
PROBLEM FORMULATION.........
.... 65 APPENDIX B:
INITIATION OF RUPTURE......
70
I. INTRODUCTION Deterministic computational simulation of earthquake ground motion, directed toward developing engineering design motion criteria, is an emerging technology in earthquake engineering.
In current practice, ground motion simulations often utilize rather sophisticated computational methods to model elastic wave propagation from source to site (for ex ample, Wiggins, et al., 1978; Swanger and Boore, 1978; Apsel, 1979).
On the other hand, the earthquake source itself is usually specified on largely intuitive grounds.
The displace ment discontinuity time history (slip function),
by means of which the seismic source is represented, is generally pre scribed without rigorous consideration of fault dynamics.
Since computed ground motion depends linearly on the prescribed fault slip, uncertainty in the slip function translates to un certainty in predicted ground motion.
Although the process of seismic faulting is complex and poorly understood in detail, some of the assumptions governing the prescription of the slip function for ground motion computation can be tested with computer simulations of the fault dynamics. The objective of the current study is to help define the sensitivity of earthquake slip functions to fault geometry, frictional strength, and prestress configura tion. The primary analytical tool for the study is the I4TRES computer code operating on the ILLIAC IV computer. I4TRES is a three-dimensional finite difference code which exploits the parallel processing capability of the ILLIAC. Three-dimensional simulations utilizing one-half to three-quarter million finite difference nodes can be performed with a few hours computing time on this system.
Six detailed finite difference earthquake simulations were specified for the study in the original statement of work.
Subsequent improvements in the finite difference code reduced 1
the computational requirements substantially, however, and several additional calculations were performed. in addition, a number of calculations performed under separate funding were available, and these have been analyzed as well.
In all, ten rather detailed three-dimensional calculations have been utilized for this study.
The first four of these calculations were directed to ward understanding the influence of fault width on fault dy namics. For these calculations, we have restricted consider ation to the case in which faulting occurs at constant rupture velocity and uniform stress-drop (prestress minus frictional stress).
Material properties were linearly elastic and the frictional behavior of the fault was such that the solutions scale linearly with the stress-drop.
The other six calculations assess the effect of the frictional strength and the prestress configuration in con trolling fault slip.
The frictional behavior is governed by a simple constitutive model which includes a frictional strength and a displacement-weakening law.
Rupture velocity is not a prescribed constant, but is derived from the dynamic solution --
that is, rupture is "spontaneous."
Off the fault plane, material behavior is constrained to be linearly elastic.
In Section II, the constant rupture velocity modeling is discussed.Section III discusses the results of the soon taneous rupture calculations using the displacement-weakening model, which is described in Appendices A and B.Section IV summarizes the results of the project and the implications for deterministic ground motion studies.
2
II.
FIXED RUPTURE VELOCITY MODELING
2.1 INTRODUCTION
Closed-form theoretical solutions describing fault mechanics are available for only the most idealized of prob lems.
Perhaps the most useful of such analytical results is the solution of Burridge and Willis (1969), giving the slip history on an elliptical shear crack which initiates at a point in a prestressed elastic wholespace and grows without bound at a constant (azimuth-dependent) rupture-propagation velocity. This solution is characterized by singularities at the crack edge in both shear stress and velocity, and the strength of these singularities grows as the crack dimension increases.
In addition, this so-called "self-similar" solu tion cannot account for the effect of rupture termination.and the subsequent arrest of sliding on the fault.
More recently, numerical methods have been applied to investigate dynamic models of faulting in which rurpture is confined to a rectangular surface (Madariaga, 1977; Archuleta and Day, 1980).
While these numerical solutions have resolved only low-frequency (less than 1 Hz) behavior of the fault slip, they demonstrate that edge effects, associated with the narrow dimension of the fault, substantially influence the slip function.
In this section, we assess the influence of fault width on fault dynamics, with emphasis on the high-frequency behavior which is of primary engineering interest. For simplicity we have restricted consideration to the case in which faulting occurs at constant rupture velocity and uniform stress-drop.
2.2 DESCRIPTION
OF THE CALCULATIONS We performed four three-dimensional finite difference earthauake simulations for this study.
The problem formulation 3
is outlined in Appendix A.
Each calculation was for a rectan gular fault surface in a uniform whole space, with rupture initiated at the center of the rectangular fault, Figure 2.1.
The prestress direction was aligned with the long dimension of the fault in each case.
The following parameters were employed for all four cal culations:
P wave velocity a=
6.0 km/sec S wave velocity
= 3.46 km/sec Density p
= 2.7 gm/cm3 Rupture velocity VR = 3.12 km/sec Tectonic shear stress a = 1,000 bars Frictional stress C = 900 bars Three of the calculations were for a fault length Z of 8 km and fault widths w of 1.5, 4, and 8 km, respectively. The fourth calculation was for a fault length of 16 km and a fault width of 4 km.
A number of physical quantities compete in the seis mological literature for the appellation "stress drop."
Here we reserve the term stress drop, denoted by da, for the quan tity a a ; the shear prestress minus the dynamic frictional stress (a -
a has also been called the "effective stress").
The shear prestress minus the final static shear stress on the fault, which may differ from ac, will be referred to as the "static stress drop."
The numerical solutions scale with the prescribed stress drop,
,a. _Equivalently, the solutions scale with average computed static stress drop, which differs very little from da in these calculations.
Results in the figures have been scaled for da = 100 bars.
4
9 of 1
\\
X I
Figure 2.1. Rupture geometry and coordinate system.
I 5
I
In each case, the medium in the vicinity of the fault plane was represented by cubic finite difference zones of 0.15 km on a side. This mesh refinement is sufficient to retain accuracy for frequency components up to about 5 Hz.
The numerical grid was large enough that no reflection from the exterior grid boundary returned to the fault zone during the calculation.
The formulation of Appendix A was simplified somewhat for these finite difference calculations by prohibiting fault slip in the direction normal to the prestress direction. This restriction does not have a significant effect on the solution.
It is known, for example, that for the self-similar, expanding elliptical crack, the slip is everywhere parallel to the shear prestress. Furthermore, Madariaga (1976) has shown that the perpendicular component of slip is very small for the case of a finite circular crack. A second simplification was that slip was not permitted to recommence at a point once a slip velocity zero first occurred at that point. This simplifica tion has little effect, since Day, et al. (1978) have shown that most of the fault plane experiences a stress relaxation, rather than reloading, following the arrest of sliding.
(This may not be true for highly nonuniform stress-drop models, and this simplification is not employed in the calculations of spontaneous rupture discussed in Section III.)
2.3 NUMERICAL RESULTS We first examine the low-to intermediate-frequency character of the earthquake slip function, which can be sum marized by the static slip and the rise time. Emphasis is on analysis of the three calculations with Z = 8 km, since the fourth calculation (w = 4 km, Z = 16 km) was not run to com pletion, i.e., fault slip velocity was not everywhere zero when the calculation terminated.
6
Figure 2.2 shows the calculated slip histories for selected points along the fault length for the square (8 km x 8 km) fault. The final static offset is greatest near the center of the fault, decreasing as the observation point approaches an edge.
The rise time of the slip function is also greatest at the center, where it equals approximately the fault length divided by the shear speed. The rise time de creases as the edges of the fault are approached. Arrest of slip occurs first at the outer edge and propagates inward.
The center of the fault is the last point to heal.
Figure 2.3 shows slip histories calculated along the fault length for the 4 km x 8 km simulation. The time his tories in this case are similar to those in Figure 2.2, but the rise time and static slip near the center of the fault are reduced. Apparently, the tise time and static slip are con trolled by the distance to the nearest edge of the fault.
Figure 2.4 shows slip histories for the 1.5 km x 8 km case. Again, it is clear that fault width controls static amplitude and rise time. From Figure 2.4, we see that beyond a distance of about one fault width from the hypocenter, the slip function for the case w = 1.5 km remains essentially constant in shape with increasing distance.
Figure 2.5 shows slip histories for the simulation with w = 4 km, 2 = 16 km.
The important point here is that, again, the slip function shape is nearly invariant with distance, beyond one fault width from the hypocenter.
Figure 2.6 shows the static slip along the fault center line for w = 4 km and w = 1.5 km (for 2 = 8 km).
The hor-i zontal lines show the static solution for an infinitely long fault (Knopoff, 1958).
Except near the end of the fault, the static slip is essentially constant along the length of the fault, and is very well predicted from Knopoff's static solution.
7
.3.0
=
0.3 2.0 0.75
- 1.
0.
2.25 3.0
- 3. 75 1.0 0
1.0 2.0 3.0 TIY.!
(sec)
Ficure 2.2.
Slip histories along the center-line o; the 8 k1 x 8 ca fault plane, at several distances from the focus.
The observation points are aligned with the direction of sliz.
08
i.e6 1.4 1.2 1.0 03
- 0. G
- 0.4
_3 0.22
- 0.
3.70 0
1.0 2.0 3.0 Figure 2.3.
S 1 i; hi st-6:or S along the center-Il-e o~the 4 k.- x 8 '-=.
fault ?!anre.
The obser~iaticn ::cints ar-e aligned with the sli; direction and the 1on-c di.ensiom of:
t-he rault.
9
ii P 4I I:fl
-fl
(
1 in tn tu
-1
'I-, O Il, 14 I~~
44t I,)~a (3.U rU 0
r(
I II) ~Ii)*
I,
-$1FCDP 1v'.
C) 14
.- ri
1.4 1.4 r
1.95 1.00 4.05 1.2 4.95 6.0 1.0 7.05 0.8_
H
/*/
Ui) 0.6
- I I
0.4 a
0.2.
0 1
2 3
TIME (sec)
Figure 2.5.
Slip histories along the center-line of the 4 km x 16 km fault plane.
The observation points are aligned with the slip direction and the long dimension of the fault.
1 4 1.2 0
w.= 4 kma o9 1.0a AG O
0.6o o0
. 4 W
1
- 1. 5 kmn 0.2 I
I I
1.0 2.0 3.0 4.0 DISTANCE (km)
Figure 2.6.
Static slip along the center-line of the fault, as a function of distance from the focus, for the 4 km x 8 km and 1.5 km x 8 kin faults.
Figure 2.7 shows the slip rise times along the center line for w = 4 km and w = 1.5 km (Z = 8 km).
Rise time was defined to be the time required to attain 90 percent of the static value of slip. The horizontal lines represent the time for a shear wave to travel from the edge to the center line. For w -
1.5 km, the rise time at first decreases with distance from the hypocenter, then approaches a constant level of w/2a. For w = 4 km, the rise time again decreases with dis tance from the center, but the effects of the end of the fault interfere to further reduce the rise time before a constant level can be established.
We turn our attention now to the high frequencies. Here it is appropriate to focus on the slip velocity function, and particularly on the peak slip velocity. Figures 2.8 through 2.11 show the slip velocity time histories along the center lines of the four fault models.
The slip velocities have been low-pass filtered to remove frequencies in excess of 5 Hz, which is close to the highest frequency that can be reliably computed in the finite difference mesh with the 150 m zone size.
It can be seen from Figures 2.8 through 2.11 that the peak slip velocity increases with hypocentral distance. From Figure 2.10 (1.5 km x 8 km) and Figure 2.11 (4 km x 16 km),
however, it is evident that the rapid growth in peak velocity ceases beyond a hypocentral distance of about one fault width.
Figure 2.12 illustrates this phenomenon more clearly. The broken curves represent the peak low-passed (5 Hz) velocities obtained along the fault center-line for the 1.5 km x 8 km and 4 km x 16 km simulations, respectively.
In both cases, the peak slip velocity apparently approaches a uniform level with increasing hypocentral distance.
Figure 2.13 shows slip velocities across the fault width, at a fixed distance of 6 km along the fault length, for the case w = 4 km, 2 = 16 km.
The figure illustrates the near uniformity of peak slip velocity across the fault 13
k-0. 5 0)w 1.4 kin 1
0 0
.04.
0ISVN U
kk 40i i
n.5 kin,.x.
0 i als
C)(V (1)
C in '-I v, 'Lj.c in~
in v
n 4
-4 (n)
(1) 0 5,
q) o 4-i in :i; Inn if)
U-nl tj
)i I
t
()
11i I.1 0
l l *r L--
4 in 0,
F-'
(U a) 0
't CV.
- 0)
[(0.
41J F) r:' a))~
in 44~
34-1
-,1
- 4) o 1:3
) 1x AJrMIR 110 -lniin
3.0 r = 3.75 3.0 2.25 00.75
\\3 I
0 0.2 0.4 0.6 0.8 TIME (sec)
Figure 2.9.
Slip velocity histories obtained along the center of the 4 km x 8 km fault, at several distances from the focus.
Observation Points are aligned with the slip direction and the long dimension of the fault.
The time histories have been low-cass filtered wi t
- a. 5 Hz cutoff and the time delays correspondin-to ructure arrival have been removed.
16
SL IP VIL.0CT1Iy (itse IiD (D
t~0 1 ti 0 Jit 0
(1)~
4 to (D
I-i H Dti iClurrt (D
I-f-(D MD HI-1b0.0 1-I -~0 1".
fl -
-C
('iD 0
H Ut 0 t-CO ul fI rL) L) 0 C
HA (UND (t U) J
- 1-1 iZ 1-
-4 rp w (D ft J)
Ura 0 fl-P) ft H 0 u
WH)I 0
'.J cl-0, 1 1ik U)t 0
f P.~
(D O f.
D
<H C.
( - U)0)
- i r
-, 0 u0 -1 u
W~C 0' to 110l(DI it)-
L1 Ibk (I
(D 0. 1-1 1-k to(D (Dfu P-S
- 0) 0 - U) t ON (D
,(
W O. 0 1-1 1 0 7 Il i Cl)(
fu 10 '] (D 0 :3 u) Ul 0
() D C)
4 4.95 6.0 70 4.05 33.0 o
r
=1.95
'I I
>4i 2
7.0 r~l.4.95
- 6.
4.0 0
3 3
r 1
fault Th tiehsoishvtenlwps I
- I S
~
.\\
fitee wit
- a.
5zcuof is4 1
1 0
0 1
2 3
TIME (sec)
Figure 2.11.
Slip velocity histories along the center of the 4 kmn x 16 km fault, at several distances from the focus.
Observation points are aligned with the slip direction and the long dimension of the fault.
The time histories have been low-cass filtered with a 5 Hz cutoff.
18
0 I
I I
I I
I 5
Circular Crack 4
>1 H
U 3!
2 4."
~~~
w = 1. 5 km I
0 I
I I
I II 0
1 2
3 4
5 6
7 8
DISTANCE (ki)
Figure 2.12.
Peak slip velocity as a function of focal distance, alonqthe center of the long dimension oC the fault.
The large-dashed curve is for the 1.5 km x 8 km case, the small-dashed curve is for the 4 km x 16 km case.
The solid curve represents the analytic solution for an expanding circular crack.
Peak velocities were determined from low-pass filtered slip histories with a 5 Ilz cutoff.
(D SLIP VELOCITY (m/sec)
U- (
-1 0 (
(0nU) 0
~0 H-C (CD (D
ft ol (D.
r0 fb (1 (1) 0 ti
()
ulUt Cl)
(D (D
DH (D
width, for hypocentral distances greater than w. The slight decrease in peak velocity very near the fault edge may be due to the slip function rise time near the edge becoming compar able to that of the 5 Hz low-pass filter.
We can interpret these results for peak slip velocity using the closed-form analytic solution of Kostrov (1964) for an expanding circular crack. The analytic solution for the velocity discontinuity A on a circular crack expanding at a uniform rupture velocity is
. a 7 + r/VR VT (T + 2r/VR) where T is the reduced time (time minus rupture arrival time),
r is the distance from the center of the crack to the obser vation point, and C is a constant which equals 0.81 for a Poisson's ratio of 0.25 and rupture velocity of 0.98 (Dahlen, 1974).
This solution is singular at the rupture arrival time (except at the hypocenter), and approaches C for T large compared to the rupture travel time.
In order to compare this to the numerical solution, we approximate the effect of a low pass filter by averaging the analytic solution over a "cutoff" period 1/f. The resulting expression for peak slip velocity is
= C 1 (2rf/V + 1)1/2 (2.2) pa R
1/
which, for >> V /r, is proportional to r /2 That is, in the absence of edge effects and nonlinearities, peak slip ve locity would increase as the square root of distance from the initiation point.
Equation (2.2) is plotted as a solid line in Figure 2.12.
From Figure 2.12 we see that edge effects do not act to modify the peak velocity in the numerical solution within a hypocentral distance of about one fault width; up to that distance, the behavior of peak slip velocity closely follows 21
that of Kostrov's expanding circular crack solution. The in fluence of the fault edges is to terminate this growth of peak velocity at a distance r = w.
In summary, then, the peak slip velocity increases with distance from the focus up to distance of approximately one fault width, then remains nearly uniform at that level over most of its length.
Thus, over most of its length,.a long, narrow fault will have a peak slip velocity, after low-pass filtering with cutoff frequencY f, of approximately Aa
/V+11/2 s
(2wf/V+
1)1 (2.3)
This.implies that peak velocity is proportional to do and proportional to the square root of fault width. For the case of VR = 0.9a, and when f >> VR/w, Equation (2.3) implies a value of approximately /
AT 7
/pS for the peak slip velocity.
2.4 APPLICATION TO GROUND MOTION ESTIMATION We wish to assess the applicability of these numerical results to strong motion modeling. First, we discuss the scaling of the solutions with stress drop.
Second, we discuss the scaling with fault width and examine the high-frequency behavior of the slip function, which is of particular impor tance for design motion estimation.
Finally, to test the plausibility of the results, we apply the proposed scaling.
to a source-receiver configuration similar to that of the Station 2 accelerograph for the 1966 Parkfield earthquake.
2.4.1. Stress Drop Scaling The solutions presented in Section 2.2 scale directly with the stress-drop Aa. The static stress drop, which is a function of position, also scales directly with Aa.
Actually, for the simple models treated thus far, static stress 22
drop is nearly equal to Aa over most of the fault. Thus, the solutions scale directly with static stress drop, a quantity which has been widely estimated for earthquakes.
However, particular attention should be given to some of our simplifying assumptions when relating the solutions to seismic estimates of static stress drop.
First of all, we have neglected the effect of finite material strength. This effect, if it has an important role for frequencies of engineering interest, could cause deviation from simple stress drop scaling.
Secondly, the model assumes a simply-connected fault plane and coherent rupture propagation. As a result, the static stress drop approximately equals da over most of the fault plane, and the static stress drop inferred from final offset and gross fault dimension is representative of the actual static stress drop (see Figure 2.6 for example).
If, on the other hand, the rupture had left unbroken patches, or had progressed irregularly, the static stress drop inferred from final offset divided by gross fault dimension would have been lower. In other words, the results should be interpreted with the understanding that static stress drop as measured seismically is not necessarily representative of the stress-drop a a=.
2.4.2 Scalina with Fault Width If the aspect ratio w/ is held fixed, then each of the four fault slip functions obtained in Section III can be scaled to correspond to different fault widths, provided the time axis is also scaled. The displacement axis scales linearly with w, the time axis scales as w, and velocity is consequently independent of w. However, the maximum frequency represented in the solution will go as 1/w, so if we scale, say, the 1.5 km x 8 km fault slip function to represent, say, 23
a 9 km x 48 km fault, the scaled numerical solution would re solve frequencies up to only 0.83 Hz.
Fortunately, we can do much better than these rigorous scaling relations suggest. First, we have noted that for Z >> w, the slip function is nearly independent of distance along the fault length for distances greater than w. We can thus reasonably apply the results for the 1.5 km x 8 km fault to any long, narrow fault. Secondly, we have noted that the shape of the initial slip-velocity peak is approximated very well by the square root singularity of Kostrov's solution, Equation (2.1), for focal distances less than the fault width, and by
/T(Z + 2w/VR) for focal distances greater than the fault width. If we truncate this function at T = w/2$, as the calculated rise times shown in Figure 2.7 suggest, we get the slip function s(T) = C A/(-
+ 2w/VR)
[H(T)
H(T -
w/25)].
This expression, in turn, provides a static slip which is 2
within five percent of the static value dow/p$
suggested by Figure 2.6 for the fault center-line, at focal distance greater than the fault width. More generally, we can define a rise time TR by T
= t y
+ x7 R
s V
R where 24
1~~
2~1w~
R 2
+
lxi) t = min (2.4) s 1
1 2 x2 1
V ~
x2
+
lyi R2 Referring to Figure 2.14, ts is a rupture arrest time given by the smaller of (1) the sum of the rupture travel time over path A plus the shear wave travel time over path A', or (2) the sum of the rupture travel time over path B plus the shear wave travel time over path B'.
Then the expression s(T) = C a
L-T4
+ 2 min (r,w)/VRj
[H()
H(T -
TR (2.5) closely approximates the rise time, static slip, and peak veloc ity of our finite difference solutions, throughout the fault surface. Equation (2.5) provides a means of scaling our re sults to long faults and high frequencies.
2.4.3 Model Plausibility The model we have treated appears to incorporate much of the physics appropriate to earthquake faulting. Rupture initiates at a point and spreads at finite rupture velocity; slip is driven by the tectonic prestress and the dynamic stress concentration developed at the crack edge; slip ceases when fault-plane tractions fall below the kinetic friction level. On the other hand, the physical problem has been highly idealized in these calculations.
In particular, rup ture growth-is coherent, stopping of rupture growth is abrupt, stress-drop is uniform, and crack-tip stress concentra tions are unconstrained.
The rupture process in actual 25
(000)
Hypocenter Figure 2.14.
Geometry for estimating stopping phase arrival time.
26
earthquakes is certainly substantially more complex than our model permits, although the details are poorly understood.
The question we address here is whether we have included suf ficient physics to obtain plausible predictions of ground motion within the frequency band of engineering interest. It ought to be noted that other investigators (Del Mar Technical Associates, 1978) have asserted that rupture incoherence and crack-tip yielding are the essential source phenomena con trolling ground motion generation over most of that frequency range.
In that view, high-frequency generation would be in dependent of stress-drop and source size. This is in contrast to the present idealized but rather rigorous treatment, which predicts that high-frequency generation scales with do and /.
As a test of the plausibility of the model, we apply Equation (2.5) to compute ground acceleration near a buried strike-slip fault. Figure 2.15 shows the fault geometry and receiver location. The fault is 15 km long, 5 km wide, buried between the depths of 2.5 and 7.5 km. The receiver is located at the surface, in the plane of the fault, so that the only non-zero component of motion is in a direction normal to the fault plane, and is predominantly shear waves. A two-dimen sional numerical integration over the fault surface was per formed, using the integral solution of Haskell (1969) with the slip function of Equation (2.5).
The integration mesh was.sufficient to retain accuracy for frequencies up to 10 Hz.
The calculation was for a uniform wholespace (with the free surface effect approximated by doubling the computed motion).
Accelerograph Station 2 is similarly situated with re spect to the 1966 Parkfield earthquake (although that event apparently involved faulting over a length of 30 km or so).
Therefore, for comparison, we scale our theoretical source to the stress-drop of 25 bars estimated by Archuleta and Day (1979) for the Parkfield event. Figure 2.16 shows the computed accelerations for three choices of rupture velocity.
27
Receiver 71-11
/ 1---I7 77 77 2.5 km Ilypocenter
- 7. 5 km 2
km 15 km Figure 2.15.
Problem geometry for accelerogram synthesis. The observer lies in the plane of the fault.
V
=.753
-. 2.
-. 3 3
2 -V,
=.8 O
0
-. 2.
-. 4 1.0.
.5 0
-. 5
-1.0 TIME (sec)
Figure 2.16.
Synthetic accelerograms computed with Equation 2.5, for three different ructure velocities VR 29
Peak acceleration increases with rupture velocity, ranging from 0.24 g to 0.90 g. This range brackets the peak acceler ation value of 0.5 g recorded at Station 2 for the Parkfield event. Apparently, then, the simplicity of the rupture model has not led to unacceptable near-field accelerations.
In particular, the coherent rupture idealization did not result in the enormous (greater than 6 g) acceleration peak reported by Del Mar Technical Associates (1979) for a coherent rupture-model. The large, but brief, acceleration peaks reported in their study may be attributable instead to the method of integration over the fault plane.
in that in tegration scheme, the function evaluations (slip function and propagation Green's function) were made at mesh points whose spacing was many times the smallest wavelength of interest, and time shifts of these functions, corresponding to smooth rupture over each mesh subdivision, were approximated.
Small inconsistencies in the smooth-rupture approximation between adjacent mesh subdivisions may have led to the observed acceler ation peaks.
Although coherent rupture propagation did not lead to unacceptable near-field peak motion, the abrupt stopping of rupture propagation, which was artificially enforced at the fault edge, was responsible for the predominant acceleration peak in each accelerogram in Figure 2.16.
The peak acceler ation in each case is clearly attributable to the first arriving stopping phase from the upper edge of the rectangular fault. The preeminence of the stopping phases may indicate.
a significant weakness of the uniform rupture velocity/stress drop model.
30
III.
SPONTANEOUS RUPTURE MODELING
3.1 INTRODUCTION
In Section II, we analyzed an earthquake model in which rupture velocity and stress-drop were uniform over the fault plane, stopping of rupture was abrupt, and large stress concentrations were permitted to develop at the rupture front.
In this section, we examine the effects of relaxing these con straints.. Rather than specifying a rupture velocity, we specify a failure criterion and permit rupture to ensue spon taneously. The failure criterion is such that shear stress in the fault plane is limited by a frictional strength, and a fixed amount of energy is dissipated per unit area encompassed by extension of the rupture.
3.2 FAULT CONSTITUTIVE MODEL A description of the fault-plane constitutive model is given in Appendix A, and we briefly overview it here.
In the displacement-weakening model, shear traction on the fault plane, which is initially aT, is limited by a frictional strength a.
When the shear traction reaches c, further stress buildup is relieved by a fault slip. The slip velocity at a point is opposed by a frictional traction which is a.
prescribed, decreasing function of the total slip which has occurred at the point. When the total slip reaches a constant, d
do, cohesion is considered destroyed, and the frictional traction is held at a fixed "dynamic friction" level a.
3.3 DESCRIPTION
OF THE CALCULATIONS Six three-dimensional finite difference earthouake simulations were analysed for this portion of the study.
These are listed in Table 3.1.
In each case, the maximum frequency accurately propagated by the finite difference grid 31
TABLE 3.1 FAULT PARAMETERS FOR SPONTANEOUS RUPTURE COMPUTATIONS PROBLEM d
(m) a (bars) a (bars) a (bars)
- 9. (km) w (km) f (Hz)
DESIGNATION o
o T
f max I
0.1 1020 1000 900 m
0 15 II 0.08 1050 1000 900 15 III 0.08 1050 1000 900 16 4
5 IV 0.1 1020 1000/900 900 c
m 5
V 0.1 1020 1000/925 900 m0 m
5 VI 0.08 1020 1000/925 900 18 6
5
is given by f max The parameters Z and w describe the rectan gular region within which rupture is permitted --
outside the Z by w rectangle, the strength increases abruptly, preventing further extension. Note that in four cases, no a priori limits are placed on rupture size; faulting in these cases either advances without bound or is retarded by heterogeneities in the prestress configuration. In all cases, the hypocenter is at the coordinate origin, the x direction is parallel to the direction of shear prestress, and symmetry is assumed about both the x and y axes. The density and P and S veloc ities, p, a, and 6, are as given in Section 2.1.
Although slip was permitted to have a component in the y direction (perpendicular to the prestress), this component was invariably far smaller than the x component, and will be neglected here.
In the first two calculations, the frictional strength and stress-drop are everywhere uniform, so that once initiated, the rupture does not stop growing. In these cases, we utilize small zoning (50 m cubes) and focus on the high frequency ( 15 Hz) character of the slip over a limited region surrounding the hypocenter.
In the third calculation, frictional strength and stress-drop are uniform over a 4 km wide and 16 km long strip, outside of which rupture is prevented.
Zone size is 130 m.
The final three calculations involve nonuniform dis tributions of prestress aT (with a0 and a. uniform).
For these cases, Table 3.1 gives the maximum and minimum of a on the fault plane. The prestress configurations for Problems IV, V, and VI are diagrammed in Figure 3.1.
Hatching delineates the regions of elevated prestress.
In Problems IV and V, rupture propagation stops spontaneously; in Problem VI, rupture reaches the strength barriers, with Z = 18 km and w = 6 km.
In all three cases, the zone size is 150 M.
33
1050 15 PROBLEM IV 1000 950Of
-T 900 1050 o
PROBLEM V S1000 U) 950 T
f 900
.p.
y 1050 PROBLEM VI 1000 950 x
900 I
I I
I I
0 2
4 6
8 x (km)
Figure 3.1.
Nonuniform stress configurations for Problems IV, V, and VI.
3.4 NUMERICAL RESULTS FOR UNIFORM STRESS Figure 3.2 shows contours of the rupture front at 0.1 second time intervals, for Problem 1. The rupture front is defined to enclose regions of the fault plane in which the shear traction has relaxed to a f. As described in Appendix B, a minimum rupture velocity of 0.56 has been enforced to initiate rupture, and the quarter-circular contour at 0.2 seconds reflects this.
Subsequent contours represent spontaneous rup ture propagation, however. Propagation is most rapid in the x direction, least rapid in the y direction, leading to a roughly elliptical rupture front.
Figure 3.3 shows the rupture velocities observed in the two principal directions, as a function of hypocentral dis tance. The flat part of the curves represents the minimum rupture velocity. As the minimum is exceeded, rapid accelera tion of the rupture front occurs.
In the y direction, this acceleration gradually decreases, and the rupture velocity smoothly approaches the shear wave velocity; at 2 km distance the rupture velocity has reached 0.95S.
In the x direction, the rupture reaches the shear wave velocity at a distance of just over 0.5 km. The rupture velocity then levels off some what, before accelerating again toward the P wave velocity, reaching 0.9a at a distance of 1.3 km.
Slip velocities along the x and y axes are shown in Figures 3.4 and 3.5, respectively. These have been low-pass filtered, with a 15 Hz cutoff frequency. As in the fixed rupture velocity study, we find that peak slip velocity in creases with focal distance, in the absence of edge effects.
The increase appears to be similar in form to the /V increase found for circular, fixed rupture velocity models.
Peak slip velocity, at a given focal distance, is higher on the y axis than on the x axis.
This azimuthal variation is qualitatively 35
0 t=0O. 2 0.4 0.6 0.8 2
0 1
2 x (kmn)
Figure 3.2.
Rupture front contours for Problem I, at 0.1 second intervals.
36
6 direction 5
0 U
3 V direction G
a 2
E eH 0
0 1
2 3
HYPOCENTRAL DISTANCE (kn)
Figure 3.3.
Rupture velocity, in x and y directions, for Problem I.
37
4 2.25 3
1.5' 0.7 I
~
2
=
0.3
- 1.
I 0
I I'
0
.2
.4
.6
.8 1.0 TIME (sec)
Figure 3.4.
Slip velocities along x axis, for Problem I. The time histories have been low-passed filtered, with a 15 Hz cutoff.
38
2.20 6
1.5 U
0.75 C.,*,Ii 0.3, 1
1..
I 2
S I*
2 0ma 0
.2
.4
.6
.8 1.0 TIME (sec)
Figure 3.5.
slip velocities along the y
- axis, for Problem I.
Time histories have been low-pass filtered with a
15 Hz cutoff g
.9 I
39
explained by the Burridge and Willis (1969) solution for the self-similar, expanding elliptical crack:
2 2
2 2
s(t) = C' t 2 x
Hi xt (3.1)
TS2 V2 V2 V2 x
y x
v where C' is a constant and Vx and V are the rupture velocities x
y in the x and y directions, respectively. This expression pre dicts that peak slip velocity at a given focal distance alcng the y axis will approach
/V times the peak slip velocity x y at the same focal distance along the x axis. At r = 2.2 km, this accounts for about half of the observed difference between the two azimuths.
The discrepancy is not surprising, consider ing that Equation (3.1) strictly applies only for constant rup ture velocities, with Vx less than the Rayleigh wave velocity and V less than a.
Problem II differs from Problem I only in the frictional strength ao and weakening parameter d.
The frictional strength for Problem II is al + 1.5 (a -
a-), compared to a, T
r
+ 1.2 (aT -
- 04) in Problem I; d is reduced for Problem II sufficiently to keep the specific fracture energy the same for both problems.
Figure 3.6 shows contours of the rupture front at 0.1 second time intervals, for Problem II.
Again, the rupture fronts are approximately elliptical, with propagation most rapid in the x direction.. As Figure 3.7 indicates, the rupture velocities in the x and y directions accelerate somewhat less rapidly than in the case of Problem I. In the x direction the leveling off of rupture velocity at about 3 km/sec is more distinct than in Problem I, and the acceleration toward the P velocity occurs later.
40
e0 t=0.2 0.46 1
- ~
0.8 2
-1.0 0
1 2
x (km)
Figure 3.6.
Rupture front contours for Problem II, at 0.1 second intervals.
41
5 direction S4 C-4 c..
3direction 2
2 0
12 3
HYPOCENTPRAL DISTANCE (kmn)
Figure 3.7.
Rupture velocity, in x and y directions, for Problem II.
42
Figures 3.8 and 3.9 show the low-passed (15 Hz) slip velocities, along the x and y axes, respectively. These are very similar to those for Problem I shown in Figures 3.4 and 3.5, but generally a few percent higher. An exception is the velocity at x = 2.25 km, which is a few percent lower than for Problem I.
In-Problem III, we revert to a coarser grid (150 m cubes) and restrict the fault size to 4 km to 16 km.
The fault parameters, Table 3.1, are otherwise as for Problem II.
Rupture front contours are shown in Figure 3.10.
Initial rupture fronts are again elliptical, and the rupture velocity accelerates to near the P velocity. Rupture acceleration is rapid out to about 3 km in the x direction, where it reaches about 4.5 km/sec, then accelerates more slowly, reaching about 5.5 km/sec at 8 km focal distance. Low-passed (5 Hz) slip velocities along the x axis are shown in Figure 3.11.
Over this distance range, the peak slip velocity has not yet settled at a constant value, as occurred for the same fault geometry at fixed rupture velocity (Figure 2.11).
This is apparently due to the fact that the rupture velocity is still increasing in Problem III, even at 8 km focal distance. Figure 3.12 demonstrates that the peak slip velocity is nearly constant across the fault width, as for the fixed rupture velocity case (see Figure 2.13).
3.5 NUMERICAL RESULTS FOR NONUNIFORM STRESS For Problem IV, the prestress exceeds the dynamic friction a-by 100 bars, over a circular region 1.43 km in radius, and is equal to a outside that region (see Figure 3.1).
The frictional strength a is uniform, 120 bars above a f* Figure 3.13 shows contours of the rupture front at 0.3 second intervals.
The rupture accelerates rapidly over the prestressed patch, then abruptly decelerates as it breaks 43
4 I
2.25 1.5**
0a 3
-Il Il E
30.7 1
cnI
\\1*
It 2
I I0.
- a 2
1Ito' r =0.3 I
1 I
0
.2
.4
.6
.8 1.0 TIME (sec)
Figure 3.8. Slip velocities along the x~ axis f or Problem II.
Time histories have been low-sass filtered with a 15 Hz cutoff.
g4
2.20 6
1.5 CIIL r =0.3 2
0 4.2
.4
.6
.8 1.0 TIME (sec)
Figure 3.9.
Slip velocities along the y axis, f or Problem II.
Time histories have been low-pass filtered with a 15 Hz cutoff.
45
0 0.9 1.2
- 1. 5 1..2.1 3
0 3
6 9
x (km)
Figure 3.10. Rupture front contours for Problem III, at 0.3 second intervals.
46
0 44 7 05 3.0 2I Cn
- 0.
4.95 T
IM 3
4.05 3.0I47
- 1 N.
- t.
I 4
f m
I'1 0
~
0 1
2 5
Hz cutoff 47*
4 S.
o ~
- 1.
- \\
.0.45 1.1 0*.
3 1.8 2
3 TI.ME (sec)
Figure 3.12.
Slip velocities across the width of the fault at x
=
6 km, for Problem III.
Time histories have been low-pass fLtered with a 5 Hz cutoff.
48
0 0.6 1.5 0.9 1
- 2.
3 4
0 1
2 3
4 x (kn)
Figure 3.13.
Rupture front contours for Problem IV, at 0.3 second intervals.
0 49
into the zero-stress-drop region.
In the y direction, deceler ation is very abrupt, and the rupture penetrates only about 150 meters into the low-stress region.
In the x direction, however, the rupture penetrates about 500 meters into the low-stress region. By 1.5 seconds, rupture growth has ceased.
Figure 3.14 shows both peak slip velocity (low-passed, 5 Hz) and rupture velocity along the x axis. We note that the two quantities roughly parallel each other. When the rupture decelerates, the peak slip velocity falls off. A maximum rup ture velocity of 4 km/sec is reached at the edge of the prestressed zone.
Problem V differs from IV only in having a prestress 25 bars above a outside the circular high-stress patch. In spite of this non-zero stress-drop, we see from Figure 3.15 that the rupture still stops spontaneously, this time at about 2.1 seconds.
Rupture again decelerates outside the high stress patch, but overshoots it by 0.6 km in the y direction and 1.5 km in the x direction.
Figure 3.16 again demonstrates the strong linkage between peak slip velocity and rupture velocity.
Rupture velocity in the x direction peaks at about 4 km/sec, drops sharply at the edge of the circular patch, then recovers to about 1.7 km/sec and smoothly decelerates to zero.
The slip velocity mirrors this behavior.
Problem VI contains several patches of 100 bars ef fective stress, surrounded by a 25 bar effective stress region, with a and a f uniform (see Figure 3.1).
Figure 3.17 shows a
r the rupture front contours at 0.1 second intervals. A fairly complex pattern of rupture develops. Along the y axis, for example, rupture stops just after 1 second. As rupture advances on the other parts of the fault plane, however, the stress intensity along the y axis increases, and rupture recommences at about 1.8 seconds. Along the x axis, rupture 50
PEAX SL.-? TTEOC:TYZ 0
3 -
u TtRz VELocI 2
0 Figure 3.14.
Peak slip velocity (low-passed, 5 Ez cutoff) and rupture velocity along the x axis, for Problem IV.
51
t0-O.3 1.5) 0.6 2.1 0.9 1 1
2 3
4 0
1 2
3 4
x (km)
Figure 3.15.
Rupture front contours for Problem V, at 0.3 second intervals.
52
405
?.STRE 1
0 0
950 --
00 DISTA-NCz
('c-0 Figure 3.16.
Peak slim velocity (low-passed, 5 Fz cutoff) and rupture velocity along the x axis, for Problem V.
53
t 1.0 1.9 0.5 1
5 2.0 2.5 t = 0.5 2.0 2.5 ---
0 3
6 9
x (km)
Figure 3.17.
Rupture. front-contours for Problem VI, at 0.1 second intervals.
54
accelerates rapidly as it breaks each high-stress patch, and decelerates between patches. At 1.1 seconds, and then again at 1.9 seconds, the rupture front "jumps," leaving unbroken areas behind, which subsequently break.
Figure 3.18 shows slip velocities and rupture velocities along the x axis. The close relationship between slip velocity and rupture velocity is especially visible here. The dashed portions of the rupture velocity curve represent regions that ruptured out of sequence as the rupture front jumped ahead to a high-stress patch.
We note that local rupture velocities in excess of the P wave velocity occur at the edges of the stress concentrations, although average rupture velocity is, of course, less than the P wave velocity. Very abrupt drops in rupture velocity occur whenever the crack propagates into a low-stress region.
The drop in rupture velocity causes the stress intensity to build up, and the rupture velocity then rebounds somewhat. As'the rupture front propagates further from the high-stress patch, we expect the stress intensity, and therefore the rupture velocity, to fall (see, for example, Kostrov, 1966; Husseini, et al., 1975), and this is what is observed in Figure 3.18.
As rupture proceeds coherently across a stress concen tration, peak slip velocity increases smoothly, as it did in the constant stress models.
Drops in rupture velocity, however, are accompanied by drops in peak slip velocity. When the rupture breaks into a new stress concentration, the peak slip velocity behaves approximately as though the new stress concentration represented an independent event.
Figure 3.19 shows slip velocity time histories observed along the x axis.
The late time level for all the points shown approaches approximately 80 percent of the average stress-drop (40 bars) divided by the impedance (;.).
PEAK SLIP VELOCITY 3
2
- 1.
0 8
RUPTURE VELOCITY 6
2 0
1050 STRESS o
1000 950 900 t 0
1 2
3 4
5 6
7 8
9 DISTANCE (km)
Fiaure 3.18.
Peak slip velocity (low-passed, H
Hz cutoff) and rupture velocity along the x axis, for Problem VI.
Dashed curves for rupture velocity indicate regions which ructured out of seauence. For example, rupture occurred at x = 5.6 km while the region between 3.0 and 5.6 km was still intact.
56
4.0 3.0 7.05 4.05 II' En U1O 3.0 05 r =1.9 2.0 -3.0 1
60 u
I'I 1?0 0
1 2
3 TIME (sec)
Figure 3.19.
Slip velocities along the x axis, for Problem VI.
Time histories have been low-pass filtered with a 5 1Iz cutoff.
3.6 DISCUSSION In this section, we examine the imvlications of these numerical results for modeling strong ground motion. First we assess the effect of finite frictional strength on slip velocity. Then we discuss the effect of a complex stress dis tribution, and the relationship between rupture velocity and peak slip velocity.
We find, from the results of Problems I, II and III, that the frictional strength has a minimal effect on peak slip velocity over the frequency range considered. There is no support in these results for the conjecture that peak slip velocity is bounded by (a -
7,)/p6. Taken as a whole, the slip velocities for Problems I, II and III shown in Figures 3.4, 3.5, 3.8, 3.9, 3.11 and 3.12 support the approximation suggested in Section II:
peak slip velocity is approximately 2wf/3 Ac/pa, where f is the high frequency cutoff, and w is the fault width (focal distance r would replace w in the absence of edge effects, as in Problems I and II).
This expression gives a reasonable estimate for both a 5 Hz cutoff frequency (Problem III) and a 15 Hz cutoff frequency (Problems I and II).
We have seen that the presence of stress concentrations leads to irregular rupture propagation. For the relatively simple prestress configuration of Problem VI, the rupture veloc ity fluctuates, goes to zero and restarts, and even locally exceeds the P wave velocity of the medium; the rupture front sometimes jumps ahead, leaving unbroken patches. One of the features of the DELTA randomness model (Del Mar Technical Associates, 1979) is that random fluctuations introduced in the rupture arrival times imply that rupture velocity sometimes locally exceeds a, while the average rupture velocity is less than S. Our numerical results indicate that this feature is physically acceptable, and the phenomenon can occur in the presence of even fairly modest stress heterogeneities.
58
The DELTA randomness model accompanies the irregular rupture advance with a spatially invariant slip function. In contrast, the current numerical results indicate that the slip function is strongly linked to the rupture velocity. In par ticular, the peak slip velocity varies roughly in parallel with the rupture velocity. Essentially, each isolated stress con centration on the fault plane acts as an independent source.
As an individual stress concentration ruptures, peak slip velocity starts at a low level, around Aa/pS, increases approx imately as V/ over regions of smooth, coherent rupture propa gation, and reverts to a low value when the stress concentration is fully encompassed and rupture slows down.
To estimate the effect of artificially coupling an irregular rupture to a uniform slip function, we use the expression derived by Richards (1973) for the first-motion approximation to the shear-wave acceleration U due to an expanding circular crack:
C dea V2 R~
us 2
2 2
2 H (t -
R/3)
(3.2)
Rp$
(1 -
V /$
sin 8) where R contains the double couple radiation pattern, e is the angle from the fault normal to the observer, and R is the hypocentral distance. We will also use the corresponding expression for an expanding, uniform circular dislocation:
2 VR Us 2
2 3/2 s (t -
R/3),
(3.3)
$R (1 V2/8 sin 2) where s(t) is the slip velocity. Equation (3.3) can be deduced from Savage (1966).
59
We have noted that, in a complex rupture, each individual stress concentration within the ccmoosite event behaves essen tially as a separate source. We will consider an observer in the plane of a fault (e = 900 ), situated on a radiation anti node, and at a hypocentral distance of 10 km from a particular isolated high-stress patch (da = 100 bars) within a complex rupture. Now we compare two estimates of the initial shear wave acceleration:
(1) that obtained if the high-stress patch is assumed to be a separate event governed by Equation 3.2, with a 100 bar stress-drop, and (2) that obtained (using Equa tion 3.3) if the stress-concentration is represented as an isolated dislocation which has jumped ahead of the rupture front, on which initial slip velocity is everywhere fixed at 8 m/sec.
The first method of estimation corresponds physically to the behavior of our dynamic model; the second corresponds to the DELTA randomness model in which rupture velocity is irregular, but peak slip velocity is uniform,.with a prescribed value of 8 m/sec, throughout the fault plane.
If we assume in both cases that the patch ruptures at a velocity of 0.98, an initial acceleration of 0.67 g would be predicted for the observer by the first method, and an initial acceleration of 2.7 g would be predicted by the second method. Thus, an anomalously high acceleration is generated by artificially imposing a uniform slip function in the presence of irregular rupture propagation.
60
IV.
SUMMARY
Numerical solutions for a simple, constant rupture velocity, constant stress-drop earthquake model show a significant influence of fault width on the slip function.
For points on the fault which are more than one fault width away from the focus, the slip function is spatially quite uniform. -The rise time is given by w/28, the static slip by 2
dow/P82, and the peak slip velocity by V/27f Aa/p8, where w is the fault width, da is the stress-drop, f the maximum frequency, a the shear speed, and p the density. The slip function can be closely approximated by a simple, closed-form expression (Equation 2.5)..
The simple physical model, applied to compute ground acceleration near a buried strike-slip fault, gives peak acceleration comparable to that observed for the Parkfield event.
A displacement-weakening constitutive model has been ap plied to simulate spontaneous rupture in the presence of both uniform and heterogeneous prestress. Rupture velocities in excess of S persistently occur in the direction of prestress, while rupture velocities in the direction normal to the pre stress are less than 8.
For uniform prestress the above expression for peak slip is still valid for the spontaneous rupture model.
For nonuniform prestress, each isolated stress concen tration behaves roughly as an independent event. Rupture velocity becomes quite irregular, and locally can even exceed the P wave velocity a.
Peak slip velocity is strongly coupled to the rupture velocity, and there is considerable risk in ignoring this coupling when formulating kinematic fault models.
In particular, a uniform dislocation fault model, combined with an irregular rupture velocity, generates large acceleration peaks which might not be justified by the fault dynamics.
S6
V. REFERENCES Andrews, D. J. (1976a), "Rupture Propagation with Finite Stress in Antiplane Strain," JGR, 81, pp. 3575.
Andrews, D. J. (1976b), "Rupture Velocity of Plane Strain Shear Cracks," JGR, 81, pp. 5679.
Apsel, R. J. (1979), "Dynamic Green's Functions for Layered Media and Applications to Boundary-Value Problems,"
Ph.D. Dissertation, University of California, San Diego.
Archuleta, R. J. and S. M. Day (1980), "Dynamic Rupture in a Layered Medium:
1966 Parkfield Earthquake," submitted to BSSA.
Burridge, R. and J. R. Willis (1969), "The Self-Similar Problem of the Expanding Elliptical Crack in an Anisotropic Solid," Proc. Cambridge Phil. Soc., 66, pp. 443-468.
Burridge, R., G. Conn and L. B. Fretfnd (1979), "The Stability of a Rapid Mode II Shear Crack with Finite Cohesive Traction," JGR, 85, pp. 2210-2222.
Cherry, J. T., E. J. Halda and K. G. Hamilton (1976), "A Deterministic Approach to the Prediction of Free Field Ground Motion and Response Spectra from Stick-Slip Earthquakes," Earthquake Engineering and Structural Dynamics, Vol. 4, pp. 315-332.
Cherry, J. T. (1977), "Users' Manual for the TRES Code,"
Systems, Science and Software Report SSS-R-77-3128, January.
Dahlen, F. A. (1974), "On the Ratio of P-Wave to S-Wave Corner Frequencies for Shallow Earthquake Sources," BSSA, 64, pp. 1159-1180.
Day, S. M., T. C. Bache, T. G. Barker and J. T. Cherry (1978),
"A Source Model for the 1975 Pocatello Valley Earth quake," Systems, Science and Software Report AFGL-TR 79-0001 submitted to the Air Force Geophysics Laboratory, December.
DELTA (Del Mar Technical Associates) (1978), "Simulation of Earthquake Ground Motions for San Onofre Nuclear Generating Station Unit 1," Final Report for Southern California Edison Company, submitted for review to the Nuclear Regulatory Commission.
62
DELTA (Del Mar Technical Associates) (1979), "Simulation of Earthquake Ground Motions for San Onofre Nuclear Generating Station Unit 1,.Supplement 1," Report for Southern California Edison Company. Submitted for review to the Nuclear Regulatory Commission.
Haskell, N. A. (1969), "Elastic Displacements in the Near Field of a Propagating Fault," BSSA, 59, pp. 865-908.
Husseini, M. I., D. B. Jovanovich, M. J. Randall and L. B.
Freund (1975), "The Fracture Energy of Earthquakes,"
Geophys. J. R. Astr. Soc., Vol. 43, pp. 367-385.
Ida, Y. (1972), "Cohesive Force Across the Tip of.a Longitudinal Shear Crack and Griffith's Specific Surface Energy,"
JGR, 77, pp. 3796-3805.
Knopoff, L. (1958), "Energy Release in Earthquakes," Geophys.
J., 1, pp. 44-52.
Kostrov, B. V. (1964), "Self-Similar Problems of Propagation of Shear Cracks," J. Appl. Math. Mech., Vol. 28, pp.
1077-1087.
Kostrov, B. V. (1966), "Unsteady Propagation of Longitudinal Shear Cracks," J. Appl. Math. Mech., 30, pp. 1241-1248.
Madariaga, R. (1976), "Dynamics.of an Expanding Circular Fault,"
BSSA, 66, pp. 639-666.
Neuber, H. (1937), Kerbspannungslehre, Springer-Verlag, Berlin; Theory of Notch Stresses: 'Principles of Exact Calcula tion of Strength with Reference to Structural Form and Material, The Office of Tech. Info., AEC-TR-4547, 1958, (English Translation).
Richards, P. (1973), "The Dynamic Field at a Growing Plane Elliptical Shear Crack," Int. J. Solids and Structures, 9, pp. 843-861.
Savage, J. C. (1966), "Radiation from a Realistic Model of Faulting," BSSA, 56, pp. 577-592.
Swanger, H. J. and D. M. Boore (1978), "Importance of Surface Waves in Strong Ground Motion in the Period Range of 1 to 10 Seconds," Proceedings of the Second International Conference on Microzonation for Safer Construction Research and Application, San Francisco, California, November 26 -
December 1, 1978, pp. 1447-1457.
63
Wiggins, R. A., G. A. Frazier, J.
Sweet and R. Apsel (1978),
"Modeling Strong Motions from Major Earthquakes,"
Proceedings of the Second International Conference on Microzonation for Safer Construction -
Research and Application, San Francisco, California, November 26 December 1, 1978, pp. 693-700.
0 64
APPENDIX A PROBLEM FORMULATION In the fault simulations reported here, we treat fault ing as a propagating stress relaxation due to shear failure on a planar surface. We specify (1) the initial state of the medium and its constitutive properties, and (2) the boundary conditions to be satisfied on the plane in which faulting is permitted. For the fault simulations at constant rupture velocity,Section II, we prescribe the evolution of the fault surface. For the simulations of spontaneous rupture,Section III, the boundary condition on the fault plane embodies a constitutive law permitting shear failure. The subsequent growth of the fault surface is not prescribed, but results from the constitutive law and the fault dynamics.
In either case, the mathematical model of faulting is then solved using a three-dimensional finite difference method (Cherry, 1977).
A.1 INITIAL CONDITIONS AND CONSTITUTIVE PROPERTIES OF THE MEDIUM For time t less than zero we assume that an equilibrium state of stress exists, with displacement and velocity every where zero.
The equilibrium configuration is such that the prospective fault plane experiences a shear traction a and a compressional normal traction N which may vary with position on the fault plane.
In this study, we assume that material behavior off the fault plane is linearly elastic. Since average stress changes associated with faulting are modest, on the order of a few hundred bars, and since faulting represents a relaxa tion of stress, except near the fault edges, linear elasticity is probably an adequate model of material behavior away from the immediate zone of faulting.
65
A.2 BOUNDARY CONDITIONS The rupture surface is assumed to occupy a prescribed plane, with unit normal vector n. Across this plane, we per mit a tangential displacement discontinuity s(x,t), and re quire continuity of the normal displacement:
lim^
^
s(xt)
(u(x + en) -
u(x en), n
- s
= 0, with e > 0, where u is the displacement vector. We also re quire continuity of traction across the fault plane. T de notes the tangential traction vector, defined in terms of the stress tensor a by
-r n
- a I -nn),
where I is the identity tensor. We also define T as the tangential traction vector on the fault plane which would be required to enforce continuity of velocity.
A.2.1 Fixed Rupture Velocity Case For the part of the study reported in Section II, we specify the geometry of the fault surface as a function of time, rather than deriving the fault surface from the model as a consequence of a failure mechanism. The rupture surface is assumed to nucleate from a point, and to expand symmetrically at a constant, prescribed rupture velocity (V R), until it reaches prescribed boundaries.
Z(t) denotes the fault surface at time t, w and Z the width and length of a rectangular fault, and x and y Cartesian coordinates in the fault plane.
Z(t) consists of all points x,y such that 2
2 2
2 x + y <VRt and 66
1x <2 02 and Iyl < w 2
Two significant features of this model of rupture are (1) the rupture velocity accelerates instantaneously to its final velocity, and (2) the rupture advance terminates ab ruptly, i.e., the rupture velocity decelerates instantaneously to zero at a prescribed boundary.
Non-zero slip is permitted only within Z(t):
s(x,t) = 0 for x not on '(t)
Within Z, we define af to be a sliding frictional traction which opposes the slip velocity s and is proportional to the normal stress:
S C7fj~
f or
>s
>0
-f 0
for 0
A A
where a is -
n*
- n, the product of the normal stress and the coefficient of dynamic friction (a is presumed to be positive).
The boundary condition of Z is I
if I I > a
-f
-f 0c c -f 67
A.2.2 Spontaneous Rupture Case For the simulations reported on in Section III, we do not specify E(t),
but permit a displacement discontinuity s(x,t) anywhere in the prescribed fault plane. The tangential traction vector is governed by a displacement-weakening failure model. In this model, a finite frictional strength a0 (which may be a function of position) is assigned to the fault plane prior to initiation of sliding. Slip commences at a point when necessary to prevent a stress concentration in excess of 00 from occurring. This relative displacement, with ampli tude denoted by s, is assumed to weaken the fault plane at that point according to some prescribed weakening function f(s), until the total slip equals the value do, at which state the fault has lost all cohesion and reaches the kinetic friction level, a.
Displacement-weakening mechanisms of this form have been investigated, in two dimensions, by Ida (1972) and Burridge, et al. (1979) using analytical methods and by Andrews (1976a, b) using a finite difference method. Cherry, et al.
(1976) generalized the mechanism to include material yielding near the fault plane.
Ida's (1972) analysis indicates that as long as f(s) is assumed to be fairly smooth, the behavior of the slip function is insensitive to the precise shape of f.
We have simply assumed f to be linear. The boundary condition on the fault plane then takes the form S
-f(s)--
if I
> f(s)
T cif
_.1 I < f (s) where 68
1 -
a
+
if 0 < s d 0 0
d 0 f-0 f(s)
=
if s>dO 69
APPENDIX B INITIATION OF RUPTURE Appendix A outlines a rupture model based on a displace ment-weakening friction law.
Once started, the rupture process proceeds spontaneously, governed by the displacement weakening model in conjunction with the fault dynamics.
However, some additional mechanism is required to initiate rup ture from the equilibrium prestress configuration. We might imagine, for example, that a relatively small region of the prospective fault plane has been weakened, and that the shear stress there falls from aT to a.
This initial crack then slides stably under a slowly increasing tectonic load.
Eventually, a situation develops analogous to the critical crack in elastic fracture mechanics, and accelerating crack growth ensues.
Quasi-static modeling of the process leading to the un stable crack growth is beyond the scope of this study.
In our fault modeling, a mechanism of this sort is approximated by enforcing a minimum rupture velocity of half the shear wave velocity within a circular region of radius a (Andrews,.1976b, used a similar method to initiation plane-strain shear cracks).
The choice of this value for minimum rupture velocity is a compromise between approximating quasi-static crack inception (favored by rupture velocity approaching.zero) and reducing computation time (favored by a high minimum rupture velocity).
The assumption of a circular incipient crack is arbitrary.
For the purpose of estimating a suitable value of the "critical" radius a, we assume that, at the onset of spon taneous growth, the circular crack initially expands uniformly, retaining circular shape. We then seek a balance between strain energy release rate and the energy dissipation rate at the crack edge, per unit increase of the crack radius.
We 70
start with Neuber's (1937) solution for the static slip on a circular shear crack in a Poisson solid:
s(r) 24 T
f r
s~r
=7w 2
a 1-p7 a
where s is the static slip, p the density, 8 the shear wave speed, a the radius, and r the distance from the crack center.
The change in strain energy U in creating the crack is 2
- 2) 8a3
_(r f
2
'a 7p82 and the work done against friction W is 3
W = 2a (a
8a) f
.T f 7 o2 so that total "available" energy, E, is
)2 8
T f
3 E =- (U +W) a2.
7 Pa2 The displacement-weakening mechanism dissipates energy Z at the rupture front in the amount
-=
Tr(a a ) d a
da 0
f 0
per unit increase in crack radius.
The desired value of a is that corresponding to a stationary value of E -
Z:
2 7-, pa2 (aO
- a f) a 7=
28(
0 f
dO (B.1) 24 (a f)
For simulations I through V, discussed in Section III, the critical radius, a, as defined by Equation (B.1),
was 71
approximately 350 n; for simulation VI, a was approximately 300 m.
72