Energy cascade in the Garrett-Munk spectrum of internal gravity waves

We study the spectral energy transfer due to wave-triad interactions in the Garrett-Munk spectrum of internal gravity waves (IGWs) based on a numerical evaluation of the collision integral in the wave kinetic equation. Our numerical evaluation builds on the reduction of the collision integral on the resonant manifold for a horizontally isotropic spectrum. We directly evaluate the downscale energy flux available for ocean mixing, whose value is in close agreement with the empirical finescale parameterization. We further decompose the energy transfer into contributions from different mechanisms, including local interactions and three types of nonlocal interactions, namely parametric subharmonic instability (PSI), elastic scattering (ES) and induced diffusion (ID). Through analysis on the role of each type of interaction, we resolve two long-standing paradoxes regarding the mechanism for forward cascade in frequency and zero ID flux for GM76 spectrum. In addition, our analysis estimates the contribution of each mechanism to the energy transfer in each spectral direction, and reveals new understanding of the importance of local interactions and ES in the energy transfer.


Introduction
Internal gravity waves (IGWs) are ubiquitous features of the ocean but are filtered out by the quasi-geostrophic description of the system.Although they generally account for only a small fraction of the kinetic energy of the overall ocean, their existence has profound physical significance: they play an important role in transferring momentum, heat and tracers across the ocean, and their breaking drives most of the turbulence that leads to the inhomogeneity of ocean mixing which in turn affects the large-scale circulation.
In the IGW continuum, energy is supplied at large scales by atmospheric and tidal forcings and is dissipated at small scales.Understanding the energy transfer across scales driven by nonlinear processes is one topic of central importance in physical oceanography.Such understanding will not only shed light on the physical interpretation of IGW spectral forms, generally considered as the Garrett-Munk spectrum (GM) and its variations (Garrett & Munk 1972, 1975;Cairns & Williams 1976), but also provide an estimate of the energy flux toward small dissipation scales (downscale energy flux) that is available for ocean mixing.The latter aspect, as one focus of our current work, is particularly important given the fact that information at small dissipation scales is difficult to obtain from both measurements and modeling.
While the downscale energy cascade of the IGW continuum can be excited by many mechanisms such as wave-eddy interactions (e.g., Watson 1985) and bottom scattering (e.g., Kunze & Llewellyn Smith 2004), in many cases, nonlinear wave interactions are considered as a major contributor to the cascade in abyssal oceans (e.g., Müller et al. 1986;Polzin & Lvov 2011).In quantification of spectral energy transfer due to nonlinear wave interactions, one critical tool is the wave kinetic equation (WKE) derived in the framework of wave turbulence theory (Zakharov et al. 1992;Nazarenko 2011).For systems with three-wave resonant interactions, the general form of WKE reads where n(p, t) is spectral action density with p being the vector of wavenumber and t the time, V is the interaction coefficient, ω is wave frequency, f p12 = n 1 n 2 − n p (n 1 + n 2 ) and f 1p2 = n p n 2 − n 1 (n p + n 2 ).The right-hand-side (RHS) of (1.1) is also referred to as the collision integral, which describes wave action evolution due to triad interactions.Other mechanisms such as generation/dissipation of IGWs and coupling with eddies/currents are not included.The existence of eddies and currents may be potentially important in nonlinear energy transfer (e.g., Kafiabad et al. 2019;Savva et al. 2021;Dong et al. 2020Dong et al. , 2023) ) but will not be the focus of the current work.For IGWs, p = (k x , k y , m), ω 2 = (N 2 k 2 + f 2 m 2 )/(k 2 + m 2 ) is the dispersion relation with k = (k 2 x + k 2 y ) 1/2 being the magnitude of horizontal wavenumbers, N the buoyancy frequency, and f the Coriolis frequency.The interaction coefficient V has been derived using various approaches in the literature (e.g., Olbers 1974;Müller & Olbers 1975;Olbers 1976;McComas & Bretherton 1977;Lvov & Tabak 2001, 2004;Lvov et al. 2010), which yield formulations shown to be equivalent on the resonant manifold (Lvov et al. 2012).
(1.2) Therefore, (1.1) provides the energy transfer rate through collections of triad interactions in the spectral space.Such WKE characterizes the spectral evolution in the kinetic (or nonlinear) time scale τ N L p and is valid only for weakly nonlinear waves whose linear time scale τ L p (= 2π/ω p , i.e., wave period) is (much) smaller than τ N L p , that is, the normalized Boltzmann rate (Nazarenko 2011;Lvov et al. 2012) The evaluation of IGW energy cascade for Garrett-Munk spectra based on (1.1) was first undertaken by McComas et al. in a series of papers (McComas & Bretherton 1977;McComas & Müller 1981a,b).A major argument made in these works is that the collision integral in (1.1) is dominated by three types of nonlocal (i.e., scale-separated in either vertical wavenumber or frequency, or both) interactions, namely parametric subharmonic instability (PSI), elastic scattering (ES) and induced diffusion (ID) (see schematic illustrations in Figure 1).McComas & Müller (1981a) further argued that the Garrett-Munk spectrum is in equilibrium with respect to ES so that the downscale energy flux can be calculated from PSI and ID contributions.This simplification allowed an analytical formulation of the downscale energy flux, which laid the foundation of finescale Figure 1.The wavenumber vectors and action balance for the three types of scale-separated interactions with p0 = p1 + p2, where the action transfer direction is based on GM spectra.For PSI and ES, +1 denotes one unit of action received by the wave mode as a sink (red dot), and −1 denotes one unit of action sent by a mode as a source (green dot).Induced diffusion can reverse direction with sinks and sources (grey dots), which is determined by the spectral slopes to be discussed in §3.2.For all three mechanisms, the action transfer regarding the highest-frequency wave p0 is always opposite to those regarding p1 and p2, with energy conservation guaranteed by ω0(∂n0/∂t) = ω1(∂n1/∂t) + ω2(∂n2/∂t).
While McComas et al.'s calculation provides an estimate of energy flux in the same order of finescale parameterization, the interaction mechanisms involved in the calculation suffer from physical inconsistencies that have never been completely resolved.As summarized in Polzin & Lvov (2011), at least two confusing paradoxes exist: (a) In frequency space, both dominant mechanisms of PSI and ID are believed to transfer energy toward low frequencies (i.e., backward cascade).This is not realistic for a balanced IGW continuum unless there is an energy injection at high frequency into the ocean, which is not known.Thus, there must exist a "missing" mechanism that moves energy to high frequencies to form a forward cascade.(b) The GM76 wave action spectrum is independent of vertical wavenumber m in the range of high m, which leads to a zero-flux state for ID (since diffusion requires gradient in m at least at the leading order).This makes obscure McComas et al.'s calculation where an artificial correction in ID has to be applied to enable cascade in the high-m high-ω range of the spectra and raises questions on what the actual mechanism is for such cascade in this range.The paradoxes (a) and (b) have been partly addressed by Dematteis & Lvov (2021) and Dematteis et al. (2022), mainly for a modified GM76 spectrum that serves as a stationary solution to (1.1) [action spectrum n(k, m) ∼ k −3.69 m 0 instead of the standard GM76 n(k, m) ∼ k −4 m 0 ].In these works, it was necessary to consider the non-rotating condition f = 0 such that (1.1) becomes scale-invariant and yields a power-law solution.For this modified power-law spectrum, it was identified that ID provides a non-zero and frequencyforward flux by considering the complete diffusion tensor (i.e., including the off-diagonal elements) and that local interactions (which had been ignored in McComas et al.) play a major role in the downscale energy cascade.By applying a combined analytical and numerical approach, the authors explicitly calculated the downscale energy flux, which is within a factor of 2 compared to the prediction by finescale parameterization (1.4).
In spite of the significant progress achieved in Dematteis & Lvov (2021) and Dematteis et al. (2022), the paradoxes (a) and (b), along with a convincing evaluation of downscale energy flux in quantitative consistency with (1.4) in the WKE framework, have never been addressed for the original problem of the GM76 spectrum.This is an important task considering that the GM76 spectrum, although not a stationary state of (1.1) as understood now, is still largely considered as a standard model for realistic IGW spectra among most physical oceanographers.We undertake this task leveraging the fast rise of computational power that has enabled a direct numerical evaluation of the complete spectral energy transfer based on (1.1).Since our approach is purely numerical and we do not seek a scale-invariant field, we are able to naturally incorporate a finite value of f that was not treated in Dematteis & Lvov (2021) and Dematteis et al. (2022).Our results show that the energy flux across the critical vertical scale of 10 m (generally considered as the scale where dissipation starts to take over linear wave dynamics) is approximately 1.5×10 −9 W kg −1 , with a factor of 1.5 greater than the finescale formula (1.4).We further decompose the cascade into different mechanisms and show that the downscale flux is mainly provided by PSI and local interactions, with ID contributing almost zero flux.This is in sharp contrast to McComas et al.'s calculation (flux based on PSI plus ID) and addresses the paradox (b).In addition, we find that there exists a clear frequency-forward cascade, supplied mainly by local interactions and ES, which addresses paradox (a).The role of ES, which was previously hypothesized to have no effect on energy cascade in GM76, is now revealed because of the adoption of finite Coriolis frequency f .

