AFCRL-TN-60-649 THE UNIVER SITY OF MICHIGAN 2886-1-T PRESSURE PULSE RECEIVED DUE TO AN EXPLOSION IN THE ATMOSPHERE AT AN ARBITRARY ALTITUDE PART I STUDIES IN RADAR CROSS SECTIONS XLI by V. H. Weston August 1960 Report No. 2886-1-T on Contract AF 19(604)-5470 Prepared for GEOPHYSICS RESEARCH DIRECTORATE AIR FORCE CAMBRIDGE RESEARCH LABORATORIES AIR RESEARCH AND DEVELOPMENT COMMAND UNITED STATES AIR FORCE BEDFORD, MASSACHUSETTS

THE UNIVERSITY OF MICHIGAN Qualified requestors may obtain copies of this report from the Armed Services Technical Information Agency, Arlington Hall Station, Arlington 12, Virginia. ASTIA Services for Department of Defense Contractors are available through the Field of Interest Register or a "need-to-know" certified by the cognizant military agency of their project or contract. ii

THE UNIVERSITY OF MICHIGAN 2886 -- 1 - T TABLE OF CONTENTS Page Introduction 1 Initial Nomenclature 4 1. Characterization of the Source of the Pressure Pulse 4 2. The Equations of Motion 6 3. Boundary Conditions 9 4. The Unperturbed Variables 10 5. The Green's Function 11 6. Eigenfunction Expansion for the Green's Function 13 7. An Approximate Technique to the Determination of the Low Frequency Free Wave Modes 19 8. Comments on the Accuracy of the Approximations 29 9. Simple Source Model 33 10. Method of Evaluation of p (a, 0, t) 34 11. Comments 46 APPENDICES A: Simple Source Representation B: Radial Equation for the Vertical Velocity 00 C: Evaluation of the Integrals F (q) = -1/2+nexp i wq - d w References iii

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T INTRODUCTION There have been many recent efforts devoted to the study of "free" waves in the atmosphere and the pressure pulse produced by large explosions such as the Krakatoa volcanic eruption of 1883, the great Siberian meteorite of 20 June 1908, and, recently, nuclear explosions. As an example of the latter, the Japanese have recorded such atmospheric waves from nuclear explosions (R. Yamamoto, Ref. 1). The pressure pulse produced by such explosions can be decomposed into two portions, a low frequency "gravity" wave train, plus high frequency "acoustic" wave trains. Theoretical analysis has been mainly devoted to the "gravity" wave portion of the pressure pulse produced by large explosions on the ground for simple models of the atmosphere. (Pekeris, Refs. 2 and 3, Scorer, Ref. 4.) Recently Dikii (Ref. 5) has discussed the gravity and acoustical type waves. Penney et al, (Ref. 6) have calculated both the gravity wave and acoustic wave portions of the pressure pulse produced by explosions on the ground. Here the interest will be mainly concerned with the gravity wave, produced by explosions not only on the ground, but at various heights in the atmosphere. At present the effect of winds is ignored, but will be considered later. They do have an influence upon the pressure pulse, although the high-frequency acoustical portion is affected the most.

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T It is assumed that the explosion is symmetric about an axis normal to the earth's surface, and that beyond a specific distance away from the center of the explosion the perturbed values of pressure etc, will be small compared to the unperturbed value. Hence beyond this distance (characterized by a surface ) the hydrodynamical equations may be linearized. The explosion can then be represented in terms of the excess pressure and normal velocity on this surface,. Unfortunately analytical solutions are not available for large explosions. An analytic solution has been found for an intense spherically symmetric explosion (J. L. Taylor, Ref. 7) but it holds down to only over-pressure of about 20 atmospheres. Hence for good source models, values of the excess pressure and normal velocity must be obtained from observational data. Hence the hydrodynamical equations are linearized with the viscosity terms neglected (actually these become important at very high altitudes, roughly around 200 km). The time dependence of the equations is removed by taking a Laplace transform of them. A second order differential equation involving the excess pressure transform is obtained. The boundary conditions that are imposed, are, that the vertical velocity must vanish at the earth's surface, and the total kinetic energy in a solid angle subtended at the earth's surface be finite. This latter condition (which is a restriction upon the behavior of the excess pressure at high altitudes) is the same condition Scorer and Penney used. A ring source Green's function is then defined. It is shown that the excess pressure may be represented in terms of an integral over the surface surrounding the source, the integral containing the Green's function together with the values of the excess pressure and normal velocity on>.

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T An eigenfunction expansion is then obtained for the Green's function involving the various modes. Only the contribution die to the modes corresponding to the "free' waves (those waves which propagate around the earth without exponential attenuation) are then retained. It is shown that for certain temperature models of the atmosphere a good approximation can be made to the "gravity" wave mode. This approximation gives good results for Scorer's model of the atmosphere. The pressure pulse for the directly received wave (as contrasted to the antipodal wave) is computed for various ranges, for two models of the atmosphere. In doing so, a simple source model is taken, namely a point source in space with a delta function dependence in time. The intensity of the explosion is given in terms of volume of gas introduced. In the calculation of the pressure pulse at various ranges from the source, emphasis is placed upon the main body of the gravity wave portion. The tail of the pulse up to the so-called cut-off point is not considered. In calculating the head of the pulse a new asymptotic technique is introduced which gives very good results, for intermediate and long ranges. -3 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T PRESSURE PULSE RECEIVED DUE TO AN EXPLOSION IN THE ATMOSPHERE AT AN ARBITRARY ALTITUDE PART I Initial Nomenclature: (1) Spherical polar coordinates (r, 0, ) (2) Radius of earth r = a (3) T = absolute temperature R = gas constant /- = ratio of specific heats p = gas pressure p = gas density c2 IRT =RT - po/po u = velocity = (ur, u0, u) (4) Unperturbed state denoted by subscripts zero, i.e. T, p, etc. (5) Perturbed state is represented in terms of the unperturbed state plus a perturbation term (which has no subscript) i. e. p + p, To + T... 1. Characterization of the Source of the Pressure Pulse A large explosion in the atmosphere is usually detected at large distances by changes in pressure. Hence the excess pressure p(r, t) will be solved for. At a sufficient distance from the source of the explosion it is known that the excess variables, pressure, density and temperature are small in comparison with the corresponding unperturbed variables for the atmos - phere. -4-.

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T Let the source of the explosion be at the point (R 0, 0) where Ro > a (the radius of the earth). It will be assumed that the explosive effects are independent of the azimuth angle 0. The source can be enclosed by a surface of revolution > with axis of revolution 0 = o (i.e. the z axis), such that outside this surface, the excess pressure, density, etc. are small in comparison with their respective unperturbed values. The surface E may be a small sphere with centre (Ro, 0, 0). However, its actual form will not be specified. Let c represent the region outside 1 and outside the earth (r = a). Hence for the region 7, the assumption that the excess variables are small in comparison with their respective unperturbed values, is valid. + The source of the pressure pulse will be characterized by boundary conditions on the surface. To specify the problem uniquely p (r, t) and n. u for r on > must be specified. n represents the unit outward normal to A. Hence set the initial conditions. u(r, t) E 0 p(r, t) E 0 t < 0 andreco p(r, t) - 0 In addition viscosity and the effects of winds will be neglected for the present. -5 -

