Low-frequency Internal Gravity Waves Are Pseudo-incompressible

The anelastic condition is a relatively simple replacement for the continuity equation that captures significant density variation in the mean properties of the fluid. For instance, in a gravitationally stratified fluid with velocity, u , and time-averaged density that varies with height, ρ₀(z), the continuity equation is replaced with

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot \left({\rho }_{0}{\boldsymbol{u}}\right)=0.\end{eqnarray} \tag{ 1 }$$

This expression can be derived from the full continuity equation,

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\boldsymbol{\nabla }}\cdot \left(\rho {\boldsymbol{u}}\right)=0,\end{eqnarray} \tag{ 2 }$$

by making two assumptions that are often appropriate for flows of low Mach number: (1) the time derivative of the mass density ρ is inconsequential and (2) the fractional fluctuations of the density around the background density are small, i.e., ∣ρ₁/ρ₀∣ ≪ 1 where ρ = ρ₀ + ρ₁. The popularity of the anelastic approximation arises from two important properties. When the continuity equation is replaced by the anelastic condition, Equation (1), sound waves are removed as a permissible solution to the fluid equations and the mass flux ρ₀ u can be written using stream functions.

Brown et al. (2012) and Vasil et al. (2013) both remarked that, when the anelastic form of the continuity equation is employed, the fluid equations are no longer energy conserving without modifications to the momentum equation. To enforce conservation of energy, an otherwise unmotivated change to the buoyancy force is required. For an inviscid fluid, the vertical momentum equation can be written in the following form:

$$\begin{eqnarray}&&{\rho }_{0}\displaystyle \frac{{Dw}}{{Dt}}=-{\rho }_{0}\displaystyle \frac{d}{{dz}}\left(\displaystyle \frac{{P}_{1}}{{\rho }_{0}}\right)+\displaystyle \frac{g{\rho }_{0}}{{c}_{p}}{s}_{1}+\displaystyle \frac{{N}^{2}}{g}{P}_{1},\end{eqnarray} \tag{ 3 }$$

with the pressure P and specific entropy density s decomposed into a steady hydrostatic background and a fluctuation about that background, P = P₀ + P₁ and s = s₀ + s₁. The vertical velocity is w, c_p is the specific heat capacity at constant pressure, and z is the height within the atmosphere with concomitant unit vector $\hat{{\boldsymbol{z}}}$ antialigned with gravity, ${\boldsymbol{g}}=-g\hat{{\boldsymbol{z}}}$. Further, the quantity ${N}^{2}={{gc}}_{p}^{-1}{{ds}}_{0}/{dz}$ is the square of the atmosphere's buoyancy or Brunt–Väisälä frequency. In Equation (3), we have ignored the density fluctuation in the inertial term on the left-hand side, subtracted the steady hydrostatic component from the force balance, and used the ideal gas law to rewrite the density fluctuation in terms of the pressure and entropy fluctuations. To ensure energy conservation, the term involving the buoyancy frequency must be discarded or be physically subdominant. In a convection zone, where efficient heat transport drives the atmosphere toward an adiabatic gradient with N² ≈ 0, this approximation is completely justified and has been coined the Lantz–Braginsky–Roberts (LBR) formulation of the anelastic approximation (Lantz 1992; Braginsky & Roberts 1995). Conversely, in a stably stratified region, the term is not small and cannot generally be self-consistently ignored.

We will examine two distinct formulations of the anelastic approximation. Both replace the continuity equation with the anelastic condition (1). One of these approximations—which we will dub the "fiducial" anelastic approximation—will make no further assumptions, leaving the momentum equation unmodified. The other formulation will be the LBR anelastic approximation as discussed above, which ensures energy conservation by excising a specific term from the momentum equation.

The pseudo-incompressible approximation as proposed by Durran (1989) modifies the continuity equation under the assumption that Eulerian fluctuations of the gas pressure can be ignored. Following Durran (2008), we start by defining the potential density ρ_* for an ideal gas,

$$\begin{eqnarray}&&{\rho }_{* }\equiv \rho \,{e}^{s/{c}_{p}}.\end{eqnarray} \tag{ 4 }$$

If we take the convective derivative of the potential density and utilize the continuity Equation (2) and the thermal energy equation,

$$\begin{eqnarray}&&\rho T\displaystyle \frac{{Ds}}{{Dt}}=Q,\end{eqnarray} \tag{ 5 }$$

we obtain a prognostic equation for the potential density:

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\rho }_{* }}\left(\displaystyle \frac{\partial {\rho }_{* }}{\partial t}+{\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}{\rho }_{* }\right)=-{\boldsymbol{\nabla }}\cdot {\boldsymbol{u}}+\displaystyle \frac{Q}{{c}_{p}\rho T},\end{eqnarray} \tag{ 6 }$$

where T is the temperature and Q represents all irreversible thermodynamic processes, such as thermal diffusion, viscous heating, radiative transfer, etc. Finally, by invoking Equation (4) and the equation of state for an ideal gas,

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\rho }_{* }}\displaystyle \frac{\partial {\rho }_{* }}{\partial t}=\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial \rho }{\partial t}+\displaystyle \frac{1}{{c}_{p}}\displaystyle \frac{\partial s}{\partial t}=\displaystyle \frac{1}{\gamma P}\displaystyle \frac{\partial P}{\partial t},\end{eqnarray} \tag{ 7 }$$

we replace the time derivative of the potential density with the time derivative of the gas pressure:

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot \left({\rho }_{* }{\boldsymbol{u}}\right)=\displaystyle \frac{{\rho }_{* }}{\rho }\left(\displaystyle \frac{Q}{{c}_{p}T}-\displaystyle \frac{1}{{c}^{2}}\displaystyle \frac{\partial P}{\partial t}\right).\end{eqnarray} \tag{ 8 }$$

In the preceding equations, γ is the gas's adiabatic exponent and c is the sound speed given by c² = γ P/ρ.

Equation (8) is an exact form of the continuity equation for which no approximation has been made other than the gas being ideal. The pseudo-incompressible approximation is achieved by assuming that the term involving the time derivative of the gas pressure is negligible:

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot \left({\rho }_{* }{\boldsymbol{u}}\right)=\displaystyle \frac{{\rho }_{* }}{\rho }\displaystyle \frac{Q}{{c}_{p}T}.\end{eqnarray} \tag{ 9 }$$

Such an approximation is valid in the limit of infinite sound speed and is consistent with slow motions of low Mach number for which a displaced parcel of fluid rapidly reaches pressure equilibration with its new surroundings. Most importantly, making this approximation removes sound waves from the fluid equations in the same way that anelasticity does. Durran's form of the pseudo-incompressible approximation (Durran 2008) involves replacing the continuity equation by the preceding equation but otherwise leaving the other fluid equations unmodified—specifically, the momentum equation remains the same. For a derivation of the pseudo-incompressible approximation valid for a nonideal, multicomponent gas, we refer the reader to Almgren et al. (2006) and Nonaka et al. (2010). Further, we note that the pseudo-incompressible approximation can be derived using a Lagrangian formulation with a constraint that the Eulerian pressure fluctuation vanishes (Vasil et al. 2013).

For isentropic motion, the pseudo-incompressible condition reduces to a form that is reminiscent of the anelastic relation:

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot \left({\rho }_{* }{\boldsymbol{u}}\right)=0,\end{eqnarray} \tag{ 10 }$$

with the mass density replaced by the potential density. However, for flows with low Mach number, thermodynamic fluctuations are small and we can safely linearize Equation (10), replacing the potential density by the potential density of the hydrostatic background atmosphere (denoted by "0" subscripts):

$$\begin{eqnarray}&&{\rho }_{* 0}\approx {\rho }_{0}{e}^{{s}_{0}/{c}_{p}}=\left(\displaystyle \frac{\hat{\rho }}{{\hat{P}}^{1/\gamma }}\right)\,{P}_{0}^{1/\gamma }.\end{eqnarray} \tag{ 11 }$$