Numerical Method
For numerical evaluation of spectral energy transfer, we use the WKE derived in Lvov & Tabak (2001, 2004) with detailed formulation of the interaction coefficient V provided in Lvov et al. (2010Lvov et al. ( , 2012) ) and Pan et al. (2020).This WKE is derived from a Hamiltonian formulation of the dynamical equation of IGWs, in which the isopycnal vertical coordinate has to be used.The isopycnal vertical wavenumber, m i , is converted from its Eulerian counterpart by m i = (g/ρN 2 ) m. Hereafter, we use notation m throughout the paper for convenience, which represents isopycnal m i in the formulation of WKE [e.g., (2.1) and the appendices] and Eulerian m in presenting the results in §3.
The WKE, in the form of (1.1), involves a collision integral over six dimensions in p 1 and p 2 .One can reduce the dimension of integration by integrating only on the resonant manifold defined by the delta functions.Since GM76 spectrum is horizontally isotropic, it is convenient to first integrate out the horizontal angle dependence, then reduce the integration by applying the delta functions.With detailed formulation provided in Appendix A, this procedure leads to an integration over only two dimensions, namely magnitudes of horizontal wavenumbers k 1 and k where functions h + , h − , g + ′ , g − ′ (which additionally depend on k, m, k 1 , k 2 ) and the roots m * + and m * − are defined in Appendix A, k = k 2 x + k 2 y ∈ R + , and m ∈ R (taking both positive and negative values).The numerical integration of (2.1) is rather straightforward, but care has to be taken in terms of the root finding for m * + and m * − with details discussed in Appendix B. Our numerical code, implemented in FORTRAN with Message Passing Interface (MPI) for parallel computation, is made available on GitHub at https://github.com/yue-cynthia-wu.
Our numerical approach, in principle, shares some level of similarity to "Method 3" in Eden et al. (2019b) regarding the evaluation of the collision integral (2.1) on the resonant manifold, but the latter is implemented for a different version of WKE.Additionally, our procedure (of using cylindrical coordinates and integrating out the horizontal angle dependence) does enforce horizontal isotropy of the IGW spectrum.This feature is beneficial for our planned subsequent work (beyond this paper) to integrate the WKE in time while exactly maintaining the horizontal isotropy as was done in Olbers et al. (2020).In addition, Eden et al. (2019b) employed other methods using broadened delta functions to compute the collision integral, but the results are more noisy and do not show clear advantage.Indeed, as understood recently in both pure mathematical derivation (e.g.Deng & Hani 2023) and numerical studies (e.g.Hrabski & Pan 2020, 2022), the WKE should be considered as a result of maintaining sufficient quasi-resonances from the dynamical equation, so using broadening in the delta function is somewhat redundant.Such broadened delta functions, on the other hand, might be physically meaningful if finite size effect is important, as in situations described in Pan & Yue (2017).
Despite the observed variability in the spectral forms of the realistic IGW fields in different seasons and at different geographical locations, we start with the GM76 model which is one of many realistic possible spectra where m * = 10 −2 m −1 is a reference vertical wavenumber, and the functions The action density spectra, which is needed in evaluating (2.1), can be calculated by n(k, m) = E(k, m)/ω = E(ω, m)(∂ω/∂k)/(2πωk) considering horizontal isotropy.
In the high-ω high-m regime of the spectrum, we have E(ω, m) ∼ ω −2 m −2 and n(k, m) ∼ k −4 m 0 (see Figure 2).In order to set our computation for a physical problem that reflects the size of the real ocean, we consider a domain with a horizontal circular radius of 42.4 km and a vertical extent of 2.1 km.To evaluate (2.1), we use 1080×1080 grids of uniform spacing in both k and m, with the smallest resolved scales to be 40 m and 2 m in the two directions, i.e.,

