Applied Acoustics 27 (1989) 163-201

Ocean Propagation Models

C. H. H a r r i s o n
YARD, ScientificHouse, 40--44 Coombe Road, New Malden, Surrey KT3 4QF, UK

A BS TRA C T
Four approaches to underwater sound propagation modelling are reviewed
including rays, normal modes, Green's function integral, and parabolic
equation. Specific programs are discussed. Particular reference is made to the
problems in running models and the applicability in various regimes,
especially deep-water environments.

1 INTRODUCTION
The wave equation and the Helmholtz equation occur in many branches of
physics, and the basic types of solution are always the same. However, the
problems in different fields vary tremendously because of emphasis on
different types of source, receiver, environment, geometry and so on.
Whereas propagation of light and radar in air is relatively straightforward,
underwater sound propagation can be extremely complex and calculations
may require several hours of computation on a CRAY-size computer. In the
underwater environment there are many different types of problem, some far
simpler than others. These include detecting distant ships or submarines,
echo sounding, short baseline location, seismic surveying and acoustic
tomography.
A typical problem of interest in this paper is to calculate intensity (and
phase) from a receiver somewhere in the water column separated by
distances of 1-100km in water of depth 100-5000m. The water is usually
assumed to be stratified (i.e. depth-dependent sound speed with constant
163
Applied Acoustics 0003-682X/89/$03-50 O 1989 ElsevierSciencePublishers Ltd, England.
Printed in Great Britain

C. H. Harrison

164

density), the sea surface is treated as a perfect pressure-release reflector, and