The last equivalency in Equation (11) arises by noting that the potential density is the density that a fluid parcel would possess if displaced adiabatically to a fiducial height in the atmosphere where ${P}_{0}=\hat{P}$, ${\rho }_{0}=\hat{\rho }$, and s₀ = 0. Like the anelastic approximation, the flow field can be expressed using stream functions when Equations (10) and (11) are valid:

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot \left({P}_{0}^{1/\gamma }{\boldsymbol{u}}\right)=0.\end{eqnarray} \tag{ 12 }$$

We remind the reader that these two equations were derived using two assumptions: (1) the advective timescales are fast compared to diffusion times, i.e., isentropic motion; and (2) thermodynamic fluctuations are small compared to the background atmosphere.

Here, we demonstrate that any viable soundproofing technique must produce an appropriate coefficient for the first-derivative term that appears in Equation (23). This term is crucial for energy conservation. To see this, consider the vertical energy flux for an acoustic-gravity wave, $F(z)=\left\langle w\,{P}_{1}\right\rangle$, where angular brackets 〈〉 indicate a temporal average over a wave period. Because the second term on the right-hand side of Equation (21) is 90° out of phase with the vertical velocity, in a time average the second term's contribution vanishes and the energy flux can be written just in terms of the Lagrangian pressure fluctuation,

$$\begin{eqnarray}&&F(z)=\left\langle w\delta P\right\rangle =\displaystyle \frac{1}{4}\left(w\,\delta {P}^{* }+{w}^{* }\,\delta P\right),\end{eqnarray} \tag{ 27 }$$

where the superscript asterisks denote complex conjugation. By employing Equation (19), one can demonstrate that this flux is inversely proportional to the mass density and proportional to the Wronskian of the Lagrangian pressure fluctuation and its complex conjugate,

$$\begin{eqnarray}&&F(z)=-\displaystyle \frac{i{\omega }^{3}}{4\alpha {\rho }_{0}}\left(\delta P\displaystyle \frac{d\,\delta {P}^{* }}{{dz}}-\delta {P}^{* }\displaystyle \frac{d\,\delta P}{{dz}}\right).\end{eqnarray} \tag{ 28 }$$

Abel's Identity tells us that, to within an unknown multiplicative constant, C, the Wronskian depends only on the coefficient of the first-derivative term in the ODE. For the ODE here, the necessary integration is trivial to perform:

$$\begin{eqnarray}&&{ \mathcal W }\left\{\delta P,\delta {P}^{* }\right\}(z)=C\,\exp \left(-\int \displaystyle \frac{{dz}}{H}\right)=C\,{\rho }_{0}.\end{eqnarray} \tag{ 29 }$$

Hence, the energy flux is constant with height even though the coefficients of the ODE are vertically variable,

$$\begin{eqnarray}&&F(z)=-\displaystyle \frac{i{\omega }^{3}C}{4\alpha }=\mathrm{constant}.\end{eqnarray} \tag{ 30 }$$

The constancy of the energy flux with height in the atmosphere is one way to characterize the conservation of energy by acoustic-gravity waves.

From this analysis, we can deduce that any approximation that incorrectly reproduces the first-derivative term, may produce wave solutions with energy fluxes that vary with height. Consequently, such approximations will fail to conserve energy. For example, if the first-derivative term is artificially set to zero, the flux will be inversely proportional to the mass density and F(z) will spuriously increase with height. This is the fundamental reason why Brown et al. (2012) and Vasil et al. (2013) found a lack of energy conservation when applying a variety of anelastic approximations to an isothermal atmosphere. Those approximations failed to correctly reproduce the first-derivative term of the ODE. Here, we show that it is a general property for any stratification, not just an isothermal one.

For a general stratification, the coefficients of the ODE (23) are functions of height and the solutions will not be sinusoidal. However, by making a change of variable that converts the ODE into standard form (i.e., a Helmholtz equation that lacks a first-derivative term), a local dispersion relation can be generated which is appropriate in a WKB framework (e.g., Bender & Orszag 1999). The required change of variable involves the square root of the mass density, $\delta P={\left(\alpha {\rho }_{0}\right)}^{1/2}\psi$. We include the constant α inside the square root purely for the sake of symmetry in later sections when we explore various soundproofing techniques. Here, its inclusion is unnecessary and only introduces a multiplicative constant that factors out of the resulting ODE:

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}\psi }{{{dz}}^{2}}+{k}_{z}^{2}\psi =0,\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{k}_{z}^{2}(z)=\displaystyle \frac{{\omega }^{2}-{\omega }_{c}^{2}}{{c}^{2}}-{k}_{h}^{2}\left(1-\displaystyle \frac{{N}^{2}}{{\omega }^{2}}\right).\end{eqnarray} \tag{ 32 }$$

In the preceding equations, k_z (z) is a local vertical wavenumber and ω_c (z) is the acoustic-cutoff frequency, which depends on the stratification through the density scale height H:

$$\begin{eqnarray}&&\displaystyle \frac{{\omega }_{c}^{2}}{{c}^{2}}\equiv \displaystyle \frac{1-2{H}^{{\prime} }}{4{H}^{2}}.\end{eqnarray} \tag{ 33 }$$

We denote vertical derivatives of atmospheric quantities using a superscript prime, i.e., the vertical derivative of the density scale height is given by ${H}^{{\prime} }\equiv {dH}/{dz}$.

From the preceding analysis, we see that acoustic-gravity waves vary over two relevant vertical spatial scales: a local vertical wavelength and an envelope scale. The wavelength is given by the local dispersion relation (32) and hence depends on the wave frequency as well as the characteristic frequencies of the atmosphere—i.e., the buoyancy frequency N, the acoustic-cutoff frequency ω_c , and the Lamb frequency k_h c. The envelope scale is associated with vertical variation of the envelope function ${\left(\alpha {\rho }_{0}\right)}^{1/2}$ that appears in the change of variable above. This function provides a local amplitude of the wave function (in a WKB sense). Because the envelope function only depends on the mass density, the envelope scale is solely determined by the atmospheric stratification through the density scale height H. For later convenience, we choose to define the envelope scale Λ as twice the scale length associated with the envelope function, such that Λ = H:

$$\begin{eqnarray}&&{{\rm{\Lambda }}}^{-1}\equiv -\displaystyle \frac{2}{{\left(\alpha {\rho }_{0}\right)}^{1/2}}\displaystyle \frac{d{\left(\alpha {\rho }_{0}\right)}^{1/2}}{{dz}}={H}^{-1}.\end{eqnarray} \tag{ 34 }$$

The nondimensional forms of the other fluid variables can be generated through Equations (18), (19), and (26), and they are related to the Lagrangian pressure fluctuation through differential operators:

$$\begin{eqnarray}&&\tilde{u}={\left(\displaystyle \frac{{k}_{h}}{g}\right)}^{1/2}u=-\displaystyle \frac{{\tilde{\rho }}_{0}^{-1}}{1-{\varepsilon }^{4}}\left(\displaystyle \frac{d}{d\zeta }+{\varepsilon }^{3}\right)\delta \tilde{P}\sim { \mathcal O }(1),\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&\tilde{w}={\left(\displaystyle \frac{{k}_{h}}{g}\right)}^{1/2}w=\displaystyle \frac{i\varepsilon \,{\tilde{\rho }}_{0}^{-1}}{1-{\varepsilon }^{4}}\left(1+\varepsilon \displaystyle \frac{d}{d\zeta }\right)\delta \tilde{P}\sim { \mathcal O }(\varepsilon ),\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}&&{\tilde{P}}_{1}=\left(\displaystyle \frac{{k}_{h}}{g\hat{\rho }}\right){P}_{1}=-\displaystyle \frac{\varepsilon }{1-{\varepsilon }^{4}}\left(\displaystyle \frac{d}{d\zeta }+{\varepsilon }^{3}\right)\delta \tilde{P}\sim { \mathcal O }(\varepsilon ),\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&{\tilde{s}}_{1}=\displaystyle \frac{{s}_{1}}{{c}_{p}}=\displaystyle \frac{{\tilde{N}}^{2}{\tilde{\rho }}_{0}^{-1}}{1-{\varepsilon }^{4}}\left(1+\varepsilon \displaystyle \frac{d}{d\zeta }\right)\delta \tilde{P}\sim { \mathcal O }(1),\end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray}&&{\tilde{\rho }}_{1}=\displaystyle \frac{{\rho }_{1}}{\hat{\rho }}=-\displaystyle \frac{1}{1-{\varepsilon }^{4}}\left({\tilde{N}}^{2}+\displaystyle \frac{\varepsilon }{\tilde{H}}\displaystyle \frac{d}{d\zeta }+\displaystyle \frac{{\varepsilon }^{4}}{{\tilde{c}}^{2}}\right)\delta \tilde{P}\,\,\sim { \mathcal O }(1),\end{eqnarray} \tag{ 47 }$$

