Membrane CHT THESIS

Waves on a Spherical Membrane
by
Carl Harlan Tape
Thesis submitted to the University of Oxford

for the Degree of Master of Science by Research
Department of Earth Sciences

and
Worcester College
University of Oxford
July 15, 2003
Advisor: Prof. John H. Woodhouse

i
Abstract
The purpose of this study is to develop a numerical model for wave propagation on a spherical
membrane and to compare the results with predictions made from ray theory. Such membrane
waves are an analogue for seismic surface waves. The two-dimensional wave equation on
a sphere is solved using a finite-difference method on a spherical grid of several thousand
hexagonal faces. The spherical grid is first tested using a homogeneous phase velocity field
with an initial Gaussian source. The resulting numerical solutions, u(, , t), agree with
analytical solutions for the homogeneous case. Next a heterogeneous phase velocity field,
c(, ), is used; the solutions reveal the variation in amplitude and phase due to the lateral
heterogeneity.
Surface wave ray-tracing equations are then used to calculate the amplitude and phase
anomalies for a set of source-receiver pairs. Exact ray theory calculates these values along
the exact ray path, which may deviate considerably from the great-circle path between the
source and receiver, which is the basis for the linearized ray theory calculations. We find that
multipathing multiple ray paths between a source and receiver is quite common for the
current resolution of long-period surface wave phase velocity maps, and it gives an indication
of the divergence of linearized ray theory from exact ray theory.
Using the results from the numerical model, we are able to examine the regimes under
which ray theory predictions are invalid, and we discuss the results in terms of the ray theory
validity condition, , where is the wavelength of the waves in the numerical model
and is the minimum scalelength of heterogeneity in the phase velocity map. We show that
exact ray theory is better than linearized ray theory at predicting both phase and amplitude
anomalies. We find that ray theory predictions of phase are quite stable. Predictions of
amplitude, however, are valid only when the validity condition is well-satisfied.
ii
Waves on a Spherical Membrane
Extended abstract
The purpose of this study is to develop a numerical model for wave propagation on a
spherical membrane and to compare the results with predictions made from ray theory. Such
membrane waves are an analogue for seismic surface waves. Ray theory has been used in
seismology to explain observations of path differences, amplitude differences, and phase dif-
ferences of surface waves with respect to a spherical reference model (e.g., Woodhouse and
Wong, 1986; Laske et al., 1994). Furthermore, it is a standard tool used in creating 3D global
models of the earth. However, ray theory is an approximate theory and cannot always be
used to explain observations. Thus, numerical models are helpful, especially as earth models
become increasingly more detailed. This study presents a two-dimensional numerical model
aimed at investigating the regimes under which ray theory predictions are invalid.
Analytical solutions
Our study is concerned with analytical and numerical solutions to the wave equation on the
sphere:
1 2u
2 u = f ,
c2 t2
where u = u(, , t) is the wave field, c = c(, ) is the frequency-independent phase velocity,
f = f (, , t) is the forcing term, and 2 u is the surface Laplacian of u on a sphere with
radius a, where , are spherical coordinates.
For the analytical solutions, we specify a constant phase velocity function, as well as the
forcing term f (, , t), which describes the source. The source we use is simply a Gaussian in
time (source-time function) multiplied by a Gaussian in space (initial-shape function).
Numerical solutions
We employed a finite-difference method (FDM) on a special spherical grid in solving the wave
equation numerically, using a method developed in the field of climate modelling (Williamson,
1968). The grid was constructed by subdividing triangular patches on the sphere into four
iii
new triangles at each level of subdivision. This process, which was carried out on grids with
up to 122,882 points, achieves a relatively uniform distribution of gridpoints. The numerical
approximation for the surface Laplacian ( 2 ) in the wave equation is based on a local FDM
which averages the function gradients at a point with its six adjacent neighbors.
We are able to test the numerical solutions by comparing them with the analytical solutions
for the homogeneous case, that is, c(, ) = const. We find good agreement between the
numerical and analytical solutions, and this gives us confidence in the numerical solutions for
the heterogeneous case. For the heterogeneous case, we use surface wave phase velocity maps
for the function c(, ) (Trampert and Woodhouse, 1995).
The solution to the wave equation is the displacement field, u(, , t), i.e., a seismogram
for a fixed surface point, or a snapshot of the distorted spherical membrane for a fixed time.
Comparing the seismograms from the homogeneous earth case with those from the hetero-
geneous case, we calculate the phase anomaly P and amplitude anomaly A, where P and
A are the values that minimize the residuals between the two seismograms for the first four
arrivals, orbits R1R4. This gives us our basis for comparison with the phase and amplitude
anomalies calculated in ray theory.
Ray theory predictions
With surface wave ray theory, the rays, confined to the spherical surface, are the orthogonal
trajectories of the wavefront of the surface wave. Ray theory can be used to calculate the
amplitude anomalies and phase anomalies for a given heterogeneous earth model (with respect
to a homogeneous earth model). We examine the results from ray theory using exact ray-
tracing as well as a linearized ray theory, in which the calculations are carried out along the
great-circle path between the source and receiver.
We can identify regimes under which ray theory predictions are inaccurate. The regime
can be changed by any combination of the following parameters: lmax, the truncation degree
of the spherical harmonic expansion of the phase velocity map c(, ); , the power of the
map; T , the period of the map; and , the path length of the ray, which has no upper limit
for the purposes of ray tracing.
We find that multipathing multiple ray paths between a source and receiver is quite
common for the current resolution of long-period surface wave phase velocity maps (e.g.,
iv
Trampert and Woodhouse, 1995). In particular, for paths beyond R1 for a map expanded in
spherical harmonics up to lmax=20, we would expect at least three paths between the source
and receiver; for R4 ray paths in an lmax=40 map, we find cases with up to 15 ray paths.
Multipathing does not occur in linearized ray theory, which assumes that the ray deviates
only slightly from great-circle path between the source and the receiver. Thus, multipathing
gives an indicating of the divergence of linearized ray theory from exact ray theory.
A comparison of three models
The numerical model provides us with our ground-truth estimate of the phase and am-
plitude anomalies, as it is the only model that calculates wave propagation. Furthermore,
it incorporates properties of a specified, albeit simplified, source. Ray theory assumes an
infinitesimally-short-wavelength wave, and it is this approximation which leads to deviations
from the actual picture of wave propagation. The accuracy of the numerical model is limited
by the number of gridpoints per wavelength (gpw): too few gpw results in unwanted numerical
dispersion.
Here we will state three conclusions:
1. The results from the numerical model phase anomalies and amplitude anomalies
correlate well with the results from ray theory, and this correlation is best when ,
where is the wavelength of the waves in the numerical model and is the minimum
scalelength of the heterogeneity in the phase velocity map. We see that exact ray-tracing
is a better predictor of phase and amplitude than linearized ray theory, based on the
comparison with lmax 8.
2. The phase anomalies predicted from ray theory are much better than the predictions of
amplitudes, which are more sensitive to shorter-scalelength structure. Even in regimes
where > , the phase predictions are well-correlated between the numerical model
and ray theory. This is one reason why phase remains the primary seismic data in
constructing velocity models of the earth.
3. In regimes of extreme multipathing (i.e., R4 orbits, T =40s, lmax=20), exact ray the-
ory predicts a large spread phase and amplitude anomalies, as each ray path carries
information on the distorted wave traveling between source and receiver. It is possible
v
that linearized ray theory is better in these extreme cases, or perhaps there is a way of
extracting one ray or combining the information from the rays get a single value. This
warrants further investigation.
We show that, with increased computing power, we would be able to input shorter and
shorter waves into the numerical model, thereby allowing us to see the waves response
to all of the structure in the more-detailed phase velocity maps. Thus, we could test
ray theory in the regime where multipathing is most prevalent, and where we expect
ray theory to fail.
vi
Acknowledgements
I would like to thank my advisor, John Woodhouse, who provided guidance throughout my
studies at Oxford. The basic idea behind the project was suggested by him, and he was instru-
mental in helping to develop the analytical solutions (Chapter 3) and the more generalized
ray-tracing equations (Appendix E). It has been a pleasure working with him over the past
two years, and I could not have asked for a better introduction into the field of seismology.
I would like to thank the Rhodes Trust at Oxford for making the experience possible from
a financial standpoint. I am also grateful for their financial support in attending meetings at
AGU 2002 in San Francisco and GSA 2001 in Boston. Worcester College provided support
for both of these meetings as well.
I have benefitted from discussions with Ana Ferreira, and I look forward to reading her
thesis in two years time (Surface Waves in a Laterally Varying Earth). She also generated
the data for the surface wave dispersion figures in Appendix F. Gareth Funning helped in
downloading the IRIS data for Figure 1.3. I would also like to thank Paul Davis for helping
me get started with wave propagation in Matlab, and my dad Walt for helping to clarify
some of the details in Appendix A. Also, I have benefited from discussions on finite-difference
methods with Nick Trefethen at Oxford, as well as from his two-semester course on scientific
computing. Finally, thanks to my MSc examiners, Jeannot Trampert at Utrecht and Barry
Parsons at Oxford, whose critical reading helped to improve the thesis.
Contents
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . i
Extended abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi
Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi
1 Introduction 1
2 The two-dimensional wave equation 8

2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Waves on a rectangular membrane . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Waves on a circular membrane . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.4 Numerical solutions of waves on a planar membrane . . . . . . . . . . . . . . 12
2.5 Waves on a spherical membrane . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3 Analytical spherical membrane waves 18

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.2 The solution for given forcing f = f 1 (, , t) . . . . . . . . . . . . . . . . . . . 21

3.3 The solution for given forcing f = t f1 (, , t) . . . . . . . . . . . . . . . . . 24
3.4 The solution for given initial shape g(, ) . . . . . . . . . . . . . . . . . . . . 25
3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4 Spherical hexagonal grids 34

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.2 Constructing the hexagonal grids . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.3 Properties of the hexagonal grids . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.4 Calculating the Laplacian on the hexagonal grids . . . . . . . . . . . . . . . . 42
4.5 Error in the numerical approximation of the Laplacian . . . . . . . . . . . . . 44
4.6 Order of accuracy for the numerical approximation . . . . . . . . . . . . . . . 46
4.7 Variations and optimizations on the hexagonal grids . . . . . . . . . . . . . . 47
4.8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5 Numerical spherical membrane waves 60

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.2 The phase velocity function c(, ) . . . . . . . . . . . . . . . . . . . . . . . . 61
5.3 Numerical iteration for wave propagation . . . . . . . . . . . . . . . . . . . . 64
5.4 Numerical solutions for homogeneous and heterogeneous phase velocity maps 64
vii
CONTENTS viii
5.5 Stability and accuracy in the numerical model . . . . . . . . . . . . . . . . . . 66
6 Ray calculations on the sphere 79

6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
6.2 Surface wave ray-tracing equations . . . . . . . . . . . . . . . . . . . . . . . . 81
6.3 Path anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
6.4 Phase anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
6.5 Amplitude anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.6 Linearized ray theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
6.7 Exact ray-tracing: Ray diagram cross-section curves . . . . . . . . . . . . . . 91
6.8 Exact ray-tracing: Variation with , lmax, , and T . . . . . . . . . . . . . . 93
7 Numerical model vs ray theory 113

7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
7.2 Exact ray theory vs linearized ray theory . . . . . . . . . . . . . . . . . . . . 115
7.3 Obtaining phase and amplitude anomalies from the numerical model . . . . . 117
7.4 Ray theory validity condition . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
7.5 Numerical model vs ray theory . . . . . . . . . . . . . . . . . . . . . . . . . . 119
7.6 Summary and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
A Functions on the sphere 138
B Spherical harmonic functions 144

B.1 Spherical harmonic conventions . . . . . . . . . . . . . . . . . . . . . . . . . . 144
B.2 Spherical harmonic conventions in Matlab . . . . . . . . . . . . . . . . . . . . 146
B.3 Spherical harmonic expansions . . . . . . . . . . . . . . . . . . . . . . . . . . 148
C Euler angles and rotations 151
D Derivations for selected equations 155

D.1 Equation (3.12) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
D.2 Equations (3.27) and (3.33) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
D.3 Equations (3.20) and (3.62) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
D.4 Some definitions and properties . . . . . . . . . . . . . . . . . . . . . . . . . . 159
E Ray-tracing equations 161

E.1 Receiver at a fixed longitude . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
E.2 Equation (6.19) ds/d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167
E.3 Equation (6.11) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
F Wave parameters and dispersion 170
G Flow chart of computer programs 172
Bibliography 174
List of Figures
1.1 Rayleigh waves and Love waves . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.2 Global Rayleigh waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3 Alaska earthquake of November 2002 . . . . . . . . . . . . . . . . . . . . . . . 7
2.1 Rectangular and circular membranes . . . . . . . . . . . . . . . . . . . . . . . 9

2.2 Waves on a rectangular membrane . . . . . . . . . . . . . . . . . . . . . . . . 15
2.3 Waves on a square membrane . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4 Waves on a circular membrane . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.1 Initial-shape function g(, ) and source-time function h(t) . . . . . . . . . . 31

3.2 The width of the Gaussian initial-shape function . . . . . . . . . . . . . . . . 32
3.3 Surface plot of Il () . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.4 Analytical solutions for wave propagation on the sphere . . . . . . . . . . . . 33
4.1 Discretizations of the sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

4.2 Sketch of integration parameters . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.3 Parameter for optimizing the gridpoint spacing . . . . . . . . . . . . . . . . . 48
4.4 Relationship between successive orders of triangular grids . . . . . . . . . . . 51
4.5 Construction of the triangular grids . . . . . . . . . . . . . . . . . . . . . . . 51
4.6 The order q=0 triangular grid . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.7 Construction of a hexagonal grid from a triangular grid . . . . . . . . . . . . 52
4.8 Spherical hexagonal grids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.9 Geometric properties of the hexagonal grids, I . . . . . . . . . . . . . . . . . . 54
4.10 Geometric properties of the hexagonal grids, II . . . . . . . . . . . . . . . . . 55
4.11 Numerical method for calculating the Laplacian at a gridpoint . . . . . . . . 55
4.12 Error in the numerical approximation of the Laplacian, I . . . . . . . . . . . . 56
4.13 Error in the numerical approximation of the Laplacian, II . . . . . . . . . . . 57
4.14 Error in the numerical approximation of the Laplacian, III . . . . . . . . . . . 58
4.15 Order of accuracy of the numerical approximation of the Laplacian . . . . . . 59
5.1 Phase velocity map: Homogeneous earth . . . . . . . . . . . . . . . . . . . . . 62

5.2 Surface wave phase velocity map (vary T ) . . . . . . . . . . . . . . . . . . . . 69
5.3 Phase velocity maps: Vary lmax . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.4 Phase velocity maps: Vary . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
5.5 Amplitude spectra for phase velocity maps . . . . . . . . . . . . . . . . . . . . 72
5.6 Phase velocity maps: TW95 vs TW96 . . . . . . . . . . . . . . . . . . . . . . 73
5.7 Seismograms for a homogeneous velocity map . . . . . . . . . . . . . . . . . . 74
ix
LIST OF FIGURES x
5.8 Seismograms for a heterogeneous velocity map . . . . . . . . . . . . . . . . . 75

5.9 Waveforms for a heterogeneous phase velocity map . . . . . . . . . . . . . . . 76
5.10 Wavefront propagating on a heterogeneous earth . . . . . . . . . . . . . . . . 77
5.11 Variation in seismograms with varying lmax . . . . . . . . . . . . . . . . . . . 78
6.1 Source and receivers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

6.2 Great-circle ray paths . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
6.3 Introduction to ray-tracing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
6.4 Phase velocity map and transverse derivatives . . . . . . . . . . . . . . . . . . 99
6.5 Ray diagram plus cross-sections . . . . . . . . . . . . . . . . . . . . . . . . . . 101
6.6 Ray diagram cross-sections . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
6.7 Path deviation from great circles. . . . . . . . . . . . . . . . . . . . . . . . . . 103
6.8 Ray diagram cross-sections: Variation with lmax, , and T . . . . . . . . . . 104
6.9 Zoom-in on Figure 6.8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
6.10 Variation in amplitude anomaly with lmax, , and T . . . . . . . . . . . . . . 106
6.11 Variation in phase anomaly with lmax, , and T . . . . . . . . . . . . . . . . 107
6.12 Ray diagram cross-sections: TW95 vs TW96 . . . . . . . . . . . . . . . . . . 108
7.1 Flowchart of seismogram generation . . . . . . . . . . . . . . . . . . . . . . . 114

7.2 Exact ray theory vs linearized ray theory: Vary lmax . . . . . . . . . . . . . 124
7.3 Exact ray theory vs linearized ray theory: Vary . . . . . . . . . . . . . . . . 125
7.4 Exact ray theory vs linearized ray theory: Vary T . . . . . . . . . . . . . . . 126
7.5 Exact ray theory vs linearized ray theory: Vary orbit . . . . . . . . . . . . . . 127
7.6 Obtaining the amplitude and phase anomalies from the numerical model . . . 128
7.7 Wavelength as a function of source width . . . . . . . . . . . . . . . . . . . . 129
7.8 Numerical results vs ray theory: lmax=6 . . . . . . . . . . . . . . . . . . . . 130
7.9 Numerical results vs ray theory: lmax=12 . . . . . . . . . . . . . . . . . . . . 131
7.10 Numerical results vs exact ray theory: Vary lmax . . . . . . . . . . . . . . . . 132
7.11 Numerical results vs exact ray theory: Vary . . . . . . . . . . . . . . . . . . 133
7.12 Numerical results vs exact ray theory: Vary T . . . . . . . . . . . . . . . . . . 134
7.13 Numerical results vs exact ray theory: Vary . . . . . . . . . . . . . . . . . . 135
7.14 Numerical results vs exact ray theory: Vary orbit . . . . . . . . . . . . . . . . 136
7.15 Shorter wavelengths versus shorter scalelengths . . . . . . . . . . . . . . . . . 137
A.1 Function mappings for spherical coordinates . . . . . . . . . . . . . . . . . . . 138

A.2 Extending ge(, ) to the - plane, I . . . . . . . . . . . . . . . . . . . . . . . 142
A.3 Extending ge(, ) to the - plane, II . . . . . . . . . . . . . . . . . . . . . . . 143
A.4 Extending ge(, ) to the - plane, III . . . . . . . . . . . . . . . . . . . . . . 143
B.1 Sample functions of Ylm (, ) . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
C.1 Euler angle rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153

C.2 Euler angle rotation (in reverse) . . . . . . . . . . . . . . . . . . . . . . . . . . 154
E.1 Sketches for ray-tracing equations . . . . . . . . . . . . . . . . . . . . . . . . . 169
F.1 Dispersion of surface waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171
G.1 Flowchart of computer programs . . . . . . . . . . . . . . . . . . . . . . . . . 173

List of Tables
4.1 Spherical hexagonal grids used in other studies . . . . . . . . . . . . . . . . . 36

4.2 The five platonic solids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.3 Geometric properties of the spherical grids . . . . . . . . . . . . . . . . . . . . 40
4.4 Properties of the hexagonal grids used in this study . . . . . . . . . . . . . . 40
4.5 Properties of the hexagonal grids of Heikes and Randall (1995a) . . . . . . . 41
4.6 Subdivided cubic grids used in other studies . . . . . . . . . . . . . . . . . . . 47
4.7 Optimizing the hexagonal grids . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.1 Properties of phase velocity maps of Trampert and Woodhouse (1995) . . . . 68

5.2 Minimum scalelength in a phase velocity map . . . . . . . . . . . . . . . . . . 68
6.1 Phase anomaly convention . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

6.2 Great-circle path lengths for global surface waves . . . . . . . . . . . . . . . . 97
6.3 Dependence of multipathing on orbit and lmax . . . . . . . . . . . . . . . . . 109
6.4 Dependence of multipathing on and T . . . . . . . . . . . . . . . . . . . . . 111
6.5 Dependence of multipathing on TW95 vs TW96 . . . . . . . . . . . . . . . . 112
F.1 Wave parameters and some related equations . . . . . . . . . . . . . . . . . . 170
xi
Chapter 1
Introduction
This is a study in global seismology, and because is leans on the theoretical side, it is important
to establish its context in the earth sciences before going into any details. The basic goal of
global seismology is to improve our physical picture of the earths interior, from the center
to the surface. Our pictures of the earths interior tomographic images are expressed
in terms of seismic velocities, which serve as proxies for physical properties of the material
through which the seismic waves travel 1 . The tomographic images (e.g., Ritsema et al., 1999)
give us a snapshot of what lies deep beneath our feet, and these snapshots provide the best
starting point for discussion of topics such as mantle convection and plate-driving mechanisms
(e.g., Albarede and van der Hilst, 2002). Thus an important goal of global seismology is to
map the earths interior with as much detail as possible, thereby establishing the images
essential for discussing and modelling the dynamic processes of the earths interior.
Surface waves provide the best constraints on the structure 2 of one of the most interesting
regions of the earth, the upper mantle, which extends to a depth of approximately 700 km.
Long-period surface waves periods approximately in the range T = 20 to 300s sample
the upper mantle, hence their label as mantle waves. These waves provide insight into
1
In particular, the expressions for the compressional velocity (vp ) and shear velocity (vs ) in an isotropic
solid are given by:
s
+ 34
r

vp = vs = ,

where is the shear modulus, is the bulk modulus, and is the density.
2
By detail of structure, we mean the detail of the seismic velocity coverage. The detail in velocity struc-
ture could be a manifestation of the thermal structure or chemical structure (or a combination). The term
structure in this thesis will generally refer to the structure of the phase velocity variation.
1
CHAPTER 1. Introduction 2
the lateral heterogeneity and anisotropy within the upper mantle. Lateral heterogeneity is
helpful in characterizing subduction zones, plumes, and the fate of subducted material, while
anisotropy in the upper mantle is one of the indicators of mantle flow direction, which is
considered to be a principal driver of the plates.
There are two types of surface waves: Rayleigh waves, which propagate with vertical
and radial displacements, and Love waves, which propagate with transverse displacements
only (Figure 1.1). After large earthquakes, such as the one in central Alaska in November
2002 (M7.9), the surface waves may be detected traveling around the earth multiple times.
Figure 1.2 shows a theoretical example of these global surface waves, and Figure 1.3 shows
that these surface waves are readily identifiable in seismograms. It is data such as these
to which theories in seismology are applied in order to construct the images of the earths
interior.
Our three-dimensional picture of the earth has improved drastically over the past 25
years. The development of the Preliminary Reference Earth Model (PREM), a 1D earth
model that remains the standard, marked a starting point for tomography (Dziewonski and
Anderson, 1981). Woodhouse and Dziewonski (1984) produced the first 3D global images of
the Earths interior, and today there exist several 3D models with greater detail of structure
(e.g., Ritsema et al., 1999). Many models of the upper mantle rely on surface wave data
(an exception is Grand (2002)), and some are constructed from surface wave data alone 3
(Boschi and Ekstr
om, 2002). The CaltechOxford model (Ritsema et al., 1999), one of the
most accepted global models today, incorporates surface waves, including higher-mode surface
wave arrivals. Improvements in the global models have come with improvements in theory
and with improvements in data coverage and data quality. In this study, we look at one of
the most commonly used tools in seismology, ray theory.
The purpose of this study is to develop a numerical model for wave propagation on a
spherical membrane and to use the results to check predictions made from ray theory. We use
this model as an analogue for surface wave propagation, in order to investigate the properties
of ray theory for wave propagation on a sphere. Ray theory has been used in seismology
to explain observations of path differences, amplitude differences, and phase differences of
surface waves with respect to a spherical reference model (e.g., Woodhouse and Wong, 1986;
3
A global model of the crustal thickness down to the Moho is crucial in these models (Bassin et al., 2000;
Mooney et al., 1998).
Laske et al., 1994). Furthermore, it is a standard tool used in creating the global models
of the earth. However, ray theory is an approximate theory and cannot always be used to
explain observations. Thus, numerical models are helpful, especially as earth models become
increasingly more detailed. This study presents a two-dimensional numerical model aimed at
investigating the regimes under which ray theory predictions fail.
Here we will emphasize some key points, and then outline the remainder of this report.
1. This study is divided into two approaches to examining wave propagation on a spherical
membrane: (1) a numerical model, in which the wave propagates on the spherical
membrane according to the wave equation, and (2) surface wave ray theory, whereby the
rays, confined to the spherical surface, are the orthogonal trajectories of the wavefront
of the surface wave. The numerical solution provides the amplitude information in
the form of a global snapshot at each timestep, while ray theory provides amplitude
and phase information for a particular source-and-receiver pair. Our goal is to use the
numerical model to examine the limitations of ray theory.
2. Our problem is mathematically a two-dimensional problem, as the term membrane

suggests. The propagation of a Rayleigh wave induces displacements down to depths of
hundreds of km, but we are only concerned with the displacement of the surface, which
can be adequately modelled by the wave equation. Thus our membrane waves are an
analogue for the surface waves. The surface wave phase velocity maps c(, ) that we
input into the numerical model and into the ray tracing include information about the
3D structure of the upper mantle.
3. The mathematics of our problem is stripped of the complexity of wave propagation

that we know to exist in response to the 3D, laterally heterogeneous, anisotropic, and
anelastic structure of the upper mantle. In nature these effects give rise to dispersion,
scattering, reflection, anisotropy, and attenuation of the surface waves 4 . Our wave
equation simply includes a forcing term f (, , t) for the source and a phase velocity
map of the earth c(, ).
4
We do, however, encounter several of these effects in our numerical model as a result of discretizing a
continuous problem; numerical dispersion is the most predictable artifact of the numerical solution.
In this thesis we first look at analytical and numerical solutions to waves on a planar
membrane (Chapter 2), which provides context for examining analytical solutions to waves
on a spherical membrane (Chapter 3). In order to calculate numerical solutions to waves on a
spherical membrane, we need a spherical grid that will allow us to numerically approximate the
Laplacian (2 ) at each time-step. After constructing and testing the grid (Chapter 4), we then
compare the numerical solutions of waves on the spherical membrane with analytical solutions
for the homogeneous case (i.e., phase velocity constant everywhere) (Chapter 5). We then find
numerical solutions for the heterogeneous case, using phase velocity maps of Trampert and
Woodhouse (1995) (Chapter 5). Lastly, we use ray theory to predict the phase and amplitude
variation of the waves on the spherical membrane (Chapter 6), and we investigate how these
predictions compare with those from the numerical model (Chapter 7).
Rayleigh Wave
vertical
radial
e
rs
s ve
an
r
t
Love Wave
vertical
radial
e
rs
s ve
an
r
t
Figure 1.1: Ground displacements for Rayleigh waves and Love waves (after Shearer , 1999).
The direction of propagation in each case is in the radial direction. Displacement of Rayleigh
waves is restricted to the vertical and radial directions; displacement of Love waves is restricted
to the transverse direction only. These surface waves travel 45 km/s see Appendix F for
more details.
Global Surface Waves it s eve

orb n
or
dd
bi
o
source
ts
o
at latitude +90
R1 (t = 3000),
receiver at
R3 (t = 13000), o
latitude -30
R5, etc.
R2 (t = 7000),
R4 (t = 17000),
R6, etc.
m
aj
rc
or
a
or
ar
min
R1 R2 R3 R4
( = 120o) ( = -240o) ( = 480o) ( = -600o)
lat = -30.00, lon = 0.00 lat = -90.00, lon = 0.00

600 600 R1+R2 R3+R4
400 400
Displacement
Displacement
R1 R2 R3 R4
200 200
0 0
-200 -200
10,000 s
-400 -400
-600 -600
0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
Time (s) x 104 Time (s) x 104
Figure 1.2: Global surface waves. Schematic diagram illustrating the nomenclature of global
surface wave orbits, depicted for a source at the North Pole and a receiver at lat=30 .
R denotes Rayleigh waves (G = Love waves), and each number corresponds to the relative
arrival of the orbit at a particular receiver (see theoretical seismogram). The time between
even or odd orbits is the time it takes for the wave to propagate around the earth (10,000
seconds, or 2.8 hours). Inset shows a seismogram for the special case of the receiver located
at the antipode, the point on the earth directly opposite the source. In that case, the arrivals
constructively interfere, and the observer experiences two pulses, rather than four.
90o R1, R180 90o o

0o 60o 120o 3 o -120o -60o 0
A
60o 60o
30o 30o
R2 , R
0o 0o
4
-30o -30o
-60o -60o
-90o -90o
0o 60o 120o 180o -120o -60o 0o
x 105 First Six Hours of Denali quake: Unfiltered, Vertical Component
5
B
Amplitude
-5
1 2 3 4 5 6
P S G1+R1 R2 R3
Time (hours)
x 104 First Six Hours of Denali quake: Filtered T=135-150s, Vertical Component
2
about 3 hours
C
Amplitude
R3 R4
-2 R2
R1
1 2 3 4 5 6
Time (hours)
Figure 1.3: The M7.9 Alaska earthquake of November 3, 2002 (lat=63.42 , lon=144.96 ).
Here we show the event recorded in Borovoye, Kazakhstan (BRVK) (lat=53.06 , lon=70.28 ).
A. Great circle path containing the source (?) and receiver (2); also shown is the pole to the
great circle (). B. Unfiltered record of the first six hours. Note the relative arrival times and
relative amplitudes of the P -wave, the S-wave, and the surface waves (R1, G1). C. Filtered
record of the first six hours. The long-period Rayleigh waves are recorded orbiting the earth
in both directions more than once (see Figure 1.2). The time for these elastic waves to circle
the earth is about three hours. The great-circle path lengths are as follows: R1 = 60.32 ,
R2 = 299.68 , R3 = 420.32 , R4 = 659.68 .
Chapter 2
The two-dimensional wave equation
Contents
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Waves on a rectangular membrane . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Waves on a circular membrane . . . . . . . . . . . . . . . . . . . . . . . . 10
2.4 Numerical solutions of waves on a planar membrane . . . . . . . . . . . 12
2.5 Waves on a spherical membrane . . . . . . . . . . . . . . . . . . . . . . . 14
2.1 Introduction
The wave equation is given by
1 2u
= 2 u , (2.1)
c2 t2
where c is the phase velocity of the wave field u. In the homogeneous case considered in this
chapter, c is a constant.
In the following sections we will examine the solutions to Equation (2.1) for waves (i)
on a rectangular membrane and (ii) on a circular membrane (Figure 2.1). The boundary
condition in these cases is that the edge of the membrane is fixed for all t 0. These cases
will provide context for studying waves on a spherical membrane. Details of the following
analytical solutions for waves on rectangular and circular membranes can be found in Kreyszig
(1988) and references therein.
8
CHAPTER 2. The two-dimensional wave equation 9
A y B y
b
Figure 2.1: Rectangular
R
x
and circular membranes for
x
a the two-dimensional wave
equation.
2.2 Waves on a rectangular membrane
The two-dimensional wave equation for waves in the xy-plane is
1 2u 2u 2u
= + 2 , (2.2)
c2 t2 x2 y
where u(x, y, t) is the displacement in the vertical direction. Our boundary condition is u = 0
for t 0 on the boundary of a length a-by-b rectangular membrane (Figure 2.1A). Our initial
conditions consist of a given initial displacement and initial velocity:
u(x, y, 0) = g(x, y) (2.3)

u
= v(x, y) . (2.4)
t t=0
The general solution is the two-dimensional Fourier sine series

X
X
u(x, y, t) = umn (x, y, t) (2.5)
m=1 n=1
X X mx ny

= (Bmn cos mn t + Bmn sin mn t) sin sin , (2.6)
m=1 n=1
a b
where
r
m2 n2
mn = c + 2 . (2.7)
a2 b
Using the initial conditionsEquations (2.3) and (2.4)in Equation (2.6), we obtain
X
X mx ny
u(x, y, 0) = Bmn sin sin = g(x, y) (2.8)
a b
m=1 n=1
X mx ny
u X

= Bmn mn sin sin = v(x, y) . (2.9)
t t=0 m=1 n=1
a b
Supposing that g(x, y) and v(x, y) can each be developed in such a series, Equations (2.8)

and (2.9) can be used to determine the coefficients B mn and Bmn in Equation (2.6) (e.g.,
Kreyszig, 1988):
Z bZ a mx ny
4
Bmn = g(x, y) sin sin dx dy (2.10)
ab 0 0 a b
Z bZ a mx ny
4
Bmn = v(x, y) sin sin dx dy . (2.11)
a b mn 0 0 a b
We can plot Equation (2.6) once we specify the initial conditions g(x, y) and v(x, y). Let
us choose an initial shape as a Gaussian at the center of the (a-by-b) rectangular membrane,
and choose an initial velocity of zero:
1
u(x, y, 0) = g(x, y) = 2 e((xa/2) + (yb/2) ) /
2 2 2
(2.12)

u
= v(x, y) = 0 , (2.13)
t t=0
where is included as a normalization constant. Plots of Equation (2.6) for three different
times are shown in Figure 2.2. If we use a square membrane instead of a rectangle, then we
see a more symmetric reflection of the wave for the initial Gaussian (Figure 2.3).
2.3 Waves on a circular membrane
We now look at the case with a circular membrane of radius R (Figure 2.1B). Using the
polar-coordinate form of the Laplacian, the wave equation is
1 2u 2u 1 u 1 2u
= + + . (2.14)
c2 t2 r 2 r r r 2 2
Considering solutions u(r, t) that are radially symmetric (i.e., no -dependence), Equation (2.14)
reduces to
1 2u 2u 1 u
2 2
= + . (2.15)
c t r 2 r r
With the membrane fixed at r = R, we have the boundary condition u(R, t) = 0 for t 0.
The initial position and velocity do not depend on and are given by
u(r, 0) = g(r) (2.16)

u
= v(r) . (2.17)
t t=0
The general solution is the series

X
u(r, t) = um (r, t) (2.18)
m=1
X r
m
= (am cos m t + bm sin m t) J0 , (2.19)
R
m=1
where m = c m /R, and m are the positive (real) roots of J0 (x), which is the Bessel function
of the first kind of order zero:
(1)j x 2j

X x2 x4 x6
J0 (x) = 2
= 1 2 + 2 2 2 2 2 2 2 + . (2.20)
(j!) 2 2 2 4 6 2 4 6 8
j=0
Using the initial conditionsEquations (2.16) and (2.17)in Equation (2.19), we obtain

X r
m
u(r, 0) = am J0 = g(r) (2.21)
R
m=1
X r
u m
= bm m J0 = v(r) . (2.22)
t t=0 R
m=1
Supposing that g(r) and v(r) can each be developed in such a series, Equations (2.21) and
(2.22) can be used to determine the coefficients a m and bm of the Fourier-Bessel series in
Equation (2.19) (Kreyszig, 1988):

