Verifications of the high-resolution numerical model and polarization relations of atmospheric acoustic-gravity waves

Comparisons of amplitudes of wave variations of atmospheric characteristics obtained using direct numerical simulation models with polarization relations given by conventional theories of linear acoustic-gravity waves (AGWs) could be helpful for testing these numerical models. In this study, we performed high-resolution numerical simulations of nonlinear AGW propagation at altitudes 0–500 km from a plane wave forcing at the Earth’s surface and compared them with analytical polarization relations of linear AGW theory. After some transition time te (increasing with altitude) subsequent to triggering the wave source, the initial wave pulse disappears and the main spectral components of the wave source dominate. The numbers of numerically simulated and analytical pairs of AGW parameters, which are equal with confidence of 95 %, are largest at altitudes 30– 60 km at t > te. At low and high altitudes and at t < te, numbers of equal pairs are smaller, because of the influence of the lower boundary conditions, strong dissipation and AGW transience making substantial inclinations from conditions, assumed in conventional theories of linear nondissipative stationary AGWs in the free atmosphere. Reasonable agreements between simulated and analytical wave parameters satisfying the scope of the limitations of the AGW theory prove the adequacy of the used wave numerical model. Significant differences between numerical and analytical AGW parameters reveal circumstances when analytical theories give substantial errors and numerical simulations of wave fields are required. In addition, direct numerical AGW simulations may be useful tools for testing simplified parameterizations of wave effects in the atmosphere.


Introduction
Observations show the frequent presence of acoustic-gravity waves (AGWs) generating at tropospheric heights and propagating to the middle and upper atmosphere (e.g., Fritts and Alexander, 2003).These AGWs can break and produce turbulence and perturbations in the atmosphere (Gavrilov and Yudin, 1992;Gavrilov and Fukao, 1999).For example, sources of AGWs could be mesoscale turbulence and convection in the troposphere (e.g., Fritts and Alexander, 2003;Fritts et al., 2006).Turbulent AGW generation may have maxima at altitudes 9-12 km in the regions of tropospheric jet streams (Medvedev and Gavrilov, 1995).
Non-hydrostatic models are useful for direct numerical simulations of waves and turbulence in the atmosphere.For example, Baker and Schubert (2000) simulated nonlinear AGWs in the atmosphere of Venus.They modeled waves in the atmospheric region with horizontal and vertical dimensions of 120 and 48 km, respectively.Fritts and Garten (1996), also Andreassen et al. (1998) and Fritts et al. (2009Fritts et al. ( , 2011)), simulated the instabilities of Kelvin and Helmholtz and turbulence produced by breaking atmospheric waves.These models simulate turbulence and waves in atmospheric regions with limited vertical and horizontal dimensions.The models exploited spectral methods and Galerkintype series for converting partial differential equations (versus time) into the ordinary differential equations for the spectral series components.Yu and Hickey (2007) and Liu et al. (2008) developed two-dimensional numerical models of atmospheric AGWs.Gavrilov and Kshevetskii (2013) described a twodimensional model for high-resolution numerical simulat-N.M. Gavrilov et al.: Verifications of the high-resolution AGW model ing nonlinear AGWs using a finite-difference scheme taking into account hydrodynamic conservation laws as described by Kshevetskii and Gavrilov (2005).This approach increases the stability of the numerical scheme and allows us to obtain non-smooth solutions of nonlinear wave equations.This permitted us to get generalized physically acceptable solutions to the equations (Lax, 1957;Richtmayer and Morton, 1967).Gavrilov and Kshevetskii (2014a) created a threedimensional version of this algorithm for simulating nonlinear AGWs in the atmosphere.They modeled waves produced by sinusoidal horizontally homogeneous wave forcing at the Earth's surface.Karpov and Kshevetskii (2014) used a similar numerical three-dimensional model to study AGW propagation from local non-stationary wave excitation at the Earth's surface.They showed that infrasound going from tropospheric sources could provide substantial mean heating in the upper atmosphere.Dissipating nonlinear AGWs can also create accelerations of the mean flows in the middle atmosphere (e.g., Fritts and Alexander, 2003).However, details of the mean heating and mean flows created by non-stationary nonlinear AGWs in the atmosphere need further studies.
Numerical models of atmospheric AGWs require verifications.For plane stationary wave components with small amplitudes, conventional linear theories (e.g., Gossard and Hooke, 1975) give the dispersion equation and polarization relations, which connect wave frequency, vertical and horizontal wave numbers and ratios of amplitudes of different wave field variations.One can expect that such relations could exist between corresponding parameters of the numerical model solutions.Therefore, theoretical polarization relations could be useful for verifying direct simulation models of atmospheric AGWs.
In this paper, using the high-resolution numerical threedimensional model by Gavrilov and Kshevetskii (2014a, b), we made comparisons of calculated ratios of amplitudes of different wave fields with polarization relations given by the conventional linear AGW theory.We considered simple AGW forcing by plane wave oscillations of vertical velocity at the surface, which is similar to the assumptions made in analytical wave theory.We found height regions of the atmosphere where numerical results agree with analytical ones, and regions of their substantial disagreement.
Theoretical dispersion equation and polarization relations are widely used for developing simplified parameterizations of AGW dynamical and thermal effects in the general circulation models of the middle atmosphere.Therefore, comparisons of numerically modeled and analytical polarization relations are useful for both verifications of numerical models, and obtaining limits of analytical relation applicability and for verifying AGW parameterizations.

