2.1. Internal waves and the GM spectrum

Garrett and Munk have observed that the internal-wave spectrum is separable in frequency-vertical wavenumbers. In other words, it can be accurately represented as a product of a function of frequency and a function of vertical wavenumber. The GM energy spectrum is therefore represented in the two-dimensional domain of vertical wavenumber $m$ with $m\in [m_{min},m_{max}]$ and frequency $\sigma$ with $\sigma \in [f,N]$ as

(2.1)\begin{gather} e_\textrm{GM}(m,\sigma) = N_0 N b^2 E A(m) B(\sigma), \end{gather}

(2.2a,b)\begin{gather}A(m)=\frac{2}{\rm \pi}\frac{m_\star}{m_\star^2+m^2},\quad B(\sigma) =\frac{2f}{\rm \pi}\frac{1}{\sigma\sqrt{\sigma^2-f^2}},\quad m_\star{=} \frac{3{\rm \pi}}{b}\frac{N}{N_0}, \end{gather}

normalized in such a way that $\int _{m_{min}}^{m_{max}}A(m)\,\textrm {d} m\simeq \int _0^\infty A(m)\,\textrm {d} m =1$, $\int _{f}^{N}B(\sigma )\,\textrm {d}\sigma \simeq \int _0^\infty B(\sigma )\,\textrm {d}\sigma =1$, and the total energy density (per unit mass) is therefore given by $N_0N b^2 E$, in units of $\textrm {J}\ \textrm {kg}^{-1}$. Here $N$ and $N_0=0.00524\ \textrm {s}^{-1}$ are, respectively, the buoyancy frequency and the reference buoyancy frequency, $f=2\times 7.3\times 10^{-5} \sin (l)\ \textrm {s}^{-1}$ is the Coriolis parameter computed at latitude $l=32.5^\circ$, $b=1300\ \textrm {m}$ is the scale height of the ocean, $E=6.3\times 10^{-5}$ is the GM specification of the non-dimensional energy level and $m_\star$ is a reference vertical wavenumber. Furthermore, $m_{min} =2{\rm \pi} (2600\ \text {m})^{-1}$, $m_{max} =2{\rm \pi} (10\ \text {m})^{-1}$ are the physical cutoffs imposed by the ocean depth and by wave breaking, respectively.

2.2. Wave-turbulence interpretation of the GM spectrum

Despite the GM far reaching combination of simplicity and descriptive power of available field measurements, its phenomenological nature does not necessarily provide an explanation to the underlying physics. Since the 1970s, the concept of nonlinear interactions has become the leitmotiv in the search for a physical interpretation of the GM spectrum starting from the primitive equations of a stratified ocean (Olbers Reference Olbers1973, Reference Olbers1976; McComas & Bretherton Reference McComas and Bretherton1977; Pelinovsky & Raevsky Reference Pelinovsky and Raevsky1977; Voronovich Reference Voronovich1979; McComas & Müller Reference McComas and Müller1981; Holloway et al. Reference Holloway, Müller, Henyey and Pomphrey1986; Müller et al. Reference Müller, Holloway, Henyey and Pomphrey1986; Caillol & Zeitlin Reference Caillol and Zeitlin2000; Lvov & Tabak Reference Lvov and Tabak2001). The quadratic nonlinearity in the primitive fluid equations and a dispersion relation allowing for three-wave interactions imply that internal waves interact through triads. In a weakly nonlinear regime, three-wave resonant interactions are responsible for slow, net energy transfers between different wavenumbers (Davis et al. Reference Davis, Jamin, Deleuze, Joubaud and Dauxois2020). This process can be described by a wave kinetic equation, the evolution equation of the action spectrum of the internal wave field (Hasselmann Reference Hasselmann1966; Zakharov, L'vov & Falkovich Reference Zakharov, L'vov and Falkovich1992; Nazarenko Reference Nazarenko2011). In the present paper, we use the three-dimensional wavenumber domain $\boldsymbol {p}=(\boldsymbol {k},m)$, where $\boldsymbol {k}$ and $m$ are the horizontal and the vertical wavenumbers, respectively. Note that $\boldsymbol {k}$ is a two-dimensional horizontal wave vector and we define its norm as $k:=|\boldsymbol {k}|$. The dispersion relation of internal gravity waves is given by $\sigma ^2_\textbf {p} = f^2 + N^2({k^2}/{m^2})$, which can be used to switch from one domain to the other, since only two of the three variables $k$, $m$ and $\sigma$ are independent. The action spectrum $n(\boldsymbol{k},m)$ is related to the energy spectrum $e(\boldsymbol{k},m)$ (now both intended as three-dimensional spectra) via $e(\boldsymbol {k},m)=\sigma n(\boldsymbol {k},m)$, where for simplicity we use the quantities in brackets to specify the domain of dependence of the quantity of interest. We assume horizontal and vertical isotropy, so that we have $e(m,\sigma )=4{\rm \pi} k e(\boldsymbol {k},m)({\textrm {d} \sigma }/{\textrm {d} k})^{-1}$, after integrating over the horizontal azimuthal angle and considering a positive definite m. Considering the scale-invariant (or non-rotating) limit $f\ll \sigma \ll N$, which yields the scale-invariant dispersion relation (here defined as the positive branch) $\sigma = N k/|m|$, (2.1) transforms into

(2.3)\begin{equation} n_\textrm{GM}(\boldsymbol{k},m) = \frac{1}{{\rm \pi}^3} \frac{E b^2 N_0\, f m_\star}{N} k^{{-}4}, \end{equation}

which represents the non-rotating limit of the GM three-dimensional action spectrum, in the horizontal wavenumber–vertical wavenumber domain.

The wave-turbulence theory for internal waves was revisited with the generalized random phase and amplitude formalism (Choi, Lvov & Nazarenko Reference Choi, Lvov and Nazarenko2004, Reference Choi, Lvov and Nazarenko2005; Nazarenko Reference Nazarenko2011) in the series of works (Lvov & Tabak Reference Lvov and Tabak2001; Lvov, Polzin & Tabak Reference Lvov, Polzin and Tabak2004; Lvov et al. Reference Lvov, Tabak, Polzin and Yokoyama2010), where the internal-wave kinetic equation was derived starting from the primitive equations of motion in hydrostatic balance by using isopycnal coordinates. A detailed review is found in the introductory paragraphs of Lvov et al. (Reference Lvov, Tabak, Polzin and Yokoyama2010) and will not be repeated here. The main steps can be schematized as follows: (i) the primitive equations of a vertically stratified ocean in hydrostatic balance and with no background rotation are rewritten in isopycnal coordinates under the Boussinesq approximation. The scale-invariant limit of the dispersion relation in the new variables reads (with no background rotation)

(2.4)\begin{equation} \sigma(\boldsymbol{k},m) =\frac{g}{\rho_0 N} \frac{k}{|m|}, \end{equation}

where $g$ is the acceleration of gravity, $\rho _0$ is the reference density and the vertical wavenumber $m$ is now an inverse density. (ii) In the isopycnal formulation, the equations of motion are reduced to Hamiltonian form for the two conjugate fields $\phi$ and $\varPi$, the velocity potential and the normalized differential layer thickness. (iii) The machinery of wave turbulence is applied by switching to Fourier space and introducing the complex canonical normal variables $c_{\boldsymbol {p}}$ and $c_{-\boldsymbol {p}}^*$, representing complex amplitudes of the normal modes of the system. Under the assumption of spatial homogeneity, the action spectral density is defined as

(2.5)\begin{equation} \langle c_{\boldsymbol{p}_1} c_{\boldsymbol{p}_2}^* \rangle = n_{\boldsymbol{p}_1}\delta_{\boldsymbol{p}_1-\boldsymbol{p}_2}, \end{equation}