Z R r
2 m
am = r g(r) J0 dr (2.23)
R (J1 (m ))2
2
0 R
Z R r
2 m
bm = r v(r) J0 dr . (2.24)
c m R (J1 (m ))2 0 R
We can plot Equation (2.19) once we specify the initial conditions g(r) and v(r). Let us
choose an initial shape as a Gaussian at the center of the circular membrane, and choose an
initial velocity of zero:
1 2 2
u(r, 0) = g(r) = 2 er / (2.25)

u
= v(r) = 0 , (2.26)
t t=0
where is included as a normalization constant. Plots of Equation (2.19) for three different
times are shown in Figure 2.4.
2.4 Numerical solutions of waves on a planar membrane
We now examine numerical solutions of the wave equation on a planar membrane. This
will serve as a basis for examining numerical solutions for waves on a spherical membrane
(Section 5.3). Here we examine the case of a rectangular membrane, and include a forcing
term in the two-dimensional wave equation:
1 2u
2 u = f , (2.27)
c2 t2
where f = f (x, y, t) is the forcing term, u = u(x, y, t) is the wave field, represented as a vertical
displacement, and 2 u = 2 u/x2 + 2 u/y 2 . The solutions are subject to the boundary
condition u = 0 on the boundary of the membrane for all t 0:
u(x, 0, t) = u(x, b, t) = 0 , 0xa (2.28)
u(0, y, t) = u(a, y, t) = 0 , 0yb ,
where a and b describe the boundary of the rectangular membrane (Figure 2.1A).
The wave equation represents a continuous problem in space and time. In order to obtain
numerical solutions, we must discretize the problem in both space and time. This involves
discretizing the differential equations (and the functions within), as well as the boundary
conditions. Here we use a finite-difference method in space and time to solve the 2D wave
equation. Using the second-order leap-frog finite-difference approximation for u
, u xx , and
uyy , Equation (2.27) becomes

= c 2 2 u + f
u (2.29)
= c(x, y)2 (uxx + uyy + f (x, y, t) )

u (2.30)
!
ui,j i,j i,j
t+1 2ut + ut1
i,j 2 ui+1,j 2ui,j i1,j
t + ut ui,j+1 2ui,j i,j1
t + ut
= c t
+ t
+ fti,j ,
dt2 dx 2
dx 2
(2.31)
where ui,j 1
t is the displacement at each gridpoint, t is the time increment index , i is the spatial
increment index in the x-direction, j is the spatial increment index in the y-direction, dt is the
time increment, and dx is the space increment 2 (e.g., Trefethen, 1996). Thus, Equation (2.30)
is the continuous form of the wave equation, and Equation (2.31) is our discrete form. For
simplicity, we shall substitute del2 for the discrete finite difference approximation of the
Laplacian operator in the xy-plane:
ui,j i,j i,j

t+1 2ut + ut1 2
= ci,j del2 ui,j
t + f i,j
t , (2.32)
dt2
which, solving for ut+1 , gives
ut+1 = 2 ut ut1 + c2 dt2 (del2(ut ) + ft ) , (2.33)
where we have eliminated the i, j superscripts for clarity.

1
Note that the subscripts in the index notation do not refer to partial differentiation.
2
The spatial increment dx is the same in both directions, even though we have a rectangular domain. We
could put a dy in the second denominator of Figure 2.31 to indicate the direction of the spatial step.
We used Equation (2.33) in a program to produce numerical solutions to the wave equation
on a rectangular membrane. First we specify the forcing term f = f (x, y, t) and the initial
conditions ut1 = ut = g(x, y). For each time-step, the iteration is:
ut1 = ut (2.34)
ut = ut+1 (2.35)
ut+1 = 2 ut ut1 + c2 dt2 (del2(ut ) + ft ) (2.36)
ut |boundary = 0. (2.37)
With initial conditions given by Equation (2.12) and no forcing term (f (x, y, t) = 0) and
c(x, y) = const, the numerical solutions for the rectangular membrane waves agreed with
those for the analytical case (i.e., Figures 2.2 and 2.3).
It should be noted that there are various ways of obtaining numerical solutions to the
wave equation. An alternative approach to using finite difference methods is to use spectral
methods (or pseudo-spectral methods). An excellent treatment of these two methods as
applied to the wave equation can be found in works by Trefethen (Trefethen, 1982, 1996,
2000). A third approach it to use finite element methods.
2.5 Waves on a spherical membrane
The work of Tanimoto (1990) was one of the motivations for this study. Tanimoto (1990)
modelled waves on a planar membrane using a second-order finite-difference method on a
region bounded by constant latitudes and longitudes. He first looked at surface wave prop-
agation through California, comparing his numerical results with ray theory predictions and
with seismograms for the Whittier Narrows earthquake of 01 October 1987.
With the California example, the area was small enough so that the curvature of the Earth
was neglected. In order to consider wave propagation on a global scale, Tanimoto (1990)
solved the wave equation on the sphere and incorporated absorbing boundary conditions
(Clayton and Engquist, 1977), which work fairly well up to about 75 [latitude] (p. 100).
Our goal was to model waves on a spherical membrane without using a Cartesian grid, which
always requires special boundary conditions for high latitudes.
time = 0.00
g (x, y)
membrane fixed
at y = 1 membrane fixed
at x = 2
membrane fixed
at x = 0 membrane fixed
at y = 0
time = 0.50
time = 1.75
Figure 2.2: Analytical solutions of waves on a rectangular membrane. Three time frames for
the initial conditions of (i) an initial shape of a Gaussian (Eq. 2.12 with = 0.1; top figure)
and (ii) an initial velocity of zero (Eq. 2.13). These are plots of Eq. 2.6 with c = 1, summing
the modes up to m = 20 and n = 20. The boundary is fixed.
time = 0.00 g (x, y)
membrane fixed
at y = 1 membrane fixed
at x = 1
membrane fixed
at x = 0 membrane fixed
at y = 0
time = 0.50
time = 1.75
Figure 2.3: Analytical solutions of waves on a square membrane. The only difference from
Figure 2.2 is that we have changed the width of the membrane to from a = 2 to a = 1. The
initial Gaussian wave has a more symmetric reflection off the square boundary here than off
the rectangular boundary in Figure 2.2, as expected.
time = 0.00
g (r, )
membrane fixed
at r = 1.0
time = 0.40
time = 1.20
Figure 2.4: Analytical solutions of waves on a circular membrane. Three time frames for the
initial conditions of (i) an initial shape of a Gaussian (Eq. 2.25 with = 0.1; top figure) and
(ii) an initial velocity of zero (Eq. 2.26). These are plots of Eq. 2.19 with c = 1, summing
the modes up to m = 20. At t = 1.20, the wave has just reflected at the boundary and is
returning on the underside. The boundary is fixed.
Chapter 3
Analytical solutions of waves on a

spherical membrane
Contents
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.2 The solution for given forcing f = f1 (, , t) . . . . . . . . . . . . . . . . . 21

3.3 The solution for given forcing f = t f1 (, , t) . . . . . . . . . . . . . . . . 24
3.4 The solution for given initial shape g(, ) . . . . . . . . . . . . . . . . . . 25
3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.1 Introduction
This chapter presents analytical solutions to the wave equation on a spherical membrane.
Some of the more important derivations are shown, while others are included in Appendix D.
Appendix A examines the issue of functions being well-defined on the sphere, and Appendix B
discusses spherical harmonic functions.
We wish to solve the following equation for u = u(, , t):
1 2u
2 u = f , (3.1)
c2 t2
where u is the wave field, c = c(, ) is the frequency-independent, constant phase velocity,
18
CHAPTER 3. Analytical spherical membrane waves 19
f is the forcing term, and

2 1 1 u 1 2u
u = 2 sin + (3.2)
a sin sin2 2
is the surface Laplacian of u on a sphere with radius a, where , are spherical coordinates
( is colatitude). In our results, we plot the wave field u as a radial displacement on the
sphere. It is assumed that the forcing term f = f (, , t) is a given function, vanishing at
(sufficiently large) negative times t. In this form, Equation (3.1) is the two-dimensional wave
equation on the sphere. Note that attenuation and dispersion are not considered.
To solve Equation (3.1) we first expand f and u in spherical harmonics :

X l
X
f = flm Ylm (, ) (3.3)
l=0 m=l

X l
X
u = ulm Ylm (, ) (3.4)
l=0 m=l
where flm = flm (t) and ulm = ulm (t) are given by
Z 2 Z
1
flm = Ylm (, ) f (, , t) sin d d (3.5)
Nl 0 0
Z 2 Z
1
ulm = Ylm (, ) u(, , t) sin d d , (3.6)
Nl 0 0
where Nl is the normalization integral (assumed independent of m), such that
Z 2 Z
0
Ylm
0 (, ) Ylm (, ) sin d d = Nl l0 l m0 m , (3.7)
0 0
where

1 if l = l0 and m = m0
l 0 l m0 m = (3.8)
0 otherwise .
Equation (3.7) is a statement of the orthogonality of spherical harmonic functions, which are
discussed in Appendix B. Substituting Equations (3.3) and (3.4) into Equation (3.1), we
obtain
l
! l
! l
1 2 X X X X X X
ulm Ylm (, ) 2
ulm Ylm (, ) = flm Ylm (, ) ,
c2 t2
l=0 m=l l=0 m=l l=0 m=l
(3.9)
or, for each spherical harmonic term individually,
1 d2 ulm m
Yl (, ) ulm 2 Ylm (, ) = flm Ylm (, )
c2 dt2
1 d2 ulm m l(l + 1) m
2 2
Yl (, ) + ulm Yl (, ) = flm Ylm (, )
c dt a2
1 d2 ulm l(l + 1)
2 2
+ ulm = flm
c dt a2
d2 ulm
+ l2 ulm = c2 flm , (3.10)
dt2
where
p
c l(l + 1)
l . (3.11)
a
The solution of Equation (3.10) can be obtained by the method of Greens functions, or by
Laplace or Fourier transforms. The solution is, assuming that u lm and dulm /dt, like flm ,
vanish at (sufficiently large) negative times (see Appendix D):
Z t
sin l (t )
ulm (t) = c2 flm ( ) d . (3.12)
l
Substituting Equation (3.12) into Equation (3.4) gives

X l
X Z t
2 sin l (t )
u = c flm ( ) d Ylm (, ) (3.13)
l
l=0 m=l
From Equation (3.5) we see that
Z 2 Z
1
flm ( ) = Ylm (, ) f (, , ) sin d d , (3.14)
Nl 0 0
and we combine Equations (3.13), (3.14), and (3.7) to obtain

Z 2 Z

X l
X 2
Z t Ylm (, ) f (, , ) sin d d
c sin l (t )
u = Z0 2 Z0 d m
Yl (, ) .
l
l=0 m=l Ylm (, ) Ylm (, ) sin d d
0 0
(3.15)
Equation (3.15) is our general expression for the wave field u(, , t) on the surface of the
sphere, given a forcing term f . In the next two sections we will examine the solutions for two
different forcings, f1 and f2 .
3.2 The solution for given forcing f = f1(, , t)
We assume that the forcing term f = f (, , t) in Equation (3.1) has the form:
f1 (, , t) = h1 (t) g(, ) , (3.16)
where
2 2
et /2
h1 (t) = (3.17)
2
is the source-time function and
2 2
e /2
g(, ) = (3.18)
2
is the initial-shape function (Figure 3.1). The denominators in Equations (3.17) and (3.18) are
included for convenience only, to facilitate the consideration of the properties of the solutions
as 0 (g ()) or 0 (h (t)) the point source solutions (see Section 3.5).
They are chosen by virtue of:
Z 2 2
!

et /2
dt = 1 (3.19)
2
Z 2 2
!

e /2
d = 1 . (3.20)
0 2
Intuitively, we would expect the wavelength of propagating the wave on the sphere to decrease
if either or decrease. (This is shown later in Figure 7.7).
It is somewhat easier to interpret changes in on g(, ) if we write as a function of
wid , the angular width of the Gaussian at a specified fraction r of its maximum amplitude.
If we set g(r , ) = r/2 and solve for in Equation (3.18), we get

= r
2 ln r

(wid , r) = wid , (3.21)
2 2 ln r
where we have set wid 2 r as shown in Figure 3.2. Thus, and wid can be thought of
interchangeably.
Because f1 (, , t) is independent of , the solution will be cylindrically symmetric; thus
flm = 0 for m 6= 0, and similarly ulm = 0 for m 6= 0. With m = 0 and no -dependence,
Equation (3.15) becomes
Z

X 2
Z t 2 Yl0 () f1 (, ) sin d
c sin l (t )
u = Z0 d 0
Yl () (3.22)
l
l=0 2 Yl0 () Yl0 () sin d
0
Writing Yl0 () = Yl0 () = Pl (cos ) and substituting the forcing term (Eq. 3.16), we get
Z 2 2

e /2
Z t Pl (cos ) h1 ( ) sin d
X 2 0 2 sin l (t )
u f1 = c Z d Pl (cos )
l
l=0 Pl (cos ) Pl (cos ) sin d
0
Z 2 2

e /2
Z t h1 ( ) Pl (cos ) sin d
X
2 2 sin l (t )
= c Z0 d Pl (cos )
2 l
l=0 (Pl (cos )) sin d
0
Z 2 2

e /2

X Pl (cos ) sin d Z t
2 0 2 sin l (t )
= c Z h1 ( ) d Pl (cos ) .
l
l=0 (Pl (cos ))2 sin d
0
(3.23)
We use the notation uf1 for this solution to emphasise that it is the solution correspond-
ing to the forcing term f1 (, , t). Clearly the result Equation (3.23) is independent of the
normalization of the Pl functions. However, using the standard normalization, we have
Z
1
(Pl (cos ))2 sin d = 1 , (3.24)
0 l+ 2
and Equation (3.23) becomes

X Z t
2 1 sin l (t )
u f1 = c (l + 2 ) Il () h1 ( ) d Pl (cos ) , (3.25)
l
l=0
where
Z 2 2
e /2
Il () = Pl (cos ) sin d . (3.26)
0 2
These integrals will need to be calculated numerically. Plots of these integrals are shown in
Figure 3.3A, where we see that Il () 1 as l 0 or 0. Figure 3.3B shows how we
select the upper limit lmax for the infinite summation in the displacement function u f .
The time integral in Equation (3.25) will give a complicated result involving error functions
of complex argument. However, the integral simplifies greatly if the upper limit is replaced
by . This will be a good approximation for all times t i.e., at times after the source
has ceased to act. Using the source-time function h 1 (t) (Eq. 3.17), we find (see Appendix D)
Z 2 2
!

e /2 sin l (t ) sin l t 2 2 /2
d = e l (t ) , (3.27)
2 l l
and Equation (3.25) becomes

X sin l t 2 2 /2
uf1 (, , t) = c2 (l + 21 ) Il () e l Pl (cos ) (t ) , (3.28)
l
l=0
where the Il () are given by Equation (3.26). Note that the width of the source-time function,
, and the width of the initial-shape function, , both influence the solution. The physical
meaning of Equation (3.28) is that it is the seismogram for a fixed surface point, or a snapshot
of the distorted spherical membrane for a fixed time.
The spherical average of the solution is calculated by integrating Equation (3.28) over the
entire sphere. However, only the l=0 term gives a non-zero contribution:
Z 2 Z
1
u
f1 (t) = uf1 (, , t) sin dd
4 0 0
Z 2 Z
1
= uf1 (, , t)|l=0 sin dd
4 0 0
Z 2 Z
1
= c2 12 I0 () t sin dd
4 0 0
1 2 1
= 2 c 4 I0 () t 2 2
1 2
= 2 c I0 () t , (3.29)
where in the third step we have used

2 1sin l t 2 2 /2
uf1 (, , t)|l=0 = lim c (l + e l
2 ) Il () Pl (cos )
l0 l
sin l t
= c2 12 I0 () lim
l0 l
2 1
= c 2 I0 () t .
Note that the spherical average increases linearly with time.

3.3 The solution for given forcing f = t
f1 (, , t)
Our forcing term thus far has been Equation (3.16). Let us now use a forcing term

f2 (, t) = f1 (, t) = h2 (t) g() , (3.30)
t
where
2 2
d t et /2
h2 (t) = h1 (t) = 2 . (3.31)
dt 2
The forcing term is again cylindrically symmetric, as we have not changed the initial-shape
function g(); the source-time function h 2 (t) is shown in Figure 3.1. In this case the expression
for displacement, Equation (3.25), is given by
Z 2 2
!
X t
e /2 sin l (t )
u f2 = c 2 (l + 21 ) Il () d Pl (cos ) , (3.32)
2 2 l
l=0
where we have simply substituted the forcing term h 2 ( ). In place of Equation (3.27) we have
(see Appendix D)
Z 2 2
!

e /2 sin l (t ) 2 2
d = cos l t el /2 (t ) , (3.33)
2 2 l
and in place of Equation (3.28) we have

X 2 2 /2
uf2 (, , t) = c2 (l + 21 ) Il () cos l t el Pl (cos ) (t ) . (3.34)
l=0
In this case, the spherical average (compare Eq. 3.29) is constant:
Z 2 Z
1
u
f2 (t) = uf2 (, , t) sin dd = 1
2 c2 I0 () . (3.35)
4 0 0
Figure 3.4 shows plots of Equation (3.34) at two particular times. As the wave pulse

spreads out over the sphere, the amplitude of the pulse varies approximately as K/ sin ,
where K is the maximum displacement of the pulse at the equator ( = /2).
3.4 The solution for given initial shape g(, )
We now look at the case with zero forcing, but with a prescribed initial shape. The solution
in the case that u(, , t) and u(, , t)/t are given at time t = 0, will be similar, except
that in place of Equation (3.10) we have
1 d2 ulm
+ l(l + 1) ulm = 0 , (3.36)
c2 dt2
subject to initial conditions that u lm and dulm /dt have prescribed values at time t = 0. Let
us assume that u(, , 0) = g(, , t) is given, and that the time derivative is zero, at t = 0.
(These are analogous to the initial conditions we set for waves on the circular membrane, e.g.,
Eqs. 2.16 and 2.17.) Then we require the solution, u lm (t) of Equation (3.36) satisfying
ulm |t=0 = glm (3.37)

dulm
= 0, (3.38)
dt t=0
where glm are the spherical harmonic coefficients of the given initial-shape function:
Z 2 Z
1
glm = Ylm (, ) g(, ) sin d d . (3.39)
Nl 0 0
In this case, the solution (to Eq. 3.36) is simply
ulm (t) = glm cos l t , (3.40)
which can be compared to Equation (3.12). Let us assume we have the initial-shape function
g(, , t) given by Equation (3.18). Again, the terms for m 6= 0 will all vanish, and we find:

X
ug (, , t) = (l + 21 ) Il () cos l t Pl (cos ) . (3.41)
l=0
We use the notation ug to emphasise that this is the solution corresponding to the given
initial shape g(, , t). Note that for t=0 the solution is simply g().
The spherical average is again constant (compare Eq. 3.35):
1
u
g (t) = 2 I0 () . (3.42)
Comparing Equations (3.29), (3.35), and (3.42), we find:
u
f2 = c2 u
g (3.43)
u
f1 f2 = t c 2 u
= tu g , (3.44)
where the fi subscript denotes the forcing term used for the displacement function u fi .
3.5 Discussion
Equations (3.28), (3.34), and (3.41) give the solutions to three problems: (i) with forcing f 1 ,
(ii) with forcing f2 = f1 /t, and (iii) with given initial shape g (and no forcing):

X
2 sin l t 2 2 /2
uf1 (, , t) = c (l + 12 ) Il () e l Pl (cos ) (t ) (3.45)
l
l=0
X
2 2 /2
uf2 (, , t) = c2 (l + 12 ) Il () cos l t el Pl (cos ) (t ) (3.46)
l=0

X
ug (, , t) = (l + 12 ) Il () cos l t Pl (cos ) , (3.47)
l=0
where Il () is given by Equation (3.26) and l is given by Equation (3.11).

In the case of the forcing term f1 , the time-derivative of the solution is given by

!
d X sin l t 2 2 /2
2 1
u f1 = c (l + 2 ) Il ()
e l Pl (cos )
dt l
l=0

X
2 1 d sin l t 2 2 /2
= c (l + 2 ) Il () e l Pl (cos )
dt l
l=0

X 2 2 /2
= c2 (l + 21 ) Il () cos l t el Pl (cos ) (t ) . (3.48)
l=0
Thus, examining Equations (3.46) and (3.48), we have the relationship
u f1 = u f2 (t ) (3.49)
u f |h = h1 (t) = uf |h = h 1 (t) (t ) , (3.50)
i.e., the time-derivative of the displacement due to a Gaussian source-time function h 1 (t)
is equal to the displacement due to a time-derivative Gaussian source-time function h 1 (t).
Equation (3.50) comes naturally by differentiating Equation (3.1) with respect to time:
1 2 u
2 u = f = g(, ) h 1 (t) . (3.51)
c2 t2
Solutions for an impulsive source ( = 0)
We now look at the case where = 0 (impulsive force). In this limit, the Gaussian source-time
function h1 (t) becomes a delta function. Equation (3.45) and Equation (3.48) become

X
2 sin l t
u f1 0 = c (l + 12 ) Il () Pl (cos ) (t 0) (3.52)
l
l=0
X
u f1 0 = c2 (l + 12 ) Il () cos l t Pl (cos ) (t 0) , (3.53)
l=0
where uf1 0 = uf1 0 (, , t) is the displacement due to an impulsive force with source-time
function h1 (t). The same process follows for the solutions due to the source-time function
h2 (t) (Eq. 3.31):

X 2 2 /2
u f2 = c2 (l + 12 ) Il () (l sin l t) el Pl (cos ) (t ) (3.54)
l=0

X
u f2 0 = c 2 (l + 12 ) Il () cos l t Pl (cos ) (t 0) (3.55)
l=0

X
u f2 0 = c2 (l + 12 ) Il () (l sin l t) Pl (cos ) (t 0) . (3.56)
l=0
Examining Equations (3.55), (3.53), and (3.47), we have the following relationship:
uf2 0 = u f1 0 = c2 ug (t 0) . (3.57)
The effect of nonvanishing is to introduce a frequency decay, arising from the convolution
of the instantaneous solution uf1 0 with the source-time function h1 (t) (Eq. 3.17). Thus, for
general , we have
u f1 = uf1 0 h1 (t) (D.33) (3.58)

d
u f1 = (uf1 0 h1 (t))
dt
= u f1 0 h1 (t) (D.31)
= c2 h1 (t) ug (3.57), (D.27)

Z t
= c2 h1 ( ) ug (, , t ) d (D.26)

Z
!
t X
2 1
= c h1 ( ) (l + 2 ) Il () cos l (t ) Pl (cos ) d (3.47)
l=0

X Z t
= c2 (l + 21 ) Il () h1 ( ) cos l (t ) d Pl (cos ) , (3.59)
l=0
where * denotes convolution (Eq. D.26), and it is understood that u g is taken to vanish for
negative times. Equation (3.59) is in agreement with the result obtainable by taking the time
derivative in Equation (3.23):
Z !
d X t
sin l (t )
u f1 = c2 (l + 12 ) Il () h1 ( ) d Pl (cos )
dt l
l=0
X Z t
2 1 d sin l (t )
= c (l + 2 ) Il () h1 ( ) d Pl (cos )
dt l
l=0
X Z t
2 1 d sin l (t )
= c (l + 2 ) Il () h1 ( ) d Pl (cos )
dt l
l=0
X Z t
= c2 (l + 21 ) Il () h1 ( ) cos l (t ) d Pl (cos ) . (3.60)
l=0
The Fourier transform of h1 (t) is (see Eq. D.25)
2 2
e
h(l ) = el /2 . (3.61)
Note that, after making the approximation t (replace the upper limit of t with in
Eq. 3.59 or Eq. 3.60), the effect of the convolution is to introduce a factor e
h(l ) in each term
of uf1 0 (Eq. 3.52) to get uf1 (Eq. 3.45). The same is true for the time-derivative, introducing
the factor e
h(l ) in each term of u f1 0 (Eq. 3.53) to get u f1 (Eq. 3.48).
Solutions for a point source ( = 0)
We now look at the case where = 0 (point source). In this limit, the Gaussian initial-shape
function g() becomes a delta function. The integrals I l () (Eq. 3.26) simplify (see Ap-
pendix D):
Z 2 2
e /2
Il (0) = lim Pl (cos ) sin d = Pl (1) = 1 . (3.62)
0 0 2
The spatial point-source solutions ( = 0) are thus (see 3.453.47):

X sin l t 2 2 /2
uf1 = c 2
(l + 21 ) e l Pl (cos ) (t ) (3.63)
l
l=0
X
2 2 /2
uf2 = c2 (l + 12 ) cos l t el Pl (cos ) (t ) (3.64)
l=0

X
ug = (l + 12 ) cos l t Pl (cos ) , (3.65)
l=0
where the * is used to distinguish the spatial point-source form of the solutions.
Lastly, we can present solutions for a point source in both time ( = 0) and space ( = 0):

X sin l t
uf1 0 = c2 (l + 21 ) Pl (cos ) (t 0) (3.66)
l
l=0

X
uf2 0 = c2 (l + 12 ) cos l t Pl (cos ) (t 0) , (3.67)
l=0
The solutions in Equations (3.65), (3.66), and (3.67), however, contain arbitrarily high fre-
quencies and wavenumbers and will probably not be useful unless a frequency filter is applied.
-3
x 10
1.8
1.6
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0
-800 -400 0 400 800
Time (s)
-6
x 10
5.0
4.0
3.0
2.0
1.0
0
-1.0
-2.0
-3.0
-4.0
-5.0
-800 -400 0 400 800
Time (s)
Figure 3.1: Initial-shape function g(, ), and two source-time functions, h 1 (t) and
h2 (t) = h1 /t (Eqs. 3.18, 3.17, and 3.31). The two forcing terms discussed in this paper
are f1 (, , t) = h1 (t) g(, ) and f2 (, , t) = h2 (t) g(, ). The initial-shape function is plot-
ted on the order q=4 triangular grid. For the examples here, the parameters governing the
width of the Gaussians are = 0.0713 and = 204.5. For the analytical solutions, of course,
the gridpoint separation d and the timepoint separation dt are not relevant.
220
200
180
160
140
120
100
Equator
Equator
80
60
40
20
0
100 80 60 40 20 0 20 40 60 80 100
Figure 3.2: The relationship between and wid for the Gaussian initial-shape function g().
In this example, g() is plotted with wid = 20 at a height of 0.05A (r = 0.05), which gives
(wid , r) = 0.0713 (Eq. 3.21).
A B
1.0
150
0.9
0.8 125
0.7
100
0.6
0.5 75
0.4
50
0.3
0.2 25
0.1
0 0
0 0 5 10 15 20 25 30 35 40
0 0.05
5
10 n,
15
20
25 0.10 au ssia
Harmonic De 30 fG
gree, L th o
Wid
Figure 3.3: Surface plot of Il () (Eq. 3.26). A. We see graphically that I l () 1 as l 0

or 0. The latter case represents the point source approximation (see Eq. 3.18). B. The
narrower the Gaussian is in the initial-shape function, the greater the number of l terms
required such that Ilmax () < tol, where tol = 103 in this plot. and wid can be thought
of interchangeably (Figure 3.2).
Frame 22
(time = 0.0 s)
800
600 736.6
Displacement
400
200
-200
,
-400
-600
-800
-90 -60 -30 0 30 60 90
Latitude
800
600
141.4
Displacement
400
200
-200
Frame 38 ,
-400
(time = 801.1 s)
-600
-800
-90 -60 -30 0 30 60 90
Latitude
800
Max Displacement
600
400
K = 109.1
,
200
0
-90 -60 -30 0 30 60 90
Latitude
Figure 3.4: Analytical solutions for wave propagation on the sphere, using the forcing term
f2 (, , t) of Figure 3.1. The wave is shown both as a 3D figure and as a 2D plot for two
different times t. The waveform in the 2D plots is Equation (3.34) for fixed t and varying
. The bottom plot shows how the maximum
values of the wave pulse varies as a function
of latitude. The dashed curve is K/ sin , where K is the maximum displacement of the
pulse at the equator and is the colatitude. As expected, the displacement is maximum at
the poles and minimum on the Equator. See also Figure 5.7.
Chapter 4
Spherical hexagonal grids
Contents
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.2 Constructing the hexagonal grids . . . . . . . . . . . . . . . . . . . . . . . 37
4.3 Properties of the hexagonal grids . . . . . . . . . . . . . . . . . . . . . . . 39
4.4 Calculating the Laplacian on the hexagonal grids . . . . . . . . . . . . . 42
4.5 Error in the numerical approximation of the Laplacian . . . . . . . . . . 44
4.6 Order of accuracy for the numerical approximation . . . . . . . . . . . . 46
4.7 Variations and optimizations on the hexagonal grids . . . . . . . . . . . 47
4.8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.1 Introduction
Every numerical computation is carried out on a discrete number of gridpoints. The choice of
the distribution of gridpoints will depend on the numerical method that is being used. The
two basic types of numerical methods are global methods and local methods. Global meth-
ods refer to spectral methods, which use all the gridpoint values at each timestep iteration;
local methods, such as finite difference methods or finite element methods, use only nearby
gridpoint values at each iteration. Historically, these methods were developed roughly in the
following periods (Trefethen, 2000): 1950s, finite difference methods; 1960s, finite element
methods; and 1970s, spectral methods.
34
CHAPTER 4. Spherical hexagonal grids 35
The use of spherical geodesic grids 1 in approaching wave propagation on a sphere was
developed in the field of climate modelling in the late 1960s. The aim of these grids was to
cover the sphere with grid cells as uniformly as possible, and in a manner that allowed for
the implementation of a finite difference method to solve equations of fluid flow. Developing
grids inspired by Buckminster Fullers geodesic dome, Williamson (1968) and Sadourny et al.
(1968) discretized the sphere into spherical triangles that were nearly uniform in area. Part of
their motivation was to develop an alternative grid to the standard latitudelongitude grids,
which have an unnecessarily high resolution of grid cells near the poles (Figure 4.1).
Although spectral methods are predominant in meteorology 2 , they are far from ideal.
Stuhne and Peltier (1996) point out three fundamental problems of the spectral approach to
fluid dynamics: (1) it over-resolves polar regions, (2) its results may be less physical in that
it has the potential to manifest anomalous long range interactions, and (3) it has poor
theoretical efficiency at high spatial resolution (p. 58). These shortcomings have motivated
a resurgence of the local methods over the past ten years, a debate which is reviewed in most
of the introductions of the post-1990 studies in Table 4.1 (e.g., Swarztrauber et al., 1997;
Ringler et al., 2000; Randall et al., 2002).
Developments in local approaches have included variations on the grids of Williamson and
Sadourny, notably those by Masuda and Ohnishi (1987) and Heikes and Randall (1995a),
among others (Table 4.1). In these studies, once the sphere is discretized, the grid is tested
using some version of equations describing fluid flow, such as the shallow-water equations.
Many authors emphasize the need to develop and test numerical models (Williamson et al.,
1992).
Instead of modelling fluid waves in the atmosphere, we shall be modelling Rayleigh wave
propagation on the surface of the earth. Our spherical geodesic grid and finite difference
calculation is based on those used in the climate studies above (Heikes and Randall, 1995a).
Various names for the grids have been used (Table 4.1); our nomenclature will be triangular
grid for a spherical grid of triangles, and hexagonal grid for a spherical grid of predomi-
nantly hexagons.
1
This nomenclature refers to the fact that the grids are characterized by arcs between gridpoints on the
unit sphere, i.e., geodesics (Williamson, 1968).
2
The method which is cited as most successful at modelling fluid flow is in fact the pseudo-spectral method
(Stuhne and Peltier , 1996).
Figure 4.1: Two discretizations of the sphere. One problem of the latitudelongitude dis-
cretization (left) is that the areas of the grid cells go to zero as the cells approach the poles.
Alternative discretizations (right) try to keep the area of the cells as uniform as possible.
Table 4.1: Spherical grids used in other studies, each which used an icosahedron as the
initial triangular grid to subdivide. The third column indicates the type of numerical method
used: FDM = finite difference method, FEM = finite element method. The fourth column
indicates whether the study employed a discrete Laplacian operator on the grid, which is of
interest to our study.
Study Name of Grid Method 2
Williamson (1968, 1970) spherical geodesic grid FDM yes

Williamson (1969) triangular grid FDM yes
Sadourny et al. (1968) icosahedralhexagonal grid FDM yes
Sadourny and Morel (1969) hexagonal grid FDM yes
Cullen (1974) icosahedral grid FEM yes
Baumgardner and Frederickson (1985) icosahedral grid FEM yes
Sword et al. (1986a,b) icosahedral grid FEM yes
Masuda and Ohnishi (1987) icosahedralhexagonal grid FDM yes
Wang and Dahlen (1995a) spherical tessellation grid N/A no
Wang et al. (1998)
Heikes and Randall (1995a,b) twisted icosahedral grid FDM yes
Stuhne and Peltier (1996, 1999) icosahedral grid FEM yes
Swarztrauber et al. (1997) icosahedral geodesic grid FDM yes
Thuburn (1997) hexagonalicosahedral grid FDM yes
Ringler et al. (2000) spherical geodesic grid FDM yes
Randall et al. (2002) FDM yes
Ringler and Randall (2002a,b)
Table 4.2: The five platonic solids, plus the buckeyball for comparison. The internal angle
for each solid is given by , with the number of faces meeting
P at each vertex in parenthesis.
The sum of the interior angles at each vertex is given by .
P
platonic solid polygonal face faces, f vertices, v edges, e
tetrahedron triangle 4 4 6 60 (3) 180

