An application of WKBJ theory for triad interactions of internal gravity waves in varying background flows

Motivated by the question of whether and how wave–wave interactions should be implemented into atmospheric gravity‐wave parametrizations, the modulation of triadic gravity‐wave interactions by a slowly varying and vertically sheared mean flow is considered for a non‐rotating Boussinesq fluid with constant stratification. An analysis using a multiple‐scale WKBJ (Wentzel–Kramers–Brillouin–Jeffreys) expansion identifies two distinct scaling regimes, a linear off‐resonance regime, and a nonlinear near‐resonance regime. Simplifying the near‐resonance interaction equations allows for the construction of a parametrization for the triadic energy exchange which has been implemented into a one‐dimensional WKBJ ray‐tracing code. Theory and numerical implementation are validated for test cases where two wave trains generate a third wave train while spectrally passing through resonance. In various settings, of interacting vertical wavenumbers, mean‐flow shear, and initial wave amplitudes, the WKBJ simulations are generally in good agreement with wave‐resolving simulations. Both stronger mean‐flow shear and smaller wave amplitudes suppress the energy exchange among a resonantly interacting triad. Experiments with mean‐flow shear as strong as in the vicinity of atmospheric jets suggest that internal gravity‐wave dynamics are dominated in such regions by wave modulation. However, triadic gravity‐wave interactions are likely to be relevant in weakly sheared regions of the atmosphere.


INTRODUCTION
spectrum is influenced to a considerable degree by modulation by a spatially and time-dependent resolved flow (Bretherton, 1966;Eckermann and Marks, 1997;Senf and Achatz, 2011). Especially at large vertical wavenumbers, the observed GW spectrum exhibits a slope, somewhat independent of time and location (Dewan and Good, 1986;Smith et al., 1987;Fritts and Vanzandt, 1993), which is reminiscent of the quasi-universal spectrum GWs are often thought to exhibit in the ocean (Garrett and Munk, 1972;1975;Polzin and Lvov, 2011). The universality of that spectrum is considered an indication of a transfer of energy in wavenumber (e.g., Olbers and Eden, 2013), usually attributed to nonlinear wave-wave interactions (Olbers, 1976;McComas and Bretherton, 1977;Pomphrey et al., 1980;Mueller et al., 1986;Lvov and Tabak, 2001;Lvov et al., 2004). Wave-turbulence theory (Hasselmann, 1962;Caillol and Zeitlin, 2000;Nazarenko, 2011;Eden et al., 2019) is a well-established tool for studies of corresponding spectra, considering statistical ensembles of GW fields, which most often focus on resonant triad interactions. In all of these the influence of mean-flow shear and varying stratification are neglected.
A complementary approach is Wentzel-Kramers-Brillouin-Jeffreys (WKBJ) theory (Bretherton, 1966;Grimshaw, 1975;Achatz et al., 2010;Achatz et al., 2017) which, instead of considering continuous wave spectra, describes the development of locally monochromatic GW fields which feature a nearly discrete spectrum. Moreover, the WKBJ approach takes into account nonlinear interactions between GWs and a spatially and temporally varying mean flow. WKBJ theory is the basis of present-day GW parametrizations, however for most applications relying on a steady-state approximation where GWs instantaneously assume an equilibrium distribution defined by the available sources and the mean flow (Kim et al., 2003;Quinn et al., 2020). While the GWs are modulated by the mean flow in this approximation, a GW impact on the mean flow is only possible once GWs dissipate, for example, by wave breaking. Non-dissipative direct GW-mean-flow interactions, relying on explicit GW transience, can only be described once the steady-state approximation is dropped. In numerical implementations this tends to lead to instabilities due to caustics (e.g., Rieper et al., 2013a) that can however be avoided when WKBJ theory is translated into a spectral formulation (Muraschko et al., 2015). Using this approach, Bölöni et al. (2016) have shown that direct, non-dissiative GW-mean-flow interactions dominate over dissipative effects in the dynamics of upward propagating GW packets and the wind induced by them. Hence it seems appropriate to generalize GW parametrizations accordingly.
Many process studies have investigated GW-GW interactions in the atmosphere (e.g., Dong and Yeh, 1988;1991;Fritts et al., 1992;Yi and Xiao, 1997;Huang et al., 2007). However, these have not alleviated the obvious deficiency of WKBJ-based GW parametrizations that they do not take such interactions into account (Kim et al., 2003). Shear effects are not of leading-order importance in the ocean (e.g., Garrett and Munk, 1972;1975;Mueller, 1976;Elipot et al., 2010), so there it seems appropriate just to supplement the spectral wave-action equation resulting from WKBJ by nonlinear scattering integrals as derived from a wave-turbulence theory (e.g., Olbers and Eden, 2013) assuming a zero or constant large-scale flow. In the atmosphere, however, it appears that the modulation of GWs by the large-scale flow is the dominant effect, so that a consistent numerical treatment of GW propagation through a sheared environment, while simultaneously undergoing wave-wave interactions, seems to be more important. Once a numerical implementation of a corresponding theory is available, one could better investigate the relevance of GW-GW interactions in the atmosphere as such. So far it seems to be unclear whether the typical lifetime of an atmospheric GW, between emission from its source and its turbulent breaking, gives nonlinear triad interactions enough room to act. If so, consecutive wave-wave interactions, that is, wave turbulence, could be an efficient mechanism for the nonlinear dissipation itself. Furthermore, an interesting question in this context is how much triad interactions are affected by wave modulation due to varying large-scale flows. Such modulation changes GW wavenumber and frequency so that a triad might be brought into and out of resonance. Hence strongly sheared environments might actually suppress nonlinear interactions, while such interactions, as described by wave-turbulence theory, might be more effective in less-sheared locations of the atmosphere.
With this motivation in mind, the work reported here builds on the study of Grimshaw (1988), who proposed a WKBJ theory for wave-wave interactions modulated by a slowly varying background flow. He considered the effect of the mean-flow shear on near-resonant triad interactions of internal GWs and outlined a possible approach to computing asymptotically the energy exchanges among the members of a triad that passes through resonance. The focus of the present study is the first (to our knowledge) implementation of such local resonant triad interactions into a numerical WKBJ model. In particular, we revisit the theory introduced by Grimshaw (1988) (Sections 2 to 4), simplify the equations to a quasi-one-dimensional setting (Section 5), and propose an interaction parametrization which allows for a straightforward application of the local interaction equations to the WKBJ modulation equations (Section 6). As an efficient tool for modelling the WKBJ equation system we use the spectral ray-tracing algorithm introduced by Muraschko et al. (2015) and expand it by a triad-interaction module. The resulting model is verified by constructing test cases of two interacting wave trains that generate a third wave train in the presence of a shear flow, and comparing the WKBJ simulations against wave-resolving simulations (Section 7). In general, wave modulation by a variable background stratification or a sheared mean flow are equally important in the atmosphere (cf. Achatz et al., 2017). However, we restrict the analysis to the case of Boussinesq dynamics with a constant background stratification and zero rotation for the sake of simplicity.

