Geophys. J. R . astr. SOC.

(1981) 64,425-446

Synthesis of seismic surface waves

N. J . Kerv* Department of Applied Mathematics and Theoretical Physics,

University of Cambridge, Silver Street, Cambridge CB3 9E W

Received 1980 April 25; in original form 1979 September 11

Summary. The reflection method is used to produce complete sets of

Rayleigh wave dispersion curves in a given phase velocity-frequency window
for horizontally stratified media, including modes of very high numerical
order, and theoretical surface wave seismograms are then synthesized.
The difficulties encountered when attempting to complete large numbers
of dispersion curves are discussed. Particular problems arise from models with
a low-velocity zone, when curves in a certain portion of the dispersion dia-
gram split into two families, the crustal and the channel modes. The reflection
method provides a convenient framework in which to examine this pheno-
menon heuristically and so devise a method to overcome the difficulties.
Seismograms are produced by mode summation and it is found that body-
wave behaviour, as well as surface-wave features can be synthesized. The
effect of truncating the mode series at a number comparable with the number
of modes used in previous studies is examined. It is found that although the
S-coda at longer ranges is not adversely affected, the arrivals attributable to
body-waves are severely distorted. T h s must call into question the validity
of synthetic seismograms generated by summation of only a few modes.

1 Introduction

The study of seismic surface waves was initiated when Rayleigh (1887) obtained a solution
t o the elastodynamic equations of motion representing energy trapped near the free surface
of a homogeneous elastic halfspace. Following this, other studies of trapped waves in simple
structures were carried out, notable examples being SH waves trapped in a uniform layer
overlying a homogeneous halfspace (Love 191 1) and P-SV waves trapped at the boundary
of two welded halfspaces (Stoneley 1924). The first discussion of the propagation of surface
waves in a more complicated multi-layered medium was by Sezawa (1927) who presented a

* Present address: AUWE, Southwell, Portland, Dorset DT5 2JS.

426 N. J. Keny
method applicable to structures with an arbitrary number of layers. The disadvantage of his
method is that for an n-layered model the condition for existence of a trapped mode (the
dispersion equation) is expressed as the vanishing of a 4 n - 2 rowed determinant for
Rayleigh waves or a 2n - 1 rowed determinant for Love waves.
No theoretical advances were made until Thomson (1950) presented a method for
obtaining the waves field in terms of a product of n 4 x 4 or 2 x 2 matrices for P-SV or SH
motion respectively. This method was corrected and cast into a form suitable for surface-wave
dispersion computations by Haskell (1953) and was first used on a computer by Dorman,
Ewing & Oliver (1960). Knopoff (1964) presented an alternative and more efficient method
which was first used for numerical work by Randall (1967). Gilbert & Backus (1966) pre-
sented the extension of the theory to arbitrary horizontal stratification by introducing the
propagator matrix, a generalization of the Haskell layer matrix.
The extension of Haskell's method to the problem of excitation of surface waves by
sources was done by Harkrider (1 964) and the computation of theoretical surface-wave
seismograms was first achieved by Knopoff, Schwab & Kausel(l973) using the fundamental
and first five higher modes. Their work has since been extended to slightly larger numbers of
modes (e.g. by Mantovani et al. 1977, using 11 modes).
The Thomson-Haskell or propagator method can be used to compute complete synthetic
seismograms, incorporating both surface and body waves. An alternative method, the reflec-
tion method, hitherto used for body-wave computations only, has been advanced by several
authors (e.g. Fuchs 1968). In this method the response is given in terms of the reflection and
transmission properties of the medium. The application of this method to coupled P-SV
waves was first presented by Kennett (1974) who used a formal framework applicable to all
other cases. The theoretical advantage of the reflection method is that results obtained in this
formalism are often capable of interpretation in a physically appealing manner, for example,
it provides a rigorous derivation of the formerly heuristically justified ray series expansion.
The numerical advantages of the method are that the quantities involved in the computations
are 2 x 2 matrices for P-SV motion and scalars for SH motion, and that the response is
given in a form which excludes the unphysical growing solutions present in the propagator
method. Kennett (1974) also derived the dispersion equation in the reflection method
formalism, when it can be interpreted as a constructive interference condition (cJ: Tolstoy &
Usdin 1953), and the formalism has been extended to obtain the response to buried sources
by Kennett & Kerry (1979).
This paper will describe the use of the reflection method to compute complete dispersion
curves in a given frequency window, including the theoretical modifications necessary t o
deal with the complications arising from the structure of the dispersion diagram when the
model contains a low-velocity channel. Although the method is applicable to arbitrary
horizontal stratification, the models used will be multi-layered. Because of the large numbers
of modes involved (about 170) it is necessary to sacrifice some of the accuracy of description
of the individual curves, and in particular those corresponding to energy trapped mainly in
the low velocity channel. However, it will be shown that this can be done without degrading
the final seismograms. The dispersion diagram will then be used t9 generate seismograms
for buried sources. Because we are using a considerably larger number of modes than has
been used in previous work, we would expect the seismograms to display high-frequency
phenomena in the early arrivals, w h c h is not the case when few modes are used.
T h ~ spaper relies heavily on the current state of the theory of the reflection method as
presented in Kennett & Kerry (1979). In order to avoid unnecessary and unproductive
repetition of material frequent references will be made to that paper. Sections and equations
in Kennett & Kerry (1979)will be denoted by the prefix KK. As in KK, the theory presented
 Synthesis of seismic surface waves 427
here is kept as general as possible so that it is valid for either P-SV or SH motion. However,
the numerical results are all for P-SV motion.
To compute seismograms it is necessary to transform from (r, Q1, z , t ) space to ( k , n , z , 0)
space by means of the transformations given in (KK2.1). The partial differential equations
of motion then become ordinary differential equations, which depend on the structure, and
on the two parameters k and w . A more convenient parameter pair is c and w or p and w,
where
1
c = - = w lks
P
c is the phase velocity and p the phase slowness.