THE UNIVERSITY OF 2886 - 1 - T p(r, t) = 0 t < f (r, t) t > n. u s 0 t < g(r, t) t > MICHIGAN 0 0 0 0 r E where f (r, t) and g (r, t) are given functions such that f (r, t) and g (r, t) both vanish when t-ooo. 2. The Equations of Motion. To begin with note that the explosive effects are independent of 0, hence all the variables will be independent of 0. From Ramsey (Ref. b ), we have the equations of motion au +(u.)u + K - at 1 _ (p+po) + (-g, 0) (P+P, ) (1) where 1: - + K = - [u+ Iu2 cot0 +, u + u -u cot cot the equation of continuity a (p + p) + (u. ) (p + o) (+ p )P 7. = 0, 8t ~ 0 0 - (2) the adiabatic energy equation a ( +a ----. )(p+pp) = c (P+Po)+(uV) (P+Po), (3) at (P+op0 ) aP -6 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T and the equation of state (P+P) = (p+pO) R(T+To) If use is made of the hydrostatic equations for the unperturbed atmosphere and second order pressure and density terms neglected, equation (1) becomes au 1 gp - + (u. 7) u + K — (,, 0). (4) at - 0- pAlso equations (2) and (3) become ap — o+ Ur P + p7 u = 0 (5) ap, P+ ur c + ur P (6) at r o za where the prime indicates differentiation with respect to r. Now, since all the variables u, p, p vanish for t<0 and r E, take the Laplace transform of the above equations, and set 00 U(r, s) = $ et u(r, t)dt O 00 P(r, s) = eSt p(r, t)dt (7) P(r, ), = e-t p(r, t)dt 0 -7 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T Hence 1 g O sU+- VP +(,, 0) = 0 (8) Po Po r 0 s~ + UrpO + pO U = O (9) sP+Ur p = c p + Ur (10) [Since u (r, t)-0 or t-,oo for all r- oCS, U (r, t) is bounded above. Hence it can be shown that for sufficiently large s, the transforms of the velocity squared terms are negligible in comparison with sU (provided U ~ 0). If the hydrostatic equation = -gP is used, equations (10) and (9) become s LP-co] = U [cP + gPo (11) s c2V. U+- P-gUr = 0 (12) Po Using equation (11), eliminate f from equation (8) to give sU p 1/2 = -L (P -1/2 p) - Ph- p-1/2 A (13) sUo r 0 where g h = _ -+ (14) s2 P 0 g (15) A - + - (15) Ic2 2 p o0 o0 -8 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T and the operator L has components i1 a L -( -- - h ar 1 a 1 r aO rsin0 a0 (16) The velocity vector U may now be eliminated from (12) and (13) giving -1/2 p)q(r) 'V. L(p -1/2 P) + -, r (p -1/2 p) 0 (17) with q(r) 2 r2 C2 L r2 + — A2 h h' 2 + A' - A-+-A h r (18) 3. Boundary Conditions: For r on the surface,, the following is obtained P(r, s) = o e -t f(r, t) dt e -t g(r, t) dt [: =E(r, s) = F (r, s).1. (19) n. U(r, s) 0 o It is also required that u. n = 0 on the earth's surface (i.e. at r = a). This condition becomes U r = 0 for r = a. Hence from equation (8) and (10) one obtains aP+ g p= 0 ar C2 r = a (20) -9 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T or ( -1/2 P) + A(p -1/2 p) 0 r = a (21) arO r = a (21) where A is given by (15). A boundary condition is required for r-oo. Here the condition that will be placed, is that the total kinetic energy in a solid angle whose vertex is the center of the earth, be finite. Hence 00oo 0O 3 r2 p (r) U 2 dr < oo and i p 1 PJ -dr oo "a ~ a 4. The Unperturbed Variables: Before proceeding further, several relations are needed involving the unperturbed variables po, po, and To. From the hydrostatic equation 0 -gP and equation of state PO = RpoTo we obtain on integration the following 1 r P (r) = Po(a)exp L — 5 (-) dr (22) R a To Po (a) g 1 PR (r) exp (23) R To (r)ex R 0 0 a 0~~ -10 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T The following relation can be obtained from (23) Po Ta_ P - T- 2 (24) p0 To c 5. The Green's Function: Define a ring source Green's Function G (r, r ) to be the solution of - -o q(r) ' (r -ro) ( 0 - 0 ). LG + G -- -r (25) r2 2 r sin 0 r2 which is integrable squared over (a, oo) and satisfies the boundary condition a G + AG 0 for r = a (26) aar Later on an eigenfunction expansion will be obtained for G (r, r ). However in the remainder of this section it will be shown how the Laplace transform of p (r, t), namely P (r, s) satisfying the boundary conditions in Section 3, may be derived from G (r, r ), Recall the equation for the excess pressure, namely. L (p -1/2 p) + (p -1/ P) q 0 (7) Multiply equations (27) and (25) by G and p / P respectively, then subtract the two. If in the resulting equation we interchange r and ro we have G (r, r) o Lo [p-l1/ (r) P (r, s - -1/2 (ro) P (r, s ~ Lo G (r, r) G~~~~~~~ (to r)0'L o- ooo- - o 2 -o-1/2 (, (r, s S(r- ) S(- - ) (28) rro o 0 LC-o>r sin 0 ro -11 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T This may be rewritten in the form 0. G L (p -/2 p) - (p-1/2p) L G = lr -/ P z f(r o) S(O ) (0 -. -77 (29) ro o 2 7r sin 0 Integrate over the domain P with variables of integration ro. We obtain using the divergence theorem, -1/Z -1/Z p) - n G L (p 1 P) - (p 1 P)L G dSo - /2(r) P(r, s) (30) S where n is the unit normal to S directed inwards in i. The surface S comprises of the surfaces, r = a, and r = oo. Since both G and P are squared integrable on (a, oo), the surface integral at infinity vanishes. The integrand for the surface r = a becomes 1 (PP) (p/Zp) G1 h G r o <r which vanishes since G and (p 1/2 P) satisfy the same boundary condition at r = a. Hence (30) becomes P(r, s) = p1/2 (r) fn.GLo(p 1/ p) - (pl/2 P) LO dSo To simplify the integral in (28), note from (13) that (32) L(po1/2 ) = -sUp 1/2 - pp 2Ah-1 (33) -12 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T thus (32) becomes P(r, s) = -sp1/ (r)J G(r, r )pl/2 (r)n. U(r, s)dS0 -Po/(rf /2 PO (r) P(ro s)() Gn. r+ n L G dS (34) Thus it is seen that if n. Ur and P are known on the surfacee, the transform of the excess pressure P (r, s ) can be found. 6. Eigenfunction Expansion for the Green's Function: Let Lr and L represent the differential operators a /r2 a V\ Lr = - + Yq(r) (a r oo ) (35) ar h a r 1 a a8 LV = _ —_-(sin -) (o007r) (36) sin 0 a0 a0 Ience the Green's function is given by (Lr+L0) G(r, ro) - (r ro) 2 (0i-no) r -0?-72 sin 0 First the operator Lr must be investigated. This can be written in the normal form for the Sturm-Liouville operators a + (r) with ) = 38) L r~l - ar[~r) ]+ 1 'q (r) with p(r) = h (38) Lr- = ar r h The parameter s may be chosen sufficiently large such that h(r) is positive over the range (a, oo), hence p(r)> 0 for (a, oo). — 13 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T We will consider the case where the parameter s is real and positive. Thus, p (r) and q (r) are real for (a, oo). The following restriction will also be placed: any model of the atmosphere that is used must have p (r) continuous. The asymptotic values of p( r) and q(r) for large r have to be investigated. For large r, c2 = KRRTo has the asymptotic value 1 c2r c2 (co) + 0(-) (39) p 1 where c2 (oo) = 'RTo (o). Equation (24) obtains ~._ 0 ( -) for large r. P0 r Thus it is seen that 1 A'JO (-) r2 and (40) 1 h- 1 + 0 ( r which together with equations (22) and (38) gives the asymptotic values for large r p(r) r2 (41) r2 S2 q (r)) - c2(o) c 2 ( coD) The equation LrY - X = 0 has two solutions w1 (r, X) and w2 (r, X) whose asymptotic values are 1 w (r, X) -erk(oo) r -14 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T 1 wZ (r, X)- e-rk(o3) r s2 where k (oo) is the positive root of k2 (oo) =.The second solution c2(oo) c 2 ( 00) w2 (r, X) is square integrable on (a, oo). Let ((r, Ai ) represent the set of solutions of Lri -XiK = 0 (42) with the property that 00 f ()) dr oo (43) a and d -- + A1 = 0 for r = a. (44) dr From Friedrichs (Refs. 9,10 ) it may be shown using (41) that the spectrum of the operator Lr is totally discrete, and if Z-Arepresents the manifold of all functions V/(x) such that 00 IJ y 2 dr< oo. a then for CYJ and the V/ such that the functions Lr r.c, the following expansion holds lru^ ai t(r i) (45) LrVY/ Xi ai/(r, Xi) (46) Therefore, the Green's functions G (r, ro) may be expanded in such a form. -15 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T 1 In order to simplify the analysis set X = - (, 2 + — ) andlet 0 (r, i) represent the solutions of 2 1 Lr0 + (i +-)0 = 0 1 4 (47) so that 00 f02) a dr oo and d r (r, ( i) do-(r, -i) +AO(r, u.)) 0 dr at r = a. (48) For the present let 0 (r,, ) be the solution of Lr0 + (2 +)0 = 0 L r + (Pt 2 +_) = 0 (49) which belongs toc. Multiply equations (47) and (49) by 0 (r, p ) and 0 (r, u ) respectively. Then subtract the resulting two equations and integrate from a to oo. Using (48) one obtains 00 + a2 0 (a, pi) d0(r,pA) J 0(r, /1) (r, i) dr = (2-) ------ + A a E2 -~ ) (h) dr a I ia (r, r=a r=a (50) Let /p = u j where i f j, then it can be seen immediately using (48) that 00 a 0(r, pi) 0(r, p.j) dr = 0 (51) Let pu approach,. i and take the limit, obtaining OD (r, ui 2 d r =ia at0 a 3 (a, I)) f 0(r,pJ)2 dr = - - A- A( )a, )) p)+ a 2t i (h) r a a Ip p=i (52) -16 -

THE UNIVERSITY OF M/lICHIIGAN 2886 - 1 - T Expand the Green's function G (r, r ) in terms of 0 (r, p i) by G(r, ro) = 2 ai 0(r, i) i (53) Substitute expansion (53) into equation (37). Using the property given by (46), one obtains 0 (r, pi) i 1+ 4 + L0 a 8 (r -ro) S( - 0o) (54) 2 r sin 0 Multiply both sides of (54) by 0 (r, p j ) and integrate from a to co, thus giving L 0 - +- a. -(uj 4 J (0-e00) 0 (ro, uj) (55) 0 0 2 7r sin 0 If (r, /j)2 dr a 0 (ro, pj) Set a. -oSJ 0(r, uj)2 dr a hence equation (55) becomes xj (0) (56) 1 a axj n — (sin 0 — ) - sin 0 a 0 a (2 1 (A i + 4 ) xi &(0 - Oo) 2 7r sin 0 (57) Solutions of the homogeneous equation corresponding to (57) are the Legendre functions P-1/2 + ij (cos 0) and P-1/2 + ip (-cos 0) (58) the first of which is finite at 0 = 0, and the second finite at 0 = 7r. -17 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T By the usual techniques it can be shown that the solution of (57) which is finite over the range (0 0! 7r) is given by P1/2 + ipj (cos 0) P-1/2 + il j (-cos 00) 0 0 0 1 X = - 1 (59) 4 cosh 7r i 1/2 + ij (cos 0) P-1/z + i/. ( -cos 0) 0 0 The Green's function expansion is now given by _ o j s ~(r,/~j)2 dr (60) G (r, ro) = I0 (r, j) (r j) Xj (0) / 0 (r, Uj)2 dr (60) a However the interest for the particular problem at hand is in evaluating the excess pressure p (r, t) on the earth's surface i.e., at r = a, at a long distance from the pressure pulse source. Hence the value of G (r, r ) for r = a, and 0 such that 0> 00 is required. In this case (60) simplifies giving G(a, ro) = fj(a,0; ro, ) (61) j 2 cosh 7u j with j 0(r0, Ij ) Ptj (cos 0) Ptj (-cos 0) _ 3~~r r a jA r /u +._A.(r, + ( )1} /r \ where for simplification tj has been set equal to i/u j - 1/2. However to evaluate the inverse transform, we require the Green's function for s complex. But for s complex, p(r) given by (38) is complex. For this case, one cannot deduce See Friedlander (Ref. 14) page 170. See Friedlander (Ref. 14) page 170. -18 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T that there exists a totally discrete set of eigenvalues Q j, hence the expansion given by (6) is questionable. For p (r) discontinuous we apply the boundary condition that the pressure and vertical component of the velocity must be continuous at any points of discontinuity of To, the orthogonality properties of the eigenfunctions are retained. Thus we may write n(s) t. G(a, 0; ro, 00) + R(p) (63) j - 2 cosh Jir.j where the summation is over the t'free" wave modes (those modes such that ij is real for s pure imaginary, or, setting s = iw, for real frequency w). The remainder R (uA ) represents the contribution due to the remaining modes. For present purposes it does not matter whether the remaining eigenvalues form a totally discrete set or not.' 7. An Approximate Technique to the Determination of the Low Frequency Free Wave Modes: We now require the values of the eigenvalues p j and the eigenfunctions corresponding to the free wave modes for real frequencies w. This requires a solution of the differential equation (47) with s replaced by io. However, this equation is difficult to solve for most temperature models of the atmosphere. Simple solutions can be found for an isothermal atmosphere (Pekeris, Refs. 2, 3) or an atmosphere composed of isothermal strata (Yamamoto, Ref. 1 Penney, Ref. 6 ). However, for the case of Scorer's model, (Ref. 4 ) (isothermal stratosphere and tropsphere with a constant temperature gradient), numerical methods have -to be used. IlThis is explained in more detail in Section 9. -19 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T For these various models, it is shown that there is a low frequency mode which we will call the gravity wave. This low frequency mode possesses a continuous spectrum in w from U = 0*to some cut-off frequency c c. The cut-off frequency depends upon the model of the atmosphere. The low frequency mode is essentially a surface wave, the energy being propagated around the earth. However above cut-off, the character of the wave changes into a wave which propagates energy upward. The temperature of the actual atmosphere increases quite rapidly above 100 km. For instance, at 140 km the temperature is given as 880~ K. For models of the atmosphere possessing this characteristic, namely that the upper atmosphere is the warmest, there exists a set of higher frequency modes with a continuous spectrum. Penney et al have shown these modes or "branches" for models comprised of isothermal layers, with the top layer the warmest. We will give here a new approach to the problem which will enable one to find good approximate analytical solutions to the gravity waves for some simple atmospheric models. The approach here is based upon an approximate technique for solving the radial equation for the radial component Ur of the velocity. Set /A2 +- _= -co2 a2 (64) 1 1/2 = iW? cZ a + - (65).... For some models of the atmosphere there is a lower frequency cut-off at a non-zero value of w. -20 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T where X -1 is the phase velocity of the wave*. From Appendix B it is shown that [r2 /2 - 2 a2 1/2 Ur (r, 0, Lw) = P-/(cos 0)/ - M (r) (66) r ill - l/2 r2 p1/M ) r2po0 where M is a solution of the equation 3 a 1, a ha M + M - (-) + — [Aa + - A' - A2 +- } = 0 (67) 4 ( a 2r r2 where -a(r) = w2 K -2a2] The best way to make use of this equation is to break the atmosphere into two+ regionsspecifiedby (i) a r b and (ii) b;r. Inthe lower region (a c r a b), the atmosphere will be composed of the troposphere and possibly the lower portion of the stratosphere. In the upper region (b A r), the temperature will be assumed to be a very slowly varying function of altitude. We This is seen using the asymptotic relations for large real wXa,- ip/cuXa, / z and P-wXa - 1/2 (cos ( ) / cos 2xa os- wXand noting / sin 0 LXa that a 0 is the arc length along the earth's surface. Hence for harmonic time dependence, X-1 is the phase velocity of the wave propagating around the earth. Note: At the interface of the two regions we will assume that the temperature is continuous (although the derivatives may or may not be continuous). -21 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T will use two different approximate techniques for finding M and hence Ur (r, w) in each of these regions. The solutions will be matched across the interface r = b, by requiring that the vertical velocity and pressure be continuous functions. The latter condition can be shown to be equivalent to the criteria that the derivative of the vertical velocity with respect to r, be continuous. Solution for the Upper Atmosphere (r >b) Here To is a slowly varying function of r, in fact we may take To constant. If the region is isothermal we will take TO = T s, and if the temperature should vary slightly, we will take the value of To at r = b to be Ts. For this region rewrite equation (26) in the form M" + M - (-) +- ( —) ) ) - 3 + R (r)| = O (68) l 4 a 2 a a r ~1/2 where we have set 3 A2 _ h (69) A + 3 and R(r) = (70) Since the temperature is slowly varying or constant we may neglect the remainder term R (r). Solving for M with R set equal to zero we have r M a -/Z~/2exp - /3 (r) dr (71) b -22 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T This is essentially the W. K. B. approximation. Of course for this to be sufficiently accurate, we require that the remainder R be sufficiently small, and this in turn holds provided that a does not vanish in, or close to the region b 4 r. In addition, at sufficient distances out there may be a turning point given approximately by the value for r for which j3 vanishes. However, we will assume* that the turning point is sufficiently far away from r = b, so that expression (71) is valid at least for the lower part of the upper atmosphere. We now have C (X) U (r, w) =- - exp - j (r)dr (72) rZ Po b ha where P = A - (73) - r which for the isothermal case has the explicit form,. <^ L)2 1. 2/4. 2) 9 2 1/2 - / gc r2 5 / cs The constant C 1 () is determined from the boundary conditions at r = b. Solution for the Lower Atmosphere (a z r - b) Because X 1 is the phase velocity of the gravitational wave we expect that there is some point in or just outside the interval (ay r b) for which This assumption breaks down when we approach the "so-called"l cut-off frequency. -23 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T (r2 / c2 _ X2 a2) vanishes. Hence the equation (67) has one or more singular points. Performing an order analysis it is seen that the term ha / r2 is small compared to the other terms for low frequencies (periods greater than three minutes). Thus a good approximation is given by the equation 3 '2 1 M" + M - 4 (-) +-(At +-) -A' A 0 (75) 4 4 / a 2 / and the general exact solution of (34) is r M = a-l/ 2e(r) C2 (X) + e -2.(x) dx (76) a r with ((r) = ~ Adx a r 2- g 1 T l2 2 t0dr (77) r \pj/2 c2 2 U(r, ) 1 et C2 (Xr))+ ae-2 (x) dx (78) PO a Determination of the Eigenvalue X = X"' (first approximation) First of all we must match the solutions for the two regions at r = b, from the boundary conditions that U, and Ur are continuous. -24 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T Using the continuity condition of Ur at r = b, we have from (72) and (78) C1(X) = e(b) C(X) + ae-2(x ) dx (79) a Before we employ the continuity condition of U at r = b we note that the corresponding derivatives have the form, 2 1 Y/g = Ur + - c2 - /3(r) b r (80) r 2 cs and U 1 K I g aei (r) U Ur L- r+ + 1 - a r b (81) r r r rz Po where in (80) we have taken the temperature constant To = Ts, and in (81) we have used the relation g 1 P a = A = + c 2 2 Po Thus for r = b we may equate the right hand sides of (80) and (81) giving U (b, o) [:(b) ) = ) r2 =b (82) r 2 - 1/2 Si 0 r = b Substitute the value for Ur (b, w) from (78), and rearrange terms to give b a(b, )e-2 (b) C )(X) - ae2-(x) dx -1 (83) Thus the values C 1 (X) and CZ (X) given by (79) and (83) are determined. However we must satisfy the remaining boundary condition namely that U must -25 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T vanish at r = a. Since the solution (78) in the region a < r - b, is uniquely determined apart from X, we require a condition upon X in order to satisfy the boundary condition. From (78) it is easily seen that for Ur to vanish at r = a we must have C2(X) = 0 (84) where we shall denote the real solutions of (84) by X. The first approximation to X* is given by equation (84). Higher Order Approximation for Ur (r, w) in a r b The problem is to find higher order approximations to the equation M" + M f(r) + g(r) = 0 (85) in the range a - r ~ b, and where for convenience of analysis we set 3 1 a f (r) = +- + Aa - A' - A2 (86) 4 a a ah g(r) = 2 (87) r The two independant solutions M 1 and M2 of M" + Mf (r) = 0 (88) are M1 = a-1/2 e(r) (89) r M2= -1/ 2 e j(r) a e-g~(x) dx (90) -26 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T As the first approximation to M we take r M1 -= -l/2 e(r ) C2 (X)+ f ae-2f(x) dx (91) a where C2 (X) is given by equation (83). To find higher order approximations of (85) we write the equation in the form M + f(r)M = -g(r)M (92) and treat the right-hand side as the inhomogeneous portion of the differential equation. Equation (92) may be solved by the method of variation of parameters giving b M = M1+ g(t) M(t) l (t) M (r)- M2 (t)M1 (r r)dt (93) r or expressing it in terms of vertical velocity we have b Ur (r, ) (r, ) +J U(r +(t, w) k(r, t) dt (94) r where h(t) 1l/2 (t) r k (r, t) = l ) e(t) + (r) a (x, >) e (x) dx (95) r2 pol/2 (r) t and Ur1 (r, L) is the first-order approximation. This is a Volterra integral equation which can be solved explicitly by successive iterations. -27 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T We may take the second order approximation b Ur2 (r, ) = Ur1 (r, w) + Ur (t, w) k(r, t) dt (96) r and obtain similar expressions for the higher order approximations. The next question that arises concerns the matching of the boundary conditions at r = b. However, it can be shown that the higher approximations, in fact, even the exact solution found by the iteration technique, are such that the boundary conditions, namely, that Ur, Ur be continuous at r = a, are automatically satisfied where for r ) b we have taken the solution given by equation (72) with C 1 (X) specified by (79). The remaining problem is to consider higher order approximations to X'. Let us consider the second order approximation only. We require that Ur vanish at r = a. From (96) and (78), we have C2(X) + f Ur1 (t, w) k(r, t) dt = 1/2 a2 Po (a) a i.e., b h ~Ut) 2 C2 (X) + t dt e t) () ) dx 2 f ( ) a a t e-2j(x) dx = 0 (97) Solving equation (97) will give the second approximation to '. From equations (B. 1), (B. 2) and (B. 3) of Appendix B, the radial component of the excess pressure P(r) is related to the vertical velocity compon ent by -28 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T Lue r - 1/2 o1/2 (r)P (r) = e r) e rr p (r) U (98) hence the second approximation to P (r) is on using (96) p-1/2(r) P(r) = ie-(r) + h ) e t) C + J ae ( ) dxdt t2 La (99) 8. Comments on the Accuracy of the Approximations: In order to show the accuracy of the approximations, we give two examples. First we compare our first and second approximations with the results obtained by Scorer. Scorer considered the case where the temperature gradient was constant in the troposphere and the stratosphere was isothermal, i. e., a r fb T = T1 (r-a) ( T)1 b r T = T5 where J. = b - a is the height of the troposphere T 1 = 286. 91~K is the ground temperature T = 229. 53~K is the stratosphere temperature The acceleration due to gravity g, was considered constant with the value 9. 806 x 10-3 km sec -2. In addition the constants,, ( are given ~/ = 1.403 = 9. 6137 km -29 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T Scorer solved the differential equation for a "modified" pressure function numerically, and hence obtained the eigenvalue K* which corresponds to our XA co. We compare the results between Scorer's theory and our approximate theory with the Table I below where X:'is in units km-1 sec. TABLE I 5000 2 First Approximation Second Approximation Scorer's _to A x ' "to x toValue of X 0 3. 1899 3.18947 3. 18950 1 3.1913 3.19148 3.19148 2 3.1929 3. 19355 3. 19356 3 3. 1945 3. 19570 3.19574 4 3. 1962 3. 19793 3. 19793 5 3. 1980 3.20025 3.20022 6 3.2000 3.20266 3.20271 7 3.2021 3.20517 3.20521 8 3.2045 3.20778 3.20781 9 3.21052 10 3.21338 3.21342 11 3.21638 12 3.21953 3.21956 13 3.22286 14 3. 22637 3. 22643 We immediately see that our first approximation gives the phase velocities (1/X\ ) accurate to three figures, and for the lower frequency modes, the error is one unit in the fourth figure. These errors really are insignificant compared to the errors introduced through the simplification of the atmospheric model, i. e., neglect of winds could produce a change in the phase velocity in the second figure. The second approximation is extremely accurate. -30 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T For the second case, we take an atmosphere whose region below 50 geopotential kilometers corresponds to measured values. Above the 50 geopotential kilometers we take the atmosphere as isothermal. We will consider the case where the measured values of temperature are given as follows (where for convenience of later reference we will call this atmosphere II): Altitude (geopotential km) surface 1.5 3.1 5.9 9.7 12.4 14.2 16.6 20 25 30 35 40 45 50 Temperature (degrees centigrade) 29 18 10 -5 -31 -53 -67 -81 -66 -55 -45 -34 -20 -5 + 7 with this temperature model, we calculate X* from the second approximation by setting equation (103) equal to zero. In Table II we compare the second approximations to X* to those calculated exactly (this is discussed below). TABLE II 5000 2 Second Exact Value Approximation to X of X * 0 3.1519 3.1528 1 3.1577 3.1582 2 3.1624 3. 1628 3 3. 1663 3. 1670 4 3.1697 3.1709 5 3. 1727 3. 1740 6 3.1753 3.1784 7 3. 1777 3. 1821 8 3. 1799 3. 1857 -31 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T It is seen that the approximate value of X* is accurate to four figures for small values of a, but for higher values the accuracy decreases to three figures. For this model of the atmosphere the accuracy is not as good as in Scorer's model. The main reason for this is that the error is proportional to (w,)4 for to / 0, where A= b - a, is the interval of integration (in this case = 49 km, in Scorer's model J = 9. 6137 kin). The approximate technique is good for low frequency waves and not too large non-isothermal intervals. However, for models of the atmosphere composed of alternating isothermal and constant temperature gradients sections, an approximate technique of this type should give good results. In the isothermal sections, analytical solutions in the form of exponentials can be found, and in the sections with constant temperature gradient, good approximate analytical solutions can be found. For large intervals and high frequencies, numerical methods are necessary. However, for very high frequencies asymptotic techniques may be used (see Friedlander, Ref. 14). To obtain a numerical solution we will return to the equation for the radial component P (r) of the excess pressure. From equation (B. 2) and (B. 4) of Appendix B we have d 121 p r d P P/Z] p + Fl(r) - cXZag p-1/Z = 0 (100) drth drwh s = where for the free waves q (r) is given by (18) with s = it. r Set = p1/ p(r) exp f Adx (101) a -32 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T with A given by (15), then (100) becomes h' 2 1 ga2 "y' + ' F-2 A- h - + h 2 J = 0 (102) 1h r 2 r2 Setting s = iw in (13) and using (101) we obtain 1/2?,' A p/ ic U, = - -h exp Adx o r h hI a thus the boundary condition at r = a, namely that Ur vanish requires Y' = 0 at r= a. and the boundary condition at a discontinuity in To (but not To) requiring that Ur, P be continuous, sets the condition that both ^T /h and V/ be continuous. The upper atmosphere was assumed isothermal and the upper boundary condition at the junction r = b, is (X-) g h (103) where j3 is given by (74) and ht, hs are the values of h for the lower and upper atmospheres at r = b, respectively. Equation (102) was solved numerically using the method of RungeKutta. In addition g was assumed constant, the term a 2/ r2 was set to unity and the term 2 / r neglected. 9. Simple Source Model: Consider the case where the surface l may be taken as a small sphere with centre ( Ro 0) and radius i. In addition we assume that locally the dis turbance produced by the source is spherically symmetric about (Ro, 0). If we take the limit as shrinks to a point, it is shown in Appendix A that we can -33 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T write in place of (34) P(a, 0, s) = p1/(a)G(a, 0;Ro, 0) I(s, Ro) where I (s, Ro) is a function of s representing the behavior of the source at (Ro, 0). 10. Method of Evaluation of p (a, 0, t): The excess pressure at (a, 0) may be evaluated from (104) by means of the inverse Laplace transform giving 1 b+i oo p(a, 0, t) =- p 1/2 (a) G(a, 0;Ro,0) I(s, Rc (104) ) ds (105) where the line Real s = b is to the right of all the poles of the integrand of (105). We shall assume that the position of the source (RO, 0) and the observer (a, 0) is such that the contribution of the modes other than the free wave modes is negligible. It can be shown that the other modes which satisfy the required boundary conditions at r = a and r = o will attenuate since the phase velocities are complex. We shall place the restriction that the observer should be in the geometrical shadow region i.e., 0 > cos-1 - iRo -34 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T In the illuminated region, the modes other than the "free" wave modes become important. In addition, appealing to known results of acoustic and electromagnetic theory, the Green's function cannot solely be represented by a residue series, but must be represented by a contour integral in the jt plane. Hence for the geometric illuminated region, other techniques than the modal expansion must be used. From (63) and (105) we have b+i oo p(a, 0, t) 2= 1 jioo etp,1/Z (a) I(s, R,) 2f (a. BRo'R Ods Z rib-io j= 2cosh(r oj) (106) where we have considered the contribution of the "free" wave modes only. Making the following decomposition 1 e_. j - e - 2 rtI j 2 cosh ( 7r w j ) cosh 7r j (107) equation (106) becomes 1 bb+i oo p(a, 0, t) = ri b-i oo n(s) est pl/2 (a) I(s, Ro) > fj e/J ds j=1 1 b+ioo 27r b-i oo n(s ) e- 2 r j eStpo1/2 (a) I(s, Ro) > f eos7r ds j = 1 cosh 7rA j (108) Using the relation Pi j - 1/2 (cos 0o) 1 for 0% = 0 (109) -35 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T and the asymptotic expression 2 Pij j-1/2 (-Coso) = sin - N[ ~~~7r sin 0.r(i j + 1/2) r( j + 1) 7 (110) cos i ij (7r - 0) --- (110) 4 for 0 < 0 < 7r together with the fact that for s in the right half-plane Real j > 0 (111) it can be shown that the second integral of (108) represents the contributions to the excess pressure of the modes that have travelled at least once completely around the earth. Hence if we consider only the contribution to the excess pressure pulse that arrives directly, we may take 1 b+i oo n(s) p(a, 0, t) = estpl/2 (a) I(s, Ro) ~ fj eT/ j ds 2-7ri b-ioo j=T We may now take b = 0 and set s = iw in (112), giving 1 +o n p(a, 0, t) = - r eiwt p01/2 (a) I(iw, RO) j fj e j dw 2 7r - ODj=l 3 (112' (113: For certain models of the atmosphere (in the sense of the variation of the atmospheric temperature To), there exists only one free mode below some critical frequency oc. For this case we obtain p(a, O, t) = 1 tc ewt- 7tr" Pol/2 (a) f'" (a, 0;Ro 0) I(iw, R ) dw c (e14) 27r - c (114) See footnote, page 171 of Reference 14. -36 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T where we specify the eigenvalue by /". For models of the atmosphere in which there exists a set of higher frequency modes with continuous spectrum in w greater than their respective cut-off frequencies, we must of course, include their integrals corresponding to the integral given by (114) for the low frequency mode. However, for the present, we will consider the excess pressure pulse produced by the low frequency gravity wave. From (64) we have Frr 1i11/2 = - i (jX a)2 +- ' i X* a for wX' a>>1 (115) From (110) we have 2 1/2 7" P/'" -1/2 (-cosO)v airo s Isin cos k'" a(7r- O) - (116) p -1/2~~ ~" aTrsin 4 hence we see that from (116) we may write p/2 (a) f* (a,R 0;R, 0) I(iu, Ro) = K(a, 0, R, u) co |I|X' a (- 0)4 (117) where from (61), (109) and (116) we have, on using the approximate relation, /o/ iwX a P., (a) Xn ) 1/'-2.^^ (RO,*) K = I(iw, R) (2o(a)) + A (118) The integral (114) can be written in the form 1 c Wc 4]iitt- ILA`t ar cos (a, 0, t) = -- eos ar WX' a (It -0) - Kdo (119) 2T -T 4 4 _W -37 —

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T Decompose K in the following manner K = Ke + Ko where Ke and Ko are the even and odd functions of w. Making use of the fact that X:" is an even function of w we can rewrite (119) in the form 1 Wc p(a, 0, t) = 7r cos[ (t - X*a7r) cos Owka ((r - 0) — Ke dw L 4 +- sJ sin [(t- X a7r) cos X* a ( 7r - 0-) - -K d Decompose p (a, 0, t) into two parts (120) p (a, 0, t) = pl (a, 0, t) + p2 (a, 0, t) where pl (a, 0, t) 1 2 c r o 0~ Ke cos t - WX a - dow eKo sin w4; a 4 Ko sin t - X':a 0- - dw L 4 *] i 27r (121) p2(a, 0, t) 1 27r wc Ke cos o J Ko sin o [t - wXka (2- - 0) +- do t - wX*a (2 - 0) +- dw L ~~~~4i i 27r (122) The physical significance of expressions pi (a, 0, t) and p2 (a, 0, t) is that p (a, 0, t) represents the pressure wave which has travelled directly

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T to the observer, whereas p2 (a, 0, t) is the pressure wave which has travelled through the antipode. For a simple point source at altitude zo = (Ro - a), the source function (given in Appendix A) is I(iw, Ro) = -iwpo1/2(Ro) V (123) where V is the volume of gas introduced. In this case Ke = 0, and (118)becomes /po (a) po (Ro) ''/2 K = Ko = i(,/) VQ(w) (124) \7r a sin 0wl J where for convenience we have set -2 0 (RO, x* ) Q (w) = a 7+AO r = a*| L h kJr= '" ' h X Xi` (125) Hence the pressure waves which have travelled both directly to the observer and by way of the antipodes are given respectively by 0~(a)Po c(Ro 12 1V WC 7r1 pl(a, 0, t) = ( - - J (wX)12 Q() cos t -wX' a +- d \ 2 7ras in 0 (7 )(L1 4J (126) o(a) po(Ro) V v we 3_ p2(a, O,,t)r = 0( a ~R ) - (wk')1/2 Q(w) cos wt-wXka(27r-0)+1dw 2 7r a sin 6 L 4 (127) The integrals in (126) and (127) may be equated to the real parts of Wc,-1/2 r1 (w(w ) Q(w) exp it t - wX*a0 +- du 0^.~~~~~~L 4 (128) -39 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T and Jc (t ) 1 /2 Q() exp i ot-wX*a (27r - 0)+ 3 d[ (129) o respectively. For large ranges (a0 and (ZTr - 0) a sufficiently large) these integrals may be evaluated by the method of stationary phase. If wo is the value of w for which t - a((w')' = 0 (130) where the prime indicates the derivative with respect to w, then the stationary phase technique gives the asymptotic approximation to (128) aI Q(Wo) X 9 cos t - o X' a0 (131) Similarly if wo is the value of U for which t - a(27r - 0) (wXX*)' = 0 (132) then the stationary phase technique gives the asymptotic approximation to (129) - - Q(Wo) - cos wo t-woX' a (27ir-0) + -(133) a(27z-0) (00)" L_ ow P However the stationary phase technique for (128) is not accurate unless (wX)" / [(CWX)' 3/24 Jf and (wo) iv / [(wX)' 2 a0. It has been shown by Penney (Ref. 6 ), that for most models of the atmosphere under consideration, this critique fails when the range is of the order of 1000 km or less. For -40 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T these cases it is best at present to perform numerical integration to evaluate (128) in the manner that is suggested in Reference (6). In addition, even at very large ranges, the stationary phase technique fails for the initial part of the received pressure pulse. We shall give here an asymptotic technique for calculating the initial part or head of the pulse. We will consider integral (128). Define the following series x'() = ( 1 2n (134) n= 0 (Xk)1/2 Q (w) = o/ Q (o) j7 dnc Z2n (135) n=0 where do = 1, and the time interval -= t - a 0. (136) By the stationary phase technique, the condition 't = 0 or t = a0Xo would indicate the beginning of the pulse. However for l|l small, the stationary phase technique fails and we shall show that r = 0 does not represent the beginning of the pulse, but a point near the peak of the first compression. The integral (128) may now be written in the form X012 Q(o) f Wl/2 dnczn 1 - iaO(w5X +d 7X3+...)+..] o n =0 exp i uwZ - c3X1a0 +-dw d (137) exp [[~o~'-~03klae +4Ida0 -41 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T Making the substitution w = O(X a0e)-1/3 (138) q = (X a) -1/3 (139) and defining the integral Fn(q) = J 1-/2+ne e i[.q -3 do (140) o expression (137) has the asymptotic form 0 Q (o) ei7/4{ Fl (q) + 1 a /3 F3(q) + — F (q) +.. \Xia&(q)+0 (XiaO)2/3 (Xla0)4/3 iaO iaO kla05/3 X2 F6 (q) - iX + d.. (141) (X "1aO )a5l/3 (1 )7/3 3 +d F8 ()+J (141 In Appendix C, Fo (q), F1 (q) and F2 (q) are expressed in terms of Airy integrals, and Fn (q) for n ) 3 are expressed in terms of Fm (q); m = 0, 1, 2. Thus expression (141) may be represented in terms of Airy integrals (which are tabulated, Ref. 13). After much algebraic manipulation, the realpart of (141) has the following asymptotic form 2Q (o) - Ai Ai C1 + Ai2 C2 -[(Ai )2 + y (Ai)2 C31 (142) 3X a0 L where Ai (y) is the Airy integral, with Ai being its derivative. The argument y is given by Y - (12Xa 0)1/3 (143) -42 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T The functions C i (r) have the forms = 1 20( 2x 8t2 - = 1+(l1 ) 1 (121a0)+ 2 (12X1 a 0) 3 3X1 ((12k1 a0)22 C1 (T) 7[3 +dlj i2 xi c (r) 1 (12 xl a e)2/3 - 7 X 4 3[3+ 6X (1ZXa 0)5/3 k I C3 () 10 (12 k1 a )4/3 5 3 + 12 a )7/3 (12 X1 a e)7/3 2 11[X3 +dl X 2 4X2 6X1 (12X1 a 0)4/3 3X1 16[X3 + dl2] 3X1 (144) For the real part of integral (129) one may derive in a similar manner the asymptotic expression 23, 1 -2Q(o) -- ----- ()Al I (y)Bi (y) + — 3v 3Xa(27r - 0) a- 2 7r lw/ (145) where the argument of the Airy integrals in this case is given by with 12 X1 a (2r -_ )]1/3 = t - a(27r - 0) X (146) (147) Because the distance around the earth by way of the antipodal route is extremely large, the higher order terms in the asymptotic expansions are unnecessary, -43 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T except for very large T'1 in which case, the stationary phase technique is valid. The directly received pulse p1 (a, 0, t) is calculated for Scorer's model of the atmosphere for explosions both on the ground and at a height of 9. 6137 km, and the pulse forms are plotted in Fig. 1. The integral representing the pulse was calculated using the stationary phase technique for r = t - a X o > 160 sec. and using the asymptotic expansion involving the Airy integrals for r' - 160 sec. Even though the use of the stationary phase technique for the case when a = 3600 km (the distance from the source) is somewhat a delicate matter, there is good matching of the results obtained by the two techniques. The complete pulse forms are not given, the dashed line indicating that the remaining portion of the pulses are not given. The transit time of the first crest (given approximately by t = a X ) is 3 hr. 11 min. and 5 hr. 19 min. for distances 3600 and 6000 km respectively. In addition, the received pressure pulse p 1 (a, 0, t) is calculated for atmosphere II for a range of 7000 km for explosions on the ground and at a height of 39 km. A table of the intensity function required in the evaluation of the pressure pulse integral, and the X's are given below. The value of the X's for the higher values of o's are not calculated, The stationary phase technique was used for > 260 secs., the other asymptotic technique for Z: A 260 secs. The main body of the pressure pulses observed at ground level at range of 7000 km from both an explosion on the ground and an explosion at 39 km are plotted in Fig. 2. The transit time of these pulses is roughly 6 hr. 8 min. -44 -

l, IAA \A AA.~~~~~~~ (a) 4 3 2 1 I \. (b) 1 Tq mr I 1 \ I I I 0 5 (Minutes) 10 5 10 (Minutes) 15 d q 1 -H 3 2 - 1 - I (c) AM AA. (d) 3 2 1 NJ I Co I co I I I I I ~ -4 0 z 1-4 ' ' I I I I 0 10 I I I 0 5 (Minutes) 5 10 (Minutes) 15 FIG. 1: THE PRESSURE PULSE AT THE GROUND AT A DISTANCE OF 3600 km (a) and (c), and 6000 km (b) and (d) FOR SCORER'S ATMOSPHERE, FOR AN EXPLOSION ON THE GROUND (a) and (b), AND AT A HEIGHT OF 9. 6137 km (c) and (d). ONE UNIT OF AMPLITUDE CORRESPONDS TO 0. 614 u BARS, PER 1 km OF GAS RELEASED.

1. AA A A^ 0 5 10 15 20 25 minutes ca rTl -1.0~o g _.02 (b) o I -, -.01 -.010 -.02 Z FIG. 2: THE HEAD OF THE PRESSURE PULSE AT GROUND LEVEL AT A DISTANCE OF 7000 km FROM AN EXPLOSION ON THE GROUND (a), AND AT A HEIGHT OF 39 km (b), FOR ATMOSPHERE MODEL II. ONE UNIT OF AMPLITUDE CORRESPONDS TO 1 p BAR, PER 1km. OF GAS RELEASED. (The tail of the pulse which extends for a very long time interval beyond the dashed lines is not given. )

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T I 5000 2 X( () 1 Q ( Ro) o (a) Q Ro) Explosions on the ground Explosions at 39 km I 0 3. 15283 10.71.1856 1 3. 15819 10. 83.1848 2 3. 16277 10.90. 1862 3 3. 16695 10. 68.1843 4 3.17091 10.42.1829 5 3. 17472 10.09.1807 6 3.17844 9.74.1783 7 3. 18209 9.37.1758 8 3.18569 8.99.1730 9 3.18923 8.61.1699 10 3.19274 8.25.1671 11 3.19620 7.86.1635 12 3.19962 7.54.1609 13 3.20301 7.19.1576 14 3.20635 6.86.1542 15 3.20965 6.51.1500 16 3.21290 6.19.1462 17 3.21611 5.89.1424 18 3.21928 5.60.1387 19 3.22239 5.30.1344 20 3.22545 5.03.1304 21 3.22846 4.77.1263 22 3.23142 4.51.1221 23 3.234 33 4.30.1187 24 3.237 18 4.09.1151 25 3.239 9; 3.87.1111 One important difference between explosions on the ground and at an altitude of 39 km is that the gravity wave portion of the pressure pulse is less attenuated for the latter case. As expected the amplitude on the ground of the pressure pulses is less for high altitude explosions. For atmosphere II, the ratio of amplitudes of the pressure pulses observed on the ground due to an explosion at 39 km and -45 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T on the ground is. 017, where the same volume of gas is released in both cases. An estimate of the ratio of the amplitudes of the gravity wave portion of the pressure pulses is given by 2 exp - g dz Z co(z) where z, z 1 are the altitude of the explosions. However this holds for idealistic explosions characterized by the sudden release of a specific amount of gas. 11. Comments: The pulse forms in the previous section have been obtained from an idealistic source model. The source representation could be more generalized by taking a different time dependence (like the ones mentioned in Ref. 6) rather than the delta function. Another question arises as to the relation between the volume of gas introduced and the total energy released by an actual explosion, for instance a nuclear explosion. Penney et al do give some estimates for explosions on the ground. However, explosions at various altitudes will require different constants of proportionality. A problem arises with regard to the upper boundary condition, namely that the total kinematic energy be finite. For an isothermal stratosphere or upper atmosphere, the excess pressure has the following functional dependence on altitude z yg 2es exp [zS-46 --46 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T where the flat earth approximation is taken for simplication by equating r to a, and taking g constant in the expression for 13 given by (74). For Scorer's model 'g / (2 c 2) is greater than 3 for the gravity wave. Hence the excess pressure increases exponentially with altitude. At high altitudes the hydrodynamics equations can no longer be linearized, as the excess variables approach and even become greater than the magnitude of the unperturbed variables. For atmosphere II, g/(2 c 2) is greater than 3 for 5000 co2 less than 13. In this case the excess pressure of the gravity wave increases exponentially only for the lower frequencies. This suggests that when taking an isothermal model for the upper atmosphere, it probably is best to start the isothermal level above 100 or 110 km, in which case the temperature may be taken very large, say 350~ K or higher. For this case one would expect that the excess pressure of the gravity wave attenuates with increasing altitude except for the extreme lower frequencies. However, at high altitudes an additional factor becomes important, namely viscosity. This should be included in the model of a very warm isothermal layer at high altitudes. Winds have an effect upon the gravity wave although not so much as on the high frequency acoustical waves. Apart from jet streams, there are good strong, steady winds at high altitudes. For instance, in the winter there are winds of magnitudes of 50 m/sec and 100 m/sec at heights 40 and 60 km respectively around latitudes 30~ to 60~. (See Ref. 15. ) These winds can cause an appreciable change in the phase and group velocities of the gravity wave. -47 -

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T In Part II, the effect of winds and better models of the upper atmosphere will be investigated. Acknowledgement The author wishes to thank D. Brown, D. R. Hodgins and L. Evans for the computational work that was performed. -48 -

THE UNIVERSITY OF MICHIGAN 2866 - 1 - T APPENDIX A Simple Source Representation: Given the inverse Laplace transform of the excess pressure in terms of an integral over a surface, namely P(a, 0, s) = -sp1/2 (a) \ G(a, 0; r, 0) n U(r ) po/2 (rO) dSo (a) 2 (r)P(r) n. i + n-L G dS - ^/2 (a) \ Po 2 (ro)P( r) - o }0 - -0 0 where the variables of integration are (ro, 0 ), we will consider the case where E is a small sphere with center (Ro, 0) and radius c. We will assume in addition that the disturbance (characterized by the transform of the normal velocity to the surface and the pressure) is spherically symmetric about (Ro, 6), in which case n U (r0) and P(r o) are functions of s and c only. We shall assume that, the velocity and the pressure satisfy the acoustic equation at the surface Z, in which case U and P are proportional to -2 -1 -2 and c respectively. Retaining only the dominant terms for small e we have P(a, 0, s), - spol/2 (a) p 1/2 (Ro) G(a, 0; Ro, 0) n. U dSo — 0 0 A-1 A - 1

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T / po/2 (a) G(a, 0; Ro, 0) I(s, Ro) For an instantaneous velocity source at r =(R 0, 0) at time t we obtain n.U (rO, s) n-U(ro) est and J n-U(ro) dSo = V is the volume of gas introduced Thus for a simple point source at time to I(s, R0) = - sp/2 (R0) V est and for to - 0+, this becomes I(s, Ro) = - spo2 (Ro) V. In the special case where the explosion is on the ground, the surfaces is a hemisphere, and V will be the volume of gas passing through it. A-2

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T APPENDIX B Radial Equation for the Vertical Velocity: For the convenience of analysis set iWp /2 U (r, 0, w) = (r) Pi _ (cos ) (B. 1) then using the relationship po1/2 0(r) Pi - Y2 (cos ) = P(r, 0, w) (B. 2) together with radial component of equation (13), we obtain 0' + AP h — - (B. 3) h But we have from Eqns. (47) and (64) d (2 d0)A + t(ro) w2 X2a p 0 (B. 4) dr eh dr/ Hence from (B. 3) substitute- J + A h-1] in place of ' h in (B. 4) and repeating the procedure after differentiation we obtain dr2 X W2 r2 Ar2 X - - A x - 02 x2 a2 + -- 0 (B.5) dr c2 B-

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T Setting r2X = -p (B. 6) (B. 7) a, 2 r -i a2 ] This reduces to V' - AI = aO (B. 8) Eliminate 0 from (B. 8) and (B. 3) to obtain TI j) a'. + A2 A + A -__ ha r2 (B.9) (B. 10) Setting p = a1/2 M we obtain 2 M" + M - 3 a?' 4 au +a, 2a + -A' - A2 a + h aCY where now from (B1), (B.6), and (B.) we have where now from (B. 1), (B. 6), and (B. 10) we have (B.11) Ur- r 2 2 2 2;1 M(r) r2 Pi/ - 1/2 (Cos 0) (B. 12) po(r) B-2

THE UNIVERSITY OF 2886 - 1 - T APPENDIX C MICHIGAN Evaluation of the Integrals: 0O Fn(q) -I O ' /2 +xp i wjq - o d2 w [. From Franz (Ref. 12) we have the relation 00 0 3 / -1 /2 - Ce a 2 e d y /3 2 7r/ 3 4 - i 7/3 4 41/3) (C.1) where 1 A(x) - 2 00 -00 d i xr df e -73] (C.2) However we may express A(x) in terms of the usual Airy integral A.(x), the relation by -1/3 - 1/3 A(x) - 3 7r A (-3 1 x) (C. 3) where o0 0 o (C. 4) Hence the integral (C. 1) may be given in the better known notation O 3 -1/2 aa - (y 2 e 3/2 5/3 -1/6 2i7r/3 () do =- TT 2 3 Ai ( Ai(y) (C. 5) C-1

THE UNIVERSITY OF MMICHIGAN 2886 - 1 - T -i7r/3 )-1/3 y = -e c(12) (C.6) where Performing the transformation i7r/3 a = qe i c = we I (C.7).i /6 Equation (C. 5) may be written in the form F (q) = 0 3/2 5/3 -1/6 -iir/12 A ( 2iir/3 Y3 e A-(y) Ai (e y) y - = - 1/3q y~~~~ (C. 8) (C. 9) with From Miller (Ref. 13) we have A (e 2i7r/3 Y) A. (ei /3 y) - 1 1 i 4r/3 2- e 2 Ai() - i Bi(y) (C. 10) with 0 Bi(Y) _ o L exp 1 3 t 3 + yt] + sin 1 t3 + 13 ytj dt (C.11) hence (C. 8) becomes 3/2 2/3 -1/6 iir/4 Fo(q) =- 2 3 e 0~~~~~ Ai(Y) - iBi(y) (C. 12) C -2

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T The other integrals Fn(q) may be evaluated from F (q) by differentiation with respect to q, ie! n d F (q) Fn(q) = i dq (C. 13) -n/3 n dF() (12) i dy In particular we obtain Fi(q) = 2 3 e 2 ei 4 2Ai(Y) Ai(y) - i Ai(y)Bi(y) 1~~~~~[ + A.(y) B'i(y which on using the Wronskian relation A. Bi AB. - 1 -1i Bi r (C.14) reduces to 3/2 -1/2 i37r/4 7r 23 e Fl(q) = Ai(y) Ai(y)' - i 'Ai(y)B (y) + (C.15) 1 ' L 2lrJJ (.15) In addition we obtain C-3

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T 1/2 F2(q) - 2(73/3) -1/3 i57r/4 (12) e (A,)2 I 2 + y(Ai) - i (YAiBi + Ai B) The remaining Fn(q) may be expressed in terms of Fo, F1, F2 by the relation (c.16) -n/3 Fn(q) - aCn(12) i FO(q) + n(12) (n i - 1) F1 (q) + 'n(12) (n - 2) i F2(q) where the functions a n(q), fn(q) and f n(q) are given in the table below n%^n~n n na n - On - n 3 2 4y 0 4 0 6 4y 5 8y 16y2 10 6 28 80y 16y 2 3 7 32y 108 + 64y 112y 8 288y 672y2 (220 + 64y ) _,... l.. C-4

THE UNIVERSITY OF MICHIGAN 2886 - 1 - T REFERENCES 1. R. Yamamoto, "A Dynamical Theory of the Microbarographic Oscillations Produced by the Explosions of Hydrogen Bombs," Journal of the Meteorological Society of Japan, Series II, 35:5, (1957). 2. C. L. Pekeris, "The Propagation of a Pulse in the Atmosphere, " Proceedings of the Royal Society, A 171, 434, (1939). 3. C. L. Pekeris, "The Propagation of a Pulse in the Atmosphere, Part II," Physical Review, 73:2, 145, (1948). 4. R. C. Scorer, "The Dispersion of a Pressure Pulse in the Atmosphere," Proceedings of the Royal Society, A 201 137, (1950). 5. L. A. Diki, "Acoustical and Gravity Vibrations in the Atmosphere," Izvestiia A.N. SSSR; ser. geofiz, (8): 1186-1194, (1959). 6. J. N. Hunt, R. Palmer, and Sir William Penney, "Atmospheric Waves Caused by Large Explosions," Phil. Trans. Roy. Soc. London, 18 February 1960. 7. J. L. Taylor, "An Exact Solution of the Spherical Blast Problem," Phil. Mag., Vol. 46, 317, (1954). 8. A. S. Ramsey, A Treatise on Hydromechanics, Part II, C. Bell and Son, (1949). 9. K. 0, Friedrichs, "Criteria for Discrete Spectra," Symposium on Electromagnetic Theory, p. 85, Interscience, (1950). 10. K. 0. Friedrichs, Courant Anniversary Volume, Interscience, New York, pp. 145-160, (1948). 11. E. C. Titchmarsh, Eigenfunction Expansion, Part I, Oxford, (1946). 12. W. Franz, "On the Green's Function of the Cylinder and Sphere," Z. fur Naturforschung, 9a, 705-714, (1954). 13. J. C. P. Miller, "The Airy Integral," British Assoc. for the Advancement of Science, Math. Tables Part - Vol. B, Univ. Press, Cambridge, (1946). 14. F. G. Friedlander, Sound Pulses, Cambridge Univ. Press, (1958). 15. R. J. Murgatroyd, "Winds and Temperatures Between 20 km and 40 km — a Review, " Quart. Journ. of Royal Met. Soc., 83, No: 358, 417-458, October, 1957.

UNIVERSITY OF MICHIGAN 3 9015 03530 0337