FLOW REGIMES, NON-DIMENSIONALIZATION AND SCALING ASSUMPTIONS
We consider the non-rotating inviscid Boussinesq equations, where v = (u, v, w) T , p, b, and N denote the velocity vector, the pressure, the buoyancy, and the buoyancy frequency associated with the background stratification, respectively. e z is the vertical unit vector. Note that we have scaled the pressure with the reference density so that it does not appear in the equations. For convenience we denote the horizontal velocity vector as u = (u, v, 0) T . The material derivative, D t , is defined by D t = t + v ⋅ ∇. We non-dimensionalize the governing equations with the help of the scaling parameters summarized in Table 1 Note that, in this scaling regime, rotation is a higher-order effect and is set to zero for simplicity. Thus, the non-dimensionalized governing equations for non-hydrostatic internal gravity waves in non-rotating Boussinesq dynamics are given by Here, the hatted variables denote the non-dimensional variables. Unless indicated otherwise, we will consider the non-dimensional variables without explicitly denoting the hat in the course of this study. Moreover we introduce a small parameter to scale the wave modulation and strength of the nonlinearities. To establish a consistent balance between modulation and nonlinearity, we follow a WKBJ approach with weak wave amplitudes of order O( ) and thus seek solutions of the form with y representing any of the fields v, p or b. The first term above is an expansion of the large-scale flow in terms of the scale-separation and wave-amplitude parameter, , while the second constitutes the wave field. The compressed coordinates, (T n , X n ), are defined by (T n , X n ) = ( n t, n x). In doing so we introduce a three-scale system where the fast scales, (T 0 , X 0 ), correspond to the wave oscillations and the slow scales, (T 2 , X 2 ), correspond to the slow variation of the mean flow which in turn causes a slow modulation of the wave fields. The intermediate scales, (T 1 , X 1 ), as explained below, are associated with the nonlinear wave-wave interaction. Choosing a wave field with leading order O( ), we balance the strength of the nonlinear terms with the modulation (Grimshaw, 1988;Glebov et al., 2005). The summation over the index represents the superposition of several wave trains in the solution. Moreover, for each wave train, we define the wave frequency, , and wave vector, k = (k , l , m ), as compressed temporal and spatial derivatives of the wave phase, , so that where the subscript indicates the scale of the derivative, that is, ∇ 2 = ( X 2 , Y 2 , Z 2 ) T . We hence construct wave solutions with slowly varying amplitudes, wavenumbers and frequencies on the compressed scales (T 2 , X 2 ). It should be noted that, for a superposition of wave trains with slowly varying amplitudes, the various harmonics may be separated and the equations may be written for the individual wave trains, only if the corresponding frequencies and wavenumbers are sufficiently separated. In particular the frequency difference of any two wave trains must be at least − ∼ O(1). A rigorous treatment may be done with the aid of the weak asymptotic method as introduced by Danilov (2001).
Following similar arguments, the quadratic nonlinear terms are only important where the conditions for a resonant triad are satisfied or nearly so. The behaviour is then analogous to the spectral passage through resonance of harmonic oscillators (Neu, 1983). In the case of an isolated triad, the quadratic nonlinearities scale with an exponential phase factor e iΔ ∕ 2 , where the phase difference is defined as Δ = ± ± − with signs depending on the various triad combinations (cf. Grimshaw, 1988). In the case of an exact and static (i.e., time-independent) resonance, one finds Δ ≡ 0 such that the phase factor becomes unity. In terms of wave vectors and frequencies that is, These are the the well-known resonance conditions of the classical interaction with constant stratification and zero background. However, if the resonance is not exactly satisfied or the phase difference is a function of time and space due to wave modulation, the phase factor enters the nonlinear interaction equations. For visualization of the local scaling, we locally expand the phase difference, Δ , in the compressed time. In particular one finds Thus the typical exponential term, e iΔΦ∕ 2 , due to quadratic nonlinearities is oscillating with the fast time-scale, T 0 , in general but becomes a function of the intermediate time-scale, T 1 , near resonance. The latter is the case as long as T 1 ∼ O(1) and hence Δ = Δ ∕ T 2 ∼ O( ). Consequently, the quadratic nonlinear terms are important only in an -neighbourhood around resonance (Grimshaw, 1988). A similar argument can be employed in all spatial dimensions such that one may obtain an analogous condition for the wavenumbers, Δk ∼ O( ).
We thus follow Grimshaw (1988) and Glebov et al. (2005), and consider two distinct regimes: the linear offresonance solution, where the nonlinear triad terms can be neglected, and the weakly nonlinear near-resonance solution.

THE LINEAR OFF-RESONANT SOLUTION
As long as the nonlinear terms do not come into play, the off-resonance solution is equivalent to the classical linear internal-gravity-wave theory, and all fields depend on the slow coordinates, (T 2 , X 2 ), only. Consequently the resulting equation hierarchy is equivalent to the well-known linear WKBJ theory for non-hydrostatic internal gravity waves (e.g., Achatz et al., 2010;Sutherland, 2010). Therefore we will only briefly review the most important results here, obtained after inserting (7) into Equations (4)-(6) and sorting in terms of powers of and the phase factor.

Leading-order mean flow evolution
In view of the assumptions of Boussinesq dynamics and weak wave amplitudes, the leading-order mean-flow velocity is purely horizontal and incompressible Furthermore it is governed by Hence the leading-order horizontal mean flow is independent of the wave field. Moreover, we find that B (0) 0 = B (1) 0 = 0, and that the leading-order mean-flow pressure, P (0) 0 , and the leading-order buoyancy, B (2) 0 , are in hydrostatic balance. We also obtain W (0) 0 = 0, and with this the evolution of the leading-order buoyancy is given by and therefore is linked to the leading-order vertical mean wind, W (4) 0 .

Dispersion and polarization relations
The internal gravity wave evolution is characterized by the following dispersion relation (e.g., Sutherland, 2010) with the intrinsic frequency,̂= − k ⋅ U (0) 0 . The polarization relations are Note that we restrict our analysis to the internal gravity wave evolution and neglect the vortical mode corresponding to the solution̂= 0. The next-order wave equations reveal that the next-order mean-flow velocities vanish, that is, U (1) 0 ≡ 0.

The eikonal equations
We use the standard definition for the group velocities corresponding to the extrinsic and intrinsic frequencies where ∇ k denotes the derivatives with respect to the corresponding wavenumbers ∇ k = ( k , l , m ) T . Using the dispersion relation (Equation 15) one may derive the evolution of the frequencies and wavenumbers -the eikonal equations. Specifically, where the explicit form of the intrinsic group velocity,ĉ g, , and the extrinsic group velocity, c g, , are given bŷ

Wave action conservation
The linear wave action conservation in standard form is where the wave action,  , is defined as the ratio of the wave energy and corresponding intrinsic frequency,

Wave impact on the mean flow and leading-order vertical winds
Exploiting higher orders one may find that the second-order horizontal mean flow, U (2) 0 , the next-order buoyancy, B (3) 0 , and the leading-order vertical wind, W (4) 0 , are directly interacting with the wave field. Moreover, the leading-order vertical mean flow is connected to the leading-order buoyancy by Equation (14) and therefore affects the hydrostatic relation. However, the leading-order vertical mean flow is three orders smaller compared to the leading-order wave amplitude, W (1) . Also, there is no feedback onto the wave field. Thus the wave impact onto the mean flow will be treated as a higher-order effect and will not be taken into account here, in accordance with the weak-wave-amplitude assumption.