Numerical model
The three-dimensional numerical AGW model calculates velocity components u, v, and w along horizontal (x, y) and vertical, z, axes, respectively.The model also calculates departures of pressure p , temperature T , and density ρ from background hydrostatic stationary fields p 0 , T 0 , and ρ 0 , respectively.Gavrilov and Kshevetskii (2014a) described the set of hydrodynamic nonlinear equations used in the model.The set includes equations of continuity, momentum and heat balance.At the upper boundary z = 500 km, the conditions involve zero vertical gradients of perturbations of temperature, pressure, density and horizontal velocity, and zero vertical velocity.At the Earth's surface, the lower boundary conditions consist of zero perturbations of temperature, pressure, density and horizontal velocity (see Gavrilov andKshevetskii, 2013, 2014a, b).In this study, we assume horizontal periodicity of wave solutions: where r denotes any of the calculated variables, and L x = mλ x , L y = nλ y are the horizontal dimensions of the considered atmospheric region, m and n are integer constants, and λ x and λ y are wavelengths along horizontal axes x and y, respectively.Variations of vertical velocity w 0 = w(x, y) at the ground z = 0 generate AGWs in the model.
The used numerical scheme is analogous to the two-dimensional algorithm described by Kshevetskii and Gavrilov (2005).It is a modification of the method by Lax and Wendroff (1960).This algorithm involves the conservation laws of momentum, mass and energy.The main differences of our scheme from the classical Lax and Wendroff (1960) algorithm are the implicit approximating equations of hydrodynamics at the first half step in time, which diminish errors of description of acoustic waves (Kshevetskii, 2001a, b, c).
We use a numerical scheme similar to the two-dimensional algorithm developed by Kshevetskii and Gavrilov (2005).Used hydrodynamic equations (see Gavrilov andKshevetskii, 2013, 2014a) can be presented in the conservation law forms where s represents any one of momentum, energy or mass per unit volume, and X, Y , and Z denote components of the respective quantity fluxes along axes x, y, and z.Additionally to the model by Kshevetskii and Gavrilov (2005), the current energy balance equation contains terms describing heating caused by viscosity.The used numerical method exploits the Lax and Wendroff (1960) scheme, which approximates Eq. ( 2) with the second-order finite-difference analog where n, i, j , k and t, x, y, z are the grid node numbers and grid spacing in t, x, y, and z, respectively.This algorithm gives possibilities to select physically appropriate solutions of the equations (Lax, 1957;Richtmayer and Morton, 1967).It keeps the numerical scheme stability and allows us consideration of non-smooth solutions of nonlinear AGW equations.In addition, we exploit a staggered grid, where temperature, pressure and density are specified at the same nodes, but for the velocity components u, v, w the mesh points are half grid spacing shifted along axes x, y, and z, respectively.To compute s n+1/2 at the first time half step, we apply the implicit equation This substantially complicates simulations, but Kshevetskii (2001a, b, c) found that such structures of finitedifference schemes do not accumulate errors caused by acoustic waves.
In this study, we employ vertical profiles of background T 0 , ρ 0 , and p 0 given by the model of standard atmosphere MSIS-90 (Hedin, 1991) for average geomagnetic activity in January.The average spacing of the height grid is about 170 m, but it varies from 12 m near the ground (because of high gradients in the boundary layer) to about 1.2 km at altitudes of about 500 km, depending on inhomogeneities of vertical temperature profiles.The horizontal grid spacing is 1/60 of horizontal wavelengths taken in the wave source Eq. ( 2).Time spacing is automatically determined to guarantee stability of the numerical algorithm and is equal to 0.14 and 0.24 s for AGWs analyzed in this study with period τ = 2 × 10 3 s and horizontal phase speeds 30 and 100 m s −1 , respectively.
The numerical model involves kinematic molecular heat conductivity and viscosity, increasing versus altitude inversely proportional to the background density.We also include background turbulent heat conductivity and viscosity, taking their vertical profiles with the maxima of 10 m 2 s −1 near the ground and at altitude of 100 km and the minimum of ∼ 0.1 m 2 s −1 in the stratosphere (Gavrilov, 2013).The model does not include some effects, for example, wave dissipation caused by ion drag and radiative heat exchange, which are less important for modeling highfrequency AGWs.