cube (hexahedron) square 6 8 12 90 (3) 270
octahedron triangle 8 6 12 60 (4) 240
dodecahedron pentagon 12 20 30 108 (3) 324
icosahedron triangle 20 12 30 60 (5) 300
buckeyball pentagon (12) 32 60 90 108 (1) 348
(q=0 hexagonal grid) hexagon (20) 120 (2)
4.2 Constructing the hexagonal grids
Each point on a grid (triangular or hexagonal) is always found on the unit sphere, which
represents the surface of the earth. The hexagonal grids are generated from the triangular
grids. Below we explain how to generate triangular grids and hexagonal grids.
Triangular grids
The algorithm for generating successively finer triangular grids is simple: take each existing
triangular patch, find the midpoints of the three sides, and then construct four new triangular
patches. This process, shown in planar view in Figure 4.4, is carried out for triangular
gridpoints on the unit sphere.
What type of spherical polyhedron do we choose for our lowest order triangular grid? It
seems reasonable to construct our triangular grid from a polyhedron whose faces are uniform
in area and whose vertices lie on the unit sphere: this limits us to five polyhedra, the platonic
solids (Table 4.2). Properties of the platonic solids are: (i) the same number of faces meet at
each vertex, (ii) each face is a regular congruent polygon, (iii) the sum of the interior angles
at each vertex is <360 , and (iv) each solid satisfies Eulers Polyhedron Theorem,
f +ve=2, (4.1)
where f is the number of faces, v is the number of vertices, and e is the number of edges
(see Table 4.2). To construct our order q=0 triangular grid, we begin with a dodecahedron,
subdividing each face into five 72 /54 /54 triangles and then subdividing each triangle into
four new triangles (Figure 4.5). Our order q=0 grid contains the vertices of both the dodeca-
hedron and icosahedron3 (Figure 4.6) and is not a platonic solid. It is clear from the criteria
for the platonic solids that the maximum number of equal-area triangles that can form a
spherical triangular grid is 20 (Table 4.2); our goal in constructing a finer mesh of points on
the sphere is to keep the area of the triangles as uniform and small as possible.
The mechanism for constructing higher order triangular grids (Figure 4.4) shows that, be-
ginning with our subdivided dodecahedron, the number of faces, f t , of the qth-order triangular
grid is
ft (q) = 60 4q , (4.2)
which can be seen in Table 4.3.
Hexagonal grids
The hexagonal grids are constructed from the triangular grids, a process shown schematically
in Figure 4.7. The corner of each hexagonal face is found by calculating the central point of
the three vertices forming the face of the corresponding triangular grid. The central point of
the vectors T1 , T2 , and T3 is given by
(T3 T1 ) (T2 T1 ) T1 T 3 + T 3 T 2 + T 2 T 1
V = = , (4.3)
|(T3 T1 ) (T2 T1 )| |T1 T3 + T3 T2 + T2 T1 |
where T1 , T2 , T3 , and V are all on the unit sphere (Augenbaum and Peskin, 1985). Note
that T1 , T2 , and T3 are ordered in a clockwise sense, looking from the exterior, so that V 2
lies in the correct hemisphere (see Figure 4.7B).
The order q=0 hexagonal grid is in fact the well-known buckeyball (C 60 molecular struc-
ture), with 12 regular pentagonal faces and 20 regular hexagonal faces (having two types of
faces, it is not a platonic solid Table 4.1). The buckeyball has the property that the dis-
3
It is quite simple to construct and visualize the icosahedron. It can be constructed from three mutually
perpendicular rectangles whose sides have the golden ratio t = (1 + 5)/2 1.62. The vertices are the corners
of the rectangles and are thus given by the points (0, t, 1), (t, 0, 1), (t, 1, 0). (Note that these points
are on the unit sphere.)
tance between one point and its neighbors is the same for all points. Higher-order hexagonal
grids are shown in Figure 4.8A; the relationship between successive orders of hexagonal grids
is shown in Figure 4.8B. The number of faces f h of the qth-order hexagonal grid is
fh (q) = 30 4q + 2 , (4.4)
which can be seen in Table 4.3.

Here we note two important aspects of the grids. First, note that the vertices of the
triangular grids are the face centers on the hexagonal grids (Figure 4.7). Similarly, the vertices
of the hexagonal grids are the face centers of the triangular grids. This gives us the equalities
fh = v t (4.5)
vh = f t . (4.6)
Function values are calculated at the center of each hexagonal face, or, in other words, at each
triangular vertex. Second, note that the triangular vertices accumulate as the grids become
finer; in other words, the vertices of the coarser triangular grids represent a subset of the
vertices of the finer grids (Figure 4.5). This is not the case with the hexagonal grids, whereby
the gridpoints of the coarser grids do not coincide exactly with gridpoints of the finer grids
(Figure 4.8).
So far we have qualitatively presented how to construct the spherical grids, but we have
not presented an explicit algorithm for indexing the vertices, which is crucial. Our code
was originally written by Woodhouse in 1995 and modified over the course of this study to
accommodate the numerical method (Section 4.4).
4.3 Properties of the hexagonal grids
Table 4.4 shows some properties of the hexagonal grids. The values are for a sphere with
earths radius:
earth radius a = 6371 km
circumference 2a = 40, 030 km
surface area 4a2 = 5.1 108 km2 .

Table 4.3: Geometric properties of the triangular grids (Figure 4.5) and the hexagonal grids
(Figure 4.8A). The order q=0 triangular grid is a subdivided dodecahedron. Each hexagonal
grid and triangular grid satisfies Equation (4.1).
triangular grid hexagonal grid

order, q faces, ft vertices, vt edges, et faces, fh vertices, vh edges, eh
0 60 32 90 32 60 90
1 240 122 360 122 240 360
2 960 482 1 440 482 960 1 440
3 3 840 1 922 5 760 1 922 3 840 5 760
4 15 360 7 682 23 040 7 682 15 360 23 040
5 61 440 30 722 92 160 30 722 61 440 92 160
6 245 760 122 882 368 640 122 882 245 760 368 640
q ft = 60 4q vt = 30 4q + 2 et = 90 4q fh = v t vh = f t eh = e t
Table 4.4: Properties of the hexagonal grids used in this study. f eq is the approximate number
of faces on the equator. d is the average angular distance between face centers. Areas and
lengths assume a spherical earth with radius a = 6371 km. Compare with Table 4.5.
area of distance between

faces (km2 ) face centers (km)
order, q faces, fh feq aave amin /amax dave (d) dmin /dmax
0 32 10 15 939 515 0.9413 4320 (38.86 ) 0.8940
1 122 19 4 180 856 0.9142 2208 (19.86 ) 0.8608
2 482 39 1 058 225 0.9070 1111 (9.99 ) 0.8520
3 1 922 78 265 382 0.8776
557 (5.00 ) 0.8498
4 7 682 155 66 397 0.8700 278 (2.50 ) 0.8492
5 61 440 310 16 603 0.8681
139 (1.25 ) 0.8491
6 245 760 620 4151 0.8676 69 (0.63 ) 0.8490
Table 4.5: Properties of the hexagonal grids (twisted icosahedron polyhedron) of Heikes
and Randall (1995a). feq is the exact number of faces on the equator. d is the average
angular distance between face centers. Areas and lengths assume a spherical earth with
radius a = 6371 km. Compare with Table 4.4.
area of distance between

faces (km2 ) face centers (km)
order, q faces, fh feq aave amin /amax dave (d) dmin /dmax
0 42 10 12 140 000 0.885 3755.5 (33.8 ) 0.881
1 162 20 3 149 000 0.774
1916.4 (17.2 ) 0.848
2 642 40 794 600 0.763 962.4 (8.7 ) 0.840
3 2 562 80 199 100 0.742
481.7 (4.3 ) 0.838
4 10 242 160 49 800 0.736 240.9 (2.2 ) 0.837
5 40 962 320 12 450 0.733
120.5 (1.1 ) 0.837
q 40 4q + 2 5 2q+1
With the hexagonal grids, each corner of a hexagon (or pentagon) is a junction of three edges.
And since each edge has exactly two (hexagonal) vertices, we have 2e = 3v (Augenbaum and
Peskin, 1985). Using Equation (4.1), which is true for an arbitrary arrangement of points on
the surface of a sphere, we see that v = 2f 4 and e = 3f 6 (Augenbaum and Peskin, 1985).
These relationships are evident in Table 4.3 for the hexagonal grids.
Area calculations for each hexagonal face are done by subdividing each face into six (or
five) triangles, and then summing the areas of the spherical triangles. The area of a spherical
triangle with sides having arc lengths a, b, and c is
s !
s sa sb sc
1
ai = 4 tan tan tan tan tan , (4.7)
2 2 2 2
where s = (a + b + c)/2. Thus, subdividing each face into spherical triangles, the area of
hexagonal face with center Ti is
N
X
Ai = ai , (4.8)
i=1
where N = 6 or N = 5.
Figures 4.9 and 4.10 show the distribution of three geometric properties of the hexagonal
grids: area of faces, distance to each neighbor, and length of hexagonal face edge. The
variation in these values is larger than one might have expected, based on our subdivision
algorithm. We see that the distribution of areas is not particularly uniform the deviation
from the mean is about 10%.
Comparison of hexagonal grids with those used in other studies
Several climate modelling studies use a subdivided icosahedron as the order q=0 triangular
grid (Table 4.1). A comparison of our order q=0 triangular grid the subdivided dodeca-
hedron with the subdivided icosahedron shows the following differences: the subdivided
dodecahedron has 60 triangular faces, each triangle having approximately 67.6 /56.2 /56.2
internal angles (Figure 4.5); the subdivided icosahedron has 80 triangular faces, 60 having
approximately 68.8 /55.6 /55.6 internal angles and 20 being equilateral triangles.
For the subdivided icosahedron, the number of faces on the triangular grid is f t (q) = 80 4q
(compare Eq. 4.2), and the number of faces on the hexagonal grid is f h (q) = 40 4q + 2 (com-
pare Eq. 4.4). Table 4.5 shows some properties of the initial-icosahedron hexagonal grids,
which can be compared with the properties of the initial-dodecahedron hexagonal grids (ours)
in Table 4.4. Based on the values in the two tables, our hexagonal grids appear to be better
at having the greatest uniformity of face area (a min /amax 1) and distance between faces
(dmin /dmax 1). Optimizations of the initial-icosahedron grids are discussed in Section 4.7.
4.4 Calculating the Laplacian on the hexagonal grids
The hexagonal grids are discretized in a way that makes it possible to calculate the Laplacian
2 u of a given function u(, ) on the surface of the sphere by using a numerical approxima-
tion. We employed the approximation used in the climate studies to calculate the Laplacian
on the sphere4 (Masuda and Ohnishi, 1987; Heikes and Randall, 1995a).
4
Generally the Laplacian 2 is one of several operators approximated on the grids, such as the Jacobian
J(, ) and the flux-divergence operator (), where and are arbitrary scalar fields (e.g., Masuda
and Ohnishi, 1987).
From Gausss theorem, we have

dn
Z I C
u
2 u dA = ds , (4.9) ds
A C n
A
where u is a scalar function and u/n is the gradient
of u, normal to a small segment ds of a curve C (Fig-
ure 4.2). Now, consider a particular face center T 0 (, )
on the hexagonal grid. The function u(T 0 ) = u0 is de- Figure 4.2: Sketch of integra-
tion parameters for Eq. 4.9. The
fined at the center of the face, which has a spherical area scalar function u is defined every-
A0 defined by Equation (4.8). where in space.
If the faces on the grid are small, then 2 u can be considered constant across each face,
and we can approximate the left-hand side of Equation (4.9) by
Z Z

u dA uT0
2 2
dA = 2 uT0 A0 . (4.10)
A A
The line integral on the right-hand side of Equation (4.9) can be approximated by summing
the u-gradients between T0 and its N = 6 (or N = 5) neighbors:
I N
X
u li
ds (ui u0 ) , (4.11)
C n Li
i=1
where Li is the arc distance to the adjacent face-center, and l i is the arc length of the
corresponding hexagonal face edge (Figure 4.11). Using these results, we find
N
1 X li
D2 (u)|T0 2 uT0 (ui u0 ) (4.12)
A0 Li
i=1
N N
!
X li X li
= ui u0 , (4.13)
Li A0 Li A0
i=1 i=1
where D2 is our notation for our discrete Laplacian. Equation (4.13) emphasizes the depen-
dence of the Laplacian on the function values u i at neighboring points. Figure 4.11 shows the
variables used in calculating the Laplacian at each hexagonal face.
4.5 Error in the numerical approximation of the Laplacian
The error in the numerical Laplacian approximation is determined by comparing the numerical
value of 2 u(, ) with the actual value, which can be the (exact) analytical value, if we
choose an appropriate test function u(, ). Spherical harmonic functions Y lm (, ) are a
convenient choice for u(, ), since they are well-defined on the sphere (Appendix A) and are
eigenfunctions of the Laplacian operator:
2 Ylm = l(l + 1)Ylm . (4.14)
Thus, with u Ylm , from Equation (4.12) we obtain [2 u]nu and from Equation (4.14) we
obtain [2 u]an . Both of these functions are defined at each gridpoint of the order q grid.
We define the error in the numerical calculation of the Laplacian at T i (, ) as

[2 u]nu [2 u]an
Ti Ti
D2error(, ) = , (4.15)
max {[2 u]an }
where max finds the maximum single value of 2 u on the sphere. Thus, for each point on
the sphere, the absolute difference between the model approximation and the actual value is
normalized to the maximum actual value on the sphere 5 .
Figure 4.12A shows the distribution of the error over the q=6 grid, and we see that the
error is not at all uniformly distributed. Although the average error at each point is small
(D2error(, ) 104 ), there is a small fraction of points having errors 1020 times greater
than the average.
Figure 4.12B-C takes a closer look to show how isolated the bad error gridpoints are.
The error at gridpoint 112 (D2error(112) = 6.807 10 3 ) is on average 24 times greater than
the error at one of its neighboring points. Furthermore, the error at 112 does not decrease
with finer grids, which means that the numerical approximation for the Laplacian operator is
zero-order accurate at this point. What makes these particular gridpoints so bad? The answer
seems to lie in the geometric design of the grids. We see in Figure 4.9 that gridpoint 112 is
a particularly anisotropic point, where the areas are relatively large in one direction and
relatively small in the perpendicular direction. Our discrete approximation of the Laplacian
5
This normalization gives the smallest value for the error. If we normalized to the analytical value at each
point, however, we would have points with near-infinite error, where Ylm (, ) = 0.
(Eq. 4.12) is not effective at these gridpoints.
Sphericallyaveraged error
We are also interested in calculating the spherical average of the error. One way of calculating
the sphericallyaveraged error is take the mean of all the error values, since the gridpoints
are nearly uniformly distributed:
D2error(l, m, q) = mean { D2error(, ) } . (4.16)
We also adopt the norm conventions of Heikes and Randall (1995b):
1 X
||fi || = Ai |fi | (4.17)
Ae
cells
" #1/2
1 X
||fi ||2 = Ai (fi )2 (4.18)
Ae
cells
||fi || = max {|fi |} , (4.19)
where fi is a function defined at all the gridpoints (i.e., cells) on the sphere and A e = 4a2
is the area of the sphere6 . We use these norms to establish three more measurements of the
sphericallyaveraged error7 :
2
[ u]nu [2 u]an
D2one(l, m, q) = (4.20)
|| [2 u]an ||
2
[ u]nu [2 u]an
2
D2two(l, m, q) = (4.21)
|| [2 u]an ||2
2
[ u]nu [2 u]an

D2inf(l, m, q) = . (4.22)
|| [2 u]an ||
Thus, given a test function u Ylm on a given grid q, we have four equations for the
sphericallyaveraged error in terms of l, m, and q. Figure 4.13 shows how three of the error
values vary with grid order q for a particular Y lm function. We immediately see that the error
converges for the one-norm and two-norm, but not for the infinity norm, indicating that there
must be at least one bad gridpoint (note Eq. 4.19).
6
To be precise, Ae is the sum of the areas of all the faces, which is equal to 4a2 to within 105 .
7
These are essentially equivalent to the three standard normalized global errors put forth by Williamson
et al. (1992).
In general, we see that the error in the numerical Laplacian increases as degree l increases,
which is what one would expect, since the higher the degree, the larger the surface gradients.
Also, we observe relative lows in error values for test functions with order m = 5 N (N =
0, 1, . . . , 8); this can be attributed to the symmetry of the test function matching the symmetry
of the grids (Figure 4.14). (The vertical color streaks for these m-values in Figure 4.15 indicate
this feature as well.) Our findings suggest that the orientation of the grids has an effect on
the errors, something that was first noted by Williamson (1968).
The error in the numerical Laplacian was also compared between (1) the grid calculations
using arc distances and spherical areas and (2) the grid calculations using chord distances
and planar areas (see Figure 4.11). There is a negligible difference in the results for the finer
grids (q 4). We used arc distances and spherical areas, since there is no increase in the
calculation time, as the geometric properties of the grids are calculated once prior to iterating
the wave equation.
4.6 Order of accuracy for the numerical approximation
As Heikes and Randall (1995b) point out, it is difficult to analytically investigate the order of
accuracy of a discrete operator on an irregular grid. Suppose our finite difference form of the
Laplacian operator is mth order for the order-q spherical grid having average grid spacing h:
[2 u]an [2 u]nu = O [ hm ] (4.23)
Suppose we have grids having two different spacings, h 1 and h2 , where h1 = 2h2 (i.e., grid 2
is finer than grid 1). The we have:
2
[ u]nu (h1 ) [2 u]an O [ (h1 )m ] O [ 2m (h2 )m ]
= = = 2m , (4.24)
|| [2 u]nu (h2 ) [2 u]an || O [ (h2 )m ] O [ (h2 )m ]
where || || is some norm (Heikes and Randall, 1995b), and m is the order of accuracy of the
Laplacian operator. [2 u]nu (h1 ) is a vector of values of the numerical Laplacian evaluated at
each gridpoint on the grid having average gridpoint spacing h 1 .
Table 4.6: Studies of wave propagation on the global scale that use subdivided cubic grids.
FDM = finite difference method; SEM = (pseudo-)spectral element method.
Study Name of Grid Numerical Method
Sadourny (1972) quasi-uniform spherical grid FDM

Rancic et al. (1996) expanded spherical cube FDM
Ronchi et al. (1996) cubed sphere FDM
Komatitsch and Tromp (1999, 2002a,b) cubed sphere SEM
Capdeville et al. (2003) cubed sphere SEM
Geometrically, the order of accuracy is the slope of the error-vs-gridpoint-spacing curve

(Figure 4.13B). We can plot the order of accuracy between successive grids q and q+1, and
plot these values for each of the three norms (Figure 4.15). We see that our discrete approx-
imation of the Laplacian operator for the finest grids (q 4) is nearly second-order accurate
for the one-norm and two-norms, but zero- or first-order accurate for the infinity norm (Fig-
ure 4.15). This can be seen in other figures as well. The leveling of the infinity-norm curve in
Figure 4.13B indicates zero-order accuracy for at least one gridpoint, since the infinity norm
measures the greatest absolute deviation of the numerical value from the analytical value
(Eq. 4.22). In Figure 4.12 we can actually identify where these particularly bad points are
they happen to be at the gridpoints of the lower-order q grids.
In short, our numerical approximation for the Laplacian is nearly second-order, except at a
few bad points on the grids. This is consistent with the findings of Heikes and Randall (1995b)
(their Figure 6), who used the same numerical approximation but used slightly different grids
and different test functions.
4.7 Variations and optimizations on the hexagonal grids
It is apparent that discretizing the sphere as uniformly as possible is a challenging problem

(e.g., Saff and Kuijlaars, 1997; Cui and Freeden, 1997; Augenbaum and Peskin, 1985). Fur-
thermore, the task of communicating among the points, as required by numerical methods,
is an additional barrier. A wealth of studies in global climate modelling have used grids and
numerical methods similar to ours (Table 4.1). However, these grids have their weaknesses, as
we have shown that there are specific points in these grids that do not converge, supporting
the findings of Heikes and Randall (1995b) and Sword et al. (1986a).
Heikes and Randall (1995b) developed an improved grid by optimizing the distribution of
gridpoints such that the value (ri )4 was minimized for a particular order q grid (Figure 4.3).
This is a rather brute-force approach that requires a global minimization to a problem whose
numerical method is local. If the finite difference approximation is bad at certain points,
why not improve the method? Nevertheless, Heikes and Randall (1995b) show that with the
optimized grid, the errors converge to zero at all gridpoints on the sphere.
li r
i
A B
Li
Figure 4.3: Schematic diagram illustrating the parameter r i , defined as the distance between
the midpoint of a face edge (A) and the midpoint of the corresponding distance-between-
face-centers (B). This is the optimization parameter used by Heikes and Randall (1995b) for
optimizing the gridpoint spacing.
There are negative consequences associated with optimizing a particular parameter of

the grid. Heikes and Randall (1995b) perturbed the gridpoints to obtain a more uniform
area coverage or, equivalently, a more uniform density of gridpoints but this was at
the expense of increasing the range of distance-between-cells on the grid (Table 4.7). Sword
et al. (1986a) found similar results when optimizing the area coverage 8 , reporting that the
regularization increased the measure of local anisotropy in some regions, while smoothing
out sharp contrasts (p. 72). They defined local anisotropy as the difference between the
maximum and minimum distance-to-neighboring-cell for each triangular vertex.
8
Sword et al. (1986a) do not explain specifically how they optimized the areas: A least-squares problem
was set up, with the goal of making the areas of the individual triangles as equal to each other as possible
(p. 67).
Table 4.7: Optimizing the hexagonal grids: results from Heikes & Randall. Ratios from
the unoptimized grids are from Heikes and Randall (1995a) and shown earlier in Table 4.5.
Ratios from the optimized grids are from Heikes and Randall (1995b) and Randall et al.
(2002). After optimization in 1995, the area ratios were better (i.e., closer to 1) but the
distance ratios were worse.
Ratio of smallest cell Ratio of smallest distance

to largest cell to largest distance
grid hex between cell centers
order faces unopt- optimized optimized unopt- optimized optimized
q fh imized 1995 2002 imized 1995 2002
0 42 0.885 0.885 0.885 0.881 0.881 0.881
1 162 0.774 0.868 0.916 0.848 0.811 0.820
2 642 0.763 0.880 0.942 0.840 0.787 0.881
3 2562 0.742 0.877 0.948 0.838 0.778 0.799
4 10492 0.736 0.870 0.951 0.837 0.775 0.790
5 40962 0.733 0.867 0.952 0.837 0.776 0.788
Other studies have tried different discretizations of the sphere, such as a subdivided cube
(Table 4.6), which has seen developments in the field of seismology (Komatitsch and Tromp,
2002b). Cui and Freeden (1997) is a comprehensive study of the equidistribution of points
on the sphere, and they discuss the convergence of the subdivided tetrahedron, octahedron,
and icosahedron, in comparison with other distributions of gridpoints. They state: To our
surprise the results are not very convincing for larger pointsets (p. 606). Nevertheless, it
seems that the distribution of gridpoints in these subdivided platonic solids makes it easier
to implement the numerical approximations, which is, after all, is the purpose for using these
grids.
4.8 Summary
The important points from this chapter are included here.
1. Numerical methods on the sphere are required to analyze global wave propagation,
particularly in the fields of climate modelling, oceanography, and seismology. Present
approaches are predominantly spectral (global methods), but finite difference and
finite element methods (local methods) pose as appealing alternatives and have seen
significant developments (Table 4.1).
2. Novel distributions of gridpoint systems on the sphere such as the subdivided icosahe-
dron (Table 4.1) or dodecahedron (ours) have been developed in order to implement
local methods of numerical approximation.
3. Triangular grids are constructed by subdividing each triangular face into four faces
and then projecting the points to the unit sphere (Figure 4.5). Hexagonal grids are
constructed from triangular grids (Figure 4.7).
4. The hexagonal and triangular grids should be considered interchangeable. The vertices
of the triangular grids are the centers of the faces on the hexagonal grids. The function
values are calculated at each hexagonal face center, or, equivalently, at each triangular
vertex (Figure 4.7).
5. Our numerical approximation for the Laplacian operator ( 2 ) the discrete Laplacian
is a finite difference method and depends on the function values at the six (or five)
neighboring points, as well as on the geometric properties, such as the area of the cell
and distance-to-neighbors (Figure 4.11).
6. The discrete Laplacian is nearly second-order accurate at almost all of the gridpoints.
However, at a few gridpoints, the discrete approximation is zero-order (Figure 4.12).
These bad points may act as scatterers when we propagate waves on these grids.
7. It is possible to improve the discrete Laplacian by optimizing the gridpoint distribution,

but we did not attempt this.
bold = order q triangles (1)

dash = order q+1 triangles (4)
solid = order q+2 triangles (16)
Figure 4.4: Relationship between successive orders of triangular grids. Each triangle is
divided into four smaller triangles. This process is applied to spherical triangular patches on
the unit sphere (see Figure 4.5).
platonic dodecahedron order q=0 triangular grid order q=1 triangular grid
Figure 4.5: Construction of the q=0 triangular grid and higher order triangular grids.
Beginning with the order q=0 grid, the midpoints between vertices are calculated and then
projected to the surface of the unit sphere (see Figure 4.4).
+ =
Icosahedron Dodecahedron New set of Order q=0

(12 vertices) (20 vertices) 32 vertices triangular grid
Figure 4.6: Geometric interpretation of the order q=0 triangular grid. The order q=0
triangular grid is the superposition of two platonic solids (see Table 4.2).
A T1
V T2
T3
triangular grid centers of triangles construct a hexagonal grid

= hexagonal corners hexagonal grid
B
T1
T3
V
T2
order q=0 centers of triangles order q=0

triangular grid = hexagonal corners hexagonal grid
order q=1 centers of triangles order q=1

triangular grid = hexagonal corners hexagonal grid
Figure 4.7: Construction of a hexagonal grid from a triangular grid. A. Hexagonal grid
corners (gray dots) are calculated from the triangular gridpoints (black dots) via Eq. 4.3.
B. The method in A applied to the order q=0 and order q=1 triangular grids. The centers
of the hexagonal grid faces are the vertices from the corresponding triangular grid.
q=0 q=1
q=2 q=3
q= 1 and q= 2
Figure 4.8: A. Hexagonal grids, orders q=0 to q=3. Each order q hexagonal grid has 12
pentagonal faces (gray) and 30 4q 10 hexagonal faces (see Table 4.3). The q=0 hexagonal
grid is a buckeyball. B. Relationship between successive orders of hexagonal grids. The
lower-order grid (q=1) contains larger faces with * at each center; the higher-order grid
(q=2) contains smaller faces with at each center.
km2
km
km
Figure 4.9: Geometric properties of the order q=6 hexagonal grid, showing the latitude-
longitude variations of the hexagonal face areas (A), the distance between hexagonal faces
(B), and the lengths of the hexagonal face edges (C). Superimposed circles denote the 122
order q=1 triangular vertices (i.e., hexagonal face centers); circles with labels denote the 32
order q=0 triangular vertices (a subset of the q=1 vertices). Figure 4.10 shows another view
of these data. Circle in A indicates gridpoint 112 (see Figure 4.12).
Areas of hex faces Distances between hex face centers Hex face edge lengths
mean = 4.151e+03 + 1.146e+02 km mean = 6.954e+01 + 3.860e+00 km mean = 4.016e+01 + 7.011e+00 km
4400
46
A
76
B 44
C
4300
74 42
Distance (km)
40
Length (km)
Area (km )
2
4200
72 38
36
4100
70 34
32
4000
68 30
28
3900
66
26
0 2 4 6 8 10 12 0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
x 104 x 105 x 105
sorted hex face number Distance Number Edge Number
q q
(122882 faces) (368640 = 90*4 ) (368640 = 90*4 )
Figure 4.10: Geometric properties of the order q=6 hexagonal grid, showing the distributions of the
hexagonal face areas (A), the distance between hexagonal faces (B), and the lengths of the hexagonal
face edges (C). The data are the same as those shown in Figure 4.9, only here they are sorted in
increasing order. Solid line is the mean; dashed lines show the standard deviation.
u6
u5 u1
l1
u0 L1
T6
T5 T1
T0 l2
T4 T2
T3 L2
u4 u2
Hexagonal Grid u3
(order q=1)
Figure 4.11: Finite-difference approximation for the Laplacian at a hexagonal face center
T0 . The function u(, ) is defined at all face centers as u(T i ) = ui , and li and Li are arc
distances on the unit sphere. The Laplacian at T 0 is calculated via Eq. 4.12 (shown here),
where A0 is the spherical area of the target face with center T 0 on the unit sphere. As shown
in Figure 4.7, the hexagonal face centers are the vertices of the triangular grids.
-3
x 10
A
Latitude
Longitude
Del2[Ylm] Error, Order q=6 grid: error box = 1.5 deg; Order q=6 grid: error box = 1.5 deg;
max = 6.807e-03 at (lat = 7.8764, lon = -72.0000) x 10-3 11
max = 6.807e-03 at (lat = 7.8764, lon = -72.0000)
B C
10
107763 113992
9 26979
107760 113987
26977 28532
Latitude
Latitude
107761
107759 113988
8 109027 115015
112 112
108261 114639
27109 28689
108259
108260 114640
7
108266 114642
27107
108263 114644
5
-75 -74 -73 -72 -71 -70 -69
Longitude Longitude
Figure 4.12: Error in the numerical Laplacian approximation for the order q=6 hexagonal
grid. A. The error is calculated for each of the 122,882 gridpoints via Eq. 4.15, and then the
values are interpolated onto a latitudelongitude grid with cells 1 by1 . Lack of coverage
near the poles in the interpolation produces the white areas. Superimposed circles denote
the 122 order q=1 triangular vertices; circles with labels denote the 32 order q=0 triangular
vertices (a subset of the q=1 vertices). Large circle indicates gridpoint 112 (see Figure 4.9).
The color scale is saturated in order to accentuate the pattern of the error on the grid.
B. Zoom-in on the maximum value of the error, which occurs at index 112, one of the q=1
vertices. Here the color scale is the default, with the range spanning twice the maximum
value of the error. C. Configuration of hexagonal gridpoints used to calculate the error in the
numerical Laplacian at 112. See Section 4.5.
Ylm: L=6, M=1; q=0 to q=6 Ylm: L=6, M=1; q=0 to q=6
101 101
A B
Error in discrete Laplacian (D2)
Error in discrete Laplacian (D2)

100 100
92
1.
=
10-1 10-1 O
is
e
op
sl
10-2 10-2
one-norm one-norm
two-norm two-norm
inf-norm inf-norm
one-norm fit
10-3 10-3
0 1 2 3 4 5 6 102 103
Grid order q Gridpoint spacing (km)
Figure 4.13: Error in the numerical approximation of the Laplacian. The error norms are
calculated via Eqs. 4.204.22 for a Y lm (, ) test function with l = 6 and m = 1. In A the
x-axis data are plotted in terms of the grid order q; in B they are plotted on a log-scale in
terms of the average gridpoint spacing for a particular order q grid (these values are shown
in the sixth column of Table 4.4). Notice that the error in the infinity norm does not appear
to converge to zero for this test function. The slope of the curves in B gives the order of
accuracy for a particular norm. Here the line is fit using the one-norm data for the gridpoint
spacing for q=3 to q=6. See Section 4.5.
degree L=0 to L=40; q=6
0.025 A
Error in discrete Laplacian (del2inf)
0.020
0.015
0.010
0.005
0
0 5 10 15 20 25 30 35 40
order M
degree L=0 to L=40; q=6
0.025 B
Error in discrete Laplacian (del2inf)
0.020
0.015
0.010
0.005
0
0 5 10 15 20 25 30 35 40
order M
Figure 4.14: Error in the numerical approximation of the Laplacian, plotted using the
infinity norm (Eq. 4.22) for the order q=6 grid for all Y lm test functions with harmonic
degree 0 l 40 and harmonic order 0 m l. A and B compare the infinity norms for
two orientations of the order q=6 grid. A. Standard orientation of the grid (right), revealing
the five-fold symmetry (with ) about the z-axis view is looking down the z-axis (order
q=0 grid). The dips in error for the test functions having order m = 5 N (N = 0, 1, . . . , 8)
are probably due to the test functions symmetry coinciding with the symmetry of the grids.
B. Rotated orientation of the grids. Here the grids have been rotated by = 90 about the
y-axis. The relative decreases in error are no longer present. See Section 4.5.
one-norm order of accuracy: two-norm order of accuracy: inf-norm order of accuracy:
q=0 to q=1 q=0 to q=1 q=0 to q=1
40 40 40 2.0
20 20 20 1.5
L
L
0 0 0 1.0
q=1 to q=2 q=1 to q=2 q=1 to q=2
40 40 40 2.0
20 20 20 1.5
L
L
0 0 0 1.0

q=2 to q=3 q=2 to q=3 q=2 to q=3
40 40 40 2.0
20 20 20 1.5
L
L
0 0 0 1.0

q=3 to q=4 q=3 to q=4 q=3 to q=4
40 40 40 2.0
20 20 20 1.5
L
0 0 0 1.0

q=4 to q=5 q=4 to q=5 q=4 to q=5
40 40 40 2.0
20 20 20 1.5
L
0 0 0 1.0

q=5 to q=6 q=5 to q=6 q=5 to q=6
40 40 40 2.0
20 20 20 1.5
L
0 0 0 1.0
0 20 40 0 20 40 0 20 40
M M M
Figure 4.15: Order of accuracy of the numerical approximation of the Laplacian. The order
of accuracy is a function of the grid order q, the Y lm test function, and the norm that is the
gauge for the error. Thus, each colored point is a comparison between the error in the discrete
Laplacian of two spherical-hexagonal grids for a particular Y lm test function. Note that the
one-norms and two-norms are nearly second-order accurate (red) for the finest grids, while
the infinity-norms are zero- or first-order accurate (blue). See Section 4.6.
Chapter 5
Numerical solutions of waves on a

spherical membrane
Contents
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.2 The phase velocity function c(, ) . . . . . . . . . . . . . . . . . . . . . . 61
5.3 Numerical iteration for wave propagation . . . . . . . . . . . . . . . . . . 64
5.4 Numerical solutions for homogeneous and heterogeneous phase veloc-

ity maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
5.5 Stability and accuracy in the numerical model . . . . . . . . . . . . . . . 66
5.1 Introduction
In the previous chapter we developed our spherical grids and presented our numerical method
for calculating the discrete approximation of the Laplacian operator on the grids. The spatial
part of the wave equation has thus been discretized, and now it is a matter of choosing a phase
velocity function c(, ) (Section 5.2), and then iterating the wave equation (Section 5.3).
However, there are several complexities that arise in the numerical problem, most of which
result from discretizing a continuous equation into discrete spatial units and temporal units
(Section 5.5). As mentioned in the introduction, our wave equation (Eq. 3.1) is quite simple,
as it is based on a non-dispersive, isotropic spherical membrane. Nevertheless, the discrete
60
CHAPTER 5. Numerical spherical membrane waves 61
form of the wave equation, including the irregular discretization of the spatial gridpoints,
induces several of these effects, notably numerical dispersion and anisotropy (e.g., Trefethen,
1982). Thus we must take care in establishing stability and accuracy conditions, in order to
establish numerical solutions that match the analytical solutions for the homogeneous case
(c(, ) = const).
5.2 The phase velocity function c(, )
Amplitude spectrum of c(, )
For the purposes of comparisons with ray theory, we can select any reasonable function to
insert into the wave equation for c(, ). It is convenient to write c(, ) as an expansion in
spherical harmonics (see Appendix B). The amplitude spectrum of c(, ), R(l), is calculated
for each degree l via
l h
X i
R(l) = (Al0 )2 + (Alm )2 + (Blm )2 , (5.1)
m=1
where Alm and Blm are the spherical harmonic coefficients describing c(, ) (e.g., Eq. B.15).
We normalize the amplitude spectra to the maximum coefficient value for all the spectra being
compared in a given plot.
Homogeneous earth
The phase velocity for the heterogeneous earth, c(, ), is expressed as a perturbation from a
homogeneous earth reference model, c hom (, ) = const:
c(, ) c(, ) chom

