Quasi-hydrostatic equations for climate models and the study on linear instability

An advanced “quasi-hydrostatic approximation” of 3-dimensional atmospheric-dynamics equations is proposed and justified with the practical goal to optimize atmospheric modelling at scales ranging from meso meteorology to global climate. For the vertically quasi-hydrostatic flow with inertial forces negligibly small compared to gravity forces, the asymptotically exact equation for vertical velocity is obtained. In the closed system of hydro/thermodynamic equations, the pressure is determined by the total air mass above, so that mass instead of pressure is considered as a dependent variable. In such a system, the 5 sound waves are filtered, though the horizontal inertia forces are taken into account in the horizontal momentum conservation equations. The major practical result is an asymptotically exact equation for vertical velocity in the quasi-hydrostatic system of the atmospheric dynamics equations. Investigation of the stability of solutions to the system in response to small shortwave perturbations has shown that solutions have the property of shortwave instability. There are situations when the growth rate of the perturbation amplitude tends to 10 infinity. 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 15 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.

sound waves are filtered, though the horizontal inertia forces are taken into account in the horizontal momentum conservation equations. The major practical result is an asymptotically exact equation for vertical velocity in the quasi-hydrostatic system of the atmospheric dynamics equations.
Investigation of the stability of solutions to the system in response to small shortwave perturbations has shown that solutions have the property of shortwave instability. There are situations when the growth rate of the perturbation amplitude tends to 10 infinity. 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 15 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 wave. However, as discussed in Orlanski (1981), the aspect ratio between vertical and horizontal scales significantly 20 influences the dispersion relation that, in turn, determines the stability of numerical calculation. Most global climate models are based on a system of dynamic equations in quasi-hydrostatic approximation White and Bromley (1995); White et al. (2005).
The atmospheric predictability problem has been discovered and studied for a long time since the early works like (Lorenz, 1969a(Lorenz, , 1982. Lorenz Edward proposed three approaches to this problem (Lorenz, 1969b), within which the dynamical ap-charov, 1982; Batchelor and Batchelor, 2000;Gill, 2016;Loitsyanskii, 1970;Sedov, 1997;Holton and Hakim, 2012;Nigmat-40 where ρ and v = (v x , v y , v z ) are density and velocity of air, p is pressure, Π ij is the viscous (turbulent) stress tensor, f = (f x , f y , f z ) is the Coriolis force vector, and g is the gravity acceleration (g ≈ 9.81m s −2 ). 50 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: (1.6) Here Q * is total heat influx into unit volume due to the turbulent heat flux, q, radiation heat absorption, Q, and latent heat consisted of evaporation (condensation) heat of droplets, J (e) l (e) , and the heat due to melting (solidification) of ice particles, J (m) l (m) . Here, it is taken into account that the horizontal turbulent heat transfer is much smaller than the vertical one, so that the horizontal heat transfer is due to advection, determined by horizontal velocities v x and v y .

60
Using the state equation(1.5), the thermodynamic energy equation (1.6) can be rewritten in terms of the following relation between the time derivatives of pressure and density: (1.7) The equations (1.1) -(1.7) form a closed system. However, for filtering the sound wave, the quasi-hydrostatic equation is used instead of the vertical momentum equation (1.4). Thus, for the new system of equations with quasi-hydrostatic approxi-65 mation, it is necessary to derive an equation to estimate the vertical velocity.
2 Asymptotics of hydrodynamic equations with vertical quasi-hydrostatic approximation

Meteorological and climate scales
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 , 70 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 .
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, as well as inertial forces, can be evaluated as: Here Ω ≈ 0, 727 × 10 −4 s −1 is Earth's angular velocity of rotation. Here A hor , A ver , and A cor are scales (characteristic values) 80 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 non-extreme climate processes with scales (2.1), the following dimensionless parameters are small Thus, the inertial forces are negligibly small compared to gravity.

85
It makes sense to consider the asymptotics of the equations of hydro-and thermodynamics of the atmosphere (1.1) − (1.5), when ε → 0, (2.5) for decreasing values of V hor , V ver and increasing values of L hor , L ver , τ . Here, the viscous forces are small and not considered.

90
For the scales (2.1) and small value of , from the momentum equations (1.2) -(1.4) we get 1 As a result, the pressure of air at a certain height is determined by the total mass above in the vertical column: (2.7) 1 Here and below O(ϕ) means, that on the surface of Earth (z = 0). In particular, one can take H ≈ 40 km, then x, y, 0). (2.8) The distribution of vertical velocity is evaluated using the continuity equation (1.1) and thermodynamic energy equation (2.9) Taking into account the quasi-hydrostatic equation for pressure (2.7), we obtain (2.10) Substituting this expression into (2.10), we get As a result, substituting (2.12) into (2.9), we obtain the equation for vertical velocity: The first two terms on the right side of the equation above are estimated as follows: (2.14) Using the mean-value theorem, we estimate the third term on the right side of (2.13): whereρ(z) andρ(z) are the mean densities of the considered integrals over vertical coordinate from z to H in the denominator and numerator, respectively. Here we haveρ(z)/ρ(z) = O(1).
To estimate the fourth and fifth terms on the right side of (2.13), the horizontal gradients of pressure is estimated from ( where Ma ≡ V hor /C is the Mach number of the horizontal motion and C = (γp/ρ) 1 /2 = 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.17) and (2.18), the terms of horizontal advection of pressure gradients are negligibly small as → 0 .
As a result, the equation (2.12) and the differential equation following from it (2.13) take the following form: (2.20)

Asymptotic system of equations in case of small inertial foeces
Given the fields of ρ, v x and v y , the fields of M andṀ can be calculated using (2.7) and (2.11): was presented in the book (Lorenz and Lorenz, 1967), but in recent years, it has hardly been mentioned. A similar but different equation for the adiabatic regime (Q * = 0) was also obtained in (Eliassen, 1949). However, these works do not show that the last two terms of (2.22) are in the order of , i.e., the inertial forces are negligibly small in comparison with gravity force.

140
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: (2.23) Here f = 2Ω sin θ is Coriolis parameter (θ is the latitude). In this system, the temperature can be calculated from the equation of state of air: If the density, ρ(x, y, z, t), and surface pressure, p(x, y, 0, t), are known, the pressure, p(x, y, z, t), can be calculated by integration of the quasi-hydrostatic equation (2.7). Then temperature T (x, y, z, t) can be obtained via the state equation (2.24).
The vertical velocity, v z , and the pressure, p, are inertialess as they are determined by the distributions of density ρ, horizontal velocities, v x , v y , at the same moment of time. Thus, for non-extreme meteorological processes in the atmosphere with scales 155 (2.1), the inertia effects only through the horizontal velocities.
Until now, we have proven the following theorem. values of the inertial forces (horizontal,vertical and Coriolis) to the gravity force, approaches to zero (see (2.5)).

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.23) 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.23) can be represented in the dimensionless form The derivative ∂ 2 T /∂z 2 , as a component of heat source Q * through eddy heat transfer, is not considered as a part of the differential operator (3.2).
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 175 scales: here ∆T is the characteristic value of temperature change in period τ , then accounting that τ Q * = ρc p ∆T , we have the following estimations for the dimensionless scales (3.3)

