## Abstract

With the aim of contributing to the improvement of subgrid-scale gravity wave (GW) parameterizations in numerical weather prediction and climate models, the comparative relevance in GW drag of direct GW–mean flow interactions and turbulent wave breakdown are investigated. Of equal interest is how well Wentzel–Kramer–Brillouin (WKB) theory can capture direct wave–mean flow interactions that are excluded by applying the steady-state approximation. WKB is implemented in a very efficient Lagrangian ray-tracing approach that considers wave-action density in phase space, thereby avoiding numerical instabilities due to caustics. It is supplemented by a simple wave-breaking scheme based on a static-instability saturation criterion. Idealized test cases of horizontally homogeneous GW packets are considered where wave-resolving large-eddy simulations (LESs) provide the reference. In all of these cases, the WKB simulations including direct GW–mean flow interactions already reproduce the LES data to a good accuracy without a wave-breaking scheme. The latter scheme provides a next-order correction that is useful for fully capturing the total energy balance between wave and mean flow. Moreover, a steady-state WKB implementation as used in present GW parameterizations where turbulence provides by the noninteraction paradigm, the only possibility to affect the mean flow, is much less able to yield reliable results. The GW energy is damped too strongly and induces an oversimplified mean flow. This argues for WKB approaches to GW parameterization that take wave transience into account.

## 1. Introduction

The parameterization of gravity waves (GWs) is of significant importance in atmospheric global circulation models (GCM), in global numerical weather prediction (NWP) models, and in ocean models. In spite of the increasing available computational power and the corresponding increase of spatial resolution of GCMs and NWP models, for the time being, an important range of GW spatial scales remains unresolved both in climate simulations and in global NWP (Alexander et al. 2010). Numerous studies indicate that a representation of GWs is necessary for a realistic description of various aspects of the middle-atmospheric circulation [e.g., the Brewer–Dobson circulation (Butchart 2014)] and hence the zonal-mean winds and temperature (Lindzen 1981; Houghton 1978), the quasi-biennial oscillation (QBO) (Holton and Lindzen 1972; Dunkerton 1997), sudden stratospheric warmings (SSW) (Richter et al. 2010; Limpasuvan et al. 2012), and—via feedback loops—the tropospheric circulation [e.g., the North Atlantic Oscillation (Scaife et al. 2005, 2012)].

Parameterizations of the gravity wave drag are indeed applied in most GCMs or NWP models (Lindzen 1981; Medvedev and Klaassen 1995; Hines 1997a,b; Lott and Miller 1997; Alexander and Dunkerton 1999; Warner and McIntyre 2001; Lott and Guez 2013). In some way or other they all use the Wentzel–Kramer–Brillouin (WKB) theory, but with some important simplifications: that is, (i) the assumption of a steady-state wave field and background flow, (ii) the neglect of the impact of horizontal large-scale flow gradients on the GWs, and (iii) one-dimensional vertical propagation. Under these conditions, the wave-dissipation or nonacceleration theorem states that GWs can deposit their momentum only where they break. In theoretical analyses of this problem in a rotating atmosphere, Bühler and McIntyre (1999, 2003, 2005) point out that the steady-state assumption can lead to the neglect of important aspects of the interaction between GWs and mean flow. By wave-resolving numerical simulations and analyses on the basis of a nonlinear Schrödinger equation, Dosser and Sutherland (2011) have demonstrated the relevance of direct GW–mean flow interactions as well. Still missing, however, is an explicit assessment of the significance of the direct interaction between transient GWs and mean flow as represented by WKB, called direct GW–mean flow interactions in the following. WKB modeling for diagnostic purposes, as by the Gravity-Wave Regional or Global Ray Tracer (GROGRAT) model (Marks and Eckermann 1995; Eckermann and Marks 1997), is a well-established tool (e.g., Eckermann and Preusse 1999), but such analyses leave out the GW impact on the large-scale flow. A semi-interactive approach to studies of the interaction between GWs and solar tides has been described by Ribstein et al. (2015), however, with a simplified treatment of the GW impact on the solar tides, using effective Rayleigh-friction and thermal-relaxation coefficients. The numerical implementation of a fully interactive WKB theory, allowing direct GW–mean flow interactions, is not a trivial task and should best be accompanied by validations against wave-resolving data. In a Boussinesq framework, the representation of direct GW–mean flow interactions by a WKB algorithm has been studied by Muraschko et al. (2015) for vertically propagating idealized wave packets with variable vertical wavenumber. WKB theory had been implemented there in a two-dimensional phase space spanned by the physical height and the vertical wavenumber. The phase-space representation (Bühler and McIntyre 1999; Hertzog et al. 2002; Broutman et al. 2004) turned out to be effective to avoid numerical instabilities due to caustics (i.e., when ray volumes representing GWs become collocated in physical space but have different vertical wavenumbers and thus group velocities). Muraschko et al. (2015) could demonstrate the validity of their approach by comparisons against wave-resolving large-eddy simulation (LES) data.

