## Abstract

Stationary gravity waves, such as mountain lee waves, are effectively described by Grimshaw’s dissipative modulation equations even in high altitudes where they become nonlinear due to their large amplitudes. In this theoretical study, a wave-Reynolds number is introduced to characterize general solutions to these modulation equations. This nondimensional number relates the vertical linear group velocity with wavenumber, pressure scale height, and kinematic molecular/eddy viscosity. It is demonstrated by analytic and numerical methods that Lindzen-type waves in the saturation region, that is, where the wave-Reynolds number is of order unity, destabilize by transient perturbations. It is proposed that this mechanism may be a generator for secondary waves due to direct wave–mean-flow interaction. By assumption, the primary waves are exactly such that altitudinal amplitude growth and viscous damping are balanced and by that the amplitude is maximized. Implications of these results on the relation between mean-flow acceleration and wave breaking heights are discussed.

## 1. Introduction

Atmospheric gravity waves generated in the lee of mountains extend over scales across which the background may change significantly. The wave field can persist throughout the layers from the troposphere to the deep atmosphere, the mesosphere and beyond (Fritts et al. 2016, 2018). On this range background temperature and therefore stratification and background density may undergo several orders of magnitude in variation. Also, dynamic viscosity and thermal conductivity cannot be considered constant on such a domain (Pitteway and Hines 1963; Zhou 1995).

Two predominant regimes for the waves can be identified: the homosphere and the heterosphere (Nicolet 1954). These two are separated by the turbopause, usually somewhere in the mesopause region. Below the turbopause molecular viscosity is negligible. Hence, if diffusion occurs, it is caused by turbulence. Due to the missing damping and the thinning background air, mountain waves amplify exponentially when extending upward. This phenomenon is also called anelastic amplification.

At certain heights the amplitudes cannot grow any further due to limiting processes. Those may be static or dynamic instabilities that act on the small scale comparable to the wavelength. For instance, Klostermeyer (1991) showed that all inviscid nonlinear Boussinesq waves are prone to parametric instabilities. The waves do not immediately disappear by the small-scale instabilities, rather the perturbations grow comparably slowly such that the waves persist in their overall structure over several more wavelengths. However, turbulence is produced. Lindzen (1970) modeled the effect of turbulence on the wave by harmonic damping with a constant kinematic eddy viscosity. The eddy viscosity is exactly such that it saturates the wave, meaning that viscous damping and anelastic amplification are balanced (Lindzen 1981; Fritts 1984; Dunkerton 1989; Becker 2012). Pitteway and Hines (1963) referred to this instance as amplitude-balanced wave.

Above the turbopause molecular viscosity takes over. It is usually modeled by a constant dynamic viscosity. In combination with the thinning background density, the kinematic viscosity increases exponentially with height resulting in a rapid decrease in amplitude.

In the process of becoming saturated the amplitude becomes considerably large such that the waves cannot be considered linear. Pioneering work on the mathematical description of nonlinear gravity waves was accomplished by Grimshaw (1972, 1974). This paper aims to extend Lindzen’s linear wave saturation theory with the aid of Grimshaw’s nonlinear wave description.

## 2. How the modulation equations solve the compressible Navier–Stokes equations—A brief review

The nonlinear governing equations for our investigations are the two-dimensional dissipative Grimshaw’s modulation equations being first introduced by Grimshaw (1974). Before presenting them, we want to give a brief review in this section how they solve the dimensionless compressible, Reynolds-averaged Navier–Stokes equations (NSE) asymptotically. Detailed derivations can be also found in Achatz et al. (2010) and Schlutow et al. (2017). Length and time are nondimensionalized via

where *L*_{r} ≈ 1...10 km denotes the reference wavelength (hence the subscript *r*) and *υ*_{r} is the reference velocity. Note that variables labeled with an asterisk denote dimensional quantities throughout the paper. To separate the hydrostatic background from the flow field associated with the wave the first ingredient necessary is a small scale separation parameter

