On Linear Instability of Atmospheric Quasi-hydrostatic Equations in Response to Small Shortwave Perturbations

A set of 3-dimensional atmospheric-dynamics equations with quasi-hydrostatic approximation is proposed and justified with the practical goal to optimize atmospheric modelling at scales ranging from meso meteorology to global climate. Sound waves are filtered by applying the quasi-hydrostatic approximation. In the closed system of hydro/thermodynamic equations, the inertial forces are negligibly small compared to gravity forces, and the asymptotically exact equation for vertical velocity is obtained. Investigation of the stability of solutions to this system in response to small shortwave perturbations has shown that solutions have the property of shortwave instability. There are situations when the increment of the perturbation amplitude tends to infinity, corresponding to absolute instability. It means that the Cauchy problem for such equations may be ill-posed. Its formulation can become conditionally correct if solutions are sought in a limited class of sufficiently smooth functions whose Fourier harmonics tend to zero reasonably quickly when the wavelengths of the perturbations approach zero. Thus, the numerical scheme for the quasi-hydrostatic equations using the finite-difference method requires an adequately selected pseudo-viscosity to eliminate the instability caused by perturbations with wavelengths of the order of the grid size. The result is useful for choosing appropriate vertical and horizontal grid sizes for modelling to avoid shortwave instability associated with the property of the system of equations. Implementation of pseudo-viscosities helps to smoothen or suppress the perturbations that occur during modelling.


INTRODUCTION
The advantage of quasi-hydrostatic equations is the efficiency of numerical calculations of atmospheric circulation by filtering the sound waves. The pseudo-incompressible approximation div v = 0 or the anelastic approximation divρ v = 0 is also used in certain set of equations for this purpose, but they violate the mass conservation equation. Also, the anelastic equation set is suitable in smallscale models, and it plays an analogous role as the quasi-hydrostatic approximation does in large-scale models [1]. Different acoustic-wave-filtered systems of equations used in the atmospheric sciences are compared [2]. Besides these physical approaches, for the purpose of sound waves filtering one can also apply numerical techniques such as the variational approach [3] and semi-implicit temporal integrators [4]. However, as discussed in [5], the aspect ratio between vertical and horizontal scales significantly influences the dispersion relation that, in turn, determines the stability of numerical calculation. As stated in [6], almost all global models are based on the vertically hydrostatic equations. The vertically quasi-hydrostatic approximation is adopted in the atmospheric models for studying global long-term non-extreme climatic processes, while vertical non-hydrostatic terms (vertical inertia force) may be significant when taking into account the small-scale processes such as cloud microphysics or phase change of water in the air [7,8].
Along numerical calculations with finite difference methods for large scale global models the smallscale numerical non-physical perturbations arise. The growth of such perturbations is associated with the grid sizes, and instability may occur during calculation. That is why long-term calculations of the global models are restricted by numerical blow-ups (instability). The atmospheric predictability problem has been discovered and studied for a long time since the early works [9,10]. Lorenz Edward proposed three approaches to this problem [11], within which the dynamical approach is adopted in the analysis of present paper. The influence of scales on the stability (or predictability) is discussed in more recent works [12,13].
We introduce a system of 3-dimensional atmospheric-dynamics equations with a quasi-hydrostatic approximation for climate modelling. Based on this system, we investigate the linear instability in response to shortwave perturbations. The present paper focuses on this problem with the practical goal to assure realistic numerical modelling, with emphasis on the scale effects. In this section, we introduce the primitive equations used in atmospheric sciences in general. The derivation of quasi-hydrostatic equations with modification [14] in the process of derivation of the vertical velocity is introduced in Section 2. In Section 3, the system of equations for small perturbations in dimensionless form and the cubic characteristic equation with complex coefficients are derived. In section 4, the stability of solutions to the system of equations with response to shortwave perturbations is analyzed through solving the cubic characteristic equation. Section 5 focuses on the elimination of instability caused by the shortwave perturbations with the aid of pseudo-viscosities.
As the total mass of the atmosphere is located in the layer with thickness, H, of order 10 km, which is negligibly small compared to the Earth's radius (R ≈ 6000 km). So, the atmospheric dynamics outside polar zones can be considered in a local Cartesian coordinate system x, y, z, t, where x and y are directed horizontally along latitude and longitude, respectively, z is the height over the surface, and t is time.
Due to the smallness of the ratio H/R, the curvature terms containing H/R can be neglected. Thus, the equations of conservation of mass and momentum for inviscid air read [15][16][17]: where ρ and v = (v x , v y , v z ) are density and velocity of air, p is the pressure, f = (f x , f y , f z ) is the Coriolis force vector, and g is the gravity acceleration (g ≈ 9.81 m s −2 ). We consider air as an ideal gas, from whose state equation for the pressure p and the specific internal energy u are determined by the specific gas constant R ≈ 286 m 2 s −2 K −1 and the adiabatic index γ ≈ 1.4: where T is the air temperature, and c v is the specific heat capacity at constant volume. The equation of conservation of energy in the form of the thermodynamic equation reads: where Q * is total heat influx into unit volume due to the turbulent heat flux, radiation heat absorption, and latent heat of wataer phase transfer. The equations (1.1)-(1.6) form a closed system. However, for modelling long-term large-scale climatic processes, the quasi-hydrostatic equation is used instead of the vertical momentum equation (1.4). Thus, for the closure of the system of equations with quasi-hydrostatic approximation, it is necessary to derive an equation to estimate the vertical velocity.

