## Abstract

This paper compares two different approaches for the efficient modeling of subgrid-scale inertia–gravity waves in a rotating compressible atmosphere. The first approach, denoted as the pseudomomentum scheme, exploits the fact that in a Lagrangian-mean reference frame the response of a large-scale flow can only be due to forcing momentum. Present-day gravity wave parameterizations follow this route. They do so, however, in an Eulerian-mean formulation. Transformation to that reference frame leads, under certain assumptions, to pseudomomentum-flux convergence by which the momentum is to be forced. It can be shown that this approach is justified if the large-scale flow is in geostrophic and hydrostatic balance. Otherwise, elastic and thermal effects might be lost. In the second approach, called the direct scheme and not relying on such assumptions, the large-scale flow is forced both in the momentum equation, by anelastic momentum-flux convergence and an additional elastic term, and in the entropy equation, via entropy-flux convergence. A budget analysis based on one-dimensional wave packets suggests that the comparison between the abovementioned two schemes should be sensitive to the following two parameters: 1) the intrinsic frequency and 2) the wave packet scale. The smaller the intrinsic frequency is, the greater their differences are. More importantly, with high-resolution wave-resolving simulations as a reference, this study shows conclusive evidence that the direct scheme is more reliable than the pseudomomentum scheme, regardless of whether one-dimensional or two-dimensional wave packets are considered. In addition, sensitivity experiments are performed to further investigate the relative importance of each term in the direct scheme, as well as the wave–mean flow interactions during the wave propagation.

## 1. Introduction

As one of the most fundamental physical modes in meteorology, gravity waves (GWs) are ubiquitous buoyancy oscillations in the atmosphere. The sources of excited gravity waves include, among others, topographic forcing (Smith 1980; Menchaca and Durran 2017), convection (Alexander et al. 1995; Lane et al. 2001), the jets (Zhang 2004; Plougonven and Zhang 2014; Hien et al. 2018), frontal systems (Snyder et al. 1993; Griffiths and Reeder 1996), and shear instability (Bühler et al. 1999; Bühler and McIntyre 1999). GW dynamics is also influenced by rotation, especially when the respective waves have relatively long horizontal wavelengths and short vertical wavelengths. These are called inertia–gravity waves (IGWs), and they are characterized by parcel oscillations that are elliptical instead of straight lines in the pure gravity wave case without the Earth rotation effect (Gill 1982; Nappo 2002). In addition to that, IGWs have relatively low intrinsic frequencies and small vertical group velocities. Many studies have documented the signals and the life cycles of IGWs associated with complex background flows, using observations (Wang and Geller 2003; Plougonven et al. 2003; Gong et al. 2012), numerical investigations of observed cases (Zhang et al. 2013, 2015), and idealized simulations (Wei and Zhang 2014, 2015; Wei et al. 2016).

GWs play a significant role in atmospheric dynamics on various spatial and temporal scales. For example, GWs can generate and modulate atmospheric turbulence (Shapiro 1980; Lane et al. 2004), and they can also initiate and organize moist convection (Zhang et al. 2001; Lane and Zhang 2011). As an important candidate, GWs can also potentially contribute significantly to building the atmospheric energy spectra (Callies et al. 2014; Bierdel et al. 2016; Sun et al. 2017). GWs are found to likely link the small-scale small-amplitude initial error to the rapid upscale error growth and thus limit mesoscale predictability (Zhang et al. 2007; Sun and Zhang 2016; Bierdel et al. 2018). The most important impact globally is because GWs can travel over large distances from their sources and transfer significant amounts of momentum and energy to high altitudes, which contributes to the forcing of the circulation and the variability of the middle atmosphere (Holton and Lindzen 1972; Houghton 1978; Lindzen 1981; Dunkerton 1997; Richter et al. 2010; Limpasuvan et al. 2012; Butchart 2014). The dynamics of the middle atmosphere can influence the tropospheric circulation by downward control (Haynes et al. 1991), and it can be very important for the forecasting of weather (Baldwin and Dunkerton 2001) and climate (Scaife et al. 2005, 2012).

Despite the increasing computational power, an important range of GW spatial scales remains unresolved in most atmospheric global circulation models (GCM) or in global numerical weather prediction (NWP) models (Alexander et al. 2010), and GW parameterizations are still applied in those models (Medvedev and Klaassen 1995; Hines 1997a,b; Lott and Miller 1997; Alexander and Dunkerton 1999; Warner and McIntyre 2001; Lott and Guez 2013). Most of the parameterization schemes are based on the Wentzel–Kramer–Brillouin (WKB) theory, but with some oversimplifications, including the assumption of a steady-state wave field and background flow, instantaneous GW propagation, and one-dimensional vertical propagation. Those abovementioned oversimplifications are often made to ensure numerical stability and model efficiency, but they can lead to the neglect of important aspects of the interaction between GWs and mean flow (Bühler and McIntyre 1999, 2003, 2005; Dosser and Sutherland 2011; Bölöni et al. 2016). The spectral phase-space representation (Bühler and McIntyre 1999; Hertzog et al. 2002) turns out to be effective to avoid numerical instabilities due to caustics. The issue of caustics occurs when by wave–mean flow interactions a wave field loses its initial locally monochromatic character (Tabaei and Akylas 2007; Rieper et al. 2013a). Based on the spectral approach, a prognostic Lagrangian WKB GW ray-tracing model has been implemented, using a multidimensional phase space spanned by physical space and the wavenumber space, and it has been successfully validated against wave-resolving large-eddy simulation (LES) data in different idealized settings, including one-dimensional (1D) vertically propagating idealized wave packets with variable vertical wavenumbers in a nonrotating Boussinesq atmosphere (Muraschko et al. 2015), 1D wave packets propagating in a nonrotating compressible atmosphere (Bölöni et al. 2016), and two-dimensional (2D) wave packets of parameterized submesoscale GWs interacting with resolved mesoscale GWs in a rotating Boussinesq atmosphere (Wilhelm et al. 2018). In addition, this approach has also been used to study the interaction between GWs and solar tides (Ribstein et al. 2015; Ribstein and Achatz 2016).