THE NONLINEAR NEAR-RESONANCE SOLUTION
The near-resonance solution hinges on the quadratic triad terms and depends on the intermediate-scale coordinates (T 1 , X 1 ). Hence the wave amplitudes and second-as well as higher-order mean flow also vary on the intermediate scales. The leading-order mean-flow contributions are assumed to depend on the slow coordinates (T 2 , X 2 ) only, as they correspond to the slowly varying background and wave modulation. We next derive the asymptotic hierarchy closely following Grimshaw (1988) and Achatz et al. (2010).

Leading-order mean-flow evolution
Similarly to the off-resonance solution, the horizontal mean flow and pressure are of order O(1), U (0) 0 = U (0) 0 (T 2 , X 2 ) and P (0) 0 = P (0) 0 (T 2 , X 2 ). They represent the slowly varying background state and are therefore assumed to be dependent on the slow scales only. Also the leading-order mean-flow buoyancy, B (2) 0 , and leading-order vertical wind, W (4) 0 , are of order O( 2 ) and O( 4 ), respectively. We thus set B (2) 0 = B (2) 0 (T 2 , X 2 ). The leading-order incompressibility criterion requires that Averaging Equation (23) over the intermediate scales, and requiring sub-linear growth of , we find that the leading and next-order mean-flow velocities, U (0) 0 and U (1) 0 , are incompressible: The evolution of the leading-order horizontal mean flow is given by Again one may average Equation (25) over the intermediate scales and obtain, via the sub-linear . One then finds that Equation (25) can be separated so that where Equation (28) represents the hydrostatic balance of the mean flow to leading order. Thus the evolution of the leading-order mean flow is equivalent to the linear regime (Equation 13). The leading-order mean-flow buoyancy evolves as Leading-order mean-flow buoyancy and vertical mean wind are therefore linked similarly to the linear regime. Even though we may neglect the small leading-order vertical mean flow, we will use this statement to obtain the formal leading-order matching conditions for the two regimes.

Dispersion, polarization, and interaction equations
Due to the small wave amplitudes, the dispersion relation as well as the polarization relations are retained from the linear evolution (cf. Equations 15 and 16). In contrast to the off-resonance solution, the wave amplitude and wave action equations comprise the nonlinear triad terms. Projecting the next-order wave evolution equations onto the normalized polarization relations (cf. Achatz et al., 2010), one arrives at the wave amplitude equation, with the interaction coefficients (cf. McEwan and Plumb, 1977) The wave action evolution then follows as with the wave action,  , being defined analogous to the linear solution (cf. Equation 22). Here, the interaction term, T (1) , is given by Thus in the near-resonance solution, wave action is conserved up to the exchange of energy between the modes. Note that, due to the cubic nonlinearities, that is T (1) ∼ W (1) W (1) W (1) * (and complex conjugates), Equation (33) is ill posed when the amplitude, W (1) , is initially zero (cf. Equation 22). Instead it is necessary to solve the amplitude equation (Equation 31) near resonance.

Energy conservation
Naturally the linear solution comprises the wave-action conservation (Equation 21). Near resonance, while the wave action of each wave train is not conserved, wave triads exchange energy such that the sum of all wave energies is conserved. Therefore one may assess the total energy balance by considering an individual triad with resonance conditions The evolution of the corresponding wave energies is then given by the three coupled equations When summing the contributions of all members of the triad one gets where A ∈ R is equal to Here we note again that all wavenumbers and intrinsic frequencies depend on the slow time and spatial scales only. Thus, in an -neighbourhood around exact resonance, the resonance conditions (Equations 35 and 36) remain valid. Applying these conditions to Equation (41) yields Thus the wave energy is conserved during the interaction of any triad. For any wave train that is not a member of a resonant triad, the interaction terms vanish due to the asymptotic scale separation. Hence we have conservation of the sum of the wave energies over all packets near resonance and write accordingly

Wave impact during interactions
While the evolution of U (1) 0 on the intermediate coordinates, (T 1 , X 1 ), is independent of the waves, it is influenced by the GW momentum flux convergences on the large-scale coordinates, (T 2 , X 2 ) (not shown). However, since the near-resonance solution is valid in an -neighbourhood in T 2 around the exact resonance, the slow change of the next-order mean flow is of the order O( 2 ) with respect to the leading order. Also the term describing the impact of U (1) 0 on the wave fields in (31), , may be interpreted as a next-order correction to the dispersion relation. We therefore neglect the next-order horizontal mean flow, U (1) 0 . Similarly, the leading-order vertical mean flow, W (4) 0 , is driven by GW fluxes. Moreover it is related to the leading-order buoyancy similar to the off-resonance solution. We thus conclude that the vertical mean flow is generally dependent on the waves near resonance. However, we neglect this effect since there is no feedback on the wave field and the largest non-zero vertical mean flow is three orders smaller than the assumed wave amplitude.

Matching the solution regimes
The prognostic and diagnostic equations for the regimes near and far from resonance were summarized above. Naturally we require that in the limit, → 0, the mean flow and wave amplitudes of the same hierarchy must match at the regime transition. To determine the regimes of validity of the solutions, one may consider the validity of the phase expansions. In particular, the near-resonance solution is valid in an -neighbourhood around the exact resonance on the slow scales, (T 2 , X 2 ) ∈ O( ). Hence, we seek conditions so that the off-resonance solution matches the near-resonant solution in the limit as the resonance manifold is approached. This corresponds to the limit (T 1 , X 1 ) → ±∞ (cf. Glebov et al., 2005).

Mean flow
The leading horizontal mean flow in both solutions are non-divergent (Equations 12 and 23) and hydrostatic (Equations 13 and 28). Moreover, the leading-order buoyancy is independent of the wave field in the linear regime (Equation 14). The evolution in the near-resonance solution, averaged over (T 1 , X 1 ), is given by Equation (30). Consequently the only matching condition for the leading-order mean flow is given by

Wave amplitudes
By assumption, the wave properties, i.e., and k , depend on the large-scale coordinates (T 2 , X 2 ) and consequently obey the eikonal equations (Equations 18 and 19) in both solutions. Moreover, the leading-order waves follow the dispersion and polarization relations (Equations 15 and 16). While wave action is conserved off-resonance (Equation 21), it is subject to nonlinear exchange on the intermediate coordinates in the near-resonance solution (Equation 33). To find the matching conditions between the two solutions, we thus seek the limit of the near-resonance solution for (T 1 , X 1 ) → ±∞. First we note that the interaction term scales with exponential functions of phase differences: By assumption, one has − T 2 Δ = Δ ∼ O(1) and ∇ 2 Δ = Δk ∼ O(1) in the limit (T 1 , X 1 ) → ±∞. Thus the nonlinear forcing term can be expanded to leading order This term is dependent on the short-scale coordinates (T 0 , X 0 ), and must therefore vanish after averaging over the large scales (cf. Danilov, 2001). Thus, in the limit (T 1 , X 1 ) → ±∞, the wave action equation becomes a conservation law similar to the linear solution (Equation 21): We conclude that in the limit (T 1 , X 1 ) → ±∞ the wave amplitudes are not driven by interaction on the intermediate scales. The formal matching condition is given by