where *H*_{θ} ≈ 10...100 km is the potential temperature scale height. This choice for the scale separation parameter was introduced by Achatz et al. (2010).

The authors considered inviscid flows. To take viscous damping into account, we need to compare inviscid and viscous terms. The (eddy) viscosity is not a constant throughout the atmosphere and among others depends on temperature. Midgley and Liemohn (1966, their Fig. 9) gave a realistic vertical profile of the effective kinematic viscosity combining eddy and molecular effects. Hodges (1969) computed the eddy diffusion by gravity waves near the mesopause to be *K*_{r(eddy)} ≈ 10^{6...7} cm^{2} s^{−1}, which compares well with Midgley and Liemohn (1966) and Lindzen (1981). In the scaling regime of Achatz et al. (2010) as well as Schlutow et al. (2017) these values correspond to Reynolds numbers of Re ≈ 10...100 when the Mach number is Ma ≈ 0.1...0.01. These numbers are also supported by a review of Fritts (1984).

Taking these arguments into account, a realistic flow regime in terms of Mach, Froude, Reynolds, and Prandtl number most suitable for internal gravity waves in the middle/upper atmosphere region is then found by assuming

where *ρ*_{r}, *p*_{r}, *κ*_{r}, and *μ*_{r} represent the reference density, pressure, thermal conductivity, and dynamic viscosity, respectively. The constant *κ* = 2/7 is the ratio of the ideal gas constant *R* to the specific heat capacity at constant pressure *c*_{p} for diatomic gases. Equations (3) provide a distinguished limit for multiple-scale asymptotic analysis. Introducing compressed coordinates for the horizontal, vertical, and time axis,

separates the slowly varying background from the fast oscillating wave fields. Under the assumption of stable stratification to the leading order, the hydrostatic background is determined by the vertical temperature profile *T*(*z*). Two dimensionless background variables, which vary on the large scale, will appear in the final modulation equations, the Brunt–Väisälä frequency or stratification measure *N* and the background density *ρ*. They have to be calculated from the temperature by solving

which originate from the hydrostatic assumption in combination with the ideal gas equation of state.

With the given scale separation of background and wave field, we can formulate the scaled two-dimensional compressible NSE,

in terms of the wave-field variables, velocity **v**, buoyancy *Nb*, and kinematic pressure *p*, where $\u2207^=\u2061(\u2202/\u2202x^,\u2202/\u2202z^)T$ and **e**_{z} the unit vector pointing in the vertical direction; $D/Dt^$ denotes the material derivative, and Λ is the total kinematic viscosity. The dots represent higher-order terms that we do not give explicitly but must not be neglected. The interested reader finds these terms in Schlutow et al. [2017, their (2.19)].

Wentzel–Kramers–Brillouin (WKB) theory provides us with a multiple-scale method to solve the scaled NSE asymptotically by the spectral ansatz

where c.c. stands for the complex conjugate and h.h. for higher harmonics. The vector **U** = (**v**, *b*, *p*)^{T} contains the prognostic variables and *ϕ* is the wave’s phase.

We want to give a remark at this point. By construction, the WKB ansatz is a nonlinear approach. In the limit *ε* → 0 the amplitudes are finite. In fact, the ansatz converges to the nonlinear plane wave of Boussinesq theory, which is known to be an analytical solution. Mean-flow interaction and Doppler shift are leading-order effects, which is in contrast to weakly nonlinear theory where they appear as higher-order corrections. In the *ε* limit, the weakly nonlinear approach converges to the linear plane wave.

The nonlinear WKB ansatz is inserted into the compressible NSE and terms are ordered with respect to powers of *ε* and harmonics.

## 3. Governing equations: Grimshaw’s dissipative modulation equations

