Spectral stability of nonlinear gravity waves in the atmosphere

We apply spectral stability theory to investigate nonlinear gravity waves in the atmosphere. These waves are determined by modulation equations that result from Wentzel-Kramers-Brillouin theory. First, we establish that plane waves, which represent exact solutions to the inviscid Boussinesq equations, are spectrally stable with respect to their nonlinear modulation equations under the same conditions as what is known as modulational stability from weakly nonlinear theory. In contrast to Boussinesq, the pseudo-incompressible regime does account for the altitudinal varying background density. Second, we show for the first time that upward-traveling wave fronts solving the inviscid modulation equations, that compare to pseudo-incompressible theory, are unconditionally unstable. Both inviscid regimes turn out to be ill-posed as the spectra allow for arbitrarily large instability growth rates. Third, a regularization is found by including dissipative effects. The corresponding traveling wave solutions have localized amplitude and blow up unconditionally by embedded eigenvalue instabilities but the instability growth rate is bounded from above. Additionally, all three types of nonlinear modulation equations are solved numerically to further investigate and illustrate the nature of the analytic stability results.


Introduction
Gravity waves have a significant impact on the dynamics of the earth's atmosphere. Since the resolution and the upper boundary of numerical weather and climate models steadily increase, a better understanding of these waves becomes necessary in order to construct more precise subgrid-scale parametrizations [Fritts and Alexander, 2003]. Most gravity waves are excited in the lower atmosphere. However, they may propagate deep into the higher atmospheric layers above the stratopause. In these regions amplitudes have become so large, due to the thin background air, that linear theory is not applicable anymore.
Pioneering work on nonlinear gravity waves was accomplished by Bretherton [1966] and Grimshaw [1972]. These authors applied Wentzel-Kramers-Brillouin theory to find leading-order asymptotic equations to the Euler equations. These asymptotic equations may be called modulation equations as they describe the evolution of the wave characteristics like amplitude, frequency and wave number [Whitham, 1974].
The stability of gravity waves is of particular interest for modelers of subgrid-scale parametrizations to predict the wave's breaking height. Spectral stability theory is an especially useful approach to study nonlinear waves [Sandstede and Scheel, 2000, Sandstede, 2002, Kapitula and Promislow, 2013.
The key idea is to linearize the equations around a given wave solution which translates the problem of stability to finding the spectrum of a differential operator L. The spectrum consists of all complex λ for which the operator L − λ is not invertible.
The paper is structured as follows. We consider two-dimensional gravity waves being horizontally homogeneous in an isothermal atmosphere. Our analysis will start with the inviscid Boussinesq regime in section 2. Here, plane waves do not only fulfill the nonlinear modulation equations analytically but also the Boussinesq equations themselves. From weakly nonlinear theory it is established in the literature that the plane wave is prone to modulational instabilities. We apply spectral stability analysis on the nonlinear modulation equations and find the very same criterion for modulational instability. The nature of the instability is also examined numerically.
The Boussinesq equations are limited in the sense that they do not cover anelastic growth of the amplitude. When a gravity wave propagates upwards its amplitude may increase due to the decreasing background density. This phenomenon influences the breaking height as the stability depends sensitively on amplitude. Therefore, we continue by investigating an extended set of modulation equations in section 3 which agrees with the inviscid pseudoincompressible regime [Durran, 1989, Achatz et al., 2010. Here, the background density is an explicit function of height. These modulation equations possess upward-traveling wave solutions that take the form of a front where the envelope velocity is always greater than the derivative of the frequency with respect to the vertical wave number, what one may call the linear group velocity. That those two velocities are not the same, is a consequence of the nonlinearity. The spectral stability analysis reveals that the traveling wave fronts are unconditionally unstable. The dynamics of the instabilities will be illustrated numerically.
The spectra of both inviscid regimes allow instabilities with arbitrarily large growth rates which is evidence for ill-posedness [Fetecau and Muraki, 2011]. In section 4, a regularization will be found by including dissipation to the pseudo-incompressible modulation equations in accordance with the compressible Navier-Stokes equations. These modulation equations exhibit both upward and downward-traveling wave solutions where the amplitude is of finite extent making them wave packets. In terms of the spectral stability analysis we find out that the wave packets are also unconditionally unstable. The waves destabilize by a superposition of exponentially growing eigenfunctions. In contrast to the inviscid waves, a maximum instability growth rate can be computed explicitly. We will solve the modulation equations numerically for the wave packets to demonstrate the evolution of these instabilities.

Inviscid Boussinesq plane waves
We start our investigations with the modulation equations for gravity waves in the Boussinesq approximation. Two spatial dimensions (X, Z), one horizontal and one vertical, are assumed and the wave shall be homogeneous in X. For the derivation of these equations see, e.g., Muraschko et al. [2015]. We may give this coupled set of non-dimensionalized nonlinear equations here in conservative flux form (1) The prognostic variables k z , A and u denote the vertical wavenumber, the wave action density and the mean-flow horizontal wind, respectively. Their domain in Z is the real line. The horizontal wavenumber K x is constant and without loss of generality positive, hereinafter. Two additional positive constant parameters, the Brunt-Väisälä frequency N and the background density ρ, account for the ambient atmosphere which is assumed to be uniformly stratified. The intrinsic frequency is defined throughout this article by the dispersion relation for non-hydrostic gravity waves as a function of k z , soω and its derivativeω which denotes the linear group velocity. Throughout this article primes denote the derivative with respect to the vertical wave number. The first equation of (1) describes the conservation of the wave's troughs and crests. Its flux is exactly the extrinsic frequency, i.e. the Dopplershifted intrinsic frequency. The second equation governs the conservation of wave action density being the ratio of wave energy density and the intrinsic frequency. The third equation accounts for the acceleration of the mean flow which is related to the Stokes drift known from water waves.
The sign ofω is ambiguous, but can be chosen without loss of generality to be positive due to the symmetry of the problem. Consequently, the wave action density has to be strictly positive sinceωA represents the wave energy density.

Plane wave solution
For simplification we introduce the "specific" wave action density a by and rewrite system (1) in vector form for our solution vector y = (k z , a, u) T with the flux function Equation (5) has a stationary solution Y = (K z , A, U ) T = const. that satisfies as F is autonomous which corresponds to a monochromatic plane wave. It is well known that plane waves are solutions to the fully nonlinear Boussinesq equations [Mied, 1976].

Spectral stability of the plane waves
In order to investigate the stability of the plane waves Y , we perturb (5) according to Y + y and linearize around the stationary solution, which yields with the Jacobian Since the Jacobian of the flux does not depend on time, (8) must have a solution of the form that gives us an eigenvalue problem involving the differential matrix operator As L Y is a differential operator and therefore operates on an infinitedimensional space, in general one cannot expect to obtain only discrete eigenvalues, λ ∈ C, like one would get for finite-dimensional matrices. Instead L Y generates a spectrum which is the set of λ's such that the operator L Y − λ is not invertible. The stability problem can be transfered to finding the spectrum of the operator L Y . We say the solution Y is spectrally stable if the spectrum of L Y is contained in the left hand side of the complex plane, so for every λ in the spectrum ℜ(λ) ≤ 0 must hold. We can motivate this definition when we consider (10). If there is a point λ in the spectrum with positive real part, then the corresponding eigenfunction for the perturbation grows exponentially in time with the real part of λ being the instability growth rate and the solution blows up. For negative or vanishing real part the eigenfunction decays or remains bounded in time, respectively, such that the solution is stable.
We consider L Y as a closed, densely defined operator on L 2 (R, R 3 ), the space of vector valued square integrable functions on the real line equipped with the norm The eigenvalue problem (11) can be reformulated as a linear ordinary differential equation (ODE) where B(λ) = −DF −1 (Y )λ is a spatially constant matrix. The ODE (14) has a general solution in terms of the Fourier transform which translates the ODE into an algebraic equation having nontrivial solution only if det(B(λ) − iµ) = 0. This constraint gives us a parametrization for the spectrum of L Y in the complex plane by three different zeros The spectrum in the complex plane is depicted in figure 1. It consists of straight lines going through the origin. As long as the discriminant in (18) and (19) is positive the spectrum is the imaginary line (panel a) and we can conclude that the plane wave is spectrally stable. If however the discriminant is negative, the spectrum splits into two straight lines that extend into the right half of the complex plane (panel b) and the wave becomes spectrally unstable. This can only happen whenω ′′ (K z ) becomes negative, which occurs if This result is well known and commonly referred to as modulational instability. Whitham [1974], Grimshaw [1977] and Sutherland [2001] derived this instability criterion from weakly nonlinear theory where the amplitude of the wave is used as the expansion parameter. Therefore, the coupling to the mean flow is a higher order effect. However, the modulation equations from the WKB theory, that are presented here, turn out to be nonlinear in the leading order as the wave-mean-flow interaction and the Doppler shifting of the frequency appear in the coupled set of leading order modulation equations.
With the stability analysis of this section we know under what circumstances a plane wave becomes modulationally unstable, but how does the perturbation of an unstable wave evolve in time?

Numerical simulation of the plane waves
We explore the nature of the modulational instability further by numerical computations of the evolution of the unstable plane wave. For this endeavor a finite-volume scheme was implemented to solve (1). A precise description of the algorithm is given in appendix A. The computational results are plotted in figure 2 for an initially unstable wave according to the modulational stability criterion of (20). The axes are labeled by the dimensional height z and time t. The length scale is fixed by a reference wave length L r and the time scale is given by a reference Brunt-Väisälä frequency N r . The small scale separation parameter ε, being the WKB-expansion parameter, measures the ratio between L r and the length scale on which the wave's envelope varies. To seed the instability, the initial conditions are perturbed by a Gaussian peak in the middle of the domain of the order of magnitude 10 −5 . The parameters are set to N = 1.2 and K x = 1.8. We observe that the perturbation, that excites instability modes, blows up quickly in all three prognostic variables. The instability possesses a small scale structure with short wavelengths.

Inviscid pseudo-incompressible waves
In this section we want to extend the model in order to incorporate anelastic effects as well which are omitted in the Boussinesq regime. The modulation equations for horizontally homogeneous waves in the pseudo-incompressible regime exposed to an isothermal inviscid medium may be written as Here, in contrast to the Boussinesq model the background density for the isothermal ambient atmosphere is a function of height, It can be derived from the hydrostatic equation and the equation of state for ideal gases. The system (21) was introduced by Grimshaw [1972] to describe nonlinear internal gravity waves in an inviscid compressible fluid in the environment of a slowly varying medium. Achatz et al. [2010] demonstrated that these are WKB solutions to the compressible Euler equations and likewise to the pseudo-incompressible equations of Durran [1989] when applying the same asymptotic scaling. In contrast to the Boussinesq equations they can capture the well known anelastic amplification of wave packets when they travel upwards. A derivation of these modulation equations is also presented in Schlutow et al. [2017]. The Euler equations are scaled in terms of a small parameter where L r denotes a reference wave length and H θ the potential temperature scale height. A distinguished limit for gravity waves is found to be i.e. assuming a small Mach and Froude number. The nonlinear wave ansatz allows that the envelope, the phase and the induced mean flow vary on the same length scale as the background given by H θ while the wave itself oscillates on a much shorter scale L r . We call the ansatz nonlinear because the mean flow interacts with the wave field already in the leading order.

Traveling wave solutions
To get rid of the varying coefficient in (21) due to the height-depending background density we can again introduce the specific wave action density via A(Z, T ) = ρ(Z)a(Z, T ). When we now write our system in vector form, an additional inhomogeneity G appears on the right hand side because of the derivative of the background density with respect to Z. The inhomogeneity is given by and is the same flux as in (5). By the substitution with the specific wave action density both the flux and the inhomogeneity become autonomous. Therefore, the modulation equation can only assume a stationary solution Y = (K z , A, U ) T if the inhomogeneity becomes zero. This is the case when A = 0 orω ′ (K z ) = 0. None of the according solutions corresponds to a gravity wave, so we omit them. Alternatively, one can ask for traveling wave solutions, i.e. solutions of the form Y (Z, T ) = Y (Z − CT ) where C is a constant that we will refer to as envelope velocity.
A traveling wave solution Y (ζ), where ζ = Z − CT , must then fulfill This system of ODEs can be be simplified to a single first-order ODE by integrating the first equation giving Integration of the third minus the second equation results in Equations (29) and (30) are diagnostic equations for the specific wave action density and the mean-flow horizontal wind that we can treat as functions of K z now. The two integration constants A 0 and U 0 must be determined by some initial conditions. The single ODE reads then where primes denote the derivative with respect to K z . Schlutow et al. [2017] showed that (31) has solutions that represent upward traveling wave fronts. These solutions converge at the infinities to constant values which are called asymptotic rest states, so The asymptotic rest states correspond to the zeros of f in (31), i.e.ω ′ (K − z ) = 0 and A(K + z ) = 0. According to (2) and (30), one finds that K − z = 0 and is a necessity to guarantee that the denominator of f does not assume a zero such that f stays bounded. Note in particular that C > 0 so the waves travel upwards. Given a K x and K + z Schlutow et al. [2017] computed a critical envelope velocity being computed in terms of (30) and (29). In optics or signal processing one would identify this traveling wave front as a "down-chirp" as it drops in frequency when observed at a fixed height. The mean-flow horizontal wind u is accelerated in positive direction in a fashion such that a persistent mean flow is induced. Remarkably, traveling wave backs as solutions to (21) cannot exist which we will show in appendix B. Also note that linear group and envelope velocity for nonlinear waves are not necessarily the same. We defined the linear group velocity as the derivative of the frequency with respect to the wave number. It actually only coincides with the envelope velocity if the wave packet is quasi-monochromatic. Quasi-monochromatic waves may be found in linear and weakly nonlinear theory. Here in contrast we allow the vertical wavenumber to vary on the same order as the amplitude and the mean flow, so the prerequisite for monochronicity is generally not given. Insteadω ′ describes the velocity of the group only locally in the vicinity of a given k z . For a detailed derivation and discussion we refer to Majda [2003, pp 86].

Stability of the inviscid traveling wave front
To assess the stability of the inviscid pseudo-incompressible traveling wave fronts we first recast (25) in translational coordinates We benefit from this coordinate transformation because the traveling wave front Y = Y (ζ) is a stationary solution in this system. Second, we linearize the transformed equation (25) around the traveling wave solution with the Jacobians and Since the Jacobians of the flux and the inhomogeneity are independent of the time variable τ by autonomy and the stationarity of Y , equation (37) must have solutions of the form which yields an eigenvalue problem (L Y − λ)y = 0. The linear operator is given by As for the Boussinesq plane waves we also translated the problem of linear stability into the task of finding the spectrum of a linear operator. However, in contrast to the plane waves this operator exhibits varying coefficients due to the spatial dependency of the traveling wave front Y = Y (ζ). In conclusion, we cannot compute the spectrum straightforwardly by a Fourier transformation as we did for the plane waves.
But fortunately, analytical progress is possible utilizing the asymptotic rest states. The idea is essentially to approximate the operator using these states by an asymptotic operator for which the spectrum is easy to compute. We can connect the asymptotic to the original operator by Fredholm theory. For a detailed introduction on this topic we refer to Sandstede [2002], Kapitula and Promislow [2013].
Let L be a closed linear operator defined on some Hilbert space X with a domain being dense in X. The operator L is called Fredholm if the dimension of its kernel and the codimension of its range are finite. The Fredholm index is then defined by In terms of the above definition we can separate the spectrum of L into two sets, the essential and the point spectrum.
• The essential spectrum of L is the set of all λ ∈ C such that one of the following is true: either L − λ is not Fredholm, or L − λ is Fredholm but ind(L − λ) = 0.
• The point spectrum of L is the set of all λ ∈ C such that L − λ is Fredholm and ind(L − λ) = 0 but L − λ is not invertible.
With the aid of these definitions one can now apply Weyl's essential spectrum theorem [Kapitula and Promislow, 2013, pp 29] which states that a well-behaved perturbation L of a Fredholm operator L 0 is Fredholm, that the Fredholm indices are the same, and that L and L 0 share the same essential spectrum.
Let us define the asymptotic operator to (41) as One can readily show that the original operator L Y is indeed a well-behaved perturbation of L ∞ in the required sense [Kapitula and Promislow, 2013, pp 39]. Therefore, we can compute the essential spectrum of L ∞ and make use of Weyl's essential spectrum theorem by identifying it with the essential spectrum of L Y . The asymptotic eigenvalue equation (L ∞ − λ)y = 0 can be recast as the ODE having a piecewise constant coefficient matrix generated by From (41) we find that The eigenvalues ν ± (λ) of B ± (λ) contain all the information we need to characterize the essential spectrum. The operators L ∞ − λ and hence L Y − λ are Fredholm if, and only if, all eigenvalues ν ± (λ) have non-zero real part, i.e. B ± (λ) is hyperbolic [Sandstede, 2002]. The Fredholm index can be written as where i ± are the Morse indices of B ± (λ). The Morse index of a hyperbolic matrix is the dimension of its unstable subspace. So, if for a λ ∈ C the operator is Fredholm and the Morse indices differ, i.e. i − (λ) = i + (λ), then λ is in the essential spectrum of L Y . We can circumnavigate the computation of the Morse indices for every λ ∈ C by investigating the boundaries of the essential spectrum directly. The essential spectrum is delimited by curves where the λ's are such that the matrices B ± (λ) fail to be hyperbolic, that is where some eigenvalues become purely imaginary, so ν ± (λ) = iµ and µ ∈ R. These λ's satisfy the dispersion relation Plugging (46) into (48) and simplifying provides We substitute the asymptotic rest states from (35) and solve for λ which results in six curves in the complex plane being parametrized by µ where we also identified and replacedω ′′ (K − z = 0) = −N/K 2 x . These curves are denoted as Fredholm borders. They are the boundaries of the regions in the complex plane where the operator L Y − λ is Fredholm. For a typical set of parameters defining a pseudo-incompressible inviscid traveling wave front, the Fredholm borders are plotted in figure 3. Arrows depict their orientations which are determined by the direction of increasing µ. As λ crosses a Fredholm border the Fredholm index may change. And indeed, one can prove that ind(L Y − λ) • increases by 1 when crossing the graph of λ + (µ) from right to left with regard to its orientation.
• decreases by 1 when crossing the graph of λ − (µ) from right to left with regard to its orientation.
An asymptotic analysis reveals that real valued λ ≫ 1 are not in the essential spectrum. This means that the Fredholm index in the regions containing these λ's is zero. Since we know the Fredholm index in one region, we can compute the other indices of the remaining regions by the rule from above. The regions with Fredholm index other than zero belong to the essential spectrum. Now, as we know the essential spectrum of the operator, we can analyze it with respect to instabilities which means to check for ℜ(λ) > 0 in the essential spectrum. We can see in figure 3 that, for this particular choice of parameters, the spectrum reaches far into the right half of the complex plane. And indeed we find from (52) that λ + 3 always has positive real part −ηω ′ (K + z ) because η < 0 as the background density decreases with height andω ′ (K + z ) > 0 which is a prerequisite for the existence of the traveling wave front. Also, it can be shown that ℜ(λ − 2 ) > 0 independent of the choice of parameters. We can therefore deduce that the traveling wave front is spectrally unstable with respect to the modulation equations for every set of parameters defining the wave and the background.

Numerical simulation of the inviscid traveling wave front
We want to compute the evolution of the essential spectrum instability of the inviscid traveling wave front numerically to learn about its dynamics. For this purpose, the modulation equations (21) are integrated in time with the numerical method presented in appendix A. The parameters are set to N = 1.7, K x = 1.5 and C = 0.8. By integrating equation (31) numerically we obtain the initial conditions.
The numerical results of the time integration are plotted in figure 4. We labeled the axes with the dimensional height z and time t where N 2 r = g/H θ and g the Earth's acceleration according to Schlutow et al. [2017]. An unstable wave packet grows rapidly directly behind the front. The simulation is terminated by the numerics being unable to handle the small-scale instability. For typical values for the mesosphere, H θ = 20 km and L r = 1 km, we can estimate that after only 15 min the wave blows up.

Dissipative Grimshaw waves
So far, we have seen in section 2 and 3 that both the inviscid Boussinesq plane wave and the inviscid pseudo-incompressible traveling wave front exhibit instabilities which grow rapidly in our numerical simulations. They also suffer from an unphysical spectrum. Both spectra allow for unstable modes with arbitrarily large instability growth rates which makes instabilities potentially catastrophic. According to Kapitula and Promislow [2013, p 61] the operators are not well-posed. This problem originates intrinsically from the nature of the Euler equations that we approximated. By omitting diffusive effects we removed the damping which is especially efficient in suppressing high-frequency perturbations.
In the Boussinesq case, Fetecau and Muraki [2011] proposed a regularization by including higher order terms to the dispersion relation. In contrast to this ad hoc method we want to introduce a regularization to the pseudoincompressible regime by incorporating dissipative, instead of dispersive, effects in an asymptotic fashion.
The modulation equations for nonlinear gravity waves including dissipation according to Grimshaw [1974] assuming horizontal homogeneity and isothermal background are given by Note that the only difference to the inviscid case of the previous section is an additional term on the right hand side that accounts for the dissipative effects. It consists of the squared norm k 2 (k z ) = K 2 x + k 2 z as function of k z , the wave action density and a variable Λ that comprises the effects of viscosity and thermal conductivity.
The dissipative modulation equations can be derived from the compressible Navier-Stokes equations in the same fashion as if one starts from the compressible Euler equations. The small parameter is likewise the ratio of reference wave length and potential temperature scale height. The enhanced distinguished limit (cf. equation 24) from the asymptotic scaling reads then So, we assume a large Reynolds number and by the Prandtl number that the viscosity is of the same order as the thermal conductivity. This scaling is realistic for two regions. Either the lower atmosphere when we treat the equations as Reynolds averaged with a simple turbulence model such that the Reynolds number represents eddy viscosity. Or for the regions of the higher atmosphere, i.e. the mesosphere and lower thermosphere, where molecular viscosity becomes large. The nonlinear WKB wave ansatz is the very same as for the inviscid case. The phase, wave envelope, mean flow and background vary on the same length scale H θ . Mean flow and wave field interact in leading order. A thorough derivation of the resulting modulation equations can be found in Grimshaw [1974] which is indeed a consistent dissipative extension of Schlutow et al. [2017] with regard to the scaling presented here.
The variable Λ is a function of altitude Z and is found to be proportional to the inverse background density which implies that the impact of dissipative effects increases with height. However, we treat Λ as a constant henceforth in order to continue with our analytical approach. On the one hand, we are interested in a regularization of the inviscid traveling wave fronts. Setting Λ constant enables us to construct an autonomous flux and inhomogeneity giving rise to traveling wave solutions with asymptotic rest states similar to the inviscid waves. It is reasonable to assume that the non-autonomy caused by a varying Λ may change the solution but the stability analysis based on the asymptotic rest states is most likely untouched though.
On the other hand, Grimshaw [1974] argues that this assumption is not a great restriction as a varying dissipation is only of importance in the vicinity of critical levels where the extrinsic frequency approaches zero such that the wave action density amplifies. In our scenario there are no critical levels present. Grimshaw [1974] also tested the influence of varying versus constant Λ numerically and found no significant difference.

The traveling wave solution
To construct traveling wave solutions to the dissipative nonlinear modulation equations (56) we exploit the transformation generated by the specific wave action density a(Z, T ) = ρ −1 (Z)A(Z, T ) where ρ is the same exponentially decreasing background density as before. This transformation removes the Z-dependencies, such that the coefficients of the resulting system become constant. Written in vector form the system reads for the solution vector y = (k z , a, u) T . The flux and the inhomogeneity are given by that states a three-dimensional ODE. We can integrate the first equation of the system, since the right hand side is zero, which yields a diagnostic equation for the mean-flow horizontal wind with an integration constant U 0 that depends on the initial conditions. Thus, we reduce the dimension to two and U is now treated as a function of K z . With a mild abuse of notation we rewrite the solution vector to Y = (K z , A) T and the system of ODEs (63) becomes The flux and inhomogeneity are updated by substituting (64), so and The ODE (65) has equilibria if the inhomogeneity vanishes. Those points are Y = (K z , 0), i.e. all points on the K z -axis in the two-dimensional phase space. So, there may be solutions with two asymptotic rest states, Y + and Y − , on the K z -axis where the specific wave action density vanishes at the infinities and the vertical wave number approaches some constant values. The solution is then the trajectory in phase space parametrized by ζ that connects those two equilibria. From a point of view of dynamical systems we can investigate the character of the equilibria by linearizing the ODE (65) around these states, which yields in order to find out if an equilibrium is an attractor or a repeller. The Jacobians in (68) can be computed as and Evaluating them at the equilibria yields and Equation (68) has solutions of the form Y (ζ) = Y e σζ giving us the characteristic equation σ 2 − pσ = 0 with the coefficient being also a root Here, we introduced ∆(K z ) =ω ′ (K z ) − C as the difference of linear group and the envelope velocity. Equilibria are attractive if p < 0 and repellent if p > 0. In particular we can identify that p(K − z ) > 0 and p(K + z ) < 0 must be true if we want to associate K ± z with the asymptotic rest states at ζ = ±∞. The sign of p depends partly on the denominator ∆.
For downward-traveling wave packets, i.e. C < 0, the sign of ∆ is positive asω ′ (K ± z ) > 0 for K ± z < 0 which we anticipate here and explain subsequently. For upward-traveling wave packets, i.e. C > 0, we must make sure that ∆ does not assume a zero. Sinceω ′ has a global maximum at −K x / √ 2, it proved useful to choose By this choice ∆ is strictly negative. It can be also established by continuity that there must be aK z between K − z and K + z such that p(K z ) = 0. When the equilibrium passes this point, p switches the sign being determined by the numerator of (73) as the denominator, ∆, does not assume a zero.
To evaluate the trajectory numerically we use the explicit fourth-order Runge-Kutta method and rearrange (65) by applying the chain rule to By matrix inversion we obtain In particular det Df (Y ) = 0 gives a singularity if which parametrizes a curve in the phase space. Horizontal isoclines, which are curves in phase space where the vectors are horizontal or formally h = (h Kz , 0) T , are derived from (75) being Three vertical isocline where h = (0, h A ) T turn out to be vertical straight lines that cross the K z -axis at Prototypical phase portraits together with the trajectories computed numerically by the Runge-Kutta method are depicted in figure 5 for an upward and a downward-traveling wave packet, respectively. In terms of the previous considerations we can exploit phase plane analysis to show that (65) is indeed solved by trajectories of the form of heteroclinic orbits. In both plots the K z -axis consists of equilibria (grey line) and the A-axis is an isocline denoted by K iso z where the vectors are vertical. This implies that a trajectory initialized in the second quadrant, i.e. where K z < 0 and A > 0, stays in it. By the symmetry of the problem and the choice of the branch ofω, it is sufficient, without loss of generality, to focus on the second quadrant only.
Let us consider the upward-traveling wave (C > 0) which is shown in panel (a) of figure 5. Besides the aforementioned vertical isocline, we find a second one and two additional branches of the horizontal isocline A iso . These four curves delimit three regions outside of the singularity where the sign of the components of h are constant. In particular the region outside of the singularity between the two branches of A iso has vectors that possess only negative horizontal and negative vertical components. So, for every K + z on the K z -axis between the two branches of A iso we already know that by p(K + z ) < 0 we obtain an attractor which implies that there is a trajectory ending in the point (K + z , 0) for ζ → +∞ that we can trace backwards in ζ. This point is connected by the trajectory with some point on the right branch of A iso if K + z is far enough to the right such that it is clear from the singularity. All vectors, h, between the right branch of A iso and the A-axis point upwards and to the left. Hence, there must be a K − z in this region which is a repeller due to p(K − z ) > 0 such that it connects the crossing point on the right branch of A iso with (K − z , 0) when we keep tracing back. The resulting concatenated trajectory is then a heteroclinic orbit that we refer to as wave packet because the specific wave action density is of finite extent in contrast to the traveling wave front in the inviscid case.
Let us now consider the case C < 0 corresponding to downward-traveling wave packets. As depicted in figure 5 on the right panel we also find the two branches of A iso . However, the character of the equilibrium points between these two changes. In contrast to the upward-traveling wave, these equilibria are repellent, so p(K − z ) > 0. The components of h possess constant sign in the region between the two branches of A iso and therefore all vectors point upwards and to the left. A trajectory that starts at K − z sufficiently close to the left branch of A iso must cross it at some point when integrating forwards. All vectors to its left point downwards and hence the trajectory must land on an equilibrium point K + z giving us a heteroclinic orbit in accordance with a downward-traveling wave packet.
In conclusion, we found both up and downward-traveling wave packets solving the dissipative Grimshaw modulation equations. An illustration of these waves will be given in section 4.4 when we study the evolution of instabilities which is presented in the following section.

Stability of the dissipative Grimshaw traveling wave packets
In this section we will investigate the stability of the dissipative traveling wave packets. We linearize the modulation equations (60) around the traveling wave solutions in the translational coordinate system presented in (36) with the Jacobian matrices and Since the traveling wave solution Y (ζ) is stationary in the translational coordinates and the system is autonomous, equation (81) must have solutions of the form which yields an eigenvalue problem (L Y − λ)y = 0. The linear operator is given by We proceed with the same method being presented in section 3.2. As the traveling wave solution converges to its asymptotic rest states at the infinities, we can construct an asymptotic operator that exhibits the same essential spectrum as L Y itself.
The new eigenvalue problem (L ∞ − λ)y = 0 is rearranged to a first-order ODE, with the piecewise constant coefficient matrix and The purely imaginary eigenvalues of B ± (λ) give us the Fredholm borders as curves in the complex plane parametrized by µ ∈ R The Fredholm borders associated with λ ± 1 and λ ± 2 consist of the imaginary axis. They are invariant for every set of parameters defining the wave solution and can not cause any unstable essential spectrum.
However, the Fredholm borders associated with λ ± 3 have a real part in which we can identify p from (73). So, we substitute and see that and hence the stability of the wave packets depends on the sign of (91). In the previous section, we figured out that Therefore, the sign of (91) is determined by the sign of p at the asymptotic rest states, so the Fredholm border • λ + 3 is found on the right and λ − 3 on the left hand side of the imaginary axis for C > 0, • λ − 3 is found on the right and λ + 3 on the left hand side of the imaginary axis for C < 0.
To obtain the Fredholm indices in the delimited regions, we compute the Morse indices i ± for real valued λ ≫ 1, i.e. far to the right of the unstable λ ± 3 . They are the sum of the algebraic multiplicities of the unstable eigenvalues of the coefficient matrices B ± (λ). For upward-traveling wave packets we find i ± = 3 and for downward-traveling i ± = 0. Thus, according to (47) the Fredholm index in this region is zero and the Fredholm indices of the remaining regions are calculated in terms of the rule presented in section 3.2 for λ's crossing the Fredholm borders. In particular, the region left of the rightmost Fredholm border λ ± 3 has ind(L Y − λ) = +1 and hence belongs to the essential spectrum which is consequentially contained on the right hand side of the complex plane. To put it in a nutshell, the dissipative traveling wave packets are unconditionally unstable.
The essential spectrum for the upward-traveling wave packet is plotted in figure 6. It consists of a vertical band in the complex plane delimited by the two Fredholm borders, λ + 3 to the right and λ − 3 to the left. The remaining four borders lie on the imaginary axis. For the downward-traveling wave packets the essential spectrum is very similar. Only the orientation of the Fredholm borders switches and λ + 3 swaps places with λ − 3 . The Fredholm indices stay untouched.
We want to point out that (91) may give us a maximum instability growth rate which was a prerequisite for well-posedness as introduced in the beginning of this section. Therefore, it is worthwhile to investigate the spectrum and, by that, the nature of the instability in more detail.

Point spectrum and embedded eigenvalues
So far, we were only concerned with the essential spectrum, but to complete our argumentation about stability of the traveling wave packets we have to take the point spectrum also into account (cf. section 3.2). We want to show that (91) yields indeed a maximum instability growth rate, i.e. there exists no point spectrum to the right of the essential spectrum.
It suffices to consider the ODE associated with the eigenvalue problem for L Y given by (85). For upward-traveling waves in the limit ζ → +∞ the ODE has the coefficient matrix B + (λ) from (89) with eigenvalues ν + having only positive real part for λ's in the regions to the right of the essential spectrum. With the help of theorem 11.2 in Hartman [2002, p 301] using the eigenvalues ν + one can show that the corresponding "eigenfunctions" diverge at plus infinity and hence are not in L 2 . For downward-traveling waves the same argument holds at minus infinity. In conclusion, for both traveling wave solutions there is no unstable point spectrum and the rightmost Fredholm border gives a bound on the instability growth rate. Hence, the system was regularized in the desired sense by including dissipative effects.
In the remainder of this section we want to discuss how positive real valued λ's in the spectrum manifest themselves as instabilities. We observe that the Fredholm index in the unstable essential spectrum is +1 for both upward and downward-traveling wave packets. By the definition (42), a positive Fredholm index implies that the kernel of L Y − λ has at least dimension one and is in particular non-empty. Hence, there exist y ∈ L 2 , y = 0 that solve L Y y = λ y. Those y's can be called eigenfunctions and the corresponding λ's eigenvalues. In conclusion, the unstable essential spectrum consists of embedded eigenvalues and the instabilities are superpositions of exponentially growing eigenfunctions.
This observation is important as it prohibits convective instabilities [Sandstede and Scheel, 2000] which are characterized by an unstable essential spectrum that can be stabilized in an exponentially weighted space (instead of the usual L 2 ). Convective instabilities are transported on the domain towards the infinities away from the traveling wave solution such that the solution can survive. Even though it is possible here to stabilize the essential spectrum in our case, the interpretation that the instability may be convective is flawed due to the embedded eigenvalues.

Numerical simulation of the dissipative traveling wave packet
In this section we will analyze the essential instabilities of the dissipative traveling wave packets numerically. The modulation equations (56) are solved with the method presented in appendix A. Let us investigate the two classes of dissipative traveling wave solutions separately starting with the upward-traveling wave packet. The parameters defining the wave are given by N = 2.1, K x = 1.5, Λ = 1.1 and C = 1.6 . The initial conditions are computed by the fourth-order Runge-Kutta method solving the ODE (75). The results of the numerical time integration are shown in figure 7. The labeling of axes is the same as in section 3.3 for the numerical simulations of the inviscid waves. An instability is excited by an initial perturbation of Gaussian shape of the order of magnitude 10 −3 added to the wave packet. The instability growth quickly but, in contrast to the two inviscid cases of the previous sections, it possesses a large scale structure. In terms of (91) we can compute the maximum instability growth rate analytically and find ℜ(λ) max = 5.1 . By computing the L 2 -norm of the difference between the numerical and a reference solution, the growth rate of the actually excited instability can be computed numerically. It is approximately 3.8 . This rate is smaller than the expected value because the initial perturbation was not optimal and the numerical error is of diffusive character. Let us continue with the dissipative downward-traveling wave packet being defined through the parameters N = 2.2, K x = 1.3, Λ = 1.1 and C = −0.7 . Likewise, it is initially perturbed by a Gaussian 10 −3 -deviation from the exact traveling wave solution. We plotted the results of the numerical time integration in figure 8 where we can observe that the wave blows up by a largescale instability. The maximum growth rate according to (91) is 8.9 . And the approximated numerical growth rate of the actually excited instability is 6.8 . In conclusion, the numerical results fit nicely into our theoretical framework from the spectral stability analysis. Note that the linear group velocity of the downward-traveling wave packet is positive everywhere whereas the envelope velocity is negative. This seemingly paradoxical situation is a logical consequence of the nonlinear nature of these waves.

Conclusion
We studied the modulation equations of two-dimensional nonlinear gravity waves in the isothermal atmosphere and their stability. The waves were assumed to be horizontally homogeneous. Two different asymptotic regimes of the Euler equations were of interest, the Boussinesq and the pseudoincompressible equations. The former does not account for altitudinal variations of the background density but possesses the convenient property of having plane wave solutions that we examined with respect to stability introducing the notion of spectral stability analysis. Here, we were able to reproduce the criterion which is known from modulational instability of the weakly nonlinear theory. Allowing waves to propagate on times scales where the background density changes, brought us to the inviscid pseudoincompressible modulation equations. They exhibit solutions of the form of upward-traveling wave fronts which destabilize unconditionally.
The unstable operator spectra of the linearized inviscid modulation equations of the two regimes were revealed to be unphysical. Mathematically speaking, the systems were ill-posed. Therefore, we deduced that they may be regularized by considering the compressible Navier-Stokes instead of the Euler equations. We found that the emerging dissipative modulation equations can have upward and downward-traveling wave packet solutions. The spectral stability analysis showed that the wave packets are unstable due to embedded eigenvalues.
We emphasized that linear group and envelope velocities can differ for nonlinear waves. Remarkably, the dissipative downward-traveling wave packet is an especially interesting case. Its phase and envelope velocities are negative, though the linear group velocity is positive everywhere.