AGW polarization relations
The comparisons considered in this paper used relations obtained from a theoretical model of monochromatic AGWs in the plain rotating atmosphere.Conventional linear theories suppose that wave components v , p , ρ , and T are small deviations from stationary background values v 0 , p 0 , ρ 0 , and T 0 .In agreement with Hines (1960), Beer (1974), and Matsuno and Shimazaki (1981), we can look for solutions to atmospheric wave equations for AGW spectral components in the following form where p 0s is the surface pressure; axis x is directed along horizontal wave phase velocity; σ , k and m are frequency, horizontal and vertical wave numbers; and U , V , W , P , R and are complex amplitudes of respective values.Assuming homogeneity of v 0 and T 0 , one can obtain (see Hines, 1960, andBeer, 1974) a dispersion equation relating frequency and wave numbers, which can be written in the form of where f is the Coriolis parameter, N is the isothermal Brunt-Vaisala frequency, c is the sound speed, and ω a is the highest frequency of acoustic waves, ω = σ − ku 0 .Beer (1974) found that Eq. ( 6) could be an appropriate approximation for slowly varying background temperature and wind if one uses the following expressions: where γ a = g/c p , g is the acceleration by gravity, H is the atmospheric scale height, and c p is the heat capacity at constant pressure.Applying the technique by Beer (1974), we can get the following polarization relations: where α = 1/(2H ); = (2 − γ )/(2γ H ); γ = c p /c v .Equation (8) does not allow calculating of wave amplitudes, but gives the opportunity to find their ratios.At f = 0, Eq. ( 8) are equivalent to the polarization relations obtained by Hines (1960).In a nondissipative atmosphere, according to Eq. ( 5), AGW amplitudes should grow with altitude, so that An important AGW characteristic is the wave momentum flux, the vertical component of which, F mz , is as follows: where sign denotes averaging over the wave period.