Energy transfer in spectral space
The energy transfer in spectral space ∂E(k, m)/∂t is calculated by multiplying ∂n(k, m)/∂t with ω, where ∂n(k, m)/∂t is obtained by numerically evaluating the WKE (2.1).In Figure 3, we plot (2πk)(mk) ∂E(k, m)/∂t, with the factor 2πk accounting for horizontal azimuth integration and mk accounting for the plot in log-log axis.More precisely, with these prefactors, the total ∂E/∂t can be conveniently computed by integrating the values over the area in Figure 3 ] (we will later show in Figure 5 that our simulation conserves the total energy so that ∂E/∂t = 0).This plotting technique is also used in Eden et al. (2019a,b), which facilitates an unbiased visualization of energy transfer.We further split Figure 3 into two panels, showing respectively the source [with ∂E(k, m)/∂t < 0, providing energy] and sink [with ∂E(k, m)/∂t > 0, receiving energy] regions.Here, the terminologies "sink" and "source" are inherited from Eden et al. (2019b) and are used to indicate the direction of energy cascade.If a stationary spectrum is assumed, they can be physically related to the generation and dissipation mechanism that has to balance the spectral energy transfer.We see that energy is transferred from waves in the frequency band [2f, 4f ] to both lower and higher † Since the interaction coefficients in the WKE used in this paper are derived under hydrostatic approximation, the results at regions of large k/m should be considered as an approximation.Olbers (1974Olbers ( , 1976) ) and Müller & Olbers (1975) provided a nonhydrostatic version of the WKE for IGWs.It is important to verify that the WKE under the above setting stays in the weakly nonlinear regime and provides a valid approximation of the dynamics.In particular, one may concern about the rapid modal evolution in the high-wavenumber high-frequency regime, as first pointed out in Holloway (1978), which may violate condition (1.3) regarding the normalized Boltzmann rate ε p .For this purpose, we check the values of ε p in the selected scale limits.As shown in Figure 4, there indeed exist large spectral regions where m > 0.1 m −1 and/or ω > 0.9N with |ε p | ∼ O(1), indicating the failure of WKE in the regime of high wavenumbers and/or high frequencies.These regions with |ε p | not much less than 1 will be treated with caution in the subsequent discussion of energy fluxes.
We further define the energy fluxes across k 0 , m 0 and ω 0 , respectively, based on energy conservation in finite domain where k max = 0.016 m −1 , m max = 0.32 m −1 , and 1 an indicator function.Equation (3.1)-(3.3)are energy fluxes only due to nonlinear wave-triad interactions within the selected scale limits, which do not include fluxes entering the IGW field from the largescale end by generation nor draining from the smallscale end by dissipation.The evaluations here are based on the conservation of total energy ∂E/∂t by the WKE.In the prefactors 4πk, 2πk comes from the integration over horizontal azimuth and 2 accounts for the vertical symmetry over ±m.The energy flux in all three directions are plotted in Figure 5 (black curves).We see that P α (α max ) ≈ 0 with α = k, m and ω, indicating energy conservation.We remark here that energy conservation is only approximately achieved by our numerical algorithm since the roots m * + and m * − in (2.1) found by the rootfinding algorithm (Appendix B) do not exactly lie on the discrete m-grid points, which breaks the symmetry when looping over three wavenumber vectors in a triad.This mainly affects the high-frequency regime of the spectrum (Figure 3 with strong sink and source where ω > 0.7N ) †. Since the hydrostatic approximation is also violated in this regime, we discard the contribution of triads with frequencies greater than the cutoff frequency ω cutoff = 0.7N in the calculation of energy fluxes.The fluxes defined in Equation (3.1)-(3.3)are energy transfer only due to nonlinear wave-wave interactions within the selected scale limits, and the values have to be zero at the boundaries (e.g., k min and k max , etc.) due to energy conservation.It is clear from Figure 5 that the GM76 spectrum, as expected, does not yield a constant energy flux in any spectral direction.While the flux across frequency is bi-directional [Figure 3 and 5(c)], the fluxes across horizontal and vertical wavenumbers are mainly in the forward direction (toward small scales).
To evaluate the downscale energy flux P d that provides energy for small-scale dissipation and mixing, i.e., to evaluate the finescale parameterization, we may follow two approaches.The first approach is to consider P d = P m (m c ) with the critical vertical wavenumber m c evaluated through ´mc Polzin et al. 2014).This approximation encapsulates the energy escaping the internal-wave field at the critical wavenumber m c ≈ 0.6 m −1 past which internal waves become unstable to shear instability.The first approach gives P d = P m (m c ) = 1.5 × 10 −9 W kg −1 , with a factor of Figure 5. Energy fluxes due to nonlinear wave-wave interactions within the selected scale limits across (a) horizontal wavenumbers, (b) vertical wavenumber, and (c) frequency for GM76.In (b), dashed and dash-dotted lines denote m cutoff = 0.2 m −1 (corresponding to the smallest scale that 90% of the waves with |εp| < 0.2; Figure 4) and the critical vertical wavenumber mc = 0.6 m −1 , respectively.In (c), solid, dashed and dash-dotted lines denote ω = 2f, 3f and 4f , respectively.Colored curves denote PSI, ES, ID and local interactions (LT) using selection criteria defined in §3.2.
1.5 greater than P finscale = 1.0×10 −9 W kg −1 from (1.4) (equation rescaled for our values of f and N ).However, as seen in Figure 4, the Boltzmann rate close to m c contains large regions of values that are not much less than 1, making the validity of WKE questionable.The second approach is to instead define a cutoff vertical wavenumber m cutoff ≈ 0.2 m −1 , below which only 10% of the computed waves violate the weak nonlinearity assumption with |ε p | > 0.2.The second approach gives P d = P m (m cutoff ) = 1.6 × 10 −9 W kg −1 .This approximation yields an upper bound of energy available for dissipation while (almost) free of uncertainties associated with the first approach.We should acknowledge that there exists a gap between the (vertical) length scales where WKE breaks down and the scale of 10 meter in the vertical.Fortunately, the energy-flux curve between the two scales (corresponding to m c ≈ 0.6 m −1 and m cutoff ≈ 0.2 m −1 ) is relatively insensitive to the vertical wavenumber (Figure 5b), making our estimate quite robust.