SUMMARY OF DIMENSIONAL EQUATIONS IN 1.5D
For the application of the above-derived equations, we revert to dimensional variables. Moreover we assume homogeneity of mean flow and wave amplitudes in the horizontal direction such that the equations become effectively one-dimensional (cf. Muraschko et al., 2015). Finally we also assume that the horizontal mean-flow velocity, U (0) 0 , as well as the horizontal wave vectors, k h , have an x-component only. Under these assumptions the incompressible, constant, and hydrostatic mean flow satisfies in both regimes The waves are governed by the eikonal equations (Equations 18 and 19), the dispersion and polarization relations (Equations 15 and 16), and the wave action or wave amplitude equations (Equations 21,31, and 33). The dimensional eikonal equations are with the dispersion relation and the group velocitieŝ The polarization relations are While wave action conservation holds for the linear off-resonant solution, that is, the near-resonance regime requires additional interaction terms such that where the interaction coefficients are given by However, Equation (57) is ill posed when the wave amplitude, w , is zero at an initial time. Instead we solve the complex wave amplitude equation given by where the second-order horizontal mean flow is neglected. The evolution of the phase functions, , along the wave characteristics is given by the definition of the wavenumber and frequency where we have rescaled the phase function such that = −2 . Hence the small parameter does not appear explicitly in the equations.
Equations (51), (52), and (56) are equivalent to Grimshaw's modulation equations (Grimshaw, 1977;1988) for weakly nonlinear non-hydrostatic internal gravity waves. Here, Equation (57) replaces Equation (56) where near-resonant triad interactions are relevant and the nonlinearities come into play. This system may be employed numerically to estimate wave-wave interactions in the context of WKBJ ray-tracing simulations as discussed below.

A SEMI-EMPIRICAL PARAMETRIZATION FOR THE INTERACTION EQUATIONS
In the previous sections we have presented a weakly nonlinear multi-wave WKBJ theory based on the non-hydrostatic Boussinesq equations. The resulting modulation equations are summarized in Equations (49)-(59) assuming horizontal homogeneity so that they are effectively one-dimensional. These equations may be solved numerically using several approaches.

6.1
Phase expansion around resonance Following Grimshaw (1988), near the manifold of exact resonance one may expand the phase functions, , to second order in time and space and project the resulting interaction equations onto a space-time direction which is perpendicular to the resonance manifold in z-t space. In the limit → 0, the exchange of energy among the members of a triad implied by the near-resonance solution then appears as a "jump" across the resonance manifold. However, in implementing this approach, we encountered certain difficulties. In particular, the projection onto the cross-resonance coordinate leads to singularities and secular growth in the interaction equations where the space-time trajectory of any triad member is parallel to the resonance manifold. Also, singularities may occur where the second-order truncation of the phase expansion becomes invalid and the equations have to be rescaled. Since both these issues do arise at rather common settings of wavenumbers and background shear strengths, we do not follow this approach.

Equivalent window method
We observe that the exponential term in the interaction equations (Equation 59) acts as an integration window limiting the interaction, depending on the spectral deviation from resonance. We thus suggest finding a spectral window function with an equivalent width.
In such a case, the interaction equations may be solved as if in exact resonance but limited in terms of spectral deviation from resonance. This approach is explained below.
The asymptotic theory presented earlier comprises two scaling regimes. While the off-resonance solution follows linear dynamics on slow time-and spatial scales, with corresponding coordinates, (T 2 , X 2 ), the near-resonance solution is characterized by the interaction of GW triads on intermediate time-and spatial scales, with a dependence on (T 1 , X 1 ). In both cases the background is assumed to vary on the slow scales only. Thus near resonance and in the asymptotic limit → 0, the characteristic length-scales of both the wave train amplitudes and the background shear are virtually infinite with respect to the interaction (i.e., intermediate) scales. Motivated by this asymptotic limit, we consider GWs in a constant background shear z u 0 ≠ 0 with infinite extent in the vertical. Similarly, the slowly varying wavenumbers, m , are assumed to be homogeneous in the vertical such that z m = 0.
Consequently, the local tendency of the wave frequencies can be expressed as where we have used the eikonal equations (51)-(52), the fact that the wave frequencies and wavenumbers are related by − z = z t = t m , as well as the above assumption z m = 0. Expanding the phase difference locally in time, one finds where the linear term vanishes when expanding around exact resonance. Moreover, in resonance one has (Δk) 0 = 0. Also, without loss of generality, we set (Δ ) 0 ≡ 0. Finally the phase difference becomes approximately For an explicit triad, the interaction equations (Equation 59) thus are where the phase difference is Δ = 2 + 3 − 1 and the interaction coefficients are given by Inserting the local phase evolution (Equation 64), one finds that the right-hand sides in Equations (65)-(67) are independent of the vertical coordinate, z. Thus the homogeneity assumption can be repeated for the wave amplitudes, w . Consequently the vertical gradients vanish at any time and height, z w ≡ 0. Hence the system Equations (65)-(67) simplifies to where the dephasing, ( t Δ̂) 0 ∕2, can be expressed in terms of the wavenumbers using Equation (62): This system of equations is equivalent to the evolution of plane waves in a background shear flow with infinite extent. This image may be useful to understand the asymptotic, i.e., local, passage through resonance.
For comparison, we set up a simplified system making use of the fact that, in a small neighbourhood around exact resonance, the resonance conditions are satisfied approximately. To balance the limited width of validity of approximately exact resonance, one may introduce a window function G(t − t 0 ) of the form where represents the Heaviside function. This represents a symmetric box with value G = 1 around the resonance time t 0 , and value G = 0 elsewhere. A simplified system then reads To evaluate G(t − t 0 ), one may compare Equations (76)-(78) to Equations (71)-(73) with a quadratic dephasing near-resonance as described above. Integrating the exponential term over the time corresponding to a local asymptotic expansion, i.e., from −∞ to ∞, yields the Fresnel integral Equating the absolute value of the real part of I F with the integral over the window, G(t − t 0 ), yields giving a first-order estimate for the window half-width, t † . However, in order to practically apply this estimate, one would need to solve for the exact resonance manifold and consecutively reconstruct the effective interaction region before integrating the model. We therefore evaluate t † in terms of an effective spectral resonance deviation which can be readily diagnosed during run time of the model integration. In particular, we define a normalized resonance function R such that Thus the resonance function at the boundaries of the interaction window G(t − t 0 ) is given by Equation (82) directly gives a leading-order estimate for the spectral resonance width. This estimate covers the dependencies on the wavenumbers and background shear but may need tuning for a global parametrization. We therefore introduce a tuning parameter, , and set the spectral window where denotes the Heaviside step function and R † = |R(t 0 ± t † )|. The parameter may then be optimized for best agreement between simulation results of Equations (71)-(73) and Equations (76)-(78). This formulation now allows for simulations that locally diagnose resonance deviations and enable interactions where necessary without solving for exact resonance manifolds.
Note that the described procedure approximates the interaction equations (65)-(67) such that the total changes of the magnitude of the complex wave amplitudes, w , are recovered. However, the parametrization introduces a modification of the phases of the complex wave amplitudes which are then bound to differ from the full solution. Analyzing interaction systems with constant phase difference, Bustamante and Kartashova (2009) found that this may modify the evolution of the wave amplitudes in terms of both energy exchange rates and maximum energy exchange. While the effective resonance interval, R † , is chosen such that the total energy exchange is recovered, deviations in the exact evolution between the full and the parametrized solution may occur.