= . (5.2)
chom chom
For our reference model, we use the values calculated from PREM (Dziewonski and Anderson,
1981), i.e., chom (, ) = cPREM . The amplitude spectrum of our reference model is, by defini-
tion1 , R(l) = 0; any different homogeneous earth model will have a nonzero l=0 component
(Figure 5.1).
1
In other words, the deviation of the reference model from itself is zero.
1.0
A B
0.8
Amplitude
0.6
0.4
0.2
0
0 2 4 6 8 10
-2.50 -1.90 -1.40 -0.85 -0.25 0.85 1.40 1.90 2.50 Harmonic Degree, L
Phase Velocity, Percentage
Figure 5.1: Homogeneous earth reference model. A. The phase velocity is constant on
a homogeneous earth. Phase velocity maps are plotted as a phase velocity perturbation
with respect to a homogeneous reference model. Thus, the light blue corresponds to a 0%
perturbation from the reference model. B. The amplitude spectrum of the reference model
is R(l) = 0. Any other homogeneous earth model would have a non-zero l=0 component
(dashed line).
Heterogeneous earth
By selecting an actual phase velocity map for c(, ), our problem becomes particularly rel-
evant to seismology, since ray theory has been used with these maps to make predictions of
phase and amplitude for surface wave arrivals in seismograms (e.g., Woodhouse and Wong,
1986). We will be examining three variations on the phase velocity maps of Trampert and
Woodhouse (1995): (1) the period T (Figure 5.2), (2) the truncation degree lmax of the
spherical harmonic expansion (Figure 5.3), and (3) a perturbation factor that controls the
power of the maps (Figure 5.4). Thus, with three parameters to test, we may think of the
phase velocity as
c(, ) = c(, , T, lmax, ) . (5.3)
Note that the color-scale in the phase velocity maps is always the same.
Vary the period T . Trampert and Woodhouse (1995) provide Rayleigh and Love wave
phase velocity maps for periods of 40, 60, 80, 100, 150, and 200 seconds. Figure 5.2 shows
three of these maps, and the corresponding amplitude spectra are shown in Figure 5.5A-B.
The PREM phase velocity for each map is shown in Table 5.1. For each map, the spherically
averaged phase velocity value differs from the PREM value, indicating the presence of a
nonzero l=0 component.

Upon examining Figure 5.2, we see that the shorter-period map (T =40s) contains a greater
amount of shorter-scalelength structure (l > 20). Table 5.1 shows the periods and corre-
sponding wavelengths for each of the maps of Trampert and Woodhouse (1995). In general,
the surface wave maps average the structure down to the depth of the corresponding wave-
length of their period2 . For example, the T =40s map is sensitive to lateral structure with
scalelength as short as 160 km, and it shows the averaged vertical structure down to a
depth of approximately 160 km.
Vary the truncation degree lmax. We can easily create smoother functions by truncat-
ing the spherical harmonic expansion at some lmax < 40. Figure 5.3 shows three maps, with
differing lmax values, and the corresponding amplitude spectra are shown in Figure 5.5C.
Table 5.2 shows how the minimum scalelength changes as a function of lmax.
Vary the power . We now include a perturbation factor in the phase velocity maps:

c(, )
, 01. (5.4)
chom
Thus, the maps with =0 represent the homogeneous earth reference model, and the maps
with =1 are exactly those of Trampert and Woodhouse (1995). Figure 5.3 shows three maps,
with differing values, and the corresponding amplitude spectra are shown in Figure 5.5D.
Other sets of phase velocity maps. A review of studies producing surface wave
phase velocity maps can be found in Dahlen and Tromp (1998, Section 16.9) (e.g., Ekstr
om
et al., 1997; Laske and Masters, 1996), and further discussion is presented in Trampert and
Woodhouse (2001). We also tested the maps of Trampert and Woodhouse (1996), which
incorporate more major arc data (i.e., R2 orbits) than those of TW95. Figure 5.6 shows a
map of TW96 in comparison with a map of TW95. We immediately notice that the map
of TW96 contain a greater amount of shorter-scalelength structure (Figure 5.6C). Trampert
and Woodhouse (1996) state that the models in TW95 are long wavelength approximations
(p. 22) of the models in TW96.
2
It is more accurate to talk of the sensitivity kernels of the Rayleigh waves (fundamental spheroidal mode).
These functions give more insight into how the propagating waves sample the structure with depth (see
Woodhouse, 1996).
5.3 Numerical iteration for wave propagation
The numerical iteration for wave propagation on the spherical membrane is essentially the
same as that of the planar membrane (Eqs. 2.342.36). The only difference is that here, for
the discrete Laplacian approximation 2 ut , we use D2 of Equation (4.12) instead of del2.
Furthermore, we do not require a boundary condition, only the initial conditions, which are
ut1 = ut = ut+1 = 0. Thus the iteration at each timestep is:
ut1 = ut (5.5)
ut = ut+1 (5.6)
ut+1 = 2 ut ut1 + c2 dt2 (D2(ut ) + ft ) . (5.7)
where ft = f (, , t) is the (discrete) forcing term, u t = u(, , t) is the displacement at time

t at all - gridpoints on the sphere, and D2 is the discrete Laplacian for the hexagonal grids.
The phase velocity function is entered as c(, ) = c PREM for the homogeneous case, or as a
phase velocity map for the heterogeneous case.
5.4 Numerical solutions for homogeneous and heterogeneous

phase velocity maps
This section contains numerical solutions to the wave equation on a spherical membrane
for both a homogeneous earth and a heterogeneous earth. The solution, u(, , t), can be
thought of as a global snapshot of the wave field (fixed time) or as a seismogram at a receiver
(fixed space). Figures 5.75.10 shows plots of the solutions 3 . For the homogeneous earth
reference model (Figure 5.1), we have chosen the PREM model for T =40s surface waves,
which corresponds to c(, ) = 3.9280 km/s (Table 5.1). For the heterogeneous earth model,
we have chosen the T =40s, lmax=12 (=1.0) Rayleigh wave phase velocity map of Trampert
and Woodhouse (1995) (middle map in Figure 5.3). The minimum scalelength of this map is
= 3202 km (Table 5.2). The forcing term for the solutions in this section is f = f 2 (, , t),
with = 0.0713 and = 204.5 (Figure 3.1).
In Figure 3.4 we introduced the analytical solutions to waves on the sphere for the forcing
3
Our set of data for these figures consists of 450 timesteps of solutions on the order q=4 grid, which has
7682 gridpoints. In other words, we have 7682 seismograms, one for each gridpoint.
term f2 (, , t). The solution, uf2 (, , t) (Eq. 3.32), provides a basis for comparison with
the numerical solutions for the homogeneous case. Figure 5.7 shows the agreement between
numerical and analytical seismograms for the homogeneous case. This agreement for the
homogeneous case gives us confidence that the numerical solutions for the heterogeneous case
are reasonably accurate.
In Figure 5.8 we use a heterogeneous phase velocity map, and the numerical solutions are
significantly different from the homogeneous numerical solutions, especially at later times,
when more interference of the waves has occurred. Figure 5.9 shows some of the waveforms
for two selected times.
Figure 5.10 shows a snapshot of the wavefront for a source at the North Pole. In the case
of the waves on the sphere, it is not obvious what the wavefront actually is. We define the
wavefront as the maximum value of the displacement for a given longitude (see Figure 3.4).
Thus we can plot the position of the wavefront and the amplitude of the wavefront, each as
a function of longitude. Because of the sinusoidal nature of the pulse, we keep track of both
the top and the bottom of the pulse (Figure 5.10).
For a homogeneous velocity field, and for a source at the North Pole, the wavefront will
always have constant latitude; on a - latitude-longitude grid, a wavefront propagating
between the North Pole and South Pole will be a horizontal line. The horizontal lines in
Figure 5.10A show the position of the homogeneous wavefront; the corresponding horizontal
lines in Figure 5.10B show that there is no variation in the amplitude of the top and bottom
of the pulse, as we would expect.
Some aspects of the position of the heterogeneous wavefront in Figure 5.10A are what we
might expect. First, we see that the heterogeneous pulse is leading the homogeneous pulse,
and this gap increases with time. This is due to a positive l=0 component in the phase velocity
map, meaning that the spherical average of the map is greater than the PREM value, causing
the wavefront to propagate faster than in PREM (e.g., Trampert and Woodhouse, 2001). (The
amplitude is not similarly affected by the l=0 component.) Second, we see that the leading
portion of the wavefront is along the longitude of 180 , corresponding to the high-velocity
anomaly of the Pacific basin. In general, it is difficult to predict the shape (or position) of
the wavefront, especially for a complicated phase velocity map and after many timesteps.
Scalelength and wavelength
In Figure 5.11 we show how the scalelength of the heterogeneity affects the numerical solutions
of wave propagation. We pick one receiver at 90 from the source, and we record the
seismogram for four phase velocity maps: lmax=4, lmax=8, lmax=12, and lmax=20; the
corresponding minimum scalelengths are = 8896, 4709, 3202, and 1953 km, respectively
(Table 5.2). We notice that the predominant shape of the arrivals associated with each orbit
changes very little for lmax > 8. This is because in that regime, the minimum scalelength
of the heterogeneity is less than the wavelength of the wave, which is 8300 km in this
example. In reality and in the numerical model we do not expect the wave to see
heterogeneity with scalelengths longer than the wavelength of the wave. Nevertheless, we see
that these relatively short-scalelength structures do affect the wave propagation (Figure 5.11).
Thus far we have looked at numerical results calculated using the order q=4 spherical grid
(7682 gridpoints). In Chapter 7 we return to the numerical solutions and present results for
the order q=6 grid (122,882 gridpoints), which allows for shorter wavelengths. As we discuss
in Sections 7.4 and 7.5, shorter wavelengths allow us to explore the validity of ray theory
predictions in shorter-scalelength earth models.
5.5 Stability and accuracy in the numerical model
An analysis of stability and accuracy is helpful for any numerical model (Trefethen, 1996).
An unstable solution is one in which, for example, the solution blows up exponentially.
An inaccurate solution is one in which, for example, dispersion arises due to numerical ef-
fects. Both of these are unwanted, and both can generally be predicted prior to running the
numerical program.
There is a general stability condition that states the relationship between the timestep
and spacestep for a numerical model (see Trefethen, 1996):
dx 1
R, (5.8)
dt c
where R is some critical value. Using dx = a d (a is earth radius), Equation (5.8) becomes
a d
dt R. (5.9)
c
Table 4.4 shows d for the hexagonal grids, orders q=0 to q=6. Note that if we increase c by
a factor of 4, then dt decreases by a factor of 4.
Given a uniform grid and a discretized wave equation, it is possible to calculate a value for

R, or a more complicated stability condition. By trialanderror, it appears that R = 1/ 2
is a good value for choosing the time-step dt in our problem. Larger R values will cause
the solutions to increase indefinitely; smaller R values are allowed, although they increase
computation time without decreasing error. In addition, error introduced at each time-step
may lead to a more inaccurate solution.
The inaccuracy of our numerical solutions is manifested primarily by numerical dispersion.
The effects of numerical dispersion are entirely predictable on regular grids (e.g., Trefethen,
1982) (in fact, this nuisance can be rather interesting). A numerical dispersion relationship
(for a particular FDM) will show the extent to which one would expect dispersion for a
wave with a particular wavelength. On irregular grids, however, it is more a matter of trial-
and-error in trying to determine the minimum wavelength that achieves a satisfyingly small
amount of numerical dispersion4 . We know that the dispersion is related to the wavelengths
of the waves propagating on the grid, and thus to our source parameters: the source-time
function h(t) and the initial-shape function g(, ). The width of the source-time function
must be scaled according to the time-step; thus we set the width of the function to be some
multiple of dt. Similarly the initial-shape function g(, ) must be greater than or equal to
some multiple of d. These scaling values are chosen by trial-and-error in order to minimize
the numerical dispersion in our problem. These aspects of numerical grids are generally
discussed in terms of the number of gridpoints per wavelength (gpw). Our values of 20 gpw
is considerably greater than those of Komatitsch and Tromp (2002a) (6 gpw) and Baig et al.
(2003) (7 gpw), two studies that model wave propagation using a pseudo-spectral method.
4
The study of Sword et al. (1986a) suggests that it might be possible to approximate the accuracy conditions
by assuming a regular planar grid of hexagons, and this would eliminate some of the trial and error.
Table 5.1: Wave parameters for the Rayleigh wave phase velocity maps of Trampert
and Woodhouse (1995) (also shown for T =200s.). The wavelengths are calculated via
= 2a/(l + 12 ); the phase velocity values are calculated using normal-mode summation with
PREM (see Appendix F). Each map is expanded up to lmax=40. Note that the values in
this table are properties of the probing surface waves used to create the maps.
PREM
Period, T (s) Harmonic Wavelength, Phase Velocity, c 0
approx. exact Degree, l (km) (km/s)
40 40.043 254 157.29 3.9280
60 60.147 167 238.99 3.9734
80 79.910 124 321.53 4.0236
100 100.393 97 410.57 4.0896
150 149.124 62 640.48 4.2950
200 200.837 43 920.23 4.5820
Table 5.2: Minimum scalelength in a phase velocity map expanded in spherical harmonics
up to degree lmax. The scalelengths are calculated via = 2a/(l + 12 ) (see Appendix F).
lmax Scalelength, (km)
4 8895.59
6 6158.49
8 4709.43
10 3812.40
12 3202.41
14 2760.70
16 2426.07
18 2163.79
20 1952.69
40 988.40
TW95, Rayleigh, Lmax=40, =1.00
T = 40s
T = 80s
T = 150s
-2.50 -1.90 -1.40 -0.85 -0.25 0.85 1.40 1.90 2.50

Figure 5.2: Phase velocity maps of Trampert and Woodhouse (1995) for Rayleigh waves,
periods T =40s, T =80s, and T =150s. Each map is expanded in spherical harmonics up to
degree lmax=40. Values are expressed as 100(c c 0 )/c0 , where c0 is from PREM (Table 5.1).
See Figure 5.5A-B.
TW95, Rayleigh, T=40s, =1.00
Lmax = 20
Lmax = 12
Lmax = 08
-2.50 -1.90 -1.40 -0.85 -0.25 0.85 1.40 1.90 2.50

Figure 5.3: Varying lmax on the phase velocity maps. Three different lmax truncations of
the T =40s, lmax=40 Rayleigh wave phase velocity map (Figure 5.2A). See Figure 5.5C.
TW95, Rayleigh, T=40s, Lmax=20
= 1.00
= 0.25
= 0.10
-2.50 -1.90 -1.40 -0.85 -0.25 0.85 1.40 1.90 2.50

Figure 5.4: Varying on the phase velocity maps (see Eq. 5.4). Three different values
for the T =40s Rayleigh wave phase velocity map, truncated at lmax=20 (Figure 5.2A). See
Figure 5.5D.
Varying T: Varying T:
A TW95, Rayleigh, Lmax=40, =1.00 B TW95, Rayleigh, Lmax=40, =1.00
1.0 100
T=40s T=40s
T=80s T=80s
T=150s T=150s
0.8
-1
10
0.6
Amplitude
Amplitude
10-2
0.4
10-3
0.2
-4
0 10
0 10 20 30 40 0 10 20 30 40
Harmonic Degree, L Harmonic Degree, L
Varying Lmax: Varying :

C TW95, Rayleigh, T=40s, =1.00 D TW95, Rayleigh, T= 40s, Lmax=20
1.0 100
= 1.00
= 0.25
= 0.10
0.8
-1
10
0.6
Amplitude
Amplitude
10-2
0.4
0.2
-3
10
0
0 5 10 15 20 0 5 10 15 20
Figure 5.5: Amplitude spectra for phase velocity maps of Trampert and Woodhouse (1995)
shown in Figures 5.25.4, which show variations in T , lmax, and . A. Variation in maps with
period T (Figure 5.2). Spectra are normalized to the maximum coefficient value for all the
maps (l=2 for T =40s). The decay for l > 20 is a primarily product of the inversion procedure
and is not indicative of the earth structure (Trampert and Woodhouse, 1996; Trampert and
Snieder , 1996). B. Same as A, with amplitude plotted on a logarithmic scale. C. Variation
in maps with truncation degree lmax (Figure 5.3). D. Variation in maps with (Figure 5.4).
Note that the amplitude is plotted on a logarithmic scale.
Rayleigh, T=80s, Lmax=40, =1.00
TW95
A
TW96
B
-2.50 -1.90 -1.40 -0.85 -0.25 0.85 1.40 1.90 2.50

0
10
C TW96 (T=80s)
TW95 (T=80s)
-1
10
Amplitude
TW96
-2
10
TW95
-3
10
-4
10
0 5 10 15 20 25 30 35 40
Harmonic Degree, L
Figure 5.6: Surface wave phase velocity maps for T =80s, lmax=40, and =1.0. A. Map
from Trampert and Woodhouse (1995). B. Map from Trampert and Woodhouse (1996).
C. Amplitude spectra of the two maps. Note the greater amount of shorter-scalelength struc-
ture in TW96, especially for l > 12. See Section 5.2.
R1 R2+R3 R4+R5
400
Global or b it s eve
n
d or 300
Rayleigh od
bi
displacement
source
ts
o 200
Waves at latitude +90
100
R1 (t = 3000), 0
R3 (t = 13000), receiver at
o
-100
R5, etc. latitude -30
u,
-200
-300 Latitude = +90

R2 (t = 7000), R4 (t = 17000), R6, etc.
-400
0 0.5 1.0 1.5 2.0
t, time (s) 4
x 10
R1 R2 R3 R4 R1 R2 R3 R4
400 400
300 300
displacement
displacement
200 200
100 100
0 0
-100 -100
u,
u,
-200 -200
-300 Latitude = +60 -300 Latitude = +30
-400 -400
0 0.5 1.0 1.5 2.0 0 0.5 1.0 1.5 2.0
t, time (s) 4 t, time (s) 4

x 10 x 10
R1 R2 R3 R4 R1 R2 R3 R4
400 400
300 300
displacement
displacement
200 200
100 100
0 0
-100 -100
u,
u,
-200 -200
-300 Latitude = 0 -300 Latitude = -30

-400 -400
0 0.5 1.0 1.5 2.0 0 0.5 1.0 1.5 2.0

x 10 x 10
R1 R2 R3 R4 R1+R2 R3+R4
400 400
300 300
displacement
displacement
200 200
100 100
0 0
-100 -100
u,
u,
-200 -200
-300 Latitude = -60 -300 Latitude = -90
-400 -400
0 0.5 1.0 1.5 2.0 0 0.5 1.0 1.5 2.0

x 10 x 10
Figure 5.7: Analytical and numerical solutions for seismograms at selected latitudes for the
wave propagation on the sphere using a homogeneous velocity map (c(, ) = 3.9280 km/s).
The dashed line is the analytical solution given by Eq. 3.34 for fixed (, ) and varying t;
the solid line is the numerical solution for the model. The descriptive sketch showing the
nomenclature is drawn for a seismogram at lat=30 (see Figure 1.2). The double-arrow
indicates the minimum time between surface wave arrivals in each seismogram. This value is
a minimum for a receiver at the source ( = 0 ; lat=90 ) and antipode ( = 180 ; lat=90 ),
and is a maximum for a receiver at = 90 (lat=0 ). See Section 5.4.
o o
Gridpoint 1340 : (lat, lon) = (-0.74 , 9.06 )
150
R1
100
R2
Displacement
50
-50
-100 R4
R3
-150
0 0.5 1.0 1.5 2.0
4
Time (s) x 10
o o
Gridpoint 1711 : (lat, lon) = (-0.74 , -81.06 )
150
R1
100 R2
Displacement
50
-50
R4
-100
R3
-150
0 0.5 1.0 1.5 2.0
4
Time (s) x 10
Figure 5.8: Seismograms for the wave propagation on the sphere using a heterogeneous phase
velocity map (T =40s, lmax=12, =1.0, Trampert and Woodhouse (1995)). The dashed line
is the numerical solution for the homogeneous case; the solid line is the numerical solution
for the heterogeneous case. The homogeneous solution is the same in both plots, since the
latitudes are the same. The agreement between the homogeneous analytical and numerical
solutions can be seen in Figure 5.7 for lat=0 . The variation here is due to the heterogeneity
of the phase velocity map. The gridpoint number is the index into the order q=4 triangular
grid (7862 gridpoints). See Section 5.4.
Frame 270 (time = 12417.3 s) Frame 350 (time = 16422.9 s)

300 300
200 200
= 8000 km
100 100
Displacement
Displacement
0 0
-100 -100
-200 -200
-300 -300
-90 -60 -30 0 30 60 90 -90 -60 -30 0 30 60 90
Latitude Latitude
Frame 270 (time = 12417.3 s): Longitude 0 Frame 350 (time = 16422.9 s): Longitude 0
300 300
200 200
100 100
Displacement
Displacement
0 0
-100 -100
-200 -200
-300 -300
-90 -60 -30 0 30 60 90 -90 -60 -30 0 30 60 90
Latitude Latitude
Frame 270 (time = 12417.3 s): Longitude 108 Frame 350 (time = 16422.9 s): Longitude 108
300 300
200 200
100 100
Displacement
Displacement
0 0
-100 -100
-200 -200
-300 -300
-90 -60 -30 0 30 60 90 -90 -60 -30 0 30 60 90
Latitude Latitude
Figure 5.9: Waveforms for a heterogeneous phase velocity map (T =40s, lmax=12, =1.0,
Trampert and Woodhouse (1995)). The top two plots show the agreement between numerical
(solid) and analytical (dashed) solutions for the homogeneous case. Each of the bottom four
plots compares the heterogeneous numerical solution (solid) with the homogeneous numerical
solution (dashed). See Section 5.4.
Position of Wavefront: Frame 352 (time = 16523.1 s)
90
A
60
30
Latitude
-30
-60
-90
-120 -60 0 60 120
Longitude
Height of Wavefront: Frame 352 (time = 16523.1 s)

400
300 B
200
100
Height
-100
-200
-300
-400
-120 -60 0 60 120
Longitude
Figure 5.10: Position and amplitude of a wavefront due to a source at the North Pole
(NP). The horizontal lines are wavefronts that are propagating through a homogeneous phase
velocity map (cPREM = 3.9280 km/s); the other curves are wavefronts propagating through
a heterogeneous phase velocity map (T =40s, lmax=12, =1.0, Trampert and Woodhouse
(1995)). A. Position of the wavefront, superimposed on the phase velocity map: the solid
line is the top of the pulse, the dashed line is the bottom of the pulse. At this time, the
wavefront is traveling northward, having travelled NPSPNPSP. B. Amplitude of the
top (solid) and bottom (dashed) of the wavefront. This view is looking at a cross-section of
the wavefront. See Section 5.4.
T=40s (TW95), Lmax=04, eps=1.00; q=4, thwid=20;
150
100
Displacement 50
50
100 Homogeneous Numerical Solution

Heterogeneous Numerical solution
150
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time (s) 4
x 10
150
100
50
Displacement
50
100
150
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time (s) 4
x 10
150
100
50
Displacement
50
100
150
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time (s) 4
x 10
150
100
50
Displacement
50
100
150
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time (s) 4
x 10
Figure 5.11: Variation in heterogeneous earth seismograms with varying lmax. The lmax-
value of each phase velocity map increases from lmax=4 at the top to lmax=20 at the bottom
(T =40s, =1.0). The wave has a wavelength of 8300 km. Each box designates the arrival
associated with each orbit, R1R4. See Section 5.4.
Chapter 6
Ray calculations on the sphere
Contents
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
6.2 Surface wave ray-tracing equations . . . . . . . . . . . . . . . . . . . . . . 81
6.3 Path anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
6.4 Phase anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
6.5 Amplitude anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.6 Linearized ray theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
6.7 Exact ray-tracing: Ray diagram cross-section curves . . . . . . . . . . . 91
6.8 Exact ray-tracing: Variation with , lmax, , and T . . . . . . . . . . . . 93
6.1 Introduction
In the previous section, we showed wavefronts that were produced from numerical solutions
to wave propagation on the sphere using a heterogeneous phase velocity map. The results
show us the position and the amplitude of the wavefront as a function of time (one time
frame shown in Figure 5.10). Ray theory whereby the rays are the orthogonal trajectories
to the wavefronts can be used to approximate the differences between heterogeneous and
homogeneous solutions for (1) paths, (2) amplitude, and (3) phase. In this chapter we present
the ray-tracing equations for waves on a sphere, as well as a linearized version of these
equations (Section 6.6). In Section 6.8 we show how the amplitude and phase anomaly
79
CHAPTER 6. Ray calculations on the sphere 80
predictions vary as we alter the path length , as well as one of three parameters of the
phase velocity maps c(, ): lmax, , and T (see Section 5.2).
Source and receivers
Orbit Arc distance For our purposes, we can essentially choose arbi-
R1 0 < 180 trary sources and receivers. For simplicity, we will
R2 180 < 360
look at one (fictitious) source, situated at (lat, lon)
R3 360 < 540
R4 540 < 720 = (0 , 0 ), and 25 receivers, situated approximately
= 90 from the source (Figure 6.1). This is done in
order to maximize the separation between arrivals of Rayleigh wave orbits (e.g., Figure 5.7),
which is particularly important for the numerical solutions, but not for the ray-tracing. Based
on the nomenclature for global surface waves (Figure 1.2), we can identify the arc-distance
range for each orbit (see table above). In Table 6.2 we have listed the arc distances for R1R4
for the 25 receivers in Figure 6.1.
Since we will be comparing the results from ray tracing with the results from the numerical
model, we must select receivers that are located exactly at the gridpoints in the numerical
model. The process of selecting the stations is discussed in Appendix C.
90o 4 90o o
0o 60o 24 120o 180o -120o -60o 0
20
12
60o 13 60o
7
15
11 1
30o 9 30o
10 6
3 14 2
0o 25 16 0o
8
23 -30o
-30o 17
22
-60o 18 -60o
21 19
5 -90o
-90o
0o 60o 120o 180o -120o -60o 0o
Figure 6.1: Source (star) and the 25 receivers used for the ray-tracing calculations. Receivers
are clustered around the great circle that is = 90 from the source (bold lines). The hexagon
indicates the antipode of the source. Path lengths for global Rayleigh wave orbits R1R4 are
listed in Table 6.2.
6.2 Surface wave ray-tracing equations
Ray theory was applied to surface waves by Woodhouse (1974), extended in Woodhouse and
Wong (1986), and analyzed in numerous studies (e.g., Tromp and Dahlen, 1992a,b; Wang
et al., 1993; Wang and Dahlen, 1994, 1995b; Laske et al., 1994). Below we present some
results of Woodhouse and Wong (1986) that are pertinent to our study.
For waves on a curved surface, the surface displacement field may be approximated by:
u(x, t) = A(x, t) ei (x,t) , (6.1)
where x = (x1 , x2 ) are the surface coordinates, t is time, A is the amplitude, and is the
phase (Woodhouse, 1974; Woodhouse and Wong, 1986). Woodhouse and Wong (1986) derive
equations for ray tracing and amplitude calculations on a spherical earth. Their surface
coordinates are given by x1 = cot and x2 = , where is colatitude and is longitude.
The ray-tracing equations are written as a second-order differential equation for = ()
(Woodhouse and Wong, 1986, Eq. 33):

d2 2
+ = ( ) ln c(, ) , (6.2)
d2 1 + 2
where d/d and c = /k is the phase velocity defined everywhere on the surface of
the sphere. The initial conditions (0) = 0 and (0) = 0 must be specified; they give the
starting colatitude, 0 (e.g., on the equator, 0 = cot(/2) = 0), and the take-off angle of the
ray, , via 0 = tan . Note that with = cot , we have 1/(1 + 2 ) = sin2 .
Homogeneous case: On a homogeneous earth, we have c(, ) = c hom , and the solutions
are great-circle ray paths. In this case, Equation (6.2) becomes
d2
+ = 0, (6.3)
d2
and the solutions are given by
() = 0 sin , (6.4)
which plot great-circle paths on a Mercator projection. Other forms of this solution are shown
in Appendix E, Equations (E.15)(E.19); Figure 6.2 shows some examples of great-circle ray
paths.
90
60
30
Latitude
-30
-60
-90
0 60 120 180 240 300 360
Longitude
Figure 6.2: Great-circle ray paths: ray paths for the homogeneous earth. Shown here are
plots of Eq. E.19 (Appendix E) for 35 paths, with take-off angles = 85 to +85 .
6.3 Path anomaly
On a homogeneous earth, the ray paths are great-circle paths. Once the heterogeneity is
added, the rays will deviate from great-circle paths. Figure 6.3 shows three rays for a hetero-
geneous phase velocity map, with a starting latitude of 0 (0 = 0), with starting longitude
0 = 0 , and with finishing longitude f = 180.0 . The homogeneous solutions are great
circles, and the heterogeneous solutions are curves that deviate from the great circles. The
rays are solutions to Equation (6.2), which is solved in Fortran using a fourth-order Runge
Kutta algorithm1 . Figure 6.3C shows that the deviation from the homogeneous ray paths
(i.e., great-circle paths) starts at 0 , increases to a maximum, and then decreases to 0 , where
the heterogeneous ray crosses the homogeneous ray, which is in the vicinity of the antipode.
It is helpful to modify the ray-tracing equations to find only those rays between a particular
source and receiver. The calculations are made easier by shifting the phase velocity maps so
1
Originally I wrote the ray-tracing equations of Woodhouse and Wong (1986) in Matlab and used the
differential equation solver ode45, which compares the results from a fourth-order algorithm and a fifth-order
algorithm at each step. However my Matlab program proved to be relatively slow, so in the end all the
ray-tracing was performed in Fortran.
that the source is on the equator (Appendix C). Woodhouse and Wong (1986) use a
coordinate system in which both the source and the receiver are situated on the equator, i.e.,
= 0 at longitudes 0 = 0 and f = . This imposes the boundary conditions
(0) = () = 0 (6.5)
on Equation (6.2). In Appendix E we show a generalized version of the equations, where

only the source is on the equator and the receiver is at a particular longitude f but not
necessarily on the equator.
The goal is to determine how the finishing latitude at the receiver varies as a function of
take-off angle, or, in this coordinate system, how (, 0 ) varies as 0 changes. With that
information, we can use a convergence method to find the rays that hit the receiver. This
gives rise to the notation (Woodhouse and Wong, 1986):

0 () = (, 0 ) (6.6)
0
0
0 () = (, 0 ) = (, 0 ) . (6.7)
0
Differentiating Equations (6.2) and (6.5) gives (Woodhouse and Wong, 1986, Eqs. 3839):

d2 0 2 0
+ 0 = 0
( ) ln c(, )
d2 1 + 2 1 + 2

2 0
+ 1 ( 2 ) ln c(, )
1 + 2 1 + 2

2
+ 1 0 ln c(, ) (6.8)
1 + 2
0 (0) = 0 (6.9)
0 (0) = 1 . (6.10)
Equations (6.2) and (6.8) can be solved numerically to obtain both () the ray path
and 0 ().
Using Newtons method (heterogeneous earth)
We use Newtons method to determine the right take-off angles for the ray-tracing, i.e., the
take-off angles of rays that converge on the receiver. If we shoot a number of rays from
a particular source (situated on the equator) and stop the rays at a particular longitude,
then we can plot as a function of 0 , which can be thought of as the finishing latitude,
f = 90 f , as a function of take-off angle , or simply f (). To obtain Newtons method
of approximating the root 0 , we expand 0 (n ) as a first-order Taylor series (see Appendix E)
f
0 (n ) 0 (f ) , (6.11)
f0
where n = (f ) is the finishing value for the next iteration and f is the finishing
value for the current iteration. Thus, using Equation (6.11) we can find the 0 of the ray
that converges on the receiver. It is important to note that Newtons method only works in
two cases: (1) for complicated curves when your initial guess of the root ( 0 ) is sufficiently
close to the actual root, and (2) for simple, quadratic curves, whereby your guess need not
be sufficiently close to the root.
Therefore, what is needed is the take-off angle of the ray that will finish sufficiently close to
the receiver, so that, using Newtons method, the next take-off angle brings the ray closer to
finishing at the receiver. Woodhouse and Wong (1986) used a linearized ray theory prediction
of initial take-off angle (Eq. 6.29), which worked well for the smooth (lmax=8) maps used in
their study. However, we use phase velocity maps that are much more detailed (lmax 40),
and in many cases the linearized ray theory prediction of take-off angle is far from the actual
take-off angle (of the exact ray path). In our study we shoot out a fan of rays in the direction
of the receiver, effectively creating a segmented f () curve, where f is the finishing latitude
and is the take-off angle. We then apply Newtons method over the dintervals where
there is a change of sign. This way we do not use a linearized prediction of take-off angle,
and we account for cases involving more than one path between the source and receiver (i.e.,
multipathing).
6.4 Phase anomaly
The phase is given by
Z

= ds , (6.12)
ray c(, )
where ds is an infinitesimal arc length of the ray, c(, ) is the phase velocity along the ray,
and is the (local instantaneous) angular frequency. Note that the phase is proportional
to the frequency and inversely proportional to the phase velocity. Following Woodhouse and
Wong (1986), we put the source and receiver on the equator ( = /2), with the source at
= 0 and the receiver at = . The phase and its derivative are defined at each = k
along the ray as
Z k
ds
(k ) = d (6.13)
0 c(, ) d
d ds
= . (6.14)
d c(, ) d
Homogeneous case: We have c(, ) = c hom , which gives
ds
= a (6.15)
d
d a
= (6.16)
d chom
Z k
a a
(k ) = d = k , (6.17)
0 chom chom
where a is the radius of the earth.

