On strongly nonlinear gravity waves in a vertically sheared atmosphere, Part I: Spectral stability of the refracted wave

We investigate strongly nonlinear stationary gravity waves which experience refraction due to a thin vertical shear layer of horizontal background wind. The velocity amplitude of the waves is of the same order of magnitude as the background flow and hence the self-induced mean flow alters the modulation properties to leading order. In this theoretical study, we show that the stability of such a refracted wave depends on the classical modulation stability criterion for each individual layer, above and below the shearing. Additionally, the stability is conditioned by novel instability criteria providing bounds on the mean-flow horizontal wind and the amplitude of the wave. A necessary condition for instability is that the mean-flow horizontal wind in the upper layer is stronger than the wind in the lower layer.


Introduction
The importance of gravity waves for the atmospheric dynamics and hence for weather and climate forecasting was first established in [6] and [17]. Gravity waves redistribute energy vertically and thereby couple the different layers of the atmosphere [9,3]. Usually excited in the troposphere, gravity waves may persist deep into the upper atmospheric layers [10,11,12] where they interact with the mean flow. They exert drag onto the horizontal mean-flow momentum, produce heat when dissipating [2], and cause increased mixing of tracer constituents such as green-house gases [21]. To this day, many questions regarding the sources, propagation and dissipation of gravity waves are still to be answered.
From a conceptional point of view, three distinct types of theories-each with its own benefits and challenges-may be employed to study these questions. First and foremost, linear wave theory provides the foundation for our understanding of gravity waves. It predicts dispersion, refraction, reflection, wave packet formation etc. as long as the amplitudes of the waves can be considered to be infinitesimally small [5].
For large-amplitude gravity waves, perturbation theory is a versatile and powerful tool leading to weakly nonlinear wave theory. The key idea of weakly nonlinear theory is to expand the linear solution to the fluid dynamical equations in an asymptotic series assuming small but finite amplitude. Corrections of the order of the amplitude squared to the linear model add new phenomena to wave theory such as wave-mean-flow interaction and modulational instabilities [27,15,24,25,26]. For internal gravity waves in the oceans, weakly nonlinear theory is almost always sufficient. However, this is often not the case for atmospheric gravity waves. When excited over mountains, tropospheric gravity waves may have amplitudes of the order of magnitude of the background wind. In the stratosphere and higher, such strong amplitudes are rather the rule than the exception. This effect is due to anelastic amplification which is absent in the ocean and caused by the compressibility of air. The majority of atmospheric gravity waves is excited in the troposphere. When they extend into higher altitudes, they encounter exponentially decreasing background density. It is then energy conservation that lets the waves grow exponentially causing strong amplitudes.
The mathematical description of gravity waves beyond weakly nonlinear theory was pioneered by Grimshaw [13,14] and further analyzed in [1] and also [22] in terms of modulation equations. In this asymptotic theory the small expansion parameter is the ratio between the scale on which the background changes and the wavelength. There is, however, no further restriction on the velocity amplitude and therefore it is allowed to be of the same order as the background flow in contrast to weakly nonlinear theory where the amplitude itself is the expansion parameter. Hence, we may call Grimshaw's modulation equations a strongly nonlinear wave theory.
The notion of "strongly nonlinear" is stipulated in the following sense. A superposition of solutions from the class of strongly nonlinear waves does not lead to small errors when inserted into the governing equations (the Euler or Navier-Stokes equations) as it is the case for the weakly nonlinear theory. Instead such a superposition generates terms by wave-wave interaction of order unity in the asymptotic limit. However, this notion does not necessarily imply that the nonlinearities in the governing equation, i.e. the advection terms, are bigger than the pressure gradient term, say, in terms of nondimensional analysis. In fact, the asymptotic expansion of Grimshaw's theory is such that the perturbation advection terms vanish due to the solenoidality of the wind field. Some theoretical investigations on strongly nonlinear effects and their implications with respect to observations and modeling were performed in [23,20].
The objective of this paper is to investigate the interaction of strongly nonlinear waves with a very thin shear layer. This renders a common situation in the real atmosphere, e.g., when a mountain wave meets the tropospheric jet and gets refracted [7]. In the setting of linear theory, this situation was studied in the seminal work [8] using a layered model. They approximated the height-dependent background by piecewise constant functions. This idea was advanced in [18] for linear wave packets. We will adopt the idea of a layered background atmosphere but for strongly nonlinear waves.
Our main results are necessary and sufficient criteria for the stability of a refracted non-hydrostatic wave in a two-layered background. If the mean-flow horizontal wind in the lower layer is stronger than the wind in the upper layer, than the wave is stable with respect to the modulation equations. If, however, the upper-layer mean-flow horizontal wind is stronger, which typically occurs for a gravity wave entering the jet, the refracted wave becomes unstable if both mean-flow horizontal winds meet particular upper bounds and if the amplitude of the wave is sufficiently large.
This work is structured as follows. In section 2, we will introduce Grimshaw's modulation equations as our governing equations. Section 3 will be dedicated to the strongly nonlinear wave refracted at a discontinuity as a particular solution to the modulation equations. The stability of the refracted wave will be analyzed in section 4. We will summarize the results in section 5 and offer some concluding remarks in section 6.
In the second part, we will investigate the wave solution and the stability results numerically by means of Large Eddy Simulations.

