Turbulent spectrum of 2D internal gravity waves

Abstract

We find the turbulent energy spectrum of weakly interacting 2D internal gravity waves without relying on the hydrostatic approximation. This spectrum is an exact solution of the kinetic equation after the shear modes are removed, and it agrees with the oceanic Garrett–Munk spectrum for frequencies large compared to the Coriolis frequency and vertical scales small compared to the depth of the ocean. Among a continuous family of solutions to the kinetic equation, we show that the turbulent solution is the special solution with non zero angular dependent radial flux. Our solution provides an interesting insight of how turbulent energy cascade proceeds in anisotropic systems - similarly to isotropic turbulence it is self-similar in scale, but its angular part is peaked along the curve of vanishing frequency.

Introduction In 1966, Hasselmann laid the groundwork by introducing perturbative equations for weakly interacting waves [1]. Since then, despite a wealth of oceanic measurements and continuous theoretical advancements [2, 3, 4, 5], including a derivation of the renowned phenomenological Garrett-Munk spectrum [6, 7], a theoretical derivation of the energy spectrum for weakly interacting internal gravity waves remains an open problem. A promising avenue lies in the kinetic approach; however, the Boussinesq equation both in 2D and 3D is an anisotropic, non-canonical Hamiltonian equation, which turns its kinetic description into an interesting and non-trivial problem since the classical wave turbulence approach is almost irrelevant. Indeed, previous attempts at weak wave turbulence analysis in three dimensions [8, 9, 10, 11] have fallen short of providing a definitive prediction for the energy spectrum. While observations emphasize the central role dispersive internal gravity waves play in natural processes like the ocean’s climate cycle [12], theoretically and experimentally decoupling these from the evolution of slow modes, degrees of freedom with vanishing frequency, proves difficult. Slow modes, such as shear and domain modes in 2D and 3D and vortical modes in 3D, are non-linearly interacting with the dispersive waves, and tend to occupy a prominent part of the energy [13, 14, 15]. Compared to 3D, the 2D problem holds a few advantages: it has no vortical modes, it is cheaper for direct numerical simulation and its, recently derived [16], kinetic equation takes a particular simple form due to the existence of a second quadratic invariant. All suggest that 2D Boussinesq can serve as fertile ground for theoretical and experimental investigation. The relevance of the 2D internal gravity wave description extends to both experimental settings, such as long water tanks, and natural phenomena like internal tides around isolated 1D topographical features like the Hawaiian ridge [17]. Moreover, understanding the 2D case can serve as a foundation for comprehending the 3D dynamics. In our previous work, [16], we derived the isotropic (integrated over the angle) part of the turbulent energy spectrum of weakly interacting 2D internal gravity waves. Here we complete the derivation of the angular part, when shear modes are removed. In a direct approach, the kinetic equation diverges around the turbulent solution we find. This, we claim, indicates a meaningful limit on the applicability range of the kinetic description due to the existence of slow modes in the original dynamical equations. We address this singularity using a regulator and discuss its relevance to practical observations. We interpret the solution in terms of a angular dependent radial flux. Along a curve of steady generalized solutions we find that the turbulent solution is the unique one with a non-zero radial flux. We believe that our methodology can be extended to 3D and to waves in plasma.

Figure 1: Logarithm of the energy spectrum $e_{k}=C_{0}K^{-3}(\omega_{k}/N)^{-2}=C_{0}K^{-3}\cos\theta_{k}^{-2}$ on the regularized $(k_{x},k_{z})$ plane, outside an angular source centered at radius $K_{f}$ and emitting a radial energy flux $\Pi^{r}(\theta_{k})=\Pi_{0}\sin(\theta_{k})\hat{k}/K_{f}$.

The two-dimensional Boussinesq equations restricted to a vertical $xz$-plane can be written as

$$\Delta\psi_{t}+\left\{\psi,\Delta\psi\right\} =-N^{2}\eta_{x}$$ (1)

