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.

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.

References 4–7 have done fundamental studies of the interaction between GWs and mean flow that have more recently been extended by Refs. 8–10 to consider a wider range of atmospheric stratification, higher GW harmonics, and to give a deepened account on the conditions for GW impacts on the mean flow. Whatever the considered GW amplitudes, we will consider GWs in interaction with a synoptic-scale flow in a hydrostatic reference atmosphere, with profiles θ̄(z),ρ̄(z),π̄(z) of potential temperature, density, and Exner pressure, respectively, that only depend on altitude 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/θ̄)dzθ̄, with g the gravitational acceleration, is f/N=O[ε(5α)/2], where α = 0, 1 denotes moderately strong or weak stratification, and ε=O(1/10) is a small parameter. Observations11,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 T00 a typical atmospheric temperature, representative horizontal and vertical length scales and a representative time scale of such GWs are
(1)
whence the wind scales are obtained in the large-amplitude case as advective wind scales (Uw, Ww) = (Lw, Hw)/Tw. 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 (Ls, Hs, Ts) = (Lw, Hw, Tw)/ɛ, i.e., ɛ is our scale-separation parameter. Note that the synoptic-scale wind scales are (Us, Ws) = (Uw, Ww) so that Us/(Lsf) = ɛ, i.e., ɛ also is the Rossby number of the synoptic-scale flow.9,10
The next step is to use the wave wind, length, and time scales, and T00, to non-dimensionalize the equations of motion
(2)
(3)
(4)
(5)
leading to
(6)
(7)
(8)
(9)
where Dt = t + v · ∇ indicates the material derivative, u and w are the horizontal and vertical components of the total wind v, respectively. cp and cV = cpR are the specific heat capacities at constant pressure and volume, respectively, with R the ideal gas constant of dry air. f0 = 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,
(10)
(11)
(12)
where Θ̄(j) and Π̄(j) are due to the reference atmosphere, all terms proportional to the phase factors exp 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 βk=βϕ/ε=βXϕ=exX+eyY+ezZϕ 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).
Inserting (10)(12) into (6)(9) and sorting by terms with equal powers in ɛ and in the phase factor exp /ɛ 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.,
(13)
(14)
(15)
with B0(0)=Θ0(0)/Θ̄(0) a non-dimensional synoptic-scale buoyancy, and ∇X,h = exX + eyY the horizontal gradient operator in the slow spatial variables.
The leading-order results for the wave field are that frequency and wave number satisfy a dispersion relation ω = Ω(k, X, T) for either geostrophic modes (GM) or gravity waves (GW), i.e.,
(16)
where the dependence on 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Θ̄(α)/Θ̄(0). The components of the wavenumber vector are defined so that k = kex + ley + mez, and kh=k2+l2 is the absolute magnitude of the horizontal wavenumber kh = kex + ley. The dispersion relations entail the eikonal equations
(17)
(18)
for the development of wave number and frequency, with cg=kΩ=exk+eyl+ezmΩ the local group velocity. The wave amplitudes satisfy the polarization relations
(19)
(20)
(21)
with ω̂=ωkU0(0) the intrinsic frequency. One also finds that for GWs Bβ(0)=0 if β > 1, and only for j > 0 one can have Bβ(j)0 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.
A central result from the next-order terms in ɛ is the wave-action conservation equation
(22)
for the GW wave-action density A=Egw/ω̂, with
(23)
the GW energy density, where R̄(0)(Z) is the leading-order reference-atmosphere density, decaying strongly from the ground to higher altitudes. The derivation of this result uses the geostrophic and hydrostatic equilibrium (14) and (15) of the synoptic-scale flow. For the leading-order contributions to the synoptic-scale flow one obtains the prognostic equations
(24)
(25)
(26)
where P̄(0)=R̄(0)Θ̄(0) is the leading-order mass-weighted potential temperature of the reference atmosphere. While there is no direct GW impact on the synoptic-scale Exner pressure, GW entropy-flux and momentum-flux convergence appear in prognostic equations for mean-flow potential temperature and horizontal momentum. With ĉg=kω̂=cgU0(0) being the intrinsic group velocity, the entropy fluxes are
(27)
and the momentum-flux tensor
(28)
has the elements
(29)
(30)
(31)
(32)
(33)
The last GW term in the horizontal momentum equation (26) is the so-called elastic term that appears only at moderately strong stratification, i.e., for α = 0. It also vanishes in the absence of rotation.
Note that the results above represent a fully nonlinear theory where the nonlinear advection terms turn out to not contribute to the leading orders of the wave equations because the solenoidality property
(34)
of the wave velocity field, following directly from the leading-order wave part of the Exner-pressure equation, eliminates self-advection of the wave field to leading order.
The eikonal equations (17) and (18) can lead to a breakdown of the local monochromaticity assumed in the derivation of the wave-action equation (22), e.g., if an initial locally monochromatic GW field has a spatial dependence of wave number so that neighboring regions have group velocities leading to caustics where rays cross. This calls for a theory describing the dynamics of the superposition of several wave fields. However, because of the nonlinear advection terms, the superposition of wave fields with strong amplitudes close to breaking, as assumed above, does not even allow the derivation of dispersion and polarization relations anymore. In the locally monochromatic case the advection terms do not contribute to the leading orders of the wave equations, because of the solenoidality property (34). In the case of a superposition of several wave fields, however, mutual advection of the different wave fields only does not contribute to the leading order of the equations if the wave amplitudes are sufficiently weak, i.e., when the wave fields in (10)(12) are multiplied by a factor ɛn with, e.g., n = 1, termed weakly nonlinear case or n = 2, the quasilinear case. One sets
(35)
(36)
(37)
where the index β 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.,
(38)
(39)
With the replacements β(k,ω̂)(kβ,ω̂β) one also obtains the polarization relations (19)(21).
In the quasilinear case n = 2 one obtains from the next-order terms in ɛ the wave-action conservation equation
(40)
for the GW wave-action density Aβ=Egw,β/ω̂β, with
(41)
the energy density of the GW field indicated by β. 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
(42)
one can show from the eikonal equations and the wave-action equations that it satisfies the conservation equation
(43)
where we have introduced the more compact notation k̇=XΩ. This equation indicates that the integral of phase-space wave-action density over the total phase-space volume is conserved. Moreover, cg = ∇kΩ and the definition of k̇=XΩ imply that the six-dimensional phase-space velocity is non-divergent,
(44)
Hence the conservation equation (43) can also be written
(45)
In contrast to wave-action density A=βAβ=d3kN in position space, the spectral wave-action density is conserved along rays in phase space satisfying dT(x,k)=(cg,k̇). We also note that the superposition (42) allows for an arbitrary number of components with arbitrary wavenumbers each, so that effectively spectral wave-action density can also be a truly continuous function of wavenumber.

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.