where ${\tilde{\rho }}_{0}={\rho }_{0}/\hat{\rho }$ is the nondimensional atmospheric density.

We can immediately see that internal gravity waves possess motions that are nearly horizontal for low frequencies. The vertical velocity component w is small by a factor of ε. Furthermore, while the Lagrangian pressure fluctuation remains of order unity in size, $\delta P\sim { \mathcal O }(1)$, the Eulerian pressure fluctuation becomes small, ${P}_{1}\sim { \mathcal O }(\varepsilon$). Both the entropy and density fluctuations remain of order unity. The fact that the Eulerian pressure fluctuation vanishes in the limit of low frequency is consistent with the pseudo-incompressible approximation and ensures that the internal gravity waves and acoustic waves decouple in that limit. However, because the mass density fluctuation does not vanish, these limits further suggest that this decoupling is not accomplished through the anelastic limit. We explore this result fully in the next subsection.

In order to make obvious the relative magnitude of terms in subsequent equations, we define alternate dimensionless variables for the vertical velocity and Eulerian pressure fluctuation:

$$\begin{eqnarray}&&\tilde{w}\equiv \varepsilon \tilde{W},\qquad {\tilde{P}}_{1}\equiv \varepsilon \tilde{{\rm{\Theta }}}.\end{eqnarray} \tag{ 48 }$$

Both $\tilde{W}$ and $\tilde{{\rm{\Theta }}}$ are of order unity because the prefactors in their definitions absorb the leading-order behavior as the frequency becomes small.

Consider the dimensional form of the continuity Equation (15), where the equation of state (16) is used to replace the density fluctuation:

$$\begin{eqnarray}&&i\omega \displaystyle \frac{{\rho }_{1}}{{\rho }_{0}}=i\omega \left(\displaystyle \frac{{P}_{1}}{{\rho }_{0}{c}^{2}}-\displaystyle \frac{{s}_{1}}{{c}_{p}}\right)={{ik}}_{h}u+\displaystyle \frac{{dw}}{{dz}}-\displaystyle \frac{w}{H}.\end{eqnarray} \tag{ 49 }$$

In order to soundproof the equation set, we need to eliminate the term involving the Eulerian pressure fluctuation. This term is responsible for producing the pressure fluctuations that generate the restoring force for acoustic oscillations.

The anelastic approximation does indeed eliminate this pressure term, but it is overkill and removes the entire left-hand side of the continuity equation above. In particular, the term involving the entropy fluctuation is also thrown away. For low-frequency internal gravity waves, this is inconsistent. If the continuity equation is nondimensionalized, it becomes obvious that the entropy term is of the same order as other terms that are retained by the anelastic approximation:

$$\begin{eqnarray}&&i{\varepsilon }^{2}\displaystyle \frac{\tilde{{\rm{\Theta }}}}{{\tilde{\rho }}_{0}{\tilde{c}}^{2}}-i\varepsilon {\tilde{s}}_{1}=\left[i\tilde{u}+\displaystyle \frac{d\tilde{W}}{d\zeta }\right]-\varepsilon \displaystyle \frac{\tilde{W}}{\tilde{H}}.\end{eqnarray} \tag{ 50 }$$

The leading-order behavior consists of the two order-unity terms that appear in square brackets on the right-hand side of Equation (50). The first correction for nonzero frequency is comprised of the two first-order terms, ${ \mathcal O }(\varepsilon );$ one of these is the aforesaid entropy term. The term involving the Eulerian pressure fluctuation is second order, ${ \mathcal O }({\varepsilon }^{2})$.

The lowest-order self-consistent approximation that one could make would be to keep just the leading-order terms, resulting in an assumption of incompressibility, ∇ · u ≈ 0. The next self-consistent approximation would be the retention of all zero-order and first-order terms. As we will show next, this approximation is equivalent to the pseudo-incompressible condition.

We demonstrate pseudo-incompressibility by using the energy Equation (14) to replace the entropy fluctuation in Equation (49) with the vertical velocity and then combining the two first-order terms using the definition of the buoyancy frequency, N² = g/H − g²/c²:

$$\begin{eqnarray}&&i\omega \displaystyle \frac{{P}_{1}}{{\rho }_{0}{c}^{2}}=\left[{{ik}}_{h}u+\displaystyle \frac{{dw}}{{dz}}\right]-\displaystyle \frac{{gw}}{{c}^{2}}.\end{eqnarray} \tag{ 51 }$$

The last term on the right-hand side is equal to the vertical velocity divided by the scale height for the potential density, i.e., the density scale height for an adiabatic density stratification:

$$\begin{eqnarray}&&\displaystyle \frac{1}{{H}_{* }}\equiv -\displaystyle \frac{1}{{\rho }_{* 0}}\displaystyle \frac{d{\rho }_{* 0}}{{dz}}=\displaystyle \frac{g}{{c}^{2}}.\end{eqnarray} \tag{ 52 }$$

Hence, the terms on the right-hand side of Equation (51) can be cleanly combined:

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \left({\rho }_{* 0}{\boldsymbol{u}}\right)=i\omega \displaystyle \frac{{\rho }_{* 0}{P}_{1}}{{\rho }_{0}{c}^{2}}\sim { \mathcal O }({\varepsilon }^{2}).\end{eqnarray} \tag{ 53 }$$

A self-consistent low-frequency approximation is to discard all second-order terms, leading to the pseudo-incompressible approximation, ${\boldsymbol{\nabla }}\cdot \left({\rho }_{* 0}{\boldsymbol{u}}\right)=0$.

When transformed into the spectral representation, the vertical component of the momentum Equation (13) is given by

$$\begin{eqnarray}&&-i\omega {\rho }_{0}w=-\displaystyle \frac{{{dP}}_{1}}{{dz}}-g{\rho }_{1},\end{eqnarray} \tag{ 54 }$$

and nondimensionalization of this equation yields

$$\begin{eqnarray}&&-i{\varepsilon }^{2}{\tilde{\rho }}_{0}\tilde{W}=-\displaystyle \frac{d\tilde{{\rm{\Theta }}}}{d\zeta }-{\tilde{\rho }}_{1}.\end{eqnarray} \tag{ 55 }$$

It is now obvious from the preceding equation that the inertial term on the left-hand side becomes the smallest term in the low-frequency limit; it is a second-order correction. The right-hand side consists solely of terms that are zero-order in the dimensionless frequency. Hence, to first order, the balance is simply the hydrostatic relation between the perturbed pressure and the perturbed density:

$$\begin{eqnarray}&&-\displaystyle \frac{{{dP}}_{1}}{{dz}}-g{\rho }_{1}\approx 0.\end{eqnarray} \tag{ 56 }$$

The pseudo-incompressible and fiducial anelastic approximations both leave the vertical momentum equation unmodified, but the LBR formulation of the anelastic approximation drops a term whose removal is formally valid only in an adiabatic (or near-adiabatic) stratification. The vertical momentum equation (54) can be rewritten in the following manner:

$$\begin{eqnarray}&&-i\omega w=-\displaystyle \frac{d}{{dz}}\left(\displaystyle \frac{{P}_{1}}{{\rho }_{0}}\right)+\displaystyle \frac{{{gs}}_{1}}{{c}_{p}}+\displaystyle \frac{{N}^{2}}{g}\displaystyle \frac{{P}_{1}}{{\rho }_{0}},\end{eqnarray} \tag{ 57 }$$