The model equations
We consider vertically modulated, horizontally periodic, non-hydrostatic, and strongly nonlinear gravity waves in the unbounded x-z-plane which we model by Grimshaw's modulation equations [14,1,22], This coupled set of equations governs the evolution of vertical wavenumber k z , wave action density ρa, and mean-flow horizontal wind u. The horizontal wavenumber k x = K x is a constant being, without loss of generality, positive. The Extrinsic frequency ω =ω + k x u is determined by the Doppler-shifted intrinsic frequency which depends on wavenumber due to the dispersion relation for non-hydrostatic waves, Its derivative with respect to the vertical wavenumber, represents the vertical linear group velocity. Two coefficients appear representing the background state of the atmosphere, the Brunt-Väisälä frequency N and the density ρ In Whitham's modulation theory [27], equation (1a) represents conservation of waves. The second equation (1b) yields conservation of wave action. Finally, (1c) describes the acceleration of mean-flow horizontal momentum due to the convergence of the flux of horizontal pseudo-momentum k x ρa.

The refracted wave solution
In this section, we construct an analytical solution to the modulation equations reminiscent of a typical mountain wave that encounters a jet at a certain height-e.g. the tropospheric or mesospheric jet-and experiences refraction. This is a classical situation from linear wave theory. In contrast to linear theory however, we will see that the wave drags the background flow by direct nonlinear interaction resulting in a self-induced mean flow and Doppler shifting.
When multiplying k x to (1b), then subtracting (1c) and integrating we find We also findŪ as an integration constant that can be identified as the horizontal background wind, the wind without the occurrence of a wave (a = 0), and u is really the mean-flow horizontal wind. Let us consider a vertically sheared background flow being piecewise constant with the discontinuity at z = 0,Ū (z) = Ū 1 , z < 0, Then a stationary solution is given by The mean-flow horizontal wind really becomes a diagnostic variable due to (5), Here and further on in this investigation, j = 1 denotes the layer below and j = 2 the layer above the discontinuity at z = 0. The constant pieces are determined by integration of (1a) and (1b),ω The integration constant in (9a) must be zero as the extrinsic frequency is readily the derivative of the phase function with respect to time. Therefore, a stationary phase requires vanishing frequency. By means of this argument and (2), (9a) can be rearranged to obtain Equation (9b) provides an interesting physical implication: the vertical wave action flux is invariant throughout the layered atmosphere. And in particular, it is constant crossing the interface of the two layers. In order to ensure internal waves we must find Otherwise the discriminant in (10) is negative corresponding to evanescent or rather external waves. We have introduced a new parameter J j which will become important in the following analysis. Please note that the mean-flow horizontal wind u has become a diagnostic variable due to (5). Let us concatenate the prognostic variables of the wave solution into a vector, A sketch of the wave which would be generated by this solution is plotted in figure  1. Due to the larger background velocity in the upper layer, the phase lines are more inclined than in the lower layer.

