Schlutowetal 2017

This draft was prepared using the LaTeX style file belonging to the Journal of Fluid Mechanics 1
Finite-amplitude gravity waves in the

atmosphere: traveling wave solutions
Mark Schlutow1 †, R. Klein1 and U. Achatz2
1
Institut für Mathematik, Freie Universität Berlin, Arnimallee 6, 14195 Berlin, Germany
2
Institut für Atmosphäre und Umwelt, Goethe-Universität Frankfurt, Altenhöferallee 1, 60438
Frankfurt am Main, Germany
(Received xx; revised xx; accepted xx)
Wentzel-Kramers-Brillouin (WKB) theory was employed by Grimshaw (Geophys. Fluid

Dyn., 6, 131–148, 1974) and Achatz et al. (JFM, 210, 120–147, 2010) to derive mod-
ulation equations for non-hydrostatic internal gravity wave packets in the atmosphere.
This theory allows for wave packet envelopes with vertical extent comparable to the
pressure scale height and for large wave amplitudes with wave-induced mean flow speeds
comparable to the local fluctuation velocities.
Two classes of exact traveling wave solutions to these nonlinear modulation equations are
derived here. The first class involves horizontally propagating wave packets superimposed
over rather general background states. In a co-moving frame of reference, examples
from this class have a structure akin to stationary mountain lee waves. Numerical
simulations corroborate the existence of close-by traveling wave solutions under the
pseudo-incompressible model and reveal better than expected convergence with respect
to the asymptotic expansion parameter.
Traveling wave solutions of the second class also feature a vertical component of their
group velocity but exist under isothermal background stratification only. These waves in-
clude an interesting nonlinear wave-meanflow interaction process: A horizontally periodic
wave packet propagates vertically while draining energy from the mean wind aloft. In
the process it decelerates the lower level wind. It is shown that the modulation equations
equally apply for hydrostatic waves in the limit of large horizontal wavelengths.
Aside from these results of direct physical interest, the new nonlinear traveling wave
solutions provide a firm basis for subsequent studies of nonlinear internal wave instability
and for the design of subtle test cases for numerical flow solvers.
Key words:
1. Introduction
The stratification of the earth’s atmosphere gives rise to internal gravity waves which
obtain kinetic energy and momentum from the mean flow when excited, transport them
vertically as well as horizontally, and redistribute them by non-linear interaction and
breaking. By these processes, gravity waves have a major impact on the large-scale
middle atmospheric circulation (Fritts 2003; Becker 2012). In fact, when wave packets
propagate upwards in a density-stratified medium they amplify and ultimately tend to
break, generate turbulence, and transfer heat and momentum to the mean flow. At the
† Email address for correspondence: mark.schlutow@fu-berlin.de
2 M. Schlutow, R. Klein and U. Achatz
same time, internal waves appear on a considerable spectrum of wavelengths, so that even
high resolution numerical models can only resolve the long-wave part of this spectrum
whereas the small-scale part together with the effects of wave breaking and turbulence
need to be parametrized (McLandress 1998; Alexander & Dunkerton 1999).
In developing such parameterizations it is observed that wave breaking and the sub-
sequently developing turbulence are the result of wave instabilities. Numerous different
instability mechanisms are known and provide wave breaking criteria in gravity wave
parametrizations. Specifically, these schemes aim to predict those positions where wave
drag is produced and affects the mean flow, e.g.: if the perturbation of the background
state due to a wave leads to unstable stratification locally, then the wave is said to be
statically unstable. Yet, even statically stable waves can generate strongly growing modes
and become unstable subsequently if their local Richardson number is less than a quarter
(Lelong & Dunkerton 1998; Achatz 2007; Liu et al. 2010).
By means of Floquet theory, which is applicable for linear partial differential equa-
tions with periodic coefficients, parametric instabilities were considered by Mied (1976),
Lombard & Riley (1996) and others. These authors exploited that monochromatic plane
waves in a uniformly stratified background are exact solutions of the Boussinesq equations
that are nonlinear in a particular way. More recently, effects of modulational instability
in a weakly nonlinear regime, in which wave-induced mean flows are of the order of the
wave-amplitude squared, were analyzed by Sutherland (2001, 2006) and Tabaei & Akylas
(2007).
The underlying basis for analytical stability studies are traveling wave solutions. They
are defined as stationary wave-like solutions in a translational coordinate system moving
with a constant relative velocity. The aforementioned plane waves are prime examples of
this kind. Once traveling wave solutions are found, linearization of the chosen fluid flow
equation system relative to these stationary solutions generates an eigenvalue problem.
If the spectrum of the resulting linear operator exhibits a positive real part, then the
traveling wave solution is linearly unstable.
All above-mentioned stability analyses assume a Boussinesq-type fluid. While this is
appropriate for small-scale flows, it is of limited value for flows covering scales comparable
to the pressure scale height or larger. Thus, Klein et al. (2010); Klein (2011); Achatz et al.
(2010) discuss the validity of various flow models for the description of internal wave
dynamics in comparison with the full compressible Euler equations. They found that,
in this sequence, the Boussinesq, anelastic, pseudo-incompressible, and full compressible
flow equations have an increasing range of validity.
In particular, Achatz et al. (2010) presented a theory that simultaneously covers (i)
finite amplitude, (ii) nonlinear wave packets with (iii) wave envelopes extending over
a pressure scale height or more in (iv) strong, non-uniform stratification. This is of
interest, e.g., in studying breaking waves in the stratosphere. The theory is derived
from the full compressible flow equations through a multiple-scales approach combined
with Wentzel-Kramer-Brillouin (WKB) theory that assumes slow variation of the wave
envelope and phase. Out of the family of sound-proof models (Boussinesq, anelastic,
pseudo-incompressible), only the latter covers this and all the less demanding regimes
and yields results in all of them that are identical to leading order with those obtained
from the full compressible model. In this sense, this theory is, to the best of the authors’
knowledge, the most general description of internal wave packets in the atmosphere
currently available, and is thus adopted as the basis for the present work. The theory
was tested favourably against numerical simulations by Rieper et al. (2013a).
The main motivation for the present work is to prepare the ground for internal wave
stability analysis of internal wave packets in a deep atmosphere covering at least a
Finite-amplitude gravity waves in the atmosphere: traveling wave solutions 3
pressure scale height and for the development of related subgrid scale parameterizations.
To this end we derive traveling wave solutions of the large scale WKB equations derived
by Achatz et al. (2010) as references for linearizations. Being exact solutions of a nonlinear
flow model valid in a physically relevant regime, however, these solutions may also be of
interest independently of the stability question: The horizontally traveling wave solutions
strongly resemble mountain lee waves in their co-moving frame of reference, so that the
present theory may provide a new analytical handle to study nonlinear properties of
the latter. The class of vertically traveling wave solutions includes narrowly supported
structures that redistribute the kinetic energy of the background flow to produce propa-
gating shock-like transitions with strong vertical shear somewhat reminiscent of the wave
patterns that induce the QBO, albeit on much larger horizontal scales (Baldwin et al.
2001). Our main focus here is on the small-scale part of the spectrum, however, so that
we defer more detailed analyses of larger scale flows by related methods to future work.
This paper is structered as follows. In §2 we will revisit the scaling assumptions of
Achatz et al. (2010) to derive the governing equations. A spectral WKB ansatz for
finite-amplitude waves will be exploited in §3. We will examine, in addition to that,
the existence of weak asymptotic solutions. This will lift the already well established
theory on a stronger mathematical basis. Solvability constraints for the ansatz functions
lead to the modulation equations predicting the evolution of the wave properties. In §4 we
will investigate the energy budget of the wave-mean-flow interaction and show that the
modulation equations conserve the sum of wave and mean-flow energy which allows to
compute an energy exchange rate. §5 contains our main results. Here, we derive traveling
wave solutions for the modulation equations by means of a translational coordinate
transformation. These traveling waves give us leading order solutions for the governing
equations. We will find two distinct classes of traveling wave solutions. For the first one,
we will present an analytic solution and for the second one we can show the structure
of the wave and for which initial conditions and set of parameters a solution exists.
The exact traveling wave solutions are asymptotic solutions to the fully nonlinear Euler
equations. We will test the validity of the solutions with respect to the Euler dynamics
numerically in §6. In §7 we will show consistency of hydrostatic gravity waves with the
modulation equations. The paper will be concluded with some final remarks in §8.
2. The governing equations

In this section we will present a derivation of equations, that describes the dynamics of
saturated gravity waves in a hydrostatic background flow, starting from the dimensionless
Euler equations.
2.1. Scaling the Euler equations for large amplitude internal waves
We consider gravity waves in a two-dimensional domain neglecting the Coriolis force.
Viscous effects as well as external heating sources are ignored, too. With respect to these
basic assumptions the time evolution of the flow is governed by the compressible Euler
equations for an ideal gas which may be written in non-dimensionalized form as
1−κ 1 1
0 = Dt v + θ∇π + 2 ez
κ Ma2 Fr
0 = Dt θ
κ
0 = Dt π + π∇·v (2.1)
1−κ
Dimensionless variable x z t v T p
Reference value Lref Lref Lref /vref vref Tref pref
Dimension in SI units m m s ms−1 K Pa
Table 1. Non-dimensionalization of variables
with (x, z) the horizontal and vertical coordinate and the corresponding unit vectors
ex and ez . Here, Dt = ∂t + v·∇ is the material derivative. The velocity vector field is
denoted by v. The Exner pressure is related to the usual pressure by
π = pκ (2.2)
where κ = 2/7 is the ratio of the ideal gas constant R and the specific heat capacity at
constant pressure cp . From the temperature T we define the potential temperature to be
θ = T /π. (2.3)
The independent and dependent variables are non-dimensionalized according√to table 1.
The Mach and the Froude number are given by Ma = vref /Cs and Fr = vref / gLref with
r
1
Cs = RTref (2.4)
1−κ
the speed of sound and g the gravitational acceleration. Let Hp , Hθ , Lref , and Nref denote
the pressure and potential temperature scale height, a typical wavelength, and reference
Brunt-Väisälä frequency, respectively. The distinguished asymptotic limit adopted by
Grimshaw (1974) and Achatz et al. (2010) is reproduced by assuming:
(i)low Mach number