Contributions of PSI, ES and ID triads
In this section, we discuss the decomposition of energy transfer into contributions from different mechanisms, namely nonlocal interactions of PSI, ES and ID, and local interactions.The nonlocal interactions exhibit scale separation in frequency or wavenumber, or both, based on which our classification method is designed.For a resonant triad, we rank the frequencies from high to low as (ω H , ω M , ω L ) and the magnitude of vertical wavenumbers as (|m H |, |m M |, |m L |).We can then naturally classify nonlocal interactions according to some threshold values of the two metrics as follows: • PSI: The above criterion characterizes PSI as scale-separated in m and halving in ω, ES as scale-separated in ω and halving in m, and ID as scale-separated in both m and ω.
In practice, we use ξ = η = 2 and ϵ = α = 0.1 for results below, but we note that the major conclusions hold for a reasonable range of parameters selected.We also remark that the above choices of ξ and η are "conservative" for local interactions, in the sense that some interactions with moderate ξ and η (say slightly greater than 2) are classified as nonlocal, instead of local, interactions.This is not a drawback for this paper, as we shall show that even for this "conservtive" choice of local interactions, their role in energy cascade is significant, in sharp contrast to McComas et al.'s early conjecture.In the following sections, we discuss the roles of each mechanism in spectral energy transfer.