The Boussinesq setting, however, leaves out the amplitude growth experienced by atmospheric GWs due to propagation into altitudes with decreasing density. That process, however, is central for the ensuing wave breaking due to static or dynamic instability at large GW amplitudes. By the wave-dissipation theorem, steady-state GW parameterizations depend on wave breaking as the mechanism leading to a large-scale GW drag. How this mechanism—the only one represented by present GW parameterizations—competes with GW drag by direct GW–mean flow interactions and how well the latter can be captured in the atmosphere by a WKB algorithm have remained mainly unanswered questions to date. These are addressed here by investigations in a non-Boussinesq atmosphere, where the WKB algorithm is supplemented by a wave-breaking scheme. Steady-state WKB simulations are considered as well, representing the GW parameterization approach in present weather and climate models. Like Muraschko et al. (2015), we consider idealized cases of upwardly propagating horizontally homogeneous GW packets. LES provide wave-resolving reference data.

Our investigations are described as follows: the theoretical background of the work is presented in section 2, while the corresponding numerical models are introduced in section 3. This is followed in section 4 by the presentation of the results. In section 5, the main findings of the work are summarized, and conclusions are drawn.

## 2. Theoretical background

We start out with the compressible two-dimensional Euler equations without rotation, which describe the evolution of the fluid in the *x*–*z* plane:

where *g* is the gravitational constant, denotes the heat capacity at constant pressure, *R* is the ideal gas constant, , is the Exner pressure with *p* pressure and surface pressure, and *u* and *w* are the velocity components in the *x*–*z* plane. We assume that the flow consists of a reference part constant in time as well as a large-scale background part and a small-scale wave part both changing in time.

### a. WKB theory

A notation for each variable in Eqs. (1)–(4) can be introduced as , where the first term denotes the reference part, with zero wind, while the other two refer to the large-scale background and the wave parts, respectively. The flow is assumed to be periodic in the *x* direction. By linearization of Eqs. (1)–(4) about the reference and large-scale background, introducing the wave pressure and using an appropriate WKB scaling (see an explanation below), one obtains the Boussinesq polarization and dispersion relations at first order:

and

where is the intrinsic frequency, and *N* denotes the Brunt–Väisälä frequency. In the polarization relation [Eq. (5)], *U*_{w}, *W*_{w}, *B*_{w}, and *P*_{w} denote WKB wave amplitudes of *u*_{w}, *w*_{w}, *b*_{w}, and *p*_{w}, where is the wave buoyancy.

By “appropriate WKB scaling” it is meant that a scale separation between the potential temperature scale height and the wavelength is assumed and that a corresponding WKB ansatz (Bretherton 1966; Grimshaw 1975; Achatz et al. 2010; Rieper et al. 2013a) is imposed

where *k* is a constant horizontal wavenumber, always assumed to be positive, and the local phase *ϕ* defines the local vertical wavenumber and the local frequency . The wave amplitude , the local frequency, and vertical wavenumber, like the large-scale background , depend only slowly on *z* and *t*.

The WKB approximation at the next order leads to the wave-action conservation equation

where is the vertical group velocity and is the wave-action density, with

the wave energy. From the definition of the vertical wavenumber and frequency via the local phase, one derives a prognostic equation,

for the vertical wavenumber. A solution method for the *field equations* [Eqs. (8) and (10)] is the ray technique, observing that along characteristic, so-called ray trajectories, defined by , wavenumber, and wave-action density satisfy the *ray equations*:

where the dispersion relation [Eq. (6)] is used to calculate the local intrinsic frequency and hence the group velocity. By definition there is a unique wavenumber and a unique frequency at each vertical location.

The system is closed by a prognostic equation for the mean flow. Based on Achatz et al. (2010) and Rieper et al. (2013a), for example, it is obtained as

with denoting the complex conjugate and thus

The problem with Eqs. (11)–(15) is that, after being initialized from some fields of *m*, , and , they very often lead to so-called caustics, where wavenumber, frequency, and wave-action density are not unique anymore. This happens when rays with different wavenumbers cross in space. Then the solution is not well defined anymore, and numerical instabilities become a serious problem (Rieper et al. 2013a) in attempts to obtain a local wavenumber by averaging the crossing rays. As demonstrated by Muraschko et al. (2015) in the Boussinesq context, however, this problem can be circumvented by considering the wave fields as a superposition of (infinitely) many WKB wave fields, characterized by a field index *β*, each having wavenumber and wave-action density and and satisfying Eqs. (11)–(15) separately. In phase space, spanned by wavenumber and position, here *m* and *z*, one introduces a wave-action density

with *δ* denoting the Dirac delta function. It can then be shown that

and since

also

In this representation, wavenumber is not a prognostic field but a coordinate. The only wave field to be predicted is . Moreover, because of Eq. (18), the phase-space flow is volume preserving so that rays cannot cross. Again, one can resort to a ray technique, now however in phase space. Defining the rays by their phase-space velocity and , with given by Eq. (10), one simply solves along these rays:

that is, one keeps the conserved phase-space wave-action density along ray trajectories. For diagnostic purposes one can also determine the superposition of constituting wave-action densities,

and the corresponding total wave energy density:

The ray equations are to be coupled to a mean-flow equation with a wave impact that is the superposition of the wave impact of each of the constituting wave fields, characterized by and ; hence,

and thus

In a nutshell, the GW field and the mean flow are coupled and have an impact on the time evolution of each other: the GW field is influenced by the mean flow via its impact on , and in turn the mean flow is modified by the GW phase-space wave-action density via Eq. (24). This direct coupling is clearly transient. It is nonlinear, but different spectral components can only interact indirectly with each other by GW–mean flow interactions. Rigorously this is only correct at small amplitudes so that one might speak of a weakly nonlinear theory. As will be seen below, however, it yields quite useful results even at large amplitudes, close to breaking. Transience also means that wave propagation is described in a prognostic manner, different from present-day steady-state GW parameterizations, where everything is instantaneous.

### b. Wave breaking

As WKB theory does not account for the possible turbulent wave breakdown at large GW amplitudes, the coupled ray and mean flow equations above have been supplemented with a saturation criterion in an attempt to add additional parameterization of this process. Comparisons between simulations with or without this “turbulence scheme” also enable an assessment of the relevance of wave breaking for the GW drag, as compared to the direct GW–mean flow interactions described by WKB.

It is assumed that saturation occurs if static instability sets in at a certain height *z* during the wave propagation (Lindzen 1981) so that somewhere within a complete wave cycle , or, after an additional multiplication by ,

Comparison with Eq. (7) shows that this occurs in a locally monochromatic GW field if , or, using Eq. (9),

We transfer this from the locally monochromatic situation to the spectral treatment represented by the phase-space approach by taking Eq. (22) into consideration, suggesting

The free parameter *α* represents well-known uncertainties of the criterion Eq. (25). Stability analyses (e.g., Lombard and Riley 1996; Achatz 2005) and direct numerical simulations (Fritts et al. 2003, 2006; Achatz 2007; Fritts et al. 2009) indicate that GWs are unstable already below the static-instability threshold, and strongly nonhydrostatic, modulationally unstable wave packets also tend to break earlier (Dosser and Sutherland 2011). Another issue is that the criterion does not account for the possibility of destructive interference of different spectral components that would retard the onset of static instability.

Once the static-instability criterion is satisfied at height *z*, turbulence is assumed to be generated that acts to damp wave-action density to an extent that the GW field becomes again statically stable. Following Lindzen (1981) and Becker (2004), the turbulent fluxes are modeled by eddy viscosity and diffusivity so that small scales are damped more strongly than larger scales. The buoyancy equation, for example, is supplemented by a diffusion term

with the turbulent eddy diffusivity coefficient . By Fourier transformation in space and integration over a short time interval , one obtains the following as change of the buoyancy amplitude:

Employing identical eddy viscosity and diffusivity, an analogous equation,

can be derived for the wave amplitude. Hence, after a saturation step

and thus the turbulent eddy diffusivity is computed as

Further details regarding the numerical implementation of this wave-breaking parameterization are discussed in section 3c.

In summary, the weakly nonlinear coupled GW–mean flow equations [Eqs. (11), (12), (20), and (24)] describe the time evolution of a *transient* GW field through a *transient* large-scale background flow in a *direct* manner. In addition, wave breaking is accounted for in the WKB models by applying the saturation criterion [Eq. (27)] and reducing the wave-action density proportionally to , as prescribed in Eq. (30), if necessary.

### c. Steady-state WKB theory

As mentioned in the introduction, current GW parameterization schemes are based on a steady-state WKB theory (Nappo 2002; Coiffier 2011; Fritts and Alexander 2003; Kim et al. 2003). The assumption of a steady wave-action density profile reduces Eq. (8) to

Hence, the pseudomomentum flux is altitude independent, and the GW drag in Eq. (15) vanishes. This is the nonacceleration paradigm. It is the reason why steady-state WKB schemes rely on wave breaking, thereby imposing a nonzero pseudomomentum-flux convergence and hence tendencies for the induced wind. To compute the equilibrium profile of the wave field, first the vertical group velocity profile is obtained via Eq. (11) from a vertical wavenumber profile:

where is the constant extrinsic frequency, with a “source” altitude where vertical wavenumber and wave-action density are prescribed. From Eq. (33), one then obtains the wave-action density profile:

Wave breaking is assumed wherever the static-instability condition Eq. (26) is fulfilled, which amounts to setting the wave-action density profile there to

using the same *α* uncertainty parameter as explained in section 2b. Notably, this approach leads to instantaneous pseudomomentum-flux profiles. Variations of the boundary conditions at the source altitude are communicated immediately throughout the whole altitude range of a model, while in a more realistic transient approach any signal propagates at the group velocity. There are various possibilities of implementations of steady-state parameterizations (Fritts and Alexander 2003; Alexander et al. 2010): for example, by allowing a spectrum via a superposition of components as just described, each with its own values of vertical wavenumber and wave-action density at the source altitude. Lott and Guez (2013), for example, suggest picking these in a stochastic manner from a random sample. However, all of these approaches are instantaneous, and they only allow GW–mean flow interactions where GWs break.

