\documentclass[a4paper,12pt]{article}

\usepackage{lmodern}
\usepackage{textcomp}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}

\usepackage[english]{babel}

\usepackage{amssymb, amsmath}
\usepackage[top=0.8in,bottom=0.8in,left=2cm,right=2cm,a4paper]{geometry}
\usepackage{graphicx}

\usepackage{natbib}
\bibliographystyle{unsrtnat}

\usepackage[usenames,dvipsnames]{color}
\usepackage[normalem]{ulem}

\usepackage{array}

\usepackage{soul}

\title{Reply to referee 2}
\author{}
\date{}

\newcommand{\diff}{\text{d}}
\newcommand{\dxdydz}{\diff x \diff y \diff z}
\newcommand{\R}{\mathcal{R}}
\newcommand{\kmax}{k_{\rm max}}
\newcommand{\eez}{\boldsymbol{e}_z}
\newcommand{\vv}{\boldsymbol{v}}
\newcommand{\kk}{\boldsymbol{k}}
\newcommand{\bnabla}{\boldsymbol{\nabla}}
\newcommand{\bomega}{\boldsymbol{\omega}}
\newcommand{\bOmega}{\boldsymbol{\Omega}}
\newcommand{\ok}{\omega_{\kk}}
\newcommand{\bv}{Brunt-V\"ais\"al\"a }
\newcommand{\eep}{\boldsymbol{e}_{p\kk}}
\newcommand{\thk}{\theta_{\kk}}

\newcommand{\hatf}{\hat{f}}
\newcommand{\hatff}{\hat{\boldsymbol{f}}}
\newcommand{\hatvt}{\hat{v}_{t}}

\newcommand{\epsK}{{\varepsilon_{\!\scriptscriptstyle K}}}
\newcommand{\epsKK}{{\varepsilon_{\!\scriptscriptstyle K 2}}}
\newcommand{\epsKKKK}{{\varepsilon_{\!\scriptscriptstyle K 4}}}

\newcommand{\Add}[1]{{\color{blue}#1}}
\newcommand{\Remove}[1]{{\color{red}\st{#1}}}
\newcommand{\Comment}[1]{{\color{green}#1}}

\newcommand{\itemit}[1]{\item{\it #1}}

\begin{document}
\maketitle

\noindent The authors thank the referee for his/her useful and constructive comments
that helped us to improve the manuscript. We realized that our Methods section was
clearly not 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
reorganized and strengthened.

Below we have responded to all questions. Corresponding corrections are made in blue in
the new draft.

\subsection*{Major comments}

\begin{enumerate}

\itemit{The methods section is vague about the choice of vertical domain size and
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:

``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 Brunt-V\"ais\"al\"a 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.}''


\itemit{The poloidal forcing is meant to excite only waves but you do not discuss its
impact on the potential vorticity (PV), which governs the vortical mode.  In fact I was
surprised that you do not discuss PV in your paper more generally, given how important
a concept it is in geophysical fluid mechanics.}

PV is known to be much less dynamically important in stratified flows than in
stratified and strongly rotating flows. The Ertel PV $\bOmega \cdot \bnabla (N^2 z + b)
= \Omega_z (N^2 + \partial_z b) + \bOmega_h \cdot \bnabla_h b$ is not quadratic and is
not equal to the Charney PV. There is no PV cascade like in quasi-geostrophic flows.
%
The forcing in poloidal velocity produces and removes Ertel PV but we do not know what
can be deduced from that.
%
We add some facts concerning Potential Vorticity (PV) and Potential Enstrophy (PE) in
the method section:

``\Add{Without dissipation and forcing, equations (7-9) 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 (7-9) if no dissipation and no 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.

``\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.}''

We also refer to the study of \cite{waite_potential_2013} in the discussion section:

\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 ~ \diff x
\diff y \diff z$ 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}''}

\itemit{I felt overwhelmed by the great number of non-dimensional parameters introduced
here; even with all these definitions one could still invent more, ie parameters based
on the hyperviscosities.  Basic physical insight is lost for me here, in particular as
I would like to see parameters based on controllable parameters such as the energy
input. It's a personal choice/failing but I lost interest in the detailed numerical
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:

\Add{We simulate forced dissipative 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
Appendix~A, 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} ~~~~ \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 hope that it is clearer that the simulations are characterized by two input
parameters ($N,\,\R_i$) and two output turbulent non-dimensional numbers ($F_h,\,\R$),
which are commonly used in other studies.

Only few of our simulations remain (weakly) affected by hyperviscosity. They correspond
to $0.45< \kmax \eta \lesssim 1$. We checked that these simulations don't change the
results presented in this study (see answer to reviewer 1) so the relevant
dimensionless numbers are only the horizontal Froude and the buoyancy Reynolds numbers
defined in equation (17).

We introduce many physical quantities in equations (21-32). This constitutes a
technical part that is rather difficult to read. Yet, all the presented quantities are
useful (and used in the manuscript) to discuss the spectral energy budget. In our
opinion, removing some of this diagnostic would make the manuscript unclear and less
precise.


\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
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:

``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}.''


\itemit{In the discussion session you say that the energy conversion rate does not
fluctuate around zero.  As only the kinetic energy is forced, shouldn't the KE->PE
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
$(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:

``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.''

\end{enumerate}

\subsection*{Minor comments}

\begin{enumerate}

\itemit{(1) why introduce a symbol $k_b$ here, which is non-standard anyway?  This
conditions is $k U_h/N << 1$, which is a standard criterion for horizontal phase
speeds, not scale separation as stated here.}

The buoyancy scale $L_{\rm b} = U_h/N$ and its inverse, the buoyancy wave vector
$k_{\rm b} = N/U_h$, are well established concepts in the study of stratified
turbulence (see e.g.
\cite{billant_self-similarity_2001,kimura_energy_2012,maffioli_vertical_2017,yokoyama_energy-based_2019}).
$L_{\rm b}$ corresponds physically to the typical vertical displacement of fluid
particles with respect to their floatability level, which is also the size of the
layers observed in stratified turbulent flows \cite{waite_stratified_2011}.
Consequently, the condition $k U_h/N = k / k_{\rm b}\ll 1$ can be seen as a
length-scale separation.


\itemit{(3) maybe it could be stated how these are derived from each other?  Also, the
folklore in oceanography is that the internal wave energy flows downscale in wavenumber
but upscale in frequency. Do the quoted WWT theories agree with that?}

We explain how the change of coordinates $(k_h, k_z) \rightarrow (\omega, k_z)$ is done
after equation (3) by adding ``\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}$.}''.

The question of energy transfers in anisotropic wave systems is a subtle question. For
example, the direction of the energy cascade depends on the considered triad for
inertial wave turbulence \cite[]{david_locality_2023}. Generally speaking, the
transfers also depend on the exponents of the energy spectrum. In the case of
stratified flows, \cite{dematteis_origins_2022} argued that the global energy flux is
direct in both $k_z$ and $\omega$ for some spectra due to the dominance of Induced
Diffusion (ID) contribution. Yet, the Parametric Subharmonic Instability (PSI)
mechanism is responsible for inverse cascade 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.

Note also that a lot of forcing mechanics 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 a upscale in frequency energy cascade.


\itemit{typo p4 hyper viscosity}

This typo is corrected.


\itemit{(20) is an interesting definition. But I don't understand (20b)}

(20a) is relevant because we remove one over two degree of freedom of the horizontal
velocity when we removing vortical modes. However, the reviewer's remark make us
realize that (20b) is wrong. Indeed, the dissipation still occurs in the three spatial
directions 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.


\itemit{Figure 4: I find it surprising that there is no “hump” in the spectrum at the
injection frequencies, as is typical for turbulent spectra.  Why is that?}

The energy tends to be accumulated at small $k_h/k$ (see 2-dimensional spectra on
Figure 8, and Figure 14), corresponding to slow internal gravity waves and/or vortical
modes. It follows that small frequencies have more energy than forced frequencies such
that no ``hump'' is observed on the spectra shown in Figure 4 (now Figure 5).


\itemit{The forcing renewal time $T_c$ is chosen as a wave period; why is that? Why not
force with white noise in time?}

We chose $T_c$ to be equal to the wave period to excite frequencies of the same order
than the linear internal gravity $\ok = N k_h/k$. In contrast, a white noise would
excite all frequencies, favoring the appearance of vortical modes and non linear
structures.

\end{enumerate}

\bibliography{../../input/main}

\end{document}
