## Abstract

As present weather forecast codes and increasingly many atmospheric climate models resolve at least part of the mesoscale flow, and hence also internal gravity waves (GWs), it is natural to ask whether even in such configurations subgrid-scale GWs might impact the resolved flow and how their effect could be taken into account. This motivates a theoretical and numerical investigation of the interactions between unresolved submesoscale and resolved mesoscale GWs, using Boussinesq dynamics for simplicity. By scaling arguments, first a subset of submesoscale GWs that can indeed influence the dynamics of mesoscale GWs is identified. Therein, hydrostatic GWs with wavelengths corresponding to the largest unresolved scales of present-day limited-area weather forecast models are an interesting example. A large-amplitude WKB theory, allowing for a mesoscale unbalanced flow, is then formulated, based on multiscale asymptotic analysis utilizing a proper scale-separation parameter. Purely vertical propagation of submesoscale GWs is found to be most important, implying inter alia that the resolved flow is only affected by the vertical flux convergence of submesoscale horizontal momentum at leading order. In turn, submesoscale GWs are refracted by mesoscale vertical wind shear while conserving their wave-action density. An efficient numerical implementation of the theory uses a phase-space ray tracer, thus handling the frequent appearance of caustics. The WKB approach and its numerical implementation are validated successfully against submesoscale-resolving simulations of the resonant radiation of mesoscale inertia GWs by a horizontally as well as vertically confined submesoscale GW packet.

## 1. Introduction

Internal gravity waves (GWs) play a significant role in atmospheric dynamics on various spatial scales (Fritts and Alexander 2003; Kim et al. 2003; Alexander et al. 2010; Plougonven and Zhang 2014). Already in the lower atmosphere GW effects are manifold. Examples include the triggering of high-impact weather (e.g., Zhang et al. 2001, 2003) and clear-air turbulence (Koch et al. 2005), as well as the effect of small-scale GWs of orographic origin on the predicted larger-scale flow (e.g., Palmer et al. 1986; Lott and Miller 1997; Scinocca and McFarlane 2000) and the GW impact on the generation of high cirrus clouds and polar stratospheric clouds (e.g., Joos et al. 2009). Even more conspicuous than in the lower atmosphere, however, are GW effects in the middle atmosphere. The general circulation in the mesosphere is basically controlled by GWs (Lindzen 1981; Holton 1982; Garcia and Solomon 1985). This also seems to be of relevance to both medium-range weather forecasts and climate modeling in the troposphere. Middle-atmosphere circulation influences the lower layers by downward control (Haynes et al. 1991), and there is evidence of the importance of the middle atmosphere for long-range forecasting of winter weather (Baldwin and Dunkerton 2001; Kidston et al. 2015; Hansen et al. 2017; Jia et al. 2017) and climate (Scaife et al. 2005, 2012) in the Northern Hemisphere.

As a substantial portion of the GW spectrum involves scales too small to describe explicitly in current-resolution climate models, accounting for such small-scale GWs poses an important parameterization problem to atmospheric dynamics. With rising computing power available, an increasing number of studies of middle-atmosphere global GW dynamics uses models that can resolve a part of the GW spectrum (Kawatani et al. 2009, 2010a,b; Brune and Becker 2013). This raises the question of whether the neglected subgrid-scale (SGS) GWs could impact the resolved flow, and, if so, how their effect could be taken into account. In regard to these issues, to the best of our knowledge, global weather forecast codes, with horizontal mesh distances of *O*(10) km, still generally use a parameterization of SGS GWs, while high-resolution local-area codes used by the weather services, with mesh distances of *O*(1) km, typically do not. Since the “effective resolution” in such codes is well above their mesh distances (Skamarock 2004; Ricard et al. 2013), one might suppose that even there a considerable portion of the GW spectrum is not captured. However, a systematic investigation of the potential impact of SGS GWs on the resolved mesoscale flow is lacking at present.

Available GW parameterizations (e.g., Lindzen 1981; Palmer et al. 1986; McFarlane 1987; Alexander and Dunkerton 1999; Warner and McIntyre 2001; Scinocca 2003; Orr et al. 2010) invariably rely on WKB theory (Bretherton 1966) for describing the interaction between scale-separated waves and (resolved) mean flow. However, the specific implications of this theory may depend on the scales involved. The classic scenario is the interaction between a resolved synoptic-scale flow and unresolved mesoscale inertia GWs. The corresponding WKB theory (Grimshaw 1975; Achatz et al. 2017) as well as the generalized Lagrangian-mean theory (Andrews and McIntyre 1978a,b; Bühler 2009) show that the wave amplitude is controlled by wave-action conservation, while the synoptic-scale flow is described by a quasigeostrophic potential vorticity that is affected by the GWs via pseudomomentum-flux convergence. For efficiency reasons, parameterizations use these theoretical results with drastic simplifications: (i) lateral GW propagation and the impact of horizontal mean flow gradients are ignored, and (ii) the time-dependent transient wave–mean flow interaction is replaced by an equilibrium picture where, because of the nonacceleration paradigm, GWs can only modify the resolved flow when they break. In the present context, the latter steady-state approximation especially may not be entirely justified. As pointed out by Bühler and McIntyre (1998, 2003, 2005), wave transience is potentially important, and recently Bölöni et al. (2016) have shown that in many cases it can attain at least an equally important role as turbulent wave breaking in mediating the impact of GWs on the resolved flow. It is also essential to keep in mind that the standard WKB approach assumes from the outset geostrophic and hydrostatic balance of the synoptic-scale flow. It is therefore not obvious that this theory can be applied to the interaction between a mesoscale resolved flow and mesoscale or submesoscale SGS GWs, which seems to be the most appropriate scenario for GW parameterizations in mesoscale-resolving models. In this setting, a modified WKB theory that allows for a mesoscale unbalanced large-scale flow would be most useful.

Of related interest is that packets of small-scale GWs are capable of radiating larger-scale GWs. This possibility was first suggested by Bretherton (1969) for two-dimensional small-scale GW packets with isotropic scaling, and more recently has been investigated further by Van den Bremer and Sutherland (2014) for wave packets of various aspect ratios. The radiation of large-scale waves hinges on a resonance mechanism, wherein the vertical phase velocity of the emitted long waves matches the vertical group velocity of the small-scale wave packet, which acts as a traveling wave source. Furthermore, the vertical wavenumber of the long wave is set by the scale of the wave packet envelope. In a related study, Tabaei and Akylas (2007) show that the longwave radiation process is especially enhanced if the small-scale wave packet is “flat” (i.e., its envelope is elongated in the horizontal relative to the vertical) so that both the horizontal and vertical envelope scales can be compatible with free, nearly steady, long GWs. All of these studies assume the small-scale GWs to be nonhydrostatic. They do not investigate, however, which small-scale GWs in general are able to interact with given mesoscale long GWs.

Moreover, no prior study examines the feasibility of a model for SGS GWs in a resolved mesoscale flow. Closest to this comes the WKB model of Tabaei and Akylas (2007). These authors, however, report numerical instability problems once the wave–mean flow interaction develops caustics where the initially locally monochromatic small-scale GW field exhibits multivalued wavenumbers. This problem, also observed by Rieper et al. (2013a), can be circumvented, however. As shown by Muraschko et al. (2015) and Bölöni et al. (2016), a spectral approach based on phase-space wave-action density yields numerically stable and fast algorithms for the efficient integration of the coupled equations of small-scale GWs in a larger-scale flow.

