Acoustic ray tracing in the atmosphere : with gravitational effect and attenuation considered

An acoustic ray tracing model is developed to take into account the impacts of gravitational field and realistic atmospheric attenuation. Ray tracing equations are deduced from the real part of the dissipative dispersion 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 frequency 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 predicted 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.


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 temperature 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 observed at infrasonic frequencies because of the low attenuation in the lower atmosphere [Bedard and Georges 2000].Consequently, natural and artificial sources could be detected and monitored through the network of infrasonic stations, such as the International Monitoring System [Le Pichon et al. 2009].In addition, as the propagation 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 propagation properties of acoustic waves.
Ray tracing is an efficient way to study the propagation 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 available, 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 applied to study the sound propagation [e.g., Le Pichon et al. 2010].HARPA is derived from the dispersion relation 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 [Garcés et al. 1998] is developed for computation of travel times as a reformulation of the original method presented by Buland and Chapman [1983].All of these three models are developed to study the atmospheric wave propagation through a windy, non-isothermal, non-dissipative atmosphere.
Theoretically, the high frequency approximation adopted by the above three models confines the ray tracing modeling to acoustic signals above 5 Hz for which the gravity influence is neglected.Moreover, the attenuation 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 were observed in the cases of large bolides [Hedlin et al. 2010], meteorites [Le Pichon et al. 2013], volcano [Kanamori et al. 1994, Evers andHaak 2005] and explosions [Stevens et al. 2002, McKenna et al. 2007].Under the influence of gravity, the propagation of acoustic waves becomes more dispersive and anisotropic as frequency decreases.For modeling the propagation of acoustic waves below 5 Hz, the high frequency approximation 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 developed by Sutherland and Bass [2004], the attenuation effects 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 instead 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 expression 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 [Sonnenschein et al. 1997].When studying the acoustic wave propagation, especially the infrasonic propagation, both the effects of attenuation and the gravitational dispersion 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 extension.
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 dissipative dispersion relation.In Section 2.1, the dissipative 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 substituted 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 andBass 2004, 2006], aiming to include the non-isothermal effect and the realistic atmospheric attenuation.This ray tracing model is validated by the comparison with FDTD solution in Section 4. Numerical examples are shown in Section 5 to study the influences of the gravitational effect and attenuation, followed by the conclusion in Section 6.

Inclusion of gravity and absorption into ray tracing equations
In order to introduce gravity influence and dissipation 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, entropy, 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 derivation, the viscous term is neglected in following calculations.
With the knowledge of the thermodynamic relationship, the entropy and temperature in the atmosphere can be given by (2) In Equations (2), c = c p /c v is the ratio between specific 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 irrotational wave, the variable With the approximation of a stratified atmosphere, Equation (1) can be simplified and linearized to yield where V 0 = (u 0 ,v 0 ,0) is the wind profile of the atmosphere, is the adiabatic sound speed with T 0 as the ambient temperature.

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=(k x ,k y ,k z ) are angular frequency and wave number components, respectively.Substitution into Equation (4) yields a set of fifth-order algebraic equations for (u 1 ,v 1 ,w 1 ,t 1 /t 0 ,p 1 /p 0 ), which can be expressed in the matrix form as (6) In the coefficient matrix of Equation ( 6), X = ~− V 0 • k, the Doppler-shifted frequency, is the coefficient of kinematic viscosity of air, H = c 2 /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 inhomogeneity 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.Subsequently, the dispersion relation Equation ( 7) can be decomposed 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 dispersion relation, g 2 (c − 1)(k 2 x + k 2 y ), plays an important role in the real part of the dispersion relation, as it represents the influence of the gravity that generates gravity dispersion.This term is discarded in the derivation of Equation ( 9), because the imaginary part of the dissipative 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 variations along a ray path.By substituting k" into Equation ( 10), the ray equations can be determined on the basis of Hamilton equations.