where $\delta (\cdot )$ is a Dirac delta, and the angular brackets denote averaging on a suitably defined statistical ensemble: under the standard assumptions of random phases and amplitudes (Zakharov et al. Reference Zakharov, L'vov and Falkovich1992; Nazarenko Reference Nazarenko2011), in the joint limit of large box and small nonlinearity the following wave kinetic equation is derived, assuming isotropy in the horizontal plane (for simplicity, here written in the non-rotating limit):

(2.6)\begin{align} \frac{\partial n_{\boldsymbol{p}}}{\partial t} &= \frac{8{\rm \pi}}{k}\int \left(f^{\boldsymbol{p}}_{12} |V^{\boldsymbol{p}}_{12}|^2\delta_{m-m_1-m_2} \delta_{\sigma_{\boldsymbol{p}}-\sigma_1-\sigma_2} \frac{k k_1 k_2}{\varDelta_{\boldsymbol{p}12}} - (0\leftrightarrow1) - (0\leftrightarrow2)\right) \nonumber\\ &\quad \times \textrm{d} k_1\,\textrm{d} k_2 \,\textrm{d} m_1\,\textrm{d} m_2, \end{align}

where $n_{\boldsymbol {p}}=n(\boldsymbol {k},m;t)$ is the three-dimensional action spectrum defined in (2.5), $f^{\boldsymbol {p}}_{12}=n_1n_2-n_{\boldsymbol {p}}(n_1+n_2)$, $V^{\boldsymbol {p}}_{12}$ is the matrix element describing the magnitude of nonlinear interactions between the triad of wavenumbers $\boldsymbol {p}$, $\boldsymbol {p}_1$ and $\boldsymbol {p}_2$, given below by (3.3). Furthermore, the two delta functions impose the conservation of vertical momentum and energy in each three-wave interaction. The $\varDelta _{\boldsymbol {p}12}$, given by (3.2), is a factor coming from integration of the horizontal momentum delta function, proportional to the area of the triangle with sides $k$, $k_1$ and $k_2$. (iv) The wave kinetic equation (2.6) with dispersion relation (2.4) is fully scale-invariant and the consequent theory of power-law spectra (Zakharov et al. Reference Zakharov, L'vov and Falkovich1992) was worked out in Lvov et al. (Reference Lvov, Tabak, Polzin and Yokoyama2010). Assuming a solution of type

(2.7)\begin{equation} n_{\boldsymbol{p}}\propto k^{{-}a} |m|^{{-}b}, \end{equation}

the stationary solution corresponding to constant energy flux, i.e. the Kolmogorov– Zakharov (KZ) spectrum, can be derived by Zakharov–Kuznetsov conformal mapping (Zakharov et al. Reference Zakharov, L'vov and Falkovich1992) yielding $a=7/2$, $b=1/2$:

(2.8)\begin{equation} n^{PR}(k,m) = k^{-({7}/{2})}m^{-({1}/{2})}. \end{equation}

Such a solution was derived in Pelinovsky & Raevsky (Reference Pelinovsky and Raevsky1977) and again in Lvov & Tabak (Reference Lvov and Tabak2001) and is known as the Pelinovski–Raevski (PR) spectrum. (v) A KZ spectrum is a valid solution of the wave kinetic equation if and only if the locality conditions are satisfied, i.e. when the collision integral on the right-hand side of the wave kinetic equation converges. It turns out that this is not the case for the PR spectrum; more precisely, in the $a-b$ power-law space Lvov et al. (Reference Lvov, Tabak, Polzin and Yokoyama2010) found that the collision integral converges only on the segment $b=0, 3.5 < a < 4$. (vi) On this convergence segment, it was shown by direct numerical integration that the collision integral is zero for $a\simeq 3.7$, locating the scale-invariant stationary solution of the wave kinetic equation at the point $a=3.7, b=0$. Since this is not far from the $a=4, b=0$ point of (2.3), such a solution has therefore been put forward as the possible theoretical explanation of the GM spectrum provided by wave turbulence.

In the present work. we use the wave-turbulence kinetic equation to analyse how the stationary scale-invariant internal-wave spectrum is formed, and we calculate the corresponding energy fluxes related to this spectrum. This quantity is modelled phenomenologically, as interpretation of the available data, by what is known as the finescale parametrization of the oceanic turbulent mixing (Polzin et al. Reference Polzin, Toole and Schmitt1995, Reference Polzin, Garabato, Huussen, Sloyan and Waterman2014; Whalen, Talley & MacKinnon Reference Whalen, Talley and MacKinnon2012; MacKinnon et al. Reference MacKinnon2017; Liang et al. Reference Liang, Shang, Qi, Chen and Yu2018), and represents a fundamental building block of the GCM.

3.1. Locality conditions

Let us consider the wave kinetic equation of internal gravity waves in a non-rotating frame, in hydrostatic balance, and in the scale-invariant limit. This is described by (2.6), with the dispersion relation (2.4), expressed in isopycnal coordinates. By integrating analytically the two remaining Dirac deltas, we simplify the collision integral, reducing it to a double integral. The wave kinetic equation thus takes the following form:

(3.1)\begin{equation} \left.\begin{gathered} \partial_t n_{\boldsymbol{p}} = \mathcal{I}(k,m;a,b):= \int_0^\infty \textrm{d} k_{1}\,\textrm{d} k_2 {\mathcal{J}}(k,k_1,k_2,m),\\ {\mathcal{J}}(k,k_1,k_2,m)= \frac{8{\rm \pi}}{k}(R^{\boldsymbol{p}}_{12}\, f^{\boldsymbol{p}}_{12} - R^1_{\boldsymbol{p}2}\,f^1_{\boldsymbol{p}2} - R^2_{\boldsymbol{p}1} \,f^2_{\boldsymbol{p}1}),\\ R^{\boldsymbol{p}}_{12} = k k_1 k_2 |V^{\boldsymbol{p}}_{12}|^2/(|{g^{\boldsymbol{p}}_{12}}'|\varDelta_{\boldsymbol{p}12}). \end{gathered}\right\} \end{equation}

Here $f^{\boldsymbol {p}}_{12} = n_1n_2 - n_{\boldsymbol {p}} (n_1 + n_2)$ and the area of the triangle of sides $k,k_1,k_2$, coming from integration over angles under the assumption of isotropy, is given by

(3.2)\begin{equation} \varDelta_{\boldsymbol{p}12} = \tfrac12 \sqrt{2(k^2 k_1^2 + k^2 k_2^2 + k_1^2 k_2^2)-k^4-k_1^4-k_2^4}.\end{equation}

The expression of the matrix elements reads (Lvov et al. Reference Lvov, Tabak, Polzin and Yokoyama2010)

(3.3)\begin{gather} V^{\boldsymbol{p}}_{\boldsymbol{p}_1 \boldsymbol{p}_2} = \sqrt{kk_1k_2} \left(\frac{k^2+k_1^2-k_2^2}{2kk_1}\sqrt{\left|\frac{m_2^\star}{mm_1^\star}\right|} +\frac{k^2+k_2^2-k_1^2}{2kk_2}\sqrt{\left|\frac{m_1^\star}{mm_2^\star}\right|} \right.\nonumber\\ \hspace{-8pc}\left.+\frac{k^2-k_1^2-k_2^2}{2k_1k_2}\sqrt{\left|\frac{m}{m_1^\star m_2^\star}\right|} \right), \end{gather}

(3.4)\begin{gather} {g^{\boldsymbol{p}}_{12}}' = \frac{\textrm{sign} (m_1^\star) k_1}{(m_1^\star)^2} - \frac{\textrm{sign} (m_2^\star) k_2}{(m_2^\star)^2}, \end{gather}

where $m_1^\star , m_2^\star$ are given by the solution of the resonance conditions, i.e. the joint conservation of momentum and energy in each triadic interaction. Thus, in the four-dimensional space spanned by $k_1$, $k_2$, $m_1$, $m_2$, the problem is now restricted to the resonant manifold, parametrized by two independent variables $k_1$ and $k_2$ as summarized in table 1.

Table 1. The six independent solutions to the resonance conditions (Lvov et al. Reference Lvov, Tabak, Polzin and Yokoyama2010).

Note the symmetries of the resonant manifold: the solution $({\rm Ia})$ is obtained from solution $({\rm Ib})$ through permutation of the indices $1\leftrightarrow 2$. We also notice that solutions $({\rm IIa})$, $({\rm IIb})$ reduce to solutions $({\rm IIIa})$, $({\rm IIIb})$, respectively, under permutation of the indices $1\leftrightarrow 2$.

Based on the condition $b=0$ in (2.7), we consider a scale invariant solution which is horizontally isotropic and independent of the vertical wavenumber:

(3.5)\begin{equation} n_{\boldsymbol{p}} = n(\boldsymbol{p}) \propto k^{{-}a}. \end{equation}

Since the collision integral is scale invariant in $k$, it is sufficient to calculate it for a fixed value (e.g. $k=1$ for simplicity), and then retrieve the solution for any value of $k$ by using the scale-invariance relation involving the homogeneity degree of the collision integral.

Integration is performed in the kinematic box, defined by the three triangular relations: $k+k_1\ge k_2$, $k_1+k_2\ge k$, $k+k_2\ge k_1$. We differentiate three different regions of the kinematic box: near-collinear region ($A_C$ and $B_C$), extreme scale-separated region (the infrared region ${\rm IR}$ and the ultraviolet region ${\rm UV}$), and the region of unclassified triads, denoted as ($A_{U}$ and $B_{U}$), as shown in figure 1. Here, $A_C$ and $B_C$ are named near-collinear regions since the resonant triads tend to the collinear limit approaching their boundary given by $k_2=|k_1-k|$, a relationship that can be fulfilled only by degenerate triangles with their sides lying on the same line. The thickness of the regions $A_C$ and $B_C$ is given by the parameter $k_\textrm {IR}$: small values of $k_\textrm {IR}$ imply that the resonant triads inside these regions are close to the collinear limit. We refer to ${\rm IR}$ and ${\rm UV}$ as the extreme scale-separated regions: for the triads in ${\rm IR}$, two wavenumbers have finite horizontal momentum, and one wavenumber has vanishing horizontal momentum. For the triads in ${\rm UV}$, two wavenumbers have very large horizontal momentum, and one wavenumber has a much smaller horizontal momentum. All the possible resonances in ${\rm IR}$ and ${\rm UV}$ constitute the so-called named triads. Finally, $A_{U}$ and $B_{U}$ include all the non-collinear, unclassified triads.Exploiting symmetries, the right-hand side of (2.6), which we denote by $\mathcal {I}(k,m;a,b)$ after introducing the ansatz (2.7), can be reorganized as follows:

(3.6)\begin{align} \mathcal{I}(k,m;a,b) = \left[\int_{A_C} + \int_{A_{U}} + 2 \left(\int_{B_C} + \int_{B_{U}} + \int_\textrm{IR} + \int_\textrm{UV} \right)\right]{\mathcal{J}}(k,k_1,k_2,m) \,\textrm{d} k_1\,\textrm{d} k_2, \end{align}

where a sum over the six solutions to the above resonance conditions is implicit. With $b=0$, the conditions on the exponent $a$ for convergence of the collision integral on the right-hand side of (3.1) come from the infrared (${\rm IR}$, red in figure 1) and the ultraviolet (${\rm UV}$, dark blue in figure 1) regions of integration. The details for the computation of the following results are given in the supplementary materials (Dematteis & Lvov Reference Dematteis and Lvov2020). Both singularities involve a first and a second cancellation between equal and oppositely signed leading terms. For the infrared contribution, we obtain

(3.7)\begin{align} \mathcal{I}_\textrm{IR} \simeq{-}16{\rm \pi} a k^{{-}2a+4} m \int_0^{k_\textrm{IR}/k} {\textrm{d}x} \int_{{-}x}^x {\textrm{d}y}\, x^{{-}a-1} \frac{y^2(y^2-x^2)}{\sqrt{x^2 - y^2}} = 2{\rm \pi}^2 \frac{a}{4-a}m k^{{-}a }k_\textrm{IR}^{{-}a+4}, \end{align}

where $k_\textrm {IR}$ is the (small) height of the red region in figure 1. The integral converges if $a<4$. Also notice that the integral is positive. For the ultraviolet contribution, we obtain

(3.8)\begin{align} \mathcal{I}_\textrm{UV} &\simeq{-}32{\rm \pi} a k^{{-}2a + 4}m\int_0^{k/k_\textrm{UV}} {\textrm{d}x} \int_0^x {\textrm{d}y}\, \frac{k^2}{x^3} x^{a-8} \left[ (x-y)^4 + x^2(x-y)^2\right]/\sqrt{(2x-y)y} \nonumber\\ &\simeq{-}14 {\rm \pi}^2 \frac{a}{a-3}k^{{-}a+1} m k_\textrm{UV}^{3-a}, \end{align}

where $k_\textrm {UV}$ is the $k_1$ coordinate of the left boundary of the ultraviolet region. The integral converges when $a>3$. Note that this contribution is negative, providing possibility for this contribution to balance the positive contribution from (3.7). This observation will later be exploited in § 4.2 to find the steady-state solution composed of a balance of infrared and ultraviolet contributions. The contribution (3.7) is given by the resonance conditions $(\textrm {Ia})$ and $(\textrm {IIa})$, the infrared ID resonances, while the contribution (3.7) is given by the resonance conditions $(\textrm {IIb})$ and $(\textrm {IIIc})$, the ultraviolet ID resonances. Both ES and PSI resonances turn out to be subleading.

Figure 1. The kinematic box is split into subregions: $A_C$ (light blue), $A_{U}$ (yellow), $B_C$ (green), $B_{U}$ (orange), ${\rm IR}$ (red), ${\rm UV}$ (blue). Here, $A_C$ and $B_C$ are the near-collinear regions, ${\rm IR}$ and ${\rm UV}$ are the extreme scale-separated regions, and $A_{U}$ and $B_{U}$ are the unclassified regions. A suitable Zakharov–Kraichnan transformation, see (4.1a,b), maps the regions $B_C$ and $B_{U}$ into $A_C$ and $A_{U}$, respectively.

3.2. Numerical solution: $a=3.69$

Straightforward numerical integration can be performed only in $A_{U}$ and $B_{U}$, since close to the boundaries the integrand contains integrable singularities. For this reason, numerical integration is performed adopting the following technique for integrable singularities (Heath Reference Heath2002). We take the leading-order singularity of the integrand, integrate it analytically and add the numerical integral of the difference between the integrand and the leading-order singularity. This way, the integrable singularities are integrated analytically rather than numerically, ensuring accurate results. Notice that a singular behaviour is found not only in the infrared and ultraviolet regions, but also in the collinear regions, due to the vanishing denominator $\varDelta _{\boldsymbol {p}12}$. The vanishing of denominator occurs because the area of a triangle with collinear sides tends to zero. The detail of the procedure for each of the five regions is found in the supplementary materials (Dematteis & Lvov Reference Dematteis and Lvov2020).

To convince the reader that the numerical integration is performed accurately, we show that the discretized integral is independent of the step size of the discretization grid for sufficiently fine grid. This is demonstrated in the supplementary materials (Dematteis & Lvov Reference Dematteis and Lvov2020).The width of the regions around $k_2=0$ is determined by the parameter $k_\textrm {IR}$, while the cut at large $k$'s is performed at $k_1=k_\textrm {UV}$. For the result to be general, it must be independent of the choice of $k_\textrm {IR}$ and $k_\textrm {UV}$, as long as they are finite numbers, $k_\textrm {IR}$ being sufficiently small and $k_\textrm {UV}$ sufficiently large. It turns out this is indeed the case in our numerics. In the supplementary materials (Dematteis & Lvov Reference Dematteis and Lvov2020), we show how convergence is reached as $k_\textrm {UV}$ increases, as the neglected contribution in ${\rm UV}$ vanishes. Independence of the result upon variations of $k_\textrm {IR}$ is even more robust.

According to our results, the stationary spectrum of internal waves is given by

(3.9)\begin{equation} n(k,m) \propto k^{{-}a_0},\quad a_0\simeq 3.69. \end{equation}

The result in figure 2, obtained with a choice $k_\textrm {UV}=1/k_\textrm {IR}=16$, confirms convergence of the collision integral for $3 < a < 4$: divergence of the integral is found both as $a\to 3^+$ and $a\to 4^-$. In figure 2, we show the contributions of each region of the kinematic box. The $a=3$ divergence is negative and due to the region ${\rm UV}$ (ultraviolet). The contribution of ${\rm UV}$ is always negative and tends to zero as $a\to 4$. On the other hand, the $a=4$ divergence is positive and due to the region ${\rm IR}$. The contribution of ${\rm IR}$ is always positive and tends to zero as $a\to 3$. At the stationary solution $a=a_0$, the contributions of ${\rm IR}$, $A_{U}$, $B_{U}$, ${\rm UV}$ are close to zero, while the contributions of $A_{C}$ (positive) and $B_{C}$ (negative) are large and cancel out. Notice that the results are obtained for $A_C$ and $B_C$ being thin slices of width $1/16$: the points in the collinear region correspond to triads of wave vectors that have angles between each other's horizontal components of $3^\circ$ or less!

Figure 2. (a) Contributions of each subregion (as split in (3.6)) to the integral for $b=0$ and varying $a$. (b) The base-$10$ logarithm of the magnitude of the integrand is shown, for the solution $a=3.69$, $b=0$. The colourmap labelled by the left colourbar indicates negative values, and the right colourbar indicates positive values. Also, we show here the schematic representation of the downscale energy transfers. The thicker arrow represents the stationary transfer between near-collinear regions, and the thinner arrow between regions with extreme scale separation. The fluxes of energy are explained in § 5.

Let us consider an arbitrary reference number $k=1$. The sign of the integrand in figure 2(b) indicates the direction of the energy transfers. We see from figure 2 that the contribution from the wavenumbers with $k_1<1$ is positive. Therefore, there is a net energy flow from wavenumbers smaller than reference number $k=1$ to the reference wavenumber $k=1$. On the other hand, the contribution from the wavenumbers with $k_1>1$ is negative. Consequently the wavenumber $k=1$ constantly pumps energy towards higher wavenumbers. Therefore, we conclude that the energy transfer is directed towards high horizontal wavenumbers. We elaborate on this further in § 5.

The outflowing energy from the small-wavenumber near-collinear triads is balanced by the inflowing energy at the large-wavenumber near-collinear triads, implying a stationary flow mediated by $k=1$, which is represented as a thick arrow in figure 2. The outflowing energy from the infrared region is balanced by the inflowing energy entering the ultraviolet region, giving a stationary energy flow between these two regions mediated by $k=1$. This energy transfer is represented as a thin directed arrow connecting the two regions. The quantitative justification of the balance is given by figure 2(a).

The value of the exponent appears to be characterized importantly by a balance of the regions $A_C$ and $B_C$. This suggests that a suitable transformation mapping one region into the other could potentially make the search for the steady self-similar spectrum amenable to analytical treatment. This task is addressed in the next section.

4.1. Collinear limit

We start by analysing the contributions in the regions $A_C$ and $B_C$ by decomposing them into separate subcontributions from the three resonance types $(\textrm {I})$, $(\textrm {II})$ and $(\textrm {III})$ (see table 1). The different contributions are plotted in figure 3. We realize that in the region $A_C$ the leading contribution is given by the resonance condition $(\textrm {I})$, while in region $B_C$ the leading contribution is given by the resonance condition $(\textrm {II})$. Let us consider only these two contributions, naming them the main collinear contributions:

(4.2)\begin{align} \mathcal{I}_{A_C}+2\mathcal{I}_{B_C} &\simeq \int_{k_\textrm{IR}}^{1-k_\textrm{IR}} \textrm{d} k_2 \int_0^{k_\textrm{IR}} {\textrm{d}x} \frac{T^k_{k_1, k-k_1}}{\sqrt{2kk_1(k-k_1)x}} \nonumber\\ &\quad - \int_{1+k_\textrm{IR}}^{k_\textrm{UV}} \textrm{d} k_2 \int_0^{k_\textrm{IR}} {\textrm{d}x}\, \frac{T^{k_1}_{k, k_1-k}}{\sqrt{2k k_1(k_1-k)x}} \nonumber\\ &\quad - \int_{1+k_\textrm{IR}}^{k_\textrm{UV}} \textrm{d} k_2 \int_0^{k_\textrm{IR}} {\textrm{d}x} \,\frac{T^{k_2}_{k, k_2-k}}{\sqrt{2k k_2(k_2-k)x}}, \end{align}

where $T^0_{12} = kk_1k_2 |V^0_{12}|^2 f^0_{12}/|{g^0_{12}}'|$. Transforming the integral into an integral in the region $A_C$, by using the symmetries of (4.1a,b) (with $k_\textrm {UV}=k^2/k_\textrm {IR}$) and the scale-invariant properties of the integrand, we obtain (renaming $\tilde k_1 \rightarrow k_1$, $\tilde k_2 \rightarrow k_2$)

(4.3)\begin{align} \mathcal{I}_{A_C}+2\mathcal{I}_{B_C} \simeq \int_0^{k_\textrm{IR}} {\textrm{d}x} \int_{k_\textrm{IR}}^{1-k_\textrm{IR}} \textrm{d} k_2 \,\frac{T^k_{k_1, k-k_1}}{\sqrt{2kk_1(k-k_1)x}} \left[ 1 - \left( \frac{k}{k_1} \right)^{r+3} - \left( \frac{k}{k_2} \right)^{r+3} \right], \end{align}

where $r$ is the degree of homogeneity of the integrand: $r=3-2a$ (same dependence on $k$ as in the computation of the PR spectrum). We used the property of the Zakharov–Kraichnan transformation for an homogeneous function of degree $r$ with respect to horizontal wave numbers:

(4.4)\begin{equation} \mathcal{J}\left(k_1,k,k_2,m\right)=\mathcal{J}\left(\frac{k}{\tilde k_1}k, \frac{k}{\tilde k_1} \tilde k_1,\frac{k}{\tilde k_1} \tilde k_2,m\right) = \left(\frac{k}{\tilde k_1}\right)^r \mathcal{J}(\tilde k,\tilde k_1,\tilde k_2,m), \end{equation}

and a factor $(k/\tilde {k_1})^3$ appeared due to the Jacobian of the coordinate change. Now, we notice that if $r+3=-1$, the integrand contains a factor $[k-k_1-k_2]$. Since we are considering the near-collinear region, horizontal momentum conservation implies that such a factor vanishes (more precisely, it would vanish in the limit $k_\textrm {IR}\to 0$): $A_C$ is a thin slice lying on the line $k_2 = k-k_1$. Therefore, we have that $\mathcal {I}_{A_C}+2\mathcal {I}_{B_C}=0$ for $6-2a=-1$, which determines the value of $a=7/2$ as the critical exponent. We have therefore found analytically the steady-state solution for the reduced kinetic equation dominated by the balance between the contribution of type $(\textrm {I})$ in the region $A_C$ and the contribution of type $(\textrm {II})$ in the region $B_C$, which are the largest contributions from the near-collinear regions. This solution is therefore given by

(4.5)\begin{equation} n(k,m) =k^{-({7}/{2})}m^0. \end{equation}

Since this solution is close to the PR spectrum (2.8), we propose to call this solution the modified PR spectrum. The difference between (4.5) and (2.8) is that the latter is the formal solution, corresponding to a non-local spectrum (i.e. implying a divergent collision integral). The former solution, on the other hand, is a physically relevant solution corresponding to a local action spectrum (i.e. whose collision integral is finite). Note, however, that a part of the resonances have been neglected.

Figure 3. The contributions from $A_C$ (a), and $B_C$ (b) (the latter multiplied by $2$ to account for its symmetric $B_C'$ by permutation $k_1\leftrightarrow k_2$) are split into their subcontributions from the three resonance types, showing that $A_C$ is dominated by the contribution $(\textrm {I})$, while $B_C$ is dominated by the contribution $(\textrm {II})$. Computing the balance between these two contributions leads to a theoretical estimate $a = 7/2$.

The result $a=7/2$ coming from the collinear limit of (4.3) involves only the two leading contributions in figure 3. This observation provides an intuition on how the exponent $a$ is determined by the kinetic equation. The sum of the subdominant contributions in figure 3 is negative and almost independent of $a$. When added to the main contribution crossing zero at $a=7/2$, this negative contribution makes the zero-crossing point shift toward the right and in figure 4(a) the total collinear contribution is shown to cross zero at $a\simeq 3.69$.

Figure 4. (a) The net contribution of the three types of triads, showing that each type balances to zero independently at the convergent stationary solution. (b) Relative error of (3.8) with respect to the fully numerically integrated contribution of region $\textrm {UV}$. Around $a=3.7$, we see that a choice of $k_\textrm {UV}=16$ ($y=4$) implies an error of approximately $2\,\%$, which we consider acceptable. Thus, we consider $k_\textrm {UV}=16$ as a reasonable delimitation of the ultraviolet region, as well as $k_\textrm {IR}=1/16$ as the delimitation of the infrared region.

4.2. Extreme scale-separated triads

In this section, we consider the contribution that comes from the extreme scale-separated triads, the sum of the infrared and the ultraviolet contributions. The former is positive and tends to $+\infty$ as $a\to 4^-$; the latter is negative and tends to $-\infty$ as $a\to 3^+$. In figure 4(a), we show the total contribution of the extreme scale-separated resonances and we observe numerically that it crosses zero around $a\simeq 3.69$, too.In § 3.1, we proposed a balance between positive infrared (3.7) and negative ultraviolet (3.8) contributions as a way to form the steady-state spectrum of internal waves. This balance hinges upon the choice $k_\textrm {UV}=k^2/k_\textrm {IR}$, as explained in § 4.1.

In figure 4, we show the relative error of the leading-order analytical expression (3.8), with respect to the numerically computed contribution, as a function of $k_\textrm {UV}$. The result is shown for $a=3.5$ and $a=3.7$. Using the expression $k_\textrm {UV}=2^y$, with $k=1$, we observe that the analytical approximation is good starting from values of $y$ between $4$ and $5$. We choose $y$ (and therefore $k_\textrm {UV}$ and $k_\textrm {IR}$) large enough for an accurate approximation with the leading-order expression (3.8) (or (3.7)), yet small enough so that all extreme scale-separated triads (3.8) and (3.7) are actually included in the ultraviolet and infrared regions. We make an arbitrary choice $y=4$ ($k_\textrm {UV}=1/k_\textrm {IR}=16$) so that the leading-order error of (3.7) and (3.8) is approximately $2\,\%$ (see figure 4). Our results are insensitive to this specific choice. The balance between the two expressions gives

(4.6)\begin{equation} \left.\begin{gathered} 2{\rm \pi}^2 \frac{a}{4-a} k_\textrm{IR}^{4-a} = 14{\rm \pi}^2\frac{a}{a-3}\frac{1}{k_\textrm{IR}^{3-a}}\\ h_y(a):=(2a-7) y \log 2 = \log 7 +\log\left(\frac{4-a}{a-3}\right)=:g(a). \end{gathered}\right\} \end{equation}

First, we notice that the presence of a factor $7$ on the right-hand side breaks the symmetry that would imply the two contributions balance out at $a=7/2$, in the middle of the convergence interval $(3,4)$. Secondly, we notice that the function $g(a)$ has an inflection point at $a=7/2$, making its Taylor expansion of first order have an error of third order. Using linear interpolation centred at $a=7/2$, $g_{7/2}(a)=\log 7 -2(2a-7)$, as an accurate approximation to $g(a)$, demanding that $h_y(a) = g_{7/2}(a)$, we obtain

(4.7)\begin{equation} a= \frac72 +\frac{\log 7}{2y\log2 + 4}. \end{equation}

If we adopt $y=4$, as chosen throughout the paper, we obtain the solution $a\simeq 3.70$. Using $y=5$ would yield a solution $a\simeq 3.68$.

In this section, we have shown how the formation of the stationary solution of the wave kinetic equation can be interpreted via two independent balances, between near-collinear triads and between the ID triads of the extreme scale-separated regions. This is consistent with recent results from direct numerical simulations where the dominating interactions were located in the ID regime but also throughout all of the interval $k_1\in [0,1.4k]$ (Pan et al. Reference Pan, Arbic, Nelson, Menemenlis, Peltier, Xu and Li2020). Despite the latter effect was investigated in Pan et al. (Reference Pan, Arbic, Nelson, Menemenlis, Peltier, Xu and Li2020) as a broadening of ID due to large nonlinearity, we find that it may be consistent with the collinear resonances depicted in figure 2(b).

5.1. Physical dimensions and energy conservation

Using the scale-invariant properties of the collision integral, the right-hand side of (3.1) can be rewritten considering the appropriate physical dimensions as

(5.1)\begin{equation} I(\boldsymbol{k},m) = |m|^{{-}2b+1} k^{{-}2a+4} (V_0A)^2 \mathcal{I}(k=1,m=1;a,b), \end{equation}

where $\mathcal {I}(k=1,m=1;a,b)$ is non-dimensional, $V_0$ is the dimensional prefactor of the matrix element defined below in (5.5) and $A$ is the prefactor of the GM spectrum defined in (2.2a,b) above. A simple way to check the dimensional consistency of the prefactor in (5.1) is to consider the contribution from the extreme scale-separated region, (4.6), which is analytically tractable:

(5.2)\begin{align} & 2{\rm \pi}^2 \frac{a}{4-a} k^{{-}a} k_\textrm{IR}^{4-a} - 14{\rm \pi}^2 \frac{a}{a-3}\frac{k^{{-}a+1}}{k_\textrm{IR}^{3-a}} \nonumber\\ &\quad = mk^{{-}2a+4} \left(2{\rm \pi}^2 \frac{a}{4-a} x_\textrm{IR}^{4-a} - 14{\rm \pi}^2 \frac{a}{a-3}x_\textrm{IR}^{a-3}\right), \end{align}

where $x_\textrm {IR}=k_\textrm {IR}/k$ and the term in brackets is a non-dimensional function of the exponent $a$ that vanishes at $a=a_0$. The factor $(V_0 A)^2$ comes from having both the matrix elements and the spectrum to the second power in the collision integral, and by introducing the appropriate dimensional constants that have been omitted so far. The dimensional properties of the collision integral are indeed the same also in the other integration regions.

Next, we compute the spectral energy fluxes, recalling that the energy density is given by

(5.3a,b)\begin{equation} e(k,m) = 2{\rm \pi} k \sigma(\boldsymbol{k},m) n(\boldsymbol{k},m),\quad \sigma(\boldsymbol{k},m)=\alpha \frac{k}{|m|},\quad \text{with}\ \alpha=\frac{g}{\rho_0 N}, \end{equation}

where the scale-invariant dispersion relation (2.4) is used. Using (5.1), the stationary-wave kinetic equation for the energy density assumes the simple form

(5.4)\begin{equation} \dot e(k,m) = 2{\rm \pi} (A V_0)^2 k^{{-}2a+6} m^{{-}2b} \mathcal{I}(k=1,m=1;a,b)=0. \end{equation}

5.2. Dimensional prefactors

The dimensional factor coming from the matrix elements in (5.1) is given by

(5.5)\begin{equation} |V_0|^2= \frac{N}{32\rho_0}. \end{equation}

The factor $A$ is the dimensional prefactor of the GM spectrum, our observational input for the oceanic wave field. The procedure to obtain $A$ is explained in the rest of this paragraph. The non-rotating limit of the GM spectrum in $k-m$ coordinates is given by (2.3). However, in isopycnal coordinates the spectrum needs to be multiplied by a factor $N^2\rho _0/g$, which gives

(5.6)\begin{equation} n_\textrm{GM}(k,m) = \frac{1}{{\rm \pi}^3 g} E b^2 N_0 N\rho_0 f m_\star k^{{-}4},\quad m_\star{=} 3{\rm \pi}\frac{N}{b N_0}, \end{equation}

where $E=6.3\times 10^{-5}$ is the non-dimensional energy level, $b=1300$ m, $\rho _0=1000$ kg m$^{-3}$, $N_0=0.00524$ s$^{-1}$ and $f = 2\cdot 7.3\times 10^{-5} \sin (l)$ (at latitude $l=32.5^\circ$). (A prefactor of $4$ instead of $3$ is known to imply a more accurate asymptotic fit in the large-wavenumber regime (Polzin & Lvov Reference Polzin and Lvov2011). However, for simplicity here we keep the factor appearing in the original 1976 GM parametrization as is. We will see that the choice does not affect the order of magnitude of the estimate of the flux.) The values given here are the ones of the standard GM parametrization. Since for the GM spectrum we have $a=4$ instead of $a=3.69$, we introduce the modified version of (5.6):

(5.7)\begin{equation} n_{\textrm{GM},c}(k,m) = \frac{1}{{\rm \pi}^3 g} E b^2 N_0 N \rho_0 f \frac{m_\star}{{k_\star}^{(1-s)/2}} k^{-(4-(1-s)/2)},\quad k_\star{=} m_\star{/}r, \end{equation}

where $s=2a-7 = 0.38$. This is a slightly less steep, dimensionally consistent version. The non-dimensional parameter $r=N/f$ quantifies the horizontal-to-vertical anisotropy (Polzin & Lvov Reference Polzin and Lvov2011). Now, the GM spectrum is normalized so that the total energy is expressed in units of J$/$kg, as usual in physical oceanography. On the other hand, in the wave-turbulence formalism, the total energy is expressed as a density per unit of volume of the physical space. Here, the physical space has units of $\textrm {kg}\ \textrm {m}^{-1}$, given by an area in the horizontal directions times a density in the vertical. It is, therefore, possible to switch from one representation to the other multiplying by the appropriate density, which in this case is the characteristic density of the isopycnal coordinates, $\varPi=g/N^2$, in units of metres (normalized differential thickness of an isopycnal layer). So, we obtain the equivalence: $n_{WKE}=n_\textrm {GM} \cdot g/N^2$, which applied to (5.7) finally gives the dimensionally consistent factor

(5.8)\begin{equation} A=\frac{1}{{\rm \pi}^3} E b^2 \rho_0\, f \frac{N_0}{N}\frac{m_\star}{k_\star^{(1-s)/2}}. \end{equation}

5.3. Dissipated power at high wavenumbers

In order to compute the energy flux towards high wavenumbers, we need to consider the physical cutoffs of the problem. Natural cutoffs are imposed on the vertical wavenumber by the depth of the ocean and by the wave breaking cutoff, and on the frequency by the inertial frequency and the buoyancy frequency,

(5.9a–d)\begin{equation} m_{min} =\frac{2{\rm \pi} g}{\rho_0 N^2}\times(2600\ \text{m})^{{-}1},\quad m_{max} =\frac{2{\rm \pi} g}{\rho_0 N^2}\times(10\ \text{m})^{{-}1},\quad \sigma_{min} = f,\quad \sigma_{max} = N. \end{equation}

These limiting values define a rectangle in $\sigma - m$ space that translates into a trapezoid in $k-m$ space, with inclined sides given by

(5.10a,b)\begin{equation} k_{min}(m) =\frac{f}{\alpha}m,\quad k_{max}(m) =\frac{N}{\alpha}m. \end{equation}

The collision integral contains all of the necessary information on the spectral energy transfers. In the following, we compute numerically the outflowing power at high wavenumbers from computation of the contribution of the resonant triads with an output wavenumber such that $m>m_{max}$ or $k>k_{max}(m)$, i.e. assuming that the production of a wave beyond the physical high wavenumber cutoff results in complete dissipation of its energy. At the same time, it is assumed that the region within the physical cutoffs is an inertial range with no sources nor sinks, where energy is transferred exclusively via resonant interactions.

Let us define the part of the collision integral which contributes to the dissipation of energy by transferring it beyond the dissipation cutoffs:

(5.11)\begin{equation} \left.\begin{gathered} I_{diss}(\boldsymbol{k},m;k_{max}): = \frac{{N}^2}{g} (V_0 A)^2 \alpha^{{-}1} |m|^{{-}2b+1} k^{{-}2a+4} \mathcal{I}_{diss}(\kappa),\\ \mathcal{I}_{diss}(\kappa) = \int_{\varOmega_{h}(\kappa)} \mathcal{J}(k=1,k_1,k_2,m=1) \textrm{d} k_1 \,\textrm{d} k_2, \end{gathered}\right\} \end{equation}

where $\kappa = k_{max}/k$ and $\varOmega _{h}(\kappa )$ is the set of triads transferring energy to output waves beyond the horizontal dissipation cutoff $k_{max}$, and is defined in table 2. Moreover, a factor ${{N}^2}/{g}$ is added to account for transition to isopycnal coordinates, the inverse of the factor in (5.8). Here, $I_{diss}(\boldsymbol {k},m;k_{max})$ quantifies the amount of wave action that wavenumber $\boldsymbol {p}$ sends via resonant interactions beyond the dissipation threshold $k_{max}$, per unit time and per unit volume of Fourier space.

Table 2. All of the non-negligible contributions for $\kappa >1$ are negative (see (b) of figures 2 and 3). The same holds for $\mu >1$. Thus, wavenumber $\boldsymbol {p}$ is always an incoming wave in the triads here considered. The table represents the conditions under which every type of resonance results in the dissipation of energy due to at least one output beyond the dissipation cutoffs. During the numerical computation of the integrals $C_{h}$ and $C_{v}$, these conditions are applied for every point of the kinematic box. In the case of a decay into two waves of which only one is dissipated, only the fraction of energy of the dissipated wave must be accounted for, multiplying the contribution by the weight shown in the table.

The power (per unit of m) dissipated beyond $k_{max}$ at fixed $m$ (now considered as the positive definite magnitude of m) is related to the integral of $I_{diss}(\boldsymbol {k},m)$ for all values of $k$:

(5.12)\begin{align} F_{diss}(m) &= \int_{k_{min}}^{k_{max}} 4{\rm \pi} k \omega(\boldsymbol{p}) I_{diss}(\boldsymbol{k},m;k_{max}) \,\textrm{d} k \nonumber\\ &= \frac{{N}^2}{g}\int_{k_{min}}^{k_{max}} 4{\rm \pi} |m|^{{-}2b+1}k^{{-}2a+5} \alpha \frac{k}{|m|} (V_0 A)^2 \alpha^{{-}1} \mathcal{I}_{diss}(\kappa)\textrm{d} k \nonumber\\ &= 4{\rm \pi} \frac{{N}^2}{g} (V_0 A)^2 k_{max}^{{-}s} C_{h},\quad \text{with}\ C_{h} := \int_1^{{N}/{f}} \kappa^{s-1} \mathcal{I}_{diss}(\kappa)\textrm{d}\kappa, \end{align}

where we recall that $s=2a-7=0.38$, and $b=0$. In a similar fashion, we can obtain an analogous relation for the energy flux dissipated vertically at vertical wavenumbers larger than $m_{max}$:

(5.13)\begin{equation} \left.\begin{gathered} G_{diss}(k) = 4{\rm \pi} \frac{{N}^2}{g} (V_0 A)^2 k^{-(1+s)}{m_{max}} C_{v},\quad C_{v} := \int_1^{{m_{max}}/{m_{min}}} \mu^{{-}2} \mathcal{K}_{diss}(\mu)\textrm{d}\mu,\\ \mathcal{K}_{diss}(\mu) = \int_{\varOmega_{v}(\kappa)} \mathcal{J}(k=1,k_1,k_2,m=1) \,\textrm{d} k_1 \,\textrm{d} k_2, \end{gathered}\right\} \end{equation}

where $\mathcal {K}_{diss}(\mu )$ ($\mu=m_{max}/m$) is the analogue of $\mathcal {I}_{diss}(\kappa )$ in the vertical direction, i.e. the collision integral restricted to the triads with output waves beyond the $m_{max}$ threshold. The conditions defining $\varOmega _{v}$ are found in table 2. These integrals are computed numerically. The computation of $C_{h}$ is quite straightforward since the kinematic box is expressed in horizontal wavenumber coordinates. The quantity $\mathcal {I}_{diss}(\kappa )$ is computed inside a loop spanning all values of $\kappa$, checking the constraint of $\varOmega _{h}(\kappa )$ for every point. Note that since the position of the right boundary depends on $m$, (5.10a,b), in the computation of $\mathcal {I}_{diss}(\kappa )$ we have to impose that $k_1>k_{max}m_1/m$ (or the same for $k_2$) for point $(k_1,m_1)$ to be past the absorbing boundary. For the computation of $C_{v}$, on the other hand, at every loop iteration we have to integrate over all of the kinematic box to compute $\mathcal {K}_{diss}$ and check point by point whether $m_1$ or $m_2$ have ‘crossed’ the boundary at $m_{max}$ (which is independent of $k$).

One important point requires particular care. In the previous sections, we relied upon the leading order of the infrared contribution of (3.7), which is obtained after the second cancellation where a negative singularity for $k_1>1$ cancels exactly with a positive singularity for $k_1<1$. However, since now in (5.11) we consider an integration region where $k_1>\kappa >1$, the positive singularity due to $k_1<1$ is no longer present and thus the use of (3.7) is not justified. In the integrand of (5.11), whose expression is found in the supplementary materials (Dematteis & Lvov Reference Dematteis and Lvov2020), the finite point singularity therefore has exponent $-a+3/2$ rather than the $-a+2$ resulting after the second cancellation (leading to (3.7)). Integrating twice according to (5.11), the exponent of the singularity of the integrand defining $C_{h}$ in (5.12) is $\beta _{h} = -a+7/2\simeq 0.19$. Numerical computation of $\mathcal {I}_{diss}(\kappa )$ as $\kappa \to 1^+$ confirms this analytical prediction and gives a best least square fit of $\mathcal {I}_{diss}(\kappa ) \sim -d_{h}(\kappa -1)^{-\beta _{h}}$, for $\kappa -1\ll 1$, with $d_{h}=127\times 8{\rm \pi}$, $\beta _{h}=0.19$. An analogous computation for the integral $C_{v}$ defined in (5.13) leads to a best-fit scaling of its integrand of $\mathcal {K}_{diss}(\mu ) \sim -d_{v}(\mu -1)^{-\beta _{v}}$, for $\mu -1\ll 1$, with $d_{v}\simeq 9.0\times 8{\rm \pi}$, $\beta _{v}\simeq 0.75$. In this case, an analytical evidence of the exponent is not simply available since the kinematic box is not expressed in vertical wavenumber coordinates. In figure 5, we show the rapidly decaying behaviour of the integrands of $C_{h}$ and $C_{v}$. Therefore, $C_{h}$ and $C_{v}$ are ‘universal constants’ which are relatively insensitive to the limits of integration in (5.12) and (5.13). The figure also shows a closer look to the scaling of the singularities as $\kappa \to 1^+$ and $\mu \to 1^+$, with the respective best-fit scalings. Therefore we obtain:

(5.14)\begin{align} C_{h} \simeq \int_1^{{N}/{f}} \kappa^{s-1} \mathcal{I}_{diss}(\kappa) &\simeq{-}\frac{d_{h} {k_\textrm{IR}}^{1-\beta_{h}}}{1-\beta_{h}} + \int_{1+k_\textrm{IR}}^{{N}/{f}} \textrm{d}\kappa \kappa^{s-1} \mathcal{I}_{diss}(\kappa) \nonumber\\ &\simeq{-}8{\rm \pi} (16.6 + 54.6 + 10.2) \simeq -8{\rm \pi} \times 81.4.\end{align}

Here, the three numbers $16.6$, $54.6$ and $10.2$ represent the partial contributions by the infrared ID interactions, near-collinear interactions and unclassified interactions, respectively. The ultraviolet contribution is negligible.Then, choosing a value $\mu _\textrm {IR}=1/16$ to delimit the infrared region for vertical wavenumbers, we have

(5.15)\begin{align} C_{v} \simeq \int_1^{{m_{max}}/{m_{min}}} \mu^{{-}2} \mathcal{K}_{diss}(\mu) &\simeq{-}\frac{\textrm{d}_{v} {\mu_\textrm{IR}}^{1-\beta_{v}}}{1-\beta_{v}} + \int_{1+\mu_\textrm{IR}}^{{m_{max}}/{m_{min}}} \textrm{d}\mu \,\mu^{{-}2}\mathcal{K}_{diss}(\mu) \nonumber\\ &\simeq{-} 8{\rm \pi}(18.0 + 14.0 + 3.1) \simeq{-}8{\rm \pi}\times 35.1. \end{align}

Again, the numbers $18.0$, $14.0$ and $3.1$ represents infrared (ID), near-collinear and unclassified contributions, respectively, with a negligible contribution from the ultraviolet regions.

Figure 5. Panels (a) and (b) show the integrands of $C_{h}$ and $C_{v}$, respectively, and the areas coloured in red represent $C_{h}$ and $C_{v}$ themselves. Both integrals have a finite point integrable singularity (the exponents $\beta _{h}$ and $\beta _{v}$ are in the interval $(0,1)$) at their left boundaries: the insets show the asymptotic power-law scaling of the singularities.

Now, we compute the power crossing the $k_{max}$ boundary, denoted as $P_{h}$ (horizontal boundary) and the power crossing the $m_{max}$ boundary, denoted as $P_{v}$ (vertical boundary), via

(5.16)\begin{equation} \left.\begin{gathered} P_{h} = \int_{m_{min}}^{m_{max}} F_{diss}(m)\,\textrm{d} m = D_{h}(1-s)^{{-}1} \left(\frac{N}{\alpha}\right)^{{-}s} (m_{max} - m_{min})^{1-s},\\ P_{v} = \int_{({f}/{\alpha})m_{max}}^{({N}/{\alpha})m_{max}} G_{diss}(k)\,\textrm{d} k = D_{v} s^{{-}1}\left[\left(\frac{f}{\alpha}\right)^{{-}s}- \left(\frac{N}{\alpha}\right)^{{-}s}\right]m_{max}^{1-s}, \end{gathered}\right\} \end{equation}

where

(5.17a,b)\begin{equation} D_{h} = 4{\rm \pi} \frac{{N}^2}{g}(V_0 A)^2 C_{h},\quad D_{v} = 4{\rm \pi} \frac{{N}^2}{g}(V_0 A)^2 C_{v} \end{equation}

Using (5.5), (5.8), (5.16) and (5.17a,b), we obtain

(5.18) \begin{equation} \left.\begin{gathered} P_{h} = \frac{\varGamma C_{h}}{1-s}\left[1- \left(\frac{\ell}{2b}\right)^s \right] f^{1+s} N^{2}E^2,\\ P_{v} = \frac{\varGamma C_{v}}{s}\left[1-\left(\frac{f}{N}\right)^s\right] f N^{1+s} E^2, \\ \varGamma = \frac{3}{4{\rm \pi}^3} \left(\frac{3\ell}{2b}\right)^s \frac{\rho_0 b^2 N_0^{1-s}}{g \ell},\quad s=0.38, \quad \ell = 10\ {\rm m}. \end{gathered}\right\} \end{equation}