The PSI mechanism
Parametric subharmonic instability represents the decay of a low-wavenumber parent wave into two nearly identical high-wavenumber child waves with half frequencies.One unit of action of the parent wave p 0 is transferred into two units of action of the child waves p 1 and p 2 (see Figure 1a).
Using the criterion defined above, we compute the spectral energy transfer due to PSI, with result shown in Figure 6a.In terms of energy cascade in frequency, we see that PSI contributes predominantly to the backward cascade, namely moving energy from frequency band [2f, 4f ] (source) to [f, 2f ] (sink).This frequency cascade is accompanied by a strong forward cascade in vertical wavenumbers, which is also clear from the figure.The physical picture of PSI revealed here is consistent with the existing understanding that PSI is effective in tidal damping (i.e., extracting energy above 2f ) and that PSI contributes significantly to downscale energy cascade in the IGW continuum (McComas & Bretherton 1977).

The ES mechanism
The energy transfer due to ES is plotted in Figure 6b, which shows a clear forward cascade in frequency.This was not understood by previous theory of McComas et al., which instead postulated that ES can be neglected for energy transfer in any vertically symmetric IGW spectrum.In fact, the previous postulation to neglect ES' contribution is a bit surprising given that the dynamics of ES is similar to a diffusion process that can be understood in analogy to ID (the previous researchers do have a better understanding of ID as will be discussed later in the paper).Consider an ES triad p 0 = p 1 + p 2 (Figure 1b), where p 2 is the near-inertial mode, and p 0 and p 1 the high-frequency modes with ω 0 ≈ ω 1 + f .Given the fact that the action spectra are red with respect to ω, i.e., most action contained in low frequencies, it is reasonable to set n 2 ≫ n 0 , n 1 and n 0 < n 1 .According to WKE, we then have where C = 4π|V (p 0 , p 1 , p 2 )| 2 is a constant for this triad.The sign of ∂n/∂t indicates p 0 is a sink while p 1 and p 2 are sources.Consumption of one unit of action of p 2 combined with one unit of action of p 1 results in generation of one unit of action of p 0 .This process can be equivalently described as diffusion from p 1 to p 0 (i.e., toward higher frequency) which in the meanwhile extracts energy from p 2 .
It is clear from the above analysis that a finite value of f is critical to enable the forward cascade in frequency, which we indeed incorporate in our calculation.The assumption of f = 0 used in previous research (either for convenience or for obtaining a scale-invariant WKE as in Dematteis & Lvov 2021;Dematteis et al. 2022) is perhaps one reason leading to the neglect of ES in energy transfer.With the above analysis, we can conclude that forward frequency cascade by ES shown in Figure 6b should essentially be expected, and we reach consistency in theory and numerical results.

