1.

W711 1W! "1ng '1'“ 61,11
1H‘NI111WI11 II 11911
”I “I1 “11$ 1:. .: '..1|1.:1"' III.

7h. WE», <
{IN-ll . fidpf-

  
   
   
 

  
    
 

    
 
 

M
111 1% 13-4,

. I' "i n
t" '1‘ ' lfig': : :3 n . A '
‘.'.l ' . ' .

1‘1} 0 3 3? b! '0 ‘g ‘. ' . J a 90
4’5 i. ". '_ .

I .1 I l‘ .1 .p; .l :‘ V “d .‘ 1.. ’r. I. L . V N
‘I‘I “:‘h‘ ' ‘ I md ( £11411“ «Moutggg‘q‘ r : ? Liar 1

‘ + ' .. ' . v I I . 0:. " 19:8). “ ..' I
‘ | 1%.; 1 '1“ "1 “lirzftliz l 9| . :I I'M in: 15-133.... C" . w ‘
2m. M I ~ . . -
1H .11'1111 [1‘

. I'r '. 23:1": :23“:
«1’ 491.«I.11Inh‘31§'u4 I; “391% |..

      
 

   
 
 
  
  

 
   
 

   
     

11 ,:11,1~,1I{1I ' 11111. 1 :31; I 1‘ ‘1 I
11.11)? ‘ :: 11:1? 1113?. II (I.
311111111“ ~ ”91

f '?.|l0!.1“1A.‘ ‘3‘ yo I
. .. ' .

I11. ‘
I.’ °.
' IL
‘A
. '1 “I.“ . 9
0 .0 ‘5‘.
1 II"
0.? IV
I
.,;. ':

   
  
   
 

  
 
 

  
 

VII!

1111s;

 

    
 
  

l '1‘ 'I 1‘: I .
' :11 l: I
.1

II" I
("'1' {11h '11I'
1 IIII"1'§,..,

III1HH1111

11‘1'1HI‘IJ III "131'?” 1
II1"|.‘_'11':'E' :
. 1‘11 W1 1’ 111 '[1I111 ’1' 115111111 “1111”” 11!
11111.1(” i [11111111H111I1Ill1‘11t1111!

I
1111.111111N1IIJ1UI’ I1

“1' 1‘;

_f_.-:§_v Am— " "3;, -
M ». _ -. 4
—-—-""-‘:.- ..‘--‘ -

> i
:"‘T“

 
  

    
  

4. ,l ‘ ,
I 11".. 5‘1?" 1r ‘
1:::: ea,

“in .- .1
1 111...». 212111. .‘

I11” ‘1”: _ .-
11} 111111.111. f1:11 ’ I"
11'1“ [$1. 11131 1 1
Il ;

$1.11.?11sz1113111111 II 1
141?”; j '1 1111"‘211.:: , d1
“1111111111111‘111” 1g 1
1 11M VI}:

111111.11.”
141.5! 11111: ”[1th .:'_;1

11‘ 1'.‘ I' I
I1 II. I, .
I 1 1111131111"

       
 

7‘4—
‘3
*‘v

1

 
     
   
  
  

   
   

   
 

: k u '1 u “ .
:9" '1‘ ‘,

11 ‘91")? lfir”; I
M 3 11" l
‘ HI

12-: g I 3.- 11311133 S1t|

I 1w 1711;131W1 .
1‘ 1n .. IIIIg-
.. «I

#13311.» 111“”
. uZ‘ [f'lPIY1111'u1i’L 1‘
I431- “é: $156.3;
" 51-.‘,:.1:11111fi(1,r.11;17%1.111 r11

- .
:03.“

   

 

      

 
 
  

  
    
 
 

   
   

  

  
 

 

“I‘—
MW“:
mf-
V51»
..;._-
:g—g
-L.“
.~
..
.;
_..
m:-"

1111 {1 113911111.

 
   

     

 

 
 

1' 11I 1.0111111111111111: 1 11 11-;

p11I .111W111'1

I1I ll J"! "1
M111” WII1I13 "t" H 11
“II IIIIiliI'IIIi'I IELI'IEII'LII1III'

 

      
    

         
 

 

I}!

     

“Willi“Milli“W“\\\\\\\\\\|\\\\\\l\\\l

135519
3 1293 104

Juana-amen? wamggyl m mas-E.“-

I

 

This is to certify that the

thesis entitled
THE EFFECTS OF SUBDUCTING SLABS ON

SEISMIC AMPLITUDES AND TRAVEL-TIMES FROM
FOREARC AND OUTER RISE EARTHQUAKES

presented by

William John Rogers, Jr.

has been accepted towards fulfillment
of the requirements for

M.S. degree in Geology

14;“ * 3/:
a

ajor pr 8301

Date 4 May, 1982

0-7639

 

 

IVIESI.) RETURNING MATERIALS:
Place in book drop to
LJBRAfiJES remove this checkout from
4—! your record. FINES will

be charged if book is
returned after the date
stamped below.

 

 

 

 

er graft? If“
my: W 3' 7

arrow i

 

 

~—- *H*—~‘ H“ .

THE EFFECTS OF SUBDUCTING SLABS ON
SEISMIC AMPLITUDES AND TRAVEL-TIMES FROM
FOREARC AND OUTER RISE EARTHQUAKES

By

William John Rogers, Jr.

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements

for the degree of
MASTER OF SCIENCE

Department of Geology

1982

// "a

K15

ABSTRACT
THE EFFECTS OF SUBDUCTING SLABS ON

SEISMIC AMPLITUDES AND TRAVEL-TIMES FROM
FOREARC AND OUTER RISE EARTHQUAKES

By

William John Rogers, Jr.

The assumption that subducting slabs are too far from earthquakes
occurring in the forearc and outer rise to cause teleseismic
mislocations is investigated for the Aleutians using seismic ray tracing.
Results indicate that effects should be seen for earthquakes as far out
as the outer rise, but these effects become lOwer in intensity and
affect only a small area. Thus, teleseismic hypocentral locations are
probably relatively accurate seaward of the trench. Theoretical slab
effects include a low seismic amplitude zone with a high amplitude
band at its edge. Observed amplitude data agrees fairly well with
theoretical calculations and indicates the model is generally correct.
The subducting slab also distorts focal mechanism solutions. Short
period waveforms for an earthquake are consistent between stations, and
focal depths determined using depth phases for forearc and outer rise

earthquakes are generally shallower than ISC determinations.

In memory of

Dr. William J. Rogers, Ed.D

ii

ACKNOWLEDGMENTS

I would like to thank my committee chairman and advisor, Kazuya
Fujita, who did everything which could be expected to see this study
through, and more. He also taught me that foreign foodstuffs can be
delicious even if you can't pronounce them. I would like to thank the
other members of my committee, F. William Cambray, Masaharu Kato, and
James W. Trow, for their help, criticism, and encouragement. Thanks
to the faculty and staff at both Northwestern University and Michigan
State University for their understanding and cooperation. Mom paid for
part of this; so did the Department of Geology of Michigan State
University, and I am grateful to both parties.

This research was supported in part by National Science Foundation

grant EAR 80-25267 and by Petroleum Research Fund grant 612366.

TABLE OF CONTENTS

Chapter 1: Introduction ........................................ 1
1.1 The nature of forearc and outer rise earthquakes ........... 5
1.2 Rheology of the lithosphere ....... ' ......................... 7
1.3 Regional stress in the lithosphere ......................... 15
1.4 Effects of subducting slabs on earthquake location ......... 16
Chapter 2: Real and Theoretical Slab Effects ................... 26
2.1 Ray tracing and model slabs ................................ 26
2.2 Ray tracing results ........................................ 29
2.3 Amplitude effects .......................................... 47
2.4 Waveform effects ........................................... 60
2.5 Effects on focal mechanism plotting ........................ 71
2.6 The practical limits of slab mislocation effects in

the Aleutians .............. 71
Conclusions ..................................................... 77
Bibliography .................................................... 79
Appendix I: Program "Raytrak" ................................... 84
Appendix II: Using "Raytrak" .................................... 107

iv

FIGURES

1. Index map of study area ....................................
2. Adak local seismological network ...........................
3. Sample Raytrak plot output .................................
4. Thermal model of slab ......................................
5. Seismic velocity model of slab .............................
6. Ray tracing results ........................................
7. Ray path sketch ............................................
8. Relative amplitude vs. zero residual line ..................
9. ISC earthquake depths vs. depth phase earthquake depth
determinations.
10. "Escalating" waveform arrivals at several
seismological observatories ..
11. "Clean" waveform arrivals at several
seismological observatories ..
12. Ray emergence points corrected for slab ray path distortion
effects vs. the same emergence points uncorrected ...........
13. Stations detecting the earthquake of February 27, 1970,
showing which ones would have shown slab effects for
different values of 6 ......................................
14. Stations detecting the earthquake of May 30, 1967, showing
which ones would have shown slab effects for different
values of 6 ...............................................
15. Percentages of stations affected for different values of 6..

61

63

65

7O

7O

72
76

TABLES

1. Earthquakes used in this study ............................. 50
2. WWSSN stations in networks-locations in
degrees/minuites/seconds/ ............. 51

vi

CHAPTER 1
INTRODUCTION

Precise earthquake location is vital to many fields of geophys—
ical study. Often the seismicity of a certain area and its associated
effects are the main or only information available for investigation
of the local tectonic processes. This is especially true of tectonic-
ally active offshore areas such as island arc and trench systems. The
distribution of observed seismicity was originally the main argument
for explaining these island arc-trench systems as the sites and surface
expression of subducting lithosphere (Isacks et al., 1968). As the
precision of hypocentral locations increased, the extreme thinness of
Benioff-Wadati zones became apparent (Isacks and Barazangi, 1977).
More recently, further refinement in earthquake location procedures
has allowed the discovery of double seismic zones in subducted slabs
(Umino and Hasegawa, 1975; Barazangi and Isacks, 1979; Fujita and
Kanamori, 1981).

This thesis investigates the accuracy of the conventional earth-
quake hypocentral location method as applied to earthquakes in the
forearc and outer rise regions of a subduction zone. It has been
generally assumed that these earthquakes are well located by teleseis-
mic means, and that no correction for the effects of the subducting
slab on signals from these earthquakes needs to be made (Sleep, 1973).
An attempt will be made to decide whether this assumption is correct.
It is hoped that it will be possible to improve the accuracy of
location of these forearc and outer rise earthquakes. This would im-

prove our knowledge of the subduction process, the state of stress at

1

FIGURE 1.

INDEX MAP OF STUDY AREA.
Circles indicate ISC epicenters for earthquakes considered in

detail in this thesis.

gou<
.
ma.

L1

.Somh_

co
0 00.

new.

moz<mmoz<

.m.

IT

Lum—

k<m_

womb.

FIGURE 1

W?"

m" "of!

 

 

 

 

I7 '32.

 

 

 

 

 

 

u rzwo
° 73
“30.10 A
:00 ’ -
C
.72 4“; 20
4mm,” ..4| .
/o/.z—.2 2.3:: o"
«o
A (3‘3
<‘/ V '
nooo Nam

 

 

FIGURE 2

ADAK LOCAL SEISMOLOGICAL NETWORK

Triangles indicate instrument locations for this network.

convergent plate boundries, the structure of the subducting slab, and
the rheology of the lithosphere.

The earthquakes studied in this thesis occurred in the central and
western Aleutian Islands, Alaska (Figure 1). Studies of this area are
made easier by the nuclear tests which were carried out on Amchitka
Island. In addition, seismicity and crustal structure are well cali-
brated at the eastern end of the study area, in the region covered by
the Adak local microearthquake network (Figure 2).

Portions of this thesis have been presented as Fujita et al.(1980).

1.1 The Nature of Forearc and Outer Rise Earthquakes

Most of the earthquakes in a subduction zone occur in the shallow
thrust region. These are the most powerful earthquakes as well as the
most frequent. They are not the only earthquakes present in subduc-
tion zones, however. Earthquakes also occur downdip of the shallow
thrust region within the subducting slab. These concern this study
only in as much as it is by their presence that the length of the
slab is generally estimated. The other region of seismicity in the
subduction zone is the forearc and outer rise.

Earthquakes can occur anywhere in the forearc or outer rise areas,
but tend to be concentrated near the trench axis (Forsyth, 1980).
These earthquakes tend to have normal fault mecnanisms, indicating
their being caused by tensional forces. Therefore, it has long been
assumed that these earthquakes are caused by the bending of the litho-
sphere before subduction (Stauder, 1968a, 1968b). There was some
thought that a portion of these earthquakes might be in the accretion-

ary prism, associated with the extreme deformation this feature is

subjected to, but there seems to be no evidence which conclusively
indicates any earthquakes occurring in the accretionary prism at all
(Chen et al., 1981). Because of this information, it appears that the
earthquakes in the forearc and outer rise are caused by the bending in
the slab as Stauder (1968a, 1968b) stated. This occurs because when
the plate is bent, having a finite thickness rather than an infinites-
imal one, only the "neutral plane" approximately midway between the
top and bottom of the slab bends without being otherwise stressed.
Since the slab is bent downward, the part of the slab above the neutral
plane is put in tension with the amount of stress depending linearly
on the distance above the neutral plane. The part of the slab below
the neutral plane is in compression, and again the amount of compres-
sion depends linearly on the distance from the neutral plane

(Johnson, 1970). The exact position of the neutral plane depends on
the regional stresses, if any.‘ A regional compressional or tensional
stress will displace the neutral plane upward or downward, respect-
ively (Johnson, 1970). The position of the neutral plane should

also depend on the properties of the material being bent. If these
properties are not uniform, some displacement of the neutral plane from
the ideal position should be seen.

It seems that far more tensional earthquakes than compressional
ones are recorded from the forearc and outer rise. If the plate being
bent is of uniform properties there should be equal numbers of
compressional earthquakes at the bottom of the slab. Several possible

reasons for this disparity exist:

[I] There are equal numbers of compressional earthquakes,
but for some unknown reason they are not detected. In the
absence of an explanation of why this should be so, it seems
unlikely.

[2] The lower parts of the lithosphere are hotter than the
upper (Toksoz et al., 1973) and may distort plastically with-
out breaking sharply to cause earthquakes.

[3] A combination of the above may occur.

1.2 Rheology of the Lithosphere

Better location of forearc and outer rise earthquakes would help
determine the rheology of the lithosphere. Even basic questions about
this, such as the thickness of the lithosphere and the definition of
what the lithosphere is are more difficult to answer adequately than
might at first appear.

The lithosphere is generally defined as the crust of the Earth and
that part of the mantle above the low (seismic) velocity asthenosphere.
The thickness of the lithosphere is thought to be about 100 km. The
exact thickness of the lithosphere varies greatly from place to place.
The lithosphere thickens as a function of cooling and therefore of age
(Caldwell and Turcotte, 1979; Baer, 1981). An extreme example of the
range in thickness occurs in marginal seas such as the Japan Sea, where
the lithosphere has a thickness of about 30 km while the nearby
Pacific Ocean lithosphere has a thickness of about 75 km (Segawa and
Tomoda, 1976).

The problem with this is the question of what method to use for
determining the thickness of the lithosphere. One way is to use seismic

methods (e.g. Jin and Herrin, 1980). These use observations of seismic

waves to determine where the lithosphere ends. However, it is not
certain that materials which behave elastically over very short seismic
time spans will also behave elastically over years or millions of years.

There are also thermal methods. The thermal method involves
assuming conditions at the midocean ridge axis, cooling rates of the
lithosphere, and from these calculating the changes in lithospheric
thermal structure with time. In the absence of conclusive data about
the composition and temperature-pressure conditions of the lithosphere
below relatively modest depths, it is unfortunately impossible to
decide whether a given temperature is in itself an absolute indicator
of rheological behavior.

Another problem is that of time scales. It is believed that rocks
behave differently over different time spans. While the behavidr of
rocks at high temperatures and pressures can be tested in the labora-
tory ( although not to the extent of duplicating conditions very deep
in the Earth) it is impossible to experimentally check the behavhdr of
rocks on the time scale of millions or even hundreds of years. Liu
(1980) points out that a common assumption of lithospheric models based
on elastic thin plate theory is that the upper part of the mantle
behaves elastically over long periods of time while the lower part does
so only over very short periods, behaving viscously over longer
periods of time. Watts et al. (1980) , using thin elastic plate
theory, claim that observations indicate that the mechanical thickness
of the lithosphere is only one half to one third it's "seismic" or
"thermal" thickness. Such assumptions are needed to fit observed data,
which require a perfectly elastic plate of around 20—30 km thickness

and relatively high strength, while seismic studies and thermal

models indicate lithospheric thickness to be much greater.

The mechanical thickness of the lithosphere can also be estimated
by observing its behavior under certain loads and bending moments.
These include the benavior of the lithOSphere at subduction zones, its
subsequent straightening, and its deformation under point loads such
as seamounts (e.g. Bodine et al., 1981; Forsyth, 1980; Caldwell et al.,
1976; Hanks, 1979). These studies have the advantage of avoiding the
problems of time dependent behavior of rocks to some extent, since the
rocks studied have been under stress for very long periods of time,
being real world examples and not laboratory samples. Results of this
type of study, however, are very model dependent. Estimates of the
mechanical thickness of the lithosphere depend largely on the assumed
overall rheology of the lithosphere. Since the results of these
studies are model dependent there has been a large degree of success
' among various authors in fitting the observed behavior of the litho-
sphere to wildly different rheological models, achieved by altering
the properties of the lithospheric materials to fit the data without
ever altering the assumed model.

Many different models for the rheology of the lithosphere have
been pr0posed. These include a simple thin elastic plate, an elastic
plate under regional compression, an elastic-perfectly plastic plate
with or without regional compression, an elastic-perfectly plastic
plate with variable strength in different layers, and an all-viscous
model (Forsyth, 1980) as well as more complex models such as that of
Beaumont (1979) which combine several distinct depth-dependent
rheologies.

The distribution of forearc and outer rise earthquakes is already

10

well enough known to disqualify some rheological models. For example,
consider the thin elastic plate model (e.g. Caldwell et al., 1977).
This is the simplest rheological model for the lithosphere generally
proposed. It can explain the geometry of the outer rise fairly well,
and can do so without requiring large amounts of regional compression
in the oceanic crust; it can explain the geometry of all subduction
zones if it is assumed that those with shorter than normal wavelengths
are under regional tensile stresses of about 10 kilobars and bending
plate stresses of about 15 kilobars (Forsyth, 1980). Forsyth (1980)
claims that it is somewhat unreasonable to assume that the lithosphere
can handle such stresses. Many authors claim that the lithosphere
either can not or does not support such stresses. For example, Watts
et al. (1980) claim that the lithosphere can withstand a much more
modest 1 kilobar stress over long periods of time. Kirby (1980)
believes on the basis of laboratory experiments that the ultimate
compressive strength of the lithosphere in its strongest layer is about
8 kilobars; since this is only for the strongest layer, the overall
average strength of the lithosphere would be much less. Lambeck (1980)
believes on the basis of observation of the crust's configuration under
loads and considerations of isostasy that the average crustal stress

is about 500 bars and the maximum stress differences between maximum
and minimum on a given piece of the lithosphere is about 1-1.5 kilo-
bars. Solomon et al. (1980) state that average differential stress in
the lithosphere is only about 200-300 bars. Liu (1980) believes the
lithosphere to be even weaker, being able only to support differential
stresses of around 50 bars. Because of these considerations, it seems

unlikely that the lithosphere behaves in a perfectly elastic manner.

11

None of this, however, is able to completely disprove the con-
tention that the lithosphere behaves as a thin elastic plate. The most
convincing argument against thin elastic plate theory comes from con-
sideration of the distribution of earthquakes in the forearc and outer
rise. If these earthquakes are caused by the bending of the plate,
they should be evenly distributed across the area of bending under thin
elastic plate theory. This is because a bending perfectly elastic
member will have an even amount of stress everywhere across the area
of its bending, and so earthquakes will be equally likely anywhere in
that area. The present knowledge of the distribution of forearc and
outer rise earthquakes is good enough to show that these earthquakes
are concentrated in the forearc region, not distributed evenly across
the outer rise, which fact effectively disproves the thin elastic plate
theory of lithospheric behavior (Forsyth, 1980).

Viscous models of the lithosphere have also been proposed. These
models can explain the topography of the trench and outer rise under
all conditions. They do not require large bending stresses; a stress
of less than one kilobar is sufficient. However, these models are not
able to explain the fact that the bending lithosphere can hold a fairly
high stress without relaxing over a period of millions of years, as in
the Puerto Rico trench (Forsyth, 1980) or in the case of the Michigan
Basin, where the stresses caused when this feature first appeared have
not relaxed significantly in 450 million years (McAdoo et al., 1978).
Therefore, it would seem that the lithosphere does not behave in a
completely viscous manner over time scales significant in regard to
plate tectonics.

The simplest mechanical model for the lithosphere which seems to

12

fit all the observed data relatively well is the elastic~perfectly
plastic model (e.g., McAdoo et al., 1978). This model suggests that
the lithosphere behaves elastically up to a critical yield stress.
Below this yield stress any deformation undergone by the lithosphere
will be recovered after the stress causing the deformation is removed.
At the yield stress the material breaks and deforms in a perfectly
plastic manner, causing permanent deformation. A model based on a
simple single-layer elastic-perfectly plastic model still requires
fairly large regional compression to explain all the observed topography
of the sea floor (McAdoo et al., 1978). Again, the distribution of
forearc and outer rise earthquakes would tell us much about whether
this rheological model of the lithosphere is acceptably close to the
real behavior of the lithosphere or not. Forsyth (1980) observes that
tensional earthquakes in the forearc and outer rise extend to at least
a depth of 25 km, even in subduction zones which are assumed to be
under strong compression, and that this is deeper than earthquakes
could occur under an elastic-perfectly plastic single layer model.
Unfortunately, this is cutting the observations very close to their
limit of accuracy. In the areas of subducting slabs, errors on the
order of 10 km in hypocentral depth are believed to occur using
teleseismic hypocentral location proceedures (Fujita et al., 1981);

if similar errors occur for earthquakes in the forearc and outer rise,
Forsyth's objection to the single layer elastic-perfectly plastic
model (1981) may not be valid. In order to make such fine distinctions
it would be useful to determine whether these hypocentral location
methods are really as accurate in the forearc and outer rise as seems

to be generally believed.

13

The rheological model which Forsyth (1980) prefers to the single
layer elastic-perfectly plastic model on the basis of this fine dis-
tinction of earthquake depths is the elastic—perfectly plastic variable
strength model. Under this rheology, the lithosphere.behaves as under
the simple elastic-perfectly plastic model except that different layers
of the lithosphere have different critical yield stresses. This model
reduces or eliminates the need for regional compressive stresses
(Chapple and Forsyth, 1979) and can explain all the observed topography
of the ocean floor at least to the level of accuracy now availible; it
can also explain the distribution of the earthquakes in the forearc and
outer rise (Forsyth, 1980).

These more or less simple rheological models for the lithosphere
are useful in predicting the behavior and seismicity of the lithosphere
under different conditions. However, it seems unreasonable to assume
that any simple model with only one or two discrete layers explains
the actual behavior of the lithosphere in every way. Aside from phase
boundaries,temperature and pressure which vary continuously with depth
should lead to continuously variable physical properties rather than
discrete boundaries.

Several authors have proposed quite involved rheological models
for the lithosphere. Beaumont (1979) has developed a rheology based
on the olivine flow laws observed in the laboratory, extrapolated to
higher temperatures and pressures. He assumes that the temperature of
the lithosphere is a linear function of depth in any given location.
This leads to a rather complex rheological cross section for the
lithosphere. In a typical area with lithosphere 100 km thick, the

crust might behave plastically via cataclastic flow down to a depth of

14

12 km, at a temperature of 180° C. Below this, to a depth of about 58
km and a temperature approaching 800° C, the crust behaves in a pseudo-
elastic manner. The rheology of this "pseudo-elastic“ area is based

on diffusional (Newtonian) creep. Given time the lithosphere will
creep out from under any stress applied to it, but since the strain
rate for this zone is slow it will take a very long time, up to a
billion years. Below this level the lithosphere behaves in an increas-
ingly plastic manner until the base of the lithosphere is reached at
about 100 km, where the temperature is around 1,300° C.

The time scales involved im Beaumont's (1979) model can be very
long. When a load is applied to the lithosphere, plastic yielding in
the lower part of the plate reduces the apparent elastic thickness of
the lithosphere from its thermal or seismic thickness to at most half

6 or 107 years. In

this value. This takes place over a period of 10
the upper layers of the lithosphere, diffusional creep will reduce the
apparent flexural rigidity over time scales as long as 109 years.
This is long enough so that the oceanic crust will appear to be elastic
over plate tectonic time scales. However, continents (which persist
longer than oceanic crust) will eventually be allowed to reach a
relatively unstressed isostatic equilibrium over most of their area,
as is observed (Beaumont, 1979).

Bodine et al. (1981) have proposed a model of the lithosphere with
three different rheologies. They use an elastic-plastic rheology based
on cataclastic flow to describe the lithosphere at shallow depths. At

greater depths the authors propose a temperature activated creep as

the factor which governs yield stress. For shear stresses less than

15

2 kilobars, power law creep dominates and is expressed as

e =7O(os)3-exp(-Qp/RT)
where c is the strain rate, as is the shear stress, T is the temperature
in degrees Kelvin, R is the universal gas constant and 0p is the activa-
tion energy. At shear stresses greater than 2 kilobars, power law creep
breaks down and the authors state that the Born law works better,
expressing the creep as

E = 5.7-1011-exp(-Qd/RT[I - oS/a.s-1o“]2)
with all notation as before, except for Qd' Qd is an activation energy
not equal to 0p, chosen to assure that both laws agree at 2 kilobars.

Whether the rheological model considered is meant to describe the

behavior of the lithosphere as a whole or actually to describe the
behavior of rocks at all the specific levels within the lithosphere,
accurate data about observed lithospheric behavior is necessary.
Sufficiently good data might make it possible to decide which simple
model works best and whether the more complicated models are generally
worth the extra trouble they cause. In any case, more accurate
location of earthquakes in general and forearc and outer rise earth-

quakes in particular will help answer these questions.

1.3 Regional Stresses in the Lithosphere

A question which is closely related to rheology is the question of
the regional field of stresses throughout the lithosphere. There are
various means of approaching this problem. Lambeck (1980) uses topo-
graphy and isostatic considerations to estimate the average stress on
the lithosphere. This method is unfortunately dependent on the

rheology assumed. Another way of approaching the problem is by

l6

investigation of intraplate earthquakes (Richardson et al., 1976,1977;
Richardson and Solomon, 1977) on the assumption that these reflect
accurately the regional stress field. Unfortunately, earthquakes in
midplate or passive margin regions may occur along pre-existing fault
systems. These may date from process on midocean ridges; they are not
necessarily well aligned with the present day stress field, but represent
a direction of weakness along which slip can easily occur. The mech-
anisms of the resulting earthquakes may relate to the present-day stress
field in only a relatively vague manner (Stein, 1978,1979). This is an
especially perplexing problem, since a better understanding of the
worldwide regional stress fields would probably help not only in
determination of lithospheric rheology but might give information on the
driving forces of plate tectonics as well (Richardson et al., 1976).

The seismicity of the forearc and outer rise is one of the best
sources of data to solve these problems. As was shown in the last
section, seismicity is one of the chief sources of data for determining
lithospheric rheology. The questions of rheology and regional stress
are so closely intertwined that to answer one is a large part of answer-

ing the other.

1.4 Effects of Subducting Slabs on Earthquake Location

One of the major sources of error in hypocentral location is the
lateral inhomogeneity of the upper mantle seismic velocity in and near
subduction zones. The effect of this inhomogeneity has long been
detected for earthquakes in the shallow thrust region of various subduc-
tion zones (e.g. Davies and McKenzie, 1969; Biswas, 1973; Suyhiro and

Sacks, 1979; Fujita et al., 1981). These mislocation effects have also

17

been observed in connection with nuclear explosions. It was observed
that teleseismic location of nuclear tests in the middle of the Pacific,
far from subduction zones, are accurate to within 5 km of the actual
focus (Sykes, 1966) while similar solutions for tests in the Aleutians
tend to be mislocated farther north and deeper than the actual foci,
presumably due to the effect of the subducting slab beneath these
islands (Sykes, 1966; Herrin and Taggert, 1968; Jacob, 1972). The
disagreement between local seismic network hypocentral locations and
teleseismic hypocentral solutions for the same events can be on the
order of 50 km (Murdock, 1969; Fujita et al., 1981).

The cause of these effects is fairly well known, although their
magnitude generally is not. The subducting slab is cooler, denser, and
more rigid than the ambient mantle around it. Because of these dif-'
ferences the seismic velocity is assumed to be higher in the slab than
near it. Seismic waves passing down the subducting slab travel faster
than they would travelling outside the slab (Engdahl et al., 1977).
Because most hypocentral location methods assume that the Earth is
spherically symmetric, a lateral inhomogeneity such as a subducting
slab causes earthquake hypocenters to be determined erroneously.
Specifically, the seismic stations downdip of the slab receive signals
early, and since the symmetric earth model assumes that a short travel
time indicates a short distance, earthquakes occurring in the shallow
regions of a subduction zone appear to be further down the slab than
they actually are. This results in a mislocation both in epicentral
location and depth of focus, often on the order of 50 - 100 km (Fujita

et al., 1981). Conversely, earthquakes deep in the subduction zone can

18

send seismic waves upward along the slab, causing them to be mislocated
(Spencer and Engdahl, 1981).

In order to improve the accuracy of hypocentral locations it is
necessary to remove this slab effect. An important means of doing
this is the use of local seismological networks such as the one in
operation on and near Adak Island in the central Aleutians (Figure 2).
These local networks are able to detect smaller earthquakes than the
teleseismic seismological stations are. More important from the point
of view of earthquake location methods, they are presumably largely
free of slab effects. Seismic waves travelling from earthquakes in the
shallow thrust zone do not travel very far if at all in the subducting
slab. This should presumably make hypocentral locations obtained from
local networks more accurate than teleseismic solutions for the same
event (Fujita et al., 1981). This provides a check on the accuracy of
the teleseismic hypocentral locations and a way of determining the slab
effects.

Unfortunately, this method may not be accurate for forearc and
outer rise earthquakes. These earthquakes are sufficiently far from
any land-based local network so that their seismic signals follow
deeper paths than those from shallow thrust zone earthquakes on their
way to the local network seismological stations. These deeper ray
paths may pass through the lower parts of the lithosphere or inhomogen-
eous upper mantle, with unpredictable effects on travel times and
amplitudes. For these reasons, local network hypocentral solutions for
forearc and outer rise earthquakes are considered relatively unreliable.
This makes it impossible to meaningfully compare the locations

obtained from local networks with those from teleseismic ones.

19

Another problem in attempting to accurately locate forearc and
outer rise earthquakes is that their magnitude is small relative to
thrust zone earthquakes. Therefore, their initial signal strength is
lower. In addition, Sleep (1973) has noted an extensive region of low
seismic amplitude arrivals at teleseismic distances downdip of the
slab. The seismic waves from forearc earthquakes which pass through
the slab generally fall into this "shadow zone" and are therefore
further weakened. This lessens the signal to noise ratio. The signals
detected at seismological stations in this shadow zone are likely to
be emergent, making accurate determination of arrival (and travel)
times difficult. This may create a tendency to pick arrival times as
being later than they actually are.

It has generally been assumed that the forearc and outer rise are
far enough from the subducting slab so that the slab does not affect
their signals. This view is strengthened by the fact that little or no
systematic bias in travel times has been observed. This may be an
illusion, due to the possible tendency to pick arrival times late and/or
the difficulty of determining errors in teleseismic hypocentral
locations for these earthquakes by comparison with local network
solutions. The first would not be a matter of much concern if the
tendency to misread arrivals exactly counteracted any slab effect, but
this is not to be expected. The error caused by reading earthquakes
differently at the seismological stations in the "shadow zone" is
probably random, and if it cancels any slab effect it does so only by
chance.

One way to examine the effect of the subducting slab on seismic

waves from forearc and outer rise earthquakes is through seismic ray

20
tracing. This is a method of calculating the path followed by a seismic
wave passing through various bodies with different seismic velocities.
The simplest version of this is Snell's Law, which predicts the change
in ray path between two discrete layers. There are various approaches
to ray tracing. Sorrels et al. (1971) used Snell's Law by tracing rays
through a medium made up of a series of discrete constant-velocity
layers separated by planar or spherical interfaces. They used this
method to examine the anomalies in teleseismic data for the nuclear
test LONGSHOT, set off on Amchitka Island in the Central Aleutians.
This method works but can only handle velocity gradients by using
large numbers of very thin layers. This leads to some lack of pre-
cisiOn and some increase in digital input.

In order to use gradients of seismic velocity rather than sharply
bounded layers, it is necessary to use differential equations. The
equations developed to predict the path of a ray in an arbitrarily
inhomogeneous three-dimensional medium arise from Fermat's principle of
stationary time (Jacob, 1970; Julian, 1970) from which, in the
simplest case, Snell's Law also arises.

This thesis examines the effects of a subducting slab on seismic
signals from forearc and outer rise earthquakes by using the ray-tracing
method of Julian (1970). This involves a mathematical solution of
Fermat's Principle in a three dimensional arbitrarily inhomogeneous
medium. The sulution consists of five simultaneous first order
differential equations which give the variation with time of ray
direction, based on seismic velocity and its spatial derivatives at the
point in question. This method will be used only to find the ray path

and travel time, hence the travel time residual. It is also possible

21

to quantitatively determine the effect of geometric spreading of the
rays on the seismic amplitude at any point, but this requires consider-
ation of ten additional differential equations and a large amount of
additional calculation (Julian, 1970). For this reason, the calcdla4
tion of seismic amplitudes based on Julian's equations is not attempted
here. Instead, amplitudes will be found qualitatively using the
observation that areas from which rays are refracted show a relatively
low density of emergence points; ray emergence point density is there-
fore associated with seismic amplitude (Sleep, 1973). If a distribution
of rays uniformly spaced in azimuth and take-off angle are traced
throughout their paths and the emergence points of these rays at the
surface of the Earth are plotted, the ray emergence points should be
concentrated in areas of higher seismic amplitude arrivals and thinned
out in any "shadow zone" (Sleep, 1973).

Even ignoring the quantitative determination of amplitudes, the
ray tracing equations of Julian (1970) require large amounts of calcula-
tion. The cost of using this procedure for large numbers of individual
earthquakes would be high. However, unless the slab is contorted or
torn, its effects are a smooth function of epicentral distance from the
center of arc curvature for any given constant focal depth (Fujita et
al., 1981). It is therefore only necessary to determine the effects of
the subducting slab at various evenly spaced points across a range of
distances and to interpolate for slab effects at intermediate points.
Once a table of slab effects has been produced in this manner for a
sufficiently accurate slab model, it should no longer be necessary to
actually trace rays at all in order to obtain the improvement in

accuracy ray tracing provides (Fujita et al., 1981).

22

FIGURE 3.
SAMPLE RAYTRAK PLOT OUTPUT .
Small crosses indicate ray emergence points for rays calculated in

previous runs. Large crosses are for the rays most recently calculated.

DIST. IS 12.60

20 RAYS.

9N6.

 

."---'.-..

23

O
b
0
a
O
\

FIGURE 3

~
0
o ‘ ‘
a ‘ ‘
.0
. C
O
o ".
O
o o '
o ‘ '
o o
o’. O
‘ o

 

 

 

24

FIGURE 4 ,
THERMAL MODEL OF SLAB (from Fujita, Engdahl, and Sleep, 1982)
Triangle at zero depth and distance of 11.4 represents location of

volcanic front. Dotted outline is the outline of the subducting slab.

FIGURE 5
SEISMIC VELOCITY MODEL OF SLAB (from Fujita, Engdahl, and Sleep, 1982)
Triangle at zero depth and distance of 11.4 represents location of

volcanic front. Dotted outline is the outline of the subducted slab.

DEPTH (km)

25

DISTANCE FROM POLE (deg)
ll " l2 l3

L

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

400—
' zoOO°c-'
500
FIGURE 4
DISTANCE FROM POLE (deg)
o no u g :2 l3 :4
l r . l 1"“ l 4:590 I
4 6.75
:00 —
a 25——-—~
E zoo ._
$5
I a
F use
B
o 300 .7: —
9.00
400 — -
9.25
9,50 kin/sec—
soc

 

FIGURE 5

CHAPTER 2
REAL AND THEORETICAL SLAB EFFECTS

2.1 Ray Tracing and Model Slabs

The ray tracing for this thesis was done using a computer program
first developed by Julian (1970) and modified thrice since (see
Appendix I). This program operates on the “shooting" system, sending
rays out from the source at known aximuths and take-Off angles to
observe the results rather than trying to home in on the ray path
between a given source point and receiver point (Julian and Gubbins,
1977). It was modified first by Sleep in 1972, then for interactive
use by Fujita in 1978, and most recently by the author for automatic
plotting of the ray emergence points using a Tektronix graphics
terminal. As modified, the program uses ready-made Tektronix routines
to create a polar projection plot showing the distance and azimuth of
each ray emergence point from the earthquake epicenter. It shows these
points as two sizes of crosses; the rays most recently ordered by the
Operator as large crosses and any previous rays as.small ones. It also
produces a simple polar coordinate grid for easy reference. The
resulting plot (Figure 3) provides a useful display of what has been
done and what remains to be done.

In order to use ray tracing, a seismic velocity nodel for the sub-
ducting slab is needed. This thesis uses a thermal model (Figure 4)
of the slab calculated using the methods and parameters of Tokst et
al.(1973). The specific model used for the Aleutian slab in this
thesis was developed by Fujita et al. (1981). A dip angle of 57.3° and
a convergence rate of 4.5 cm/yr. were assumed. The model slab is 300

KM long and 80 km thick.

26

27

Fujita et al (1981) constructed a velocity model based on this
thermal model (Figure 5). This was done by assuming that the
difference in seismic velocity between the slab and the normal mantle
velocity distribution depends linearly on the temperature difference
between them. The seismic velocity at any point is given by Sleep
(1973) as:

V(x,z) = Vo(z) + aV/aTE T°(z) - T(x,z)] +

z aV/mam [mo(z) - m(x.z)]

where:
V = seismic velocity,
T = calculated temperature,

o(subscript) indicates an unperturbed quantity,
x = horizontal coordinate,

vertical coordinate, and

Z .

m the amount of a given phase present.

However, Fujita et al. (1981) point out that the final term of
this equation ( which deals with velocity differences due to phase
changes between slab and ambient mantle) is small with respect to the
others, and is probably too small to be significant considering the
data currently available. Because of this, the last term of the
equation above is dropped, and it becomes:

V(x,z) = V°(z) + aV/aT[ T°(z) - T(x,z)]
with the notation as before. i

The value for the derivative of velocity with respect to temp-
erature is assumed to be constant; the value used for aV/aT was

-0.009 km/sec-°C. This is the value used by Fujita et al. (1981).

The resulting velocity model of the slab is shown in Figure 5.

28

There is some question of how long the slab should be, since this
depends on at least two poorly known factors. These are the age of
the lithosphere (which determines its thickness) and the duration of
subduction. In addition, it is not known whether the slab ends where
deep seismicity ends, as most authors seem to believe, or if it con-
tinues far beyond this point (Sleep, 1973; Jordan, 1977). The slab
model used here is 300 km long. This is longer than the slab model
used by Sleep (1973) for the Aleutians, in-which paper the theoretical
subducting slab had little or no effect on seismic signals from the
forearc and outer rise. It is far shorter than Sleep's (1973) slab
model for the Tonga arc which he felt did not fit the Observed data for
the Aleutians as well as the shorter slab model. It is somewhat less
than the length of slab indicated by a recent velocity inversion
study of the Aleutians, where the data seemed to favor a slab of about
300 km in length, with the best single fit to the data occurring for a
slab 386 km long (Spencer and Engdahl, 1981). It has been observed
that deep seismicity in the Aleutians extends to a depth of somewhat
over 250 km (Engdahl, 1977). It would seem that 300 km is a reasonable
length for the slab model based on the data currently availible.

The slab model used here is 80 km thick. This thickness may be
doubtful. The oceanic lithosphere subducting near Adak Island in the
central Aleutians is in the neighborhood of 65 million years Old
(Pitman et al., 1974). Oceanic lithosphere of this age should be
around 60 km in thickness (Yoshii et al., 1976) rather than the
assumed 80 km.

Assuming that the subducting slab forms a symmetrical arc about

its center of curvature, which is a reasonable approximation in most

29

areas (Fujita et al., 1981), the only geographical variable which bears
on whether or not seismic signals from a given earthquake will be
affected by the slab is the distance between the earthquake and the
slab. This can be given in terms of the distance between the earth-
quake and the center of arc curvature, defined as a geocentric angle.
This was done by Fujita et al. (1981), who designated this distance as
6 ; it will be so desingated in this thesis. The center of curvature
of the arc is taken as 63.314° N, 178.059° W (Engdahl, 1973), and the
bottom hinge point of the bent slab is set at 6 = 12.5 (Fujita, 1979)

for the slab model.

2.2 Ray Tracing Results

Rays were traced through the slab model from a fixed source depth
of 25 km and variable distances from the center of arc curvature. In
the Aleutians, the trench axis falls at around 6 = 13.0°; the range of
6 covered in this thesis was from 6 = 12.4’to 6 = 14.6°. This covers
the forearc, trench, and outer rise. Ray tracing was done across this
area at intervals of O.2° of 6, except for the distance 6 = 13.1°,
which was checked in order to be certain that the spacing used was
adequately close for accurate interpolation of intermediate values.
It was found that the slab effect varies smoothly both as a function of
azimuth and of 6, confirming the identical Observation of Fujita et
al. (1981). Based on this, it appears that the O.2° spacing for values
of 6 is adequately close in this range of 6.

Results of the computer ray tracing are shown in figures 6A-6L.
On the right side of each figure is a contoured plot of the theoretical

travel time residuals indicated by the computer study. On the left are

30

FIGURE 6
RAY TRACING RESULTS

Polar distance plots showing slab effects for different values of
6. The radius of each plot is 100° of geocentric angular distance.
Left side of each plot is a contour of travel time residual in seconds
early, while the right side shows where the emergence points of rays of
a given take-off angle (italic numbers) and azimuth ( straight numbers)
would fall. Dark circles are added to these points in order to
emphasize areas of amplitude effects except where this would obscure

detail.

31

 

 

 

 

 

 

A = 12.4

FIGURE 6-A

32

 

 

 

 

 

 

‘—

 

12.6

8

FIGURE 6-B

33

 

 

 

 

 

8 - 12.8

FIGURE 6-C

34

 

 

 

 

 

uh

 

13.0

.8

-D

FIGURE 6

35

 

 

 

 

13.2

A

FIGURE 6-E

36

 

 

 

 

8 - 13.4

FIGURE 6-F

37

 

 

 

 

8 13.6

FIGURE 6-G

38

 

 

8 - 13.8

FIGURE 6-H

39

 

 

 

8 = 14.0

FIGURE 6-I

4O

 

 

 

 

14.2

8

FIGURE 6-J

 

41

 

 

 

 

 

 

8 - 14.4

FIGURE 6-K

43

shown lines of equal take-off angle. Running across these are another
set of lines, generally radial; these are lines of equal take-Off
azimuth. The azimuth lines are spaced 5° apart, the take-off angle
lines are spaced 2°. Radius on these plots marks the epicentral
distance ( A ) from the source, in degrees. Zero degrees azimuth
points in the downdip direction of the slab, which is by definition
toward the center of curvature of the island arc. As it has been
assumed that the slab is symmetric about its center of curvature,

only half each of the full-circular travel time and ray emergence
point plots are shown. The other half could be obtained by reflection
across the zero-azimuth line.

Figures 6A-6L, taken individually, show how slab effects vary with
azimuth and A; taken as a whole they show how slab effects vary with 6.
It can be seen on most Of the figures that although a smaller (i.e.near
vertical) take-off angle generally causes the ray to emerge farther
from the source, this is not always so. Over a narrow range an
increase in takeoff angle causes a decrease in A. This area marks the
edge of the slab effects. Between this narrow band and the seismic
source, the rays are going through a relatively great part Of the
subducted slab. Outside of this band, the rays are either going
through a relatively small amount of the subducted slab or are missing
it altogether, depending on the location of the earthquake (Figure 7).
Rays which pass down the length of the slab are refracted into paths
which take them much farther from the source than rays of their take-
off angle would normally go. The narrow band between the slab-affected
rays and the relatively unaffected rays contains those rays which

barely graze the bottom of the slab; in this area the slab effects

44

FIGURE 7
RAY PATH SKETCH (not to scale)

Dashed lines are a sketch of ray paths from an earthquake close to
the subducting slab. Dotted lines are ray paths from an earthquake

farther from the slab.

45

 

 

FIGURE 7

 

 

46

rapidly die away, so that for a small range Of take-Off angle a smaller
take-off angle causes the ray to hit less of the slab and emerge closer
to the source rather than farther away.

Slab effects also change greatly with 6.. The reason for these
changes can be seen fairly easily from the geometry of the situation
(Figure 7). If an earthquake occurs above the subducting slab, even
the rays with shallow take-Off angles will strike the slab and be
affected by it. They will strike the surface of the slab at a
relatively shallow angle, refract into it and continue within the slab
for a relatively long time. The effect of the slab on these rays will
be correspondingly great. As the seismic source is moved seaward, the
situation changes. Increasingly, rays with small take-off angles will
miss the slab completely and continue to their destinations unmolested.
Rays at some intermediate take-Off angle will strike the underside of,
the slab and refract down it. Since their angle of incidence at the
slab surface will be greater than that for rays from more arcward
seismic sources, these rays will follow a path which tends more to cut
across the slab than follow it. They will still be refracted further
from the source than rays with this take-Off angle would normally go,
but not as far as the rays from arcward sources with shallow take-off
angles would. Because of this, the area in which seismic signals are
affected by the slab contracts. In addition, not being so strongly
refracted, these rays are not geometrically spread as much, and so the
emergence point density within the area of slab effects is greater for
seaward sources than for arcward ones. With increasing values of 6

this process continues until only signals to seismic stations almost

47

directly above the slab would be affected by it; at this point the slab

effect has become a station effect rather than a source effect.

2.3 Amplitude Effects

The slab can have effects on travel time due to its anamalous
seismic velocity. It can also affect seismic amplitude, mainly by
geometric spreading of rays. The first Of these questions to be con-
sidered is that of predicted seismic amplitude. This is shown by the
density of ray emergence points at the surface of the Earth. These
points, plotted in figures 6A-6L, are the emergence points of rays
which were sent Out in an evenly spaced distribution of azimuths and
takeeoff angles, so that the ray emergence point density reflects amp-
litude (Sleep, 1973).

Were there no slab effect, the lines of equal take-Off angle would
describe concentric circles about the origin. They do not. Although
all points considered in this thesis are relatively far from the
subducting slab, there are still noticible slab effects over most of
the range in 6 considered. At 6 = 12.4°, for example, the emergence
pOints for rays sent out at a take-off angle of 33° emerge at values
of A which vary between 85° and 34° across the range Of azimuth (with
respect to the arc normal direction) of 0° to 90°. The azimuth of the
emergence point can vary from the take-off azimuth by as much as 30°.
The low emergence point density caused by this ray path distortion is
prominent. This shadow zone, reported by Sleep (1973) for earthquakes
closer to the island arc, is still present and even strong at this value
of 6 ; the associated low amplitude seismic signals should be seen at

seismological stations in this shadow zone as well. The ray tracing

48

studies also show a previously mentioned area of high density of ray
emergence points in a narrow band at the edge of this shadow zone. A
narrow band of unusually high amplitude seismic signals should be
associated with this, provided that rays travelling to this area along
different routes arrive in phase. As will be seen later, they should
do so.

Depending on the location of the earthquake both along the arc and
in 6, this area of high amplitude should be detected at stations in.
western North America or northern Europe. Due to the narrowness of the
high amplitude band, the alignment of earthquake, slab, and seismologic-
al station is critical to its detection.

The shadow zone loses intensity as the earthquake location is
moved away from the slab. At any value ofISwhere it was possible to
make a comparison, the emergence point density in the shadow zone is
less than that at the same value Of A but outside the shadow zone. At
the same time, the emergence point density in the shadow zone is higher
than it would be for the shadow zone observed at a lower value of 6.
As an illustration of this difference in density, consider the relative
ease of directing a ray so that its emergence point falls within a
desired 10° by 10° area. At 6 = 12.4° this can be done without having
to use increments Of take-Off angle or azimuth smaller than O.25°. At
6 = 13.0° a 1.0° increment was adequate for this purpose, while at 6
= 14.0° the necessary increment is never less than 2.0°. As distance
from the slab increases, it is becoming progressively easier to obtain
rays emerging in the shadow zone. Therefore the ray emergence point
density in the shadow zone increases accordingly, and with it the seis;

mic amplitude (Sleep, 1973). In other words, as the distance between

49

TABLE 1
EARTHQUAKES USED IN THIS STUDY

50

:Pumppsm own soc; mgaamc can meowumuop .mmspp

 

 

mm». ._.uop. ..mm~..oqm 4-44 .ufi. mH~.¢~H ...um~JHm .neumaumo mnmflqofi ma<
mm» :uon mad o.m m.mfi mm. umH.m- zm¢.om omnofiusfi eumfi.ma ma<
mm» uop omm m.m o.mH mm zom.mmfi zm¢.om .meummuma ema~.- cm:
mm» uop Neg e.m m.m~ we zmo.m- zo~.om nonmaumfi H~m~.o~ own
was uop emfi ~.m o.m~ m zom.mufi zmm.om NmuHmumo Han.m~ L52
o: oo_ m¢~ ~.m m.m~ mm 3 w.o- zmm.om moumouma oumfl.m poo
o: ass moH o.m u--- we mem.-~ zm~.~m mmnmoumo o~m~.e~ xpza
mm» nan mom m.m ---- m mmm.mm~ . zem.Hm mmummumm oumH.mH co:
o: uop mmH m.m ~.m~ m :~o.m- zmo.om emnomnmfi omafi.m~ cm:
mm» :uon “mm 0.0 m.¢H N.“ 3mm.m- zmH.om nmnnonuo oum~.u~ 3mm
mm» nan mm“ m.m m.¢~ mm umc.mnfi zmH.~m mfiumeuwo mom~.e ca<
o: guoa mfifl o.m m.~H wm, woe.nua zum.om maummufia momfl.m~ poo
mm» cop no“ o.m u--- on 3 m.o- z ~.om mmuemumo Noofi.om. Aw:
mm» ucmucmumcoc in- ~.m m.- we mmm.mufi z~o.om cmuumumo mmmH.~ «can
mm» :uoa mum «.0 c.5H m www.mma zNo.om Noummnmo mom~.~ poo
o: span av“ m.m m.~H om www.mhfi zmo.om qfiueeueo momH.m~ mcza
mm» :uon mam m.m m.m~ om mmm.~nfi zmm.om mounmumo momfi.om co:
mm» asm Hm” m.m u--- me z¢¢.mufi. zem.Hm ofiunflumo mmm~.n no;
o: cuoa fimfi ¢.m m.m~ mm www.mufi z-.om mfinmeumfi com~.m~ uamm
~:o_uumccou .zumz. _..um: .mno a. az.i._mzm.n. Ex.umm.n .uuauwmcog. muauwumg_ .mewh mama

 

,>o=hm wmzp 2H cum: mmx<=ozpa<u ”H ~3m<h

51

 

TABLE 2
WWSSN Stations in Networks -- locations in degrees/minuites/seconds
Code Location Latitude Longitude

 

Location Network:

 

 

COP COpenhagen, Denmark 55°41' N 12°26' E
CTA Charters Towers, Queensland, Australia 20 OS 18 S 146 15 16 E
DAG Danmarks Havn, Greenland 76 46 12 N 18 46 12 E
ESK Eskdalemuir, Scotland, U.K. 55 19 00 N 3 12 08 W
GDH Godhavn, Greenland 69 15 N 53 32 W
KBS Kingsbay, Svalbard 78 55 03 N 11 55 26 E
KEV Kevo, Finland 69 45 19 N 27 00 24 E
KIP Kipapa, Hawaii, USA 21 25 24 N 158 00 54 W
KON Kongsberg, Norway 59 38 57 N 9 35 54 E
LAB Lahore, Pakistan 31 33 N 74 20 E
MAT Matsushiro, Honshu, Japan 36 32 30 N 138 12 24 E
NIL Nilore, Pakistan 33 39 00 N 73 15 06 E
NOR Nord, Greenland 81 36 N 16 41 W
STU Stuttgart, Germany 48 46 19 N 9 11 42 E
UME Umea, Sweeden 63 48 54 N 20 14 12 E
Both networks:

ALQ Albuquerque, New Mexico, USA 34 S6 33 N 106 27 27 W
DAL Dallas, Texas, USA 32 SO 46 N 96 47 02 W
FLO Florissant, Missouri, USA 38 48 06 N 90 22 12 W
LON Longmire, Washington, USA 46 45 00 N 121 48 36 W
Amplitude Network:

802 Bozeman, Montana, USA 45 36 N 118 38 W
COR Corvallis, Oregon, USA 44 35 09 N 123 18 12 W
FVM French Village, Missouri, USA 37 59 02 N 90 25 34 W
GOL Golden, Colorado, USA. 39 42 01 N 105 22 16 W
GSC Goldstone, California, USA 39 42 01 N 116 48 17 W
JCT Junction City, Texas, USA 30 28 46 N 99 48 08 W
MNN Minneapolis, Minnesota, USA 44 54 52 N 93 11 24 W
MSO Missoula, Montana, USA 46 49 as N 113 56 26 w
TUC Tucson, Arizona, USA 32 18 35 N 110 46 56 W

data from ISO Catalog

52

the earthquake and the slab is increased, the amplitudes in the shadow
zone are moving toward their normal nO-slab values. In addition to
this effect, the more obvious shrinking of the shadow zone with
increasing 6 can be seen on figures 6A - 6L.

The high amplitude band is relatively thin and well-defined at
6 = 12.4°. As 6 increases, this band broadens, becoming relatively
unfocused. The expected effect of this would be that this high amp-
litude band would follow the behavior of its low amplitude counterpart.
As 6 increases, the high amplitude zone indicated by the area of high
ray emergence point density should disperse and become less apparent.
The shape of the high amplitude zone also changes, contracting toward
the seismic source. It should be noticed that while the band of high
amplitudes becomes less sharp beyond 6 = 12.4, it can not become any
sharper as 6 decreases beyond this point; the ray emergence points
which make up the high amplitude zone are essentially all along a
single line already (figure 6A).

Seismic data was checked to see if these theoretical amplitude
effects actually existed. WWSSN records for nineteen earthquakes in
the central and west-central Aleutians were examined (Table 1).
Depending on the location of the earthquake, one or both of two sets
of seismic stations were considered (Table 2). These two sets of seis-

mic stations were the "amplitude network" and the "location network".

The amplitude network was meant to check for the narrow band Of
unusually high amplitudes at the edge of the low amplitude "shadow
zone" (Sleep, 1973). It consists of WWSSN seismological stations in
the western United States, where this band was expected to fall. The

location network was meant to circle the source area as well as possible.

53

in order to collect travel time and waveform data for earthquake re-
location and depth phase studies. Due to the Arctic and Pacific Oceans
and their geographic relationship to the Aleutians, there are large
gaps in the desired circle. This situation is worsened by the fact
that what stations there are in the central and western Pacific

usually produce illegible records.

Amplitudes were Obtained from short period vertical component
records. The amplitude data recorded was processed tO remove apparent
amplitude differences based on station magnification. All amplitudes
were normalized to what they should have been if all the seismological
stations used a magnification of 100,000x. In addition, the data was
corrected to remove the amplitude effects of the focal mechanism when
possible, using the formulae for focal mechanism amplitude effects
given by Aki and Richards (1980). Some of the focal mechanisms used
were taken from Stauder (1968b) and the rest were derived from ISC
data. Some focal mechanisms were not well constrained by the data and
amplitudes for these earthquakes could not be completely corrected. A
focal mechanism amplitude correction requires one to know the position
of the ray going to a given station with respect to both focal planes.
_ Some earthquakes studied in this thesis are of relatively low magnitude
and the data did not seem sufficient to constrain both focal planes; in
this case an estimate based on the one known focal plane was made. In
some cases, for various reasons, the focal mechanism for a given
earthquake was unobtainable. These cases are indicated in Table 1
and on the relevant figures.

The corrected data was then normalized again, this time relative

to the amplitude of the strongest arrival recorded. This was given

54

FIGURE 8
RELATIVE AMPLITUDE VS. "ZERO RESIDUAL LINE"
Amplitudes are normalized so that the strongest is 10 and the

rest are in proportion to this value.

55

flop-u

«.Op 3

vs: .24. :33:

 

aha. .sudo.

1.6— In

1.0- an

'5.—

00¢:

6.6:.

 

upcuoo

 

o.N—Iu

$00. .00 30...:

 

050— .O— .0:-

FIGURE 8

56

poo-«00.500 ‘..¢I£0.‘ O:

woo—.3 .3:

on... .« 2....

 

80:00-500 '02-'50... GB

0.0:: 050— 0— 33¢...

 

0.0...

can. .v :3. .

.320230 [Se-cool 2.

0.5. a. I... .0“ £00

 

.6. coop .s .00.

as A
saw

0.

FIGURE 8 (cont.)

57

CO-«OOquO ecu—CQSOQE O:

V524

 

. comp .mp 0::—

co:oo:oo .EmEasooE on

ad. a 3.36 .80

O.

 

FIGURE 8 (cont.)

58

the arbitrary amplitude value of 10, with weaker signals becoming nines,
sights, and so on. All values which were not exactly an integer number
of tenths of the strongest amplitude were rounded to the nearest tenth
except for those which would have been rounded to zero. These were
arbitrarily set as amplitude 1 instead; therefore any station with an
amplitude listed as 0 did not recedve any noticeable signal at all.
These amplitude values were plotted on a polar graph of similar format
to that of figure 6, in order to see whether there was a band of high
amplitude and if so, how well it corresponded with the theoretical one.

Figure 8 shows the results of this examination. The curved line
on each figure indicates the line of zero travel time anomaly for that
value of 6 nearest the earthquake for which theoretical slab effects
were calculated. At this line and slightly to either side of it the
slab wave should arrive at the same time as the mantle wave. Being in
phase, the two should reinforce each other leading to higher amplitudes;
notice that this line coincides fairly well with the zone of high
ray emergence point density (Figure 6).

All values of 6 for these earthquakes are calculated using ISC
‘hypocentral locations. If these are in error due to a slab mislocation
effect, the theoretical high amplitude line should not quite pass
through the observed high amplitude area, if any.

At the smaller values of 6 represented, the correlation between
the "zero line" and the high amplitude area in the Observed data is
poor. For the earthquake of August 10, 1975 ( 6 = 12.8°), for
example, it looks as if there is a high amplitude band in the data,
but it falls inside the “zero line" for this distance. Moving out to

about 6 = 13.0 the correlation between data and "zero line“ becomes

59

good. For example, consider the earthquakes of June'2, 1966, April

4, 1969, and October 25, 1968. For all of these events the correlation
between zero line and high amplitudes is good, although not perfect.
This correlation seems to persist out to 6 = 13.2° or so, but beyond
this point the correlation between theoretical "zero line" and high
amplitude signals becomes very poor.

These results seem to support the idea that at least some forearc
and outer rise earthquakes may be mislocated due to slab effects. If
earthquakes closer to the slab (the ones placed by the ISC at 6 =
12.8°, for example) are mislocated, their actual location would be
farther from the slab. Their true “zero lines" would fall inside the
zero lines as presently plotted, and in that case the correlation
between zero line and high amplitude data would be good. At 6 = 13.0°
as indicated by the ISC, the slab mislocation effects have decreased,
and therefore the "zero lines" plotted on these figures already
correlate with areas of high amplitude.

At higher values Of 6, approxomately 6 = 13.4° and out, there seems
to be much scatter of high amplitude signals away from the zero line.
Many high amplitude stations do not fall on or near the line,
including some of the highest amplitude ones. This scatter is probably
associated with the fading of both the shadow zone and the high
amplitude band at its edge. Neither of these are as intense as they
were for earthquakes closer to the slab. The amplitude anomalies they
cause do not contrast as sharply with the general range of amplitudes
as before. The results are therefore more Open to distortion by
receiver structure or error; also possibly due to maladjustment of

the receivers themselves (Stein, personal communication).

60

2.4 Waveform Effects

It was proposed in Chapter 1 that a possible reason for the fact
that slab effects had not been observed for forearc and outer rise
earthquakes is that the arrivals from these earthquakes at seismological
stations in the shadow zone are systematically picked late. It was
thought that this tendency might be defeated by matching waveforms
of signals received at different stations and attempting to correlate
them. At the same time it was possible to check the accuracy Of ISC
hypocentral depth determinations by study of depth phases, if any.

Location network stations were used for this (Table 2). As with
the amplitude study, short period vertical component records were used
when available. Where short period records were unreadable, arrival
times could sometimes be obtained from the long period vertical records
for the same event. Where practical, amplitude network events were
used for this study as well.

Results of the attempt to check the depth of focus using depth
phases are given in Table 1. Where depth phases can be read, they seem
to indicate that these earthquakes fall in the depth range of 7 to 20
km below mean sea level, assuming that rocks in this area have seismic
velocities of around 6 km/sec. In general, focal depths indicated by
depth phases were shallower than those indicated in the ISC Bulletin,
whenever there was any significant difference between them. Some of
the ISC solutions were very shallow, e.g. about 3 - 5 km below msl, and
in these cases the depth phases seemed to indicate a deeper source.

The average depth as determined by the ISC was 26.3 km. Maximum was 48
km, minimum was 3 km, with a standard deviation of 14.6 km. Taking

6 km/sec as the average seismic velocity between the hypocenter and the

Distance from Pole (deg.)

61

 

 

 

 

 

 

'4'...
,—
0 :3-0
~00
-O
0).. r
1'" 031:9.
.'o
N .—
.1-
9
F—
' 1-
CD
" 1 1 1 L_ l L
o o O 0
o O O
,. N a)
(uni) uideo
FIGURE 9

ISC EARTHQUAKE DEPTHS VS. DEPTH PHASE EARTHQUAKE DEPTH DETERMINATIONS

62

FIGURE 10
"ESCALATING" WAVEFORM ARRIVALS AT SEVERAL SEISMIC OBSERVATORIES.

‘63

 

ESK

 

CTA

 

 

P pP DWD’

FIGURE 10

 

WW "2.12.0.

TUC

 

 

 

 

P pP pwps

September 25, 1964

64

FIGURE 11
"CLEAN" WAVEFORM ARRIVALS AT SEVERAL SEISMIC OBSERVATORIES.

 

65

STU

DAG

KBS

UME

 

 

P pP pwp
March 27, 1974

FIGURE 11

66

sediment-water interface, depth phases give an average hypocentral depth
of 17.6 km. Maximum depth using this method was 27 km, minimum was
13.5, with a standard deviation of 3.23 km. If the picks of depth
phases were correct, this would seem to indicate that seismicity in the
trench and outer rise is at a relatively uniform level. This would
also indicate that the depth of earthquakes in the forearc and outer
rise is not as great as had been generally believed, removing one
Objection to the single layer elastic-perfectly plastic rheology as a
model for the lithosphere (Forsyth, 1980). A comparison of ISC and
pwP depths is shown in Figure 9.
1 An attempt was made to determine if, as suggested earlier, earth-
quakes in the forearc and outer rise might be systematically picked
late at European seismological Observatories due to their low magnitude.
The method for doing this was suggested by the observation that seismic
signal waveforms on short period seismograms are relatively consistent
for a given earthquake no matter where received (Hong and Fujita,1981).
It was hOped that by comparing waveforms of emergent arrivals with the
waveform of a clear, impulsive arrival at a different station, it would
be possible to more precisely determine when the sighal actually
arrived at the station with the emergent signal. This was partially
successful. The waveform matching method does work. This is because
the short period waveform seems to be, aside from some high frequency
attenuation, amazingly consistent from station to station (Figures
10 and 11). This similarity seems to continue regardless of geographic
location or focal mechanism, and may indicate that most of the short
period waveform is determined by source structure and/or the prop-

agation of the fracture along the fault as the earthquake occurs.

67

However, the utility of this sort of "wave form correction" is
limited by the type of waveform involved. The two types of waveform
involved can be called "clean" and "escalating" waveforms. Figure 10
is a trace of an escalating waveform, while figure 13 is a clean wave-
form. In a clean signal the amplitude of the first arrival is almost
as large as that of the strongest wave in the wave train. The first
arrival may in fact be the largest-amplitude wave itself. In this case
a wave form correction is not called for, since if the largest swing of
the signal is seen the first swing will almost always be seen as well.
This was confirmed by some wave form matching done with this type
Of signal. On the other hand, the escalating signal calls for a wave
form correction in many cases. In this type of arrival, the signal has
low amplitude phases just before the first major impulse arrives. It
was suggested (Fujita, personal communication) that these "predecessors"
might be a wave which travelled in the subducting slab when the main
wave was a mantle wave, and this may be so for some stations and some
situations. However, for the earthquakes checked here, the relation—
ship between the arrival time of the weak arrival and that of the main
signal is constant, regardless of station location or azimuth from the
source. It would seem that a slab effect should be highly azimuth
dependent. Since this is not, it is probably not a slab effect.

These escalating signal earthquakes may be one of two things.
They may be a "slow“ earthquake, or they may be a multiple shock event
with the first shock being very much smaller than the second. There
may be other explanations, but these are the two most obvious.

In these escalating arrivals, the weak wave arrives a second or

two before the main signal. It is about a tenth the amplitude of the

68

main signal, more or less. It is easy to see that with this waveform
a low amplitude arrival may lead to the weak wave escaping notice.

"Wave form corrections" were applied to escalating signal earth-
quakes as needed. A copy was made of a clear example of the short
period signal, and its waveform was matched with that of the emergent
arrival in question to see where the weak wave would have been if it
could be seen. Usually, with this prOcedUre, the weak wave could be
found on the "emergent" signal as well as the strong one, but if it
was not known to be there it would probably escape notice. These weak
waves appeared under emergent conditions as a very small signal or in
some cases as a change in the character of the background noise.

It should be noted that the utility of these wave form corrections
is severely limited. They are a certain amount of trouble, and should
not be used if not really needed. They do not help clarify all earth-
quakes, only those with escalating arrivals. Finally, only escalating
waveform earthquakes within a certain range of magnitudes will benefit
from this correction. If the magnitude is too low the weak wave will
be seen almost nowhere and if the amplitude is too high it will be
seen almost everywhere. In either case the arrival time picks will be
consistent and no wave form correction will be needed. A reasonable
guess would be that an escalating earthquake of Mb 6.0 or greater
would show its weak_wave almost everywhere the earthquake was recorded,
but this depends strongly on the specific ratio of amplitudes between
weak wave and main signal. It is difficult to determine, even as a
guess, what the lower limit of the useful range for the wave form

correction would be.

69

FIGURE 12
RAY EMERGENCE POINTS CORRECTED FOR SLAB RAY PATH DISTORTION VS. THE SAME
EMERGENCE POINTS UNCORRECTED.

7O

“VG—005500

 

cozofioa 525.9%.
f; p .0“ mmmimowo

60—00550023

 

FIGURE 12

71

2.5 Effects on Focal Mechanism Plotting

Another theoretical effect of the subducting slab was noticed.
This is not an effect on seismic signals but on seismic practice. As
previously mentioned, at 6 = 12.4° rays leaving the source at a take-
off angle of 33° can emerge at epicentral distances of between 34°
and 85° from the source over about a 90° range in azimuth. According
to standard angle of incidence tables (e.g. Pho and Behe, 1972) a ray
arriving at a station 85° from the source should leave that source at
a take-Off angle of around 17.6°. If a focal mechanism is plotted
using standard practice, the assumption will be that the 17.6° take-off
angle is the correct one, when in fact the ray left the source with a
take-off angle of 33°. This can lead to more or less serious errors
in focal mechanism plots, the degree of error depending largely on the
distance between the earthquake and the slab. It is not clear whether
the theoretical ray paths for this study are sufficiently accurate
to allow correction Of focal mechanisms for this effect. In any case,
some focal mechanisms were plotted both with and without this theor-
etical ray path correction, and it is Observed that the "corrected"
mechanisms seem to show less intermixing of compressional and

dilitational arrivals compared to uncorrected ones (Figure 12).

2.6 The Practical Limits of Slab Mislocation Effects in the Aleutians
The problem remains that, although theoretically present, slab
effects have not been detected for forearc and outer rise earthquakes.
Examination of arc normal slab effects diagrams (Figure 6) shows that
earthquakes across most of the study area should show noticeable slab
effects. However, these slab effects will not be generally noticed

unless a statistically significant percentage of seismological stations

72

 

 

February 27. 1970 Mb 6.0. 327 Observations
FIGURE 13

STATIONS DETECTING THE EARTHQUAKE OF FEB. 27, 1970, SHOWING WHICH
ONES WOULD HAVE SHOWN SLAB EFFECTS FOR DIFFERENT VALUES OF 6.

73

 

 

 

May 30. 1967 Mb 5.0. 107 observations

FIGURE 14

STATIONS DETECTING THE EARTHQUAKE OF MAY 30, 1967, SHOWING WHICH ONES
WOULD HAVE SHOWN SLAB EFFECTS FOR DIFFERENT VALUES OF 6.

74

are affected by them for a given earthquake. Using the arc normal slab
effects diagrams, it is possible to investigate this question.

Figures 13 and 14 show the location of all ISC seismological
stations which reported receiving signals from two of the earthquakes
studied in this thesis: that of February 27, 1970, and that of May 30,
1967 respectively. On both plots the earthquake epicenter is at the
origin. It can be seen that the question of statistical significance'
of slab effects is complicated by uneven distribution of seismometers.
Large empty areas exist north of the epicenters, due to the Arctic Ocean
and Bering Sea; an even larger empty area is the Pacific Ocean, south
of the epicenters. From the distribution of seismological stations
shown on these figures, it appears that for Aleutian earthquakes slab
effects will be detected at a statistically significant number of
stations if the area of slab effects extends into North America or
Europe.

Also shown on Figures 13 and 14 are lines marking the limits of
slab effects for different values of 6. These lines are the lines of
zero travel time anomalies; they mark the absolute limit of possible
slab effects. Under actual operating conditions the areas in which
travel time effects would be noticed would be smaller, as arrival times
can be determined only to within about 0.3 seconds and any deviation
smaller than this could not be detected.

At 6 . 12.4 most of the reporting seismological staions for both
earthquakes could theoretically be affected by slab effects. The -O.3
second travel time residual contour line extends well into northern
Europe. Stations farther north, in Scandinavia and Greenland, should

show even greater travel time effects. As regards amplitude anomalies,

75

most stations in northern Europe or farther north fall into a
relatively pronounced “shadow zone" (Sleep, 1973) and so should show
noticeable amplitude loss as well.

The trench axis in the Aleutians falls around 6 = 13.0°. Moving
beyond this to 6 = 13.2°, the travel time residual at the distance of
Stuttgart has been reduced essentially to zero, and that near Copen-
hagen and Eskdalemuir has fallen to about 0.3 seconds, about the limit
of readability and just a tenth of what it was at 6 = 12.4°. At this
level it is doubtful that the residual would be noticed. Travel time
effects and amplitude effects should still be noticed in the
Scandinavian and Greenland areas, although at this point it becomes -
doubtful whether the number of stations affected is still statistically
significant. Moving out onto the outer rise, at 6 = 13.4° the shadow
zone has retreated so that no slab effects should be detected south of
northern Finland and Greenland. At this point slab effects would
probably go unnoticed simply because no teleseismic stations would be
in a position to record them. Beyond this point the slab effects,
although still seen on the arc normal slab effects diagrams, rapidly
decline to insignificance.

Figure 15 shows the maximum possible percentages of stations
detecting slab effects as a function of 6. It can be seen that the
statistical significance of the number of stations affected drops
rapidly with increasing 6. Slab effects cause errors in teleseismic-
ally determined hypocentral locations, but as the number of affected
staions decreases the probable mislocation decreases as well. It
would seem that for earthquakes much seaward of the trench the

teleseismic hypocentral location would be fairly accurate. There are

76

 

 

HDC>'-'

~\

\“
\s‘ ---------- Feb. 27, 1970(327)
80 "_ ‘s
\\ ------ May 30, 1967(107)
“o \\
2 x
c: \
3 60 — \‘
“" \
<t \
\\
2 x
.2 x
.. _ \
£3 44) \\
(I) \
\\

s
o\° 20 — \\

O 1 1 1 ........ 1

l2.4 12.8 13.2 l3.6 14.0

Distance from Pole (deg)

FIGURE 15
PERCENTAGES OF STATIONS AFFECTED FOR DIFFERENT VALUES OF 6.

77

two problems with this, however. First, as shown in Figure 15,
earthquakes of smaller magnitude show a higher percentage of seismo-
logical stations in areas which should show slab effects. This is
because smaller earthquakes tend to register only at relatively close
stations, which then are more Often in slab effect areas than distant
stations are. Second, earthquakes are far more common in and around
the trench (aside from the shallow thrust zone) than elsewhere
(Forsyth, 1980) which puts them on a statistical borderline as far as
teleseismic mislocations are concerned.

In conclusion, it seems that in the Aleutians larger earthquakes
(Mb 6) at 6 = 13.0° or more should be relatively well located by tele-
seismic means without slab corrections. Inside this line, i.e. inside
the trench axis, slab corrections will become increasingly important.
Smaller earthquakes, however, tend to be mislocated more easily than
large ones and may require slab mislocation corrections well out onto
the outer rise. It should be noticed that this statistical conclusion

agrees with the conclusion drawn from the Observed amplitude data.

CONCLUSIONS
[1] Slab effects in the forearc and outer rise include amplitude
effects, travel time anomalies of up to -4 seconds, and incidentally,
distortion of ray paths sufficient to obscure focal mechanisms some-
what. All effects are less pronounced and occur over a smaller area
the farther the seismic source is from the slab. This makes these
effects less distinct and harder to detect as distance from the slab

increases.

78

[2] Theoretical ray tracing indicates that slab effects should be seen
for some earthquakes on the forearc and outer rise, being detected
mainly at stations in Europe and the western United States due to the
geography and geometry involved. Observed amplitude effects seem to
correlate well with theoretical ones. The correlation is somewhat
faulty closer to the slab, probably due to slab-induced mislocation

of the eathquakes. Correlation also worsens with increasing distance
beyond about 6 = 13.2°, probably due to the fact that slab effects
should weaken and become less distinct as the seismic source is moved
further from the slab.

[3] Statistical studies of the number of stations detecting slab
effects for earthquakes of different magnitudes and at different
distances from the slab indicate that for Mb = 6 or so, slab effect
mislocation corrections can probably be dispensed with beyond the
trench axis in the Aleutians, but for smaller events the slab correct-
ions may need to be carried out onto the outer rise. The rapid decline
in the number of stations affected with increasing distance from the
slab, along with the relatively small magnitude Of slab effects at
this distance, probably explain why these effects have not been

previously observed.

BIBLIOGRAPHY

79

BIBLIOGRAPHY

Aki, K., and Richards, P.G., 1980. [Quantitative Seismology: Theory
and Methods. W.H. Freeman, San Francisco. 932 p.

 

Baer, A.J., 1981. Geotherms, evolution of the lithosphere, and plate
tectonics: Tectonophysics, v. 72, p. 203 - 227.

Barazangi, M., and Isacks, B.L., 1979. A comparison of the spatial
distribution of mantle earthquakes determined from data provided
by local and by teleseismic networks for the Japan and Aleutian
arcs: Bull. Seismol. Soc. Am., v. 69, p. 1763 - 1770.

Beaumont, C., 1979. On rheological zonation of the lithosphere during
flexure; Tectonophysics, v. 59, p. 347 - 365.

Biswas, N. N., 1973. P-wave travel time anomalies: Aleutian - Alaska
region: Tectonophysics, v. 19, p. 361 - 367.

Bodine, J.H., Steckler, M.S., and Watts, A.B., 1981. Observations
of flexure and the rheology of the oceanic lithosphere: J.
Geophys. Res. v. 86, p. 3695 - 3707.

Caldwell, J.G., Haxby, W.F., Karig, D.E., and Turcotte, D.L., 1976.
On the applicability of a universal elastic trench profile: Earth
Planet. Sci. Lett., v. 31, p. 239 - 246.

Caldwell, J.G., Turcotte, D.L., Haxby, W.F., and Karig, D.E., 1977.
Thin elastic plate analysis of outer rises (abstract): In
Talwani, M., and Pitman, W.C. III (eds.), Island Arcs, DEEp Sea
Trenches and Back—Arc Basins. Am. Geophys. Union, Washington,
DIC. p. 467.

Chapple, W.M., and Forsyth, D.W., 1979. Earthquakes and bending of
plates at trenches: J. Geophys. Res., v. 84, p. 6729 - 6749.

Chen, A.T., Frolich, C., and Latham, G.V., 1981. The seismicity Of the
forearc marginal wedge (accretionary prism): submitted to J.
Geophys. Res.

Cleary, J., 1967. Azimuthal variation of the Longshot source term:
Earth Planet. Sci. Lett., v. 3, p. 29 - 37.

Davies, 0., and McKenzie, D.P., 1969. Seismic travel-time residuals and
plates: Geophys. J. R. astr. Soc., v. 18, p. 51 - 63.

80

Engdahl, E.R., 1973. Relocation of intermrdiate depth earthquakes in
the central Aleutians by seismic ray tracing: Nature Phys. Sci.,
v.245, p 23 - 25.

Engdahl, E.R., 1977. Seismicity and plate structure in the central
Aleutians. Ig_Talwani, M., and Pitman, W.C. III (eds.), Island
Arcs, deep Sea Trenches and Back-Arc Basins. Am. Geophys. Union,
Washington, 0.0. p. 259 - 271.

Engdahl, E.R., Sleep, N.H., and Lin, M.-T., 1977. Plate effects in
North Pacific subduction zones: Tectonophysics, v. 37, p. 95-116.

Forsyth, D.W., 1980. Comparison of mechanical models of the oceanic
lithosphere: J. Geophys. Res., v. 85, p. 6364 - 6368.

Fujita, K., 1979. Tectonics Of divergent and convergent plate margins:
Ph.D. thesis, Northwestern University, 300 p.

Fujita, K., Engdahl, E.R., and Sleep, N.H., 1981. Subduction zone
calibration and teleseismic relocation of thrust zone events in the
central Aleutian Islands: Bull. Seismol. Soc. Am., v. 71,

p. 1805 - 1828.‘

Fujita, K., and Kanamori, H., 1981. Double seismic zones and stresses
of intermediate depth earthquakes. Geophys. J. R. astr. Soc
v.66, p.131-156.

Fujita, K., Rogers, W.J. Jr., and Sleep, N.H., 1980. Enhanced arrivals
due to subducting slabs (abstract): Eos (Trans. Am. Geophys.
Union) v. 61, p. 1045.

Hanks, T.C., 1979. Deviatoric stresses and earthquake occurance at the
outer rise: J. Geophys. Res. v. 84, p. 2343 - 2347.

Herrin, E., and Taggart, J., 1968. Source bias in epicenter
determinations: Bull. Seismol. Soc. Am., v. 58, p. 1791 - 1796.

Hong, T.-L., and Fujita, K., 1981. Modelling of depth phases and
source processes of some central Aleutian earthquakes: Earth
Planet. Sci. Lett., v. 53, p. 333 - 342.

Isacks, B.L., and Barazangi, M., 1977. Geometry of Benioff zones:
lateral segmentation and downwards bending of the subducted
lithosphere. Ig_Talwani, M., and Pitman, W.C. III (eds.),
Island Arcs, Deep Sea Trenches and Back-Arc Basins. Am.
Geophys. Union, Washington, 0.6., p. 99 - 114.

Isacks, 8., Oliver, J., and Sykes, L.R., 1968. Seismology and the
new global tectonics: J. Geophys Res., v. 73, p. 5855 - 5899.

Jacob, K.H., 1970. Three-dimensional seismic ray tracing in a
laterally heterogeneous spherical earth: J. GeOphys. Res., v. 75,
p. 6675 - 6689.

81

Jacob, K.H., 1972. Global tectonic implications of anomalous seismic
P-travel times from the nuclear explosion Longshot: J. Geophys.
Res.., v. 77, p. 2556 - 2573.

Jin, D.J., and Herrin, E., 1980. Surface wave studies of the Bering

S a ° '
2e44and Alaska area. Bull. Seismol. Soc. Am., v. 70, p. 2117 -

Johnson, A., 1970. Physical Processes in Geology. Freeman, Cooper,
and Co., San Francisco. 577p.

Jordan, T.H., 1977. Lithospheric slab penetration into the lower mantle
beneath the Sea of Okhotsk: J. Geophys., v. 43, p. 473 - 496.

Julian, B.R., 1970. Ray tracing in arbirtarily heterogeneous media:
Lincoln Lab. Tech. Note 1970-45.

Julian, B.R., and Gubbins, D., 1975. Two-point seismic ray tracing: a
comparison of methods (abstract): Eos (Trans Am. Geophys. Union)
v. 56, p. 1027.

Kirby, S.H., 1970. Tectonic stress in the lithosphere: constrants
provided by the experimental deformation of rocks: J. Geophys.
Res., v. 85, p. 6353 - 6363.

Lambeck, K., 1980. Estimates of stress differences in the crust from
isostatic considerations: J. Geophys. Res., v. 85, p. 6367 - 6402.

Liu, H.-P., 1980. The structure of the Kurile Trench - Hokkaido Rise
system computed by an elastic time-dependent plastic plate model
incorporating rock deformation data: J. Geophys. Res., v. 85,
p. 901 - 912.

McAdOO, D.C., Caldwell, J.G., and Turcotte, D.L., 1978. On the
elastic-perfectly plastic bending of the lithosphere under
generalized loading with application to the Kuril Trench:
GeOphys. J. R. astr. Soc., v. 54, p. 11 - 26.

Murdock, J.N., 1969. Short-term seismic activity in the central
Aleutian region: Bull. Seismol. Soc. Am., v. 59, p. 789 - 797.

Pho, H.-T., and Behe, L., 1972. Extended distances and angles of
incidence of P waves: Bull. Seismol. Soc. Am., v. 62, p. 885 - 902.

Pitman, W.C. III, Larson, R.L., and Herron, E.M., 1974. The age of the
ocean basins: Geol. Soc. Am. map.

Richardson, R.M., and Solomon, S.C., 1977. Apparent stress and stress
drop for intraplate earthquakes and tectonic stress in the plates:
Pure Appl. Geophys., v. 115, p. 317 - 331.

82

Richardson, R.M., Solomon, S.C., and Sleep, N.H., 1976. Intraplate
stress as an indicator of plate tectonic driving forces: J.
GeOphys. Res. , v. 81, p. 1847 — 1856.

Richardson, R.M., Solomon, S.C., and Sleep, N.H., 1977. Tectonic
stress in the plates: Rev. Geophys. Space Phys., v. 17, p. 981 -
1019.

Segawa, J., and Tomoda, Y., 1976. Gravity measurements near Japan
and study of the upper mantle beneath the oceanic trench-marginal
sea transition zones. Ig_$utton, G.H., Manghnani, M.H., and
Moberly, R., editors: The Geophysics of the Pacific Ocean Basin
and its Margin (A.G.U. Monograph # 19), American Geophysical
Union, Washington, D.C., p. 35-52.

Sleep, N.H., 1973. Teleseismic P-wave transmission through slabs:
Bull. Seismol. Soc. Am., v. 63, p. 1349 - 1373.

Solomon, S.C., Richardson, R.M., and Bergman, E.A., 1980. Tectonic
stress: models and magnitudes: J. GeOphys. Res., v. 85,
p. 6086 - 6092.

Sorrells, G.G., Crowley, J.G., and Veith,-K.F., 1971. Methods for
computing ray paths in complex geological structures. Bull.
Seismol. Soc. Am., v. 61, p. 27 - 53.

Spencer, C.P., and Engdahl, E.R., 1981. Inversion for hypocenters and
structure beneath the central Aleutian island arc (abstract):
Eos (Trans. Am. Geophys. Union) v. 62, p. 335.

Stauder, W., 1968a. Mechanism of the Rat Island earthquake sequence of
February 4, 1965, with relation to island arcs and sea floor
spreading: J. Geophys. Res. v. 73, p. 3847 - 3858.

Stauder, W., 19686. Tensional character of earthquake foci beneath
the Aleutian trench with relation to sea-floor spreading: J.
Geophys. Res., v. 73, p. 7693 - 7701.

Stein, 5., 1978. An earthquake swarm on the Chagos-Laccadive ridge and
its tectonic implications: Geophys J. R. astr. Soc., v.55,
p. 577 - 588.

Stein, 5., 1979. Intraplate seismicity on bathymetric features - the
1968 Emperor Trough earthquake: Jour. GeOphys. Res., v. 84,
p. 4763 - 4768.

Suyehiro, K., and Sacks, I.S., 1979. P- and S-wave velocity anomalies
associated with the subducting lithosphere determined from travel-
time residuals in the Japan region: Bull. Seismol. Soc. Am.,

v. 69, p. 97 - 114.

Sykes, L.R., 1966. The seismicity and deep structure of island arcs:
J. Geophys. Res. v. 71, p. 2981 - 3006.

83

Tokst, M.N., Sleep, N.H., and Smith, A.T., 1973. Evolution of the
downgoing lithosphere and the mechanism Of deep focus earthquakes:
GeOphys. J. R. astr. Soc., v. 35, p. 285 - 310.

Umino, H., and Hasegawa, A., 1975. On the two-layered structure of
deep seismic plane in northeastern Japan arc (in Japanese):
Zisin (ii), v. 27, p. 129 - 135.

Watts, A.B., Bodine, J.B., and Streckler, M.S., 1980. Observations of
flexure and the state of stress in the oceanic lithosphere:
J. Geophys. Res., v. 85, p. 6369 - 6376.

Yoshii, T., Kono, Y., and Ito, K., 1976. Thickness of the oceanic
lithosphere. 13 Sutton, G.N., Manghanani, M.H., Moberly, R.,
and McAffe, E.V. (eds.), The Geophysics of the Pacific Ocean
Basin and its Margins (GeOphys. Monograph 19). Am. Geophys.
Union, Washington, D.C. p. 423 - 430.

 

APPENDICES

APPENDIX I

PROGRAM "RAYTRAK"

PROGRAM RAYTRAK(INPUTyOUTPUT.TAPE4=INPUTpTAPES=OUTPUT.
+TAPE5.TAPE7:TAPEB)
C

C....ROUTINE TO TRACE RAYS THROUGH UPPER MANTLE VELOCITY STRUCTURE
C....BY JULIAN AND SLEEP (1973). MODIFIED FOR INTERACTIVE USE .
C....ON CDC 6800 BY FUJITA (197B). VERSION 4.1 31 OCTOBER. 1978 .
C....FREEFORM INPUT/OUTPUT. MODIFIED FOR BACKUP LOAD IN EVENT OF DUMP.

C....PROGRAM MODIFIED FOR PLDTTING OF RAY EMERGENCE POINTS USING
C....TEKTRONIX PLOTTING ROUTINEs AND THE TEKTRONIX MODEL 4010-1
C....GRAPHICS TERMINAL BY W.J.ROGERS JR..(1981).
DIMENSION PLDTAZ (400)
DIMENSION PLOTRAD (400)
DIMENSION TINE<Iooi 24(1001 zsrroo: PHIR(100).TDAR(100).AZRII
+00). THETAR(100). "(3). DELAYIIoo). TIMIIOO). AZR71100). 0(100)
DIMENSION KTYPEIS). KTYPEIIZ). ICOUNTS(2). RAYDATIS). FMT(10)
REAL,IDENT.LAT1.LAT2