To leading order of the nonlinear WKB analysis of the NSE with the turbulence model of Lindzen one obtains Grimshaw’s dissipative modulation equations supported by a mean-flow horizontal kinematic pressure gradient [Grimshaw 1974, their (4.1)–(4.6)]. Note that the leading-order WKB analysis of the dissipative pseudoincompressible equations (Durran 1989) lead to the same modulation equations. We present them in slightly different notation as Grimshaw and with the prerequisite that the wave field is horizontally periodic (0 < *k*_{x} = const) and only modulated in the *z* direction,

The modulation equations govern the evolution of vertical wavenumber *k*_{z}, wave action density *ρa*, and mean-flow horizontal wind *u*. Equations (9a)–(9c) are closed by

where **k**, *ω*, and *ω*′ represent the wavenumber vector, extrinsic frequency, and vertical linear group velocity, respectively. Note that primes denote derivative with respect to the vertical wavenumber throughout this paper. Extrinsic frequency is defined by the sum of intrinsic frequency and Doppler shift. It is linked to the wavenumber vector by the dispersion relation for nonhydrostatic gravity waves. It was shown in Schlutow et al. (2017) that the modulation equations equally hold for hydrostatic waves, where the horizontal wavelength is much larger than the vertical, if the dispersion relation is replaced by *ω* = *Nk*_{x}/|*k*_{z}| + *k*_{x}*u*. The prognostic variables determine the asymptotic solution as described in the previous section and explained in Schlutow et al. (2017).

The first equation, (9a), essentially describes the evolution of the phase in (8) as by definition ∇*ϕ* ≡ **k** and −∂*ϕ*/∂*t* ≡ *ω*. Cross differentiating these expressions, they add up to zero. The second equation, (9b), governs the conservation of wave action density being the ratio of wave energy density and the intrinsic frequency. In terms of the polarization relation the leading-order first harmonics are computable via

with $B=\u2061(\omega \u2212kxu)a/2$ the wave’s buoyancy amplitude. The third equation, (9c), accounts for the acceleration of the mean flow. With a slight abuse of notation (we drop the indices) we obtain the zero harmonics,

The variable *p* corresponding to the mean-flow horizontal kinematic pressure is unknown and needs additional investigation to close the system (9).

The governing equations, (9), can be reformulated in vector form:

with a flux **F** and an inhomogeneity **G** where **y** = (*k*_{z}, *a*, *u*)^{T} is the prognostic vector.

## 4. General stationary solutions of the modulation equations

In this section we will explore general stationary solutions before we focus on particular solutions for which we will present stability analysis in the upcoming sections.

In the inviscid limit (i.e., Λ → 0), the modulation equations assume stationary solutions where ∂*p*/∂*x* = 0, which can be computed analytically by a formula of Schlutow et al. [2017, their (5.20)]. When we multiply (9b) by *k*_{x} and subtract (9c), we obtain

Thus, to be consistent with the inviscid limit, the dissipative modulation equations assume stationary solutions only if the right-hand side of (13) vanishes, which provides eventually a closure for the mean-flow horizontal kinematic pressure gradient,

This result implicates that the mean-flow horizontal kinematic pressure gradient balances the viscous forces acting on the horizontal mean-flow wind. Such a flow configuration is referred to as antitriptic flow in the literature (Jeffreys 1922).

Then a general stationary solution **Y**(*z*) = (*K*_{z}, *A*, *U*)^{T} depicting a mountain lee wave fulfills

with $|K|2=Kx2+Kz2$ and Ω′ = *ω*′(*K*_{z}). Note that we label the stationary solution by capital letters. Given any well-behaved, slowly varying mean-flow horizontal wind *U*(*z*), the remaining two variables of the general stationary solution are given explicitly by

with $F$ being the vertical wave action flux at *z* = 0.

## 5. Lindzen-type mountain lee wave and the wave-Reynolds number

So far, we considered all background variables of our governing PDE as functions of *z*. In this section we will prescribe these functions in a piecewise fashion in order to construct a typical mountain lee wave that gets saturated by some small-scale instability process comparable to Lindzen (1981). The complete solution is divided into an unsaturated as well as a saturated middle-atmospheric part and a deep-atmosphere part.