THE SYSTEM OF QUASI-HYDROSTATIC EQUATIONS
We consider non-extreme hydrodynamic and thermodynamic processes in the Earth's atmosphere with climatic scales, i.e., time scale, τ , horizontal length scale, L x , L y ∼ L hor , vertical length scale, L z ∼ L ver , horizontal velocity scale, V x , V y ∼ V hor , and vertical velocity scale, V z ∼ V ver , with the following values: τ ≥ 10 2 s, L hor ≥ 10 3 m, L ver ≥ 10 2 m, V hor ∼ 10 1 m/s, V ver ∼ 10 0 m/s. (2.1) With the use of these scales, we introduce dimensionless values of time, t, horizontal coordinates, x, y, and horizontal velocities, v x , v y , vertical coordinate, z, and vertical velocity, v z : As a result, the horizontal, vertical, and Coriolis accelerations can be evaluated as dv Here Ω ≈ 0.727 × 10 −4 s −1 is Earth's angular velocity of rotation; A hor , A ver , and A cor are scales (characteristic values) of the horizontal, vertical, and Coriolis accelerations, respectively. According to (2.2), the dimensionless values (denoted with over-line above) have the order of unity. Thus, for nonextreme climate processes with scales (2.1), the following dimensionless parameters are small Thus, the inertial forces are negligibly small compared to gravity, i.e., ε → 0. For the scales (2.1) and small value of ε, the quasi-hydrostatic approximation can be adopted, then the pressure of air at a certain height is determined by the total mass above in the vertical column 1) : Taking into account the quasi-hydrostatic equation for pressure (2.4), we obtain Substituting this expression into (2.5), we get dp dt Then the vertical velocity is evaluated using the continuity equation (1.1) and thermodynamic energy equation (1.6) following [14] ∂v z ∂z The horizontal gradients of pressure is estimated from (2.3) and horizontal momentum equa- then the fourth and fifth terms on the right side of (2.7) can be estimated where Ma ≡ V hor /C is the Mach number of the horizontal motion and C = γp/ρ = 300−350 m/s is the sound speed. The magnitude of Ma 2 is small under the scales of (2.1). Furthermore, we show that the asymptotics Ma 2 → 0 as ε → 0 corresponds to the following relation: Here we consider that the magnitude of L hor does not exceed the magnitude of C 2 /g ∼ 10 4 m by orders. Thus, according to (2.8) and (2.9), the terms of horizontal advection of pressure gradients are negligibly small as ε → 0. Then by neglecting the small terms of ε in equations (2.6) and (2.7), we have dp dt = gṀ, (2.10) The equation for the vertical velocity v z in the form (2.7) was presented in [18,19], but in recent years, it has hardly been mentioned. A similar but different equation for the adiabatic regime (Q * = 0) was also obtained in [20]. However, these works do not show that the terms (2.8) are on the order of ε, i.e., negligibly small when the inertial forces are negligibly small in comparison with gravity force.
In total, we have the following system of hydro-and thermodynamic equations for the inviscid air in the field of gravity with the quasi-hydrostatic approximation along the vertical coordinate [14]: Here f = 2Ω sin θ is Coriolis parameter (θ is the latitude). In this system, the temperature is calculated from the state equation of air (1.5). In [21] a different diagnostic equation for v z is used in the set of equations, and temperature T is directly obtained from the thermodynamic equation, which serves as a prognostic equation. Although the thermodynamic equation (1.6) is not in our set of equations, it is utilized to derive the diagnostic equation for the vertical velocity. Following (1.1)-(1.4) and (1.6) for inviscid air and accounting that the velocity is perpendicular to the Coriolis force ( f · v = 0), the energy conservation law writes For the scales (2.1) and climatic processes, ΔT ∼ 10 K, then the term Q * dominates on the right side of (2.13). Thus, the omission of v z in such scales and processes does not affect energy conservation during the long-term numerical calculation.