In figure 6, a graphical interpretation of the computations in the above paragraph is shown. The physical cutoffs define a rectangular box (figure 6a) within which energy is transferred through resonant interactions toward high vertical wavenumbers and high horizontal wavenumbers. If energy is provided by large-scale forcing at low frequency, for simplicity represented as a mesoscale eddy in the figure, it ends up being dissipated at high frequencies or high wavenumbers, represented as wave breaking with generation of turbulent vortices. In $k-m$ , space this inertial box translates into a trapezoid whose lateral sides are straight lines defined by the dispersion relation. The integrals in (5.16) quantify the contribution to the dissipated power along the dissipative sides of the box. In the figure, the red areas quantify how dissipation is distributed along these dissipative sides of the inertial box.

Figure 6. (a) The inertial box in $\sigma -m$ space. (b) The inertial box in $k-m$ space. The low frequency energy forcing is represented in the picture as a mesoscale eddy, while wave breaking (with generation of turbulent vortices) represents dissipation at high wavenumbers. The distribution of energy dissipation along the absorbing boundaries is depicted in the two insets showing $F_{diss}(m)$ and $G_{diss}(k)$, whose integrals (red area) give, respectively, $P_{h}$ and $P_{v}$, (5.16). In these two plots, $k$ is in units of $\textrm {m}^{-1}$, $m$ is in units of $\textrm {m}^3\ \textrm {kg}^{-1}$, $G_{diss}(k)$ is in $\rm {W}\ \textrm {m}\ \textrm {kg}^{-1}$ and $F_{diss}(m)$ is in $\rm {W}\ \textrm {m}^{-3}$, so that both integrals give a power per unit mass.