DIMENSION 2(5)
COMMON /CMMDL/ R0.DR.NR.LAT1.DLAT.LAT2.NLAT.RADIUSIDRSGIDLSO.DRDL:

+LON1.LON2:H:DM2.NUMP.F
COMMON ICOMID/ ITHMAX

COMMON /CMRAY/ RFOCITHFOCIRNIN'DTIEPS'IRAY!IRAYPL!I"VFITHAX
COMMON /CMAMP/ 115'ISLOFFINPCPPIISIIPRAY
DATA KIYPE1/2.2/.KTYPE/1.1.1.1.1l'
DATA ICOUNTS/1y1/yRAYDAT/12.0.0.1:0.1:1.:20./
CALL INITT130)
CALL ANMODE
PI 8 3.1415926
IRAYPL = 0
PRINT 9250
9250 FDRMATI! TRACK RAYS>*)
READ(4:§I IPRAY
F 8 57.29579

c.... READ IN MODEL

10 CONTINUE
REHIND B
REWIND 5
NNN=O
READ (5.9030) KTABLE.NUMP.NIERMX:NRAD.ITHMAX:IRAYIKREADvIPUNCleDL
1T0:NOERR:NPCP:IIS
WRITE (8.9030) KTABLE:NUMP.NIERMX.NRAD:ITHMAX.IRAYyKREADyIPUNCH.ID
1LTO.NOERR:NPCP:IIB

84

85

WRITE (7.9030) KTABLE.NUMP.NIERMX.NRAD.ITHMAX.IRAY.KREAD.IPUNCH.ID
1LTO.NOERR.NPCP.I18

IDLTMX = IDLTO+KTABLE-1

READ (5.9020) EPS.DTEST.POLON.TMAX.DT.FLAT.TOAEPS

WRITE (8.9010) EPS.DTEST.POLON.TMAX.DT.FLAT.TOAEPS

WRITE (7.9010) EP5.DTEST.POLON.TMAX.DT.FLAT.TOAEPS

C
C.... INPUT‘OCCURS IN MODL5
C
CALL MODL5
IF (NRAD.ED.1) GO TO 20
C
C.... CONVERT TO RADIANS
C
TOAEPS = TOAEPS/F
20 CONTINUE
READ (5.9080) IRAD.JREAD.KSKIP
READ (5.9090) DMAX.DFOC.BOT
WRITE (8.9040) DFOC
C

C.... MAIN LOOP
C.... GET DELAY TABLE FOR ORDINARY RAY

.IF (KSKIP.NE.0) GO TO 50
IF (KREAD.E8.0) GO TO 30
READ (5.9050) FMT
READ (5.FMT) (TIM(K).K=1.KTABLE)
GO TO,‘0
30 CONTINUE.
40 CONTINUE
C
C.... A NEW EPICENTER AND SET OF RAYS IS CALCULATED ON EACH CYCLE OF LOO
C
50 CONTINUE
C
C.... READ IN EPICENTER PARAMETERS
C .
WRITE(8.9240)
9240 FORMAT(* ANG.FOCAL DX. FROM ARC CURVATURE CENTER FOR THIS SET?*/)
READ(4.*) SMDEL
80 CONTINUE
70 WRITE (8.9080)
READ(4.§) ITOA.IPHI
IGRAB = ITOAiIPHI
IF (IGRAB.LE.10) GO TO 80
WRITE (8.9070) IGRAB
READ (4.9210) SET
IF (SET.NE.3HYES) GO TO 70
80 WRITE (7.9080) IRAD.JREAD.KSKIP
WRITE (7.9010) DMAX.DFOC.BOT

WRITE (8.9100)
READ(4.*) DELTOA.DELPHI.TOASTR.AZRATS

TOASTR = TOASTR/F-
AZRATS = AZRATS/F

(Inn “000

86

THFOCD-= SMDEL

THFOC = SMDEL/F

DELTOA = DELTOA/F

DELPHI = DELPHI/F

BOT3 BOT/F ,

RMIN = RADIUS-OMAN

IF (DMAX.ED.0.) RMIN = 0.
RFOC = RADIUS-DFOC
ISLOFF = 0

BEGIN TO MAKE TABLE OF CALCULATED RAYS
N = JREAD+1

ITDA*IPHI+JREAD.LE.N74

‘ IF IIPHI*ITOA.EO.O) GO TO 130

CI...
COCO).