either by linearizing Equation (3) or by dividing the vertical momentum Equation (54) by the mass density and pulling the density into the gradient operator that appears in the pressure force by use of the chain rule. The LBR formulation of the anelastic approximation removes the term involving the buoyancy frequency, even in stable stratifications where the buoyancy frequency is not small. We demonstrate that this approximation is inconsistent with low-frequency gravity waves by considering the nondimensional form of Equation (57):

$$\begin{eqnarray}&&-i{\varepsilon }^{2}\tilde{W}=-\displaystyle \frac{d}{d\zeta }\left(\displaystyle \frac{\tilde{{\rm{\Theta }}}}{{\tilde{\rho }}_{0}}\right)+{\tilde{s}}_{1}+\varepsilon {\tilde{N}}^{2}\displaystyle \frac{\tilde{{\rm{\Theta }}}}{{\tilde{\rho }}_{0}}.\end{eqnarray} \tag{ 58 }$$

The LBR approximation inconsistently ignores the first-order term while retaining the inertial term (which is second-order).

Because the pseudo-incompressible approximation is self-consistent in its treatment of the continuity equation and correct to first order in the frequency, we expect that this approximation should produce low-frequency internal gravity waves that are correct to first order in the dispersion relation and in the wave functions. To demonstrate that this expectation is true, we rederive the wave equation for internal gravity waves but with the continuity Equation (15) replaced by the pseudo-incompressible condition, ${\boldsymbol{\nabla }}\cdot \left({\rho }_{* 0}{\boldsymbol{u}}\right)=0$. We simply present the result:

$$\begin{eqnarray}&&\begin{array}{l}\left\{\displaystyle \frac{{d}^{2}}{{{dz}}^{2}}+\left(\displaystyle \frac{1}{H}+\displaystyle \frac{{\omega }^{2}}{g}{\theta }_{\mathrm{PI}}\right)\displaystyle \frac{d}{{dz}}-{k}_{h}^{2}\left(1-\displaystyle \frac{{N}^{2}}{{\omega }^{2}}\right)\right.\\ \qquad \qquad \qquad \qquad \left.+\left[\displaystyle \frac{{N}^{2}}{{c}^{2}}+\displaystyle \frac{{\omega }^{4}}{{g}^{2}}{\theta }_{\mathrm{PI}}\right]\right\}\delta P=0.\end{array}\end{eqnarray} \tag{ 61 }$$

In this expression, θ_PI(z) is a dimensionless function that depends on the temporal frequency ω, the horizontal wavenumber k_h , and the potential density ρ_*0 through the following definitions:

$$\begin{eqnarray}&&{\alpha }_{\mathrm{PI}}(z)\equiv {g}^{2}{k}_{h}^{2}-{\omega }^{4}+{\omega }^{2}\displaystyle \frac{g}{{H}_{* }},\end{eqnarray} \tag{ 62 }$$

$$\begin{eqnarray}&&{\theta }_{\mathrm{PI}}(z)\equiv -\displaystyle \frac{g}{{\omega }^{2}}\displaystyle \frac{{\alpha }_{\mathrm{PI}}^{{\prime} }}{{\alpha }_{\mathrm{PI}}}=\displaystyle \frac{{g}^{2}}{{\alpha }_{\mathrm{PI}}}\,\displaystyle \frac{{H}_{* }^{{\prime} }}{{H}_{* }^{2}}.\end{eqnarray} \tag{ 63 }$$

Compared to the fully compressible equations, the quantity α_PI has been augmented by ω² g/H_* and is therefore no longer a constant function of height.

A direct comparison of Equation (61) with the wave equation for a fully compressible fluid (23) reveals that there are three spurious terms: both of the terms involving θ_PI, as well as the term (N²/c²)δ P. To demonstrate that all of these spurious terms are small in magnitude and can be safely ignored in the low-frequency limit, we nondimensionalize Equation (61),

$$\begin{eqnarray}&&\begin{array}{l}\left\{\displaystyle \frac{{d}^{2}}{d{\zeta }^{2}}+\left(\displaystyle \frac{\varepsilon }{\tilde{H}}+{\varepsilon }^{3}{\theta }_{\mathrm{PI}}\right)\displaystyle \frac{d}{d\zeta }+{\tilde{N}}^{2}-{\varepsilon }^{2}\right.\\ \qquad \qquad \qquad \qquad +\left.\left[{\varepsilon }^{2}\displaystyle \frac{{\tilde{N}}^{2}}{{\tilde{c}}^{2}}+{\varepsilon }^{6}{\theta }_{\mathrm{PI}}\right]\right\}\delta \tilde{P}=0,\end{array}\end{eqnarray} \tag{ 64 }$$

and we recognize that the function θ_PI is an order-unity quantity for low frequencies:

$$\begin{eqnarray}&&{\theta }_{\mathrm{PI}}=\displaystyle \frac{1}{1-{\varepsilon }^{4}+{\varepsilon }^{2}{\tilde{H}}_{* }^{-1}}\,\displaystyle \frac{{H}_{* }^{{\prime} }}{{\tilde{H}}_{* }^{2}}\quad \sim { \mathcal O }(1).\end{eqnarray} \tag{ 65 }$$

Thus, all of the spurious terms are second-order or higher in the dimensionless frequency ε and the Lagrangian pressure fluctuation that is generated by Equation (61) is correct to first order.

Based on this result, we should expect the local dispersion relation to also be correct to first order, and this is indeed the case. The transformation that converts the ODE into a Helmholtz equation has the same form as we found for the fully compressible equations,

$$\begin{eqnarray}&&\delta P={\left({\alpha }_{\mathrm{PI}}\,{\rho }_{0}\right)}^{1/2}\psi ,\end{eqnarray} \tag{ 66 }$$

but now the function α = α_PI(z) varies with height. This change of variable leads to the following local dispersion relation,

$$\begin{eqnarray}&&\begin{array}{l}{k}_{z}^{2}(z)={k}_{h}^{2}\left(\displaystyle \frac{{N}^{2}}{{\omega }^{2}}-1\right)-\displaystyle \frac{{\omega }_{c}^{2}}{{c}^{2}}\\ \,+\,\left[\displaystyle \frac{{N}^{2}}{{c}^{2}}-\displaystyle \frac{{\omega }^{2}}{2g}\left({\theta }_{\mathrm{PI}}^{{\prime} }+\displaystyle \frac{{\theta }_{\mathrm{PI}}}{H}\right)+\displaystyle \frac{{\omega }^{4}}{{g}^{2}}\left({\theta }_{\mathrm{PI}}-\displaystyle \frac{{\theta }_{\mathrm{PI}}^{2}}{4}\right)\right],\end{array}\end{eqnarray} \tag{ 67 }$$

with a nondimensional form given by

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{k}_{z}^{2}(z)}{{k}_{h}^{2}}={\varepsilon }^{-2}{\tilde{N}}^{2}-\left(1+{\tilde{k}}_{c}^{2}\right)\\ +\,\left[\displaystyle \frac{{\tilde{N}}^{2}}{{\tilde{c}}^{2}}-\displaystyle \frac{{\varepsilon }^{2}}{2}\left(\displaystyle \frac{{\theta }_{\mathrm{PI}}^{{\prime} }}{{k}_{h}}+\displaystyle \frac{{\theta }_{\mathrm{PI}}}{\tilde{H}}\right)+{\varepsilon }^{4}\left({\theta }_{\mathrm{PI}}-\displaystyle \frac{{\theta }_{\mathrm{PI}}^{2}}{4}\right)\right].\end{array}\end{eqnarray} \tag{ 68 }$$

All of the terms contained by the error term, E_PI, in the preceding equations are spurious and do not appear in the local dispersion relation for a fully compressible fluid. However, all spurious terms appear as a correction that is smaller than the leading-order behavior by a factor of ε² or smaller. Hence, the pseudo-incompressible approximation leads to a local dispersion relation that is correct to first order.

Finally, the envelope scale can be read directly from the coefficient of the first-derivative term in the ODE, Λ⁻¹ = H⁻¹ + ω² θ_PI/g. To first order in the frequency, the envelope scale is simply the density scale height.

For the fiducial anelastic approximation, where the only modification to the fully compressible fluid equations is made to the continuity equation, the resulting ODE for the Lagrangian pressure fluctuation is as follows:

$$\begin{eqnarray}&&\begin{array}{l}\left\{\displaystyle \frac{{d}^{2}}{{{dz}}^{2}}+\left(\displaystyle \frac{1}{{H}_{* }}+\displaystyle \frac{{\omega }^{2}}{g}{\theta }_{\mathrm{FA}}\right)\displaystyle \frac{d}{{dz}}-{k}_{h}^{2}\left(1-\displaystyle \frac{{N}^{2}}{{\omega }^{2}}\right)\right.\\ \qquad \qquad \ \left.+\,\left[\left(\displaystyle \frac{{\omega }^{2}}{{c}^{2}}+{k}_{h}^{2}\right){\theta }_{\mathrm{FA}}-\displaystyle \frac{{H}_{* }^{{\prime} }}{{H}_{* }^{2}}\right]\right\}\delta P=0,\end{array}\end{eqnarray} \tag{ 69 }$$

where the α and θ functions take on subtly but crucially different forms,

$$\begin{eqnarray}&&{\alpha }_{\mathrm{FA}}(z)\equiv {g}^{2}{k}_{h}^{2}-{\omega }^{4}+{\omega }^{2}\displaystyle \frac{g}{H},\end{eqnarray} \tag{ 70 }$$

$$\begin{eqnarray}&&{\theta }_{\mathrm{FA}}(z)\equiv -\displaystyle \frac{g}{{\omega }^{2}}\displaystyle \frac{{\alpha }_{\mathrm{FA}}^{{\prime} }}{{\alpha }_{\mathrm{FA}}}=\displaystyle \frac{{g}^{2}}{{\alpha }_{\mathrm{FA}}}\,\displaystyle \frac{{H}^{{\prime} }}{{H}^{2}}.\end{eqnarray} \tag{ 71 }$$

Here, α_FA and θ_FA differ from the pseudo-incompressible case, Equations (62) and (63), by the appearance of H instead of H_*.

A direct comparison of Equation (69) with the ODE (23) appropriate for a fully compressible fluid illustrates that fiducial anelastic generates a variety of spurious and incorrect terms. Specifically, the terms in the square brackets are spurious and the entire coefficient of the first-derivative term is incorrect. To ascertain the magnitude of these mistakes, we nondimensionalize:

$$\begin{eqnarray}&&\begin{array}{l}\left\{\displaystyle \frac{{d}^{2}}{d{\zeta }^{2}}+\left(\displaystyle \frac{\varepsilon }{{\tilde{H}}_{* }}+{\varepsilon }^{3}{\theta }_{\mathrm{FA}}\right)\displaystyle \frac{d}{d\zeta }+{\tilde{N}}^{2}-{\varepsilon }^{2}\right.\\ \qquad \qquad +\,\left.{\varepsilon }^{2}\left({\theta }_{\mathrm{FA}}-\displaystyle \frac{{H}_{* }^{{\prime} }}{{\tilde{H}}_{* }^{2}}\right)+{\varepsilon }^{4}\displaystyle \frac{{\theta }_{\mathrm{FA}}}{{\tilde{c}}^{2}}\right\}\delta \tilde{P}=0.\end{array}\end{eqnarray} \tag{ 72 }$$

Fiducial anelastic performs rather poorly in reproducing the behavior of low-frequency internal gravity waves. The ODE is correct only to leading order in ε, with inconsistencies appearing at first order in the coefficient of the first derivative. The first term in this coefficient contains the reciprocal of the scale height of the potential density, where it should instead possess the reciprocal of the density scale height—see Equation (38).

Interestingly, conversion of the ODE to standard form—via the change of variable $\delta P={\left({\alpha }_{\mathrm{FA}}\,{\rho }_{* }\right)}^{1/2}\psi$—results in a local dispersion relation that is correct to first order:

$$\begin{eqnarray}\begin{array}{rcl}{k}_{z}^{2}(z) & = & {k}_{h}^{2}\left(\displaystyle \frac{{N}^{2}}{{\omega }^{2}}-1\right)-\displaystyle \frac{1+2{H}_{* }^{{\prime} }}{4{H}_{* }^{2}}\\ & & +\,\left[{k}_{h}^{2}{\theta }_{\mathrm{FA}}-\displaystyle \frac{{\omega }^{2}}{2g}\left({\theta }_{\mathrm{FA}}^{{\prime} }-\displaystyle \frac{{\theta }_{\mathrm{FA}}}{{H}_{* }}\right)-\displaystyle \frac{{\omega }^{4}}{4{g}^{2}}{\theta }_{\mathrm{FA}}^{2}\right],\end{array}\end{eqnarray} \tag{ 73 }$$

or

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{k}_{z}^{2}(z)}{{k}_{h}^{2}} & = & {\varepsilon }^{-2}{\tilde{N}}^{2}-\left(1+\displaystyle \frac{1+2{H}_{* }^{{\prime} }}{4{\tilde{H}}_{* }^{2}}\right)\\ & & +\,\left[{\theta }_{\mathrm{FA}}+\displaystyle \frac{{\varepsilon }^{2}}{2}\left({k}_{h}^{-1}{\theta }_{\mathrm{FA}}^{{\prime} }-\displaystyle \frac{{\theta }_{\mathrm{FA}}}{{\tilde{H}}_{* }}\right)-\displaystyle \frac{{\varepsilon }^{4}}{4}{\theta }_{\mathrm{FA}}^{2}\right].\end{array}\end{eqnarray} \tag{ 74 }$$

In addition to all of the spurious terms that appear in the square brackets, the acoustic-cutoff frequency is incorrect:

$$\begin{eqnarray}&&\displaystyle \frac{1+2{H}_{* }^{{\prime} }}{4{\tilde{H}}_{* }^{2}}\ne {\tilde{k}}_{c}^{2}=\displaystyle \frac{1-2{H}^{{\prime} }}{4{\tilde{H}}^{2}}.\end{eqnarray} \tag{ 75 }$$

For ease of comparison, in Table 1 we have reworked the right-hand side of Equation (73) to extract the correct form of the acoustic-cutoff frequency. Despite these issues, the errors all appear at second order or higher in the dimensionless frequency ε, meaning that the erroneous terms divided by the leading-order behavior are small by a factor of ε². The fact that the ODE itself is incorrect at first order manifests in the envelope function, ${\left({\alpha }_{\mathrm{FA}}\,{\rho }_{* }\right)}^{1/2}$, which is wrong at all orders. As we will see in a subsequent section this results in first-order errors in the wave functions even though the dispersion relation is correct to first order.

In the framework of the LBR anelastic approximation, in addition to the anelastic treatment of the continuity equation, i.e., ${\boldsymbol{\nabla }}\cdot \left({\rho }_{0}{\boldsymbol{u}}\right)\approx 0$, a term in the vertical momentum equation is removed. When these two modifications to the fluid equations are adopted, the resulting ODE that describes internal gravity waves becomes

$$\begin{eqnarray}&&\begin{array}{l}\left\{\displaystyle \frac{{d}^{2}}{{{dz}}^{2}}+\left(\displaystyle \frac{1}{H}+\displaystyle \frac{{\omega }^{2}}{g}{\theta }_{\mathrm{LBR}}\right)\displaystyle \frac{d}{{dz}}-{k}_{h}^{2}\left(1-\displaystyle \frac{{N}^{2}}{{\omega }^{2}}\right)\right.\\ +\,\left.\left[\left({k}_{h}^{2}+\displaystyle \frac{{\omega }^{2}}{{gH}}\right){\theta }_{\mathrm{LBR}}-\displaystyle \frac{{H}^{{\prime} }}{{H}^{2}}\right]\right\}\delta P=0,\end{array}\end{eqnarray} \tag{ 76 }$$

where α and θ are now

$$\begin{eqnarray}&&{\alpha }_{\mathrm{LBR}}(z)\equiv {g}^{2}{k}_{h}^{2}-{\omega }^{4}+{\omega }^{2}\left({N}^{2}+\displaystyle \frac{g}{H}\right),\end{eqnarray} \tag{ 77 }$$