$$\eta_{t}+\left\{\psi,\eta\right\} =\psi_{x}.$$

Here $z$ is the vertical and $x$ is the horizontal coordinate with corresponding velocities $w$ and $u$, $\psi$ is a stream function such that $(\psi_{x},\psi_{z})=(w,-u)$ and $-\Delta\psi$ is the vorticity, $\eta$ is the vertical displacement, $N$ the constant buoyancy frequency and $\left\{g,f\right\}=\partial_{x}g\partial_{z}f-\partial_{z}g\partial_{x}f$. The vertical buoyancy force $b=-N^{2}\eta$ opposes vertical displacements and derives from a consideration of potential energy in the presence of gravity and non-uniform density. (1) has two quadratic invariants: the total energy $E=\int d\mathbf{x}(-\psi\Delta\psi+N^{2}\eta^{2})$ and the pseudomomentum $P=\int d\mathbf{x}\,\eta\Delta\psi$.

Wave mode expansion and the kinetic equation. We consider a periodic domain $\mathbf{x}\in\left[0,L\right]^{2}$ and expand the fields $(\psi,\eta)=\sum_{\alpha}Z_{\alpha}(t)g_{\alpha}(\mathbf{x})$ in terms of linear wave modes, where $Z_{\alpha}(t)$ are complex scalar wave amplitudes and the $g_{\alpha}(\mathbf{x})$ are eigenvector functions for the linear part of (1),

$$-i\omega_{\alpha}\begin{pmatrix}-\Delta&0\\ 0&N^{2}\end{pmatrix}g_{\alpha} ={N^{2}}\begin{pmatrix}0&1\\ 1&0\end{pmatrix}\partial_{x}g_{\alpha}.$$ (2)

The multi-index $\alpha=(\sigma,\mathbf{k})$ combines branch choice $\sigma=0,\pm 1$ and wave number $\mathbf{k}\in(2\pi\mathbb{Z}/L)^{2}$. In polar coordinates, $\mathbf{k}=K(\cos\theta_{k},\sin\theta_{k})$, the dispersion relation is

$$\omega_{\alpha=(\sigma,\mathbf{k})}=\sigma N\cos\theta_{k}.$$ (3)

The choice $\sigma=\pm 1$ corresponds to horizontally right-going or left-going waves and $\theta_{k}=\pm\pi/2$, where $\omega_{k}=0$, characterizes the shear modes. The expansion diagonalizes the energy $E=\sum_{\alpha}E_{\alpha}=\sum_{\alpha}Z_{\alpha}Z_{\alpha}^{*}$ and yields the exact modal wave energy equation

$$\dot{E_{\alpha}}=\sum_{\beta,\gamma}V_{\alpha}^{\beta\gamma}\operatorname{Re}\left(Z^{*}_{\alpha}Z^{*}_{\beta}Z^{*}_{\gamma}\right),$$ (4)

with the interaction coefficients $V_{\alpha}^{\beta\gamma}$ given in [16]. The wave expansion diagonalizes the pseudo-momentum as well, $P=\sum_{\alpha}\sigma_{\alpha}N^{-1}KE_{\alpha}$, so horizontally right-going waves carry positive pseudo-momentum. The kinetic equation that describes the slow evolution of the averaged energy density $e_{\alpha}=\overline{E_{\alpha}}$, where the overbar denotes averaging over an initial Gaussian statistical ensemble,

$$\overline{Z_{\beta}^{*}(0)Z_{\alpha}(0)}=\delta_{\alpha\beta}e_{\alpha}(0),$$ (5)

was derived in [16] and is given by

$$\dot{e}_{\alpha} =St_{\alpha}(e_{\alpha})$$ (6)

$$={\pi}\!\!\int\limits_{\omega_{\alpha\beta\gamma}}\!\!\omega_{\alpha}\,\Gamma_{\alpha\beta\gamma}^{2}(\omega_{\alpha}e_{\beta}e_{\gamma}+\omega_{\beta}e_{\alpha}e_{\gamma}+\omega_{\gamma}e_{\alpha}e_{\beta}).$$