90

Cl...

C

100

110

120

130.

LOOP TO CALCULATE RAY PARAMETERS
INITIALIZE TAKE-OFF AZIMUTH

AZR(N) = AZRATS
I8 8 0
IO = IG+1

INITIALIZE TAKE-OFF ANGLE

.TOAR(N) = TOASTR

IT 3 0

IT 3 IT+1
GO TO 140
CONTINUE

IN 8 N+I

IF (NIER.GT.1) GO TO 120
TOAR(N) = TOAR(N-1)-DELTOA
AZR(N) = AZR(N-I)

IF (IT.LT.ITOA) GO TO 100
CONTINUE

AZR(N) = AZR(N-1)+DELPHI
IF (IG.LT.IPHI) GO TO 90
N = N-I

CONTINUE

GO TO 210

140 CONTINUE

C.... NOW READY TO CALCULATE A RAY

C

C

150

NIER = 0 '
CALL RAYS (TOAR(N).AZR(N).TIME(N).Z.IER

C.... CHECK IF ERROR OCCURRED DURING CALCULATION

C

87

IF (IER.EG.0.0R.NOERR.NE.O) GO TO 180
NIER 8 NIER+1
IF (NIER.GT.NIERMX) GO TO 190
C
C.... RAY WAS TRAPPED IN LOW VELOCITY ZONE AND DID NOT EMERGE
C
TOAR(N) 8 TOAR(N)+EPS
C
C.... VARY RAY TO GET IT UNTRAPPED
C
GO TO 150
180 CONTINUE .
2(1) 8 RADIUS-2(1)
C .
C.... CHECK IF RAY ENDED AT CORE BOUNDARY
C
IF (2(1).GT.FLAT) GO TO 180
C .
8.... CHECK IF RAY ENDED AT SURFACE PROPERLY
C
IF (ABS(Z(1)).LT.POLON) GO TO 170
C .
C.... RAY DID NOT EMERGE CORRECTLY
C
WRITE (8.9110) N.Z(1)
WRITE (7.9110) N.Z(1)
170 CONTINUE
C .
C.... SAVE RAY PARAMETERS FOR LATER OUTPUT
C
THETAR(N) 8 2(2)
PHIR(N) 8 2(3)
Z4(N) =_Z(4)
.IZSLN) 8.215)
GO TO 110
180 NIER 8 NIER+1
C
C.... RAY HAS HIT CORE
C
IF (NIER.GT.NIERMX) GO TO 200
C
C.... VARY RAY TO AVOID CORE
C
TOAR(N) = TOAR(N)+BOT
GO TO 150
190 CONTINUE
WRITE (8.9120)
WRITE (7.9120)