Comparisons of the numerical model and polarization relations
In this study, using the high-resolution nonlinear numerical model described in Sect.2, we simulated hydrodynamic fields produced by spectral AGW components and compared ratios of their amplitudes with those predicted by the analytical polarization relations Eqs. ( 7) and ( 8).To make simulations matching the linear AGW theory (see Eq. 5), we used nonlinear AGWs with forms of plane waves and suppose horizontally periodical distributions of vertical velocity at the Earth's surface moving along axis x of the form of where k = 2π/λ x and c x are horizontal wave-number and phase speed along the horizontal axis x in the direction of the wave propagation; W 0 is the amplitude.Equation (11) represents the plane wave of vertical velocity at the lower boundary, which may correspond to spectral components of convective and turbulent AGW sources (Townsend, 1965(Townsend, , 1966)).Medvedev and Gavrilov (1995) studied AGW generation caused by nonlinear interactions in meteorological and turbulent atmospheric processes.They found a variety of wavelengths, amplitudes and other parameters of created AGWs.
In this paper, we describe simulations for wave modes with c x = 30 m s −1 and c x = 100 m s −1 with unchanged period τ = 2×10 3 s and amplitudes W 0 = 0.3 cm s −1 .The modeling was performed beginning from the MSIS initial state (zero wave fields) and the windless background flow at t = 0, when the wave source Eq. ( 11) was triggered at the lower boundary.Gavrilov and Kshevetskii (2014a, b) demonstrated that, after triggering the wave source at t = 0, fast acoustic and very long gravity wave modes would quickly reach very high altitudes.Simulations demonstrate that, in the horizontally periodic approximation of Eq. ( 1), these initial pulses can reach altitudes of 100 km and higher in a few minutes and form quasi-vertical wave fronts analogous to those in Fig. 1a, b, c of the paper by Gavrilov andKshevetskii (2013, 2014a).These initial waves dissipate because of molecular viscosity and heat conduction.When time increases, more and more of the waves with longer vertical wavelengths are taken away by dissipation; therefore, vertical wavelengths should decrease in time at a given height in the middle atmosphere (Heale et al., 2014).After some transition time, initial AGW wave modes disappear and wave vertical structure matches the main spectral component of the wave source (Eq.11) with horizontal wave number k and phase speed c x .
To estimate AGW amplitudes in the numerical model solution, we calculated standard deviations of the corresponding wave fields over all nodes of the horizontal grid at the considered altitude.For the sinusoidal wave component, this standard deviation is equal to a half AGW amplitude.Therefore, ratios of amplitudes of horizontally homogeneous stationary sinusoidal AGWs should be equal to the ratios of the corresponding standard deviations.Simulated standard deviations of wave fields in horizontal planes located at different heights grow in time throughout transition intervals after activating the wave forcing and then tend to constant values different at each height (see Gavrilov and Kshevetskii, 2014b).In the horizontally periodical approximation of Eq. ( 1), these standard deviations are approximately equal to half wave amplitudes at large t, when the AGW process tends to become quasi-stationary.For a plane spectral AGW component with vertical wavelength λ z , the vertical group velocity is c gz ≈ λ z /τ , and the time of its energy arriving at altitude z is t e = z/c gz .For the considered main spectral components of the wave source (Eq.11) with τ = 2 × 10 3 s and average λ z ∼ 10 km for c x = 30 m s −1 , λ z ∼ 35 km for c x = 100 m s −1 .Therefore, one can get t e /τ = z/λ z ∼ 1, 6, 10 and t e /τ ∼ 0.3, 1.7, 2.9 at heights 10, 60, and 100 km, respectively, for both c x .Thus, lengths of the transition intervals are longer for smaller c x .These intervals grow with altitude and may be longer than ten wave periods at a height of 100 km.
Table 1 represents standard deviations at different altitudes calculated with the numerical model and with analytical polarization relations and their ratios for AGW with c x = 30 m s −1 .Table 1 contains simulated SDs at each altitude averaged over n model outputs during the initial transient interval t < t e (bottom part of Table 1) and for quasistationary waves t > t e (upper part of Table 1).Respective data numbers n for each altitude are presented in Table 1.Respective values obtained from analytical linear AGW theory (see Sect. 3) are calculated using average background values and are placed in the columns labeled "Lin" at each altitude in Table 1.Consideration of Fig. 5 of the paper by Gavrilov and Kshevetskii (2014b) shows that standard deviations of wave fields simulated with the numerical model vary in time due to definite variations and irregular perturbations.Standard deviations of each average numerically simulated parameter are given in Table 1.
For comparisons of numerically simulated values with analytical ones in Table 1, we use a standard t test giving the probability of the null hypothesis about equity of averages 6.6 ± 1.0 1.12 1.9 ± 0.3 1.51 3.6 ± 0.7 2.22 7.1 ± 3.9 1.7 10.1 ± 5.1 2.51 P /U , 10 −3 s m −1 3.1 ± 0.4 0.36 2.0 ± 0.3 0.47 3.6 ± 0.9 0.46 9.0 ± 5.4 0.42 12.6 ± 7.0 0.536 F mz , 10 −5 kg m −2 s −1 0.07 ± 0.04 0.32 0.12 ± 0.03 0.14 0.03 ± 0.01 0.03 0.05 ± 0.02 0.03 0.04 ± 0.01 0.03 of two irregular quantities (Rice, 2006).Approximately, the probability of equity of two respective average values in Table 1 is larger 95 %, if the difference between them is less than 1.96 multiplied by the standard deviation of the average value (Rice, 2006).In this study, we considered only cases where the standard deviations in  1, which corresponds to the initial transition time interval.Gavrilov and Kshevetskii (2014b) showed that vertical structures of transient waves are different from those predicted by the linear AGW theory during the transition interval after activating the surface wave source Eq. ( 11).The bottom part of Table 1 shows that numerically simulated wave amplitude W is smaller than that predicted by AGW theory at high altitudes, because these values refer to small t < t e , when the energy of the main wave component does not yet reach the considered altitude.Numerical and analytical amplitude ratios are also substantially different in the bottom part of Table 1 for t < t e .
In the upper part of Table 1 for quasi-stationary AGWs at t > t e , the numerically simulated AGW amplitudes W are slightly smaller than the analytical values at altitudes up to 60 km.This can be caused by small AGW dissipation at low altitudes and by partial reflections of the wave energy from inhomogeneities of background atmospheric fields in the numerical model.Wave dissipation becomes larger at altitude 100 km due to growths in cinematic viscosity and heat conductivity; therefore, simulated amplitude W in the upper part of Table 1 becomes much smaller than that predicted by nondissipative AGW theory.In addition, one can see substantial differences in numerically simulated and analytical ratios of some AGW amplitudes, which can be due to influences of dissipative effects.At low altitudes, differences in simulated and analytical ratios of AGW amplitudes can reflect the influence of lower boundary conditions.In particular, the condition u = 0 at the Earth's surface makes AGW amplitudes of horizontal velocity at low altitudes smaller than that predicted by the AGW theory for the free atmosphere.The upper part of  5.9 ± 0.4 3.24 9.4 ± 0.9 6.6 12.4 ± 0.9 16.2 7.9 ± 0.4 8.77 15.0 ± 0.8 30 P /W , 10 −3 s m −1 7.0 ± 0.5 3.59 13.1 ± 1.4 6.89 8.5 ± 0.5 11.1 6.5 ± 0.9 7.4 11.9 ± 0.4 13.8 R/ 1  1 reveals the numerically simulated AGW momentum fluxes F mz Eq. ( 10) calculated as ρ 0 u w averaged over horizontal planes at fixed altitudes and over respective time intervals.For comparisons, Table 1 also contains momentum fluxes F mz given by Eq. ( 10) and calculated from numerically simulated amplitudes W and U .The upper part of Table 1 shows that, at t > t e , wave momentum flux F mz is almost constant at altitudes 10-60 km due to relatively small dissipation and reflection of wave energy.At an altitude of 100 km, wave dissipation increases and F mz decreases, producing strong wave accelerations of the mean flow, which are proportional to the vertical gradient of F mz .In the bottom part of Table 1 for t < t e , values of F mz are much smaller than respective F mz values for t > t e , because during the initial transition interval, the energy of the main AGW modes of the wave source (Eq. 11) does not yet reach high altitudes.
Table 2 is the same as Table 1, but for AGW components with c x = 100 m s −1 , which has a longer vertical wavelength.In the upper part of Table 2 for t > t e , we have a smaller number of pairs equal with confidence 95 % (marked with bold font) than that in the upper part of Table 1.This may be connected to the stronger influence of vertical inhomogeneities of background temperature profiles on a faster AGW with a longer vertical wave number and with larger partial reflection of faster AGW energy.Stronger reflections lead to smaller amplitudes W at altitudes below 100 km in the upper part of Therefore, waves with longer vertical wavelengths can propagate faster from the surface to the upper layers and dissipate less in the middle atmosphere, where they can have amplitudes larger than those for shorter vertical wavelength AGWs (see Gavrilov and Kshevetskii, 2014b).Similar to Table 1, we have larger amounts of equal (with 95 % confidence) numerically simulated and analytical AGW parameters at altitudes 30 and 60 km.At low and high altitudes and at t < t e (in the bottom part of Table 2), numbers of equal pairs are smaller due to the influence of the lower boundary conditions, larger dissipation and AGW transience, respectively.
Tables 1 and 2 contain comparisons of the numerical results and linear polarization relations at altitudes below 100 km, where considered AGW modes are quasi-linear and almost nondissipative.At higher altitudes, growing wave amplitudes and molecular viscosity and heat conduction lead to fast growing wave-induced mean flows, which violate assumptions of conventional AGW theories and change ratios of wave amplitudes of different hydrodynamic fields.Therefore, we found poor agreement between numerical and analytical wave results above altitude 100 km and do not include them in Tables 1 and 2. These disagreements become larger with increases in amplitudes of the lower boundary wave sources due to higher nonlinear effects and faster growths in the wave-induced jet streams above 100 km.To get better agreements, improved analytical AGW theories taking into account transient processes, high wave dissipation and fast changes in background fields are required.
In the areas of Tables 1 and 2, where numerical and analytical parameters are close, one can use analytical formulae for descriptions and estimations of the wave fields.Opposite to that, areas of substantial differences between numerical and analytical AGW parameters in Tables 1 and 2 reveal regions where numerical simulations are required.
Relations of linear AGW theory are frequently used for simplified parameterizations of AGW dynamical and thermal effects for their use in the numerical models of atmospheric general circulations (e.g., Lindzen, 1981;Holton, 1983;Gavrilov and Yudin, 1992, etc.).Similar parameterizations are also being developed for highly dissipative AGWs in the upper atmosphere (e.g., Vadas and Fritts, 2005;Yigit et al., 2008).Sometimes, different parameterizations give different results.Direct numerical simulation models of atmospheric AGWs may be useful tools for testing and verifying simplified parameterizations of wave effects.