### a. The unsaturated middle-atmospheric solution

First, we assume that the background atmosphere is piecewise isothermal. From (5), this assumption implies that *N* = const, some typical value for the middle atmosphere. Then (6) gives

where *H* = *κN*^{−2} = const denotes the dimensionless (local) pressure scale height. Second, we assume that the mean-flow horizontal wind is piecewise constant as well, so *U* = const. These assumptions can be weakened, which we will discuss in the concluding section 7. It follows immediately by (18) and horizontal periodicity that *K*_{z} = const making it a plane wave.

In the middle atmosphere viscosity is negligible. Therefore, below the breaking height *z*_{break} the integral in (19) vanishes and the amplitude of the wave grows exponentially with the inverse density,

Because of the vanishing dissipation, the mean-flow horizontal kinematic pressure gradient also vanishes according to (14).

### b. The saturated middle-atmospheric solution

Above *z*_{break}, the wave saturates by some small-scale instability process producing turbulence, which balances the anelastic amplification. The exact kinematic eddy viscosity that keeps the amplitude leveled in (19), such that *A* = const, is then given by

### c. The deep-atmosphere solution

The saturated middle-atmospheric solution is valid below the turbopause at *z*_{turbo} as in the heterosphere the molecular viscosity dominates, which is modeled by a constant dynamic viscosity *μ*_{mol} implying

In the deep-atmospheric region the integral in (19) has also an analytic solution for *A*, which decays quickly for *z* → +∞. Here, typical values *N* and *U* for the deep atmosphere are used.

### d. The wave-Reynolds number

In Fig. 1, such a prototypical mountain wave, which saturates, is illustrated. By inspection of (22), its behavior in the different regions may be characterized by a “wave-Reynolds” number. It can be written as

When reintroducing the dimensional variables by reversing the nondimensionalization of Schlutow et al. (2017), so

and using the definitions (3), the wave-Reynolds number reads

The dimensional linear vertical group velocity is denoted by *C*_{gz} and represents the velocity scale for this type of Reynolds number. The term $D=Hp*\u22121|K|*\u22122$ defines its length scale. Estimates of the wave-Reynolds number in the different regions are depicted in Fig. 1.

## 6. Stability of the saturated wave

In this section we will investigate the saturated mountain wave with respect to stability. In particular, we are interested in the region where Re_{wave} = *O*(1) as here the amplitude is at its maximum and one can expect most likely nonlinear behavior. We assess stability by analyzing the evolution of small perturbations. Due to their smallness, we can linearize the governing equations in vector form (12) around the stationary solution **Y**. Applying the ansatz

for the perturbation, results in an eigenvalue problem (EVP):

where we dropped the hat over **y**. The Jacobian matrices of the flux and inhomogeneity evaluated at **Y** are given by

Note that the second and third row of each Jacobian are linearly dependent. Hence, we can solve the third equation of (29) for

which reduces the dimension of the system, so with a slight abuse of notation, **y** = (*k*_{z}, *a*)^{T} and

If one assumes that the perturbation decays sufficiently fast toward the edges of the region where Re_{wave} = *O*(1) (in Fig. 1 at 70 and 110 km), in other words, that the edges have no influence on the perturbation, then one can extend the domain to the infinities. Consequently, we consider the differential operator associated with the EVP as closed and densely defined on *L*^{2}, the space of vector valued square integrable functions on the real line equipped with the norm

In conclusion, we translated the problem of stability to finding the spectrum of a linear operator of an EVP. Broadly speaking, the spectrum is the set of all $\lambda \u2208C$ for which the EVP has admissible solutions.

The EVP is solved by the Fourier transform

which yields an algebraic equation

It has nontrivial solutions only if the coefficient matrix is singular [det(*M*) = 0], which happens if

This characteristic polynomial has two zeros,

which determine the spectrum of the linear operator of the EVP as curves in the complex plane parameterized by the spatial eigenvalue *μ*, which can also be interpreted as the vertical wavenumber of the perturbation. The real part of *λ* represents the instability growth rate, if it is positive, and the imaginary part is a frequency. Thus, (39) provides also a dispersion relation linking the perturbation’s wavenumber with its frequency.