However, there is still no corresponding detailed validation of an IGW parameterization in a rotating compressible atmosphere against data from idealized wave-resolving simulations. This gives one of the motivations for the current study. More importantly, there are actually two available approaches for IGW parameterization, here called pseudomomentum scheme and direct scheme, respectively. The pseudomomentum scheme exploits the fact that in a Lagrangian-mean reference frame the effect of GWs on the large-scale flow only appears in the momentum equation (Andrews and McIntyre 1978). Application of this theory to an Eulerian-mean reference frame leads to a pseudomomentum-flux convergence by which the large-scale momentum is to be forced (Andrews and McIntyre 1976, 1978). As will be shown below this is at least justified if the large-scale flow is in geostrophic and hydrostatic balance. The direct scheme does not rely on any balance assumption with regard to the large-scale flow, and the large-scale flow is forced by anelastic momentum-flux convergence in the momentum equation, an elastic term also in the momentum equation, and entropy-flux convergence in the entropy equation, as given by Grimshaw (1975) and Achatz et al. (2017). All present-day operational IGW parameterizations represent, one way or other, simplified versions of the pseudomomentum approach, where the vertical gradient of pseudomomentum-flux convergence forces the resolved flow, when wave dissipation occurs (Fritts and Alexander 2003; Kim et al. 2003), and neither elastic nor thermal effects are taken into account. In contrast, to the best of our knowledge, there is still no numerical application of the direct scheme. This further motivates us to understand the differences between these two available schemes, to investigate which scheme performs better in validations against data from idealized wave-resolving simulations, and to verify how trustworthy the pseudomomentum approach is when the large-scale flow is unbalanced.

This article is arranged as follows. The theory part related to the two available approaches for the IGW parameterization will be shown in section 2, followed in section 3 by a brief introduction to the numerical models, and in section 4 by a presentation of the various cases used. Section 5 will present the budget analysis of the wave-induced forcing terms from the two different schemes in different idealized wave packet profiles. Based on the prognostic Lagrangian WKB GW ray-tracing model, both schemes will be verified against the data from wave-resolving simulations in section 6. In section 7, sensitivity experiments are performed to further investigate the relative importance of each wave-induced forcing term, as well as the wave–mean flow interactions during the wave propagation. Section 8 contains a summary.

## 2. Theory

For an explanation of the theoretical underpinnings of the two respective approaches we follow the presentation of Achatz et al. (2017) where, expanding on previous work by Grimshaw (1975), the theory is discussed mostly in nondimensional form. We translate the essentials into dimensional form and choose, for easier tractability, a heuristic formulation. For all mathematical details, the reader is referred back to Achatz et al. (2017).

Starting point are the compressible Euler equations on an *f* plane, (e.g., Durran 1989), with Coriolis parameter *f*, without external sources or sinks,

where **u** and *w* are the horizontal and vertical components of the total wind **v**, respectively; $cp$ and $cV=cp\u2212R$ are the specific heat capacities at constant volume and pressure, respectively, with *R* the ideal gas constant of dry air; *θ* is potential temperature; *π* is the Exner pressure; and *g* is the gravitational acceleration.

Within this setting we consider a superposition of an exclusively altitude-dependent hydrostatic reference atmosphere at rest, a rather general synoptic-scale flow, and a locally monochromatic small-scale wave field:

Here $\theta \xaf\u2061(z)$ and $\pi \xaf\u2061(z)$ are potential temperature and Exner pressure of a horizontally symmetric and steady reference atmosphere. All other fields depend on all three spatial coordinates and on time. Using a $\alpha =0,1$ we discriminate between a weakly stratified atmosphere, as in the troposphere, indicated by $\alpha =1$, and a moderately strongly stratified case $(\alpha =0)$, characterizing the stratosphere, so that

to be understood so that $\theta \xaf0$ is a constant in the second case, and $\theta \xaf1/\theta \xaf0=O\u2061(\epsilon )$, where $\epsilon \u226a1$ is the Rossby number. Capital letters indicate the leading-order and next-order [with a tilde, so that $|V\u02dc|/|V|=O\u2061(\epsilon )$, for example] synoptic-scale flow, and hatted quantities the leading-order amplitudes of a locally monochromatic wave field, with local phase *ϕ*, yielding local wavenumber $k=\u2207\varphi $ and frequency $\omega =\u2212\u2202\varphi /\u2202t$. All additional contributions are next-order terms, in an expansion of all fields in terms of *ε*, and higher harmonics of the wave, induced by nonlinear interactions.

The synoptic-scale flow is assumed to vary on horizontal and vertical scales $Ls$ and $Hs$, respectively, and on the time scale $Ts$. Following standard considerations one can show that for a rotating atmospheric layer with inertia *f* and characteristic temperature $T00$ these synoptic scales are represented by

The various dynamical fields scale as

The wave field is assumed to vary on spatial and temporal scales $Lw$, $Hw$, and $Tw$ shorter than those of the synoptic-scale flow,

and its dynamical fields are those of a large-amplitude wave close to the margin of static instability, leading to