2 The Rayleigh function

2.1 SURFACE WAVES

Surface waves are solutions of the elastodynamic equations of motion for a semi-infinite
medium bounded above by a free surface in the absence of any sources. The condition for
such a solution to exist is expressed as the vanishing of the Rayleigh function, and using
the reflection method formalism was derived in (KK 5.24) as:

A = det [(ND t NUR~~)(G~)-']

= 0. (2.1)
The matrices N U , ND are the stress transformation matrices for upgoing and downgoing
waves respectively (cf: KK 2.17). The matrices R and T are the reflection and transmission
coefficient matrices (cf: KK4.6). Subscripts U and D will be used to indicate that the coef-
ficients are generated respectively by upward or downward travelling waves incident on some
elastic medium which may be thought of as removed from the structure and sandwiched
between two homogeneous halfspaces. An (unordered) pair of superscripts will be used to
denote the medium of interest: superscript L will denote that the medium terminates at
zL, beneath which the structure is conventionally assumed homogeneous. Superscript F will
denote that the medium terminates at the free surface and 0 that it terminates just below
the free surface. Thus R g L and REL are different quantities, the second taking into account
free surface reflections and multiples. Quantities such as RLL have no physical meaning and
do not appear, but RCo, the free surface reflection matrix arises naturally and is denoted i.
The Rayleigh function, A, may be alternatively expressed by the introduction of a notional
receiver at depth z = Z R (cf: KK 5.27) as:

A = A' det [ ( I - RERREL)(TEL)-' 1,

where A' = det [(ND + N u R E ~ ) ( T ~ ~ ) -is' ]the Rayleigh function for the structure between
the free surface and z = z R . Because the existence and hence non-singularity of TEL and the
condition A ' # O imply the existence of REL and RER (cf Section 2.2) we may derive a
modified Rayleigh function:
= det ( I - RER REL). (2.3)
Note that (KK 5.25) may be derived from (2.3) by setting Z R = Z O .
It would appear at first glance that (2.3) is a better choice than (2.1) or (2.2) for the Ray-
leigh function, as its evaluation certainly involves less computational effort. However, it
suffers from three serious drawbacks. First, it is not possible to choose a normalization for
428 N. J. Keny
the transformation matrix, D (cJ: Kennet, Kerry & Woodhouse 1978) so that r\ is real,
wheras with non-attenuative media the energy normalized form of D given in -2.16)
ensures that A as defined by (2.1) is a complex function of constant argument and so may
be made real byJrivia1 scaling. Secondly, for a particular value of zR there may be values of
w , c such that A does not exist (cf Section 2.2). We cannot get round thls by a different
choice of zR as is not independent of zR whereas A may be evaluated from (2.2) and is
therefore clearly independent of z R . This indicates the need for the alternative expressions.
Thirdly, and most seriously, for any choice of Z R , near some roots 2\ is a slowly varying
function of w , c and near others it is a rapidly varying function. So much so that any numerical
estimate of it will fail to exhibit some of the roots. Thus we abandon (2.3) as a suitable form
for the Rayleigh function.

2.2 REFLECTION COEFFICIENTS AND TRAPPED MODES

We must now turn attention to those circumstances under which the reflection and trans-
mission coefficients for a given medium are not well defined. Consider the medium z1< z < z2
as shown in Fig. 1. From (KK4.4) the wave vectors at z l , and z2 are related by

V(zi) = QGi, ~ 2V(zd,

) (2.4)
and the reflection and transmission coefficients are given by (KK4.7) and are well defined
unless QZ2 is singular. In this case we may choose to be the null eigenvector of QZ2and
construct a wavefield from

V(z) =: Q(z, z2)(0, vg))T. (2.5)

Then the wavefield at z1 is given by

so that the wavefield represents a trapped mode for the whole space generated by sandwiching
the given structure between two halfspaces as indicated in the figure. The consequence of
this is that unless for a given structure the phase velocity and frequency are chosen so as to
admit a trapped mode, then the reflection and transmission coefficients exist. On energy
conservation grounds this will be the case unless waves of at least one type are evanescent
on both sides of the structure, and in the absence of a proof or counter-example we will
assume that trapped modes cannot exist unless waves of all types are evanescent on both
sides of the structure.
I
I

=1

Figure 1. A slab of stratified elastic material containing a low-velocity channel.

 Synthesis of seismic surface waves 429
In an analogous way it may be shown that the reflection coefficients from a structure
bounded above by a free surface also exist unless the given structure admits a trapped (surface
wave) mode.

2.3 COMPUTER ALGORITHM

Initially we attempt to produce dispersion curves using (2.1) to define the secular function.
To evaluate the reflection and transmission coefficients needed to compute the Rayleigh
function we use the iterative procedure described in KK Section 4.5. In this procedure we
build up the complete structure an interface at a time, starting at the deepest interface. We
use the iterative relations (KK 4.21) which we reproduce here for convenience

Note that in the first equation of (2.7) Rh3 which is the 'old' value of the overall reflection
coefficient occurs together with Rb2 etc, the coefficients for the layer or interface being
added, to produce Rk3, the 'new' value of the overall reflection coefficient. Thus RD , the
reflection coefficient for the current structure can be kept as a running variable during the
iterative process. Det(TD) may also be found directly since the second equation of (2.7)gives
det (Th3) = det (TA2) det (Th3)/det (I - R g R g ) , (2.8)
so det(TD) may be evaluated in a similar fashion without the need to evaluate the entire
matrix Tl, first.
The strategy for finding roots of the dispersion equation will be discussed fully in Section 3,
but will clearly involve some sort of iterative procedure for refining the root. This brings us
to the question of root accuracy, i.e. what is a suitable criterion for terminating the iteration?
Several different answers to this problem have been suggested. Our method is appealing
because it is extremely cheap computationally compared to, say, testing for equality of
kinetic and potential energy which would involve an additional integration through the
structure.
Consider the problem of determining the eigenfunction once a root has been found.
Let the wave vector at the surface be given by