Conclusions
In this study, we performed high-resolution numerical simulations of nonlinear AGW propagation to the middle and upper atmosphere from a plane wave forcing at the Earth's surface and compared them with analytical polarization relations of linear AGW theory.Such comparisons may be used for verifications of numerical models of atmospheric AGWs.Numerical simulations show that, after triggering the wave source Eq. ( 11) at t = 0, fast acoustic and very long gravity wave modes would quickly reach very high heights.After some transition time t e (increasing with altitude), initial AGW wave modes disappear and wave vertical structure matches the main spectral component of the wave source Eq. ( 11) with horizontal wave number k and phase speed c x .The numbers of numerically simulated and analytical pairs of AGW parameters, which are equal with confidence 95 %, are largest at altitudes 30 and 60 km at t > t e .At low and high altitudes and at t < t e , numbers of equal pairs are smaller, because of the influence of the lower boundary conditions, larger dissipation and AGW transience, which can produce substantial inclinations from conditions, assumed in conventional theories of linear nondissipative stationary AGWs in the free atmosphere.
Reasonable agreements between numerically simulated and analytical wave parameters in atmospheric regions, which correspond to the scope of the limitations of the AGW theory, may be considered as evidence of adequate descriptions of wave processes by the used nonlinear numerical model.Areas of substantial differences between numerical and analytical AGW parameters reveal atmospheric regions, where analytical theories give substantial errors and numerical simulation of wave fields is required.Direct numerical simulation models of atmospheric AGWs may be useful tools for testing and verifying simplified parameterizations of wave effects.