## 3. Test cases and numerical models

Simulations have been done for a set of idealized test cases, where horizontally homogeneous quasi-monochromatic GW packets are initialized in an isothermal background with a reference temperature resulting in a constant buoyancy frequency . This implies a reference density profile

where is the density scale height. Some of the test cases involve a prescribed background jet as an initial mean flow with a half-cosine wave shape:

where is the maximal magnitude of the jet initialized at height , and is the width (i.e., vertical extent) of the half-cosine shape. In these cases, the wave-induced mean flow is diagnosed as : that is, the initial mean wind is subtracted from the total mean wind to get the one induced by the GW. We remark in this context that integrating Eq. (17) in wavenumber space, assuming a vanishing wave-action density flux at the boundaries, and multiplying the result by the constant horizontal wavenumber yields, without saturation scheme,

Therefore, comparing with Eq. (24), one obtains

so that is in the absence of wave breaking the residual between , often termed the wave-induced wind, and its initial value.

The GW packets are initialized with a Gaussian- or a cosine-shaped buoyancy amplitude envelope in the vertical direction: that is,

where is the height of the wave envelope maximum, is the initial vertical wavenumber, and is the initial amplitude factor implying static instability if . The parameter *σ* defines the vertical size of the GW packet : namely, for the Gaussian wave packet and for the cosine-shaped wave packet. The envelope of the cosine-shaped wave packets is limited to the interval [i.e., outside this vertical range]. In the horizontal *x* direction, the wave packet is initialized with infinite extent and a constant wavenumber *k*. To initialize the idealized wave packet in the wave-resolving LES, the following perturbations are prescribed at initial time :

In the transient WKB simulations (we introduce this terminology for the nonsteady-state WKB simulations) the GW packets are initialized via the corresponding monochromatic phase-space wave-action density : that is,

As a numerical representation of Eq. (45), the initial phase-space wave-action density is set as

where is a narrow initial wavenumber width of the wave packet. A typical value of the initial ratio in our numerical experiments is .

Seven idealized test cases have been investigated. Three cases elaborate the refraction and the reflection of hydrostatic GW packets from a background jet, while four other cases aim to study static and modulational instability of hydrostatic and nonhydrostatic wave packets, including the process of wave breaking. The initial wave packet characteristics , and vary from case to case as well as the magnitude and height of the jet. In all cases, the negative frequency branch of Eq. (6) has been used so that positive vertical wavenumbers correspond to upward-directed group velocities. For the specific settings for each case see Tables 1 and 2 and the corresponding explanations in section 4. The LES resolution is , *dz* ≈ λ_{z0}/30, while the WKB simulations have been done at a vertical resolution of *dz* ≈ λ_{z0}/10: that is, at a resolution 3 times coarser than the reference LES (see further details in Tables 1 and 2). Both LES and WKB simulations with an increased resolution have been performed without observing significant changes in the results, which confirms that a convergence in the numerical results has been reached with the resolution described above.

### a. Reference LES model

The reference LES model called PincFloit solves the pseudo-incompressible equations [i.e., a soundproof approximation of the Euler equations [Eqs. (1)–(4)] (Durran 1989)]. A third-order Runge–Kutta time scheme and a finite-volume spatial discretization are applied, which involve an adaptive local deconvolution model (ALDM) (Hickel et al. 2006) for turbulence parameterization. Alternatively, the Monotonic Upstream-Centered Scheme for Conservation Laws (MUSCL; van Leer 1979) can also be used in the finite-volume scheme. Tests using both schemes for our cases did not show an important sensitivity. It is important to mention that, in contrast to the WKB simulations, the reference LES is fully nonlinear and enables the description of wave–wave interactions as well as turbulent wave dissipation, which, with a high resolution, implies a realistic description of compressible flows. The PincFloit model has been described and tested in detail by Rieper et al. (2013b).

### b. Eulerian WKB model

The Eulerian implementation of the WKB equations solves the flux form [Eq. (17)] of the phase-space wave-action density equation using the MUSCL finite-volume discretization on the *z*–*m* plane, with an equidistant staggered grid in both the *z* and *m* direction. As a consequence, the phase-space wave-action density and the derivatives and are defined at cell centers, while and are defined at cell edges, as are the rest of the variables, including the wave energy . In addition, the mean flow Eq. (24) is solved using simple centered differences on the vertical part of the same staggered grid. A fourth-order Runge–Kutta time scheme is used to evolve all the prognostic variables. A detailed description of the model, there in Boussinesq approximation, is given by Muraschko et al. (2015).

### c. Lagrangian WKB model