v(zD>=(vU, vDIT, (2.9)

Then the two boundary conditions to be satisfied are:
Only downgoing waves at infinity
VU = R ~ ~ v D . (2.10)
No stress at the free surface
NU vu -k N D VD = 0. (2.1 1)

We use this form rather than VD =

-
Vu to avoid the singularity in R when considering the
fundamental mode at high frequency when its phase velocity asymptotes to the Rayleigh
wave speed appropriate to the material at the surface. The amplitude VD is related to the
430 N. J. Keny
energy density near the surface. If Vu is calculated from (2.10) then the value yielded by
the left hand side of (2.1 1) is the surface stress.
From (2.10) and (2.1 1) we obtain:

(2.12)

so that VD is the null eigenvector of ( N D + N u R g L ) . Thus it would seem reasonable to

terminate the iteration process when two estimates of the root, corresponding to opposite
signs of the Rayleigh function differ by less than some predetermined small number and the
smallest eigenvalue of (ND t N U RgL), evaluated at one of the two estimates is less than some
other predetermined small number. This will be the 'surface stress per unit near-surface
energy density' in the surface wave mode.
Two models that have been used for computational work are shown in Fig. 2 . Note that
both these models have low-velocity channels. It transpires that when the above procedure is
used to compute dispersion curves, for those values of phase velocity such that there are
propagating waves in the surface zone and in the low-velocity zone, but not in the lid, there
are some roots of the dispersion equation for which apparently neither of the eigenvalues
of the matrix in (2.12) become small as the root is further refined. The reason for this is
that for these modes the energy density near the surface is small, most energy residing in the
low-velocit y channel.

2.4 CRUSTAL WAVES AND C H A N N E L WAVES

It has been noted (e.g. see Panza, Schwab & Knopoff 1972; Frantsuzova, Levshin &
Shkadinskaya 1972) that for a structure with a low-velocity zone the dispersion curves in
a particular portion of w-c space can be divided into two 'families', the crustal modes and

1000-
Figure 2. Models used for numerical work in this paper, derived by applying the earth-flattening
approximation to spherical models and producing a layered approximation. (a) Model T7, (b) model
1066B.
 Synthesis of seismic surface waves 43 1

I I / /

P(Z)

Figure 3. Schematic representation of a stratified elastic halfspace containing a low-velocity zone, showing
the four zones into which it is divided when there are three S-wave turning points.

the channel modes. This phenomenon is related to the problem described at the end of the
previous section. Consider the structure shown in Fig. 3. If
blvz <c <hid 9 (2.13)
the structure can be naturally divided into four zones by the three turning points forS-waves.
Zone 1 is the surface zone with S-waves propagating, zone 2 the lid with S-waves evanescent,
zone 3 the low-velocity channel with S-waves propagating, and zone 4 the substructure
with S-waves evanescent. With the ‘receiver’ at Z R in the channel as indicated, we may write
the dispersion equation (2.2) as
FR RL
det (ND t N ~ Rdet ~(I -~RU) R D ) = 0, (2.14)
since the other factors are non-zero by the assumptions of Section 2.2. Since the lid is a
region of evanescence, the transmission coefficients through any portion of it will become
exponentially small as w + 00. Choosing zA at the ‘peak’ of the lid as indicated it is reasonable
to assume by virtue of the iterative relations (2.7) that

(2.15)

Thus we see that the lid acts as a screen, so that the material on one side does not ‘see’ the
material on the other, and we may approximate the dispersion equation as
AR RL
det (ND t N ~ Rdet ~(I -~RU) R D ) = 0. (2.16)
which is the product of the dispersion equation for surface waves on the structure above
z=zA

det ( ND t N ~ R= 0~, ~ ) (2.17)

432 N. J. Keny

f r e e surface free surface

Figure 4. Schematic representation of the splitting of the structure into the surface zone and low-velocity
channel at high frequency.

and the dispersion equation for channel waves on the structure below z = z,4
det (I - R e R R g L )= 0, (2.18)
which may be expressed pictorially as in Fig. 4. This provides an explanation for the existence
of these two families, the crustal modes being those for which (2.17) is approximately
satisfied and the channel modes being those for which (2.18) is approximately satisfied.
It also suggests an explanation of the apparent paradox of our inability to solve (2.12).
Since R g R 2 R g A it is tempting to assume that also R g L z REA. Although this is true
for most values of w , c, if we chose values corresponding to a root of (2.18), i.e. corresponding
to a trapped mode in the low-velocity channel then we cannot use the iterative relations
(2.7) to yield the above approximation owing to the singularity then implied in T k L .However,
it may be seen from (2.7) that the region over w h c h the approximation R g L R g A breaks
down is given by
R L) 5 6,
A RR D
det ( I - RU (2.19)

where 6 is of the order of the decay factor for a double pass through the lid by S-waves.
Since the elements of the matrix ( T - R ( j R R E L ) are of order unity, accurate evaluation of
this determinant is impossible on a computer when 6 is comparable with the machine
accuracy. Thus for sufficiently high frequencies accurate numerical evaluation of R E L is
impossible near channel wave frequencies and it is this which causes the problem.
It must now be considered whether this situation is due to a drawback of the method,
or whether it is unavoidable when using a finite accuracy machine. The author believes
the latter to be the case, as this method is the only one presented to date in which the
exponentially growing solutions are completely excluded. The fine details of the crustal-
channel interaction are due to extremely weak coupling across the lid, and these exponentially
small effects are lost when they are swamped by the order unity terms. When, however, a
method is used which does not exclude the growing solutions, these will swamp the weak
coupling effects at much lower frequencies. Thus all other methods will be at least as bad,
and usually worse. than this one.

2.5 THE MODIFIED COMPUTFR ALGORITHM

In view of the remarks in the previous sections we modify the termination criterion of the
iterative procedure described in Section 2.3. If condition (2.13) is satisfied so that the
possibility of channel waves arises, we pick zli to lie at the top of the low-velocity channel,
 Synthesis of seismic surface waves 43 3
immediately beneath the turning point zp2.We then perform a single pass through the struc-
ture to calculate the quantities REL, RD”L,d e t ( g L ) , d e t ( e R ) , det(TEL) and hence evaluate
the Rayleigh function from (2.1). We define

q E det ( I - R g R R s L ) , (2.20)
which by virtue of (2.7) may be evaluated as
77 = det (TgR) det (TgL)/det (TgL). (2.21)
If 7) is larger than some predetermined small number we proceed as in Section 2.3. If 7) is
smaller we perform a second pass downwards through the structure from zo to zR to
compute RER. The criterion for termination of the iteration in this case is that the smallest
eigenvalue of ( I - R ER EL)
R must be less than the small number chosen for termination in
Section 2.3.
By this method we can classify the roots a priori as channel or crustal according as to
whether 77 is small or not.

3 Dispersion curve computation

3.1 STRUCTURE REDUCTION

In view of the remarks made in previous sections that an evanescent region acts as a screen
so that the structure beyond it normally has little effect on propagation properties, it would
seem reasonable to suppose that the fine details of the ‘deep’ structure of a stratified half-
space do not significantly affect the dispersion curves. If this is the case it would clearly
be advantageous to replace the structure below some suitable level z = ztmnc by a simpler
structure if we can thus effcct a significant saving in computer time with negligible loss of
accuracy. We replace the structure below zmnCby the simplest possible: a uniform halfspace
whose material properties are those of the original structure at ztrunc.
We determine ztrunc as follows: if zpL is the deepest turning point in the structure, so
that below zpL all waves are evanescent, we define for z > zpL

(3.1)
VL
Then ztruncis defined by the condition

u(Ztrunc) urnax. (3 3
This condition was first used by Wiggins (1 968) (equation 12).
The effect of choice of urn= on root accuracy may be determined as follows. For a
sample of roots (mi, ci)of the dispersion equation we consider the structures with y 1 = 1 , 2 , . . .
layers retained beneath the deepest turning point. We calculate u from (3.1) and an estimate
of root accuracy from

(3.3)

where A, is the Rayleigh function for the reduced structure and A, the Rayleigh function
for the entire structure. Fig. 5 gives a plot of loglo(A m)against u.
In subsequent numerical work we use urnax = 10, which will be adequate when seeking
roots to an accuracy of
 - 2.0

- 4.0

- -6.0

I
-14.0 -- i x

- 16.0

- 18.0 i1 x
-20.0 ---" " ' ' I 1 ' 1 ' ' 1 1 ' 1 1 1 ' ' 1 ' 1 ' ' 1 1 ' ' l ' ' " ' ' ' t ' l

00 5.0 10.0 15.0 20.0 25 0 300

urnox

Figure 5. Plot of the error introduced into the roots of the dispersion equation by truncating the structure
against structure truncation parameter, urnax.

3.2 DISPERSIONCURVE GENERATINGPROGRAM

The basic requirement of any dispersion curve generating program is simple to state: given
a real function A of two variables w and c we seek zeros of that function in a predetermined
portion of a - c space. An obvious method, which we shall in fact use, is to fix one of the
variables and vary the other, searching until we find two points at which the function takes
opposite signs. We then use some iterative process to obtain an accurate estimate of the root.
Guided by Schwab & Knopoff (1972) we use quadratic interpolation to refine the root.
The first problem is the choice of which parameter to fix while searching for roots. We
choose to fix c and vary a.This has three advantages: first, since the models we have used
for numerical work consist of uniform layers we can utilize the fact that for a given structure
the interface coefficients are independent of frequency, (cf KK48) noting that the matrix
D is independent of frequency), so that considerable computational effort may be saved;
secondly, the zeros of A are approximately evenly spaced in w while the spacing in c is
highly irregular; and thirdly, the procedure can be adapted directly to spherical geometry
where the analogous parameter to c only takes a set of discrete values.
The second problem is to ensure that all roots have been found for a particular value of
c. This can in general be done by keeping track of each particular mode curve as c varies,
although the details of a computer program to do this are necessarily complicated. It must
also be borne in mind that if two mode curves osculate closely it may not be possible to
'split' them without inordinate expenditure of computer time.
The third problem is to choose a suitable new value for the parameter c having found all
the roots at one particular value, c o , say. A reasonable requirement in accurate determination
of the dispersion curves is that Am, the frequency jump on a given mode curve between
adjacent values of c , is small.
 Synthesis of seismic surface waves 43 5
Asymptotic theory indicates that, for fixed c , the slope of the dispersion curves, dwldc,
is proportional to w, so that this condition might best be expressed as
la4 < 4 , (3.4)
where ( is some small number. Thus we set c = c,, + 6c where

1 dw
w dc

which should ensure that (3.4) is satisfied. We may evaluate dw/dc by means of the implicit
function formula

Unfortunately, this procedure is not entirely satisfactory, as we find that if c -+ plVz + then
dw/dc I -+ some finite limit, but do/dc I channel -+ - m. Thus we find that for values of c

near flkz accurate determination of the channel modes involves exceedingly fine spacing of
the c values compared to that required for accurate determination of the crustal modes.
This problem may be overcome by reconsideration of the eventual aim of the program,
which must be to produce surface wave seismograms for either a surface or a buried source,
but with always a surface receiver. It is then reasonable to suppose that the channel modes
which cause the above problem will not make a noticeable contribution to the seismograms
since they represent energy trapped at depth. Thus we settle for considerably less accuracy
in the channel wave branches of the dispersion diagram by choosing

(3.7)

Recall that we can detect whether a mode is crustal or channel a priori by the criterion given
in Section 2.5.
Having chosen a new value of c and having determined the corresponding roots of the
dispersion equation we now have to decide whether the picture we obtain by joining points
on the dispersion curves by straight line segments will give a sufficiently accurate picture.
The problems this can involve are illustrated by consideration of Fig. 6. The roots corres-
ponding to crustal modes are indicated by crosses, and the roots corresponding to channel
modes by triangles. The estimated gradients from (3.6) of the curves are also indicated by
arrows. The 'actual' dispersion curves are indicated in Fig. 6(a), the result of joining the
roots by mode number in Fig. 6(b), and a third possibility in Fig. 6(c). This is the result of
joining the points so that there is the minimum discrepancy between the slopes as given by
(3.6) and the slope of the straight line segment. It is clear that this is a much more satis-
factory picture than Fig. 6(b).
Having formed an assignment such as that indicated we have to decide whether it gives a
satisfactory description of the dispersion curves. There is a clear trade-off between accuracy
on the one hand and computer time taken on the other. It is quite possible than no significant
improvement in description would be attained until many tens of further c values between
c,, and c, + 6c were sampled, and we will discuss whether t h s would affect the final seismo-
grams in the next section.
We therefore perform two simple checks: first that (3.4) is in fact satisfied for all crustal
modes, where modes are considered crustal if either end point is a crustal solution of the
436 N. J. K e n y

I I i \ \

Figure 6. Roots of the dispersion equation at adjacent values of phase velocity. Crustal roots are indicated
by crosses, channel roots by triangles. (a) The true dispersion curves. @) The approximate curves corre-
sponding to joining the roots by mode number. (c) The approximate curves given by joining the roots so
as to minimize the kink in each curve (see text for details). Note that the c-scale is orders of magnitude
larger than the w-scale.

dispersion equation; and secondly that the ‘kink’ in each of the curves, both crustal and
channel, defined by
6 e = 10, - e,l t ie3 - 021, (3.8)
is less than some predetermined small number, where 01, 0 2 , O3 are the three estimates of
the slope of the mode curve given by (3.6) at the bottom end point, the slope of the straight
line segment, and (3.6) at the top end point respectively, as indicated in Fig. 6(b).
If these tests are both successful we pick a new value of c from (3.7) and then proceed as
before. If, however, either of these tests fails then we halve 6c and again proceed as before,
halving 6c until either a satisfactory picture is obtained or some minimum c-step is reached.
Sets of dispersion curves for the two models shown in Fig. 2 generated by this process
can be seen in Fig. 7. They will be discussed in Section 4.4. A fuller dispersion diagram for
model T7 is shown in Fig. 9. The parameters used to generate this were { = 0.05 (cJ: equation
3.4) and 6c (cf. equation 3.5) constrained to lie in the range 0.001 < 6c < 0.03 km s-’. Of the
order of 60000 roots of the dispersion equation are plotted. Approximately 1 hr of com-
puter time was used on an IBM 370/165 machme using double precision.
 Synthesis o f seismic surface waves 43 7

5.0

111
\

E
1
4.5
u

TJ
-
0,
u
Q
10 20 3.0 4.0 50 60
m
u
2 5.0
I
a

I I
Frequency w rad/s
Figure 7. Sections of the Rayleigh wave dispersion diagrams for the two models shown in F i g . 2.

Comparison of the time taken per layer per iteration in computing the Rayleigh function
(cfi Schwab & b o p o f f 1972)indicate that the method is comparable with Knopoffs method
for the non-attenuative media (as in this case), but will suffer no increase due to the introduc-
tion of complex wave velocities (attenuative media), thus introducing a saving of a factor
of 9 over existing methods.

4 Synthetic surface-wave seismogram generation

4.1 T H E R E S P O N S E O F A HALFSPACE TO A POINT S O U R C E

The response of a halfspace to a point source was discussed by Kennett & Kerry (1979)
(Section 5) where two expressions were derived for the disturbance at the receiver, depend-
ing on whether the receiver is above the source (KK 5.17) or below the source (KK 5.20).
In this paper we are concerned with the surface wave portion of the response to a buried
source with a receiver at the surface, so we would expect to use (KK 5.17):

~ ( Z R )= ( M E -t MBk)(Z - RBSR")-'TBS(I- R ~ L R ~ s ) - l ( R- ~xu),

L~~ (4.1)
and we note (cJ: Section 2) that the surface wave singularity is incorporated via the term
(Z - RgL R E S ) - ' . For the reasons discussed in Section 2.3 a numerical estimate of this
term may not accurately exhibit the singularity when a buried source is considered. Also
the residue at the surface wave pole will depend on source depth and so will require recompu-
tation for each source depth considered. Thus (4.1) suffers from both inefficiency and
inaccuracy and is therefore unsuitable for computing the response. A more convenient form
may be derived as follows. Let

v ( z R ) = ( V R U ) VRDlT, (4.2)
then

V(ZR)= Q(zR, zs)V(zs-> (4.3)

438 N. J. Keny
and

From (KK 5.10)

V(Z,t) = V(z,-) + 2,
so that

U = (uu, UD)* = V(ZL)- Q ( Z L , Z S ) ~ . (4.7)

Thus, since from the boundary conditions V(ZL) represents only downgoing waves, we may
obtain, using (KK4.9)

Uu = ( @ - ' { R g L Z ~ - Xu), (4.8)

and (4.6) may be inverted in a similar manner to obtain

(4.10)

(cf KK4.17), and since the receiver is above the source, the condition of vanishing stress at
the surface is

vRU + ~ f vKD
2 = 0, (4.1 1)

and the displacement at the receiver is

W(ZR) = MEV R U + MEVRD. (4.12)

Now we may eliminate Uu from (4.8) and (4.9) and solve the resulting equation with (4.1 1)
to obtain

(4.13)

where we have also made use of the iterative relations (2.7). This may be shown to be
equivalent t o (KK 5.1 7) by use of the iterative relations and the fact that (KK4.18):

R -1 FR -RFR (4.14)
- (F22) ( 21) - u .
Equation (4.13) is better than (4.1) as the surface wave singularity is now incorporated via
the term (FF2 t FFl R:L)-', which for a surface receiver is (ND + N U R g L ) - ' . This will
exhibit the singularity more satisfactorily, and the residue at the surface wave poles may be
computed independent of source depth.
 Synthesis of seismic surface waves 439
4.2 THE S U R F A C E W A V E RESPONSE

The source used for computational work is a vertical dip-slip dislocation. This source is
equivalent to a double-couple without moment so that the force is given by (KK Al) with

Fi = 0,
Mxz = M z x = Mo,
Mij = 0 otherwise. (4.15)

Then (cf KK A2 and KK A3) the only angular orders involved are m = f 1. The seismograms
are obtained by inverting the transforms (KK 2.1)

wz(r,@, t ) =
277
1m

-! dw exp (- iwt) 1
-m
W

m=-m
exp (irn@)I0 m

dkWm(kr)U(k,m, a). (4.16)

We will work with the vertical component only, as this will not involve the Love waves. As
we are only interested in the surface wave portion of the seismogram we may distort the
contour in the k-plane to pick up the contributions from the poles and disregard the rernain-
ing contour integral which corresponds t o the body wave portion of the seismogram and
contains arrivals outside the phase velocity window we are concerned with.
Using the symmetry properties of U (Hudson 1969) we obtain

w,(T,@, t ) = i cos @ 1 dw exp (- iwt) ki(w)HJ1)(ki(w)T) U, (4.17)

modes -m

where we have used the fact (KK A2) that U(k, -1, w ) = - U(k, 1, w ) and where
ReskilwU= l h ((K - ki(w))U), (4.18)
K+ki(w)
or
af 1
R e s k i I w l r =I l a k l ,wheref--
W U'
Note that

(4.19)

where cg = l/[dk(w)/do], the group velocity.

The seismogram may now be evaluated as

w,(T, @, t ) = i cos @ 1
modes sy
-c.J
d o exp (-iwt) ki(w)ff$')(ki(U)r) (-1
c
- -1
1
cg.
Res, I ciu (4.20)

cg may be determined from dw/dc by

(4.2 1)

We prefer (4.20) to (4.17) to evaluate the seismograms as the residues are evaluated with
c held constants so that we can take advantage of the fact that the interface coefficients are
independent of frequency to evaluate the integrand more efficiently.
440 N. J. Keny
To evaluate the integrals (4.20) for each of the modes we will use the fast Fourier trans-
form algorithm (Cooley & Tukey 1965). It is necessary to use some form of interpolation to
provide values of the integrand at small frequency steps when computing seismograms at
large epicentral distances. Suppose we wish to compute the seismogram at a given range R with
a time window -T < t Q T and with sampling every 6T. Then by the usual relationships for
the discrete Fourier transform the frequency spectrum will be -52 < w < R , R = n/6T with
frequency step 652 = n/T. For accurate evaluation of the integral the argument of the Hankel
function must not vary too much between frequency points, so that if Omax is the maxi-
mum permissible change in argument we must have
dk R
6kRZ6QR -=66-<<m,, (4.22)
dw cg
whence 652 < (cg/R)Bmax for all points on the dispersion curve, so that

Cgminernax
6!2 = (4.23)
R '
where cgmin is the minimum value of the group velocity on the given dispersion curve. If
'gmin
is the minimum value of the group velocity on all dispersion curves the spectra can
b e combined before inversion. Independent of the initial spacing of points on the dispersion
curves it will be necessary to interpolate values to ensure that (4.23) is satisfied at large epi-
central distances.

4.3 MODE EXCITATION

Equation (4.20) essentially involves integration along the dispersion curves. We must therefore
consider the effect on these integrals of the approximations introduced into the dispersion
curves by the algorithm described in Section 3.2. We will first examine the relative amplitudes
of the responses clue to crustal and channel modes. Equation (4.1 3 ) indicates that the ampli-
tude will be given by
OL -1
ResuIc(ND NURD ) . (4.24)
By use of the iterative relations (2.7) we obtain

OR -1
(NDtNURD ) (4.25)
The only singular term on the right hand side is (Z - RER R :")-', the dispersion function
appropriate to channel modes, which is multiplied by the terms fiR gR.
and Thus the
surface displacement due to channel modes will be less than that due to crustal modes by
the decay factor due to a double pass through the evanescent lid zone.
Consider now Fig. 8 which indicates the integrand in (4.20) for two adjacent modes near
an osculation point. This is analogous to Fig. 6, and indicates the integrand as approximated
by joining up the points on the dispersion curves in different ways. Fig. 8(a) indicates the
actual integrand, Fig. 8(b) that given by joining up the points by mode number as in Fig. 6(b),
and Fig. 8(c) indicates the integrand as given by the alternative method of joining the points,
as in Fig. 6(c). It is clear that this third alternative will give a considerably more accurate
estimate of the integral. It is also worth pointing out that the loss of accuracy in the descrip-
tion of the channel-wave dispersion curves is unimportant in this instance - the decision we
took in Section 3.2 to settle for considerable loss of detail in the channel waves has only
discarded unimportant information.
 Synthesis of seismic surface waves 44 1

W
(c)

Figure 8. Schematic representation of the magnitude of the integrand in equation (4.20) at points on
adjacent mode curves near an osculation point. (a) The true integrand. (b) and (c) The approximate
integrand implied by the choices in Fig. 6(b) and (c) respectively.

3.0
1
tui L
'
~ ~ L LAU iI -
L
~~

LLI LLLULLLUL
~

,~d
00 1.o 2.0 3.0 4.0 5.0 6.0
Frequency w rad/s

Figure 9. The Rayleigh wave dispersion diagram for model T7.

442 N. J. Keny
4.4 NUMERICAL RESULTS

A comparison of corresponding portions of the dispersion diagrams in the frequency range

0.063 < w < 6.3, i.e. a period range of 1-100 s for the two models illustrated in Fig. 2, is given

expected in the light of Section 2.4.

180
-
in Fig. 7. Both sets of dispersion diagrams clearly exhibit the general structure that is to be

5
t
m

pj120
L
0
m
'
I
150

90

60

30
T I
i

0
1000 1200 1LO0 1600
R km
Source Depth 10 km

1 u .
t
i
+
200 100 600 800 1000
1200 lL00 1600
R km
Source Depth 60 km
(b)
Figure 10. Synthetic surface-wave seismograms for model T7 generated by a vertical dip-slip source with
smoothed 6-function time dependence. R is epicentral distance and the seismograms are plotted in
reduced time as shown. (a) Source depth 10 km, (b) source depth 60 km, (c) source depth 200 km,
(d) source depth 60 km, showing comparison with seismograms produced by including only 15 higher
modes.
 Synthesis of seismic surface waves 443
For values of phase velocity in the range 4.45 < c < 4.77 for model 1066B and 4.52 < c <
4.67 for model T7 there are three turning points for S-waves in the medium and as expected
the dispersion curves in these phase velocity ranges form two families. The problem that this
causes can be clearly seen in the flatness of the channel mode curves compared to the crustal
mode curves. Other differences in the two dispersion diagrams can be explained in terms of
the differences between the two models. First, there are many more modes at each value of
phase velocity for model 1066B than for model T7 - at c = 5.1 there are about 70 higher
modes in the period range 1-100s for model 1066B, but only about 45 higher modes for

120 -

90 -
Ln

" -
Ln
h
0
Lo 6 0 -
I
5 -
I
t
30 -

O L LU
200 LOO 600 800 1000 1200 1400 1600
R km
Source Depth 200 km
(C)

0 200 LOO 600 800 1001' 1200 1400 1600

R km
Source Depth 60 km
Figure 10 (d)
444 N. J. Keny
model T7. This is due to the fact that the turning point for a given phase velocity occurs at
greater depth in model 1066B than in model T7 so that asymptotic theory predicts more
modes for model 1066B.
Another feature of interest is that the osculation points between crustal and channel
modes are much ‘tighter’ for model 1066B than for model T7. This is due to the fact that
t h e lid in model 1066B is much thicker, extending from about 20 to 285 km in depth as
compared to the lid in model T7 which extends from about 55 to 155km and also to
t h e fact that the lid is much deeper. Thus the cross-lid coupling is much weaker for model
1066B than for model T7, making the osculation points tighter.
Since reflections from velocity gradients occur as low-order effects in asymptotic theory,
and both the models have strong velocity gradients and interfaces, we would expect to see
features of the dispersion diagrams not explained by first-order theory. One such feature is
that the crustal modes can be traced apparently continuing, as aligned ‘kinks’ in the dispersion
curves, into the region of the dispersion diagram with phase velocities greater than the shear
wave velocities in the lid. This is not predicted by first-order theory and is caused by the
strong velocity jump at the Moho. This feature of the dispersion curves will necessarily
affect the synthetic seismograms and will give rise to strong near surface reverberations.
Fig. 9 displays the complete dispersion diagram for model T7 in the period range 1-100 s,
and the phase velocity range up to 6.8 km sC1, corresponding to surface waves penetrating
t o the 700 km transition zone. A feature of this diagram, also not predicted by asymptotic
theory, is the perceptible upward kink in all the dispersion curves in the phase velocity
range 5.4 < c < 5.7. This velocity range corresponds to the S-wave velocity change across the
400 km transition zone which gives rise to this feature.
Synthetic surface wave seismograms for model T7 computed from equation (4.20) are
shown in Fig. 10. As discussed in Section 4.2, the source is a dip-slip dislocation with a vertical
fault-plane, and the time dependence is a smoothed &-function.The model receiver has a
broadband frequency response between the -3dB points at 1.5 and 66 s.
As we are summing all modes penetrating to 700 km we would expect, in view of the work
of Nolet & Kennett (1978), to synthesize body-wave type arrivals for S-waves turning above
7 0 0 km.
The seismograms are plotted without distance-related amplification so that for each source
depth the true relative amplitudes of the seismograms are preserved.
For source depth 1Okm. (Fig. lO(a)) the tail of the S-wave train is well represented as a
higher mode surface-wave train with no clear phases, and the fundamental mode is very clear.
However, there is a precursor developing at longer ranges which is a numerical arrival caused
b y phase velocity windowing.
For source depth 6 0 km (Fig. 1O(b)) the S coda can be seen developing at longer ranges,
with a very much reduced amplitude in the fundamental mode, the amplitude being down
b y a factor of 10 compared with the previous case. Also clearly visible are the multiple surface
reverberations of the pulse-type arrivals. Thus as we anticipated we have successfully synthe-
sized body waves by mode summation. The seismogram at 50 km range shows a long period
onset which is also caused by phase velocity windowing.
For source depth 200 km (Fig. 1O(c)) no arrivals are apparent until the horizontal range is
comparable with the depth of the source. This is again due to windowing effects as the direct
ray from the source to the receiver would have ray parameter corresponding to a turning
point below 700km and hence the surface wave modes required to synthesize this ray would
have phase velocity greater than the maximum included in Fig. 9. At longer ranges there is
a good representation of pure body-wave arrivals and the wave trains show very little surface
wave character.
 Synthesis of seismic surface waves 44 5
Fig. 10(d) gives a comparison of the seismograms generated by summation of all the
modes (as in Fig. 1O(b)) with those generated by including only 15 modes in the summation.
The differences are very clear. Only the surface wave type arrivals at long times are well
represented in the truncated sum, the open triangle indicating the point after which the
waveforms agree. The earlier arrivals, and the entire seismograms at shorter ranges where
the surface-wave train has not developed, exhibit a lack of high-frequency energy. This is to
be expected as the effect of truncating the mode series is to apply a low-pass fitter with a
cut-off frequency reducing at higher phase velocities. This will have more effect on the
earlier arrivals, and because of this the pulse shapes of the body-wave type arrivals for the
truncated series have become corrupted and are non-causal with ill-defined onsets.

Acknowledgments
The author would like to thank Dr B. L. N. Kennett for many useful discussions and continual
encouragement during the course of this research. The work was carried out while the author
was supported on a NERC studentship.

References
Cooley, J. W. & Tukey, S. W., 1965. An algorithm for the machine calculation of complete Fourier series,
Math. Comp., 19, 297-301.
Dorman, J., Ewing, M. & Oliver, J., 1960. Study of shear-velocity distribution in the upper mantle by
mantle Rayleigh waves, Bull. seism. SOC.Am., 50, 87-115.
Frantsuzova, V. I., Levshin, A. L. & Shkadinskaya, G. V., 1972. Higher modes of Rayleigh waves and
upper mantle structure, in Computational Seismology, ed. Keilis-Borok, V. I., Consultants Bureau,
New York.
Fuchs, K., 1968. The reflection of spherical waves from transition zones with arbitrary depth-dependent
elastic moduli and density, J. Phys. Earth, 16, 27-41 (special issue).
Gilbert, F. & Backus, G . E., 1966. Propagator matrices in elastic wave and vibration theory, Geophysics,
31,326-332.
Harkrider, D. G., 1964. Surface waves in multilayered elastic media: I-Rayleigh and Love waves from
buried sources in a multilayered elastic halfspace, Bull. seism. Soc. Am., 54,627-679.
Haskell, N. A., 1953. The dispersion of surface waves on multilayered media, Bull. seism. Soc. Am., 43,
17-34.
Hudson, J. A., 1969. A quantitative evaluation of seismic signals at teleseismic distances - 11. Body waves
and surface waves from an extended source, Geophys. J. R. astr. Soc., 18, 353-370.
Kennett, B. L. N., 1974. Reflections, rays and reverberations, Bull. seism. SOC.Am., 64, 1685-1696.
Kennett, B. L. N. & Kerry, N. J., 1979. Seismic waves in a stratified halfspace, Geophys. J. R. astr. Soc.,
57,557-583.
Kennett, B. L. N., Kerry, N. J. & Woodhouse, J. H., 1978. Symmetries in the reflection and transmission
of elastic waves, Geophys. J. R. astr. Soc., 52, 215-229.
Knopoff, L., 1964. A matrix method for elastic wave problems, Bull. seism. Soc. Am., 54,431-438.
Knopoff, L., Schwab, F. A. & Kausel, E., 1973. Interpretation of Lg, Geophys. J. R. astr. Soc., 33,
389-404.
Love, A. E. H., 1911. Some Problems of Geodynarnics, Cambridge University Press.
Mantovani, E., Schwab, F. A., Liao, H. & Knopoff, L., 1977. Generation of complete theoretical seismo-
grams for SH - 11, Geophys. J. R. astr. Soc., 48, 531-536.
Nolet, G. & Kennett, B. L. N., 1978. Normal-mode representations of multiple-ray reflections in a
spherical earth, Geophys. J. R. astr. Soc., 53, 219-226.
Panza, G. F., Schwab, F. A. & Knopoff, L., 1972. Channel and crustal Rayleigh waves, Geoplzys. J. R.
astr. Soc., 30, 273-280.
Randall, M. J., 1967. Fast programs for layered halfspace problems, Bull. seism. Soc. Am., 57, 1299-
1316.
446 N. J. Keriy
Rayleigh, Lord, 1887. On waves propagated along the plane surface of an elastic solid,Proc. Lond. math.
SOC., 17,4-11.
Schwab, F. A. & Knopoff, L., 1972. Fast surface wave and free mode computations, in Methods of
Cotnputational Physics, vol. 11, ed. Bolt, B. A., Academic Press, New York.
Sezawa, K., 1927. Dispersion of elastic waves propagated on the surface of stratified bodies and on curved
surfaces, Bull. Earthq. Res. Inst. Tokyo, 3, 1-18.
Stoneley, R., 1924. Elastic waves at the surface of separation of two solids,&oc. R . Soc. Lond. A , 106,
416-428.
Thomson, W. T., 1950. Transmission of elastic waves through a stratified soIid medium, J. appl. Phys.,
21, 89-93.
Tolstoy, I. & Usdin, E., 1953. Dispersive properties of stratified elastic and liquid media: a ray theory,
Geophysics, 18, 844-870.
Wiggins, R. A., 1968. Terrestrial variational tables for the periods and attenuation of the free oscillations,
Pliys. Earth platzet. Iriteriors, I , 201 -266.