The attenuation and growth of acoustic waves
As acoustic waves propagate in the atmosphere, the wave amplitude usually drops due to the attenuation effect.Growth in the wave amplitude may occur because of the atmospheric density gradient in the vertical direction.Thus, in the dissipative stratified atmosphere, the imaginary wave vector can be considered as a combination of omnidirectional attenuation and growth in the vertical direction.Assuming the attenuation efficient and growth factor as a and b respectively, we obtain that (11) By substituting Equation (11) into Equation ( 9) and ( 10 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 atmosphere, 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.

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)g 2 /c 2 , ~g is the isothermal Vaisala-Brunt frequency; and ~a = c/(2H) is the acoustic cut-off frequency; the term is the contribution of attenuation.If both gravitational effect and attenuation effect are neglected, the dispersion will become 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 frequency 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 obtain 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 relation.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 relation, 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 ~2 a can be taken as the contribution of the inhomogeneity of atmosphere density.
Both limiting cases demonstrate that the inclusion of the gravity and attenuation are reasonable enrichment for the classic acoustic dispersion relation.The effect caused by gravity and attenuation can be studied by the group velocity which can be derived as where u 0 and v 0 represent zonal wind and meridional wind respectively.If both gravitational effect and attenuation are neglected, Equation ( 20) would be sim- (20) plified as the group velocity derived from the classic acoustic dispersion relation.In Equation ( 20), we can find that the terms which mainly controlled 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.

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 where r is a point on the ray path, ∇ r and ∇ k are the gradient operators in r and k space, respectively, and H is the dispersion relation [Yeh and Liu 1974].By substituting Equation ( 16) into ( 21), the ray equations can be derived as In Equations ( 22) and ( 23 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 attenuation considered, the group velocity would be modified.What's more, the gradient of kinematic viscosity also can enhance the gradient of the group velocity.Both terms demonstrate that the acoustic trajectory is influenced by effect attenuation in a dissipative atmosphere.
For given atmospheric parameters including the wind profile, temperature field and the profile of viscosity g, the ray trajectory in the dissipative atmosphere can be determined by numerically solving Equations (22-24).

Parameter modifications in the realistic atmosphere
Derivation of the dissipative dispersion is very complex, because a fifth-order determinant has to be solved.In order to simplify the derivation, two assumptions are employed.Firstly, the atmosphere is assumed isothermal 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 momentum equation is considered.In the realistic atmosphere, 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 cutoff 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- sphere and thermosphere, respectively.The prediction of acoustic attenuation in the realistic 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 coefficient 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 dispersion relation.
According to the algorithms developed by Sutherland and Bass [2004Bass [ , 2006]], which can precisely estimate the acoustic absorption up to 160 km, the net attenuation coefficient a t (Np/m), is written in the form of linear summation (26) where a cr is the attenuation coefficient associated with the combination of classical and rotational relaxation losses, a diff is the attenuation coefficient due to diffusion loss, and A max,i is the maximum loss per wavelength for the ith molecular vibration relaxation component with a relaxation frequency f vib,i .In the realistic atmosphere, four primary gas components, O 2 , N 2 , CO 2 , and O 3 are considered.Typical profiles of attenuation coefficient at frequencies of 0.05, 0.5, 5 and 8Hz are shown in Figure 2. Clearly, the total attenuation coefficients for all considered frequencies increase significantly in the thermosphere.
From Equation ( 14), in the windless conditions, the parameter g that describes the acoustic attenuation can be given by ( 26) with Equation ( 27), we understand the parameter g as an artificial viscosity, which is determined as a function of frequency, temperature, 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.

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 emitted eastward with shooting angle of 15° at the height of 10 km.The pressure perturbation of the packet has the following form: (28) where p c = 100 Pa; x c = 8 km; z c = 10 km; k x and k z are the horizontal and vertical component of the wave number vector respectively; and v x and v z are the half width of the wave packet in the horizontal and vertical direction respectively.In the simulation, v x = v z = 1.6 km.According to the background profiles, the parameter k x and k z are initially set by using the dispersion Equation ( 16).The other initial quantities such as u 1 , w 1 and  are derived from the polarization equations of acoustic waves in the atmosphere [Beer 1974].

sin p e e
The background profiles of density, temperature shown by Figure 3A are provided by the MSISE-00 [Picone et al. 2002], and the wind profile shown by Figure 3B is provided by HWM-93 [Hedin et al. 1996].The attenuation coefficient profile is derived by the algorithms founded by Sutherland andBass [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 position 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 obtained from the FDTD solution and the ray path predicted 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.

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  numerical results are presented.The first is calculated using the acoustic ray tracing model for the realistic atmosphere, as described in Section 2, in which both gravity dispersion and acoustic dissipation are considered.The second is calculated using the lossless ray tracing model that considers only the gravity dispersion, 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 neglected [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).

Trajectories
In this calculation, the input parameters of MSISE-00 and HWM-93 are set as same as the validation experiment 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 infrasonic 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 tracing model for the realistic atmosphere, the acoustic gravity wave ray tracing model, and the classical ray tracing model, respectively.Two channels can be identified from these modeled trajectories.One is the reflection from the height less than 50 km, which is called "stratospheric ducting", and the other is "thermospheric 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 dotted curves, and the effect due to the gravity can be evaluated 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 coefficient 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 dispersion relation would be reduced to the traditional dispersion of acoustic gravity wave when the attenuation could be neglected.As the altitude increases, the attenuation 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 ducting.For thermospheric ducting, in the height range 100 km < z <140 km, the turning points of ray trajectories with the realistic attenuation are at an altitude lower than that of rays in the lossless atmosphere.Obviously, the refraction is enhanced by the strong acoustic attenuation in the thermosphere.
By comparing the dotted and dashed lines, discrepancies occur in both stratospheric ducting and thermospheric 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 dotted lines are above the dashed lines.However, this phenomenon 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 comparisons between the trajectories are complicated, we will describe the discrepancies of the three kinds of trajectories via some concrete parameters such as turning point and propagation distance.

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- fluences associated with the gravity dispersion and dissipation 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 elevation 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 gravitational effect, the rays with shoot angles 20° and 21° no longer belong to the stratospheric ducting but the thermospheric ducting, suggesting that the gravitational effect 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 propagation channel of thermospheric ducting.As indicated in Figure 5, little difference is observed between solid line and dotted lines in stratospheric ducting in both Figure 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 propagation distance between solid line and dotted lines increases 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 deviations between solid line and the other lines shown in Figure 5 indicate that the strong attenuation is a dominant 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.

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 nonisothermal 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 comparing 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 effect 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 thermospheric ducting in the realistic atmosphere.For thermospheric ducting, the refraction is enhanced by the strong absorption.This results in the dominance of dissipative effect over gravitational effect in thermospheric ducting.
To conclude, by taking the gravity into consideration, the acoustic frequency in our model is no longer constrained by the high frequency approximation.Instead, it is available to accurately predict the propagation of infrasound with frequencies lower than 5 Hz.Consideration of the dissipation effect in our model can also improve the computation accuracy of ray trajectories and the efficiency of wave amplitude evaluation.Numerical simulation of finite amplitude wave propagation in air using a realistic atmospheric absorption 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.
), the imaginary part and real part of the dissipative dispersion relation can be rewritten

Figure 2 .
Figure 2. Total attenuation coefficient as a function of altitude for frequencies of 0.05, 0.5, 5 and 8Hz (from left to right).

Figure 4 .
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 .
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 continuous capture every 30s.

Figure 5 .
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 dissipation.Shooting angles of 5 Hz acoustic waves vary from 3 to 43°w ith a step of 8°.

Figure 6 .
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 gravity dispersion and acoustic dissipation.Shooting angles vary from 3~43° with a step of 1°.A B