Building on the above brief review of related prior literature, the goals of the present paper are (i) a systematic investigation of which smaller-scale GWs are able to interact resonantly with given typical mesoscale GWs, (ii) the development of a WKB theory for the efficient description of this interaction, (iii) the implementation of a numerical algorithm for this theory, and finally (iv) the validation of the WKB theory and its numerical implementation against submesoscale-resolving simulations of the radiation of mesoscale GWs by horizontally and vertically confined submesoscale GW packets as considered earlier by Bretherton (1969), Van den Bremer and Sutherland (2014), and Tabaei and Akylas (2007).

The paper is structured as follows. The GW scales of interest are identified in sections 2a and 2b; these scales form the basis for nondimensionalizing the Boussinesq equations and formulating an appropriate multiscale asymptotic ansatz in section 2c. Leading- and next-order WKB approximations are discussed in section 2d, which eventually yield (section 2e) a coupled energy-conserving equation system of linear Boussinesq equations for the mesoscale, and one-dimensional ray equations for the submesoscale dynamics. Subsequently, our numerical models are described in section 3; the initial conditions of the numerical experiments are motivated in section 4a, and section 4b briefly discusses the postprocessing of the model output data. In section 4c, a kinematic analysis similar to ship wake theory is used to predict the geometry of the induced mesoscale wave disturbance, on the assumption of steady state forcing by a propagating submesoscale GW packet. In section 4d, the simulation results of different test cases are presented and compared against the theoretical predictions. Finally, the article concludes with a summary and discussion of the main findings in section 5.

## 2. Theory: Basics and formalism

For simplicity the interaction between mesoscale and submesoscale GWs is studied in a rotating, incompressible, and inviscid Boussinesq atmosphere with height-dependent background stratification, characterized by Coriolis parameter *f* and Brunt–Väisälä frequency *N*(*z*). Under these flow conditions, the governing equations are

where *D*/*Dt* = ∂/∂*t* + **v** ⋅ ∇ is the material derivative, **e**_{z} denotes the unit vector pointing upward, and **v** is the full and **u** the horizontal velocity vector, while *w* stands for the vertical velocity component. Furthermore, *p* and *b* are the density-weighted pressure deviation and the buoyancy deviation, respectively, from a reference atmosphere with stratification *N*^{2}(*z*), generally slowly varying in the vertical. These equations can capture essential aspects of local dynamics at various scales, as long as the vertical length scale of the waves is smaller than the atmospheric density scale height.

### a. Wave scaling

In the following a (resolved) mesoscale flow interacting with (unresolved) smaller-scale motions is considered, termed “submesoscale” for simplicity. The question arises which, if any, submesoscale motions are able to leave an impact on the mesoscale flow. Here this issue is addressed by considering possible interactions between a mesoscale and a submesoscale GW.

The mesoscale GW (subscript *m*) is taken to have horizontal and vertical length scales *L*_{m} and *H*_{m}, respectively, with an aspect ratio of the order

where *N*_{*} is a characteristic value of the Brunt–Väisälä frequency *N*, so that the mesoscale intrinsic frequency is equally affected by rotation and stratification, as follows from the general GW dispersion relation

Here, *N*^{2} = *O*() is meant to be the local value of the stratification, and *k* and *l* are the horizontal and *m* the vertical wave vector components, so that *L*_{m} = *O*[1/()^{1/2}] and *H*_{m} = *O*(1/*m*_{m}). These scaling assumptions are met, for instance, in the case where the horizontal and vertical scale are smaller by a synoptic-scale Rossby number than the Rossby deformation radius and the vertical scale height, respectively; in such situations one can take (*H*_{m}, *L*_{m}) = (1, 100) km (Achatz et al. 2017). The submesoscale wave (subscript *w*) has shorter vertical and horizontal scales *H*_{w} and *L*_{w}, respectively, satisfying

with *η* ≪ 1 a small parameter and *p* > 0, so that the submesoscale aspect ratio is of the order

The respective wave amplitudes are chosen as large as possible while keeping the analysis tractable. Hence, the submesoscale wave is assumed to be at the margin of static instability, where *db*/*dz* = *O*(*N*^{2}), so that the buoyancy-amplitude scale is

implying vertical displacements on the order of the submesoscale vertical scale. As for the mesoscale flow, if a marginally statically stable buoyancy amplitude is assumed, it turns out that the submesoscale-wave frequency *ω* = + **k**_{wh} ⋅ **u**_{m}—with **k**_{wh} indicating the horizontal submesoscale wave vector and **u**_{m} representing the mesoscale horizontal wind—is dominated by the Doppler term due to the strong mesoscale-flow horizontal winds, while the intrinsic frequency is relatively small; as a result, submesoscale motions are mainly transported by the mesoscale flow in such a regime. For this reason, the mesoscale buoyancy-amplitude scale is restricted to satisfy

which actually agrees with the submesoscale buoyancy-amplitude scale.

The remaining scales of interest follow from the buoyancy equation, (2), and the continuity equation, (3), with the material derivative scaling with the intrinsic frequency, which provides the inverse time scale. Specifically, the buoyancy equation yields a mesoscale vertical-wind scale:

From the continuity equation follows a mesoscale horizontal-wind scale

and one also notes for later reference that the mesoscale time scale is *T*_{m} = 1/*f*. Likewise, one obtains for the submesoscale horizontal- and vertical-velocity scales

where Ω_{w} is the scale-dependent submesoscale intrinsic-frequency scale. From the dispersion relation, (5), to a good approximation,

Of these, the last, strongly nonhydrostatic regime is generally modulationally unstable (Sutherland 2001) and hence not considered here.

### b. Regimes of interaction between mesoscale and submesoscale motions

Decomposing the total flow into mesoscale and smaller-scale submesoscale motions, (**v**, *b*) = (**v**_{m}, *b*_{m}) + (**v**_{w}, *b*_{w}), it is expected that the latter can only influence the former via flux-convergence terms. For such an interaction to be possible, these terms must be of the same magnitude as (or larger than) the leading mesoscale terms in the governing equations. To meet this condition, submesoscale wave fields are considered, with scaling as introduced above, that are spatially modulated on the mesoscale, in response to the two-way interaction between mesoscale and submesoscale flow. As explained in the appendix, it is then possible to identify a sufficiently scale-separated regime, where submesoscale motions may interact significantly with mesoscale GWs. In this regime the small-scale GWs lie in the midfrequency range *f*/*N*_{*} < *a*_{w} ≤ 1, and their scale separation is obtained by setting *p* = 2.

Specifically, incorporating this finding in (7) leads to a submesoscale aspect ratio

For the midfrequency range *f*/*N*_{*} < *a*_{w} ≤ 1 this implies, under the requirement of a sufficiently strong scale separation *η* ≪ 1,

According to (6), this scaling bears a stronger scale separation in the horizontal than in the vertical

While, asymptotically, *q* ≥ 0 is a free parameter, for atmospheric applications where *f*/*N*_{*}= *O*(10^{−2}), a sufficiently scale-separated general scaling regime can be identified in the finite range 0 ≤ *q* ≤ 1. In view of (17) and (18), two characteristic limit cases thus arise:

The first (nonhydrostatic) limit case is the one also discussed by Tabaei and Akylas (2007). Here the scale separation is quite large and the mesoscale-wave amplitude that can be affected is rather small. The most interesting case for atmospheric applications is the second (hydrostatic) limit, for which the mesoscale-wave impact is the strongest. For instance, taking (*H*_{m}, *L*_{m}) = (1, 100) km and *f*/*N*_{*}= 10^{−2}, from (20) it is then found that (*H*_{w}, *L*_{w}) = (0.1, 1) km. Notably, this scale estimate is in good agreement with present-day local-area weather-forecast-code mesh distances (see section 1).

For later reference, Table 1 provides an overview of the scales deduced in this section. It is worth noting that the pressure scales follow from ∂*p*_{m,w}/∂*z* = *O*(*b*_{m,w}), which hold both in the hydrostatic and nonhydrostatic regime, so that

as can be also verified from the GW polarization relations. Finally, the submesoscale time scale is the inverse intrinsic-frequency scale

### c. Nondimensional equations, multiscale asymptotics, and WKB ansatz

In the next step, the scaling for the submesoscale GWs derived above (see Table 1) is used to nondimensionalize the governing equations, (1)–(3). After substituting

where the subscript *h* denotes the horizontal components (here of the position vector **x**), the dimensionless equation system reads

Next, “compressed” variables are introduced to describe the slow variations of the resolved mesoscale flow, as compared to those of the submesoscale flow,

and then a WKB ansatz to describe a locally monochromatic submesoscale wave with slowly varying amplitude, wavenumber, and frequency is used. For a generic variable *ξ* it reads

where (**X**, *T*) indicates the (slowly varying) amplitude and *ϕ*(**X**, *T*)*η*^{−2} the (rapidly varying) phase. Following Tabaei and Akylas (2007), the latter is defined as *ϕ*(**X**, *T*) = *ϕ*_{0}(**X**_{h}) + *ηϕ*_{1}(**X**_{h}, *Z*, *T*), so that the local horizontal wavenumber, vertical wavenumber, and frequency are

where ∇_{X} denotes compressed spatial derivatives. Finally, all fields are expanded in the small scale-separation parameter *η*, taking into account the scaling derived above, and using subscripts 0 and 1 for the mesoscale and submesoscale parts, respectively:

Note that in this ansatz both mesoscale-field and submesoscale-wave amplitude as well as spatial and temporal scaling are all given in terms of the scale separation parameter *η*. There is no separate amplitude parameter. Higher harmonics of the submesoscale waves are neglected, as they can be shown to not contribute at leading order due to the dispersive GW dispersion relation (Achatz et al. 2017).

### d. Order analysis

After inserting the multiscale asymptotic ansatz, (34), into the nondimensional equations, (26)–(28), all terms are sorted by equal powers of *η* and the phase factor . Terms without phase factor describe mesoscale dynamics, while those proportional to the phase factor yield information on submesoscale wave dynamics; other harmonics, like and higher, are not considered, as noted above. The following will be kept concise, as the procedure is standard (e.g., Achatz et al. 2010, 2017).

#### 1) Leading-order results

The leading order of the vertical momentum equation establishes that the mesoscale flow is hydrostatic,

while the leading-order submesoscale terms in the equations can be summarized as , where is the antihermitian coefficient matrix

with the well-known Kronecker delta *δ*_{q,0}; the intrinsic frequency , that is, the frequency relative to the mesoscale velocity; and

the vector of leading-order submesoscale wave amplitudes. Nontrivial submesoscale wave amplitudes require a vanishing determinant of , leading to either the balanced solution or the GW dispersion relation

The corresponding null vector yields the submesoscale wave-amplitude polarization relations

where defines the buoyancy amplitude relative to the margin of static instability.

While the leading-order horizontal wavenumber does not develop in time, vertical wavenumber and frequency do. From their definition and the dispersion relation

one obtains their prognostic eikonal equations

where *c*_{gz} = ∂Ω/∂*m* is the (intrinsic) vertical group velocity. No horizontal group velocities appear since, at the submesoscales considered, to leading order, energy is transported only vertically.

#### 2) Higher-order results

The next-to-leading (“second”) orders in *η* yield the following mesoscale-flow equations:

while the “third” order of the mesoscale part of the continuity equation reads

with . As expected, a submesoscale-wave impact only exists in the horizontal momentum equation.

With the definition

using the shortcut , the submesoscale wave terms of the equations give for the next-to-leading orders in *η* the equation set , where

contains the next-order wave amplitudes. The matrix has a nonvanishing null space and thus **r**_{q} may not project onto it. This amounts to , with for the complex conjugate transpose, yielding

where * denotes the complex conjugate, for the energy density

Using the dispersion relation, (38), and polarization relations, (39), the energy flux and the shear-production term are expressed as

Thus, the wave-action conservation equation is obtained from (48), yielding

where is the wave-action density.

### e. Redimensionalization

The essential meso- and submesoscale equations derived in section 2d are finally transformed back into the original coordinates and redimensionalized by applying the substitutions

with , and their energy density is given by

with the corresponding wave-action density .

The so-called ray equations, consisting of the eikonal equations, (41) and (42), as well as the wave-action density equation, (52), describe completely the submesoscale dynamics, while the mesoscale dynamics is governed by the buoyancy, continuity, and momentum equations, (35) and (43)–(45). After back-transformation to the original coordinates and redimensionalization, a coupled equation system for the interaction between meso- and submesoscale GWs is obtained. Specifically, the mesoscale prognostic equations are

where **V** = **U** + *W***e**_{z}. These equations are linear in the mesoscale variables since the corresponding wave amplitude is sufficiently low. All mesoscale-flow nonlinearities disappear in the asymptotic limit of small *η*. The submesoscale-wave forcing acts via the convergence of the vertical flux of pseudomomentum **k**_{h}, as one finds from (44) and (51). The submesoscale dynamics is given by

Vertical propagation of the waves is too fast for lateral propagation effects to matter within the time elapsing during the vertical propagation over a mesoscale vertical length scale. Note that (65) follows from (66) and the dispersion relation, (58), and hence is not an independent equation in the submesoscale model. The vertical group velocity *c*_{gz} may be calculated with the help of the current vertical wavenumber *m*.

Finally, it is to be mentioned that the coupled equation system (61)–(67) is energy conserving. For the mesoscale energy density *ε*_{m} = 1/2(|**U**|^{2} + *B*^{2}/*N*^{2}), one finds from (61)–(64) that its tendency obeys

while it can be seen from the dimensional version of (48), (50), and (51) that the evolution of submesoscale energy density is governed by

Therefore, the local total energy density *ε*_{t} = *ε*_{m} + *ε*_{w} satisfies

As a result, if there is no mesoscale pressure flux at all boundaries, and no submesoscale energy flux *ωc*_{gz} at the vertical boundary of the domain, or if periodic boundary conditions hold, the spatially integrated total energy density *E*_{t} is conserved:

## 3. Description of the numerical models

In this section the numerical code used for validation tests is described: the WKB code PincFloit–WKB is an implementation of the theory presented above. PincFloit without WKB submesoscale-wave model, but instead in a setting explicitly resolving the submesoscale waves, is used for large-eddy simulations (LESs) to provide data against which to validate the WKB model as well as its underlying theory. Since the GW dynamics is invariant with regard to rotation of the horizontal coordinate system, both codes have been used in 2D mode.

### a. The PincFloit–WKB model

#### 1) Submesoscale flow: The Lagrangian WKB model

The numerical implementation of the interaction between submesoscale and mesoscale flow is achieved by coupling a Lagrangian phase-space ray tracer (Muraschko et al. 2015) to the mesoscale-resolving model PincFloit. As can be read directly from the submesoscale GW equations, (66) and (67), along rays satisfying

vertical wavenumber and wave-action density develop according to

while the frequency can be obtained from the wavenumber and local mesoscale flow by the dimensional version of the full dispersion relation, (40). In (73) the replacement *δ*_{q,0} → 1 has been done, since in the case *q* = 1, one has *k*^{2} ≪ *m*^{2} anyway. Direct implementation of this model leads, however, to serious numerical instabilities (Tabaei and Akylas 2007; Rieper et al. 2013a). These are due to caustic situations where rays cross, leading to multivalued wavenumbers and wave-action densities. This can be avoided by taking a spectral perspective (Muraschko et al. 2015), where (72) and (73) indicate movement through a phase space spanned by vertical position and wavenumber. Wave-action density (**x**, *t*) is then replaced by the spectral phase-space wave-action density (**x**, *m*, *t*) (e.g., Dewar 1970; Olbers 1976; Bühler and McIntyre 1999; Hertzog et al. 2002; Muraschko et al. 2015), developing along the rays according to

with *d*_{r}/*dt* = ∂/∂*t* + *c*_{gz}∂/∂*z* + ∂/∂*m* as the phase-space material derivative. The physical space wave-action density can be retrieved from it by the wavenumber integral (**x**, *t*) = (**x**, *m*, *t*), and the submesoscale momentum flux, for example, is obtained from

This requires reconstructing the full phase-space dependence of from infinitely many rays. In a first discretization step 1 therefore collects rays carrying nonzero wave action in a number of rectangular ray volumes with constant . In principle the individual ray velocities will deform these ray volumes arbitrarily strongly. In a second discretization step, this deformation is simplified by prescribing the ray volumes to keep a rectangular shape. More details on this and the corresponding momentum-flux reconstruction are given by Muraschko et al. (2015) and Bölöni et al. (2016).

The submesoscale momentum flux as well as energy density are smoothed after regridding over a window of 3 × 3 PincFloit finite-volume cell equivalents [see sections 3a(2) and 3a(3) below] in order to avoid artificial peaks resulting from sampling problems due to the ray discretization. Moreover, in the simulations the effect of submesoscale horizontal group velocity is indeed small compared with that of the vertical group velocity, as found in theory, but still noticeable in the comparisons with the LES. Therefore, this effect has been incorporated by allowing the ray volumes to also propagate in the horizontal direction for several simulations,

implying as well a generalized phase-space material derivative *d*_{r}/*dt* = ∂/∂*t* + **c**_{g} ⋅ ∇ + ∂/∂*m*, with **c**_{g} the 3D group velocity. Finally, the ray tracer has been supplemented by a simple saturation scheme (Bölöni et al. 2016) to account for turbulent wave breaking. The wave-action density of the submesoscale GW packet is locally reduced, when its amplitudes reach the upper limit of static stability. Results show that the saturation is important for the total energy budget. In the present Boussinesq context, however, it has not contributed significantly to the instantaneous wave field distribution and simulation results discussed below.

#### 2) Mesoscale flow: PincFloit

The pseudoincompressible flow solver with implicit turbulence modeling (PincFloit), originally developed by Rieper et al. (2013b) to solve the pseudoincompressible equations of Durran (1989), modified appropriately to integrate the Boussinesq equations, (1)–(3), has been used at mesoscale resolution to simulate the resolved mesoscale flow. To account for the impact of the unresolved submesoscale waves, the momentum equation has been supplemented by the corresponding convergence of horizontal pseudomomentum flux, as indicated by (64) and (76):

The latter is provided by the Lagrangian WKB code described above. As adumbrated by (78), in PincFloit the leading-order mesoscale dynamics is identified with the full resolved mesoscale nonlinear flow. Technically, PincFloit uses a finite-volume discretization with a staggered grid. Time integration is performed by an adaptive third-order Runge–Kutta scheme with a CFL criterion. Pressure is computed, using the nondivergence constraint, (61), by solving the corresponding Poisson equation. The latter is done using a biconjugate gradient stabilized method (BiCGSTAB) method (van der Vorst 1992). More details can be found in Rieper et al. (2013b).

#### 3) The coupled PincFloit–WKB model

PincFloit and the Lagrangian WKB model are coupled interactively, so as to simulate the transient interaction processes of resolved mesoscale and unresolved submesoscale GWs, as derived in section 2. At every Runge–Kutta substep, information is exchanged between the mesoscale and the submesoscale dynamics. The Lagrangian WKB model determines the momentum flux convergence of the submesoscale waves via the discretization of (76) and updates the mesoscale wind field, which is then delivered to PincFloit. Hereafter, the latter integrates the Boussinesq equations, (1)–(3), at mesoscale resolution. After that, the new wind and background values are provided to the Lagrangian WKB model, which solves the ray-tracing equations, (72), (73), and (75), [as well as (77), if intended], yielding an updated submesoscale wave momentum flux and thus closing the circle.

Another remark is that this coupled system conserves the sum of mean flow and wave energy too. From the Boussinesq equations, (1)–(3), it can be derived for the mesoscale-flow energy density *ε*_{m} = 1/2(|**V**|^{2} +*B*^{2}/*N*^{2}) that

while Bölöni et al. (2016) have shown for the submesoscale wave energy density that

so that total energy ∫*ε*_{t} *d *^{3}*x* = ∫(*ε*_{m} + *ε*_{w}) *d *^{3}*x* is conserved under suitable boundary conditions (zero or periodic) for the respective fluxes.

### b. PincFloit–LES

In LES mode PincFloit is used to integrate the fully nonlinear Boussinesq equations, (1)–(3), with the above WKB submesoscale wave model switched off. Its resolution is chosen fine enough that the initial submesoscale wave field is completely resolved, and that it captures wave–wave interactions and interactions between all waves and the larger-scale turbulent eddies. Motivated by the results from corresponding benchmark tests (Remmler et al. 2015), small-scale turbulence is not parameterized by the implicit adaptive local deconvolution method (ALDM; see, e.g., Hickel et al. 2006), as originally implemented into PincFloit, but by a dynamic Smagorinsky method (Germano et al. 1991). The corresponding Smagorinsky coefficient is averaged over a local spatial window of 5 × 5 finite-volume cells so as to stabilize the scheme.

## 4. Numerical experiments

PincFloit–WKB and PincFloit–LES were used to simulate the propagation of a spatially confined wave packet in a uniformly stratified (*N* = 0.02 s^{−1}) atmosphere on an *f* plane (*f* = 10^{−4} s^{−1}, except in the test case COR, where *f* = 0) with zero initial ambient flow. It is well known from longwave–shortwave interaction theory (Tabaei and Akylas 2007; Van den Bremer and Sutherland 2014) that such a packet of small-scale waves is able to generate a mean flow consisting of mesoscale wave structures connected to a resonance phenomenon (see also section 4c). In turn, the mesoscale waves may have an influence on the propagation of the wave packet. All simulations, investigating the resonant behavior of various submesoscale wave packets, are two-dimensional. The horizontal *x* axis is chosen to point into the direction of **k**_{h} = *k***e**_{x} and all initial fields are only dependent on *x* and *z*, as is then also the case henceforth. All models use periodic boundary conditions in *x* and *z*.

### a. Initialization

Consequently, a locally monochromatic wave packet with horizontal and vertical wavenumbers *k* and *m*_{0}, respectively, is initialized. It is vertically as well as horizontally confined, with a Gaussian envelope amplitude characterized by the standard deviations *σ*_{x} and *σ*_{z}. Its buoyancy field is thus

Herein the amplitude parameter is chosen so that at infinite *σ*_{z} static stability is given at *x* = *x*_{0}, that is, *N*^{2} + ∂*b*/∂*z* > 0 , for values . The initialization fields for the LES model are determined from the full rotational GW polarization relations (see, e.g., Bühler (2009), section 8.2)