The ID mechanism
The energy flux due to ID is plotted in Figure 6c, which shows a very weak transfer compared to those from other interaction mechanisms.One could further expect that ID contributes insignificantly to the downscale energy cascade for GM76 spectrum.This is in strong disagreement with results in McComas & Müller (1981a,b) that ID contributes nearly 30% of the total downscale energy cascade.The result from McComas et al. roots from the hypothesis that GM76 spectrum yields a constant downscale energy flux which relies on a logarithmic correction to the ID flux.It is clear from our analysis (Figure 5) that the constant flux hypothesis is incorrect, and thus the logarithmic correction has no meaningful ground.
The ID mechanism for GM76 or more general IGW spectra can be conveniently understood from a diffusion equation in the high-m high-ω regime: ∂n/∂t = ∂/∂p i (D ij ∂n/∂p j ) with p i = (k x , k y , m) for i = (1, 2, 3) and D ij as the diffusion coefficient matrix.This equation, including the detailed formulation of D ij , can be derived by taking dominant terms in the WKE or from a WKB approximation of the dynamic equation, which was first done in McComas & Müller (1981b) and re-derived by Lvov & Polzin (2022) for IGWs [also see derivations in other physical contexts such as MHD turbulence (Nazarenko et al. 2001), Rossby waves (Connaughton et al. 2015) and surface gravity waves (Korotkevich et al. 2023)].With D 33 being the dominant element in D ij for IGWs, the leading order effect of the diffusion takes place in the vertical wavenumber direction and can be approximated by a one-dimensional diffusion equation ∂n/∂t = ∂/∂m(D 33 ∂n/∂m).Along with this approximation, one can see that in the high-m high-ω regime the direction of diffusion is fully determined by the dependence of wave action spectrum on m, i.e., action diffuses down gradient in vertical wavenumber direction, according to n(k 0 , m) at a given k 0 .
The physical picture of ID revealed above through the diffusion equation can be alternatively explained directly via WKE using a similar argument as our previous one for ES.Take an ID triad p 0 = p 1 + p 2 (Figure 1c) with p 1 the low-m low-ω mode and consider n 1 ≫ n 0 , n 2 , we obtain from WKE that We see that the exchange of action between the two high-m high-ω modes, p 2 and p 0 , depends on the relative magnitudes between n 2 and n 0 .The one that is greater between n 2 and n 0 serves as the source of the diffusion and the other as the sink.If we further assume that p 0 and p 2 have the same horizontal wavenumber k 0 , i.e., p 1 is (almost) vertical, then the diffusion direction is again determined by the dependence of n(k 0 , m) on m.Now let us consider a general power-law spectrum in the high-m high-ω regime: E(ω, m) ∼ ω −2−p m −2−p−q , equivalent to n(k, m) ∼ k −4−p m −q .Based on the above analysis (either from the diffusion equation or WKE), the leading-order diffusion direction is controlled only by the parameter q.For the GM76 spectrum with p = q = 0, it is expected that the leading-order diffusion vanishes, i.e., the GM76 spectrum is indeed approximately a zero-flux state for ID.The energy transfer in Figure 6c, in fact, mainly comes from the (off-diagonal) sub-elements in D ij (other than D 33 ).The effects of the sub-diffusion is analyzed analytically in Dematteis et al. (2022) for the scale-invariant case.For our case with finite f (which breaks the scale invariance), an analytical treatment is generally much difficult.Nevertheless, our direct numerical calculation shows that the transfer generated by the sub-diffusion is weak compared to other interaction mechanisms.
We can further deduce the ID dynamics under conditions q > 0 and q < 0. For q > 0 and q < 0 respectively, the action diffuses toward higher and lower m (with the same k), indicating a backward and a forward cascade in frequency.We can leverage our numerical tools to verify these inferences.For the former with q > 0, we consider a GM75 spectrum with p = 0 and q = 0.5.The energy transfer due to ID is computed as shown in Figure 7a, where we do see a dominating backward cascade in frequency.For the latter, we consider a realistic spectrum E(ω, m) ∼ ω −2 m −1.8 with p = 0 and q = −0.2taken from field measurements in the Southeast Subtropical North Pacific (Polzin & Lvov 2011) reporting E(ω, m) ∼ ω −2 m −1.9∼−1.75 .The ID energy transfer, as shown in Figure 7b, indeed exhibits a dominating forward cascade.In summary, the ID mechanism can lead this assumption, it is likely related to some simple calculation regarding the interaction coefficient V and red IGW action spectra.However, a comprehensive understanding of the problem also relies on other factors, such as the number of triads participating in energy transfer and the specific form of the GM76 action spectrum.Our direct calculation of (1.1) encapsulating all factors clearly demonstrates the paramount importance of local interactions for the GM76 spectrum.The result here also echos those in Dematteis & Lvov (2021) and Dematteis et al. (2022) for a modified GM76 spectrum in the scale-invariant case.