Stability of small shortwave perturbation to the original solution
We consider some solution to the system (3.2): which is named after the original solution. And exists another solution U (d) to (3.2) with a small perturbation U to U: Substituting (3.7) into the system (3.2) and subtracting equations (3.2) for the original 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 original solution (unperturbed) to the system (3.2), and the column matrix F is determined by the perturbation of the heat flux Q (see Appendix B).

190
We consider small harmonic perturbations as the solution to the system (3.8), 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 original solution U vary within distances of L hor and L ver , whch are much larger than the wavelengths of perturbations l x , l y , l z . Thus, these matrices of 200 coefficients can be regarded as constants. Then system (3.8) 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: Substituting (3.11) into the system of equations (3.8), 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 det D = 0 reduces to a cubic characteristic equation with respect to ω * : where b * , c * , d * are complex parameters (see Appendix C) defined by the original solution to system (3.6).
Obviously, if among the roots of the characteristic equation (3.13) exists such a root ω * = ω + iω * * that increment ω * * > 0, the amplitude of perturbation in (3.11) increases with time. This indicates the instability of the investigated solution (3.6). If at the same time then such a solution holds absolute instability, and the Cauchy problem of the system (3.2) is ill-posed (Godunov et al., 1976).
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.6) requires the study of the nonlinear behavior of the system (2.23).
4 Linear instability of original solutions 220

