The interaction between small-scale waves and a larger-scale flow can be described by a multi-scale theory that forms the basis for a new class of parameterizations of subgrid-scale gravity waves (GW) in weather and climate models. The development of this theory is reviewed here. It applies to all interesting regimes of atmospheric stratification, i.e., also to moderately strong stratification as occurring in the middle atmosphere, and thereby extends classic assumption for the derivation of quasi-geostrophic theory. At strong wave amplitudes a fully nonlinear theory arises that is complemented by a quasilinear theory for weak GW amplitude. The latter allows the extension to a spectral description that forms the basis of numerical implementations that avoid instabilities due to caustics, e.g., from GW reflection. Conservation properties are discussed, for energy and potential vorticity, as well as conditions under which a GW impact on the larger-scale flow is possible. The numerical implementation of the theory for GW parameterizations in atmospheric models is described, and the consequences of the approach are discussed, as compared to classic GW parameterizations. Although more costly than the latter, it exhibits significantly enhanced realism, while being considerably more efficient than an approach where all relevant GWs are to be resolved. The reported theory and its implementation might be of interest also for the efficient and conceptually insightful description of other wave-mean interactions, including those where the formation of caustics presents a special challenge.

## I. INTRODUCTION

With horizontal and vertical wavelengths down to at most 1 km and 100 m respectively, mesoscale atmospheric waves such as internal gravity waves (GWs) will not all be simulated explicitly by operational climate models within the foreseeable future. However, without taking their influence into account, climate models miss essential circulation aspects even on the planetary scale.^{1–3} Hence they must be parameterized, requiring a solid theory for the interaction between the waves and a mean flow. Corresponding studies have led in recent years to the emergence of a new class of GW parameterizations (GWP). The present review is to give an overview on these developments, from theoretical investigations to the implementation of a new GWP into a state-of-the-art climate model. For this purpose we first review the dynamics of large-amplitude, locally monochromatic GWs, then discuss spectra of weak-amplitude GWs, next describe the GW impact in general on the mean flow, touch on conservation properties, and finally give a sketch of the hence resulting numerical developments.

## II. LARGE-AMPLITUDE LOCALLY MONOCHROMATIC WAVES

*z*. We also assume an

*f*-plane with a constant Coriolis parameter. One can assume that the ratio between Coriolis frequency

*f*and Brunt–Vaisala frequency $N=(g/\theta \u0304)dz\theta \u0304$, with

*g*the gravitational acceleration, is $f/N=O[\epsilon (5\u2212\alpha )/2]$, where

*α*= 0, 1 denotes moderately strong or weak stratification, and $\epsilon =O(1/10)$ is a small parameter. Observations

^{11,12}indicate that, in the spectrum of GWs, most energy is carried by the waves that are in scale just below the synoptic scale, which we assume to be resolved by climate models of interest. A careful analysis then turns out that, with

*T*

_{00}a typical atmospheric temperature, representative horizontal and vertical length scales and a representative time scale of such GWs are

*U*

_{w},

*W*

_{w}) = (

*L*

_{w},

*H*

_{w})/

*T*

_{w}. Waves of such amplitude are close to locally inducing static instability, i.e., negative vertical derivatives of total potential temperature, which would cause them to break. The corresponding mean-flow synoptic scales are (

*L*

_{s},

*H*

_{s},

*T*

_{s}) = (

*L*

_{w},

*H*

_{w},

*T*

_{w})/

*ɛ*, i.e.,

*ɛ*is our scale-separation parameter. Note that the synoptic-scale wind scales are (

*U*

_{s},

*W*

_{s}) = (

*U*

_{w},

*W*

_{w}) so that

*U*

_{s}/(

*L*

_{s}

*f*) =

*ɛ*, i.e.,

*ɛ*also is the Rossby number of the synoptic-scale flow.

^{9,10}

*T*

_{00}, to non-dimensionalize the equations of motion

*D*

_{t}=

*∂*

_{t}+

**v**· ∇ indicates the material derivative,

**u**and

*w*are the horizontal and vertical components of the total wind