Dimensionless Form of Vertical Quasi-Hydrostatic Equations
In addition to (2.2), we introduce the following dimensionless variables: where ρ 0 is the characteristic density at the surface of the Earth (z = 0).
The system of equations (2.12) is non-linear or quasilinear due to the convection terms in continuity equation and horizontal momentum equations. With the use of (3.1) the system (2.12) can be represented in the dimensionless form where U is a column matrix of dependent variables (superscript T denotes transposed matrix), B t , B x , B y , B z , B are matrices of coefficients defined by dependent variables in U (see Appendix A.1).
In the dimensionless system of equations the following dimensionless parameters are defined: they determine the external forcings, i.e., heat flux, gravity, and Coriolis force, respectively. If we use the following values for scales: here ΔT is the characteristic value of temperature change in period τ , then the characteristic value of total heat flux can be evaluated by τ Q * = ρc p ΔT . Therefore, we have the following estimations for the dimensionless scales (3.3) 3.2. Stability of Small Shortwave Perturbation to the Basic Solution We consider some solution to the system (3.2): (3.5) and name it the basic solution. Assume that there exists another solution U (d) to (3.2) with a small perturbation U to U: Substituting (3.6) into the system (3.2) and subtracting equations (3.2) for the basic solution U and neglecting the nonlinear terms of the second and third-order of δ, we obtain a system of quasi-linear differential equations for the perturbation U : where B as B t , B x , B y , B z is a matrix, which consists of parameters depending on the basic solution (unperturbed) to the system (3.2), and the column matrix F is determined by the perturbation of the heat flux Q (see Appendix A.2).
We consider small harmonic perturbations as the solution to the system (3.7), which can be represented in the following form with the help of imaginary unit i = √ −1 and complex variables: We assume that the wavelengths of shortwave perturbations are much shorter than the characteristic length scale, and the heat flux perturbation is zero (Q = 0): For such perturbations, we have F = 0. The coefficients B t , B x , B y , B z , B with parameters of the basic solution U vary within distances of L hor and L ver , which are much larger than the wavelengths of perturbations l x , l y , l z . Thus, these matrices of coefficients can be regarded as constants. Then system (3.7) becomes a homogeneous system of linear equations with constant coefficients, and among its solutions for t > 0, there exist solutions in the form of traveling wave: Here the real part ω is the frequency of the perturbation, and the imaginary part ω * * is the increment of the perturbation amplitude. Substituting (3.8) into the system of equations (3.7), we obtain a homogeneous linear system of algebraic equations Therefore, for the homogeneous system of algebraic equations DA * = 0 having a nonzero solution, the condition detD = 0 reduces to a cubic characteristic equation in ω * : where b * , c * , d * are complex parameters (see Appendix A.3) defined by the basic solution to system (3.5).
Obviously, if among the roots of the characteristic equation (3.10) there exists such a root ω * = ω + iω * * that the increment ω * * > 0, then the amplitude of perturbation in (3.8) increases with time. This indicates the instability of the investigated solution (3.5). If at the same time then such a solution holds absolute instability, and the Cauchy problem of the system (3.2) is illposed [22]. The explanation to the energy principle of the absolute instability (or blow-up) is: 1) the perturbations in the form of (3.8) occupies the whole space, but the blow-up occurs at a certain spacetime point; thus, the infinite amplitude of the perturbation at this point does not indicate the infinite total energy of perturbation; 2) such shortwave perturbations occur due to numerical errors. In contrast, the negative increment ω * * < 0 ensures stability. In the case of neutral stability of the linear approximation (ω * * = 0), the conclusion about the instability of solution (3.5) requires the study of the nonlinear behavior of the system (2.12).