Estimating the tuning parameter
In order to estimate the tuning parameter, , we compare numerical results of the reference system given by Equations (71)-(73) and the system Equations (76)-(78) with (83) for a range of specific resonances. With a length-scaleL = 5 km∕2 , we consider triads with non-dimensional horizontal wavenumbersk 1 = 0.2 and k 2 =k 3 = 0.1. Whilem 1 andm 3 follow from the resonance conditionsm 2 +m 3 −m 1 = 0 and̂2 +̂3 −̂1 = 0, the non-dimensional vertical wavenumber of the second GW is chosen from the intervalm 2 ∈ [0.56, 10]. Here the lower boundary corresponds to the lowest wavenumber at which Equations (71)-(73) can be solved while avoiding passage through a second distinct resonance. The upper boundary is chosen in accordance with the asymptotic assumption m 2 = O(1). We thus consider dimensionful horizontal wavelengths of the order ∼10 km and vertical wavelengths of the order ∼1 km. In general, atmospheric GWs cover a broad range of spatial scales (Callies et al., 2014). Here we choose rather small wavelengths for GWs in order to be consistent with the simplification of non-rotating Boussinesq dynamics. A critical reader will notice that most of the cases considered imply |m ∕k | ≫ 1, not quite consistent with the original non-dimensionalization of the Boussinesq equations, where equal horizontal and vertical length-scales have been assumed. However, this is only an apparent violation of the assumptions. Redoing the scale-asymptotic analysis with anisotropic scaling leads to a limiting form of the presented formulae with k 2 << m 2 in Equations (53) and (58). All other results stay intact. However, while our formulation is quite general, we present results for relatively anisotropic test cases so that the interaction coefficients (Equation 58) are larger and strong energy exchanges may be observed.
As for wave and shear amplitudes, in accordance with the weakly nonlinear theory, the initial amplitude of the wave trains relative to the corresponding static instability criterion are varied between 10 −3 and 10 −1 . The chosen background shear strength is equal to a value corresponding to the maximum shear in our reference simulations which are introduced in the later part of this study. In particular z u = 2 ∕40,000 s −1 ≈ 1.6 × 10 −4 s −1 . This value corresponds to relatively weak sheared background winds in the atmosphere where jet strengths U = O(10 m⋅s −1 ) and a tropopause height H = O(10 km) imply z u = O(10 −3 s −1 ). However, as shown in the discussion, there are large areas, for instance in the midlatitudes and polar regions in spring and autumn after the breakdown of the polar night jet, that are well represented by a shear strength z u = O(10 −4 s −1 ). Moreover, it allows for a strong wave modulation (order O(1)) on time-scales approximately two orders of magnitude longer than a typical wave period (∼10 3 s). The chosen background shear is therefore consistent with both the scaling assumptions of the asymptotic theory as well as observed atmospheric conditions. Stronger shears will be discussed below as well, and it will be shown there that the associated wave modulation by the mean flow partially suppresses the nonlinear interactions.
Optimal values for the parameter are then found for 189 central wavenumbers m 2 and 21 different amplitudes within the given intervals. In particular we use a Nelder-Mead procedure to find the least-square deviation between the the two model results (Nelder and Mead, 1965). We find that the optimal value for is approximately constant for all central wavenumbers and in the limit of small amplitudes (Figure 1). For amplitudes near 10 −1 the optimal tuning parameter decreases with minimum values as low as 0.04. In this regime, the characteristic time-scales of the exactly resonant system are comparable to the time-scales given by the dephasing (cf. Equation 74). This causes a potential systematic bias to this method at large amplitudes. The median of all optimal values for 189×21 evaluations, spanning the intervalsm 2 ∈ [0.56, 10] and ∈ [10 −3 , 10 −1 ], is equal to = 0.9969 ≈ 1. The approximately constant shows that the form of the effective spectral interaction threshold, R † (Equation 82), covers the dependency on the wavenumbers with high accuracy. This strongly suggests that a globally constant may describe the spectral interaction threshold across a wide range of wavenumber scales. Moreover the constant optimal tuning parameter, ≈ 1, for small amplitudes suggests that the derived effective spectral interaction threshold, R † , in combination with a global tuning parameter, , are appropriate to parametrize the spectral passage through a resonance across all scales covered by the asymptotic theory. However, for amplitudes approaching the limit of static instability at O(1), the parametrization, understandably, may not be as accurate.
For a qualitative error estimate, we visualize the difference between numerical solutions to the simplified system (Equations 71-73) and the parametrized system (Equations 76-78) in Figure 2. Here we distinguish wave amplitudes solving the simplified system, w (s) j , and wave amplitudes solving the parametrized system, w (p) j . The parameters for the shown simulation are z u = 2 ∕40,000 s −1 ,k 1 = 0.2,k 2 =k 3 = 0.1,m 2 = 5, = 1, and = 10 −2 . Considering the trajectories of the solutions on the complex plane (cf. Figure 2a-c), it is evident that the parametrized solution (blue) cannot reproduce the evolution of the phases of the complex wave amplitudes of the simplified solution (red). Despite this obvious shortcoming, the parametrization predicts the total change of the absolute value of the wave amplitudes with small errors (cf. Figure 2d-f).
To highlight the quantitative error, we compute solutions to both systems for the previously used parameter space and evaluate the ratios of the resulting energies after the interaction. In particular we choose relative wave amplitudes ∈ [10 −3 , 10 −1 ], central wavenumberŝ m 2 ∈ [0.56, 10], a background shear z u = 2 ∕40,000 s −1 , and a parametrization constant ≡ 1. As a diagnostic we j and w (s) j , relative to their initial value. Each column shows one of the three triad members calculate the ratio of the wave energies of the two solutions after the interaction for each triad member, E We find that the deviation is generally smaller than 8% with a mean and median value <2% for all three triad members over the whole tested parameter space (Figure 3). Systematic biases arise at initial relative amplitudes > 0.04. Here the energy transferred to the generated triad member may be either under-or overestimated ( Figure 3a). Strongest deviations occur near = 0.1 where the parametrization leads to overestimates of the transferred energy while overestimating and underestimating the energy of either one of the generating triad members (Figure 3a-c).