Constituents of energy flux
The energy flux in directions k, m and ω due to each mechanism are plotted in Figure 5.For forward cascade in k, we see that majority of the cascade is provided by local interactions, with ES contributing a relatively small fraction.For forward cascade in m, local interactions and PSI each contribute approximately half to the total flux.The energy transfer in frequency exhibits a bi-directional flux.The backward flux, moving energy from [2f, 4f ] to lower frequencies, is supplied by both PSI and local interactions with similar magnitudes.The forward cascade results from local interactions, ES and ID (mainly from sub-diffusion process) with descending contributions.Among all directional cascades, local interactions are the only mechanism participating significantly in all of them, which was neglected in McComas et al.'s early works.
We are eventually in a good position to state our new understanding regarding paradoxes

Conclusions and discussions
Through direct evaluation of the collision integral in WKE of IGWs, we study the spectral energy transfer for the GM76 spectrum.Our calculation of the downscale energy flux, through its maximum value over all vertical wavenumbers, provides an estimate in close agreement with the finescale parameterization.We also analyse different interaction mechanisms, resolving some long-standing paradoxes in the field.Our new understanding includes: (3) ID can provide cascade toward different directions depending on spectral slopes of the IGW continuum.For GM76, the leading-order flux by ID vanishes, with sub-diffusion process providing a weak forward frequency cascade.(4) The ES mechanism provides a forward frequency cascade (but no cascade in wavenumbers) in vertically symmetric IGW fields, which was not investigated in previous works.Our capability of numerical evaluation of WKE opens a new door to an advanced understanding of oceanic IGW cascade and mixing.Among all possible directions of future study, an immediately fruitful one is to evaluate the flux properties and magnitudes for various IGW spectral forms.As revealed in field measurements (e.g., Polzin & Lvov 2011) and wave turbulence theory (e.g., Lvov et al. 2010 in terms of stationary solutions to the WKE), the oceanic IGW spectrum exhibits large variability, deviating from the GM76 model.Under this circumstance, it is not clear whether the current finescale formula (1.4), developed mainly referencing to the GM76 spectrum, is sufficiently robust for all IGW spectral forms.Our numerical method provides direct approach through which this problem can be studied.In addition, combining our WKE predictions with recent highresolution regional ocean simulations (e.g.Nelson et al. 2020;Thakur et al. 2022;Skitka et al. 2023) may bring new insights to the field.
We also would like to mention that there is an active debate on the relative importance between the nonlinear wave-wave interactions and wave-eddy interactions in IGW energy cascade, and our current paper clearly does not consider the latter.In order to consider both processes, it may be beneficial, as a first-order approximation, to include an additional eddy-diffusion term (Kafiabad et al. 2019;Dong et al. 2020) in the WKE and study the full equation (note that this eddy diffusion is linear assuming stationary eddy field).While these are future directions we plan to consider, the current work perhaps already has shed some light on the problem in terms of understanding the importance of local interactions that is only part of the wave-wave interactions.In the inside integral, for given k, m, k 1 and k 2 , where the subscript m + (m − ) of the interaction coefficient V and quadratic function f p12 (f 1p2 ) denotes projection on m = m 1 + m 2 (m 1 = m + m 2 ).Finally, CL is reduced to an integration over two dimensions k 1 and k 2 by integrating out delta function in frequency where m * + 1 is the root of g + (m 1 ) = 0 for summation interactions, and m * − 1 is the root of g − (m 1 ) = 0 for reduction interactions.The denominators g + ′ (m 1 ) and g − ′ (m 1 ) are the m 1 -derivatives of functions g + (m 1 ) and g − (m 1 ), respectively [also see Eden et al. 2019b, eq. (8)] We note that (A 4) and (A 6) involve a singularity (first-order pole) when p, p 1 and p 2 lie on the same vertical plane, i.e., k = k 1 + k 2 or k 1 = k + k 2 , leading to S △ = 0 and corresponding to the collinear triad interactions described in Dematteis & Lvov (2021).However, this is an integrable singularity since it is not involved in (A 1) before the coordinate transform and angle integration.Therefore, given sufficiently fine grid resolution of k 1 and k 2 , the singularity point can be neglected in the integration of (A 8).For finite grid resolution, we can approximately retrieve the contribution from the nearby region of the singularity point using (A 1) (and its counterpart for reduction interactions).Specifically, for given k and m, contributions from vicinity of all singularities in the integral of (A 8) can be accounted for by treating (A 1) with the following procedures: set (due to isotropy) k = (k, 0); integrate out the delta function in k with respect to the integration over k 2x and k 2y ; consider k 1 in the same direction as k by setting k 1y = 0 and convert the integration along k 1y into the integral multiplying by the grid size ∆k; change the dummy variable k 1x into k 1 ; integrate out the ω and m delta functions with respect to the integration over m 1 and m 2 .These procedures lead to