To derive the kinetic equation, which is first non trivial closure to modal energy evolution (4) the joint kinetic limits of big box and long nonlinear times, $L\rightarrow\infty$ and $t\omega\rightarrow\infty$, are taken. So the discrete sum in (4) is replaced by an integral over the resonant manifold:

$$\int\limits_{\omega_{\alpha\beta\gamma}}\!\!\!\!\!=\!\int d\beta d\gamma\,\delta(\omega_{\alpha}+\omega_{\beta}+\omega_{\gamma})\delta\left(\mathbf{k}_{\alpha}+\mathbf{k}_{\beta}+\mathbf{k}_{\gamma}\right),$$ (7)

where $\int\!\!d\alpha=\sum_{\sigma=\pm 1}\int\!\!d\mathbf{k}$. The interaction coefficients are proportional to the frequencies, i.e.,$V^{\beta\gamma}_{\alpha}/\omega_{\alpha}=V_{\beta}^{\alpha\gamma}/\omega_{\beta}=V_{\gamma}^{\alpha\beta}/\omega_{\gamma}=\Gamma_{\alpha\beta\gamma}$, where

$$\Gamma_{\alpha\beta\gamma}=\frac{1}{\sqrt{8}}\sum_{\chi=\alpha,\beta,\gamma}\!\!\!\!\sin\theta_{\chi}\!\!\!\sum_{\chi=\alpha,\beta,\gamma}\!\!\!\sigma_{\chi}K_{\chi},$$ (8)

which allows the relatively compact and simple kinetic description given by (6). In the vicinity of small frequencies the kinetic equation must be interpreted carefully. In the discrete sum on the lattice, (4), shear modes are well separated from waves with non zero frequency, however as $L\rightarrow\infty$ the collision integral (6) includes integration arbitrarily close to slow modes. While, slow modes cannot be created by resonant interactions, as $\omega_{\alpha}\rightarrow 0$ the kinetic equation needs to include off diagonal correlators apart from (5), such as $\overline{Z_{\alpha}^{2}}$, which would yield a very complicated description. Such correlators oscillate with frequency $\sim 2\omega_{\alpha}$ and over long times can be neglected as long as $\omega_{\alpha}\!>\!\epsilon\!>\!0$ [16]. From that perspective, anisotropic systems with dispersion relation that vanish on a curve rather than on a point, pose a special problem for kinetic description. This does not affect the derivation in [16] of the homogeneous part of the energy spectrum which depends on the homogeneity degree of the interaction coefficient, frequency and dimension. However, as we continue towards the derivation of the angular dependence of the energy spectrum we anticipate trouble along the $k_{z}$ axis where $\omega_{k}=0$. To overcome this obstacle, without trading the simplicity of the kinetic equation in a cumbersome version, we consider the kinetic equation on the $(k_{x},k_{z})$ plane with a small sector of width $2\delta$ around the $k_{z}$ axis removed and then investigate the limit as $\delta\rightarrow 0$ numerically; see Figure 1. To obtain our results we numerically integrate the collision integral in Mathematica, our notebook is attached to the supplementary material.