the sea bed may be treated as a reflector or multiple-layered refracting
medium with arbitrary density, sometimes demanding inclusion of shear
wave effects and surface scattering. A number of models are capable of
tackling range-dependent environments in which the bathymetry and sound
speed structure are allowed to vary relatively slowly in range. A few can
handle three-dimensional effects to a limited extent.
The aim of this paper is to provide a review of some of these programs for
the non-specialist in underwater acoustics. The examples are taken from
various underwater environments, and are woven into the discussion which
emphasises some of the practical problems in running the models. Some
rules-of-thumb are suggested for deciding which model is appropriate.
The review also refers to some comparisons, tests and a few unusual
applications. The models considered are some of the well known ones
developed with physical understanding in mind and minimal approximation
or empiricism. These are: GRASS (ray); SNAP and SUPERSNAP (normal
mode); SAFARI and FFP (Green's function); and PAREQ and IFD
(parabolic equation).

2 REGIMES A N D E N V I R O N M E N T S
It is difficult to appreciate the complexity of propagation modelling without
understanding the many environments and propagation regimes. These are
distinguished partly according to the various mathematical approaches and
partly according to oceanography and geophysics; many references discuss
the subject.l - 5 On top of this there are variants of these environments and
other applications of the usual models such as propagation from air to
water,6-a propagation in air, 9 seismic modelling, 1 ultrasound, 1 and radio
propagation over undulating terrain (the multiple-knife-edge diffraction
problem which is often tackled fairly crudely l~). Despite the complications
of diffraction and wave treatment the phenomena that need to be modelled
can be understood very simply by considering the corresponding ray paths.
2.1 Shallow water

The term 'shallow water' is often taken as synonymous with the Continental
Shelf where it has special strategic significance. For shipping-line frequencies of less than a few hundred hertz, where the wavelength is a
significant proportion of the water depth, shallow water behaves like a
waveguide in which propagation is dominated by bottom reflection. Ray
paths may be relatively steep, and consequently relatively straight, in the

Ocean propagation models

165

water column so that the influence of the water column velocity profile is
minimal. On the other hand, the propagation may be strongly influenced by
the structure of the sediment and the underlying rock since sudden changes
in velocity and density and the inevitable gradual increase in velocity with
depth provide strong reflections, particularly if there is a critical angle, and
upward refracted paths with potentially low losses.
Absorption in the sediment increases with frequency and typical values 12
are between 0.05 and 0.5dB/m/kHz, i.e. at 100Hz the loss might reach
l0 dB after a 1 km path length. By contrast, the losses in water are extremely
low, although they still rise with frequency. At 100 Hz a 10 dB loss is reached
after l0 000 km. The corresponding ranges for 10 kHz, though, are 10 m for
sediment and 10 km for water.
An additional factor in shallow water is the possibility of a solid rock sea
bed which supports shear waves. Under certain conditions the boundary
between the solid and fluid layer (whether water or sediment) may exhibit
interface waves, the effects of which are significant for sources and receivers
near the bottom. 13 A common phenomenon is the result of the fact that the
shear critical (grazing) angle, if it exists, is always smaller than the
compression wave critical angle, so that there is a tendency for energy to leak
into the bottom by generation of downward-propagating shear waves. This
effect can be screened out or toned down by an intermediate layer of fluid
sediment between the rock and water.
2.2 Deep water

In deep water propagation can be very complex because there are so many
different types of path. The sound speed increases roughly linearly with
depth, temperature and salinity. The salinity in the oceans usually plays a
minor role, although in the polar regions, in landlocked seas such as the Baltic,
and in the vicinity of some ocean fronts and eddies, there are significant
effects. In the North Atlantic in winter there is typically a well mixed warm
surface layer, several hundreds of metres thick which, when superimposed
on the depth effect, causes a minimum in sound speed at about 1000 m with a
maximum at several hundred metres. Above this maximum there is a surface
duct where rays from a shallow source can be trapped by upward refraction
and surface reflection.
Below the maximum is a deep sound channel centred on 1000 m, bounded
entirely by refraction. However, to make use of this regime both source and
receiver need to be at depths below the maximum. Steeper rays from a nearsurface source travel more deeply, and the roughly linear increase of velocity
with depth causes upward refraction. There is often a slight focusing effect so
that the rays that nearly hit the bottom bunch up near the surface at around

166

C. H, Harrison

60 km. This region of sudden high intensity is known as a 'convergence zone'

and it repeats at multiples of about 60 km.
Two extremely important phenomena resulting from refraction are
caustics, where initially adjacent rays converge, and shadow zones, where
adjacent rays diverge. In the velocity profile just described a shallow source
might form a wedge-shaped shadow extending from the velocity maximum,
which defines the bottom of the surface duct, to the downward-refracted
ray below. Conversely, a shallow receiver would not be able to see a source in
this shadow zone via the direct refracted path. However, a downwardlooking or omnidirectional beam receiver would also be sensitive to bottom
reflections (bottom bounce) or bottom refracted paths which would fill the
shadow in. Consequently, the nature of the sea bed can be extremely
important for near-surface sources and receivers.
The Atlantic in summer has a more gradual change in temperature with
depth, so that the profile is closer to a parabola. In all seasons very shallow
surface ducts (30 m or so) may form by mixing during the day, and these may
have significant effects above a few kHz where the duct thickness is greater
than a few wavelengths. From the modelling point of view it should be
remembered that the surface mixing that gives rise to the isothermal layer
constituting the duct is caused by wind and wave action, which necessarily
makes the surface rough. At frequencies of several kHz (wavelength less than
1 m), scattering effects are certainly significant and difficult to model from
first principles since the mechanisms are a complicated mix of surface shape
effects and surface-modulated volume scattering by entrained bubbles. 1'~4
A topic of increasing interest is Arctic and under-ice propagation. 15'16
Here the profile is virtually linear, giving rise to an upward-refracting deep
surface duct. Two of the complications added by ice cover are the roughness
of the lower surface and the effects of shear waves propagating in the ice.

2.3 Range dependence

An improvement over the assumption of a stratified medium is made by
including changes of velocity profile and water depth (bathymetry) with
range. Slow changes in bathymetry may deflect bottom-reflected paths into
the sound channel and vice versa. Relatively abrupt changes are
encountered with ocean fronts where there may be significant changes over a
horizontal distance of 100 m. As will be seen, there are a number of models
that can handle range dependence, but not many that can cope with abrupt
changes. In any case, it is extremely hard to validate range-dependent
models against real data, and it is also difficult to devise analytical
benchmarks.

Ocean propagation models

167

2.4 Three-dimensional effects

Clearly, refraction effects are also possible in the horizontal plane, especially
in the vicinity of ocean fronts. Increasing computer power has made it
possible, for instance, to investigate the horizontal effects of an eddy passing
across a fixed source 17 (made tangible by cine film of a very large number of
runs showing intensity at a fixed depth in the horizontal plane). A less
obvious horizontal bending, shadowing and focusing effect is caused by
propagation across a slope, each bottom reflection causing a slight
horizontal deflection.18-2o In fact, the horizontal curvature can be viewed
as stemming from the differing phase velocities of the vertical normal modes
in the varying-depth waveguide. 21
Although there are clearly three-dimensional effects which are important
for some of the time, propagation modelling in this domain is still in its
infancy. A lot of recent effort has gone into deriving exact solutions for the
wedge, 22 and a conical seamount, 2a and confirming these experimentally. 24
There is no truly general three-dimensional propagation model although
there are some close approaches. 25
3 MODELLING TREATMENTS
The two most widely used concepts in understanding propagation problems
are rays and normal modes. Naturally, the approaches to modelling have
followed these concepts, but emphasis in the past 20 years has been on
computer-intensive wave treatments which are broader in scope than the
normal mode approach. Nevertheless, in reconciling results from different
types of model one still has to resort to simple ideas and calculations using
rays, modes, images, flux and so on.
It is natural to use whichever approach is computationally most efficient
for the circumstances, and so in shallow water or in a duct where there are
vast numbers of reflected rays it is most useful to think in terms of modes. 26
Conversely, at short range or in deep water there may be only a few
'eigenrays' connecting the source and receiver but a very large number of
modes, so it is advantageous to think in terms of rays.
One can group propagation models into seven types as follows, although
the first four are in more c o m m o n usage than the subsequent two, and the
last is much more special-purpose.
(1)
(2)
(3)

ray tracing;
normal mode;
Green's function types;

C. H, Harrison

168

(4)
(5)
(6)
(7)

parabolic equation types;

coupled-mode types;
finite element types;
empirical.

The first six types have been reviewed by De Santo 27 and the first four by
Jensen. 2s This review will restrict itself to the first four types although some
references are given for coupled mode, 29-33 finite element 27'34'35 and
empirical 1,36 - 3 8 models.
Examples of computer codes are given below: a more complete list is given
in references 39 - 4 2 :
ray,
mode,
Green's function,
PE,

GRASS, 43'.4 PLRAy45;

SNAP, 46 SUPERSNAp47;
SAFARI, 4s-s FFP51,s2;
PAREQ, 53.54 IFD. 55'56

There are other models such as FACT 57 and R A Y M O D E ss that do use ray
concepts, although transmission loss is calculated without going through the
motion of tracing the rays. These models are designed to be used 'hands-off"
and so mode calculations are inserted automatically as appropriate, source
and receiver may be swapped, etc. Therefore for the research user it is often
difficult to control what calculations are executed.
A brief description of the four approaches is given below, but more detail
can be found in Refs 2-4, 27-28 and 59-61. All four can be derived from the
inhomogeneous Helmholtz equation, i.e. the wave equation with a point
harmonic source,
V % + k2(r)4, = - ~(r)
(1)
where k =

w/c is the

local wavenumber.

3.1 Ray tracing

In the limit of high frequencies the Helmholtz equation can be reduced to an
Eikonal equation which describes the path swept out by rays in three
dimensions. 6 In an inhomogeneous medium the changes in direction cosine
(0q fl, y) for an incremental step, ds, in the ray are given by39:

(1 dfl=
dy=

~-~[(1
ds[(1

oc
fl2)~y 2 ~c

Uy-y

oc1

Oc - o~flOc]
Ox[

)'fl-~z

Oc fl ~c7

(2)

Ocean propagation models

169

A convenient assumption in a marine environment (and several others,

such as air) is that the medium is horizontally stratified, and then one can
directly calculate range, r, and delay time, t, from
r = S cot 0 dz

t = S (c sin 0)- 1 dz

(3,4)

using Snelrs law, that csec 0 is constant for each ray. Also, the radius of
curvature of the ray, R, can be written in terms of velocity, its gradient and
the grazing angle at one depth as
R = c see O/(dc/dz)

(5)

The intensity, /, relative to the intensity at unit distance Io is usually

calculated from the horizontal spread of the rays. With an initial ray
elevation angle and sound speed 0o, Co and a local angle and sound speed 0, c,
the formula is
I = Ioco/(Cr sin 0 IOr/O0o I)
(6)
With the exception of the vicinity of focuses and caustics this gives
adequate representation of the loss. In these regions corrections can be
applied. '.3
There are a number of other useful ray or hybrid concepts such as ray
invariants,18- 21 and flux. 26'62'63 These are related to the WKB formula for
mode number 64 and the adiabatic approximation, 21'32 and are useful when
the medium is range-dependent. Although these approaches have not
formed the basis of computer codes they can provide insight through
analytical or numerical calculations. Simple calculations may also be based
on the method of images 65 in some circumstances. The concept of 'fuzzy'
rays has also been developed recently. 66
Frequency-dependent effects are sometimes included in ray-tracing
models, but without detailed corrections to the phase of bottom reflections,
ray tracing cannot compete with normal-mode models in shallow water. A
rule of thumb is that ray tracing can be used when the water depth or the
thickness of a refracting duct is many wavelengths. One disadvantage of ray
tracing is that many rays need to be computed to provide a reasonable
spread at the receiver end. Another is that ray computations need to be made
at all ranges out to the specified receiver. This is not so with normal mode
models. An advantage is that rays can easily be traced through varying
velocity profiles and over undulating sea bed. Directionality of the source
can be inserted, in principle, by selecting particular initial ray angles and
weighting accordingly.
3.2 Normal mode
The Helmholtz equation can be separated in range and depth assuming
cylindrical symmetry and vertical stratification. The solution of the resulting

('. H. Harrison

170

one-dimensional equation in depth can be expressed as the sum of the

discrete normal modes (i.e. the solution of a homogeneous Helmholtz
equation) and one or more branch cut integrals. 27 The acoustic pressure is
given by
i ~ , q~.(zs)~.(zr)H~o1)(K.r)

(7)

where q~, are the normal modes evaluated at the source and receiver depths zs
and z,, and K. are the eigenvalues. The discrete sum represents loss free
modes equivalent to up- and down-going rays travelling at a well defined
angle given by the horizontal wavenumber K, = k(z) cos O(z). The branch cut
integrals represent lossy modes travelling at lower wavenumbers (steeper
angles), and the contributions that are only effective at short range where the
up- and down-going waves do not balance properly yet.
The normal modes are calculated by choosing trial values of K. and using
numerical techniques such as Runge-Kutta iteratively, until the boundary
conditions are matched at depths zero and infinity with the correct number
of zero crossings. 6v These methods are adequate for shallow water in which
there are a small number of modes, and they are the basis of the model
SNAP. In deep water the far larger number of modes is a problem in itself,
but the close packing in wavenumber coupled with the refracting ducts
means that low-order modes are likely to contain upward- and downwarddecaying portions which lead to stability problems with standard shooting
methods such as Runge-Kutta. Porter and Reiss 4v'6a have used a finite
difference approach which is adopted in SUPERSNAP to formulate the
mode calculation as an algebraic eigenvalue problem which is more stable.
Ferla et al. report some high-frequency normal mode calculations in deep
w a t e r . 69

Strictly speaking, because the normal mode method is based on

separation of range and depth variables, it cannot be used unless the medium
is horizontally stratified. Nevertheless, if the range variations in velocity
profile and depth are only slight, the energy in each mode remains constant
as the mode propagates (the adiabatic approximationS2'ss). In fact, the rays
corresponding to a particular mode become steeper when the water becomes
shallower because of the reflections from the sloping bottom. These steeper
rays in shallower water correspond to exactly the same mode as the original
one in deep water, although the mode shape will have squashed vertically to
fit into the water column. Weak variations of velocity profile and water
depth are incorporated into S U P E R S N A P via the adiabatic approximation.
Various other effects have been included in normal-mode models. The
effect of different bottom types may be incorporated through the velocity

Ocean propagation models

171

profile, density and absorption. 67 Shear wave effects, 7 and surface or

bottom roughness 7a can also be treated.
Other variants of the normal-mode programs include those which
compute the complete spectrum of modes. This is necessary under certain
circumstances, for instance when the last mode has passed cut-off (discretemode programs would give zero result) or when there are downwardpropagating shear waves, weak reflections or when the receiver is close to the
source. The extra computation involves evaluation of a branch cut
integral.27.72-74 It is probably more straightforward to use the Fast Field
Program in this case since it routinely evaluates the whole spectrum of
modes.

3.3 Green's function solutions

The complete solution of the inhomogeneous Helmholtz equation for a
stratified medium can be written as the Hankel transform of the vertical
Green's function

dp(r,z,) = fo G(K, z,, zs) Jo(Kr)KdK

(8)

This is the basis of the Fast Field Program (FFP) 51'52 and SAFARI. 48 - s0
In both programs the bulk of the computation time is taken up with
calculation of the Green's function G (as a function of horizontal
wavenumber K). The Hankel transform is approximated by a Fourier
transform which is implemented by F F T algorithm. The original F F P used
the Thomson-Haskell method at each K to evaluate G by matching the
solutions at each layer boundary and the source. S A F A R I starts with the
same equations but uses a global matrix scheme to solve for all K and all
layers simultaneously.
The Green's function, as a function of K for fixed source and receiver
depth, is an intermediate output of the programs, and it is extremely useful
not only for checking performance but for providing insight into the physics.
The discrete and virtual modes are shown in Fig. l(a), and the corresponding
transmission loss is shown in Fig. l(b). The distinction between the Green's
function and the normal modes can be seen by analogy with a violin string
stretched across the water column. The normal modes are the shapes of the
violin string at its many resonances (the violin frequency corresponds to the
horizontal wavenumber). Between resonances the amplitude is zero, and so
the sequence is discrete. The Green's functions are the shapes of the violin
string when driven by a harmonic source at some particular position along
the length (i.e. the source depth). These exist for any frequency of the
harmonic source, although there will be a large amplitude at each resonance.

C. H. Harrison

172

SAFARI

20.0
F

* 15.0

50.OH~

SD =

50.OH

RD:

IO0.OH

10.0
C3
Z
I%
(._b
LLI
F--

5.0

z
'~

0,0
314

4734.2

9154.2

13574.2

17994.2

2241

.2

( 1 O* *-- 5)
(a)

50

SAFARI

-TLRA

613 0

70

u'l

80

O
J

90

1013

1113

13
0.0

15.O
RANGE

313.O

45.0

60.0

tKM;

(b)
Fig. 1. (a) SAFARI discrete and continuous mode spectrum. (b) SAFARI transmission loss.

The effect o f including bottom losses or any other losses is to broaden and
shorten the resonance peaks in a calculable way. When the broadening is
severe the modes are referred to as 'virtual' modes 75-77 and the modes
decay exponentially with range (as is apparent from considering the Fourier
transform o f an exponential). Similarly, the method can handle mode cutoffs, shear waves in solid layers (sediment, seabed, or ice) and near field as

Ocean propagation models

173

opposed to far field. S A F A R I (and FFP) is a very powerful program for

research because it gives an accurate solution of the wave equation for mixed
liquids and solids (including compression, shear and interface waves13), at
any range including the near field, above, below or near mode cut-offs.
S A F A R I has been implemented by Schmidt on an array processor, making
calculations with multiple frequencies, sources and receivers feasible. This
has led to the ability to produce synthetic seismograms and contour plots of
ultrasonic beams? Its two shortcomings are that it can (so far) only deal
with horizontally stratified media, and its computation time is rather long
despite the efficiency of the code.

3.4 Parabolic equation

The parabolic equation takes a different starting point from SAFARI, F F P
and the normal-mode programs. For waves travelling predominantly within
a small range of angles (not necessarily horizontal) one can approximate the
elliptic wave equation to a parabolic equation by taking out the main
oscillating part of the solution in a function S(r). Thus the velocity potential
~b can be written in terms of a slowly varying function, , of range and
depth2S,54:

q~ = (r, z)S(r)

(9)

and the parabolic equation is

02
Oz---f + 2ik o ~r + k2(n 2 - 1) = 0

(10)

where the wavenumber k = kon has been written in terms of an arbitrarily

chosen constant ko (the 'reference wavenumber', related to the 'reference
sound speed' by k o = co~c) and the refractive index n = Co/C which is assumed
to be slowly varying. It is evident from eqn (10) that if~, is completely defined
over a vertical line at some given range then and d2/Oz 2 are known, and
consequently de~Or is known. Therefore, ~k can be calculated for all z at the
next range step r + dr, and so on. This 'marching' solution allows one to
start with a known wavefront at some point and follow its horizontal
progress as it diffracts through the medium and around the humps in the sea
bed. Clearly, this approach can cope with horizontal variations in velocity
profile and variations in depth. 78 A limitation is that it needs to be started up
some way away from the source. For this purpose a normal-mode program
is often used. An alternative is to assume the initial amplitude distribution to
be a Gaussian function of depth, centred on the true source location. 2s'55
Equation (I0) can be converted into manageable form by taking vertical

174

C. H. Harrison

Fourier transforms to give the split-step-Fourier algorithm s4 which is used

in PAREQ:
~O(r+ Ar, z) = exp (iko(n 2 - 1)Ar/2) ~" -x {exp (-iArs2/2ko)

~,~(~(r,z))} (11)

The process of recalculating ~ at successive vertical planes is almost

identical to the propagation of Huygen's wavelets, and in fact the above
formula can be derived from the Kirchhoffapproximation treating the plane
at r as a diffraction screen. The various diffraction terms can be identified
here simply by taking the Fourier transform of both sides of eqn (11) and
writing the right-hand side as a convolution. Since the inverse FT of
exp ( - iArs2/2ko) is (ko/iAr2n) 1/2 exp (ikoz2/2Ar), the result is

ql(r + Ar, z)= (\iA---r~)

k

1/2

exp(iko(n 2 - 1)Ar/2)

x fo~ ~,(r,z')exp{[iko(z
J--

z')2]/2Ar} dz'

oO

Given the field on the vertical plane at r the second term in the integral is
precisely the Fresnel diffraction term in passing from an arbitrary point (r, z')
on the first screen to the observation point on the second at (r + Ar, z). The
exponential term outside the integral corrects the phase already included in
S(r) (i.e. kor, see eqn (9)) by adding in the phase due to the slowly varying
refractive index but taking out kor. The result, koAr(n-1) is equal to
koAr(n 2 - 1)/2 to first order since n is always very close to unity. The
multiplier turns into the usual ko/2niAr of the Kirchhoff approximation
when allowance is made for the other dimension of the diffraction screen
(horizontal, out of the plane of propagation) which supplies an extra
(ko/2niAr) x/2 term.
The method has been extended to wider angles (40-60 ) and stronger
bottom interaction by using finite difference methods, and one such model is
Implicit Finite Difference (IFD). ~9 There have been developments in which
density is treated more comprehensively.8 A number of operator techniques
for manipulating the parabolic equation have been presented at Yale. al
Recent extensions include calculating the field in a uniform a2 or sheared a3
current, and close approaches to three-dimensional modelling. 25'a4- 86

3.5 Model comparisons

There seems to be a law of nature that says that modellers do not go to sea
and experimentalists avoid using advanced models. Therefore, empirical
propagation laws tend to be extremely simple, whereas advanced models go
largely unverified against trials data. There are many good reasons for this.

Ocean propagation models

175

One is that transmission loss is very sensitive to environmental parameters,

and measuring all the quantities in the sea trial demanded by the model and
also guaranteeing to fulfil all the conditions of validity of the model in the
trial is almost always impossible exactly. Therefore honest comparisons are
difficult, although by varying some of the unknown parameters one can
probably obtain a reasonable fit. Nevertheless, an area that does have scope
for controlled tests is the reconciliation of existing models with themselves
and other analytical solutions. Indeed, it is important to reach a consensus
and map out regions of validity before attempting to make comparisons
with trials data.
Inter-model comparisons are occasionally made in the literature. Jensen
e t a / . 2a'sT'S8 have compared SNAP, FFP and PAREQ. Gilbert et aL 29 have
compared a normal mode model (COMODE) with a coupled mode model
(CUPYL) and IFD. Stickler 73 has made comparisons between FFP and
normal-mode solutions that include various branch cut integrals. A number
of parabolic equation corrections have been compared in Refs 89-93. Some
comparisons have been made by Tolstoy 92 between the parabolic equation
and ray tracing in a focusing environment. A methodology for comparison
and choice in models is given by De Santo 27 (pp. 121-34).
There have also been a number of workshops and specialc0nference
sessions based on model comparisons. The AESD workshop 94 compared
some wave treatment models, including FFP, discrete and continuous
normal-mode models and the parabolic equation. Another workshop at
N O R D A 95 concentrated on parabolic equation methods. More recently,
there have been the ASA sessions organised by Felsen at Anaheim 96 on
benchmarks, and on tests against some range-dependent and threedimensional analytical solutions. 9~

4 PRACTICAL LIMITATIONS A N D T U N I N G
Having briefly reviewed the mode of operation of the models it is useful to
look at some of the more practical limitations. These depend on validity in
the given regime or environment, computation time and very often a lot of
fine tuning and user experience. As well as the oceanographic and geoacoustic inputs, which are themselves often difficult to define, most models
have numerical inputs which need to be chosen by experience or trial and
error. For instance, ray-tracing models require a ray density at the source;
parabolic equation models require a reference sound speed and an
understanding of the angle limits and constraints on medium variability;
S A F A R I requires various compromises to be made in numerical
integration. One can usually devise simple formulae as an aid to choosing a

176

C. H. Harrison

starting point, but there is no substitute for repeating runs with altered
numerical inputs until convergence is reached.
The following discussions are restricted to GRASS (ray), SNAP and
S U P E R S N A P (normal mode), S A F A R I (Green's function) and P A R E Q
and I F D (parabolic equation). The comments are aimed at deep-water
environments with source and receiver within a few hundred metres of the
surface.
4.1 GRASS
It is obvious that the chosen ray density at the source must provide adequate
ray coverage at distant points of interest. However, it is not always easy to
predict what that density should be without a trial run. Ray traces with
coarse and fine ray densities for an identical shallow surface duct are shown
in Figs 2(a) and (b). The apparent size of the wedge-shaped shadow zone is
quite different in the two cases. To a certain extent the worries are illusory
because these rays (surface-reflected) do not contain much energy, but the
incoherently added intensity (Fig. 2(c)) shows steps corresponding to the
wide ray spacing in range, even for the high-density example. This kind of
problem arises when the ray location at a distance is very sensitive to initial
ray angle, i.e. lar/aOol is large, but this is exactly the condition that the
intensity contribution is low. Whether or not the result is important depends
on whether or not there are stronger contributions from elsewhere.
The obvious solution is to carry on increasing the ray density, but one
penalty is increased numerical error, and another more practical limitation
is computation time, which is roughly proportional to total range covered
and number of rays.
Although GRASS is essentially a high-frequency approximation, the
'coherent' option sums rays with regard to their phase and is capable of
showing some frequency dependence other than simple absorption effects.
In principle, a ray treatment can handle bottom reflections and bottomrefracted paths, but G R A S S cannot handle both correctly. One choice is to
use a compromise reflection loss curve which attempts to cover both cases,
but this cannot possibly treat the horizontal offset associated with the
refracted path correctly. Another is only to use the model with the true
reflection loss curve in the case where refracted paths are non-existent or
extremely weak. The converse case where refracted paths dominate over
reflected paths cannot be handled, despite the fact that refracted rays can be
traced. This is because the sediment layer would have to be treated as part of
the water column, and G R A S S does not have the facility to include
absorbing layers in the water.
In a duct ray treatments are valid as long as the duct supports a

Ocean propagation models

Ray t r a c e

(a)

177

ZS = 3 8 m

O,

20-

. ~.60-

801

100

,
1

,
3

Ray trace

(b)

,
4

. . . .

5
6
Range (km)

10

10

ZS:38m

.C:
1
0

Range (kin)

(c)

50"

607o--

8o j

110

\
,

210

'

4.'0

610

8"0

10-0

Fig. 2. (a) Coarse ray trace: 1 ray/degree; 16 spread. (b) Fine ray trace: 100 ray/degree;
1 spread. (c) Transmission loss showing steps due to finite but small ray spacing.

C. H. Harrison

178

substantial number of modes (i.e. definitely greater than one, but not
necessarily an extremely large number). Thus, treating the duct crudely as a
rectangular well of height h, the mode number is related to ray angle by
2H sin 0 = n2. For weak ducts with velocity contrast, Ac,
sin 0 = (1

c2/(c -[- At)2) 1/2 "~ (2Ac/c) 1/2

(12)

and a rough formula for the cut-off frequency (i.e. n = 1) is

f 2 = c3/8h2Ac

(13)

For instance, ray treatment in a surface duct with Ac = 4 m/s and h = 500 m
requires frequencies to be well above 20 Hz. In fact, since the duct crosssection is closer to linear than rectangular, the values of h and kc should be
reduced to compensate, so that the cut-off frequency is two or three times
higher. Better approximations to the formula may be made using the WKB
phase integral 1 for n. The cut-off frequency for the complete water column
is usually far lower, and for typical frequencies ray treatments are valid in
deep water (provided that some allowance is made for Lloyd's mirror, the
bottom dipole effect, etc. as appropriate). The cut-off frequency is now
roughly
f2 = [4H2(co 2 _ c12)] -1
(14)
where H is the total depth and the larger velocity spread is defined by the
limits c o and c 1. For Co -- 1500 m/s, c 1 = 2000 m/s and H = 4000 m the cut-off
frequency is 0.3 Hz.
Above the mode cut-off the dependence of intensity on mode number is
weak, because although the number of modes increases in proportion to
frequency the wavenumber in the denominator of the mode sum formula
(originating from the square of the Hankel function Ho(Knr ) in eqn (7)) also
increases in proportion to frequency. So, apart from some saw-tooth effects
with very low mode number the response is flat. A more important effect for
sources and receivers removed from the centre of the duct is the frequencydependent reduction in mode amplitude near the boundaries. Crudely, the
intensity for one mode is proportional to sin 2 ~.zs sin 2 ~.z,, and with both
source and receiver near the surface this reduces to ~,4.z~z,Z, where (in terms of
the horizontal wavenumber, K,, and the wavenumber in the medium, k)
V. = (k2 - Kz.) 1/2 = (2rtf/c)sin 0.

(15)

Thus intensity will rise with the square of frequency until saturation when
sin 2 ?,z approaches 1 for most modes. This can only be modelled with
G R A S S by a very large number of coherent rays.
A clear indication of the importance of modal effects at frequencies as
high as 5 kHz (in typical shallow surface ducts) is shown by the comparison
of Fig. 3(a) (an I F D intensity contour run on a C R A Y II) with Fig. 3(b) (ray

Ocean propagation models

o r

24

5000.OOHZ,

I~.~o

SO.

5.00M

~~~--~'--~-

~
.~

179

CONDR

..,

"

.,

.,
.o,I

,.,.
,
.

.'.

~
.

OI

;:~*':
o

.:

.:
o

(B)

.;

48

, v 0

~ 0 0

72

I--

o~
Ld
CD

~~.,,~I~I~,
II~'~
.~,~
~o".~.
~,"'~,''
~,.

..

o.~

.~

96

~ .',~
,
,~ ',o

.-~

<3

:o

120
0.0

2.0

.i .0

RANGE

6.0

8.0

I0

(KM)

(a)
RAY

t'~/

TRACE

PANCE

ZS-SM

(KM)

(b)
Fig.

3.

(a) I F D at 5 kHz in a surface duct, showing leakage in the form of downward beams
emerging at every cycle. (b) Corresponding ray trace.

180

C. H. Harrison

plot) for a surface duct. There are strong similarities in the shadow zone
boundaries, but the sharp-edged shadow envisaged by ray exponents is
clearly a myth. This area is partially covered by duct leakage which clearly
includes downward-propagating beams repeating at the cycle distance and
emerging from the duct just before the ray turning points at ranges of 1"3, 4.6
and 7.6 km.
4.2 S N A P and S U P E R S N A P
Normal-mode models are relatively robust, and they have the potential
advantage of being able to calculate intensity at any range or depth without
making computations at all intermediate ranges. Computation time consists
of the 'overhead' of calculating the modes, which increases in proportion to
the number of depth points, the number of modes and the uncertain but
usually small number of iterations, and then the mode summation for each
chosen range and depth point.
For most purposes the discrete mode sum is adequate. However, there are
several weaknesses in deep water apart from the large number of modes and
the convergence problems already mentioned. As a further extension of the
analogy with the violin string which was introduced in section 3.3, a normal
mode is a vertical standing wave caused by upward- and downward-going
waves interfering. The amplitudes of the two waves need to be comparable
to form a true standing wave, and this cannot be so until there have been at
least a few reflections or refractions to turn the down-going rays around.
Therefore the normal mode solution is incomplete for ranges shorter than
the cycle distance 9a for each mode (or ray family).
There are a number of other important effects that cannot be handled by
discrete normal mode models, essentially because the effects require lossy
modes to be included. Lossy modes form a continuous, rather than discrete
spectrum, and this requires evaluation of a branch line integral as well as the
sum. An important return at low frequencies in deep water at ranges before
the first bottom-refracted arrival is bottom reflection. Reflections from a
density discontinuity where velocity is more or less continuous are
necessarily lossy because the Rayleigh reflection coefficient is a constant
independent of angle, and there is no critical angle. The reflections therefore
do not feature in the usual mode sum. Nevertheless, these reflections may be
quite significant, as seen in the comparison of S U P E R S N A P with S A F A R I
(Fig. 4). This example has a water depth of about 1000 m. At ranges beyond
about 10km there are many deep bottom-refracted paths, and the good
agreement implies that a discrete sum is perfectly adequate. At shorter
ranges S A F A R I shows many rapidly interfering b o t t o m reflections which
are not seen in S U P E R S N A P .

Ocean propagation models

50.0

181

SUPERSNAP
I'LRAN,COH
SD30.OOM
RD40.OOM
F
50.OOHZ

SS024

&
n

70.0 i
'JI
0
._I

90.0

110.0
0.0

15 .0
RANGE

30.0

45.0

GO.O

(KM)

(a)

SAFARI-TLRA
SD30.OOM
RD40.00M
F
50.00HZ

50 0
60

70 0
Q

U')
O
.-1

80
90

I00
110

0
0.0

15.0
30.0
RANGE {KM)
(b)

45.0

60.0

Fig. 4. (a) SUPERSNAP and (b) S A F A R I showing the neglected short-range reflections
due to density differences at the sea bed.

182

C. H. Harrison

Shear effects have been included in normal mode models, 67 but in many
environments of interest lossy modes are again important, and these can be
handled more faithfully by SAFARI. For instance, the shear critical angle is
always smaller than the compression critical angle (since shear speed is less
than compression speed) so that the spectrum includes a set of decaying
modes, which is neglected by the discrete sum, but which may nonetheless be
significant.
The adiabatic approximation can be used for 'slowly varying' environments, but it is sometimes difficult to see what this means in any other way
than the statement that there is no mode coupling. In discussing applications
of ray invariants, the author of Ref. 21 looked at the equivalence between
adiabatic modes and ray angles and their joint conditions for validity in a
reflecting duct. Rays become steeper in shallow water in a reversible way,
provided that the reflecting surfaces are smooth and their shapes change
slowly. This is equivalent to the adiabatic approximation. It is clear that in
crossing a ridge, say, the process is not reversible if the critical angle is passed
at any stage; it is also clear that undulations in between bottom bounces will
produce a virtually random additional angle to the ray. Therefore the
condition is roughly that ray angles (at the sea bed) should always be
considerably greater than the slope of the sea bed, so that there are many
bounces within a small change of environment.
4.3 SAFARI

SAFARI is an exact solution of the Helmholtz equation (in a stratified

medium) for ranges greater that a wavelength or so. Compression and shear
velocities with separate absorptions and stepped density profiles can be
handled. Limitations are not fundamental since they stem from the problem
of resolution in wavenumber which can be alleviated (without changing the
code) given enough computation time and numerical accuracy.
The number of points in the Fourier transform, N, constrains the range,
re, and the usable portion of this is r, = rF. The usual Fourier transform
relation is

2nN=AKrF

N =fc-~ (1 c _ ~ )

(16, 17)

where AK, f, c 1 and c 2 are, respectively, the wavenumber resolution,

frequency, and lower and upper phase velocity assumed in calculating the
Green's function (as a function of horizontal wavenumber K). The CPU time
to carry out the FFT is relatively small, but calculation of the N Green's
function values for entry into the FFT increases in proportion to N.
Calculations become much more efficient for large numbers of receivers, 1

Ocean propagation models

183

because the Green's function has been matched at all the layer boundaries
and its form is known in each layer so that it is easily evaluated anywhere
without starting from scratch. Without increasing N the only way of
improving the resolution in deep water where there is a large number of lossfree modes (each corresponding to a spike in wavenumber, i.e. a pole in the
complex wavenumber plane) is to offset the line of integration from the real
wavenumber axis by adding a small imaginary part to K. This has the effect
of moving the integration line further from the poles (which are on or very
close to the real axis) so that the integral appears smoother; consequently
resolution with the same Nis improved. A criterion has been given 4s for the
offset which can be written in terms of dB per wavelength, J, as
J = 60cl/frr

(18)

Substituting for the maximum range this reduces to

(19)
A rule of thumb for having achieved a well sampled K-plot is that the
function should look like a modulated sine wave with a smooth envelope.
This is nearly always true for low K, but for high K a jagged envelope
demonstrates poor sampling.
This effect and the consequences for transmission loss are shown in Figs
5(a)--(f). Firstly a small offset results in jagged peaks (a) and a 'noisy'
transmission loss plot (b). The optimum offset case (c) shows a much neater
K-plot with a smooth envelope and the low order modes have been lumped
together. The loss plot (d) is relatively smooth but still rather noisy at long
range. A large offset results in a very smooth K-plot (e), but the penalty is that
the intensity begins to increase with range (f). In effect, the artificial insertion
of a small imaginary part to the wavenumber e damps the modal resonances,
but to retain the low loss of the original modes the formal mathematics 48
needs to compensate by amplifying the result by an exponential exp (~r) at
each range r. This delicate balance is upset if ~ becomes too large.
A complementary point is that the Green's function amplitude must tail
off to zero for high and low K to avoid superimposing Fourier transform
noise. This appears at first sight to be a reasonable approximation in the first
case, but as the function broadens with increasing offset the left-hand side
rises (relatively) so as to be appreciably 'chopped' by the end of the Fourier
transform. In fact, what has happened is that the area under the left-hand
side has remained constant and so has the area under the right-hand side
(low order modes), but the original spikiness of the modes gave the mistaken
impression that the low K Green's function amplitude was negligible.
The earlier equation for N makes it clear that for fixed computation time

C. H. Harrison

184

(a)

SAFARI

20.Q

SO=

40.0H7
50,OH

RD =

IOO.OM

15.0

C3

IO,O
z
c~
-5
Ld

5.O

O.O
4188

6938 8

9688.8

(b)

12438.8

15188.8

I ?938 8

(10"*-5)
SAFAR]-TLRA

50.0
60.0
70.0
80.0
90.0

I00.0
110.0
0.0

15.0
RANGE

Fig. 5.

30.0

45.0

GO.O

IKMI

The effect of altering the integration line offset from the real K-axis: (a), (b) small
offset; (c), (d) optimum offset; (e), (f) large offset.

(i.e. fixed N) the user may trade range for velocity contrast or frequency. At a
fixed frequency this usually means that he or she has to open up the phase
velocity limits Cl and c2 to accommodate the real environment. This puts a
restriction on the maximum range. When there is an abrupt density change
at the sea bed the wavenumber plot (integrand) often has significant
amplitudes for phase velocities well above the maximum velocity in the
sediment (i.e. wavenumbers much lower than 2rcf/cB), and a suitable phase
velocity must either be found by trial and error or set to infinity.
An interesting check for SAFARI in a non-trivial environment is given by
an analytical image calculation. Imagining source images in the multiple
surface and bottom image planes at depths -I-2nil +_ z, one can group the

Ocean propagation models

185

SAFARI

40.0

(c)
0

* 30.0

0
20.0
E3
Z

13~ 1 0 . 0

t----

0.0

4188.8

6938.8

9688.8

12438.8

15188.8

?938.8

(10"*-5)
(d)

SAFARI

50.0

SDRDF -

60.0

-TLRA

50.OOM
IO0.OOM
40.00HZ

70.0
o
ul

80.0

o
_3

90.0

100.0

110.0

.0

51.0
RANGE

30.0

45.0

60.0

(KM}

Fig. 5.--contd
sources according to the number of bottom interactions (in isovelocity
water). The true source and one image (in the surface) have no bottom
interactions and produce a Lloyd's mirror effect. The next four images (or
four rays) have one bottom interaction (but different numbers of surface
interactions). The next four have two bottom interactions, and so on. It is
assumed that the angle differences between each order of reflection will be so
great that very rapid oscillations will result, but that they can be neglected
because usually one order of reflection dominates for most ranges. Thus the
main point o f interest is the remaining relatively slow beats in the spatial
pattern. The transmission loss for each reflection order is easily shown to be
TL = - 1 0 l o g P
P = 16(R"/r') 2

[sin 2 (a,)sin 2 (at) + sin 2 (a~r)(1 -- sin2 (a,) -- sin 2 (at))]

(20)

186

C. H. Harrison
SAFARI

40.0

(e)
0

30.0

F =

40.OHZ

SD=

50.OM

RD=

~00. i

20.0
L--%
Z
El". I0.0
0

W
Z.

0.0

4~88.a

6938

9688.8

(0

12438.8

15188.8

7938.8

(I0"*-5)

SAFAR]

50.0

SDRDF

60.0

--rLRA
50.OOM
lO0.OOM
aO.OOHZ

70.0

,4

ul
c3

80.0

96.G

lO0.G

110.0
.C

15.0

30.0

RANGE

45.0

60 .0

IKMI

Fig. 5.--contd.
for n > 1, where
r '2

= (2nil) 2 + r 2

a, = 2knHz,/r'

(21)

a s = 2knHz,/r'
as, = k z : , / r '

If the bottom interface is marked only by a change in density p, the Rayleigh

reflection coefficient is
R = (p - 1)/(p + 1)
(22)
This is a constant for all angles, and Fig. 6 compares numerical evaluation
of eqn (20) for n = 1 and 2 and Lloyd's mirror with S A F A R I in the case
where R = 0.1. Agreement is extremely good, and the residual interference

Ocean propagation models

(a)

187

o
i

N.0

SD-

50,,

N-t

RD-

100-,,

[]

N-2

-4 0 0 0 ~

- lO.OdB

u3
If}
0
J

50H~

0
O3
,

vo

40

20

RANGE

(b)

SAFARI

50.0

G O .O

m
1:3

70.O

80.0
O
J

60

(KH)

SORDF

-TLRA

50.OOM
IO0.OOM
50.OOHZ

9oo

lO0.O
I|O.O
O.O

5.0

RANGE

Fig. 6.

30.0

45.0

60.0

(KMI

An image calculation (a) used as a diagnostic tool to compare with SAFARI (b).

can at least partly be accounted for by the neglected beating between paths
with different numbers of bottom reflections. Other examples have shown
that the image calculation is extremely useful for filling in the short range
reflections which are omitted by other models.
4.4 PAREQ and IFD
There is a considerable literature on the parabolic equation and its
shortcomings, s1.99-12 The PE has a fundamental elevation angle

188

C. H. Harrison

c~

_c)

oO.t~0

12.00

6.00

lB,O0

RANGE

~4.00

30.00

36.00

42,00

(KM)

(a)
PAREQ-CONDR
0.0

50.OOHZ,

SO-

(B)

15.00M

....."~.'.! "';'"

'

i'"" - "

'.~
I000.0

2000.0

3000.0

W
0

4000.0

5000.0
0.0

8 .0

12 .O

18 .0

RANGE

t
24 .0

J
30.0

J
36 .0

l
42 .0

(KMI

(b)
Fig. 7. Angle or phase errors in the parabolic equation, illustrated by the discrepancies in
range to the first bottom refracted arrival at the sea surface. (a) GRASS ray trace; (b) PAREQ.

Ocean propagation models

0.0~F-

50.OOHZ, SO.

189

IS.OOM

CONOR
"

800 .ol ~,~ '"-,9

.0

~4o00

"

":.,~

'

-~L\ ~ .,, .
L- 1600

(B)

"

i,
~.~

,''! ~-

0.0

9.0
RANGE

"
18.0
(KM)

27.0

36.0

45.0

(c)
0.O

50.OOHZ, SO-

IS.OOM

'~ ~ ~

.....

CONDR

:.

:('

(B)

1600.0
"l- 2400.0
p-

8_
W
0

4000~0

0.0

S.0
RANGE

18.0

27.0

36.

45.0

IKMI

(d)
Fig. 7.--comd

(c) I F D (N); (d) IFD (W). A strongly refracting sediment layer between 4000
and 5000 m, is assumed.

restriction regardless of its implementation because it is an approximation

to the Helmholtz equation. I F D has two implementations, 56 one due to
Tappert 53 (narrow angle) and one due to Claerbout la (wide angle), referred
to here as IFD (N) and I F D (W). PAREQ and IFD (N) are often quoted as
having a limit of order 20 , whereas IFD (W) has a limit of 40 . A weakness in
these models is that there is no sudden change in the output at these angles,
and contour and loss plots continue to look realistic to the uninitiated. This
is because the wide angle returns are effectively mapped into narrow angles

190

C. H. Harrison

rather than being ignored.I4 An environment with a sediment velocity of

1600 m/s and low absorption is all that is required to produce 20 ~ rays. After
a long enough range the steep rays will die out, of course, but at low
frequencies there are often significant returns remaining within the first
convergence zone range.
A ray trace in which rays are transmitted through the sea bed at 4000 m
and refracted from the sediment back to the water is shown in Fig. 7(a). A
repeated pattern of rays penetrating to the bottom of the sediment (5000 m)
is seen. At 50Hz, where the wavelength is 30m, this plot is reasonably
representative except that it says nothing about the intensity. Equivalent
intensity contour plots for PAREQ, I F D (N) and I F D (W) are shown in Figs
7(b), (c) and (d). Since the loss is quite low (0.15 dB/wavelength) the repeated
pattern still exists in each case. What is striking is the varying distances to the
first bottom-refracted arrival at the surface; these are 10, 16, 21 and 16 km
respectively. The grazing angle of the rays at the sea bed is 54 (with a
sediment velocity gradient of 1 m/s/m) so it is not surprising that it is difficult
to predict performance once the rough angle limit has been exceeded.
T h o m s o n & Wood14 have suggested a practical method for realising De
Santo's correction. 27 This consists of separating out each horizontal
wavenumber component by Fourier transforming in range and performing
a weighting and mapping to a new wavenumber before transforming back to
range. In effect, the large-angle components are separately set to larger
angles while the small-angle components are left untouched. An alternative
approach is taken by Tolstoy, 92 where a similar effect is achieved by altering
the velocity profile progressively away from the velocity minimum. Provided
that the velocity contrast is small it is possible to make corrections by
choosing the reference sound speed carefully; a 'natural' choice is described
in Ref. 105. Using the WKB approach this can be converted into a very
simple practical formula.106
There is an additional effect (which may be seen in Figs 7(b), (c) and (d) at
very short ranges when source and receiver are at different depths. The
marching solution results in zero intensity for Izs - z,[ > r tan a, where ct is a
constant angle depending on the algorithm and the range step length.
The standard versions of I F D and P A R E Q cannot handle density
changes in the sediment exactly. Instead PAREQ, for instance, emulates the
change by altering the velocity profile. However, it can easily be seen that a
boundary with a density discontinuity but continuous velocity gives a
constant reflection coefficient, independent of angle. Thus at short range
there will be many reflections at steep angles right up to 90 . This can only be
emulated if the two velocities are close in value and consequently the
reflection coefficient is very low. More rigorous work on inclusion of density
has been done by St Mary. 8

Ocean propagation models

191

A practical limitation to the PE is the cost of computation time. This rises

in proportion to the number of calculation points in depth and range. The
usual step sizes of 2/4 and 2/2 result in computation time being proportional
to f2. However, this is offset to a certain extent by the fact that the
calculation necessarily includes many receiver depths and is therefore
directly suitable for contouring. In the vertical there is, in principle, some
scope for increasing the step sizes since the spatial variation cannot be more
rapid than the vertical wavelength of the highest order mode (at least for
calculations starting beyond the cycle distance). Experimenting with step
size is often made difficult by the assumption elsewhere in the code (e.g. in the
Gaussian initialisation) that the step size is no greater than a quarter of a
wavelength.
5 REGIONS OF APPLICABILITY
Each model has a range of input parameters for which the model outputs are
reliable. There are two independent questions that need to be answered in
any application. One is, what are these limits for a particular model? The
other is, what are the required values of the inputs for the environment in
question? These 'fuzzy' areas may not overlap, and there may even be areas
where no model works. The input parameters of all models break down into
operational parameters (i.e. frequency, source depth, receiver depth and
range), environmental parameters (i.e. velocity profile, water depth, sediment
depth, absorption and velocity and their gradients in the sediment, density
etc.) and numerical tuning factors (e.g. reference sound speed, depth and
range increments in IFD; integration offset, velocity contrast versus FFT
size trade-off in SAFARI; and range and angle increments in GRASS).

5.1 Effective angles dictated by the environment

A clearer picture of the conditions under which the model has to operate is
necessary in order to proceed further, and this can be seen by attempting to
limit the number of parameters. The approach here is to take only three
parameters: range and frequency, with the (horizontally stratified)
environment drastically reduced to one parameter. These parameters define
a volume that must be covered by the model.
The single parameter used to describe the environment is the 'effective
angle', as described below. At a particular range the total energy from an
omnidirectional source is spread into a range of angles; Snelrs law and the
reflection coefficients dictate the relation between angle and the complete
velocity and density profile, but attenuation effects will eventually narrow
the angle down at long range. The range of angles is related to the number of

192

C. H. Harrison

remaining modes that have non-zero amplitudes at source and receiver

depth (as can be seen from the WKB solution) and is related to Weston's
invariants and flux formulation 26 which are useful concepts for rangedependent media.
The initial energy split between surface duct, water column and sediment
is determined by the angle at the source for the limiting ray in each duct, i.e.
0 = cos- 1 (cs/c)

(23)

where cs is the velocity at the source and c corresponds to the maximum

velocity in either the surface duct, the water column or the entire profile
including the sediment. From simple flux arguments or WKB-mode
formulae the distribution of energy in the ducts is in the proportions
O~. (0,,, - 0~). (Oh - Ow)

lid"

Hw

"(lib + Hw)

where the subscripts refer to the surface duct (d), water column (w) and
sediment (b), and H is the appropriate layer thickness.
The 'effective angle' is defined, rather loosely, as the steepest angle that
succeeds in making a significant contribution for the given source-receiver
combination at the given range. It is especially useful in defining the acoustic
effects of the environment because it can be compared directly with angle
limits of models or translated into wavenumber limits as appropriate. There
are, of course, many complications, particularly at low frequency, such as
near-boundary effects, which need to be borne in mind as well.
At a single frequency there is a curve in the r-O plane which defines the
angles 0elf that must be catered for at each range. Examples are shown in
Figs 8(a) and (b) for low and high frequency respectively. The region that
must be covered is above and to the left of the curve, i.e. for a given range r,
all 0 < [01 < Ocff(r).
With a dense sediment there will be reflections at angles up to 90 , but at
extremely short range these will be overpowered by the direct path. The
effective angle will therefore be smaller than 90 depending on the detailed
geometry, but steep angles will become significant at ranges comparable
with the water depth. F r o m then on the effective angle will fall to the critical
angle and level out after a range of about the cycle distance 9s (2Hcot 0c). If
the energy is spread more or less uniformly in angle each contribution will
drop out at a range where its total loss, RL, exceeds a certain level. The
relation between R L and range r is easily calculated by considering the
absorption in transitting the sediment, and for small grazing angles, 0,

R L = rfaO2/Hc '
(24)
where f, a, H and c' are frequency, absorption in dB per wavelength, water

Ocean propagation models

193

Range

~'J . . . . . . . . . . . . . . . . . . . . . . . .
o l a f

....

n/2[..."""'~"=2H cotO~

T-d
fu

(a)
o

Range

Oduct

0~
rJ2

(b)
Fig. 8. The regime defined by the environment in terms of a single parameter, the "effective
angle', versus range for (a) low and (b) high frequencies.

depth and velocity gradient in the sediment respectively. If the drop-out level
is taken as a constant for all contributions then the range is
r w. H c ' / f ~ O 2

(25)

and has already been calculated for mode stripping. 17 If there is a surface
duct (above cut-off) the effective angle reduces finally to a value defined by
the velocity contrast and depth of the duct.
Naturally, energy partitioning between surface duct, water column and
sediment-refracted arrival must also b e considered. At low frequencies
where the surface duct cannot support any modes its energy contribution is
zero so that the surface duct can be ignored. At higher frequencies, despite
the low spread of angles, there may be relatively high energy density in the
duct because of its limited depth. In extreme circumstances this may make it
possible to ignore the bottom returns despite their high angles. It is difficult
to invent a general rule, though, because the relative strengths of the various
returns change in an arbitrary manner. The best way to find out is to perform
a pilot run to see which returns dominate and where, and then to use a ray
trace to estimate the appropriate angles.

5.2 Anglescoveredby GRASS

GRASS can be run in two ways, one having a reflection loss table strictly
covering reflections and ignoring bottom-refracted paths, and the other

C. H. Harrison

194

Range
S3
S~

Range

S 3 Si

. ........... Parabolic
Equation

Sa

I
I

_~

S~

S2

Normal
Mode

..-"

Si

Equation I

Normal
Mode

GRASS

7
Fig. 9.

vr ~ { 1 -2c o s 0 ) ~l
(a)

~2

fl

(b)
Thc areas of coverage given by G R A S S {ray),normal-mode parabolicequation and
SAFARI models for (a) low and (b) high frequencies.

having a compromise table covering all upward-turning rays. Neither case

features in the low-frequency graph, Fig. 9(a), but the latter case occupies the
whole of the high-frequency graph, the only limitation being the practical
consideration of computation time and accuracy at the right-hand edge. The
strict reflection loss case which is shown in Fig. 9(b) is only valid for the case
in which bottom-reflected paths (total loss R L = rLO/2H) are stronger than
refracted paths (eqn (24)). This condition is
0 >>c'L/2f~

(26)

5.3 Angles covered by SNAP and SUPERSNAP

The discrete normal-mode solution is valid at ranges beyond the cycle
distance and frequencies above the cut-off of the entire water-sediment
column. Shooting methods such as SNAP suffer not only from excessive
computation time when there are many modes but also from difficulties with
convergence in the mode shapes (particularly low-order modes). Thus they
run into difficulties at high frequencies, but may also have trouble in deep
water at low frequencies. The convergence problems are alleviated by
SUPERSNAP although computation time will still limit performance.
There is no particular angle limit in the normal-mode approach (other than,
by definition, the critical angle) but, for a fixed computation time, frequency
could be increased by artificially restricting the number of modes (in
principle, the restriction could be to any given set, not necessarily the lowest
order modes). This is shown by the horizontal line in Fig. 9.

5.4 Angles covered by SAFARI

The coverage given by SAFARI cannot be shown once and for all in Fig. 9,
since there is a choice to be made by the user between range and velocity

Ocean propagation models

195

contrast (effective angle). The computation time is roughly proportional to

the number of points in the FFT, N and, putting cos 0~ff = c~/c2, the earlier
formula (eqn (17)) corresponds to
C P U ~: N

N : f c l 1rE(1 -- cos Oeff)

(27)

For fixed C P U the 0eff(r ) curve is therefore

0etf = cos-1 [1 - Ncl/frF]

(28)

This is shown by the dotted line in Fig. 9. The user can choose any point on
this dotted line to define 0ef f and rr, and so he has a choice of the rectangular
boxes of validity defined by S 1, S 2, $3 etc. However, he does not have a free
hand because the answers will be incorrect or suspect if the effective angle in
Fig. 9 does not cover the required angle shown in Fig. 8, so there is always a
limit on range.
A possible, although tedious, way ofextending the range of SAFARI is to
calculate loss for short extensions in range from non-zero starting points in a
piecewise fashion. Another way is to calculate short-range TL from a high
velocity contrast run, and then calculate long-range TL from a separate low
velocity contrast run. The overlap is usually well behaved.
5.5 Angles covered by PAREQ and IFD
The parabolic equation is valid for large ranges and angles below some limit
for the implementation regardless of frequency. The effect is a flat but hazy
cut-off, as shown by the solid line in Fig. 9. At short ranges there may be
phase errors, and at very short ranges there are guaranteed angle violations
unless the receiver is at the same depth as the source. The practical limitation
of computation time (and ultimately numerical accuracy) is shown by the
dotted line on the right.
5.6 General comments on coverage
At short range the only wave treatment to give reliable coverage is SAFARI.
This is complemented by and overlaps with SUPERSNAP at longer ranges.
At frequencies where SAFARI becomes too expensive to run there is no other
model that can handle short-range returns, and simple image calculations
may be able to fill the gap. In a strongly range-dependent environment where
neither SAFARI nor S U P E R S N A P is applicable there is a strong possibility
of a hole in coverage at high frequencies, which becomes a certainty at low
frequencies. One hopes that, in practice, the environment only requires small
angles that either GRASS or PE can handle.

196

c. H. Harrison

6 CONCLUSIONS
A review of four types of propagation model has been given including
G R A S S (rays), S N A P and S U P E R S N A P (normal mode), S A F A R I and
F F P (Green's function), and P A R E Q and I F D (parabolic equation).
References have been m a d e to various model comparisons. Brief
descriptions of the models and the environments to which they apply have
led to a discussion of their limitations in terms of validity, computation time
and tuning. Finally, an attempt has been m a d e to map out regions of
applicability and to present a way of choosing models in practice.

ACKNOWLEDGEMENT
A substantial part of the work on which this review is based was carried out
under contract to the Admiralty Research Establishment (UK).

REFERENCES
1. Urick, R. J., Principles of Underwater Sound. McGraw-Hill, New York, 1967.
2. Tolstoy, I. & Clay, C. S., Ocean Acoustics: Theory and Experiment in
Underwater Sound. McGraw-Hill, New York, 1966.
3. Clay, C. S. & Medwin, H., Acoustical Oceanography. John Wiley, New York,
1977.
4. Officer, C. B., Introduction to the Theory of Sound Transmission. McGraw-Hill,
New York, 1958.
5. Kinsler, L. E,, Frey, A. R., Coppens, A. B. & Sanders, J. V., Fundamentals of
Acoustics. John Wiley, New York, 1982.
6. Weinstein, M. S. & Henney, A. G., Wave solution for air-to-water sound
transmission. J. Acoust. Soc. Am., 37 (1965)899-901.
7. Urick, R. J., Noise signature of aircraft in level flight over a hydrophone in the
sea. J. Acoust. Soc. Am., 52 (1972) 993-9.
8. Chapman, D. M. F., The effects of seabed interaction on air-to-water sound
transmissions. J. Acoust. Soc. Am. 82 Suppl. 1 (1987) DDD2.
9. Harrison, C. H., Propagation in air using SAFARI (unpublished).
10. Schmidt, H., Modelling of pulse propagation in layered media using a new fast
field program. In Hybrid Formulation of Wave Propagation and Scattering, ed.
L. B. Feisen. Martinus Nijhoff, Dordrecht, 1984, pp. 337-56.
11. Meeks, M. L., Radar Propagation at Low Altitudes. Artech House, MA, 1982.
12. Hamilton, E. L., Geoacoustic modeling of the sea floor. J. Acoust. Soc. Am., 68
(1980) 1313--40.
13. Rauch, D., Experimental and theoretical studies of seismic interface waves in
coastal waters. In Bottom Interacting Ocean Acoustics, ed. W. A. Kuperman &
F. B. Jensen. Plenum PresS, New York, 1980, pp. 307-27.
14. Crowther, P. A., Acoustical scattering from near surface bubble layers. In

Ocean propagation models

15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.

30.
31.
32.
33.
34.
35.

197

Cavitation and Inhomogeneities in Underwater Acoustics, ed. W, Lauterborn.

Springer-Verlag, Berlin, 1980, pp. 194-204.
Arctic acoustics. J. Acoust. Soc. Am., 82, Suppl. 1 (1987) C,M.
Under-ice propagation. Z Aeoust. Soc. Am., 80, Suppl. 1 (1986) Y.
Kuperman, W. A., Porter, M. B. & Perkins, J. S., Three-dimensional
oceanographic acoustic modeling of complex environments. J. Acoust. Soc.
Am., 82, Suppl. ! (1987) S1.
Weston, D. E., Horizontal refraction in a three dimensional medium of
variable stratification. Proc. Phys. Soc. Lond., 78 (1961) 46-52.
Weston, D. E. Guided propagation in a slowly varying medium. Proc. Phys.
Soc. Lond., 73 (1959) 365-84.
Harrison, C. H. Three-dimensional ray paths in basins, troughs and near
seamounts by use of ray invariants. J. Acoust. Soc. Am., 62 (1977) 1382-8.
Harrison, C. H., Acoustic shadow zones in the horizontal plane. J. Acoust. Soc.
Am., 65 (1979) 56--61.
Buckingham, M. J., Theory of three-dimensional acoustic propagation in a
wedgelike ocean with a penetrable bottom. J. Acoust. Soc. Am., 82 (1987)
198-210.
Buckingham, M. J., Theory of acoustic propagation around a conical
seamount. J. Acoust. Soc. Am., 80 (1986) 265-77.
Glegg, S. A. L. & Yoon, J. R., Range-dependent cut-on frequency in a threedimensional ocean environment. J. Acoust. Soc. Am., Suppl. 1, 82 (1987) CCI.
Siegmann, W. L., Kriegsmann, G. A. & Lee, D., A wide-angle threedimensional parabolic wave-equation. J. Acoust. Soc. Am., 78 (1985) 659-64.
Weston, D. E., Rays, modes and flux. In Hybrid Formulation of Wave
Propagation and Scattering, ed. L. B. Felsen. Martinus Nijhoff, Dordrecht,
1984, pp. 47-60.
De Santo, J. A. Ocean Acoustics. Springer-Verlag, New York, 1979.
Jensen, F. B., Numerical models in underwater acoustics. In Hybrid
Formulation of Wave Propagation and Scattering, ed. L. B. Felsen. Martinus
Nijhoff, Dordrecht, 1984, pp. 295-335.
Gilbert, K. E., Evans, R. B., Chin-Bing, S. A., White, D., & Kuperman, W. A.,
Some new models for sound propagation in bottom-limited ocean
environments. In A coustics and the Sea Bed, ed. N. Pace. Bath University Press,
UK, 1983, pp. 243-50.
Evans, R. B., A coupled mode solution for acoustic propagation in a
waveguide with stepwise depth variations of a penetrable bottom. J. Acoust.
Soc. Am., 74 (1983) 188-95.
Rutherford, S. R. & Hawker, K. E., Consistent coupled mode theory of sound
propagation for a class of nonseparable problems. J. Acoust. Soc. Am., 70
(1981) 554-64.
Pierce, A. D. Extension of the method of normal modes to sound propagation
in an almost-stratified medium. J. Acoust. Soc. Am., 37 (1965) 19-27.
Milder, D. M., Ray and wave invariants for SOFAR channel propagation. J.
Acoust. Soc. Am., 46 (1969) 1259-63.
Kalinowski, A. J., Survey of finite element related techniques as applied to
acoustic propagation in the ocean--l. Finite element method and related
techniques. Shock and Vibration Digest, 11 (1979) 9-16.
Murphy, J. E. & Chin-Bing, S. A., The finite element method applied to ocean
acoustic propagation. J. Acoust. Soe. Am., 81, Suppl. 1 (1987) D1.

198

C. H. Harrison

36. Marsh, H. W. & Schulkin, M., Shallow water transmission. J. Acoust. Soc. Am.,
34 (1962) 863.
37. Lindop, P. M., The deduction of approximate values of sediment depth from
propagation loss measurements. In Ocean Seismo-Acoustics, ed. T. Akal
& J. M. Berkson. Plenum Press, New York, 1986, pp. 243-52.
38. Harrison, C. H. & Cousins, P. L., A study of propagation loss dependence on
sediment layer thickness using the fast field program. In Ocean SeismoAcoustics, ed. T. Akal & J. M. Berkson. Plenum Press, New York, 1986,
pp. 139-48.
39. Urick, R. J. Sound propagation in the sea. Peninsula Publishing, Los Altos, 1982.
40. Etter, P. C., Underwater acoustic modelling techniques. Shock and Vibration
Digest, 13 (1981) 11-20.
41. Etter, P. C., Underwater acoustic modelling techniques. Shock and Vibration
Digest, 16 (1984) 17-23.
42. Etter, P. C., Numerical modelling techniques in underwater acoustics, aT.
Acoust. Soc. Am., 82, Suppl. 1 (1987) UU1.
43. Cornyn, J. J., GRASS: a digital-computer ray-tracing and transmission-lossprediction system. Volume I--overall description. Nav. Res. Lab., Report No.
7621, 1973.
44. Cornyn, J, J., GRASS: a digital-computer ray-tracing and transmission-lossprediction system. Volume II--user's manual. Nav. Res. Lab., Report No.
7642, 1973.
45. Bartburger, C. L., P L R A Y - - a ray propagation loss program. Nav. Air. Dev.
Center, Report No. NADC-77296-30, 1978.
46. Jensen, F. B. & Ferla, M. C., SNAP: The SACLANTCEN normal-mode
acoustic propagation model. SACLANT Memo SM~121, 1979.
47. Porter, M. & Reiss, E. L., A numerical method for ocean-acoustic normal
modes. J. Acoust. Soc. Am., 76 (1984) 244-52.
48. Schmidt, H., S A F A R I User's Guide. SACLANT ASW Research Centre, La
Spezia, Italy, 1987.
49. Schmidt, H. & Jensen, F. B., A full wave solution for propagation in
multilayered viscoelastic media with application to Gaussian beam reflection
and fluid-solid interfaces. J. Acoust. Soc. Am., 77 (1985) 813-25.
50. Schmidt, H. & Glattetre, J., A fast field model for three-dimensional wave
propagation in stratified environments based on the global matrix method. J.
Acoust. Soc. Am., 78 (1985) 2105-14.
51. Kutschale, H. W., Rapid computation by wave theory of propagation loss in
the Arctic Ocean. Report No. CU-8-73. Columbia University, Palisades, New
York, 1973.
52. Di Napoli, F. R., Fast field program for multilayered media. NUSC Tech.
Report 4103, Naval Underwater Systems Center, 1971.
53. Tappert, F. D., The parabolic approximation method. In Wave Propagation
and Underwater Acoustics, ed. J. B. Keller & J. S. Papadakis, Springer-Verlag,
Berlin, 1977, pp. 224-87.
54. Jensen, F. B. & Krol, H. R., The use of the parabolic equation method in sound
propagation modelling. Report No. SM-72, SACLANT ASW Centre, La
Spezia, Italy, 1975.
55. Lee, D. & Botseas, G., IFD: an implicit finite-difference computer model for
solving the parabolic equation. NUSC Tech. Report No. 6659, New London,
CT, 1982.

Ocean propagation models

199

56. Botseas, G. & Lee, D., IFD: wide angle capability. NUSC Tech. Report No.
6905, New London, CT, 1983.
57. Spofford, C. W., The FACT model, Vol. 1. AESD, Office of Naval Research,
Report No. 109, 1974.
58. Yarger, D. F., The user's guide for the RAYMODE propagation loss program.
NUSC Tech. Memo. No. 222-10-76, 1976.
59. Keller, J. B. & Papadakis, J. S., Wave Propagation and Underwater Acoustics.
Springer-Verlag, Berlin, 1977.
60. Brekhovskikh, L. M. & Lysanov, Yu., Fundamentals of Ocean Acoustics.
Springer-Verlag, Berlin, 1982.
61. Brekhovskikh, L. M., Waves in Layered Media. Academic Press, New York,
1980.
62. Weston, D. E., Acoustic flux methods for oceanic guided waves. J. Acoust. Soc.
Am., 68 (1980) 287-96.
63. Weston, D. E., Acoustic flux formulas for range-dependent ocean ducts. J.
Acoust. Soc. Am., 68 (1980) 269-81.
64. Harrison, C. H. & Ainslie, M. A., The suitability of the WKB method to fast
numerical acoustic pressure calculations. In Acoustics and Ocean Bottom, ed.
A. Lara, C. Ranz & R. Carbo, Federation of the Acoustic Societies of Europe,
Imprenta Moderna, 1987.
65. Biot, M. A. & Tolstoy, I., Formulation of wave propagation in infinite media by
normal coordinates with an application to diffraction. J. Acoust. Soc. Am., 29
(1957) 381-91.
66. Porter, M. B. & Bucker, H. P., Gaussian beam tracing for computing ocean
acoustic fields. J. Acoust. Soc. Am., 82 (1987) 1349-59.
67. Miller, J. F. & Ingenito, F., Normal mode FORTRAN programs for
calculating sound propagation in the ocean. NRL Memo. Report No. 3071,
1975.
68. Porter, M. B. & Reiss, E. L., A numerical method for bottom interacting ocean
acoustic normal modes. J. Acoust. Soc. Am., 77 (1985) 1760--7.
69. Ferla, M. C., Jensen, F. B. & Kuperman, W. A., High frequency normal-mode
calculations in deep water. J. Acoust. Soc. Am., 72 (1982) 505-9.
70. lngenito, F. & Wolf, S., Acoustic propagation in shallow water overlying a
consolidated bottom. J. Acoust. Soc. Am., 60 (1976) 611-7.
71. Kuperman, W. A. & lngenito, F., Attenuation of the coherent component of
sound propagating in shallow water with rough boundaries. J. Acoust. Soc.
Am., 61 (1977) 1178-87.
72. Labianca, F. M., Normal modes, virtual modes, and alternative representations in the theory of surface-duct sound propagation. J. Acoust. Soc. Am., 53
(1973) 1137-47.
73. Stickler, D. C., Normal-mode program with both discrete and branch line
contributions. J. Acoust. Soc. Am., 57 (1974) 856-61.
74. Bartberger, C., Comparison of two normal-mode solutions based on different
branch-cuts. J. Acoust. Soc. Am., 61 (1977) 1643.
75. Williams, A. O. Pseudoresonances and virtual modes in underwater sound
propagation. J. Acoust. Soc. Am., 64 (1978) 1487-91.
76. Tindle, C. T., Virtual modes and mode amplitudes near cut off. J. Acoust. Soc.
Am., 65 (1979) 1423-8.
77. Tindle, C. T., Stamp, A. & Guthrie, K. M., Virtual modes and the surface
boundary condition in underwater acoustics. J. Sound Vibr., 49 (1976) 231--40.

200

C. H. Harrison

78. Jensen, F. B. & Kuperman, W. A., Sound propagation in a wedge-shaped

ocean with a penetrable bottom. J. Acoust. Soc. Am., 67 (1980) 1564-6,
79. St Mary, D. F. & Lee, D. Accurate computation of the wide-angle wave
equation. In Computational Acoustics: Wave Propagation, ed. D. Lee, R. L.
Sternberg & M. H. Schultz. International Association for Mathematics and
Computers in Simulation. North-Holland, New York, 1988.
80. Knightly, G. H. & St Mary, D. F., Derivation of higher-order parabolic
equations which include density variations. J. Acoust. Soc. Am., 82 Suppl. 1
(1987) UU2.
81. Lee, D., Sternberg, R. L. & Schultz, M. H. (eds), Computational Acoustics:
Wave Propagation. International Association for Mathematics and Computers in Simulation. North-Holland, New York, 1988.
82. Brent, R. I., Siegmann, W. L. & Jacobson, M. J., Effects of uniform horizontal
currents in the parabolic approximation method. J. A coust. Soc. Am., 82 (1987)
545-58.
83. Robertson, J. S., Siegmann, W. L. & Jacobson, M. J., Acoustical effects of
ocean current shear structures in the parabolic approximation. J. Acoust. Soc.
Am., 82 (1987) 559-73.
84. Baer, R. N., Perkins, J. S., Wright, E. B. & Adams, B. B., An acoustic model for
bathymetric scattering with low-frequency applications. In Ocean SeismoAcoustics, ed. T. Akal & J. M. Berkson. Plenum Press, New York, 1986, pp.
345-53.
85. Baer, R. N., Propagation through a three-dimensional eddy including effects
on an array. J. Acoust. Soc. Am., 69 (1981) 70-5.
86. Perkins, J. S. & Baer, R. N., An approximation to the three-dimensional
parabolic-equation method for acoustic propagation. J. Acoust. Soc. Am., 72
(1982) 515-22.
87. Jensen, F. B. & Kuperman, W. A., Environmental acoustic modelling at
SACLANTCEN. Report SR-34, SACLANT ASW Research Centre, La
Spezia, Italy, 1979.
88. Jensen, F. B. & Kuperman, W. A., Consistency tests of acoustic propagation
models. Report SM-157, SACLANT ASW Research Centre, La Spezia, Italy,
1982.
89. McDaniel, S. T. & Lee, D., A finite-difference treatment of interface conditions
for the parabolic wave equation: the horizontal interface. J. Acoust. Soc. Am.,
71 (1982) 855-8.
90. Thompson, D. J. & Chapman, N. R., A wide-angle split-step algorithm for the
parabolic equation. J. Acoust. Soc. Am., 74 (1983) 1848-54.
91. Pierce, A. D., Augmented adiabatic mode theory for upslope propagation
from a point source in variable-depth shallow water overlying a fluid bottom.
J. Acoust. Soc. Am., 74 (1983) 1837-47.
92. Tolstoy, A., Berman, D. H. & Franchi, E. R., Ray theory versus the parabolic
equation in a long-range ducted environment. J. Acoust. Soc. Am., 78 (1985)
176-89.
93. Brock, H. K., Buchal, R. N. & Spofford, C. W., Modifying the sound-speed
profile to improve the accuracy of the parabolic-equation technique. J. A coust.
Soc. Am., 62 (1977) 543-52.
94. Spofford, C. W., A synopsis of the AESD Workshop on Acoustic Propagation
Modeling by Non-Ray-Tracing Techniques. Report TN-73-05, Office of
Naval Research, Washington DC, 1973.

Ocean propagation models

201

95. Davis, J. A., White, D. & Cavanagh, R. C., NORDA Parabolic Equation
Workshop. Report TN-143, Naval Ocean Research and Development
Activity, NSTL Station, MS, 1982.
96. Quality assessment of numerical codes. Special session in: J. Acoust. Soc. Am.,
80, Suppl. 1 (1986) I, R.
97. Numerical solutions of two benchmark problems. Special session in: J. Acoust.
Soc. Am., 81, Suppl 1 (1987) Q.
98. Tindle, C. T. & Weston, D. E., Connection of acoustic beam displacement,
cycle distances and attenuations for rays and normal modes. J. Acoust. Soc.
Am., 67 (1980) 1614-22.
99. Fitzgerald, R. M., Helmholtz equation as an initial value problem with
application to acoustic propagation. J. Acoust. Soc. Am., 57 (1975) 839-42.
100. McDaniel, S. T., Propagation of normal modes in the parabolic approximation. J. Acoust. Soc. Am., 57 (1975) 307-11.
101. McDaniel, S. T., Parabolic approximations for underwater sound propagation. J. Acoust. Soc. Am., 58 (1975) 1178-85.
102. Palmer, D. R., Eikonal approximation and the parabolic equation. J. Acoust.
Soc. Am., 60 (1976) 343-54.
103. Claerbout, J. F., Fundamentals of Geophysical Data Processing. McGraw-Hill,
New York (1976), pp. 206-7.
104. Thompson, D. J. & Wood, D. H., A post processing method for removing
phase errors in the parabolic equation. J. Acoust. Soc. Am., 82 (1987) 224-32.
105. Pierce, A. D., The natural reference wavenumber for parabolic approximations in ocean acoustics. Comput. Maths. Applics, 11 (1985) 831-41.
106. Harrison, C. H., A formula for reference sound speed (unpublished).
107. Weston, D. E., Intensity-range relations in oceanographic acoustics. J. Sound
Vibr., 18 (1971) 271-87.