Equation (5.18) is the main result of our paper. The outgoing power toward large wavenumbers is the quantity that is modelled by the finescale parametrization formula of ocean mixing (Polzin et al. Reference Polzin, Toole and Schmitt1995) as

(5.19)\begin{equation} P_{fs} = (1.9 \times 10^7\ \textrm{m}^2)\cosh^{{-}1}\left(\frac{N}{f}\right) f N^2 E^2, \end{equation}

which is in agreement with (5.18) regarding the dependence upon the main physical parameters, if the lower-order corrections $(f/N)^{s}$ and $\cosh ^{-1}(N/f)$ are neglected. The predicted intensities can also be compared directly, using the standard parameters of the GM spectrum (5.6) with $N=N_0$, which yields

(5.20)\begin{equation} \left.\begin{gathered} P_{h} \simeq{-}0.5 \times 10^{-8}\ \textrm{W kg}^{{-}1},\\ P_{v} \simeq{-}1.5 \times 10^{-8}\ \textrm{W kg}^{{-}1}. \end{gathered}\right\} \end{equation}

This amounts to a total dissipated power

(5.21)\begin{equation} P_{tot} = P_{h} + P_{v} \simeq{-}2.0 \times 10^{-8}\ \textrm{W kg}^{{-}1}, \end{equation}