The Lagrangian implementation of the WKB equations solves the advective form [Eq. (19)] of the phase-space wave-action density equation. This is done using a ray technique in phase space. Defining the ray velocities by and , with given by Eq. (10), one simply solves Eq. (20) along these: that is, one keeps the conserved phase-space wave-action density. This procedure is discretized numerically by gathering rays in finite ray volumes around a characteristic carrier ray each, with uniform phase-space wave-action density (see Fig. 1). By Eq. (20) that uniformity is conserved. Because the phase-space velocity is divergence free, each ray volume moreover preserves its volume content in phase space, but arbitrary shape deformations are possible. In a second discretization step, we constrain each ray volume, however, to keep a rectangular shape, responding nonetheless, in a volume-preserving manner, to local stretching and squeezing. The prognostic equation for the evolution of the ray-volume edge length in *m* is given by

with and being the rectangle edges in the *m* direction. Because of the conserved area *a*, the evolution is given through Eq. (47) as well. The wave energy and the wave-induced mean wind are computed on an equidistant vertical grid, which is staggered in a consistent manner with the Eulerian WKB model regarding the background variables, the momentum flux, and wave energy. The prediction of the mean flow is done as in the Eulerian model. Given that the distribution of the rays might get sparse in the vertical during the time evolution of the wave field, the projection of the momentum flux and the energy from the rays to the vertical grid is supplemented by a smoothing. This consists of computing the average over a certain number of neighboring vertical layers (the corresponding values are indicated for all numerical experiments in Tables 1 and 2). See Muraschko et al. (2015) for further specifications of the model. In comparison to that study, the reference density profile equation [Eq. (37)] has been implemented and used in the mean-flow equations [Eqs. (23) and (24)]. In line with this modification, the originally periodic bottom and upper boundary conditions have been changed to allow for a free outflow through these boundaries. This change was necessary, given the realistic growth of wave amplitudes, which are due to the quasi-realistic density profile, and which are essentially nonperiodic in vertical.

The wave-breaking parameterization has been implemented only in the Lagrangian model since it is much more efficient than its Eulerian counterpart (see section 4a), and thus this is the WKB variant intended for future numerical studies. In the Lagrangian framework, the analytical criterion equation [Eq. (27)] is to be rewritten as

Here, is the number of ray volumes overlapping with layer *i* of the Lagrangian WKB model, with and with height , *nz* being the number of vertical levels, and the vertical resolution (layer depth), with the height of the model top. The quantities and are the wavenumber and squared wave buoyancy amplitude, respectively, of the carrier ray of the *j*th ray volume relevant to layer *i*, and is the Brunt–Väisälä frequency value in layer *i*. The term is computed analytically from

where and are the edges of the *j*th ray volume in the *m* direction so that . The factor denotes the fraction of the *j*th ray being in the *i*th vertical layer, as ray volumes might overlap with several vertical layers (see Fig. 1). Using the dispersion relation [Eq. (6)], one obtains the following from Eq. (49):

### d. Steady-state WKB model

A numerical model based on the steady-state WKB theory (section 2c) has been implemented as well, in order to enable a comparison with the transient WKB simulations and thus an assessment of present-day GW parameterizations. For optimal correspondence between transient and steady-state simulations, the time-dependent boundary values for wavenumber and wave-action density at the “source” altitude in the steady-state simulation have been set to the corresponding values diagnosed from the transient simulation.

## 4. Results

We first give a comparative validation of the Eulerian and the Lagrangian WKB models, demonstrating the superiority of the latter. We then discuss a case where no wave breaking is active but where the negligence of direct GW–mean flow interactions would make a fundamental difference. Finally, we compare the relative importance of direct GW–mean flow interactions and of wave breaking in cases where both are active. There we also demonstrate the limitations of a steady-state approach with wave-breaking scheme.

### a. Comparative validation of the Eulerian and Lagrangian WKB models

The refraction of a hydrostatic GW packet () by a jet has been studied using a conventional WKB ray tracer in physical space by Rieper et al. (2013a). This WKB model failed to reproduce the LES simulation [see Figs. 11a and 11b in Rieper et al. (2013a)] and crashed as a result of numerical instabilities. These numerical instabilities were due to caustics (i.e., to rays crossing in physical space). A deliberate application of the phase-space representation to avoid numerical instabilities due to caustics was the main innovation of Muraschko et al. (2015); however, in their study a Boussinesq reference atmosphere has been used, which did not allow an investigation of the case referred to in Rieper et al. (2013a). Thus, as a first proof of concept for the phase-space approach in an atmosphere-like configuration with a variable reference density, the very case described in Rieper et al. (2013a) has been revisited via LES and the transient phase-space WKB model simulations. A Gaussian-shaped hydrostatic GW packet is initialized at an altitude of 10 km, propagating upward and refracted by a low-speed jet () with its maximum at 25 km (see case REFR in Table 1). The time evolution of the induced wind profiles is shown in Fig. 2, which reveals a good correspondence between the transient WKB simulations and the LES (shaded colors). This proves that the phase-space approach applied in both WKB models helps in avoiding numerical instabilities due to caustics and producing a realistic induced mean flow. Wave energy diagnostics of the transient WKB simulations compare to LES very well too (not shown).