Both in the quasilinear and in the weakly nonlinear case the prognostic equations for the mean flow are still (24)(26), but now the contributing fluxes are due to the full spectrum, i.e.,
(46)
and the momentum-flux tensor has the elements
(47)
(48)
(49)
(50)
(51)
Note that these fluxes only contribute to leading order in the large-amplitude case. It is known, however, that their effects accumulate over longer times so that they cannot be ignored. The cleanest way to handle this would be to introduce a correspondingly slow time scale. For simplicity we avoid this step and just keep the wave impacts on the mean flow as outlined above. In the following we will use (46)(51) also for the monochromatic large-amplitude case, where it is understood that n = 0 and that N(X,k,T)=A(X,T)δkk(X,T), with k as an argument of N not depending on X and T anymore.
Because of the geostrophic and hydrostatic equilibrium (14) and (15) the synoptic-scale mean flow satisfies an extended quasigeostrophic theory. Using these equilibrium conditions, one can derive from the prognostic equations (24)(26) a single prognostic equation
(52)
for the quasigeostrophic potential vorticity (QGPV)
(53)
with G=d3kĉgkN and H=d3kĉglN the fluxes of the zonal and meridional components of pseudomomentum ph=d3kkhN. Inverting QGPV yields the leading-order synoptic-scale Exner-pressure fluctuations Π0(0) whence one can obtain, using geostrophic and hydrostatic equilibrium, the leading-order horizontal wind and potential-temperature fluctuations of the synoptic-scale flow. QGPV is forced by the vertical curl of the pseudomomentum-flux convergences, and in the absence of waves it is conserved.
The practitioner needs the results above in their dimensional form. These are as follows: A re-dimensionalization of the GW dispersion relation in (70), by the substitutions
(54)
(55)
(56)
(57)
(58)
(59)
leads to the dimensional GW dispersion relation
(60)
with N2=(g/θ̄)dzθ̄, and the trivial reformulation
(61)
for the weakly nonlinear and quasilinear case with a superposition of spectral components. Substituting
(62)
one obtains from (19)(21) the dimensional polarization relations
(63)
(64)
(65)
where bβ=gθβ/θ̄ is the dimensional buoyancy of the βth GW component. The corresponding results for the locally monochromatic large-amplitude case are obtained from (63)(65) by dropping the β-index.
Likewise, we obtain for the locally monochromatic case and for the quasilinear spectral case, respectively, the dimensional GW wave-action equations
(66)
where (cg, cg,β) = (∇kω, ∇kωβ) is the GW group velocity, and (A,Aβ)=(Ew/ω̂,Ew,β/ω̂β) the GW wave-action density, with
(67)
the wave-energy density. ρ̄ is the reference-atmosphere density.
The wave-action equations together with the dimensional forms
(68)
(69)
of the eikonal equations (17) and (38), where
(70)
leads to the spectral wave-action equations
(71)
(72)
for the spectral wave-action density
(73)
in the quasilinear case and N(x,k,t)=A(x,t)δkk(x,t) in the locally monochromatic large-amplitude case. Here as well we have defined k̇=xΩ. In deriving (72) from (71) one again exploits the non-divergence
(74)
of the phase-space velocity.
The GW impact on the synoptic-scale mean flow, indicated by angle brackets ⟨⋯⟩, is captured within synoptic scaling by supplementing the entropy equation by GW entropy-flux convergence,
(75)
and the horizontal-momentum equation by GW momentum-flux convergence and the elastic term
(76)
where the entropy flux can be obtained from the spectral wave-action density via
(77)
and the mass-specific momentum-flux tensor has the elements
(78)
(79)
(80)
(81)
(82)
In Eqs. (75) and (76) the mean flow is understood to be its full expansion
(83)
(84)
(85)
so that (75) is in the two leading orders in ɛ consistent with (13) and (25), and that (76) is in the two leading orders consistent with (15) and (26).
The fluxes (77)(82) are to be coded into weather-forecast and climate models, but the leading-order synoptic-scale mean-flow dynamics can also be expressed by the QG potential-vorticity equation
(86)
with G=d3kĉgkN and H=d3kĉglN the fluxes of the zonal and meridional components of GW pseudomomentum ph=d3kkhN,
(87)
and ψ=cpθ̄0δπ/f the streamfunction, where δπ=ε1+αΠ0(0) is the leading-order synoptic-scale Exner-pressure fluctuations, and θ̄0=T00Θ̄(0) is the leading-order reference-atmosphere potential temperature. The latter depends on 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,
(88)
and the leading-order synoptic-scale potential temperature fluctuations δθ=ε1+αT00Θ0(0), via hydrostatic equilibrium,
(89)

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.

