Modulational stability of nonlinear saturated gravity waves

TL;DR

This paper investigates the stability of nonlinear saturated gravity waves using a theoretical framework, introducing a wave-Reynolds number to analyze wave behavior and potential secondary wave generation.

Contribution

It introduces a wave-Reynolds number to analyze the stability of nonlinear gravity waves and demonstrates their destabilization in the saturation region through analytic and numerical methods.

Findings

Lindzen-type waves destabilize when wave-Reynolds number is around unity.

Secondary waves may be generated via wave-mean-flow interactions.

Implications for wave breaking heights and mean-flow acceleration are discussed.

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 non-dimensional 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, i.e. 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.

Equations

$$(x^{\ast },z^{\ast },t^{\ast })=\left(L_{\mathrm{r}}\hat{x},L_{\mathrm{r}}\hat{z},\frac{L_{\mathrm{r}}}{{\upsilon }_{\mathrm{r}}}\hat{t}\right)$$

$$\epsilon =L_{\mathrm{r}}/H_{\theta }\ll 1$$

$$\frac{{\upsilon }_{\mathrm{r}}}{\sqrt{(1-\kappa )p_{\mathrm{r}}/{\rho }_{\mathrm{r}}}}\equiv \mathrm{Ma}$$

$$\frac{{\upsilon }_{\mathrm{r}}}{\sqrt{g{L}_{\mathrm{r}}}}\equiv \mathrm{Fr}$$

$$\frac{{\rho }_{\mathrm{r}}{\upsilon }_{\mathrm{r}}{L}_{\mathrm{r}}}{{\mu }_{\mathrm{r}}}\equiv \mathrm{Re}$$

$$\frac{c_p{\mu }_{\mathrm{r}}}{{\kappa }_{\mathrm{r}}}\equiv \mathrm{Pr}$$

$$(x,z,t)=(\epsilon \hat{x},\epsilon \hat{z},\epsilon \hat{t})$$

$$N^2$$

$$\frac{1}{\rho }\frac{\mathrm{d}\rho }{\mathrm{d}z}$$

$$\frac{\mathrm{D}\mathbf{v}}{\mathrm{D}\hat{t}}+\hat{\nabla }p-Nb{\mathbf{e}}_z$$

$$\frac{\mathrm{D}b}{\mathrm{D}\hat{t}}+Nw$$

$$\hat{\nabla }\cdot \mathbf{v}$$

$$U(\hat{x},\hat{z},\hat{t};\epsilon )=U_{0,0}(x,z,t)+\left(U_{0,1}(x,z,t)e^{i\varphi (x,z,t)/\epsilon }+c.c.\right)+h.h.+O(\epsilon )$$

$$\frac{\partial k_z}{\partial t}+\frac{\partial \omega }{\partial z}$$

$$\rho \frac{\partial a}{\partial t}+\frac{\partial {\omega }^{\prime }\rho a}{\partial z}$$

$$\rho \frac{\partial u}{\partial t}+\frac{\partial {\omega }^{\prime }{k}_x\rho a}{\partial z}$$

$$|\mathbf{k}{|}^2$$

$$\omega$$

$${\omega }^{\prime }$$

$$\frac{\partial \mathbf{y}}{\partial t}+\frac{\partial \mathbf{F}(\mathbf{y})}{\partial z}=\mathbf{G}(\mathbf{y})$$

$$\rho \frac{\partial }{\partial t}(k_{x}a-u)=\frac{\partial p}{\partial x}-k_x\Lambda |\mathbf{k}{|}^2\rho a$$

$$\frac{\partial p}{\partial x}=k_x\Lambda |\mathbf{k}{|}^2\rho a.$$

$$\omega (K_z,U)$$

$$\frac{\partial {\Omega }^{\prime }\rho A}{\partial z}$$

$$U$$

$$K_z$$

$$A$$

$$\rho (z)={\rho }_0{e}^{-z/H}$$

$$A(z)$$