If the jet blows in the appropriate direction, in case of the negative frequency branch, by increasing the speed of the jet, the refraction of the wave packet turns into reflection because of the strong vertical wind shear and thus the intensive tendency in the vertical wavenumber [see Eq.(10)]. A linear estimate for the jet speed threshold for reflection is [see, e.g., the work of Muraschko et al. (2015)]

which, in the current case, gives . To achieve a reflection with full certainty and also to validate the models under strong-gradient conditions, was chosen for our next case, REFL (see details in Table 1). Note that reflection also implies caustics and thus could not be properly handled with a conventional ray tracer. Also, reflection is a great challenge for the WKB theory in itself because it implies a change of sign of the wavenumber and thus an increase of the vertical wavelength far over the envelope scale, which breaks the scale separation assumption. In spite of these challenges, the Lagrangian transient WKB simulation provides results with a good agreement with the LES, which is demonstrated in Figs. 3a and 3b by plotting the time evolution of the wave energy profile. It is apparent, however, in Fig. 3c that the Eulerian transient WKB simulation does not produce a satisfactory reflection. This turned out to be due to the strong vertical wind shear, which made the Eulerian WKB code too diffusive, and thus almost fully dissipating the GW packet near the reflection level. By increasing the resolution of the Eulerian model by a factor of 10 in both the *z* and *m* directions, a good agreement with the LES can be achieved (not shown); however, at the same time the computational time gets exhaustive (i.e., 2–3 times larger than the computing time of the LES). It is to be mentioned here that, even by using the same vertical resolution in the Eulerian and the Lagrangian transient WKB simulations, the latter is more efficient computationally than the former, by a factor of 10–100, depending on the number of rays used in the Lagrangian model. The better efficiency of the Lagrangian model is due to the fact that, in this framework, (i) there is no necessity to span the whole phase-space volume, including all its regions where the wave-action density is negligibly small, and (ii) the prediction of the latter is done by solving the trivial conservation equation [Eq. (20)]. In contrast, in the Eulerian model, the prognostic equation [Eq. (17)] is solved using the MUSCL finite-volume scheme, which requires a relatively expensive reconstruction of the fluxes on the cell edges. These findings regarding the efficiency of the transient WKB simulations motivated the use of the Lagrangian model in all our further studies.

### b. Role of direct GW–mean flow coupling

Initializing a hydrostatic GW packet with somewhat shorter horizontal wavelength () and somewhat larger vertical wavelength () results in a reflection threshold based on Eq. (52). An interesting experiment is to see whether the reflection happens if a jet speed maximum of is set (see case PREFL in Table 1). As shown in Fig. 4a, the GW packet in the LES is only partially reflected and a part of the wave packet is just refracted by the jet. This is due to the fact that Eq. (52) is a linear estimate for while the wave–mean flow interaction is a nonlinear process, as described by the LES. This shows that the linear approximation equation [Eq. (52)] is underestimating the jet speed threshold of reflection in a realistic nonlinear flow. The partial reflection (i.e., the behavior of the LES) is qualitatively reproduced by the Lagrangian WKB model (Fig. 4b), which suggests that the weakly nonlinear WKB dynamics are successful in capturing important aspects of the nonlinear interactions between the wave and the mean flow. What happens is that the GW packet based on Eq. (23) induces a wind, which is comparable to the jet speed but blowing in the opposite direction () and thus reduces the net velocity of the jet. Indeed, running the Lagrangian WKB model in a decoupled mode [i.e., ignoring Eq. (23) and thus not permitting the GW packet to modify the mean flow] the wave packet is fully reflected as predicted by the linear estimate [Eq. (52)] (see Fig. 4c). An additional experiment points out the importance of the variable reference density profile: the Lagrangian WKB model, in its coupled mode, predicts full reflection if a Boussinesq reference medium is used (Fig. 4d). In this case, the weakly nonlinear wave–mean flow interaction is active, but, in the absence of amplitude growth due to the density effect, it is too weak to induce winds that are strong enough to alter the jet significantly. All this shows that, during the propagation of GWs in an atmosphere-like medium from their source till their breaking levels, they continuously interact with their background flow and modify it significantly. In addition to these wave-induced changes in the background flow, the waves themselves can switch to inherently different regimes (e.g., changing the direction of propagation). No wave breaking occurs in this case. Because of the nonacceleration paradigm, steady-state WKB simulations would therefore not be able to reproduce the observed dynamics.

### c. Role of wave breaking

Finally, a set of cases with unstable GW packets or with GW packets turning into unstable regimes has also been studied. Comparisons between LES and transient WKB simulations with wave-breaking parameterizations serve to validate the latter. More importantly, however, these cases are to provide an assessment of the comparative importance of direct GW–mean flow interactions, as represented by the WKB model without wave-breaking parameterization, and the wave-breaking process itself, as represented by that latter turbulence scheme.