C.... LOOP STEP ABANDONED ALL RAYS WERE TRAPPED

GO TO 210'
200 CONTINUE

 

('3

CI...
CI...

210

Cl...

C...-

CI...

220
230

c...-

U

240

. LOOP STEP ABANDONED RATS HIT

88

IT;

I".

WRITE (8.9130)
WRITE (7.9130)

LOOP STEP RAY CALCULATION FINISHED
PREPARE TO OUTPUT RAY TABLE

CONTINUE

INITIALIZE THE DISTANCE ROUTINE

0(1) a DIST1(O..THFOC.PHIR(1).THETAR(1)18F

CALCULATE DISTANCE AND AZIMUTH OF EACH RAY RELATIVE TO EVENT

DO 240 I8 8 1.N ,.
D(IG) 8 DIST(PHIR(IG).THETAR(IO)) m
AZR7(IG) 8 ASIN(SIN(PHIR(IG))ISIN(D(IG))*SIN(THETAR(IG)))
ZETA 8 DIST(AZR7(IG).D(IG))

GET CORRECT BRANCH 0F AZIMUTH
IF (AOSIZETA-THETARIIO)1.0T.DTESTi‘A2R7IIO) = PI-AzR7IIOi
CONvERT TO DEGREES

AZR7(IG) 8 AZR7(IO)*F
D(IG) 8 D(IG)*F

CALCULATE TRAVEL TIME DELAY RELATIVE TO TRAVEL TIME TABLE

IDLT 8 D(IG) . . .

IF (IDLT.LE.IDLT0.0R.IDLT.GE.IDLTMX) GO TO 220
.AF.8MD(IG)-IDLT.

IDLT 8 IDLT-IDLTO+1 . .

DELAY(IG) 8 TIME(IG)-(TIM(IDLT)*(1.-AF)+TIM(IDLT+1)*AF)
GO TO 230

SET DELAY T0 0.0 IF DELTA IS NOT PART OF TRAVEL TIME TABLE

DELAYlIO) . o.o
CONTINUE