Heterogeneous case: We specify a phase velocity map c(, ) or c(, ):
c(, ) chom
c = c(, ) = 1 + c chom . (6.18)
chom
In this case, we have

1/2
ds ()2 1
= a + (6.19)
d [1 + 2 () ]2 1 + 2 ()
d ds
= (6.20)
d c(, ) d
Z k
ds
(k ) = d . (6.21)
0 c(, ) d
Equation (6.19) is derived in Appendix E, Equation (E.34).

The phase anomaly P along the ray is defined with respect to a homogeneous earth model:
P (k ) (k ) het (k ) hom (k ) . (6.22)

Substituting Equations (6.21) and (6.17), we obtain
P (k ) het (k ) hom (k )
Z k Z k
ds a
= d d
0 c(, ) d 0 chom
Z k
ds a
= d
0 c(, ) d chom
Z k
chom ds
= a d
chom 0 c(, ) d
Z k " 1/2 ! #
a chom ()2 1
= + 1 d , (6.23)
chom 0 c(, ) [1 + 2 () ]2 1 + 2 ()
where a k is the arc distance of the hom-ray, and ds/d is given by Equation (6.19). Equa-
tion (6.23) is presented as Equation 42 in Woodhouse and Wong (1986).
In essence it is immaterial for our purpose what we choose for c hom and , as long as they
are the same for both the ray-tracing and the wave propagation on the hexagonal grids. In
our study, the homogeneous surface wave phase velocity is calculated using the PREM model
(chom = cPREM (T ); see Appendix F)2 . We let = 1, which means that the phase anomaly
is simply a travel-time difference and has units in seconds. Values for c hom = c0 are listed in
Table 5.1.
2
A different choice for the reference model could be the mean value of the particular phase velocity map
(chom = cmean ()).
Note that with this convention, a positive phase anomaly at the receiver arises when the
ray has traveled through a slower region of phase velocity (late arrival); a negative phase
anomaly arises when the ray has traveled through a faster region (early arrival) (Table 6.1).
It is important to keep in mind that the phase is essentially averaging the phase velocity
map along the ray path. The phase at a particular point depends on the path the ray has
taken from the source. Phase velocity variation (in other words, structure) in the vicinity
of the receiver does not have as great an effect on the phase anomaly as it does on amplitude
anomaly (Um and Dahlen, 1992).
Table 6.1: Phase anomaly convention. P het hom , and 1/c.
Ray travels through a. . . phase velocity phase phase anomaly relative arrival
faster (blue) region chet > chom het < hom P <0 early
homogeneous earth chet = chom het = hom P =0 on time
slower (red) region chet < chom het > hom P >0 late
6.5 Amplitude anomaly
The amplitude anomaly A is defined as the square root of the ratio of geometrical spreading
in the homogeneous case to the geometrical spreading in the heterogeneous case:

d 1/2
d hom

A = d
, (6.24)
d het

where d is an infinitesimal take-off angle width at the source, and d is the ray-tube width at
some distance from the source (Appendix E, Figure E.1B). In other words, when the rays
passing through a heterogeneous earth are converging relative to the rays in the homogeneous
earth, the amplitude anomaly increases; when the rays are diverging, the amplitude anomaly
decreases. In Appendix E (Eq. E.28) we derive the following expression for the amplitude
anomaly:
1/4 1/4
d 1/2 1 sin2 f cos2 f 1 + f2 sin2 f
d hom

A = d
= 1/2 1/2 , (6.25)

d het |sin f | f0 1 + 02
which is a more general case than Equation 41 of Woodhouse and Wong (1986) (Eq. E.29).
Equation (6.25) is a ray-theory prediction for the amplitude anomaly due the heterogeneity
of the phase velocity maps.
6.6 Linearized ray theory
Linearized ray theory is an approximation for the (exact) ray-tracing equations in the case of
a slight heterogeneity c = c(, ) c 0 , where c0 = chom . The calculations for the anomalies
are made not along the exact ray path but along the great circle connecting the source and
receiver (i.e., = /2). Hence its alternative name, Path Integral Approximation.
The linearized ray theory offers a physical interpretation of the phase velocity maps, in
terms of phase, path, and amplitude. In particular, the linearized ray theory states that
the phase anomaly will be proportional to the integral of slowness, 1/c(, ), the path will
be proportional to c(, ), and the amplitude anomaly will be proportional to c(, )
(see Eqs. 6.28, 6.29, 6.30, and 6.31 below). Thus, the maps of the derivatives of phase
velocity in Figure 6.4 are relevant to the calculation of arrival-angle and amplitude anomalies
for the case of a source and a receiver situated on the equator. The basic idea is that the
ray path is affected only by gradients perpendicular to its propagation direction: the c and
c derivatives push and pull the ray off the great-circle path.
Phase anomaly
In the case of slight heterogeneity we can use a first-order approximation:

!
1 1 1 1 1 c
= = 1 . (6.26)
c(, ) c0 + c c0 1 + c
c0
c0 c0
In terms of c the phase can thus be written as

Z Z

= ds = ds
path c(, ) path c0 + c
Z
c
1 ds
c0 path c0
Z
a c
= 1 d , (6.27)
c0 0 c0
where we have substituted ds = a d.

Fermats principle states that the travel time (or phase) between the source and receiver
is stationary, and thus can be approximated by replacing the integral along the exact path by
the integral along the unperturbed or great-circle path. Thus, we substitute Equations (6.27)
and (6.17) into Equation (6.22) in order to get the approximation for the phase anomaly:
Plinear het hom

Z Z
a c a
= 1 d d
c0 0 c0 0 c0
Z
a c(, )
= d . (6.28)
c0 0 c0
This is Equation 43 of Woodhouse and Wong (1986). Keep in mind that in this equation the
integration path is the great circle between source and receiver, whereas in Equation (6.23)
it is over the actual ray path.
Path anomaly
The path anomaly refers to the deviation of rays from great-circle paths. The quantity gives
the slope of the ray along the ray path and may be calculated using linearized ray theory, in
order to estimate the ray path between the source and receiver. At the source and receiver,
the linearized ray theory provides estimates for the initial take-off angle and the off-azimuth
arrival angle. These equations are given by

Z
d

h i
0 = = csc sin( ) c1
0 c , d (6.29)
d =0 0 2
Z
d

h i
= = csc sin c1
0 c , d . (6.30)
d = 0 2
We should make an important distinction about the linearized ray theory prediction for
the path. The calculation of the path is performed along the great-circle (equatorial) path,
as indicated by the term c(/2, ) in the two equations above. So even though the linearized
ray theory prediction of the ray path is not a great-circle path, the calculation of it is based
on information (i.e., the transverse derivative of the phase velocity) that lies exactly on the
great-circle path. It is not accurate to say that linearized ray theory predicts a great-circle
path, but it is true that the great-circle path is the only region of the phase velocity map
that is incorporated into the linearized calculations 3 .
Laske has done extensive work in measuring and predicting off-azimuth arrival angles,
and she has used these data in refining maps of surface wave phase velocity (Laske et al.,
1994; Laske, 1995; Laske and Masters, 1996). In Figure 7 of Laske (1995), she compares the
off-azimuth arrival angle between exact ray tracing and the linearized ray tracing, noting that
the linear approximation is accurate enough to model the observed angles reliably (p. 253).
However, it is clear from the figure that the correlation decreases significantly for higher
orbits (up to R4). Furthermore, the phase velocity map used (Wong, 1989) is expanded
up to lmax=12; thus she was probably not concerned with multipathing (see Table 6.3).
An additional finding of these studies is that Laske has identified several stations whose
azimuths have been misaligned by as much as 6 , and this misalignment has been confirmed
and corrected in a couple cases.
We do not use the ray theory predictions (linear or exact) of take-off angle or arrival
angle in our comparison with the numerical results, mainly because this would be difficult to
measure with the case of the numerical model. We instead focus on ray theory predictions of
phase anomalies and amplitude anomalies.
3
Of course, for a lowlmax map with broad structures, the transverse derivatives c and c evaluated on
the great-circle path will provide accurate information over a broader region than they would for a highlmax
map with short-scalelength structures.
Amplitude anomaly
The linearized amplitude anomaly is given by (Woodhouse and Wong, 1986, Eq. 52):
Z
ln Alinear = 1
2 csc sin( ) c1
0 sin 2
cos c , d . (6.31)
0 2
The derivations and figures in Woodhouse and Wong (1986) show the behavior of the lin-
earized amplitude anomaly as a function of orbit. The anomalies of the even orbits form a
line whose slope is opposite to that of the line formed by the anomalies of the odd orbits.
6.7 Exact ray-tracing: Ray diagram cross-section curves
This section examines the variation in path, phase, and amplitude anomalies based on calcula-
tions from exact ray-tracing. Figure 6.5A-B shows a surface wave ray diagram superimposed
on a T =40s, lmax=20 (=1.0) phase velocity map. The figure shows both the homogeneous
rays (i.e., rays traced through a homogeneous earth) and the heterogeneous rays, which de-
viate from the great-circle paths. In the remainder of this section we will be looking at
cross-section curves of ray diagrams, as shown in Figure 6.5C-D. By cross-section, we mean
this: take a source and receiver that are rotated to the equator with the source at 0 = 0 .
Shoot a batch of rays over a range of take-off angles in the direction of the receiver, which
is situated at f , where f will depend on the orbit (Table 6.2). Then stop the rays at any
|| |f | (typically we stop at f ) and plot the finishing latitude as a function of take-off
angle at the longitudinal meridian . We denote the cross-section curves as f (), where
is the take-off angle and f is the finishing latitude. Similarly, we have cross-section curves
for phase anomaly and amplitude anomaly given by P f () and Af (), respectively.
The first thing that we notice in Figure 6.5 is that there is substantial deviation of the
heterogeneous rays from the great-circle paths even more than 25 with some of the
finishing rays in Figure 6.5D. The reference cross-section curve is the (red) sloping straight
line formed by the finishing homogeneous rays. Secondly, we see that multipathing occurs,
that is, multiple ray paths existing between source and receiver (again, assuming the receiver
is on the equator). The number of ray paths finishing at a particular latitude can be found
by drawing a horizontal line (f = const) through the cross-section curve and then counting
the number of intersections; we denote these points with the green vertical lines for the
equatorial finishing points with f = 0 .

Each ray path carries information about the amplitude anomaly and phase anomaly.
Thus we can plot amplitude and phase anomaly values along the cross-section curves, as in
Figure 6.5E-H. Here we note a few observations based on this example. First, notice that the
amplitude anomalies (A) approach infinity where the slope of the cross-section curve is zero,
and that A 0 in areas where the slope is steep. This is consistent with our definition of
amplitude anomaly (Eq. 6.24). Steep slopes in the cross-section curves correspond to regions
of large geometrical spreading, and thus low amplitude anomalies; shallow slopes correspond
to regions of small geometrical spreading, and thus high amplitude anomalies 4 . The regions
of infinite-amplitude anomaly are caustics.
Secondly, we note that at none of the equatorial multipathing points do we observe
near-infinite amplitudes. In other words, the multipathing rays in Figure 6.5F (vertical lines
intersecting amplitude anomalies) do not fall in the regions of near-infinite A. This observation
highlights the subtlety of the definition of caustics. As stated in the first sentence of Wang
et al. (1993), The geometrical theory of surface wave propagation on a laterally heterogeneous
earth model diverges at caustics, where neighboring rays cross. In general, multipathing
presents crossing ray paths, but not crossing neighboring paths. In other words, a caustic is
a case of multipathing, but the reverse is not necessarily true.
Thirdly, it appears that the phase anomalies (P ) exhibit less variation than do the am-
plitude anomalies. Furthermore, there appears to be a correlation between the phase and
the finishing latitude, which is particularly noticeable when the rays have undergone a rel-
atively large amount of focusing and defocusing (highlmax maps combined with long path
lengths) (Figures 6.5D,F). The ray diagram suggests why this might be occurring. Notice
that at f = 640.5 many of the rays are tending to the high-velocity anomaly (blue)
beneath the eastern Indian Ocean. To be precise, we should say that the phase anomaly
correlates to the phase velocity c(, ), which we showed in Equation (6.28). It appears that
after significant focusing, the rays have been channeled into a high velocity zone, and thus
the finishing latitudes are correlated to the finishing phase velocities and thus to the finishing
phase anomalies P . This relationship between finishing latitude and phase was only recently
4
Based on the definition of A (Eq. 6.24), we should get exact agreement between a scatter plot of A( f ())
d
versus 1/ d [ f () ]. This shows that the amplitude anomaly can be calculated directly from the cross-section
curve.
observed and warrants further study 5 .

To emphasize the importance and variation of the ray diagram cross-sections, we present
a figure of cross-section curves for the first 15 receivers (Figure 6.6). The variability in these
diagrams shows why it is necessary to calculate a range of take-off angles if one hopes to
find the ray paths between the source and receiver. Also, note that the there are, perhaps
unexpectedly, very sharp changes in the cross-section curves, which, upon closer examination,
appear as smooth curves (Figure 6.9). In Figure 6.7 we plot the difference between the
finishing latitudes of the heterogeneous and homogenous rays; this difference is a measure of
the deviation of the exact ray path from the great-circle path, which is the integration path for
linearized ray theory. We see that at the finishing (R4) longitude f , there are heterogeneous
rays that deviate from the great-circle paths by more than 30 .
Figures 6.56.7 show the exact ray-tracing results for one particular phase velocity map.
Now we present results for 15 different phase velocity maps, each one an alteration of five
maps of Trampert and Woodhouse (1995). The results, which we will discuss later, are in
Figures 6.86.11.
6.8 Exact ray-tracing: Variation with , lmax, , and T
Now we will examine how the exact ray-tracing results vary in terms of four parameters:
the path-distance6 , the truncation level of the phase velocity map lmax, the power of the
phase velocity map , and the period of the phase velocity map T . These are the parameters
discussed and presented in figures in Section 5.2. It is particularly helpful to refer to the
amplitude spectra of the phase velocity maps (Figures 5.5 and 5.6C) when analyzing these
results.
5
Intuitively, one might expect the paths which are deviating the most from the great circle to have travelled
a longer path (in distance) and therefore arrive at a later time. However, the plots of deviation versus
phase anomaly show no correlation, while plots of finishing latitude versus phase anomaly show a moderate
correlation, suggesting that the influencing factor may be a quality of the map itself.
6
is the ray-path distance between the source and the receiver, whether it be a great-circle path, a
linearized ray theory path, or the actual ray path. We discuss the variation in path length in terms of f ,
which is the finishing longitude of the rays. For the equatorial great-circle path, f = .
Variation with
Figure 6.5 shows the effect on the path, amplitude, and phase anomalies of increasing the path-
length . Deviation from the great-circle paths (i.e., Figure 6.7) increases with path-length
. Intuitively, longer paths sample a greater amount of heterogeneity, and therefore the rays
have more room for focusing and defocusing. For R4 orbits in our lmax=8 model, we find
maximum great-circle path deviations of approximately 15 . For path lengths approximately
twice as long, Woodhouse and Wong (1986) report deviations of nearly 25 in their lmax=8
model.
One way of discussing this is in terms of the variation with orbit, which was the motivation
for Woodhouse and Wong (1986). Keep in mind, however, that odd and even orbits represent
fundamentally different sets of rays, since they take off in opposite directions. Woodhouse
and Wong (1986) provide a good a summary of some of these results:
The character of amplitude and phase variations in mantle waves has been shown
to be diagnostic of substantial lateral refraction. Sometimes this is so large that
the great circle path can appear to be slow for orbits of one sense and fast for those
of the opposite sense. There is a tendency for even and odd orbit paths to polarize
and thus to sample somewhat different areas of the globe. Amplitude anomalies
are predicted to have systematic and opposite trends for even and odd orbits and
this is very commonly observed in the data. The slope of these trends is sensitive to
the great circle integral of the second transverse derivative of phase velocity, with
a range dependent kernel which also depends upon epicentral distance. (p. 772)
Although these observations were based on a comparison of actual data with ray theory
predictions (exact and linear) using relatively broad phase velocity maps (M84C of Woodhouse
and Dziewonski (1984), expanded up to lmax=8), it is possible that the general ideas will
hold for maps with a greater amount of shorter-scalelength structure. We did not make a
systematic study of the regions sampled by even and odd orbits, but we did analyze the effect
of increasing the orbit on the ray theory predictions (see Chapter 7).
Variation with lmax
The left column of Figures 6.86.11 shows the variation in the path, amplitude, and phase
with lmax. The value lmax controls the minimum-scalelength structure in the phase velocity
maps. As an infinite-frequency approximation, ray theory will see all the structure in the
maps. Thus if we decrease lmax for a given map, the deviations of the path/amplitude/phase
anomalies should decrease. This is exactly what we observe in the figures, for lmax =
40, 20, 12, 8, 4.
Variation with
The motivation for including as a testing parameter is that it allows us to investigate the
convergence of the results of exact ray theory to those of the linearized ray theory, as 0.
The middle column of Figures 6.86.11 shows the variation in the path, amplitude, and phase
with . The value controls the power of the phase velocity map. Thus if we decrease for a
given map, the deviations of the path, amplitude, and phase anomalies should decrease. This
is exactly what we observe in the figures, for = 1.0, 0.8, 0.6, 0.4, 0.2.
Variation with T
Here we show some results for phase velocity maps corresponding to different periods. The
right column of Figures 6.86.11 shows the variation in the path, amplitude, and phase with T .
As discussed in Section 5.2, each phase velocity map has different power in each part of the
amplitude spectrum (Figure 5.5A-B). The longer-period maps contain less shorter-scalelength
structure, and thus the deviation from linearized theory is less. We observe that the variation
in path/amplitude/phase anomalies decreases with increasing T = 40, 60, 80, 150, 200 s.
TW95 vs TW96
We present a comparison between the maps of Trampert and Woodhouse (1995) and Trampert
and Woodhouse (1996) in Figure 5.6. Both of the maps are expanded up to lmax=40; however,
the maps of TW96 contain more shorter-scalelength structure (i.e., greater power) for l > 12,
as shown in Figure 5.6C. As expected, this gives rise to exact ray paths deviating considerably
more from the great-circle paths for the maps of TW96, as shown in Figure 6.12. The
correlation between the finishing latitudes of TW95 and TW96 diminishes for lmax > 12.
Also note that the cross-section curves in Figure 6.12 are for the R2 orbits, not the R4 orbits,
as shown earlier (Figures 6.6 and 6.8).
Multipathing
We are interested in documenting the influence of the aforementioned parameters on multi-

pathing. Tables 6.3, 6.4, and 6.5 show some regimes under which multipathing occurs. The
trends in Table 6.3 are consistent with what we might expect: multipathing is more common
for longer path-lengths (i.e., R4) and for maps with more shorter-scalelength structure, that
is, maps with greater lmax. It is clear from this table that the current resolution of our phase
velocity maps, lmax 20 poses a significant barrier to using ray theory for orbits beyond R1.
Table 6.4 also shows what we would expect as well: as 0, multipathing ceases to occur.
The occurrence of multipathing raises an important question. If there is more than one
ray between the source and receiver, and if each ray carries separate information about the
phase and amplitude anomaly, what can we do? This question will be discussed in the next
chapter, when we introduce the numerical results for comparison with the results from ray
theory.
Table 6.2: Receivers used in the ray-tracing calculations and the numerical model (Sec-
tion 6.1, Figure 6.1). The arc distance and the orbit distances depend on the source, which
we have set at (0, 0). The receiver11 boxed value refers to the scenario for the rays plot-
ted in Figures 6.86.11; the receiver5 boxed value refers to the scenario for rays plotted in
Figure 6.12.
Receiver Arc distance Path length f

# lat lon R1 R2 R3 R4
1 39.13 -74.60 78.11 78.11 281.89 438.11 641.89
2 9.11 79.45 79.58 79.58 280.42 439.58 640.42
3 9.11 79.45 79.58 79.58 280.42 439.58 640.42
4 88.53 0.00 88.53 88.53
271.47
448.53 631.47
5 88.53 180.00 91.47 91.47 268.53 451.47 628.53
6 22.33 81.16 81.83 81.83 278.17 441.83 638.17
7 53.88 84.56 86.80 86.80 273.20 446.80 633.20
8 13.20 78.77 79.07 79.07 280.93 439.07 640.93
9 26.95 83.86 84.53 84.53 275.47 444.53 635.47
10 22.49 86.23 86.52 86.52 273.48 446.52
633.48
11 44.66 97.71 95.47 95.47 264.53 455.47 624.53
12 67.44 84.78 88.00 88.00 272.00 448.00 632.00
13 58.36 96.65 93.48 93.48 266.52 453.48 626.52
14 13.48 91.30 90.00 90.00 270.00 450.00 630.00
15 49.35 95.35 93.48 93.48 266.52 453.48 626.52
16 0.00 90.00 90.00 90.00 270.00 450.00 630.00
17 31.51 93.22 92.74 92.74 267.26 454.74 627.26
18 62.77 98.74 93.99 93.99 266.01 450.00 630.00
19 72.00 90.00 90.00 90.00 270.00 450.00 630.00
20 72.00 90.00 90.00 90.00 270.00 450.00 630.00
21 72.00 90.00 90.00 90.00 270.00 450.00 630.00
22 49.47 91.94 91.26 91.26 268.74 451.26 628.74
23 26.97 85.53 86.01 86.01 273.99 446.01 633.99
24 80.91 85.29 89.26 89.26 270.74 449.26 630.74
25 0.00 90.00 90.00 90.00 270.00 450.00 630.00
8
B
African Rift
Homogeneous Earth
Heterogeneous Earth
6
4
Latitude (deg)
-2
-4
-6
-8
0 20 40 60 80 100 120 140 160 180
Longitude, (deg)
8
C
y
6 ra
p
Deviation (deg)
to
ay
4 d le r
mid
bottom ray
2
0
0 20 40 60 80 100 120 140 160 180
Longitude, (deg)
Figure 6.3: Introductory figure for ray-tracing. Three surface wave rays superimposed
on the T =40s, lmax=20 (=1.0) phase velocity map (see Figure 5.3A). Here the source is
at (lat, lon)=(0 , 0 ), and we have shot three rays eastward with take-off angles of =
5 , 0 , +5 . The homogeneous rays (red) are great-circle paths, and they all converge at the
antipode (lat, lon)=(0 , 180 ). The heterogeneous rays (green) deviate from the great-circle
paths. A. Ray paths plotted to scale, superimposed on the outline of the land masses (Africa,
Indonesia). B. Ray paths plotted on a vertically exaggerated scale. The black line indicates
the maximum latitudinal deviation of the heterogeneous ray path from the homogeneous ray
path. This occurs for the top ray at lon=125.20 and has a value of 7.42 . C. Magnitude of
the deviation of the heterogeneous rays from the homogeneous rays. See Section 6.3.
TW95, T=40s and Lmax=12
-2.50 -1.90 -1.40 -0.85 -0.25 0.85 1.40 1.90 2.50
Percentage
-25.0 -19.0 -14.0 -8.5 -2.5 8.5 14.0 19.0 25.0
Percentage
-200 -152 -112 -68 -20 68 112 152 200
Percentage
Figure 6.4: A phase velocity map c, and its first and second derivative with respect to the
colatitude . Note the different color scales. See Section 6.6.
50
A
40
30
20
Latitude (deg)
10
-10
-20
-30
-40
-50
-600 -540 -480 -420 -360 -300 -240 -180 -120 -60 0
Longitude, (deg)
o o
R4: 100 rays, take-off angles = 150 to 210 (Lmax = 20)
50
B Homogeneous Earth
40 Heterogeneous Earth
30
20
Latitude (deg)
10
o
150
0
-10
-20
-30
-40
-50
-600 -540 -480 -420 -360 -300 -240 -180 -120 -60 0
Longitude, (deg)
PATH ( f = 105.0 ): Lmax = 20 PATH ( f = 640.5 ): Lmax = 20
50 50
40 C 40 D
30 30
Finishing Latitude (deg)

20 20
10 10
0 0
10 10
20 20
30 30
40 40
50 50
150 160 170 180 190 200 210 150 160 170 180 190 200 210
Takeoff Angle Takeoff Angle
AMPLITUDE ( f = 105.0 ): Lmax = 20 AMPLITUDE ( f = 640.5 ): Lmax = 20

10 10
E F
8 8
Amplitude Anomaly
Amplitude Anomaly
6 6
4 4
2 2
0 0
150 160 170 180 190 200 210 150 160 170 180 190 200 210
PHASE ( f = 105.0 ): Lmax = 20 PHASE ( f = 640.5 ): Lmax = 20
100 G 100 H
0 0
Phase Anomaly (s)
Phase Anomaly (s)
100 100
200 200
300 300
400 400
150 160 170 180 190 200 210 150 160 170 180 190 200 210
Figure 6.5: Ray theory predictions of path, amplitude, and phase. AB show a ray diagram
for Rayleigh waves; CH show the finishing latitudes, amplitude anomalies, and phase anoma-
lies at two specified finishing longitudes (1000 rays per plot). A. Ray diagram for Rayleigh
waves (100 rays, d = 0.4 ) superimposed on the phase velocity map of Figure 5.3A, which is
expanded up to lmax=20. The source is at (0, 0), the intersection of the Greenwich meridian
and the equator; the finishing longitudevalue for each ray is f = 640.5 , equivalent to
the R4 rays for a receiver on the equator at = +79.5 . Take-off angles range from 150
to 210 , as measured counter-clockwise from the equator. B. Same rays as in A. Red rays
show the great circles that we would expect for the first and last ray of the homogeneous
case. The two f values indicate the positions of the ray diagram cross-sections in C and
D. C. Cross-section of the ray diagram at f = 105.0 . The diagram indicates that one
distinct ray path crosses the equator at the longitude of the receiver (vertical green line).
D. Cross-section of the ray diagram at f = 640.5 . The diagram indicates that nine dis-
tinct ray paths cross the equator (vertical green lines). In other words, ray theory predicts
nine R4 rays between the source and receiver. E-F. Amplitude anomaly at f = 105.0
and f = 640.5 . G-H. Phase anomaly at f = 105.0 and f = 640.5 .
Rec 1, R4 orbit (Lmax=20) Rec 2, R4 orbit (Lmax=20) Rec 3, R4 orbit (Lmax=20)
60 60 60

40 40 40
20 20 20
0 0 0
20 20 20
40 40 40
60 60 60
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210
60 60 60

40 40 40
20 20 20
0 0 0
20 20 20
40 40 40
60 60 60
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210
60 60 60

40 40 40
20 20 20
0 0 0
20 20 20
40 40 40
60 60 60
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210
60 60 60
40 40 40
20 20 20
0 0 0
20 20 20
40 40 40
60 60 60
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210
60 60 60
40 40 40
20 20 20
0 0 0
20 20 20
40 40 40
60 60 60
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210
Takeoff Angle (deg) Takeoff Angle (deg) Takeoff Angle (deg)
Figure 6.6: Ray diagram cross-sections. A cross-section is a plot of the finishing latitude at
a specified longitudinal meridian f as a function of take-off angle (145 215 ): red is
for a homogeneous earth; blue is for a heterogeneous earth, which is indicated by the phase
velocity map (T =40s, lmax=20, =1.0). Shown here are the R4 ray diagram cross-sections
for the first 15 receivers (see Figure 6.1). Vertical green lines show the rays that converge on
the receiver, which had been rotated to the equator (finishing latitude is 0 ). See Figure 6.7
and Section 6.7.
30 30 30
20 20 20
Deviation (deg)
Deviation (deg)
Deviation (deg)
10 10 10
0 0 0
10 10 10
20 20 20
30 30 30
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210
30 30 30
20 20 20
Deviation (deg)
Deviation (deg)
Deviation (deg)
10 10 10
0 0 0
10 10 10
20 20 20
30 30 30
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210
30 30 30
20 20 20
Deviation (deg)
Deviation (deg)
Deviation (deg)
10 10 10
0 0 0
10 10 10
20 20 20
30 30 30
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210
Rec 10, R4 orbit (Lmax=20) Rec 11, R4 orbit (Lmax=20) Rec 12, R4 orbit (Lmax=20); dev = 28.29
30 30 30
20 20 20
Deviation (deg)
Deviation (deg)
Deviation (deg)
10 10 10
0 0 0
10 10 10
20 20 20
30 30 30
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210
Rec 13, R4 orbit (Lmax=20) Rec 14, R4 orbit (Lmax=20); dev = 30.92 Rec 15, R4 orbit (Lmax=20)
30 30 30
20 20 20
Deviation (deg)
Deviation (deg)
Deviation (deg)
10 10 10
0 0 0
10 10 10
20 20 20
30 30 30
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210
Figure 6.7: Path deviation from great circles. Each data point shows the finishing latitude
for a ray in the heterogeneous model minus the finishing latitude for a ray in the homogeneous
model having the same take-off angle at the source. These are obtained by subtracting the
homogeneous cross-section curves from the heterogeneous cross-section curves in Figure 6.6.
In each plot we are seeing the deviation of the R4 rays for the T =40s, lmax=20, =1.0 phase
velocity map. The maximum positive and negative deviation from great circle paths for the
15,000 ray paths in this figure can be found in the plots for receivers 14 and 12, respectively.
T=40s, Lmax=40, eps=1.00; T=40s, Lmax=20, eps=1.00; T=40s, Lmax=20, eps=1.00;
50 50 50
40 40 40

30 30 30
20 20 20
10 10 10
0 0 0
10 10 10
20 20 20
30 30 30
40 40 40
50 50 50
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210

50 50 50
40 40 40

30 30 30
20 20 20
10 10 10
0 0 0
10 10 10
20 20 20
30 30 30
40 40 40
50 50 50
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210

50 50 50
40 40 40

30 30 30
20 20 20
10 10 10
0 0 0
10 10 10
20 20 20
30 30 30
40 40 40
50 50 50
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210

50 50 50
40 40 40
30 30 30
20 20 20
10 10 10
0 0 0
10 10 10
20 20 20
30 30 30
40 40 40
50 50 50
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210

50 50 50
40 40 40
30 30 30
20 20 20
10 10 10
0 0 0
10 10 10
20 20 20
30 30 30
40 40 40
50 50 50
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210
Figure 6.8: Ray diagram cross-sections: Variation with lmax, and , and T . Each plot is for
the R4 rays between the source and receiver 11; the great-circle path length is f = 624.53 .
Left column shows variation with lmax; middle column shows variation with ; right column
shows variation with T . Thus each figure represents exact ray-tracing results for a different
phase velocity map. See Figures 6.6 and 6.9 and Section 6.8.
ZOOMIN figure: Orbit 4, Rec 11, Lmax=40, T=40s, =1.0

takeoff angle inc = 0.0701;
10
4
10
194 194.5 195 195.5 196 196.5 197
Takeoff Angle (deg)
Figure 6.9: Zoom-in on the box shown in the upper-left plot of Figure 6.8, which is the
phase velocity map with the greatest amount of short-scalelength structure (T =40s, lmax=40,
=1.0). The zoom-in box is shown to show that the kinks in the upper plots of Figure 6.8
are in fact smooth curves, and could be made smoother with more rays. There are 1000
points in each curve (i.e., 1000 rays, take-off angles, or finishing latitudes), so the take-off
angle increment is d = 70 /1000 = 0.070 . The vertical green lines indicate the point closest
to the zero and are only for graphing purposes; the actual converging rays are calculated by
using a convergence method and not by using these data.
10 10 10
8 8 8
Amplitude Anomaly
Amplitude Anomaly
Amplitude Anomaly
6 6 6
4 4 4
2 2 2
0 0 0
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210

10 10 10
8 8 8
Amplitude Anomaly
Amplitude Anomaly
Amplitude Anomaly
6 6 6
4 4 4
2 2 2
0 0 0
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210

10 10 10
8 8 8
Amplitude Anomaly
Amplitude Anomaly
Amplitude Anomaly
6 6 6
4 4 4
2 2 2
0 0 0
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210

10 10 10
8 8 8
Amplitude Anomaly
Amplitude Anomaly
Amplitude Anomaly
6 6 6
4 4 4
2 2 2
0 0 0
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210

10 10 10
8 8 8
Amplitude Anomaly
Amplitude Anomaly
Amplitude Anomaly
6 6 6
4 4 4
2 2 2
0 0 0
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210
Takeoff Angle Takeoff Angle Takeoff Angle
Figure 6.10: Same as Figure 6.8, only for amplitude anomaly. Left column shows variation
with lmax; middle column shows variation with ; right column shows variation with T . Ver-
tical green lines indicate the amplitude values for the converging rays identified in Figure 6.8.
100 100 100
Phase Anomaly (s)
Phase Anomaly (s)
Phase Anomaly (s)

0 0 0
100 100 100
200 200 200
300 300 300
400 400 400
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210
100 100 100

Phase Anomaly (s)
Phase Anomaly (s)
Phase Anomaly (s)

0 0 0
100 100 100
200 200 200
300 300 300
400 400 400
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210
100 100 100

Phase Anomaly (s)
Phase Anomaly (s)
Phase Anomaly (s)

0 0 0
100 100 100
200 200 200
300 300 300
400 400 400
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210
100 100 100

Phase Anomaly (s)
Phase Anomaly (s)
Phase Anomaly (s)
0 0 0
100 100 100
200 200 200
300 300 300
400 400 400
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210
100 100 100

Phase Anomaly (s)
Phase Anomaly (s)
Phase Anomaly (s)
0 0 0
100 100 100
200 200 200
300 300 300
400 400 400
150 160 170 180 190 200 210 150 160 170 180 190 200 210 150 160 170 180 190 200 210
Takeoff Angle Takeoff Angle Takeoff Angle
Figure 6.11: Same as Figure 6.8, only for phase anomaly. Left column shows variation with
lmax; middle column shows variation with ; right column shows variation with T . Note
the correlation with the path anomalies in Figure 6.8, particularly for the higher lmax cases.
Vertical green lines indicate the phase values for the converging rays identified in Figure 6.8.
Rec 5, R2, T=80s (TW95), Lmax=40 Rec 5, R2, T=80s (TW96), Lmax=40
50 50
30 30

10 10
10 10
30 30
50 50
150 160 170 180 190 200 210 150 160 170 180 190 200 210
Takeoff Angle (deg) Takeoff Angle (deg)
50 50
30 30