r
κ
Ma = ε 1 , (2.5)
1−κ
(ii)leading-order stable stratification
Hp dΘ RTref
= O(1) (ε → 0) where Hp = ≈ 10 km, (2.6)
Θ dz g
(iii)non-hydrostatic motions
xref = zref ≡ Lref ≈ 4 km, (2.7)
(iv)large internal wave Froude number
s
vref g dΘ
Frint = = O(1) (ε → 0) where Nref = ≈ 0.02 Hz. (2.8)
Nref Lref Θ dz
On account of (2.5),(2.6), and (2.8) we find that Lref /Hp = O(ε) and we let
Lref = εHθ where Hθ = Hp /κ ≈ 40 km (2.9)
for specificity. Combining the scaling assumptions (i)–(iv), we obtain the distinguished
limit Fr = ε1/2 which links the Froude to the Mach number, so that ε is the only
parameter appearing in (2.1).
2.2. Asymptotic scalings
Here we prepare the ground for the subsequent summary of the WKB analysis in
section 3 by anticipating the essential spatio-temporal and amplitude scalings of the
dependent variables from Achatz et al. (2010). Let
(X, Z, T ) = (εx, εz, εt) (2.10)
denote compressed coordinates that resolve spatial scales comparable with the pressure
scale height and the associated advection time scale. Furthermore, for any variable V ,
let
V0 (Z) and Vε(i) (x, z, t; ε) (2.11)
denote its slowly varying leading-order horizontal background and the collection of its
ith and higher order contributions in ε with wave-induced variations on the fast and slow
scales. Then the asymptotic scaling of the WKB solutions to be discussed below can be
summarized as
(0)
 
  uε (x, z, t; ε)
v (0)
wε (x, z, t; ε)
 
θ  =  (1) .

(2.12)
 θ0 (Z) + εθε (x, z, t; ε) 
π (2)
π0 (Z) + ε2 πε (x, z, t; ε)
In the velocity field the wave appears at O(1), whereas potential temperature and Exner
pressure decompose into leading-order hydrostatic background variations θ0 , π0 satisfying
dπ0 1
=− , (2.13)
dZ θ0
and higher-order fluctuations. Potential temperature perturbations are of order O(ε),
whereas Exner pressure deviations appear first at O(ε2 ). See (Achatz et al. 2010) and
section 3 below for more detail.
For later reference, we use the ideal gas equation of state to express the leading order
density and temperature as
(1−κ)/κ
π0
ρ0 = (2.14)
θ0
and
T0 = θ0 π0 , (2.15)
respectively. In addition to that we also define the “logarithmic background derivative”

of any f0 ∈ {θ0 , π0 , ρ0 , T0 } by
1 df0 d ln(f0 )
ηf = = . (2.16)
f0 dZ dZ
Taking the derivative w.r.t. Z of (2.14) and substituting with (2.16) we find that
1−κ
ηπ = ηρ + ηθ . (2.17)
κ
We can identify ηθ = N 2 to be the local Brunt-Väisälä frequency squared, −ηρ−1 may be
interpreted as the local density scale height. Combining (2.13) – (2.17) it turns out that
all background variables can be expressed by the background temperature,
1 1
ηρ = −ηT − and N 2 = ηT + . (2.18)
κT0 T0
So, once the vertical background temperature profile is known, all profiles of the remain-
ing background variables are determined. If we now inject the distinguished limit (2.9) as
well as the multiple-scale representation from (2.12) into the Euler equations (2.1) and
substitute by (2.17), we arrive at a system of evolution equations
Dt u + ∂x P = −ε N B∂x P
Dt w + ∂z P − N B = −ε (N B∂z P − N 2 P ) + ε2 N 3 BP
(2.19)
Dt B + N w = −ε (N 2 + ηN )wB
∂x u + ∂z w = −ε (N 2 + ηρ )w + C
for the new set of prognostic variables
!T
(1)
T 1 θε
U = (u, w, B, P ) := u(0)
ε , wε(0) , , θ0 πε(2) (2.20)
N θ0
where B and P relate to the buoyancy and the kinetic pressure, respectively.
The expression

κ
C = ε2 (N 2 + ηρ ) Dt P + P (∂x u + ∂z w) − ε3 (N 2 + ηρ )N 2 wP (2.21)
1−κ
in (2.19) denotes higher-order terms which do not contribute to the leading-order WKB
solutions because either they are literally of higher order and just do not appear, or
they do appear in the expanded equations up to the order considered but cancel each
other. We note in passing that (2.19) without the C-terms constitutes Durran’s pseudo-
incompressible model equations, and this corroborates the conclusion by Achatz et al.
(2010) that the full compressible and pseudo-incompressible models share the WKB-wave
packet solutions in the considered flow regime. We will refer to system (2.19) henceforth
as scaled Euler equations.
3. WKB theory for finite-amplitude gravity waves with non-uniform

stratification
In this section we will apply a weak asymptotic scheme to derive the modulation
equations whose traveling wave solutions are the main objective of this survey.
3.1. A spectral WKB expansion

The scaled Euler equations are a set of nonlinear, first-order PDEs whose coefficients
depend only on the stretched coordinate Z due to the background stratification. To
account for the nonlinearities we exploit an asymptotic spectral expansion approach
M
X X −1
U= εn Ûn,m eimε Φ
+ o(εM ) (3.1)
n=0 |m|6n+1
for the vector (2.20). The subscript n keeps track of the order in ε while the m-subscript
carries the order of the harmonics. So, resonant coupling of the wave with the mean
flow, i.e. zeroth harmonics, occurs already in leading order (Chu & Mei 1970; Miura
∗
& Kruskal 1974). We restrict Ûn,−m = Ûn,m for real-valuedness, where ∗ indicates the
complex conjugate of a variable. To account for the slowly varying background state we
apply the WKB assumption that the coefficients and the phase are functions of only the
slow variables (X, Z, T ), such that
(X, Z, T ) 7→ Ûn,m (X, Z, T ), Φ(X, Z, T ). (3.2)
The differential operators become
∂t = ε∂T and ∇ = ε∇ε = ε(∂X , ∂Z )T . (3.3)
The phase should be also expanded in ε as well, like
Φ = Φ(0) + ε2 Φ(2) + O(ε4 ). (3.4)
We need only the even gauge functions because effects by any Φ(1) are already captured
by arg(Û0,1 ). Similarly, higher order odd phases are as well included in the corresponding
coefficients. Since we will only analyze the system up to M = 1, these higher order effects
do not matter in this work, so we keep simply Φ ≡ Φ(0) . By means of (3.3) we can define
the local wave vector and the frequency,
ω = −∂T Φ, (3.5)
k = ∇ε Φ. (3.6)
By inserting the ansatz (3.1) into the governing equations (2.19), we obtain an equation
of the form
M
X X −1
εk Ck,l (X, Z, T )eilε Φ(X,Z,T )
+ o(εM ) = 0. (3.7)
k=0 |l|6k+2
The coefficients Ck,l represent the terms which contain the ansatz functions Ûn,m , Φ,
their derivatives, and the background variables. One can now establish that in the
limit ε → 0 the coefficients must vanish independently of the order in ε and the order
of the harmonics, i.e. Ck,l = 0 for all possible k and l. The result is an equation
hierarchy providing constraints on the ansatz functions which will give us the modulation
equations. In order to achieve independence with regard to the harmonics Whitham
(1965); Grimshaw (1974); Miura & Kruskal (1974); Tabaei & Akylas (2007) assumed a
priori that the solution vector U is additionally to its spatio-temporal dependence also
a periodic function of the fast phase ε−1 Φ, so
U (ε−1 Φ, X, Z, T ; ε) = U (ε−1 Φ + 2π, X, Z, T ; ε). (3.8)
If one now averages (3.7) with respect to the fast phase over a period 2π, then one can
exploit that the harmonics are mutually orthogonal. In terms of this argument it was
justified that the claim of independently vanishing coefficients must hold with regard to
the average.
We encounter two issues, when we want to apply this technique to our approach in order
to establish independently vanishing coefficients: (i) The fast phase depends on ε and
(ii) it also depends itself on space and time. Due to these properties the assumption of
periodicity must be justified in terms of the asymptotic series expansion of U . Miura
& Kruskal (1974), who applied the nonlinear WKB method to the Korteweg-DeVries
equation, wrote in this regard: “it remains unsettled whether or in what sense the series U
[...] approximates the solution of the original problem.” However, we show in appendix A
an alternative approach where periodicity is not needed by introducing the notion of
weak asymptotic solutions based on Danilov et al. (2003). It requires only that the ansatz
functions are twice differentiable and it is hence more general than the averaging method.
The key ideas to overcome the aforementioned issues are (i) to formulate (3.7) in a weak
sense and (ii) to substitute the Cartesian coordinates by a map onto some phase following
coordinates.
3.2. Leading-order analysis

In the following sections we present the equation hierarchy obtained by the weak
asymptotic analysis. From C0,0 = 0 we find solenoidality of the velocity field of the
primary harmonics,
ikx û0,1 + ikz ŵ0,1 = 0. (3.9)
And constraints for some mean-flow components (the zeroth harmonics)
B̂0,0 = ŵ0,0 = 0. (3.10)
From C0,1 = 0 we get a linear homogeneous system of equations for the primary harmonics
M(−iω̂, ik)Û0,1 = 0 (3.11)

with the anti-Hermitian coefficient matrix
 
−iω̂ 0 0 ikx
 0 −iω̂ −N ikz 
M(−iω̂, ik) =  . (3.12)
 0 N −iω̂ 0 
ikx ikz 0 0
Here, ω̂ = ω − û0,0 kx represents the intrinsic frequency that is the frequency one observes
in a reference frame moving with the leading order mean flow. It differs from ω only by
a Doppler-shift term. Equations C0,2 = 0 hold by the same solenoidality (3.9).
3.3. First-order analysis

From C1,0 = 0 we obtain next order constraints for the mean flow
∂T û0,0 + ∂X P̂0,0 = −û∗0,1 ∂X û0,1 − ŵ0,1

∗
∂Z û0,1
∗ ∗
−ikz û0,1 ŵ1,1 − ikz û1,1 ŵ0,1
∗
−iN kx B̂0,1 P̂0,1 + c.c., (3.13)
2
(∂Z − N )P̂0,0 − N B̂1,0 = −û∗0,1 ∂X ŵ0,1 − ∗
ŵ0,1 ∂Z ŵ0,1
∗ ∗
−ikx û1,1 ŵ0,1 − ikx û0,1 ŵ1,1
∗
−iN kz P̂0,1 B̂0,1 + c.c., (3.14)
N ŵ1,0 = −û∗0,1 ∂X B̂0,1 − ∗
ŵ0,1 ∂Z B̂0,1
2 ∗
−(N + ηN )ŵ0,1 B̂0,1 + c.c., (3.15)
∂X û0,0 = 0, (3.16)
where c.c. stands for the complex conjugate of the previous terms. Primary harmonics
are fixed by C1,1 = 0: it follows a linear inhomogeneous system for Û1,1 ,
M(−iω̂, ik)Û1,1 = −R(Û1,2 , ik)Û0,−1 − L(δT , ∇ε )Û0,1 (3.17)