Instability of resting-state solution
We consider the regime when there is no motion and heat influx in the original solution, and the density is uniform horizontally: (4.1) The horizontal k hor and vertical wavenumbers k ver are defined With resting-state conditions (4.1), the cubic characteristic equation (3.13) is simplified as Here N is the dimensionless value of the Brent-Väisälä frequency N , C is sound speed. For scales (3.4) and standard atmosphere (Atmosphere, 1976) we have the following estimates: (4.4) The dependency of N 2 − G 2 on height for the standard atmosphere is presented in Figure 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.
For the cubic equation (3.13) despite the root ω (1) * = 0, expressions of the other two roots are shown in Appendix D when k x , k y , k z 1. Among these roots, there always exists one with a positive increment (ω * * > 0), which indicates the instability 235 of the resting-state solution.
The small Coriolis force can take place at near-equatorial latitudes (θ 1) and a height not close to 4 km. In such situation for k hor → ∞, From the general formulae of roots (Appendix D), we obtain the expressions of nonzero roots of (3.13) for small Coriolis force: When k hor → ∞ and k ver → ∞, there are various asymptotics, including If, in contrast to (4.5), the Coriolis force prevails: which can take place at high latitudes or altitudes near 4 km, then according to Appendix D we get As in (4.8) for k hor → ∞ and k ver → ∞, exist different asymptotics, including (4.11) From (4.8) and (4.11) it follows that the resting-state of the system of hydrodynamic equations with the vertical quasihydrostatic approximation (2.23) 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 absolute unstable. For perturbations when 255 the vertical wavelength is many times shorter than the horizontal wavelength (l ver l hor , k ver k hor 1), the resting-state tends to be neutral 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 /2k hor are shown in Figure 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 per-260 turbation 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 by (Arakawa and Konor, 2009). However, in the original system the following equation (thereinafter, the Arakawa inexact approx- (4.14) If instead of equation (2.19) or (4.13) the Arakawa inexact approximation (4.12) is adopted, then instead of (3.13) 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.7) and 280 (4.10). In the case of the predominance of Coriolis force, this difference is fundamental, namely, instead of absolute instability in (4.11), neutral stability follows from (4.17). In addition to the inexact equation (4.12), the conclusion about shortwave instability in (Arakawa and Konor, 2009) is made from the values of the real part of the root ω, not the imaginary part ω * * .
Such conclusion is inaccurate, as can be seen from (4.16).
In the case of small Coriolis force, the Arakawa inexact approximation gives absolute instability, with following asymptotics  The root with positive increment ω * * when f < Gk hor /k ver can lead to absolute instability (3.14). In particular, or small Coriolis force, The Holton inexact vertical quasi-hydrostatic approximation (Holton and Hakim, 2012) with equation of the constant local pressure (thereinafter, the Holton inexact approximation) can also be used: (4.23) Using this approximation instead of equation (2.19) or its equivalence (4.13), 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 (ω i.e., the density increases vertically, such original solution of resting-state is unstable for small Coriolis force f < N κ 4 , and analogously to (4.22), absolutely unstable.
One can also adopt the vertical quasi-hydrostatic approximation with the quasi-incompressibility dρ dt = 0, (4.28) from which together with the continuity equation, instead of (2.19) the vertical velocity v z is obtained by the following equation (4.29) For such a system , the characteristic equation at resting-state (4.1) has the same form (4.24), as for the Holton inexact approximation (4.23).
The vertical quasi-hydrostatic approximation with constant local density is proposed in (Marchuk, 1974)  The shortwave instability at resting-state is also studied by (Moore, 1985) for a two-dimensional system (v y = v y = 0) with approximation (4.29), 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.24).

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.23) is simplified as (4.34) its root is For a standard atmosphere, ∂ρ/∂z < 0, one always has ∂ρ/∂z M ∼ −10 −2 , (4.36) 350 from which it follows that ω * * ≈ − ∂v z ∂z . (4.37) 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 355 height (∂v z /∂z < 0), and such state is unstable (ω * * > 0).

Instability of solution with motion
Using the roots of the cubic equation (3.13), 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.13), we set the parameters of the original solution with motion as 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 C that the values k hor /k ver strongly change the coefficients of the cubic equation (3.13). For an increase of k ver , the growth of positive increment ω * * greatly weakens.
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 (4.39) 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 ).
5 Use of pseudo-viscosities to eliminate linear instability 375

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. The boundary conditions in the framework of the boundary value problem can significantly limit the growth of the perturbation. Nevertheless, the shortwave instability, of course, is one of the reasons of the possible instability of finite-difference calculations.

380
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 with turbulent viscosity in the horizontal momentum equations are relatively small in comparison with other terms.

385
However, the pseudo-viscosities, which are much greater than real turbulent 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.12) through where µ and λ are coefficients of shear and bulk pseudo-viscosities, respectively. We set the dimensionless values of these 390 coefficients growing with increasing wavenumbers: µ = c µ µ a k x k y k z n , λ = c λ µ a k x k y k z m µ a = µ a ρV ver L ver µ a = 1.8 × 10 −5 kg/(m · s) (5.3) Here µ a is the viscosity of the air.

Influence of pseudo-viscosities for resting-state solution
Similar as (3.13) the characteristic equation for a resting-state solution, taking into account the shear pseudo-viscosity in accordance with (5.2), 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.3)) 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 damps and disappears eventually.

Influence of pseudo-viscosities for the solution with motion
where b 1 (µ) and c 1 (µ) are the complex parameters determined by the original solution and shear pseudo-viscosity µ. 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 Figure 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 425 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.30) and (4.31) in comparison with the exact approximation.

435
It should be noted that for large values k j * * 1, there are large gradients of perturbation parameters (like impulse), and this leads an absolute instability (when k j * * → ∞) in response to relatively longer wavelength range (k hor < 30, for fixed k ver = 15). In addition, the amplitude of initial perturbation (5.9) increases for x, y, z < 0.

Conclusions
In this paper, an equation for vertical velocity refines and simplifies the equation obtained by Lorenz and Lorenz (1967).

440
Different from the approximations used in other works such as (Holton and Hakim, 2012;Marchuk, 1974), 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 quasi-hydrostatical approximation

460
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 relates 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 465 perturbation, so that numerical solutions to the asymptotically exact quasi-hydrostatic 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 470 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 quasi-hydrostatic equations of the present paper will be beneficial, in particular, for Appendix C: Complex Coefficients of the Cubic Equation (3.13): The following coefficients b * , c * and d * are expressed for k x , k y , k z 1.