SONG_finalonline:Layout 6 ANNALS OF GEOPHYSICS, 57, 5, 2014, A0542; doi:10.4401/ag-6474 A0542 Acoustic ray tracing in the atmosphere: with gravitational effect and attenuation considered Yang Song*, Yuannong Zhang, Chen Zhou, Zhengyu Zhao Wuhan University, Department of Space Physics, School of Electronic Information, Wuhan, China ABSTRACT An acoustic ray tracing model is developed to take into account the im- pacts of gravitational field and realistic atmospheric attenuation. Ray tracing equations are deduced from the real part of the dissipative dis- persion relation, while the acoustic attenuation coefficient and growth rate in a stratified moving atmosphere are deduced from the imaginary part of the dispersion relation. To account for the non-isothermal effect and realistic attenuation, the buoyancy frequency and the cut-off fre- quency are substituted by the values in the slowly varying atmosphere, and the attenuation coefficient is corrected by the realistic absorption. In the validation by numerical experiment, the ray trajectory obtained by this ray tracing model agrees well with the result calculated by the FDTD method. It is shown that the acoustic trajectory can be accurately pre- dicted by this ray tracing model. The numerical results for 5 Hz acoustic waves show that in the stratospheric ducting the gravitational effect plays a leading role while the attenuation effect could be neglected. But for the thermospheric ducting, the contribution of the absorption becomes more important. 1. Introduction It is widely considered that the acoustic wave is one of the essential forms of energy transportation in the atmosphere. Depending on the profiles of the tem- perature and wind, waves could be ducted through the troposphere, stratosphere, and thermosphere [Drob et al. 2003]. This phenomenon, called long range sound propagation (LRSP) [Kulichkov 2004], is frequently ob- served at infrasonic frequencies because of the low at- tenuation in the lower atmosphere [Bedard and Georges 2000]. Consequently, natural and artificial sources could be detected and monitored through the network of in- frasonic stations, such as the International Monitoring System [Le Pichon et al. 2009]. In addition, as the prop- agation is controlled by the background profiles, it is also possible to study the properties of the atmosphere by analyzing infrasonic signals, such as probing the high altitude wind [Le Pichon et al. 2005]. Therefore, it is im- portant to comprehensively understand the propaga- tion properties of acoustic waves. Ray tracing is an efficient way to study the propa- gation of acoustic waves in the atmosphere. It can be used to calculate the trajectories of acoustic waves as well as a number of propagation parameters including travel time and arrival angle. With the ray path avail- able, the location and intensity of the fluctuations can be predicted [Pilger and Bittner 2009]. Since the early model developed by Barry [1963], enormous progresses have been achieved in ray tracing studies. There are three classical ray tracing models, namely HARPA, WASP-3D and Tau-P, all of which have been widely ap- plied to study the sound propagation [e.g., Le Pichon et al. 2010]. HARPA is derived from the dispersion rela- tion in a stratified moving atmosphere on the basis of Hamilton equations [Jones et al. 1986]. WASP-3D [Dessa et al. 2005] is also a Hamiltonian solution extended from 2D Cartesian ones [Virieux et al. 2004], but it is deduced from the Eikonal equation. Tau-P model [Gar- cés et al. 1998] is developed for computation of travel times as a reformulation of the original method pre- sented by Buland and Chapman [1983]. All of these three models are developed to study the atmospheric wave propagation through a windy, non-isothermal, non-dis- sipative atmosphere. Theoretically, the high frequency approximation adopted by the above three models confines the ray trac- ing modeling to acoustic signals above 5 Hz for which the gravity influence is neglected. Moreover, the at- tenuation of acoustic wave is not considered well in the three models. However, in the atmosphere acoustic waves can extend down to much lower frequencies ~0.02 Hz [Drob et al. 2003]. In addition, there are many types of signals in the 0.02-5 Hz frequency band, due to the radiation by some natural or artificial sources. For instance, signals with frequencies from 0.5 to 5 Hz Article history Received December 10, 2013; accepted September 2, 2014. Subject classification: Acoustic ray tracing, Dissipative dispersion relation, Gravity effect, Attenuation. were observed in the cases of large bolides [Hedlin et al. 2010], meteorites [Le Pichon et al. 2013], volcano [Kanamori et al. 1994, Evers and Haak 2005] and explo- sions [Stevens et al. 2002, McKenna et al. 2007]. Under the influence of gravity, the propagation of acoustic waves becomes more dispersive and anisotropic as fre- quency decreases. For modeling the propagation of acoustic waves below 5 Hz, the high frequency ap- proximation is no longer valid and the effect of gravity needs to be taken into account. Damping of acoustic waves in the realistic atmosphere is another important ingredient of acoustic propagation modeling [de Groot- Hedlin et al. 2011], which however is often ignored in ray tracing modeling. Based on the calculation of the speed and attenuation at altitudes up to 160 km devel- oped by Sutherland and Bass [2004], the attenuation ef- fects on the propagation of infrasound is studied without considering gravitational effect [Bass et al. 2007]. The damping effect of low frequency waves is not significant near the ground, but it increases remarkably at higher altitudes. In addition, the accumulation of attenuation effects along the long-range ray path is nonnegligible. Therefore, in order to model the ray path of acoustic waves precisely, the wave attenuation also needs to be included. Ray theory for lossy media has been discussed in radio wave propagation by Jones [1970]. To take losses into account, the complex phase refractive index in- stead of the real index defined in the ideal media is used, so that the ray path can be calculated using the extended Snell’s law or Haselgrove’s equations. The ex- pression of such a path is in complex space, so only those rays ended in real space are physically significant. Another way is developed by extending Hamilton’s equations for lossy media under the constraint that the space, time and group velocity should be real [Sonnen- schein et al. 1997]. When studying the acoustic wave propagation, especially the infrasonic propagation, both the effects of attenuation and the gravitational dis- persion should be taken into account. Other than using the complex index, it could be a better way to obtain the ray equation through the Hamilton’s equations by employing the dissipative dispersion relation as an ex- tension. In the present study, an acoustic ray tracing model in the realistic atmosphere is developed on the basis of analyzing both imaginary part and real part of the dis- sipative dispersion relation. In Section 2.1, the dissipa- tive dispersion relation is derived from the governing equations of an isothermal dissipative atmosphere. In Sections 2.2~2.4, the acoustic attenuation coefficient and growth rate in the moving atmosphere are deduced from the imaginary part of the dispersion relation, while ray tracing equations are obtained from the real part of the dissipative dispersion relation. Modifications for the parameters are presented in Section 3. The buoyancy frequency and the cut-off frequency are sub- stituted by the values in the slowly varying atmosphere. The acoustic attenuation coefficient is corrected by the theory of absorption in the atmosphere up to 160 km [Sutherland and Bass 2004, 2006], aiming to include the non-isothermal effect and the realistic atmospheric at- tenuation. This ray tracing model is validated by the comparison with FDTD solution in Section 4. Numer- ical examples are shown in Section 5 to study the influ- ences of the gravitational effect and attenuation, followed by the conclusion in Section 6. 2. Inclusion of gravity and absorption into ray trac- ing equations In order to introduce gravity influence and dissi- pation effects into acoustic ray tracing theory, we start with the compressible fluid equations with molecular viscosity and thermal diffusivity in three dimensional Cartesian systems [Beer 1974]. (1) where t, p, s and T are particle density, pressure, en- tropy, and temperature respectively, V = (u,v,w) is the wind velocity, is the Stokes derivative, g is the gravitational acceleration, h and p are shear and bulk viscosity. For the convenience of formula deriva- tion, the viscous term is neglected in following calcula- tions. With the knowledge of the thermodynamic rela- tionship, the entropy and temperature in the atmos- phere can be given by (2) In Equations (2), c = cp/cv is the ratio between spe- cific heat at constant pressure and constant volume, and R is the universal gas constant. It is assumed that the basic quantities such as t, p and V can be expressed as the sum of the ambient value (with subscript 0) and the perturbation caused by the passage of a sound wave (with subscript 1), i.e. (3) As the acoustic wave in the atmosphere is an irro- tational wave, the variable V satisfies ∇∇·V1=∇ 2V1. VDt D t $2 2 d= + V V g VV t Dt D p T Dt D s 3 0 0 2 $ $ 2 2 d d dd d t t t t p h t h + - + = + = = + ^ h Z [ \ ] ]] ] ]] ` j const ,lns c p T R p v yt t = + =c m , , V V Vp p p0 1 0 1 0 1t t t= + = + = + SONG ET AL. 2 3 With the approximation of a stratified atmosphere, Equation (1) can be simplified and linearized to yield where V0 = (u0,v0,0) is the wind profile of the atmosphere, is the adiabatic sound speed with T0 as the ambient temperature. 2.1. Dissipative dispersion relation Assume that the acoustic wave has a plane wave solution in the form (5) where P, R, X, Y and Z are polarization terms, ~ and k=(kx,ky,kz) are angular frequency and wave number components, respectively. Substitution into Equation (4) yields a set of fifth-order algebraic equations for (u1,v1,w1,t1/t0,p1/p0), which can be expressed in the matrix form as (6) In the coefficient matrix of Equation (6), X = ~ − V0 · k, the Doppler-shifted frequency, is the coefficient of kinematic viscosity of air, H = c2/cg is the scale height. According to the Cramer’s rule, for a non-trivial solution of Equation (6), the determinant of the coefficient matrix must vanish. Therefore, the dissipative dispersion relation becomes (7) Equation (7) is complex because the inhomogene- ity and dissipation of the atmosphere results in complex wave number components. Obviously, if the absorption is neglected, Equation (7) reduces to the dispersion re- lation of acoustic gravity waves [Yeh and Liu 1974]. Since the wave vector k is complex, it can be written as (8) where k' = (k'x,k'y,k'z) is the real part of the wave vector, known as the wave vector in an ideal medium, and k" = (k"x,k"y,k"z) is the imaginary part of the wave vector that can characterize the attenuation and growth of wave amplitude in the inhomogeneous atmosphere. Subse- quently, the dispersion relation Equation (7) can be de- composed into two equations, i.e., the imaginary part, Equation (9) and the real part, Equation (10). (9) (10) The first term on the left side of the dissipative dis- persion relation, g2 (c − 1)(k2x + k 2 y), plays an important role in the real part of the dispersion relation, as it rep- resents the influence of the gravity that generates grav- ity dispersion. This term is discarded in the derivation of Equation (9), because the imaginary part of the dis- sipative dispersion relation is mainly associated with the attenuation relation between frequency and imaginary wave vector k". However imaginary wave vector k" represents the changes of wave amplitude, we keep it in the derivation of Equation (10). The imaginary wave vector k" can be derived from Equation (9) in order to calculate the amplitude varia- tions along a ray path. By substituting k" into Equation (10), the ray equations can be determined on the basis of Hamilton equations. 2.2. The attenuation and growth of acoustic waves As acoustic waves propagate in the atmosphere, the wave amplitude usually drops due to the attenua- tion effect. Growth in the wave amplitude may occur because of the atmospheric density gradient in the ver- tical direction. Thus, in the dissipative stratified atmos- phere, the imaginary wave vector can be considered as a combination of omnidirectional attenuation and growth in the vertical direction. Assuming the attenu- ation efficient and growth factor as a and b respectively, we obtain that (11) By substituting Equation (11) into Equation (9) and (10), the imaginary part and real part of the dissipative dispersion relation can be rewritten as 0 Dt D u x p Dt D p Dt D p g Dt D x z w Dt D p c Dt D z p c z w u v y v w z w u y v z w 3 4 3 4 3 4 0 0 0 1 1 0 0 1 1 2 1 0 0 1 1 1 2 1 0 1 0 1 0 1 0 1 0 1 0 1 0 2 0 1 0 2 0 1 0 2 12 2 2 2 d 2 2 d 2 2 2 2 2 2 2 2 2 2 2 2 dt p h t p h t t p h t t t t t t t + = + + = + + + = + + + + + = - + - = Z [ \ ] ] ] ] ]] ] ] ] ] ]] ` ` ` c j j j m exp p P p R X u Y v Z w A i t k y k y k zx y z 0 1 0 1 1 1 1 t t ~ = = = = = = - - -^ h6 @ ,VDt D t 0 0 $2 d= + 3 4 0g p h t= +` j 1 0 g k ik i k c k ik i gk yx z 2 22 2 2 2 2 2 c g g c X X X X - - + + - + - = ^ ^ ^ ^ h h h h k k ki= +l m 2 0k k kc gkz0 2 2$ g cX+ - =l m l l yx1 2 2 0k k k k k k g k gk k c z 2 22 0 2 2 2 2 $ $c g g c X X X - + + - + - - - = l l l m l m l m m ^ ^ ^ ^ h h h h6 @ , ,k k k k k k k k x x y z z a a a b= = = + yk m l l m l l m l l c RT0 0c= ACOUSTIC RAY TRACING IN THE ATMOSPHERE (4) 1 i k ik i k ik i k ik H H g i i igHk igHk igHk g i u v w p p 0 0 0 0 0 0 0 0 1 0 0 0 0 x x y y z z 2 2 2 2 1 0 1 g g g c c t t X X X X X X + - + - + - - - - - - - = p r q q q q q q qq p r q q q q q q q q t v u u u u u u uu t v u u u u u u u u (12) The attenuation coefficient and growth rate can then be solved from Equation (12) as (14) (15) In Equation (14), minus sign means the decrease of wave amplitude; Doppler-shifted frequency X represents the impact of wind on the attenuation by changing the propagation distance in a unit time. In the windless at- mosphere, X reduces to ~, and this coefficient becomes the classical absorption coefficient [Beer 1974]. Equation (15) indicates that an exponential amplification or decay with height can be caused by the inhomogeneity of the atmosphere in the vertical direction. 2.3. Final form of dissipative dispersion relation Substituting the solved attenuation coefficient a and growth factor b into Equation (13), the real part of the dissipative dispersion relation becomes the final form of dispersion, written as (16) the term , which is the contribution of the gravitational effect, where ~2g = (c− 1)g2/c2, ~g is the isothermal Vaisala-Brunt frequency; and ~a = c/(2H) is the acoustic cut-off frequency; the term is the con- tribution of attenuation. If both gravitational effect and attenuation effect are neglected, the dispersion will be- come the classic acoustic dispersion relation (17) According to the theory of atmospheric waves, the waves whose frequency is greater than acoustic cut-off frequency are identified as acoustic waves, and the waves whose frequency is less than Vaisala-Brunt fre- quency are the internal gravity waves. The dispersion relation Equation (16) can be employed to describe both acoustic waves and internal gravity waves, but this work mainly focuses on the propagation of the acoustic wave, the internal gravity waves will not be discussed. Actually, when the dispersion is employed to ob- tain the ray tracing equations, the accuracy of the ray tracing equations depend on the dispersion relation equation. So the gravitational term Gr and attenuation term At are of great importance to the dispersion rela- tion. It can be discussed in two limiting cases that are low frequency case and high frequency case. For the acoustic wave with low frequency which is close to the acoustic cut-off frequency (~ → ~a), the dissipative effect is so small that could be neglected, meanwhile the gravitational effect is enhanced, so the dispersion relation is transformed as (18) Equation (18) is the classic dispersion relation of acoustic gravity wave. As to the acoustic waves with high frequency (~ >> ~a), the dispersion equation becomes (19) Different from the classic acoustic dispersion rela- tion, Equation (17), the attenuation term At is taken into consideration, because as the frequency becomes higher, the attenuation becomes stronger. In the high frequency approximation, the gravitational dispersion term , but the scale height of a mean density variation still exists. The term ~2a can be taken as the contribution of the inhomogeneity of at- mosphere density. Both limiting cases demonstrate that the inclusion of the gravity and attenuation are reasonable enrich- ment for the classic acoustic dispersion relation. The ef- fect caused by gravity and attenuation can be studied by the group velocity which can be derived as where u0 and v0 represent zonal wind and meridional wind respectively. If both gravitational effect and at- tenuation are neglected, Equation (20) would be sim- 2c k k k gk 0z0 2 2 2a b g cX+ + - =l l l l^ h x y 2 g k k c k k k g k k c k k k 1 4 2 0 z z z 2 2 2 4 0 2 2 2 2 2 2 3 2 2 c a b ab c a b a b X X X X X - + + - - - - - + + - + =- - l l l l l l l l l ^ ^ ^ ^ h h h h c c m m c k 2 0 2a g X =- l c g H2 2 1 0 2b c = = c k Gr At 04 2 2 2X X- + + =l x yGr c k kg a 2 2 2 2 2 2~ ~X= + -l l^ h 7 2At k k c c kH k 4 z2 2 2 2 2 2 2 23g gX X X X= - -+l l l l^ h` j 0c k2 2 2X - =l c k Gr 02 2 2X - + =l 0c k Ata 2 2 2 2~X - - + =l x yc k k 0g 2 2 2 2 2 "~ X+l l^ h x y x y x y k c k k At c k c k At u k c k k At c k c k At v k c k k At c k At 4 2 1 2 1 4 2 1 2 1 4 2 1 2 1 x g x g x k y g y g y k z g z k 3 2 2 2 2 2 2 2 2 0 3 2 2 2 2 2 2 2 2 0 3 2 2 2 2 2 x y z 2 2 2 2 2 2 2 2 2 2 2 2 ~ ~ ~ ~ ~ ~ ~ ~ X X X X X X X X =- - + + - + + + =- - + + - + + + =- - + + - + X X X l l l l l l l l l l l l l l ^ ^ ^ h h h Z [ \ ] ] ] ] ] ]] ] ] ] ] ] ]] SONG ET AL. 4 (13) (20) 5 plified as the group velocity derived from the classic acoustic dispersion relation. In Equation (20), we can find that the terms which mainly con- trolled by attenuation modify the valve of the group velocity. The gravity term Gr brings in the gravitational dispersion and makes the sound propagation become inhomogeneous. Taking the gravity and attenuation into consideration, the profile of group velocity in the atmosphere can be predicted more accurately. 2.4. Ray tracing equations Equation (16) is a real variable equation. On the basis of Equation (16), ray equations can be derived from the Hamilton equations, expressed as (21) where r is a point on the ray path, ∇r and ∇k are the gra- dient operators in r and k space, respectively, and H is the dispersion relation [Yeh and Liu 1974].By substitut- ing Equation (16) into (21), the ray equations can be de- rived as (22.a) (22.b) For the components of the wave number, (23) In Equations (22) and (23), (24.a) (24.b) (24.c) From Equations (22-24), we can observe that the parameter g and play important roles in the ray tracing equation. As mentioned in Section 2.3, with at- tenuation considered, the group velocity would be modified. What’s more, the gradient of kinematic vis- cosity also can enhance the gradient of the group velocity. Both terms demonstrate that the acoustic tra- jectory is influenced by the effect due to the attenua- tion in a dissipative atmosphere. For given atmospheric parameters including the wind profile, temperature field and the profile of vis- cosity g, the ray trajectory in the dissipative atmosphere can be determined by numerically solving Equations (22-24). 3. Parameter modifications in the realistic atmosphere Derivation of the dissipative dispersion is very com- plex, because a fifth-order determinant has to be solved. In order to simplify the derivation, two assumptions are employed. Firstly, the atmosphere is assumed isother- mal so that the gradients of pressure and temperature can be described in a simple form in terms of the scale height. Secondly, only the viscous dissipation in the mo- mentum equation is considered. In the realistic atmos- phere, the temperature varies with height, and other dissipation mechanisms, such as heat conduction and molecular relaxation losses, also play a role. Provided that the atmospheric temperature varies very slowly and the high order derivatives of scale height can be ignored, modifications of ~a and ~g can be obtained as (25) In Equation (25), ~an and ~B are the acoustic cut- off frequency and Vaisala-Brunt frequency in the non-isothermal atmosphere, respectively [Beer 1974, Matsumura et al. 2012]. By using the temperature and wind profiles given by the MSISE-00 [Picone et al. 2002] and HWM-93 [Hedin et al. 1996] empirical reference models, the typical values of acoustic cut-off frequency and Vaisala-Brunt frequency in the realistic atmosphere are shown in Figure 1. Varying with altitude, the acoustic cut-off frequency reaches the extreme value of about 0.0039 Hz and 0.005 Hz at the height of the strato- andAt Atk2 2Xl ,dt dr H dt d Hkk rd d= = z d d z x Ck H g k k Ak Bu 2 2z x 3 2 2 0 g c g X X = + - + + l l l l ^ h z d d z y Ck H g k k Ak Bv 2 2z y 3 2 2 0 g c g X X = + - + + l l l l ^ h z d dk k d d 0, dz z z Ck H g k k D E 2 2 x y x z 3 2 2g c gX X = = = = - + + + dk l l l^ h 2 4 2 A c c k g k 2 7 g 2 2 2 2 2 4 2 2 2 2 ~ g g gc X X X X = - - + + -l l ^ h 4 2 2 7 2 c k c k k k c k B H 6 a 2 2 2 2 2 2 2 2 2 2 3 3 4 ~ g g g X X X X X X = - - - + + + - l l l l l^ h 2 2 7 4C c c k 2 2 2 2 2 2 2 4g gX X X=- - + l yx 4 D z c k z z k c k k u k v c k2 a z a g x z y 2 2 2 2 0 2 2 2 2 2 0 2 2 2 2 2 2 2 2 2~ ~ ~ X X X = + - + + - + - + l l l l ^ ^ ^ ^ ^ ^ h h h h h h ; 6 E @ 7 7 2 2 E c z c k z k k c k z z c k k u k c k k H k H k c k z H k c k H k H c H v 8 7 6 2 4 x z z z z z y z 3 2 4 2 2 2 2 2 2 2 2 0 0 2 2 2 3 2 4 2 2 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 g g g g g g g g X X X X X X X X X X = + - + - + + + - + - - - - + - - l l l l l l l l l l l l ^ ^ h h ` ` c j j m z2 g z2 g ,T g z T T g z T an a B g 2 2 2 2 2 2 ~ ~ ~ ~= + = + ACOUSTIC RAY TRACING IN THE ATMOSPHERE (24.d) (24.e) sphere and thermosphere, respectively. The prediction of acoustic attenuation in the real- istic atmosphere should consider the contribution of different dissipation mechanisms including viscosity, diffusion and molecular relaxation [Brown and Hall 1978]. The coefficients for different kinds of dissipation mechanism can be combined as one net attenuation co- efficient in terms of linear superposition. In this study, kinematic viscosity g is modified by the value of the net attenuation coefficient, so that we can account for the realistic attenuation effects into dissipative disper- sion relation. According to the algorithms developed by Suther- land and Bass [2004, 2006], which can precisely estimate the acoustic absorption up to 160 km, the net attenua- tion coefficient at (Np/m), is written in the form of lin- ear summation (26) where acr is the attenuation coefficient associated with the combination of classical and rotational relaxation losses, adiff is the attenuation coefficient due to diffusion loss, and Amax,i is the maximum loss per wavelength for the ith molecular vibration relaxation component with a relaxation frequency fvib,i. In the realistic atmosphere, four primary gas components, O2, N2, CO2, and O3 are considered. Typical profiles of attenuation coefficient at frequencies of 0.05, 0.5, 5 and 8Hz are shown in Fig- ure 2. Clearly, the total attenuation coefficients for all considered frequencies increase significantly in the ther- mosphere. From Equation (14), in the windless conditions, the parameter g that describes the acoustic attenuation can be given by (27) Combining Equation (26) with Equation (27), we understand the parameter g as an artificial viscosity, which is determined as a function of frequency, temper- ature, pressure and relative humidity. Using the modified values of ~a, ~g, and g, the ray tracing model can compute the trajectory of acoustic waves more accurately, since the contributions of both temperature gradient and acoustic attenuation are taken into account. 4. Comparison with FDTD solution To validate this ray tracing model, we simulate the propagation of a Gaussian infrasound wave packet by FDTD method, and compare the trajectory obtained by FDTD solution with the result calculated by the ray tracing model. By using a finite-difference time-domain method [Wochner et al. 2005] consisting a dispersion relation preserving scheme in space and a Runge-Kutta scheme in time, this FDTD model is built on the basis of the governing equations of acoustic waves, Equation (4). In the FDTD experiment, the Gaussian infrasound wave packet at the frequency of 0.5 Hz is set to be emit- ted eastward with shooting angle of 15° at the height of 10 km. The pressure perturbation of the packet has the following form: (28) where pc = 100 Pa; xc = 8 km; zc = 10 km; kx and kz are the horizontal and vertical component of the wave num- ber vector respectively; and vx and vz are the half width of the wave packet in the horizontal and vertical direc- tion respectively. In the simulation, vx = vz = 1.6 km. According to the background profiles, the parameter kx and kz are initially set by using the dispersion Equation (16). The other initial quantities such as u1, w1 and t1 / / / A c f f f f12 ,maxt cr diff ii vib vib 22 # # a a a= + + +^ ^^h h h 6 6 @ @ " , / 2 .k c t 2 g ~ a = - ,p x z k x x k z zc x x z z x c z c 2 2x c z c 2 2 2 2 = - + -v v - - - sinp e e ^ ^ ^ ^ ^ h h h h h 6 @ SONG ET AL. 6 Figure 1. Profiles of acoustic cut-off frequency (solid) and Vaisala- Brunt frequency (dashed), calculated using the MSISE-00 and HWM-93 empirical model. Figure 2. Total attenuation coefficient as a function of altitude for frequencies of 0.05, 0.5, 5 and 8Hz (from left to right). 7 are derived from the polarization equations of acoustic waves in the atmosphere [Beer 1974]. The background profiles of density, temperature shown by Figure 3A are provided by the MSISE-00 [Pi- cone et al. 2002], and the wind profile shown by Figure 3B is provided by HWM-93 [Hedin et al. 1996]. The at- tenuation coefficient profile is derived by the algorithms founded by Sutherland and Bass [2004, 2006] which is already exhibited in Figure 2. For the MSISE and HWM models, the time we set is the 300th day of the year at 9:00 UT; while the geographic location of the source is (39.92°, 116.46°) and the parameter F10.7 index is 150. The simulation results are represented by the acoustic intensity of the packet. Continuous capture of the infrasound packet every 30s during the propagation are shown by Figure 3C. As seen in Figure 3C, the po- sition and altitude of the packet are changed due to the refraction controlled by the atmospheric profile. Because of the effects of attenuation and geometric spreading, the intensity of the infrasound packet becomes weaker as the packet gets higher. In this case, the trajectory of the infrasound is the wave energy propagation path which can be represented by the positions of the wave packet’s energy center. The propagation path of the packet ob- tained from the FDTD solution and the ray path pre- dicted by the ray tracing model are compared in Figure 4. The comparison shows that the ray path predicted by the ray tracing model is equal to the propagation path of the infrasound packet obtained from the FDTD solution within the computational precision. Based on the results shown in Figure 4, we can conclude that this ray tracing can be proven accurate and reliable. 5. Numerical results To investigate the impacts of gravity dispersion and acoustic dissipation, based on the background models of the MSISE-00 and HWM-93, three cases of ACOUSTIC RAY TRACING IN THE ATMOSPHERE Figure 4. Comparison between the trajectory obtained by the FDTD solution and ray solution. The FDTD solution is denoted by ‘+’, while the ray solution is represented by the solid line. Figure 3. (A) background profiles of density and temperature; (B) background profiles of zonal wind and meridional wind; (C) the acoustic intensity of the infrasound packet in the form of continu- ous capture every 30s. A B C numerical results are presented. The first is calculated using the acoustic ray tracing model for the realistic at- mosphere, as described in Section 2, in which both gravity dispersion and acoustic dissipation are consid- ered. The second is calculated using the lossless ray tracing model that considers only the gravity disper- sion, which can also be called “the acoustic gravity wave ray tracing model” [Yeh and Liu 1974]. The last is obtained from the classical ray tracing model in which either gravity dispersion or acoustic dissipation is neg- lected [Jones et al. 1986]. The numerical results for the above three cases are compared with each other under the condition of long range sound propagation (LRSP). 5.1. Trajectories In this calculation, the input parameters of MSISE- 00 and HWM-93 are set as same as the validation ex- periment in Section 4. Under the background profiles of density, temperature and wind shown by Figure 3A and Figure 3B, We assume that the acoustic wave with a frequency of 5 Hz is emitted eastward by an infra- sonic source on the ground. Shooting angles of the rays vary from 3° to 43°, with an increment of 8°. Figure 5 shows the wave trajectories for the three cases. The solid, dotted, and dashed lines denote the wave trajectories obtained using the acoustic ray trac- ing model for the realistic atmosphere, the acoustic gravity wave ray tracing model, and the classical ray tracing model, respectively. Two channels can be iden- tified from these modeled trajectories. One is the re- flection from the height less than 50 km, which is called “stratospheric ducting”, and the other is “thermos- pheric ducting” that occurs at higher altitudes from the height of thermosphere. The effect of acoustic dissipa- tion can be registered by comparing the solid and dot- ted curves, and the effect due to the gravity can be eval- uated by comparing the dotted and dashed curves. We first discuss the effect of acoustic dissipation and then the gravitational effect. For the rays with shooting angles of 3°-19°, which are in the channel of stratospheric ducting, the two sets of trajectories coincide well. But for the other rays with shooting angles of 27°-43° that are in the channel of thermospheric ducting, there are obvious differences between the two sets of trajectories. It indicates that the dissipation of acoustic waves has an effect on the trajectory of 5 Hz acoustic waves. The influence due to the dissipation is quite significant for thermospheric ducting but less important to stratospheric ducting. It is mainly because the attenuation coefficient varies with altitude. As shown in Figure 2, the net attenuation co- efficient for 5 Hz waves reaches a peak valve (~0.03 dB/km) at the altitude of 10km in the troposphere and stratosphere. As expressed by Equation (18), the dis- persion relation would be reduced to the traditional dis- persion of acoustic gravity wave when the attenuation could be neglected. As the altitude increases, the atten- uation coefficient increases, reaching ~400 dB/km in the thermosphere. Because the attenuation coefficient is very small (< 0.03 dB/km), the effect of dissipation has little influence on the rays in the stratospheric duct- ing. For thermospheric ducting, in the height range 100 km < z <140 km, the turning points of ray trajec- tories with the realistic attenuation are at an altitude lower than that of rays in the lossless atmosphere. Ob- viously, the refraction is enhanced by the strong acoustic attenuation in the thermosphere. By comparing the dotted and dashed lines, dis- crepancies occur in both stratospheric ducting and ther- mospheric ducting. For the trajectory with shooting angle 3°, due to the gravitational effect, the dotted lines are below the dashed line in the region z <160 km. After the intersection of the two types of lines, the dot- ted lines are above the dashed lines. However, this phe- nomenon does not occur for all rays. Take the trajectory with shooting angle 19° for an example. The dotted lines are always above the dashed lines. Since compar- isons between the trajectories are complicated, we will describe the discrepancies of the three kinds of trajec- tories via some concrete parameters such as turning point and propagation distance. 5.2. Turning point and propagation distance Turning point and propagation distance of the three kinds of trajectories as a function of elevation are shown in Figure 6. With the same range of shooting angle but a smaller step of 1°, a clearer view of the in- SONG ET AL. 8 Figure 5. Ray paths estimated based on three ray tracing models. The solid and dotted lines show the results calculated using our acoustic ray tracing model and the acoustic gravity wave ray tracing model respectively, which both consider the effect of gravity dispersion. The dashed lines represent the results calculated using the classical ray racing model that ignores both gravity dispersion and acoustic dis- sipation. Shooting angles of 5 Hz acoustic waves vary from 3 to 43° with a step of 8°. 9 fluences associated with the gravity dispersion and dis- sipation can be obtained. A turning point is an altitude that characterizes the transformation of wave rays between the two channels, i.e., stratospheric ducting and thermospheric ducting. In Figure 6A, for the solid and dotted lines, the eleva- tion range of stratosphere ducting is from 3° to 19°, while the elevation range of stratosphere ducting is 20°~43°. For the dashed lines, the elevation range of stratospheric ducting and thermospheric ducting are 3°~21° and 22°~43°, respectively. Due to the gravita- tional effect, the rays with shoot angles 20° and 21° no longer belong to the stratospheric ducting but the ther- mospheric ducting, suggesting that the gravitational ef- fect play an important role in determining the channel of acoustic ray trace. By comparing the dotted line and the dashed line in Figure 6A, the calculation deviation of the turning point is about 1.5 km in stratospheric ducting, and it will rise to about 8 km (at the initial elevation 31°) in the propa- gation channel of thermospheric ducting. As indicated in Figure 5, little difference is observed between solid line and dotted lines in stratospheric ducting in both Fig- ure 6A and Figure 6B. This difference illustrates that the deviation between the solid line and dashed line in stratospheric ducting is caused by the gravity. In thermospheric ducting, deviation of the turning point between solid line and dotted lines could reach about 15 km, as indicate in Figure 6A. Controlled by the deviation of the turning point, deviation of the propa- gation distance between solid line and dotted lines in- creases along with the initial elevation, and reaches the peak value of about 160 km, then decreases to the level of approximately 60 km (Figure 6B). These great devia- tions between solid line and the other lines shown in Figure 5 indicate that the strong attenuation is a domi- nant factor for the ducts in the thermosphere. As seen in Figure 2, the profile of acoustic attenuation varies with frequency. The attenuation decreases when frequency goes down, so for the propagation with the frequency smaller than 5 Hz, the influence due to attenuation in the thermospheric ducting will become weaker. 6. Conclusions An acoustic ray tracing model including the effects of gravity and attenuation is developed in the present study. The propagation of acoustic waves in non- isothermal lossy atmosphere can be described by this model. By comparing the ray path estimated by this model with the FDTD solution, the ray tracing model has been proven accurate and reliable. The influences due to absorption and gravity are investigated by com- paring the numerical results obtained using this acoustic ray tracing model with the other two models, i.e., acoustic gravity wave ray tracing model and the classic ray tracing model that ignores the two effects. The results indicate that in stratospheric ducting, the gravitational effect plays a leading role while the ef- fect associated with the dissipation is weak. Because of the gravity, some rays ducted in the stratosphere as a lossless medium can merge into the channel of ther- mospheric ducting in the realistic atmosphere. For ther- mospheric ducting, the refraction is enhanced by the strong absorption. This results in the dominance of dis- sipative effect over gravitational effect in thermospheric ducting. To conclude, by taking the gravity into considera- tion, the acoustic frequency in our model is no longer constrained by the high frequency approximation. In- stead, it is available to accurately predict the propaga- tion of infrasound with frequencies lower than 5 Hz. Consideration of the dissipation effect in our model can also improve the computation accuracy of ray trajec- tories and the efficiency of wave amplitude evaluation. ACOUSTIC RAY TRACING IN THE ATMOSPHERE Figure 6. Comparison of the turning height (A) and propagation distance (B) as a function of elevation angle. The solid and dotted lines show the results calculated using our acoustic ray tracing model and the acoustic gravity wave ray tracing model respectively, which both consider the effect of gravity dispersion. The dashed lines represent the results calculated using the classical ray racing model that ignores both grav- ity dispersion and acoustic dissipation. Shooting angles vary from 3~43° with a step of 1°. A B References Barry, G. (1963). Ray tracings of acoustic waves in the upper atmosphere, J. Atmos. Terr. Phys., 25 (11), 621-629. Bass, H.E., C.H. Hetzer and R. Raspet (2007). On the speed of sound in the atmosphere as a function of altitude and frequency, J. Geophys. Res., 112, D15110; doi:10.1029/2006JD007806. Bedard, A.J., and T.M. Georges (2000). Atmospheric In- frasound, Phys. Today, 53 (3), 32-37. Beer, T. (1974). Atmospheric waves, Adam Hilger, London. Brown, E.H., and F.F. Hall (1978). Advances in Atmos- pheric Acoustics, Rev. Geophys., 16 (1), 47-110. Buland, R., and C. Chapman (1983). The computation of seismic travel times, J. Acoust Soc. Am., 73, 1271- 1302. de Groot-Hedlin, C., M.A.H. Hedlin and K. Walker (2011). Finite difference synthesis of infrasound propagation through a windy, viscous atmosphere, Geophys. J. Int., 185 (1), 305-320. Dessa, J.X., J. Virieux and S. Lambotte (2005). Infra- sound modeling in a spherical heterogeneous at- mosphere, Geophys. Res. Lett., 32, L12808; doi:10. 1029/2005GL022867. Drob, D.P., J.M. Picone and M. Garcés (2003). Global morphology of infrasound propagation, J. Geophys. Res., 108 (D21), 4680; doi:10.1029/2002JD003307. Evers, L.G., and H.W. Haak (2005). The detectability of infrasound in The Netherlands from the Italian vol- cano Mt. Etna, J. Atmos. Terr. Phys., 67 (3), 259-268. Garcés, M.A., R.A. Hansen and K.G. Lindquist (1998). Traveltimes for infrasonic waves propagating in stratified atmosphere, Geophys. J. Int., 135, 255-263. Hedin, A.E., E.L. Fleming, A.H. Manson, F.J. Schmidlin, S.K. Avery, R.R. Clark, S.J. Franke, G.J. Fraser, T. Tsuda, F. Vial and R.A. Vincent (1996). Empirical wind model for the upper, middle and lower at- mosphere, J. Atmos. Terr. Phys., 58, 1421-1447. Hedlin, M.A.H., D. Drob, K. Walker and C. de Groot- Hedlin (2010). A study of acoustic propagation from a large bolide in the atmosphere with a dense seis- mic network, J. Geophys. Res., 115, B11312; doi:10. 1029/2010JB007669. Jones, R.M. (1970). Ray theory for lossy media, Radio Sci., 5, 793-801. Jones, R.M., J.P. Riley and T.M. Georges (1986). HARPA - a versatile three-dimensional Hamiltonian ray-trac- ing program for acoustic waves in the atmosphere above irregular terrain, NOAA Special Report; http: //cires.colorado.edu/~mjones/raytracing. Kanamori, H., J. Mori and D.G. Harkrider (1994). Exci- tation of atmospheric oscillations by volcanic erup- tions, J. Geophys. Res., 99 (B11), 21947-21961. Kulichkov, S.N. (2004). Long-range propagation and scattering of low-frequency sound pulses in the mid- dle atmosphere, Meteorol. Atmos. Phys., 85, 47-60. Le Pichon, A., E. Blanc and D. Drob (2005). Probing high-altitude winds using infrasound, J. Geophys. Res., 110, D20104; doi:10.1029/2005JD006020. Le Pichon, A., J. Vergoz, E. Blanc, J. Guilbert, L. Cer- anna, L. Evers and N. Brachet (2009). Assessing the performance of the International Monitoring Sys- tem’s infrasound network: Geographical coverage and temporal variabilities, J. Geophys. Res., 114, D08112; doi:10.1029/2008JD010907. Le Pichon, A., E. Blanc and A. Hauchecorne (2010). Infrasound Monitoring for Atmospheric Studies, Springer, New York. Le Pichon, A., L. Ceranna, C. Pilger, P. Mialle, D. Brown, P. Herry and N. Brachet (2013). The 2013 Russian fireball largest ever detected by CTBTO infrasound sensors, Geophys. Res. Lett., 40, 3732-3737; doi:10. 1002/grl.50619. Matsumura, M., H. Shinagawa and T. Iyemori (2012). Horizontal extension of acoustic resonance between the ground and the lower atmosphere, J. Atmos. Terr. Phys., 75-76, 127-132. McKenna, M.H., B.W. Stump, S. Hayek, J.R. McKenna and T.R. Stanton (2007). Tele-infrasonic studies of hard-rock mining explosions, J. Acoust. Soc. Am., 122 (1), 97-106. Picone, J.M., A.E. Hedin, D.P. Drob and A.C. Aikin (2002). NRLMSISE-00 - Empirical model of the at- mosphere: statistical comparisons and scientific is- sues, J. Geophys. Res., 107, 1468; doi:10.1029/2002 JA009430. Pilger, C., and M. Bittner (2009). Infrasound from tro- pospheric sources: Impact on mesopause tempera- ture?, J. Atmos. Terr. Phys., 71, 816-822. Sonnenschein, E., D. Censor, I. Rutkevich and J.A. Ben- nett (1997). Ray trajectories in an absorbing iono- sphere, J. Atmos. Terr. Phys., 59 (16), 2101-2110. Stevens, J.L., I.I. Divnov, D.A. Adams, J.R. Murphy and V.N. Bourchik (2002). Constraints on infrasound scaling and attenuation relations from Soviet explo- sion data, Pure. Appl. Geophys., 159, 1045-1062. Sutherland, L.C., and H.E. Bass (2004). Atmospheric absorption in the atmosphere up to 160 km, J. Acoust. Soc. Am., 115 (3), 1012-1032. Sutherland, L.C., and H.E. Bass (2006). Erratum: “At- mospheric absorption in the atmosphere up to 160 km”, J. Acoust. Soc. Am., 120 (5), 2985. Virieux, J., N. Garnier, E. Blanc and J.X. Dessa (2004). Paraxial ray tracing for atmospheric wave propaga- tion, Geophys. Res. Lett., 31, L20106; doi:10.1029/ SONG ET AL. 10 11 2004GL020514. Wochner, M.S., A.A. Atchley and V.W. Sparrow (2005). Numerical simulation of finite amplitude wave prop- agation in air using a realistic atmospheric absorp- tion model, J. Acoust. Soc. Am., 118 (5), 2891-2898. Yeh, K.C., and C.H. Liu (1974). Acoustic-gravity waves in the upper atmosphere, Rev. Geophys., 12 (2), 193- 216. *Corresponding author: Yang Song, Wuhan University, Department of Space Physics, School of Electronic Information, Wuhan, China; email: phsongyang@gmail.com. © 2014 by the Istituto Nazionale di Geofisica e Vulcanologia. All rights reserved. ACOUSTIC RAY TRACING IN THE ATMOSPHERE << /ASCII85EncodePages false /AllowTransparency false /AutoPositionEPSFiles false /AutoRotatePages /None /Binding /Left /CalGrayProfile (Dot Gain 20%) /CalRGBProfile (sRGB IEC61966-2.1) /CalCMYKProfile (U.S. Web Coated \050SWOP\051 v2) /sRGBProfile (sRGB IEC61966-2.1) /CannotEmbedFontPolicy /Warning /CompatibilityLevel 1.3 /CompressObjects /Tags /CompressPages true /ConvertImagesToIndexed true /PassThroughJPEGImages true /CreateJobTicket false /DefaultRenderingIntent /Default /DetectBlends true /DetectCurves 0.1000 /ColorConversionStrategy /LeaveColorUnchanged /DoThumbnails false /EmbedAllFonts true /EmbedOpenType false /ParseICCProfilesInComments true /EmbedJobOptions true /DSCReportingLevel 0 /EmitDSCWarnings false /EndPage -1 /ImageMemory 1048576 /LockDistillerParams true /MaxSubsetPct 100 /Optimize false /OPM 1 /ParseDSCComments true /ParseDSCCommentsForDocInfo true /PreserveCopyPage true /PreserveDICMYKValues true /PreserveEPSInfo true /PreserveFlatness true /PreserveHalftoneInfo false /PreserveOPIComments false /PreserveOverprintSettings true /StartPage 1 /SubsetFonts true /TransferFunctionInfo /Apply /UCRandBGInfo /Preserve /UsePrologue false /ColorSettingsFile (None) /AlwaysEmbed [ true /AndaleMono /Apple-Chancery /Arial-Black /Arial-BoldItalicMT /Arial-BoldMT /Arial-ItalicMT /ArialMT /CapitalsRegular /Charcoal /Chicago /ComicSansMS /ComicSansMS-Bold /Courier /Courier-Bold /CourierNewPS-BoldItalicMT /CourierNewPS-BoldMT /CourierNewPS-ItalicMT /CourierNewPSMT /GadgetRegular /Geneva /Georgia /Georgia-Bold /Georgia-BoldItalic /Georgia-Italic /Helvetica /Helvetica-Bold /HelveticaInserat-Roman /HoeflerText-Black /HoeflerText-BlackItalic /HoeflerText-Italic /HoeflerText-Ornaments /HoeflerText-Regular /Impact /Monaco /NewYork /Palatino-Bold /Palatino-BoldItalic /Palatino-Italic /Palatino-Roman /SandRegular /Skia-Regular /Symbol /TechnoRegular /TextileRegular /Times-Bold /Times-BoldItalic /Times-Italic /Times-Roman /TimesNewRomanPS-BoldItalicMT /TimesNewRomanPS-BoldMT /TimesNewRomanPS-ItalicMT /TimesNewRomanPSMT /Trebuchet-BoldItalic /TrebuchetMS /TrebuchetMS-Bold /TrebuchetMS-Italic /Verdana /Verdana-Bold /Verdana-BoldItalic /Verdana-Italic /Webdings ] /NeverEmbed [ true ] /AntiAliasColorImages false /CropColorImages true /ColorImageMinResolution 150 /ColorImageMinResolutionPolicy /OK /DownsampleColorImages true /ColorImageDownsampleType /Bicubic /ColorImageResolution 300 /ColorImageDepth -1 /ColorImageMinDownsampleDepth 1 /ColorImageDownsampleThreshold 1.10000 /EncodeColorImages true /ColorImageFilter /DCTEncode /AutoFilterColorImages true /ColorImageAutoFilterStrategy /JPEG /ColorACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /ColorImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000ColorACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000ColorImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasGrayImages false /CropGrayImages true /GrayImageMinResolution 150 /GrayImageMinResolutionPolicy /OK /DownsampleGrayImages true /GrayImageDownsampleType /Bicubic /GrayImageResolution 300 /GrayImageDepth -1 /GrayImageMinDownsampleDepth 2 /GrayImageDownsampleThreshold 1.10000 /EncodeGrayImages true /GrayImageFilter /DCTEncode /AutoFilterGrayImages true /GrayImageAutoFilterStrategy /JPEG /GrayACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /GrayImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000GrayACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000GrayImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasMonoImages false /CropMonoImages true /MonoImageMinResolution 1200 /MonoImageMinResolutionPolicy /OK /DownsampleMonoImages true /MonoImageDownsampleType /Bicubic /MonoImageResolution 1200 /MonoImageDepth -1 /MonoImageDownsampleThreshold 1.08250 /EncodeMonoImages true /MonoImageFilter /CCITTFaxEncode /MonoImageDict << /K -1 >> /AllowPSXObjects false /CheckCompliance [ /None ] /PDFX1aCheck false /PDFX3Check false /PDFXCompliantPDFOnly false /PDFXNoTrimBoxError true /PDFXTrimBoxToMediaBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXSetBleedBoxToMediaBox true /PDFXBleedBoxToTrimBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXOutputIntentProfile (None) /PDFXOutputConditionIdentifier () /PDFXOutputCondition () /PDFXRegistryName (http://www.color.org) /PDFXTrapped /Unknown /CreateJDFFile false /SyntheticBoldness 1.000000 /Description << /ENU (Use these settings to create PDF documents with higher image resolution for high quality pre-press printing. The PDF documents can be opened with Acrobat and Reader 5.0 and later. These settings require font embedding.) /JPN /FRA /DEU /PTB /DAN /NLD /ESP /SUO /NOR /SVE /KOR /CHS /CHT /ITA >> >> setdistillerparams << /HWResolution [2400 2400] /PageSize [595.000 842.000] >> setpagedevice