3.1. Energy transfer in spectral space

The energy transfer in spectral space $\partial E(k,m)/\partial t$ is calculated by multiplying $\partial n(k,m)/\partial t$ with $\omega$, where $\partial n(k,m)/\partial t$ is obtained by numerically evaluating the WKE (2.1). In figure 3, we plot $(2{\rm \pi} k)(mk)\,\partial E(k,m)/\partial t$, with the factor $2{\rm \pi} k$ accounting for horizontal azimuth integration, and $mk$ accounting for the plot in the log-log axis. More precisely, with these prefactors, the total $\partial E/\partial t$ can be computed conveniently by integrating the values over the area in figure 3, i.e. $\partial E/\partial t = \iint (2{\rm \pi} k)(mk)\,(\partial E(k,m)/\partial t) \,\textrm {d}[\log (m)]\,\textrm {d}[\log (k)]$. (We will show later, in figure 5, that our simulation conserves the total energy so that $\partial E/\partial t=0$.) This plotting technique is also used in Eden, Chouksey & Olbers (Reference Eden, Chouksey and Olbers2019a) and Eden et al. (Reference Eden, Pollmann and Olbers2019b), which facilitates an unbiased visualization of energy transfer. Figure 3(a) shows the source regions (with $\partial E(k,m)/\partial t<0$, providing energy), and figure 3(a) shows the sink (with $\partial E(k,m)/\partial t>0$, receiving energy) regions. Here, the terminologies ‘sink’ and ‘source’ are inherited from Eden et al. (Reference Eden, Pollmann and Olbers2019b) and are used to indicate the direction of energy cascade. If a stationary spectrum is assumed, then they can be related physically 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 frequencies. The component of forward frequency cascade is also seen in the results of Eden et al. (Reference Eden, Pollmann and Olbers2019b). We will further investigate the mechanisms leading to this cascade in § 3.2, addressing paradox (a) in the theory of McComas et al.