Boundary conditions and steady solutions of the kinetic equation. The frequency and the interaction coefficients (8) are symmetric under dilation transformation. That is, for $\lambda>0$, $\omega_{\left(\sigma,\lambda\mathbf{k}\right)}=\omega_{\left(\sigma,\mathbf{k}\right)}$ and $\Gamma\left(\lambda\mathbf{k}_{\alpha},\lambda\mathbf{k}_{\beta},\lambda\mathbf{k}_{\gamma}\right)=\lambda\Gamma\left(\mathbf{k}_{\alpha},\mathbf{k}_{\beta},\mathbf{k}_{\gamma}\right)$. Assuming the steady spectrum $e_{\alpha}$ is symmetric with respect to dilation as well, fixes the homogeneous separable form $e_{\alpha}=C_{0}K^{w}e_{\alpha}^{\theta}\left(\theta_{k}\right)$, where $C_{0}>0$ is a constant, usually referred to as the Kolmogorov constant. We consider the case that energy is symmetrically distributed between horizontally left-going and right-going waves, $e_{(-,\mathbf{k})}=e_{(+,\mathbf{k})}$ and then the pseoudomomentum is zero on average $\overline{P}=0$, so that $e_{\alpha}^{\theta}=e_{\theta}$. Concentrating further on solutions that correspond to a symmetric forcing about the $k_{x}$ and $k_{z}$ axes, we can equivalently represent the solution in the form motivated by $\omega_{k}/N\sigma=\cos\theta_{k}$

$$e_{\alpha}=C_{0}K^{w}(\cos\theta_{k})^{2f_{w}(\theta_{k})},$$ (9)

where $\cos\theta_{k}^{2f_{w}}\equiv\left(\left(\cos\theta_{k}\right)^{2}\right)^{f_{w}}$, so the spectrum is non negative. The homogeneous spectrum turns the collision integral (the RHS of the kinetic equation (6)) to a homogeneous function and hence separable in angle and wave number amplitude

$$St_{\alpha}\equiv K^{2w-2w_{0}-2}C_{0}^{2}St_{\theta}\left(w,f_{w},\theta_{k}\right),$$ (10)

where $w_{0}=-3$ was found in a recent work by the authors [16]. The explicit form for the angular part of the collision integral, $St_{\theta}$, is found by parametrization of the resonant manifold and is presented in the supplementary material. The separability of the collision integral suggests the existence of a one dimensional continuous family of solutions parameterized by $w$. To choose the relevant physical solution we need to assign boundary conditions for the steady equation, which we do through the energy flux. Energy is an integral of motion so the kinetic equation (6) can be written formally as a continuity equation

$$\dot{e_{\alpha}}\left(t\right)+\text{div}\Pi_{\alpha}\left(t\right)=0,$$ (11)

where $\Pi_{\alpha}\left(t\right)=\left(\Pi^{r}_{\alpha}(t),\Pi^{\theta}_{\alpha}(t)\right)$ is a two dimensional flux, its divergence equals the minus of the collision integral

$$\text{div}\Pi_{\alpha}(t)=-St_{\alpha}(e_{\alpha}),$$ (12)

so a steady solution has zero divergence flux. Due to the gauge freedom $\Pi_{\alpha}+\nabla\times A_{\alpha}\hat{z}$, with a scalar $A_{\alpha}$, any steady solution corresponds to infinitely many choices of fluxes. However, fixing a flux fixes a unique solution. In steady turbulence existing between forcing and dissipation that are widely separated in scale and symmetric about the $k_{x}$ and $k_{z}$ axes, far away from the source we do not expect existence of an angular flux. This also aligns with the fact that there are no resonant interaction among triads $(\mathbf{k},\mathbf{p},\mathbf{q})$ on the same spherical shell $K=K_{p}=K_{q}$ but with different angles $(\theta_{k},\theta_{p},\theta_{q})$. Having only radial flux (12), reads

$$\text{div}\Pi_{\alpha} =K^{-1}\partial_{K}(K\Pi^{r}_{\alpha})\equiv-St_{\alpha}(e_{\alpha}).$$ (13)

multiplying by $K$ and integrating $\int_{K}^{\infty}dK$, the homogeneous collision integral (10) assuming $w<w_{0}$ yields the flux

$$K\Pi^{r}_{\alpha}=-\frac{K^{2\left(w_{0}-w\right)}}{2\left(w_{0}-w\right)}C_{0}^{2}St_{\theta}\left(w,f_{w},\theta_{k}\right).$$ (14)