with R being a matrix
 
ikx û1,2 − ikz ŵ1,2 2ikz û1,2 0 0
 2ikx ŵ1,2 ikz ŵ1,2 − ikx û1,2 0 0
R(Û1,2 , ik) =   (3.18)
 2ikx B̂1,2 2ikz B̂1,2 −ikx û1,2 − ikz ŵ1,2 0
0 0 0 0
and L being a linear differential operator
 
δT ∂Z û0,0 0 ∂X
0 δT 0 ∂Z − N 2 
L(δT , ∇ε ) = 
0
 (3.19)
0 δT 0 
∂X ∂Z + N 2 + ηρ 0 0
with δT = ∂T +û0,0 ∂X +ikx û1,0 +ikz ŵ1,0 . Equations C1,2 = C1,3 = 0 yield the solenoidality,
ikx û1,2 + ikz ŵ1,2 = 0, (3.20)
for the velocity field of the secondary harmonics.
3.4. Second-order analysis

From C2,0 = 0 we obtain among others
∂X û1,0 + (∂Z + N 2 + ηρ )ŵ1,0 = 0. (3.21)
Here, we have already used the solenoidality (3.9). At this point we want to stop the
iteration as we found enough constraints to construct leading order solutions.
3.5. Derivation of the modulation equations

The objective of this section is to exploit solvability constraints from the asymptotic
analysis (3.9)–(3.21) in order to derive a closed set of equations which we will denote
as modulation equations. Equation (3.11) has non-trivial solution only if det(M) = 0
yielding the dispersion relation
N 2 kx2
ω̂ 2 = . (3.22)
kx2 + kz2
Note that sgn(ω̂) = ±1 indicates two different branches of the solution. By definition
(3.5), the dispersion relation is essentially a prognostic equation for Φ which is by
definition (3.6) of the form of a Hamilton-Jacobi equation. We find a particular vector
U ∈ ker(M) providing the polarization relation
T
kz ω̂ 2

kz ω̂ ω̂
U(ω̂, k) = −i , i , 1, − i 2 (3.23)
kx N N kx N
which is exactly the same polarization one would obtain in Boussinesq theory. For a
detailed discussion on polarization relations in different sound-proof approximations we
refer to Davies et al. (2003). We deduce that rank(M) = 3. Since hence dim(ker(M)) = 1,
we get a general solution for (3.11) by
Û0,1 = A U(ω̂, k) where Û0,1 = (û0,1 , ŵ0,1 , B̂0,1 , P̂0,1 )T (3.24)
with a yet free complex valued scalar A, that we later refer to as wave amplitude. Note
in particular that B̂0,1 ≡ A. By means of solenoidality of the several harmonics in the
velocity, (3.9) and (3.20), we can show that R(Û1,2 , k)Û0,−1 = 0 which simplifies (3.17).
A necessary condition for the solvability of (3.17) is that any vector LÛ0,1 ∈ img(M) has
to be orthogonal to every U ∈ ker(M) by the Rank-nullity theorem:
0 = U ∗T LÛ0,1 . (3.25)
If equation (3.25) holds, then
∗T
0 = Û0,1 LÛ0,1 + c.c. (3.26)
∗T
0 = Û0,1 LÛ0,1 − c.c. (3.27)
must be true independently, if A 6= 0, which will serve to fix the degrees of freedom in

A = |A| exp i arg(A) . (3.28)
This argument essentially means that the real and the imaginary part of the rhs of (3.25)
must vanish. Equation (3.26) is true if
kx ∂ ω̂
∂T |A|2 + (∇ε + ηρ ez ) · (cg |A|2 ) + |A|2 ∂Z û0,0 = 0 (3.29)
ω̂ ∂kz
which is derived using (3.16). This equation of motion yields essentially the evolution of
the wave energy E e ≡ 2ρ0 |A|2 that will be extensively discussed in §4. With the help of
the consistency relations, ∂T k + ∇ε ω = 0 and ∂X kz = ∂Z kx , that we obtain by cross
differentiating the definitions (3.5) and (3.6), we find
∂ ω̂
(∂T + cg · ∇ε )ω̂ = −kx ∂Z û0,0 (3.30)
∂kz
with the group velocity
∂ ω̂
cg = + û0,0 ex . (3.31)
∂k
Combining ρ0 ×(3.29) with (3.30) and exploiting that ρ0 (∂Z + ηρ ) · = ∂Z (ρ0 · ) result in
a conservation law
∂T A + ∇ε · (cg A) = 0, (3.32)
for the wave action density A = E/ω̂.

e
Let us return to the orthogonality of the kernel vector and the image vector (3.25).
Exploiting (3.27) gives an equation of motion for the argument of A being basically some
slow phase in addition to the fast phase ε−1 Φ, i.e.
(∂T + cg · ∇ε ) arg(A) = −kx û1,0 (3.33)
which contains the first order mean-flow horizontal wind. With the aid of the derivations
from above, we can handle the mean-flow equations (3.13)-(3.16), and (3.21). Using once
more solenoidality from (3.9) and (3.20) it is possible to settle all the terms depending
on Û1,1 , so
∂T û0,0 + ∂X P̂0,0 = −∂X |û0,1 |2 − (∂Z + ηρ )(û∗0,1 ŵ0,1 ) + c.c., (3.34)

ŵ1,0 = 0, (3.35)
∂X û0,0 = 0, (3.36)
∂X û1,0 = 0. (3.37)
If we plug (3.24) into (3.34), we obtain an equation of motion for the leading order
mean-flow horizontal wind

ρ0 ∂T û0,0 + ∂X P̂0,0 = −∇ε · (cg − û0,0 ex )kx A . (3.38)
The six equations (3.22), (3.32), (3.33), (3.36), (3.37), and (3.38), that we derived in this
section, build a coupled system of equations for the six unknowns
Φ, |A|, arg(A), û0,0 , P̂0,0 , and û1,0 .
These are the aforementioned modulation equations. Once they are solved, the leading
order, weak asymptotic solution of the scaled Euler equations has been found. It is
worth mentioning that further solvability constraints for the first order mean flow will
show up when we keep iterating the asymptotic scheme. The equations are in line with
the calculations of Achatz et al. (2010). However in addition, we provide explicitly a
prognostic equation for the slow phase arg(A) and a constraint for the first order mean-
flow horizontal wind since they contribute to the leading order solution. Equation (3.14)
containing B̂1,0 can be considered to be uncoupled as it does not affect neither the leading
order solution nor does it appear in one of the other modulation equations.
4. Energy balance
An important aspect of wave dynamics is the energy budget. As our theory allows for a
wave-mean-flow interaction, we want to investigate the energy exchange. We can rewrite
e = 2ρ0 |A|2 , so
(3.29) in terms of the wave energy density, which we may define as E

e + ∇ ε · cg E
∂T E e = −cgz kx A ∂Z û0,0 . (4.1)
On the lhs we find a conservative flux, that transports energy with the group velocity, and
on the rhs an energy source term that is associated with the vertical shear of the mean-
flow horizontal wind. Let us now consider the mean-flow equation (3.38). Multiplying it
by û0,0 , using the chain rule and (3.16), give us an evolution equation for the mean-flow
kinetic energy density E = ρ0 û20,0 /2, so

∂T E + ∇ε · ρ0 û0,0 P̂0,0 ex + cgz kx Aû0,0 ez = cgz kx A ∂Z û0,0 . (4.2)
On the lhs we find also a conservative flux and on the rhs the same energy source term
as in (4.1) but with the opposite sign making it an energy exchange rate. Because, if we
add (4.1) and (4.2), we obtain a conservation law for the total energy density E = Ee + E.
We will refer to the energy exchange rate, i.e. the rhs of (4.2), as q. In conclusion, the
total energy is conserved by the modulation equations. The sign of the exchange rate
determines the direction of energy flow. E.g., an upward propagating wave packet with
downward-facing phase speed vector transfers energy to the mean flow if the vertical
shear is positive and vice versa.
5. Traveling wave solutions of the modulation equations

A key aspect for derivations of stability properties are traveling wave solutions. These
are stationary, i.e. time-independent, solutions in a translational coordinate system
moving with relative velocity C. Performing a coordinate transformation of the governing
equations and linearizing at the traveling wave solution results in an eigenvalue problem
(due to the stationarity in the translational coordinates) for the time evolution of a small
perturbation. The eigenvalues tell us whether the perturbation is bounded or whether
it grows in time which implies instability. If we find traveling waves for the modulation
equations, we automatically obtain weak asymptotic traveling wave solutions for the
scaled Euler equations. Consider the translational coordinate system Γ = (ξ, ζ, τ ) defined
by the transformation (X, Z, T ) 7→ Γ (X, Z, T ):

ξ X
= − CT and τ = T, (5.1)
ζ Z
with C = (Cx , Cz )T being the constant relative velocity. Note that the unit vectors are
preserved and the mapping is one-to-one. The nabla operator and the time derivative
become
∇ε = ex ∂ξ + ez ∂ζ (5.2)
and
∂T = ∂τ − C · ∇ε , (5.3)
respectively. The modulation equations in the new coordinates may be written as
−(∂τ − C · ∇ε )Φ = ω̂(k) + û0,0 kx , (5.4)
(∂τ − C · ∇ε )a + (∇ε + ηρ ez ) · (cg a) = 0, (5.5)

∂τ + (cg − C) · ∇ε arg(A) = − kx û1,0 , (5.6)