LINEAR INSTABILITY OF BASIC SOLUTIONS 4.1. Instability of Resting-State Solution
We consider the regime when there is no motion and heat influx in the basic solution, and the density is uniform horizontally: The horizontal k hor and vertical wavenumbers k ver are defined With resting-state conditions (4.1), the cubic characteristic equation (3.10) is simplified as Here N is the dimensionless value of the buoyancy frequency N , C is sound speed. For scales (3.4) and standard atmosphere [23] we have the following estimates: The dependency of N 2 − G 2 on height for the standard atmosphere is presented in Fig. 1. It is shown that N 2 − G 2 ∼ 10 −4 s −2 , and at height z ≈ 4 km, the sign of N 2 − G 2 changes. The small Coriolis force can take place at near-equatorial latitudes (θ 1) and a height not close to 4 km. In such situation for From the general formulae of roots (Appendix A.4), we obtain the expressions of nonzero roots of (3.10) for small Coriolis force: When k hor → ∞ and k ver → ∞, there are various asymptotics, including If, in contrast to (4.4), the Coriolis force prevails: which can take place at high latitudes or altitudes near 4 km, then according to Appendix A.4 we get As in (4.6) for k hor → ∞ and k ver → ∞, there exist different asymptotics, including From (4.6) and (4.8) it follows that the resting-state of the system of hydrodynamic equations with the vertical quasi-hydrostatic approximation (2.12) is unstable, in particular, in case of shortwave perturbations with vertical wavelengths l ver = 2π/k ver longer than horizontal (l hor = 2π/k hor , k ver , k hor 1, and k ver < k hor ), is absolutely unstable. For perturbations when the vertical wavelength is many times smaller than the horizontal wavelength (l ver l hor , k ver k hor 1), the resting-state tends to be neutrally stable.
The positive values of the increment of the amplitude of perturbation ω * * at resting-state for k z = 15 and k z = 150 with the horizontal wavelengths k x = k y = √ 2 2 k hor are shown in Fig. 2 by the lines with indicator number 0.
The main conclusion is that the shorter the horizontal wavelengths (the larger k hor 1), the faster the amplitude of perturbation at resting-state develops. Moreover, the shorter the vertical wavelengths (the larger k ver 1), the more slowly the perturbation at resting-state grows.