Integrating $\int_{0}^{K}dK$ would yield the same expression for $w>w_{0}$. Here it is clear the pair $(w_{0},f_{\omega_{0}})$ corresponds to a special solution - the unique solution that gives a non zero radial flux with zero divergence. As $w\rightarrow w_{0}$, due to the diverging denominator in (14), the radial flux is given by

$$\Pi^{r}_{\alpha}=-\frac{1}{2K}\frac{d}{dw}St_{\theta}\left(w,f_{w_{0}},\theta_{k}\right)\mid_{w=w_{0}}.$$ (15)

Note here that the $w$ derivative is taken after the angular part $f_{w}=f_{w_{0}}$ has been fixed in the angular collision integral. Isotropic wave turbulence [18] is a particular case of (15) where $St_{\theta}\left(w,f_{w_{0}},\theta_{k}\right)=St_{\theta}\left(w\right)$ is independent of the angle. The integral over the angle gives the constant non-zero total radial flux $\int_{0}^{2\pi}K\Pi^{r}_{\alpha}\left(\theta\right)d\theta=\Pi_{0}$. Having the explicit power law of the radial part of the turbulent solution, $w_{0}=-3$ from previous work, we now turn to find the power of the angular part of the solution, $f_{w_{0}}$.

Obtaining the turbulent solution on the regularized plane. While the radial part of the energy spectrum is insensitive to the validity range of the kinetic equation as the frequency $\omega_{k}\rightarrow 0$, the angular part is. To obtain the angular part of turbulent energy spectrum we regularize the angular part of the collision integral by removing an angular sector of width $2\delta$. That is, the angular integration over the incoming angles on the resonant manifold, (7), is performed over $\theta_{\beta},\theta_{\gamma}\in\left[-\pi/2+\delta,\pi/2-\delta\right]$. In terms of frequencies, this restricts $\cos\theta_{\beta},\cos\theta_{\gamma}\in\left[-1,-\epsilon\right]\cup\left[\epsilon,1\right]$, for small positive $\epsilon$, $\epsilon\sim\delta$, so we use $\delta$ in both contexts. Denote the regularized angular collision integral by $St_{\theta}\left(w,f_{w},\theta_{k},\delta\right)$, we look for the angular function $f_{w_{0}}$ that is the zero of the regularized functional in the limit $\delta\rightarrow 0$:

$$\lim_{\delta\rightarrow 0}St_{\theta}\left(w_{0},f_{w_{0}},\theta_{k},\delta\right)=0.$$ (16)

The Ansatz (9) for the energy spectrum makes the main term in the collision integral (6) proportional to

$$\omega_{\beta}e_{\alpha}e_{\gamma}+\omega_{\gamma}e_{\beta}e_{\alpha}+\omega_{\alpha}e_{\beta}e_{\gamma}\propto\!\!\!\!\!\sum_{\chi=\alpha,\beta,\gamma}\!\!\!\!\!\sigma_{\chi}K_{\chi}^{-w}\left(\cos\theta_{\chi}\right)^{1-2f_{w}},$$ (17)

suggesting the curve $w=2f_{w}-1$ is special. Indeed, we find numerically that $\lim_{\delta\rightarrow 0}St_{\theta}(-3,-1,\theta_{k},\delta)=0$. Note that (16) is just the principal value of the collision integral around divergences coming from the curve where the frequency vanishes. These divergences cancel each other and yield a zero of the collision integral, which we refer to as a generalized zero. So (16) is equivalent to integrating the collision integral over the regular domain with the energy spectrum $e_{\alpha}=K^{-3}(\omega_{k}+i\delta)^{2f_{w}}$ and then taking $\delta\rightarrow 0$.

The turbulent energy spectrum and its flux. The energy spectrum of weakly interacting 2D internal gravity waves in the inertial range between large scale forcing and small scale dissipation, when shear modes are removed and the forcing is symmetric about the $(k_{x},k_{z})$ axes and does not produce pseudo-momentum on average, is