(∂τ − Cz ∂ζ )û0,0 + ∂ξ P̂0,0 = − (∇ε + ηρ ez ) · (cg − û0,0 ex )kx a , (5.7)
∂ξ û0,0 = 0, (5.8)
∂ξ û1,0 = 0. (5.9)
For the sake of the following arguments, we replaced A by the “mass-specific” wave
action density a = ρ−10 A. These equations depend on the background variables N and
ηρ which are functions of Z = ζ + Cz τ due to the transformation, so they depend now in
particular on time. In other words, the background in the translational reference frame is
not hydrostatic anymore since the frame moves relatively to the background. We conclude
that stationary solutions can only exist in two disjoint cases:
C1 The reference frame does not move vertically, such that Cz = 0. In this case we find
that Z = ζ and hence the background is again time-independent. We will refer to this
case as horizontally traveling wave solution. It will be examined in §5.1.
C2 The background variables are constant in the first place and in general Cz 6= 0. If we
combine the definition of the background temperature (2.15) with (2.17), we obtain
κ−1 1
= ηρ + N 2 . (5.10)
κ T0
From this brief calculation we can deduce that the background temperature must be
constant if the rhs is constant. By this point we will refer to this case as isothermal
traveling wave solution which will be investigated in §5.2.
In both cases, a stationary solution in the co-moving reference frame to (5.4) is given by
Z ζ
Φ(ξ, ζ) = Kx ξ + kz (ζ 0 )dζ 0 (5.11)
ζ0
Where Kx is a constant being, without loss of generality, positive and kz the root of the
4th order polynomial derived from the dispersion relation
2 N 2 Kx2
Cz kz + (Cx − û00 )Kx = . (5.12)
Kx2 + kz2
a b c
Cx
√
N2 xobs x u ρ0
Figure 1. Illustration of the horizontally traveling wave solution without mean flow experiencing
non-uniform stratification. Panel a: Squared Brunt-Väisälä frequency profile. Panel b: Horizontal
wind component multiplied by the square root of the background density in the stationary
coordinate system (x, z). Also, the translational coordinate system (ξ, ζ) that moves with the
constant velocity Cx to the right. Panel c: Profile that one would observe at xobs .
The point ζ0 denotes a reference point for the boundary conditions and may be chosen
as ζ0 = 0, i.e. the bottom.
5.1. The horizontally traveling wave solution