One can readily show that either *λ*_{1} or *λ*_{2} of (39) has positive real part for reasonable wave solutions implying that they are unconditionally unstable. We want to point out that when *H* → ∞ and Λ → 0, (39) reduces to the spectrum of the inviscid nonlinear Boussinesq plane waves (Schlutow et al. 2019) having the classical modulational stability criterion Ω″ > 0 (Grimshaw 1977).

### a. Transient (in)stability

In this subsection we want to investigate the characteristics of the instability that is presented in the previous section. We put the “in” of the section title into parentheses because there is the possibility that the primary wave “survives” the instability. One particular type of those harmless instabilities is called transient instability. Despite the fact that the instability’s norm grows exponentially in time, the instability vanishes at each given point in space if one waits long enough.

To show such characteristics a prerequisite for mathematical rigor must be fulfilled: the linear differential operator of the EVP must be well-posed, which means its spectrum has a maximum real part or, in other words, the instability growth rate is bounded. The reader finds this tedious verification in the appendix. We want to give an interesting remark at this point. The “well-posedness” depends on a criterion, Ω″ > 0, which turns out to be equivalent to the modulational stability criterion of plane waves in Boussinesq theory.

The key idea for transient instabilities (Kapitula and Promislow 2013) is to ask for solutions of the EVP, (29), in a weighted space $L\alpha 2$ with $\alpha \u2208R$ an exponential weight, such that

Doing so, we get the same curves as in (39) but with *iμ* → *iμ* − *α*, so

When we choose

provided Ω″ > 0 the spectrum is stabilized, so $Re\u2061(\lambda 1,2\alpha )\u22640$ for all $\mu \u2208R$. The negative exponential weight penalizes solutions at −∞. To converge in the weighted norm, the perturbation in *L*^{2} must therefore decay exponentially at −∞ for all times. But simultaneously, its *L*^{2} norm grows exponentially in time. This seeming paradox resolves when the perturbation propagates sufficiently fast toward +∞ (i.e., upward). Perturbations of this behavior are called transient instabilities (Sandstede and Scheel 2000). For illustrative purpose, the unweighted and weighted spectrum of an example wave, which we will discuss in more detail in the following section, are plotted in Fig. 2.

### b. Numerical investigation of the transient instabilities

We use the finite-volume numerical solver presented in Schlutow et al. (2019) for the governing equations, (9), to compute the evolution of a tiny Gaussian initial perturbation of the saturated wave in the region where Re_{wave} = *O*(1). The results are shown in Fig. 3. The simulation is set up by *N* = 1 and *K*_{x} = 1. We discretize the equations on 2000 grid points in *z* and integrate in *t* over 6000 time steps. As can be seen in the figure, the perturbation amplifies exponentially and propagates to the right (i.e., upward), as theory suggests. We can also observe that the perturbation’s wavelength compares with the scale height. Or in other words the instability varies on the large scale. The corresponding spectra to this case as computed by (39) and (41) are plotted in Fig. 2. The amplitude and the vertical wavenumber undergo strong modulations due to the exponentially growing perturbation. Furthermore, the initially constant mean-flow horizontal wind experiences acceleration. The analytic maximum instability growth rate according to (A6) is 2.2 for this particular case. In terms of an approximated *L*_{2} norm that we compute numerically, the actual growth rate of the Gaussian perturbation is found to be 1.6. The difference between theoretical and observed rate occurs because the perturbation is not optimal.

## 7. Summary and discussion

In this paper we investigated nonlinear mountain waves, which are governed by Grimshaw’s dissipative modulation equations being asymptotically consistent with the compressible Navier–Stokes equations and the dissipative pseudoincompressible equations likewise.