$$e_{\left(\sigma,\mathbf{k}\right)}=C_{0}K^{-3}(\cos\theta_{k})^{-2}=C_{0}K^{-1}k_{x}^{-2}.$$ (18)

Actually, for a finite total energy flux $\Pi_{0}$ the Kolmogorov constant $C_{0}$ is given by

$$C_{0}=\sqrt{\Pi_{0}/2\pi}\,\delta^{3/2},$$ (19)

so as $\delta\to 0$ a finite flux $\Pi_{0}$ can be maintained with a smaller and smaller energy spectrum. Conversely, at fixed $C_{0}$ the flux diverges as $\Pi_{0}\sim\delta^{-3}$, which is a neat indication for the breakdown of the kinetic description given by (6) for interactions with and among slow modes. The spectrum is plotted on the regularized domain in Figure 11. This spectrum is a density over $(k_{x},k_{z})$ so in order to compare it with the oceanic Garrett–Munk spectrum [7] we use the full, non-hydrostatic dispersion relation (3) in order to rewrite (18) as a spectral density over positive $(\omega,k_{z})$. This yields

$$K^{-1}k_{x}^{-2}\,\mathrm{d}k_{x}\mathrm{d}k_{z} =N\omega^{-2}k_{z}^{-2}\,\mathrm{d}\omega\mathrm{d}k_{z}.$$ (20)

This agrees with the shape of the Garrett-Munk spectrum in the limit of large $k_{z}$ and $\omega$. The practical relevance of (18) for experiments and numerical observation with $\delta\neq 0$ is discussed in the appendix. We find numerically the angular dependent turbulent flux (15) normalized by the total flux to be,

$$\lim_{\delta\rightarrow 0}K\Pi_{\alpha}/\Pi_{0}=\sin\theta_{\alpha},$$ (21)

which implies that the energy transferred within an angular sector is linear in frequency

$$\lim_{\delta\rightarrow 0}\int_{0}^{\theta}d\theta^{\prime}K\Pi_{\alpha^{\prime}}/\Pi_{0}=\cos\theta_{k}=\omega_{k}/\sigma N.$$ (22)

The contribution to the energy flux at frequency $\alpha=(\sigma_{\alpha},\mathbf{k_{\alpha}})$ can be divided into three classes of triad $(\alpha,\beta,\gamma)$ interactions wrt to the corresponding frequency branches $(\sigma_{\alpha},\sigma_{\beta},\sigma_{\gamma})$. Fixing $\sigma_{\alpha}=+1$, we find that interactions among waves belonging to the same branch with $\sigma_{\alpha}=\sigma_{\beta}=\sigma_{\gamma}$ transfer energy backwards in wave number $K$. That is, the contribution to the flux $\Pi_{\alpha}$ from interactions among the same branch is negative: $\left.{\Pi_{\alpha}}\right|_{(+,+,+)}\!<\!0$ for all $\omega_{k}$, the contribution of triads $\sigma_{\alpha}\neq\sigma_{\beta}=\sigma_{\gamma}$ is negligible and the contribution of the triads $\sigma_{\alpha}=\sigma_{\beta}\neq\sigma_{\gamma}$ and $\sigma_{\alpha}=\sigma_{\gamma}\neq\sigma_{\beta}$ is positive to a forward cascade, so that $\left.{\Pi_{\alpha}}\right|_{(+,-,+)}\!\!=\!\!\left.{\Pi_{\alpha}}\right|_{(+,+,-)}\!\!>0$. Furthermore, the fraction of energy transported backward is half of the energy transported forward, $\left.{\Pi_{\alpha}}\right|_{(+,+,+)}/\left.{\Pi_{\alpha}}\right|_{(+,-,+)}=-\frac{1}{2}$, and is transported by the set of waves belonging to the same branch on which the pseudo-momentum is sign definite. This suggests an existence of an inverse energy cascade of internal gravity waves when interactions among waves of different branches is suppressed, for example in the neighborhood of a source emitting horizontally unidirectional internal waves, such as the Hawaiian ridge [17]. The total angular flux $\Pi_{\alpha}$, and the fluxes restricted to the different triads $\left.{\Pi_{\alpha}}\right|_{(+,+,+)}$, $\left.{\Pi_{\alpha}}\right|_{(+,-,+)}$ and $\left.{\Pi_{\alpha}}\right|_{(+,+,-)}$ are plotted in Figure 4.