$$\Lambda ={\Omega }^{\prime }{H}^{-1}|\mathbf{K}{|}^{-2}=\mathrm{const}\text{in }[z_{\mathrm{break}},z_{\mathrm{turbo}}).$$

$$\frac{\partial P}{\partial x}=K_x{H}^{-1}{\Omega }^{\prime }A\rho \text{in }[z_{\mathrm{break}},z_{\mathrm{turbo}}).$$

$$\Lambda (z)=\frac{{\mu }_{mol}}{\rho (z)}\text{in }[z_{\mathrm{turbo}},+\infty ).$$

$$\frac{{\Omega }^{\prime }{H}^{-1}|\mathbf{K}{|}^{-2}}{\Lambda }⪌1.$$

$$H=\frac{\epsilon }{L_{\mathrm{r}}}{H}_p^{\ast },\Lambda =\frac{{\rho }_{\mathrm{r}}}{{\mu }_{\mathrm{r}}}{\nu }^{\ast },N=\frac{L_{\mathrm{r}}}{{\upsilon }_{\mathrm{r}}}{N}^{\ast },\mathbf{K}=L_{\mathrm{r}}{\mathbf{K}}^{\ast },$$

$${\mathrm{Re}}_{\mathrm{wave}}\equiv \frac{C_{gz}D}{{\nu }^{\ast }}.$$

$$\mathbf{y}(z,t)=\hat{\mathbf{y}}(z)e^{\lambda t}$$

$$\lambda \mathbf{y}+\frac{\partial }{\partial z}\left({\frac{\partial \mathbf{F}}{\partial \mathbf{y}}|}_{\mathbf{Y}}\mathbf{y}\right)={\frac{\partial \mathbf{G}}{\partial \mathbf{y}}|}_{\mathbf{Y}}\mathbf{y}$$

Full text

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. 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 upwards. 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 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 Bretherton (1966) and Grimshaw (1972). This paper aims to extent 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 review in this section how they solve the dimensionless compressible, Reynolds-averaged Navier-Stokes equations (NSE) asymptotically. Length and time are nondimensionalized via

[TABLE]

where $L_{\mathrm{r}}\approx 1\dots 10$ km denotes the reference wavelength (hence the subscript r) and $\upsilon_{\mathrm{r}}$ is the reference velocity. To separate the hydrostatic background from the flow field associated with the wave the first ingredient necessary is a small scale separation parameter

[TABLE]

where $H_{\theta}\approx 10\dots 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. In order 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, 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_{\mathrm{{r(eddy)}}}\approx 10^{6\dots 7}\,\mathrm{cm}^{2}\mathrm{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 $\mathrm{Re}\approx 10\dots 100$ when the Mach number is $\mathrm{Ma}\approx 0.1\dots 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

[TABLE]

where $\rho_{\mathrm{r}},p_{\mathrm{r}},\kappa_{\mathrm{r}}$ and $\mu_{\mathrm{r}}$ represent the reference density, pressure, thermal conductivity and dynamic viscosity, respectively. The constant $\kappa=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

[TABLE]

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, that 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 $\rho$. They have to be calculated from the temperature by solving

[TABLE]

that 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,

[TABLE]

in terms of the wave field variables, velocity $\mathbf{v}$, buoyancy $Nb$, and kinematic pressure $p$ where $\hat{\nabla}=(\partial/\partial\hat{x},\,\partial/\partial\hat{z})^{\mathrm{T}}$ and $\mathbf{e}_{z}$ the unit vector pointing in the vertical direction. $\mathrm{D}/\mathrm{D}\hat{t}$ denotes the material derivative and $\Lambda$ is the total kinematic viscosity. The dots represent higher-order terms that we do not give explicitly but must not be neglected.

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

[TABLE]

where c.c. stands for the complex conjugate and h.h. for higher harmonics. The vector $U=(\mathbf{v},\,b,\,p)^{\mathrm{T}}$ contains the prognostic variables and $\phi$ is the wave’s phase. This ansatz is inserted into the compressible NSE and terms are ordered with respect to powers of $\varepsilon$ 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, Eqs. (4.1 – 4.6)). Note that the leading-order WKB analysis of the dissipative pseudo-incompressible 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}=\mathrm{const}$) and only modulated in the $z$-direction,