We begin with energy. As can be shown in the derivation of the wave-action equation, GW energy density
(90)
satisfies
(91)
with
(92)
the wave-energy flux. The last term describes the exchange with the synoptic-scale mean flow. The latter has an energy density
(93)
where the first part is the kinetic energy density and the second part is the density of available potential energy. As can be derived from the QGPV equation (86), it obeys
(94)
Hence the prognostic equation for the total energy density Es + Ew,
(95)
contains only flux terms, so that under suitable boundary conditions the volume-integrated total energy is conserved.
For the derivation of a potential-vorticity conservation property one needs a prognostic equation for pseudomomentum ph=d3kkhN. As can be derived with the help of the wave-action-density equation (71), this prognostic equation is
(96)
The vertical curl of this yields
(97)
Comparing with the QGPV equation (86) one obtains the conservation equation
(98)
for an extension Π of quasigeostrophic potential vorticity that contains contributions from the synoptic-scale flow, that are linear in the synoptic-scale streamfunction, and the negative of the gravity-wave pseudovorticity ez×ph/ρ̄, that is nonlinear in the gravity-wave amplitudes. This conservation property can only be broken by non-conservative effects, e.g., GW sources or dissipative GW breaking.
Of considerable consequence for the properties of GW parameterizations is the direct consequence
(99)
of (98). Hence the leading-order synoptic-scale mean flow is not influenced by GWs if GW pseudovorticity is a Lagrangian invariant of the flow. The latter is the case, e.g., if
  • 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.