With these scaling assumptions one obtains the following main results: (i) Frequency and wavenumber are connected by the dispersion relations of either geostrophic modes, not to be pursued here any further, or IGWs,

where *k*, *l*, and *m* are the wavenumber components in *x*, *y*, and *z* directions, respectively, and $N2=\u2061(g/\theta \xaf)d\theta \xaf/dz$ is the squared Brunt–Väisälä frequency. Hence, with

and its wavenumber gradient $cg=\u2207k\Omega $ is the local group velocity, the prognostic equation for the local wavenumber is given by

(ii) The IGW amplitudes satisfy the polarization relations, with $kh=kex+\u0142ey$ the horizontal wavenumber and $b^=g\theta ^/\theta \xaf0$ the wave buoyancy amplitude,

(iii) Higher harmonics are negligible to leading order. (iv) The leading-order amplitudes satisfy wave-action conservation:

with $A=Ew/\omega ^$ the wave action, where

is wave energy, and $\rho \xaf$ the leading-order reference-atmosphere density.

Of central interest here are the following results. (v) To leading order the synoptic-scale flow is horizontal,

in geostrophic equilibrium,

with $\Psi =cp\theta \xaf0\Pi /f$ the streamfunction given by the leading-order synoptic-scale pressure fluctuations, and the leading-order synoptic-scale flow is hydrostatic,

so that the thermal wind relation

holds. Here the first term on the right-hand side is the classic result for the weakly stratified troposphere, while the second term contributes under more strongly stratified conditions (e.g., in the stratosphere). (vi) The wave impact on the synoptic-scale flow is exerted by wave-entropy flux-convergence in the entropy equation

where the asterisk indicates complex conjugation, and in the momentum equation

by anelastic wave-momentum flux-convergence $\u2212\u2061(1/\rho \xaf)\u2207\u22c5\u2061(1/2)R\u2061(\rho \xafv^u^*)$ and an elastic term $(f/g)ez\xd7\u2061(1/2)R\u2061(u^b^*)$. As usual, only the next-order ageostrophic fields appear in the Coriolis and in the pressure-gradient term while the leading-order geostrophic fields cancel each other. (vii) If the leading-order synoptic-scale flow is in geostrophic and hydrostatic equilibrium, it is described completely by its streamfunction. A prognostic equation for the latter can be obtained from the entropy equation, (25); the horizontal momentum equation, (26); and the leading-order Exner-pressure equation

using geostrophy, (22); the hydrostatic equilibrium, (23); the leading-order horizontal character of the flow, (21); and the IGW polarization relations. The result is a generalized quasigeostrophic potential vorticity (PV) equation

for quasigeostrophic PV

Here $c^g=cg\u2212U$ is the intrinsic group velocity. One sees that quasigeostrophic PV is driven by the vertical curl of the convergence of IGW pseudomomentum flux $c^gkhA$, and this determines completely the time development of the leading-order synoptic-scale flow.

Of considerable importance is the following observation, following a similar argument given in a simplified Boussinesq context by Bühler (2009): *if* (i) in the entropy equation, (25), the entropy fluxes are neglected and (ii) the anelastic momentum-flux convergence and the elastic term in the horizontal-momentum equation, (26), are replaced by pseudomomentum-flux convergence, that is,

then one obtains the same PV equation, (28), *provided the synoptic-scale flow is horizontal, geostrophic, and hydrostatic*. Hence the dynamics will obey (28) and the synoptic-scale flow will develop in a realistic manner. The corresponding calculations are nearly the same, only a bit simpler, as those leading from the original equations shown in (25) and (26)–(28), described by Achatz et al. (2017), and are therefore not detailed here any further. Therefore a recipe for an IGW parameterization can be to implement it exclusively into the horizontal-momentum equation, via pseudomomentum-flux convergence (Andrews and McIntyre 1976). In a single-column approximation, where all horizontal flux gradients are neglected, this is what present-day operational IGW parameterizations do (e.g., Alexander and Dunkerton 1999; Warner and McIntyre 2001; Scinocca 2002, 2003; Orr et al. 2010), and this pseudomomentum approach guarantees that a geostrophically and hydrostatically balanced flow is affected correctly by parameterized IGWs. However, the basis for this approach is even more fundamental. In their seminal work Andrews and McIntyre (1978) show that in a Lagrangian-mean reference frame the effect of arbitrary-amplitude GWs on a large-scale flow indeed only occurs in the momentum equation, by an appropriate forcing. Nonetheless an issue remains how to transfer this result to the Eulerian-mean reference frame that atmospheric models use, both for weather and climate applications. Andrews and McIntyre (1976), in the derivation of their generalized Eliassen–Palm-flux convergence for the zonal-mean zonal wind equation, use a scaling where weak wave amplitudes are assumed, and where hence approximate geostrophy and hydrostaticity of the mean flow allow the derivation of an approximate thermal-wind balance that can be used to avoid GW entropy-flux convergence and its application in a mean-flow entropy equation. Because of the weak-amplitude assumption advection of zonal momentum by the residual mean flow turns out to be negligible as well, so that one is left with a direct acceleration of the zonal-mean zonal wind by a the convergence of pseudomomentum (or generalized Eliassen–Palm) flux. The reliance of the pseudomomentum-flux result of Achatz et al. (2017) on large-scale balance appears to be directly related to this, however without the need for weak wave amplitudes. The following question may be asked: what happens if the required balance conditions are not fulfilled anymore? Although the theory of Achatz et al. (2017) exploits synoptic scaling, and hence derives geostrophic and hydrostatic balance, one can see (25) and (26) as intermediate results not explicitly relying on the large-scale flow being balanced, while the derivation of (28) makes heavy use of large-scale flow balance. It might therefore be more realistic to use in a more direct approach the unmodified entropy-flux convergence, anelastic momentum-flux convergence and elastic term. This could be the case, for example, if large-scale flow at different scales are of interest, or also ageostrophic flow such as the residual circulation. A central goal of the present work is a comparison of the two approaches in validations against wave-resolving simulations, and it will be shown that the direct approach is considerably more reliable.