The test case definition
As a generic test case, we use the generation of a third wave train through the resonant interaction of two given wave trains. A predefined sinusoidal background shear modulates the given wave trains to enable spectral passage through resonance. In particular we use a periodic domain of height H = 40 km (so that H∕L = 16 ) and a background shear flow defined as u = u 0 sin(2 z∕H). In the tests below, u 0 ∈ [0, 1 m s −1 ]. The maximum shear is therefore approximately z u 0 ≈ 1.6 × 10 −4 s −1 . Based on this mean-flow profile, we derive initial conditions for the wave trains such that the modulated wave-wave interaction takes place near the centre of the domain. Note that the vertical scale H is much larger than the considered vertical wavelengths. Thus this structure of the mean flow is compliant with both the periodic boundary conditions of the simulations and the slow modulation assumption. As initial conditions for the given wave trains we consider Gaussian wave packets with vertical velocities The centre of the wave trains, z j , and initial wavenumbers, m j , are chosen such that the wave trains are in resonance and overlap in the centre of the domain after approximately half the integration time, as described below. The latter is chosen such that at the final time the energy exchange is negligible. The initial wave amplitudes are taken to be a given fraction, , with respect to the static F I G U R E 3 Ratios of wave energies corresponding of the solutions the simplified system (Equations 71-73) and the parametrized system (76-78) after the interaction. Here, the background shear is set to z u = 2 ∕40,000 s −1 , and the parametrization constant is ≡ 1. (a), (b), and (c) are associated with j = 1, j = 2, and j = 3, respectively instability criterion (e.g., Achatz et al., 2017). In particular All other initial fields are determined through the polarization relations (Equation 55). Wave amplitudes considered below are in the range ∈ [10 −2 , 10 −0.4 ]. The initial packet width is constant with = 2 km, the horizontal wavelengths are set to 2 = 3 = 50 km such that 1 = 25 km. These wavelengths correspond to non-dimensional wavenumbersk 1 = 0.2 andk 2 =k 3 = 0.1. They are chosen such that the interaction coefficients (Equation 32) are large enough for the resonance conditions to permit a wide range of wavenumbersm 2 for exact resonances. Validations against cases with larger horizontal wavenumbers have been done as well (not shown), with qualitatively similar results to those reported here.

WKBJ validation against wave resolving simulations
For qualitative and quantitative comparisons, we run wave-resolving simulations as well as WKBJ ray-tracing simulations with equivalent initial conditions. In particular we employ, in Boussinesq mode, the code PincFloit with a second-order MUSCL (Monotonic Upstream-centred Scheme for Conservation Laws) scheme utilizing an MC (monotonized central) flux limiter (Rieper et al., 2013b;Wilhelm et al., 2018). The time integration is realized through a third-order Runge-Kutta scheme. To circumvent numerical attenuation, it is run with a resolution that permits for at least 12 points per wavelength in the initial conditions. For comparison we use the spectral WKBJ ray-tracing code Note: The corresponding horizontal wavenumbers are (k 1 ,k 2 ,k 3 ) = (0.2, 0.1, 0.1). By convention, the first wavenumber corresponds to the generated wave train such thatk 1 =k 2 +k 3 andm 1 =m 2 +m 3 introduced by Muraschko et al. (2015), augmented by an interaction module corresponding to the parametrized solution method presented in Section 6. A brief technical description of the implementation can be found in the Appendix. The wave-resolving simulations are carried out in two dimensions with a periodic horizontal domain. In particular we set the domain such that it has a horizontal extent equal to a multiple of the horizontal wavelengths of the wave, and we employ periodic boundary conditions. We choose three distinct resonances with different resonant wavenumber triads around which we construct the simulation. In particular the central triads are characterized by the non-dimensional vertical wavenumbers,m , as summarized in Table 2. To define appropriate initial conditions for the test cases, the ray tracer was employed with the interaction scheme disabled and using a negative time step, thus integrating backward in time. The initial conditions were chosen to be Gaussian wave packets in exact resonance, and the reverse integration time was set to half the desired model integration time for the corresponding test case. The maximum amplitude position and mean wavenumber of the resulting wave trains were used to set z and m in the initial conditions for the test cases (cf. Equation 84). The initial amplitudes are set to = 0.1 where the background shear is varied and the background velocity is set to u 0 = 1 m⋅s −1 where the amplitudes are varied. In general the amplitudes do not only vary during triadic interactions but also elsewhere due to the wave modulation and the wave action conservation (Equations 21 and 22). Thus the initial amplitudes are not equivalent to the amplitude at the transition to the near-resonance regime. However, in cases of small background velocities, u 0 ≤ 1 m⋅s −1 , the effect is small and the initial amplitude remains a good estimate as to how strong the waves are relative to static instability near resonance. For stronger background velocities, where the modulation effect dominates the experiments, the amplitude modulation influences the strength of the interactions. It is therefore balanced by modifying the initial amplitude such that the amplitudes are comparable in the interaction regime.
For the analysis and visualization, the fields from the wave-resolving simulations are first Fourier-transformed in both spatial directions, then separated by wavenumber and projected onto the polarization relations (Equation 55). After this filtering of each of the three contributing waves, the corresponding fields are used to determine the wave energy densities of the individual wave trains in physical space. This procedure is similar to the one used by Borchert et al. (2014), however without a separate local Fourier transform around each grid cell. To estimate energy exchanges, the energy densities were integrated in the vertical and compared among the individual triad members. For convenience we normalize the vertically integrated wave energy densities by the sum of all triad components.

Energetics of the interacting wave trains
An example of the total wave energy density corresponding to the case withm 2 = 10 and u 0 = 1 m⋅s −1 is shown in Figure 4. Note that the energy of the background flow is filtered and not shown due to the projection onto the wave modes. Naturally, the WKBJ simulations do not account for the variation on the scale of the wavelengths and lack structure with respect to the wave-resolving simulations (Figure 4b). As a result, interference patterns with high wave energy densities are not present in the WKBJ simulations. However, the evolution of the wave train amplitudes are reproduced both qualitatively and quantitatively.
Separating the fields into the distinct wave trains and integrating the wave energy densities in the vertical as explained above, we find a generally good agreement in the temporal evolution of the individual wave energy densities ( Figure 5). For early times, t ≤ 24 hr, and late times, t ≥ 48 hr, there is approximately no wave-wave

F I G U R E 5
Integrated wave energy densities of the separated wave trains in the simulations withm 2 = 10 and u 0 = 1 m⋅s −1 , showing (a) the absolute wave energy density, (b) the relative wave energy density, and (c) the mean absolute wavenumber, with wave-resolving (solid) and WKBJ simulation (dashed) results. The three colours represent the two initial wave trains (blue and yellow) as well as the generated wave train with zero inital energy (turquoise). The black lines in (a) depict the sum of the individual wave energy densities. The solid lines are broken where the wavelengths of the different wave trains have similar absolute values and therefore cannot be separated interaction and the evolution of the wave trains is dominated by the wave modulation through the background shear flow (Figure 5a). Naturally the weakly nonlinear multi-scale WKBJ theory is an approximation to the fully nonlinear dynamics. As an example, the waves in the wave-resolving simulations may be modulated not only due to the prescribed mean-flow shear but also due to the self-induced mean flow (Sutherland, 2006a). However, this effect is neglected in the WKBJ simulations, as it appears in the theory only at higher orders in the scale separation parameter, . Moreover the scale separation assumption where → 0 and the weak wave amplitude assumption (Equation 7) are potential sources for errors. These systematic biases or other higher-order effects associated with wave modulation and not accounted for cause a mismatch in the evolution of the wave energy densities (Figure 5a), particularly for t ≥ 48 hr. To account for the effect of the modulation on the wave energy, we normalize the individual wave energy densities by the total wave energy density, that is, the sum of the individual wave energy densities (Figure 5b). These relative wave energy densities show a qualitatively good agreement throughout the whole simulation -including the dynamics during the interaction at times from approximately 24 hr to approximately 48 hr. However, a closer look reveals that the wave energy of the generated wave differs by up to 30 % relative to the LES (Figure 5a,b). Possible reasons for this mismatch could be a systematic bias in the phases of the complex wave amplitudes introduced by the parametrization (cf. Section 6), inaccuracies of the parametrization parameter, , near = 0.1 (cf. Figure 1), or higher-order modulation effects as discussed above. Despite this shortcoming, the total energy exchange is well reproduced for various settings as discussed below. If, for example, one would consider an experiment with the same initial conditions but neglected wave-wave interaction, the superimposing wave trains would be modulated with constant wave action densities and not exchange any energy. In contrast, an experiment where the wave-mean-flow interaction is negligible, but the near-resonant wave-wave interaction is included, would exhibit constant wave energies, but also no energy exchange, as it is the modulation which brings the waves into resonance. Note that where the absolute value of the wavelengths of the wave trains are too close (cf. Figure 5c), a separation of the wave trains is numerically difficult and therefore omitted (cf. gaps in Figure 5a,b).
Values of relative wave energy densities at the final time of the simulation serve as benchmark for comparisons in the further analysis of interaction simulations under varying conditions. In particular we consider varying background velocities as well as wave amplitudes (Figures 6-8).