**v**, respectively.

*c*

_{p}and

*c*

_{V}=

*c*

_{p}−

*R*are the specific heat capacities at constant pressure and volume, respectively, with

*R*the ideal gas constant of dry air.

*f*

_{0}= 1 is a non-dimensional placeholder for the Coriolis frequency. We then introduce slow variables (

**X**,

*T*) =

*ɛ*(

**x**,

*t*), and insert into the non-dimensional equations the WKB expansions for a superposition of a hydrostatic reference atmosphere, a synoptic-scale flow, and a wave field with its higher harmonics,

*iβϕ*/

*ɛ*are contributions from the wave (subscript

*β*= 1 for the basic wave, and

*β*≥ 2 for its

*β*th higher harmonic), and the rest constitutes the synoptic-scale part (subscript 0). Both the wave amplitudes and the synoptic-scale flow are only slowly varying in space and time, as are the local wavenumbers $\beta k=\beta \u2207\varphi /\epsilon =\beta \u2207X\varphi =ex\u2202X+ey\u2202Y+ez\u2202Z\varphi $ and frequencies

*βω*= −

*β∂*

_{t}

*ϕ*/

*ɛ*= −

*β∂*

_{T}

*ϕ*. The leading orders of the expansions follow from a dimensional analysis of the basic equations. For the synoptic-scale flow they agree in the weak-stratification case

*α*= 1 with standard quasigeostrophic scaling (e.g., Ref. 13).

*ɛ*and in the phase factor exp

*iϕ*/

*ɛ*one finds from the leading orders in

*ɛ*that the mean flow is to leading order horizontal, and that it is in hydrostatic and geostrophic balance, i.e.,

_{X,h}=

**e**

_{x}

*∂*

_{X}+

**e**

_{y}

*∂*

_{Y}the horizontal gradient operator in the slow spatial variables.

*ω*= Ω(

**k**,

**X**,

*T*) for either geostrophic modes (GM) or gravity waves (GW), i.e.,

**X**and

*T*enters via that of the mean-flow leading-order horizontal wind $U0(0)$ and the non-dimensional Brunt–Vaisala frequency $N02=dZ\Theta \u0304(\alpha )/\Theta \u0304(0)$. The components of the wavenumber vector are defined so that

**k**=

*k*

**e**

_{x}+

*l*

**e**

_{y}+

*m*

**e**

_{z}, and $kh=k2+l2$ is the absolute magnitude of the horizontal wavenumber

**k**

_{h}=

*k*

**e**

_{x}+

*l*

**e**

_{y}. The dispersion relations entail the eikonal equations

*β*> 1, and only for

*j*> 0 one can have $B\beta (j)\u22600$ if

*β*> 1. Hence GW higher harmonics do not contribute to the leading order, but potentially only to the next but leading order. This is not the case for the GM higher harmonics, but in the following we ignore this option and only consider wave fields without GM contribution.

*ɛ*is the wave-action conservation equation

*α*= 0. It also vanishes in the absence of rotation.

## III. A SPECTRUM OF WEAK-AMPLITUDE WAVES

*ɛ*

^{n}with, e.g.,

*n*= 1, termed

*weakly nonlinear*case or

*n*= 2, the

*quasilinear*case. One sets

*β*indicates different wave fields. Those could also include higher harmonics of another wave field, but due to the weak wave amplitudes these do not contribute to all relevant orders. Take first a locally monochromatic wave field, where the wave amplitudes are non-zero only for a single

*β*. Inserting (35)–(37) into the equations of motion (6)–(9), sorting by terms with equal powers in

*ɛ*, and separating wave part and mean-flow part by averaging over scales sufficiently longer than the wave scales one retrieves the mean-flow results (13)–(15). From the leading-order wave contributions one again obtains the dispersion relations

*ω*

_{β}= Ω(

**k**

_{β},

**X**,

*T*) for either geostrophic modes or gravity waves, entailing also the eikonal equations, i.e.,

*quasilinear case*

*n*= 2 one obtains from the next-order terms in