10 10
10 10
30 30
50 50
150 160 170 180 190 200 210 150 160 170 180 190 200 210
50 50
30 30
10 10
10 10
30 30
50 50
150 160 170 180 190 200 210 150 160 170 180 190 200 210
50 50
30 30
10 10
10 10
30 30
50 50
150 160 170 180 190 200 210 150 160 170 180 190 200 210
Figure 6.12: Ray diagram cross-sections (see Figure 6.6): TW95 vs TW96. Each plot is for
the T =80s, =1.0 map (Figure 5.6) for variable lmax (40, 20, 16, 12); left column is TW95,
right column is TW96. Vertical green lines indicate distinct ray paths between the source
and the receiver, which has been rotated to the equator. See Table 6.5 and Section 6.8.
Table 6.3: Dependence of multipathing on orbit and lmax (T =40s, =1.0). Each number
indicates the number of distinct ray paths for a particular source-receiver pair; a - indicates
a single ray (i.e., no multipathing). Source location and receiver locations are shown in
Figure 6.1; path lengths are listed in Table 6.2. This page shows the results for receivers 115;
the following page shows the results for receivers 1625. The horizontal circled numbers refer
to the plots in Figure 6.6; vertical circled numbers refer to the plots in the left column of
Figure 6.8. See Section 6.8.
Orbit lmax Receiver Number
01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
R1 4 - - - - - - - - - - - - - - -
6 - - - - - - - - - - - - - - -
8 - - - - - - - - - - - - - - -
10 - - - - - - - - - - - - - - -
12 - - - - - - - - - - - - - - -
14 - - - - - - - - - - - - - - -
16 - - - - - - - - - - - - - - -
18 - - - - - - - - - - - - - - -
20 - - - - - - - - - - - - - - -
40 - - - - - - - - - - - - - - -
R2 4 - - - - - - - - - - - - - - -
6 - - - - - - - - - - - - - - -
8 - - - - - - - - - - - - - - -
10 - - - - - - - - - - 3 - - - -
12 - - - - - - - - - - 3 - - - -
14 - - - - - - - - - - 3 - - - -
16 - - - - - - - - - - 3 - - - -
18 - - - - - - - - - - 3 - - - -
20 - - - - - - - - - - 3 - - - -
40 - - - - - - - - - - 3 3 - 3 -
R3 4 - - - - - - - - - - - - - - -
6 - - - - - - - - - - - - - - -
8 - - - - - - - - - - - - - - -
10 - - - - - - - - - - - - - - -
12 - - - - - - - - - - - - - - -
14 3 - - - - - - - - - - - - - -
16 3 - - 3 - - - - - - - - - - -
18 3 - - 3 - - - - - - - - - 3 -
20 3 3 - - - 3 - - - 5 - - - - -
40 3 - - 3 - 5 - 3 3 3 - 3 - 5 3

R4 4 - - - - - - - - - - - - - - -
6 - - - - - - - - - - - - - - -

8 - - - - - - - - - - - - - - -
10 - - - - - - - - - -
3 - - - -
12 - - - - 3 - 3 - - - 3 - 3 - -
14 - - - - 3 - 3 - - - 3 - 3 - -

16 - - - - 3 3 - - - - - - 3 - -
18 - - -
-
5 5 5 - -
-
5
3 3 - -

20 - - 3 3 3 5 3 - 3 3 5
5 3 - -
40 - 3 9 5 -5 7 7 - 5 3 9 11 3 5 -
Orbit lmax Receiver Number

16 17 18 19 20 21 22 23 24 25
R1 4 - - - - - - - - - -
6 - - - - - - - - - -
8 - - - - - - - - - -
10 - - - - - - - - - -
12 - - - - - - - - - -
14 - - - - - - - - - -
16 - - - - - - - - - -
18 - - - - - - - - - -
20 - - - - - - - - - -
40 - - - - - - - - - -
R2 4 - - - - - - - - - -
6 - - - - - - - - - -
8 - - - - - - - - - -
10 - - - - - - - - - -
12 - - - - - - - - - -
14 - - - - - - - - - -
16 - - - - - - - - - -
18 - - - - - - - - - -
20 - - 3 - - 3 - - - -
40 - 3 5 - - - - 3 - 3
R3 4 - - - - - - - - - -
6 - - - - - - - - - -
8 - - - - - - - - - -
10 - - - - - - - - - -
12 - - - - - - - 3 - -
14 - - - - - - - - - -
16 - - - - 3 - - - - -
18 - - 3 - - - - - - -
20 - - 3 - 5 - - - - -
40 3 3 5 5 - - - 9 5 -
R4 4 - - - - - - - - - -
6 - - - - - - - - - -
8 - - - - - - - - - -
10 - - - - - - - - - -
12 - - - - - - - - - -
14 - - - - - - - - - -
16 - - - - - 3 3 - - -
18 - - 3 - - 3 - - - 3
20 - 3 3 3 - 3 - 3 - 5
40 3 5 13 7 3 9 5 11 - 15
Table 6.4: Dependence of multipathing on , a parameter that adjusts the power of the phase
velocity map, and T , the period of the map. Each number indicates the number of distinct ray
paths for a particular source-receiver pair; a - indicates a single ray (i.e., no multipathing).
Source location and receiver locations are shown in Figure 6.1. Shown here is the number of
paths between the source and each receiver for the R4 orbit. The circled numbers refer to
the plots in the middle column () and right column (T ) of Figure 6.8. See Section 6.8.
Receiver Vary Vary T

(R4) (lmax=20, T =40s) (lmax=20, =1.0)
1.0 0.8 0.6 0.4 0.2 40 60 80 100 150
1 - - - - - - - - - -
2 - - - - - - - - - -
3 3 - - - - 3 - - - -
4 3 - - - - 3 3 3 - -
5 3 3 - - - 3 3 3 3 -
6 5 - - - - 5 - - - -
7 3 - - - - 3 - - - -
8 - - - - - - - - - -
9 3 - - - - 3 3 - - -
10
3
3 3 - -
3 3 - - -

11 5 3 - - - 5 3 - - -
12 5 5 - - - 5 5 3 - -
13 3 3 - - - 3 3 3 3 -
14 - - - - - - - - - -
15 - - - - - - - - - -
16 - - - - - - - - - -
17 3 - - - - 3 3 - - -
18 3 3 - - - 3 3 - - -
19 3 - - - - 3 - - - -
20 - - - - - - 3 - - -
21 3 3 3 - - 3 5 3 - -
22 - - - - - - - - - -
23 3 - - - - 3 3 - - -
24 - - - - - - - - - -
25 5 3 - - - 5 5 - - 3
Table 6.5: Dependence of multipathing on tw95 vs tw96, i.e., the phase velocity of maps
from Trampert and Woodhouse (1995) and Trampert and Woodhouse (1996) for the same T ,
lmax, and . Each number indicates the number of distinct ray paths for a particular source-
receiver pair; a - indicates a single ray (i.e., no multipathing). Source location and receiver
locations are shown in Figure 6.1. Shown here is the number of paths between the source and
each receiver for the R2 orbit. For the R1 orbit there are no cases of multipathing for the 25
source-receiver pairs. The circled numbers refer to the plots in Figure 6.12. See Section 6.8.
Receiver Vary lmax (T =80s, =1.0)

(R2) lmax=08 lmax=12 lmax=16 lmax=20 lmax=40
tw95 tw96 tw95 tw96 tw95 tw96 tw95 tw96 tw95 tw96
1 - - - - - - - - - 5
2 - - - - - - - - - 5
3 - - - - - - - - - -
4 - - - - - - -
- -
-

5 - - -
-
-
-
-
3
-
7

6 - - - - - 3 - - - 3
7 - - - - - 3 - 3 - 3
8 - - - - - - - 3 - 3
9 - - - - - - - - - 3
10 - - - - - - - - - 3
11 - - - - - - - - - -
12 - - - - - - - - - 3
13 - - - - - - - - - 3
14 - - - - - - - - - 5
15 - - - - - - - - - 1
16 - - - - - - - - - 5
17 - - - - - - - - - 5
18 - - - - - - - - - -
19 - - - - - - - - - 7
20 - - - - - - - - - -
21 - - - - - - - 3 - 3
22 - - - - - - - - - 3
23 - - - - - - - - - -
24 - - - - - - - - - 3
25 - - - - - - - - - 9
Chapter 7
Comparison of results from

numerical model with results from
ray theory
Contents
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
7.2 Exact ray theory vs linearized ray theory . . . . . . . . . . . . . . . . . . 115
7.3 Obtaining phase and amplitude anomalies from the numerical model . 117
7.4 Ray theory validity condition . . . . . . . . . . . . . . . . . . . . . . . . . 117
7.5 Numerical model vs ray theory . . . . . . . . . . . . . . . . . . . . . . . . 119
7.6 Summary and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
7.1 Introduction
In the previous chapter we illustrated how according to ray theory the amount of deviation
from great-circle paths depends on the distance of propagation and the properties of the
phase velocity map (lmax, T , and ). In this chapter we introduce an additional parame-
ter that is extremely important when considering the results from the numerical model: the
wavelength of the wave (Section 7.4). We compare some of the results from the ray theory
113
CHAPTER 7. Numerical model vs ray theory 114
Numerical Model
Homogeneous Earth
Seismogram
Numerical Exact Linearized

Model Ray-Tracing Ray-Tracing
Heterogeneous Earth Heterogeneous Earth Heterogeneous Earth

Seismogram Seismogram Seismogram
(Pnum & Anum) (Pexact & Aexact) (Plinear & Alinear)
Scatterplots of P, A
Figure 7.1: Flow chart for creating scatterplots comparing seismograms from three models:
numerical model, exact ray theory, and linearized ray theory. Each heterogeneous earth
seismogram is specified in terms of a phase anomaly and amplitude anomaly. The calculation
of Pnum and Anum is shown in Figure 7.6. (See Appendix G for more details.)
calculations both exact ray theory and linearized ray theory with results from the nu-
merical model, and we discuss them in term of the ray theory validity condition (Section 7.5).
The results are plotted as scatterplots correlating the phase or amplitude anomalies pre-
dicted by one theory with those predicted by another theory. Thus, each point in the scat-
terplot represents a correlation between a seismogram predicted by one method versus a
seismogram predicted by another method (Figure 7.1). In each scatterplot we show a mini-
mum of 100 data points, representing R1R4 for 25 receivers; the number of points will exceed
100 in cases of multipathing occurrences in exact ray theory predictions. The units of the
phase anomaly are seconds1 . The limits on phase anomaly in the scatterplots are always the
same; the limits on the amplitude anomaly vary from figure to figure.
In Section 7.2 we compare the results between exact ray theory and linearized ray theory.
In Section 7.5 we compare the results between the numerical model and exact ray theory.
1
Normally, phase is expressed in units of radians. With our single-frequency numerical calculations, the
traveltime difference could be converted to a phase shift by linear scaling, knowing the wavelength in the
numerical model. Any scaling, of course, would not affect the appearance of the scatterplot correlations.
7.2 Exact ray theory vs linearized ray theory
Here we will examine the agreement (and disagreement) between linearized ray theory pre-
dictions and exact ray theory predictions. In general these results should be correlated, and
in some limiting cases, they should match exactly. However, as we will see, they are poorly
correlated under many circumstances. As we did in Section 6.8, we will discuss how the
correlation varies in terms of lmax, , T , and . These are the parameters discussed and
presented in figures in Section 5.2.
Linearized ray theory is based on slight perturbations of the ray path from the great-
circle path. Thus, any effect that causes an increase in path deviation will make the results
of linearized ray theory less valid, and hence the correlation with the exact ray-tracing results
worse. In essence, the figures in this section give a different perspective on the material
presented in Section 6.8. One of the key points here is that exact ray theory may (and often
does) predict multiple ray paths, while linearized ray theory always predicts one ray path.
Variation with lmax
Figure 7.2 shows plots of exact ray theory predictions of amplitude anomalies and phase
anomalies versus predictions from linearized ray theory. With increasing lmax, the deviations
from the great-circle paths increase (Figure 6.8), and the correlation between exact and linear
ray theory worsens (Figure 7.2).
We should note three points from Figure 7.2: (1) the correlation fails at very different
lmax levels for amplitude (lmax 8) and phase (lmax 20); (2) linearized ray theory
predictions of amplitude anomalies are greater than those of exact ray theory, especially for
lmax 12; (3) exact ray theory predicts a very wide spread in phase anomalies in cases of
extreme multipathing (blue dots). It is quite remarkable how steady the correlation is for
phase; the correlation for lmax=20 is not much worse than for lmax=4. Also note that the
range of phase anomalies remains the same for each value of lmax. This suggests that the
phase anomaly is primarily sensitive to low-degree structure.
Variation with
The motivation for including as a testing parameter is that it should provide the precise
scenario in which the results of exact ray theory converge to those of the linearized ray theory,
i.e., linearized ray theory is guaranteed to be a good approximation for small enough values
of . As 0, the ray paths of the exact ray-tracing approach the great-circle paths used
with linearized ray theory, as we observed in Figure 6.8 (middle column).
Figure 7.3 confirms our hypothesis nicely: as 0, both the amplitude and phase anoma-
lies approach agreement for the two ray theories. For the phase velocity map in Figure 7.3,
the exact ray theory predictions begin to deviate from the linearized predictions when 0.4.
Variation with T
Finally, we show some results for different period phase velocity maps to illustrate the effects of
the varying degree of structure of these maps. As discussed in Section 5.2, each phase velocity
map has different power in each part of the amplitude spectrum (Figure 5.5A-B). Figure 7.4
shows how the correlation between linearized ray theory and exact ray theory increases with
increasing period. This is because the longer-period maps contain less shorter-scalelength
structure, and thus the deviation from linearized ray theory is less.
Variation with
The deviation from great-circle paths increases with path length , as we showed in Figure 6.5.
Thus, we would expect the discrepancy between the linearized calculations and the exact
calculations to increase with path length as well. Figure 7.5A-P shows the influence of path
length on the exact-vs-linear ray theory correlation. It is clear that for longer path lengths (i.e.,
higher orbits), the correlation is poorer (and is worse for higher lmax values). In particular,
for the T =40s map of Trampert and Woodhouse (1995), linearized ray theory predictions
of phase agree with exact ray theory up to orbits of R3, even for lmax=40 (Figure 7.5O).
Linearized ray theory predictions of amplitude agree with exact ray theory only for R1 for
lmax=40 (Figure 7.5I) and up to R2 for lmax=8 (Figure 7.5B).
The average phase anomaly becomes increasingly negative as the orbit increases (Fig-
ure 7.5E-H). This is a result of the non-zero (and positive) l=0 component of the phase veloc-
ity map, meaning that the heterogeneous model is on average faster than the PREM model.
As a result, the rays arrive earlier for longer paths, indicated by the increasingly negative
phase anomalies for increasing orbit (see Table 6.1). The only positive phase anomalies (i.e.,
late arrivals) occur for R1 orbits (Figure 7.5E,M). It would be helpful as an additional test to
remove the l=0 component, which would remove the bias of the phase anomaly correlation
(e.g., Trampert and Woodhouse, 2001).
7.3 Obtaining phase and amplitude anomalies from the nu-

merical model
The numerical solutions u(, , t) from the wave propagation on the hexagonal grids give both
the homogeneous case (hom) and the heterogeneous case (het). The numerical solutions at
the 25 receivers provide 25 pairs of hom-het seismograms, each containing the R1R4 arrivals.
In simplest terms, we obtain the phase and amplitude anomalies by shifting and amplifying
hom to match het; this is described further in Figure 7.6.
In general, there is some discrepancy between the het seismogram and the ray-theory-
corrected hom seismogram (Figure 7.6B). These discrepancies, generally very small, are ef-
fectively not considered in the comparison of numerical model results with ray theory re-
sults, since we specify the heterogeneous earth seismogram in terms of only two parameters,
P and A. The discrepancies are minimal when the wavelength of the input wave is much less
than the minimum scalelength of the heterogeneity (Figure 5.11), in other words, when the
ray theory validity condition is satisfied.
7.4 Ray theory validity condition
In the next section we compare phase and amplitude predictions from the numerical model
and ray theory; here we note where we expect ray theory predictions to fail. Ray theory
predictions are valid when the wavelength of the wave, , is much less than the scalelength
of the heterogeneity, , in the earth model, i.e.,
, (7.1)
which we will refer to as the validity condition.

From a practical standpoint, we would like to know how much less the wavelength needs to
be in order to obtain a good correlation between our ground-truth model (numerical model)
and exact ray theory. Furthermore, we want to examine the sensitivity of amplitude anomalies
and phase anomalies to the validity condition. These questions can be examined by approach-
ing the validity condition from either side, by varying the wavelength or the scalelength.
Changing the scalelength
The easiest way of testing the validity condition is by varying the scalelength of the phase
velocity maps. Table 5.2 shows the relationship between the lmax of a phase velocity map
and the minimum possible scalelength, min , which we will simply call the scalelength .
Changing the wavelength
We can also test the validity condition by varying the input wavelength in the numerical
model. As we mentioned in Section 5.5, the limiting factor on using shorter wavelengths in
our numerical model is the gridpoint spacing d. And based on the construction of our grids,
this value will halve for every increase in grid-order q, which will results in a factor of four
times as many gridpoints. Furthermore, we require twice as many time-steps if we halve
the gridpoint spacing, so in total we have a substantial increase in computation time when
increasing the grid order.
Figure 7.7 shows a plot of how the wavelength of the wave in the numerical model varies
as a function of the source parameters ( wid ) and . The range of values for (, ) shows
the approximate range of wavelengths that our numerical grids, orders q = 4, 5, 6 can sustain
while managing accurate solutions. For example, the order q=6 grid, containing 122,882
points, can support the shortest wavelengths: 2050 km. This is the value used in all of
the numerical results, with the exception of the top two plots in Figure 7.11.
The validity condition of Wang and Dahlen (1995b)
Wang and Dahlen (1995b) extend the notion of the validity condition in Equation (7.1) by
incorporating Fresnel areas around the ray path. They express the original validity condition,
Equation (7.1), as
s
1, (7.2)
l
where s is a calculation of the average scalelength of an earth model (Wang and Dahlen,
1995b) and l is the angular degree corresponding to the wavelength of surface waves. They
derive a new validity condition,
s
1, (7.3)
4l
which is more stringent than the condition in Equation (7.2). Wang and Dahlen (1995b)
emphasize that this condition does not depend on path length.
7.5 Numerical model vs ray theory
Finally we have reached the stage where we can compare results of our numerical model
with results from ray theory (exact or linear). The numerical model can be thought of as
our ground-truth model in the sense that it provides the most complete picture of wave
propagation, and we have demonstrated good agreement between numerical and analytical
solutions (e.g., Figure 5.7). It is the numerical model which we set as the standard to which
ray theory predictions are compared.
Figure 7.8 shows results from all three models for a very broad, lmax=6 phase velocity
map. The results from the numerical model correlate with the results from ray theory, both
exact and linearized. The phase anomalies are extremely well correlated, and the agreement
is significantly better than the amplitude anomalies. This is in support of our findings in
Chapter 6 that phase anomalies are better behaved than amplitude anomalies. Thirdly, and
very importantly, we see that the amplitude results from exact ray theory are better correlated
with the numerical model than are the results from linearized ray theory.
Figure 7.9 shows the same results as Figure 7.8, only for a map with more shorter-
scalelength structure (lmax=12). The presence of colored points indicates a regime of multi-
pathing. In the plots comparing linearized ray theory and the numerical model, the colored
points indicate the corresponding multipathing rays predicted by exact ray theory (of course,
both linearized ray theory and the numerical model calculate only a single arrival associated
with each orbit). Although in general the exact ray theory predictions are better than the
linearized ray theory predictions, we see that in some cases for example, the multipathing
rays the exact ray theory predictions are worse (Figure 7.9). The remaining scatterplots
will compare predictions from the numerical model and exact ray theory.
Variation with lmax ()
Figure 7.10 shows the systematic worsening of the correlation of amplitude anomalies and
phase anomalies between the numerical model and ray theory for increasing lmax. Beginning
with lmax=4, we have good correlation for both amplitude and phase. As we increase lmax,
thereby increasing the amount of short-scalelength structure, the correlation decreases. We
see that the ray theory predictions for amplitude break down seriously, while the results for
phase remain quite good. Thus, ray theory predictions of phase based on an lmax=20 map
are probably accurate enough for most purposes (especially for R1 and R2).
In Figure 7.10 the correlation breaks down at lmax 10 for amplitude and lmax 20
for phase. Thus, in terms of the validity condition (Eq. 7.1), the ray theory predictions of
amplitude are accurate down to scalelengths of 3812 km with 2050 km, and the
predictions of phase are accurate even when (Table 5.2).
Variation with
By using the order q=6 spherical grids, we can introduce shorter wavelengths. Based the
validity condition ( ), we would expect shorter wavelengths to increase the correlation
between the numerical results and the ray theory results (both exact and linear). Indeed, this
is what we see in Figure 7.11. The only difference between the two scenarios is the wavelength
of the input wave, and the correlation is clearly better for the shorter wavelength.
Variation with T
Exact ray theory predictions of phase and amplitude are better for longer-period phase veloc-
ity maps, which contain less shorter-scalelength structure. This is illustrated in Figure 7.12,
which shows the variation with period for lmax=12 phase velocity maps.
Variation with
Exact ray theory predictions are better for phase velocity maps with low power (Figure 7.13).
The reduced-power maps diminish the effects of focusing and defocusing, thereby keeping
the exact ray paths closer to the great-circle paths, along which the linearized ray theory
predictions are calculated.
Variation with
Figure 7.14A-X shows the dependence of the exact ray theory predictions on orbit for R1R4,
using the T =40s (=1.0) phase velocity map. Ray theory predictions of amplitude and phase
are worse for higher orbits, and this is most noticeable for higher-lmax maps, where the
validity condition is less well-satisfied (Figure 7.14Q-X). In particular, the lmax=12 map is
a regime in which it is not well-satisfied, considering the wavelengths in our numerical model
( 2050 km, 3202 km). Wang and Dahlen (1995b) emphasize that their more stringent
validity condition (Eq. 7.3) does not depend on orbit number. However, even some of their
figures show a worsening of exact ray theory predictions with increasing orbit, especially for
higher degree models2 (lmax > 20).
7.6 Summary and conclusions
Our results are based on global surface wave Rayleigh wave orbits up to R4, i.e., path lengths
up to 640 (Table 6.2). Here results from three models are summarized: numerical
model, exact ray-tracing, linearized ray-tracing. We also raise some questions and point out
directions of future work.
1. Numerical model
2
See Figures 9 and 11 in Wang and Dahlen (1995b). The error of the ray theory predictions increases with
orbit, and is especially noticeable for Love waves, G1G4. Their ground-truth calculations are based on a
coupled-mode method described by Um and Dahlen (1992).
(a) With a non-standard distribution of gridpoints on a sphere (Section 4.2), and with
a relatively simple finite-difference method (Figure 4.11), we are able to produce
numerical solutions to the wave equation on a spherical membrane. The numerical
model offers an appealing approach to examining surface wave propagation, and
for comparing its results with predictions from ray theory.
(b) Our numerical solutions to wave propagation on a spherical membrane agree with
analytical solutions for the homogeneous case (Figure 5.7). This gives us confidence
in the numerical solutions for the heterogeneous case.
(c) The minimum number of gridpoints per wavelength on our spherical grids could
be improved by optimizing the distribution of gridpoints, or perhaps by using a
better finite difference approximation to the Laplacian. This would allow us to use
shorter wavelengths in our grids without going to finer grids.
2. Ray theory
(a) Linearized ray theory results i.e., phase and amplitude anomalies are cal-
culated along the great-circle path between the source and receiver. Exact ray
theory results are calculated along the actual ray path, and thus when the actual
path deviates from the great-circle path, the exact results deviate from the lin-
earized results. Deviation from linearized ray theory increases with: (1) increasing
path-length , (2) increasing lmax, (2) increasing , and (4) decreasing period T .
(b) Multipathing multiple ray paths between the source and receiver tends to
be the norm for path lengths longer than R1 and for our current resolution of the
phase velocity maps (lmax=20) (Tables 6.36.5). Multipathing is an indication
of the divergence between linearized ray theory and exact ray theory; in a regime
with no multipathing, linearized ray theory predictions agree with exact ray theory
predictions.
(c) Ray diagram cross-section curves, f (), show the variation of the finishing latitude
as a function of take-off angle. The variability of these curves is surprisingly sharp
and almost kinked (Figure 6.9). It is not obvious what controls the sharpness
of the cross-section curves. Certainly the path length and the scalelength of the
structure (lmax) play a role, but this topic warrants further analysis.
3. Combined results
(a) Exact ray theory is better than linearized ray theory at predicting amplitude and
phase anomalies for lmax < 10 phase velocity maps (Figure 7.8). In this regime
multipathing is rare.
(b) Linearized ray theory may be better than exact ray theory at predicting amplitude
and phase anomalies in regimes of extreme multipathing (e.g., lmax 20 phase
velocity maps), where exact ray theory predicts seemingly implausible results (e.g.,
Figure 7.2). Based on the phase velocity maps used here, it appears that in the
regime < 0.4, linearized ray theory predictions of amplitude and phase agree with
those from exact ray theory (Figure 7.3).
(c) We check the validity of the predictions of ray theory by altering the two parameters
of the validity condition, , where is the wavelength in the numerical model
and is the scalelength of the phase velocity map (Figures 7.10 and 7.11). Our
results show that the validity condition is more stringent for amplitudes than for
phase (Figure 7.10), something we might have expected based on the results from
linearized ray theory (Figure 7.2).
Physically, in the regimes approaching the validity condition ( ) or outside
the validity condition ( > ), the very-long-period waves in the numerical model
effectively smooth over some of the heterogeneity in the earth model, but with
ray theory no such smoothing occurs (Figure 7.15).
(d) In other words, ray theory is particularly sensitive to the smallest scalelengths in
the model, whereas physical waves are not they tend to see only elements
of the model comparable with or larger than their wavelength. It may be that
a smoothed or filtered version of the model (i.e., not simply a truncated lmax
version) could yield more stable and accurate results in ray theory, a conjecture
that remains for future work.
4. Extending this study
(a) An increase in computing power would allow us to compare results from ray-
tracing with those from the numerical model for higher lmax phase velocity maps.
With more computing power, we can use finer grids and thus propagate shorter-
wavelength waves in the numerical model. These waves need to be sufficiently
shorter than the scalelength structure, in order to approach the infinitely-short
wavelength approximation employed by ray theory.
(b) With finer grids, we would enter the regimes of multipathing (lmax > 12) and
perhaps find a way to utilize the information from the multiple ray paths. It could
be that, with a fine enough numerical grid, we could identify the surface wave
arrivals corresponding to each ray path.
(c) Real surface waves, such as those in Figure 1.3, include a range of frequencies. The
wave propagation in our numerical model is non-dispersive (i.e., phase velocity, c,
is independent of frequency). It would be desirable to introduce dispersion into
the numerical model, in order to model surface wave properties more realistically.
(d) It is also possible to use ray theory to create more realistic waveforms. For exam-
ple, we could combine the phase and amplitude anomaly results from all periods
to make a single waveform. Realistic waveforms from the numerical model or
waveforms from ray theory could be compared with observed seismograms (e.g.,
Woodhouse and Wong, 1986; Laske, 1995).
Amplitude Anomaly (RMS = 0.8989) Phase Anomaly (RMS = 3.6982)

T=40s (TW95), Lmax=20, eps=1.00; R1 to R4; T=40s (TW95), Lmax=20, eps=1.00; R1 to R4;
20
0
1
Exact [P (rad)]
Exact [Ln(A)]
20
0
40
1
2 60
2 1 0 1 2 60 40 20 0 20
Linear [Ln(A)] Linear [P (rad)]

20
0
1
Exact [P (rad)]
Exact [Ln(A)]
20
0
40
1
2 60
2 1 0 1 2 60 40 20 0 20

20
0
1
Exact [P (rad)]
Exact [Ln(A)]
20
0
40
1
2 60
2 1 0 1 2 60 40 20 0 20
Figure 7.2: Vary lmax. Correlation between linear and exact ray theory for varying the lmax
value of the phase velocity maps (T =40s, =1.0). Each plot is of the exact ray theory predic-
tions vs the linearized ray theory predictions; left column is for amplitude anomalies, right is
for phase anomalies. The lmax value decreases from lmax=40 at the top to lmax=4 at the
bottom. The correlation increases as lmax decreases. Note the multipath rays for the exact
ray tracing: black = 1 ray, red = 3 rays, green = 5 rays, blue = 7 rays. See Section 7.2.
0.5 20
0
0.25
Exact [P (rad)]
Exact [Ln(A)]
20
0
40
0.25
60
0.5
0.5 0.25 0 0.25 0.5 60 40 20 0 20

0.5 20
0
0.25
Exact [P (rad)]
Exact [Ln(A)]
20
0
40
0.25
60
0.5
0.5 0.25 0 0.25 0.5 60 40 20 0 20

0.5 20
0
0.25
Exact [P (rad)]
Exact [Ln(A)]
20
0
40
0.25
60
0.5
0.5 0.25 0 0.25 0.5 60 40 20 0 20
Figure 7.3: Vary . Correlation between linear and exact ray theory for varying the value of
the phase velocity maps (T =40s, lmax=6). Each plot is of the exact ray theory predictions vs
the linearized ray theory predictions; left column is for amplitude anomalies, right is for phase
anomalies. The value increases from =1.0 at the top to =0.1 at the bottom. As expected,
the correlation increases as 0. There is no multipathing here, due to the relatively low
lmax=6 of the map. See Section 7.2.

1.5 20
1
0
0.5
Exact [P (rad)]
Exact [Ln(A)]
20
0
0.5 40
1
60
1.5
1.5 1 0.5 0 0.5 1 1.5 60 40 20 0 20

1.5 20
1
0
0.5
Exact [P (rad)]
Exact [Ln(A)]
20
0
0.5 40
1
60
1.5
1.5 1 0.5 0 0.5 1 1.5 60 40 20 0 20

1.5 20
1
0
0.5
Exact [P (rad)]
Exact [Ln(A)]
20
0
0.5 40
1
60
1.5
1.5 1 0.5 0 0.5 1 1.5 60 40 20 0 20
Figure 7.4: Vary T . Correlation between linear and exact ray theory for varying the period T
of the phase velocity maps (lmax=12, =1.0). Each plot is of the exact ray theory predictions
vs the linearized ray theory predictions; left column is for amplitude anomalies, right is for
phase anomalies. The period T increases from T =40s at the top to T =150s at the bottom.
The correlation increases as T increases. Note the multipath rays for the exact ray tracing:
black = 1 ray, red = 3 rays. See Section 7.2.
Amplitude: Lmax=08, R1 Amplitude: Lmax=08, R2 Amplitude: Lmax=08, R3 Amplitude: Lmax=08, R4
1.5 A 1.5 B 1.5 C 1.5 D

1.25 1.25 1.25 1.25
RAYexact
RAYexact
RAYexact
RAYexact
1 1 1 1
0.75 0.75 0.75 0.75
0.5 0.5 0.5 0.5
0.5 0.75 1 1.25 1.5 0.5 0.75 1 1.25 1.5 0.5 0.75 1 1.25 1.5 0.5 0.75 1 1.25 1.5
RAYlinear RAYlinear RAYlinear RAYlinear
Phase: Lmax=08, R1 Phase: Lmax=08, R2 Phase: Lmax=08, R3 Phase: Lmax=08, R4

100 100 100 100
E F G H
0 0 0 0
RAYexact
RAYexact
RAYexact
RAYexact
100 100 100 100
200 200 200 200
300 300 300 300
400 400 400 400
400 300 200 100 0 100 400 300 200 100 0 100 400 300 200 100 0 100 400 300 200 100 0 100
4 I 4 J 4 K 4 L
3 3 3 3
RAYexact
RAYexact
RAYexact
RAYexact
2 2 2 2
1 1 1 1
0 0 0 0
0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4

100 100 100 100
M N O P
0 0 0 0
RAYexact
RAYexact
RAYexact
RAYexact
100 100 100 100
200 200 200 200
300 300 300 300
400 400 400 400
400 300 200 100 0 100 400 300 200 100 0 100 400 300 200 100 0 100 400 300 200 100 0 100
Figure 7.5: Vary orbit (and lmax). Comparison of amplitude anomalies (rows 1 and 3) and
phase anomalies (rows 2 and 4), as predicted from exact ray theory (RAYexact) and from
linearized ray theory (RAYlinear). Each column shows the results for a particular orbit, R1
(left column) to R4. Each phase velocity map is T =40s (=1.0) for two lmax values. Note
the multipath rays for the exact ray tracing: black = 1 ray, red = 3 rays, green = 5 rays,
blue = 7 rays. Each scatterplot contains a minimum of 25 points. A-H. Variation
with orbit for lmax=8 for amplitude anomalies (A-D) and for phase anomalies (E-H).
I-P. Same as A-H, only for lmax=40. See Section 7.2.
A 150
Homogeneous Numerical Solution

R1
Heterogeneous Numerical Solution
100
R2
50
Displacement
-50
R4
-100
R3
-150
0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
4
x 10
Time (s)
B R4 arrivals: HOM = red, HET = blue (fixed)

Initial Data, RMS = 25.1975;
R4 arrivals: HOM = red, HET = blue (fixed)
Final Data, RMS = 3.3818;
150 150
P P = -176.80s; A = 0.8257
100 A 100
Displacement
Displacement
50 50
0 0
-50 -50
-100 -100
-150 -150
1.6 1.7 1.8 1.9 1.6 1.7 1.8 1.9
4 4
Time (s) x 10 Time (s) x 10
Figure 7.6: Obtaining the amplitude and phase anomalies from the numerical model.
A. Seismogram at a receiver approximately 90 from the source. Boxes show the time in-
tervals picked to designate the impulse associated with each orbit. In this example, the
extremely-long period wave propagated through the T =40s, lmax=8 (=1.0) phase velocity
map. B. Method for determining the phase anomaly and amplitude anomaly in the numerical
model, shown here for the R4 arrival in A. A non-linear least-squares program calculates the
phase shift P (seconds) and amplitude change A that best fits the homogeneous earth seismo-
gram (red) to the heterogeneous earth seismogram (blue). Minor discrepancy remains between
the fitted seismogram and the heterogeneous numerical solution, as the two parameters P and
A do not completely describe the heterogeneous numerical solution (see Section 7.3).
= 9576 km
10000
8000