In a first comparison between the pseudomomentum and the direct approach, one can easily point out that the pseudomomentum approach does not include any elastic or heating effect by IGWs while, mathematically speaking, each momentum-flux component

with *i* being *x* or *y*, and *j* being *x*, *y*, or *z*, can be related to its corresponding pseudomomentum-flux component

by

with

One sees that the absolute value of the momentum flux is always greater than or equal to the absolute value of the pseudomomentum flux. The difference becomes notable in the low-frequency limit where $\omega ^=O\u2061(f)$ if

that is, the ratio between horizontal and vertical wavelength must be of the order or larger than $N/f$. As the whole goes with the square of this ratio, corresponding factors get rapidly quite large. For example, for typical values of $N/f=102$, wavelength ratios of 10^{2} or 2 × 10^{2} lead to a factor 2 or 5, respectively. Also of interest is that in some cases certain contributions to the pseudomomentum flux vanish, that is, $PMFxx=0$ if $k=0$ and $PMFyy=0$ if $l=0$, while the corresponding momentum flux is nonzero. As we will show below this can have consequences for the structure of the mean-flow response. Finally we note that all $\gamma ij=1$ in the nonrotational case. Also the elastic terms and the GW entropy fluxes vanish then. Hence, the difference between the approaches exists only outside the tropics, and it is negligible for midfrequency and high-frequency GWs as well.

## 3. Description of the numerical models

In this section, the numerical code used for the validation cases is described, used either in a wave-resolving high-resolution mode with a subgrid-scale (SGS) turbulence parameterization to provide reference data or in a low-resolution mode with a WKB module switched on for validations of the two parameterization approaches.

### a. PincFloit LES

The Pseudoincompressible Flow solver with implicit turbulence modeling (PincFloit) solves on an *f* plane the pseudoincompressible equations of Durran (1989) in flux form (Klein 2009; Rieper et al. 2013b):

where *ρ* and $\rho \xaf$ are total density and reference-atmosphere density, respectively, and $\theta \u2032=\theta \u2212\theta \xaf$ the deviation of potential temperature from its reference-atmosphere value. The Σ is the viscous stress tensor with elements

where *η* is the dynamic shear viscosity coefficient, $\nu t$ is the turbulent viscosity (see below) and Pr is the turbulent Prandtl number, here taken to be 1/2, while *Q* is the total heating, here zero in all wave-resolving simulations. The elliptic problem for the pressure is solved using the hypre package (Falgout et al. 2006).

In the original numerical implementation of Rieper et al. (2013b) the code uses an implicit scheme for the parameterization of SGS turbulence. More recently Remmler et al. (2015) have shown in benchmark validations against direct numerical simulations of GW breaking that although this approach works quite well, it is surpassed in performance in some cases by the dynamic Smagorinsky method (Germano et al. 1991). The present implementation therefore calculates the resolved fluxes by an upstream-centered scheme for conservation laws [Monotonic Upwind Scheme for Conservation Laws (MUSCL)] following van Leer et al. (1991) with a modern total variation diminishing (TVD) limiter (Toro 1999; Kemm 2010), and parameterizes turbulence effects by the dynamic Smagorinsky approach. Following corresponding rules (Germano et al. 1991; Remmler et al. 2015) a turbulent Smagorinsky coefficient $CS2$ is determined, then averaged for stability reasons over a local five-point smoothing window along all the available spatial directions, and finally a turbulent kinematic viscosity is determined via

where

and

The turbulent viscosity is then used in the dynamic shear viscosity coefficient via

where $\eta m$ is the molecular dynamic shear viscosity.

### b. PincFloit/MS-GWaM

The numerical implementation of the transient interaction processes between IGW and large-scale flow is achieved by coupling the Multi-Scale Gravity-Wave Model (MS-GWaM), a Lagrangian phase-space WKB ray tracer, to the PincFloit code in a coarse-resolution setting. MS-GWaM is based on the methods described by Muraschko et al. (2015), Bölöni et al. (2016), and Wilhelm et al. (2018). It predicts the spectral distribution of wave action by a phase-space wave-action density $N$ so that wave-action density and energy density are

Its development is predicted by

with $k\u02d9$ given by (15). In the Lagrangian approach wavenumber and phase-space wave-action density are followed along rays parallel to the local **c**_{g} and $k\u02d9$. These respond to changes in the resolved horizontal flow **u** and stratification *N*^{2}, where from here on we denote the total resolved flow by lowercase letters. If coupled to MS-GWaM, PincFloit is used without molecular and turbulent viscosities, since IGWs should lead to the leading SGS effects.

In the *direct approach* the IGW impact on the resolved flow is mediated by anelastic momentum-flux convergence, the elastic term, and by entropy-flux convergence. Since in pseudoincompressible dynamics the nondiffusive continuity equation, that is, (41) with $\nu t=0$, is equivalent, together with the divergence constraint, (42), and the equation of state, (43), to the entropy equation

Here the subscript GW indicates IGW fields and the angle brackets stand for a local average over typical IGW wavelengths. Following the considerations above, the momentum fluxes are obtained from