*ɛ*the wave-action conservation equation

*β*. One then observes that leading and next-order wave terms giving the results above are strictly linear in the wave fields. Hence the superposition (35)–(37) is a solution as well, with each component satisfying its own set of eikonal equations (38) and (39) and wave-action equation (40). Defining for this superposition the spectral wave-action density

**c**

_{g}= ∇

_{k}Ω and the definition of $k\u0307=\u2212\u2207X\Omega $ imply that the six-dimensional phase-space velocity is non-divergent,

In the *weakly nonlinear case* *n* = 1 the wave-action equation is to be supplemented by the effect of wave-wave interactions, of either GWs with GWs or GWs with GMs. To the best of our knowledge, the corresponding scattering integrals have only been worked out fully for Boussinesq dynamics without mean flow,^{14,15} and first steps for Boussinesq dynamics with non-vanishing mean flows have been taken by Ref. 16. In the atmospheric context they have so far been simply ignored.

*n*= 0 and that $N(X,k,T)=A(X,T)\delta k\u2212k(X,T)$, with

**k**as an argument of $N$ not depending on

**X**and

*T*anymore.

## IV. WAVE IMPACT ON THE BALANCED MEAN FLOW

## V. SUMMARY OF THE RESULTS IN DIMENSIONAL FORM

*β*th GW component. The corresponding results for the locally monochromatic large-amplitude case are obtained from (63)–(65) by dropping the

*β*-index.

**c**

_{g},

**c**

_{g,β}) = (∇

_{k}

*ω*, ∇

_{k}

*ω*

_{β}) is the GW group velocity, and $(A,A\beta )=(Ew/\omega \u0302,Ew,\beta /\omega \u0302\beta )$ the GW wave-action density, with

*ɛ*consistent with (13) and (25), and that (76) is in the two leading orders consistent with (15) and (26).

*z*only in the case with moderately strong stratification (

*α*= 0), while it is a constant in the weakly stratified case (

*α*= 1). The streamfunction also yields the leading-order synoptic-scale horizontal wind via geostrophic equilibrium,

## VI. CONSERVATION PROPERTIES

Neither GW energy is conserved nor is QGPV. The interaction between GWs and synoptic-scale mean flow leads to an exchange between the two so that the corresponding conserved quantity comprises contributions from both components.

### A. Energy

*E*

_{s}+

*E*

_{w},

### B. Potential vorticity

### C. Non-acceleration

GW amplitudes and wavenumbers are steady,

GW pseudovorticity does not vary horizontally, and

the GWs are not affected by sources or sinks.

Classic GW parameterizations assume steady-state GW fields and they do not take horizontal variations of the GWs fields into account. Hence they rely exclusively on GW sources and GW breaking as processes leading to a GW impact on the resolved flow.

## VII. CONSEQUENCES FOR GW PARAMETERIZATIONS

### A. Numerical implementation

^{16–19}one sub-divides the part of phase space with non-zero wave-action density into rectangular so-called ray volumes. Wave-action density is conserved along trajectories (called rays) satisfying

*x*and

*k*,

*y*and

*l*, or

*z*and

*m*, i.e.,

*x*Δ

*k*, Δ

*y*Δ

*l*and Δ

*z*Δ

*m*are conserved, with Δ

*x*, Δ

*y*, Δ

*z*, Δ

*k*, Δ

*l*, and Δ

*m*the ray volume extent in the six respective phase-space directions. Figure 1 illustrates this for the space spanned by

*z*and

*m*. Volume deformations away from an original rectangular shape are to be represented by splitting a corresponding volume into a sufficiently large number of rectangular ray volumes.

GW dissipation is simulated by a classic saturation approach based on Ref. 20 and adapted by Refs. 18 and 21 to spectral wave-action dynamics. A static-instability breaking threshold, with a tuning factor close to 1, is determined by the criterion that, within a resolved-flow volume cell, the constructive interference of all ray volumes can lead to a total GW signal with a negative vertical derivative in potential temperature larger than the positive derivative given by the stratification of the resolved flow. Once this threshold is exceeded turbulence is invoked, causing turbulent viscosity and diffusivity that provide a dissipative right-hand side to the spectral wave-action equation (72) so that the GW field is kept at the threshold of static instability.

