Scattering of Elastic Waves in Layered Media:
A Boundary Integral—Normal Mode Method
by
R onald T. K essel
B .A .S c., U niv ersity o f W aterloo, 1984 M .S c., U niversity o f W aterloo, 1989
A D issertation Subm itted in Partial Fulfillm ent o f the R equirem ents for the D egree of
Doctor of Philosophy
in the D epartm ent o f Physics W e accept this dissertation as conform ingto the required standard
D r. T rev o r.W . D aw so n , S u p erv iso r (A djunct P rofessor, D ept, o f P h ysics and A stronom y)
Dr. John. W eaver, Supervisor (Dean of Science, Dept, of Physics and Astronomy)
Dr. Robert W . Stewart, Departmental Member (Adjunct Professor, Dept, of Physics and A stronom y)
Dr. Pauline Van Den Driessche, Outside-M enjfej^fDepartment of M athematics)
Dr. Henrik Schmidt^ External Exam iner '(DefjM Jf Ocean Engineering, Massachusetts Institute o f Technology)
© R onald T. K essel, 1996 U niversity o f V ictoria
All rights reserved. This dissertation may not be reproduced in whole or in part, by photocopying or other means without the permission of the author.
ii
Supervisor: Dr. Trevor W. Dawson
Abstract
Elastic waves are used in geoacoustics to identify rem ote objects and events. C om puter models for such applications are being pressed to handle propagation in three dim ensions, and scattering from penetrable objects o f arbitrary shape. One approach approxim ates earth and ocean as horizontally stratified media into which the features o f interest are embedded. Here the Boundary Integral Equation (BIE) method for harmonic elastic wave scattering is applied to layered media with penetrable fluid or solid inclusions. A new indirect method is used to ev alu ate sin g u lar and poorly convergent BIE co e fficien ts, the num erically troublesom e coefficients being inferred from free-field solutions o f the integral equation in the absence o f scattering. The new model is very flexible. M ultiple inclusions that are penetrable or impenetrable, passive or actively vibrating, are permitted. Inclusions may have edges, and may pass through interfaces between layers. A combined equation method is used to over-determ ine the solution in the event o f num erical instability. C om prehensive num erical tests recom m ended for all BIE scattering models are described. C entral to the m odel is a new norm al m ode m odel o f propagation, now packaged u n d er the nam e SAM PLE, an acronym for S eism o-A coustic M ode P rogram for L ayered E n v iro n m e n ts. Designed for the rigorous demands o f the BIE method, SAM PLE com putes the total elastic field (displacem ent vector and stress tensor) for general point forces and force couples (with and without momen*), using an com plete mode series that is valid in both the near field and far. Modes are found using a stable scattering matrix method together with singular value decom position (SVD), in a robust root-finding routine, The search includes all mode types: P-SV and SH; propagating and evanescent; proper and im proper (on any Riem ann sheet); and interface and duct modes. C lose mode pairs (double roots) in m ultiple channels are easily identified and resolved during the search. Little-known properties of modes that were discovered when using SAMPLE are also reported.
j
Examiners: _ _ _______________________________
Dr. Trevor W . D aw son, Supervisor
(Adjunct Prof., Dept, o f Physics and Astronom y)
/ ________________
Dr. John. W eaver, Supervisor Dr. Robert W. Stewart, Departmental Member (D ean o f S cience, Dept, o f Physics and Astronom y) (Adjunct Prof., Dept, o fjd iy sic sjp i^ A stronom y)
Dr. Pauline y a n Den D ricsschc, OutsldcTvIcmbcr ^ D r. Henrilc Schnjidtf^xtecffef^xam iner (Department o f M athematics) (Dept, o f Ocean Engineering,
C o n t e n t s
iii L is t o f F ig u r e s v ii L is t o f T a b le s X A c k n o w le d g m e n ts xi D e d ic a t io n x ii 1 I n t r o d u c t i o n 1 1.1 M otivation ... ... , | 1.1.1 Need for stratified m edia...1.1.2 Need for solid m e d i a ... 1.1.3 Need for three-dim ensional m o d e l s ...
1.2 T w o cardinal difficulties in s c a tte r in g ...
...
6
1.3 B ackground for the B IE m ethod ... 1.3.1 Review of recent d e v e lo p m e n ts ...
...
7
1.3.2 B IE m eth o d in perspective ... 1.3.3
A new indirect B IE m e t h o d ...
1.3.4A
flexible com bined B IE m e t h o d ... ...II
1.1A com plete G reen ’s function for layered m e d i a ...
...
LI
1.4.1 H istory ... 1.4.2 A new norm al m ode program : S A M P L E ... ... 161.5 M odel verification ... ... 18
1.6 S u m m a r y ... ... 11)
2 W a v e s in e l a s ti c m e d ia 21 2.1 F u n d am en tal e q u a tio n s ... ... 21
2.2 P la n e waves an d energy a b s o r p ti o n ... ... 25
2.3 Field contin u ity and boundary c o n d i t i o n s... ... 27
2.3.1 Solid-solid interface ... ... 27
2.3.2 An interface with a f l u i d... ... 28
2.3.3 Im p en etrab le b o u n d a rie s... ... 29
2.3.4 Infinite dom ains ... ... 30 2.4 Som igliana’s I d e n t i t y... ...
C O N T E N T S iv
2.6 T h e boundary integral p e r s p e c t i v e ... 38
2.7 D eterm ining th e boundary fields ... 39
2.7.1 A com bined equation m ethod... 41
3 T h e i n d i r e c t B I E m e t h o d 43 3.1 Numerical approxim ation of th e boundary i n t e g r a l ... 43
3.2 A flexible com bined equation m e t h o d ... 46
3.2.1 Surface-held equations ... . . . . 47
3.2.2 Null-field e q u a tio n s ... 47
3.2.3 Com bined equation m e t h o d ... 48
3.2.4 C om pressing th e B IE coefficient m a t r i x ... 50
3.2.5 A m odular approach to th e com bined eq uation m ethod ... 51
3.2.6 Specifying th e scatterin g problem in a flexible w a y ... 52
3.3 A consistency test using H uygens’ p r i n c i p l e ... 53
3.4 T roublesom e integral coefficients ... 54
3.5 Using H uygens’ principle to evaluate troublesom e integral coefficients in d i rectly ... 55
3.5.1 Several troublesom e e l e m e n t s ... 57
3.5.2 T h e p a rtic u la r s o l u t i o n s ... 58
3.5.3 T h e im p o rtan ce of the residue term ... 58
3.5.4 N onunique indirect c o e f f ic ie n ts ... 59
3.6 Using H uygens’ principle to te st th e entire B IE m ethod: the free-field te st . 60 3.7 T esting th e b o u ndary field: the null-field t e s t ... 61
3.8 M easuring error in num erical t e s t s ... 61
3.9 Exam ple: P lane wave scatterin g by s p h e r e s ... 62
3.9.1 P lane wave sca tte rin g from a rigid sphere in a f l u i d ... 63
3.9.2 Fluid sphere in a f l u i d ... 68
3.9.3 Solid sphere in a s o l i d ... 68
4 T h e f ie ld f r o m a p o i n t s o u r c e in la y e r e d m e d ia 76 4.1 T h e tw o-point b o u ndary value problem ... 77
4.1.1 T ransform ation to stress-displacem ent v e c to r s ... 79
4.1.2 T ransform ation to wave v e c to r s ... 81
4.1.3 T h e solution to th e hom ogeneous e q u a tio n ... 83
4.2 T h e sca tte rin g m a trix m eth o d ... 85
4.2.1 S catterin g m a trix for a hom ogeneous l a y e r ... 85
4.2.2 S cattering m a trix for an interface betw een l a y e r s ... 87
4 2.3 S cattering m a trix for im penetrable b o u n d a r i e s ... 91
4.2.4 S catterin g m a trix for an infinite h alf space an d th e im po rtan ce of R iem ann s h e e t s ... 94
4.2.5 P o in t Sources ... 95
1.3 L ad d er diagram for th e global m a trix m eth o d ... 98
4.4 A single source and r e c e i v e r ... 101
C O N T E N T S v
5 M o d e s in la y e r e d m e d ia 112
5.1 T h e discrete and continuous s p e c t r u m ... 113
5.2 T h e condition for inodes ... 117
5.2.1 Singular value decom position ( S V D ) ... I IS 5.2.2 T he channel m a trix m e t h o d ... 119
5.3 T h e search for m o d e s ... 120
5.3.1 R estricting tlve search area... 120
5.3.2 Using SVD to search for m o d e s ... 122
5.3.3 W eakly coupled channels ... 12-1 5.3.4 Searching for SH m o d e s ... 125
5.4 Verifying th e m ode s e a r c h ... 125
5.4.1 Node c o u n t i n g ... 127
5.4.2 C ontour I n t e g r a t i o n ... 127
5.4.3 C ontinuity and b o u ndary c o n d i t i o n s ... 128
5.4.4 R ayleigh’s P r i n c i p l e ... 128
5.5 F ive properties of m o d e s ... 129
5.5.1 R ayleigh’s principle ... 130
5.5.2 M ode ortho g o n ality ... 130
5.5.3 G roup v e l o c i t y ... . 131
5.5.4 M ode ex citation ... 132
5.5.5 M ode n o r m a l i z a ti o n ... 134
5.6 V ertical energy flow in a m ode ... 134
5.7 T h e G reen’s function for bounded waveguides by m odal s u m m a tio n ... 137
5.7.1 T he m ode series ... 137 5.7.2 Series convergence ... 139 5.7.3 T he G reen ’s function s i n g u l a r i t y ... 139 5.8 M odel verification ... 141 5.8.1 A rctic ocean m o d e l ... 141 5.8.2 G u ten b u rg e a rth m o d e l ... 150 6 C a ta lo g u e o f n o r m a l m o d e s 159 6.1 M odes as cylindrical waves ... 160
6.1.1 P ro p ag a tin g m o d e s ... 161
6.1.2 E vanescent m o d e s ... 161
6.1.3 P ro p er and im proper m o d e s ... 163
6.1.4 V ertical energy flux ... 168
6.2 M ode trap p in g m e c h a n i s m s ... 171
6.2.1 M odes tra p p e d by constructive reflection ... 171
6.2.2 M odes in bounded waveguides ... 178
6.2.3 Interface m o d e s ... 181
6.3 Identifying channels in th e channel m atrix m e t h o d ... 186
6.4 Deep bounded waveguides—T h e o r y ... 187
6.4.1 R ad iatio n m odes close to th e E JP branch l i n e s ... 187
6.4.2 T h e spacing of radiation m o d e s ... 190
C O N T E N T S vi
6.4.4 C orrespondence w ith th e unbounded w a v e g u id e ... 192
6.5 W eakly coupled sound channels— T h e o r y ... 195
7 B I E m e t h o d fo r la y e r e d m e d ia 2 0 4 7.0.1 C onstraints due to elem ent s i z e ... 204
7.0.2 C onstraints due to m ode p r o l i f e r a t i o n ... 205
7.1 T h e p artic u la r solutions for layered m e d i a ... 206
7.2 Interfaces and e d g e s ... 207
7.3 D em onstration of th e B IE m ethod for layered m e d i a ... 208
7.3.1 H em ispherical cavity in an ice p la te ... 209
7.3.2 Rigid sphere half buried in sand ... 229
7.3.3 C ircular ice d o m e ... 241
8 C o n c lu s io n s 255 8.1 Conclusions regarding layered m e d i a ... 256
8.1.1 Schem atic lad d er d i a g r a m ... 256
8.1.2 N orm al m odes an d S V D ... 257
8.1.3 T he m ode search program S A M P L E ... 257
8.1.4 T h e properties of m o d e s ... 259
8.2 T h e indirect B IE m e t h o d ... 261
8.2.1 F u tu re r e s e a r c h ... 262
8.3 C losing R e m a r k s ... 263
List o f F igures
2.1 T h e volum e and surface of integration in B e tti’s reciprocity relatio n ... 2.2 H uygens’ principle for wavefront re co n stru ctio n ... 3.1 Q u a d rila teral b o u ndary elem ent and n odes... 3.2 B oundary elem ents for th e rigid sphere... 3.3 P o la r plot o f th e sca tte red field due to the rigid sp h ere... 3.4 B oundary elem ents for the fluid and solid spheres... 3.5 P o la r plot of th e sca tte red field due to th e fluid sphere... 3.6 P o la r plot of th e scattered field due to the solid sphere... 3.7 P lo t of th e to ta l field surrounding the solid sphere, k r — 2.0... 3.8 P lo t of th e to ta l field inside and outside the solid sphere, k r = 2.0... 4.1 E J P b ra n ch cuts for vertical slowness... 4.2 T w o-port scatterin g m a trix netw o rk ... 4.3 S catterin g m a trix netw ork for a hom ogeneous lay er... 4.4 S catterin g m a trix netw ork for an i n t e r f a c e ... 4.5 S catterin g netw ork for an im penetrable plane boundary... 4.6 S um m ation o p erato r for a p oint source... 4.7 T h e schem atic ladder diag ram ... 4.8 T h e global m a trix ... 4.9 C om ponents in th e schem atic diagram can be combined in different ways. . 4.10 T w o sca tte rin g networks connected to g eth er... 4.11 T h e waveguide reduced to its sim plest form ... 4.12 A single source and receiver... ... 5.1 T h e slowness integral in the Fourier-Bessel tran sfo rm ... 5.2 “L”-sliaped region covering the E J P branch c u ts... 5.3 T h e search surface w \ resembles a conical hole close to a m ode... 5.4 Searching for a close m ode p a ir... 5.5 T h e search function for modes in th e A rctic waveguide... 5.6 P-SV m ode locations for the A rctic waveguide... 5.7 V ertical m ode functions for P-SV modes in th e A rctic w a v e g u id e ... 5.8 V ertical dependence of integrand in contour integral representation of rgg , 5.9 T h e norm al stress rzz as a function of range in th e A rctic w aveguide... 5.10 T h e vertical displacem ent uz as a function of range in the, ice . , ,
VII 32 36 45 65 67 70 71 73 74 75 84 86 88 90 92 98 102 103 1 0 4 105
107
108 115121
123 126 144 145 146 147 148 149L I S T O F F IG U R E S viii
5.11 W ave speed and density profile for tlie G uten b u rg e a rth m odel... 152
5.12 Search function for m odes in th e G uten b u rg e a rth m o d e l... 153
5.13 P-SV mode locations for th e G u ten b u rg e a rth m odel... 154
5.14 Vertical m ode functions for th e G u ten b u rg e a rth m odel... 155
5.15 V ertical dependence of in teg ran d in contour integral re p resen tatio n of u z . . 156
5.16 T h e displacem ent uz as a function of range in the G u ten b u rg e a rth m odel. . 157
5.17 D ispersion of the R ayle'gh interface wave in th e G u ten b u rg e a rth m odel. . . 158
6.1 P-SV and SH m odes for a free ice p la te ... 164
6.2 Evanescent P-SV m odes in th e A rctic w areguide... 165
6.3 P-SV m odes in a 3 in thick floating ice p la te a t 750 H z... 166
6.4 U pw ard energy flow in an infinite basem ent lay er... 170
6.5 P-SV m ode locations for th e tw o-layer leaky m ode en v iro n m en t... 176
6.6 V ertical m ode functions for p ro p er S-leaky m odes . 177
6.7 M ode locations in the unbounded and deep bounded w aveguides... 180
6.8 P ro p er leaky m ode in th e u nbounded waveguide and its ra d ia tio n co u n terp art in th e b o u n d e d ... 182
6.9 P-SV m odes for the two-layer m odel of A rvelo et a l... 185
6.10 T h e reduced schem atic lad d er diagram for a deep basem ent la y e r... 189
6.11 N orm al stress tzz as a function of range, in the infinitely deep waveguide and th e corresponding deep bounded w a v e g u id e ... 196
6.12 T h e reduced schem atic ladder diagram for a waveguide su p p o rtin g tw o weakly coupled channels a t o nce... 197
6.13 R esonant m odes in th e alm ost sym m etric w aveguide... 201
6.14 T ransfer of sound energy betw een weakly couped m odes... 202
7.1 B oundary elem ents for the hem ispherical cavity... 211
7.2 428 sources generating p a rtic u la r solutions outside th e hem ispheral cavity. . 212
7.3 T h e relative stren g th of m odes in th e m ode series... 213
7.4 Free-field te s t for th e ice cavity w ith incident P-SV w aves... 215
7.5 Null-field te st for th e ice cavity... 217
7.6 Incident field generated by a vertical p o in t force a t th e surface of th e ice. . 218
7.7 S cattered field due to th e cavity in the ice w ith incident P-SV w aves 219 7.8 P o la r plot of th e sca ttered field in a horizontal plane: incident P-SV waves. 220 7.9 D isplacem ent field along a line over the cavity: incident P-SV waves 221 7.10 Free-field te st for th e ice cavity w ith incident SH w aves... 222
7.11 Incident field generated by a a horizontal tw isting source in th e ice... 223
7.12 S cattered field due to th e cavity in the ice w ith incident SH w a v e s ... 224
7.13 P o la r plot of th e sca tte red field in a horizontal plane: incident SII waves. . 225
7.14 D isplacem ent field along a line over th e cavity: incident SH waves... 226
7.15 T o tal field a t the surface of the ice... 227
7.16 Vertical com ponent of th e sc a tte re d field uz a t th e ice surface . 228 7.17 Schem atic diagram of half-buried rigid sp h ere... 231
L I S T O F F IG U R E S ix
7.19 P o in t source locations generating th e p articu lar solutions for the half-buried
sphere... 233
7.20 P-SV m odes for th e shallow ocean with a sand layer... 23*1 7.21 SH modes for the sand la y e r... 235
7.22 Free-held te st for th e half-buried sphere... 230
7.23 S cattered field due to th e half-buried sphere in the vortical p lan e... 237
7.24 T o tal field in a. h orizontal plane a t the top of the sand lay er... 238
7.25 S cattered field in a horizontal plane a t the top of the sand la y e r... 239
7.26 P lo t of the sca tte red field a,t th e top of the sand layer... 240
7.27 Schem atic diagram of an ice dom e in a 3m ice p la te ... 2*1,2 7.28 B oundary elem ents for th e ice d om e... 244
7.29 SH m odes for th e A rctic waveguide a t 20 H z... 245
7.30 In terio r null-field p o in t locations for the ice dom e... 2*16 7.31 Sources g enerating p a rtic u la r solutions for perspective 0 (w aveguide) 247 7.32 Sources g enerating p a rtic u la r solutions for perspective 1 (ice dom e)...248
7.33 Free-field te s t for perspective 0 (w aveguide)... 249
7.34 Free-field te st for perspective 1 (ice dom e)... 250
7.35 Incident and to ta l field in a horizontal plane a t the ice surface... 251.
7.36 S cattered field and its vertical com ponent u~ in a horizontal plane a t the ice surface... 252
7.37 R adial u T an d az im u th al u,j, com ponents of the scattered field... 253
X
List o f Tables
3.1 S etup for scatterin g from a rigid sphere... 64
3.2 Free-field test sum m ary for rigid sphere b o u ndary... 66
3.3 Setup for scattering from a fluid sphere... 69
3.4 S etup for plane wave scatterin g from a solid sphere... 72
5.1 A rctic profile with floating ice layer. . . 143
5.2 T h e G u ten b u rg earth m odel... 1 5 1. 6.1 10 m sand layer on a lim estone sea floor... 171
6.2 Tw o-layer leaky mode en vironm ent... 178
6.3 A sim ple deep bounded waveguide. . . 179
6.4 T h e two-layer Jeaky m ode environm ent of Arvelo et a l... 184
6.5 A lm ost sym m etric fluid waveguide... 203
7.1 S etup for scattering from a hem ispherical cavity in an ice p la te ... 210
7.2 Setup for scattering from a rigid sphere half buried in a sand la y e r... 229
Acknowledgments
I thank Dr. Trevor Dawson for exemplary guidance through so many mathematical and computational difficulties. I thank the Canadian Departm ent o f National Defence for financial support throughout this work, through the Defence Research Establishm ent Pacific (DREP) and through the Esquirnalt Defence Research D etachm ent (EDRD), and I thank the scientists there who helped through many fruitful discussions and administrative support, especially Dr. John Ozard, Dr, D ave T homson, and Dr. Gary Brooke. I thank the members o f my exam ining committee, Dr. JohnV/eaver, Dr. Pauline van den Driessche, Dr. Robert Stewart, and Dr. H enrik Schmidt for their helpful criticism and
xii
I dedicate this publication to my parents, in the ye a r o f their forty-first wedding anniversary.
C h a p t e r 1
In tro d u ctio n
The ocean is an extrem ely complicated acoustic m edium .
T h a t is how Brekovskikli an d Lysanov begin th eir in tro d u cto ry tex t lor ocean acoustics [17]. I t is perhaps as m uch a w arning to the confident acoustician as to the beginner, for although the law s of acoustic wave m otion can be expressed rath er simply, their application to realistic environm ents takes ou d aunting com plexity, mainly because the ocean is highly variable over th e distances of interest. T h e w arning is app ro p riate here as well, w here the goal is to p redict th e way sound in tera cts w ith solid o bjects in the ocean, such as features of b o tto m topography, sub-b o tto m stru ctu re , or a surface ice canopy. Solid media com plicate the analysis because elastic waves in solids m ust be calculated using vector fields ra th e r th a n scalar fields, an d because th e actual properties of th e sea bed and ice canopy im p o rta n t p aram eters for th e analysis— can only be determ ined with difficulty.1 When ocean acoustics includes solid m edia th is way, it is often called geoacoustics to distinguish it from th e tra d itio n a l strictly fluid case, while em phasizing its affinity with geophysics.2
W hen analytically in tra c ta b le , th e theory of elastic wave motion can nevertheless
'F o r exam ple, see C lay and Medwin [28] who review m ethods of m easuring the properties of th e sea bed; H am ilton [54] who re p o rts the elastic p ro p erties of m any different kinds of sedim ents; and Brooke and O/.ard [18] who rep o rt the elastic properties of sea ice.
2T he aflininty w ith geophysics is widely accepted. Tolstoy and Clay [128, p,20(J] conclude th a t “Ocean acoustics is, in the final analysis, a branch o f geophysics.” Jensen e t al, [66, p.adj define a yeoacowitic model as “a m odel of th e real seafloor w ith em phasis on m easured, ex trap o lated , and predicted values of those m aterial properties im p o rta n t for th e m odelling of sound transm ission.”
C H A P T E R 1. I N T R O D U C T I O N
be applied to realistic geoacoustic environm ents using com puter m odels whose com putations im plem ent the governing equations. T h e m odels now in use are th e result of continual efforts to improve realism an d com putation speed.3 T he tw o goals are usually incom patible because greater realism engenders g reater m odel complexity, and therefore g reater execution tim e, except where m atched by m arked im provem ents in com puter h ard w are and softw are. T he added com plication of solid m edia is ju s t one exam ple of this. W hen we speak of realism , it need hardly be m entioned th a t all m odels are for th e m ost p a r t ab stra ctio n s of reality, in which a m yriad unaccountable physical properties are idealized by just a few dom inant ones. M oreover, it is the judicious sim plification of the real world th a t m akes a m odel especially useful and insightful. T hus, th e realism of a m odel can be posed m ost concisely in reverse, by citing a m odel’s few abstractions ra th e r th an its long list of omissions.
T h e goal of this thesis, th en , is to include th ree im p o rta n t p ro perties of th e ocean: horizontal stratifica tio n , solid features, an d wave p ropagation in th ree dim ensions. E ach of these have been m odeled to some degree independently, b u t only recently have researchers a tte m p te d to unite all three in a full-wave scatterin g m odel. T h a t is m y objective here.
By “full-wave,” I m ean th a t th e m odel u ndertakes to solve th e full scattering problem , posed m ath em atically as a b o undary value problem , w ith o u t m aking theoretically m otivated approxim ations from the o u tse t, such as the high frequency approxim ation of ray theory [66, ch a p t. 3], or th e horizontal forw ard propag atio n of th e parabolic equation m ethod [66, chapt. 6]. In principle, th e aim of a full-wave m odel is to solve th e b o undary value problem exactly, though in practice, of course, a num erical solution always entails approxim ation and erro r to som e ex te n t, and at tim es these m ay be disastrous. Progress in full-wave scatterin g m odels has advanced along tw o m ain fro n ts, using th e F in ite Differ ence (F D ) and B oundary Integral E q u atio n (B IE ) m ethods, whose relative advantages and disadvantages will be com pared below. I have chosen th e B IE m eth o d because it is favored for scatterin g in large dom ains [88] [27] such as th e ocean and e a rth .
T hree-dim ensional scatterin g m odels of any kind ten d to be difficult to use, and
3Many applications of ocean acoustics th a t use co m p u ter models are reviewed by Clay and M edwin [28], including echo ranging, m onitoring biological life, and m easurem ent of sea bed properties. More recently, an acoustic m ethod for d etec tin g changes in th e average te m p e ra tu re o f th e ocean has been developed to address concerns a b o u t global w arm ing [93]. A ctive and passive sonar are of g reat im p o rtan ce in undersea w arfare [59] [120],
C H A P T E R 1. I N T R O D U C T I O N 3
one always hopes th e re m ay be a sim pler approach based on some efficacious approxima* B u t a rigorous full-wave m odel is in m any ways a step tow ards im proved approxim ate models because it p erm its freehanded experim entation th a t is rarely, if ever, possible with physical experim ents. It p erm its, for exam ple, exam ination of the wave field everywhere in the problem to identify w here helpful approxim ations m ight be m ade, and it perm its variation of the problem itself to identify its m ost im p o rta n t features. A full-wave model also provides a benchm ark to judge th e accuracy of faster bu t approxim ate “sh o rt-c u t” m ethods.
T h e sca tte rin g of geoacoustic waves is a vast topic and one cannot hope to touch its m any sides in a single p ro ject. My tre a tm e n t of scattering is therefore restricted in three im p o rta n t ways. F irs t, I only consider detei’m inistic scattering, in which the shape and elastic p roperties of th e sca tte rin g features are com pletely specified; in c o n tra st to statistical scattering, arguably of equal im portance, in which small-scale random roughness is handled statistically. B ut ju s t as determ inistic m odels have served as the basis for a statistical tre a tm e n t of im pen etrab le roughness [37] [102], my own model might serve for penetrable roughness, though I do n o t explore th e possibility here. Secondly, I only consider harm onic (con stan t frequency) wave scatterin g because the B IE m ethod solves ju s t one frequency at a tim e, an d th a t in itself is a considerable undertaking. In principle, however, the model could be extended to tra n sie n t waves using F our'*r synthesis in the frequency dom ain as others have done [67][50]. T hirdly, the em phasis will be on shallow oceans, in which the sound is likely to in te ra c t significantly w ith the sea bed, a t acoustic frequencies less than I kHz, for which sea w a te r is essentially tra n sp a re n t,'1 though some geophysical applications a t much lower frequency are also included. Am ong the problem s I a tte m p t here are
• norm al m ode com p u tatio n s in an A rctic ocean model, and the G u ten b u rg earth model from geophysics;
• sca tte rin g of plane waves from rigid, fluid, and solid spheres for which analytic solu tions are available;
• sca tte rin g (750 Hz) from a hem ispherical cavity on the floor of a 3 m ice plate a t 750 Hz;
C H A P T E R 1. I N T R O D U C T I O N 4
• sca tte rin g (200 Hz) from a 1.193 m rigid sphere half buried in sand in 50 m of w ater;
• sca tte rin g (20 Hz) from an ice dom e in a floating ice layer in th e A rctic ocean.
H ut the new m odel is likely to find applications w herever the sca tte rin g of waves in layered m edia is o f in te re st, in geophysical prospecting5, n o ndestructive testin g using ultrasonics, and electrom agnetics for exam ple.
G eoacoustic m odeling has grow n ou t of th ree great fields of study: elastic wave propagation, ocean acoustics, and geophysics. It is im possible to survey its m anifold history here. C hin-B ing et al. [26], for exam ple, catalogue th irty -six research-oriented com puter models in ocean acoustics alone, and these are classified according to seven m ain m odel types,6 whose underlying theory is given in an excellent com panion p u blication by Jensen e t al. [66]. D etailed reviews of roughly th e sam e nu m b er of geophysical m odels have been e d 'te d by D oornbos [39], Bolt [14], and Chin et al. [25]. M ore to our purpose, in this in tro d u cto ry ch ap ter I will review th e m otivation for th e m ain features of th e m odel, then introduce th e m ain difficulties facing every full-wave sca tte rin g m odel, th e n survey recent developm ents concerning the B IE m ethod an d layered m edia, a n d th en outline how th e new model will be verified.
1.1
M o tiv a tio n
1.1.1
N eed for stratified m edia
G eoacoustic waves depend on th e distu rb an ce th a t causes th em and th e elastic p aram eters of the m edium they traverse. In th e ocean, th e elastic p aram eters vary contin uously alm ost everyw here according to th e w ater salinity, te m p e ra tu re , an d pressure [17], and they m ay ju m p suddenly as a t th e tran sitio n from fluid to solid a t th e ocean floor. In th e absence of d istin c t scatterin g features, the length scale of th e horizontal variation is
‘ Seismic waves have been used by oil geologists in the search for d istin ct porous s< elementary stru ctu res, such as a P innacle reef, th a t can signify an oil reservoir [117, sect. 23.6],
flT hc seven m ain m odel types used by C hin-B ing e t al. are th e ray, p arab o lic eq u atio n , norm al mode, contour integral, coupled m ode, finite elem ent and finite difference models. O ddly enough, th e BIG m ethod for scatterin g w as not included. Presum ably th is is because B IE models for geoacoustics are fairly recent, and so far only th eir inventors use them w ith confidence. Jensen e t al. [66] include a discussion of B IE m ethods for tw o-dim ensional fluid media.
C H A P T E R 1. I N T R O D U C T I O N 5
often m uch g reater th a n th a t of the vertical, and greater th an the acoustic wavelengths of in terest, so it is usual to tre a t th e ocean as a horizontally stratified waveguide m ade up of hom ogeneous layers; th e elastic param eters of each layer approxim ating the vertical varia tio n in stepw ise fashion [16] [44] [66] [70] [128]. T h e vertical reverberation th a t dom inates geoacoustic d a ta can b e realistically m odeled this way. F eatures of bo tto m topography or surface ice can th en be m odeled as p en etrab le elastic inclusions em bedded in the layered waveguide. O f course, no t all m odels assum e horizontally invarian t m edia, such as the fi n ite elem ent, finite difference, parabolic equation and geom etric ray models, for exam ple. N evertheless, a locally stratified m odel is alm ost always used in practice, if only because th e horizontal variation of th e ocean is rarely known with certainty.
1.1.2
N eed for solid m edia
E arly m odels in ocean acoustics consisted of strictly fluid layers, with solids in cluded ap p roxim ately simply as fluids, by o m ittin g shear stresses alto g eth er.7 Increasingly, however, a tte n tio n has been given to low frequency propagation and shallow ocean envi ronm ents, for which th e im portance of shear waves in solid m edia to transm ission loss in shallow oceans, reverberation, an d scatterin g is now well established [17] [66] [77] [42]. To om it th em often precludes im p o rta n t resonance effects and energy loss m echanisms, and at tim es th e d o m inant p a r t of th e physical reality.8
1.1.3
N e ed for th ree-dim en sion al m odels
A tw o-dim ensional propagation m odel assum es a high degree of sym m etry in the wave field so th a t th e field in th ree dim ensions is com pletely determ ined by its values in a single v ertic al plane. In planar sym m etry, th e m edia, source, and wave field are all assum ed to be inv arian t in one horizontal direction; the represen tativ e plane being any plane p erpendicular to th e direction of invariance. T h e planar approach is often inadequate
7T he review of earlier num erical models for ocean acoustics by DiNapoli and D avenport (1970) [08] deals entirely w ith fluid m edia, although th ey m ention th a t a solid bottom could be included by term in atin g the w ater colum n below using th e equivalent reflection coefficient for solid sedim ent layers.
8Scholte interface waves on th e sea floor, for exam ple, do not occur in strictly fluid m edia (see S ection(6.2.3)).
C H A P T E R 1. I N T R O D U C T I O N 6
because 1) the elem entary source is a uniform line source lying parallel to th e direction of invariance, which is n o t typical of th e com pact sources encountered in reality; 2) because spherical spreading for body waves, and cylindrical spreading for tra p p e d waves (norm al m odes) ca n n o t occur; and 3) because out-of-plane scatterin g is precluded by th e assum ed sym m etry. T hese lim itatio n s have been overcome to some degree using “tw o-and-a-half”- dirnensional models, in which th e field due to a single p o in t source is co n stru cted using a Fourier transform of th e p la n a r field, while th e elastic m edia rem ains in v arian t in one direction. Faw cett an d Daw son [46] applied the m eth o d to a p lan ar B IE m odel of acoustic wave sca tte rin g in fluid waveguides, an d Schm idt has since applied it to solid layers [110]. S catterin g from com pact three-dim ensional inclusions cannot be tre a te d th is way.
In axial sym m etry, th e m edia, source, an d wave field are ro tatio n ally sym m etric ab o u t a com m on vertical axis, and th e represen tativ e plane is a half plane w ith one edge along th e vertical axis. C ylindrical m odels are som etim es used as “building blocks” to assem ble larger lion-sym m etric dom ains in a piecewise m an n er [45] [66, sect. 5.10], b u t the extension to solid m edia has yet to be m ade. H ere again, an im p o rta n t in stan ce of out-of- plane s c a tte r, m ode conversion betw een vertically (P-SV ) and horizontally (SH) polarized wave m otion, is precluded by th e assum ed sym m etry.
In general, th e sca tte rin g problem does n o t suit eith er th e p lan ar or axial ideal ization an d a full three-dim ensional tre a tm e n t is required.
1.2
T w o card in al d ifficu lties in sc a tte r in g
A geoacoustic sca tte rin g m odel faces tw o fun d am en tal difficulties. T h e first is th a t of extrem e length scales. To represent th e in teractio n of waves w ith th e shape and sudden c o n tra st of th e inclusion, th e m odel m u st replicate th e equation of m otion and boundary conditions in its vicinity oil a scale m uch less th a n th e inclusion’s dim ensions and the geoacoustic w avelength. A t th e sam e tim e, th e m odel m u st p ro p ag ate waves over very large distances, w hether to d ista n t receivers or from d ista n t sources. H andling both scales a t once is ra th e r like m easuring very small and very large distances w ith th e sam e m eter stick; the stick m ay be w ell-suhcd to m easuring one length scale o r th e o th er, b u t
C H A P T E R 1. I N T R O D U C T I O N 7
n o t to b o th a t once. T h e disparity can be very severe in ocean acoustics when the inclusion lies in unconsolidated sedim ents where the shear (S) w avelength— the length of the m eter stick— can b e very sm all. I am no t aw are of any full-wave scatterin g model (including my own developed here) th a t p u rp o rts to include shear- wave scatterin g in such extrem e cases.
T h e second fun d am en tal difficulty is th a t of com pleteness, for the model m ust have access to all nine com ponents of th e elastic wave f ie ld - th r e e com ponents of displacement, and six unique com ponents of stress— to m atch boundary conditions on an arb itra ry in terface betw een solid m edia in three dim ensions. Such rigor is uncom m on in propagation modeling. Usually only p a rt of th e field is required, perhaps the norm al stress in ocean acoustics, or th e displacem ent vector in seismology.
1.3
B a ck g ro u n d for th e B IE m e th o d
1.3.1 R eview o f recent developm ents
A m ong th e recent developm ents in th e B IE m ethod for ocean acoustics a,re con trib u tio n s by S chuster and Sm ith (1985) [112] whose B IE m ethod for rigid o bjects in a tw o-dim ensional fluid m edia includes re v erb eran t interactions betw een the scatterin g inclu sion and th e layered m edia approxim ately using a Born series. Seybei't and Casey (1988) [115] were possibly th e first to apply the B IE m ethod to pen etrab le dom ains with a view tow ards ocean acoustics, although they considered unbounded homogeneous dom ains rath er th a n layered. S eybert an d W u (1988) [116] applied the B IE m ethod to a homogeneous fluid halfspace. Lu (1989) [86] used his hybrid ray-m ode m ethod for layered m edia in a two- dim ensional B IE m eth o d for strictly fluid m edia including p enetrable scattercrs. Dawson an d Faw cett (1990) [31] considered scatterin g in two-dim ensions by im penetrable deform a tions in a strictly fluid waveguide, which they la te r (1990) extended to “tw o-and-one-half” dim ensions [46] as cited earlier, an d which Dawson (1.991) [32] extended to long repeated boundary deform ations using sca tte rin g m atrices. Dawson (1991) [33] then developed a three-dim ensional B IE m ethod for im penetrable inclusions in a strictly fluid waveguide. G erstoft an d S chm idt (1991) [50] developed the first tw o-dim ensional B IE m ethod for a r b itrarily layered elastic m edia following th e m ethod of K aw ase (1988) [67] in geophysics,
C H A P T E R 1. I N T R O D U C T I O N
who com puted th e tim e dom ain response of a canyon in a hom ogeneous solid h alf space in two dim ensions, and Schm idt (1993) [110] la te r extended th eir m odel to “tw o-and-one-half” dimensions as cited earlier. X u and Yan (1993) [135] considered tw o-dim ensional scatterin g in the elem entary oceanic waveguide using norm al m odes,9 which th ey la te r used for source localization in a shallow ocean w ith a large rigid inclusion (1994) [136]. W u (1993) [132] considered three-dim ensional scatterin g in th e elem entary oceanic w aveguide using b o th the m ethod of images and norm al m odes. Eliseevnin and Tuzhilkin (1995) [41] apply Kirch- hoff’s app ro x im atio n to th e B IE m ethod for rectangular vertical screens in the elem entary oceanic waveguide. I t rem ains to develop a three-dim ensional B IE m eth o d for layered m edia involving solids.
1.3.2
B IE m eth o d in p ersp ectiv e
P erh a p s th e best-know n num erical m ethods for b o u n d ary value problem s are the F inite Difference (F D ) and F in ite E lem ent (F E ) m ethods. In each case th e elastic dom ain is subdivided into volum e elem ents whose dim ensions are sm all com pared w ith b o th the scattering inclusion a n d th e wavelength, an d whose corners are nodes a t which th e unknow n field variables are to be com puted. T h e field a t each node is re la ted to its neighbors using a num erical ap p roxim ation to th e equation of m otion. T his gives a sparse b anded system of equations, whose solution yields th e field variables a t every node th ro u g h o u t th e volum e. T h e m eth o d s are very flexible because th e scattering body can have any shape an d the elastic p a ra m e te rs can vary alm ost arbitrarily.
T h e m ain disadvantage of th e FD and F E m ethod is th a t a very large n um ber of nodes are needed to span a large dom ain, especially in th ree dim ensions. In geoacoustics th e dom ain is reasonably assum ed to be infinitely large, b u t th e grid of nodes m u st be term in ated som ew here, by a b o u ndary contrived to m inim ize any effect on the wave p ro p agation, th ereb y m aking a seamless connection to th e o m itted u n b o unded dom ain. Being im perfect, these false boundaries m ust be kept as far a p a rt as possible, hence th e dom ain spanned by th e grid of nodes m ust be as large as possible. T h e F D and F E m eth o d s have
“T he elem entary oceanic waveguide consists of a single hom ogeneous fluid layer bounded above and below by free ami rigid boundaries, respectively.
C H A P T E R 1. I N T R O D U C T I O N 9
therefore been lim ited to tw o-dim ensional m odels [40] [66] [47], or to sh o rt range (less than te n w avelengths) three-dim ensional models [20] [19]. In light of th e first cardinal difficulty of scattering, th a t of handling incongruous distances, it would app ear th a t the FD .aid FE m ethods are w ell-suited for short range propagation in the vicinity of the scattering inclu sion, h u t they face considerable difficulties in the far field due to com puter lim itations. T he second difficulty, th a t of com pleteness, m ust be resolved by th e finite difference or element scheme relatin g each node to its neighbors. Interested readers can consult the references.
T h e B oundary In teg ral E quation (B IE ) m ethod is an a ltern ativ e approach in which only the boundary of th e scatterin g body need be subdivided into small surface elem ents, ra th e r th a n th e e n tire dom ain [8]. A ssociated w ith each elem ent are nodes, also on the boundary where th e field is to be com puted. T he field at each node is related to th a t a t all th e others thro u g h a b o u ndary integral representation of the field, giving a full system of equations whose solution yields th e field a t all nodes sim ultaneously. T h e field can then bo com puted a t points off th e boundary, once again using its integral representation. As we will see, th e m eth o d is a realization of H uygens’ principle for wave reconstruction, in which th e field w ithin a dom ain is represented as th e sum of wavelets radiated from secondary sources on a boundary. T he field ra d ia ted by each wavelet is the G reen’s function for a suitably chosen dom ain in th e absence of sca tte rin g .10 In an unbounded dom ain, the G reen’s function propag ates a wavelet over large distances in a single evaluation, w ithout an intervening grid of nodes as in the FD and F E m ethods. T his is why th e B IE m ethod can accom m odate unbounded dom ains. T he tw o cardinal difficulties of scatterin g are therefore relegated to th e G re en ’s function; for it m u st propagate a wavelet over a rb itra ry distances between nearby nodes on the boundary, or from th e boundary to d ista n t receivers and it m u st be com plete, re tu rn in g all com ponents of th e elastic field for a wavelet.
T h e m ain adv an tag e over the F D and F E m ethods, then, is a considerable reduc tion in th e num ber of nodes, from a grid spanning a large three-dim ensional dom ain, to a grid spanning a finite tw o-dim ensional dom ain of th e boundary of th e scatterin g inclusion.
10A G reen ’s function is th e field due to a fundam ental po in t source. T h e fundam ental source in th e HIM m ethod for elastic waves is a point force, which m akes th e G reen ’s function a second rank tensor (see Section (2.4)). As we will see, th e field due to a wavelet is th e superposition of two fields originating from the sam e point: one field going as th e G reen’s displacem ent tensor (sim ple), and th e o th e r as th e th ird -ran k G reen’s stress tensor (com plex) (see equation (2.62)).
C H A P T E R 1. I N T R O D U C T I O N 10
U nfortunately, this gain is bought a t th e price of added com plexity, for n o t only does th e m ethod entail the inversion of a full m a trix system , but ce rtain B IE coefficients in th a t m atrix are difficult to com pute due to th e singularity of th e G reen ’s function. To illus tra te , im agine com puting the wavelet superposition a t a p rin t directly on the surface where those wavelet sources are continuously d istrib u te d . Very close to th e co m p u tatio n point the G reen’s function of the wavelets becomes infinitely large, b u t th e b o u ndary integral over the wavelets, defined in the C auchy P rincipal Value (C P V ) sense, rem ains bounded. T h e situ atio n is fu rth e r com plicated in layered m edia because its G reen’s function m u st be expressed in term s of integral transform s th a t are com putationally intensive; m ore so when evaluated close to th e source th a n far away. Indeed, a good p o rtio n of th is thesis is dedi cated to com puting th e com plete G reen’s function for layered m ed ia in b o th th e n ear and far field of the source. In the final analysis, th en , th e B IE m eth o d is no t a short cu t around th e FD and F E m ethods, b u t it ranks w ith them as an altern ativ e num erical m ethod th a t is especially suited to large three-dim ensional dom ains, and this is im p o rta n t for geoacoustics.
To evaluate troublesom e singular m atrix coefficients in the B IE m eth o d , some modelers have ex tra cted the singularity from th e integral an d in teg rated it analytically, leaving a regular com ponent to be in teg ra te d num erically w ith o u t difficulty [31] [33] [135], while o thers have perfected a direct num erical in teg ratio n [52] [53], b u t sim ilar m ethods have yet to be accom plished for three-dim ensional elastic waves in layered elastic m edia due to its form idable G reen ’s function. I follow a ra th e r different course.
1.3.3
A new indirect BIE m eth o d
One of the principal innovations in this thesis is th e in d irec t com p u tatio n of tro u blesome B IE coefficients, by inferring th e ir values from a fam ily of know n solutions to the integral equation. T h e m ethod is sim ilar to one proposed by Niku an d B rebbia [94], in which p articu lar solutions to th e boundary integral eq uation are used to infer all b o u ndary integral coefficients, w hether num erically troublesom e or not. T heirs was a general discus sion, removed from th e details of any given b o undary value problem and w ithout num erical exam ples, b u t they an ticip ated th a t th e indirect com putations m ay be num erically u n s ta ble. T h a t is w hat I found when trying to infer all coefficients in early trials w ith strictly
C H A P T E R . 1. I N T R O D U C T I O N II
fluid tw o-dim ensional problem s. I therefore modified the m ethod to com pute only the tro u blesome coefficients indirectly, while com puting all straightforw ard coefficients numerically. This gave good re su lts.11
T h e indirect B IE m eth o d can be applied m ore generally, to com pute any HIE coefficients th a t are num erically troublesom e, no t ju st those troubled by singularity. This m akes it ideal for layered m edia, where convergence problems in the G reen’s function can lead to additional troublesom e coefficients.
1.3.4
A flexible com bined BIE m eth o d
T h e B IE m eth o d described here is very flexible. It perm its several scattering objects a t once, th e ir boundaries being disconnected, in co n tact, or em bedded one inside another. T hey can be p en etrab le fluids or solids, or im penetrable with linear boundary conditions, w hether passive an d actively v ib ratin g (to model a hydrophone for exam ple). T hey can have edges an d corners. In layered m edia, they can pass through J lie interfaces betw een layers, and ex ten d above or below its horizontal plane boundaries.
All of this m ay ap p e ar excessive, b u t alm ost all of these features are required for layered m edia, even in relatively sim ple problem s. For in layered m edia, the scattering object m ay be em bedded p a rtly in solid layers, and partly in fluid, calling for different boundary conditions in each. A nd where the o b ject passes through a solid-fluid interface, the boundary field m ay be discontinuous, calling for tre a tm e n t much as if the boundary had an edge. If the o b ject con stitu tes a deform ation of an im penetrable plane boundary (a ridge in a floating ice p la te , for exam ple, or a m o u n tain on otherw ise level te rra in ), then it may extend outside th e layered m edia, and its boundary may be penetrable in some places, but im penetrable in o thers. A ctively v ib ratin g boundaries should be included in all BIE models, if only to im plem ent th e com prehensive free-field test described in this thesis. P erhaps the only ex tra v ag an t featu re, th en , is to p erm it more than one sca tte rin g o b ject a t a tim e, B ut as we will see, provision for several o b jects is only a m inor extension of the provision for a single o b ject th a t is penetrable.
u U nfortunately, J did n o t come across Niku and B reb b ia’s work until after prelim inary trials Jliowed th a t instability was a problem when inferring all coefficients indirectly, and after I modified th e indirect m ethod for troublesom e coefficients only.
C H A P T E R 1. I N T R O D U C T I O N 12
A general. B IE mode! is com plicated by th e fact th a t th e in teg ra n d of th e integral representation of th e field generally involves th e boundary values of b o th th e field and its normal derivatives. Hence, changes to the boundary conditions m ay change th e form of the integral equation being sol ved, which in tu rn m ay have drastic consequences for num erical stability. For exam ple, if the norm al derivatives of the represented held variable are pre scribed on the boundary (th e N eum ann boundary conditions in scalar wave th eo ry ), then the integral represen tatio n yields an integral equation of th e second kind, which ordinarily have excellent num erical stability [4].12 B u t if the represented field is itself prescribed on the boundary (th e DirichJct boundary conditions in scalar wave theory), th e n an eq uation of the first kind results, which are notorious for being num erically u n stab le [4],13 T h en again, for penetrable sca tte rin g objects, the integral rep resentation rem ains as a m ixed eq u atio n —of the second kind in th e represented field variable, and of th e first in its norm al derivatives. T he general properties of the mixed equation have apparen tly n o t been studied, a t least not in th e geoacoustic literature. It is n o t known, for instance, w hether the occasional non-uniqueness th a t affects equations of th e second kind, and the in h ere n t in stab ility of the first kind, m ay not b o th affect the solution of a. m ixed equation. All of this m eans th a t a general B IE model, though proven successful for one or an o th er problem , m ay nevertheless become unstable and fail for o thers, possibly if only the boundary conditions are changed. P recautions against num erical instab ility are i Imrcfore an im p o rta n t p a r t of a general BIE model.
In my ow n m odel, the th re a t of num erical instability is countered in th e following way. F irst of all, it has several built-in num erical tests— consistency an d com prehensive tests th a t can verify different stages in th e model for every application. T hese tests can be applied in any B IE m ethod for sca tte rin g problem s, b u t they h ave a p p a ren tly received
l2T hc integral representation for the elastic displacem ent field (Som igliana’s id en tity (2.48)) reduces to an equation of the second kind in displacem ent when zero-traction b oundary conditions are enforced from th e o u tset, in scatterin g from a hollow t cnch or cavity in th e e a rth , for exam ple. In acoustic (scalar) wave scatterin g , th e H elm holtz integral equation f .r the pressure sim ilarly reduces to an equation of th e second kind for im p en etrab le rigid objects [73] [22],
''’iSomigliana's identity (2.48) reduces to an equation of th e first kind in tra c tio n when zero-displacem ent bou n d ary conditions are enforced, in scatterin g from an im penetrable rigid body in welded c o n ta c t with a solid, for exam ple. In acoustic (scalar) m odels, a different integral representation of th e field is ordinarily used for th e N eum ann an d Dirichlet problem s, to give an equation th e second kind in each [22] [31] .
C H A P T E R 1. I N T R O D U C T I O N 13
little m ention in th e lite ra tu re , probably because they are little known.
Secondly, in principle I follow the com bined equation m ethod developed by Schonck [106] for th e N eum ann problem for acoustic waves. T he combined m ethod exploits the fact th a t a b o undary integral represen tatio n of the held also takes different forms according to th e location of th e point where the field is represented: w hether inside, outside, or on the boundary of in teg ratio n . Schenck began w ith the Helm holtz integral representation for the acoustic pressure, applied a t points on th e boundary (th e surface-Jkld equation), giving an equation of the second kind in the unknow n pressure field. Using the theory of integral equations, he d em o n strate d th a t its solution is not unique at isolated critical frequencies, where ruinous num erical instability occurs. To force a unique solution, ho over-determ ined the problem by com bining the integral representation at points outside the boundary (the null-field equation, an equation of the first kind) with the surface-field equation in the num erical m ethod. T he extension now to elastic waves and penetrable objects pushes th e com bined m eth o d far beyond Schenck’s original analysis: from the Helm holtz integral equation for acoustic (scalar) fields, to Sornigliana’s identity for elastic (vector) fields; from a surface-fteld eq uation strictly of the second kind for the Neum ann problem , to one th a t is generally m ixed for m any different problems: an d from one im penetrable object, which requires ju s t one integral represen tatio n for the sole penetrable dornai n , to several penetrable o bjects, which calls for a different integral representation in each penetrable dom ain. It is im possible to p redict th e consequences of the. extension. At the very least, we can expect to solve th e equivalent N eum ann problem for clastic waves (traction-free boundaries when Som igliana’s id en tity is used) with m uch the same success Schenck had w ith acoustic waves, As I show in the exam ples in this thesis, however, the m ethod works for a much wider iuuge of problem s, as intended from th e o u tset.
T h e new m odel’s broader success is due in p a rt to its flexibility, T he surface and null-field equations can be com bined in any proportion so'ccted by the user, so th a t instabilities, if th ey occur, m ight be overcom e by changing the proportions one way or another. To this end the user can choose any points on the boundary of the scatterin g o bject to apply th e surface-field eq uation, ra th e r than always using the elem ent nodes as D IE m ethods ordinarily do. T h e num ber of points m ay actually exceed the num ber of
C H A P T E R 1. I N T R O D U C T I O N 14
nodes on the boundary, over-determ ining the solution to th e surface-field eq u atio n to a limited degree, even before it is com bined w ith th e null-field equation. (T he freedom to choose boundary points also m akes it possible to avoid nodes on edges of th e boundary, or along the interface betw een fluid layers intersected by th e b o u ndary in layered m edia, where th e surface-field equation has a m ore com plicated specialized form .) In a sim ilar way, the user can choose any points in the e x u rio r dom ain to apply th e null-field eq u atio n . In my experience, I have on occasion faced w hat a t first appeared to be disastrous num erical instability, but I have yet to encounter a problem where it could n o t be overcom e using this flexible com bined eq uation m ethod.
As with m ost m odels, th e range of problem s th a t can be solved is lim ited by the com puter. As a rule, B IE m ethods are biased tow ards low an d m id-frequency sca tte rin g , in which the dim ensions of th e o bject are no t large com pared to th e w avelength. For as the frequency increases, sm aller boundary elem ents are needed to sam ple th e small- scale variations in th e field, thereby increasing the n um ber of b o u ndary elem ents, as well as the com putation tim e an d m em ory dem ands. A t this stage, the indirect coefficient com putations ap p e ar lim ited by th e n um ber of troublesom e b o u n d ary elem ents encountered a t once. For not only does th e com putational effort increase w ith th e n u m b er of troublesom e elem ents, but the indirect com putations become unstable. O rdinarily, however, less th an six troublesom e elem ents would be encountered a t a tim e, which th e indirect m e th o d now m anages w ithout difficulty.
1.4
A c o m p le te G r e e n ’s fu n ctio n for la y ered m ed ia
T h e B IE m ethod requires all nine com ponents of th e G reen’s function in b o th the near and far field of a point source. N um erical m ethods for com puting th e G reen’s function in layered m edia have developed independently, w ithout reference to th e B IE m ethod.
1.4.1 H istory
Thom son (1950) [124] and Haskell (1953) [56] took th e first steps tow ards com puting th e G reen’s function with th eir p ro p ag ato r m a trix m ethod for p lan e waves travelling
C H A P T E R 1. I N T R O D U C T I O N 15
in layered m edia. M any im provem ents have been m ade since: Haskell (1964) [57] and H arkrider (1964) [55] added source excitation; G ilbert and Backus (1966) [49] form ulated th e fu ndam ental th eo ry behind th e p ro p a g ato r m a trix m ethod; H udson (1969) [60] extended th e m eth o d to han d le p ropagation in th ree dimensions; and K en n ett (1979) [69], whose work features prom inently here, recast th e m a trix problem using a num erically stable recursive scheme. T here are m any text-book tre a tm e n ts of the m atrix m ethod an d its variations such as Ew ing, Ja rd e tsk y an d Press (1957) [44], Takeuchi and Saito (1972) [122], Aki and R ichards (1980) [1], Brekliovskikh (1980) [16], an d K en n ett (1983) [71], and Jensen ct al. (1994) [66].
P erh ap s th e m ost n o tab le developm ent in ocean acoustics has been the global m atrix m eth o d due to Schm idt and Jensen (1985) [107], which is th e basis of the SAFARI m odel for sound p ropagation [109]. SA FA RI now ranks am ong the most, widely-used models in ocean acoustics, b u t its wave propagation m odule is nevertheless incom plete for a three- dim ensional BJE m eth o d , oecause it only com putes p art of the field (th e norm al stress, and vertical a n d h orizontal displacem ent for P-SV waves excited by axially sym m etric sources), and because it becom es in accu rate w ithin a few wavelengths of a point source (due to asy m p to tic tre a tm e n ts of Bessel functions in its F ast Fourier Transform integration scheme). D r. S chm idt has recently released an upgraded version of SAFARI called OASES [ I I 1], built in p a rt on th e work of Kim [74], th a t com putes all field variables due to explosive point sources an d point forces, w ith an option to use a full integration schem e th a t is accurate in th e near field. To m y knowledge, OASES is the first general purpose model adequate for a three-dim ensional B IE m ethod.
B o th SA FA RI and OASES are w avenum ber (contour) in tegration models inasm uch as the field is represented as th e superposition of cylindrical waves whose wavenumbers span th e positive real axis. I follow a ra th e r different approach, based ori a norm al mode model o f wave propag atio n , in th e hope th a t it will be m ore com putationally efficient in the three- dim ensional B IE m eth o d where the G reen ’s function m u st be evaluated several million tim es. A norm al m ode m ethod m ay prove more efficient since it represents the field as a series of d istin ct waves, th e norm al m odes, ra th e r th an as an integral over a continuous spectrum of waves. T h e characteristic vibrations of each m ode can be com puted beforehand,
C H A P T E R 1. IN T R O D U C T IO N 16
and saved, to speed up subsequent evaluations. A m ong th e existing m ode m odels for ocean acoustics is SNAP by Jensen and F erla [65], K R A K EN by P o rte r [99], an d m ore recently ORCA by W estw ood [131]. T hough proven successful for a wide range of ocean propagation models, none of these is adequate for th e B IE m eth o d , m ainly because th ey do n o t com pute all com ponents of th e elastic field, an d because th ey concen trate on only p a r t of th e m ode series— th e d o m inant P-SV m odes v ib ratin g in th e w ater colum n— while o m ittin g short- range, strongly evanescent m odes and SH m odes. A m ore com plete n orm al m ode m odel is therefore required.
1.4.2
A new normal m o d e program: S A M PL E
To com pute the G reen ’s function by m odal sum m ation, I develop a scattering m a trix tre a tm e n t of layers th a t in effect converts K e n n e tt’s recursive schem e to a global m a trix m ethod. A t the sam e tim e, I in troduce a schem atic re p resen ta tio n of vertical propagation through layered media,, m uch like th e diagram s an engineer m ight use to analyze a linear system . T his aids th e discussion because it helps us visualize how m a trix m ethods for lay ered m edia work (or fail). T here are m any m atrix m ethods for layered m edia, an d th ey can all be illu stra te d using diagram s of th is kind.
For sim plicity, m ost norm al m ode program s search for a m ode a t a single dep th , usually in the w ater colum n, w ith th e effect of layers above an d below being represented by an up p er and low er reflection m atrix . B u t as we will see, searching a t a single d e p th is likely to miss m odes whose vibrations do n o t extend thro u g h th a t d ep th . To ensure th a t all m odes are d etected , I te st for singularity a t m any d epths sim ultaneously by com puting the singular values of the global m atrix . T h e m ethod has been im plem ented in a robust norm al m ode program called SA M PLE— an acronym for Seism o-A coustic M o d e P ro g ra m for L ayered E nvironm ents. T h e m odel is new inasm uch as it:
• u ndertakes to find all p ropagating m odes, plus a significant p o rtio n of th e long series of evanescent m odes, needed for near-field propagation;
• com putes all nine com ponents of th e three-dim ensional elastic field due to a p oint source;
C H A P T E R 1. IN T R O D U C T IO N 17
• su p p o rts a general class of seismic sources, including general point forces and couples w ith an d w ith o u t m o m en t;14
• uses a m ode search m ethod based on singular value decom position th a t can a u to m ati cally identify an d resolve uncom m only close m ode pairs (double roots) th a t may occur w hen th e layered m edia supports weakly coupled sound channels.15
SA M PL E has also been a valuable research tool, leading to a com prehensive un d erstan d in g of m ode behavior, such as
• upw ard energy flow in an infinitely deep basem ent lay er;16
• w ell-behaved leaky m odes on th e p ro p er R iem ann sh eet;17
• a Rayleigh m ode “of th e second kind” for a solid layer in welded co n tact with a rigid b o u n d a ry;18
• th e correspondence betw een a deep bounded waveguide and its u nbounded c o u n te rp a rt. 10
A s we will see, a com plete norm al m ode representation of the G reen’s function is only possible w hen th e layered m edia is bounded above an d below by perfectly reflecting plane boundaries. T h u s it will be assum ed, as it often is [17], th a t th e ocean-atm osphere interface is a traction-free plane boundary, an d th a t th e sea bed is term inated a t d ep th by a reflecting in terface.20 A com prehensive norm al m ode approach for arb itra rily elastic m edia has ap p a re n tly never been rep o rted , so th e technique developed here, and the properties of m odes it has revealed, are of in tere st in th eir own right.
I n th e course of this work it becam e ap p aren t th a t th e mode sum very close to th e source required so m any (evanescent) m odes th a t the preparatory m ode com putations
14Section (4.2.5). ls Section(5.3.3) and (6.5). 16Section (6.1.4). 17Section (6.2.1). I#Section(6,2.3). I9Section(6.4).
20T he correspondence betw een an infinitely deep and a corresponding deep bounded waveguide is the su b ject of Section (6.4).