[TABLE]

where $\mathbf{k}$, $\omega$, and $\omega^{\prime}$ 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. The prognostic variables determine the asymptotic solution as described in the previous section and explained in Schlutow et al. (2017).

The first equation (9a) describes the evolution of the phase in (8) as by definition $\nabla\phi\equiv\mathbf{k}$ and $-\partial\phi/\partial t\equiv\omega$. 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 $U_{0,1}$ are computable from the wave action density. The third equation (9c) accounts for the acceleration of the horizontal mean-flow wind $u_{0,0}$. The variable $p$ corresponding to the mean-flow horizontal kinematic pressure $p_{0,0}$ is unknown and needs additional investigation to close the system. Note that we dropped the indeces for the mean-flow variables.

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

[TABLE]

with a flux $\mathbf{F}$ and an inhomogeneity $\mathbf{G}$ where $\mathbf{y}=(k_{z},a,u)^{\mathrm{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. $\Lambda\rightarrow 0$, the modulation equations assume stationary solutions where ${\partial p/\partial x=0}$ which can be computed analytically by a formula of Schlutow et al. (2017, Eq. 5.20). When we multiply (9b) by $k_{x}$ and subtract (9c), we obtain

[TABLE]

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

[TABLE]

This result implicates that the mean-flow horizontal kinematic pressure gradient balances the viscous forces acting on the horizontal mean-flow wind.

Then, a general stationary solution $\mathbf{Y}(z)=(K_{z},\,A,\,U)^{\mathrm{T}}$ depicting a mountain lee wave fulfills

[TABLE]

with $|\mathbf{K}|^{2}=K_{x}^{2}+K_{z}^{2}$ and $\Omega^{\prime}=\omega^{\prime}(K_{z})$. Note that we label the stationary solution by capital letters. Given any reasonable mean-flow horizontal wind $U(z)$, the remaining two variables of the general stationary solution are given explicitly by

[TABLE]

with $\mathcal{F}$ an integration constant.

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.

5.1 The unsaturated middle atmospheric solution

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

[TABLE]

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

In the middle atmosphere viscosity is negligible. Therefore, below the breaking height $z_{\mathrm{break}}$ the integral in (17) vanishes and the amplitude of the wave grows exponentially with the inverse density,

[TABLE]

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

5.2 The saturated middle atmospheric solution

Above $z_{\mathrm{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 (17), such that $A=\mathrm{const}$ , is then given by

[TABLE]

In this region, the mean-flow horizontal kinematic pressure gradient depends on $z$ only by the background density and can be computed inserting (20) into (12),

[TABLE]

In particular, it is positive while the mean-flow horizontal wind is negative which can be seen from (9e) and (13).

5.3 The deep atmosphere solution

The saturated middle atmospheric solution is valid below the turbopause at $z_{\mathrm{turbo}}$ as in the heterosphere the molecular viscosity dominates which is modeled by a constant dynamic viscosity $\mu_{mol}$ implying

[TABLE]

In the deep atmospheric region the integral in (17) has also an analytic solution for $A$ which decays quickly for $z\rightarrow+\infty$. Here, typical values $N$ and $U$ for the deep atmosphere are used.

5.4 The wave-Reynolds number

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

[TABLE]

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

[TABLE]

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

[TABLE]

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=H_{p}^{\ast-1}|\mathbf{K}|^{\ast-2}$ 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 $\mathrm{Re}_{\mathrm{wave}}=\mathit{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 (10) around the stationary solution $\mathbf{Y}$. Applying the ansatz

[TABLE]

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

[TABLE]

where we dropped the hat over $\mathbf{y}$. The Jacobian matrices of the flux and inhomogeneity evaluated at $\mathbf{Y}$ are given by

[TABLE]

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

[TABLE]

which reduces the dimension of the system, so with a slight abuse of notation, ${\mathbf{y}=(k_{z},a)^{\mathrm{T}}}$ and

[TABLE]

If one assumes that the perturbation decays sufficiently fast towards the edges of the region where $\mathrm{Re}_{\mathrm{wave}}=\mathit{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 extent 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

[TABLE]

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\in\mathbb{C}$ for which the EVP has solutions.

The EVP is solved by the Fourier transform

[TABLE]

which yields an algebraic equation

[TABLE]

It has nontrivial solutions only if the coefficient matrix is singular which happens if

[TABLE]

This characteristic polynomial has two zeros

[TABLE]

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

One can readily show that either $\lambda_{1}$ or $\lambda_{2}$ of (37) has positive real part for reasonable wave solutions implying that they are unconditionally unstable. We want to point out that when $H\rightarrow\infty$ and $\Lambda\rightarrow 0$, (37) reduces to the spectrum of the inviscid nonlinear Boussinesq plane waves (Schlutow et al., 2018) having the classical modulational stability criterion $\Omega^{\prime\prime}>0$ (Grimshaw, 1977).

6.1 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 its norm grows exponentially in time, it vanishes at each given point.

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, $\Omega^{\prime\prime}>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 (27) in a weighted space $L^{2}_{\alpha}$ with $\alpha\in\mathbb{R}$ an exponential weight, such that

[TABLE]

Doing so, we get the same curves as in (37) but with $i\mu\rightarrow i\mu-\alpha$, so

[TABLE]

When we choose

[TABLE]

provided $\Omega^{\prime\prime}>0$ the spectrum is stabilized, so $\Re(\lambda_{1,2}^{\alpha})\leq 0$ for all $\mu\in\mathbb{R}$. The negative exponential weight penalizes solutions at $-\infty$. To converge in the weighted norm, the perturbation in $L^{2}$ must therefore decay exponentially at $-\infty$ for all times. But simultaneously, its $L^{2}$-norm grows exponentially in time. This seeming paradox resolves when the perturbation propagates sufficiently fast towards $+\infty$, i.e. upwards. Perturbations of this behavior are called transient instabilities (Sandstede and Scheel, 2000). For illustrative purpose, the unweighted and weighted spectrum of an example wave, that we will discuss in more detail in the following section, are plotted in Fig. 2.

6.2 Numerical investigation of the transient instabilities

We use the finite-volume numerical solver presented in Schlutow et al. (2018) for the governing equations (9) to compute the evolution of a tiny Gaussian initial perturbation of the saturated wave in the region where ${\mathrm{Re}_{\mathrm{wave}}=\mathit{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 and integrate over 6000 time steps. As can be seen in the figure, the perturbation amplifies exponentially and propagates to the right, i.e. upwards, 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 (37) and (39) 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 (46) 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 pseudo-incompressible 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 upwards. We tested this analysis solving the modulation equations numerically and found excellent agreement with the theory.

In the framework of saturated nonlinear wave theory, that 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 upwards while growing and transferring kinetic energy to the mean flow.

Our investigations have two implications being of interest for gravity wave parametrizations in numerical weather prediction and climate modeling. First, the induced mean-flow behaves wave-like. 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); 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., 2018). In conclusion, our results are valid in a much more realistic atmosphere.

Acknowledgments

The author thanks Prof. Ulrich Achatz and Prof. Rupert Klein 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 A Well-posedness: Are the instability growth rates finite?

In this appendix we demonstrate that the linear differential operator of the EVP (27) is well-posed. This property is essential for existence, uniquenes 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 $\lambda_{1,2}$ of (37) by $D_{Y}(\mu)$. The real part being the instability growth rate may be expressed by

[TABLE]

where

[TABLE]

The function $\Re(\lambda_{1,2}(\mu))$ has no poles. Thus, we are only concerned with its behavior at the infinities. We have

[TABLE]

Therefore,

[TABLE]

In conclusion the operator is well-posed (Kapitula and Promislow, 2013) if $\Omega^{\prime\prime}>0$ and ill-posed otherwise.