6000
Wavelength,
4000
2000
0.06 = 1846 km
0.04 200
Sp 150
at
ial
W 0.02 100
idt
h,
50 id th,
lW
0 pora
0 Tem
Figure 7.7: Wavelength (, ) as a function of spatial source width and temporal source
width . Each point represents the wavelength of the wave pulse when it reaches the equator,
having originated at the North Pole. It is calculated by taking the Fourier transform of the
pulse for the analytical solution u f2 (, , t) (Eq. 3.46); the corresponds to the frequency
at the maximum of the frequency spectrum. In the point source limit, when 0 and
0, we see that 0, as we would expect. The values shown here cover the range of
wavelengths used in the numerical calculations. See Section 7.4.
Amplitude Anomaly: HEX and RAY Phase Anomaly: HEX and RAY
T=40s (TW95), Lmax=06, eps=1.00; q=6, thwid=5; R1 to R4; T=40s (TW95), Lmax=06, eps=1.00; q=6, thwid=5; R1 to R4;
1.7
100
0
1.35
100
RAYexact
RAYexact
1
200
0.65 300
400
0.3
0.3 0.65 1 1.35 1.7 400 300 200 100 0 100
HEX HEX
1.7
100
0
1.35
RAYlinear
RAYlinear 100
1
200
0.65 300
400
0.3
0.3 0.65 1 1.35 1.7 400 300 200 100 0 100
HEX HEX
Amplitude Anomaly: RAY Phase Anomaly: RAY
1.7
100
0
1.35
100
RAYexact
RAYexact
1
200
0.65 300
400
0.3
0.3 0.65 1 1.35 1.7 400 300 200 100 0 100
RAYlinear RAYlinear
Figure 7.8: Comparison of phase and amplitude anomalies, as predicted from the numerical
model (HEX) and from ray tracing, including both exact ray theory (RAYexact) and lin-
earized ray theory (RAYlinear). Each plot includes calculations carried out using an lmax=6
phase velocity map (T =40s, =1.0). Note that the exact ray theory results are a better predic-
tor of amplitude anomalies and phase anomalies than linearized ray theory. See Section 7.5.
2.5
100
2 0
100
RAYexact
RAYexact
1.5
200
1
300
0.5
400
0
0 0.5 1 1.5 2 2.5 400 300 200 100 0 100
HEX HEX
2.5
100
2 0
100
RAYlinear
RAYlinear
1.5
200
1
300
0.5
400
0
0 0.5 1 1.5 2 2.5 400 300 200 100 0 100
HEX HEX
Amplitude Anomaly: RAY Phase Anomaly: RAY
2.5
100
2 0
100
RAYexact
RAYexact
1.5
200
1
300
0.5
400
0
0 0.5 1 1.5 2 2.5 400 300 200 100 0 100
RAYlinear RAYlinear
Figure 7.9: Same as Figure 7.8, only for an lmax=12 phase velocity map (T =40s, =1.0).
Note the multipath rays for the exact ray tracing: black = 1 ray, red = 3 rays. See Sec-
tion 7.5.
3
Lmax = 20 ( = 1953 km ) 100 Lmax = 20 ( = 1953 km )
2
100
RAYexact
RAYexact
200
1
300
400
0
0 1 2 3 400 300 200 100 0 100
HEX HEX
3
Lmax = 12 ( = 3202 km ) 100 Lmax = 12 ( = 3202 km )
2
100
RAYexact
RAYexact
200
1
300
400
0
0 1 2 3 400 300 200 100 0 100
HEX HEX
3
Lmax = 8 ( = 4709 km ) 100 Lmax = 8 ( = 4709 km )
2
100
RAYexact
RAYexact
200
1
300
400
0
0 1 2 3 400 300 200 100 0 100
HEX HEX
3
Lmax = 4 ( = 8896 km ) 100 Lmax = 4 ( = 8896 km )
2
100
RAYexact
RAYexact
200
1
300
400
0
0 1 2 3 400 300 200 100 0 100
HEX HEX
Figure 7.10: Vary lmax. Comparison of amplitude anomalies (left column) and phase
anomalies (right column), as predicted from the numerical model (HEX) and from exact ray
tracing (RAYexact). The lmax value of each phase velocity map varies from lmax=20 at
the top to lmax=4 at the bottom (T =40s, =1.0), and the corresponding scalelength ()
is shown. The input wave into the numerical model remains fixed ( 2050 km), and the
correlation worsens with increasing lmax. Note the multipath rays for the exact ray tracing:
black = 1 ray, red = 3 rays, green = 5 rays. See Section 7.5.
= 8300 km 100 = 8300 km

1.5
0
1.25
100
RAYexact
RAYexact
1
200
0.75
300
0.5 400
0.5 0.75 1 1.25 1.5 400 300 200 100 0 100

HEX HEX
= 2050 km 100 = 2050 km

1.5
0
1.25
100
RAYexact
RAYexact
1
200
0.75
300
0.5 400
0.5 0.75 1 1.25 1.5 400 300 200 100 0 100

HEX HEX
Figure 7.11: Vary . Comparison of amplitude anomalies (left column) and phase anomalies
(right column), as predicted from the numerical model (HEX) and from exact ray tracing
(RAYexact). The input wavelength into the numerical model is 8300 km at the top
(order q=4 numerical grids) and 2050 km at the bottom (order q=6 grids). The phase
velocity map is the same: the lmax=6 map (T =40s, =1.0), which contains a minimum
scalelength of 6158 km. As we go to finer grids and use shorter wavelength waves in the
numerical model, the correlation between exact ray theory and the numerical results increases.
No multipathing occurs, due to the low lmax=6 of the map.
3 100
0
2
RAYexact
RAYexact
100
200
1
300
400
0
0 1 2 3 400 300 200 100 0 100
HEX HEX
3 100
0
2
RAYexact
RAYexact
100
200
1
300
400
0
0 1 2 3 400 300 200 100 0 100
HEX HEX
3 100
0
2
RAYexact
RAYexact
100
200
1
300
400
0
0 1 2 3 400 300 200 100 0 100
HEX HEX
3 100
0
2
RAYexact
RAYexact
100
200
1
300
400
0
0 1 2 3 400 300 200 100 0 100
HEX HEX
Figure 7.12: Vary T . Comparison of amplitude anomalies (left column) and phase anomalies
(RAYexact). The period of each phase velocity map varies from T =40s at the top to T =150s
at the bottom (lmax=12, =1.0). The input wave into the numerical model remains fixed
( 2050 km), and the correlation worsens with decreasing T . Note the multipath rays for
the exact ray tracing: black = 1 ray, red = 3 rays. See Section 7.5.
3 100
0
2
RAYexact
RAYexact
100
200
1
300
400
0
0 1 2 3 400 300 200 100 0 100
HEX HEX
3 100
0
2
RAYexact
RAYexact
100
200
1
300
400
0
0 1 2 3 400 300 200 100 0 100
HEX HEX
3 100
0
2
RAYexact
RAYexact
100
200
1
300
400
0
0 1 2 3 400 300 200 100 0 100
HEX HEX
3 100
0
2
RAYexact
RAYexact
100
200
1
300
400
0
0 1 2 3 400 300 200 100 0 100
HEX HEX
Figure 7.13: Vary . Comparison of amplitude anomalies (left column) and phase anomalies
(RAYexact). The value of each phase velocity map varies from =1.0 at the top to =0.1
at the bottom (lmax=12, T =40s). The input wave into the numerical model remains fixed
( 2050 km), and the correlation worsens with increasing . Note the multipath rays for
the exact ray tracing: black = 1 ray, red = 3 rays. See Section 7.5.

1.2 1.2 1.2 1.2
A B C D
1.1 1.1 1.1 1.1
RAYexact
RAYexact
RAYexact
RAYexact
1 1 1 1
0.9 0.9 0.9 0.9
0.8 0.8 0.8 0.8

0.8 0.9 1 1.1 1.2 0.8 0.9 1 1.1 1.2 0.8 0.9 1 1.1 1.2 0.8 0.9 1 1.1 1.2
E F G H
100 100 100 100
0 0 0 0
RAYexact
RAYexact
RAYexact
RAYexact
100 100 100 100
200 200 200 200
300 300 300 300
400 400 400 400
400 300 200 100 0 100 400 300 200 100 0 100 400 300 200 100 0 100 400 300 200 100 0 100
1.5
I 1.5
J 1.5
K 1.5
L
1.25 1.25 1.25 1.25
RAYexact
RAYexact
RAYexact
RAYexact
1 1 1 1
0.75 0.75 0.75 0.75
0.5 0.5 0.5 0.5
0.5 0.75 1 1.25 1.5 0.5 0.75 1 1.25 1.5 0.5 0.75 1 1.25 1.5 0.5 0.75 1 1.25 1.5
M N O P
100 100 100 100
0 0 0 0
RAYexact
RAYexact
RAYexact
RAYexact
100 100 100 100
200 200 200 200
300 300 300 300
400 400 400 400
400 300 200 100 0 100 400 300 200 100 0 100 400 300 200 100 0 100 400 300 200 100 0 100
2 Q 2 R 2 S 2 T
1.5 1.5 1.5 1.5
RAYexact
RAYexact
RAYexact
RAYexact
1 1 1 1
0.5 0.5 0.5 0.5
0 0 0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2 0 0.5 1 1.5 2 0 0.5 1 1.5 2
U V W X
100 100 100 100
0 0 0 0
RAYexact
RAYexact
RAYexact
RAYexact
100 100 100 100
200 200 200 200
300 300 300 300
400 400 400 400
400 300 200 100 0 100 400 300 200 100 0 100 400 300 200 100 0 100 400 300 200 100 0 100
HEX HEX HEX HEX
Figure 7.14: Vary orbit (and lmax). Comparison of amplitude anomalies (rows 1, 3, 5)
and phase anomalies (rows 2, 4, 6), as predicted from the numerical model (HEX) and from
exact ray tracing (RAYexact). Each column shows the results for a particular orbit, R1 (left
column) to R4. Each phase velocity map is T =40s (=1.0) for three lmax values (4, 8, 12).
The input wave into the numerical model remains fixed ( 2050 km). Note the multipath
rays for the exact ray tracing: black = 1 ray, red = 3 rays. A-H. Variation with orbit for
lmax=4 for amplitude anomalies (A-D) and for phase anomalies (E-H). I-P. Same as A-H,
only for lmax=8. Q-X. Same as I-P, only for lmax=12. See Section 7.5.
10000
L=2
L=4 A 1.0 B
used in the q=4 grids
(km)
8000
0.8
Amplitude
L=6 L=0
6000
L=4
Scalelength,
0.6
q=5
L=8
4000
L=12 0.4
q=6
2000 0.2
0 0
0 5 10 15 20 25 30 35 40
0 5 10 15 20 25 30 35 40
Figure 7.15: The relationship between the detail of a phase velocity map and the wavelength
of the probing surface waves. A. Scalelength as a function of spherical harmonic degree
l: = 2a/(l + 12 ), a is the earths radius (Appendix F). A phase velocity map of degree
lmax=40 may have structure with scalelength as short as 1000 km; a map of degree lmax=4
may have structure with scalelength as short as 9000 km (Table 5.2). Also shown are the
approximate wavelengths of the waves in our numerical model, where q indicates the order of
the spherical grid (Figure 7.7). This is shown in order to facilitate a discussion in terms of
the ray theory validity condition, (Sections 7.4 and 7.5). B. Amplitude spectrum of
the T =40s phase velocity map in Figure 5.2A (lmax=40, =1.0). The power at each degree l
is the most important factor, not simply the lmax truncation of the map. In this map, most
of the structure is within the range l < 20.
Appendix A
Functions on the sphere
In this appendix we examine what it means to be a well-defined function on the sphere.

When deriving the analytical solutions to the wave equation on the sphere (Chapter 3), we
must ensure that all of our functions on the sphere satisfy the invariance conditions presented
in Equations (A.11)(A.14). Spherical harmonic functions, the subject of Appendix B, all
satisfy the invariance conditions for the sphere.
A well-defined function on the sphere
Let F : R2 S 2 be the spherical coordinate map-

ping, i.e.,
F (, ) = (sin cos , sin sin , cos ) , (A.1)
where is the colatitude (polar angle) and is the

Figure A.1: Function mappings de- azimuthal angle.
fined in text. S 2 is the set of all
Definition. Let Z be any set and let u : R 2 Z.
points on the unit sphere. R 2 is the
set of all points in the - plane. Z Then u is well-defined on the sphere S 2 if there is a
is a set of real numbers.
: S 2 Z such that u = u
function u F . (It is really
u
that we care about.) These mappings are shown schematically in Figure A.1.
Theorem. u is well-defined on S 2 iff
F (x1 ) = F (x2 ) implies u(x1 ) = u(x2 ) (x1 , x2 R2 ) . (A.2)
138
APPENDIX A. Functions on the sphere 139
. Let y S 2 . Then y = F (x1 ) for some x1 R2 . Define u

Idea: Try to define u (y) =
u(x1 ). But what if also y = F (x2 )? If u
is to be a function, then you need u(x 1 ) = u(x2 ). In
other words, Equation (A.2) is exactly the condition that you need.
Basic lemma from trigonometry:
r1 cos 1 = r2 cos 2 (A.3)
r1 sin 1 = r2 sin 2 (A.4)
if and only if one of the three conditions is met:
(1) r 1 = r2 = 0 (A.5)
(2) r 1 = r2 and 2 = 1 + 2n (A.6)
(3) r1 = r2 and 2 = 1 + + 2n . (A.7)
Now it is just a matter of using this lemma to see when F (x 1 ) = F (x2 ). Let us suppose
that F (1 , 1 ) = F (2 , 2 ). Using Equation (A.1), we have
sin 1 cos 1 = sin 2 cos 2 (A.8)
sin 1 sin 1 = sin 2 sin 2 (A.9)
cos 1 = cos 2 . (A.10)
This is our basic lemma (Eqs. A.3A.4), with r 1 = sin 1 and r2 = sin 2 . It is true if and
only if at least one of the three conditions (Eqs. A.5A.7) is met.
Equation (A.5) gives us
(1) sin 1 = sin 2 = 0 (and cos 1 = cos 2 )
1 = 2n and 2 = 2m
or 1 = + 2n and 2 = + 2m .
(2) sin 1 = sin 2 and 2 = 1 + 2n
2 = 1 + 2m .
(3) sin 1 = sin 2 and 2 = 1 + + 2n
sin 1 = sin(2 )
cos 1 = cos 2 = cos(2 )
1 = 2 + 2m
2 = 1 + 2m .
Putting this all together, we can state the following:
A function u is well-defined on the sphere if and only if
(1) u(, + ) = u(, ) (A.11)
(2) u( + 2, ) = u(, ) (A.12)
(3) u(0, ) = const1 uNP (A.13)
(4) u(, ) = const2 uSP . (A.14)
Equations (A.11)(A.14) are the four invariance conditions.
Extending a function to the - plane
Let us now say we have a function ge(, ) defined only in the domain R 2 D : 0 ,
0 2, whose function values there are Z = ge(D). If ge(, ) meets the four invariance
conditions1 within D (Eqs. A.11A.14), then we can extend this function to all - points
(R2 ) by applying two fundamental operations and , which can be seen upon examination
of Figure A.2:
(, ) = ( , + ) (, + ) = (, + 2)
(A.15)
(, ) = ( + 2, ),
where represents a glide reflection and represents a translation. In the -direction, we

have a glide reflection of = , where is the translational part. Two successive applications
of the glide reflection () equates to one simple translation of = 2 in the -direction. In
1
If we are within D, we can speak of the invariance conditions as the boundary conditions.
the -direction, we have a simple translation of = 2. Figures A.3 and A.4 show examples of
extending a function from D to the entire - plane. In Figure A.3 the function does not meet
the invariance conditions, and it is only shown to help illustrate the mapping. In Figure A.4
the function does meet the invariance condition.
Here is a basic algorithm for extending a function ge(, ) to the - plane:
1. Begin with a function defined on D : 0 , 0 2.
2. Ensure that the function meets the boundary conditions on D, which can be determined
from the invariance conditions (Eqs. A.11A.14).
3. Extend the function to the - plane using the two operations and (Eq. A.15).
4. The new function g(, ), defined for the entire - plane, will satisfy the invariance
conditions, and thus will be a well-defined function on the sphere.
If we begin with a function ge(, ) and obtain g(, ) in this manner, we cannot be sure
that the derivatives of g(, ) are continuous at the boundaries between the sub-domains of
the - plane. In order to meet this requirement, we could set boundary conditions on the
derivatives in step 2. However, the simplest thing is to start with a function g(, ) that is
defined in the entire - plane and that meets the invariance conditions.
(2, +2) (, +2) (+2, +2)

2
(, +2) (+, +2) (+3, +2)

South Pole
Greenwich
Meridian
(2, ) (, ) (+2, )
North Pole
2 2 3 4
(, ) (+, ) (+3, )
(2, 2) (, 2) (+2, 2)
B (+, +)
z
: (, ) = (+, ) = (+2, ) = ...