We introduced a wave-Reynolds number characterizing the stationary solution. When this dimensionless quantity is of order unity, the wave amplitude saturates by small-scale instabilities and assumes a maximum as anelastic amplification by the thinning background air is exactly balanced by the turbulent damping. We analyzed this regime with respect to modulational stability as nonlinearities dominate for large-amplitude waves. It turned out that transient instabilities emerge that propagate upward. We tested this analysis solving the modulation equations numerically and found excellent agreement with the theory.

In the framework of saturated nonlinear wave theory, which we presented in this paper, the saturated mountain wave does initially not accelerate the mean-flow horizontal wind. Instead, a mean-flow horizontal kinematic pressure gradient emerges that keeps the horizontal wind constant by balancing the viscous forces. The wave persists structurally and loses energy directly to turbulence, which in turn damps altitudinal amplification. Eventually, the mean flow is accelerated by a transient instability in the saturation zone propagating upward while growing and transferring kinetic energy to the mean flow.

Our investigations have two implications being of interest for gravity wave parameterizations in numerical weather prediction and climate modeling. First, the induced mean flow behaves wavelike. Its evolution is governed by a dispersion relation for the linear perturbation. In conclusion, it may be interpreted as an upward-traveling secondary wave of larger scale than the primary wave. Secondary waves with wavelengths comparable to the primary modulation scale were investigated by Vadas and Fritts (2002), Vadas et al. (2003), and Becker and Vadas (2018). The authors propose a generating mechanism based on body forces produced by the dissipating primary wave. In contrast to this model, our secondary wave is generated by direct wave–mean-flow interaction that compares to Wilhelm et al. (2018).

Second, this novel picture of saturated mountain waves may explain the bias between the onset of small-scale instability and the actually observed mean-flow acceleration (Achatz 2007). We give an extension to the established picture where waves become unstable at the breaking height and deposit their momentum and energy onto the mean flow at this level. As it turns out, only in combination with modulational instability does an initially saturated nonlinear wave induce a mean flow. This modulational instability has its own transient velocity and growth rate, which we can compute explicitly. Then the mean-flow acceleration depends on these two quantities.

In our derivations we assumed piecewise analytic solutions being matched. This assumption can be weakened to almost arbitrary background, fulfilling hydrostatics and the ideal gas law, as well as almost any mean-flow horizontal wind. However, these functions of height must be restricted to converge to constant values at the infinities. The resulting spectrum describing the temporal evolution of the perturbation would be much more complicated but still analytically assessable by Fredholm operator theory (Schlutow et al. 2019). In conclusion, our results are valid in a much more realistic atmosphere.

## Acknowledgments

The author thanks Prof. Ulrich Achatz, Prof. Rupert Klein, and Prof. Erik Wahlén for many helpful discussions during the preparation of this article. The research was supported by the German Research Foundation (DFG) through Grants KL 611/25-2 of the Research Unit FOR1898 and Research Fellowship SCHL 2195/1-1.

### APPENDIX

#### Well-Posedness: Are the Instability Growth Rates Finite?

In this appendix we demonstrate that the linear differential operator of the EVP, (29), is well-posed. This property is essential for existence, uniqueness and continuous dependency of the solution for the perturbation on the initial data. In section 6 we showed that the spectrum of the operator extents into the right half of the complex plane rendering the saturated wave unconditionally unstable. But is the spectrum bounded from above on the real axis? Let us denote the discriminant of the square root appearing in *λ*_{1,2} of (39) by *D*_{Y}(*μ*). The real part being the instability growth rate may be expressed by

where

The function Re(*λ*_{1,2}(*μ*)) has no poles. Thus, we are only concerned with its behavior at the infinities. We have

Therefore,

In conclusion the operator is well-posed (Kapitula and Promislow 2013) if Ω″ > 0 and ill-posed otherwise.

## REFERENCES

*The Earth as a Planet*, G. P. Kuiper, Ed., University of Chicago Press, 644–712.

## Footnotes

This article is included in the Multi-Scale Dynamics of Gravity Waves (MS-GWaves) Special Collection.

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).