which lead, by the way, to the same initialization for a 2D wave packet with *l* = 0 as if one would use the leading-order polarization relations, (59), except for the meridional wind component *ν*, yielding

where denotes the intrinsic frequency of the initially monochromatic wave packet with wavenumbers *k* and *m*_{0}, including rotational effects. In line with the hydrostatic submesoscale-wave scaling regime, (20), defined above, the choice falls on *λ*_{x} = 2π/*k* = 1000 m and *λ*_{z} = −2π/*m*_{0} = 100 m, based on typical mesoscales like 2*σ*_{x} = *L*_{m} = 100 km and 2*σ*_{z} = *H*_{m} = 1 km (except for the case SCALE in section 4d). A typical mesoscale time scale for the wave packet test case is given by *T*_{m} = *f *^{−1} = 10 000 s.

In PincFloit–WKB the initial fields of the wave packet are to be defined for the Lagrangian WKB model. Following Muraschko et al. (2015), a certain set of ray volumes, each carrying specific wave properties, is “placed” in a rectangular 2.5*σ* environment around the initial center of the wave packet in both *x* and *z* directions, covering more than 98.7% of wave energy density, as shown in Fig. 1. The initial phase-space wave-action density of the quasi-monochromatic wave packet is assumed to be

where *m*_{0,1} = *m*_{0} − Δ*m*_{0}/2 and *m*_{0,2} = *m*_{0} + Δ*m*_{0}/2. It is then discretized by rectangles in physical space, and in wavenumber space by two wavenumber intervals each centered at wavenumbers *m* = *m*_{0,1} + Δ*m*_{0}/4 and *m* = *m*_{0,2} − Δ*m*_{0}/4, constituting together phase-space ray volumes within which is taken to be constant. After initialization the ray volumes propagate through the phase space in accordance with the ray equations, (72), (73), and (77), changing in response to the mesoscale flow physical location and wavenumber. Hence, an initially quasi-monochromatic wave packet can develop quite complex spectra, as also discernible from Fig. 6 below.

Resolution is chosen as follows: in PincFloit–LES around 33 grid points in the *x* direction and 20 grid points in the *z* direction per initial submesoscale wavelength are set. For PincFloit–WKB, it turned out that 10 grid points per typical mesoscale length scale are sufficient in most of the test cases (see Table 3). Parameters describing the general setup are listed in Table 2.

### b. Postprocessing

While PincFloit–WKB outputs directly the submesoscale GW momentum flux, their energy density, and mesoscale wind and buoyancy fields, the output fields of PincFloit–LES contain both meso- and submesoscale information, which have to be separated. Instead of a wavenumber filter as applied by Van den Bremer and Sutherland (2014), a running mean over two initial submesoscale wavelengths in each spatial direction filters submesoscale contributions from the mesoscale ones. The former are then obtained by subtracting the mesoscale part from the full field and subsequently, momentum flux and wave energy density are calculated via

at each grid point (*x*_{j}, *z*_{i}) for each output time *t*_{k}, where the overbar indicates an additional running mean over two wavelengths. There is no significant sensitivity to the choice of the average interval when comparing the results for an average over one up to three (and more) initial submesoscale wavelengths.

### c. Resonant interaction of mesoscale with submesoscale GWs

It is known that energy exchange between long and short GWs is particularly strong, when a packet of small-scale GWs interacts with its induced large-scale mean flow under resonance conditions. Grimshaw (1977) found such resonant behavior for a modulated wave train propagating along a horizontal channel, whereas Sutherland (2001), Tabaei and Akylas (2007), and Van den Bremer and Sutherland (2014) investigated wave-mean flow interactions for spatially localized wave packets propagating vertically through an unbounded stratified Boussinesq fluid. Tabaei and Akylas (2007), in particular, pointed out that flat wave packets, characterized by a stronger modulation in the vertical than in the horizontal, can lead to resonant forcing of large-scale waves. This resonance arises when the wave packet modulation scales are compatible with those of free, hydrostatic inertia GWs that are natural mode solutions of the Boussinesq system. Furthermore, in the small-amplitude limit, where to leading order the wave packet envelope propagates vertically as a wave of permanent form, this resonance singles out inertia GWs whose vertical phase velocity matches the vertical group velocity of the wave packet. Tabaei and Akylas (2007) further speculated that this mechanism might be responsible for the generation of inertia GWs in the real atmosphere.

As noted in section 2b, the scaling regime considered by Tabaei and Akylas (2007) corresponds to the nonhydrostatic limit, (19), which, under atmospheric conditions, exhibits a rather strong scale separation *L*_{w}/*L*_{m} = *O*(10^{−4}) and *H*_{w}/*H*_{m} = *O*(10^{−2}); as a result, the associated submesoscale waves would be quite short: *H*_{w} = *L*_{w} = *O*(10) m. The hydrostatic regime, identified in (20), however, features submesoscale GWs with length scales about one order of magnitude smaller than typical gridpoint distances of nowadays global NWP models and of the same magnitude as the distances in regional limited-area models (see section 4a). Consequently, in the latter such submesoscale waves reside on the largest unresolved scale, and the long–short GW resonance of Tabaei and Akylas (2007) is relevant to this regime.

The resonance condition noted above for small-amplitude wave packets can be written in the present notation as

where is the (intrinsic) vertical group velocity of the submesoscale wave packet, and is the vertical phase velocity of the induced, mesoscale inertia GWs. Van den Bremer and Sutherland (2014) used this condition to estimate the phase line tilt of the generated large-scale waves in a nonrotational atmosphere. Here, based on (89), the geometry of the induced mesoscale inertia GW disturbance in the presence of rotation is discussed, following a kinematic approach analogous to that of the classical Kelvin ship wave pattern (see Whitham 1974, section 12.4).

The resonance condition, (89), combined with the dispersion relation, (5), implies that the horizontal wavenumber *k*_{m} can be expressed in terms of the vertical wavenumber *m*_{m} of the induced mesoscale GWs [which are hydrostatic according to (4)]:

As noted earlier, this assumes that (i) the submesoscale wave packet envelope, which acts as forcing, propagates as a wave of permanent form at constant speed , and (ii) the induced mesoscale waves are steady in the frame of the moving source. The validity of these assumptions will be tested numerically in section 4d.

In addition, one generally has ∂*m*_{m}/∂*x* = ∂*k*_{m}/∂*z*, which combined with (90) yields

This equation for *m*_{m} can be treated by the methods of characteristics. In the far field, all characteristics (straight lines) originate from the source (which appears as a point, *x* = *z* = 0, say, in the moving frame, with *z* > 0 behind the source), and it is found that

Equation (92) determines *m*_{m} for given *z*/*x*, and the corresponding *k*_{m} follows from (90). It should be noted that |*z*/*x*| → ∞ as |*m*_{m}| → *f* /; also, when |*m*_{m}| → ∞, one has |*z*/*x*| ≈ 2|*m*_{m}|/*N* → ∞. This suggests that there must be limiting characteristics (*z*/*x*)_{lim}, where