All of the cases are discussed in terms of the simple energetics that arises in WKB theory without rotation for horizontally homogeneous GW packets. Taking the time derivative of Eq. (22) and using the flux-form wave-action Eq. (17), one obtains

Inserting from Eq. (10) and assuming a vanishing wave-action density flux at the boundaries yields

From Eq. (24), however, one finds that the time derivative of the mean-flow kinetic-energy density

is

so that

and hence the sum of wave energy and mean-flow kinetic energy is conserved if the fluxes vanish at the vertical model boundaries. It thus makes sense to consider the vertically integrated energy densities:

and their sum . To do so, we visualize normalized values

In the first case, a Gaussian-shaped hydrostatic wave packet () travels upward and becomes statically unstable during the course of its evolution (see case STIH in Table 2). The results for this case are shown in Figs. 5a–d in terms of normalized integrated energy. The wave-breaking effect can be recognized in Fig. 5a at a decay of the total energy that is not visible in the results from the transient WKB simulation in Fig. 5b. Switching on the saturation scheme with (Fig. 5c), the GW energy, and hence also the total energy, gets reduced earlier than in the LES and also results in too weak an induced mean flow in the end. The total energy and the mean flow energy can be brought into better agreement with the LES by using the saturation scheme with (see Fig. 5d). The value suggests that the static-instability criterion as applied in this study is too strict (i.e., mimics wave breaking too early/strongly). This is presumably because of the neglect of the wave phase in the saturation scheme. Finally, by looking at Fig. 5a, it is apparent that, except for the uppermost 10 km in the induced wind, the transient WKB simulation reproduces the LES vertical structures already relatively well even without the wave saturation parameterization.

The evolution of a Gaussian-shaped nonhydrostatic GW packet () is discussed next, which evolves quickly into a statically unstable regime, because of its high initial amplitude factor (case STINH in Table 2). The normalized energies of the LES in Fig. 6a imply a decay of total energy that saturates by . This is not reproduced completely by the transient WKB simulation (Fig. 6b). By switching on the wave-breaking parameterization, however, with the LES results are met rather well (Fig. 6c). As a reference, the results on the energetics from the transient WKB simulation with saturation scheme and are also plotted in Fig. 6c (dashed curves). Again, the dissipation is somewhat too strong, indicating that the static-instability criterion is too strict for best diagnosing saturation in the real atmosphere. The reduced optimal value of *α* compared to the previous case might be due to the added inclination of nonhydrostatic wave packets to become modulationally unstable. By looking at the vertical structures of the wave energy and the induced mean wind in Figs. 6a–l, one finds again that the wave-breaking parameterization provides only small corrections on top of the relatively good results provided by the transient WKB simulations without the saturation scheme.

The next case involves a cosine-shaped nonhydrostatic GW packet (), which becomes modulationally unstable during its evolution: that is, its vertical wavelength grows beyond its horizontal wavelength so that (see case MI in Table 2). In this regime, the wave-induced mean flow accelerates the trailing edge and decelerates the leading edge of the wave envelope while the amplitude grows, leading to the collapse of the wave packet because of local static instability (Sutherland 2006, 2010; Dosser and Sutherland 2011). This same case has also been studied by Rieper et al. (2013a), where their conventional WKB code broke down because of numerical instabilities related to caustics, as explained in section 4a. In contrast to that study, the Lagrangian phase-space WKB model remains numerically stable and reproduces the LES results relatively well (Figs. 7a–i), with or without the wave-breaking parameterization. The integrated-energy plots in Figs. 7a–c suggest, however, that without the saturation scheme the induced mean flow is overestimated by the WKB model, and a best fit to the LES is found if is used in the saturation scheme. The results with (Fig. 7c, dashed curves), however, are also quite acceptable.

Finally, a hydrostatic GW packet () reaching a critical layer is studied. An easterly jet with a maximum of is prescribed at so that without wave impact on the mean flow the intrinsic phase velocity would vanish at around 22–23-km height (see case CL in Table 2). The transient WKB simulation without saturation scheme seems to slightly overestimate the mean flow energy compared to the LES, (Figs. 8a,b), which can be removed by switching on the saturation scheme with (Fig. 8c). This value of *α* suggests that, in case of a critical layer, GWs tend to break as predicted by classic static-instability criteria.

Our results suggest that wave breaking is of secondary importance in comparison with the direct GW–mean flow interactions even for large-amplitude GWs. Since present-day GW-drag parameterizations exclusively rely on wave breaking, the question arises what results they would yield in the cases considered. Therefore, the cases STIH, STINH, and MI have been also been simulated using the steady-state WKB model based on sections 2c and 3d. Figure 9 shows the corresponding results for case STINH, which are to be compared with Fig. 6, where the transient WKB and the LES results are shown. The integrated energy shows that the wave energy is overdamped and that the kinetic energy of the mean flow is strongly underestimated in the steady-state WKB model. The former is also observed in the Hovmöller diagram of the wave energy. The Hovmöller diagram of the induced mean wind shows that the magnitude of the GW drag is too small in the steady-state WKB model, and also its structure is very different from that of the LES and the transient WKB simulation. One should, of course, restrict the comparison between the models to the vertical region above the source (), but there as well the results from the steady-state WKB simulation show an unrealistic structure in the induced mean flow, seemingly fully determined by wave breaking. Again, this demonstrates the dominant role of direct GW–mean flow interactions as compared to wave breaking, and it also points to limitations of present-day GW parameterizations.