which is negative since it is lost by the wave system considered. With the same parameters and the same sign convention, the finescale parametrization formula predicts a total dissipated power

(5.22)\begin{equation} P_{fs} ={-} 7.9 \times 10^{{-}10}\ \textrm{W kg}^{{-}1}. \end{equation}

These two predictions are in qualitative agreement, but with a difference of about an order of magnitude. (Had we used a factor of $4$ in place of $3$ in the definition of $m_\star$, (5.6), the fluxes would have to be multiplied by a factor of $(4/3)^{1+s}\simeq 1.5$, which does not affect the order of magnitude of the estimate. Here, we use the original choice of $m_\star$ from Cairns & Williams (Reference Cairns and Williams1976) for simplicity. We hope to consider oceanographic implications of this choice in details in future publications.) We elucidate on possible reasons in the conclusions. Moreover, the following remarks are important for the sake of clarity.

Above, we have used the terminology dissipated power to indicate the amount of energy that is absorbed by the boundary sink at small scales. As far as the internal-wave kinetic equation is concerned, this sink is acting as an actual dissipation term. In the finescale parametrization picture to which we refer (see Polzin et al. Reference Polzin, Garabato, Huussen, Sloyan and Waterman2014, § 2), this amount of energy is converted into turbulent kinetic energy at small scales, and therefore roughly equals the turbulent energy production. In turn, the latter contributes separately both to diapycnal mixing and to heat production, in proportions that can be quantified by closure assumptions beyond the scope of the present paper. We stress that the use of the word dissipation in this paper does not refer to dissipation of energy as heat, but to the role played by the sink that takes energy out of the wave system considered.