Figure 3. Energy transfer $\log [(2{\rm \pi} k)(mk)\,\partial E/\partial t]$ for GM76 with (a) sinks ($\partial E/\partial t>0$), and (b) sources ($\partial E/\partial t<0$). We plot only 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.

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 be concerned about the rapid modal evolution in the high-wavenumber, high-frequency regime, as first pointed out in Holloway (Reference Holloway1978), which may violate condition (1.3) regarding the normalized Boltzmann rate $\varepsilon _{\boldsymbol {p}}$. For this purpose, we check the values of $\varepsilon _{\boldsymbol {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 $\omega >0.9N$ with $|\varepsilon _{\boldsymbol {p}}| \sim O(1)$, indicating the failure of the WKE in the regime of high wavenumbers and/or high frequencies. These regions with $|\varepsilon _{\boldsymbol {p}}|$ not much less than 1 will be treated with caution in the subsequent discussion of energy fluxes.

Figure 4. Normalized Boltzmann rate $\varepsilon _{\boldsymbol {p}}=2{\rm \pi} (\partial n_{\boldsymbol {p}}/\partial t)/(\omega _{\boldsymbol {p}} n_{\boldsymbol {p}})$ for GM76 in the $(k,m)$ domain. The white solid, dashed, dash-dotted and dotted lines along the diagonal denote frequencies $2f$, $3f$, $4f$ and $45f$ ($=0.9N$), respectively. Regions where $m>0.1$ m$^{-1}$ and/or $\omega >0.9N$ indicate violation of the weak nonlinearity assumption in wave turbulent theory, including regions of $\varepsilon _{\boldsymbol {p}}>1$ (white) and $\varepsilon _{\boldsymbol {p}}<-1$ (grey) with obvious violation.

We further define the energy fluxes across $k_0$, $m_0$ and $\omega _0$, respectively, based on energy conservation in the finite domain:

(3.1)\begin{gather} \mathcal{P}^k(k_0)={-}\int_0^{m_{max}} \int_0^{k_0} 4{\rm \pi} k\,\frac{\partial E(k,m)}{\partial t}\,\textrm{d} k \,\textrm{d} m, \end{gather}

(3.2)\begin{gather}\mathcal{P}^m(m_0)={-}\int_0^{m_0} \int_0^{k_{max}} 4{\rm \pi} k\,\frac{\partial E(k,m)}{\partial t}\,\textrm{d} k \,\textrm{d} m, \end{gather}

(3.3)\begin{gather}\mathcal{P}^\omega(\omega_0)={-}\int_{0}^{m_{max}} \int_0^{k_{max}} 4{\rm \pi} k\,\frac{\partial E(k,m)}{\partial t}\,\mathbb{1}_{\omega\le\omega_0}\,\textrm{d} k \,\textrm{d} m, \end{gather}

where $k_{max}=0.016$ m$^{-1}$, $m_{max}=0.32$ m$^{-1}$, and $\mathbb{1}$ is an indicator function. Equations (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 large-scale end by generation nor draining from the small-scale end by dissipation. The evaluations here are based on the conservation of total energy $\partial E/\partial t$ by the WKE. In the prefactors $4{\rm \pi} k$, $2{\rm \pi} k$ comes from the integration over horizontal azimuth, and $2$ accounts for the vertical symmetry over ${\pm }m$. The energy fluxes in all three directions are plotted in figure 5 (black curves). We see that $\mathcal {P}^\alpha (\alpha _{max}) \approx 0$ with $\alpha =k$, $m$ and $\omega$, indicating energy conservation. We remark here that energy conservation is achieved only approximately by our numerical algorithm since the roots $m^{*+}$ and $m^{*-}$ in (2.1) found by the root-finding algorithm (Appendix B) do not lie exactly on the discrete $m$-grid points, which breaks the symmetry when looping over three wavenumber vectors in a triad. This affects mainly the high-frequency regime of the spectrum (figure 3 with strong sink and source where $\omega >0.7N$). (The region where $\omega \sim N$ is problematic for two reasons: (1) if $m$ is not much greater than $k$, then the hydrostatic approximation in the WKE breaks down, making the evaluation of $\partial n/\partial t$ (and $\partial E/\partial t$) untrustworthy; (2) if $m \in \mathbb {R}^+$ is small, then the root finding for the reduction interactions ($m^*_{1,{right}} \in (0,m)$) deviates from its true value due to limited resolution between 0 and $m$; see Appendix B and figure 8.) Since the hydrostatic approximation is also violated in this regime, we discard the contribution of triads with frequencies greater than the cutoff frequency $\omega _{cutoff}=0.7N$ in the calculation of energy fluxes.

Figure 5. Energy fluxes due to nonlinear wave–wave interactions within the selected scale limits across (a) horizontal wavenumber, (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 for 90 % of the waves with $|\varepsilon _{\boldsymbol {p}}|<0.2$; figure 4) and the critical vertical wavenumber $m_c=0.6$ m$^{-1}$, respectively. In (c), solid, dashed and dash-dotted lines denote $\omega =2f$, $3f$ and $4f$, respectively. Coloured curves denote PSI, ES, ID and local interactions (LT) using selection criteria defined in § 3.2.

The fluxes defined in (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 (figures 3 and 5c), the fluxes across horizontal and vertical wavenumbers are mainly in the forward direction (toward small scales).

To evaluate the downscale energy flux $\mathcal {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 $\mathcal {P}_d=\mathcal {P}^m(m_c)$ with the critical vertical wavenumber $m_c$ evaluated through $\int _0^{m_c} 2m^2\,E_k(m)\,\textrm {d} m=2{\rm \pi} N^2/10$ (Polzin et al. Reference Polzin, Naveira Garabato, Huussen, Sloyan and Waterman2014). This approximation encapsulates the energy escaping the internal-wave field at the critical wavenumber $m_c\approx 0.6$ m$^{-1}$ past which internal waves become unstable to shear instability. The first approach gives $\mathcal {P}_d=\mathcal {P}^m(m_c)=1.5 \times 10^{-9}$ W kg$^{-1}$, with a factor 1.5 greater than $P_{finescale}=1.0 \times 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 the WKE questionable. The second approach is to instead define a cutoff vertical wavenumber $m_{cutoff}\approx 0.2$ m$^{-1}$, below which only 10 % of the computed waves violate the weak nonlinearity assumption with $|\varepsilon _{\boldsymbol {p}}|>0.2$. The second approach gives $\mathcal {P}_d=\mathcal {P}^m(m_{cutoff})=1.6 \times 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 the WKE breaks down and the scale of 10 m in the vertical. Fortunately, the energy-flux curve between the two scales (corresponding to $m_c\approx 0.6$ m$^{-1}$ and $m_{cutoff}\approx 0.2$ m$^{-1}$) is relatively insensitive to the vertical wavenumber (figure 5b), making our estimate quite robust.

3.2. Contributions of PSI, ES and ID triads

In this subsection, we discuss the decomposition of energy transfer into contributions from different mechanisms, namely non-local interactions of PSI, ES and ID, and local interactions. The non-local 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 $(\omega ^H,\omega ^M,\omega ^L)$, and the magnitudes of vertical wavenumbers as $(|m^H|,|m^M|,|m^L|)$. We can then classify non-local interactions naturally according to some threshold values of the two metrics, as follows.

1. (i) PSI: $|m^M|/|m^L|>\eta$, $1/2 \le \omega ^M/\omega ^H <1/2+\epsilon /2$.

2. (ii) ES: $\omega ^M/\omega ^L>\xi$, $1/2 \le |m^M|/|m^H| <1/2+\alpha /2$.

3. (iii) ID: $\omega ^M/\omega ^L>\xi$, $|m^M|/|m^L|>\eta$.

The above criterion characterizes PSI as scale-separated in $m$ and halving in $\omega$, ES as scale-separated in $\omega$ and halving in $m$, and ID as scale-separated in both $m$ and $\omega$. In practice, we use $\xi =\eta =2$ and $\epsilon =\alpha =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 $\xi$ and $\eta$ are ‘conservative’ for local interactions, in the sense that some interactions with moderate $\xi$ and $\eta$ (say slightly greater than 2) are classified as non-local, instead of local, interactions. This is not a drawback for this paper, as we will show that even for this ‘conservative’ choice of local interactions, their role in energy cascade is significant, in sharp contrast to the early conjecture of McComas et al. In the following subsubsections, we discuss the roles of each mechanism in spectral energy transfer.

3.2.1. 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 $\boldsymbol {p}_0$ is transferred into two units of action of the child waves $\boldsymbol {p}_1$ and $\boldsymbol {p}_2$ (see figure 1a).

Using the criterion defined above, we compute the spectral energy transfer due to PSI, with results shown in figure 6(a). 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 Reference McComas and Bretherton1977).

Figure 6. As figure 3 but for different mechanisms: (a) PSI, (b) ES, (c) ID, (d) LT.

3.2.2. The ES mechanism

The energy transfer due to ES is plotted in figure 6(b), 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 the ES contribution is a little 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 $\boldsymbol {p}_0=\boldsymbol {p}_1+\boldsymbol {p}_2$ (figure 1b), where $\boldsymbol {p}_2$ is the near-inertial mode, and $\boldsymbol {p}_0$ and $\boldsymbol {p}_1$ are the high-frequency modes with $\omega _0 \approx \omega _1 + f$. Given the fact that the action spectra are red with respect to $\omega$, i.e. most action contained in low frequencies, it is reasonable to set $n_2\gg n_0,n_1$ and $n_0< n_1$. According to the WKE, we then have

(3.4)\begin{equation} \left.\begin{gathered} \partial n_0/\partial t= C f_{012} = C(n_1n_2 -n_0n_1-n_0n_2) \approx Cn_2(n_1-n_0) >0, \\ \partial n_1/\partial t={-}C f_{012} ={-}C(n_1n_2 -n_0n_1-n_0n_2) \approx{-}Cn_2(n_1-n_0) <0, \\ \partial n_2/\partial t={-}C f_{021} ={-}C(n_1n_2 -n_0n_1-n_0n_2) \approx{-}Cn_2(n_1-n_0) <0, \end{gathered}\right\} \end{equation}

where $C=4{\rm \pi} \,|V(\boldsymbol {p}_0,\boldsymbol {p}_1,\boldsymbol {p}_2)|^2$ is a constant for this triad. The sign of $\partial n/\partial t$ indicates that $\boldsymbol {p}_0$ is a sink, while $\boldsymbol {p}_1$ and $\boldsymbol {p}_2$ are sources. Consumption of one unit of action of $\boldsymbol {p}_2$ combined with one unit of action of $\boldsymbol {p}_1$ results in generation of one unit of action of $\boldsymbol {p}_0$. This process can be described equivalently as diffusion from $\boldsymbol {p}_1$ to $\boldsymbol {p}_0$ (i.e. towards higher frequency), which in the meanwhile extracts energy from $\boldsymbol {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 $f=0$ used in previous research (either for convenience or for obtaining a scale-invariant WKE, as in Dematteis & Lvov Reference Dematteis and Lvov2021; Dematteis et al. Reference Dematteis, Polzin and Lvov2022) is perhaps one reason leading to the neglect of ES in energy transfer. With the above analysis, we can conclude that the forward frequency cascade by ES shown in figure 6(b) should essentially be expected, and we reach consistency in theory and numerical results.

3.2.3. The ID mechanism

The energy flux due to ID is plotted in figure 6(c), 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 the GM76 spectrum. This is in strong disagreement with results in McComas & Müller (Reference McComas and Müller1981a,Reference McComas and Müllerb) that ID contributes nearly 30 % of the total downscale energy cascade. The result from McComas et al. originates from the hypothesis that the GM76 spectrum yields a constant downscale energy flux that relies on a logarithmic correction to the ID flux. It is clear from our analysis (figure 5) that the constant flux hypothesis is incorrect, thus the logarithmic correction has no meaningful ground.

The ID mechanism for GM76 or more general IGW spectra can be understood conveniently from a diffusion equation in the high-$m$, high-$\omega$ regime: $\partial n/\partial t=$ $(\partial /\partial p_i)(\mathcal {D}_{ij}\,\partial n/\partial p_j)$, with $p_i=(k_x, k_y, m)$ for $i=1,2,3$, and $\mathcal {D}_{ij}$ as the diffusion coefficient matrix. This equation, including the detailed formulation of $\mathcal {D}_{ij}$, can be derived by taking dominant terms in the WKE or from a Wentzel–Kramers–Brillouin (WKB) approximation of the dynamic equation, which was first done in McComas & Müller (Reference McComas and Müller1981b) and re-derived by Lvov & Polzin (Reference Lvov and Polzin2022) for IGWs (see also derivations in other physical contexts, such as magnetohydrodynamic (MHD) turbulence (Nazarenko, Newell & Galtier Reference Nazarenko, Newell and Galtier2001), Rossby waves (Connaughton, Nazarenko & Quinn Reference Connaughton, Nazarenko and Quinn2015) and surface gravity waves (Korotkevich et al. Reference Korotkevich, Nazarenko, Pan and Shatah2023)). With $D_{33}$ being the dominant element in $\mathcal {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 $\partial n/\partial t=(\partial /\partial m) (\mathcal {D}_{33}\,\partial n/\partial m)$. Along with this approximation, one can see that in the high-$m$, high-$\omega$ regime, the direction of diffusion is determined fully by the dependence of the wave action spectrum on $m$, i.e. action diffuses down gradient in the vertical wavenumber direction, according to $n(k_0, m)$ at a given $k_0$.

Alternatively, the physical picture of ID revealed above through the diffusion equation can be explained directly via the WKE using an argument similar to our previous one for ES. Take an ID triad $\boldsymbol {p}_0=\boldsymbol {p}_1+\boldsymbol {p}_2$ (figure 1c), with $\boldsymbol {p}_1$ the low-$m$, low-$\omega$ mode, and consider $n_1\gg n_0,n_2$. Then we obtain from WKE that

(3.5)\begin{equation} \left.\begin{gathered} \partial n_0/\partial t= C f_{012}= C(n_1n_2 -n_0n_1-n_0n_2) \approx C n_1(n_2-n_0), \\ \partial n_1/\partial t={-}C f_{012}={-}C(n_1n_2 +n_0n_1+n_0n_2) \approx{-}C n_1(n_2-n_0), \\ \partial n_2/\partial t={-}C f_{021}={-}C(n_1n_2 +n_0n_1+n_0n_2) \approx{-}C n_1(n_2-n_0). \end{gathered}\right\} \end{equation}

We see that the exchange of action between the two high-$m$, high-$\omega$ modes, $\boldsymbol {p}_2$ and $\boldsymbol {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 $\boldsymbol {p}_0$ and $\boldsymbol {p}_2$ have the same horizontal wavenumber $k_0$, i.e. $\boldsymbol {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-$\omega$ regime: $E(\omega,m)\sim \omega ^{-2-p} m^{-2-p-q}$, equivalent to $n(k,m)\sim k^{-4-p} m^{-q}$. Based on the above analysis (from either the diffusion equation or the 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 6(c), in fact, comes mainly from the (off-diagonal) sub-elements in $\mathcal {D}_{ij}$ (other than $\mathcal {D}_{33}$). The effects of the sub-diffusion is analysed analytically in Dematteis et al. (Reference Dematteis, Polzin and Lvov2022) 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 towards 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 7(a), where we do see a dominating backward cascade in frequency. For the latter, we consider a realistic spectrum $E(\omega,m)\sim \omega ^{-2}m^{-1.8}$ with $p=0$ and $q=-0.2$ taken from field measurements in the Southeast Subtropical North Pacific (Polzin & Lvov Reference Polzin and Lvov2011) reporting $E(\omega,m)\sim \omega ^{-2}m^{-1.9}$ to $E(\omega,m)\sim \omega ^{-2}m^{-1.75}$. The ID energy transfer, as shown in figure 7(b), indeed exhibits a dominating forward cascade. In summary, the ID mechanism can lead to different directions of energy cascade depending on the spectral slopes of the IGW continuum and should be understood with respect to the specific spectrum of interest.

Figure 7. As figure 3 but for the ID mechanism: (a,b) GM75 spectrum with $E(\omega,m)\sim \omega ^{-2} m^{-2.5}$, and (c,d) a modified GM76 spectrum with $E(\omega,m)\sim \omega ^{-2} m^{-1.8}$, for (a,c) sinks, and (b,d) sources.

3.3. Constituents of energy flux

The energy fluxes in directions $k$, $m$ and $\omega$ 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 the 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 the early works of McComas et al.

We are finally in a good position to state our new understanding regarding paradoxes (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 the works of McComas et al. For (b), the ID mechanism in GM76 is indeed approximately a zero-flux state, except that it forms a weak forward cascade in frequency through the sub-diffusion process. The McComas et al. argument about ID providing a significant portion of the downscale flux should be replaced by local interactions.