A numerical implementation of spectral wave-action dynamics could either be done in a Eulerian finite-volume formulation, based on (71), or in a Lagrangian approach starting from (72). In the latter16–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
(100)
Hence, following each point of a ray volume along the ray passing through it one obtains its propagation as a time-dependent volume within which wave-action density N is conserved. Because of the non-divergence (74) of the phase-space (ray) velocity the volume content of this volume is conserved as well. In the numerical discretization one assumes that each ray volume keeps a rectangular shape, but it is allowed to be stretched and squeezed in a volume-preserving manner. Because (74) holds separately in the two-dimensional spaces spanned by x and k, y and l, or z and m, i.e.,
(101)
(102)
(103)
this is done so that the areas Δ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.
FIG. 1.

Illustration of the Lagrangian discretization of the propagation of a phase-space volume, here in the sub-space spanned only by m and z. It is subdivided into small rectangular ray volumes. Each ray volume propagates with a mean velocity averaged from the phase-space velocities at the center of the upper and lower edge of the ray volume. Different group velocities cgz at the edges lead to a stretching or squeezing of the ray-volume extent Δz in z-direction. The extent Δm in m direction is then adjusted so that the area ΔzΔm is conserved. The procedure in the xk and yl sub-spaces is the same. Reproduced with permission from Muraschko et al., Q. J. R. Metereol. Soc. 141, 676 (2015). Copyright 2015 Wiley Publishers.

FIG. 1.

Illustration of the Lagrangian discretization of the propagation of a phase-space volume, here in the sub-space spanned only by m and z. It is subdivided into small rectangular ray volumes. Each ray volume propagates with a mean velocity averaged from the phase-space velocities at the center of the upper and lower edge of the ray volume. Different group velocities cgz at the edges lead to a stretching or squeezing of the ray-volume extent Δz in z-direction. The extent Δm in m direction is then adjusted so that the area ΔzΔm is conserved. The procedure in the xk and yl sub-spaces is the same. Reproduced with permission from Muraschko et al., Q. J. R. Metereol. Soc. 141, 676 (2015). Copyright 2015 Wiley Publishers.

Close modal

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 cg and to the wavenumber velocity k̇ 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.