Figure 2: The normalized radial, angular dependent but scale independent, flux $\Pi_{\alpha}$, and the fluxes restricted to the different triads $\Pi_{\alpha}\mid_{(+,+,+)}$ and $\Pi_{\alpha}\mid_{(+,-,+)}+\Pi_{\alpha}\mid_{(+,+,-)}$, evaluated at $\epsilon=10^{-5}$. The linear dashed green curve presents the normalized flux contained in angular section $\theta$ given by (22).

Curve of solutions. In the vicinity of the turbulent solution $(w_{0},f_{w_{0}})$, we find the curve $(w,f_{w}=w+1)$ suggested by (17) is a curve of generalized zeroes of the collision integral for $3.1\lesssim w\lesssim 2$, see Figure 4. Notably, there exists a trivial solution of the form (9), with $w_{T}=0$ and $f_{w_{T}}=0$, which is the equilibrium solution of equipartition of energy. On the plane $(w,f_{w})$ this solution is isolated from the curve passing through the turbulent solution $(w_{0},f_{w_{0}})$.

Figure 3: Curve of solutions $(w,f_{w})$ in the vicinity of the turbulent point $(w_{0},f_{w_{0}})=(-3,-1)$ evaluated at $\delta=10^{-5}$. Different colors correspond to different angles, in the region $3.1\lesssim w\lesssim 2$ these collapse on the linear curve $f_{w}=w+1$.

Previously known solutions to anisotropic kinetic equations. Previous kinetic description of 3D internal gravity waves considered a narrow spectral range where the horizontal wave number component is small compared to the vertical $k_{x}\ll k_{z}$, also known as the hydrostatic approximation. In this approximation the frequency and interaction coefficients are bi-homogeneous functions of $k_{x}$ and $k_{z}$ and hence one can look for bi-homogeneous solutions $e_{\alpha}=C_{0}k_{x}^{2w_{x}}k_{z}^{2w_{z}}$. One pair of power-laws $(2w_{x},2w_{z})$ can be derived using Zakharov—Kuznetsov transformation [19, 18] for the kinetic equation in the hydrostatic approximation, with $2w_{x}=-w_{\Gamma_{x}}-1-\frac{w_{\omega_{x}}}{2}$ and $2w_{z}=-w_{\Gamma_{z}}-1-\frac{w_{\omega_{z}}}{2}$. Here $w_{\Gamma_{x/z}}$ and $w_{\omega_{x/z}}$ are the homogeneous degrees with respect to to $k_{x}/k_{z}$ of the interaction coefficient and frequency. Our spectrum (18) is peaked near $k_{x}=0$, which suggests that the hydrostatic approximation is accurate and in this approximation our spectrum is $e_{\alpha}\sim k_{x}^{-2}\left|k_{z}\right|^{-1}$. Notably, this is different from the bi-homogeneous spectrum obtained using Zakharov transformation in the hydrostatic approximation, which is $e_{\alpha}\sim\left|k_{x}\right|^{-1/2}\left|k_{z}\right|^{-5/2}$. While both spectra have the same isotropic radial part, the angular part differs.