Table 1 .
Standard deviations and their ratios for AGW with c x = 30 m s −1 calculated with the numerical model and with analytical polarization relations (labeled as Lin) at different altitudes averaged over the initial transient interval t < t e and for quasi-stationary waves t > t e .Bold font shows the data pairs equal to probabilities larger than 95 %.
Table 1 are smaller than 0.15 of the respective average values.Pairs of AGW parameters, which we can consider equal with confidence larger than 95 %, are marked with bold font in Table 1.The numbers of those pairs are largest in the upper part of Table1at altitudes 30 and 60 km, which correspond to quasi-stationary AGWs in the free atmosphere considered in conventional AGW theory described in Sect.3.Reasonable agreements between simulated and analytical wave parameters in atmospheric regions, which correspond to the scope of the limitations of the nondissipative linear AGW theory, may be considered as evidence of adequate descriptions of wave processes by the used nonlinear numerical model.Many numerically simulated AGW parameters do not match the respective analytical values in Table1.No matches are in the bottom part of Table Table 1 for t > t e shows that the best agree- www.geosci-model-dev.net/8/1831/2015/Geosci.Model Dev., 8, 1831-1838, 2015

Table 2 .
Same as Table 1 but for AGW with c x = 100 m s −1 .
Table 2 compared to that in Table 1.On the other hand, W at altitude 100 km in the upper part of Table 2 is larger than that in Table 1 due to smaller dissipation of longer AGWs.