CONvERT TO DEGREES

TOAR(IO) 8 TOAR(IG)*F
AZR(IG) 8 AZR(IG)*F
THETAR(IG) 8 THETAR(IG)*F
PHIR(IG) 8 PHIR(IG)*F
24(18) 8 Z4(IG)*F

25(IG) 8 25(IG)*F
CONTINUE

Cl...
C

.,250

270
550
551

555
558

557
5571

5572'
559

559
580
581

280

444

445

89

OUTPUT RAY PARAMETER TABLE

'WRITE (7.9140) DFOC.THFOCD

WRITE (7.9150) N

WRITE (7.9170)

WRITE (8.9180)

DO 250 K 8 1.N
WRITE (7.9180) (K.TDAR(K).AZR(K).THETAR(K).PHIR(K).TIME(K).Z4(K)
.25(K).D(K).AZR7(K).DELAY(K))

WRITE (8.9190) (K.TDAR(K).AZR(K).D(K).AZR7(K).DELAY(K))
CONTINUE

LLL8NNN+1 '

NNN8NNN+N _

DO 270 K8LLL.NNN

MMM8(K-LLL)+1

PLOTAZ(K)8AZR7(MMM)

PLOTRAD(K)8D(MMM)

CONTINUE _

WRITE (8.550). .