The approach described above on the simulation of the interaction between subgrid-scale GWs and a resolved flow, in climate models but also in weather-forecast codes, differs in various regards from the approach presently applied by the weather and climate centers: First, we point out that present-day GW parameterizations use the observation that the GW impact on the resolved-flow QGPV in (86) could also be obtained by dropping the GW entropy-flux convergence from the mean-flow entropy equation (75), by removing the elastic term from the mean-flow momentum equation (76), and by replacing the momentum flux there by the pseudo-momentum flux, i.e., by using
(104)
and
(105)
and then starting from there the asymptotic analysis of synoptic-scale flow dynamics. Therefore the common approach is to only force the resolved-flow momentum equation, by GW pseudomomentum-flux convergence instead of GW momentum-flux convergence, and to not have the GW parameterization acting on the thermodynamics. It is computationally simpler and also has a certain conceptual attractiveness because all of the GW impact can be attributed to the pseudomomentum-flux divergence. From the theory, this could be a viable approach, but it assumes that the resolved flow is geostrophically and hydrostatically balanced. Given the general trend to increasingly finer model grids, with also larger-scale unbalanced GWs more and more being part of the resolved flow, this is however an issue that has been investigated by Ref. 22. Indeed they show, in comparisons between idealized GW resolving simulations and coarse-grid simulations with a GW parameterization either using the general (direct) approach or the classic (pseudomomentum) approach, that the former is much more reliable in the simulation of the GW-mean-flow interaction.
The second and third aspect relate to simplifications in the implementation that are mostly due to considerations of computational efficiency: Because code parallelization typically works in vertical columns, i.e., different processors are allotted different horizontal locations, it is computationally less expensive to use parameterizations for subgrid-scale processes that do not couple different horizontal positions. In the context here this amounts to ignoring the impact of horizontal mean-flow gradients on the horizontal wavenumbers, i.e., assuming k̇=l̇=0, and to ignoring all horizontal group-velocity components, and hence to replacing the spectral wave-action equations (71) and (72) by
(106)
(107)
This neglects horizontal GW propagation, and it also neglects the response of the horizontal wave number to horizontal variations of the resolved flow. Moreover, in this so-called single-column approximation one also neglects horizontal GW fluxes in the GW impact on the large-scale flow, i.e., one approximates (105) by
(108)
Finally, in a further attempt to gain in efficiency, one neglects, in the so-called steady-state approximation, the time dependence of the wave-action density, i.e., one uses instead of (106)
(109)
or, integrating in vertical wavenumber
(110)
Given a lower boundary condition, obtained from some description of the GW sources the latter equation is integrated vertically to obtain a profile of wave-action density that agrees with what one would get, from a steady lower boundary condition, after a period sufficiently long for rays to propagate from the lower boundary to the model top. Effectively this is then an equilibrium profile in agreement with instantaneous GW propagation from the lower boundary to the model top. In practice, one always has a discrete sum of spectral components that are emitted from some parameterized source. Hence one recurs to the representation (73) and solves the single-column and steady-state simplification
(111)
of (66) for each spectral component.

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.

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 models26 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 theories14 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.

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.

The authors have no conflicts to disclose.

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 sharing is not applicable to this article as no new data were created or analyzed in this study.