From the polarization relations one also obtains

In the *pseudomomentum approach* only the momentum equation is forced by the IGWs, via pseudomomentum-flux convergence leading to the replacements

At every Runge–Kutta step, information is exchanged between IGWs and large-scale flow dynamics. MS-GWaM determines the wave-induced fluxes based on either direct or pseudomomentum scheme, and updates the resolved large-scale flow in PincFloit. After integrating the PincFloit model, the new resolved large-scale flow information is given to MS-GWaM, which solves the ray-tracing equations and yields updated parameterized wave-induced fluxes, thus closing the circle. Within the last step we also use a saturation scheme for mimicking turbulent wave breakdown. Whenever the local wave amplitudes are large enough to make the wave field statically unstable, a turbulent viscosity is switched on that reduces the wave amplitudes within one time step to the saturation level. More details on the methods used by MS-GWaM can be found in Muraschko et al. (2015), Bölöni et al. (2016), and Wilhelm et al. (2018).

## 4. Wave packet validation cases

In this section, we introduce the various wave packet cases used for the validation of the two approaches at stake. In all of these an isothermal (*T* = 240 K) atmosphere at rest is initially superposed to a locally monochromatic Gaussian IGW packet. Therein,

where

and $V\u02dc$ can be obtained from $B\u02dc$ using the polarization relations (16) and (17); $m<0$ is the initial vertical wavenumber and the horizontal wavenumber components are $(k,l)=kh\u2061(sin\alpha ,\u2009cos\alpha )$, with $kh>0$ and $0\u2264\alpha \u2264\pi /2$. The factor $a0$ indicates the wave packet amplitude with regard to the static-instability threshold, reached when $a0=1$, where the isentropes in the wave packet begin to overturn. The positive dispersion-relation branch is taken, so that $m<0$ corresponds to a local group velocity pointing upward. The wave amplitudes are modulated in the horizontal in *y* direction to be reminiscent of a zonally symmetric amplitude distribution inducing a zonal wind, as is the paradigmatic picture in parameterized IGW–mean flow interactions. Four cases are considered: 1) A horizontally symmetric IGW field with $1/\sigma y=0$ and $1/\sigma x=0$ (1DWP), 2) and 3) IGW fields modulated in *y* direction, with $\sigma y<\u221e$ and $1/\sigma x=0$, where the angle *α* between the horizontal wavenumber and the horizontal direction of modulation of the wave packet is either 0 (2DWP00) or *π*/2 (2DWP90). In both of these cases the vertical group velocity is so small that a simulation of a wave packet from initialization to the point where it becomes statically unstable would be computationally very demanding. Since the IGW group velocity is about linearly increasing with the Coriolis parameter, we have therefore also decided to follow in a 4) further case (2DWP90-HAMP-HF) the approach of Lelong and Dunkerton (1998a,b) of considering an atmosphere with a Coriolis parameter an order of magnitude larger. To keep the same ratio $(kh/m)N/f$ as above, the horizontal wavelength in that case is one order of magnitude smaller. The launch amplitude has there also been chosen stronger than in all other cases, in order to further reduce the required computational time. The details of the settings are given in Table 1. Particularly, in each wave-resolving PincFloit LES, there are at least 16 grid points per wavelength of the initial IGW along each available horizontal direction, and 10 grid points per vertical wavelength of the initial IGW. In each PincFloit/MS-GWaM simulation with IGW parameterization, the resolution along each available direction is close to the wavelength of the initial parameterized IGW along the corresponding direction, and the large-scale flow induced by the parameterized IGWs can still be resolved in the coarse-resolution PincFloit/MS-GWaM simulation. The resolution in the wave-resolving simulations, using PincFloit LES is too coarse to really justify these as LES. It is, however, sufficient for our present purposes, as we are not addressing turbulent wave breaking itself.

## 5. Incipient budget analysis

In the following we first discuss the various cases in terms of the initial contributions,

to the momentum-flux convergence or pseudomomentum-flux convergence, respectively, in terms of the elastic terms,

and in terms of the wave-induced heating,

In the case of 1DWP all horizontal derivatives in the various flux convergences vanish. Hence, in the pseudomomentum scheme the vertical gradient of the pseudomomentum flux in the *x*-momentum equation (referred to as PMFC_{xz}) is the only nonzero forcing term. In the direct scheme, the vertical gradient of the momentum flux and the elastic term in the *x*-momentum equation (referred to as MFC_{xz} and ET_{x}, respectively) are the only two nonzero forcing terms. The comparison of the abovementioned three forcing terms (i.e., PMFC_{xz} in the pseudomomentum scheme and MFC_{xz} and ET_{x} in the direct scheme) is given in Fig. 1a. In the case of 1DWP, the signal of MFC_{xz} is apparently stronger than that of PMFC_{xz}. Compared with MFC_{xz}, ET_{x} is relatively small but noticeable.