The effect of the wave amplitudes
Varying the initial GW amplitudes, , i.e., varying the strength of the nonlinearities relative to the modulation F I G U R E 6 Integrated relative wave energies per wave train at the end time of the simulations with non-dimensional vertical wavenumberm 2 = 5 and (a) varying background velocity and (b) varying amplitude. While solid lines depict the wave-resolving simulations, dashed lines represent the corresponding WKBJ simulations. The three colours and corresponding symbols represent the two initial wave trains (blue/crosses and yellow/triangles) as well as the generated wave train (turquoise/x-symbols). The black line (dots) represents the total relative wave energy F I G U R E 7 As Figure 6, but for simulations with the non-dimensional vertical wavenumberm 2 = 10 of the wave trains, we find good agreement for relative amplitudes smaller than = 10 −1 (Figures 6b, 7b, and 8b). As expected, the energy exchange increases with increasing amplitude as the nonlinearities are growing stronger. For amplitudes > 0.1, i.e., closer to the criterion of static instability, the WKBJ simulations fail to reproduce the wave-resolving simulations qualitatively and overestimate the triadic energy exchange. In general larger amplitudes may be subject to stronger nonlinear effects like self-acceleration, modulational instabilities, or F I G U R E 8 As Figure 6, but for simulations with the non-dimensional vertical wavenumberm 2 = 15 overturning and turbulence (Sutherland, 2006b;Dosser and Sutherland, 2011;Bölöni et al., 2016). The wave-mean flow interactions, including the self-acceleration, are higher-order effects in the WKBJ expansion and may be taken into account to improve the method presented here in the limit of larger amplitudes. As expected a WKBJ theory may not cover the breakdown of the wave trains at amplitudes near the static instability criterion unless using a suitable parametrization (Lindzen, 1981;Bölöni et al., 2016). Also, the employed parametrization with a globally constant tuning parameter, , may lead to systematic biases at larger amplitudes (cf. Section 6).
In general, wave amplitudes are modified through both triadic interactions and wave modulation (cf. Equations 22,51,and 56). For small background flows u 0 ≤ 1 m⋅s −1 , the wave train evolution is dominated by triadic wave interactions and the amplitude variation due to wave modulation are comparatively small. Thus we do not correct the initial amplitudes with respect to the modulation while studying the influence of the wave amplitude and modulation on triadic interactions for u 0 ≤ 1 m⋅s −1 . The wave amplitudes in Figures 6-8 refer to the wave amplitudes in the initial conditions (cf. Equation 84). For runs with strong background flows, i.e., where the evolution is dominated by the wave modulation, we correct the initial amplitude such that the wave amplitudes are comparable near resonance.

The effect of the background shear strength
The wave modulation by the imposed background shear leads to a continuous spectral shift of the wavenumbers and frequencies. However, effective energy exchange is only possible near exact resonance. The resulting passage through the resonance conditions therefore leads to a more localized interaction and limits the energy exchange. A stronger background shear is associated with an increased modulation and thus generally leads to a reduction of the energy exchange between the triad members. This effect is well reproduced in the simulations presented here for m 2 = 5,m 2 = 10, andm 2 = 15 (Figures 6a, 7a, and 8a) with a mismatch for large wavenumbers (m 2 = 15) and small background velocities (u 0 ≤ 0.4 m⋅s −1 ) (Figure 8a).
For background velocity amplitudes smaller than u 0 = 0.5 m⋅s −1 , the wave-induced mean flow (e.g., Sutherland, 2006a) in the wave-resolving simulations may lead to background shear strengths comparable to the imposed shear. This shear consequently modulates the wave triads leading to a shift in wavenumbers and frequencies relative to the WKBJ simulation where the effect is not included, as explained above. This effect, albeit small at small amplitudes, may lead to significant differences between the wave-resolving simulations and the WKBJ simulation under conditions where the spectral width of the triad resonance is small. Correspondingly, the mismatch between the simulations for background flows, u 0 < 0.5 m⋅s −1 , at larger vertical wavenumbers,Lm 2 = 15, is believed to be caused by neglecting the wave-mean-flow interactions in the present WKBJ theory (cf. Figure 8a). In particular the self-induced wave modulation perturbs the near-resonant interaction such that, even at small amplitudes, the energy exchange is limited due to the frequency deviations. For comparison, we repeat the experiment with varying background flows form 2 = 15 but decrease the amplitude to = 0.05 (Figure 9). Naturally at smaller amplitudes not only the triadic wave interaction but also the self-induced wave modulation of the wave trains is reduced. Consequently the associated frequency deviation from exact resonance is smaller. At the same time we find that the triad interaction is significantly stronger in the wave-resolving simulation despite the reduced nonlinearities at small imposed mean flows, u 0 < 0.5 m⋅s −1 . This qualitatively changed behaviour of the wave-resolving simulation agrees well with the WKBJ prediction ( Figure 9).
Herein also lies a qualitative argument as to which nonlinear effects dominate the wave evolution. At amplitudes near the static instability threshold, all nonlinear effects, like the wave-wave interaction considered here, the wave modulation by the mean-flow shear, or the self acceleration excluded here, are predicted to be equally important (Achatz et al., 2017). However, reducing the amplitude changes the picture. Using a weakly nonlinear theory, we find a regime where the wave modulation F I G U R E 9 As Figure 8a, but for simulations with the non-dimensional vertical wavenumberm 2 = 15 and = 0.05 by the mean-flow shear and the wave-wave interactions dominate the dynamics while the self acceleration effect becomes a small correction (cf. Sections 2 and 5) which could be included in the theory by introducing an additional time and spatial scale (not shown). Also we have shown that at amplitudes ≤ 10 −1 the employed parametrization is valid, but may produce systematic biases at larger amplitudes (cf. Section 6). Additionally we find the qualitative importance of various effects to depend on the wave properties considered above. Consequently a quantitative mapping of the importance of the different nonlinear effects is dependent on many variables and is thus beyond the scope of this work.