FORMAT (* DO YOU WANT TO PLOT NOW? ENTER N FOR NO.*)

READ (4.551) PLOTNO

FORMAT (A1)

IF(PLOTNO.EG.1HN) GO TO 581

WRITE (8.558) . .

FORMAT(* YOU ARE GOING INTO A PLOTTING ROUTINE. TO*)

WRITE (8.557) - . . . ._

FORMAT(¥ CHANGE THE PAGE AND PROCEED WITH THE PROGRAM YOU *)

WRITE (8.5571) . -

FORMAT(* MUST_ENTER THE LETTER 'Y'. THERE IS NO PROMPT*)

"WRITE (8.5572). A

FDRMAT(* FOR.THIS.WHEN YOU HAVE UNDERSTOOD THIS ENTER TY'*)

READ (4.558) FOOL

FORMAT(A1) .

IF (FOOL.EG.1HY) GO TO 559

GO TO 555

CONTINUE

CALL PLOTZ (PLOTAZ.PLOTRAD.SMDEL.NNN.LLL)

WRITE.(8.580) A

FORMAT (8 VERY GOOD! CONGRATULATIONS! WE PROCEED.*)

BOT 8 BOTiF

WRITE (8.9200)

READ (4.9210) SET.

IF (SET.EG.3HNO ) GO TO 280

GO TO 80

CONTINUE

WRITE (8.9220)

READ (4. 9230) SETZ

IF (SET2.EG.3HYES) GO TO 10

WRITE (8.444)

FORMAT (* CATALOG TAPE 7 BEFORE LOGGING OUT OR IT WILL BE LOST!)
WRITE (8.445)

FORMAT (* AND YOU WILL LOSE MOST OF YOUR INFORMATION WITH IT.*)

90

CALL FINITT(150..90.)
C
C.... END OF
C

MAIN LOOP

9010 FORMAT
9020LFORMAT
9030 FORMAT
9040 FORMAT
9050 FORMAT
9080 FORMAT
9070 FORMAT
9080 FORMAT
9090 FORMAT
9100 FORMAT
1./)
9110 FORMAT
9120 FORMAT
9130 FORMAT
9140 FORMAT
9150 FORMAT
9180 FORMAT
9170 FORMAT

1AR.2X.20HDELTA

9180 FORMAT
1F7.3)
9190,FORMAT

(8810.2)
(8810.2)

(5X.1515)