$$\begin{eqnarray}&&{\theta }_{\mathrm{LBR}}(z)\equiv \displaystyle \frac{{g}^{2}}{{\alpha }_{\mathrm{LBR}}}\left(\displaystyle \frac{{H}^{{\prime} }}{{H}^{2}}-\displaystyle \frac{1}{g}\displaystyle \frac{{{dN}}^{2}}{{dz}}\right).\end{eqnarray} \tag{ 78 }$$

The nondimensional form of the ODE becomes

$$\begin{eqnarray}&&\begin{array}{l}\left\{\displaystyle \frac{{d}^{2}}{d{\zeta }^{2}}+\left(\displaystyle \frac{\varepsilon }{\tilde{H}}+{\varepsilon }^{3}{\theta }_{\mathrm{LBR}}\right)\displaystyle \frac{d}{d\zeta }+{\tilde{N}}^{2}-{\varepsilon }^{2}\right.\\ +\,\left.\left[{\varepsilon }^{2}\left({\theta }_{\mathrm{LBR}}-\displaystyle \frac{{H}^{{\prime} }}{{\tilde{H}}^{2}}\right)+{\varepsilon }^{4}\displaystyle \frac{{\theta }_{\mathrm{LBR}}}{\tilde{H}}\right]\right\}\delta \tilde{P}=0.\end{array}\end{eqnarray} \tag{ 79 }$$

Despite the inconsistent treatment of the vertical momentum equation, the LBR form of the anelastic approximation generates an ODE that is correct to first order in ε. The spurious terms that appear in the square brackets are second-order or higher and the coefficient of the first derivative is correct to first order. As expected, the local dispersion relation—once again achieved by the change of variable $\delta P={\left({\alpha }_{\mathrm{LBR}}\,{\rho }_{0}\right)}^{1/2}\psi$—is correct to first order:

$$\begin{eqnarray}\begin{array}{rcl}{k}_{z}^{2}(z) & = & {k}_{h}^{2}\left(\displaystyle \frac{{N}^{2}}{{\omega }^{2}}-1\right)-\displaystyle \frac{{\omega }_{c}^{2}}{{c}^{2}}-\displaystyle \frac{{H}^{{\prime} }}{{H}^{2}}\\ & & +\,{k}_{h}^{2}{\theta }_{\mathrm{LBR}}-\displaystyle \frac{{\omega }^{2}}{2g}\left({\theta }_{\mathrm{LBR}}^{{\prime} }-\displaystyle \frac{{\theta }_{\mathrm{LBR}}}{H}\right)-\displaystyle \frac{{\omega }^{4}}{4{g}^{2}}{\theta }_{\mathrm{LBR}}^{2},\end{array}\end{eqnarray} \tag{ 80 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{k}_{z}^{2}(z)}{{k}_{h}^{2}} & = & {\varepsilon }^{-2}{\tilde{N}}^{2}-\left(1+{\tilde{k}}_{c}^{2}+\displaystyle \frac{{H}^{{\prime} }}{{\tilde{H}}^{2}}\right)\\ & & +\,{\theta }_{\mathrm{LBR}}-\displaystyle \frac{{\varepsilon }^{2}}{2}\left({k}_{h}^{-1}{\theta }_{\mathrm{LBR}}^{{\prime} }-\displaystyle \frac{{\theta }_{\mathrm{LBR}}}{\tilde{H}}\right)-\displaystyle \frac{{\varepsilon }^{4}}{4}{\theta }_{\mathrm{LBR}}^{2}.\end{array}\end{eqnarray} \tag{ 81 }$$

In the previous subsections, we demonstrated that the three approximations generate errors to the vertical wavelength of internal gravity waves that are of second order in the dimensionless frequency ε. Hence, if the only test of fidelity was to reproduce the local dispersion relation, all of the soundproofing treatments would fare equally well. This is born out by a comparison of the vertical wavenumber that is achieved in an isothermal atmosphere by each treatment. This type of atmosphere is one of the most lenient of all potential atmospheres, because all of the characteristic frequencies, i.e., N, ω_c , and k_h c, become constant functions of height, as do the scale heights H and H_*. As a consequence, the vertical wavenumber k_z becomes a constant and the quantity θ vanishes identically for all approximations. When θ is zero, many of the spurious terms disappear from the local dispersion relations.

Figure 1 shows the performance of the three approximations in an isothermal atmosphere. The leftmost panel illustrates the isocontours of the vertical wavenumber achieved in a fully compressible fluid as a function of horizontal wavenumber k_h and temporal frequency ω. The remaining three panels provide the same isocontours for the soundproofing treatment indicated at the top of the panel. The solid black contours in each panel are for the fully compressible fluid, while the dashed red curves show the same contours under the relevant approximation. The value of each contour is marked in the leftmost panel. In each panel, four isocontours of the nondimensional frequency $\varepsilon =\omega /\sqrt{{{gk}}_{h}}$ are overlaid for reference and appear as dotted orange curves. To see how well an approximation reproduces the correct behavior, one should compare the red and black curves within a panel. We would expect that the differences should be small for low values of the nondimensional frequency, i.e., in the lower-right portion of the diagram, and large for high values of ε (upper left). From the four panels, it is clear that all three approximations reproduce the vertical wavenumber well—as long as the dimensionless frequency is small, i.e., ε ≲ 0.3.

Zoom In Zoom Out Reset image size

Given that an isothermal atmosphere is so special (because many of the spurious terms in the dispersion relation vanish), it is wise to examine the behavior of the local dispersion relation in a more complicated atmosphere. We have chosen to examine the vertical wavenumber in a polytropically stratified atmosphere. Such atmospheres have thermodynamic profiles that are power laws in the depth:

$$\begin{eqnarray}&&{\rho }_{0}\left(z\right)=A{(-z)}^{m},\quad {P}_{0}(z)=\displaystyle \frac{{Ag}}{m+1}{\left(-z\right)}^{m+1},\end{eqnarray} \tag{ 82 }$$

$$\begin{eqnarray}&&H(z)=\displaystyle \frac{(-z)}{m},\quad {N}^{2}(z)=\displaystyle \frac{m(\gamma -1)-1}{\gamma }\displaystyle \frac{g}{(-z)},\end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray}&&{H}_{* }(z)=\displaystyle \frac{m+1}{m}(-z),\quad {c}^{2}(z)=\displaystyle \frac{\gamma g}{m+1}(-z).\end{eqnarray} \tag{ 84 }$$

In the expressions above, A is an arbitrary constant and m is the polytropic index. Polytropes can be stably or unstably stratified, depending on the values of the adiabatic index γ and the polytropic index m; if m > (γ − 1)⁻¹, the atmosphere is stable to convective overturning.

A convenient feature of a polytropic atmosphere is that it is self-similar, lacking an intrinsic spatial scale (see Hindman & Jain 2022). Therefore, the local dispersion relation becomes independent of the horizontal wavenumber if we express the frequency in terms of our nondimensional frequency $\varepsilon =\omega /\sqrt{{{gk}}_{h}}$ and we write all of the atmospheric profiles using a nondimensional depth −k_h z. Because of this property, we can generate a single dispersion diagram that illustrates the vertical wavenumber as a function of dimensionless depth and frequency that is valid for all horizontal wavenumbers. Figure 2 presents the resulting dispersion diagram for each treatment of the fluid equations for a polytropic atmosphere with a polytropic index of m = 3 (which is stably stratified for an adiabatic index of γ = 5/3). The leftmost panel is for a fully compressible fluid, and the right three panels are for the three soundproofing formalisms. The blue region in each diagram corresponds to those depths in the atmosphere where a wave of the given frequency is vertically evanescent and the black and red contours have the same meaning as in Figure 1. The upper panels show a range of dimensionless frequency that is wide enough to contain both the branch of low-frequency internal gravity waves and the branch of high-frequency acoustic waves (if present). The lower panels show a zoomed-in view at low frequencies that focuses on the gravity waves. It should be noted that, at a given frequency, there are two turning points where the local vertical wavenumber vanishes. Hence, the internal gravity waves are vertically trapped in a wave cavity for g modes. The turning points are indicated by the thick curves. Similarly, in the fully compressible fluid, the acoustic waves are trapped in a p-mode cavity.

Zoom In Zoom Out Reset image size

