Mixed triad wave–wave interactions between Rossby and gravity waves are analytically derived using the kinetic equation for models of different complexity. Two examples are considered: initially vanishing linear gravity wave energy in the presence of a fully developed Rossby wave field and the reversed case of initially vanishing linear Rossby wave energy in the presence of a realistic gravity wave field. The kinetic equation in both cases is numerically evaluated, for which energy is conserved within numerical precision. The results are validated by a corresponding ensemble of numerical model simulations supporting the validity of the weak-interaction assumption necessary to derive the kinetic equation. Since they are generated by nonresonant interactions only, the energy transfers toward the respective linear wave mode with vanishing energy are small in both cases. The total generation of energy of the linear gravity wave mode in the first case scales to leading order as the square of the Rossby number in agreement with independent estimates from laboratory experiments, although a part of the linear gravity wave mode is slaved to the Rossby wave mode without wavelike temporal behavior.
The dynamics of rotating stratified flows are often described using linear Rossby and inertia–gravity waves (hereinafter called gravity waves), which have well-known properties. The interactions between these linear waves generated by the nonlinear terms in the equations of motions, however, are less well known but can profoundly affect the evolution of the flow field. The nonlinear wave–wave interactions within the gravity and Rossby wave components along with the mixed interactions between the two wave modes are the subject of this paper. We discuss the energy transfers by weak gravity, Rossby, and mixed wave–wave interactions in models of different complexity.
Rossby waves evolve at temporal scales much longer than the inertial period. Wave–wave interactions of Rossby waves are typically characterized by motions with spatial scales on the order of the deformation radius and larger, in both the atmosphere and the ocean, and the interactions are known to transport energy predominantly to larger scales (e.g., Rhines 1975; Salmon 1998). Internal gravity waves, on the other hand, evolve at much shorter temporal scales and also often have small spatial scales down to a few meters, which renders them hard to observe as well as to resolve in models. Olbers (1976), Pomphrey et al. (1980), Müller et al. (1986), and Henyey et al. (1986), for example, have discussed how wave–wave interactions of gravity waves in the ocean can lead to an energy transfer within the gravity wave spectrum to smaller scales through processes such as induced diffusion (e.g., Müller et al. 1986) and parametric subharmonic instability (e.g., Sutherland 2010).
Rossby wave–wave interaction is often discussed and understood as an energy flux from one wavenumber to its neighboring wavenumber, that is, local in wavenumber space, generated by the nonlinearities in the model under investigation. Similar to the equivalent energy flux in the theory of isotropic turbulence by Kolmogorov (1941), scaling arguments can then lead to an inference of the magnitude and wavenumber dependency of the flux and the energy spectrum of Rossby waves (Salmon 1998). Such an energy flux can also be diagnosed from observations or models (Frisch 1995) and shows indeed a direction toward large scales (Scott and Wang 2005). On the other hand, gravity wave–wave interactions are usually discussed in the framework of resonant triad interactions, where three waves of different wavenumbers and frequencies—satisfying a certain relation given below—exchange energy. Such triad interactions can generate therefore an energy transfer between rather different wavenumbers, that is, nonlocal in wavenumber space. Although the nonlocal triad interaction and the local energy flux framework are in principle equivalent, the former presents more information since all waves can participate in the triad interactions forming a complicated network. The net effect of these complex triad interactions can result in an energy transfer in either upscale or downscale direction, which is then diagnosed as the local energy flux in wavenumber space.
Understanding the triad of interactions between gravity waves, Rossby waves, and their mixed interactions therefore can contribute to better understanding of the dynamics of atmosphere and ocean. Wave–wave interactions have been discussed for the case of Rossby waves (e.g., Kenyon 1964; Longuet-Higgins et al. 1967; Connaughton et al. 2015) and gravity waves (e.g., Olbers 1976; Nazarenko 2011) based on the weak-interaction assumption and the kinetic equation formulated for geophysical applications first by Hasselmann (1962). The wave–wave interaction between the two principal wave modes in this framework, however, has attracted less interest so far, an exception being Lelong and Riley (1991). In a general derivation for models of different complexity, we here also allow for the resonant and nonresonant interactions of all wave modes present, that is, mixed triad interactions between (internal) gravity waves and Rossby waves, referred to as Rossby–gravity mixed interactions.
Although usually separated by vastly different time scales, it is known that both wave modes can interact with each other (e.g., Müller 1977; Lighthill 1978; Ford et al. 2000). A general framework for the interaction between long and short waves has been set by Benney (1977), investigating nonlinear differential equations that allow, in the linearized case, for long and short wave solutions. The theory can be applied to the resonant interaction of two short internal gravity waves with a long Rossby wave. The resonance condition of such a triad of waves is an approximate equality of the long wave phase velocity and the short wave group velocity. More recently, spontaneous emission or loss of balance (e.g., Vanneste 2013) has attracted considerable attention and refers to the process of gravity wave generation by a balanced flow spontaneously in the absence of any external forcing, which often happens during baroclinic instability (Plougonven and Snyder 2007). Spontaneous generation of internal gravity waves from the balanced flow has been studied in laboratory (e.g., Williams et al. 2005, 2008) and numerical studies (e.g., Molemaker et al. 2005; Sugimoto and Plougonven 2016; Chouksey et al. 2018), where it was shown that this process can transfer energy from the balanced flow to gravity waves. It was speculated that this process may represent a significant energy sink for the balanced flow in the ocean (e.g., Williams et al. 2008; Brüggemann and Eden 2015), which would modify our view of the energy cycle in the ocean and which would need consideration in parameterizations for ocean models. Here we also discuss the opposite scenario, the generation of Rossby waves or balanced flow by gravity waves, which might have importance for the energetics of the wave field.
This study is structured into four sections. Section 2 presents a theoretical layout of an evolution equation for the wave energies in spectral space—the so-called kinetic equation or scattering integral—that is derived and discussed in Fourier space for a quasigeostrophic layer model, for a single-layer reduced-gravity model, and for vertically resolved models using the weak-interaction assumption. The derivation is detailed mainly for the reduced-gravity model but also holds for the other models. From the general kinetic equation, Rossby wave, gravity wave, and mixed Rossby–gravity wave interactions are derived. Section 3 demonstrates numerically the generation of linear gravity waves from balanced flow and vice versa and validates the results with a numerical model. Section 4 presents a summary and discussion of the results.
2. The kinetic equation or scattering integral
We use as an algebraically simple model a reduced-gravity model for a single layer. All variables and also the differential operators are scaled with a velocity scale U, length scale L, and time scale 1/Ω, where Ω denotes the magnitude of the constant Coriolis parameter. This yields
where Ro = U/(ΩL) is the Rossby number, u denotes1 the scaled layer velocity, and h is the layer thickness perturbation scaled by g′/(ΩLU), with g′ being reduced gravity. Furthermore,
is the scaled phase velocity of gravity waves without rotation, is the mean layer thickness, and
denotes the Froude number. The scaled Coriolis parameter f = 1 was kept for reference. The model is considered on a double-periodic domain.
For Ro ~ Fr ≪ 1 the quasigeostrophic approximation is valid and yields in this scaling
where potential vorticity and the geostrophic streamfunction is defined by , and where the Rossby radius LR = c/f = 1 is kept for reference. Here, a small variation of the Coriolis parameter Ω on the same order of magnitude as Ro and Fr was introduced, represented by the reference parameter β ~ Ro, but Ω is kept constant in the reduced-gravity model.2
As a more complex and realistic model we also consider the primitive equations, for which Boussinesq and hydrostatic approximation are applied to the full system. Using as time scale, horizontal and vertical length scale 1/Ω, L, and H, respectively, and scaling horizontal velocity u as U, vertical velocity w from the continuity equation as UH/L, pressure p from geostrophic balance as ΩUL, and buoyancy b from hydrostatic balance as ΩUL/H yields
with the diagnostic relations ∂zp = b and ∇ ⋅ u + ∂zw = 0; the (scaled) stability frequency N = Ro/Fr = LR/L, which is kept constant in the vertical; the Froude number ; and the Rossby radius . The unscaled gravity wave speed without rotation becomes in this model , where denotes the unscaled stability frequency. Here, f = 1 was again kept for reference. For the limit of Ro ≪ 1 (but fixed relation Ro ~ Fr, LR ~ L) and introducing again a small variation in the Coriolis parameter, the quasigeostrophic approximation of Eq. (3) also yields Eq. (2) but with the three-dimensional definition of potential vorticity given by
a. The quasigeostrophic layer model in Fourier space
Before turning to the more complex models, the quasigeostrophic layer model is discussed first because of its simpler algebraic form. The approach taken here, however, is analogous to the other models. A Fourier ansatz given by
with the notation and the Rossby wave frequency . Multiplying with exp(−ik0 ⋅ x) and integrating over x yields
with the delta function
The right-hand side of Eq. (6) contains the effect of the nonlinear terms on the right-hand side of Eq. (2). The wave amplitude still contains the regular wave oscillation in time, which can be removed by setting
(frequently called “interaction representation”; see, e.g., Nazarenko 2011). Here, represents the frequency of the Rossby wave for the wavenumber vector k0 and a0(t) = a(k0, t) is the amplitude of the wave at k0. The wave amplitude is allowed to change in time, and the scaling of the amplitude a0 by f is introduced to match a0 with total wave energy in the reduced-gravity and primitive equation models below.
Substituting the interaction representation in Eq. (6) yields
with an = a(kn) and . Equation (7) is a tendency equation for the amplitude a0 of the Rossby wave at wavenumber k0. Only if the right-hand side of Eq. (7) is nonzero, that is, if there are wave–wave (triad) interactions by the nonlinear terms in Eq. (2), will the amplitude a0 of the Rossby wave change in time. By renaming the integration variables, Eq. (7) can be rewritten as
with the Rossby wave interaction coefficients
b. The reduced-gravity model in Fourier space
An equation analogous to Eq. (8) can be derived for the reduced-gravity model, but for both Rossby and gravity wave amplitudes, instead of Rossby wave amplitudes only as in Eq. (8). We apply a Fourier ansatz to the reduced-gravity wave model in Eq. (1) given by
with complex wave amplitudes with , and similar for h. This yields in Eq. (1) after multiplication with exp(−ik0 ⋅ x) and integration over x
with the system matrix (k0), the state vector , and the vector function , which are given in appendix A. The term contains all linear terms in Eq. (1) while the integral involving the vector function contains all nonlinear terms. Equation (11) is analogous to Eq. (6) of the Rossby wave case but is now a vector equation. Standard linear algebra methods, described below, are used to cast it into a scalar equation similar to Eq. (6).
The matrix has three eigenvalues ω(0) = 0 and ω(±) = ±(f2 + c2k2)1/2. Two of them, ω(±), correspond to gravity waves and the other one, ω(0), corresponds to Rossby waves, but here with vanishing frequency instead of ωR since we have assumed a constant Coriolis parameter f. It is shown in appendix B that by including a small variation in f the vanishing eigenvalue of the system gets finite and ω(0) = ωR. But the magnitude of ω(0) remains still much smaller than |ω(±)| for Ro ≪ 1. Related to the eigenvalues, the matrix has left and right eigenvectors Ps and Qs [discussed also by, e.g., Hasselmann (1970), Leith (1980), and Olbers et al. (2012)] given by
Any state vector can be expressed by the eigenvectors as
Using the decomposition equation [Eq. (13)] in Eq. (11) and scalar multiplication with Ps(k) yields a scalar evolution equation for gs(k) similar to Eq. (6) of the Rossby wave case. But before doing so, the wave oscillation in time is removed by setting
to obtain an evolution equation for the wave amplitude as as before for the Rossby wave case. It is given by
with s0 = 0, ±, the notation , , and the exchange coefficients C′, which follow using the decomposition [Eq. (13)] in the vector function and scalar multiplication with Ps(k0). The interaction coefficients can be made symmetric in the last two indices by renaming the integration and summation variables as before going from Eq. (7) to Eq. (8) in the Rossby wave case. The symmetric interaction coefficients are given in appendix A, replacing . Since
c. The vertically resolved models in Fourier space
The derivation of a corresponding amplitude equation for the primitive equation model Eq. (3) and its quasigeostrophic version is very similar to the single-layer models. For the vertically resolved quasigeostrophic model the wave ansatz
is applied, where K = (k, m) denotes the three-dimensional wavenumber vector and m is the vertical wavenumber. This yields an amplitude equation for the vertically resolved quasigeostrophic model
with the Rossby wave frequency ωR = −βkx/(k2 + f2m2/N2) and with the Rossby wave interaction coefficients
For the primitive equation model, a wave ansatz analogous to Eq. (16) is applied to u, p, and b. This leads after the same manipulations as explained above to an amplitude equation analogous to Eq. (14), with the only difference being that the integrals are taken over K instead of k. The analogous definitions of , , and , together with eigenvalues of and the corresponding interaction coefficients , are given in appendix C.
d. Wave energy and weak-interaction assumption
The sum of kinetic energy and potential energy of the waves (averaged over one wave period) in the reduced-gravity model is related to
where ns results from the normalization of the eigenvectors and is defined in appendix A. The square of the wave amplitudes as is thus related to wave energy. This definition of energy is the same in the primitive equation model, but with potential energy given by in this case. The derivation of the energy equation below is completely analogous for the reduced-gravity model and the primitive equation model, as well as for the quasigeostrophic counterparts. We will sketch the derivation here for the reduced-gravity model only.
To obtain an evolution equation for the square of the amplitudes we multiply Eq. (14) by , we multiply the complex conjugate of Eq. (14) for k0 → kn (but keeping s0) by , and we add both, which yields
For kn = k0, Eq. (20) describes the evolution of the total energy of the wave at k0—besides a scaling factor —and for kn ≠ k0, the equation describes the evolution of covariances between the waves of type s0 at different wavenumbers. However, to evaluate this evolution, triple covariances of the amplitudes need to be known, which enter on the right-hand side of Eq. (20). Formulating equations for such triple covariances would involve quadruple covariances and so forth, resulting in the typical closure problem of nonlinear systems or turbulence. By invoking the weak-interaction assumption as described by, for example, Hasselmann (1962, 1966) and Nazarenko (2011), however, the triple covariances on the right-hand side of Eq. (20) can be cast into a form involving the energies of the other waves of different wavenumber and/or mode. The result is a closed system of equations involving the energies of all modes that is called the kinetic equation or scattering integral.
The lengthy derivation starting from Eq. (20) is only sketched here: the amplitudes are decomposed into a power series in the small parameter Ro. The zero-order variables in this power series are assumed to be time independent; thus it is assumed that the amplitudes are only weakly changing. Furthermore, it is assumed that the zero-order variables are statistically independent and can be described by Gaussian statistics. Note that this second assumption is used by Hasselmann (1966) but can also be replaced by a weaker assumption on the statistics of the amplitudes (Nazarenko 2011). The first assumption on the weakly changing amplitudes is, however, essential to the procedure. Using then Eq. (14), the first-order amplitudes can be expressed to first order by integrals of products of two zero-order amplitudes. Using this, the triple covariances in Eq. (20)—which vanish for the zero-order amplitudes—can then be cast to first order into quadruple covariances of the zero-order variables. Because of the assumption of Gaussian statistics, the quadruple covariances can in turn be written as products of two second-order covariances, which are related to the total wave energy .
The resulting kinetic equation is given as
is the total energy of the wave of mode s at wavenumber kn, with k2 = k0 − k1, and where c.c. denotes the complex conjugate of the proceeding terms. For the vertically resolved model, replace kn with Kn. The time-dependent function Δ(ω, t) is given by
If the condition is met, there is a resonant triad interaction. Such resonant interactions are most important in the long run, but nonresonant interactions can also become important for short time intervals. The last term in Eq. (21) is never resonant and is related to inertial oscillations. For k0 = 0, a further nonresonant contribution related to inertial oscillations has to be added, which is not given here.
e. Energy conservation
Integrating the kinetic equation Eq. (21) over k0 and summing over s0 yields total energy conservation. Renaming one of the integration variables and changing the sign of another, using ω−s(−k) = −ωs(k), , and , yields
The term in brackets in the second line vanishes for the coefficients of the primitive equation model and the layered and vertically resolved quasigeostrophic models, such that total energy of the models is indeed conserved. Here only the first two lines of Eq. (21) have been used, the last term related to inertial oscillations balances with the term at k0 = 0.
However, the coefficients of the reduced-gravity model do not vanish in the term in brackets in the second line of Eq. (23). The reason is that the total energy of this model is not a second-order quantity as for the primitive equation model, but involves cubic terms. Only for a linearized version of the reduced-gravity model the energy gets a second-order quantity. Thus, wave energy and full energy differ in the reduced-gravity model, and Eq. (23) does not describe the total energy of the full reduced-gravity model as is the case for the primitive equation and the quasigeostrophic models. This nonconservative effect is, however, only a few percent at maximum in the numerical evaluations below.
f. Rossby wave interaction
Restricting the sums in Eq. (21) to s0 = s1 = s2 = 0—that is, to Rossby waves only—energy is already conserved and the kinetic energy equation for single-layer Rossby waves is recovered (Kenyon 1964; Longuet-Higgins et al. 1967):
with k2 = k0 − k1, with , the coefficients , and using . The function denotes taking the real part of its argument. The interactions become resonant if Δω = 0. If β = 0, all interactions are resonant, but for the quasigeostrophic models, only certain interactions become resonant.
For the vertically resolved Rossby waves, the kinetic equation is given by replacing kn with Kn in Eq. (24) with the symmetric interactions coefficients
and by replacing L−2 with , , and in the first, second, and third set of parentheses, respectively, in the fraction.
g. Gravity wave interaction
Energy is also conserved restricting the sum in Eq. (21) to s0 = s1 = s2 = ±, that is, to the gravity waves only. After several manipulations, resolving the double sum, and using general properties of the energy, the frequency, and the coefficients, the kinetic equation can be rewritten in simplified form (three terms instead of eight) as
The corresponding expression for is just the complex conjugate of the terms explicitly given on the right-hand side. For the generalization to three dimensions, just exchange k0 with K0, and correspondingly for the other wavenumber vectors. The last term cannot become resonant and is usually not considered; the first term and the middle term are called sum and difference interaction, respectively. It can be shown that for the reduced-gravity model the resonance conditions for gravity waves cannot be satisfied3 but that resonant sum and difference interactions are possible for the three-dimensional (internal) gravity waves. Similar forms of the kinetic equation for nonhydrostatic internal gravity waves have been derived by, for example, Olbers (1976), and a recent overview is given in Lvov et al. (2012). The corresponding interaction coefficients C in Eq. (26) for a nonhydrostatic model (not given) can be related to the ones given in Olbers (1974, 1976) and Müller and Olbers (1975) for the resonant interactions but differ in general for the nonresonant interactions.
h. Mixed Rossby–gravity wave interaction
As a consequence of the individual energy conservation of Rossby and gravity wave modes, the same is true for the remaining terms, that is, the mixed components of the sum describing the interaction of Rossby and gravity waves. Considering Eq. (21) for a state with initially E±|t=0 = 0 for all k—that is, initially without gravity waves—and energy only in the balanced mode E0 yields for s0 = ± the following:
This term cannot become resonant since max(|ωR|) = βLR/2 ≪ 1 and min(ω0) = f = 1 for k0 = 0. Only if β ~ Ro gets larger can the term become near resonant. However, a finite amount of net gravity wave energy is generated during an initial adjustment, as demonstrated in the next section, which we interpret as generation of the linear gravity wave mode by the balanced flow.
The reversed case is given if for all k in the initial state and energy only in the gravity wave mode E±. Then we find for s0 = 0 the following:
The first (sum interaction) term and the last term within the square brackets can hardly become resonant, but the second (difference interaction) term can become resonant. It is possible to show that holds for the coefficients of the reduced-gravity model. The gravity wave mode is then completely separated, but this is not the case for the primitive equation model. Internal gravity waves can thus also generate balanced flow.
The difference interaction term in Eq. (28) is analogous to the mechanism discussed by Benney (1977); here, however, on the level of the kinetic equation rather than of the governing differential equations. In Benney’s discussion, the gravity wave partners 1 and 2 of the resonantly interacting triad are short, the generated Rossby wave 0 is a long wave. A one-dimensional sketch of the resonance condition is shown in Fig. 1 of Williams et al. (2003). The gravity wave energy, propagating with the appropriate group velocity, moves with the phase speed of the long Rossby wave and energy can be transferred. Note, however, that we here do not restrict the triad interaction to be resonant nor do we confine the discussion to scale-separated interacting triads.
3. Numerical evaluations
a. Generation of linear gravity waves by balanced flow
To demonstrate the generation of linear gravity waves from the balanced mode, we consider a certain given spectral shape of the energy of the single-layer Rossby waves and initially no gravity waves; that is,
and E± = 0. The constant ER determines the total energy. The kinetic equation is evaluated numerically on an equidistant grid in k (or K) space with an odd number of grid points in all (positive and negative) directions including the point k = 0. It is found that the wavenumber grid needs to be of this shape to allow for exactly matching triads, for which all three wavenumbers can be exchanged. If that exchange is not possible by using a different grid, we were not able to conserve energy in the calculation of the kinetic equation. Energy conservation is checked by evaluating and is given within the numerical precision error for all examples that we checked for single-layer Rossby waves, the reduced-gravity model, and the vertically resolved model, except for the gravity waves in the reduced–gravity wave model because of the abovementioned complication of total energy that does not match wave energy in that model. [Using a modified form of the nonlinear terms in Eq. (1), that is,
yields total energy identical to the total energy of the linear waves. Only for nondivergent flow with ∇ ⋅ u = 0 do the nonlinear terms in Eq. (29) become identical to the ones in the full system of Eq. (1). The corresponding interaction coefficients derived from Eq. (29) (not shown) then indeed satisfy Eq. (23) up to numerical precision. These coefficients were used as a consistency check for the numerical algorithm for the kinetic equation of the reduced-gravity model.]
A power law with p = 2 yields a so-called stationary spectrum for E0 with ∂tE0 = 0. Energy becomes stationary since the kinetic equation [Eq. (24)] for the Rossby waves can be shown to vanish identically for p = 2 (Kenyon 1964), which is also true for the numerical evaluation of Eq. (24) within numerical precision. In a forced dissipative model, however, the spectrum cannot be stationary, and scaling laws suggest p = 3 and an inverse energy cascade (e.g., Salmon 1998). The numerical evaluation of the kinetic equation [Eq. (24)] for p = 3 and L = 1 on a 40 × 40 spatial domain for different grid resolutions is shown in Fig. 1. The Rossby wave frequency ωR is set in this case to zero, which results in a symmetry of ∂tEs with respect to the wavenumber angle. Therefore, we show ∂tEs integrated over wavenumber angle in Fig. 1 and the following figures in this section.
All interactions become resonant and the kinetic equation [Eq. (24)] simply becomes a linear function of time. Figure 1 shows energy gain (loss) at small (large) k indicative of an inverse energy cascade for all grid resolutions. The magnitude of the energy gain at low k depends on the grid resolution; for the grid with 1272 grid points it is roughly half as large as for the grid with 10232 grid points. At highest resolution, however, the magnitude of the maximal energy gain appears to converge. The region of energy loss extends farther toward larger k for higher grid resolution, but the energy gain for k > 9 for the grid with 1272 grid points (and k > 18 for 2552 grid points) is an artifact of the numerical evaluation on a grid with finite k, since it moves to larger k using higher resolution. Note that these dependencies on grid resolution are similar to those of numerical models (see below). The effect of nonvanishing ωR is an asymmetry of ∂tE0 with respect to the wavenumber angle (not shown). The region of energy gain is shifted toward smaller zonal wavenumber, while its meridional wavenumber range stays similar to the region of energy gain in Fig. 1a). This means that a zonalization of the energy by the β effect can be seen, as first described by Rhines (1975).
To validate the results of the kinetic equation, we use an ensemble of numerical model simulations with the same parameters as used for the evaluation of the kinetic equation. Each ensemble member is initialized with the same energy distribution
and E±|t=0 = 0 but with a random phase. Specifically, we are using as initial condition for the model the Fourier transform of the state vector , where ζ(k) is a complex normal number with 0 mean and a variance of 1. Different realizations of ζ are used for each ensemble member using a random number generator. For consistency, we use eigenvectors Qs (and Ps) appropriate to the discrete model equations given in appendix D instead of the analytical version of the eigenvectors Qs (and Ps) given in appendix A. The discretization of the nonlinear terms follows the energy-conserving scheme by Sadourny (1975) in the momentum equation, and we use second-order differences for the thickness advection. Energy density Es is diagnosed from the model ensemble using the eigenvector Ps, and the energy transfer ∂tEs is calculated from finite differences of Es.
Damping in the model is necessary for a stable integration and given by biharmonic friction and diffusion and by the time stepping scheme which is chosen as a quasi-second-order Adam–Bashforth interpolation with adjusted weights to allow for stable simulations of gravity waves of highest frequency. Both damping effects decrease the energy proportional to Es since both effects are linear. To eliminate the damping effect from the diagnosed energy transfers in the model, we subtract from ∂tE0 diagnosed from the model ensemble, the corresponding ∂tE0 diagnosed from an identical model ensemble but without nonlinear terms, that is, setting Ro = 0 in Eq. (1). Since the damping effect is linear and acts in a very similar way in the ensemble of the linear and the full model, we can isolate the energy transfers due to the nonlinear terms in this way and can compare them with the predictions from the kinetic equation.
Figure 2 shows that the diagnosed ∂tE0 in the model ensemble for three different times has indeed the expected shape of an inverse energy cascade. As for the prediction by the kinetic equation, ∂tE0 increases linearly in time, and also at a similar rate. The energy transfer depends in the model, however, on the details of the discretization of the nonlinear terms and in particular on model resolution. At small wavenumbers, ∂tE0 is lower than predicted by the kinetic equation using the same grid, but ∂tE0 in the model increases by increasing the model resolution. The region of energy loss is also shifted toward smaller k relative to the prediction by the kinetic equation but appears to converge to the prediction using higher resolution.
With the kinetic equation [Eq. (21)], we can also describe the generation of linear gravity wave energy by nonresonant interactions due to the mixed components in the first sum of Eq. (21). The energy gain of the gravity mode E± is given for this case by Eq. (27), and similar for the balanced mode E0. Figures 3a and 3b show the energy gain of E0 and E± by the mixed components as a function of time and wavenumber k. Both functions oscillate in time with smaller frequencies (but larger amplitude) for smaller k. There is also a net total energy gain integrated over k and time for the gravity wave model E± and a corresponding loss of energy for the balanced mode E0. The largest energy gain of E± shows up near k = 1, that is, close to the Rossby radius. The effect also shows up for different spectral laws and even for the stationary spectrum p = 2. For the same total energy level in E0, it is getting stronger (weaker) for p < 3 (p > 3) and is for p = 2 approximately 2 times as large as in Fig. 3c. To test whether the nonconservation feature for total energy in the reduced-gravity model modifies the gravity wave generation shown here, we repeat the evaluation in Fig. 3 with the model given by Eq. (29), in which energy is conserved exactly. This does not change the results. Including the β effect also does not change the results in terms of the gravity wave generation.
Figure 4 shows ∂tE± and the increase of E± integrated over k in the model ensemble. The time evolution of ∂tE± is very similar to the prediction of the kinetic equation, although the amplitudes are larger such that the total increase in E± is also larger than predicted. Different from the energy transfer of the balanced mode, ∂tE± depends only weakly on grid resolution in the model but strongly on details of the discretization (not shown). It is therefore likely that the mismatch in the amplitude of the wave generation between model and prediction is due to the model numerics.
b. The role of the slaved gravity wave mode
Part of the generated energy in the total linear gravity wave mode in the experiments discussed in the previous section is related to the so-called slaved gravity wave mode (e.g., Leith 1980; Vanneste 2013). To show this, consider Eq. (14) without applying the interaction representation. The state vector projection is then governed by
The gravity waves are initially vanishing and then—according to Hasselmann’s assumption—weakly growing, such that g± = O(Ro). The time derivative in Eq. (30) includes the fast oscillations with frequency and the slow growth of the amplitude. Kafiabad and Bartello (2018) introduce therefore a fast and slow time scale with T = Rot* and acting on gs. This yields to first order in Ro for the gravity mode s0 = ± as follows,
since the interaction coefficients C are also of O(Ro). The nonzero right-hand side will generate an oscillating fast mode g± even if g±|t=0 = 0. To suppress this oscillation it is necessary that or
The part of the state vector (or its inverse Fourier transform) is also called the ageostrophic balanced or slaved wave mode and is identical to the correction in the first iteration of the nonlinear balancing method by Machenhauer (1977) or Baer and Tribbia (1977). Since it is given by the coefficients g0, which vary on the slow time scale T only, the same holds for the slaved mode f±. It is therefore not fast but slow and should be considered to be part of the balanced mode. When the model is initialized with the fast oscillation is eliminated to first order, which makes the nonlinear balancing method useful for, for example, initializing numerical weather forecast models with observed states. In Kafiabad and Bartello (2016, 2018) and Chouksey et al. (2018), the method is used to diagnose gravity wave activity in model simulations.
An exemplary time series of the inverse Fourier transform of the actual wave mode in terms of the model’s height field h in one ensemble member, together with the same quantity related to the slaved mode fs, is shown in Fig. 4c. Our method to calculate the slaved mode follows Chouksey et al. (2018), except that the eigenvectors that we use are the ones for the discrete system from appendix D and that the nonlinear function is evaluated by integrating the model one time step, in order to be fully consistent with the discretization in the model. While h related to the slaved mode fs is only slowly varying, h related to the total linear wave mode g± shows fast oscillations on top of the slow variation, but both are on the same order of magnitude. The total energy related to f± (Fig. 4b) is about one-half of the total energy related to g±. It is likely that the same holds for the prediction by the kinetic equation, but here it appears not possible to differentiate between slaved mode and wave mode.
c. Generation of balanced flow by gravity waves
To demonstrate the generation of balanced flow by the gravity waves by resonant and nonresonant interactions, we consider the vertically resolved model since here only resonant interactions of gravity waves are possible (to first order). Such interactions are believed to generate the so-called Garrett–Munk spectrum given in dimensional form by
(Cairns and Williams 1976; Munk 1981), where the functions A and B yield power laws m−2 and ω−2 for large arguments and where the parameter and c* = 0.5 m s−1 correspond to a wavelength of about 630 m for and small ω. The functions A and B are normalized such that the total energy of the wave field is given by Egm = 3 × 10−3 m2 s−2. The spectrum, the energy transfers within the gravity wave field inferred from the kinetic equation, and its dependency on parameters are discussed in detail in a companion study (Eden et al. 2019) for the nonhydrostatic case, whereas we concentrate here on the interaction of the wave field with the balanced flow.
Figures 5a and 5b show the content (or variance preserving) spectra of the energy transfer ∂tE+ by the resonant interactions, excluding the nonresonant interactions as a function of vertical wavenumber m and horizontal wavenumber amplitude k. Note that since E+ is symmetric in m and isotropic in k, the figure shows ∂tE+ only for positive m and integrated over horizontal wavenumber angle. The method we use to calculate the resonant interactions only are described in detail in Eden et al. (2019) (their method 3); it uses the δ function in the wavenumbers in Eq. (26) to resolve the integrals in (k2, m2) and the δ function in the frequencies to resolve the integral over m1. Only the first (sum interaction) term and the second (difference interaction) term in Eq. (26) then contribute to the transfer. For the numerical calculation we use a grid of 8012 × 1601 grid points in wavenumber space and a domain size of 200-km horizontal extent and 10-km vertical extent, Ω = 10−4 s−1, and . Note that wavenumbers and energies are shown in the figure in their dimensional form, corresponding to Ro = 0.1 and Fr = 2 × 10−3. The transfers are similar to the ones shown and discussed in detail in Eden et al. (2019) for nonhydrostatic waves; differences are only seen for .
For the generation of balanced flow by the wave field given by Eq. (28), only the second term can become resonant, and it is shown in Fig. 5c for a finite Rossby wave frequency. The generation of balanced flow is orders of magnitude smaller than ∂tE+. It turns out to be proportional to β; that is, it vanishes on an f plane where ω0 = 0. Thus, similar to the generation of gravity waves by the balanced flow, resonant interactions play also no significant role for the generation of balanced flow by gravity waves. As before, this changes again for the nonresonant interactions in Eq. (28), as demonstrated next.
Figure 6a shows the nonresonant energy transfer ∂tE0 given by Eq. (28) as a function of time t. Here, method 1 by Eden et al. (2019) is used; that is, all resonant and nonresonant interactions on the three-dimensional wavenumber grid are calculated. Since this method is numerically more expensive than for the resonant interactions, we use a smaller grid of 2013 grid points and a smaller domain size with 100-km horizontal extent and 10-km vertical extent. Similar to the generation of gravity waves by the balanced flow shown in Fig. 3a, there is initially a strong generation of balanced flow by the nonresonant, mixed terms given by Eq. (28), while after a fraction of a day the generation of balanced flow ceases and oscillates around zero. For t → ∞ the state of Fig. 5c will be reached. Figure 6b shows the cumulative time integral of the generation, which becomes stationary after about one day. For comparison, the figure also shows the energy transfers within the gravity wave field only in terms of |∂tE+| related to these interactions (the integral of ∂tE+ vanishes because of energy conservation), which do not decrease and lead to a continuous increase of the time integral.
4. Summary and discussion
Using the weak-interaction assumption by Hasselmann (1966)—that is, assuming that the wave amplitudes are changing only slowly—kinetic equations for the nonlinear wave–wave interaction are derived that are closed with respect to wave energy in wavenumber space. The kinetic equations apply to a simple two-dimensional layered model and the primitive equations, in both their full versions and their quasigeostrophic versions. The result is identical to previous studies with respect to Rossby and internal gravity wave interaction (e.g., Kenyon 1964; Longuet-Higgins et al. 1967; Olbers 1976; Nazarenko 2011), but here we also allow for the interaction of all waves present in the models, that is, (internal) gravity waves and Rossby waves. The latter degenerate to the geostrophic mode with vanishing frequency on the f plane.
Since the interactions only within the gravity wave field or only within the Rossby wave field conserve energy by themselves, the same is true for the mixed interactions between linear gravity and Rossby waves in the kinetic equation. We consider here the case of initially vanishing linear gravity wave energy in the presence of fully developed Rossby wave field with a spectral power law of E0 ~ k−3 in the layered model and the reversed case of initially vanishing linear Rossby wave energy in the presence of a fully developed gravity wave field in the primitive equation model with a spectral power law of E± ~ k−2. The kinetic equations for the two cases are numerically evaluated, for which energy is conserved within numerical precision.
In the first case, the kinetic equation predicts for the Rossby wave field an energy transfer indicative of an inverse energy cascade. A numerical model ensemble with identical parameters and initialized with the same energy as used for the kinetic equation shows a similar magnitude and shape of the energy transfer, and, as in the kinetic equation, also a linear increase in time of the energy transfers. This validation supports the validity of the weak-interaction assumption necessary for the derivation of the kinetic equation for the Rossby wave case.
Only nonresonant interactions can generate linear gravity wave energy from balanced flow, which are less important in the long run than the resonant interactions between the Rossby waves but which can still generate finite energy transfers. The evaluation of the kinetic equation of the generation of linear gravity waves by the nonresonant interaction from a fully developed Rossby wave field shows energy transfers that are maximal near the Rossby radius and that are oscillating in time. Integrating the energy transfers in time yields a small but positive generation of linear gravity wave energy for large t. A similar time evolution of the energy transfers related to the generation of linear gravity waves is seen in the model ensemble, supporting the validity of the weak-interaction assumption also for this case. Both in the model ensemble and for the kinetic equation, the gravity wave generation scales as Ro2, which can be seen in Fig. 7. This can be easily understood by rewriting Eq. (27) as
with a coupling coefficient αI(k0) derived from the dimensionless integrals in Eq. (27) and assuming that E0 ≈ const and E± ≪ E0 for t < Δt. The integral on the right-hand side of the expression for the integrated generated energy only depends on c (and time Δt), which is set constant in all experiments.
The reversed case is evaluated using the primitive equations and a wave spectrum following E± ~ k−2, that is, a Garrett–Munk spectrum of internal gravity waves. In this case resonant interactions between the gravity waves and the balanced flow for initially vanishing balanced flow are in principle possible but turn out to be orders of magnitude smaller than the resonant interactions within the gravity wave field in the numerical evaluation. Similar to the reversed case, however, nonresonant interactions can generate a small amount of energy in the balanced flow, which seems to follow the same scaling as for the reversed case of the generation of gravity waves from balanced flow. The validity of the weak-interaction assumption for internal gravity wave interactions is discussed in Eden et al. (2019). It turns out that only for the high-frequency nonhydrostatic waves are the turnover time scales for those interactions getting smaller than the gravity wave periods.
Gravity wave generation by balanced flow with wavelength close to the Rossby radius is also found by Williams et al. (2008) in laboratory experiments of baroclinic instability and in many numerical model simulations (e.g., Plougonven and Snyder 2007). Williams et al. (2008) estimate a relation of the gravity wave energy generation to the square of the Rossby number Ro, identical as we found, while other studies suggest more complicated relations (Vanneste and Yavneh 2004). Brüggemann and Eden (2015) find by a subjective fit increasing small-scale dissipation for larger Ro in a numerical model of baroclinic instability. In a similar model setup, Chouksey et al. (2018) find, by using the nonlinear gravity wave diagnostic tool described above, a dependency close to Ro2 for the fraction of generated gravity wave energy, which agrees with the dependency found by Williams et al. (2008). Such a relationship of the wave energy emission on environmental parameters like the Rossby number can be used to physically constrain the coupling of eddy and wave energy reservoirs in the framework of consistent ocean model parameterizations, as suggested by Eden et al. (2014).
Here we obtain also a dependency close to Ro2 of the generation of linear gravity wave energy by balanced flow. However, the comparison to the above cited studies is hampered by the generation of the so-called slaved wave mode. This mode consists of gravity wave modes that need to be added to the balanced modes in order to suppress the generation of fast modes by the nonlinear terms to first order (in amplitude) in Ro. In the model ensemble, it can be shown that about one-half of the generated linear wave energy can be attributed to the slaved modes. Figure 7 shows that both the generation of the slaved mode and the remaining wave energy also scale with Ro2. In agreement, Kafiabad and Bartello (2018) find to first order (in amplitude) a dependency close to Ro2 for the fraction of ageostrophic (or gravity wave) energy generated by spontaneous generation in nonhydrostatic model simulations.
To resolve the difficulty of the slaved modes, it is possible to initialize the model ensemble with the balanced state , instead of using only , analogous to a nonlinear balancing method. Here, g0 and fs are the amplitudes of the balanced mode and the slaved mode, respectively. Doing so, much less linear gravity wave energy is generated in addition to the that related to , and this small energy generation scales with Ro4 as shown in Fig. 7. It is likely that this scaling behavior continues going to higher orders in the balancing. Using the energy related to in the kinetic equation, however, does not change the results much since |f±| ≪ |g0|, and Eq. (27) still dominates the evaluation of the kinetic equation.
In conclusion, the kinetic equation is able to predict the triad interactions leading to linear gravity wave (and Rossby wave) energy generation correctly up to second order in Ro but not beyond, but this is not surprising since it corresponds also to the error level in the derivation of the kinetic equation. In the setup discussed here, the slaved mode generation, which scales also with Ro2, complicates the predictions of gravity wave generation. Other situations such as, for example, stimulated loss of balance, with an energetic wave field coexisting with balanced flow, will be investigated in later studies with the kinetic equation that was derived here.
This study is a contribution to the Collaborative Research Centre TRR 181 “Energy Transfer in Atmosphere and Ocean” funded by the German Research Foundation. The authors are grateful to two anonymous reviewers for helpful comments.
Reduced-Gravity Model Properties
The linear system matrix (k), the state vector , and the vector function for the reduced-gravity model in Eq. (1) are given by
The eigenvectors of the system matrix are given by
for s = 0, ± and the horizontal components
Orthonormality holds (i.e., ) with
For the special case of pure inertial waves for k = 0 and ω± = ±f, the eigenvectors are given by
The symmetric interaction coefficients for the reduced-gravity model in Eq. (1) are given by
Effects of Variations in f
To obtain the first-order correction by the variation of the Coriolis parameter on the vanishing eigenvalue ω0 in the reduced-gravity model, we rewrite the linear version of Eq. (1) as
with the divergence δ = ∇ ⋅ u and the vorticity . Here, a small variation of the Coriolis parameter was introduced, that is, f → f + βy, before taking vorticity and divergence, but we need to take f constant again in Eq. (B1), in order to apply the Fourier transform. Equation (B1) becomes
after Fourier transform and using the state vector . The eigenvalues ωs of ′ are identical to the ones for , in particular ω0 = 0, for which we want to find the first-order (nonzero) correction. The eigenvectors of ′ are related to Qs and Ps from Eq. (A2) by a linear transformation T with and are denoted by Rs = T ⋅ Qs and Ss = Ps ⋅ T−1. To find the eigenvalues and eigenvectors of the full matrix (′ + Ro), we expand them in the small parameter Ro (or β) using ωs and Rs as zero-order and unknown first-order and , which yields
To zero order in Ro, we find the eigenvalue problem for ′ again and to first order . Multiplying with Ss yields
The first-order correction to the zero-order eigenvalues ωs can now be found with the eigenvectors of ′ given by
For s = 0, this yields in Eq. (B4) indeed the familiar Rossby wave dispersion relation with . The corrections to ω± are not relevant here, since they are small when compared with the zero order. The same holds for the eigenvectors.
Primitive Equation Properties
For the primitive equation model the linear system matrix (k, m), the state vector , and the vector function are given by
with . The second terms in the parentheses in result from the vertical advection and vanish for m1 = 0, and for m = 0 the third component of vanishes entirely. The eigenvalues of are given by ω(0) = 0 and , and the eigenvectors are the same as for the reduced-gravity model but depend via ω on m, along with the Rossby radius LR = cm/f.
The (asymmetric) interaction coefficients are given by
For m0 = 0 and m1 = 0, the corresponding terms vanish; the asterisk denotes complex conjugate. The normalization ns is the same as in Eq. (A4) but for LR = cm/f and ω as defined in this section.
Discrete Model Properties
The linear discrete equations for the reduced-gravity model are given by
with the finite-differencing operators and , where Δx is the grid spacing in the x direction, and the averaging operators and , and similar for the y direction. Dissipative and nonlinear terms and the time discretization are omitted in Eq. (D1) but are included in the numerical model. With ,where x = (iΔx, jΔy)T and similar for the other variables, the Fourier transform of Eq. (D1) yields
with , , , , and similar for , , and so forth. Note that and for Δx → 0. The system matrix then becomes
with real eigenvalues
and eigenvectors given by and with horizontal components
and the normalization
for s = 0, ±. The eigenvectors and eigenvalues converge to their continuous counterpart for Δx, Δy → 0.
The vector denotes anticlockwise rotation of the vector u by π/2; that is, = (−υ, u) for u = (u, υ).
The reason is that the double-periodic boundary conditions for which β ≠ 0 are possible for the quasigeostrophic model but not for the reduced-gravity model. Closed meridional boundaries could be used instead, which would put an additional constraint on the wave amplitudes in the next section.