(10H DEPTH 8
(10A4)

(* {TAKEOFF ANGLES. i AZIMUTHS (SEPARATE BY COMMAS)*)
(21H YOU ARE CALCULATING .I5.20H RAYS. ARE YOU SURE?)
(15X.15.5X.2I5)

(2E10.2.50X.E10.2)

(49H DELTOA.DELPHI.TOASTR.AZRATZ (SEPARATE BY COMMAS)

.F8.2.3H KM)

(1H0.I5.18H BAD RAY DEPTH IS.F12.3)

(42H ROUTINE DOES NOT CONVERGE DUE TO TRAPPING)

(31H ERROR RAYS KEEP HITTING BOTTOM)

(18H0 DEPTH OF FOCUS. F12. 5. 10H KM THFOC.F12.5.4H DEG)

(1H .I5. 21H RAYS WERE CALCULATED. I)

(4H NO..4X. 3HTOA. 4X. SHAZR. 2X.5HDELTA.3X. 4HAZIM. 3X. SHRESID)
(4H NO..4X. 3HTOA. 4XL15HAZR THETA PHI. 4X. 4HTIME. 5X.7HI AZ
AZIM RESID)
(I4.2(1X.F8.2).2(1X.F5.1).1X.F7.2.2(1X.F5.1).2(1X.F8.2).IX.

(I4.2(1X.F8.2).2(1X.F8.2).1X.F7.3)

9200 FORMAT (49H DO YOU WANT ANOTHER SET OF ANGLES IN THIS MODEL?)
9210 FORMAT (A3)
9220 FORMAT (32H DO YOU WANT TO TRY A NEW MODEL?)
9230 FORMAT (A3)
END.
SUBROUTINE VDER (R.LAT.LON.V.DVDR.DVDTH.DVDPH.DVDRDR.DVDRDT.DVDTDT
1)
COMMON IBLYTHE/ VEL(54.41).V1(294)

C
C.... CALCULATES VELOCITY AND ITS SPATIAL DERIVATIVES AS FUNCTIONS OF
C.... POSITION IN LATERALLY HETEROGENEOUS EARTH MODEL
C
REAL LAT.LAT1.LAT2.LON.LON1.LONZ
DIMENSION VV(4)
DIMENSION AV(4). AD(4). AD2(4)
COMMON ICMMDL/ R0. DR. NR.LAT1. DLAT.LAT2. NLAT. RADIUS. DRSG. DLSG. DRDL.
ILON1.LON2. H. DMZ
COMMON ICMAMP/ I15.ISLOFF
C
C....
C

LATERALLY VARYING PART OF MODEL

COMMON ICMDL1/ DR1.R1.NR1.DRISO
Y 8 (R-RO)/DR
IR 8 Y
IF (ISLOFF.EG.1) GO TO 20
IF (IR.LE.0.0R.IR.GT.NR) GO TO 30
Y 8 Y-FLOAT(IR)
8 (LAT*LATI)/DLAT

10

20
30

CIIIO

40

91

ILAT 8 X
IF (ILAT.LE.O.DR.ILAT.GT.NLAT-3) GO TO 30
X 8 X-FLOAT(ILAT)

JL 8 ILAT-1
DO 10 K 8 1.4
J 8 JL+K

CALL SPLN3 (VEL(IR.J).Y.AV(K).AD(K).AD2(K))
CALL SPLN3 (AV.X.V.DVDTH.DVDTDT)
DVDTH 8 -DVDTH/DLAT
CALL SPLN3 (AD.X.DVDR.DVDRDT.DMY)
DVDR 8 DVDR/DR

ADVDPH 8 0.

GO TO 40
CONTINUE

V 8 0.

DVDR 8 0.
DVDTH 8 0.
DVDPH 8 0.
DVDRDR O.
DVDTDT 0.'
DVDRDT 0.
GO TO 90

HAS LONGITUDE EXCEEDED THE END OF THE ISLAND ARC

IF (LON.GT.LONI) OD T0 50
z = ((LON-LONlllHliaz
GO TO so

50 IF”(LON:LT.LON2) GO TO 90

80

z = ((LON-LON2)/H)**2
IF (2.6T.40.) DO TO 70
F = EXP(-Z)

GO TO 80

70 F 8 0.

80

CI...

CONTINUE

FP_= -DM2aZ*F
DVDPH 8 V8FP

V = ViF

DVDR 8 DVDRiF
DVDTH 8 DVDTHRF

SPHERICALLY SYMMETRIC PART OF MODEL

.DMI 8 (R-RI)/DR1

I 8 DMI

X DM1-FLOAT(I)

I MINO(I.NR1)

VV(I) 8 V1(1)

VV(2) 8 V1(I+1)

VV(3) 8 V1(I+2)

VV(4) 8 V1(I+3)

IF (I.LE.0) VV(I) 8 V1(1)
IF (I.LE.-1) VV(Z) 8 V1(1)
IF (I.LE.-2) VV(3) 8 V1(1)

c...-

CI...

CC...

CC...

92

IF (I.LE.-3) VV(4) 8 V1(1)

CALL SPLN3 (VV.(.V2.DVDR1.DVDR2)
V 8 V+V2

DVDR 8 DVDR+DVDR1IDRI

RETURN

END

SUBROUTINE SPLN3 (A.X.F.DER.DER2)
DIMENSION A(4)

XHI = X’I a
XPI 3 X+1.
DMI AII)-A(4)+3.*(A(3)‘A(2))

DMZ A(1)-2.*A(2)+A(3)

F 8 0.5*X*XM1*(A(1)-X*DM1)~XM1*XP1*A(2)+0.5*X*XP1*A(3)
DER 8 X*(DM2-(1.5*X-1.)iDM1)+0.5*(A(3)-A(1))

DERZ 8 DM2-(3.*X-1.)*DMI

RETURN

END

SUBROUTINE MODL5

INITIALIZE ARRAYS

DIMENSION KYPE(2)

REAL LATI.LAT2.LON1.LON2.LADV(2)

COMMON IBLYTHE/ V(54.41).V1(294)

DIMENSION T(100). BOT(10). TP(100.10). N(10). NN(10). DVDPH(10)
DIMENSION FMT(10), .

COMMON ICMMDL/ R0.DR.NR.LAT1.DLAT.LAT2.NLAT.RADIUS.DRSO.DLSO.DRDL.

1LON1.LON2.H.DM2.NUMP.F

COMMON [COMID/ I25G,

COMMON_/CNDL1/ DR1.R1.NR1.DRISG
DATA_KYPE[1.1/.LADV/75.54.-.0009/.NR/51/
IFIX 8 252

CREATE VELOCITY OR DENSITY DISTRIBUTION TO BE USED IN MAIN PROGRAM

SKIP IF THIS PART READ IN PREVIOUSLY
IF (IZSG.NE.0) GO TO 30

SET UP SPHERICALLY SYMMETRIC PART OF MODEL

READ (5.9040) RADIUS.DR1.NR1.)RAD.ITMP.IBACK.IGRAV
WRITE (8.9040) RADIUS.DR1.NR1.NRAD.ITMP.IBACK.IGRAV
WRITE (7.9040) RADIUS.DRI.NR1.NRAD.ITMP.IBACK.IGRAV
DRISG 8 DR1**2

R1 8 RADIUS-FLOAT(NR1-1)*DR1

READ IN SPHERICAL BASIC VELOCITY OR DENSITY DISTRIBUTION

READ (5.9030) FMT

READ (5.FMT) (V1(I).I81.NR1)
WRITE (8.9070) RADIUS.R1.DR1
WRITE (7.9070) RADIUS.R1.DR1
WRITE (8.9050)

READ (4.9080) SETS

93

IF (SET3.NE.3HYES) GO TO 10
WRITE (7.9080) (I.V1(I).I81.NR1)
10 CONTINUE

C.... SET DUMMY POINTS ABOVE SURFACE NEEDED BY NUMERICAL RAY CALCULATION

DO 20 I 8 1.3
20 V1(NR1+I) 8 V1(NR1)
C.... SET UP LATERALLY VARYING PART OF MODEL.
C.... LATITUDE IS RELATIVE TO POLE OF ISLAND ARC
C.... LONGITUDE IS RELATIVE TO POLE OF ARC AND EVENT

READ (5.9100) DR.DLAT.NLAT.LON1.LON2.H
WRITE (6.9100) DR.DLAT.NLAT.LON1.LON2.H
WRITE (7.9100) DR.DLAT.NLAT.LON1.LON2.H
WRITE,(6.9090)
READI4.*) NR
WRITE (3.9110)
READ(4.*) LAT1.DVDT
DMZ = 2.19
DRSO = DR**2
R0 ; RADIUS-FLOATINR-I).DR
LATZ = LAT1+FLOATINLAT-1)*DLAT
NRITE (8.9120) R0.DR.LAT1.DLAT.LAT2.DVDT
WRITE (7.9120) R0.DR.LAT1.DLAT.LAT2.DVDT
30 CONTINUE
C
C.... END OF SKIPPED REGION
C.... READ IN LATERALLY VARYING PART OF DENSITY 0R VELOCITY
C
READ (5.9030) FMT
IF (ITMP.EO.1) GO TO 50
C a
C.... VELOCITY OR DENSITY IS READ IN DIRECTLY
C
DO 40 II = 1.NR
I = NR-II+1
READ (5.FMT) (V(I.J).J=1.NLAT)
40 CONTINUE
c A
c.... END OF MODL5 INPUT FOR THIS BRANCH
C
GO TO 170
so CONTINUE
C
C.... READ IN TEMPERATURE MODEL AND COHPUTE VELOCITY OR DENSITY
c
READ (5.FMT) ((V(I.K).K=1.NLAT).I=I.NR)
C
C.... SKIPPED IF ALREADY READ IN
c , .
IF (IZSG.NE.0) GO TO 70
C

94

C.... READ IN PARAMETERS FOR P-T DEPENDENT PHASE CHANGES
C
IF (NUMP.E0.0) GO TO 70

C
C.... EXECUTE LOOP FOR EACH PHASE
C
00 80 NP 8 1.NUMP
C .
C.... READ PHASE PARAMETERS
C
READ (5.9130) N(NP).NN(NP).BOT(NP).DVDPH(NP)
WRITE (8.9130) N(NP).NN(NP).BOT(NP).DVDPH(NP)
WRITE (7.9130) N(NP).NN(NP).BOT(NP).DVDPH(NP)
C
C.... READ CLAPEYRON CURVE
C
READ (5.9030) FMT
READ (5.FMT)_(TP(J.NP).J81.NR)
WRITE (8.FMT) (TP(J.NP).J81.NR)
WRITE (7.FMT) (TP(J.NP).J81.NR)
C
C.... END OF LOOP TO READ PHASE PARAMETERS
C
80 CONTINUE
C

C.... END OF SKIPPED REGION
C.... END OF MODL5 INPUT

70 CONTINUE
C.... CONVERT MINEAR TEMPERATURE INTO SEISMIC VELOCITY OR DENSITY

DO 180 I 8 1.NR

B4 8 V(I.1)
DO 80 K 8 1.NLAT
J 8 K
C
°C.... REVERSE HORIZONTAL INDEX IN ARRAY
C
IF (IBACK.GE.1) J 8 NLAT+1-K
C
8.... TEMPORARILY SAVE TEMPERATURE
C
T(J) 8 V(I.K)
BO CONTINUE
DO 90 K 8 1.NLAT
C
C.... COMPUTE ANOMALY WITH RESPECT TO END POINT
C
V(I.K) 8 (T(K)-B4)*DVDT
IF (IGRAV.E0.0) GO TO 90
C
C.... CALCULATE DENSITY 8 TiRHO*ALPHA

90

Cl...

95

V(I.K) 8 V(I.K)+V1(I+IFIX)
CONTINUE
IF (NUMP.EG.0) GO TO 150
COMPUTE ANOMALIES DUE TO PHASE CHANGES
DO 140 L 8 1.NUMP
ONLY CHECK OF PHASE CHANGES WHERE NECESSARY

IF (I.LT.N(L).OR.I.GT.NN(L)) GO TO 140
DO 120 K 8 1.NLAT

FIND DIFFERRANCE BETWEEN TEMPERATURE AT POINT AND P-T CURVE
ATOP'é T(KiLTPII.L)°
FIND WHICH PHASES ARE PRESENT

IF (ATOP.GE.0) GO TO 100
IF (ATOP+BOT(L).LE.O.) GO TO 110

MIXED PHASE REGION

RI 8 -ATOP/BOT(L)
GO TO‘120

ONLY HIGH TEMPERATURE PHASE PRESENT

RI 8 0.0
GO TO 120

ONLY LOW TEMPERATURE PHASE PRESENT
RI = 1.0
mWAmmeLmPM$MEmmE
T(K) = RI
COMPUTE ANOMALOUS AMOINT OF PHASE PRESENT RELATIVE TO END POIN

B 8 T(1)
DO 130 K 8 1.NLAT

ADD PHASE EFFECT TO VELOCITY OR DENSITY

V(I.K) 8 V(I.K)+DVDPH(L)*(T(K)-B)
CONTINUE

END OF PHASE CHANGE CALCULATIONS

C

96

150 CONTINUE
180 CONTINUE

C.... END OF VELOCITY CONVERSION
C

nnn

man

2

7

C
9

170 CONTINUE
.... SET DUMMY POINTS ABOVE SURFACE NEEDED BY NUMERICAL RAY CALCULATION

DO 180 I 8 1.3
DO 180 J 8 1.NLAT
180 V(NR+I.J) 8 0.
IF (NRAD.EO.1) GO TO 190

.... CONVERT TO RADIANS

_LAT1
LAT2
LON1
LON2
DLAT

LATI/F
LATZ/F
LON1/F
LON2/F
DLAT/F.
190 DLSO DLAT882
DRDL DR*DLAT
WRITE (8.9140)
READ.(4.9150) SET4
IF (SET4.EO.3HNO ) GO TO 210
DO 200 I81.54
WRITE (7.9180) (V(I.J).J84.22)
00 CONTINUE
DO 721 I81.54
, WRITE(7.9180) (V(I.J).J823.40)
21 .CONTINUE.
210 CONTINUE
RETURN

030 FORMAT(20A4)

9040 FORMAT (2F10.2.8I10)

9050 FORMAT (34H DO YOU WANT VELOCITY PRINTED OUT?)

9080 FORMAT (A3)

9070 FORMAT (SHORADIUS 8.F8.1.5X.5HR1 8.F8.1.5X.5HDR1 8.F8.1/)

9080 FORMAT (1H .5(I4.F8.3)) .

9090 FORMAT (28H NUMBER VERTICAL GRID POINTS/35H 51 FOR LONG SLAB 41 PO
1R SHORT SLAB./)

9100 FORMAT (F10.5.15X.F10.5.15.10X.3F10.5)

9110 FORMAT (11H LAT1.DV/DT./)

9120 FORMAT (5H2R0 8.F8.1.5X.4HDR 8.F8.1.5X.8HLAT1 8.F7.2.5X.8HDLAT 8.F
18.2.I7H LATZ 8.F7.2.5X.7HDV/DT =.E15.8)

9130_FORMAT (2I5.2E10.2)

9140 FORMAT (44H DO YOU WANT VELOCITY STRUCTURE PRINTED OUT?)

9150 FORMAT (A3)

9180 FORMAT (1X.19F4.2)
END
REAL FUNCTION DERPAR (Y.X.A.M)

C
CO...
C

rurara
I
I

10

20

=CIIII
CI...

97

CALCULATES DY/DX AT CONSTANT A

IMENSION A(300): X(300). Y(300): M(3)

‘ M(1)

N(Z)

MIG)
A(I)*Y(J)-A(I)*Y(K)-A(K)*Y(J)-A(J)*Y(I)+A(J)*Y(K)+A(K)*Y(I)
P A(I)*X(J)-A(I)*X(K)-A(K)*X(J)-A(J)*X(I)+A(J)*X(K)+A(K)*X(I)
DERPAR 8 S/P

RETURN

END

SUBROUTINE RAYFIX (SPHI:STHETA:PHI:THETA)

D
I
K
J
S

lllllllll

THIS IS MADE NECESSARY BY THE INANE WAY THE AMOD WORKS

PHI PHI-SPHI

PHI - AMOD(PHIrPIXZ)

IF (PHI.LT.-PI) PHI 8 PHI+PIX2
IF (PHI.GT.PI) PHI 8 PHI-PIXZ
DPHI 8 PHI

PHI 8 PHI+SPHI

IF (DPHI,GT.PIHF) GO TO 10

IF (DPHI.LT.-PIHF) GO TO 20

GO TOP4O

PHI =“PHI+PIHF

GO TO 30.

CONTINUE

PHI 8 PHI-PIHF

THETA 8 -THETA

CONTINUE

DTHETA 8 STHETA-THETA

THETA 8 AMODIDTHETAyPI)

IF (THETA.LT.-PIHF) THETA 8 THETA+PI
IF (THETA.GT.PIHF) THETA 8 THETA-PI
THETA 8 STHETA-THETA

RETURN

ENTRY RFSET

PI 8 3.1415926

PIHF 8 PI*.5

PIXZ 8 PI*Z.

RETURN

END

SUBROUTINE DERIV (T:Z:DZDT)

I I!

DEFINES DIFFERENTIAL EQUATIONS FOR RAY TRACING THROUGH
THREE DIMENSIONAL STRUCTURE

DIMENSION 2(5). 0207(5)

CALL VDER (2(1):1.5707SS-Z(2)12(3):VrDVDRyDVDTH:DVDPH:DVDR2:DVDRDT

1vDVDT2)

ETA 8 Z(1)/V
ETAI 8 1./ETA

CG...
CI...
CC...

C...-
Cl...

C

CO...
CO...

CI...

98

SINTH 8 SIN(Z(2))
COSTH 8 COS(Z(2))
COTTH 8 COSTH/SINTH
SINIR 8 SIN(Z(4))
COSIR 8 COS(Z(4))
SINZT 8 SIN(Z(SI)
COSZT 8 COS(Z(5))
DZDT(1) 8 V*COSIR
DRDT 8 DZDT(1)

DM1 8 SINIR/ETA
DM3 8 SINZT/SINTH

DZDT(2) 8 DMIiCOSZT

DZDT(3) 8 DM1*DM3

DZDT(4) 8 SINIR*(DVDR-ETAI)-COSIR/Z(1I*(COSZT*DVDTH+DM3*DVDPH)
DZDT(5) 8 ((DVDTH/SINIR-V*SINIR*COTTH)*SINZT-DVDPH*COSZTISINIRISIN
1TH1/Z(1)

RETURN

END

CALCULATES DISTANCE ON SPHERE
DISTANCE IS BETWEEN PHII:THETA1 AND PH151THETA5
PHII:THETA1 ARE SAVED FOR FUTURE CALLS TO DIST

REAL FUNCTION DIST1 (PHIlyTHETAI:PHIS:THETA5)
COMMON /HARVEY/ STORE:CT:ST
THETA2.8.THETA5

PHIZ 8 PHIS

CT 8 COS(THETA1)

ST 8 SIN(THETA1)

STORE 8 PHI1 .

B 8_PHI1-PHIZ

D 8 COSITHETAZ)*CT+ST*SIN(THETA2)*COS(B)
DIST1 8 ACOS(D)

RETURN

END 7

REAL FUNCTION DIST (PH13pTHETA3)

CALCULATES DISTANCE 0N SPHERE
DISTANCE IS BETWEEN PHILTHETAI AND PHIS'THETAS

COMMON IHARVEY/ STORE:CT.ST

PHIZ 8 PHI3

THETAZ 8 THETAB

B 8 STORE-PHIZ

D 8 COS(THETAZ)*CT+ST*SIN(THETA2)*COS(B)
DIST 8 ACOSID)

RETURN

END

SUBROUTINE DIFSYS (F:N:H:X:Y:EPS:S:NEWH)

F IS THE NAME OF A SUBROUTINE CALLED BY *CALL F(X:Z:DZ)* WHICH
STORES IN THE VECTOR DZ THE N COMPONENTS OF THE DERIVATIVE
DZ/DX ACCORDING TO THE DIFFERENTIAL EQUATION WHICH IS BEING

99

C.... SOLVED. DZ/DX8F(X.Z)

C.... sz. AND DZ MUST BE OF REAL TYPE

C.... N IS THE ORDER OF THE SYSTEM OF DIFFERENTIAL EQUATIONS.

C.... N MUST BE NO GREATER THAN MAXORD WHICH IS SET BELOW

C.... ITEST HAS TO BE SET EQUAL TO ONE IN THE CALLING PROGRAM IN ORDER T
C.... TEST THE ERROR RETURN OPTION

C.... ERROR IS SET TO 1 IF THE FUNCTION DOES NOT CONVERGW

C.... H IS THE_BASIC STEP SIZE . . .

C.... X AND Y (VECTOR) ARE THE INITIAL VALUES

C.... EPS AND SIVECTOR) ARE THE ERROR BOUNDS.

C.... ABSIEPS) SHOULD BE NO SMALLER THAN 1.0E-OS.

C.... NENH IS A FLAG WHICH IS SET EQUAL TO .TRUE. IF THE STEP SIZE USED
C.... BY DIFSYS IS DIFFERENT FROM THE STEP SIZE H GIVEN IN THE

C.... PARAMETER LIST. NEWH IS SET EQUAL TO .FALSE. OTHERWISE.

INTEGER ERROR
LOGICAL NENM
LOGICAL EPSERR
LOGICAL KONV.BO.BH
INTEGER R.SR I
DIMENSION ((5). 5(5)
DIMENSION 0(7) .
DIMENSIDN_YA(5). YL(5). YM(5). DY(5). 02(5). DT(S.7)
DIMENSION YG(S.5). YH(8.5)
COMMON /STUPID/ NEWRAY.
COMMON IDIF/ ITEST.ERROR
DATA EPSERR/.FALSE./
DATA MAXORD/ISI
ICNT . o
ERROR = 0
IF (NEWRAY.EQ.0) GO TO 20
DO 10 12 = 1.5
DO )0 J2 = 1.7
DT(IZ.JZ) = 0.0
10 CONTINUE
NEWRAY = o
20 CONTINUE
XMAS = X+H
IF (N.GT.O.AND.N.LE.MAXORD) GO TO 30
NRITE (6.9010)
STOP
so E = ABS(EPS)
IF (E.GE.I.0E-OS) GO TO so
IF (EPSERR) GO TO 40
EPSERR = .TRUE.
WRITE (6.9020)
40 E = 1.0E-os
C
0.... EACH CALL OF DIFSYS PERFORMS ONE INTEGRATION STEP OF THE
C.... EOOATION DY/DX=F(X.Y) ACCORDING TO THE METHOD 0; R. BULRIRSCH AND
C.... J. STDER (Munggxscug MATHEMATIK. IN PRESS). THE STEP SIZE HILL
C.... DE LESS THAN OR EOOAL In H. THE PROGRAM TAKES THE FIRST OF THE
C.... NUMBERS M.H/2.N/4..... AS STEP SIZE FOR WHICH NO MORE THAN 9

C....
C....
C....
C....
C....
C....
C....
C....
C....
C....
C....
C....
C....
C....
C....
C....
C....
C....

50

SO
70

80

90

100

EXTRAPOLATION STEPS ARE NEEDED TO OBTAIN A SUFFICIENTLY ACCURATE
RESULT. IF THE STEP SIZE USED IS DIFFERENT THAN THE STEP SIZE
GIVEN IN THE PARAMETER LIST. THEN THE LOGICAL FLAG NEWH WILL BE
SET EQUAL TO .TRUE.. OTHERHISE IT WILL BE SET EQUAL TO .FALSE..

X AND Y ARE THE INITIAL VALUES FOR THE STEP TO BE COMPUTED. AFTER
LEAVING THE SUBROUTINE. THE ORIGINAL VALUES OF THE PARAMETERS X
AND Y WILL HAVE BEEN REPLACED BY X+H§ AND Y(X+H*). RESPECTIVELY.
WHERE Hi IS THE STEP SIZE ACTUALLY USED. IN ADDITION THE STEP SIZ
WILL HAVE BEEN CHANGED AUTOMATICALLY TO AN ESTIMATED OPTIMAL

STEP SIZE FOR THE NEXT INTEGRATION STEP. THE ARRAY S AND THE
CONSTANT EPS ARE USED TO CONTROL THE ACCURACY OF THE COMPUTED
VALUES. THE SUBROUTINE IS LEFT. IF ALL I81: 2: . ... N TWO
SUCCESIVE VALUES FOR YII) DIFFER AT MOST BY AN AMOUNT EPS*S(I) .
EPS SHOULD NOT BE SMALLER THAN 1.0E-OG. FOR THE FIRST INTEGRATION
STEP IT IS_ADVISABLE TO SET S(I)80.0. BEFORE RETURN TO THE .
CALLING PROGRAM. THE ARRAY S WILL HAVE HAD ITS CONTENTS MODIFIED _
SO THAT S(I)8MAX(S(I):ABS(Y(I.X)). WHERE THE MAXIMUM IS TAKEN OVER
THE INTEGRATION INTERVAL (X.X+H*).

CALL F (X'YIDZI
BH 8 .FALSE..
NEHH 8 .FALSE.
DO SO I 8 1.N
YA(I) 8 YII)
A 8 XtH
ICNT 8 ICNT+1
IF (ICNT.GT.10) GO TO 230
PC 8 1.5
BO 8 .FALSE.

J 8 Jl-l

0(2) 8 2.25

IF (BO) DIZ) 8 4.0/D(2)
0(4) - 4.0*D(2)

0(61 4.0*D(4)

KDNV J.GT.2

IF (J.LE.8) GO TO 80

L 8 6

0(7) 8 64.

FC 8 0.B*FC

GO TO 90

L 8 J

D(L+l) 8 M*M

M 28M

6 H/FLOAT(M)

B 2.0*G

IF (BH.AND.J.LT.8) GO TO 140
KK 8 (M82)/2

M 8 M-1

DO 100 I = 1.N

100

110

120
130

140

150
160

170
180

190
200

210

101

YL(I) YA(I)
YM(I) YA(I)+G*DZ(I)
IF (M.LE.0) GO TO 160
D0 130 K 8 1.M
.CALL F (X+G*FLOAT(K).YM.DY)
DO 110 I 8 1.N
U 8 YL(I)+B*DY(I)
YL(I) 8 YM(I)
YMII) 8 U
U 8 ABS(U)
SII) 8,AMAX1(U.S(I)) .
IF (K.HE.KK.OR.K.EQ.2) GO TO 130
JJ 8 JJ+1
DO 120 I 8 1.
YH(JJ+1.I)
YG(JJ+1.I)
CONTINUE
GO TO 160
D0 150 I 8 1.N
YM(I) 8 YHIJ+1.I)
YL(I) 8 YG(J+1:I)
CALL F (AvYMoDY)
DO 200 I 8 1.N
V 8 DT(I.1)
DTII.1) 8 0.5*(YM(I)+YL(I)+G*DY(I))
C 8 DTII:1I
TA 8 C
IF (L.LE.0) GO TO 190
DO 180 K 8 1.L
1 8 D(K+1)*V
BI-C
V
F (B.EQ.0.0) GO TO 170
(C-VIIB
C88
8188
DTII.K+1)
DT(I:K+1) 8 U
TA 8 U+TA
IF (ABS(Y(I)-TAI.GT.E*ABS(S(I)I) KONV 8 .FALSE.
Y(I) 8 TA
IF (“ON”) GO TO 220

YMII)
YLII)