In this section we want to investigate the horizontally traveling wave solution (C1 ),
because it is possible to find analytic, explicit solutions in this case, even for arbitrary
stratification. These solutions possess already some interesting properties such as non-
linear phase and amplitude, as well as arbitrary mean flow horizontal wind.
If Cz = 0, then û0,0 disappears in (5.7), and hence (5.8) is solved by any real function
û0,0 (ζ) = uf 0 (ζ). (5.13)
Although we did not show the equation, an evolution equation for û1,0 emerges at higher
order in the asymptotic iteration. But fortunately, the same argument holds due to (5.9),
so it turns out that û1,0 (ζ) = uf 1 (ζ) is also an arbitrary function. Equation (5.12) reduces
to a quadratic polynomial that is solved by
s
N 2 (ζ)
kz (ζ) = ± 2 − Kx2 . (5.14)
Cx − uf 0 (ζ)
It is worth pointing out that a negative discriminant results in an exponential growth or
decay of the wave in the vertical direction corresponding to a reflective layer. To prevent
this, we constrain N/Kx > |Cz − uf 0 | at any point. Also, we must not have uf 0 = Cx at
any point as it possibly causes singularities which corresponds to a critical layer.
To solve (5.5), we drop the time derivative and rewrite it as
cgx − Cx
∂ζ G + ∂ξ G = 0 (5.15)
cgz
which is an equation for the vertical wave action flux G = ρ0 cgz a. Its validity requires
that cgz 6= 0 which is fulfilled if kz 6= 0 and bounded. This manipulation is feasible
as cg depends only on ζ due to (5.13) and (5.14). p As soon as we find a solution for
G, the amplitude can then be computed by |A| = Gω̂/(2ρ0 cgz ). By the method of
characteristics a solution is given by
G(ξ, ζ) = G0 (Ψ ) (5.16)
with
ζ
cgx (ζ 0 ) − Cx 0
Z
Ψ (ξ, ζ) = Kx ξ − Kx dζ . (5.17)
ζ0 cgz (ζ 0 )
The function G0 can be determined by a boundary condition at ζ0 . Equation (5.17) can
be further simplified if we plug in the definition (3.31) and use the dispersion relation
(5.4),
Z ζ
Kx2
Ψ (ξ, ζ) = Kx ξ − 0)
dζ 0 . (5.18)
k
ζ0 z (ζ
Note that the characteristic curves given through Ψ are locally orthogonal on every
line with constant phase Φ. This remarkable property will be investigated in detail in
appendix B and used to construct the substitution that was postulated in order to obtain
weak asymptotic solutions in §3 and appendix A. Taking advantage of (5.9) we can
compute a solution to (5.6) by separation of variables, so
Z ζ
uf 1 (ζ 0 ) 0
arg(A)(ξ, ζ) = LΨ Ψ (ξ, ζ) − Kx 0
dζ (5.19)
ζ0 cgz (ζ )
with LΨ some constant which may be interpreted as a slow wavenumber. A general

solution to (5.7) is found when inserting (5.5)

P0,0 (ξ, ζ) = Kx uf 0 (ζ) − Cx a(ξ, ζ) + Pf (ζ) (5.20)
with Pf some function.
Before we continue to derive an explicit, analytic solution, let us consider the energy
budget. With the aid of (5.16) we can compute the energy exchange rate as presented in
§4 as
q = Kx G∂Z uf 0 . (5.21)
If G is locally confined the energy exchange between wave and mean flow happens there-
fore on the characteristic curves defined by Ψ . And is strongest where the characteristic
curve of the peak of G meets the level of strongest vertical shear.
We were able to derive a complete general solution in the case of a horizontally
traveling wave. However, the solution is still not explicit since integrals emerged for Φ, Ψ
and arg(A). They may be solved numerically. Though, for some particular background
stratification profiles and horizontal mean-flow winds these integrals may have analytic
solutions.
Consider a background with constant mean-flow, such that
uf 0 = Uf = const. and uf 1 = 0, (5.22)
and a stratification profile
2 !−2
2 2 2 2 ζ − ζ0
N (ζ) = Nmin + (Nmax − Nmin ) 1+ (5.23)
∆ζ
2
which represents a bell function with a peak of Nmax centered at ζ0 and width ∆ζ. For
2
ζ → ±∞ it converges to Nmin . When we insert these assumptions into (5.14), we obtain
p
2 2
2 !−1
Nmax − Nmin ζ − ζ0
kz (ζ) = ±Kx 1+ . (5.24)
Nmin ∆ζ
Here, we chose Cx = Nmin /Kx + Uf . It can be integrated by means of standard

antiderivatives, so
p
2 2

Nmax − Nmin ζ − ζ0
Φ(ξ, ζ) = Kx ξ ± Kx ∆ζ arctan (5.25)
Nmin ∆ζ
and
3 !
|Nmin | ζ − ζ0 1 ζ − ζ0
Ψ (ξ, ζ) = Kx ξ ∓ Kx ∆ζ p + . (5.26)
2
Nmax 2
− Nmin ∆ζ 3 ∆ζ
The slow phase becomes simply arg(A) = LΨ Ψ .

In conclusion, for this particular stratification we found an analytic traveling wave
solution. It is plotted in figure 1 for some choice of O(1)-parameters and G0 a bell
function. This traveling wave is characterized by spatial modulation of phase as well as
amplitude which is completely controlled by the stratification (see panel a). The structure
is that of a locally confined wave packet. For ∆ζ → ∞ the background tends to uniformity
and the traveling wave merges to a plane wave. As the wave travels to the right, an
observer at position xobs (see panel b) would experience an upward propagating wave
packet with downward traveling wave troughs and crests. As we decided for a constant
mean-flow horizontal wind, there is no shear. So, there is no energy exchange between
the wave and the mean-flow.
If Cx = 0, the result resembles a stationary wave in the original (X, Z)-coordinates.
Such a wave can be interpreted as a mountain lee wave. The orography then may
fix the boundary condition for G0 . It is also possible to construct analytic traveling
wave solutions in the same fashion, where N = const., corresponding to an isothermal
background, and uf 0 (ζ) a prescribed profile. E.g., choosing uf 0 (ζ) = u0 + Λζ linear
provides also standard antiderivatives. Here, the energy exchange rate q is controlled by
the constant shear Λ. This could be used to study the energy transport due to orography
in circulation models.
5.2. The isothermal, horizontally periodic traveling wave

The previous section was dedicated to the horizontally traveling wave solutions C1
where we found exact analytic solutions to the modulation equations. This section
examines case C2, that is the isothermal traveling wave solutions, where the local Brunt-
Väisälä frequency, N , and the logarithmic background derivative of the density, ηρ , are
constant. In this case a complete leading order solution of the scaled Euler equations
cannot be achieved as we did not provide the explicit evolution equation for û1,0 , which
is mandatory to compute arg(A) in (5.6). This is because the argument in (5.13), where
û1,0 disappeared in the evolution equation, does not hold since in general we allow for a
non-zero vertical relative velocity, Cz 6= 0. However, the remaining modulation equations
are uncoupled, so they may possess some insightful solution.
Even though the 4th order polynomial in (5.12) has an analytic solution, its calculation
does not provide any further enlightenment. We differentiate (5.4) instead w.r.t. ζ, which
avoids lengthy terms and gives us an evolution equation for the vertical wave number.
Closure of the system as well as substantial simplification can be gained if we restrict
ourselves to horizontally periodic traveling waves, such that the horizontal derivatives in
the modulation equations drop out. So, we may rewrite the modulation equations like
(∂τ − Cz ∂ζ )kz + ∂ζ (ω̂ + Kx û0,0 ) = 0, (5.27)
0
(∂τ − Cz ∂ζ )a + (∂ζ + ηρ )(ω̂ a) = 0, (5.28)
0
(∂τ − Cz ∂ζ )û0,0 = −(∂ζ + ηρ )(ω̂ Kx a). (5.29)
The prime indicates the derivative w.r.t. kz , such that ω̂ 0 (kz ) = cgz (kz ).
This system has been studied numerically by Rieper et al. (2013a). A similiar system
has been investigated by Dosser & Sutherland (2011). If we combine (5.28) with (5.29)
and integrate, we obtain
û0,0 = Kx a + u0 (5.30)
with u0 = u0 (ζ + Cz τ ) being fixed by some boundary conditions. Injecting into (5.27),
so
∂τ kz + ∂ζ (ω̂ − Cz kz + Kx2 a + Kx u0 ) = 0, (5.31)
∂τ a + ∂ζ (ω̂ 0 − Cz )a = −ηρ ω̂ 0 a,

(5.32)
reduces the number of equations. As we seek for stationary solutions in the translational
reference frame, we drop the time derivatives and constrain u0 = const., making it some
mean-flow horizontal wind offset, to get rid of its τ -dependency. By that we can integrate
(5.31) and solve for
ω̂ − Cz kz
a= + a0 . (5.33)
−Kx2
The constant a0 absorbed u0 . Hereinafter, we handle a = a(kz ) as function of kz ,
so it becomes a diagnostic quantity. Inserting it into (5.32) yields a single, nonlinear,
autonomous ODE for kz , that can be written after some basic manipulation as
2ηρ ω̂ 0 (kz )a(kz )
∂ζ kz = f (kz ) with f (kz ) = . (5.34)
Kx2 a2 00 (kz )
The denominator of f occurred after differentiating (5.33) and using the identity
00
a2 (kz )
= a02 (kz ) + a(kz )a00 (kz ). (5.35)
2
The ODE (5.34) has equilibrium solutions if the numerator of f vanishes, which happens
either for ω̂ 0 (kz− ) = 0 or for a(kz+ ) = 0. By definition (3.31) it turns out that kz− = 0.
Instead of solving (5.33) for kz+ , which means again finding the root of a quartic function,
we can use it to fix its constant
ω̂(kz+ ) − Cz kz+
a0 = . (5.36)
Kx2
Both equilibria seem uninteresting at first glance: for kz+ the amplitude vanishes and for
kz− there is no oscillation in the vertical, such that wave vector is purely horizontal. But
maybe it exists some solution connecting asymptotically both equilibria which can then
be considered as asymptotic rest states, such that
lim kz = kz± . (5.37)
ζ→±∞
Such a solution may be called wave front. Consider, without loss of generality, kz+ < 0
and an initial condition kz (0) = kz0 ∈ I = (kz+ , kz− ). We want to settle two questions: first,
for which initial condition and set of parameters does a solution exist? And second, to
which of the two asymptotic rest states does the solution converge? That is to determine
the sign of f in (5.34). As we seek for upward propagating wave packets, we shall set
sgn(ω̂) = +1, such that ω̂ 0 > 0 on I. To guarantee that |A| ∈ R+ , we must constrain
a > 0 on I. In particular, we require a(kz− ) = a− > 0 which we can use to settle
ω̂(kz+ ) − N − Kx2 a−
Cz = . (5.38)
kz+
It emerges that Cz > 0, because we can estimate ω̂(kz+ ) < N for all kz+ < 0. Our
wave therefore always travels upwards. Note that by (2.18) and ηT = 0 (isothermal
background), we get ηρ = −κ−1 N 2 . So, the numerator of f is negative on I. In order to
obtain the sign of the denominator, we must consider (5.35) where
ω̂ 0 (kz ) − Cz ω̂ 00 (kz )
a0 (kz ) = and a00 (kz ) = . (5.39)
−Kx2 −Kx2
If we evaluate (5.35) at kz− = 0, we find that
00
a2 (kz− ) C2 a− N
= z4 + > 0. (5.40)
2 Kx Kx4
As the denominator must not assume a root, a cannot have a saddle point on I. Neither,
00
is a allowed to exhibit a local maximum, because a2 would switch sign and would
particularly assume a root as it is continuous. The continuity of a prohibits also a
local minimum, since it would imply a local maximum which we already ruled out. In
conclusion, a must not have any critical points which yields the restriction a0 > 0 on I,
i.e. a is strictly increasing in kz . This is fulfilled regarding (5.39) if ω̂ 0 < Cz on I. It turns
out that

Kx
arg max ω̂ 0 (kz ) = max kz+ , − √ . (5.41)
kz ∈I 2
This gives us a critical relative velocity for the existence of bounded traveling wave
solutions
( √ √
crit + ω̂ 0 (−Kx / 2), if kz+ < −Kx / 2
Cz (kz ) = √ (5.42)
ω̂ 0 (kz+ ), if − Kx / 2 6 kz+ < 0
where
√ 2 N
ω̂ 0 (−Kx / 2) = √ . (5.43)
3 3 Kx
Demanding Cz > Czcrit (kz+ ) generates by (5.38) a lower bound on
ω̂(kz+ ) − kz+ Czcrit (kz+ ) − N
a− > acrit +
− (kz ) = . (5.44)
Kx2
a− 12
bounded
z/Lref
acrit +
− (kz )
Hθ
C
unbounded
√ 0
−Kx / 2 0 kz+ 0 x/Lref 12
Figure 2. Critical asymptotic mass-spe- Figure 3. Illustration of the isothermal
cific wave action density. The isothermal traveling wave front in two dimensions
traveling wave solutions cannot exist if the propagating with velocity C. Horizontal
asymptotic rest states kz+ and a− belong to wind u∗ as given by (5.45). And the
the gray zone. potential temperature scale height Hθ for
comparison.
a b c
20
Cz Cz Cz
z/Lref
Hθ
0
kz+ 0 kz 0 a− a u0 0 u∗
Figure 4. Illustration of the isothermal traveling wave front in one dimension. (panel a)
Vertical wave number. (panel b) Mass-specific wave action density. (panel c) Horizontal wind as
given by (5.45) and the potential temperature scale height Hθ for comparison.
The critical asymptotic mass-specific wave action density acrit

− is illustrated in figure 2.
00
Given that the previous constraints hold: by continuity a2 is positive on I and the sign
of f is therefore negative. This means, for any kz0 ∈ I the solution converges to the
asymptotic rest state kz+ as ζ → +∞. Vice versa, for ζ → −∞ the solution tends to
kz− = 0. So, as long as for some kz+ we choose a− , such that it is in the white zone in
figure 2, the solution converges. As soon as a− < acrit
− (grey zone), the solution will run
into a singularity and grow without bounds.
We want to point out that Cz is also the vertical component of the phase velocity.
Multiplying by kz gives us therefore the frequency in the stationary coordinate system.
We may interpret (5.42) in terms of a critical frequency for a front-like traveling wave.
In particular, for a lower boundary condition on kz , that is an initial condition in the
stationary system, this leads to a minimal frequency for the existence of traveling wave
solutions.
Lref /km Hθ /km ε N Kx a− kz+ Cz

3.053 30.53 0.1 1 1 1 -1.3 1.1
Table 2. O(1)-parameters defining an isothermal, horizontally periodic traveling wave
The ODE (5.34) can be solved with the explicit Euler method. For some O(1)-
parameters (table 2), that fulfill the requirements, a particular solution illustrating
the traveling wave front’s properties is plotted in figure 4. Panel a shows the directly
computed solution of the vertical wavenumber kz , panel b has the mass-specific wave
action density a diagnosed from (5.33), and panel c the leading order horizontal wind
−1
u∗ = û0,0 + (û0,1 eiΦε + c.c.), (5.45)
which is also shown in figure 3 from a two-dimensional perspective. We chose u0 = −Kx a−

for (5.30), so the mean-flow horizontal wind drops to zero after the front propagated
through. To compute Φ from (5.11), we applied the very same Euler method. To compute
û0,1 according to (3.24), we set arg(A) = const. because we lack an equation for it. The
error produced by this oversimplification occurs in the wake of the front when kz tends to
zero, such that the slow phase arg(A) becomes vertically dominant and maybe produces
some long wavelength vertical oscillation. From the two-dimensional point of view, the
leading term in Φ is then still Kx ξ. However, this will not change the overall picture
as the induced mean-flow regarding (5.30) is not affected by the slow phase. Again this
traveling wave is characterized by spatial modulation of amplitude as well as phase.
The mass-specific wave action density and therefore the amplitude remain constant
in the wake of the front. So at first glance, it seems like the wave induces wave energy
into the mean flow. Because one would usually expect increasing amplitude with height
due to the thinning air when the wave propagates upwards and due to the fact that
the wave energy scales with the background density. This “extra” energy should then
accelerate the mean-flow. This interpretation turns out to be incorrect, which becomes
clear if we compute the wave-mean-flow energy exchange rate from §4 in the stationary
(X, Z)-coordinate system,
q = Kx2 ρ0 a(kz )ω̂ 0 (kz )a0 (kz )∂Z kz < 0. (5.46)
Here, we used (5.30) and the estimates from this section. We can already conclude that
energy is exchanged directly in the front where the gradient is located. The rate is always
negative. On the contrary, the wave gains energy from the mean flow and the mean-flow
horizontal wind decelerates. This surprising behavior dissolves if we consider also the
wave energy flux from (4.1). It tells us that the wave energy is transported with the
group velocity ω̂ 0 . But the front moves with Cz which is always greater than ω̂ 0 as we
suppressed any critical points for a.
In summary, the fast traveling wave front only sustains if it obtains energy from the
mean flow. One may ask if inverted traveling wave fronts possibly exist as solutions to the
isothermal, horizontally periodic modulation equations. These would be solutions where
the asymptotic rest states switched. So, kz− 6= 0 and a− = 0 but kz+ = 0 and a+ 6= 0.
These solutions would be referred to as traveling wave backs. It turns out that backs
cannot be solutions under no circumstances. We state this claim here unproved as the
argumentation goes beyond the scope of this survey.
Lref /km lx /Lref lz /Lref Nmin Nmax ∆ζ Kx LΨ ∆Ψ Gmax sgn(kz )

3.053 20 20 0.5 1.5 0.7 1.5 -0.3 1.0 3×10−5 -1
Table 3. Fixed parameters that define the traveling wave and experimental setup
6. Numerical validation of the traveling wave solutions under the full

Euler dynamics
In the previous sections we found that traveling wave solutions for the modulation
equations appear in two distinct cases. For the first case (C1 ), when the atmosphere
is non-uniformly stratified, we were able to derive exact analytic solutions. The wave
vectors can point in any direction but vertical. Though, the entire wave travels only in
the horizontal direction. For the second case (C2 ), when the atmosphere is uniformly
stratified (isothermal), we managed to show existence and the structure of the traveling
wave solutions: assuming horizontal periodicity it occurs to be a front-like solution. These
traveling waves propagate both horizontally and vertically. We found exact solutions
to the modulation equations. But they are, however, asymptotic solutions to the fully
nonlinear Euler equations. In this section we want to test the consistency of these
solutions. The traveling wave solutions of the C1 -case will be challenged numerically,
so we can check the range of validity as the traveling wave solutions could be unstable
intrinsically or only sustain for very small, even unrealistic ε. For this purpose, we will
use a conservative solver of the Euler equations, the pincFloit model.
PincFloit stands for Pseudo-Incompressible Flow Solver with Implicit Turbulence and
was introduced by Rieper et al. (2013b). This finite-volume flow solver conserves energy,
momentum as well as mass and suits therefore the purpose well. Within the framework of
the model it is possible to simulate with almost unconstrained background stratifications
apart from the uniform, weakly stratified atmosphere. The fidelity of the model was
tested against standard cases. Rieper et al. (2013a) applied pincFloit to check the range
of validity of the extended WKB theory for gravity waves of Achatz et al. (2010) which
provides the basis for our traveling wave solutions. Recently, Muraschko et al. (2015);
Bölöni et al. (2016) derived a numerical solver for the modulation equations based
on a Lagrangian ray-tracing approach for horizontally homogeneous GW packets that
was also tested with the aid of pincFloit. It discretizes the governing equations on an
equidistant staggered C-grid. From the alternative schemes, that the model provides, it
proved practical to apply the third order, low storage Runge-Kutta time stepping scheme,
which is total variation diminishing. This is combined with the MUSCL scheme for spatial
discretization and a monotonized central flux limiter. As the name of the model implies,
it is capable of solving the compressible Navier-Stokes equations. But since this survey
is only interested in the Euler equations, any diffusion or dissipation are switched off.
We validate the C1 traveling wave solution with non-uniform stratification as given by
the equations (5.22)-(5.26). They are not compactly supported on the two-dimensional
domain. Inevitably, the wave will touch the boundary of the computational domain.
In order to prevent reflection, we prescribe the analytic solution on the boundary cells
for all times. We choose a particularly challenging test case such that the background
stratification and the amplitude varies like O(1) on the whole domain. The mean-flow
horizontal wind is chosen in a way that the horizontal translational, relative velocity Cx
vanishes, so
Uf = −Nmin /Kx . (6.1)
By this choice the traveling wave is stationary in the fix coordinate system, such that
(ξ, ζ) = (X, Z), which makes it convenient to compare the initial condition with the final
state. The wave amplitude is set to be of the same order as the mean flow to provoke
nonlinear coupling. So, the setup can be viewed as a wave that wants inherently to migrate
towards the right side with a constant horizontal group velocity but faces an opposing
wind of the exact same velocity. If the model manages to balance these tendencies the
outcome stands at rest.
The parameters that define the traveling wave and the background are presented in
table 3. They were chosen in a way that they resemble a typically unresolved wave in
weather and climate models. Here, lx and lz are the domain size in x- and z-direction,
respectively. The initial wave action flux is given by a bell shaped function
4 !−2
Ψ (X, Z) − Ψ (Xpos , Zpos )
G(X, Z) = Gmax 1 + (6.2)
∆Ψ
where Gmax is the maximum flux at the peak.

To prescribe the stratification in the model, we have to provide the background density,
Exner pressure, and potential temperature. The latter can be computed by means of
standard antiderivatives from (5.23). The background Exner pressure is defined through
the hydrostatic equation (2.13) which we must integrate numerically. This is achieved by
a simple Euler forward scheme with very high resolution. It looks costly at first glance
but it only needs to be computed once at start.
Deviations between simulation and theory are measured by the weighted norms of
the numerically simulated field (num) minus the analytically expected (true) field.
Hereinafter, we consider the horizontal wind field as we did not find any significant
differences from the other fields in terms of errors. We define the weighted error norms
as

I |utrue − unum |
l1 =
I |utrue |
!1/2
I |utrue − unum |2
l2 =
I |utrue |2

maxx,z |utrue − unum |
l∞ = (6.3)
maxx,z |utrue |
where I denotes the numerical approximation to the integral
Z lx Z lz
I(u) ≈ u(x, z)dzdx. (6.4)
0 0
The numerical quadrature exploits the trapezoidal rule.

In order to test the consistency of the traveling wave solutions, the scale separation
parameter ε is altered whereas all other parameters, as given in table 3 determining the
properties of the traveling wave, are kept constant. The analytical asymptotic solution
in non-dimensional form can be written as
−1

Utrue = Û0,0 + Û0,1 eiε Φ + c.c. + O(ε) (6.5)
for ε → 0 and (X, Z, T ) fixed. We call the solution is consistent if we can reproduce the
analytically expected order of convergence from (3.1) of M = 1 numerically for a broad
Hθ /km ε t∗end /s
45.57 0.067 1204
37.23 0.082 889
30.53 0.100 660
25.03 0.122 490
20.35 0.150 360
Table 4. Varying parameters for the experiments testing the consistency of the traveling
wave.
range of ε. In §2 we defined
Lref
ε= . (6.6)
Hθ
As we vary ε, we keep Lref fixed, such that the spatial resolution per wave period remains
constant. This choice provides the benefit that we do not need to adapt the resolution
for each value of ε. If we did so, then we could not distinguish effects on the error of the
varying ε from the varying resolution. All experiments were performed with a resolution
of 256 × 512 grid cells, which corresponds to approximately 18 grid cells per wave period
which is a reasonable value to resolve all important features.
As we change ε from one experiment to another, we must constrain that T = εt = O(1),
which fulfills a necessity of the asymptotic ansatz, and hence the dimensionful overall
simulation time must change like

tref p
t∗end = O = O ε−2/3 since tref = Hθ /g = O(ε−1/2 ). (6.7)
ε
Table 4 presents the parameters defining each of the five experiments. In particular,
the simulation time t∗end is shown in the rightmost column. The potential temperature
scale height varies from roughly 20 to 46 km throughout the simulations which gives us
a considerable range for the scale separation.
The simulation results are plotted in figure 5. The stratification, horizontal wind, and
its relative deviation from the analytic solution are depicted for three of the experiments.
For smallest ε = 0.067 (top row) the stratification is moderate and the amplitude
envelope, which is given by the wave action flux G, is also rather broad in width. The
horizontal wind oscillates around a mean flow of about 12 ms−1 with an amplitude of
3 ms−1 .
The lines of constant phase are nearly straight, so are the characteristic lines Ψ , which
define the shape of the envelope. With even smaller ε the wave will get closer to the
uniformly stratified regime with infinitesimal amplitude which can be captured by the
Boussinesq theory. The relative error after the simulation is as expected very small,
much less then 1 % at the peaks. Remarkably, the error is almost perfectly in phase
with the initial wave at least in the middle of the domain. This means that errors are not
generated by non-stationary effects. The wave does not move at all, the peaks are slightly
smeared out, probably, by some numerical dissipation. By increasing ε, the width of the
stratification peak narrows. Also, the envelope confines and its maximum increases. The
mean-flow horizontal wind gets stronger and we can expect nonlinear effects.
For ε = 0.1 (second row) the mean flow approaches 15 ms−1 and the horizontal wind
varies with an amplitude of about ±10 ms−1 which produces some substantial shear.
In this case, we observe clearly some error sources at the upper boundary. These are
Stratification Initial condition Initial minus final condition

20
ε =0.067
z/Lref
0
20
ε =0.100
z/Lref
0
20
ε =0.150
z/Lref
0
0 1 0 20 0 20
N 2 /Nref
2
x/Lref x/Lref
-25 -20 -15 -10 -5 -3 0 3

Figure 5. Computational corroboration of the existence of horizontally traveling wave solutions.
(left column) Stratification represented by the local squared Brunt-Väisälä frequency. (middle)
Initial total horizontal wind utrue in ms−1 . (right) Deviation between the numerical and the
approximate asymptotic solutions after simulation time t∗end (see table 4). i.e. utrue − unum .
(middle column, bottom row) The filled white curve depicts a stream line being inflated by a
factor 64.
numerical errors. As we prescribe the expected analytic solution at the boundaries, even
small deviations of the numerical solution in the domain can cause discontinuities directly
at the boundary. These shocks than propagate into the domain.
For even higher ε, nonlinear effects start to dominate the error production. In the
experiment where ε = 0.15 (third row), the mean-flow is about 18 ms−1 and the horizontal
wind oscillates with an amplitude of about 40 ms−1 which is not captured by the color
scale of figure 5 anymore. Here, we find some quadrupole structure in the relative error
field that indicates phase shifts by some instability process which are triggered most
likely by the boundary discontinuities. The local error exceeds 3 ms−1 at some points.
For this particular case we also plotted a stream line at the bottom in figure 5 which is
l1 l2 l∞
0.1
weighted error
M = 2.93 ± 0.65 M = 3.20 ± 0.61 M = 3.09 ± 0.60
0.01
0.001
0.067 0.100 0.150 0.067 0.100 0.150 0.067 0.100 0.150
ε ε ε
Figure 6. Convergence of the error norms with regard to ε. (solid line) Linear regression on a
log-log scale applying the least square method reveals an order of convergence M and its standard
deviation. (long dashed) first order, (medium dashed) second order, and (short dashed) third
order lines.
amplified for better visibility. This line represents just as well a mountain ridge of about
200 m height that could have excited the wave.
For all five experiments, we computed the error norms according to (6.3). The results
are shown in figure 6 where we plotted the error norms against ε on a log-log scale.
A least square fit gives us the order of convergence M assuming that the error norms
converge like
ln = CεM for n = 1, 2, ∞. (6.8)
We also computed its standard deviation which are depicted aside the fitted lines in
figure 6. For all three norms we obtain convergence like O(ε3 ). Even after subtracting
the rounded error margin of 1 from the order, it is apparently still included in the
expected O(ε)-convergence with a comfortable cushion. For larger ε the model probably
overestimates the norms due to the numerical errors at the boundaries propagating into
the domain and triggering instabilities. We tested the influence of this effect by changing
the integration domain of the norms like
Z lx −δ Z lz −δ
I(u) ≈ u(x, z)dzdx (6.9)
δ δ
and δ some positive fixed number. So, we focus on a smaller box and ignore cells that
are potentially influenced by the boundaries. For some values of δ we recalculated the
error norms for every ε. There was no significant effect on the order of convergence. So
either the hypothesis of boundary discontinuities propagating into the domain causing
instabilities is wrong or they already contaminate the whole domain after the simulation
time as they propagate faster than expected or instabilities are excited by other processes.
This remains an open question.
However, for small ε the norms are presumably also overrated as the amplitudes become
so small that nonlinear coupling between the wave and the mean flow becomes inefficient
and numerical damping takes over. Taking these effects into account, with an order of
convergence of about three, we can safely conclude that the traveling waves are consistent
for a significant parameter range.
7. Hydrostatic gravity waves
In the previous section, we showed that the horizontally traveling wave solutions
(C1 ) resemble stationary mountain lee waves at least qualitatively. In the numerical
experiments though, we used rather small wavelengths of about 3 km both vertical
and horizontal. As the overall motivation of this survey is to lay the foundation for
more realistic instability criteria in gravity wave parametrizations, the natural choice
to test the asymptotic consistency of the traveling wave solutions was to use typically
unresolved small wavelengths. But mountain waves in nature are mostly hydrostatic
waves, i.e. they exhibit much longer horizontal than vertical wavelengths, such that their
dynamics are anisotropic. For the derivation of the modulation equations we assumed a
non-hydrostatic, isotropic scaling. This section is concerned with the question whether
our modulation equations and their traveling wave solutions still apply in the hydrostatic
limit.
In this regard Achatz et al. (2010, 2017) showed that the consistency between the
scale asymptotics of the Euler equations and the pseudo-incompressible equations also
holds for hydrostatic gravity waves. Bölöni et al. (2016) investigated the non-hydrostatic
modulation equations numerically applying a ray-tracer method. The model was val-
idated, among others, by initializing hydrostatic waves with a horizontal to vertical
wavelength ratio of 10. The results were compared with reference simulations computed
with pincFloit (Rieper et al. 2013b). Both models were in excellent agreement.
Based on these considerations we will rescale the Euler equations (2.19) for hydrostatic
gravity waves and solve them asymptotically with the spectral WKB method presented
in §3. We assume a horizontal to vertical wavelength ratio to be O(ε) whereas the same
distinguished limit provided in §2 shall hold. In particular, the vertical wavelength is
small in comparison to the scale height and the background varies slowly with respect to
the fast oscillations. In order to capture the anisotropic dynamics, we introduce rescaled
coordinates
(X, Z, T ) = (ε2 x, εz, ε2 t) (7.1)
and consider the following ansatz for the WKB approximation

(0)
 
  uε (x, z, t; ε)
v (1)
εwε (x, z, t; ε)
 
θ  =  (1) .

(7.2)
π  θ 0 (Z) + εθε (x, z, t; ε) 
2 (2)
π0 (Z) + ε πε (x, z, t; ε)
Note that, in contrast to (7.2), the vertical wind occurs at O(ε) (cf. Achatz et al. 2010)
whereby all the other symbols have their usual meaning. If we insert (7.2) into the Euler
equations (2.1) and substitute by
!T
(1)
T 1 θε
U = (u, w, B, P ) := u(0)
ε , εwε(1) , , θ0 πε(2) , (7.3)
N θ0
then we can exploit the series expansion (3.1) combined with the WKB assumption
(X, Z, T ) 7→ Ûn,m (X, Z, T ), Φ(X, Z, T ). (7.4)
By the weak asymptotic approach of appendix A we can order terms in powers of ε and
thereby obtain to leading order a homogenous linear system for the primary harmonics
M Û0,1 = 0 with the coefficient matrix
 
−iω̂h 0 0 ikx
 0 0 −N ikz 
M(−iω̂h , ik) = 
 . (7.5)
0 N −iω̂h 0 
ikx ikz 0 0
The intrinsic frequency for the hydrostatic waves is found to be by ω̂h = ω − kx û0,0
where ω and k are the time and spatial derivatives of the phase Φ, respectively.
The
homogenous system for Û0,1 has non-trivial solution only if det M(−iω̂h , ik) = 0, which
yields the dispersion relation
N 2 kx2
ω̂h2 = . (7.6)
kz2
The kernel of M provides the polarization relation U ∈ ker (M) which turns out to be the
same as the non-hydrostatic (3.23), so that the solution can be written as
Û0,1 = A U(ω̂h , k) (7.7)
with A ∈ C being the scalar amplitude. To leading order, we also get constraints for the
mean flow
ŵ0,0 = 0, B̂0,0 = 0, and ∂X û0,0 = 0. (7.8)
At the next order, we obtain an inhomogeneous system for the primary harmonics
being of the same form as (3.17) with only a change in the linear differential operator
 
δT ∂Z û0,0 0 ∂X
0 0 0 ∂Z − N 2 
L(δT , ∇ε ) = 
0
. (7.9)
0 δT 0 
∂X ∂Z + N 2 + ηρ 0 0
Multiplying the inhomogeneous equation from the left by the conjugate transpose of the
kernel vector and adding its complex conjugate, we obtain a conservation law,

∂T A + ∇ε · cgh A = 0 , (7.10)
for the wave action density A = ρ0 |A|2 /ω̂h which constitutes an equation of motion
for the absolute wave amplitude |A|. By subtracting its complex conjugate, we get an
equation of motion for the argument of the amplitude
(∂T + cgh · ∇ε ) arg(A) = −kx û1,0 . (7.11)
The hydrostatic group velocity is given by
∂ ω̂h
cgh = + û0,0 ex . (7.12)
∂k
Also at the next order, we find the very same equation of motion for the non-hydrostatic
mean-flow horizontal wind (3.34) which can be recast with the aid of the hydrostatic
polarization (7.7) to

ρ0 ∂T û0,0 + ∂X P̂0,0 = −∇ε · (cgh − û0,0 ex )kx A . (7.13)
In conclusion, we find hydrostatic modulation equations (7.10), (7.11), and (7.13) that
differ from the non-hydrostatic ones (3.32), (3.33), and (3.38) only due to the dispersion
relation (7.6). If we perform a Taylor expansion of the non-hydrostatic dispersion relation
(3.22) for |kx | |kz | in the variable kx /kz , i.e. assuming steep wave vectors (Bühler 2009,
p 121), we obtain
ω̂ = ω̂h + O(|kx /kz |2 ). (7.14)
Therefore, we deduce that the non-hydrostatic modulation equations also hold in the
hydrostatic limit. This statement also implies that the traveling wave solutions, which
we found for the non-hydrostatic case, are equally valid for large horizontal wavelength
or small Kx , accordingly. The result is of particular interest for the C1 -class since it
demonstrates that they indeed describe typical mountain lee waves.
8. Conclusion
We started our derivations from the dimensionless, fully compressible Euler equations.
By assuming low Mach number, leading order stable stratification, non-hydrostatic
motion, and a large internal Froude number, we found a distinguished limit similar to
Achatz et al. (2010). The scale separation parameter ε, which is now the only parameter
appearing in the equations, is the ratio of the dominant wavelength and the scale height.
The scaled Euler equations were derived from a multiple scale ansatz that divided the
wave from a hydrostatic background. We combined a spectral expansion approach with
the WKB assumption of slow modulation of amplitude and phase. So, the frequency and
wave numbers are defined locally as the temporal and spatial derivatives of the phase,
respectively. Since the system is nonlinear and we allowed for wave-mean-flow interaction,
even the leading order equations are contaminated by higher harmonics.
Instead of an averaging technique (e.g. Grimshaw 1974; Tabaei & Akylas 2007) that
assumes a priori periodicity of the asymptotic solution with respect to the fast phase, we
applied the more general approach presented in Danilov et al. (2003) of weak asymptotic
solutions.
Investigating the results of the asymptotic analysis leads to the modulation equations:
equations of motion for the phase, wave action density, and the leading order mean-
flow variables. These equations are in line with the findings of Achatz et al. (2010). But
in addition to that we found an evolution equation for the argument of the complex
valued amplitude and a constraint for the higher-order mean-flow horizontal wind which
becomes necessary to close the equation system. Computing the energy budget revealed
that these equations conserve the sum of wave and mean-flow energy.
We found two distinct classes of traveling wave solutions for the modulation equations:
the horizontally traveling wave solution and the isothermal traveling wave solution. For
the former, we were able to construct analytic solutions. The horizontally traveling wave
solution with Cx = 0 may be applied to model mountain lee waves for arbitrary strati-
fication, background horizontal wind, and orography. Typically, mountain lee waves are
hydrostatic. Although the modulation equations were derived for non-hydrostatic gravity
waves, it was shown that they also hold in the hydrostatic limit where the wave vector
becomes steep. With the additional assumption that the wave is horizontally periodic,
we investigated the existence and the shape of the isothermal traveling wave solution. It
resembles an upward propagating wave front. Quite surprisingly, the isothermal traveling
wave front decelerates the mean flow in order to sustain its constant amplitude.
The horizontally traveling waves were validated numerically against the fully nonlinear
Euler equations. Exact solutions of the modulation equations fulfill the Euler dynamics
asymptotically. The purpose of the numerical simulations was to clarify the consistency of
the solutions, i.e. range of validity of the asymptotic approximation. The traveling wave
solutions performed surprisingly well: they exceeded the analytically expected order of
convergence of one, as the scale separation parameter tends to zero. The numerically
computed order turned out to be around three. This result holds on a broad range of
scale separation parameters which proves the consistency of the traveling wave solutions
presented in this survey.
In a companion paper we will analyze the traveling waves, that we found here, with
respect to stability by analytical and numerical methods.
The authors thank the German Research Foundation (DFG) for partial support
through the research unit Multiscale Dynamics of Gravity Waves (MS-GWaves) and
through grants AC 71/8-1, AC 71/9-1, AC 71/10-1, and KL 611/25-1. M.S. thanks
Gergely Bölöni for vital support with pincFloit.
Appendix A
This appendix introduces the weak asymptotic scheme which allows to successively
iterate the asymptotic solution in order to achieve an equation hierarchy. Our objective
is to establish that the coefficients Ck,l from (3.7) vanish independently of k and l.
This is usually achieved by an orthogonality argument that the harmonics are mutually
orthogonal with respect to the inner product given by the average over the interval
[0, 2π] w.r.t. the fast phase ε−1 Φ (Miura & Kruskal 1974; Grimshaw 1974; Tabaei &
Akylas 2007). There are two issues with this approach: First, the fast phase does depend
on ε, such that the phase lines squeeze together when ε tends to zero. Second and most
importantly, Φ is actually a function of space and time, therefore solving the integral
requires an integration by substitution. We can circumnavigate these issues if we seek
for asymptotic solutions in a weak sense (Danilov et al. 2003). Consider a differentiable,
compactly supported test function S. We want to multiply S together with some test
harmonic exp(−imε−1 Φ) to (3.7) and integrate over the entire real line with regard
to the phase Φ. The aforementioned substitution must be an injective, continuously
differentiable mapping, say
(X, Z, T ) 7→ F (X, Z, T ) (A 1)
which is independent of ε with F = (Φ, Ψ, τ ). We will present a particular mapping of
this kind in appendix B. The weak formulation of (3.7) reads then
M Z
X X −1
εk S(Φ)(Ck,l ◦ F −1 )(Φ, Ψ, τ )ei(l−m)ε Φ dΦ + o(εM ) = 0. (A 2)
k=0 |l|6k+2 R
Holding τ and Ψ constant, we call (3.1) a weak asymptotic solution if it fulfills (A 2) for
all m ∈ N and all well-behaved S. We may rewrite (A 2) for the forthcoming argument
as
M Z
X X −1
εk S Ck,m ◦ F −1 + (Ck,l ◦ F −1 )ei(l−m)ε Φ dΦ + o(εM ) = 0. (A 3)
k=0 R |l|6k+2
l6=m
One can show via integration by parts that the sum over l in (A 3) vanishes like O(ε)
when ε tends to zero, as
∂ S(Ck,l ◦ F −1 ) i(l−m)ε−1 Φ
Z Z
−1 i(l−m)ε−1 Φ ε
S(Ck,l ◦ F )e dΦ = − e dΦ = O(ε).
R i(l − m) R ∂Φ
(A 4)
Here, we exploited that S has compact support. The existence of the integral requires that
Ck,l is differentiable which translates to the ansatz functions and background variables
that their second derivative must also exist. With (A 4) we can safely apply the limit
ε → 0 for (A 3) and iterate the solution: The sum over l in (A 3), i.e. the second
summation, does not contribute to the leading order due to (A 4), such that C0,m = 0
for all m and S. In the next order these terms appear, but we can insert the leading
order solution and find that they cancel out, such that also C1,m = 0 for all m and S.
We can successively proceed to order o(εM ). To put it in a nutshell, the scaled Euler
equations (2.19) posses weak asymptotic solutions in the limit ε → 0 of the form (3.1)
only if Ck,m = 0 holds for all possible k and m .
Appendix B
This appendix investigates the existence of weak asymptotic solutions as introduced in
§3 and appendix A. We already stated in §5.1 that the characteristic curves (5.18) are lo-
cally orthogonal on every line with constant phase Φ. This motivates to introduce a curvi-
linear, translational, orthogonal coordinate system Θ = (Φ, Ψ, τ ) via (ξ, ζ, τ ) 7→ Θ(ξ, ζ, τ ).
The orthogonality can be proven when we write down the contravariant local basis
vectors,
Kx kz kz Kx
eΦ = ex + ez and eΨ = ex − ez , (B 1)
kkk kkk kkk kkk
that we obtain from computing the Jacobian of the coordinate transformation Θ. It
becomes now obvious that eΦ · eΨ = 0 which proves the orthogonality. Since kz 6= 0
is restricted, its sign cannot change and hence the integrals in (5.11) and (5.18), which
determine Φ and Ψ , respectively, are monotone, continuous, and especially injective. The
mapping F = Θ ◦ Γ (X, Z, T ) fulfills therefore the requirements for the substitution
presented in appendix A. The existence of such a mapping was at this point only
postulated. As (5.11) holds also for the isothermal traveling wave solution (C2 ), we can
construct the curvilinear transformation there as well. So, this brief appendix justifies in
hindsight the validity of the weak asymptotic scheme.
REFERENCES
Achatz, U. 2007 Gravity-wave breaking: Linear and primary nonlinear dynamics. Adv. Sp. Res.
40, 719–733.
Achatz, U., Klein, R. & Senf, F. 2010 Gravity waves, scale asymptotics and the pseudo-
incompressible equations. J. Fluid Mech. 663, 120–147.
Achatz, U., Ribstein, B., Senf, F. & Klein, R. 2017 The interaction between synoptic-scale
balanced flow and a finite-amplitude mesoscale wave field throughout all atmospheric
layers: weak and moderately strong stratification. Q. J. R. Meteorol. Soc. 143 (702),
342–361.
Alexander, M. J. & Dunkerton, T. J. 1999 A spectral parameterization of mean-flow forcing
due to breaking gravity waves. J. Atmos. Sci. 56 (24), 4167–4182.
Baldwin, M.P., Gray, L.J., Dunkerton, T.J., Hamilton, K., Haynes, P.H., Randel,
W.J., Holton, J.R., Alexander, M.J., Hirota, I., Horinouchi, T., Jones, D.B.A.,
Kinnersley, J.S., Marquardt, C., Sato, K. & Takahashi, M. 2001 The Quasi-
Biennial Oscillation. Rev. Geophys. 39 (2), 179–229.
Becker, E. 2012 Dynamical control of the middle atmosphere. Space Sci. Rev. 168 (1-4),
283–314.
Bölöni, G., Ribstein, B., Muraschko, J., Sgoff, C., Wei, J. & Achatz, U. 2016
The interaction between atmospheric gravity waves and large-scale flows: an efficient
description beyond the non-acceleration paradigm. J. Atmos. Sci. 73 (12), 4833–4852.
Bühler, O. 2009 Waves and mean flow , 1st edn. New York: Cambridge University Press.
Chu, V. H. & Mei, C. C. 1970 On slowly-varying Stokes waves. J. Fluid Mech. 41, 873–887.
Danilov, V. G., Omelyanov, G. A. & Shelkovich, V. M. 2003 Weak asymptotics method
and interaction of nonlinear waves. In Asymptot. methods wave quantum Probl., 208th
edn. (ed. M. V. Karasev), pp. 33–163. Amer. Math. Soc., Providence, RI.
Davies, T., Staniforth, A., Wood, N. & Thuburn, J. 2003 Validity of anelastic and other
equation sets as inferred from normal-mode analysis. Q. J. R. Meteorol. Soc. 129, 2761–
2775.
Dosser, H. V. & Sutherland, B. R. 2011 Weakly nonlinear non-Boussinesq internal gravity
wavepackets. Phys. D Nonlinear Phenom. 240 (3), 346–356.
Fritts, D. C. 2003 Gravity wave dynamics and effects in the middle atmosphere. Rev. Geophys.
41, 1–64.
Grimshaw, R. 1974 Internal gravity waves in a slowly varying, dissipative medium. Geophys.
Fluid Dyn. 6, 131–148.
Klein, R. 2011 On the Regime of Validity of Sound-Proof Model Equations for
Atmospheric Flows. In ECMWF Workshop on Non-Hydrostatic Modelling, November
2010 . http://www.ecmwf.int/publications/library/do/references/list/201010.
Klein, R., Achatz, U., Bresch, D., Knio, O. M. & Smolarkiewicz, P. K. 2010 Regime of
validity of soundproof atmospheric flow models. J. Atmos. Sci. 67 (10), 3226–3237.
Lelong, M.-P. & Dunkerton, T. J. 1998 Inertia-gravity wave breaking in three dimensions.
Part I: convectively stable waves. J. Atmos. Sci. 55 (1997), 2473–2488.
Liu, W., Bretherton, F. P., Liu, Z., Smith, L., Lu, H. & Rutland, C. J. 2010 Breaking
of progressive internal gravity waves: convective instability and shear instability. J. Phys.
Oceanogr. 40, 2243–2263.
Lombard, P. N. & Riley, J. J. 1996 Instability and breakdown of internal gravity waves. I.
Linear stability analysis. Phys. Fluids 8 (12), 3271–3287.
McLandress, C. 1998 On the importance of gravity waves in the middle atmosphere and their
parameterization in general circulation models. J. Atmos. Solar-Terrestrial Phys. 60 (14),
1357–1383.
Mied, R. P. 1976 The occurrence of parametric instabilities in finite-amplitude internal gravity
waves. J. Fluid Mech. 78, 763–784.
Miura, R. M. & Kruskal, M. D. 1974 Application of a nonlinear WKB method to the
Korteweg-DeVries equation. SIAM J. Appl. Math. 26 (2), 376–395.
Muraschko, J., Fruman, M. D., Achatz, U., Hickel, S. & Toledo, Y. 2015 On the
application of Wentzel-Kramer-Brillouin theory for the simulation of the weakly nonlinear
dynamics of gravity waves. Q. J. R. Meteorol. Soc. 141 (688), 676–697.
Rieper, F., Achatz, U. & Klein, R. 2013a Range of validity of an extended WKB theory
for atmospheric gravity waves: one-dimensional and two-dimensional case. J. Fluid Mech.
729, 330–363.
Rieper, F., Hickel, S. & Achatz, U. 2013b A conservative integration of the pseudo-
incompressible equations with implicit turbulence parameterization. Mon. Weather Rev.
141 (3), 861–886.
Sutherland, B. R. 2001 Finite-amplitude internal wavepacket dispersion and breaking. J.
Fluid Mech. 429, 343–380.
Sutherland, B. R. 2006 Weakly nonlinear internal gravity wavepackets. J. Fluid Mech. 569,
249–258.
Tabaei, A. & Akylas, T. R. 2007 Resonant long-short wave interactions in an unbounded
rotating stratified fluid. Stud. Appl. Math. 119 (3), 271–296.
Whitham, G. B. 1965 A general approach to linear and non-linear dispersive waves using a
Lagrangian. J. Fluid Mech. 22 (2), 273–283.

Schlutowetal 2017

Uploaded by

Copyright:

Available Formats

Schlutowetal 2017

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Schlutowetal 2017

Uploaded by

Copyright:

Available Formats

This draft was prepared using the LaTeX style file belonging to the Journal of Fluid Mechanics 1

Finite-amplitude gravity waves in the

(Received xx; revised xx; accepted xx)

Wentzel-Kramers-Brillouin (WKB) theory was employed by Grimshaw (Geophys. Fluid

2. The governing equations

(i)low Mach number

(X, Z, T ) = (εx, εz, εt) (2.10)

V0 (Z) and Vε(i) (x, z, t; ε) (2.11)

respectively. In addition to that we also define the “logarithmic background derivative”

3. WKB theory for finite-amplitude gravity waves with non-uniform

3.1. A spectral WKB expansion

3.2. Leading-order analysis

M(−iω̂, ik)Û0,1 = 0 (3.11)

3.3. First-order analysis

∂T û0,0 + ∂X P̂0,0 = −û∗0,1 ∂X û0,1 − ŵ0,1

M(−iω̂, ik)Û1,1 = −R(Û1,2 , ik)Û0,−1 − L(δT , ∇ε )Û0,1 (3.17)

3.4. Second-order analysis

3.5. Derivation of the modulation equations

for the wave action density A = E/ω̂.

∂T û0,0 + ∂X P̂0,0 = −∂X |û0,1 |2 − (∂Z + ηρ )(û∗0,1 ŵ0,1 ) + c.c., (3.34)

5. Traveling wave solutions of the modulation equations

5.1. The horizontally traveling wave solution

with LΨ some constant which may be interpreted as a slow wavenumber. A general

Here, we chose Cx = Nmin /Kx + Uf . It can be integrated by means of standard

The slow phase becomes simply arg(A) = LΨ Ψ .

5.2. The isothermal, horizontally periodic traveling wave

The critical asymptotic mass-specific wave action density acrit

Lref /km Hθ /km ε N Kx a− kz+ Cz

which is also shown in figure 3 from a two-dimensional perspective. We chose u0 = −Kx a−

q = Kx2 ρ0 a(kz )ω̂ 0 (kz )a0 (kz )∂Z kz < 0. (5.46)

Lref /km lx /Lref lz /Lref Nmin Nmax ∆ζ Kx LΨ ∆Ψ Gmax sgn(kz )

6. Numerical validation of the traveling wave solutions under the full

where Gmax is the maximum flux at the peak.

The numerical quadrature exploits the trapezoidal rule.

Stratification Initial condition Initial minus final condition

-25 -20 -15 -10 -5 -3 0 3

M = 2.93 ± 0.65 M = 3.20 ± 0.61 M = 3.09 ± 0.60

ln = CεM for n = 1, 2, ∞. (6.8)

(X, Z, T ) = (ε2 x, εz, ε2 t) (7.1)

and consider the following ansatz for the WKB approximation

(X, Z, T ) 7→ Ûn,m (X, Z, T ), Φ(X, Z, T ). (7.4)

You might also like