1.
D.
Fritts
and
M.
Alexander
,
Rev. Geophys.
41
,
1003
, (
2003
).
2.
Y.-J.
Kim
,
S. D.
Eckermann
, and
H.-Y.
Chun
,
Atmos.-Ocean
41
,
65
(
2003
).
3.
M. J.
Alexander
,
M.
Geller
,
C.
McLandress
,
S.
Polavarapu
,
P.
Preusse
,
F.
Sassi
,
K.
Sato
,
S.
Eckermann
,
M.
Ern
,
A.
Hertzog
,
Y.
Kawatani
,
M.
Pulido
,
T. A.
Shaw
,
M.
Sigmond
,
R.
Vincent
, and
S.
Watanabe
,
Q. J. R. Metereol. Soc.
136
,
1103
(
2010
).
4.
F.
Bretherton
,
Q. J. R. Metereol. Soc.
92
,
466
(
1966
).
6.
D.
Andrews
and
M.
McIntyre
,
J. Fluid Mech.
89
,
609
(
1978
).
7.
D.
Andrews
and
M.
McIntyre
,
J. Fluid Mech.
89
,
647
(
1978
).
8.
U.
Achatz
,
R.
Klein
, and
F.
Senf
,
J. Fluid Mech.
663
,
120
(
2010
).
9.
U.
Achatz
,
B.
Ribstein
,
F.
Senf
, and
R.
Klein
,
Q. J. R. Metereol. Soc.
143
,
342
(
2017
).
10.
U.
Achatz
,
Atmospheric Dynamics
(
Springer
,
Berlin, Heidelberg
,
2022
).
11.
12.
J.
Callies
,
R.
Ferrari
, and
O.
Bühler
,
Proc. Natl. Acad. Sci. U. S. A.
111
(
48
),
17033
(
2014
).
13.
J.
Pedlosky
,
Geophysical Fluid Dynamics
(
Springer
,
1987
).
14.
15.
C.
Eden
,
M.
Chouksey
, and
D.
Olbers
,
J. Phys. Oceanogr.
49
,
291
(
2019
).
16.
G. S.
Völker
,
T. R.
Akylas
, and
U.
Achatz
,
Q. J. R. Meteorol. Soc.
147
,
1112
(
2021
).
17.
J.
Muraschko
,
M.
Fruman
,
U.
Achatz
,
S.
Hickel
, and
Y.
Toledo
,
Q. J. R. Metereol. Soc.
141
,
676
(
2015
).
18.
G.
Bölöni
,
B.
Ribstein
,
J.
Muraschko
,
C.
Sgoff
,
J.
Wei
, and
U.
Achatz
,
J. Atmos. Sci.
73
,
4833
(
2016
).
19.
Y.-H.
Kim
,
G.
Bölöni
,
S.
Borchert
,
H.-Y.
Chun
, and
U.
Achatz
,
J. Atmos. Sci.
78
,
1339
(
2021
).
20.
R. S.
Lindzen
,
J. Geophys. Res.
86
,
9707
, (
1981
).
21.
G.
Bölöni
,
Y.-H.
Kim
,
S.
Borchert
, and
U.
Achatz
,
J. Atmos. Sci.
78
,
1317
(
2021
).
22.
J.
Wei
,
G.
Bölöni
, and
U.
Achatz
,
J. Atmos. Sci.
76
,
2715
(
2019
).
23.
M. P.
Baldwin
,
L. J.
Gray
,
T. J.
Dunkerton
,
K.
Hamilton
,
P. H.
Haynes
,
W. J.
Randel
,
J. R.
Holton
,
M. J.
Alexander
,
I.
Hirota
,
T.
Horinouchi
,
D. B. A.
Jones
,
J. S.
Kinnersley
,
C.
Marquardt
,
K.
Sato
, and
M.
Takahashi
,
Rev. Geophys.
39
,
179
, (
2001
).
24.
C. C.
Stephan
,
H.
Schmidt
,
C.
Zülicke
, and
V.
Matthias
,
J. Geophys. Res.: Atmos.
125
,
e2019JD031528
, (
2020
).
25.
C.
Stephan
,
J.
Duras
,
L.
Harris
,
D.
Klocke
,
W.
Putman
,
M.
Taylor
,
N.
Wedi
,
N.
Zagar
, and
F.
Ziemen
,
Tellus A
74
,
280
299
(
2022
).
26.
J.
Slingo
,
P.
Bates
,
P.
Bauer
,
S.
Belcher
,
T.
Palmer
,
G.
Stephens
,
B.
Stevens
,
T.
Stocker
, and
G.
Teutsch
,
Nat. Clim. Change
12
,
499
(
2022
).
27.
I. M.
Held
,
Bull. Am. Meteorol. Soc.
86
,
1609
(
2005
).
28.
Y.-H.
Kim
and
U.
Achatz
,
Geophys. Res. Lett.
48
,
e2021GL095226
, (
2021
).
29.
J. K.
Hunter
and
J. B.
Keller
,
Wave Motion
6
,
79
(
1984
).
30.
L.
Gosse
and
N. J.
Mauser
,
J. Comput. Phys.
211
,
326
(
2006
).
31.
Y.
Manor
,
J. Phys. A: Math. Gen.
10
,
765
(
1977
).
32.
G. V.
Pereverzev
,
Nucl. Fusion
32
,
1091
(
1992
).