7.6
Energy exchange at strong background shear flows As explained above, background flows with stronger shear may lead to an increase in wave modulation and consequently a decrease in triadic wave interactions. To include values with shear strengths typical for atmospheric jet regions we augment our findings with WKBJ simulations for mean flow strengths up to u 0 = 9 m⋅s −1 . Those simulations confirm that for strong shear the simulation is dominated by wave modulation with small energy exchange due to triadic interactions ( Figure 10). While the evolution of the wave action density shows virtually no variations (Figure 10a), the wave energy densities show strong variability due to the wave modulation (Figure 10b). To estimate the strength of the wave interaction, we repeat the simulations with a disabled interaction scheme and then compare the two simulations. This allows us to compute the fraction of energy transferred during the interaction. As a result we find that the transferred energy decreases systematically from values as high as 19 % to well below 5 % when increasing the background amplitude from u 0 = 1 m⋅s −1 to 9 m⋅s −1 (Figure 11).

F I G U R E 10
Wave action and energy densities of WKBJ simulations with u 0 = 9 m⋅s −1 as a function of time. While the evolution of the wave action densities in (a) are approximately constant, the wave energy densities in (b) are dominated by the wave modulation. The three colours correspond to the two initial wave trains (blue dots and yellow dashes) and the generated wave train (turquise dashes and dots). The solid black curve represents the total wave action in panel (a) and total wave energy density in (b)

F I G U R E 11
Fraction of wave energy density of the two initial wave trains (depicted in solid blue and dashed yellow) transferred to a newly generated wave train as a function of background flow amplitude, u 0 For a representative comparison with the atmosphere, we consider the zonally averaged zonal velocity from the UARS (Upper Atmosphere Research Satellite) Reference Atmosphere Project (URAP) climatology (Swinbank and Ortland, 2003) and its vertical shear. In particular during spring and autumn, after the breakdown of the polar night jet, we find large areas of relatively unsheared zonal velocity in midlatitude to sub-polar regions (Figure 12a,c). For comparison, we depict different vertical shear strengths as grey overlays that correspond to energy loss rates of >19% (black), 10-19% (grey), and 5-10% (light grey) in the idealized simulations (Figures 11 and 12b,d). These reveal that large areas correspond to background flow shear strengths that permit triadic wave-wave interactions without suppressing the spectral energy transfer due to wave modulation. We conclude that, depending on the region and the F I G U R E 12 Zonally and monthly averaged zonal velocities from the URAP dataset for (a) March and (c) September. For comparison, we overlay in (b, d) the vertical shear of the zonal velocity relative to the scaling of the idealized simulations, (H∕2 u 0 ) z u, with the reference values H = 40 km, and u 0 = 1 m⋅s −1 . The thresholds for the relative vertical shear correspond to energy transfer rates of 5% (5.9, hatched light grey), 10% (2.5, hatched grey), and 19% (1, hatched black) in the idealized simulations (cf. Figure 11) season, the GW dynamics in the atmosphere are likely to be impacted by triadic wave-wave interactions.

SUMMARY AND CONCLUSIONS
We have presented and applied a weakly nonlinear, Boussinesq theory of non-hydrostatic internal gravity waves (GWs) in a varying mean flow with constant stratification, extending previous work by Grimshaw (1988). The theory comprises a superposition of wave trains whose amplitudes are modulated by a slowly varying background. There are three well-separated scales: the GW period and wavelength define the fast and short scales, the spatial scale of the mean flow represents the longest scales, and nonlinear GW-GW interactions act on intermediate scales. Away from resonance, GWs follow linear WKBJ dynamics, characterized by short and long scales, and wave action is conserved. In resonance, GWs depend on all three scales, and energy is exchanged between GWs in triadic interactions. Wave amplitudes are weak so that GWs are still well-defined, including dispersion and polarization relations, and there is no leading-order impact on the mean flow. The modulation by the mean flow permanently changes the GW wavenumbers so that they are brought by this process into and out of resonance. This adds an additional source of spectral and spatial variability not accounted for in theories without a varying mean flow, like wave-turbulence approaches (e.g., Nazarenko, 2011). For a numerical implementation of this theory, we have supplemented a spectral ray-tracing code (Muraschko et al., 2015) by a wave-wave interaction module. Consistent with the two scaling regimes, a parametrization for an effective spectral resonance width has been developed, allowing for fully resonant interaction within a spectral resonance window. Beyond the corresponding spectral resonance threshold, the wave triad members stop interacting and follow linear WKBJ dynamics. The universal resonance threshold is adaptive to changes in triad wavenumbers as well as background shear strengths and therefore is applicable for a wide range of wavenumbers and background shear. Only for wave amplitudes near the threshold of static instability, systematic biases may occur, possibly exacerbated by the so-far neglect of a direct, non-dissipative, transient GW impact on the mean flow, which has been shown by Bölöni et al. (2016) to potentially be as important as impacts by wave dissipation. We believe this is the first implementation of interacting internal GW triads into a WKBJ ray tracer taking into account the modulation of the waves by a slowly varying background.
The supplemented WKBJ code is validated against simulations from a wave-resolving model. In all cases considered, wave amplitudes and mean flow have been assumed to be horizontally homogeneous. Two wave packets are considered that generate a third one while spectrally passing through resonance for a range of vertical scales. Comparing the WKBJ ray-tracing simulations with corresponding wave-resolving simulations, we generally find good qualitative and quantitative agreement for the wave modulation and triadic interactions, provided the wave amplitude with respect to static instability does not exceed 0.1. However, it is clear that beyond this limit direct, non-dissipative GW-mean-flow interactions are not negligible anymore.
Depending on the strength of the mean-flow shear, two interesting regime limits emerge. On the one hand, in weak shear and at large vertical wavenumbers, nonlinear effects become visible that lead to differences between wave-resolving and WKBJ simulations. It is very well possible that this is due to the self-induced mean wind that the present WKBJ implementation does not take into account. On the other hand, in strong background shears, wave-wave interactions seem to become partially suppressed by the wave modulation on account of the mean flow. This is due to both the corresponding strong changes in wave energy density and the more rapid development of the GW wavenumbers, so that the time window for resonant triad interactions is narrower. Hence an eventual outcome of further studies of GW-GW interactions in the atmosphere might well be that wave modulation dominates the evolution of the GW spectrum in strongly sheared regions, such as jet streams, while triadic interactions dominantly shape the GW spectrum in more weakly sheared regions, such as the midlatitudes and polar regions during spring after the breakdown of the polar night jet. However, there the modulation of GW-GW interactions by the self-induced mean wind could be non-negligible, a process not taken into account by wave turbulence theory.
Obviously there still is some way to go until this picture is confirmed. Both varying stratification and rotation will have to be included in theory and numerical implementation. Compressibility effects should be considered, but most pressing seems to be an inclusion of the GW impact on the mean flow. While atmospheric winds are clearly horizontally inhomogeneous, GW parametrizations are typically single-column implementations and do not take into account the lateral propagation of realistic internal GW packets. Yet, a three-dimensional implementation of a WKBJ ray tracer could potentially carry over the concepts for the wave-wave interaction applied here with an accordingly adapted strategy for the ray-tracing geometry. It would also be of interest to consider clusters of GW-GW interactions with common triad members (cf. Walsh and Bustamante, 2020), and validate the approach against these. Finally, a challenge will be continuous GW spectra. Our theory still assumes that the GW spectrum has distinct peaks that are sufficiently separated to allow for the discretely polychromatic GW fields considered here. Smoother spectra will need further theoretical developments, which would seem to be worth the effort.