Figure 2 .
Figure 2. GM76 model of the (left) energy and (right) action density spectra, with E(ω, m) ∼ ω −2 m −2 and n(k, m) ∼ k −4 m 0 in the high-frequency, high-wavenumber regime.Curves in each figure represent the spectrum with one variable taking fixed values.Red vertical lines denote the GM76-referenced vertical wavenumber m * = 10 −2 m −1 .
which bounds the IGW frequency between f and N , while the hydrostatic version leads to (non-physical) super-buoyancy IGWs in large spectral area of the selected (k, m)-domain where k/m is large †.With above setting, our numerical calculations are performed on the Great Lakes clusters at University of Michigan with 2 nodes of 72 CPUs, and the simulation takes 6-8 hours to calculate all O(10 12 ) triad interactions.

Figure 3 .
Figure 3. Energy transfer log[(2πk)(mk)∂E/∂t] for GM76 with sinks (∂E/∂t > 0) on the left and sources (∂E/∂t < 0) on the right.We only plot the range with positive m since the spectrum at negative m is completely symmetric.The white solid, dashed, dash-dotted and dotted lines along the diagonal denote frequencies 2f , 3f , 4f and 35f (= 0.7N ), respectively.

Figure 6 .
Figure 6.As Figure 3 but for different mechanisms.
(a) and (b).For (a), we now understand that the frequency cascade is bi-directional, with the forward flux provided by local interactions, ES and ID, all elements ignored in McComas et al.'s works.For (b), the ID mechanism in GM76 is indeed approximately a zero-flux state, except forming a weak forward cascade in frequency through the subdiffusion process.McComas et al.'s argument about ID providing significant portion of downscale flux should be replaced by local interactions.
(1) Local interactions are important for energy cascade in all spectral directions, which were completely neglected in early works by McComas et al.. (2) The downscale energy flux (toward high vertical wavenumbers) is supplied by PSI and local interactions, rather than PSI and ID as understood in McComas et al.'s works.

Figure 8 .
Figure 8. Graphical interpretation on the root finding of m * 1+ and m * 1− for an example case given k, k1, k2 and m.