Conclusion and perspectives. Our work suggests that a divergence of the collision integral on the set of vanishing frequency, which is generically large for anisotropic dispersive equations, is meaningful and indicates limited applicability of the kinetic description due to the non-linear coupling of the dynamical equations to slow modes. From that perspective, the kinetic equation of weakly interacting 2D internal gravity waves has a practically relevant generalized solution $e_{\left(\sigma,\mathbf{k}\right)}=C_{0}K^{-3}(N^{-1}\omega_{k})^{-2}$. This spectrum holds in the inertial range between large-scale forcing and small-scale dissipation when shear modes are excluded, and the forcing is symmetric about the $(k_{x},k_{z})$ axes without producing pseudo-momentum on average. For vertical scales short compared to the depth of the ocean and frequencies large compared to the Coriolis frequency it aligns with the oceanic Garrett–Munk spectrum. Apparently, this is the first time that a theoretical internal wave spectrum agrees with the Garrett–Munk spectrum without relying on the hydrostatic approximation. The simplicity of our numerical results, such as (18) and (15), implies potential for analytical derivation. A natural progression of this research would be adding rotation (Coriolis force) to the equations, which would serve as a physical regulator by decoupling slow modes from waves due to the dispersion relation gap. We expect the spectrum to adapt its form to the modified dispersion relation $\omega_{(\sigma,\mathbf{k})}=\sigma N\cos\theta_{k}\sqrt{1+f^{2}N^{-2}\tan^{2}\theta_{k}}$, where $f$ is the Coriolis parameter. We believe our methodology can be extended to waves in plasma and to 3D internal gravity waves.

Acknowledgments. We thank Gregory Falkovich for useful discussions. Special thanks are owed to Vincent Labbare for his invaluable insights and extensive feedback, which significantly contributed to the development of this manuscript. This work was supported by the Simons Foundation and the Simons Collaboration on Wave Turbulence. OB acknowledges additional financial support under ONR grant N00014-19-1-2407 and NSF grant DMS-2108225. MS acknowledges additional financial support from the Schmidt Futures Foundation.

Appendix A - Practical relevance for experiments and numerical observation. For a finite $\delta$, $f_{w_{0}}=f_{w_{0}}\left(\theta_{k},\delta\right)$ depends on $\delta$ and on the angle $\theta_{k}$. For $\delta\ll 1$ we find numerically

$$f_{w_{0}}(\theta_{k},\delta)=-1-m(\theta_{k})\delta+O(\delta^{2})$$ (23)

is approaching $-1$ linearly in $\delta$ from the left with a minimal slope at $\theta_{k}=\pi/4$ so that $m(\theta)\sim-13.51+60.9(\cos\theta_{k}-\cos(\pi/4))^{2}+...$. Figure 4 demonstrates (23) for various angles and shows the slope $m(\theta)$ in an inset. The practical relevance of (23) is the broadening in angle of the measured energy spectrum, with respect to the pure theoretical prediction, in finite size systems or in DNS on a periodic domain where shear modes are removed. While the self similar slope $w_{0}=-3$ is insensitive to a small regulator, the angular spectrum will deviate from $\cos\theta^{-2}$ according to (23). In a DNS of the dynamical equations $\eqref{eq:2D Boussinesq}$ on a periodic domain, with minimal wave numbers $k_{min}=1$ and maximal wave number $k_{max}$ set by resolution $M$, one can consider the continuous approximation of the regulator $\delta\sim\arcsin(1/M)$. Then the measured energy spectrum normalized by the theoretical solution $e_{\alpha}C_{0}^{-1}K^{3}\cos\theta_{k}^{2}\sim\omega^{-2\delta m(\theta_{k})}$ will have an angular width of $1-\omega^{-2\delta m(\theta_{k}^{0})}$, with $\cos\theta_{k}^{0}$ small but in the inertial range. So for $M=4096$ and $\cos\theta_{k}^{0}=0.1$ , $1-\omega^{-2\delta m(\theta_{k}^{0})}\sim 0.06$.

Figure 4: The root $f_{w_{0}}(\theta_{k},\delta)$ of the angular collision integral $St_{\theta}(w_{0}=-3,f_{w_{0}},\theta_{k},\delta)$ for various angles. For small $\delta$, the inset shows the linear in $\delta$ correction to $f_{w}=-1-m(\theta)\delta+O(\delta^{2})$.