For further illustration, Fig. 2 shows the dependence of the abovementioned three forcing terms on zonal wavelength (top panel) and vertical wave packet scale $\sigma z$ (bottom panel). In each sensitivity test, all the other parameters are kept the same as those in 1DWP. For waves with relatively short horizontal wavelengths (hence intrinsic frequencies well above the inertia frequency), the signal of PMFC_{xz} may be still close to that of MFC_{xz}, and ET_{x} is negligible. As the horizontal wavelength increases, the signal of PMFC_{xz} becomes weaker, while MFC_{xz} is getting stronger instead. A large separation is seen between PMFC_{xz} and MFC_{xz} for waves with relatively long horizontal wavelengths, and ET_{x} becomes more noticeable but still secondary compared with MFC_{xz}. In the sensitivity test to the changes of $\sigma z$ only, the pattern of PMFC_{xz} is the same as MFC_{xz}, but the signal of MFC_{xz} is approximately 3 times as strong as that of PMFC_{xz}. The signals of both positive branch (i.e., upper branch) and negative branch (i.e., lower branch) of MFC_{xz} are getting weaker as $\sigma z$ increases. However, the negative branch is apparently much more sensitive to the change of $\sigma z$ compared with the positive branch. This asymmetric behavior is caused by the vertically decreasing density. The center value of ET_{x} is not sensitive to the change of $\sigma z$, and $\sigma z$ only changes the vertical width of ET_{x}. Generally speaking, in the current study, ET_{x} appears to be secondary compared with MFC_{xz}. However, with larger $\sigma z$, ET_{x} gets more important. Hence, it is not only the IGW aspect ratio between vertical and horizontal wavelength that determines, as described above, the relative magnitude of the dynamical terms in the direct and pseudomomentum approach, but also the IGW amplitude distribution, as here exemplified using the wave packet width.

In the case of 2DWP00, the nonzero IGW-induced forcing terms in the direct approach only include MFC_{yy}, MFC_{yz}, and ET_{y} in the *y*-momentum equation, and there is no forcing in the *x*-momentum equation and thermodynamic equation. Together with the wave-energy density initial vertical cross sections of all nonzero terms are shown in Fig. 3. There is a dipole pattern at the altitude of wave packet center for MFC_{yy}, causing a horizontally asymmetric forcing effect. Instead, the dipole structure in MFC_{yz} results in a vertically asymmetric forcing effect. The ET_{y} has consistent negative signals with minimum at the wave packet center. It is worth mentioning that the forcing from MFC_{yy} is comparable with that from MFC_{yz}, and it is very important in this case. Similar to 1DWP, ET_{y} appears to be secondary but not negligible in this case. To facilitate a comparison between pseudomomentum and direct approach, the corresponding $\gamma ij$ are shown in the upper-left corner of the subplots of the contributors to the momentum-flux convergence, indicating how much stronger the respective momentum flux is than the corresponding pseudomomentum flux. Both $\gamma yy$ and $\gamma yz$ are as large as 3.2, which indicates much weaker signals of both PMFC_{yy} and PMFC_{yz} in the pseudomomentum approach, compared to MFC_{yy} and MFC_{yz} in the direct approach. This again suggests that the forcing from MFC_{yy}, which is neglected in operational IGW schemes, is comparable with that from MFC_{yz} in this case.

Case 2DWP90 is interesting for several reasons. Figure 4 shows the initial distribution of all nonzero IGW-induced forcing terms in the direct approach. These include MFC_{yy} in the *y*-momentum equation, MFC_{xz} and ET_{x} in the *x*-momentum equation, and the IGW-induced heating term in the thermodynamic equation. In contrast to these, PMFC_{xz} is the only nonzero term in the pseudomomentum approach; that is, only *x* momentum is forced. Both MFC_{yy} and ET_{x} are due to rotational effects. Among these, only MFC_{yy} and the wave-induced heating term are horizontally asymmetric, while in the pseudomomentum approach no horizontal asymmetry is forced initially.

## 6. Numerical simulations of different idealized wave packet profiles prior to turbulent dissipation

In this section, the direct and pseudomomentum approach will be validated against wave-resolving simulations with unbalanced large-scale flows. The initial condition settings of the respective numerical experiments are given by the three idealized wave packet profiles described in the previous section (i.e., 1DWP, 2DWP00, and 2DWP90). In each wave-resolving PincFloit LES, there are 16 or ~17 grid points per wavelength of the initial IGW along each available horizontal direction, and 10 grid points per vertical wavelength of the initial IGW. As we have convinced ourselves from 2D simulations with increased resolutions (not shown), these resolutions are sufficient for the laminar phase of the occurring IGW–mean flow interactions, while the later turbulent IGW breaking would require a finer resolution. In each PincFloit/MS-GWaM simulation, the resolution along each available direction is close to the wavelength of the initial IGW in the corresponding direction, and the large-scale flow induced by the parameterized IGWs can still be resolved in the coarse-resolution PincFloit/MS-GWaM simulation. In the PincFloit LESs, the wave packets of the resolved IGWs are initialized in two-dimensional domains for 1DWP and 2DWP00, and three-dimensional domain settings are employed for 2DWP90. Instead, in the PincFloit/MS-GWaM simulations, 1DWP uses a single vertical column, and two-dimensional domain settings are employed for all two-dimensional wave packet cases. The details of the model designs are given in Tables 1 and 2.

Figure 5 compares the large-scale *x*-velocity component and the large-scale *y*-velocity component between the PincFloit LES, and the PincFloit/MS-GWaM simulations using either direct or pseudomomentum approach. Similar to the procedure in Bölöni et al. (2016) and Wilhelm et al. (2018), the large-scale wind component in PincFloit LES is obtained by a running mean over two wavelengths of the initial IGW. The selected time for this comparison in Fig. 5 is at *t* = 1500 min, which is approximately 1.5 inertial periods. Since the wave energy of the 1DWP case does not depend on the horizontal direction, the induced large-scale zonal velocity component and meridional velocity component do not change along the horizontal direction, and the large-scale vertical velocity component is zero everywhere in the domain. As depicted in Fig. 5, the direct scheme is able to reproduce the signals in the data from the wave-resolving simulation with very small errors. The pseudomomentum scheme, on the contrary, gives a much weaker large-scale horizontal wind compared with the wave-resolving data, and thus large model errors could be expected.