Modulational stability of the refracted wave
Linearizing the governing equations reduced by (5) around the stationary solution P P P combined with an ansatz for the perturbation p p p(z, t) =p p p(z)e λt , p p p(z, t) = (k z , a) T yields an eigenvalue problem (L − λ)p p p = 0 for a linear differential operator L which can be reformulated as an ordinary differential equation dp p p(z) dz = B(z, λ)p p p(z).
The coefficient matrix contains the refracted wave solution and is therefore piecewise constant in z, where To decide about the stability of our refracted wave solution, we will investigate the spectrum of the linear operator posed by the linearized governing equations. Roughly speaking, the spectrum is the set of all λ ∈ C for which our operator is not invertible. If the spectrum is contained on the left hand side of the complex plane, i.e. all ℜ(λ) ≤ 0, then we can say the solution is spectrally stable and unstable otherwise. The spectrum consists of two qualitatively different subsets defined by its Fredholm properties: the essential (continuous) spectrum and the point (matrix-like) spectrum. We elaborated on this definition in terms of Fredholm operators in the context of strongly nonlinear gravity waves in [23]. For the sake of brevity, the details are not repeated in this paper. We refer to the book of [16] for a detailed introduction and precise definitions.

The essential spectrum
The essential spectrum is given by the hyperbolicity of the matrix (16) [19,16]. A matrix is said to be hyperbolic if all its eigenvalues have non-zero real part. The boundaries of the open regions in the complex plane which comprise the essential spectrum coincide exactly with (16) being not hyperbolic. Therefore, the ansatz leads to which represents four parameterized curves in the complex plane enclosing the essential spectrum. In particular, we will find unstable essential spectrum if one of the boundary curves enters into the right half of the complex plane as either left or right from the curve we will find essential spectrum. This occurs when H j < 0 which corresponds to the classical modulational instability criterion [15,23,20]. The instability growth rate is then K x A j |H j ||µ| which happens to be unbounded in µ or in other words arbitrarily large µ results in arbitrarily large instability growth which would render the problem unphysical and mathematically ill-posed. The classical modulational instability criterion can be rewritten by insertion of (10) into (18) as providing us with a new lower bound greater than the criterion for non-evanescence in (11).

The point spectrum
In the following, we assume a stable essential spectrum, so H 1 > 0 and H 2 > 0. The point spectrum consists of eigenvalues to which the corresponding eigenfunctions belong. Eigenfunctions must be of the form p p p = p p p 1 , z < 0 p p p 2 , z > 0 ,p p p j =p p p j e σ j z + cc (22) with cc denoting the complex conjugate. We obtain the spatial eigenvalues σ j by the solvability condition of the emerging algebraic system when inserting the ansatz (22) into (14), Please note that ± indicates two linearly independent solutions. Special care has to be taken for ±K x A j H j = C j . In this very particular case, λ = 0 is the only eigenvalue. Associated basis vectors to the eigenvalues (23) may be written as In terms of the basis vectors, the eigenvectors are given by superposition, i.e. linear combination,p To be an eigenfunction at least one of the coefficients b ± j must differ from zero. Because if all coefficients vanish, then the solution is trivial and not an eigenfunction. The eigenfunctions must obey three constraints. 0 z Figure 2: Sketch of the perturbation 1. In the limit z → −∞, the eigenfunctions must vanish.
3. At the discontinuity z = 0, the eigenfunctions need to be continuous.
These three constraints constitute perturbations of finite energy as well as continuous pressure and vertical wind velocity at the discontinuity of the backround horizontal wind. Figure 2 presents a sketch of the form of the eigenfunctions in line with the constraints 1-3. They oscillate and decay away from the jump.
Let us fix in the following ℜ(λ) > 0 and by this we ask for unstable eigenvalues. We consider (23) and distinguish the four different cases that may occur regarding its sign.
• If K x √ A 2 H 2 < C 2 , then ℜ(σ + 2 ) < 0 and ℜ(σ − 2 ) < 0. Hence, constraint 2 allows for b + 2 = 0 and b − 2 = 0. Next, we investigate the continuity constraint 3, i.e.p p p 1 =p p p 2 at z = 0 which yields Note that we already anticipated that in all cases b − 1 = 0. Given b + 1 we can solve (26) for We read immediately from (27) according to constraint 1. The latter renders therefore a sufficient condition for stability as, if it is true, there exist only trivial solution for any ℜ(λ) > 0 and hence it cannot be an eigenvalue. Furthermore, b + 2 = 0 fulfilling constraint 2 is only possible if b + 1 = 0 due to (27a), but then (27b) yields also b − 2 = 0. And again we are left with the trivial solution. By the same argument, K x √ A 2 H 2 > C 2 is hence a sufficient condition for stability. In conclusion, there are nontrivial solutions, which happen to be unstable, if