Instability of Vertical Quasi-Hydrostatic System with Different Approximations
Shortwave perturbation at resting-state of the atmosphere with a vertical quasi-hydrostatic approximation is also studied in [24]. However, in the original system the following equation (thereinafter, the Arakawa inexact approximationt) is used instead of the asymptotically exact equation (2.10), which is equivalent to ∂M ∂t =Ṁ + ρv z . (4.10) Note that the equation for perturbation at resting-state corresponding to (4.9) has the form If instead of equation (2.10) or (4.10) the Arakawa inexact approximation (4.9) is adopted, then instead of (3.10) the characteristic equation becomes: Nonzero roots of this equation are determined by the following formulae: It is clear that these expressions corresponding to the usage of the Arakawa inexact approximation, differ from (4.5) and (4.7). In the case of the predominance of Coriolis force, this difference is fundamental, namely, instead of absolute instability in (4.8), neutral stability follows from (4.11).
In the case of small Coriolis force, the Arakawa inexact approximation gives absolute instability, with following asymptotics Now we consider other approximations instead of the asymptotically exact equation for the vertical velocity v z (2.11). The simplest approximation for v z is zero vertical velocity, in particular, for the perturbation of vertical velocity, v z = 0. (4.13) The corresponding characteristic equation at resting-state (4.1) has the form The root with positive increment ω * * when f < Gk hor /k ver can lead to absolute instability (3.11). In particular, when the Coriolis force is small and can be neglected, (4.14) The Holton inexact vertical quasi-hydrostatic approximation [16] with equation of the constant local pressure (thereinafter, the Holton inexact approximation) can also be used dp dt Using this approximation instead of equation (2.10) or its equivalent (4.10), the characteristic equation becomes and its roots take the following form: It follows that the solution of the resting-state to the system of equations with the Holton inexact vertical quasi-hydrostatic approximation is neutrally stable (ω Otherwise, when i.e., the density increases vertically, such basic solution of resting-state is unstable for small Coriolis force f < N κ 4 , and analogously to (4.14), absolutely unstable.
One can also adopt the vertical quasi-hydrostatic approximation with the quasi-incompressibility dρ dt = 0, from which together with the continuity equation, instead of (2.11) the vertical velocity v z is obtained by the following equation For such a system, the characteristic equation at resting-state (4.1) has the same form (4.16), as for the Holton inexact approximation (4.15).
The vertical quasi-hydrostatic approximation with constant local density is proposed in [25] (thereinafter, the Marchuk inexact approximation): Together with the continuity equation, the equation for the vertical velocity v z is obtained as following The characteristic equation for this equation at resting-state (4.1) has the form which corresponds to the neutral stability.
The Marchuk inexact approximation (4.18) and (4.19), the Arakawa approximation (4.9), the zero velocity approximation (4.13) and the Holton inexact approximation of constant local pressure (4.15) are valid only when the vertical velocity of air is comparatively small. Only the Marchuk inexact approximation leads to a neutrally stable solution at resting-state for all wavenumbers of perturbations.
The shortwave instability at resting-state is also studied in [26] for a two-dimensional system (v y = v y = 0) with approximation (4.17), and taking into account of the perturbation of the heat flux Q = 0, which is assumed to be proportional to the perturbation of vertical velocity v z . The characteristic equation is close to a specific occasion (v y = v y = 0, k y = 0) of equation (4.16).
From (4.6), (4.8), (4.12) and (4.14), it can be seen that the asymptotics for the increment of perturbations for various approximations is determined by the coefficients of the relations between k hor and k ver , namely

Instability of One-Dimensional Vertical Motion
A one-dimensional vertical atmospheric model is used to analyze physical processes in the formation of climate. In such a case, the distribution of vertical velocity is defined by the heat flux Q and horizontal inflowṀ , which should be set parametrically, and the perturbation spreads only in the vertical direction: Thus, the system of equations (2.12) is simplified as HereṀ (t, z) is regarded as known variable. In this case, the cubic equation (3.10) degenerates into a linear one, and its root is For a standard atmosphere, ∂ρ/∂z < 0, one always has ∂ρ/∂z M ∼ −10 −2 , from which it follows that ω * * ≈ − ∂vz ∂z . Thus the shortwave stability for a one-dimensional atmosphere (ω * * < 0) takes place when the vertical velocity increases with height (∂v z /∂z > 0), this happens with the presence of heating (Q > 0) and horizontal outflow above the given height (Ṁ < 0). Otherwise, when cooling occurs (Q < 0) and there is inflow above (Ṁ > 0), then the vertical velocity decreases with height (∂v z /∂z < 0), and such state is unstable (ω * * > 0).

Instability of Solution with Motion
Using the roots of the cubic equation (3.10), we study the stability of three-dimensional perturbation to a solution with motion with equal horizontal wavenumbers. Then according to (4.2) To calculate the roots of (3.10), we set the parameters of the basic solution with motion as v the values of ρ, M , ∂ρ/∂z are chosen according to the standard atmosphere [23] at height z = 10 km.
Numerical values of increments ω (j) * (j = 1, 2, 3) depending on k hor are shown in Fig. 3, they are calculated from the cubic equation (3.10) for fixed vertical wavelengths k ver = 15, 50, 150, and parameters given by the basic solution above. Numerical indicators correspond to the number of roots j = 1, 2, 3. It is clear that the positive increment ω * * grows with the increase of k hor . It follows from Appendix A.3 that the values k hor /k ver strongly change the coefficients of the cubic equation (3.10). For an increase of k ver , the growth of positive increment ω * * greatly weakens. In Fig. 2, in addition to the resting-state (numerical indicator 0), the positive values of increment ω * * depending on k horr , k ver for two sets of horizontal velocity are also shown. The lines with indicators 1 and 2 refer to the state with horizontal velocities given in the beginning of this subsection and the state with horizontal velocities twice as larger, respectively. The dimensionless horizontal velocity in the order of unity affects the increment for wavenumbers k hor ∼ 10 1 , yet the horizontal velocity does not strongly affect the increment ω * * for shortwave perturbations (k hor > 10 2 ).

Shear and Bulk Pseudo-Viscosities
This analysis does not consider the influence of boundary conditions. The boundaries can prohibit the energy influx from outside, which contributes to the kinetic energy of perturbations and the growth of perturbation amplitude. Nevertheless, the shortwave instability, of course, is one of the reasons of the possible instability of finite-difference calculations.
While modelling (3.2) by the finite-difference method, small shortwave perturbations can be generated, the minimum wavelength of which l min i is defined by the grid size Δ i : where N x , N y , N z are the numbers of nodes at characteristic length L x , L y , L z , in directions x, y, z, respectively. The terms of viscosity in the horizontal momentum equations are relatively small in comparison with other terms. However, the pseudo-viscosities, which are much greater than real viscosity of air, can be added to eliminate the instability due to numerical shortwave perturbations. They influence the horizontal momentum equations of the perturbations, and then change the form of algebraic equations (3.9) through where μ and λ are coefficients of shear and bulk pseudo-viscosities, respectively. We set the dimensionless values of these coefficients growing with increasing wavenumbers:  Here μ a is the viscosity of the air.

Influence of Pseudo-Viscosities for Resting-State Solution
Similar as (3.10) the characteristic equation for a resting-state solution, taking into account the shear pseudo-viscosity in accordance with (5.1), turns out to be A similar cubic equation for a restiing-state solution with bulk pseudo-viscosity λ has the following form It is clear that for m > 0 (see (5.2)) and k i 1 (i = x, y, z), we have However, an increase in pseudo-viscosity can not achieve the point that all three roots (ω (j) * (j = 1, 2, 3)) of the characteristic equation have negative increments ω (j) * * < 0 (j = 1, 2, 3) so that the shortwave perturbation is damped and disappears eventually. One of the roots always has a positive increment ω (3) * * > 0.

Influence of Pseudo-Viscosities for the Solution with Motion
For moving air, the pseudo-viscosities can not suppress the global shortwave perturbations of solutions to the asymptotically exact system of equations (2.12) with vertical quasi-hydrostatic approximation. This can be explained by the fact that the perturbations occupy the entire space, i.e., the perturbations are global and have infinite energy. Real perturbations from the numerical calculations on the difference grids are local and have finite energy. Such shortwave perturbations can be suppressed by pseudo-viscosities (see Section 5.4). Now we consider the influence of pseudo-viscosities (here we consider shear pseudo-viscosity) on the stability for a basic solution with motion using the Marchuk inexact approximation. Similarly to (3.10), we obtain the following characteristic equation: and c 1 (μ) are the complex parameters determined by the basic solution and shear pseudoviscosity μ. The first root corresponds to neutral linear stability (ω (1) * * ), and the rest two roots correspond to linear stability, i.e., for k i → +∞(i = x, y, z), it follows ω (2,3) * * → −∞ (see Fig. 5). That is to say, the solution to the system with the Marchuk inexact approximation for moving air can be unstable, but for shortwave perturbations with the existence of pseudo-viscosity, it turns to be stable. It means that the shear pseudo-viscosity μ (by analogy, as well as λ) can suppress the global shortwave perturbations, and this is the advantage of the Marchuk inexact approximation (4.18) and (4.19) in comparison with the exact approximation.

Influence of Local Perturbations on Instability
Instead of perturbations (3.8), whose amplitude is identical in space, we consider such perturbations that are damped in space for x, y, z > 0 with complex wavenumbers k j * = k j + ik j * * (j = x, y, z): In Fig. 6 the values of increment ω j * * (j = 1, 2, 3) for resting-state solution (4.1) are shown for k x * * = k y * * = k z * * = 0 and 20, with k ver = 15, n = 1.5, c μ = 1. The thin lines correspond to k j * * = 0, and thick lines k j * * = 20. Numerical pointors indicate the root number j = 1, 2, 3. It is shown that the shear pseudo-viscosity can stabilize the local shortwave perturbations (k x , k y , k z 30) for x, y, z > 0, since all roots have negative increments ω j * * < 0 (j = 1, 2, 3). It should be noted that for large values k j * * 1, there are large gradients of perturbation parameters (like impulse), and this leads to an absolute instability (when k j * * → ∞) in response to relatively longer wavelength range (k hor < 30, for fixed k ver = 15).

CONCLUSIONS
In this paper, an equation for vertical velocity refines and simplifies the equation obtained by [18]. Different from the approximations used in other works such as [16,25], this equation is asymptotically exact for the system of hydrodynamic equations with small inertia forces compared to the gravity force. The advantage of this system is the absence of sound waves. Indeed, when the density field, ρ(x, y, z, t), is given, the pressure field, p, can be calculated without time stepping via the quasi-hydrostatic equation. Then having ρ and p, the temperature field, T , is calculated through the equation of state. The absence of sound waves allows for much larger time steps in modelling due to the absence of the Courant-number restriction in the vertical dimension.
In our analysis the influence of boundary conditions on the instability problem is excluded. Though the boundary conditions in the framework of the boundary value problem can significantly limit the growth of the perturbation, the shortwave instability is one of the reasons of the possible instability of finite-difference calculations.
As a result, the solution to the system of atmospheric dynamics equations with vertical quasihydrostatical approximation and without viscosity is known to be unstable for shortwave perturbations occupying the entire space. Even the resting-state solution owns the shortwave instability. The increment number of the perturbation growth is related to the ratio of horizontal to vertical wavenumber of the shortwave perturbations. Such shortwave instability eventually results in the instability introduced by numerical truncation and finite-difference approximation of the derivatives error. It shows that not only usual solutions presenting three-dimensional flow patterns but even the resting-state solution is unstable under shortwave perturbations. Moreover, the larger the ratio of horizontal wavenumber k hor to vertical wavenumber k ver (or the larger the ratio of vertical grid size Δ ver to horizontal grid size Δ hor ), the larger the increment growth of the amplitude of perturbations. Particularly, an infinitely large ratio of vertical grid size to horizontal grid size leads to the absolute instability. Thus, when decreasing the horizontal grid size to achieve better accuracy, one should also decrease the vertical grid size to keep small the ratio κ 2  For the one-dimensional vertical motion (one-column model) in the standard atmosphere where ∂ρ/∂z < 0, the shortwave stability depends on the vertical velocity profile, v z (z) or, more specifically, on the sign of ∂v z /∂z, which is related to the heat input, Q, and horizontal mass inflow,Ṁ . Then, due to the heating and horizontal outflow above the given position the condition of the shortwave stability is ∂v z /∂z > 0.
The pseudo-viscosities taken proportional to the wavenumber of perturbations reduce the increment in the amplitudes of perturbation, so that numerical solutions to the asymptotically exact quasihydrostatic system become more stable. In this context, global perturbations have non-negative increment causing numerical instability. However, in the case of only local perturbations, implementation of pseudo-viscosities allows for making all increments negative, thus yielding practically sound stable solutions.
At the large ratio of the vertical to horizontal grid size, the solution with motion is also unstable for other (inexact) vertical quasi-hydrostatic approximations, namely, those with constant local density (Marchuk, ∂ρ/∂t = 0), constant local pressure (Holton, ∂p/∂t = 0), or quasi-incompressibility (dρ/dt = 0). In contrast to other known differential operators, the inexact vertical quasi-hydrostatic system with constant local density (Marchuk, ∂ρ/∂t = 0) assures the neutral stability for the restingstate solution. And introducing pseudo-viscosities can suppress global perturbations for the solution with motion.
The result about shortwave instability can guide the choice of vertical and horizontal grid sizes for modelling, and proper usage of pseudo-viscosity can reduce the predictability problem caused by shortwave perturbations during calculation.
The code for generating the figures in this paper can be downloaded by link https://doi.org/10.5281/ zenodo.3831455.

Matrices of Coefficients of System of Equations for Perturbation in Dimensionless Form (3.7)
3. Complex Coefficients of the Cubic Equation (3.10)