The interaction between the parameterized GWs and the resolved large-scale flow takes two directions. The contribution of the large-scale winds ⟨**u**⟩ to the group velocity **c**_{g} and to the wavenumber velocity $k\u0307$ influences the development of the wave-action density. In the other direction, GWs influence the resolved flow in its thermodynamics via the divergence of the GW entropy flux (77). They also change its momentum via the divergence of GW momentum flux (78)–(82), and via the elastic term, i.e., the last term in (76). That term can be obtained from the GW entropy-flux convergence as well. The wavenumber integrals in the fluxes are presently estimated to first-order accuracy, by evaluating the integrand at the center of the ray volume and multiplying it by its wavenumber volume. So far implementations have been into finite-volume dynamical model cores for resolved-flow dynamics. In this context the fluxes are then projected onto the finite-volume cell faces and finally used there.

It would seem attractive to avoid the projection of the fluxes from the Lagrangian GW model onto the resolved-flow finite volume cells, and also the interpolation of the resolved-flow winds to the location of each ray volume, by a straightforward finite-volume implementation of the spectral wave-action equation (71) in flux form. This alternative has been studied by Ref. 18. It turns out that the Lagrangian approach is computationally much more efficient. Two factors contribute to this: Firstly, in a finite-volume approach, one must span a six-dimensional phase space, with often a substantial fraction of cells not contributing essentially to the GW fluxes. Secondly, GW refraction and reflection are only captured well in the finite-volume approach if the resolution in wavenumber space is excessively high.

### B. Comparison with classic GWP

As follows from the non-acceleration result discussed in Subsection VI C the classic approach using the single-column and steady-state approximations, the first neglecting GW horizontal propagation and horizontal GW fluxes, and the second neglecting GW transience, relies exclusively on GW dissipation as a process allowing a GW impact on the resolved flow. The relevance of GW transience has been investigated by Refs. 18, 19, and 21. As shown by Ref. 18, the resolved-flow impact induced by the breaking of a single GW packet is often described quite incorrectly if GW transience is not taken into account. The consequence of this for GW parameterization in a global climate model has been studied by Refs. 19 and 21. The result is that the statistics of simulated GW fluxes, exhibiting in measurements distributions with long tails of strong fluxes, are affected significantly by a steady-state assumption, to the point that the distributions are distorted significantly away from the observational findings.

The consequences of the single-column approximation have been investigated in the mean time as well (Voelker *et al.* and Kim *et al.*, both to be submitted elsewhere). As with regard to the neglect of wave transience, they appear to be of leading order as well. The horizontal distribution of GW fluxes in the middle atmosphere (above 15 km altitude) is very different between the single-column and the general parameterization approach so that also the large-scale winds differ significantly. In the tropics this leads to a significantly changed variability, with the quasi-biennial oscillation (e.g., Ref. 23) differing in its period conspicuously between simulations using the two approaches. Hence it appears that the general approach is considerably more realistic.

## VIII. FINAL DISCUSSION

The multi-scale theory outlined above has led to the development of an extended approach, in weather and climate models, for the parameterization of subgrid-scale GWs. As is clear by now it is considerably more realistic than classic GW parameterizations, where the GW impact has been simplified by focusing on GW pseudomomentum, where GW transience is neglected, and where the effects of horizontal GW propagation and of horizontal GW fluxes are neglected. This comes at a computational prize: Simulations with the new approach are slower, by about an order of magnitude at global-model resolutions tested so far, than those using a classic GW parameterization. Yet, they are more realistic, and if one wanted to resolve the GWs that turn out to matter, simulations would be and are (e.g., Refs. 24 and 25) by many orders of magnitude slower than even this (Ref. 21). Hence the approach discussed here seems to be a reasonable compromise between the realism of GW permitting simulations, certainly of their own value as a data source for many studies, and the efficiency of climate simulations using classic GW parameterizations. Moreover, in view of the increasing complexity of ever higher resolved weather and climate models^{26} there is an ever-increasing need for a hierarchy of models that help us gain conceptual insight into the complex atmosphere.^{27} In this hierarchy the more general approach for GW parameterizations should have its place.