## 5. Summary and conclusions

The steady-state approximation to WKB theory used nowadays in GW-drag parameterizations implies that the only GW forcing on the mean flow is due to wave breaking. Transient GW–mean flow interactions can, however, act as another important coupling mechanism. This study provides an assessment of the comparative importance of these processes in typical atmospheric situations, albeit idealized. Focusing on single-column scenarios for the time being, considered GW packets are horizontally homogeneous, and the mean flow has only a vertical spatial dependence. The wave scales and amplitudes, however, are representative, although not of inertia–gravity waves affected by rotation. Fully interactive transient WKB simulations are used to describe the simultaneous development of GWs and mean flow. All of these simulations are validated against wave-resolving LES, thereby assessing the reliability of the methods employed.

The WKB algorithms used allow the simulation of transient GW development. In both variants, Eulerian or Lagrangian, the mean flow is fully coupled to the wave field. This is enabled by a spectral approach, employing wave-action density in position-wavenumber phase space, the key to avoiding otherwise detrimental numerical instabilities due to caustics. The Eulerian approach spans the whole phase space. It thus quickly tends to be expensive, often more so than LES. The Lagrangian ray-tracing approach, however, focuses on regions of phase space with nonnegligible wave action. This makes it considerably more efficient, by orders of magnitude, than the wave-resolving simulations. Certainly this might change in situations where broad spectra develop. So far, however, we have not met with such a case.

A systematic investigation of the comparative relevance of wave breaking, as compared to direct GW–mean flow interactions, has been enabled by the implementation of a simple turbulence scheme. Turbulence is invoked whenever the wave field has the possibility to become statically unstable. A flux-gradient parameterization of turbulent fluxes is used, by way of eddy viscosity and diffusivity. The ensuing damping of the GW field is hence scale selective so that small scales are damped more strongly. Generally, it is found, by comparison against the LES data, that the static-stability criterion tends to generate turbulence too quickly. This might be explained by phase cancellations between different spectral components so that higher amplitudes are required to really lead to the onset of turbulence. Nonetheless, the turbulence scheme works quite well if validated against the LES simulations.

Finally also a steady-state WKB model has been implemented, representing the approach in current GW-drag parameterizations. Caustics are not an issue in such a context so that a spectral formulation is not necessary. The simulations discussed here consider locally monochromatic GW fields that are damped to the saturation limit once static instability is diagnosed. Spectral extensions have been considered as well (not shown), but these did not yield substantially different results.

The pattern observed in our simulations is quite clear: whenever a wave impact on the mean flow is observed, the direct GW–mean flow interactions dominate over the wave-breaking effect. It is important that these interactions depend on wave transience. Without the latter they would not be possible—because of the nonacceleration paradigm—without onset of turbulence. Interesting effects arise that would be missed by a steady-state scheme. An example is partial reflection from a jet that would not occur within linear theory. The wave-induced wind contributes sufficiently so that part of the GW packet substantially changes its propagation. In the various turbulent cases considered, be it apparent wave breaking by direct static instability or triggered by modulational instability, we see a dominant impact from direct GW–mean flow interactions as well. Even when the turbulence scheme is suppressed, the results between transient WKB and LES agree to leading order, at least as far as the spatial distribution of wave energy and mean wind are concerned. Turbulence acts to next order and ensures the correct dissipation of total energy.

Turbulence without direct GW–mean flow interactions, however, fails to explain the LES data: the steady-state WKB simulations exhibit way too strong damping of the GW energy, and also the mean flow is underestimated significantly. This argues for a serious attempt at including GW-transience effects into GW-drag parameterizations. Undoubtedly, the implementation of transient WKB into climate and weather models will face considerable efficiency issues. The strong role we see for direct GW–mean flow interactions could provide motivation, however, to overcome these.

A final caveat might be that, notwithstanding the apparent success of transient WKB seen here, there are plenty of cases where this approach might also find its limitations. Not only have we neglected lateral GW propagation, an effect known to also be potentially important (Bühler and McIntyre 2003; Senf and Achatz 2011; Ribstein et al. 2015), but the limits of WKB as a whole will be reached where strong wave–wave interactions come into play, or where significant small-scale structures are present (Fritts et al. 2013).

## Acknowledgments

U.A. and B.R. thank the German Federal Ministry of Education and Research (BMBF) for partial support through the program Role of the Middle Atmosphere in Climate (ROMIC) and through Grant 01LG1220A. U.A. and J.W. 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, and AC 71/10-1.

## REFERENCES

*Internal Gravity Waves*. Cambridge University Press, 377 pp.

## Footnotes

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