# HG changeset patch
# User vlabarre <vincent.labarre@oca.eu>
# Date 1698750479 -3600
#      Tue Oct 31 12:07:59 2023 +0100
# Node ID 0c5cd11125dd937949e48fb8cb6106b5a586dfce
# Parent  885f53d75aeceb74a6bd3bb18a2634d9e912c6d5
after Sergey's corrections

diff --git a/2022strat_polo_proj/input/main.tex b/2022strat_polo_proj/input/main.tex
--- a/2022strat_polo_proj/input/main.tex
+++ b/2022strat_polo_proj/input/main.tex
@@ -143,12 +143,6 @@
 \setlength\parindent{0pt}
 
 
-% For review
-\newcommand{\Add}[1]{{\color{blue}#1}}
-\newcommand{\Remove}[1]{{\color{red}\sout{#1}}}
-\newcommand{\Comment}[1]{{\color{green}#1}}
-
-
 \begin{document}
 
 \title{Internal gravity waves in stratified flows with and without vortical modes}
@@ -347,9 +341,9 @@
 \begin{equation}
 E(k_h, k_z) \sim k_h^{-3/2} k_z^{-3/2} ~~~~ \text{and} ~~~~ E(\omega, k_z) \sim \omega^{-3/2} k_z^{-2},
 \end{equation}
-\Add{where the change of coordinates is defined by the dispersion relation $E(\omega,
-k_z) = E(k_h,k_z) \left( \partial \ok / \partial k_h \right)^{-1}$.} Yet, it was noted
-that this candidate does not satisfy the requirement of locality, i.e. the collision
+where the change of coordinates is defined by the dispersion relation $E(\omega, k_z) =
+E(k_h,k_z) \left( \partial \ok / \partial k_h \right)^{-1}$. Yet, it was noted that
+this candidate does not satisfy the requirement of locality, i.e. the collision
 integral diverges on it. In the other words, it is not a valid mathematical solution of
 the kinetic equation. Later, it was shown that power law solutions $E(k_h,k_z) \sim
 k_h^{-\alpha_h} ~ k_z^{-\alpha_z}$ have convergent collision integral's contributions
@@ -384,8 +378,8 @@
 \label{eq:Isotropic}
 E_{\rm 1D}(k) \sim \epsK^{2/3} k^{-5/3}
 \end{equation}
-is expected. It happens when $\kb < \ko \ll k \ll k_{\eta}$ where $\ko$ is the
-Ozmidov wave-vector and $k_{\eta}$ is the Kolmogorov wave-vector.
+is expected. It happens when $\kb < \ko \ll k \ll k_{\eta}$ where $\ko$ is the Ozmidov
+wave-vector and $k_{\eta}$ is the Kolmogorov wave-vector.
 
 
 The present study first deals with the existence and properties of a weak internal
@@ -422,8 +416,7 @@
 \section{Methods}
 \label{sec:methods}
 
-We use the 3D Navier-Stokes equations under the Boussinesq approximation: \Remove{with
-added hyperviscosity and hyperdiffusion:}
+We use the 3D Navier-Stokes equations under the Boussinesq approximation:
 \begin{align}
 \label{eq:Continuity}
 \bnabla \cdot \vv &= 0 \\
@@ -436,81 +429,74 @@
 where $(x,y,z)$ will represent the three spatial coordinates in the cartesian frame
 $(O, \eex, \eey, \eez)$, $\eez$ is the stratification axis, $\vv=(v_x, v_y, v_z)$ is
 the velocity, $b$ the buoyancy, $p$ the total kinematic pressure, $N$ the \bv
-frequency, $\nu$ the viscosity, $\kappa$ the diffusivity\Remove{, $\nu_4$ the
-hyperviscosity, $\kappa_4$ the hyperdiffusivity,} and $\ff$ is the velocity forcing.
-\Add{The buoyancy is defined as $b = -g \rho' / \rho_0$, where $g$ is the acceleration
-due to gravity, $\rho_0$ is the average density of the fluid at $z=0$, and $\rho'$ is
-the density perturbation with respect to the average linear density profile
-$\bar{\rho}(z) = \rho_0 + (\mathrm{d}\bar{\rho}/\mathrm{d}z)z$.} The Schmidt number $Sc
+frequency, $\nu$ the viscosity, $\kappa$ the diffusivity and $\ff$ is the velocity
+forcing. The buoyancy is defined as $b = -g \rho' / \rho_0$, where $g$ is the
+acceleration due to gravity, $\rho_0$ is the average density of the fluid at $z=0$, and
+$\rho'$ is the density perturbation with respect to the average linear density profile
+$\bar{\rho}(z) = \rho_0 + (\mathrm{d}\bar{\rho}/\mathrm{d}z)z$. The Schmidt number $Sc
 = \nu/\kappa$ is fixed to one.
 
-\Add{We simulate forced dissipative flows} with the pseudo-spectral solver
-\texttt{ns3d.strat} from the FluidSim software \cite{mohanan_fluidsim_2019} (an
-open-source Python package of the FluidDyn project \cite{fluiddyn} using Fluidfft
-\cite{fluidfft} to compute the Fast Fourier Transforms). \Add{The forcing, which will
-be described in details at the end of this section and in
-Appendix~\ref{appendix:forcing}, is computed in spectral space such that the kinetic
-energy injection rate $P_K$ is constant and equal to unity. The physical input
-parameters are the \bv frequency $N$ and the diffusive coefficients $\nu=\kappa$, but
-in practice, we identify our simulations with the couple $(N,\,\R_i)$, where $\R_i
-\equiv P_K / (\nu N^2)$ is the input buoyancy Reynolds number.} The turbulent
-non-dimensional numbers characterizing the statistically stationarity flow are the
-horizontal turbulent Froude number \Remove{, the Reynolds number,} and the buoyancy
-Reynolds number \cite[]{brethouwer_scaling_2007} that are respectively
+We simulate forced-dissipated flows with the pseudo-spectral solver \texttt{ns3d.strat}
+from the FluidSim software \cite{mohanan_fluidsim_2019} (an open-source Python package
+of the FluidDyn project \cite{fluiddyn} using Fluidfft \cite{fluidfft} to compute the
+Fast Fourier Transforms). The forcing, which will be described in details at the end of
+this section and in Appendix~\ref{appendix:forcing}, is computed in spectral space such
+that the kinetic energy injection rate $P_K$ is constant and equal to unity. The
+physical input parameters are the \bv frequency $N$ and the diffusive coefficients
+$\nu=\kappa$, but in practice, we identify our simulations with the couple
+$(N,\,\R_i)$, where $\R_i \equiv P_K / (\nu N^2)$ is the input buoyancy Reynolds
+number. The turbulent non-dimensional numbers characterizing the statistically
+stationarity flow are the horizontal turbulent Froude number and the buoyancy Reynolds
+number \cite[]{brethouwer_scaling_2007} that are respectively
 \begin{equation}
-\label{eq:FhR} F_h = \frac{\epsK}{{U_h}^2 N} \Remove{, ~~~~ Re =
-\frac{U_h^4}{\nu \epsK}} ~~~~ \text{and} ~~~~ \R = \frac{\epsK}{\nu N^2},
+\label{eq:FhR} F_h = \frac{\epsK}{{U_h}^2 N} ~~~~ \text{and} ~~~~ \R = \frac{\epsK}{\nu N^2},
 \end{equation}
-\Add{where $\epsK$ is the kinetic energy dissipation rate and $U_h$ the rms of the
-horizontal velocity.} \Add{Note that the turbulent Reynolds number is given by $Re = \R
-/ F_h^2$.} We also compute the buoyancy and the Ozmidov wave-vectors $\kb \equiv N/U_h$
+where $\epsK$ is the kinetic energy dissipation rate and $U_h$ the rms of the
+horizontal velocity. Note that the turbulent Reynolds number is given by $Re = \R /
+F_h^2$. We also compute the buoyancy and the Ozmidov wave-vectors $\kb \equiv N/U_h$
 and $\ko \equiv \sqrt{N^3 / P_K}$. All the quantities presented in this manuscript are
 computed from averaging when stationarity is reached. A list of the simulations, with
 relevant parameters and physical quantities, is given in
 Appendix~\ref{table-better-simuls}.
 
 In this study, we consider a periodic domain of horizontal size $L_x = L_y = L_h = 3$
-\Add{and vertical size} \Remove{of the domain,} $L_z$ \Remove{, is varied depending on
-the value of the \bv frequency}. We note $(n_x, n_y, n_z)$ the numbers of collocations
-points in the three spatial directions, with $n_x = n_y \equiv n_h$. \Add{We chose
-$n_z$ in order to have an isotropic mesh in physical space, i.e. $L_z/n_z = L_h/n_h$.}
-\Add{We decrease $L_z$ with $N$. Typically, $L_z \propto 1/N$, while $L_h$ is kept
-constant for all simulations. More precisely, the aspect ratio $L_z/L_h$ is $1/2$ for
-$N \leq 20$, $1/4$ for $N \leq 60$, and $1/8$ for $N \geq 80$. This choice is motivated
-by the fact that the unforced and undissipated Boussinesq equations are self-similar in
-the limit $F_h \rightarrow 0$, with similarity variable $zN/U$, where $U$ is the
-typical velocity \cite{billant_self-similarity_2001}. In that way, we simulate few
-layers (of height $L_b = U/N$) for all our simulations.}
+and vertical size $L_z$. We note $(n_x, n_y, n_z)$ the numbers of collocations points
+in the three spatial directions, with $n_x = n_y \equiv n_h$. We chose $n_z$ in order
+to have an isotropic mesh in physical space, i.e. $L_z/n_z = L_h/n_h$. We decrease
+$L_z$ with $N$. Typically, $L_z \propto 1/N$, while $L_h$ is kept constant for all
+simulations. More precisely, the aspect ratio $L_z/L_h$ is $1/2$ for $N \leq 20$, $1/4$
+for $N \leq 60$, and $1/8$ for $N \geq 80$. This choice is motivated by the fact that
+the unforced and undissipated Boussinesq equations are self-similar in the limit $F_h
+\rightarrow 0$, with similarity variable $zN/U$, where $U$ is the typical velocity
+\cite{billant_self-similarity_2001}. In that way, we simulate few layers (of height
+$L_b = U/N$) for all our simulations.
 
-\Add{To reduce the computational costs, we first simulate the transient state with a
-coarse resolution with $n_h = 320$.}
-\Add{For these simulations, two hyperdiffusive terms $-\nu_4 \nabla^4 \vv$ and
-$-\kappa_4\nabla^4{b}$ are added to (\ref{eq:Impulsion}) and (\ref{eq:Buoyancy}),
-respectively, in order to keep the dissipative range in the simulated scales and avoid
-thermalisation at small scales.}
-\Add{Once a statistically steady-state is reached, we increase the resolution of the
-simulation while decreasing hyper-viscosity. Then the simulation is run until reaching
-a new statistically steady state. The previous step is repeated until reaching the
-highest resolution.}
-We measure the turbulent kinetic dissipation rates $\epsKK$ and $\epsKKKK$ based on
-both viscosities, and the total kinetic energy dissipation rate $\epsK = \epsKK +
-\epsKKKK$. The product of the maximal wave-vector $\kmax$ with the Kolmogorov scale
-$\eta \equiv (\nu^3 / \epsK)^{1/4}$ is computed to quantify how close our simulations
-are from true Direct Numerical Simulations (DNS). \Add{In practice, it is common to
-consider that simulations are proper DNS with well-resolved small scales when
-$\kmax\eta > 1$ \cite[]{deBruynKops1998,brethouwer_scaling_2007}.} \Add{For the
-statistically stationarity state, the time average of the total energy dissipation rate
-is equal to the injection rate $P_K$ so that the product $\kmax\eta$ depends mostly on
-$\nu$ and $n_h$.} \Add{For most couples $(N,\,\R_i)$, the resolution of the larger
-simulation is fine enough ($\kmax\eta \gtrsim 1$) so that the hyperdiffusion is zero or
-negligible.} \Add{For example, the two simulations analyzed in details in the next
-section (Figures~\ref{fig:buoyancy_fields} to \ref{fig:ratioEwaves_khkz}) are proper
-DNS with $\kmax\eta$ equal to 0.99 and 1.05, respectively (see
-table~\ref{table-better-simuls}).} \Add{There are also few simulations with $0.45 <
-\kmax\eta < 1$ (19 over 78 simulations), which remain slightly under-resolved and
-affected by hyper-viscosity. In that case, small-scales and flow statistics should be
-analyzed carefuly. We checked that these simulations do not change the results
-presented here.}
+To reduce the computational costs, we first simulate the transient state with a coarse
+resolution with $n_h = 320$. For these simulations, two hyperdiffusive terms $-\nu_4
+\nabla^4 \vv$ and $-\kappa_4\nabla^4{b}$ are added to (\ref{eq:Impulsion}) and
+(\ref{eq:Buoyancy}), respectively, in order to keep the dissipative range in the
+simulated scales and avoid thermalisation at small scales. Once a statistically
+steady-state is reached, we increase the resolution of the simulation while decreasing
+hyper-viscosity. Then the simulation is run until reaching a new statistically steady
+state. The previous step is repeated until reaching the highest resolution. We measure
+the turbulent kinetic dissipation rates $\epsKK$ and $\epsKKKK$ based on both
+viscosities, and the total kinetic energy dissipation rate $\epsK = \epsKK + \epsKKKK$.
+The product of the maximal wave-vector $\kmax$ with the Kolmogorov scale $\eta \equiv
+(\nu^3 / \epsK)^{1/4}$ is computed to quantify how close our simulations are from true
+Direct Numerical Simulations (DNS). In practice, it is common to consider that
+simulations are proper DNS with well-resolved small scales when $\kmax\eta > 1$
+\cite[]{deBruynKops1998,brethouwer_scaling_2007}. For the statistically stationarity
+state, the time average of the total energy dissipation rate is equal to the injection
+rate $P_K$ so that the product $\kmax\eta$ depends mostly on $\nu$ and $n_h$. For most
+couples $(N,\,\R_i)$, the resolution of the larger simulation is fine enough
+($\kmax\eta \gtrsim 1$) so that the hyperdiffusion is zero or negligible. For example,
+the two simulations analyzed in details in the next section
+(Figures~\ref{fig:buoyancy_fields} to \ref{fig:ratioEwaves_khkz}) are proper DNS with
+$\kmax\eta$ equal to 0.99 and 1.05, respectively (see table~\ref{table-better-simuls}).
+There are also few simulations with $0.45 < \kmax\eta < 1$ (19 out of 78 simulations),
+which remain slightly under-resolved and affected by hyper-viscosity. In that case,
+small-scales and flow statistics should be analyzed carefuly. We checked that these
+simulations do not change the results presented here.
 
 The Fourier transform of the velocity field $\hatvv = (\hat{v}_x, \hat{v}_y,
 \hat{v}_z)$ can be written using the poloidal-toroidal-shear decomposition (see e.g.
@@ -531,10 +517,10 @@
 $\hatvp$ is the poloidal component, $\hatvt$ the toroidal component, $\vvs$ the shear
 modes component, $\kk=(k_x, k_y, k_z)$ denotes the wave-vector, $k =|\kk| = \sqrt{k_x^2
 + k_y^2 + k_z^2}$ is its modulus, and $k_h = \sqrt{k_x^2 + k_y^2}$ is the modulus of
-the horizontal component of the wave-vector (Figure~\ref{fig:poloidal-toroidal}).
-\Add{Since the toroidal component $\hatvt$ corresponds only to vertical vorticity
-($\hat{\Omega}_z = i k \hatvt \sin \thk$, with $\bOmega = \bnabla \times \vv$ the
-vorticity), we also denote it as the ``vortical" velocity.}
+the horizontal component of the wave-vector (Figure~\ref{fig:poloidal-toroidal}). Since
+the toroidal component $\hatvt$ corresponds only to the vertical vorticity
+($\hat{\Omega}_z = i k \hatvt \sin \thk$, with $\bOmega = \bnabla \times \vv$ being the
+vorticity), we also denote it as the ``vortical" velocity.
 
 \begin{figure}
 \centering
@@ -586,7 +572,7 @@
 \end{equation}
 are the waves modes and the pulsation of the waves, $\thk$ is the angle between $\kk$
 and the stratification axis $\eez$, and $\hata^{(0)} \propto \hatvt$ are the vortical
-modes\Remove{(corresponding to vertical vorticity)}. Equations
+modes. Equations
 (\ref{eq:StratifiedSpectralPoloidalToroidalLinear}-\ref{eq:StratifiedSpectralShearModesLinear})
 show that both shear modes and vortical modes have zero frequency. This means that
 linear waves can only exist in the poloidal velocity and the buoyancy, but not in the
@@ -598,44 +584,42 @@
 forced at large spatial scales $ \left\{\kk ~ | ~ 5 \leq k/\Delta k_h \leq 20 \right\}$
 and small angle $\left\{\kk ~ | ~ |\ok /N - \sin \theta_f| \leq 0.05 \right\}$ where
 $\sin \theta_f = 0.3$, meaning that relatively slow internal waves are forced. The
-forcing scheme is described in Appendix~\ref{appendix:forcing}. \Add{It is neither
-harmonic nor given by a stochastic differential equation. Instead, a time correlated
-forcing is computed via generations of pseudo random numbers and time interpolations.}
-Its correlation time is equal to the period of the forced waves $T_c = 2\pi /(N \sin
-\theta_f)$. \Add{The forcing is normalized such that the kinetic injection rate $P_K$
-is always equal to 1.} \Add{Forcing slow waves is motivated by oceanic applications,
-where waves are generated, among other processes, by slow tides
-\cite{mackinnon_climate_2017, nikurashin_legg_mechanism_2011}. Low frequency forcing is
-also used in order to have a scale separation between forced frequencies and the \bv
-frequency so that one can potentially reproduce features of the oceanic temporal
-spectra close to $N$.}
+forcing scheme is described in Appendix~\ref{appendix:forcing}. It is neither harmonic
+nor given by a stochastic differential equation. Instead, a time correlated forcing is
+computed via generations of pseudo random numbers and time interpolations. Its
+correlation time is equal to the period of the forced waves $T_c = 2\pi /(N \sin
+\theta_f)$. The forcing is normalized such that the kinetic injection rate $P_K$ is
+always equal to 1. Forcing slow waves is motivated by oceanic applications, where waves
+are generated, among other processes, by slow tides \cite{mackinnon_climate_2017,
+nikurashin_legg_mechanism_2011}. Low frequency forcing is also used in order to have a
+scale separation between forced frequencies and the \bv frequency so that one can
+potentially reproduce features of the oceanic temporal spectra close to $N$.
 
 The time advancement is performed using the $4^{th}$ order Runge-Kutta scheme. All
 modes with wave-number modulus larger than $\kmax = 0.8 (n_h/2) \Delta k_h$ are
-truncated to limit aliasing. \Add{We checked that this 0.8 spherical truncation is a
-good compromize consistent with our other numerical choices.} Shear modes and
-vertically invariant vertical velocity (internal waves at $\omega = N$), which are
-absent in flows bounded by walls, are also removed in our simulations \Add{by fixing
-nonlinear transfers to these modes to zero}. \Add{For the simulations without vortical
-modes, the toroidal projection of the velocity and of the nonlinear transfers are set
-to zero. \citet{smith_generation_2002} used a similar procedure in order to distangle
-the roles of potential vorticity modes and waves modes in rotating stratified
-turbulence.}
+truncated to limit aliasing. We checked that this 0.8 spherical truncation is a good
+compromize consistent with our other numerical choices. Shear modes and vertically
+invariant vertical velocity (internal waves at $\omega = N$), which are absent in flows
+bounded by walls, are also removed in our simulations by fixing nonlinear transfers to
+these modes to zero. For the simulations without vortical modes, the toroidal
+projection of the velocity and of the nonlinear transfers are set to zero.
+\citet{smith_generation_2002} used a similar procedure in order to distangle the roles
+of potential vorticity modes and waves modes in rotating stratified turbulence.
 
-\Add{Without dissipation and forcing, equations (\ref{eq:Continuity}-\ref{eq:Buoyancy})
+Without dissipation and forcing, equations (\ref{eq:Continuity}-\ref{eq:Buoyancy})
 conserve the total energy $E = \int ~ \left[ \vv^2 / 2 + b^2 / (2N^2) \right] ~ \diff x
 \diff y \diff z$ and the potential vorticity $\Pi = \bOmega \cdot \left( N^2 \eez +
 \bnabla b \right)$ is a Lagrangian invariant \cite{bartello_geostrophic_1995}. It
 follows that the spatial average of any function of $\Pi$ is conserved. As a special
-case, the potential enstrophy}
+case, the potential enstrophy
 \begin{align}
 V &\equiv \frac{1}{2} \int ~ \Pi^2 ~ \dxdydz \\ &=
 \frac{1}{2} \int ~ N^4 \Omega_z^2 ~ \dxdydz + \int ~ N^2 \Omega_z
 \bOmega \cdot \bnabla b ~ \dxdydz + \frac{1}{2} \int ~ \left( \bOmega
 \cdot \bnabla b \right)^2 ~ \dxdydz \\ &\equiv V_2 + V_3 + V_4
 \end{align}
-\Add{is an invariant of equations (\ref{eq:Continuity}-\ref{eq:Buoyancy}) if no
-dissipation and no forcing. For a flow without vertical vorticity $V_2 = V_3 =0$.}
+is an invariant of equations (\ref{eq:Continuity}-\ref{eq:Buoyancy}) in absence of
+dissipation and forcing. For a flow without vertical vorticity $V_2 = V_3 =0$.
 
 \section{Results}
 \label{sec:results}
@@ -664,20 +648,17 @@
 anisotropic. When vortical modes are removed, we rather use
 \begin{equation}
 \label{eq:IcoeffsProj}
-	\Ivelo  = \frac{2 E_{\rm kin,z}}{E_{\rm kin}} \Remove{, ~~~~ \text{and}
-~~~~ \Idiss = \frac{1 - \varepsilon_{\rm kin,z}/\varepsilon_{\rm kin}}{(1 - 1/2)}}
+	\Ivelo  = \frac{2 E_{\rm kin,z}}{E_{\rm kin}}
 \end{equation}
-because it is expected that only half of the kinetic energy is contained in
-\Remove{(respectively dissipated by)} the vertical velocity field \Remove{(respectively
-the vertical gradients) if the energy spectra are independent of the angles in this
-case}. Figure~\ref{fig:regimes} shows the variations of $\Ivelo$ and $\Idiss$ with
-$(F_h,\R)$. Small points correspond to strong anisotropic energy dissipation, while
-large points correspond to isotropic energy dissipation. Dark points correspond to
-strong large-scale isotropy, while light points correspond to isotropic flows. The
-combination of $\Ivelo$ and $\Idiss$ allows to distinguish between four regimes
-(Figure~\ref{fig:regimes}$\rm (a)$): the passive scalar regime where $\Ivelo \simeq 1$,
-the weakly stratified regime where $0.5 \lesssim \Ivelo \lesssim 1$, and the strongly
-stratified regimes where $\Ivelo \leq 0.5$. As explained in
+because it is expected that only half of the kinetic energy is contained in the
+vertical velocity field. Figure~\ref{fig:regimes} shows the variations of $\Ivelo$ and
+$\Idiss$ with $(F_h,\R)$. Small points correspond to strong anisotropic energy
+dissipation, while large points correspond to isotropic energy dissipation. Dark points
+correspond to strong large-scale isotropy, while light points correspond to isotropic
+flows. The combination of $\Ivelo$ and $\Idiss$ allows to distinguish between four
+regimes (Figure~\ref{fig:regimes}$\rm (a)$): the passive scalar regime where $\Ivelo
+\simeq 1$, the weakly stratified regime where $0.5 \lesssim \Ivelo \lesssim 1$, and the
+strongly stratified regimes where $\Ivelo \leq 0.5$. As explained in
 \cite{brethouwer_scaling_2007}, the strongly stratified flows fall in two regimes: The
 Layered Anisotropic Stratified Turbulence (LAST) regime where the dissipation is
 isotropic because a 3D turbulence range can develop, typically when $\R \geq 10$; The
@@ -726,16 +707,16 @@
 stratification by removing vortical modes. In the two next subsections, we perform a
 spatiotemporal analysis of a couple of stongly stratified turbulent simulations with
 $(N,\,\R_i\equiv P_K/\nu N^2)=(40,20)$ and aspect ratio $L_z/L_h=1/4$.
-\Add{Figure~\ref{fig:buoyancy_fields} shows the buoyancy fields for these two
-simulations with the same color scale. We observe that the flow is layered in the
-vertical direction and that overturning (horizontal vorticity) is present with or
-without vortical modes. It is a standard feature of strongly stratified turbulence
+Figure~\ref{fig:buoyancy_fields} shows the buoyancy fields for these two simulations
+with the same color scale. We observe that the flow is layered in the vertical
+direction and that overturning (horizontal vorticity) is present with or without
+vortical modes. It is a standard feature of strongly stratified turbulence
 \cite{laval_forced_2003, lindborg_energy_2006, brethouwer_scaling_2007,
 waite_stratified_2011}. With vortical modes, the vertical vorticity is not zero
 (Figure~\ref{fig:buoyancy_fields}(a)) so the buoyancy has a different structure than
 when vortical modes are absent (Figure~\ref{fig:buoyancy_fields}(b)). Without vortical
 modes the dynamics in the horizontal direction is irrotational and the buoyancy field
-has a larger amplitude.}
+has a larger amplitude.
 
 \begin{figure}
 \includegraphics[width=1.0\textwidth]{Figure3}
@@ -748,8 +729,8 @@
 
 \begin{figure}
 \includegraphics[width=1.0\textwidth]{figure4}
-\caption{\Add{Snapshots of the buoyancy fields for simulations $(N,\R_i)=(40,20)$ with
-$\rm (a)$ and without $\rm (b)$ vortical modes.} \label{fig:buoyancy_fields}}
+\caption{Snapshots of the buoyancy fields for simulations $(N,\R_i)=(40,20)$ with $\rm
+(a)$ and without $\rm (b)$ vortical modes. \label{fig:buoyancy_fields}}
 \end{figure}
 
 \subsection{Energy budget in the strongly stratified regime}
@@ -797,7 +778,7 @@
 are respectively the kinetic energy injection rate, the kinetic energy transfer, the
 potential energy transfer, the conversion of kinetic energy to potential energy, the
 kinetic energy dissipation, and the potential energy dissipation. In the last
-equations, $(\cdot)^*$ denotes the complex conjugate, \Add{$\Re(\cdot)$ the real part},
+equations, $(\cdot)^*$ denotes the complex conjugate, $\Re(\cdot)$ the real part,
 $\left\langle\cdot \right\rangle$ stands for the averaging operator, and
 $\bar{\bar{P}}_{\kk} = \mathbb{I} - \eek \otimes \eek$ is the projector onto the plane
 orthogonal to $\kk$.  In this section, we study azimuthal average of the energy budget.
@@ -968,11 +949,11 @@
 
 
 Figure~\ref{fig:spectra_khkz} shows the $(k_h,k_z)$ spectra. Obviously, the simulation
-without vortical modes has no energy in the toroidal velocity. \Add{When vortical modes
-are present, an important part of the energy is contained in one vortical mode with
-$(k_h,k_z) = (\Delta k_h, 2 \Delta k_z)$, corresponding to large, nearly vertical
+without vortical modes has no energy in the toroidal velocity. When vortical modes are
+present, an important part of the energy is contained in one vortical mode with
+$(k_h,k_z) = (\Delta k_h, 2 \Delta k_z)$, corresponding to large, nearly vertically
 stacked shear layers (Figure~\ref{fig:spectra_khkz}$\rm (a)$). Energy then tends to be
-accumulated at the smallest horizontal wave vectors, close to shear modes.} When
+accumulated at the smallest horizontal wave vectors, close to shear modes. When
 vortical modes are removed, energy is sill concentrated in the same wave-vectors, but
 in the form of poloidal an potential energy (Figure~\ref{fig:spectra_khkz}$\rm
 (d)$-$\rm (f)$). Except for this qualitative difference, the toroidal, poloidal, and
@@ -1047,9 +1028,8 @@
 \caption{Normalized conversion to potential energy $\tilde{\mathcal{B}}$
 (\ref{eq:transfers}) for the simulations at $(N,\R_i)=(40,20)$, with $\rm (a)$ and
 without $\rm (b)$ vortical modes. The magenta dotted line corresponds to
-$\gamma_{\kk}=1$, \Add{and} the green dotted line to the dissipative scale \Remove{,
-and the continuous black line to $k_z=k_h$}. The orange box corresponds to the forcing
-region. \label{fig:conversion}}
+$\gamma_{\kk}=1$, and the green dotted line to the dissipative scale. The orange box
+corresponds to the forcing region. \label{fig:conversion}}
 \end{figure}
 
 \subsection{Spatiotemporal analysis in the buoyancy range}
@@ -1331,14 +1311,13 @@
 We observed that removing vortical modes helps to have a better overall balance between
 poloidal kinetic energy and potential energy. However, the spatial spectra appear to
 behave very similarly with or without vortical modes in our simulations. The spectral
-energy budget reveals that \Add{there is no range in $(k_h,k_z)$ for which the}
-conversion between kinetic energy and potential energy \Remove{do not show
-fluctuations} \Add{fluctuates} around zero, as we would expect for a system of
-statistically stationary waves. Additionally, the conversion between potential energy
-and kinetic energy becomes small where the wave dissipation parameter $\gamma_{\kk}$
-(\ref{eq:DissipationParameter}) is larger than unity. Due to the anisotropy of
-stratified flows, this means that it exists a range of scales where waves are
-efficiently dissipated by viscosity, but not necessarily the vortices.
+energy budget reveals that there is no range in $(k_h,k_z)$ for which the conversion
+between kinetic energy and potential energy fluctuates around zero, as we would expect
+for a system of statistically stationary waves. Additionally, the conversion between
+potential energy and kinetic energy becomes small where the wave dissipation parameter
+$\gamma_{\kk}$ (\ref{eq:DissipationParameter}) is larger than unity. Due to the
+anisotropy of stratified flows, this means that it exists a range of scales where waves
+are efficiently dissipated by viscosity, but not necessarily the vortices.
 
 A spatiotemporal analysis in the buoyancy range showed that increasing stratification
 is necessary to decrease the nonlinear broadening to the level required for a weak wave
@@ -1408,7 +1387,7 @@
 dimensionless factors, which are set to one for simplicity. On
 Figure~\ref{fig:discussion}, we plot our simulations in the $(F_h,\R)$ plane, as well
 as the simulations of
-\Add{\cite{brethouwer_scaling_2007,waite_potential_2013,reun_parametric_2018,lam_energy_2021}}
+\cite{brethouwer_scaling_2007,waite_potential_2013,reun_parametric_2018,lam_energy_2021}
 and some recent experiments \cite{rodda_experimental_2022}. For the simulations of
 \cite{reun_parametric_2018}, the Froude number is computed using our definition
 (\ref{eq:FhR}), and the buoyancy Reynolds number as $\R = Re F_h^2$ in order to compare
@@ -1418,13 +1397,13 @@
 simulations do not lie in the region corresponding to (\ref{eq:DissInBuoyancy}),
 (\ref{eq:LargeRe}), and (\ref{eq:SmallChi}). On the contrary,
 \cite{reun_parametric_2018} attained such a region in their simulations, which provides
-a tangible explanation for why they obtained better signatures of WWT. \Add{It is worth
+a tangible explanation for why they obtained better signatures of WWT. It is worth
 mentioning that \cite{reun_parametric_2018} did not have to remove shear nor vortical
 modes in their simulations to observe signatures of internal wave turbulence. It
 indicates that there must be some threshold below which both shear and vortical modes
 do not grow if not directly forced. Such a situation would be analogous to rotating
 flows where there is a threshold below which geostrophic modes do not grow
-\cite{reun_experimental_2019}.} \Add{We observe that \cite{waite_potential_2013} also
+\cite{reun_experimental_2019}. We observe that \cite{waite_potential_2013} also
 attained very small $\R$. In this study, the authors forced vortical modes so their
 simulations are not well suited for WWT. Yet, they showed that potential enstrophy
 tends to be quadradic (i.e. $V \simeq V_2$) for $F_h, \R \ll 1$ and that $V_2 \propto
@@ -1434,12 +1413,12 @@
 keeping $F_h$ or $\R$ constant in \cite{lam_energy_2021}. To our knowledge, additional
 studies are needed to confirm if the wave energy ratio is for $\R\ll 1$ and $F_h \ll 1$
 a decreasing function of $\R$, and if a threshold below which both shear and vortical
-modes are stable exists.}
+modes are stable exists.
 
 \begin{figure}
 \includegraphics[width=0.66\textwidth]{figure17.png}
 \caption{Simulations of the present study and
-\Add{\cite{brethouwer_scaling_2007,waite_potential_2013,reun_parametric_2018,lam_energy_2021}},
+\cite{brethouwer_scaling_2007,waite_potential_2013,reun_parametric_2018,lam_energy_2021},
 and the experiments \cite{rodda_experimental_2022} in the $(F_h, \R)$ parameters space.
 The colored full lines corresponding to conditions (\ref{eq:DissInBuoyancy}),
 (\ref{eq:LargeRe}), and (\ref{eq:SmallChi}) (see legend). The blue dashed line
@@ -1458,11 +1437,11 @@
 
 
 \begin{acknowledgments}
-\Add{We thank two anonymous reviewers for their constructive feedback.} This project
-was supported by the Simons Foundation through the Simons collaboration on wave
-turbulence. Part of the computations have been done on the ``Mesocentre SIGAMM''
-machine, hosted by Observatoire de la Cote d'Azur. The authors are grateful to the OPAL
-infrastructure from Université Côte d'Azur and the Université Côte d’Azur's Center for
+We thank two anonymous reviewers for their constructive feedback. This project was
+supported by the Simons Foundation through the Simons collaboration on wave turbulence.
+Part of the computations have been done on the ``Mesocentre SIGAMM'' machine, hosted by
+Observatoire de la Cote d'Azur. The authors are grateful to the OPAL infrastructure
+from Université Côte d'Azur and the Université Côte d’Azur's Center for
 High-Performance Computing for providing resources and support. This work was granted
 access to the HPC/AI resources of IDRIS under the allocation 2022-A0122A13417 made by
 GENCI.
@@ -1473,16 +1452,16 @@
 \section{Forcing scheme}
 \label{appendix:forcing}
 
-\Add{Our forcing is designed to excite waves and computed via generation of pseudo
-random numbers with uniform distribution and time interpolation. Namely, we use the
-following algorithm:}
+Our forcing is designed to excite waves and computed via generation of pseudo random
+numbers with uniform distribution and time interpolation. Namely, we use the following
+algorithm:
 
 \begin{algorithm}[H]
 $t = 0$; $t_0 = 0$; \\
 Generate two complex random fields $\hat f_0(\kk)$ and $\hat f_1(\kk)$; \\
 $\hatff = \hat f_0 \, \eep$;\\
 Normalize $\hatff$ to ensure
-\Add{$P_K(t) = \sum\limits_{\kk} ~ \Re \left[ \hatff \cdot \hatvv^* + \frac{\Delta t}{2} |\hatff|^2 \right] =  1$};\\
+$P_K(t) = \sum\limits_{\kk} ~ \Re \left[ \hatff \cdot \hatvv^* + \frac{\Delta t}{2} |\hatff|^2 \right] =  1$;\\
 \While{$t \leq T$}{
 	$t = t + \Delta t$; \\
 	\If{$t - t_0 \geq T_c$}{
@@ -1495,9 +1474,9 @@
 }
 %\caption{}
 \end{algorithm}
-where $\Delta t$ is the time increment at each time step \Add{and $T$ is the final time
-of the simulation}. \Add{The random complex fields are built such that their inverse
-Fourier transform is real and they are null for unforced wavenumbers.}
+where $\Delta t$ is the time increment at each time step and $T$ is the final time of
+the simulation. The random complex fields are built such that their inverse Fourier
+transform is real and they are null for unforced wavenumbers.
 
 \section{List of simulations}
 
diff --git a/2022strat_polo_proj/input/main_changes.tex b/2022strat_polo_proj/input/main_changes.tex
new file mode 100755
--- /dev/null
+++ b/2022strat_polo_proj/input/main_changes.tex
@@ -0,0 +1,1520 @@
+\documentclass[
+aps,
+prb,
+superscriptaddress,
+reprint,
+onecolumn,
+amsfonts,
+amssymb,
+amsmath,
+]{revtex4-2}
+
+% The result seems nicer with revtex4-2, but
+% Debian 9 (Stretch) does not have revtex4-2
+
+\usepackage[utf8]{inputenc}
+\usepackage{hyperref} % Required for customising links and the PDF*
+\hypersetup{pdfpagemode={UseOutlines},
+	bookmarksopen=true,
+	bookmarksopenlevel=0,
+	hypertexnames=false,
+	colorlinks=true, % Set to false to disable coloring links
+	citecolor=blue, % The color of citations
+	linkcolor=red, % The color of references to document elements (sections, figures, etc)
+	urlcolor=black, % The color of hyperlinks (URLs)
+	pdfstartview={FitV},
+	unicode,
+	breaklinks=true,
+}
+\usepackage{graphicx}
+\usepackage{grffile}
+\usepackage{color}
+\usepackage{array}
+\usepackage{hhline}
+\usepackage[]{algorithm2e}
+\usepackage{booktabs}
+\usepackage{ulem}
+\usepackage{siunitx}
+\sisetup{
+	inter-unit-product = \ensuremath{{}\!\cdot\!{}},
+	detect-all,
+	separate-uncertainty = true,
+	exponent-product = \times,
+	space-before-unit = true,
+	output-decimal-marker = {,},
+	multi-part-units = brackets,
+	range-phrase = --,
+	% allow-number-unit-breaks,
+	list-final-separator = { et },
+	list-pair-separator = { et },
+	abbreviations
+}
+
+\usepackage{float}
+
+\linespread{1.05}
+
+
+\setlength{\tabcolsep}{7pt}
+
+% could be commented when the paper is accepted
+\graphicspath{{../tmp/}}
+
+\newcommand{\cor}[1]{\textcolor{red}{#1}}
+\newcommand{\todo}[1]{\textcolor{red}{TODO: #1}}
+
+\newlength{\figwidth}
+\setlength{\figwidth}{120mm}
+% \setlength{\figwidth}{0.7\textwidth}  % useful in single column
+
+
+\newcommand{\R}{\mathcal{R}}
+
+\newcommand{\eps}{\varepsilon}
+\newcommand{\epsK}{{\varepsilon_{\!\scriptscriptstyle K}}}
+\newcommand{\epsKK}{{\varepsilon_{\!\scriptscriptstyle K 2}}}
+\newcommand{\epsKKKK}{{\varepsilon_{\!\scriptscriptstyle K 4}}}
+\newcommand{\epsA}{{\varepsilon_{\!\scriptscriptstyle P}}}
+
+
+\newcommand{\xx}{\boldsymbol{x}}
+\newcommand{\rr}{\boldsymbol{r}}
+\newcommand{\kk}{\boldsymbol{k}}
+\newcommand{\eek}{\boldsymbol{e}_{\boldsymbol{k}}}
+\newcommand{\eeh}{\boldsymbol{e}_h}
+\newcommand{\eep}{\boldsymbol{e}_{p\kk}}
+\newcommand{\eet}{\boldsymbol{e}_{t\kk}}
+\newcommand{\eex}{\boldsymbol{e}_x}
+\newcommand{\eey}{\boldsymbol{e}_y}
+\newcommand{\eez}{\boldsymbol{e}_z}
+\newcommand{\cc}{\boldsymbol{c}}
+\newcommand{\uu}{\boldsymbol{u}}
+\newcommand{\vv}{\boldsymbol{v}}
+\newcommand{\hatvv}{\hat{\boldsymbol{v}}}
+\newcommand{\hatvp}{\hat{v}_{p}}
+\newcommand{\hatvt}{\hat{v}_{t}}
+\newcommand{\hatb}{\hat{b}}
+\newcommand{\hata}{\hat{a}}
+\newcommand{\vvs}{\hat{\boldsymbol{v}}_{s}}
+\newcommand{\ff}{\boldsymbol{f}}
+\newcommand{\hatf}{\hat{f}}
+\newcommand{\hatff}{\hat{\boldsymbol{f}}}
+\newcommand{\bomega}{\boldsymbol{\omega}}
+\newcommand{\bnabla}{\boldsymbol{\nabla}}
+\newcommand{\Dt}{\mbox{D}_t}
+\newcommand{\p}{\partial}
+\newcommand{\mean}[1]{\langle #1 \rangle}
+\newcommand{\epsP}{\varepsilon_{\!\scriptscriptstyle P}}
+\newcommand{\epsm}{\varepsilon_{\!\scriptscriptstyle m}}
+\newcommand{\CKA}{C_{K\rightarrow A}}
+\newcommand{\D}{\mbox{D}}
+\newcommand{\diff}{\text{d}}
+\newcommand{\dxdydz}{\diff x \diff y \diff z}
+\newcommand{\bv}{Brunt-V\"ais\"al\"a }
+\newcommand{\kmax}{k_{\max}}
+\newcommand{\thk}{\theta_{\kk}}
+\newcommand{\phk}{\varphi_{\kk}}
+\newcommand{\thf}{\theta_f}
+\newcommand{\ok}{\omega_{\kk}}
+
+\newcommand{\bOmega}{\boldsymbol{\Omega}}
+
+\newcommand{\oemp}{\omega_{\text{emp}}}
+\newcommand{\odoppler}{\delta \omega_{\text{doppler}}}
+
+\newcommand{\Ivelo}{I_{\rm kin}}
+\newcommand{\Idiss}{I_{\text{diss}}}
+
+\newcommand{\Pikin}{\Pi_{\rm kin}}
+\newcommand{\Pipot}{\Pi_{\text{pot}}}
+
+\newcommand{\Etoro}{E_{\rm toro} }
+\newcommand{\Epolo}{E_{\rm polo} }
+\newcommand{\Epot}{E_{\rm pot} }
+\newcommand{\Eequi}{E_{\rm equi} }
+\newcommand{\Ediff}{E_{\text{diff}}}
+\newcommand{\Ewave}{E_{\rm wave}}
+\newcommand{\Ewaver}{\tilde{E}_{\rm wave}}
+
+\newcommand{\kb}{k_{\rm b}}
+\newcommand{\ko}{k_{\rm O}}
+\newcommand{\kd}{k_{\rm d}}
+
+\setlength\parindent{0pt}
+
+
+% For review
+\newcommand{\Add}[1]{{\color{blue}#1}}
+\newcommand{\Remove}[1]{{\color{red}\sout{#1}}}
+\newcommand{\Comment}[1]{{\color{green}#1}}
+
+
+\begin{document}
+
+\title{Internal gravity waves in stratified flows with and without vortical modes}
+
+\author{Vincent Labarre}
+\email[]{vincent.labarre@oca.eu}
+\affiliation{Universit\'{e} C\^{o}te d'Azur, Observatoire de la C\^{o}te d'Azur, CNRS,
+Laboratoire Lagrange, Nice, France}
+
+\author{Pierre Augier}
+\email[]{pierre.augier@univ-grenoble-alpes.fr}
+\affiliation{Laboratoire des Ecoulements G\'eophysiques et Industriels, Universit\'e
+Grenoble Alpes, CNRS, Grenoble-INP, F-38000 Grenoble, France}
+
+\author{Giorgio Krstulovic}
+\email[]{giorgio.krstulovic@oca.eu}
+\affiliation{Universit\'{e} C\^{o}te d'Azur, Observatoire de la C\^{o}te d'Azur, CNRS,
+Laboratoire Lagrange, Nice, France}
+
+\author{Sergey Nazarenko}
+\email[]{sergey.nazarenko@unice.fr}
+\affiliation{Universit\'{e} C\^{o}te d'Azur, CNRS, Institut de Physique de Nice -
+INPHYNI, Nice, France}
+
+\begin{abstract}
+The comprehension of stratified flows is important for geophysical and astrophysical
+applications. The Weak Wave Turbulence theory aims to provide a statistical description
+of internal gravity waves propagating in the bulk of such flows. However, internal
+gravity waves are usually perturbed by other structures present in stratified flow,
+namely the shear modes and the vortical modes. In order to check whether a weak
+internal gravity wave turbulence regime can occur, we perform direct numerical
+simulations of stratified turbulence without shear modes, and with or without vortical
+modes at various Froude and buoyancy Reynolds numbers. We observe that removing
+vortical modes naturally helps to have a better overall balance between poloidal
+kinetic energy, involved in internal gravity waves, and potential energy. However,
+conversion between kinetic energy and potential energy does not necessarily show
+fluctuations around zero in our simulations, as we would expect for a system of weak
+waves. A spatiotemporal analysis reveals that removing vortical modes helps to
+concentrate the energy around the wave frequency, but it is not enough to observe a
+weak wave turbulence regime. Yet, we observe that internal gravity waves whose
+frequency are large compared to the eddy turnover time are present, and we also find
+evidences for slow internal gravity waves interacting by Triadic Resonance
+Instabilities in our strongly stratified flows simulations. Finally, we propose
+conditions that should be fulfilled in order to observe a weak internal gravity waves
+turbulence regime in real flows.
+\end{abstract}
+
+%----------------------------------------------------------------------------------------
+
+% Print the title
+\maketitle
+
+%----------------------------------------------------------------------------------------
+%	ARTICLE CONTENTS
+%----------------------------------------------------------------------------------------
+
+\section{Introduction}
+\label{sec:introduction}
+
+% WWT
+As eddies in classical hydrodynamic turbulence, waves in nonlinear systems interact and
+transfer conserved quantities along scales in a cascade process. The Weak-Wave
+Turbulence (WWT) theory aims to provide a statistical description of the system when
+the nonlinearity is small \cite{zakharov_kolmogorov_1992,nazarenko_wave_2011,
+nazarenko_wave_2015}. The applications of this theory encompass capillary-gravity waves
+\cite{falcon_experiments_2022}, gravito-inertial waves in rotating and stratified
+fluids \cite{caillol_kinetic_2000, galtier_weak_2003, medvedev_turbulence_2007}, 2D
+acoustic waves \cite{griffin_energy_2022}, elastic plates \cite{during_weak_2006},
+Alfvén waves in magnetohydrodynamics (MHD) \cite{galtier_weak_2000}, Kelvin waves in
+superfluids \cite{lvov_weak_2010}, density waves in Bose-Einstein condensates
+\cite{dyachenko_optical_1992}, and gravitational waves \cite{galtier_turbulence_2017}.
+
+% Weak nonlinearity and small nonlinear broadening
+Three key hypotheses are used in weak-wave turbulence (WWT) theory. The first one is
+weak nonlinearity of the dynamical equations. This condition is often translated in
+terms of spatial scale separation: the considered scale on one side and the saturation
+scale on the other. The saturation scale is here defined as the scale at which the weak
+nonlinearity no longer holds, inducing wave breaking. This hypothesis is quite similar
+to the separation between the integral scale and the dissipative scale in the classical
+picture of 3D hydrodynamic turbulence, which is necessary for observing an inertial
+range. The second hypothesis corresponds to a separation of time scales. It requires
+that the linear time, given by the wave period $\tau_{\rm L}$, is much shorter than the
+nonlinear time of interactions between waves  $\tau_{\rm NL}$
+\cite{nazarenko_wave_2011}. The third hypothesis is that the nonlinear broadening
+$\delta \omega \sim 1/\tau_{\rm NL}$, which measures the frequency of  nonlinear
+interactions, must be larger than the frequency gap between wave modes $\Delta \omega$
+in the discrete Fourier space. This last condition is necessary to permit enough
+nonlinear interactions among waves such that we can consider Fourier space as being
+continuous \cite{lvov_discrete_2010}.
+
+For many physical situations, at least one of the three hypotheses of WWT is broken in
+some range of scales, which reduces the validity of the theory
+\cite{biven_breakdown_2001, lvov_discrete_2010}. Yet, when scale separation is observed
+in space and time, a WWT range can emerge in wave energy spectra. Besides the practical
+difficulties for obtaining scale separation, testing WWT in isotropic systems is
+conceptually simpler. For this reason, main progress have been made in the experimental
+and numerical verification of WWT for elastic plates \cite{miquel_nonstationary_2011,
+yokoyama_identification_2014}, capillary-gravity waves \cite{pan_direct_2014,
+falcon_experiments_2022}, density waves in Bose-Einstein condensates
+\cite{zhu_testing_2022}, and 2D acoustic waves \cite{griffin_energy_2022}.
+
+Anisotropic turbulence is, generally speaking, more difficult to investigate than
+isotropic turbulence because of the reduced number of symmetries of the considered
+system. For example, it has been shown that studying 2D spatial spectra instead of 1D
+integrated spectra is essential to investigate stratified turbulence
+\cite{yokoyama_energy-based_2019}. Anisotropy makes the problem multidimensional in
+Fourier space, which makes the notion of spectral energy fluxes, and other relevant
+quantities, more difficult to define and calculate than in isotropic turbulence
+\cite{yokoyama_energy-flux_2021}. An additional difficulty is that scale separations
+required by WWT can also be anisotropic. For linearly stratified flows, the weak
+nonlinearity of the Navier-Stokes equations requires that the spatial scale separation
+\begin{equation}
+k/\kb \ll 1
+\end{equation}
+is satisfied,  where $\kb = N/U_h$ is the buoyancy wave-vector, $k$ is the wave-vector
+modulus, $U_h$ is the rms of the horizontal velocity, and $N$ is the \bv frequency.
+This condition should be fulfilled in order to avoid wave-breaking
+\cite{waite_stratified_2011}. The temporal scale separation between linear waves and
+eddies leads to the different condition
+\begin{equation}
+\tau_{\rm L}/\tau_{\rm NL} = \frac{(\epsK k^2)^{1/3}}{N k_h/k} \ll 1,
+\end{equation}
+where $\epsK$ is the kinetic energy dissipation rate, $k_h$ is the horizontal
+wave-vector modulus, $\tau_{\rm L} = 2\pi k /(N k_h)$ is the period of internal gravity
+waves, and $\tau_{\rm NL} = 2\pi/(k^2 \epsK)^{1/3}$ is the eddy turnover time
+\cite{yokoyama_energy-based_2019}. Physically, when $\tau_{\rm L}/\tau_{\rm NL} \ll 1$,
+the waves are faster than the typical time of their nonlinear interactions. Due to the
+anisotropic dispersion relation of internal gravity waves, time scale separation is
+less valid for small $k_h$, and even impossible for modes with $k_h=0$. Consequently,
+separation of times scales can be violated even if a large separation of spatial scales
+is observed (i.e. $k/\kb \ll 1$). A similar situation takes place in plasmas under the
+effect of a strong magnetic field, or strongly rotating and stratified flows. It is
+therefore more difficult to observe signatures of WWT in anisotropic systems.
+
+Stratified turbulence is not only an interesting conceptual problem, but it is also
+essential for understanding geophysical flows. \cite{staquet_internal_2002,
+vallis_atmospheric_2017}. In particular, sub-grid parameterizations in climate models
+require an understanding of the role of waves in mixing \cite{mackinnon_climate_2017,
+gregg_mixing_2018}. It is therefore not surprising that stratified turbulence received
+a particular attention from both the ``strong'' turbulence community
+\cite{billant_self-similarity_2001, waite_stratified_2004, waite_stratified_2006,
+lindborg_energy_2006, brethouwer_scaling_2007, waite_stratified_2011,
+kimura_energy_2012, bartello_sensitivity_2013, brunner-suzuki_upscale_2014,
+augier_stratified_2015, maffioli_vertical_2017}, and the WWT community
+\cite{caillol_kinetic_2000, lvov_hamiltonian_2001, lvov_weak_2010,
+dematteis_downscale_2021, dematteis_origins_2022}. It turns out that internal waves are
+effectively important in the dynamics and mixing in stratified flows
+\cite{maffioli_signature_2020, lam_partitioning_2020, lam_energy_2021}, and that
+three-wave resonant interactions are responsible for slow energy transfers between
+different scales \cite{brouzet_internal_2016, davis_succession_2020,
+rodda_experimental_2022}. Internal waves can be excited, for example, by tides
+\cite{reun_parametric_2018} or by linearly unstable wave attractors
+\cite{brouzet_internal_2016}. However, many questions and issues remain about the
+applicability of WWT to stratified flows. In particular, it is not yet clear under
+which conditions a weak wave turbulence regime could occur.
+
+It has been showed that a wave dominated region should lie in the spectral region
+defined by $\tau_{\rm L} / \tau_{\rm NL} < 1/3$ \cite{yokoyama_energy-based_2019}, in
+agreement with observations made for MHD \cite{meyrand_direct_2016}. However, the
+observation of a system of weakly interacting internal gravity waves is still delicate
+due to non-wave structures such as shear modes (purely vertical shear) and vortical
+modes (vertical vorticity). The same problem arises in rotating flows, in which the
+geostrophic modes play the role of the non-propagative structures. It has been
+ingeniously shown that tidal forcing is very efficient at triggering weakly nonlinear
+internal gravity waves \cite{reun_parametric_2018}, in a way reminiscent of the
+inertial waves' excitation by libration in rotating fluids \cite{le_reun_inertial_2017,
+reun_experimental_2019}. Experimentalists overcame this difficulty by using honeycomb
+grids at the top and bottom boundaries of a rotating tank to dissipate geostrophic
+modes efficiently, allowing them to observe weak inertial wave turbulence
+\cite{brunet_shortcut_2020, monsalve_quantitative_2020}, in a similar fashion to the
+numerical study of \cite{le_reun_inertial_2017}. More recently, it has been observed
+that resonance with the container modes also prevent to observe weak internal gravity
+wave turbulence, and that the introduction of slightly tilted panels at the top and at
+the bottom of the fluid domain allows to inhibit the emergence of these modes
+\cite{lanchon_internal_2023}. All these works point out that some wave systems can
+generate non-wave motions that can severely affect the wave dynamics and should be
+suppressed in experiments aiming to observe wave turbulence. In stratified flows these
+modes correspond to vertical shear, vertical vorticity, an eventually the container
+modes.
+
+The most commonly used prediction of WWT is the scale-invariant stationary solution to
+the kinetic equation. This solution gives the expected spatial energy spectra of waves
+in the statistically steady state. There are two types of solutions: the thermodynamic
+equilibrium solution (Rayleigh-Jeans spectra), and the non-equilibrium solutions that
+are linked to the cascade of the dynamical invariants of the system along scales
+(Kolmogorov-Zakharov spectra) \cite{nazarenko_wave_2011}. It is important to note that
+the Kolmogorov-Zakharov spectra are obtained after the Zakharov transformation, which
+can add spurious solutions. Therefore, these spectra can be considered as valid only if
+the collision integral in the original wave kinetic equation converges. The importance
+of the nonlinear interactions among internal gravity waves was recognized early,
+leading to several derivations of waves kinetic equations (see
+\cite{muller_nonlinear_1986, lvov_resonant_2012}). Derivation of the
+Kolmogorov-Zakharov spectra in the limit of large vertical wave-numbers $k_z \simeq k
+\gg k_h$ can be found in \cite{caillol_kinetic_2000, lvov_hamiltonian_2001}. This
+candidate solution corresponds to an energy cascade, and is given by the energy spectra
+\begin{equation}
+E(k_h, k_z) \sim k_h^{-3/2} k_z^{-3/2} ~~~~ \text{and} ~~~~ E(\omega, k_z) \sim \omega^{-3/2} k_z^{-2},
+\end{equation}
+\Add{where the change of coordinates is defined by the dispersion relation $E(\omega,
+k_z) = E(k_h,k_z) \left( \partial \ok / \partial k_h \right)^{-1}$.} Yet, it was noted
+that this candidate does not satisfy the requirement of locality, i.e. the collision
+integral diverges on it. In the other words, it is not a valid mathematical solution of
+the kinetic equation. Later, it was shown that power law solutions $E(k_h,k_z) \sim
+k_h^{-\alpha_h} ~ k_z^{-\alpha_z}$ have convergent collision integral's contributions
+only on the segment $\alpha_h \in ]2, 3[, \alpha_z=1$ \cite{lvov_oceanic_2010,
+dematteis_downscale_2021}. The collision integral was then computed numerically on this
+segment, and it was deduced that the only scale invariant stationary solution to the
+kinetic equation was close to
+\begin{equation}
+\label{eq:WWTpredictions}
+E(k_h, k_z) \sim k_h^{-1.69} k_z^{-1} ~~~~ \text{and} ~~~~ E(\omega, k_z) \sim \omega^{-1.69} k_z^{-1.69}.
+\end{equation}
+It was also shown that the dominant contributions to the collision integral corresponds
+to non-local transfers first identified by McComas \cite{mccomas_resonant_1977},
+notably the Parametric Subharmonic Instability (PSI) \cite{mccomas_resonant_1977,
+muller_nonlinear_1986} observed in oceans \cite{mackinnon_parametric_2013} and
+experiments \cite{rodda_experimental_2022}, consisting in the resonant interaction of a
+primary wave and two smaller-scale waves of nearly half the frequency.
+
+
+When $k/\kb \gg 1$, nonlinearity is not small and WWT is not valid. Based on the idea
+that there exists a range of spatial scales where buoyancy force has the same order of
+magnitude than inertia, the following 1D integrated energy spectra
+\begin{equation}
+\label{eq:CriticalBalance}
+E_{\rm 1D}(k_h) \sim \epsK^{2/3} k_h^{-5/3} ~~~~ \text{and} ~~~~ E_{\rm 1D}(k_z) \sim N^2 k_z^{-3}
+\end{equation}
+where predicted \cite{lindborg_energy_2006}. In the strong wave turbulence context,
+these predictions can also be obtained using critical balance arguments
+\cite{nazarenko_critical_2011, nazarenko_wave_2011}. At even smaller scales, where
+stratification is negligible, an isotropic range with the energy spectra
+\begin{equation}
+\label{eq:Isotropic}
+E_{\rm 1D}(k) \sim \epsK^{2/3} k^{-5/3}
+\end{equation}
+is expected. It happens when $\kb < \ko \ll k \ll k_{\eta}$ where $\ko$ is the Ozmidov
+wave-vector and $k_{\eta}$ is the Kolmogorov wave-vector.
+
+
+The present study first deals with the existence and properties of a weak internal
+gravity wave turbulence regime. To this end, we perform numerical simulations of
+stratified turbulence at various \bv frequency and viscosity. In the same spirit as
+simulations and experiments of rotating flows \cite{le_reun_inertial_2017,
+brunet_shortcut_2020, monsalve_quantitative_2020}, we remove shear modes in all of our
+simulations. For each values of the control parameters, we perform two ``twin''
+simulations: one where the vortical modes remain, and one where vortical modes are
+removed from the dynamics by a projection in spectral space
+\cite{craya_contribution_1957}. It is worth noting that numerical simulations of
+reduced dynamical equations (i.e., without non-wave structures) of stratified rotating
+flows in the hydrostatic balance have already been conducted to remove non-wave
+structures \cite{lvov_nonlinear_2009}. Despite the simplifications, these simulations
+reproduced some key features of oceanic internal-wave spectra, such as the accumulation
+of energy at near-inertial waves and realistic frequency and horizontal wave-number
+dependencies of spatiotemporal spectra. In the present work, we do not account for
+rotation. However, our simulations allows to investigate the role of vortical modes on
+the dynamics of stratified flows outside the hydrostatic balance approximation.
+
+
+The manuscript is organized as follows. In section \ref{sec:methods}, we present our
+methodology including a presentation of the code and the simulations. Our results are
+presented in section \ref{sec:results}. Subsection \ref{subsec:global} is devoted to
+the study of flow regimes in the control parameter space. It shows, as expected, that
+WWT is naturally more likely to occur at high stratification and without vortical
+modes. Subsections \ref{subsec:khkz} and \ref{subsec:khkzomega} deal with the
+spatiotemporal analysis of a couple of strongly stratified simulations to investigate
+further the presence of linear waves in spatial scales. In the last subsection
+\ref{subsec:waves_energy}, we propose a simple diagnostic to evaluate the predominance
+of waves in the parameter space. We give concluding remarks in section
+\ref{sec:conclusions}.
+
+\section{Methods}
+\label{sec:methods}
+
+We use the 3D Navier-Stokes equations under the Boussinesq approximation: \Remove{with
+added hyperviscosity and hyperdiffusion:}
+\begin{align}
+\label{eq:Continuity}
+\bnabla \cdot \vv &= 0 \\
+\label{eq:Impulsion}
+\p_t\vv + \vv \cdot \bnabla \vv &= b ~ \boldsymbol{e}_z - \bnabla p +
+\nu \nabla^2\vv + \ff, \\
+\label{eq:Buoyancy}
+\p_t{b} + \vv \cdot \bnabla b &= -N^2 v_z + \kappa\nabla^2{b} ,
+\end{align}
+where $(x,y,z)$ will represent the three spatial coordinates in the cartesian frame
+$(O, \eex, \eey, \eez)$, $\eez$ is the stratification axis, $\vv=(v_x, v_y, v_z)$ is
+the velocity, $b$ the buoyancy, $p$ the total kinematic pressure, $N$ the \bv
+frequency, $\nu$ the viscosity, $\kappa$ the diffusivity\Remove{, $\nu_4$ the
+hyperviscosity, $\kappa_4$ the hyperdiffusivity,} and $\ff$ is the velocity forcing.
+\Add{The buoyancy is defined as $b = -g \rho' / \rho_0$, where $g$ is the acceleration
+due to gravity, $\rho_0$ is the average density of the fluid at $z=0$, and $\rho'$ is
+the density perturbation with respect to the average linear density profile
+$\bar{\rho}(z) = \rho_0 + (\mathrm{d}\bar{\rho}/\mathrm{d}z)z$.} The Schmidt number $Sc
+= \nu/\kappa$ is fixed to one.
+
+\Add{We simulate forced-dissipated flows} with the pseudo-spectral solver
+\texttt{ns3d.strat} from the FluidSim software \cite{mohanan_fluidsim_2019} (an
+open-source Python package of the FluidDyn project \cite{fluiddyn} using Fluidfft
+\cite{fluidfft} to compute the Fast Fourier Transforms). \Add{The forcing, which will
+be described in details at the end of this section and in
+Appendix~\ref{appendix:forcing}, is computed in spectral space such that the kinetic
+energy injection rate $P_K$ is constant and equal to unity. The physical input
+parameters are the \bv frequency $N$ and the diffusive coefficients $\nu=\kappa$, but
+in practice, we identify our simulations with the couple $(N,\,\R_i)$, where $\R_i
+\equiv P_K / (\nu N^2)$ is the input buoyancy Reynolds number.} The turbulent
+non-dimensional numbers characterizing the statistically stationarity flow are the
+horizontal turbulent Froude number \Remove{, the Reynolds number,} and the buoyancy
+Reynolds number \cite[]{brethouwer_scaling_2007} that are respectively
+\begin{equation}
+\label{eq:FhR} F_h = \frac{\epsK}{{U_h}^2 N} \Remove{, ~~~~ Re =
+\frac{U_h^4}{\nu \epsK}} ~~~~ \text{and} ~~~~ \R = \frac{\epsK}{\nu N^2},
+\end{equation}
+\Add{where $\epsK$ is the kinetic energy dissipation rate and $U_h$ the rms of the
+horizontal velocity.} \Add{Note that the turbulent Reynolds number is given by $Re = \R
+/ F_h^2$.} We also compute the buoyancy and the Ozmidov wave-vectors $\kb \equiv N/U_h$
+and $\ko \equiv \sqrt{N^3 / P_K}$. All the quantities presented in this manuscript are
+computed from averaging when stationarity is reached. A list of the simulations, with
+relevant parameters and physical quantities, is given in
+Appendix~\ref{table-better-simuls}.
+
+In this study, we consider a periodic domain of horizontal size $L_x = L_y = L_h = 3$
+\Add{and vertical size} \Remove{of the domain,} $L_z$ \Remove{, is varied depending on
+the value of the \bv frequency}. We note $(n_x, n_y, n_z)$ the numbers of collocations
+points in the three spatial directions, with $n_x = n_y \equiv n_h$. \Add{We chose
+$n_z$ in order to have an isotropic mesh in physical space, i.e. $L_z/n_z = L_h/n_h$.}
+\Add{We decrease $L_z$ with $N$. Typically, $L_z \propto 1/N$, while $L_h$ is kept
+constant for all simulations. More precisely, the aspect ratio $L_z/L_h$ is $1/2$ for
+$N \leq 20$, $1/4$ for $N \leq 60$, and $1/8$ for $N \geq 80$. This choice is motivated
+by the fact that the unforced and undissipated Boussinesq equations are self-similar in
+the limit $F_h \rightarrow 0$, with similarity variable $zN/U$, where $U$ is the
+typical velocity \cite{billant_self-similarity_2001}. In that way, we simulate few
+layers (of height $L_b = U/N$) for all our simulations.}
+
+\Add{To reduce the computational costs, we first simulate the transient state with a
+coarse resolution with $n_h = 320$.}
+\Add{For these simulations, two hyperdiffusive terms $-\nu_4 \nabla^4 \vv$ and
+$-\kappa_4\nabla^4{b}$ are added to (\ref{eq:Impulsion}) and (\ref{eq:Buoyancy}),
+respectively, in order to keep the dissipative range in the simulated scales and avoid
+thermalisation at small scales.}
+\Add{Once a statistically steady-state is reached, we increase the resolution of the
+simulation while decreasing hyper-viscosity. Then the simulation is run until reaching
+a new statistically steady state. The previous step is repeated until reaching the
+highest resolution.}
+We measure the turbulent kinetic dissipation rates $\epsKK$ and $\epsKKKK$ based on
+both viscosities, and the total kinetic energy dissipation rate $\epsK = \epsKK +
+\epsKKKK$. The product of the maximal wave-vector $\kmax$ with the Kolmogorov scale
+$\eta \equiv (\nu^3 / \epsK)^{1/4}$ is computed to quantify how close our simulations
+are from true Direct Numerical Simulations (DNS). \Add{In practice, it is common to
+consider that simulations are proper DNS with well-resolved small scales when
+$\kmax\eta > 1$ \cite[]{deBruynKops1998,brethouwer_scaling_2007}.} \Add{For the
+statistically stationarity state, the time average of the total energy dissipation rate
+is equal to the injection rate $P_K$ so that the product $\kmax\eta$ depends mostly on
+$\nu$ and $n_h$.} \Add{For most couples $(N,\,\R_i)$, the resolution of the larger
+simulation is fine enough ($\kmax\eta \gtrsim 1$) so that the hyperdiffusion is zero or
+negligible.} \Add{For example, the two simulations analyzed in details in the next
+section (Figures~\ref{fig:buoyancy_fields} to \ref{fig:ratioEwaves_khkz}) are proper
+DNS with $\kmax\eta$ equal to 0.99 and 1.05, respectively (see
+table~\ref{table-better-simuls}).} \Add{There are also few simulations with $0.45 <
+\kmax\eta < 1$ (19 out of 78 simulations), which remain slightly under-resolved and
+affected by hyper-viscosity. In that case, small-scales and flow statistics should be
+analyzed carefuly. We checked that these simulations do not change the results
+presented here.}
+
+The Fourier transform of the velocity field $\hatvv = (\hat{v}_x, \hat{v}_y,
+\hat{v}_z)$ can be written using the poloidal-toroidal-shear decomposition (see e.g.
+\cite{craya_contribution_1957, smith_generation_2002, laval_forced_2003,
+godeferd_toroidalpoloidal_2010, kimura_energy_2012, maffioli_vertical_2017})
+\begin{equation}
+\hatvv = \begin{cases} \hatvp ~ \eep + \hatvt ~ \eet ~~~~ \text{if} ~ k_h \neq 0, \\
+	\vvs = \hat{v}_x ~ \eex + \hat{v}_y ~ \eey ~~~~ \text{if} ~ k_h = 0,
+\end{cases}
+\end{equation}
+where
+\begin{equation}
+\label{eq:poloidal-toroidal}
+\eek = \frac{\kk}{k}, ~~~~
+\eep = \frac{\kk \times (\kk \times \eez)}{|\kk \times (\kk \times \eez)|}, ~~~~
+\eet = \frac{\eez \times \kk}{|\eez \times \kk|}.
+\end{equation}
+$\hatvp$ is the poloidal component, $\hatvt$ the toroidal component, $\vvs$ the shear
+modes component, $\kk=(k_x, k_y, k_z)$ denotes the wave-vector, $k =|\kk| = \sqrt{k_x^2
++ k_y^2 + k_z^2}$ is its modulus, and $k_h = \sqrt{k_x^2 + k_y^2}$ is the modulus of
+the horizontal component of the wave-vector (Figure~\ref{fig:poloidal-toroidal}).
+\Add{Since the toroidal component $\hatvt$ corresponds only to the vertical vorticity
+($\hat{\Omega}_z = i k \hatvt \sin \thk$, with $\bOmega = \bnabla \times \vv$ being the
+vorticity), we also denote it as the ``vortical" velocity.}
+
+\begin{figure}
+\centering
+\includegraphics[width=0.4\textwidth]{Figure1}
+\caption{Illustration of the poloidal-toroidal basis $(\eek, \eep, \eet)$ defined by
+equations (\ref{eq:poloidal-toroidal}). $\thk$ is the angle between $\eez$ and $\eek$.
+$\phk$ is the angle between the horizontal projection of $\kk$ and $\eex$.
+\label{fig:poloidal-toroidal}}
+\end{figure}
+
+% Description of linear modes and shear modes dynamics
+In spectral space, once projected, the equations of motion
+(\ref{eq:Continuity}-\ref{eq:Buoyancy}) read
+\begin{align}
+\label{eq:StratifiedSpectralPoloidalToroidal}
+\begin{cases}
+	\dot{\hat{v}}_p &= - (\widehat{\vv \cdot \bnabla \vv}) \cdot \eep - \hatb \sin \thk - \nu k^2 \hatvp + \hatff \cdot \eep \\
+	\dot{\hat{v}}_t &= - (\widehat{\vv \cdot \bnabla \vv}) \cdot \eet - \nu k^2 \hatvt  + \hatff \cdot \eet \\
+	\dot{\hatb} &= - \widehat{\vv \cdot \bnabla b} + N^2 \hatvp \sin \thk - \kappa k^2 \hatb
+\end{cases}
+~~~~~~~~ \text{for }~ k_h \neq 0
+\end{align}
+and
+\begin{align}
+\label{eq:StratifiedSpectralShearModes}
+\begin{cases}
+	\dot{\hat{\vv}}_{s} &= - \widehat{\vv \cdot \bnabla \vv_h} - \nu k^2 \vvs, \\
+	\dot{\hatb} &= - \widehat{\vv \cdot \bnabla b} - \kappa k^2 \hatb
+\end{cases}
+~~~~~~~~ \text{for }~ k_h = 0,
+\end{align}
+where $\vv_h = (v_x,v_y,0)$ is the horizontal velocity and $\widehat{\cdot}$ denotes
+the Fourier transform. We will also note $k_i \in \Delta k_i ~ \mathbb{Z}$ with $\Delta
+k_i = 2\pi /L_i$ for $i =x,y,z$, and hereafter assume  $\Delta k_x = \Delta k_y \equiv
+\Delta k_h$. Without forcing and dissipation, the linearized equations
+(\ref{eq:StratifiedSpectralPoloidalToroidal}-\ref{eq:StratifiedSpectralShearModes}) can
+be written as
+
+\begin{align}
+\label{eq:StratifiedSpectralPoloidalToroidalLinear}
+\dot{\hata} &= -i \ok \hata, ~~~~ \dot{\hata}^{(0)} = 0, ~~~~ \text{if} ~~ k_h \neq 0 \\
+\label{eq:StratifiedSpectralShearModesLinear}
+\dot{\hatvv}_{s} &= \mathbf{0}, ~~~~
+\dot{\hatb} = 0 ~~~~ \text{if} ~~ k_h = 0,
+\end{align}
+where
+\begin{equation}
+\hata = \frac{\hatvp - i  \frac{\hatb}{N}}{\sqrt{2\ok}}, ~~~~ \text{and} ~~~~ \ok = N \frac{k_h}{k} = N \sin \thk
+\end{equation}
+are the waves modes and the pulsation of the waves, $\thk$ is the angle between $\kk$
+and the stratification axis $\eez$, and $\hata^{(0)} \propto \hatvt$ are the vortical
+modes\Remove{(corresponding to vertical vorticity)}. Equations
+(\ref{eq:StratifiedSpectralPoloidalToroidalLinear}-\ref{eq:StratifiedSpectralShearModesLinear})
+show that both shear modes and vortical modes have zero frequency. This means that
+linear waves can only exist in the poloidal velocity and the buoyancy, but not in the
+toroidal velocity nor the shear velocity.
+
+We are motivated by forcing internal gravity waves, which only involve the poloidal
+part of the velocity field and have an anisotropic dispersion relation. Therefore, we
+use an anisotropic, poloidal velocity forcing $\hatff = \hatf \, \eep$. The flow is
+forced at large spatial scales $ \left\{\kk ~ | ~ 5 \leq k/\Delta k_h \leq 20 \right\}$
+and small angle $\left\{\kk ~ | ~ |\ok /N - \sin \theta_f| \leq 0.05 \right\}$ where
+$\sin \theta_f = 0.3$, meaning that relatively slow internal waves are forced. The
+forcing scheme is described in Appendix~\ref{appendix:forcing}. \Add{It is neither
+harmonic nor given by a stochastic differential equation. Instead, a time correlated
+forcing is computed via generations of pseudo random numbers and time interpolations.}
+Its correlation time is equal to the period of the forced waves $T_c = 2\pi /(N \sin
+\theta_f)$. \Add{The forcing is normalized such that the kinetic injection rate $P_K$
+is always equal to 1.} \Add{Forcing slow waves is motivated by oceanic applications,
+where waves are generated, among other processes, by slow tides
+\cite{mackinnon_climate_2017, nikurashin_legg_mechanism_2011}. Low frequency forcing is
+also used in order to have a scale separation between forced frequencies and the \bv
+frequency so that one can potentially reproduce features of the oceanic temporal
+spectra close to $N$.}
+
+The time advancement is performed using the $4^{th}$ order Runge-Kutta scheme. All
+modes with wave-number modulus larger than $\kmax = 0.8 (n_h/2) \Delta k_h$ are
+truncated to limit aliasing. \Add{We checked that this 0.8 spherical truncation is a
+good compromize consistent with our other numerical choices.} Shear modes and
+vertically invariant vertical velocity (internal waves at $\omega = N$), which are
+absent in flows bounded by walls, are also removed in our simulations \Add{by fixing
+nonlinear transfers to these modes to zero}. \Add{For the simulations without vortical
+modes, the toroidal projection of the velocity and of the nonlinear transfers are set
+to zero. \citet{smith_generation_2002} used a similar procedure in order to distangle
+the roles of potential vorticity modes and waves modes in rotating stratified
+turbulence.}
+
+\Add{Without dissipation and forcing, equations (\ref{eq:Continuity}-\ref{eq:Buoyancy})
+conserve the total energy $E = \int ~ \left[ \vv^2 / 2 + b^2 / (2N^2) \right] ~ \diff x
+\diff y \diff z$ and the potential vorticity $\Pi = \bOmega \cdot \left( N^2 \eez +
+\bnabla b \right)$ is a Lagrangian invariant \cite{bartello_geostrophic_1995}. It
+follows that the spatial average of any function of $\Pi$ is conserved. As a special
+case, the potential enstrophy}
+\begin{align}
+V &\equiv \frac{1}{2} \int ~ \Pi^2 ~ \dxdydz \\ &=
+\frac{1}{2} \int ~ N^4 \Omega_z^2 ~ \dxdydz + \int ~ N^2 \Omega_z
+\bOmega \cdot \bnabla b ~ \dxdydz + \frac{1}{2} \int ~ \left( \bOmega
+\cdot \bnabla b \right)^2 ~ \dxdydz \\ &\equiv V_2 + V_3 + V_4
+\end{align}
+\Add{is an invariant of equations (\ref{eq:Continuity}-\ref{eq:Buoyancy}) in absence of
+dissipation and forcing. For a flow without vertical vorticity $V_2 = V_3 =0$.}
+
+\section{Results}
+\label{sec:results}
+
+\subsection{Regimes and global energy distribution}
+\label{subsec:global}
+
+In this subsection, we study the anisotropy of the flow and global energy budget in the
+$(F_h, \R)$ plane. To discuss stratified turbulence regimes, it is useful to introduce
+a large and a small scale isotropy coefficients \cite{linares_numerical_2020}. The
+former is based on the kinetic energy components
+\begin{equation}
+\label{eq:Ivelo}
+\Ivelo = \frac{3 E_{\rm kin,z}}{E_{\rm kin}}
+\end{equation}
+where $E_{\rm kin,z}$ is the vertical velocity energy and $E_{\rm kin}$ is the total
+kinetic energy. The small scale isotropy coefficient is computed using the kinetic
+energy dissipation rates, namely
+\begin{equation}
+\label{eq:Idiss}
+\Idiss = \frac{1 - \varepsilon_{\rm kin,z}/\varepsilon_{\rm kin}}{(1 - 1/3)}
+\end{equation}
+where $\varepsilon_{\rm kin,z}$ is the energy dissipation rate due to vertical
+gradients. Both $\Ivelo$ and $\Idiss$ are equal to unity for an isotropic flow.
+Conversely, $\Ivelo$ and $\Idiss$ should be close to zero if the flow is strongly
+anisotropic. When vortical modes are removed, we rather use
+\begin{equation}
+\label{eq:IcoeffsProj}
+	\Ivelo  = \frac{2 E_{\rm kin,z}}{E_{\rm kin}} \Remove{, ~~~~ \text{and}
+~~~~ \Idiss = \frac{1 - \varepsilon_{\rm kin,z}/\varepsilon_{\rm kin}}{(1 - 1/2)}}
+\end{equation}
+because it is expected that only half of the kinetic energy is contained in
+\Remove{(respectively dissipated by)} the vertical velocity field \Remove{(respectively
+the vertical gradients) if the energy spectra are independent of the angles in this
+case}. Figure~\ref{fig:regimes} shows the variations of $\Ivelo$ and $\Idiss$ with
+$(F_h,\R)$. Small points correspond to strong anisotropic energy dissipation, while
+large points correspond to isotropic energy dissipation. Dark points correspond to
+strong large-scale isotropy, while light points correspond to isotropic flows. The
+combination of $\Ivelo$ and $\Idiss$ allows to distinguish between four regimes
+(Figure~\ref{fig:regimes}$\rm (a)$): the passive scalar regime where $\Ivelo \simeq 1$,
+the weakly stratified regime where $0.5 \lesssim \Ivelo \lesssim 1$, and the strongly
+stratified regimes where $\Ivelo \leq 0.5$. As explained in
+\cite{brethouwer_scaling_2007}, the strongly stratified flows fall in two regimes: The
+Layered Anisotropic Stratified Turbulence (LAST) regime where the dissipation is
+isotropic because a 3D turbulence range can develop, typically when $\R \geq 10$; The
+viscosity affected regime where the small dissipative scales remain affected by the
+anisotropy. WWT is foreseen at small $F_h$ at some unknown values of $\R$. We observe
+that removing vortical modes does not modify this picture (Figure~\ref{fig:regimes}$\rm
+(b)$). However, for a given $(F_h,\R)$, the isotropy coefficients tend to have bigger
+values when the vortical modes are removed.
+
+\begin{figure}
+\includegraphics[width=1.0\textwidth]{Figure2}
+\caption{Classification of regimes in our simulations using the isotropy coefficients.
+The large scale isotropy coefficient $\Ivelo$ is given by the color-scale, and the
+small scale isotropy coefficient $\Idiss$ by the size of the symbols. For simulations
+with vortical modes $\rm (a)$, definitions (\ref{eq:Ivelo}-\ref{eq:Idiss}) are used.
+For simulations without vortical modes $\rm (b)$, definitions (\ref{eq:IcoeffsProj})
+are used. The dotted blue lines correspond to $\R = 10$, $F_h = 0.14$, and $F_h=1$.
+\label{fig:regimes}}
+\end{figure}
+
+We denote
+\begin{equation}
+\label{eq:energies} E_{\rm pot}  = \frac{1}{2N^2} \sum\limits_{\kk}
+|\hatb|^2, ~~~~ E_{\rm polo}  = \frac{1}{2} \sum\limits_{\kk} |\hatvp|^2, ~~~~ E_{\rm
+toro}  = \frac{1}{2} \sum\limits_{\kk} |\hatvt|^2, ~~~~ \tilde{\mathcal{D}} =
+\frac{E_{\rm polo}  - E_{\rm pot} }{E_{\rm polo}  + E_{\rm pot} }
+\end{equation}
+which are, respectively, the potential energy, the poloidal kinetic energy, the
+vortical modes (toroidal) kinetic energy, and the relative difference between poloidal
+and potential energy. We also note $E = E_{\rm pot}  + E_{\rm polo}  + E_{\rm toro} $
+the total energy. Figure~\ref{fig:global-energy} shows the vortical modes energy ratio
+$E_{\rm toro} /E$ and $\tilde{\mathcal{D}}$  as a function of $F_h$ and $\R$. For
+waves, we expect to observe an equipartition between the poloidal kinetic energy and
+potential energy. Consequently, both $E_{\rm toro} /E$ and $\tilde{\mathcal{D}}$ should
+be close to zero for a system mainly composed by internal gravity waves. We observe
+that vortical modes energy becomes dominant at high stratification (low $F_h$) if they
+are not removed from the dynamics in these simulations with $\R \geq 0.1$
+(Figure~\ref{fig:global-energy} (a)). The ratio $E_{\rm toro} /E$ is the lowest at
+intermediate stratification $F_h \simeq 0.1-0.3$ and low values of $\R$. However, the
+same simulations are marked by a predominance of potential energy over poloidal energy,
+$\tilde{\mathcal{D}} < 0$ (Figure~\ref{fig:global-energy} (b)), meaning that these
+weakly stratified flows do contain other structures than waves. When vortical modes are
+removed, we can obtain flows with a global balance between poloidal and potential
+energies ($\tilde{\mathcal{D}} \simeq 0$) at high stratification. Consequently, a flow
+governed by weak nonlinear interactions between waves may be obtained at high
+stratification by removing vortical modes. In the two next subsections, we perform a
+spatiotemporal analysis of a couple of stongly stratified turbulent simulations with
+$(N,\,\R_i\equiv P_K/\nu N^2)=(40,20)$ and aspect ratio $L_z/L_h=1/4$.
+\Add{Figure~\ref{fig:buoyancy_fields} shows the buoyancy fields for these two
+simulations with the same color scale. We observe that the flow is layered in the
+vertical direction and that overturning (horizontal vorticity) is present with or
+without vortical modes. It is a standard feature of strongly stratified turbulence
+\cite{laval_forced_2003, lindborg_energy_2006, brethouwer_scaling_2007,
+waite_stratified_2011}. With vortical modes, the vertical vorticity is not zero
+(Figure~\ref{fig:buoyancy_fields}(a)) so the buoyancy has a different structure than
+when vortical modes are absent (Figure~\ref{fig:buoyancy_fields}(b)). Without vortical
+modes the dynamics in the horizontal direction is irrotational and the buoyancy field
+has a larger amplitude.}
+
+\begin{figure}
+\includegraphics[width=1.0\textwidth]{Figure3}
+\caption{Vortical modes to total energy ratio (a) and relative difference between
+poloidal and potential energies $\tilde{\mathcal{D}} = (E_{\rm polo}  - E_{\rm
+pot})/(E_{\rm polo}  + E_{\rm pot})$ (b) vs $F_h$ for simulations with or without
+vortical modes. The red boxes indicate the simulations with $(N,\R_i)=(40,20)$
+investigated in the next subsections. \label{fig:global-energy}}
+\end{figure}
+
+\begin{figure}
+\includegraphics[width=1.0\textwidth]{figure4}
+\caption{\Add{Snapshots of the buoyancy fields for simulations $(N,\R_i)=(40,20)$ with
+$\rm (a)$ and without $\rm (b)$ vortical modes.} \label{fig:buoyancy_fields}}
+\end{figure}
+
+\subsection{Energy budget in the strongly stratified regime}
+\label{subsec:khkz}
+
+It is expected that the spatial energy budget will depend on the ratio of temporal
+timescales. In particular the nonlinearity parameter used for several physical systems,
+including stratified and rotating turbulence \cite{nazarenko_critical_2011,
+yokoyama_energy-based_2019} and MHD \cite{meyrand_direct_2016, cerri_turbulent_2022} is
+\begin{equation}
+\label{eq:NonLinearityParameter} \chi_{\kk} \equiv \frac{\tau_{\rm L}}{\tau_{\rm NL}} =
+\frac{(\epsK k^2)^{1/3}}{N k_h / k} =   \frac{k}{k_h} \left(\frac{k}{\ko} \right)^{2/3}
+= \frac{1}{\sin \thk} \left(\frac{k}{\ko} \right)^{2/3}.
+\end{equation}
+It represents the ratio between the period of the linear wave, $\tau_{\rm L} = 2 \pi k
+/ (N k_h)$, and the eddy turnover time $\tau_{\rm NL} = 2 \pi / (\varepsilon_{\rm kin}
+k^2)^{1/3}$ for a given wave-vector. In a similar way, we introduce the wave
+dissipation parameter
+\begin{equation}
+\label{eq:DissipationParameter} \gamma_{\kk} \equiv \frac{\tau_{\rm
+L}}{\tau_{\nu}} = \frac{\nu k^2}{N k_h / k} =   \frac{k}{k_h} \left(\frac{k}{\kd}
+\right)^{2} = \frac{1}{\sin \thk} \left(\frac{k}{\kd} \right)^{2}
+\end{equation}
+which represents the ratio between the period of the linear wave and the dissipation
+time $\tau_{\nu} = 2 \pi / (\nu k^2)$. In the last equation, we have introduced the
+wave dissipation wave-vector $\kd \equiv \sqrt{N/\nu}$.  The kinetic energy and
+potential energy budgets for one Fourier mode read
+\begin{equation}
+\label{eq:seb}
+	\left\langle\frac{1}{2}\frac{\mathrm{d} |\hatvv|^2}{\mathrm{d}t} \right \rangle =
+\mathcal{I}_{\kk} + \mathcal{T}_{\rm kin,\kk} - \mathcal{B}_{\kk} - \varepsilon_{\rm
+kin,\kk},  ~~~~ \text{and} ~~~~
+\left\langle\frac{1}{2N^2}\frac{\mathrm{d}|\hatb|^2}{\mathrm{d}t} \right\rangle=
+\mathcal{T}_{\text{pot},\kk} + \mathcal{B}_{\kk} - \varepsilon_{\text{pot},\kk},
+\end{equation}
+where
+\begin{align}
+\nonumber &\mathcal{I}_{\kk} = \left\langle \Re
+\left(\hatff  \cdot \hatvv^* \right)\right \rangle, ~~~~ \mathcal{T}_{\rm kin, \kk} =
+- \left\langle \Re \left(\hatvv^* \cdot \left[\bar{\bar{P}}_{\kk}  \cdot (\widehat{\vv \cdot \bnabla \vv}) \right]\right)\right\rangle, \\
+\nonumber &\mathcal{T}_{\text{pot},\kk} = - \left\langle \Re \left(\hatb^* (\widehat{\vv \cdot \bnabla b}) \right)/ N^2 \right\rangle,
+~~~~ \mathcal{B}_{\kk} = - \left\langle \Re \left(\hat{v}_{z}^* \hatb \right) \right\rangle, \\
+&\varepsilon_{\rm kin, \kk} =  (\nu k^2 + \nu_4 k^4) \frac{|\hatvv|^2}{2}, ~~~~ \varepsilon_{\text{pot},\kk} = (\kappa k^2 + \kappa_4 k^4) \frac{|\hatb|^2}{2N^2},
+\end{align}
+are respectively the kinetic energy injection rate, the kinetic energy transfer, the
+potential energy transfer, the conversion of kinetic energy to potential energy, the
+kinetic energy dissipation, and the potential energy dissipation. In the last
+equations, $(\cdot)^*$ denotes the complex conjugate, \Add{$\Re(\cdot)$ the real part},
+$\left\langle\cdot \right\rangle$ stands for the averaging operator, and
+$\bar{\bar{P}}_{\kk} = \mathbb{I} - \eek \otimes \eek$ is the projector onto the plane
+orthogonal to $\kk$.  In this section, we study azimuthal average of the energy budget.
+Namely, we computed quantities like:
+
+\begin{equation}
+F(k_h, k_z) = \frac{1}{\Delta k_h ~ \Delta k_z} \mathop{\sum \sum}_{\substack{k_h \leq k_h'< k_h + \Delta k_h \\ k_z \leq |k_z'| < k_z + \Delta k_z}} ~ F_{\kk'},
+\end{equation}
+where $F$ can be $E_{\rm pot} $, $E_{\rm polo} $, $E_{\rm toro} $, $E_{\rm kin} =
+E_{\rm polo}  + E_{\rm toro} $, $E = E_{\rm kin} + \Epot$, $\mathcal{T}_{\text{pot}}$,
+$\mathcal{T}_{\rm kin}$, $\mathcal{B}$, $\mathcal{I}$, $\epsK$, or $\epsA$.
+
+Following \cite{yokoyama_energy-based_2019}, we use the energy ratios like
+\begin{equation}
+\frac{E_{\rm toro} (k_h, k_z)}{E(k_h, k_z)} ~~~~ \text{and}  ~~~~ \tilde{\mathcal{D}}(k_h,k_z) = \frac{\Epolo(k_h,k_z)-\Epot(k_h,k_z)}{\Epolo(k_h,k_z)+\Epot(k_h,k_z)},
+\end{equation}
+in order to quantify the energy content scale by scale. In the same spirit, we also
+introduce
+\begin{equation}
+\label{eq:transfers}
+\tilde{\mathcal{T}}_{\rm kin}(k_h,k_z) = \frac{\mathcal{T}_{\rm kin}(k_h,k_z)}{\mathcal{T}_{\text{tot}}(k_h,k_z)}, ~~~~
+\tilde{\mathcal{T}}_{\text{pot}}(k_h,k_z) = \frac{\mathcal{T}_{\text{pot}}(k_h,k_z)}{ \mathcal{T}_{\text{tot}}(k_h,k_z)} ~~~~
+\text{and} ~~~~
+\tilde{\mathcal{B}}(k_h,k_z) = \frac{\mathcal{B}(k_h,k_z)}{\mathcal{T}_{\text{tot}}(k_h,k_z)},
+\end{equation}
+where
+\begin{equation}
+\mathcal{T}_{\text{tot}}(k_h,k_z) \equiv |\mathcal{T}_{\rm kin}(k_h,k_z)| + |\mathcal{T}_{\text{pot}}(k_h,k_z)| + |\mathcal{B}(k_h,k_z)| + \epsK(k_h,k_z) + \epsA(k_h,k_z).
+\end{equation}
+These quantities are useful for tracking the energy pathways scale by scale. By
+construction, $\tilde{\mathcal{T}}_{\text{pot}}$, $\tilde{\mathcal{T}}_{\rm kin}$, and
+$\tilde{\mathcal{B}}$ vary between $-1$ and $1$ depending on the amplitude and
+direction of the energy transfer or conversion ($\tilde{\mathcal{B}}>0$ corresponds to
+conversion of kinetic energy to potential energy). In a statistically stationary state
+($\langle\cdot  \rangle = 0$) and in the inertial range ($\mathcal{I}(k_h,k_z) = 0$,
+$\epsK(k_h,k_z) \simeq 0$, and $\epsA(k_h,k_z) \simeq 0$), we should have
+$\mathcal{T}_{\rm kin}(k_h,k_z) \simeq - \mathcal{T}_{\text{pot}}(k_h,k_z) \simeq
+\mathcal{B}(k_h,k_z)$. For this reason, we will present only
+$\tilde{\mathcal{B}}(k_h,k_z)$.
+
+\subsubsection{Energy spectra}
+
+
+Figure~\ref{fig:spectra_omega} displays the temporal energy spectra for simulations
+with $(N,\R_i)=(40,20)$. In the case where vortical modes are present
+(Figure~\ref{fig:spectra_omega}$\rm (a)$), we observe that potential and poloidal
+energies dominate over toroidal energy for frequencies smaller than the \bv frequency,
+i.e. $\omega/N <1$. For $\omega/N \geq 1$, the three components of the energy have very
+similar spectra. When vortical modes are removed (Figure~\ref{fig:spectra_omega}$\rm
+(b)$), potential energy is almost in perfect equipartition with poloidal energy at a
+given frequency. The potential and poloidal energy spectra behave very similarly wether
+vortical modes are present or not. This might be explained by the fact that vortical
+modes energy is smaller than potential an poloidal energy at large temporal scales
+$\omega \leq N$. Interestingly, in Figure~\ref{fig:spectra_omega}$\rm (b)$ temporal
+spectra get closer to the Kolmogorov-Zakharov spectra $E(\omega) \sim \omega^{-3/2}$
+\cite{caillol_kinetic_2000, lvov_hamiltonian_2001} than the the high frequency limit of
+the Garrett-Munk spectra $E(\omega) \sim \omega^{-2}$ \cite{garrett_internal_1979}. The
+inertial range is larger when vortical modes are removed.
+
+\begin{figure}
+\includegraphics[width=1.0\textwidth]{figure5}
+\caption{Temporal energy spectra for simulations $(N,\R_i)=(40,20)$ with $\rm (a)$ and
+without $\rm (b)$ vortical modes. The orange dotted lines correspond to the minimal and
+maximal frequencies of the linear internal gravity waves in the forcing region.
+\label{fig:spectra_omega}}
+\end{figure}
+
+
+Figure~\ref{fig:spectra_1D} shows the 1D integrated spatial spectra, and the normalized
+(i.e. divided by the total energy dissipation $\varepsilon \equiv \epsK + \epsA$)
+energy fluxes along $k_h$ for the same simulations. We first observe that removing
+vortical modes does not change the behavior of the poloidal an potential energy spectra
+(Figure~\ref{fig:spectra_1D}$\rm (a)$-$\rm (b)$). Up to the the buoyancy wave-vector,
+the poloidal and potential energies dominate with spectra that are slightly steeper
+than $k_h^{-5/3}$. After the buoyancy scale (i.e. $k=\kb$), the toroidal energy
+spectrum starts to be of the same order than the poloidal energy spectrum, which are
+shallower than $k_h^{-5/3}$. Such changes in the spectral slope around the buoyancy
+scale were reported in earlier studies \cite{kimura_energy_2012,
+bartello_sensitivity_2013, augier_stratified_2015}. They could be due to shear
+instabilities occurring at those scales, leading to non-local energy transfers toward
+small scales \cite{brethouwer_scaling_2007, waite_stratified_2011,
+augier_stratified_2015}. The vertical spectra has a spectral slope between $-3$ and
+$-2$, which is similar to earlier simulations \cite{maffioli_vertical_2017}. Therefore,
+these 1D spectra do not show strongly stratified turbulence scaling $~k_h^{-5/3}$ and
+$k_z^{-3}$ for $\kb \leq k_h,k_z \leq \ko$ \cite{billant_self-similarity_2001,
+lindborg_energy_2006}. This might be due to the fact that these simulations do not have
+a sufficient scale separation between $\kb$ and $\ko$, which can be attained only at
+very small $F_h$ \cite{lindborg_energy_2006, bartello_sensitivity_2013}, or when
+considering only the largest horizontal scales \cite{maffioli_vertical_2017}. Indeed,
+we observe that vertical spectra tends to be steeper as $F_h$ is decreased, and are
+steeper than the scaling $k_z^{-2}$ observed is 2D numerical simulations at comparable
+$F_h$ \cite{linares_numerical_2020}. The kinetic energy flux $\Pikin(k_h)$ and the
+potential energy flux $\Pipot(k_h)$ start to show a plateau over almost a decade in
+these simulations (Figure~\ref{fig:spectra_1D}$\rm (c)$-$\rm (d)$). We observe that the
+dissipation starts to be important at the Ozmidov scale, meaning that these simulations
+lie between the LAST regime and the viscosity affected regime
+\cite{brethouwer_scaling_2007}. This explains why we do not observe an isotropic
+turbulence range (with energy spectra $\sim k^{-5/3}$) in these simulations. The main
+effect of removing vortical modes is to make $\Pipot(k_h)$ larger than $\Pikin(k_h)$.
+Also, the kinetic energy dissipation $\epsK(k_h)$ is almost equal to the potential
+energy dissipation $\epsA(k_h)$ when vortical modes are removed, as we expect for a
+system of internal gravity waves with unit Schmidt number \cite{reun_parametric_2018}.
+
+\begin{figure}
+\includegraphics[width=1.0\textwidth]{figure6}
+\caption{Compensated 1D spatial energy spectra for simulations $(N,\R_i)=(40,20)$ with
+$\rm (a)$ and without $\rm (b)$ vortical modes. Normalized $k_h$ energy fluxes for the
+same simulations with $\rm (c)$ and without $\rm (d)$ vortical modes. The orange dotted
+line corresponds to the maximal wave-vector modulus of the forcing region, the black
+dotted line to $\kb$, the black dashed line to $\ko$. \label{fig:spectra_1D}}
+\end{figure}
+
+
+Figure~\ref{fig:spectra_khkz_kin} shows slices of the $(k_h, k_z)$ kinetic energy
+spectrum. The potential and kinetic energy spectra have the same trends with respect to
+$k_h$ and $k_z$, so the potential energy spectrum is not presented. As in
+\cite{yokoyama_energy-based_2019}, we observe different scaling laws depending on the
+region in the $(k_h,k_z)$ plane. Namely, important differences between small $k_z$ and
+large $k_z$ for a given $k_h$ (Figure~\ref{fig:spectra_khkz_kin}$\rm (a)$-$\rm (b)$):
+\begin{itemize}
+\item For $k_h \ll \kb$, energy is accumulated at the lowest available $k_h$ and small
+$k_z$, as usual in stratified turbulence \cite{smith_generation_2002,
+laval_forced_2003, herbert_waves_2016}. At $k_z \ll \kb$ the spectra are close to a
+$k_h^{-2}$ dependency, which is consistent with integrated energy spectra reported in
+earlier studies \cite{waite_stratified_2011, kimura_energy_2012}, but not with the WWT
+predictions (\ref{eq:WWTpredictions}). For a fixed $k_z \gg \kb$, we rather observe
+$E_{\rm kin} \sim k_h^1$ which would correspond to an equipartition of energy in
+horizontal scales.
+
+\item For $k_h \gg \kb$, the spectrum starts to depend less and less on $k_z$, and the
+$E_{\rm kin}(k_h,k_z)$ slices eventually merge around $k_h \sim \ko$, when the eddy
+turnover time becomes less than the period of any linear waves, i.e. $\chi_{\kk} > 1$.
+\end{itemize}
+
+Important differences can also be noted when looking at $k_z$ slices
+(Figure~\ref{fig:spectra_khkz_kin}$\rm (c)$-$\rm (d)$):
+\begin{itemize}
+\item For $k_z \ll \kb$, the spectrum is close to $\sim k_z^{0}$, indicating that
+horizontal layers of height larger than the buoyancy scale are in equipartition of
+energy in these strongly stratified simulations.
+
+\item For $\kb \ll k_z \ll \ko$, the spectrum is very steep at small horizontal
+wave-vectors, while it remains flat at large horizontal wave-vectors.
+
+\item For $k_z \gg \ko$, the spectral slope starts to depend less and less on $k_h$.
+Yet, the $E_{\rm kin}(k_h,k_z)$ does not merge around $k_h \sim \ko$ since small
+horizontal scales are more energetic compared to the large horizontal scales.
+\end{itemize}
+
+Due to the reduced aspect ratio of our simulations, we do not observe a net separation
+between the forcing and the buoyancy scales. Still, the kinetic energy spectrum remains
+very similar to the one presented in \cite{yokoyama_energy-based_2019}. Interestingly,
+the trends of the kinetic energy spectra appear to be very similar, regardless of
+whether vortical modes are present or not. This is a first indication that the presence
+of vortical modes is not the only obstacle to observing internal gravity wave
+turbulence in realistic flows.
+
+\begin{figure}
+\includegraphics[width=1.0\textwidth]{figure7}
+\caption{Slices of the kinetic energy spectrum $E_{\rm kin}(k_h,k_z)$ for simulations
+$(N,\R_i)=(40,20)$ with and without vortical modes. The orange dotted line corresponds
+to the maximal wave-vector modulus of the forcing region, the black dotted line to
+$\kb$, the black dashed line to $\ko$, and the green dashed line to the dissipative
+wave-vector. (a) $E_{\rm kin}$ vs $k_h$ with vortical modes. (b) $E_{\rm kin}$ vs $k_h$
+without vortical modes. (c) $E_{\rm kin}$ vs $k_z$ with vortical modes. (d) $E_{\rm
+kin}$ vs $k_z$ without vortical modes. \label{fig:spectra_khkz_kin}}
+\end{figure}
+
+
+Figure~\ref{fig:spectra_khkz} shows the $(k_h,k_z)$ spectra. Obviously, the simulation
+without vortical modes has no energy in the toroidal velocity. \Add{When vortical modes
+are present, an important part of the energy is contained in one vortical mode with
+$(k_h,k_z) = (\Delta k_h, 2 \Delta k_z)$, corresponding to large, nearly vertically
+stacked shear layers (Figure~\ref{fig:spectra_khkz}$\rm (a)$). Energy then tends to be
+accumulated at the smallest horizontal wave vectors, close to shear modes.} When
+vortical modes are removed, energy is sill concentrated in the same wave-vectors, but
+in the form of poloidal an potential energy (Figure~\ref{fig:spectra_khkz}$\rm
+(d)$-$\rm (f)$). Except for this qualitative difference, the toroidal, poloidal, and
+potential energy spectra show the same trends.
+
+\begin{figure}
+\includegraphics[width=1.0\textwidth]{figure8}
+\caption{$(k_h,k_z)$ spectra for the simulations at $(N,\R_i)=(40,20)$, with and
+without vortical modes. The cyan dotted lines correspond to $\chi_{\kk}=1/3$ and
+$\chi_{\kk}=3$, the cyan dashed line to $k=\kb$, the magenta dotted line to
+$\gamma_{\kk}=1$, and the green dotted line to the dissipative scale. The orange box
+corresponds to the forcing region. $\rm (a)$ $\Etoro$ with vortical modes. $\rm (b)$
+$\Etoro$ without vortical modes. $\rm (c)$ $\Epolo$ with vortical modes. $\rm (d)$
+$\Epolo$ without vortical modes. $\rm (e)$ $\Epot$ with vortical modes. $\rm (f)$
+$\Epot$ without vortical modes. \label{fig:spectra_khkz}}
+\end{figure}
+
+
+In order to have a finer analysis of the distribution of energy in the $(k_h, k_z)$
+plane, it is useful to look at energy ratios $\Etoro/E$ and $\tilde{\mathcal{D}}$
+(\ref{eq:energies}) shown on Figure~\ref{fig:spectra_ratio_khkz}. One clearly observes
+that $\Etoro$ is never negligible away from the forcing region, and it is very
+important for small wave-vectors when the vortical modes are not removed
+(Figure~\ref{fig:spectra_ratio_khkz}$\rm (a)$-$\rm (b)$). Apart from this important
+difference, the spectral energy budgets of the two simulations appear to share striking
+similarities. Outside the forcing region, the poloidal energy is dominant
+($\tilde{\mathcal{D}} > 0$) between the cyan dotted lines $\chi_{\kk} = 1/3$ and
+$\chi_{\kk} = 3$ (\ref{eq:NonLinearityParameter}), corresponding to a region where the
+critical balance condition $\chi_{\kk} \sim 1$ is fulfilled
+(Figure~\ref{fig:spectra_ratio_khkz}$\rm (c)$-$\rm (d)$). The potential energy is
+dominant ($\tilde{\mathcal{D}} < 0$) when $\chi_{\kk} > 3$, i.e. in a region where
+eddies are faster than waves. We observe a good equipartition between potential and
+poloidal energy ($\tilde{\mathcal{D}} \simeq 0$) when $\chi_{\kk} < 1/3$, i.e. in the
+region where waves are faster than eddies. This indicates that the ratio of temporal
+scales $\chi_{\kk}$ is important when identifying ranges in anisotropic turbulence, as
+already explained in \cite{yokoyama_energy-based_2019}. In particular, it shows that a
+wave dominated range cannot lie above $\chi_{\kk} > 1/3$.
+
+\begin{figure}
+\includegraphics[width=1.0\textwidth]{figure9}
+\caption{Ratios of energy for the simulations at $(N,\R_i)=(40,20)$, with and without
+vortical modes. The cyan dotted lines correspond to $\chi_{\kk}=1/3$ and
+$\chi_{\kk}=3$, the cyan dashed line to $k=\kb$, and the green dotted line to the
+dissipative scale. The orange box corresponds to the forcing region. $\rm (a)$
+$\Etoro/E$ with vortical modes. $\rm (b)$ $\Etoro/E$ without vortical modes. $\rm (c)$
+$\tilde{\mathcal{D}} = (\Epolo -\Epot)/(\Epolo+\Epot)$ with vortical modes. $\rm (d)$
+$\tilde{\mathcal{D}}$ without vortical modes. A wave dominated region is expected for
+the simulation without vortical where $\chi_{\kk} < 1/3$, for which $\Etoro = 0$ and
+$\tilde{\mathcal{D}} \ll 1$. \label{fig:spectra_ratio_khkz}}
+\end{figure}
+
+\subsubsection{Conversion of kinetic to potential energy}
+
+
+We can observe, in Figure~\ref{fig:conversion}, that the conversion between potential
+energy and kinetic energy $\tilde{\mathcal{B}}$ is very similar with or without
+vortical modes. This is not surprising since the vertical velocity is fully contained
+in waves modes, and not in vortical modes. Naturally, the amplitude of
+$\tilde{\mathcal{B}}$ is important only when $\gamma_{\kk} < 1$, meaning that waves
+need to not be damped too strongly by viscosity or diffusivity in order to convert
+potential energy to kinetic energy (or conversely). Yet, we observe that
+$\tilde{\mathcal{B}}$ is non-zero and has a definite sign when $\gamma_{\kk} > 1$.
+Consequently, if waves certainly exist in this range, they cannot persist because their
+kinetic energy is converted into potential energy or vice-versa. Therefore, these
+simulations are unlikely to correspond to a WWT regime, even in the buoyancy range and
+without vortical modes. Nevertheless, we observe that $\tilde{\mathcal{B}}$ fluctuates
+in time in the buoyancy range ($k\leq \kb$).
+
+
+\begin{figure}
+\includegraphics[width=1.0\textwidth]{figure10}
+\caption{Normalized conversion to potential energy $\tilde{\mathcal{B}}$
+(\ref{eq:transfers}) for the simulations at $(N,\R_i)=(40,20)$, with $\rm (a)$ and
+without $\rm (b)$ vortical modes. The magenta dotted line corresponds to
+$\gamma_{\kk}=1$, \Add{and} the green dotted line to the dissipative scale \Remove{,
+and the continuous black line to $k_z=k_h$}. The orange box corresponds to the forcing
+region. \label{fig:conversion}}
+\end{figure}
+
+\subsection{Spatiotemporal analysis in the buoyancy range}
+\label{subsec:khkzomega}
+
+From the previous subsection, it is expected that waves can dominate (but cannot
+necessarily sustain) not too far from the buoyancy range ($k \leq \kb$). To further
+assess the presence and degree of nonlinearity of waves, we performed a spatiotemporal
+analysis for small k. As explained in the previous subsection, waves are marked by an
+equipartition between poloidal kinetic energy and potential energy. This motivates us
+to introduce the equipartition energy as
+\begin{equation}
+\Eequi(k_h, k_z) = 2 \min \left\{\Epolo(k_h, k_z), \Epot(k_h, k_z)\right\},
+\end{equation}
+in order to track the presence of waves in the $(k_h, k_z)$ plane. Indeed, $\Eequi$
+corresponds to the energy contained in $\Epot$ and $\Epolo$ that is in
+potential-kinetic equipartition. Consequently, $\Eequi$ encompass the energy of the
+waves. However, it is important to note that $\Eequi$ can also contain the energy of
+structures that are not linear waves. A more careful separation between waves and
+non-wave structures requires a 4D spatiotemporal filtering \cite{lam_partitioning_2020,
+lam_energy_2021}, which is expensive in term of computational time and data storage.
+Yet, the ratio $\Eequi(k_h,k_z) /E(k_h,k_z)$ can be used to track waves in the spectral
+space at a lower cost.
+
+
+Figure~\ref{fig:omega-k} shows slices of the temporal spectrum of the equipartition
+energy $\Eequi(k_h,k_z,\omega)$ at large spatial scales, namely $k_h = 25.1$ and
+$k_z=25.1$. We observe that the energy tends to be concentrated around the linear wave
+frequency only at large spatial scales, but is dispersed at smaller scales even in the
+simulations without vortical modes. It has been shown that this dispersion could be
+quantified by considering a Doppler effect due to the advection of internal waves by
+shear modes \cite{maffioli_signature_2020}. Despite the fact that shear modes are
+removed in our simulations, we can try to test this prediction by considering that a
+Doppler shift is due to a horizontal mean flow $\boldsymbol{U}$ of amplitude $U_h$
+\begin{equation}
+\label{eq:Doppler}
+\odoppler = \max\limits_{\phk} ~ \boldsymbol{U} \cdot\kk = k U_h \sin \thk = \frac{k}{\kb} \ok
+\end{equation}
+where $\phk$ is the angle between the horizontal projection of $\kk$ and $\eex$
+(Figure~\ref{fig:poloidal-toroidal}). This gives a reasonable explanation of the
+important dispersion of the energy in temporal scale at large horizontal wave-vectors
+(Figure~\ref{fig:omega-k}).
+
+\begin{figure}
+\includegraphics[width=1.0\textwidth]{figure11}
+\caption{Slices of $\Eequi(k_h,k_z,\omega)$ for the simulations at $(N,\R_i)=(40,20)$,
+with and without vortical modes. (a) $k_z=25.1$ with vortical modes. (b) $k_z=25.1$
+without vortical modes. (c) $k_h=25.1$ with vortical modes. (b) $k_h=25.1$ without
+vortical modes. The linear dispersion relation $\omega_{\boldsymbol{k}} = Nk_h/k$ is
+plotted in black, the lines $\ok \pm \delta \omega/2$ (\ref{eq:delta-omega}) in blue,
+and the yellow dashed lines correspond to the lines $\ok \pm \odoppler$
+(\ref{eq:Doppler}). \label{fig:omega-k}}
+\end{figure}
+
+
+To observe more precisely the dispersion in temporal scales, we represent the ratio
+$\Eequi(k_h, k_z, \omega)/ \max\limits_{\omega} E(k_h, k_z, \omega)$ for some
+$(k_h,k_z)$ on Figure~\ref{fig:omega-L}, which indicates how the energy is distributed
+among temporal scales for a given spatial scale. For a system of weakly nonlinear
+internal gravity waves, we should observe a peak of $\Eequi/\max\limits_{\omega} E$
+around $\omega = \ok$ (for the energy to be concentrated around the linear frequency),
+whose maxima should be close to unity (for the energy to be concentrated in waves
+modes). We observe that the concentration of energy around the linear frequency mostly
+depends on the ratio $k/\kb$. As expected, the pics of energy around $\omega = \ok$ is
+more pronounced when $k \ll \kb$. This is an indication that linear waves are important
+if $k \ll \kb$, with or without vortical modes. When $k \simeq \kb$, the energy spreads
+over a broader range of temporal scales. We note that removing vortical modes from the
+dynamics helps to concentrate $\Eequi$ in temporal scales since temporal spectra tend
+to be more sharp. Simulations without vortical modes naturally have bigger
+$\Eequi/\max\limits_{\omega} E$ (Figure~\ref{fig:omega-L}$\rm (b)$) when compared to
+simulations with vortical modes (Figure~\ref{fig:omega-L}$\rm (a)$).
+
+\begin{figure}
+\includegraphics[width=1.0\textwidth]{figure12}
+\caption{Spatiotemporal analysis of the simulations with $(N,\R_i)=(40,20)$.
+$\Eequi(k_h, k_z, \omega)/\max\limits_\omega E(k_h, k_z, \omega)$ for the simulation
+with vortical modes $\rm (a)$ and for the simulation without vortical modes $\rm (b)$.
+Only some couples $(k_h, k_z)$ are shown. Line colors corresponds to different values
+of $k/\kb$. The vertical black lines corresponds to $\omega = \ok$.
+\label{fig:omega-L}}
+\end{figure}
+
+
+A way to measure the dispersion of energy in temporal scales is to compute the
+deviation from linear waves frequency $\delta \omega$. We estimate it using the
+equipartition energy spectrum and defining the following measure,
+
+\begin{equation}
+\label{eq:delta-omega}
+\delta \omega(k_h, k_z) = \sqrt{ \frac{\sum\limits_{\omega} ~ (\omega - \ok)^2 ~ E_{\rm equi}(k_h, k_z, \omega)}{\sum\limits_{\omega} ~ E_{\rm equi}(k_h, k_z, \omega)}}.
+\end{equation}
+
+The quantity $\delta \omega / \ok$ is another way to estimate the strength of nonlinear
+interactions. Unlike the nonlinearity parameter $\chi_{\kk}$, which is defined by
+dimensional analysis, $\delta \omega/ \ok$ requires knowledge of the spatio-temporal
+spectra. $\delta \omega/ \ok$ is called nonlinear broadening, and is particularly
+important in the context of WWT: for $\delta \omega/ \ok \gg 1$, waves' dynamics is
+strongly affected by nonlinear interactions, while when $\delta \omega/ \ok \ll 1$,
+waves can propagate with only weak nonlinear perturbations and the theory can hold
+\cite{nazarenko_wave_2011}. For a system of non-interacting linear waves, $\delta
+\omega/\ok$ should be zero. Figure~\ref{fig:nonlinear} shows this quantity in the
+$(k_h, k_z)$ plane. It shows that the nonlinear broadening is small if $k \ll \kb$ and
+$\chi_{\kk} \leq 1/3$, which is consistent with previous results
+\cite{yokoyama_energy-based_2019}. Conversely, $\delta \omega/\ok$ is large when $k \gg
+\kb$ or $\chi_{\kk} > 3$. We observed that removing vortical tends to decrease slightly
+$\delta \omega/\ok$ at small $k$ and large $\thk$. This might be explained by the fact
+that the simulations without vortical modes tend to have a larger buoyancy range when
+compared to the simulations with vortical modes.
+
+\begin{figure}
+\includegraphics[width=1.0\textwidth]{figure13}
+\caption{$\delta \omega/\ok$ for simulations with $(N,\R_i) = (40,20)$ with vortical
+modes (a) and without vortical modes (b). Dotted lines correspond to $\chi_{\kk} = 1/3$
+and $\chi_{\kk} = 3$, the dashed line to $k=\kb$, and the orange box to the forcing
+region.  \label{fig:nonlinear}}
+\end{figure}
+
+
+In order to allow a continuum of interaction between waves, it is also important that
+the nonlinear broadening remains much larger than the frequency gap between discrete
+modes in the Fourier space \cite{lvov_discrete_2010}. This leads to the condition
+$\Delta \omega \equiv \bnabla_{\kk} \ok \cdot \Delta \kk \ll \delta \omega$ where
+$\Delta \kk = (\Delta k_x, \Delta k_y, \Delta k_z)$. Taking $\Delta k_x = \Delta k_y =
+\Delta k_z = \Delta k$ leads to $|\cos \thk| \Delta k/ k \ll \delta \omega /N$, which
+is satisfied in most of our simulations for sufficiently large $k$.
+
+\subsection{Wave energy in the parameters space}
+\label{subsec:waves_energy}
+
+The previous subsections presented the dispersion in temporal scales for two
+simulations, but we did not quantify the amount of waves in the $(k_h,k_z)$ plane, nor
+in the $(F_h,\R)$ plane so far. This is the goal of this subsection.
+
+To discuss the effect of stratification, we show the integrated spatiotemporal total
+energy spectra as a function of $\thk/N = \sin \thk$ and $\omega$ for different
+simulations in Figure~\ref{fig:omega_vs_omegak}. At low $N$, we observe no big
+difference in $E(\ok, \omega)$ between simulations with and without vortical modes
+(Figure~\ref{fig:omega_vs_omegak}$\rm (a)$-$\rm (b)$). We observe that, for these
+weakly stratified flows, the energy is mostly contained in small frequency modes, but
+present over a large range of temporal scales for a given angle of the wave-vector
+$\thk$. At larger $N$, energy concentrate around the linear dispersion relation $\omega
+= \ok$ and in slow modes $\omega = 0$ for simulation with vortical modes
+(Figure~\ref{fig:omega_vs_omegak}$\rm (c)$). When vortical modes are absent, energy
+tends to accumulate around $\omega = \ok$ only (Figure~\ref{fig:omega_vs_omegak}$\rm
+(d)$). At the highest $N$, the two branches $\omega = \ok$ and $\omega = 0$ are visible
+for the simulation with vortical modes (Figure~\ref{fig:omega_vs_omegak}$\rm (e)$)
+while only the branches $\omega = \ok$ is clearly observed in the simulation without
+vortical modes (Figure~\ref{fig:omega_vs_omegak}$\rm (f)$).
+
+\begin{figure}
+\includegraphics[width=1\textwidth]{figure14}
+\caption{Integrated spatiotemporal spectra of the total energy as a function of $\ok/N
+= \sin \thk$ and $\omega$ for different simulations at the same viscosity $\nu = 1/(N^2
+\R_i)$. $\rm (a)$ $(N,\R_i)=(20,160)$ with vortical modes, $\rm (b)$
+$(N,\R_i)=(20,160)$ without vortical modes, $\rm (c)$ $(N,\R_i)=(40,40)$ with vortical
+modes, $\rm (d)$ $(N,\R_i)=(40,40)$ without vortical modes, $\rm (e)$
+$(N,\R_i)=(80,10)$ with vortical modes, and $\rm (f)$ $(N,\R_i)=(80,10)$ without
+vortical modes. The white dotted line represents $\omega  = \ok$, i.e. the linear
+dispersion relation. \label{fig:omega_vs_omegak}}
+\end{figure}
+
+In order to quantify wave energy in a flow, we propose to simply filter the
+spatiotemporal spectra $\Eequi$ by a Gaussian weight with mean $\ok$ and
+root-mean-square  $\epsilon \ok$, where $\epsilon$ is an arbitrary but small parameter,
+such that the wave energy at a given $(k_h,k_z,\omega)$ is
+\begin{equation}
+\label{eq:Ewave}
+\Ewave(k_h,k_z,\omega) =  \Eequi(k_h,k_z,\omega) ~ \exp \left[-\frac{1}{2} \left(\frac{\omega - \ok(k_h,k_z)}{ \epsilon \ok(k_h,k_z)} \right)^2 \right].
+\end{equation}
+We observed that choosing $\epsilon=0.1$ gives a sufficiently narrow window to keep
+only the spatiotemporal region corresponding to internal gravity waves, while keeping
+the window sufficiently wide to avoid binning effects. Therefore, we kept this value
+for the analysis. This window with variable width allows us to get rid of the energy of
+non-wave structures at low frequency. We first look at the distribution of the wave
+energy in the $(k_h,k_z)$ plane. In order to allow a fair comparison between
+simulations with and without vortical modes, it is necessary to introduce the wave
+energy ratio
+\begin{equation}
+\label{eq:WaveRatio}
+\Ewaver(k_h,k_z) =  \frac{\Ewave(k_h,k_z)}{\Epolo(k_h,k_z) + \Epot(k_h,k_z)} ~~~~ \text{where} ~~~~ \Ewave(k_h,k_z) = \sum\limits_{\omega} ~ \Ewave(k_h,k_z,\omega).
+\end{equation}
+
+Normalizing the wave energy by the total energy (i.e. including vortical modes energy)
+would have mask the presence of waves in the simulation with vortical modes, since they
+represent an important part of the energy (Figure~\ref{fig:spectra_ratio_khkz}$\rm
+(a)$). The wave energy ratio $\Ewaver(k_h,k_z)$ is depicted on
+Figure~\ref{fig:ratioEwaves_khkz}$\rm (a)$-$\rm (b)$. We observe that $\Ewaver$ tends
+to be higher for simulations without vortical modes, but removing vortical modes does
+not change the variations of $\Ewaver$ in the $(k_h,k_z)$ plane for these simulations.
+As we could have expected from previous works \cite{yokoyama_energy-based_2019} and the
+previous subsections, $\Ewaver$ is higher in the buoyancy range, in particular where
+$\chi_{\kk} \leq 1/3$. Yet, $\Ewaver$ tends to be bigger at very specific angles,
+smaller than the forcing angle $\theta_f$. This accumulation of internal gravity wave
+energy at specific propagation angles has been reported in earlier studies (see e.g.
+\cite{maffioli_signature_2020}). In the present simulations, we can explain this by
+invoking Triadic Resonance Instabilities (TRI) between internal gravity waves
+\cite{brouzet_internal_2016}, in particular the PSI
+\cite{mccomas_resonant_1977,muller_nonlinear_1986}. Considering that a primary wave of
+frequency $\omega_f = \omega_0^* = N \sin \theta_f$ is excited by the forcing, and then
+decays into two daughter waves of the same frequency $\omega_1^*$ according to the PSI
+mechanism, we can deduce $\omega_1^*$ and the associated propagation angle
+$\theta_1^*$:
+\begin{equation}
+\label{eq:harmonic1}
+\omega_1^* = \frac{\omega_0^*}{2} = 0.15~N ~~~~ \Rightarrow
+~~~~\theta_1^* = \arcsin \left(\frac{\omega_1^*}{N} \right)\simeq 0.15.
+\end{equation}
+Invoking, for a second time, the PSI mechanism for a wave of frequency $\omega_1^*$, we
+can deduce a second harmonic
+\begin{equation}
+\label{eq:harmonic2}
+\omega_2^* = \frac{\omega_1^*}{2} = 0.075~N ~~~~ \Rightarrow
+~~~~\theta_2^* = \arcsin \left(\frac{\omega_2^*}{N} \right)\simeq 0.075.
+\end{equation}
+A third harmonic is obtained by considering that two waves of frequencies $\omega_1^*$
+and $\omega_2^*$ interact to give a third wave with frequency
+\begin{equation}
+\label{eq:harmonic3}
+\omega_3^* = \omega_1^* + \omega_2^* = 0.225~N ~~~~ \Rightarrow
+~~~~\theta_3^* = \arcsin \left(\frac{\omega_3^*}{N} \right)\simeq 0.225.
+\end{equation}
+$\theta_f$, $\theta_1^*$, $\theta_2^*$, and $\theta_3^*$ are reported on
+Figure~\ref{fig:ratioEwaves_khkz}$\rm (a)$-$\rm (b)$. We see that they reproduce well
+with the higher values of $\Ewaver$ observed at small angles. This agreement suggests
+that waves at small frequencies are excited by TRI in our strongly stratified
+simulations.
+
+\begin{figure}
+\includegraphics[width=1.0\textwidth]{figure15}
+\caption{Wave energy ratio $\Ewaver(k_h,k_z)$ (\ref{eq:WaveRatio}) for simulations with
+$(N,\R_i)=(40,20)$ with vortical modes $\rm (a)$, and without vortical modes $\rm (b)$.
+The cyan dotted lines correspond to $\chi_{\kk}=1/3$, and the cyan dashed line to
+$k=\kb$. The orange box corresponds to the forcing region. The forcing angle $\theta_f$
+is represented by a full orange line, while first harmonics
+(\ref{eq:harmonic1}-\ref{eq:harmonic3}) excited by the TRI are indicated by dashed
+orange lines. \label{fig:ratioEwaves_khkz}}
+\end{figure}
+
+The spatiotemporal analysis performed in this subsection allows to give a more precise
+measure of the dominance of waves in the spectral space at a given $(F_h,\R)$. The
+total wave energy ratio is defined as
+\begin{equation}
+\Ewaver = \frac{\sum\limits_{k_h,k_z,\omega} ~ \Ewave(k_h,k_z,\omega)}{\sum\limits_{k_h,k_z} ~ \Epolo(k_h,k_z) + \Epot(k_h,k_z)}.
+\end{equation}
+
+Figure~\ref{fig:ratioEwaves_Fh} shows $\Ewaver$ in the plane $(F_h, \R)$. Consistently
+with the previous subsections, we observe that removing the vortical modes helps to get
+higher values of $\Ewaver$, and thus, in principle, to get closer to a WWT regime. Yet,
+unstratified flows ($F_h \geq 1$) have almost no energy in waves modes (i.e. $\Ewave
+\ll 1$), and increasing the stratification (decreasing $F_h$) is required to increase
+the wave energy ratio (Figures~\ref{fig:ratioEwaves_Fh}$\rm (a)$ and
+\ref{fig:omega_vs_omegak}), but it also tends to increases the vortical mode energy
+(Figure~\ref{fig:global-energy}$\rm (a)$) when they are not artificially removed from
+the dynamics. Figure~\ref{fig:ratioEwaves_Fh}$\rm (b)$ shows that $\Ewaver$ tends to
+decrease with $\R$ for sufficiently low $F_h$. This is consistent with the numerical
+simulations of stratified flows forced by tides \cite{reun_parametric_2018}, which show
+convincing signatures of a WWT regime at relatively small $\R$ when compared to more
+usual strongly stratified simulations. Some of our simulations at low $F_h$ and high
+$\R$ remain affected by hyperviscosity, so $\R$ may not accurately quantify the effect
+of dissipation in this case. However, discarding these simulations from the analysis
+still leads to the same conclusions.
+
+\begin{figure}
+\includegraphics[width=1\textwidth]{figure16}
+\caption{$\Ewaver$ as a function of $F_h$ $\rm (a)$ and $\R$ $\rm (b)$ for all our
+simulations with or without vortical modes. The red boxes indicate the simulations with
+$(N,\R_i)=(40,20)$ investigated in the previous subsections.
+\label{fig:ratioEwaves_Fh}}
+\end{figure}
+
+\section{Discussions and conclusions}
+\label{sec:conclusions}
+
+In order to investigate the conditions under which a weak internal gravity wave
+turbulence regime could occur, we performed direct numerical simulations of stratified
+turbulence without shear modes, and with or without vortical modes, at various Froude
+and buoyancy Reynolds numbers.
+
+We observed that removing vortical modes helps to have a better overall balance between
+poloidal kinetic energy and potential energy. However, the spatial spectra appear to
+behave very similarly with or without vortical modes in our simulations. The spectral
+energy budget reveals that \Add{there is no range in $(k_h,k_z)$ for which the}
+conversion between kinetic energy and potential energy \Remove{do not show
+fluctuations} \Add{fluctuates} around zero, as we would expect for a system of
+statistically stationary waves. Additionally, the conversion between potential energy
+and kinetic energy becomes small where the wave dissipation parameter $\gamma_{\kk}$
+(\ref{eq:DissipationParameter}) is larger than unity. Due to the anisotropy of
+stratified flows, this means that it exists a range of scales where waves are
+efficiently dissipated by viscosity, but not necessarily the vortices.
+
+A spatiotemporal analysis in the buoyancy range showed that increasing stratification
+is necessary to decrease the nonlinear broadening to the level required for a weak wave
+turbulence regime. As in \cite{yokoyama_energy-based_2019}, we observed that waves are
+present in a region delimited by the nonlinearity parameter $\chi_{\kk}$
+(\ref{eq:NonLinearityParameter}). More precisely, waves are more likely to dominate
+where $\chi_{\kk} < 1/3$, if the vortical modes are removed. We also observed evidences
+of the presence of slow waves ($k_z \gg k_h$) subject to TRI where $\chi_{\kk} \geq
+1/3$ in simulations with or without vortical modes. However, the nonlinear broadening
+is shown to be large for $k_z \gg k_h$, so slow waves are less susceptible to be
+described by the WWT theory. Although removing vortical modes is not enough to observe
+a weak wave turbulence regime, simulations without vortical modes produce $E(\ok,
+\omega)$ plots (Figure~\ref{fig:omega_vs_omegak}) characterized by a better
+concentration of energy in wave modes. Additionally, temporal spectra of strongly
+stratified simulations without vortical modes exhibit a larger inertial range in
+temporal scales (Figure~\ref{fig:spectra_omega}). The spectra are also compatible with
+the Kolmogorov-Zakharov spectrum offered by WWT. However, this spectrum is known to be
+mathematically unrealizable, so the agreement between this prediction and our
+simulations remains to be explained.
+
+Using a simple diagnostic to quantify wave energy contained in a stratified flow, we
+showed that that the buoyancy Reynolds number should not be too large in order to
+observe a WWT regime. To understand this, we propose the following conditions in order
+to limit the generation of non wave structures and to ensure the weak nonlinearity for
+the evolution of all energetic modes:
+\begin{itemize}
+\item The waves and eddies should be dissipated in the buoyancy range which leads to
+the conditions
+\begin{equation}
+\label{eq:DissInBuoyancy}
+\kd \leq \kb ~~~~ \Rightarrow ~~~~ \R \leq F_h ~~~~ \text{and} ~~~~ k_{\eta} \leq \kb ~~~~ \Rightarrow ~~~~ \R \leq F_h^{2/3}.
+\end{equation}
+if we assume that dissipation of eddies occurs at the Kolmogorov wave-vector $k_{\eta}
+= (\epsK/\nu^3)^{1/4}$. These conditions are necessary to avoid wave breaking and the
+development of 3D small-scale eddies.
+
+\item The nonlinearity parameter $\chi_{\kk}$ should remain small for all energetic
+modes. If we assume that dissipation occurs at $k \lesssim k_\eta$, this leads to the
+condition
+\begin{equation}
+\chi_{\rm max} \equiv \max\limits_{\substack{\kk \\ k \leq k_{\eta}}} \chi_{\kk}\leq 1 ~~~~ \Rightarrow ~~~~ \max\limits_{\substack{\kk \\ k \leq k_{\eta}}} \frac{1}{\sin \thk} \left(\frac{k}{\ko} \right)^{2/3} = \frac{k_\eta}{k_{\rm h, min}} \left(\frac{k_\eta}{\ko} \right)^{2/3} \leq 1
+\end{equation}
+where $k_{\rm h, min}$ represents the minimal horizontal wave-vector. If we consider a
+flow of finite size with no shear modes, $k_{\rm h, min} = 2\pi / L_h$ such that
+\begin{equation}
+\label{eq:smallchi}
+\chi_{\rm max} = \frac{k_\eta L_h}{2 \pi} \left(\frac{k_\eta}{\ko} \right)^{2/3} \leq 1 ~~~~ \Rightarrow ~~~~  \left(\frac{\epsK L_h^4}{\nu^3} \right)^{1/4} \R^{1/2} \lesssim 1.
+\end{equation}
+\end{itemize}
+Finally, we require the flow to be in a fully turbulent regime such that
+\begin{equation}
+\label{eq:LargeRe}
+Re = \R F_h^{-2} \gg 1 ~~~~ \text{and} ~~~~ \epsK \sim U_h^3 / L_h.
+\end{equation}
+Then, substituting the definition of $Re$ from equation (\ref{eq:FhR}) and condition
+(\ref{eq:LargeRe}) into equation (\ref{eq:smallchi}) yields
+\begin{equation}
+\label{eq:SmallChi}
+\R \lesssim F_h^{6/5}.
+\end{equation}
+
+
+To summarize, we propose that a weak internal gravity waves turbulence regime could
+occur only if $F_h$ and $\R$ are such that conditions (\ref{eq:DissInBuoyancy}),
+(\ref{eq:LargeRe}), and (\ref{eq:SmallChi}) are fulfilled. Note that the last
+thresholds are obtained by dimensional analysis, so they are defined up to
+dimensionless factors, which are set to one for simplicity. On
+Figure~\ref{fig:discussion}, we plot our simulations in the $(F_h,\R)$ plane, as well
+as the simulations of
+\Add{\cite{brethouwer_scaling_2007,waite_potential_2013,reun_parametric_2018,lam_energy_2021}}
+and some recent experiments \cite{rodda_experimental_2022}. For the simulations of
+\cite{reun_parametric_2018}, the Froude number is computed using our definition
+(\ref{eq:FhR}), and the buoyancy Reynolds number as $\R = Re F_h^2$ in order to compare
+to our work. Namely, we have used Table~1 in \cite{reun_parametric_2018} with $U_h =
+u_{\rm rms}$, $\epsK = \varepsilon_k$, and $Re = Re_0$. This leads to different values
+of $F_h$ and $\R$ than the ones presented by the authors. We observe that our
+simulations do not lie in the region corresponding to (\ref{eq:DissInBuoyancy}),
+(\ref{eq:LargeRe}), and (\ref{eq:SmallChi}). On the contrary,
+\cite{reun_parametric_2018} attained such a region in their simulations, which provides
+a tangible explanation for why they obtained better signatures of WWT. \Add{It is worth
+mentioning that \cite{reun_parametric_2018} did not have to remove shear nor vortical
+modes in their simulations to observe signatures of internal wave turbulence. It
+indicates that there must be some threshold below which both shear and vortical modes
+do not grow if not directly forced. Such a situation would be analogous to rotating
+flows where there is a threshold below which geostrophic modes do not grow
+\cite{reun_experimental_2019}.} \Add{We observe that \cite{waite_potential_2013} also
+attained very small $\R$. In this study, the authors forced vortical modes so their
+simulations are not well suited for WWT. Yet, they showed that potential enstrophy
+tends to be quadradic (i.e. $V \simeq V_2$) for $F_h, \R \ll 1$ and that $V_2 \propto
+\int ~ \Omega_z^2 ~ \dxdydz$ increases with $\R$ in their simulations. It suggests that
+vortical modes energy increases with $\R$ when $\R \lesssim 1$, as explained by
+\cite{lam_energy_2021}. Yet, the study of the wave energy ratio is not done while
+keeping $F_h$ or $\R$ constant in \cite{lam_energy_2021}. To our knowledge, additional
+studies are needed to confirm if the wave energy ratio is for $\R\ll 1$ and $F_h \ll 1$
+a decreasing function of $\R$, and if a threshold below which both shear and vortical
+modes are stable exists.}
+
+\begin{figure}
+\includegraphics[width=0.66\textwidth]{figure17.png}
+\caption{Simulations of the present study and
+\Add{\cite{brethouwer_scaling_2007,waite_potential_2013,reun_parametric_2018,lam_energy_2021}},
+and the experiments \cite{rodda_experimental_2022} in the $(F_h, \R)$ parameters space.
+The colored full lines corresponding to conditions (\ref{eq:DissInBuoyancy}),
+(\ref{eq:LargeRe}), and (\ref{eq:SmallChi}) (see legend). The blue dashed line
+corresponds to $Re = 500$. The colored region is where a weak wave turbulence regime is
+expected. \label{fig:discussion}}
+\end{figure}
+
+If the latter claims are true, they constitute very prohibitive conditions for
+observing WWT in numerical simulations and experiments. To illustrate this point, let
+us consider a turbulent flow at $Re = 10^3$ that typically requires a $\sim 1024^3$
+resolution. Then (\ref{eq:SmallChi}) stipulates that the flow need to attain $F_h \leq
+Re^{-5/4} = 1.78 \times 10^{-4}$ for the nonlinearity parameter to be small for all
+energetic modes. Yet, the kinetic time for a system of internal gravity waves is
+expected to grow as $F_h^{-2}$, meaning that very long simulations are required to
+attain the steady state, even for relatively small domains.
+
+
+\begin{acknowledgments}
+\Add{We thank two anonymous reviewers for their constructive feedback.} This project
+was supported by the Simons Foundation through the Simons collaboration on wave
+turbulence. Part of the computations have been done on the ``Mesocentre SIGAMM''
+machine, hosted by Observatoire de la Cote d'Azur. The authors are grateful to the OPAL
+infrastructure from Université Côte d'Azur and the Université Côte d’Azur's Center for
+High-Performance Computing for providing resources and support. This work was granted
+access to the HPC/AI resources of IDRIS under the allocation 2022-A0122A13417 made by
+GENCI.
+\end{acknowledgments}
+
+\appendix
+
+\section{Forcing scheme}
+\label{appendix:forcing}
+
+\Add{Our forcing is designed to excite waves and computed via generation of pseudo
+random numbers with uniform distribution and time interpolation. Namely, we use the
+following algorithm:}
+
+\begin{algorithm}[H]
+$t = 0$; $t_0 = 0$; \\
+Generate two complex random fields $\hat f_0(\kk)$ and $\hat f_1(\kk)$; \\
+$\hatff = \hat f_0 \, \eep$;\\
+Normalize $\hatff$ to ensure
+\Add{$P_K(t) = \sum\limits_{\kk} ~ \Re \left[ \hatff \cdot \hatvv^* + \frac{\Delta t}{2} |\hatff|^2 \right] =  1$};\\
+\While{$t \leq T$}{
+	$t = t + \Delta t$; \\
+	\If{$t - t_0 \geq T_c$}{
+		$t_0 = t$; \\
+		$\hat f_0 = \hat f_1$; \\
+		Generate $\hat f_1$;
+	}{}
+	$\hatff = \left\{ \hat f_0 - \dfrac{(\hat f_1 - \hat f_0)}{2} \left[\cos\left(\dfrac{\pi(t-t_0)}{T_c} \right)+ 1 \right]\right\}\, \eep$; \\
+	Normalize $\hatff$;
+}
+%\caption{}
+\end{algorithm}
+where $\Delta t$ is the time increment at each time step \Add{and $T$ is the final time
+of the simulation}. \Add{The random complex fields are built such that their inverse
+Fourier transform is real and they are null for unforced wavenumbers.}
+
+\section{List of simulations}
+
+\begin{table}
+\caption{Overview of the numerical and physical parameters used in the simulations.
+$\R_i = P_K/ \nu N^2$. $F_h$ and $\R$ are the turbulent horizontal Froude number and
+buoyancy Reynolds number, respectively, defined in (\ref{eq:FhR}).}
+\begin{tabular}{m{4cm}m{5.5cm}m{5.5cm}}
+\textbf{Control parameters} & \textbf{With vortical modes} & \textbf{Without vortical modes} \\
+\end{tabular}
+\input{../tmp/table_params.tex}
+\input{../tmp/table_better_simuls.tex}
+\input{../tmp/table_better_simuls_proj.tex}
+\label{table-better-simuls}
+\end{table}
+
+\newpage
+
+\bibliography{main}
+\end{document}
diff --git a/2022strat_polo_proj/review0/input/coverletter.tex b/2022strat_polo_proj/review0/input/coverletter.tex
--- a/2022strat_polo_proj/review0/input/coverletter.tex
+++ b/2022strat_polo_proj/review0/input/coverletter.tex
@@ -36,9 +36,9 @@
 \vspace{0.5cm}
 
 \paragraph{}{We have now completed the revision of our manuscript FF10204, entitled
-"Internal gravity waves in stratified flows with and without vortical modes." We have
-attached to this letter the revised version of the manuscript, a version with the
-changes highlighted, and a separate point-by-point response to each referee.}
+"Internal gravity waves in stratified flows with and without vortical modes." Attached
+to this letter is the revised version of the manuscript, a version with the changes
+highlighted, and a separate point-by-point response to each referee.}
 
 \paragraph{}{We have corrected the manuscript by improving the presentation of our
 methods. The corresponding section has been reorganized and strengthened. In
diff --git a/2022strat_polo_proj/review0/input/rebut1.tex b/2022strat_polo_proj/review0/input/rebut1.tex
--- a/2022strat_polo_proj/review0/input/rebut1.tex
+++ b/2022strat_polo_proj/review0/input/rebut1.tex
@@ -57,15 +57,15 @@
 \maketitle
 
 \noindent We thank the referee for their critical comments, that helped us to
-significantly improve our manuscript. We realized that our Methods section was not
-completely comprehensive. This is partly due to the fact that this paper was written to
-cite another article describing with another point of view the first dataset used in
-this study (without projection). Finally, we submitted this manuscript first and did
-not correctly complete the Methods section. We tried to correct the manuscript by
-better presenting our methods, and the corresponding section has been deeply
+significantly improve our manuscript. We have realized that our Methods section was not
+completely comprehensive. This is partly due to the fact that this paper was intended
+to cite another article describing, with another point of view, the first dataset used
+in this study (without projection). Finally, we submitted this manuscript first and did
+not correctly complete the Methods section. Now we have corrected the manuscript by
+presenting our methods better, and the corresponding section has been deeply
 reorganized and strengthened.
 
-The answer on all the comments are listed below. Corresponding corrections are made in
+The answers to all the comments are listed below. Corresponding corrections are made in
 blue in the new draft.
 
 \begin{enumerate}
@@ -76,25 +76,25 @@
 Instead of the exact 2/3 rule for cubic truncation, we use a 0.8 spherical truncation.
 
 
-In a non published work with Jason Reneuve, we performed an extensive study on
+In an unpublished work with Jason Reneuve, we have performed an extensive study on
 dealiasing methods (also including phase shifting). One notable result is that for
 simulations with $\kmax \eta$ slightly smaller than one, we obtain better results to
 reproduce known flows with a 0.8 spherical truncation rather than with the exact 2/3
-cubic truncation. First, a spherical truncation is better than a cubic truncation
+cubic truncation. First, the spherical truncation is better than the cubic truncation
 because it conserves the rotational symmetry. Second, for a constant resolution,
 increasing a bit the dealiasing coefficient gives more wavenumbers to represent the
 flow (for $\kmax \eta$ close to one, more wavenumbers to represent the dissipative
 range, in which the energy spectra decrease very quickly). On the one hand, for
-spherical truncation, 2/3 removes too much modes (Only $16\% = (4/3) \pi ((2/3)0.5)^3
+spherical truncation, 2/3 removes too many modes (Only $16\% = (4/3) \pi ((2/3)0.5)^3
 \%$ of the grid points are kept). On the other hand, the coefficient needs to be
-smaller than 0.94 to kill the double aliasing, which is the most annoying ones. A value
-of 0.8 seems to be a good compromize for which most of the aliasing is removed. These
-arguments are valid without hyperdiffusion, which was not used for our study with Jason
-Reneuve. However, using a bit of hyperdiffusion of course helps to limit the effect of
-the aliasing. Note also than anyway most simulations of our datasets are proper DNS
-with $\kmax\eta$ larger or close to one.
+smaller than 0.94 to kill the double aliasing, which is the most annoying one. The
+value of 0.8 seems to be a good compromize for which most of the aliasing is removed.
+These arguments are valid without hyperdiffusion, which was not used for our study with
+Jason Reneuve. However, using a bit of hyperdiffusion, of course, helps to limit the
+effect of the aliasing. Note also that, anyway, most simulations of our datasets are
+proper DNS with $\kmax\eta$ larger or close to one.
 
-Since these ideas will be published elsewhere (with serious checks) and are about a
+Since these ideas will be published elsewhere (with thorough checks) and present only a
 small detail of our numerical methods, we do not think that it makes sense to include
 them in the manuscript so we just write in the Methods section:
 
@@ -120,8 +120,8 @@
 Low frequency wave forcing is also used in order to allow a scale separation between
 forced frequencies and the Brünt-Väisälä frequency.
 
-We now motivate the low frequency forcing in the new version of the manuscript and
-included:
+We now motivate the low frequency forcing in the new version of the manuscript by
+adding the following text:
 
 \addtoman{\Add{Forcing slow waves is motivated by oceanic applications, where waves are
 generated, among other processes, by slow tides [25, 68]. Low frequency forcing is also
@@ -140,7 +140,7 @@
 as DNS I assume? How important it is in the simulations presented here, would they all
 blow up without it?}
 
-As already mentioned, the presentation of our methods was a strong weakness of our
+As already mentioned, the presentation of our methods was rather incomplete in our
 manuscript. In particular, we did not explain clearly how hyperdiffusivity is used. As
 a result the readers were lead to think that all our simulations are not proper DNS. We
 tried to better explain this aspect of our numerical methods. We first present the
@@ -170,7 +170,7 @@
 hyperdiffusion is zero or negligible.} \Add{For example, the two simulations analyzed
 in details in the next section (Figures~4 to 15) are proper DNS with $\kmax\eta$ equal
 to 0.99 and 1.05, respectively (see table~I).} \Add{There are also few simulations with
-$0.45 < \kmax\eta < 1$ (19 over 78 simulations), which remain slightly under-resolved
+$0.45 < \kmax\eta < 1$ (19 out of 78 simulations), which remain slightly under-resolved
 and affected by hyper-viscosity. In that case, small-scales and flow statistics should
 be analyzed carefuly. We checked that these simulations do not change the results
 presented here.}}
@@ -187,9 +187,9 @@
 simulations are not DNS. Secondly, keeping only DNS does not changes the results
 presented in the present study.
 
-Moreover, we decided to discard for this article few simulations for which $\kmax\eta<
+Moreover, we decided to discard few simulations in this article for which $\kmax\eta<
 0.45$, so that all simulations used for this study are either well-resolved or slightly
-under-resolved, with $\kmax$ close to $1/\eta$ and in the diffusive range linked to
+under-resolved, with $\kmax$ close to $1/\eta$ and in the diffusive range linked to the
 standard dissipation.
 
 Regarding the question about the importance of hyperviscosity and the potential ``blow
@@ -202,10 +202,10 @@
 Is it done via a spectrally localised damping term or simply put to zero at every time
 steps?}
 
-Shear modes are set to zero and not forced. For more, the shear modes part of the
-nonlinear term is removed when de-aliasing is applied at each time step. In that way,
-energy transfer to shear modes is forbidden. We change a sentence in the method section
-to make this point clearer:
+Shear modes are set to zero and not forced. Furthermore, the shear-mode part of the
+nonlinear term is removed when de-aliasing is applied at each time step. This way,
+energy transfer to the shear modes is forbidden. We have changed a sentence in the
+method section to make this point clearer:
 
 \addtoman{Shear modes and vertically invariant vertical velocity (internal waves at
 $\omega = N$), which are absent in flows bounded by walls, are also removed in our
@@ -224,9 +224,9 @@
 \caption{Time dependence of the forcing.} \label{fig:forcing}
 \end{figure}
 
-Figure~\ref{fig:forcing} shows how the forcing depends in time. We plan to include this
+Figure~\ref{fig:forcing} shows how the forcing depends on time. We plan to include this
 figure and a discussion about this time dependence in a future article studying
-specifically the dataset composed of simulations with vortical modes.
+specifically the dataset composed of simulations with the vortical modes.
 
 As described in Appendix A, the forcing is not harmonic, and not given by a stochastic
 differential equation either. The phase and the amplitude of the forced wavenumbers are
@@ -272,13 +272,13 @@
 \includegraphics[width=0.8\linewidth]{vt_filter}
 \caption{Dominant vortical mode.}
 \end{figure}
-It makes us remarks that, contrary to what we originally said, the dominant vortical
-mode consists of nearly vertical stacked shear layers with $(k_h, k_z) = (\Delta k_h, 2
-\Delta k_z)$. We have corrected this part:
+It makes us observe that, contrary to what we originally said, the dominant vortical
+mode consists of nearly vertically stacked shear layers with $(k_h, k_z) = (\Delta k_h,
+2 \Delta k_z)$. We have corrected this part:
 
 \addtoman{\Add{When vortical modes are present, an important part of the energy is
 contained in one vortical mode with $(k_h,k_z) = (\Delta k_h, 2 \Delta k_z)$,
-corresponding to large, nearly vertical stacked shear layers (Figure~8$\rm (a)$).
+corresponding to large, nearly vertically stacked shear layers (Figure~8$\rm (a)$).
 Energy then tends to be accumulated at the smallest horizontal wave vectors, close to
 shear modes.}}
 
@@ -307,7 +307,7 @@
 
 The purpose of the $k_z = k_h$ line is to delimitate regions where $k_z > k_h$ and $k_z
 < k_h$. In the end, it is not useful and burden figures so we decide to remove it. The
-figure's caption is adapted.
+figure caption is adapted accordingly.
 
 
 \itemit{Figure 16: this is an important plot in my opinion, it might be useful to
@@ -324,7 +324,7 @@
 the authors were forced to use larger viscosity at that time
 \cite[]{waite_stratified_2004, waite_stratified_2006, lindborg_energy_2006,
 waite_stratified_2011}. Unfortunately, it not always easy to correctly quantify these
-numbers from the data given in the articles, or to compare it with our simulations.
+numbers from the data given in the articles, or to compare them with our simulations.
 \cite{waite_stratified_2004} forced vortical modes so these it is not well suited for
 observing internal gravity wave turbulence. We do not see how to extract $\R$ from
 \cite{waite_stratified_2004} and \cite{waite_stratified_2006} tables. In
@@ -353,9 +353,9 @@
 or suppressed. Do the authors expect surprising behaviours with the shear modes but
 without the vortical modes?}
 
-We know thanks to a study on 2D turbulence \cite[]{linares_numerical_2020} that shear
+We know, thanks to a study on 2D turbulence \cite[]{linares_numerical_2020}, that shear
 modes can grow without toroidal modes. Therefore we can anticipate accumulation of
-energy in shear modes, which should become very strong and distore the waves, at least
+energy in shear modes, which should become very strong and distort the waves, at least
 for large $\R$.
 
 
diff --git a/2022strat_polo_proj/review0/input/rebut2.tex b/2022strat_polo_proj/review0/input/rebut2.tex
--- a/2022strat_polo_proj/review0/input/rebut2.tex
+++ b/2022strat_polo_proj/review0/input/rebut2.tex
@@ -59,14 +59,14 @@
 \begin{document}
 \maketitle
 
-\noindent The authors thank the referee for their useful and constructive comments that
+\noindent The authors thank the referee for useful and constructive comments that
 helped us to improve the manuscript. We realized that our Methods section was not
-completely comprehensive. It is partly due to the fact that this paper was written to
-cite another article describing with another point of view the first dataset used in
-this study (without projection). Finally, we submitted this manuscript first and did
-not correctly complete the Methods section. We tried to correct the manuscript by
-better presenting our methods, and the corresponding section has been deeply
-reorganized and strengthened.
+completely comprehensive. It is partly due to the fact that this paper was intended to
+cite another article describing, from another point of view, the first dataset used in
+this study (without projection). However, we submitted this manuscript first and
+therefore we need to detail and complete the Methods section of the present manuscript.
+We have now corrected the manuscript by presenting our methods better, and the
+corresponding section has been deeply reorganized and strengthened.
 
 Below we have responded to all questions. Corresponding corrections are made in blue in
 the new draft.
@@ -79,7 +79,7 @@
 vertical resolution. It should be stated clearly here how $L_z$ is varied with $N$ and
 why, and how $n_z$ is chosen.}
 
-We add the following explanations in the method section:
+We add the following explanations in the methods section:
 
 \addtoman{In this study, we consider a periodic domain of horizontal size $L_x = L_y =
 L_h = 3$ \Add{and vertical size} \Remove{of the domain,} $L_z$ \Remove{, is varied
@@ -118,14 +118,15 @@
 \begin{align} V &\equiv \frac{1}{2} \int ~ \Pi^2 ~ \dxdydz \\ &= \frac{1}{2} \int ~
 N^4 \Omega_z^2 ~ \dxdydz + \int ~ N^2 \Omega_z \bOmega \cdot \bnabla b ~ \dxdydz +
 \frac{1}{2} \int ~ \left( \bOmega \cdot \bnabla b \right)^2 ~ \dxdydz \\ &\equiv
-V_2 + V_3 + V_4 \end{align} \Add{is an invariant of equations (7-9) if no dissipation
-and no forcing. For a flow without vertical vorticity $V_2 = V_3 =0$.}}
+V_2 + V_3 + V_4 \end{align} \Add{is an invariant of equations (7-9) in absence of
+dissipation and forcing. For a flow without vertical vorticity $V_2 = V_3 =0$.}}
 
-We also give the relation between vortical modes and vertical vorticity in the text.
+We also give the relation between vortical modes and vertical vorticity in the
+following text,
 
-\addtoman{\Add{Since the toroidal component $\hatvt$ corresponds only to vertical
+\addtoman{\Add{Since the toroidal component $\hatvt$ corresponds only to the vertical
 vorticity ($\hat{\Omega}_z = i k \hatvt \sin \thk$, with $\bOmega = \bnabla \times \vv$
-the vorticity), we also denote it as the ``vortical" velocity.}}
+being the vorticity), we also denote it as the ``vortical" velocity.}}
 
 We also refer to the study of \cite{waite_potential_2013} in the discussion section:
 
@@ -144,13 +145,13 @@
 diagnostics once this parameter zoo unfolded.}
 
 As already said, we have completely reorganized and strengthened the section on our
-numerical methods. In particular, we tried to better explain how hyper viscosity is
-used and which non-dimensional numbers are important. More generally, we tried to
-simplify the equations by using simpler notations. We also removed the definition of
-the Reynolds number from equation (17). Regarding the introduction of the
-non-dimensional numbers, we now write:
+numerical methods. In particular, we now explain better how hyper-viscosity is used and
+which non-dimensional numbers are important. More generally, we have simplified the
+equations by using simpler notations. We have also removed the definition of the
+Reynolds number from equation (17). Regarding the introduction of the non-dimensional
+numbers, we now write:
 
-\addtoman{\Add{We simulate forced dissipative flows} with the pseudo-spectral solver
+\addtoman{\Add{We simulate forced-dissipated flows} with the pseudo-spectral solver
 \texttt{ns3d.strat} from the FluidSim software [61] (an open-source Python package of
 the FluidDyn project [62] using Fluidfft [63] to compute the Fast Fourier Transforms).
 \Add{The forcing, which will be described in details at the end of this section and in
@@ -185,7 +186,7 @@
 \itemit{You talk of “removing vortical modes”.  How is that done numerically?  Is it a
 similar process as in Holmes-Cerfon et al. in JFM 2013?}
 
-The mentioned paper account for rotation (Coriolis parameter $f \neq 0$) and employ a
+The mentioned paper accounts for rotation (Coriolis parameter $f \neq 0$) and employ a
 damping of geostrophic modes. Here, we project the nonlinear terms in the velocity
 equation on the poloidal manifold, forbidding energy transfer to vortical modes. Shear
 modes are removed in the same way. We modified the text as follows:
@@ -199,8 +200,8 @@
 conversion rate be nonzero a priori?}
 
 The reviewer is right: the \ul{spatially integrated} KE->PE conversion rate is nonzero
-because the energy is injected under the form of kinetic energy. Yet, we refer here to
-the KE->PE conversion spectra, shown in Figure 10. If WWT apply, we expect a range of
+because the energy is injected in the form of kinetic energy. Yet, we refer here to the
+KE->PE conversion spectra, shown in Figure 10. If WWT apply, we expect a range of
 $(k_h,k_z)$ for which the KE->PE conversion is zero in average, even if the spatially
 integrated KE->PE conversion rate is not. We have modify this sentence in order to be
 clarify this point:
@@ -251,9 +252,9 @@
 Diffusion (ID) contribution. Yet, the Parametric Subharmonic Instability (PSI)
 mechanism is responsible for upscale transfer in $\omega$ as explained in section III.D
 (see Figure 15). To our opinion, a complete prediction for energy transfers (based on
-WWT or a more general theory) remains to be find.
+WWT or a more general theory) remains to be found.
 
-Note also that a lot of forcing mechanics known to be important (internal tide,
+Note also that a lot of forcing mechanisms known to be important (internal tide,
 inertial waves, spontaneous generation from slow vortices) force rather slow waves. It
 would be surprising if the $\omega^{-2}$ spectrum down to $\omega = N$ could be
 explained only with an upscale-in-frequency energy cascade.
@@ -269,9 +270,9 @@
 (20a) is relevant because we remove one over two degree of freedom of the horizontal
 velocity when we removing vortical modes. However, this remark make us to realize that
 (20b) is wrong. Indeed, the dissipation still occurs in the three spatial directions if
-we remove vortical modes, so definition (19) does not need to adapted. We have
-suppressed (20b) and corresponding text, and corrected Figure 2. It changes $I_{\rm
-diss}$ by a factor 4/3 in Figure 2~b, so the conclusions remain the same.
+we remove vortical modes, so definition (19) does not need to be adapted. We have
+deleted (20b) and corresponding text, and corrected Figure 2. It changes $I_{\rm diss}$
+by factor 4/3 in Figure 2~b, so the conclusions remain the same.
 
 
 \itemit{Figure 4: I find it surprising that there is no “hump” in the spectrum at the