Summary
In this section, we summarize and reformulate the results of the previous sections. We may rewrite our criterion from the point spectrum in terms of the saturation amplitudes of the waves in each layer α j which relates to the wave action density by Waves become unstable due to static instabilities if α j ≥ 1 [4]. Thus, by inserting (29) and (11) into the criteria (28) and rearranging, the refracted wave is stationary stable but unstable due to perturbations from the point spectrum if To have unity as the upper bound for the right hand side in (30a) implies that J 1 ∈ (2, ∞). It represents a particularly narrowing bound on the amplitude in the lower layer as, The minimum is assumed at J 1 = 3, so at best 1 > α 1 > 8/9 ≈ 0.94. Zero as the lower bound for the left hand side in (30b) is met if J > 3/2 which is the same bound from the stable essential spectrum. Utilizing the fact that the vertical wave action flux is invariant across the interface (cf. (9b)), we find C 1 A 1 = C 2 A 2 . We combine our finding with (29) and obtain which we substitute in (30). The resulting inequality reads being the same as g(J 2 ) > g(J 1 ). Notice that the function g is continuous as well as monotonically decreasing in (3/2, ∞) and hence which implies due to (11) that |U 1 | < |U 2 |.
In conclusion, a refracted gravity wave, that is statically stable, non-evanescent, nonhydrostatic and features a stable essential spectrum of the linearized modulation equations, becomes unstable due to perturbations from the point spectrum if and

Discussion
On the one hand, the constraint on the saturation amplitude in the lower layer for an emerging instability is extreme. It is at best 94 % or higher. On the other hand, the conditions for the mean-flow horizontal wind favorable for instability are fairly common. A typical situation where this type of instability may occur is the lower edge of the atmospheric jets where a mountain wave enters a strong shear zone from below and gets refracted into the jet. It is not unlikely that such a wave has large amplitude due to anelastic amplification. The perturbation affects not only the wave properties but also the mean-flow horizontal wind as a diagnostic variable due to (5). Also we observe that perturbations decay exponentially away from the discontinuity of the background wind (cf. figure 2 for illustration). We conclude therefore that a point-spectrum instability is localized at the jump. Due to these observations we can suspect an attenuation of the sharp gradient in the background horizontal wind by an instability from the point spectrum. In conclusion, strongly nonlinear gravity waves may serve as a mitigation of strongly sheared flows.
Moreover, Grimshaw's modulation equations, which provide the basis for our analysis, were derived from the Euler equations using WKB theory. A crucial assumption in the derivation is that the mean-flow varies on a larger scale than a typical wavelength. By assuming vertically piecewise constant mean-flow horizontal wind, we potentially infringe on this WKB constraint. The range of validity of our results will be tested by means of numerical simulations of the Euler equations in a successive paper.