Figure 6 shows corresponding results, again at *t* = 1500 min, but now for case 2DWP00. For the two-dimensional wave packets in the current study, in addition to the induced large-scale horizontal wind on both directions, the large-scale vertical velocity component is nonzero. Again, the direct scheme successfully captures the signals in the data from the wave-resolving simulation. Even though the pattern from the results based on the pseudomomentum scheme is still similar to the wave-resolving simulation for all the three wind components, its signals are much weaker. Note that the fact that the pseudomomentum scheme generates much weaker large-scale wind components is consistent with the budget analysis shown in the previous section.

Figure 7 shows corresponding results for the 2DWP90 case. It again demonstrates that the direct scheme is more reliable than the pseudomomentum scheme, with the data from the wave-resolving simulation as a reference. Two further things are worth noting here. 1) Arguably the best agreement between pseudomomentum approach and wave-resolving simulation is visible in the *x*-wind component. This could be expected since this is the horizontal wind component that can also contain balanced contributions, given by the *y* derivative of the balanced part of the pressure. The other side of the coin is, however, that the unbalanced response in the meridional circulation is met well at best in structure but definitely not in strength. 2) All three large-scale wind components produced by the pseudomomentum scheme for this particular case are highly horizontally symmetrical, which is quite different from the results from the direct scheme and the wave-resolving data. This suggests that the horizontally asymmetrical forcing is apparently negligible in the pseudomomentum scheme in this case. As suggested in Fig. 4 from the previous section, this is again consistent with the budget analysis in the pseudomomentum scheme for the 2DWP90 case. Indeed here *l* = 0 and hence PMF_{yy} = 0, while MF_{yy} is nonzero and explains the asymmetry of the response.

## 7. Sensitivity experiments

Various sensitivity tests have been done. In the first of these we have convinced ourselves that our findings even hold for IGW amplitudes very close to overturning instability. To this purpose we have considered a high-amplitude two-dimensional wave packet case (2DWP90-HAMP). Here, the initial setting is the same as that in 2DWP90, except that the initial wave amplitude at the wave packet center reaches the limit of convective instability. Therefore, a stronger wave–mean flow interaction is expected. To be compared to Fig. 7 for 2DWP90, Fig. 8 shows the comparisons of the large-scale mean wind among three model codes at the same time for 2DWP90-HAMP (see also Table 3). Similar large-scale wind patterns are found in cases 2DWP90 (Fig. 7) and 2DWP90-HAMP (Fig. 8), except that the signals in 2DWP90-HAMP are approximately 4 times as strong as those in 2DWP90. With a stronger wave–mean flow interaction, this case continues to support the statement that the direct approach is generally more reliable than the pseudomomentum approach.

Next, expanding on the budget analysis in section 5, we have investigated the relative importance of the elastic effect and the heating term for the results from the direct approach. Figure 9 shows results from three sensitivity experiments, where either the elastic term has been set to zero, or the heating term, or both. These are to be compared to the results from the direct scheme for the 2DWP90-HAMP case in Fig. 8. These sensitivity experiments are denoted by 2DWP90-HAMP-NoET, 2DWP90-HAMP-NoHeat, and 2DWP90-HAMP-NoET-NoHeat, respectively. With PincFloit LES as a reference, all of the abovementioned sensitivity experiments demonstrate a much stronger negative large-scale zonal mean wind at the domain center, especially when the heating effect is switched off. This implies that it is important to take into account the heating term in order to get the correct large-scale flow response. The elastic term, however, is, at least in this study, noticeable but still of secondary importance.

In all of the cases considered so far, the vertical group velocity is so small that a simulation of a wave packet from initialization to the point where it becomes statically unstable by propagation in increasingly attenuated altitudes, hence growing in amplitude, would be computationally very demanding. Since the IGW group velocity is about linearly increasing with the Coriolis parameter we have therefore also decided to follow in a further case (2DWP90-HAMP-HF) the approach of Lelong and Dunkerton (1998a,b) of considering an atmosphere with a Coriolis parameter an order of magnitude larger. To keep the same ratio $(kh/m)N/f$ as above, the horizontal wavelength in that case is one order of magnitude smaller. The comparison among the three model codes for the 2DWP90-HAMP-HF case at a late stage of the simulation is shown in Fig. 10. Although the large-scale *x* wind, with its balanced components, is relatively reasonable in the simulations with the pseudomomentum scheme, the direct scheme is the only one able to reproduce the response of the *y*–*z* circulation in its full strength. It is also able to match the horizontal asymmetry of the response, while the pseudomomentum approach misses this feature.

Finally, Fig. 11 illustrates the importance of the wave–mean flow interaction, as well as the role of the saturation scheme, by comparing for the same case 2DWP90-HAMP-HF the large-scale *y*–*z* circulation and the wave energy among two sensitivity experiments and their corresponding control simulation, all using PincFloit/MS-GWaM with the direct approach, but at a much later stage, where we do not have PincFloit LES data available. Here, in the first sensitivity experiment the coupling process between the large-scale mean wind and the gravity waves has been turned off by assuming that the large-scale mean wind has no impact on changing the wavenumber of the gravity waves via the ray-tracing equation. In the second sensitivity experiment the saturation scheme has been turned off by assuming that the wave action density associated with each ray is always conserved during its vertical propagation. Differences are easily identified between each abovementioned sensitivity experiment and its corresponding control simulation. The findings of Bölöni et al. (2016) are confirmed that the direct wave–mean flow interaction dominates the dynamics to leading order, while turbulent wave breaking leads to next-order modifications.