From this one obtains the critical wavenumber |*m*_{m,lim}| = (3/2*f*^{2}/)^{1/2}, which along with (92) gives (*z*/*x*)_{lim} = ∓2^{3/2}*f* /*N*. Thus, a Λ-like wave pattern behind the moving wave packet is expected as in the case of a ship wake. The opening angle *α* defined by the limiting characteristics is given by

As can be seen from (19), (20), (24), and (29), in compressed coordinates the opening angle is 2tan^{−1}[1/2^{3/2}] ≈ 38.9°, independently of the scaling regime. Remarkably, this angle is identical to that of the classical Kelvin ship wave pattern on deep water (Whitham 1974).

Solving the characteristic equation, (92), for gives

Generally, according to (95) combined with (90), there are four different possible modes: two shortwave [corresponding to the plus sign in (95)] and two longwave [corresponding to the negative sign] modes. The former are not relevant here, as their length scale is not in agreement with the mesoscale derived in our theory. Moreover, one of the two longwave modes (where the phase contributions *k*_{m}*x* and *m*_{m}*z* have opposite signs) does not obey the radiation condition, so that wave energy would not be radiated away from the source. Thus, only one mode (where *k*_{m}*x* and *m*_{m}*z* have the same sign) is expected to prevail in the simulations (see Fig. 2), consistent also with Fig. 1 in Tabaei and Akylas (2007).

### d. Test case studies

Theory and its numerical realization by the PincFloit–WKB model are validated in several test case studies: a reference simulation REF has been performed first, investigating the propagation of a relatively high-amplitude wave packet () in a uniformly stratified atmosphere initially at rest. After that, the initial amplitude is varied to both lower and higher amplitudes (test cases AMPx) in order to demonstrate the robustness of the Lagrangian WKB model for a wide range of amplitudes of submesoscale waves. Then a wave packet pushing the limits of the asymptotic scaling (see Table 1) as well as the WKB (see section 2c) assumptions is initialized, thus establishing the wide applicability of the PincFloit–WKB model and the range of validity of the theory (test case SCALE). The test case COR examines the sensitivity on the ratio *f*/*N* by setting *f* = 0; strictly speaking, this regime is not represented by the theory as *η* = 0. Physically, it is to be understood as a limit case for a tropical background. Last, the test case PSINC is used to compare LESs of the REF wave packet in a pseudoincompressible atmosphere to the chosen Boussinesq dynamics. PSINC will be referred to in the discussion in section 5. Table 3 provides an overview of the essential test case-specific model setup and namelist parameters. According to (20), the chosen background parameters *N* = 0.02 s^{−1} and *f* = 10^{−4} s^{−1} imply *α* = *η* ≈ 0.07 = *O*(0.1).

Wave packets have been initialized in accordance with (82) and (84)–(86). Except for the SCALE test case, the spatial extent of the initialization box (roughly associated with the wave packet size) in PincFloit–WKB is 250 km × 2.5 km. The quantities shown in the following figures are in dimensional units. The spatial coordinates are given in kilometers, whereas the time axis is scaled by the inverse of the inertial frequency *f*. The 2D LES fields are visualized after mapping them to a 512 × 512 grid, with corresponding domain dimensions 500 km × 10 km, except for the case SCALE with 1000 km × 5 km.

Figure 3 shows the initial condition in the REF test case; the distribution of wave energy density [calculated via (88)] is displayed, where one can easily see that most of the submesoscale energy of the wave packet is contained in a range of roughly 1.5*X* × 1.5*Z * 150 km × 1.5 km with the compressed coordinates *X* and *Z* from (29). Since the wave packet amplitude is (and one may lose or gain factors as, e.g., 2π in the scaling procedure), the wave packet energy density is somewhat smaller than theoretically assumed.

For the test case REF, three slightly different PincFloit–WKB model configurations are set up. First, a single-column ray tracer with 1D spatial ray propagation, that is, only in the vertical, is initialized and no feedback of the resolved flow onto the submesoscale wave packet is allowed [i.e., the right-hand sides of both (73) and (77) equal zero; shortcut: PincFloit–WKB–1DNF]. Second, a single-column ray tracer is used that accounts for two-way scale interactions [i.e., only the right-hand side of (77) equals zero; shortcut: PincFloit–WKB–1D]. Note that, although ray-volume propagation is 1D, the ray-volume amplitudes in are not [see (86)], so that the submesoscale wave packet induces a 2D mesoscale response. Third, horizontal ray volume propagation is also allowed (shortcut: PincFloit–WKB–1.5D).

After a bit less than one inertial period (*t* ≈ 17.3 h), one has a couple of findings. In the case of a non-energy-conserving system with no coupled interaction (Fig. 4a), one observes horizontally more prolonged mesoscale wave structures than in the validating LES (Fig. 4d). The colored contours represent the wave-packet-induced mesoscale horizontal winds. Furthermore, there is no wave packet deformation with PincFloit–WKB–1DNF, which is a result of the feedback process (to be compared with Fig. 4b). Horizontal propagation of the wave packet is rather weak compared to the vertical in terms of the characteristic mesoscales. Incorporating (77), with PincFloit–WKB–1.5D (Fig. 4c), a picture, which is qualitatively as well as quantitatively very close to the LES, is obtained.

The overlay of the theoretically derived wave structures of Fig. 2 on Fig. 4a shows that indeed, the resonance condition examined in section 4c is able to explain the lateral confinement of the induced waves: the mesoscale wave patterns are reminiscent of modes that develop in the wake of a ship, and the spatial extent of them is in good agreement to the area bounded by the limiting characteristics—except close to the wave packet, as it is no point source as assumed in the theoretical derivation. Interestingly, self-acceleration effects cause a lateral constriction of the induced structures close to the wave packet, so that the feedback-allowing simulations resemble the prediction even more with regard to the limiting characteristic. The vertical wavelengths are predicted by (95) to be somewhat longer than observed in these simulations and the phase line tilt is expected to be less pronounced; these small differences are explainable by the application of the far field approximation, the disregarded feedback from the induced waves and hence the requirement of a constant group velocity of the wave packet in section 4c, as underscored by the longer wavelengths in Fig. 4a. Nevertheless, the equations derived there seem to be a very good indicator for the prediction of the induced wave mode.

Furthermore, the induced structures are evocative of the ones induced by the small-amplitude nonhydrostatic wave packet of Tabaei and Akylas (2007): they resemble plain downward-propagating inertia GWs. Their period—extracted from the Hovmöller diagram in Fig. 5—is very close to the inertial period 2π/*f*.

It shall be stressed again that—in contrast to the assumption of Van den Bremer and Sutherland (2014)—within the first inertial period, the wave packet energy density distribution varies already significantly, as the group velocity undergoes partially strong and rapid changes: the peak values of energy density are found in the lower sector of the wave packet (not shown). When mesoscale wind and vertical wind shear build up, the vertical wavenumber—and due to (72) vertical group velocity too—of a part of the submesoscale waves undergoes considerable changes according to (73), as found by investigating the time evolution of the discrete wavenumber spectrum (Fig. 6). The characteristics and phase lines of the wave mode predicted in section 4c might hence slightly differ from Fig. 2 or Fig. 4a, respectively, when abandoning the requirement of a constant group velocity. Notably, in contrast to Tabaei and Akylas (2007), it can be reported that PincFloit–WKB is stable all the time, although caustics are ubiquitous in physical space as can be seen in Fig. 6b. Similar to Bölöni et al. (2016), it is also found (not shown) that the model conserves total energy very well (e.g., not more than 2% variation in REF at the end of the simulation time), as predicted in sections 2e and 3a(3).

The AMPx test cases substantiate the finding that the refractive feedback onto the wave packet becomes larger, the larger its (initial) amplitude is, as one can clearly see in the Figs. 7 and 8. In other words, as one moves toward small amplitudes, meso- and submesoscales interact only weakly nonlinearly and the wave packet propagation is very well approximated by its initial group velocity. Consequently, the difference between the results of PincFloit–WKB–1DNF (Fig. 7a) and PincFloit–WKB–1D (Fig. 7b) is considerably smaller than in the REF test case.

Induced mesoscale momentum amplitudes, however, are about one order smaller in the small-amplitude case AMP1 compared to the REF simulations, well in proportion to the reduced momentum forcing by the submesoscale wave packet. While the frequency of the induced mesoscale waves remains close to the inertial period independently of wave packet amplitude, their vertical extent evolving during the wave packet propagation becomes larger the less refraction appears, as in the case of a small-amplitude quasi-steady propagating disturbance, where the mesoscale wind is too weak to influence the wave packet.

In the case SCALE, a very flat wave packet of hydrostatic submesoscale GWs is initialized, whose ratio of vertical to horizontal amplitude-variation length *σ*_{z}/*σ*_{x} is decreased by one order of magnitude, but whose total energy is kept constant compared to REF. Vertical amplitude variation is thus much stronger and horizontal amplitude variation weaker than originally assumed in theory. Beyond that, the vertical amplitude variation scale is now close to the vertical submesoscale wavelength. Even though the horizontal scale separation is in accordance with the general regime, this is a case at the limit of the underlying theory.

Figure 9 shows again the wave packet energy density and the induced mesoscale horizontal momentum fields as in Fig. 4. Apart from the partly bandlike representation of the wave packet due to a limited set of ray volumes—allowing a passable computing time—in PincFloit–WKB, which experiences strong wind shear changes, and a slight overestimation of the mesoscale amplitudes by PincFloit–WKB, the results from the parameterized model agree well with the LES model. An inertia GW is generated whose resonant amplification is broken off after approximately one inertial period simulation time because of wave packet energy spreading. The elongation of the induced mesoscale wave is hence smaller in the vertical than, for example, in AMP1 after that time, like in the higher-amplitude cases REF or AMP2 after approximately two inertial periods, when the wave packet has experienced strong refraction effects.

The test case COR, where *f* = 0 in contrast to the case REF, investigates the sensitivity of the implementation in PincFloit–WKB on the value of *f* /*N* (Fig. 10). Energy is radiated stronger and further laterally than in the rotational REF case, which is reflected by the induced wind structures. This finding has been carved out in detail by Van den Bremer and Sutherland (2014) and Tabaei and Akylas (2007), and underscores again the wide applicability of PincFloit–WKB.

## 5. Résumé, discussion, and outlook

The first question in our study was whether there are any submesoscale GWs that can modify mesoscale dynamics at typical scales. This has been investigated within Boussinesq dynamics for typical midlatitude tropospheric values of *f*/*N*, using scale analysis comparing the magnitude of large-amplitude submesoscale GW flux convergences with the self-consistent acceleration and heating in statically stable mesoscale GWs. A range of scales of high- and midfrequency submesoscale GWs has been identified where an impact on the mesoscale waves is possible. It encompasses the case of very small-scale high-frequency GWs investigated by Tabaei and Akylas (2007), but even more interesting is the appearance of hydrostatic GWs with scales that are just below the resolution of present-day limited-area numerical-weather-prediction codes. For these waves the vertical-scale separation is *η* = (*f*/*N*)^{1/2}, while the horizontal-scale separation is *η *^{2}.

Using multiscale asymptotics, a large-amplitude WKB theory for the interaction between locally monochromatic submesoscale GWs and a mesoscale flow has been derived. All nondimensional fields are expanded in terms of powers of *η*, and then the distinguished limit of small *η* is taken. This leads to separate but coupled equation systems, one describing explicitly resolved mesoscale dynamics, and the other depicting submesoscale dynamics, consisting of vertical ray equations for the submesoscale wave properties. Direct coupling is established on the one hand by a modification of the mesoscale-momentum equation through the vertical convergence of submesoscale pseudomomentum flux; and on the other hand by a change of submesoscale vertical wavenumber and frequency by mesoscale wind shear regions, or, though unaffected by submesoscale fluxes, the vertical variation of background stratification. An interesting and useful result might seem to be that, other than in the case of the interaction between synoptic-scale flow and mesoscale inertia GWs, neither is horizontal propagation of the small-scale GWs a leading-order effect, nor are small-scale GW horizontal flux convergences. Hence, single-column approaches to submesoscale GW modeling in mesoscale-resolving models appear well justified at first sight. This seems to be in contrast to previous findings about the relevance of lateral propagation of mesoscale waves by Senf and Achatz (2011), Sato et al. (2012), and Plougonven et al. (2017), suggesting certain limitations of single-column approaches. Indeed, caution is at place to not misinterpret results. To see this, note that the ratio between vertical and horizontal group velocity is |*c*_{gz}/|**c**_{gh}|| = ||**k**_{h}|/*m*| = *H*_{w}/*L*_{w}. Hence, the horizontal distance *L* covered by a wave packet propagating in the vertical over a distance *H* is *L* = *HL*_{w}/*H*_{w} = *H*(*N*_{*}*/f* )^{1/2}, where the scalings (4) and (20) are used. The secondary importance of lateral propagation seen here is due to the fact that, if *H* = *H*_{m} is taken then *L*/*L*_{m} = *H*/*L*_{m}(*N*_{*}*/f* )^{1/2} = (*f* /*N*_{*})^{1/2} ≪ 1, again using (4). Hence, the submesoscale wave packet travels over a considerably less horizontal distance than the scale characterizing mesoscale variations. On the other hand, the ultimate justification for single-column implementations would be that for all possible vertical distances coverable, at most the vertical model extent *H* = *H*_{top}, the horizontal distance *L* covered should be less than a horizontal mesh distance Δ*x*_{h}. This would imply *L*/Δ*x*_{h} = *H*_{top}/Δ*x*_{h}(*N*_{*}*/f* )^{1/2} ≪ 1. Typically, this condition cannot be met. It hence appears safer to take lateral propagation into account, and our simulations also show improved results if one does so.

The validity of the theory has been examined by its implementation into a mesoscale-resolving Boussinesq model. The WKB fields of submesoscale GW amplitudes and wavenumbers are discretized and predicted by a Lagrangian ray tracer (Muraschko et al. 2015; Bölöni et al. 2016) that uses a spectral phase-space representation, thereby avoiding numerical instabilities due to caustics (Tabaei and Akylas 2007; Rieper et al. 2013a) and also allowing the potential development of a spectral submesoscale GW field from initially locally monochromatic conditions. The approach is validated against simulations by a submesoscale-resolving large-eddy code.

Test cases launching a 2D Gaussian wave packet of hydrostatic submesoscale GWs furnish evidence of the applicability of the numerical approach and its underlying theory. A resonance effect occurs, which has been found in previous studies of Tabaei and Akylas (2007) and Van den Bremer and Sutherland (2014), and leads to the generation of mesoscale inertia GWs. Their characteristics correspond to free modes of the Boussinesq dynamics, which can be explained by a theoretical study of the resonance condition. The inertia GWs can have a strong impact onto submesoscale wave packets by refraction, depending on wave amplitude, and eventually cause a saturation or breakup of the effective resonant energy transfer. Further case studies show that the approach, although designed for large-amplitude submesoscale GWs, also works for low-amplitude submesoscale GWs, and that it also performs reasonably well when the scale separation between mesoscale and submesoscale is weaker than assumed in the theoretical derivations. Beyond that, the code seems to be usable for nonrotating cases as well.

The Boussinesq setting of our analysis does not say that non-Boussinesq effects are irrelevant. The vertical decrease of ambient density would play an important role in operational weather forecast and climate models. This is left out here for the mere sake of simplicity, but corresponding generalizations as in Bölöni et al. (2016) seem straightforward. Moreover, in an analysis based on fully compressible dynamics, Achatz et al. (2017) identify compressible and elastic effects in the interaction between near-inertial mesoscale waves and synoptic-scale flow. One of these are synoptic-scale pressure fluctuations that matter in a strongly stratified atmosphere. The other is an elastic mesoscale-wave term appearing in the synoptic-scale momentum equations. That term is most relevant for near-inertial waves but loses importance in the noninertial frequency range. The most interesting submesoscale waves, which are identified in the present study, are in the latter range. Hence, a fully compressible treatment seems a necessary extension of the present investigations based on Boussinesq theory; however, it is not clear, whether it will identify, at the scales of interest here, relevant non-Boussinesq effects, beyond those resulting from the ambient-density vertical dependence. For a first hint the REF case has been simulated with PincFloit–LES, but using it in the pseudoincompressible mode (Durran 1989; Rieper et al. 2013b) instead of the Boussinesq mode (test case PSINC). pseudoincompressible dynamics captures the elastic effects arising in the analysis of Achatz et al. (2017) in the momentum equation. Figure 11 shows a snapshot of the simulation, to be compared to Fig. 4. Apparently, at least in this case, elastic dynamics does not seem to be of leading-order importance.

Our analysis shows that submesoscale GWs can significantly influence a mesoscale flow, provided their amplitudes are large enough. It also derives a theory and its numerical implementation for the efficient representation of such effects in mesoscale-resolving models. Whether submesoscale GWs at the required scales and amplitudes are indeed present to a sufficient degree in the atmosphere, and by which processes they can be generated, certainly is another question that should be investigated in the future. Mesoscale wind-field spectra determined from aircraft data by Callies et al. (2014) and Bierdel et al. (2016) do not exhibit dissipation at the smallest observable scales so that one could imagine a continuation of these well into the submesoscale range. A further indication is that GWs at very small scales seem to be a relevant issue in the stable planetary boundary layer (Sun et al. 2015). Should corresponding studies yield further support for the relevance of submesoscale GWs, an unavoidable next step would have to be extending the Boussinesq theory to an analysis of the compressible Euler equations.

## Acknowledgments

U.A., R.K. and J. Wei thank the German Research Foundation (DFG) for partial support through the research unit Multiscale Dynamics of Gravity Waves (MS-GWaves) and through Grants AC 71/8-1, AC 71/9-1, AC 71/10-1, and KL 611/25-1. U.A. and B.R. thank the German Federal Ministry of Education and Research (BMBF) for partial support through the program Role of the Middle Atmosphere in Climate (ROMIC) and through Grant 01LG1220A. B.R. acknowledges funding support from the CEA and the European project ARISE 2 (Horizon 2020, GAN653980). U.A., T.R.A. and J. Wilhelm thank MISTI (MIT-Germany Seed Fund) for making their collaboration possible.

### APPENDIX

#### Derivation of a General Scaling Regime

As motivated at the beginning of section 2b, in order to identify a general scaling regime for the interaction between meso- and submesoscale GWs, the submesoscale flux convergence terms must be comparable to the leading mesoscale terms in the governing equations, (1)–(3). Specifically, the relevant convergence terms scale as

where one has additionally made use of the GW polarization relations as given in (83). For example, for the buoyancy-flux estimates, because of the phase shift between momentum and buoyancy, the fluxes obey the relations = *O*(*fL*_{w}*B*_{w}) = *O*(*f*/Ω_{w}*U*_{w}*B*_{w}) and = 0, with the tilde indicating the submesoscale wave amplitudes and the asterisk denoting the complex conjugate. The convergence terms above are to be compared to the horizontal acceleration in the horizontal momentum equation, the buoyancy time derivative in the buoyancy equation, and the buoyancy in the vertical momentum equation for the mesoscale flow. These terms scale as

Using these scalings, the following ratios for the horizontal momentum equation are found [always neglecting the modulationally unstable strongly nonhydrostatic submesoscale regime in (14)]:

The comparisons for the vertical momentum equation yield

and

and for the buoyancy equation one obtains

One sees that the submesoscale wave field can impact the mesoscale flow only via the horizontal momentum equation. Two options exist: the first is the impact of the low-frequency submesoscale waves via horizontal momentum-flux convergence [see (A10)]. Equality between this flux term and the mesoscale-flow horizontal acceleration is reached when, using (6), 1 = (1/*η*)(/) = *η*^{2p−1}. Hence, in this regime *p* = 1/2. For an appreciable scale separation in the horizontal [see (6)], one would expect, say, *η*^{1/2} = *O*(10^{−1}). This, however, implies very small *η*, so that the mesoscale-wave amplitude, (9), that can be affected is very low. Without showing this in any detail, it is also mentioned that in this case the submesoscale-wave frequency is dominated by the intrinsic part, while the Doppler term is small. This means that the mesoscale-flow impact on these waves is rather weak.

The second and more interesting option involves the impact of the midfrequency submesoscale waves via vertical momentum-flux convergence [see (A11)]. Equality between flux convergence and the mesoscale-flow horizontal acceleration is reached here, using (4) and (6), when 1 = *ηa*_{w}*N*_{*}/*f* = *η*^{2−p}. Hence in this regime *p* = 2, and the horizontal length scale separation *η*^{2} and the vertical length scale separation *η* can be small with *η* = *O*(10^{−1}), say, so that the mesoscale-wave amplitude that can be affected is stronger than in the first option. The possible range of *η* can be determined from the condition that the submesoscale waves shall be in the midfrequency range so that *f*/*N*_{*} < *a*_{w} ≤ 1. Together with (7) and the requirement of a sufficiently strong scale separation *η* ≪ 1, this implies (16).

## REFERENCES

*Waves and Mean Flows.*Cambridge Monographs on Mechanics, Cambridge University Press, 341 pp.

*Linear and Nonlinear Waves.*John Wiley and Sons, Inc., 636 pp.

## Footnotes

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

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