The pseudo-incompressible condition does minimal damage to the g-mode cavity, see Figure 2(b). The boundaries move only slightly even for the highest-frequency waves. Further, the vertical wavenumbers within the cavity are weakly affected even for frequencies that are large enough that we might suspect that the low-frequency limit is invalid. The anelastic models fare poorly, however. The fiducial anelastic approximation does a horrendous job of reproducing the wave cavity boundaries. In fact, there appears to be a residual of the acoustic cavity that is highly distorted and appears at frequencies halfway between the acoustic and gravity wave branches. While the LBR form of the anelastic approximation does not have spurious wave cavities at high frequency, it fails to reproduce the boundaries of the g-mode cavity with fidelity. The highest-frequency gravity waves that are vertically propagating have frequencies that are too high by a factor of about one-third. Further, errors in the vertical wavenumber become noticeably large for relatively low values of the dimensionless frequency, ε > 0.1.

In Sections 5.1–5.3, we found that the pseudo-incompressible approximation and the LBR formulation of the anelastic approximation both introduced errors in the Lagrangian pressure fluctuation that appeared at second order. The fiducial anelastic approximation produced errors at first order. So, at first glance, the LBR approximation seems to fare well. However, as we shall soon see, when we consider other wave variables, such as the fluid velocity components, the pseudo-incompressible approximation becomes the clear winner.

In the same manner that one derives Equations (18) and (19) for a fully compressible fluid, similar equations can be derived for each of the approximations. When pseudo-incompressibility is adopted, one obtains the following:

$$\begin{eqnarray}&&{\rho }_{0}u=-\displaystyle \frac{\omega {{gk}}_{h}}{{\alpha }_{\mathrm{PI}}}\left(\displaystyle \frac{d}{{dz}}+\displaystyle \frac{{\omega }^{2}}{g}\right)\delta P,\end{eqnarray} \tag{ 85 }$$

$$\begin{eqnarray}&&{\rho }_{0}w=\displaystyle \frac{i{\omega }^{3}}{{\alpha }_{\mathrm{PI}}}\left(\displaystyle \frac{d}{{dz}}+\displaystyle \frac{{{gk}}_{h}^{2}}{{\omega }^{2}}+\displaystyle \frac{1}{{H}_{* }}\right)\delta P.\end{eqnarray} \tag{ 86 }$$

To see the magnitude of the spurious terms, we nondimensionalize:

$$\begin{eqnarray}&&\tilde{u}=-\displaystyle \frac{{\tilde{\rho }}_{0}^{-1}}{1-{\varepsilon }^{4}+{\varepsilon }^{2}{\tilde{H}}_{* }^{-1}}\left(\displaystyle \frac{d}{d\zeta }+{\varepsilon }^{3}\right)\delta \tilde{P},\end{eqnarray} \tag{ 87 }$$

$$\begin{eqnarray}&&\tilde{w}=\displaystyle \frac{i\varepsilon \,{\tilde{\rho }}_{0}^{-1}}{1-{\varepsilon }^{4}+{\varepsilon }^{2}{\tilde{H}}_{* }^{-1}}\left(1+\varepsilon \displaystyle \frac{d}{d\zeta }+\displaystyle \frac{{\varepsilon }^{2}}{{\tilde{H}}_{* }}\right)\delta \tilde{P}.\end{eqnarray} \tag{ 88 }$$

If one compares these expressions with Equations (43) and (44), it is clear that all spurious terms appear at second order in the dimensionless frequency. Because all of the other fluid variables (i.e., ρ₁, P₁, and s₁) are linear combinations of the two velocity components (see Equation (26)), the wave functions for all of the fluid variables are correct to first order when the pseudo-incompressible approximation is utilized.

Both of the anelastic approximations falter. For the fiducial anelastic approximation, the nondimensional forms for the two velocity components are

$$\begin{eqnarray}&&\tilde{u}=-\displaystyle \frac{{\tilde{\rho }}_{0}^{-1}}{1-{\varepsilon }^{4}+{\varepsilon }^{2}{\tilde{H}}^{-1}}\left(\displaystyle \frac{d}{d\zeta }-\varepsilon {\tilde{N}}^{2}+{\varepsilon }^{3}\right)\delta \tilde{P},\end{eqnarray} \tag{ 89 }$$

$$\begin{eqnarray}&&\tilde{w}=\displaystyle \frac{i\varepsilon \,{\tilde{\rho }}_{0}^{-1}}{1-{\varepsilon }^{4}+{\varepsilon }^{2}{\tilde{H}}^{-1}}\left(1+\varepsilon \displaystyle \frac{d}{d\zeta }+\displaystyle \frac{{\varepsilon }^{2}}{{\tilde{H}}_{* }}\right)\delta \tilde{P},\end{eqnarray} \tag{ 90 }$$

and for the LBR anelastic formulation we obtain

$$\begin{eqnarray}&&\tilde{u}=-{\tilde{\rho }}_{0}^{-1}{\tilde{\alpha }}_{\mathrm{LBR}}^{-1}\left(\displaystyle \frac{d}{d\zeta }-\varepsilon {\tilde{N}}^{2}+{\varepsilon }^{3}\right)\delta \tilde{P},\end{eqnarray} \tag{ 91 }$$

$$\begin{eqnarray}&&\tilde{w}=i\varepsilon \,{\tilde{\rho }}_{0}^{-1}{\tilde{\alpha }}_{\mathrm{LBR}}^{-1}\left(1+\varepsilon \displaystyle \frac{d}{d\zeta }+\displaystyle \frac{{\varepsilon }^{2}}{\tilde{H}}\right)\delta \tilde{P},\end{eqnarray} \tag{ 92 }$$

with ${\tilde{\alpha }}_{\mathrm{LBR}}\equiv 1-{\varepsilon }^{4}+{\varepsilon }^{2}{\tilde{N}}^{2}+{\varepsilon }^{2}{\tilde{H}}^{-1}$. Both have errors in the horizontal velocity that appear at first order (i.e., the term involving $\varepsilon {\tilde{N}}^{2}$). The fiducial anelastic approximation has the added shame that the Lagrangian pressure fluctuation itself is only correct to zero order and hence all fluid variables suffer from the same deficiency. For the LBR approximation, the first-order error in the horizontal velocity u propagates to errors of similar size in the fluctuations of the Eulerian pressure P₁ and density ρ₁.

We can explore conservation of energy under each approximation by computing the vertical energy flux F(z). Using Abel's Identity, as we did for a fully compressible fluid in Section 3.1, we find a general expression for the energy flux that is valid for all three soundproofing treatments:

$$\begin{eqnarray}&&F(z)=-\displaystyle \frac{i{\omega }^{3}}{4{\rho }_{0}\alpha }{ \mathcal W }\left\{\delta P,\delta {P}^{* }\right\}(z).\end{eqnarray} \tag{ 93 }$$

Each approximation generates a distinct form for α and has a different Wronskian because the coefficients of the first-derivative term in the respective ODEs differ.

For the pseudo-incompressible equations, using Equations (62) and (63), we find that the vertical energy flux is a constant function of height:

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal W }\left\{\delta P,\delta {P}^{* }\right\}(z) & = & C\exp \left\{\displaystyle \int \displaystyle \frac{1}{{\alpha }_{\mathrm{PI}}{\rho }_{0}}\displaystyle \frac{d\left({\alpha }_{\mathrm{PI}}{\rho }_{0}\right)}{{dz}}{dz}\right\}\\ & = & C{\alpha }_{\mathrm{PI}}{\rho }_{0},\end{array}\end{eqnarray} \tag{ 94 }$$

$$\begin{eqnarray}&&{F}_{\mathrm{PI}}(z)=-\displaystyle \frac{i{\omega }^{3}\,C}{4}=\mathrm{constant}.\end{eqnarray} \tag{ 95 }$$

Hence, energy is conserved. It is interesting to note that we have not utilized the small parameter in this derivation of the energy flux. So, energy is conserved even when the low-frequency expansions have questionable validity, because the dimensionless frequency is not small.

Performing the same calculations for the two anelastic treatments reveals that the LBR formulation conserves energy (for the same reasons that the pseudo-incompressible equations do) and the fiducial anelastic equations lack energy conservation:

$$\begin{eqnarray}&&{F}_{\mathrm{FA}}(z)=-\displaystyle \frac{i{\omega }^{3}\,C}{4}\displaystyle \frac{{\rho }_{* 0}}{{\rho }_{0}}=\displaystyle \frac{i{\omega }^{3}\,C}{4}{e}^{{s}_{0}(z)/{c}_{p}},\end{eqnarray} \tag{ 96 }$$

$$\begin{eqnarray}&&{F}_{\mathrm{LBR}}(z)=-\displaystyle \frac{i{\omega }^{3}\,C}{4}=\mathrm{constant}.\end{eqnarray} \tag{ 97 }$$

The vertical energy flux F_FA derived from the fiducial anelastic equations depends on the atmosphere's specific entropy density, and thus, in an atmosphere without adiabatic stratification, the wave will deposit or extract energy as it travels.

In numerical simulations, it is hard to overstate the utility in converting the continuity equation from a parabolic prognostic equation to an elliptic PDE constraint, as is accomplished by both the anelastic and pseudo-incompressible approximations:

$$\begin{eqnarray*}&&\displaystyle \frac{\partial \rho }{\partial t}=-{\boldsymbol{\nabla }}\cdot \left(\rho {\boldsymbol{u}}\right)\longrightarrow \left\{\begin{array}{c}{\rm{anelasticity}}\\ {\boldsymbol{\nabla }}\cdot \left({\rho }_{0}{\boldsymbol{u}}\right)=0,\\ \\ {\rm{pseudo-incompressibility}}\\ {\boldsymbol{\nabla }}\cdot \left({P}_{0}^{1/\gamma }{\boldsymbol{u}}\right)=0.\end{array}\right.\end{eqnarray*}$$

In addition to removing sound waves and hence unthrottling the simulation's time step, the imposition of constraints with this form allows the fluid velocity to be expressed using stream functions. Of course, this reduces the number of variables that must be evolved from one time step to the next. However, this is done at the expense of increasing the spatial order of the now-reformulated momentum equations in a stream function form that is now devoid of any elliptic constraints. This may demand auxiliary boundary conditions on the stream functions that are not readily available. Moreover, if linear coupling in the system is treated as explicit in numerical time-stepping algorithms, it is known, specifically for spectral schemes, that the numerical accuracy of the scheme can be degraded at high resolutions. Fortunately, recent advances have shown that this degradation is avoided if linear couplings remain implicit at the expense of using fully coupled implicit time-stepping schemes (Julien & Watson 2009; Marti et al. 2016; Burns et al. 2020; Miquel 2021).

In the derivation of the pseudo-incompressible condition above, two related assumptions are made. First, the Mach number, Ma, of the flows is small, such that the advection timescale is much longer than a sound-crossing time for a typical flow structure. This ensures that fluid motions are in a constant state of pressure equilibration—i.e., the Eulerian pressure fluctuation is small. Second, we have assumed that fluctuations in the potential density are small compared to that of the background state. This later assumption is self-consistent with low-Mach number flows. Notably, unlike the anelastic approximation discussed below, it does not restrict density fluctuations to be small compared to that of the background state. Finally, since we have ignored diffusive effects in the derivation of the pseudo-incompressible constraint, i.e., we have ignored Q in Equation (8), we have made the further assumption that the Péclet number is large, Pe ≫ 1, such that the thermal diffusion timescale is long compared to the advective timescale. To summarize, for the pseudo-incompressible constraint to be valid, we must have the following ordering of timescales:

$$\begin{eqnarray}&&{\tau }_{\mathrm{sound}}\ll {\tau }_{\mathrm{adv}}\ll {\tau }_{\mathrm{diff}},\end{eqnarray} \tag{ 98 }$$

or equivalently, in terms of nondimensional numbers,

$$\begin{eqnarray}&&\mathrm{Ma}\equiv \displaystyle \frac{{\tau }_{\mathrm{sound}}}{{\tau }_{\mathrm{adv}}}=\displaystyle \frac{U}{c}\ll 1,\end{eqnarray} \tag{ 99 }$$

$$\begin{eqnarray}&&\mathrm{Pe}\equiv \displaystyle \frac{{\tau }_{\mathrm{diff}}}{{\tau }_{\mathrm{adv}}}=\displaystyle \frac{{UL}}{\kappa }\gg 1,\end{eqnarray} \tag{ 100 }$$

where U is a typical flow speed, L is a typical length scale, and κ is the thermal diffusivity.

The validity of the anelastic constraint requires the same assumption of low Mach number, Ma ≪ 1, but makes a different stricture on the effectiveness of thermal diffusion. Because we must ignore Eulerian fluctuations of the mass density in the continuity equation, the equation of state dictates that, in addition to small pressure fluctuations, we must have small entropy or temperature fluctuations. In the convection zone of a star or planet, where the stratification is essentially adiabatic, entropy fluctuations are naturally small, anelasticity holds, and the anelastic and pseudo-incompressible conditions are equivalent. However, in a region of stable stratification, the only way that the entropy or temperature fluctuations can remain small is if temperature homogeneity is diffusively maintained across flow structures (see Bannon 1996). This requires that the thermal diffusion time is short compared to the advective timescale. Summarizing, anelasticity requires

$$\begin{eqnarray}&&{\tau }_{\mathrm{sound}},{\tau }_{\mathrm{diff}}\ll {\tau }_{\mathrm{adv}},\end{eqnarray} \tag{ 101 }$$

or equivalently

$$\begin{eqnarray}&&\mathrm{Ma}\ll 1,\qquad \mathrm{Pe}\ll 1.\end{eqnarray} \tag{ 102 }$$

The limitation of low Mach number is easily met in many astrophysical and geophysical applications. Convection is sedate in the Jovian planets, in the Earth's interior, and in the deep layers of low-mass stars. Wave motions and circulations in the stably stratified regions of stars and planets similarly often have low Mach numbers. The requirements on the Péclet number are usually the more restrictive of the two assumptions. For example, the thermal diffusion time in the Sun is typically millions of years; using the solar radius as the length scale, L = R_⊙ ≈ 700 Mm, and a thermal diffusivity appropriate for photon diffusion, κ ∼ 10⁷ cm² s⁻¹, we obtain τ_diff ∼ 16 Myr. If we consider the meridional circulation at the base of the Sun's convection zone and adopt a typical flow speed of 1 m s⁻¹, we obtain an advective timescale of 20 yr, leading to a Péclet number of Pe ∼ 10⁶. Clearly, these motions are not anelastic; thermal diffusion cannot act rapidly enough to eliminate the temperature fluctuations generated by advection. However, the motions do satisfy both of the requirements for pseudo-incompressibility: Ma ≪ 1 and Pe ≫ 1. Although large Péclet numbers are often true from an astrophysical perspective, numerical simulations are often performed in regimes where $\mathrm{Pe}\sim { \mathcal O }(1)$. The anelastic approximation offers no resolution to this problem, but the pseudo-incompressible equations do. The restriction on the Péclet number Pe can be relaxed if the irreversible thermodynamic terms are retained,

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot \left({\rho }_{* 0}{\boldsymbol{u}}\right)=\displaystyle \frac{{\rho }_{* 0}}{{\rho }_{0}}\displaystyle \frac{Q}{{c}_{p}{T}_{0}},\end{eqnarray} \tag{ 103 }$$

where ρ_*0 is the potential density of the background state. Of course, the retention of Q will usually render a stream function formalism without the requirement of an elliptic constraint impossible. Further, we caution that substantial errors can arise relative to the fully compressible system depending on the treatment of the equation of state and the thermal diffusion (e.g., diffusion of entropy versus temperature) in the pseudo-incompressible approximation (Lecoanet et al. 2014).

Finally, we wish to note a final advantage of the pseudo-incompressible approximation over anelasticity. While both soundproofing schemes are well justified in a convection zone where the stratification is nearly adiabatic, if one wishes to simulate both stable and unstable regions in the same computational domain, the pseudo-incompressible approximation allows one to do so smoothly with a uniform treatment. The anelastic approximation will result in flows that violate the underlying assumptions of the approximation.