The energy fluxes at the boundaries, (5.12) and (5.13), show a non-trivial dependence on the variables $k$ and $m$. In non-isotropic systems, rather than one single stationary state, it is common to have a family of stationary states (Nazarenko Reference Nazarenko2011, § 9.2.5) whose energy fluxes depend on $k$ and $m$ in different ways. The generalized KZ spectrum (that is the PR spectrum, for internal waves) is only one among the members of the family. Then, the locality conditions imply that for the internal waves only one solution is physically meaningful, i.e. $(a,b)=(3.69,0)$. Indeed, since the solution is stationary, for any enclosed region in $k-m$ space (contained in the inertial range where no sources or sinks are present), the incoming energy flux must be equal to the outgoing energy flux (and opposite in sign). Equivalently, as long as both ends of a boundary are fixed, the total outgoing flux through the boundary must be independent of the path of the boundary in $k-m$ space. For instance, one could show that the flux through the dissipative boundary equals the flux through the forcing boundary (respectively, the red and the black boundaries in figure 6).

Finally, we acknowledge the fact that the solutions with $b=0$ are known in the literature as no flux solutions to the Fokker–Planck (diffusive) approximation to the kinetic equation (3.1). How this approximation is achieved is shown for instance in Lvov et al. (Reference Lvov, Tabak, Polzin and Yokoyama2010) by use of a straightforward leading-order approximation of the infrared ID interactions. An exact leading-order balance makes the infrared flux apparently vanish, whose trace in this paper is found in the two exact cancellations in the infrared region. However, our computations leading to (3.7) show how the next leading-order terms carry a small but non-zero flux toward high horizontal wavenumbers. Moreover, a net energy flux is indeed due to interactions outside the infrared region, such as the collinear interactions. Therefore, there is no mathematical contradiction in having $b=0$ and a non-zero downscale energy flux. Indeed, the considerations in the present paper go way beyond the diffusive approximation to the internal-wave kinetic equation, where the concept of no flux solutions was concocted.