JOURNAL OF GEOPHYSICAL RESEARCH VOLUME 67. No. 10 SEPTEMBER 1962
Propagation of Acoustic-Gravity Waves
in the Atmosphere
FRANK PRESS AND DAVID HARKRIDER
Seismological Laboratory
California Institute of Technology, Pasadena
Abstract. ]Komogeneous wave guide theory is used to derive dispersion curves, vertical
pressure distributions, and synthetic baregrams for atmospheric waves. A complicated mode
structure is found involving both gravity and acoustic waves. Various models of the atmos-
phere are studied to explore seasonal and geographic effects on pulse propagation. The influ-
ence of different zones in the atmosphere on the character of the baregrams is studied. It is
found that the first arriving waves are controlled by the properties of the lower atmospheric
channel. Comparison of theoretical results and experimental data from large thermonuclear
explosions is made in the time an,d frequency domains, and the following conclusions are
reached: (1) The major features on baregrams are due to dispersion; (2) superposition of
severM modes is needed to explain observed features; (3) scatter of d•t• outside the range
permitted by extreme atmospheric models shows the influence of winds for Ax; wind effects
and higher modes are less important for A• waves. A discussion is included on atmospheric
terminations and how these affect dispersion curves.
INTRODUCTION
Interest in the problem of the propagation of
a pulse in the atmosphere began when world-
wide pressure disturbances were observed in
connection with the explo,sion of the volcano
Krakatoa in 1883. The great Siberian meteorite
of 1908 provided additional data which were used
in attempts to correlate observations with the-
ories of pulse propagation. These studies were
prompted by a desire to account for the velocity
of the pulse, to explain its peculiar signature,
and to see if evidence could be found from pro-
gressive waves which would shed some light on
the existence of free oscillations of the atmos-
phere. A mode of free oscillations having a pe-
riod close to 12 solar hours is required by the
resonance theory of the solar atmospheric fide.
The atmospheric pulse was of further interest in
that it could provide information concerning
the structure of the atmosphere. (For a sum-
mary of early papers see Wilkes [1949].)
Interest in this problem was renewed with
the detonation of thermonuclear bombs in the
atmosphere. These 'megaton'-class explosions ex-
cited long atmospheric waves and provided data
from a world-wide net of sensitive baregraphs
[Yamomoto, 1956, 1957; Hunt, Palmer, and
x Contribution 1086, Division of Geological Sci-
ences, California Institute of Technology.
Penney, 1960; Oksman and Katajo, 1961; Car-
penter, Harwood, and Whiteside, 1961; Donn
and Ewing, 1961; Wexler and Haas, 1962]. Al-
though the earlier theoretical studies [Scorer,
1950; Pekeris, 1948; Yamomoto, 1957] provided
much insight into the nature of wave propaga-
tion, they were of limited use in analyzing the
observations because these investigators were
forced to assume oversimplified atmospheres in
order to obtain solutions.
With the advent of the high-speed digital com-
puter it became possible to obtain numerical
solutions for a more realistic atmospheric model.
In addition, the structure of the atmosphere is
sufficiently well known from rocket soundings
and satellite observations so that it is now not
a significant variable of the problem. Major
emphasis is therefore no longer placed on de-
ducing the structure of the atmosphere, but on
using a reasonably well known structure to ex-
plain the significant features observed on the
baregrams.
In this paper we present numerical solutions
for the homogeneous problem of wave propa-
gation in which the atmosphere is considered as
a two-dimensional wave guide. Phase and group
velocity dispersion curves a. nd vertical pressure
distributions are numerically evaluated for a
number of modes and are discussed in terms of
atmospheric structure. A comparison is made
3889
3890
between observations a•d theory in both the
frequency and the time domain. The inhomo-
geneous problem, in which the source and at-
mospheric excitation functions are also included,
will be treated in a following paper by the sec-
ond author.
Our procedure will be to represent the com-
plex vertical temperature structure of the at-
mosphere by a large number of isothermal lay-
ers. The solution to the equation of motion for
each layer takes a particularly simple form.
Boundary conditions at each interface and the
characteristic equation for the multilayered
wave guide are cast in a matrix formulation
suggested by Haskell [1953] which is particu-
larly suited for programming on a digital com-
puter. In practice, 20 to 40 layers are sufficient
to obtain an adequate approximation of the
real atmosphere. This approach is similar to that
of Pfeifer and Zarichny [1962], but our con-
clusions differ somewhat from theirs. Yamomoto
[1957] and Hunt, Palmer, and Penney [1960]
also used an isothermal layer representation,
but they limited themselves to only a few layers.
Although the latter authors were in error in
their formulation of the interface boundary
conditions, much of their discussion is still use-
ful in elucidating the nature of the atmospheric
pulse.
TI-IEOR¾ AND •U1VIERICAL METI-IODS
The linearized equations of motion for a con-
stant velocity layer in a horizontally stratified
atmosphere are given by [Pekeris, 1948] 2
_•_ •2 2 2
- • • + •(•- •) x• = o (•)
gmk Otto -- gm"Y• )Xm (2)
•,•p,•(z) -- i, op2(z)
2. 2 2
g,•C•m X,• + ('ygm •' (3)
ß - - ,,, "m )Xm]
•,• -- g,•'k 2 -- w 4 (4)
where the dots denote differentiation with re-
spect to the vertical coordinate z. z is taken
•. For definitions of symbols not defined in the
text, see appendix.
PRESS AND HARKRIDER
positive in the upward direction and a space and
time dependence Jo(kr) exp (io•t) is included.
Assuming azimuthal symmetry, X• is the first
time derivative of the dilatation and is given in
cylindrical coordinates by
Ow 10(ru)
x = •; + (•) r Or
where u, w, and p denote perturbations from
equilibrium of horizontal velocity, vertical ve-
locity, and pressure, respectively.
In addition, we have for the equilibrium state
in each layer
dP ø o o
dz -- -- gP P -- RKøpø
o
a -- 'y- 'yRK ø (6)
p
From (6) we obtain for layer m
o
Pm (Z) --' pmO(Zm_l)e -2xm( ...... )
0 --2Xra(Z--Zra--x/2)
'-- pine
where Xm = 'Ygm/2 a• •, p? = pmO(Zm--•/2), and
z•_•/2 is the altitude of the midpoint of the m
layer. Since p?(z•_•) = p•_ P(z•_•), we have
pmO(Zm--1) Ogm_12pmO(Zm_l) Ogm_l 2
= • o - • (7)
Om-•ø(Zm-•) "• Pm-• (Zm-•)
We assume •all motio• and impose the
bound•w conditions of continuity of vertic•
pa•icle velocity •d total pressure across the
disturbed interfaces. Retaining o•y first-order
terms, we find that the change in total pressure
of a small p•rcel which is displaced • vertical
distance V from its static equilibrium position
z is the pressure perturbation, p(z), at the zero
displacement position plus 8p -- -- gpø(z)v.
Now defining p• (z) -- p (z) + 8p, using (3)
and w -- i•v, we obtain
p•(z) i o2(z) • (s)
• • •m Xm
At the l•yer interfaces we will now require
p•_•(Zm-•) =
in order to guarantee continuity of pressure.
It is interesting to note that when there is no
temperature or gritty discontinuity across an
interface, one can use as a boundary condition
ACOUSTIC-GRAVITY WAVES IN THE ATMOSPHERE 3891
continuity of p(z) since for that particular case
p•..(zm-•) = p•-._•(zm_•) is equivalent to
Pm(Zm-O = Pm-•(Zm-O. This relation was used
by Pekeris [1948] in his model of the atmosphere.
The general solution of (1) is given by
Xm = eX•Z[Am'e-i•r•' • - Am"ei•r•Z]Jo(kr)eiWt
(•)
where
kr • •n
(lO)
for (c, k) such that (kr•.) •' > O, and
kr•m - --i k 2 1--•.,• -- • 1
for (c, k) such that (kr•.) • ( O. Here •m
(2 • _ •/•)•m ( •m for all • • 1, and
is the Brunt resonant angular frequency for the
constant velocity layer m and is given by
•.m • am • -- 1/•m.
Substituting (9) into (2) and (S), evaluating
at Zm and z•_•, and eliminating the constants
•m' and •", we obtain the following matrix
relation
[w•(z•)•= [-(am)• (am)•••Wm(Zm--•)• (11)
p•(Zm): _(a2•1 (am)•:•p•,,,(Z•-l):
where
(a2• = e
ß c os rm + • k• -- (•r • •) J
(a21• = i(•c) •
•m k • + •m•(•r•.)• •in •m
' •mø•m • •m (•r•2 (1
P m ø •m sin P m
(am)• = i
( •c) • ( •r•.)
-kmdm
(am)• = •
cos '•
•m (•r•.)J
and Pm= (kr•) din. In equations 12 we see
that mat•x elements (am)i• are real or imaginary
for (c, k) real and for j • k equal to even or
odd integers, respectively. Therefore, the ele-
ments of a matrix resulting from the matrix
multiplication of any number of layer matrices
will be real or imaginary in the same sense as
the individual matrices.
For gm : 0 the am matrix reduces to a form
equivalent to the nongravitating liquid-layer
matrix given by Dorman [1962] in his discussion
of elastic wave propagation in layered wave
guides.
The conditions for continuity of w and p•. at
interfaces and the connection between layers
described by (11) together enable us to write
the following matrix relation
Wn-t(Z•-D 1 -- A I wø(0) 1 (13)
•_•(z•_•): •p•o(O)
where A = a•_•...a•.
At z = 0, layer I is in contact with a fiat
rigid boundary where we require w0(0) - 0 and
thus peo(0) = p0(0) ---- p0. Equation 13 becomes
Ip•n-i(Zn-i)l '-- AIpl (14)
p•--i(Zn--,)-J 0
Before deriving the period equation for an
atmosphere terminated by an isothermal half-
space extending from z•_, to infinity, we consider
two special cases.
The first is an atmosphere bounded at z•_, by a
free surface. For this case we have p,•_, (z•_,) - 0,
and (14) reduces to
0 po
which in turn yields A•.•. po = 0. Therefore, the
period equation for an atmosphere with a free
surface at z._, is
A2• = 0 (16)
The second case is an atmosphere bounded by
a rigid surface at z•_•. Here we have w,•_•(z•,_•) =
0, and (14) reduces to
0
which yields the following period equation:
Axe* ---- 0 (17)
where iAi•* = Ai•, and Aiz:* is real for j • k
equal to an odd integer.
3892
For the case of an atmosphere bounded by an
isothermal half-space, we require that the nth
layer coefficient Ad' - 0. For (kr,•) • > 0 this
is equivalent to requiring that there be no
radiation from infinity into the wave guideß
For (kr,n) 2 < 0 this requirement guarantees that
the kinetic energy integrated over a column of
atmosphere will be finite.
Setting Af - 0 in the solution for w•(z) and
pv.(z) evaluated at z•_• we find that
,(zn_l)J Lp,n(Z•_,)
where
•nZn--!
b• e
-- e
PRESS AND HARKRIDER
A•:* q- (kc) •
Lb•,•J
(19)
ß g•5- -- -- it•. •(kr.,,) Jo(kr) e i•'
•knZn-- x
ß e -ikranzn_•(kc)2
pnO(Zn-1)Oln 2 i• t
ß
Substituting (18) into (14) and eliminating
from the two resulting linear equations, we
obtain as the period equation for an n layered
hMf-space
bin
A,•A• = 0
b2•
or
A•* + (kc) a g'•k? -- zo•. (kr,•.)
p.ø(z•_,)a.•' •
ß A• = 0 (20)
All quantities in (20) are Mways real for all
(c, k) real except for i(kr•), wMch is real or
imagMary dependMg on the wlues of the real
(c, k). Since in tMs paper we are interested in
undamped propagation, we now make the
requirement that (kr•): be negative. This ex-
cludes leaMng or complex modes of propagation.
Under tMs condition, (20) is real and takes the
form used in numerical calculation of dispersion
curves:
O/n
0, 2
ß A• = 0 (21)
The condition that (kr,•) • be negative now
prescribes a cutoff region in the (c, k) plane
defined by (kr,•) = 0. The boundary of this
region is obtained in terms of c and period T by
setting (10) equal to zero. This yields
2
T TBa(•nn • 1)1/2/( 3 -- 1) •/• = c _
where Tr• is the Brunt resonant period of the
half-space and is given by Tr• = 2•r/ar•. From
(22) we have the following asymptotic values of
the boundaries of the cutoff region:
T--•(lS•iTra for
\O•n/
T = 0 for
No cutoff region for
T --• • for
T = T.• for
No cutoff region for
c• co
C--Ol n
•s• _< c _• • (23)
C'--'•n
c=0
In deriving (13) we have as a by-product the
following matrix relation
Wm(Zm)l •-- Amlp ]
Pm(Zra) J 0
where A• = a,,•...ax
Furthermore, at the layer midpoints, z.•_•/a, it
can easily be shown that
Pra(Zra--!/2) LP0J
(25)
where A•_•/: = am_l/aa,•_i ''' al and a,•_l/:
is of the same form as a• with d• replaced by
d,•_•/: = d.•/2 and all other quantities remain
unchanged.
Rewriting (25) we obtain
Po
= ( A•-•/2)1•* (26)
where w = iw* and p•,.(Zr•-l/a)/po =
ACOUSTIC-GRAVITY WAVES IN THE ATMOSPHERE 3893
From the definition of pp• we have
Pm(Zm-1/2) __ (Am_1/2)22 'Jr- DinOgre *
po - (kc) ( Am-1/2)12
Now defining the r-independent part of p,,(z, r)
as Pro(z) or p,,(z, r) - •,,(z) Jo(kr), and since
u(z, =
I Op,•(z, r)
we obtain
am$(Zm--l/2) 1
--
•o cp,•
I o
P•
ß (Am-1/2)•2 .qt_ (kc) (Am-1/2)12* (28)
where u,,,(z, r) = i a,,*(z) J•(kr).
The normalized particle velocities and pres-
sures are given in (26) to (28) at the layer mid-
points rather than at the layer boundaries in or-
der to make them correspond more closely to
the smooth distribution in nature.
The matrix formulation described above is
very convenient for numerical calculations; it
was first used for earthquake surface waves by
Haskell [1953]. For our computations we pro-
grammed the dispersion calculation for an IBM
7090 computer. This program was written in
Fortran and used for production runs. As an
independent check on the program, we also
wrote a program for our small computer, the
Bendix G-15D. The Fortran program has the
option of calculating the three models for ter-
minating the atmosphere discussed above for
any given layering of the atmosphere.
The general computational procedure is to
find the zeros of a function F of phase velocity,
wave number, and the physical constants of the
layers. For an atmosphere bounded by an iso-
thermal half-space, the F function is defined as
the left-hand side of (21); for an atmosphere
bounded by a free surface, the left-hand side
of (16); and for an atmosphere bounded by a
rigid surface, the left-hand side of (17).
The flow of the program is similar to that in
the programs described by Press, Harkrider, and
Sea/eldt [1961] and Harkrider and Anderson
[1962]. The zeros of F are determined by ini-
tially specifying the phase velocity c and a trial
value of the wave number k. The elements of the
an matrix are formed at each layer and then
multiplied by the matrix of the layer above it,
starting with the layer at the surface. After the
matrix product for all layers has been calculated,
the program then combines these numerical
quantities to obtain a value for F. New trial
values of k (of increasing or decreasing size
depending on the sign of the initial F value)
specified by an input Ak are used to calculate
new F values until the root is bracketed by a
change of sign in F. Linear interpolation and
extrapolation are then repeatedly used to find
small F values until k's of different F sign are
within the precision interval desired. The re-
sulting interpolated value of k is the output
value given as the root for the input c.
The program has an additional feature in
that, as an input option, the first or second
roots (two smallest k roots) will be found for
a given c. This is accomplished by starting at
the smallest k outside the cutoff region (F com-
plex) and finding either the first or second sign
change of F. The roots associated with either
mode are then computed. For all values of c
in the free and rigid surface models and for
c _• fin in the isothermal half-space model, this
initial k is zero. For c • fin in the isothermal
half-space model, the initial k is determined by
(22).
In all models calculated we found that each
continuous dispersion curve or mode was a
monotonic decreasing function of c versus k or
T and always had the same sign change in F
through the root region. This made it easy to
track all the roots of a preselected mode. To
save computer time and to keep from jumping
modes, the k root for the previous c is used as
the starting point for the new c. Further de-
tails about root hunting procedures can be
found in papers by Press et al. [1961] and Hark-
rider and Anderson [1962].
Once a root is found by the computer, the
velocity and pressure ratios given by (26), (27),
and (28) are calculated at the midpoint in
each layer. The vertical distribution of these ra-
tios is generally diagnostic of the particular
mode and provides a check against mode jump-
ing. The group velocity is computed by nu-
merical differentiation of the phase velocity
values. Ac/Ak is obtained by perturbing c
slightly and then finding a new k root.
The program has two options for input o.f
layer constants. The first reads d• and K• ø from