## 8. Concluding remarks and discussion

This article investigates and compares two different approaches for the efficient modeling of subgrid-scale inertia–gravity waves (IGWs) in a rotating compressible atmosphere. The first approach, denoted as a pseudomomentum scheme, is closely related to the fundamental result of Andrews and McIntyre (1978) that in a Lagrangian-mean reference frame GW effects on the large-scale flow only occur in the momentum equation. It can be shown that this holds also in the Eulerian-mean reference frame used by atmospheric models, at least if the large-scale flow is in hydrostatic and geostrophic balance. The GW forcing can then be expressed by a pseudomomentum-flux convergence in the large-scale momentum equation. Present-day gravity wave parameterizations follow this route, exclusively applying a corresponding forcing, and leaving out elastic and thermal effects that may arise in the absence of large-scale-flow balance. In the second approach, called a direct scheme, the large-scale flow is forced both in the momentum equation by anelastic momentum-flux convergence and in the entropy equation via entropy-flux convergence. In addition, also an elastic term is taken into account in the momentum equation. This second approach does not rely on any balance assumption with regard to the large-scale flow, and it differs from the pseudomomentum approach wherever rotational effects matter, that is, outside of the tropics and for low-frequency IGWs.

It has been the purpose of the current study to understand the difference between the pseudomomentum and the direct approaches. To this end wave-resolving simulations have been done, to provide reference data against which to validate coarse-resolution model simulations with either of the two approaches providing the basis for a prognostic gravity wave parameterization. All simulations have been done with the code PincFloit, solving the pseudoincompressible equations by a finite-volume method. At high resolutions PincFloit has been used together with a dynamical scheme for the parameterization of the largest turbulent eddies (PincFloit LES). In the coarse-resolution simulations it has been coupled to the prognostic gravity wave model MS-GWaM, solving the coupled WKB equations using a Lagrangian approach. This model has been used in either pseudomomentum mode or direct mode.

The dynamics of 1D IGW packets, with an amplitude only depending on time and the vertical direction, already shows that the comparative performance of the two approaches is sensitive to the following two parameters: 1) the intrinsic frequency and 2) the wave packet scale. The smaller the intrinsic frequency is, the larger the discrepancy between pseudomomentum-flux convergence and momentum-flux convergence is, and, hence, the greater the differences between two schemes are. Pseudomomentum-flux convergence is then much weaker than momentum-flux convergence. The elastic term is also nonzero in that case, but it generally appears to be of comparatively smaller importance, and only begins exhibiting some effect if the wave packet scale is rather large.

Further 2D case studies, with IGW envelope depending on time and *y* and *z*, meant to be reminiscent of a zonally symmetric IGW distribution in interaction with a zonally symmetric large-scale flow, confirm the expectations from the 1D analyses. They show that only the direct approach is able to capture the response in the unbalanced components of the large-scale flow. This response is strongly underestimated by the pseudomomentum approach. It will be interesting to investigate, for example, what consequences this has for the correct description of the IGW impact on the residual circulation. Even in the balanced response, however, some structural asymmetries can only be captured by the direct approach, while the more symmetric pseudomomentum-flux convergence is not able to do so. In these 2D cases it is especially the horizontal entropy-flux convergence that seems to matter, while the elastic term does not seem to be that important. It should be stressed that these differences are only relevant for low-frequency IGWs, and outside the tropics, while for high-frequency GWs and in the tropics pseudomomentum and direct approach are approximately equivalent. Finally, in the present settings we reconfirm the findings of Bölöni et al. (2016) that the direct wave–mean flow interaction takes a leading role, as compared to turbulence effects.

In summary, our results show that the direct approach is more reliable than its pseudomomentum counterpart. In some cases, the pseudomomentum approach yields a balanced response that is of the right strength, however, not quite with the right structure. Even more conspicuous, however, appears the difference in the reproduction of the potentially most unbalanced part of the response. The meridional-circulation response from the pseudomomentum scheme is generally too weak. As the direct scheme is not much more expensive than the pseudomomentum scheme, it seems advisable to use it instead for GW parameterizations in atmospheric models, at least outside the tropics. Nonetheless, our analyses have not yet considered global applications where irreversible GW effects due to wave breaking will weigh in more strongly. For these the comparison between the two approaches, beyond the scope of the present study, is still to be done. Our results give considerable motivation for this next step.

## Acknowledgments

We benefited from the insightful comments from two anonymous reviewers on an earlier version of the manuscript. The authors also 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-2, AC 71/9-2, AC 71/10-1, AC 71/10-2, AC 71/11-2, AC 71/12-2, BO 5071/2-2, and BO 5071/1-2. Calculations for this research were conducted on the LOEWE-CSC high-performance computer of the Goethe Universität Frankfurt.

## REFERENCES

*Waves and Mean Flows*, Cambridge University Press, 182–183.

*Numerical Solution of Partial Differential Equations on Parallel Computers*, Springer, 267–294.

*Atmosphere–Ocean Dynamics*. 1st ed. Academic Press, 662 pp.

*An Introduction to Atmospheric Gravity Waves*. 1st ed. Academic Press, 279 pp.

*Geophys. Res. Let.*,

**32**, L18715, https://doi.org/10.1029/2005GL023226.

*Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction*. 2nd ed. Springer, 724 pp.

*10th Computational Fluid Dynamics Conf.*, Honolulu, HI, American Institute of Aeronautics and Astronautics, 260–282, https://doi.org/10.2514/6.1991-1552.

## Footnotes

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

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).