: (, ) = (, "2N)

y

x
Figure A.2: A. Extending the function ge(, ) to the entire - plane (R 2 ). The two basic
operations are given by Eq. A.15. The first operation for each case gives g and g ; repeating
these operations will cover the entire plane. The L-shape is shown to assist in illustrating
the transformations. The -axis can be thought of as the North Pole; the -axis can be
thought of as the Greenwich Meridian. u SP and uNP are the values of u at the poles; they
are constant for a fixed time t. B. Spherical coordinate system (S 2 ), illustrating the glide
reflection (, ) = (, + ) and the translation (, ) = ( 2N , ).
3 3
A B
2 2

(latitude)
(latitude)
0 0
- -
-2 -2
-3 -3
-3 -2 - 0 2 3 -3 -2 - 0 2 3
(longitude) (longitude)
2 2 2
Figure A.3: A. Example function ge(, ) = e ((0 ) +(0 ) )/ , defined in domain D of the
- plane. B. Our new function g(, ), found by extending ge(, ) to the - plane using
the method in Figure A.2A. In this example, ge(, ) does not satisfy the invariance conditions
given by Eqs. A.11A.14: here the North Pole and South Pole are multiple-valued.
3 3
A B
2 2

(latitude)
(latitude)
0 0
- -
-2 -2
-3 -3
-3 -2 - 0 2 3 -3 -2 - 0 2 3
(longitude) (longitude)
2 2
Figure A.4: A. The function ge(, ) = e / , defined in domain D of the - plane. B. Our
new function g(, ), after extending u(, ) to the - plane using the method in Figure A.2A.
Here ge(, ) satisfies the invariance conditions given by Eqs. A.11A.14.
Appendix B
Spherical harmonic functions
Spherical harmonic functions, Ylm (, ), are smooth, well-defined functions on the sphere
(see Appendix A). Spherical harmonic conventions have been presented and discussed in Ed-
monds (1960); Phinney and Burridge (1973); Gilbert and Dziewonski (1975); Dziewonski and
Woodhouse (1983); Dahlen and Tromp (1998), among others. The purpose of this appendix
is to present the standard conventions used in seismology those of Gilbert and Dziewonski
(1975) and then to show how Matlabs built-in functions can be modified to match the
conventions used in seismology. Then we look at expanding functions in terms of spherical
harmonic functions.
B.1 Spherical harmonic conventions
Adopting the convention of Gilbert and Dziewonski (1975), we define complex surface spherical
harmonics of degree l and order m as
Ylm (, ) = Xlm () eim (B.1)
= Xlm () cos m + i Xlm () sin m , (B.2)
where Xlm () are the generalized Legendre functions, defined as

h i1/2

2l+1 1/2 (lm)!
(1)m

4 (l+m)! Plm (cos ) if 0 m l
Xlm () = (B.3)

X
l(m) () if l m 0 ,
144
APPENDIX B. Spherical harmonic functions 145
or, equivalently,
1/2 1/2
|m| 2l + 1 (l |m|)!
Xlm () = (1) Pl|m| (cos ) (l m l) . (B.4)
4 (l + |m|)!
where we have used the definition of Dziewonski and Woodhouse (1983) to explicitly show
the case for m < 0. For the special case of the sphericallysymmetric function with l=0 and
m=0, we have X00 () = (4)1/2 .
The associated Legendre functions of degree l and order m are defined as
l+m
1 d
Plm () = l (1 2 )m/2 (2 1)l (m 0) (B.5)
2 l! d
where = cos . (With = cos , we always meet the required condition 1 1.) The
Legendre polynomial of degree l is given by the associated Legendre function of order m = 0:
l
1 d
Pl () = l (2 1)l . (B.6)
2 l! d
The associated Legendre functions (Eq. B.5) of positive order, m > 0, can be written in terms
of the Legendre polynomials Pl ():
m
2 m/2 d
Plm () = (1 ) Pl () (m > 0) (B.7)
d
Dahlen and Tromp (1998) define real surface spherical harmonics by

2Xl|m| () cos m if l m < 0

Ylm (, ) = Xl0 () if m = 0 (B.8)

2X () sin m if 0 < m l .

lm
L=8, M=0 L=8, M=1 L=8, M=2
L=8, M=6 L=8, M=7 L=8, M=8
Figure B.1: Sample functions of Ylm (, ). The functions are plots of Eq. B.14 for l = 8 with
varying m. The functions are plotted on the order q=3 triangular grids.
B.2 Spherical harmonic conventions in Matlab
Matlab has the built-in function legendre. Two types of commands are possible. First,
Plm() = legendre(l, ) are the associated Legendre functions of degree l and order m =
0, 1, . . . , l, defined as
m
m 2 m/2 d
Plm() = (1) (1 ) Pl () (m 0) (B.9)
d
= (1)m Plm () (m 0) , (B.10)
where Plm () is given by Equation (B.7). The input = cos is a vector of N values, and the
output is a (l+1)-by-N matrix, with the mth row corresponding to the associated Legendre
function of order m-1 (and degree l). For example say l = 5 and the input vector is 100
elements of = [0, ] and = cos . Then the Matlab command legendre(l, ) gives a
6-by-100 matrix of the associated Legendre functions of degree l = 5, orders m = 0 (first row)
to m = 5. Second, Xlm() = legendre(l, , 0 sch0 ) is the Schmidt semi-normalized associated
Legendre functions, defined as

Plm(cos ) if m = 0
Xlm() = (B.11)

h i1/2
(1)m 2 (lm)! Plm(cos ) if m > 0 .
(l+m)!
Substituting Equation (B.10) (with = cos ) gives

Plm (cos ) if m = 0
Xlm() = (B.12)

h i1/2
2 (lm)! Plm (cos ) if m > 0 ,
(l+m)!
where Plm ()is given by Equation (B.7). In order to adopt the convention of Dziewonski
and Woodhouse (1983), we have to alter the Matlab function, changing Equations (B.12) and
(B.4) for m 0. We alter the Matlab output for X lm () by (i) multiplying all values by
p
A = (2l + 1)/(4), (ii) changing the sign of every other order (multiply by (1) m ), and

(iii) dividing all non-(m=0) values by 2.

A Xlm() if m = 0
Xlm () = (B.13)

(1)m A2 Xlm() if m 0 .
Eq. B.13 gives the expression for converting the Matlab function Xlm() = legendre(l, , 0 sch0 )
into the Xlm () functions used by Dziewonski and Woodhouse (1983). We use real surface
spherical harmonics Ylm(, ) to use as test functions for testing the accuracy of our numerical
approximation for the Laplacian operator (Section 4.5). We define these functions as

Xl0 () if m = 0
Ylm(, ) = (B.14)

Xlm () sin m if 0 < m l ,
where Xlm () is given by Equation (B.4). Note that for each m there exists a cos m Y lm
function and a sin m Ylm function (Eq. B.8). We chose to use the test functions only for
the case with m > 0. Plots of Equation (B.14) for l = 8 are shown in Figure B.1.
B.3 Spherical harmonic expansions
Functions on the sphere can be expanded in terms of the spherical harmonic functions,
Ylm (, ). Phase velocity maps are generally expressed as a set of spherical harmonic co-
efficients. The spherical harmonic expansion for a phase velocity map c(, ) 1 is written
as
lmax
X l
X
c(, ) = (Alm cos m + Blm sin m) Xlm () , (B.15)
l=0 m=0
where Alm and Blm are the spherical harmonic coefficients up to maximum degree l = lmax,
and Xlm () is given by Equation (B.4). In the ray-tracing calculations it is necessary to know
the derivatives of c(, ):
lmax
X l
X
dc(, )
d c = (m Alm sin m + m Blm cos m) Xlm () (B.16)
d
l=0 m=0
lmax
X X l
dc(, )
d c = (Alm cos m + Blm sin m) d [Xlm ()] (B.17)
d
l=0 m=0
lmax
X X l
d2 c(, )
d c = (m Alm sin m + m Blm cos m) d [Xlm ()] (B.18)
dd
l=0 m=0
lmax
X X l
d2 c(, )
d c = (Alm cos m + Blm sin m) d [Xlm ()] . (B.19)
d 2
l=0 m=0
Using Equation (B.4) with m 0 and
1/2 1/2
m 2l + 1 (l m)!
klm = (1) , (B.20)
4 (l + m)!
we find:
Xlm () = klm Plm (cos ) (B.21)
d [Xlm ()] = klm d [Plm (cos )] (B.22)
d [Xlm ()] = klm d [Plm (cos )] . (B.23)
What is needed are the derivatives d Plm (cos ) and d Plm (cos ). The first derivatives
1
Typically the maps are written as c(, )/c0 , not c(, ).
are calculated using the recursion relation (Arfken and Weber , 1995, Eq. 12.87)
1/2 d
2 1 x2 [Plm (x)] = Pl(m+1) (x) (l + m)(l m + 1)Pl(m1) (x) , (B.24)
dx
where x = cos . (Eq. B.24 shows why the ray-tracing calculations do not work exactly at
the poles, where = 0 and = , and thus 1 x 2 = 0.) Equation (B.24) can be written in
matrix form:
T
1 Pl(m+1) (x)
1
dx Plm (x) = 2 1/2 . (B.25)
2(1 x )
(l + m)(l m + 1) Pl(m1) (x)
Making the substitution x = cos , we get:

dPlm (cos ) dPlm (cos ) d 1
dx Plm (x) = = = d Plm (cos ) (B.26)
d(cos ) d d cos sin
d Plm (cos ) = sin dx Plm (x) . (B.27)
Using Equation (B.25) with x = cos , we can write Equation (B.27) as:
T
1 Pl(m+1) (cos )
1
d Plm (cos ) = . (B.28)
2
(l + m)(l m + 1) Pl(m1) (cos )
When m = 0 for l 1 in Equation (B.25), we have to use the relationship (Arfken and Weber ,
1995, Eq. 12.81a ):
(l m)!
Pl(m) (x) = (1)m Plm (x) . (B.29)
(l + m)!
When m = 0, Equations (B.25), (B.29), and (B.28) give:
1 Pl(1) (x)
dx Pl0 (x) = 2 1/2
P l(1) (x) + P l(1) (x) = (B.30)
2(1 x ) 2(1 x2 )1/2
1

d Pl0 (cos ) = 2 Pl(m+1) (cos ) + Pl(1) (cos ) = Pl(1) (cos ). (B.31)
The associated Legendre equation is given by (Arfken and Weber , 1995, Eqs. 12.71 and 12.72):

1 d dPlm (cos ) m2
sin + l(l + 1) Plm (cos ) = 0 (B.32)
sin d d sin2

2
2 d d2 m2
1x Plm (x) 2x 2 Plm (x) + l(l + 1) = 0. (B.33)
dx2 dx 1 x2
Equation (B.32) can be solved for the second derivatives d Plm (cos ):

1 m2
[cos d Plm (cos ) + sin d Plm (cos )] + l(l + 1) Plm (cos ) = 0
sin sin2

cos m2
d Plm (cos ) = d Plm (cos ) l(l + 1) Plm (cos ) . (B.34)
sin sin2
Equations (B.28) and (B.34) are used to calculate the first and second derivatives of the
Legendre functions, which are used in the X lm functions and derivatives (Equations (B.21)
(B.23)), which are in turn used in the ray-tracing calculations.
An important property of the spherical harmonic functions is that the Y lm (, ) functions
are eigenfunctions of the Laplacian operator:

2 l(l + 1)
Ylm = Ylm , (B.35)
a2
where the Ylm (, ) functions are plotted on a sphere with radius a. This was the reason for
choosing the Ylm (, ) functions as our test functions for our numerical approximation for the
Laplacian (Section 4.5).
Appendix C
Euler angles and rotations
Rotations in three dimensions were used throughout this study. Any rotation of a rigid
body (e.g., a sphere) may be represented by three real numbers. The Euler angles (, , )
represent a useful convention for describing a rotation in 3D (Edmonds, 1960). Using the
Euler-angle convention, each 3D rotation consists of three separate rotations. Each separate
rotation is pictured in Figure C.1 and has an associated matrix, shown in Equation (C.1).
The Euler rotation D(, , ) on the point (x, y, z) consists of three successive rotations
and is given by:

cos sin cos 0 sin
0 cos sin 0
x

sin cos 0 0 1 0
sin cos 0 y . (C.1)

0 0 1 sin 0 cos 0 0 1 z
The reverse of the operation D(, , ) is given by
D1 (, , ) = D( , , ) . (C.2)
Figure C.1 illustrates the forward rotations, and the reverse operation is shown in Figure C.2.
They key concept is that the Euler angles specify rotation angles for the axes, while the
rotation matrix D is applied to the points (Edmonds, 1960). Thus, D is a 3-by-3 matrix whose
entries are determined by (, , ) via Equation (C.1).
151
APPENDIX C. Euler angles and rotations 152
Here we will illustrate an example of how we use the Euler rotations. Suppose we wish to
compare the phase anomalies and amplitude anomalies calculated using the numerical model
(i.e., spherical grids) with those calculated using ray tracing. Here are the four steps in
selecting the receivers.
1. Specify a target source (lat, lon) and a set of target receivers (lats, lons). The locations
are arbitrary if we were comparing our results (numerical or ray theory) with actual
data, then we would use the actual source and receiver locations.
2. Determine the Euler angles (, , ) that rotate the target source to the North Pole,
which is the position of the source in the reference frame of the numerical calculations.
Rotate the target source and receivers by D(, , ).
3. Find the nearest gridpoint to each target receiver. These are your receiver points (lats,
lons) in the reference frame for the numerical calculations. Carry out the numerical
calculations on the spherical grids and obtain the hex phase and amplitude anomalies.
4. Rotate these receivers (and the source) back into the original reference frame, by using
the reverse of the initial rotation, D 1 (, , ). The source ends up back at the target
source, and the receivers end up close to the target receivers. These are the actual
lat-lon points of the receivers that you will use in the ray-tracing calculations.
For the ray-tracing calculations, the Euler angles are specified such that the source is
on the equator at (0, 0) and the receiver is at (0, ). Thus, each source-receiver pair has
a different set of Euler angles to rotate them into the ray-tracing reference frame. Once a
source-receiver pair is rotated to the equator, then we calculate the phase and amplitude
anomalies associated with the ray path.
Initial points: Now rotate x axis by Now rotate z' axis about y' axis by
= -15o = -30o
z z' z''
3 3
8 8
1 7 1 7
9 9 5
5
6 6
13
12 13 12
10 11 11
10
2
z 0 15
2 y z 0 15 y'
4 14 4
1 14 1
x' x x'
-1 -1
-1 0 0
x -1 x
-1 1 -1 1
0 1 0 1
y y
Now rotate x'' axis about z'' axis by Final points

= 60o
z'' z'''
1 1 3
3 8 8
5
7 5
7
9 13 9
13
12 6
6 2
12
z 0 15
y'' z 0 2
y'''
11 11
x''' 10 15
10
x'' 4 x'''
4
14 1 14
1
-1 -1
0 0
-1 x -1 x
-1 1 -1 1
0 1 0 1
y y
Figure C.1: Euler angle rotation. The set of sample points consists of 15 points in the
quadrant x, y, z > 0. Each Euler angle (, , ) specifies a rotation of the coordinate axes.
The corresponding rotation of the sample points is given by one of the matrices in Eq. C.1.
The -rotation is a rotation in the xyplane. The -rotation is a rotation of the in the x 0 z 0
plane. The -rotation is a rotation of the in the x 00 z 00 plane. The final reference frame S 000 is
the same as the initial reference frame S.
Initial points: Now rotate x axis by Now rotate z' axis about y' axis by
= 120o = -30o
z
z' z''
8 3
1 3 1
8 5 7
9
7 5
9 13
6
13 12
6
11
12 x' 2
10
2 15
z 0 11 y z 0 y'
10 4
15 14
x 1
4 x'
1 14
-1 -1
0 0
-1 x -1 x
-1 1 -1 1
0 0 1
1
y y
Now rotate x'' axis about z'' axis by Final points

= 195o
z'' z'''
73 3
1 98 1 7
8
6
5 9
5
11 10 6
12
13
x''' 12 13
10 11
14 1
4
15
2
z 0 z 0 2 y'''
y'' 4
15
1 14
x'' x'''
-1 -1
0 0
-1 x -1 x
-1 1 -1 1
0 1 0 1
y y
Figure C.2: Euler angle rotation (in reverse). Here we reverse the operations in Figure C.1;
the final points here are the initial points in Figure C.1. The net rotation matrix is given by
Eq. C.2.
Appendix D
Derivations for selected equations
In this appendix we include derivations for several of the equations presented in Chapter 3.
Section D.4 below contains various definitions and properties used in the derivations.
D.1 Equation (3.12)
The method of Greens functions in obtaining Equation (3.12) is that the solution to Equa-
tion (3.10) is the convolution of the source-time function f lm (t) with the Greens function,
i.e., the solution to Equation (3.12) for an impulsive force. Thus, using the definitions in
Section D.4, the solution to Equation (3.12) is given by
Z t
ulm (t) = G(t) c2 flm (t) = c2 G(t ) flm ( ) d , (D.1)

where
sin l t
G(t) = . (D.2)
l
155
APPENDIX D. Derivations for selected equations 156
D.2 Equations (3.27) and (3.33)
The derivation of Equation (3.27) is as follows. Using the identity
sin(A B) = sin A cos B cos A sin B
sin(l t l ) = sin l t cos l cos l t sin l ,
we expand the left side of Equation (3.27):
Z 2 2
!

e /2 sin l t cos l cos l t sin l
d (D.3)
2 l
Z 2 2
! Z 2 2
!

e /2 sin l t e /2 cos l t
= cos l d sin l d (D.4)
2 l 2 l
Z 2 2
! Z 2 2
!

sin l t e /2 cos l t e /2
= cos l d sin l d (D.5)
l 2 l 2
"Z 2 2
! #

sin l t e /2 cos l t
= R ei l d 0 (D.6)
l 2 l
" 2 2
!#
sin l t e /2
= R F (D.21) (D.7)
l 2
sin l t h 2 2 /2 i
= R e l (D.23) (D.8)
l
sin l t 2 2 /2
= e l , (D.9)
l
R
where in obtaining Equation (D.6) we have used Heven ( ) Hodd ( ) d = 0.
The derivation of Equation (3.33) is similar to that above, where instead of Equation (D.3)
we have
Z 2 2
!

e /2 sin l t cos l cos l t sin l
d (D.10)
2 2 l
Z 2 2
!

e /2 sin l t
= 2
cos l d
2 l t
Z 2 2
!
e /2 cos l t
2
sin l d (D.11)
2 l t
Z 2 2
!

sin l t e /2
= cos l d
l 2 2
Z 2 2
!
cos l t e /2
+ ( sin l ) d (D.12)
l 2 2
Z
sin l t cos l t il
= 0 + I h2 ( ) e d (D.13)
l l

cos l t d
= I F (h1 ( )) (D.21), (3.31) (D.14)
l d
cos l t h 2 2
i
= I i l el /2 (D.24) (D.15)
l
cos l t 2 2
= l el /2 (D.16)
l
2 2 /2
= cos l t el , (D.17)
R
where in obtaining Equation (D.13) we have used Hodd ( ) Heven ( ) Heven ( ) d = 0.
D.3 Equations (3.20) and (3.62)
For Equation (3.20), we want to simplify the following expression:

Z 2 2
! Z 2

e /2 ew/2 dw
d =
0 2 0 2 2
Z
1 2
= ew/2 dw
22 0

2
= ew/2 = (0 1) = 1 , (D.18)
0
where we have made the substitution w 2 . This also shows how we selected the normal-
izations on g().
For Equation (3.62), we begin with
Z 2 2
e /2
Il (0) = lim Pl (cos ) sin d . (D.19)
0 0 2
In the limit as 0, the function in the integral becomes a delta function () centered at
= 0. Thus, as 0, 0 (a spike at the North Pole), and we can use the first-order
approximation sin . Furthermore, cos 1 and P l (1) = 1. Putting this together we
have
Z 2 2 Z 2 2
e /2 e /2
Il (0) = Pl (1) d Pl (1) d = 1 . (D.20)
0 2 0 2
The upper integration limit of can be replaced by , because g() 0 for [, ]. The
last step follows from the results in Equation (D.18). (This result can also be seen graphically
in Figure 3.3.)
D.4 Some definitions and properties
Here we list some equations used in the Chapter 3. Spherical harmonic functions are discussed
in Appendix B.
Fourier transforms
We define forward and inverse Fourier transforms:

Z
F [h(t)] = e
h() = h(t) eit dt (D.21)

h i Z
1 e 1 e
F h() = h(t) = h() eit d . (D.22)
2
Some properties of Fourier transforms that we will use are as follows:

h i Z
2
a2 t2 a2 t2
F ke = ke dt e 4a2 (D.23)

0
F (f (t)) = i F(f (t)) (D.24)
For example, using Equation (D.23), the Fourier transform of the Gaussian function h 1 (t) in

Equation (3.17) is, with k = 1/( 2 ) and a = 1/( 2 ) :
2 2 /2
F [h1 (t)] = e . (D.25)
Convolution
The convolution of f (t) with g(t) is given by the formula 1
Z
f (t) g(t) = f (t ) g( ) d . (D.26)

Properties of convolutions that we will use are as follows:
f (t) g(t) = g(t) f (t) ; (D.27)

1
Arfken and Weber (1995, p. 863) have slightly different definition of convolution.
Lets say we have u = u(x, t) and h = h(t). Then the following is true:
Z
Fa (t, x) = u h = u(t , x) h( ) d (D.28)
Z

Fb (t, x) = h u = h(t ) u(, x) d (D.29)

Z
Fa
= (u h) = u(t
, x) h( ) d = u h (D.30)
t t
Z

Fb ) u(, x) d = h u
= (h u) = h(t (D.31)
t t
Greens functions
The solution x(t) to the second-order ODE
d2 x dx
2
+ k + 02 x = h(t) (D.32)
dt dt
is given by
Z
x(t) = G(t) h(t) = G(t ) h( ) d , (D.33)

where G(t) is Greens function and satisfies Equation (D.32) for an impulsive source, i.e.,
d2 G dG
2
+ k + 02 x = (t) . (D.34)
dt dt
By definition, G(t) = 0 for t < 0, which implies that G(t ) = 0 for > t. We can split
x(t) into two integrals and simplify:
Z t Z
x(t) = G(t ) h( ) d + G(t ) h( ) d
t
Z t
= G(t ) h( ) d . (D.35)

Thus, this expression is true for t > , where is the time value within Equation (D.32).
This places the condition that the function x(t) and its derivatives must vanish at (sufficiently
large) negative times.
Appendix E
Ray-tracing equations
E.1 Receiver at a fixed longitude
The equations for path, amplitude, and phase anomaly in Woodhouse and Wong (1986) are
specified for a source and receiver situated on the equator. Here we show a slightly more
generalized version of the equations, where only the source is on the equator and the receiver
is at a particular longitude f but not necessarily on the equator. Thus, the source is on the
Greenwich meridian at (0 , 0 ) = ( 2 , 0), and the receiver is at (f , f ). Figure E.1A shows
the setup for the ray-tracing.
Variables and notation
Variable definitions follow from Woodhouse and Wong (1986):
= cot
d d cot d d
= = = csc2
d d d d

0 =
0
0
0 = = .
0
161
APPENDIX E. Ray-tracing equations 162
Referring to Figure E.1, we define the azimuth along the ray, , as follows:
h1 d 1 d
tan = = , (E.1)
h2 d sin d
or, in terms of ,

1 d d 1 1
tan = = () = sin . (E.2)
sin d d sin csc2
Note that the azimuth is measured in a direction counter-clockwise from eastward at a given
point on the ray. The initial azimuth is the take-off angle,
0 , (E.3)
and the final azimuth is the off-azimuth arrival angle, f , which can be expressed as
tan f = f sin f , (E.4)
where f is the finishing latitude of the ray at the finishing longitude f .

Initial conditions are (see Figure E.1):

0 = (E.5)
2
0 = 0 (E.6)

1 d d
tan = tan 0 = = . (E.7)
sin 0 d =0 d =0
These lead to the following expressions:
0 = cot 0 = 0 (E.8)

2 d
0 = csc 0 = tan . (E.9)
d =0
In the remainder of this section, we develop equations for the homogeneous ray paths. Then
we derive an expression for the amplitude anomaly due to a heterogeneous earth, for a receiver
at a fixed longitude f .
Homogeneous ray paths
The solutions to Equation (6.3) have the form
() = A cos + B sin . (E.10)
Differentiating Equation (E.10) we have, in two different forms,
d
= A sin + B cos
d
= A sin B cos .
Our initial conditions specify the initial position of the source and the take-off angle :
(0) = 0 (E.11)

d
0 = tan (E.12)
d =0
Using the initial conditions, we have
(0) = A = 0 (E.13)

d
= B = 0 = tan . (E.14)
d =0
Our solution can be expressed in several different forms:
() = 0 sin (E.15)
() = 0 cos (E.16)
cot = tan sin (E.17)
() = cot1 (tan sin ) (E.18)

() = tan1 (tan sin ) , (E.19)
where is colatitude and is latitude ( =

2 ). Plots of Equation (E.19) are shown in
Figure 6.2.
Amplitude anomaly
The amplitude anomaly is defined as the square root of the ratio of geometrical spreading in
the homogeneous case to the geometrical spreading in the heterogeneous case:

d 1/2
d hom

A = d
. (E.20)
d het

Figure E.1B shows the starting point for the derivation. We consider two rays emanating
from the source with an infinitesimal take-off angle difference, . The two rays form the ray
tube, whose width is measured by d.
Homogeneous Case. With the homogeneous earth, the ray paths are great-circle paths,
and the geometrical spreading on the sphere in the homogeneous case is simply given by

d
= a sin , (E.21)
d hom
where is the arc distance between the source and the receiver and a is the radius of the
earth (Woodhouse and Wong, 1986, Eq. 40b). Physically, this means that the geometrical
spreading is a minimum at the source ( = 0) and the antipode ( = ), where all the

great-circle paths converge, and it is a maximum at a distance = 2 from the source (and
antipode), where the rays are most spread out.
In our case, the source is at (0 , 0 ) and the receiver is at (f , f ). The dot-product gives
us an expression for the arc distance, which we show here for a source at ( 2 , 0):
cos = rs rr
= sin 0 cos 0 sin f cos f + sin 0 sin 0 sin f sin f + cos 0 cos f
= sin f cos f . (E.22)
Equation (E.21) then can be written as

d 1/2 1/2
= a 1 cos2 = a 1 sin2 f cos2 f , (E.23)
d hom

If our receiver is on the equator (f = 2 ), then f = , and Equation (E.23) reduces to
Equation (E.21).
Heterogeneous Case. The geometrical spreading in the heterogeneous case can be

expanded using the relationships above (Eq. E.9):

d d d tan d
= = sec2 . (E.24)
d het d tan d d0
From Figure E.1B we see that we have the relationship cos = d/d, or d = d cos .
Differentiating = cot , we have d = csc 2 d, or d = d sin2 . These two expressions
can be used in Equation (E.24):

d d
= sec2
d het d0
d cos
= sec2
d0
d sin2 cos
= sec2
d0
d
= cos sec2 sin2 . (E.25)
d0
Using trigonometric substitutions, we have
!
d d 1
= p 1 + tan2 sin2 f , (E.26)
d het d0 1 + tan2 f
where we have entered the f -subscripts to indicate that we are at the finishing longitude f .
Substituting Equations (E.4) and (E.9) into Equation (E.26), we get

d 1/2
= a f0 1 + f2 sin2 f 1 + 02 sin2 f , (E.27)
d het
where we have also substituted 0 d/d0 and inserted the earth radius, a, which is simply a
scaling factor for the geometric spreading. Note that Equation (E.27) reduces to Equation 36

of Woodhouse and Wong (1986) when the receiver is placed on the equator ( f = 2 and
f = ).
We now substitute Equations (E.27) and (E.23) into Equation (E.20) to get an expression
for the amplitude anomaly:

d 1/2
d hom

A = d
d het

1/2

a 1 sin2 f cos2 f
1/2

= 1/2

a sin2 f 1 + 02 1 + f2 sin2 f f0
1/2
1 1/4 1/4
2 2 2 2
= 2 1 sin f cos f 1 + f sin f
sin f 1 + 02 f0
1/4 1/4
1 sin2 f cos2 f 1 + f2 sin2 f
= 1/2 1/2 . (E.28)
0
|sin f | f 1 + 02
Note the inclusion of the factor f in three places, which results from allowing the receiver
to be situated at non-zero latitudes. For the case when the receiver is placed on the equator

(f = 2 and f = ), Equation (E.28) reduces to

sin 1/2 1/4

A = 0 2
1 + 2 () , (E.29)
f 1 + 0
which is exactly Equation 41 of Woodhouse and Wong (1986). Plots of Equation (E.28) are
shown in Figure 6.10, which showing the variation of amplitude anomaly as a function of
take-off angle along a finishing longitude f .
E.2 Equation (6.19) ds/d
Beginning with Eq. 12 of Woodhouse and Wong (1986) (Dahlen and Tromp, 1998, Eq. A.100),
we have the expression for the infinitesimal distance ds between a point (x 1 , x2 ) and (x1 + dx1 ,
x2 + dx2 ):
ds = (g x x )1/2 dt
= (g dx dx )1/2 , (E.30)
where x1 = cot , x2 = , and g are the covariant components of the metric tensor,
given by
a2
g11 = g12 = 0
(1 + 2 )2
a2
g21 = 0 g22 = , (E.31)
(1 + 2 )
and the contravariant components (Woodhouse and Wong, 1986, Eq. 32) are found via
[g ] = [g ]1 :
2
g 11 = a2 1 + 2 g 12 = 0

g 21 = 0 g 22 = a2 1 + 2 . (E.32)
(Note that with = cot , we have 1/(1 + 2 ) = sin2 .) Substituting Equation (E.31) into
Equation (E.30), we have
ds = [g dx dx ]1/2
1/2
= g11 (dx1 )2 + g22 (dx2 )2
1/2
= g11 d 2 + g22 d2
1/2
= g11 2 + g22 d
1/2
a2 2 a2
= + d (E.33)
(1 + 2 )2 (1 + 2 )
1/2
ds 2 1
= a + , (E.34)
d (1 + 2 )2 (1 + 2 )
where we have used x1 = , x2 = , and = d/d. From the latter follows d 2 = 2 d2 .

E.3 Equation (6.11)
We use Newtons method to help find the take-off angles of the ray(s) that converges on a
particular receiver for the ray-shooting. We expand 0 (n ) as a first-order Taylor series, where
n = () is the finishing value for the next iteration and f is the finishing value for the
current iteration:
0
0 (n ) = 0 (f ) + (n f ) (f ) +

0
0 (f ) + (0 f ) (f )

f
= 0 (f )
0 (f )
f
= 0 (f ) 0 , (E.35)
f
where in the second step we set n = 0 and in the last step we substituted Equation (6.6).
Thus, Equation (E.35) allows a method of picking the 0 s so that the rays converge on the
receiver. See Section 6.3.
R
y
S
x
ray
wa
wa
vef
vefr
r
o
on
nt
Figure E.1: Variables and setup for the ray-tracing equations on the sphere (Section E.1). The source
is situated on the equator at (0 , 0 ) = ( 2 , 0), and the receiver is at (f , f ). The take-off angle is
0 . A. The homogeneous ray path is a great-circle ray path (dashed, ); the heterogeneous ray
path (bold, het ) deviates from the great-circle path (het > ). B. Exaggerated view of the ray in
the vicinity of the receiver, showing the relationship cos = d/d, where is the azimuth of the ray.
The geometrical spreading is defined as d/d. The quote-marks highlight the fact that the earth
model is rotated into the ray-tracing reference frame (Appendix C), so that the actual North Pole,
etc., are not where we have labeled them (except in the case where the earth model is unrotated).
C. Coordinate system in the vicinity of the receiver.
Appendix F
Wave parameters and dispersion
Table F.1: Wave parameters and some related equations. The equations in the last column
incorporate l and a and are approximate for l 1 (see Henson and Dahlen, 1986; Woodhouse,
1996, Eq. 2.110).
earth radius a 6371 km
earth circumference 2a 40030 km
harmonic degree l
1 2 2a
period T T = f = = c T = c(l+ 12 )
1 c c(l+ 21 )
frequency f f= 2 = T = f= 2a
c 2 2a
wavelength = f = cT = k = l+ 21
2 2f l+ 12
wavenumber k k= c = = c k= a
2 c(l+ 21 )
angular frequency = T = ck = a
2a
velocity c c= T = f = k c= T (l+ 21 )
170
APPENDIX F. Wave parameters and dispersion 171
10.0 10.0
A B
9.0 9.0 Love
Group Velocity (km/s)

Phase Velocity (km/s)
8.0 8.0
7.0 L=3 Love 7.0

L=2
Rayleigh
6.0 6.0
L=3
L=2
5.0 5.0
Rayleigh
4.0 4.0
3.0 3.0
0 1000 2000 3000 0 1000 2000 3000
Period, T (s) Period, T (s)
5.5 5.5
C D
Group Velocity (km/s)
Phase Velocity (km/s)
5.0 5.0
Love
4.5 4.5 Love
Rayleigh
4.0 4.0
Rayleigh
3.5 3.5
0 100 200 300 0 100 200 300
Period, T (s) Period, T (s)
Figure F.1: Dispersion of surface waves, plotted as phase velocity c(T ) and group velocity
U (T ). The curves are calculated using normal-modes summation with the 1D earth model
PREM (Dziewonski and Anderson, 1981); each point corresponds to a normal mode with
harmonic degree l. A couple points are labeled in A. These curves are for the fundamen-
tal spheroidal (Rayleigh) and toroidal (Love) modes (n = 0); dispersion curves of overtone
branches (n 1) can be seen in Figures 11.5 and 11.6 of Dahlen and Tromp (1998). A. Phase
velocity dispersion c(T ). Gray box shows region of C. B. Group velocity dispersion U (T ).
Gray box shows region of D. C-D. Same as A-B, only with different axes ranges. With pe-
riods T < 300s, Love waves travel faster than Rayleigh waves. Calculations courtesy of Ana
Ferreira.
Appendix G
Flow chart of computer programs
Figure 7.1 shows a very schematic flow chart of the method to generate scatterplot com-
parisons among the three different methods: numerical model, exact ray theory, linearized
ray theory. Figure G.1 shows a more accurate depiction of the actual method, including the
naming convention for the data files.
172
APPENDIX G. Flow chart of computer programs 173
Ray-tracing (Fortran) Numerical model (Matlab)
Lats-Lons of Order of 'Earthquake' Source

Phase velocity map: source (0,0) and spherical grid: parameters:
lmax, T, 25 receivers q ,
focusdzeros.f trigrid.m
subroutine foctrc spherewaver.m
swanalytic.m
Exact ray tracing program. swnumeric.m
Shoots out 1000 rays and finds the
converging rays. Output the take-off Numerical and analytical solutions to
angles of the converging rays. waves on the spherical membrane
for both the heterogeneous and
homogeneous case. Output three
sets of 25 seismograms.
ztkoff95T040L20e10.dat
focusd.f
subroutine foctrc
subroutine foclin
Ray tracing program.

Calculate the amplitude anomalies srcCEN25recTW95T040L20e10.mat
and phase anomalies for the
exact ray theory and for the
linearized ray theory.
mA95T040L20e10.dat
mB95T040L20e10.dat
rayanoms.m hexanoms.m
seisHEX.m
Read data files and organize them
genfit.m/theoryHEX.m
in a manner suitable for
comparison with the results from
Calculate the amplitude and phase
the numerical model.
anomalies from the seismograms
from the numerical model.
ray2003TW95T040L20e10.dat hexPATW95T040L20e10q6th05.mat
rayhexanoms.m
Creates the scatter plots
Figure G.1: Flowchart of computer programs used in calculating the phase and amplitude
anomalies from ray theory (left) and from the numerical model (right). See Figure 7.1 for a
more conceptual view of the process.
Bibliography
Albarede, F., and R. D. van der Hilst, Zoned mantle convection, Phil. Trans. R. Soc. A., 360,
25692592, 2002.
Arfken, G. B., and H. J. Weber, Mathematical Methods for Physicists, 4 ed., Academic Press,
San Diego, Calif., 1995.
Augenbaum, J. M., and C. S. Peskin, On the construction of the Voronoi mesh on a sphere,
J. Comp. Phys., 14, 177192, 1985.
Baig, A. M., F. A. Dahlen, and S.-H. Hung, Traveltimes of waves in three-dimensional random
media, Geophys. J. Int., 153, 467482, 2003.
Bassin, C., G. Laske, and G. Masters, The current limits of resolution for surface wave
tomography in North America, in EOS Trans. Am. Geophys. Union, vol. F897, p. 81, 2000.
Baumgardner, J. R., and P. O. Frederickson, Icosahedral discretization of the two-sphere,

SIAM J. Num. Anal., 22, 11071115, 1985.
Boschi, L., and G. Ekstrom, New images of the Earths upper mantle from measure-
ments of surface-wave phase velocity anomalies, J. Geophys. Res., 107 (B4), 2059,
10.1029/2000JB000,059, 2002.
Capdeville, Y., E. Chaljub, J. P. Vilotte, and J. P. Montagner, Coupling the spectral ele-
ment method with a modal solution for elastic wave propagation in global earth models,
Geophys. J. Int., 152, 3467, 2003.
Clayton, R., and B. Engquist, Absorbing boundary conditions for acoustic and elastic wave
equations, Bull. Seismol. Soc. Am., 67, 15291540, 1977.
Cui, J., and W. Freeden, Equidistribution on the sphere, SIAM J. Sci. Comp., 18 (2), 595609,
1997.
Cullen, M. J. P., Integration of the primitive equations on a sphere using the finite element
method, Q. J. R. Meteorol. Soc., 100, 555562, 1974.
Dahlen, F. A., and J. Tromp, Theoretical Global Seismology, Princeton U. Press, New Jersey,
1998.
Dziewonski, A., and D. Anderson, Preliminary reference earth model, Phys. Earth Planet. In-
ter., 25, 297356, 1981.
174
BIBLIOGRAPHY 175
Dziewonski, A., and J. H. Woodhouse, Studies of the seismic source using normal-mode theory,
in Earthquakes: observation, theory and interpretation: notes from the International School
of Physics Enrico Fermi (1982: Varenna, Italy), vol. LXXXV, edited by H. Kanamori
and E. Boschi, pp. 45137, North-Holland Pub. Co., Amsterdam, 1983.
Edmonds, A. R., Angular Momentum in Quantum Mechanics, 2 ed., Princeton U. Press,

Princeton, New Jersey, 1960.
Ekstrom, G., J. Tromp, and E. W. F. Larson, Measurements and global models of surface
wave propagation, J. Geophys. Res., 102, 81378157, 1997.
Gilbert, F., and A. M. Dziewonski, An application of normal mode theory to the re-
trieval of structural parameters and source mechanisms from seismic spectra, Phil. Trans.
R. Soc. Lond., 278, 187269, 1975.
Grand, S. P., Mantle shear-wave tomography and the fate of subducted slabs, Phil. Trans.
R. Soc. A., 360, 25452567, 2002.
Heikes, R., and D. A. Randall, Numerical integration of the shallow-water equations on a

twisted icosahedral grid. Part I: Basic design and results of tests, Monthly Weather Review,
123, 18621880, 1995a.
Heikes, R., and D. A. Randall, Numerical integration of the shallow-water equations on a

twisted icosahedral grid. Part II: A detailed description of the grid and an analysis of
numerical accuracy, Monthly Weather Review, 123, 18811887, 1995b.
Henson, I. H., and F. A. Dahlen, Asymptotic normal modes of a laterally heterogeneous earth
2. Further results, J. Geophys. Res., 91, 12,46712,481, 1986.
Komatitsch, D., and J. Tromp, Introduction to the spectral element method for three-
dimensional seismic wave propagation, Geophys. J. Int., 139, 806822, 1999.
Komatitsch, D., and J. Tromp, Spectral-element simulations of global seismic wave

propagationI. Validation, Geophys. J. Int., 149, 390412, 2002a.
Komatitsch, D., and J. Tromp, Spectral-element simulations of global seismic wave

propagationII. Three-dimensional models, oceans, rotation and self-gravitation, Geo-
phys. J. Int., 150, 308318, 2002b.
Kreyszig, E., Advanced Engineering Mathematics, Wiley, New York, 1988.
Laske, G., Global observation of off-great-circle propagation of long-period surface waves,

Geophys. J. Int., 123, 245259, 1995.
Laske, G., and G. Masters, Constraints on global phase velocity maps from long-period po-
larization data, J. Geophys. Res., 101, 16,05916,075, 1996.
Laske, G., G. Masters, and W. Zurn, Frequency-dependent polarization measurements of long-

period surface waves and their implications for global phase-velocity maps, Phys. Earth
Planet. Inter., 84, 111137, 1994.
BIBLIOGRAPHY 176
Masuda, Y., and H. Ohnishi, An integration scheme of the primitive equation model with
an icosahedral-hexagonal grid system and its application to the shallow-water equations,
in Short- and Medium-Range Numerical Weather Prediction, edited by T. Matsumo, pp.
317326, Japan Meteorological Society, 1987.
Mooney, W. D., G. Laske, and G. Masters, CRUST-5.1: A global crustal model at 5 5 ,
J. Geophys. Res., 103, 727747, 1998.
Phinney, R. A., and R. Burridge, Representation of the elastic-gravitational excitation of a
spherical earth model by generalized spherical harmonics, Geophys. J. R. Astron. Soc., 34,
451487, 1973.
Rancic, M., R. J. Purser, and F. Mesinger, A global shallow-water model using an expanded
spherical cube: Gnomic versus conformal coordinates, Q. J. R. Meteorol. Soc., 122, 959
982, 1996.
Randall, D. A., T. D. Ringler, R. P. Heikes, P. Jones, and J. Baumgardner, Climate modeling
with spherical geodesic grids, Comp. Sci. Eng., 3, 3241, 2002.
Ringler, T. D., and D. A. Randall, A potential enstrophy and energy conserving numerical
scheme for solution of the shallow-water equations on a geodesic grid, Monthly Weather
Review, 130, 13971410, 2002a.
Ringler, T. D., and D. A. Randall, The ZM grid: An alternative to the Z grid, Monthly
Weather Review, 130, 14111422, 2002b.
Ringler, T. D., R. P. Heikes, and D. A. Randall, Modeling the atmospheric general circulation
using a spherical geodesic grid: A new class of dynamical cores, Monthly Weather Review,
128, 24712489, 2000.
Ritsema, J., H. J. van Heijst, and J. H. Woodhouse, Complex shear velocity structure imaged
beneath Africa and Iceland, Science, 286, 19251928, 1999.
Ronchi, C., R. Iacono, and P. S. Paolucci, The cubed sphere: A new method for the solution
of partial differential equations in spherical geometry, J. Comp. Phys., 124, 93114, 1996.
Sadourny, R., Conservative finite-difference approximations of the primitive equations on
quasi-uniform spherical grids, Monthly Weather Review, 100, 136144, 1972.
Sadourny, R., and P. Morel, A finite-difference approximation of the primitive equations for
a hexagonal grid on a plane, Monthly Weather Review, 100, 439445, 1969.
Sadourny, R., A. Arakawa, and Y. Mintz, Integration of the non-divergent barotropic vorticity
equation with an icosahedral-hexagonal grid for the sphere, Monthly Weather Review, 96,
351356, 1968.
Saff, E. B., and A. B. J. Kuijlaars, Distributing many points on a sphere, Math. Intelligencer,
19 (1), 511, 1997.
Shearer, P. M., Introduction to Seismology, Cambridge U. Press, Cambridge, 1999.
Stuhne, G. R., and W. R. Peltier, Vortex erosion and amalgamation in a new model of large
scale flow on the sphere, J. Comp. Phys., 128, 5881, 1996.
BIBLIOGRAPHY 177
Stuhne, G. R., and W. R. Peltier, New icosahedral grid-point discretizations of the shallow
water equations on the sphere, J. Comp. Phys., 148, 2358, 1999.
Swarztrauber, P. N., D. L. Williamson, and J. B. Drake, The Cartesian method for solving
partial differential equations in spherical geometry, Dynamics Atmos. Oceans, 27, 679706,
1997.
Sword, C. H., J. F. Claerbout, and N. H. Sleep, Finite-element propagation of acoustic waves

on a spherical shell, in Stanford Exploration Project (SEP), vol. 50, pp. 4378, 1986a.
Sword, C. H., J. F. Claerbout, and N. H. Sleep, Modeling global surface waves by a finite-
element method, in EOS Trans. Am. Geophys. Union, vol. 67, pp. 10941095, 1986b.
Tanimoto, T., Modelling curved surface wave paths: membrane surface wave synthetics,
Geophys. J. Int., 102, 89100, 1990.
Thuburn, J., A PV-based shallow-water model on a hexagonal-icosahedral grid, Monthly

Weather Review, 125, 23282347, 1997.
Trampert, J., and R. Snieder, Model estimations biased by truncated expansions: Possible
artifacts in seismic tomography, Science, 271, 12571260, 1996.
Trampert, J., and J. H. Woodhouse, Global phase velocity maps of Love and Rayleigh waves
between 40 and 150 seconds, Geophys. J. Int., 122, 675690, 1995.
Trampert, J., and J. H. Woodhouse, High resolution global phase velocity distribution, Geo-
phys. Res. Lett., 23, 2124, 1996.
Trampert, J., and J. H. Woodhouse, Assessment of global phase velocity models, Geo-
phys. J. Int., 144, 165174, 2001.
Trefethen, L. N., Group velocity in finite difference schemes, SIAM Review, 24, 113136, 1982.
Trefethen, L. N., Finite difference and spectral methods for ordinary and partial differ-
ential equations, available at http://web.comlab.ox.ac.uk/oucl/work/nick.trefethen/ pde-
text.html, 1996.
Trefethen, L. N., Spectral Methods in Matlab, SIAM, Philadelphia, Penn., 2000.
Tromp, J., and F. A. Dahlen, Variational principles for surface wave propagation on a later-
ally heterogeneous EarthI. Time-domain JWKB theory, Geophys. J. Int., 109, 581598,
1992a.
Tromp, J., and F. A. Dahlen, Variational principles for surface wave propagation on a laterally
heterogeneous EarthII. Frequency-domain JWKB theory, Geophys. J. Int., 109, 599619,
1992b.
Um, J., and F. A. Dahlen, Source phase and amplitude anomalies of long-period surface
waves, Geophys. Res. Lett., 19, 15751578, 1992.
Wang, Z., and F. A. Dahlen, JKWB surface-wave seismograms on a laterally heterogeneous

earth, Geophys. J. Int., 119, 381401, 1994.
BIBLIOGRAPHY 178
Wang, Z., and F. A. Dahlen, Spherical-spline parameterization of three-dimensional Earth,

Geophys. Res. Lett., 22, 30993102, 1995a.
Wang, Z., and F. A. Dahlen, Validity of surface-wave ray theory on a laterally heterogeneous
Earth, Geophys. J. Int., 123, 757773, 1995b.
Wang, Z., F. A. Dahlen, and J. Tromp, Surface wave caustics, Geophys. J. Int., 114, 311324,
1993.
Wang, Z., J. Tromp, and G. Ekstrom, Global and regional surface-wave inversion: A spherical-
spline parameterization, Geophys. Res. Lett., 25, 207210, 1998.
Williamson, D. L., Integration of the barotropic vorticity equation on a spherical geodesic

grid, Tellus, 20, 642653, 1968.
Williamson, D. L., Numerical integration of fluid flow over triangular grids, Monthly Weather
Review, 97, 885895, 1969.
Williamson, D. L., Integration of the primitive barotropic model over a spherical geodesic
grid, Monthly Weather Review, 98, 512520, 1970.
Williamson, D. L., J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber, A standard

test for numerical approximations to the shallow water equations in spherical geometry,
J. Comp. Phys., 102, 211224, 1992.
Wong, Y. K., Upper Mantle Heterogeneity from Phase and Amplitude Data of Mantle Waves,
Ph.D. thesis, Harvard University, 1989.
Woodhouse, J. H., Surface waves in a laterally varying layered structure, Geophys. J. R. As-
tron. Soc., 37, 461490, 1974.
Woodhouse, J. H., Long period seismology and the Earths free oscillations, in Seismic Mod-
elling of Earth Structure, edited by E. Boschi, G. Ekstrom, and A. Morelli, pp. 3180,
Editrice Compositori, Bologna, Italy, 1996.
Woodhouse, J. H., and A. M. Dziewonski, Mapping the upper mantle: Three-dimensional

modeling of Earth structure by inversion of seismic waveforms, J. Geophys. Res., 89, 5953
5986, 1984.
Woodhouse, J. H., and Y. K. Wong, Amplitude, phase and path anomalies of mantle waves,
Geophys. J. R. Astron. Soc., 87, 753773, 1986.

Membrane CHT THESIS

Uploaded by

Copyright:

Available Formats

Membrane CHT THESIS

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Membrane CHT THESIS

Uploaded by

Copyright:

Available Formats

Waves on a Spherical Membrane

Carl Harlan Tape

Thesis submitted to the University of Oxford

Department of Earth Sciences

Advisor: Prof. John H. Woodhouse

Waves on a Spherical Membrane

Ray theory predictions

A comparison of three models

2 The two-dimensional wave equation 8

3 Analytical spherical membrane waves 18

4 Spherical hexagonal grids 34

5 Numerical spherical membrane waves 60

5.5 Stability and accuracy in the numerical model . . . . . . . . . . . . . . . . . . 66

6 Ray calculations on the sphere 79

7 Numerical model vs ray theory 113

A Functions on the sphere 138

B Spherical harmonic functions 144

C Euler angles and rotations 151

D Derivations for selected equations 155

E Ray-tracing equations 161

F Wave parameters and dispersion 170

G Flow chart of computer programs 172

1.1 Rayleigh waves and Love waves . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.1 Rectangular and circular membranes . . . . . . . . . . . . . . . . . . . . . . . 9

3.1 Initial-shape function g(, ) and source-time function h(t) . . . . . . . . . . 31

4.1 Discretizations of the sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

5.1 Phase velocity map: Homogeneous earth . . . . . . . . . . . . . . . . . . . . . 62

5.8 Seismograms for a heterogeneous velocity map . . . . . . . . . . . . . . . . . 75

6.1 Source and receivers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

7.1 Flowchart of seismogram generation . . . . . . . . . . . . . . . . . . . . . . . 114

A.1 Function mappings for spherical coordinates . . . . . . . . . . . . . . . . . . . 138

B.1 Sample functions of Ylm (, ) . . . . . . . . . . . . . . . . . . . . . . . . . . . 146

C.1 Euler angle rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153

E.1 Sketches for ray-tracing equations . . . . . . . . . . . . . . . . . . . . . . . . . 169

F.1 Dispersion of surface waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171

G.1 Flowchart of computer programs . . . . . . . . . . . . . . . . . . . . . . . . . 173

4.1 Spherical hexagonal grids used in other studies . . . . . . . . . . . . . . . . . 36

5.1 Properties of phase velocity maps of Trampert and Woodhouse (1995) . . . . 68

6.1 Phase anomaly convention . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

F.1 Wave parameters and some related equations . . . . . . . . . . . . . . . . . . 170

2. Our problem is mathematically a two-dimensional problem, as the term membrane

3. The mathematics of our problem is stripped of the complexity of wave propagation

Global Surface Waves it s eve

lat = -30.00, lon = 0.00 lat = -90.00, lon = 0.00

Time (s) x 104 Time (s) x 104

90o R1, R180 90o o

x 105 First Six Hours of Denali quake: Unfiltered, Vertical Component

The two-dimensional wave equation

2.2 Waves on a rectangular membrane . . . . . . . . . . . . . . . . . . . . . . 9

2.3 Waves on a circular membrane . . . . . . . . . . . . . . . . . . . . . . . . 10

2.4 Numerical solutions of waves on a planar membrane . . . . . . . . . . . 12

2.5 Waves on a spherical membrane . . . . . . . . . . . . . . . . . . . . . . . 14

The wave equation is given by

2.2 Waves on a rectangular membrane

The two-dimensional wave equation for waves in the xy-plane is

u(x, y, 0) = g(x, y) (2.3)

The general solution is the two-dimensional Fourier sine series

2.3 Waves on a circular membrane

u f |h = h1 (t) = uf |h = h 1 (t) (t ) , (3.50)