ll "2

B
B
U
I
B
U
C
V

0(3) = 4.0
0(5) = 16.0
80 = .NOT.SO
M = R
R = SR
SR = 2*M
BH = .NOT.BH
NENN = .TRUE.

H 8 0.5*H

rlrlr)

250

9010
9020
9030
9040

Congo
Coo-l
c...-

C
CI...
C

C
CO...
C

102

THE FOLLOWING TEST IS TO PREVENT INFINITE LOOPONG
THE NUMBER 1.E-7 COULD BE CHANGED IF SO DESIRED

IF (ABS(H/XMAG).LT.1.0E-7) GO TO 290
GO TO 70,

H 8 FC*H

X 8 A

GO TO 250

IF (ITEST.NE.1) GO TO 240
ERROR 8 1

HRITE.(G.9030)

GO TO 250

WRITE (6.9090)

WRITE (6.9040)

STOP

RETURN

FORMAT (ZBHOORDER TOO LARGE FOR DIFSYS.)

FORMAT (SIHOERROR LIMIT TOO SMALL FOR DIFSYS. WE USE 1.0E-OG.)
FORMAT (1H0.4SHMESSAGE FROM DIFSYS.FUNCTION DOES NOT CONVERSE)
FORMAT (1H0.25HPROGRAM IS FORCED TO STOP)

END .

SUBROUTINE RAYS (TOA.AZ.TS.ZS.IER)

TRADES RAY THROUGH IND-DIMENSIONAL VELOCITY STRUCTURE
2(I)=R 6.. Z(2)8THETA H 2(3)8PHI 2(4)=I
2(5)=zETA Z(6)8DR/DIO Z(7)8DTHETA/DIO Z(8)8DPHI/DIO

REAL IDENT.LATI.LAT2
LOGICAL NENN
EXTERNAL DERIV

A NEW STEP TO OUTPUT RAY PATHS TO CARDS FOR LATER PLOTTING HAS BEE

DIMENSION 2(5). DZDT(5). 25(5). 5(5)

DIMENSION IVVV(100). JVVV(100)

DIMENSION CS(12). SN(IZ)

COMMON ICMMDL/ RO.DR.NR.LAT1.DLAT.LAT2.NLAT.RADIUS.DRSO.DLSO.DRDL.
1LON1.LON2.HH.DM2 .

COMMON /CMRAY/ RFOC.THFOC.RMIN.DT.EPS.IRAY.IRAYPL.IWVF.TMAX

COMMON ICMAMP/ IIS.ISLOFF.NPCP.I16.IPRAY

COMMON ISTUPID/ NEHRAY -
DATA CS.SN/.866..5.0..-.5.-.866.-1..-.866.-.5.0...5..866.1...5..86
16:1.IABSSP.510.v'.51'.9667’1.t‘.8669’.510o/

INITIALIZE PARAMETERS'

NEWRAY 8 1

M 8 0

JLAT 8 (NLAT-1)*10
WLAT 8 3.141593/2.-LAT1
IER 8 0

N 8 5

10
20

30
40

CI...
50

60

C...-

90

103

ITER 8,0

T 8 0.

H 8 OT

2(1) 8 RFOC

2(2) 8 THFOC

2(3) 8 0.

2(4) 8 3.141593-TOA

2(5) 3.141593-AZ
IF (IRAY.EQ.0) GO TO 10
CALL DERIV (T.Z.DZDT)
WRITE (3.9020) T.(Z(I).I81.N).(DZDT(I).I81.N)
CONTINUE
DO 20 I 8 1.N
9(1) 8 AMAX1(ABS(Z(I)).1.0)
DO 40 I 1.N
ZS(I) Z(I)
TS 8 T
GO TO 60

REUSE ORIGINAL TIME HHEN PLOTTING

H = Hz

.IF (H.LE.0.01) H = DT

ERP = EPS

OH = H

INTEGRATE ONE STEP

CALL DIFSYS (DERIU.N.H.TS.ZS.ERP.S.NEHH)
H2,= T+OH-TS

IF (ITER.NE.0) GO TO 130

CHECK FOR RAYS HITTING CORE

IF (ZS(1).LT.RMIN) GO TO 110

CHECK IF RAY HAS RETURNED TO SURFACE

IF (25(1).GT.RADIUS) GO TO 100

CHECK IF EXCESS TRAVEL TIME HAS OCCURRED NITHOUT RAY EMERGING

IF (TS.LT.TMAX) GO TO 70

SET ERROR FLAG

IER 8 1

RETURN

DO 80 I 31:"
Z(I) 8 ZS(I)

T = TS

IF (IRAY.EQ.0) GO TO 90

CALL DERIV (T.Z.DZDT)

WRITE (6.9020) T.(Z(I).I81.N).(DZDT(I).I81.N)
CONTINUE

104

IF (I1S.EQ.0) GO TO 60

C
C.... SAVE RAY PATH FOR LATER PUNCH OUT
C .
155 = (RADIUS-Z(I))IDRilo.
LGG 8 (“LAT-2(2))IDLAT*10.
IF (11S.LT.0) LCD 8 (JLAT-LOG)
M 8 M41
IVVV(M) 8 LGG+1
JVVV(M) 8 IGG+1
C .
C.... IS BUFFER FULL
C .
IF (M.GE.100) GO TO 140
GO TO 50
C I
C.... RAY ENDS AT SURFACE
C
100 RF 8 RADIUS
GO TO 120
C
C.... RAY ENDS AT BOTTOM OF MODEL
C
110 RF 8 RMIN
C
C.... DO NOT NASTE TIME WITH RAYS WHICH HIT CORE
C
IF (NPCP.EQ.0) GO TO 140
C .
C.... ITERATE TO FIND END OF RAY
C
120 ITER 8 1
CALL DERIV (T.Z.DZDT)
H 8 (RF-2(1))IDZDT(1)
GO TO 30 .
130 D 8 RF-ZS(1)
IF (ADS(D).LT.EPS*S(1)) GO TO 140
IF (ITER.GT.4) GO TO 140
ITER 8 ITER+1
CALL DERIV (TS.ZS.DZDT)
H 8 OH+DIDZDT(1)
GO TO 30
140 CONTINUE

IF(IPRAY.EQ.1) PRINT 9025. M
9025 FORMAT(1X.15) .
IF(IPRAY.EQ.1) PRINT 9030. ((IVVV(I).JVVV(I)).I81.M)
9030 FORMAT(1X.1515)
. RETURN
C
9010 FORMAT (1H1.20A4/13X.4HTIME.7X.1HR.5X.5HTHETA.4X.3HPHI.GX.4HZETA.5
1X.1H1.9X.5HDR/DT.BX.9HDTHETA/DT.4X.7HDPHI/DT.SX.BHDZETA/DT.5X.5HDI
2/DTl)
9020 FORMAT (IP9E12.4)
END

nnnnnn

101
110

105

SUBROUTINE PLOTZ (PLOTAZ.PLOTRAD.SMDEL.NNN.LLL)

THIS SUBROUTINE PLOTS THE DATA ON A TEKTRONIX TERMINAL.
ALMOST EXCLUSIVELY CANNED TEKTRONIX ROUTINES AND SO YOU SHOULD
REFER TO AN APPROPRIATE MANUAL IF YOU WANT TO UNDERSTAND ANY

OF THIS STUFF.

DIMENSION PLOTAZ (400)
DIMENSION PLOTRAD (400)

CALL ERASE

CALL DWINDO (0..100..0..190.)
CALL SWINDO (10.1013.10.1013)
CALL POLTRN (0..190..0.)

CALL MOVEA (100..0.)-

CALL DRAWA (100..180.)

CALL MOVEA (0..0.) .

CALL DRAWA (5..90.)

‘CALL MOVEA (100..0,)

CALL DRAHSA (100..180..14)
CALL.MOVEA (BQ..0.)IL
CALL DASHSA_(80..190..14)
CALL MOVEA (GO..0.) _
CALL DASHSA (SO..180..14)

.CALL.MOVEAN(40..0.)

CALL DASHSA (40..190..14)

CALL MOVEA (20.:0.)1 I .

CALL DASHSA (20..190..14)

CALL MOVEA (0..0.)“

CALL DASHA_(100..90..14)
CALL,MOVEAT(0..0.)
CALL.DASHA-(100.EGO..14)
CALL_MOVEA (0..0.)

CALL DASHA (100..90..14)

CALL MOVEA (0..0.) .

CALL DASHA,(100..120..14)

CALL MOVEA (0..0.)

CALL DASHA (100..150..14)

L8LLL-1

IF(L.EQ.0) GO TO 110

DO 101 I81.L
IF(PLOTAZ(I).GE.190.00) PLOTAZ(I)80.
CALL POINTA (PLOTRAD(I).PLOTAZ(I))
CALL MOVREL (-4.0)

CALL DRHREL (8.0)

CALL MOVREL (-4.-4)

CALL DRWREL (0.8)

CONTINUE

CONTINUE

DO 111 I8LLL.NNN
IF(PLOTAZ(I).GE.190.00) PLOTAZ(I)80.
CALL POINTA (PLOTRAD(I).PLOTAZ(I))
CALL MOVREL (-B.0)

CALL DRWREL (16.0)

(‘3
F.)

223
224
225

333

106

CALL MOVREL (-8.-B)

CALL DRHREL (0.15}
CONTINUE

CALL MOVEA (0..0.)

CALL MOVREL (0.740)

CALL ANMODE

NRITE (6.222) SMDEL
FORMAT (* ANG. DIST. 19*. 1FS.2)
WRITE (9.223) NNN

FORMAT (I5.* RAYS.*)

READ (4.225) FLAG

FORMAT (A1)
IF(FLAG.EQ.1HY) GO TO 333
GO TO 224

RETURN

END

APPENDIX 11

USING "RAYTRAK"

APPENDIX II

USING "RAYTRAK"

As mentioned in the text, Raytrak is a computer ray tracing
program using the ray tracing method developed by Julian (1970).
It is based on Fermat's principle of stationary time; that is, a
ray travels from one point to another via the path that takes the
minimum time (or the maximum, but that does not concern us here).
The program was first written by Julian and Sleep in 1973. It was
modified for interactive use by Fujita in 1978. In 1981 it was
modified by the present author to plot ray emergence points on a
Tektronix 4010-1 graphics terminal. These modifications make the
program much faster and easier to operate than it would be
otherwise.
First, a few warnings are in order.
1. The program calculates ray emergence points and stores
their locations for plotting. Since it plots all calculated
points each time (not just the ones most recently calculated)
these points can be numerous. The arrays for plot radius and
plot azimuth are set up for four hundred points, so you should
never calculate more points than this before ending the run
and either considering another distance from the slab, or
stopping completely. You will probably never calculate anything
near this number of points on a single plot, but you should know

about this limitation.
107

108

2. As set up, the program calculates all the rays first and

prints out the results afterward. This means that should the

time limit be reached before the end of the program, all your

results will be lost. If you are going to run a lot of points

it would be better to reset the time limit before you start

in order to avoid this problem.

3. The plot produced by the Tektronix is useful, but not

detailed. All of the detailed results are on a local file

called Tape 7. This should be cataloged at least temporarily

before you sign off, or it will be lost forever. Log on again

on some hard-copy terminal and list the file to get your

permanent capy.

When you run RAYTRAK, it sends you messages and prompts. You
should respond to these as follows.

TRACK RAYS>

If you respond to this message with the number 1, the computer
will print out the coordinates of the ray as it moves through the
velocity model. Its path can then be plotted. If you respond with
the number 0 these coordinates will not be printed out.

DO YOU WANT VELOCITY PRINTED OUT?

The usual answer for this is No, as the printout is extensive
otherwise.

NUMBER VERTICAL GRID POINTS

51 FOR LONG SLAB 41 FOR SHORT SLAB

This is more or less self explanitory. The long slab is 350 km
deep at the maximum, and this is the model used for this thesis. To

get this model, answer with the number 51.

109

DO YOU WANT VELOCITY STRUCTURE PRINTED OUT?

If you answer "yes", you will get extensive printout showing the
velocity model you are using. You don't usually need this, so your
usual answer should be "no".

ANGULAR FOCAL DX. FROM ARC CURVATURE CENTER FOR THIS SET?

This is a distance measured in geocentric angular distance, as in
the text of this thesis. Put in whatever distance you are interested
in, out to 14.0. The trench falls at about 12.4, the outer rise at
13.0 and the volcanic front at 11.4 in the Aleutians.

# TAKEOFF ANGLES, #AZUMUTHS

(SEPARATE BY COMMAS)

These are the number of takeoff angles and azimuths for which
rays will be calculated. The machine takes a given azimuth and keeping
that fixed calculates the results for however many takeoff angles you
order, then goes on to the next azimuth. Takeoff angle becomes
smaller for each iteration (which makes the rays generally emerge
farther from the source). Azimuth moves away from the direct downdip
direction with each iteration. How far they move each time is
decided by your response to a later question.

YOU ARE CALCULATING ____ RAYS, ARE YOU SURE?

This message only appears when you order more than 10 rays. It
serves as a safeguard. If you respond with "yes" you get the rays
calculated for you. If you answer "no" you go back to the previous
question and have to choose the number of take-off angles and azimuths

again.

110

DELTOA, DELPHI, TOASTR, AZRATZ

(SEPARATE BY COMMAS)

Respond with four floating point numbers separated from each
other by commas. Deltoa is the step in take-off angle between adjacent
rays, delphi is the step in azimuth between them, toastr is starting
take-off angle, and azratz is starting azimuth. Given, say, the
response

2., 5., 32., 0.
the computer will send the first ray off at a take-off angle of 32°
and an azimuth with respect to the downdip direction of the slab
of 0° (or in other words, directly down the slab). If you chose
to order more than one take-off angle earlier, the second ray
calculated will be of the same azimuth (0.) but will have a
take-off angle 2° SMALLER than the first one. The computer will
calculate rays in this manner until it has done as many take-off angles
as you ordered. Then, if you ordered more than one azimuth, it will
switch to an azimuth 5° greater than the first one and start over with
the original take-off angle of 32°.

If you should want to send a single ray, respond to the question
#TAKEOFF ANGLES, # AZIMUTHS with 1,1. Then when the computer asks
for DELTOA, DELPHI, TOASTR, AZRATZ you can put in anything you wish
for deltoa and delphi. In this case they are meaningless. Put in
the takeoff angle and azimuth you want for toastr and azratz.

Hopefully you will start to get results at this point. You
will probably first get a lot of messages like
19 BAD RAY DEPTH 15 12.9

which means, in the example above, that ray number 19 did not converge

111

to the surface as it was supposed to. Instead it converged to a
depth of 12.9 km. You have to correct your travel times for this
error by hand. Nevertheless, this is not a fatal problem and if the
"bad ray depth" is small the corrected solution is still fairly
accurate. Large "bad ray depths" caUse serious errors in travel time
calculations and ray emergence point location as well.

If you get

ROUTINE DOES NOT CONVERGE DUE TO TRAPPING
or

ERROR RAYS KEEP HITTING BOTTOM
you're in trouble. Either the rays are trapped in some layer and never
reach the surface or they are hitting the Earth's core. In either case
the set is useless to you, especially since the computer abandons the
calculation at this point and shuts your program down. This forces
you to start at the beginning again.

Assuming you avoid this problem, you get a list of output. You
will lose much of this on the Tektronix screen due to over-writing and
when you "change pages" to clear the screen. Don't worry about this.
The printout on your screen is just to tell you that the program is
working properly. There is a copy in the local file Tape 7 .

Your next message is

DO YOU WANT TO PLOT NOW? ENTER N FOR NO

The no message is specified because if you respond with anything
else the computer goes ahead and plots your data. Therefore enter “N"
for no, and anything else for yes. A "no" answer skips you over the
plotting routines. A "yes" answer results in your getting the following

message.

112

YOU ARE GOING INTO A PLOTTING ROUTINE. TO CHANGE THE PAGE

AND PROCEED WITH THE PROGRAM YOU MUST ENTER THE LETTER "Y".

THERE IS NO PROMPT FOR THIS. WHEN YOU HAVE UNDERSTOOO THIS

ENTER "Y".

This refers to the end of the plotting routine itself, after the
plot is finished. To avoid writing messages all over your plot, the
program is set up to just sit there and do nothing until it gets that
"y" it is waiting for. That is the only way to get out of the plotting
routine, and you can punch anything else you can think of with no
effect. This is very frustrating if you happen to forget that you're
still in the plotting routine.

So, having read the message above, enter "y". The computer will
plot a pretty picture on the screen of your terminal. When it has
finished, and when you have made a hard copy of the plot if you wish
(and if you can), enter "y" again to exit the plot routine. The
computer responds with

VERY GOOD! CONGRATULATIONS! WE PROCEED.

DO YOU WANT ANOTHER SET OF ANGLES IN THIS MODEL?

If you answer "yes", you return to the point of entering the
number of take-off angles and azimuths desired. The distance from
the slab and the slab model itself remain unchanged.

If you answer "no", you get

DO YOU WANT TO TRY A NEW MODEL?

If you answer "yes" you go back to the point of setting slab

length and so on. You proceed from there as before.

113

If you answer "no“ you get

CATALOG TAPE 7 BEFORE LOGGING OUT OR IT WILL BE LOST AND YOU

WILL LOSE MOST OF YOUR INFORMATION WITH IT.
This means what it says and I suggest you take heed of it.

The program ends immediately after this, automatically. Make
a hard COpy of tape 7 if you wish. If you want, you can use this to
make a plot of your ray emergence points by hand; you can do a
slightly better job than the Tektronix can. Other than that, you're

through.

MICHIGAN STATE UNIV. LIBRARIES
[N"WWW(“WIWll(NW)llW(IWINIIWIHWI
31293104372234