Open issues remain, both on the theoretical and on the applied side. One aspect that seems to deserve consideration is the interaction between GWs and the geostrophic vortical mode that is the other constituting component of mesoscale atmospheric dynamics (e.g., Refs. 9 and 12). Present approaches for the description of this interaction, using tools of wave-turbulence theory,^{15} do not take the presence of a leading-order mean flow into account. In the ocean context this is possible, but in the atmosphere mean winds enter to leading order. This also holds for another aspect of relevance, the so-far neglect of GW–GW interactions. It is not clear how relevant this process will be in the end, and to the best of our knowledge this has not been investigated yet. Here as well available theories^{14} suffer from the neglect of a mean flow. It would be interesting to obtain a corresponding extension of the theory, e.g., following the route indicated by Ref. 16. Such work would close a conceptual gap that we still seem to have in the weakly nonlinear regime between the quasilinear regime (allowing a spectral approach) and the large-amplitude regime (where only locally monochromatic GW fields are possible). Finally, the handling of GW breaking by the saturation approach is very crude and a theory encompassing GW-turbulence interaction still needs to be derived from the basic equations. Similar considerations hold for the emission of GWs by various processes, be it flow over orography, convection, or emission by jets and fronts, to just name the most often discussed candidates. In all of these instances closed mathematical descriptions are still to be derived from multi-scale approaches, that hopefully would be able to replace present schemes. Finally, we think that in the numerical implementation of the general approach the last word has not been spoken. It would be helpful if experts in numerical mathematics gave it a closer look, especially with regard to accuracy and efficiency.

Finally we also want to point out that it might be of interest to consider and extend the tools and techniques outlined here for other challenging problems of wave-mean-flow-interaction theory. The quasi-biennial oscillation of the equatorial zonal-mean zonal winds,^{23} e.g., is not only due to GWs but also to larger-scale tropical Kelvin and Rossby–Gravity waves. Moreover, Kelvin waves also interact with GWs in the tropics.^{28} The interaction between mesoscale GWs, larger-scale tropical waves and the zonal-mean flow could be an interesting problem to be studied using these methods. Another field of application could be caustics in general, e.g., in nonlinear acoustics,^{29} quantum mechanics,^{30} general relativity,^{31} or plasma physics,^{32} where the spectral approach outlined here could help the development of closed and comparatively simple treatments.

## ACKNOWLEDGMENTS

U.A. thanks the German Research Foundation (DFG) for partial support through the research unit “Multiscale Dynamics of Gravity Waves” (MS-GWaves, Grant Nos. AC 71/8-2, AC 71/9-2, and AC 71/12-2) and CRC 301 “TPChange” (Project No. 428312742, Projects B06 “Impact of small-scale dynamics on UTLS transport and mixing” and B07 “Impact of cirrus clouds on tropopause structure”). Y.-H.K. and U.A. thank the German Federal Ministry of Education and Research (BMBF) for partial support through the program Role of the Middle Atmosphere in Climate (ROMIC II: QUBICC) and through Grant No. 01LG1905B. U.A. and G.S.V. thank the German Research Foundation (DFG) for partial support through the CRC 181 “Energy transfers in Atmosphere an Ocean” (Project No. 274762653, Projects W01 “Gravity-wave parameterization for the atmosphere” and S02 “Improved Parameterizations and Numerics in Climate Models.”) U.A. is furthermore grateful for support by Eric and Wendy Schmidt through the Schmidt Futures VESRI “DataWave” project.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**U. Achatz**: Conceptualization (equal); Funding acquisition (equal); Writing – original draft (equal); Writing – review & editing (equal). **Y.-H. Kim**: Writing – review & editing (equal). **G. S. Voelker**: Writing – review & editing (equal).

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

## REFERENCES

*Atmospheric Dynamics*