\documentclass[a4paper,12pt]{article}

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

\usepackage{natbib}
\bibliographystyle{abbrv}

\usepackage[english]{babel}

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

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

\usepackage{array}

\usepackage{soul}

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

\newcommand{\diff}{\text{d}}
\newcommand{\R}{\mathcal{R}}
\newcommand{\kmax}{k_{\rm max}}
\newcommand{\eez}{\boldsymbol{e}_z}
\newcommand{\vv}{\boldsymbol{v}}
\newcommand{\bnabla}{\boldsymbol{\nabla}}
\newcommand{\bomega}{\boldsymbol{\omega}}

\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 We thank the referee for his/her critical comments, that helped us to
significantly improve our manuscript. The answer on all the comments are listed below.
Corresponding corrections are made in blue in the new draft.

\begin{enumerate}

\itemit{p.5: what is the origin of the 0.8 factor instead of the more classical 2/3
rule for dealiasing? Is it related to the use of hyperviscosity and hyperdiffusion?}

We performed with Jason Reneuve an extensive study on dealiasing methods (also with
phase shifting). It is not yet published but a notable result is that ... TODO...


\itemit{p.5: what is the motivation behind forcing low frequency waves only? Wouldn't
it be better to increase the wave frequency to better satisfy the timescale separation
required by WWT?}

Forcing low frequency is mainly motivated by geophysical applications: In oceans, waves
are generated by tidal flow over topography leading to internal tides
\cite{mackinnon_climate_2017}. This forcing is particularly important when compared to
other mechanisms, and motivate many studies on stratified flows. For example, radiation
and dissipation of internal waves generated by the barotropic $M_2$ (semidiurnal) tide
over rough topography are explored in \cite{nikurashin_legg_mechanism_2011}. In this
study, the Brünt-Väisälä frequency is typically $N=10^{-3} ~ s^{-1}$ to represent the
main thermocline, while the barotropic $M_2$ tide is forced at $\omega = 1.4\times
10^{-4}~s^{-1}$ by adding a body force to momentum equations. 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 motivate the low frequency forcing in the new
version of the manuscript. Yet, the referee raise in important point: It would be
preferable to force waves at higher frequencies in order to satisfy the timescale
separation required by WWT. In our opinion, this would require additional studies.

\itemit{Why is hyperviscosity added to the regular viscosity? Is it just for
regularisation purposes, but in that case, the simulations presented cannot be labelled
as DNS I assume? How important it is in the simulations presented here, would they all
blow up without it?}

Each simulation of these datasets is constructed in the following way:
\begin{enumerate} \item We start from a simulation at small resolution to simulate the transient state.
Generally, hyper-viscosity is required for regularization purposes. \item Once a
statistically steady-state is reached, we increase the resolution of the simulation
while decreasing hyper-viscosity. Then the simulation is runned until reaching a new
statistically steady state. \item The previous step is repeated until reaching highest
resolutions.\end{enumerate}
This strategy is employed because reaching a steady state in strongly stratified flows
requires prohibitively long simulations if starting at high resolution from the initial
time. We plan to explain this methodology in more details in a study focusing on the
dataset without vortical modes.

When the highest resolutions are reached (see table in Appendix B), most of the
simulations are proper DNS with $\kmax\eta$ larger or close to one. For example, for
the simulations $(N, \R_i) (40, 20)$ (which are widely used in the present study)
$\kmax\eta = 0.99$ (with vortical modes) and $1.05$ (without vortical modes). For such
simulations, a blow up is not expected if hyper-viscosity is removed.

Few simulations, with the largest $N$ or $R_i$, are not DNS. However, they can be to
some extend characterized with the standard Reynolds number. These simulations could
blow up if hyper-viscosity is removed (in particular if $\kmax\eta \ll 1$). Yet, this
does not change the conclusions of our study. As an example, re-plotting Figure 3 of
our manuscript while keeping only DNS ($\kmax\eta \geq 1$) gives the following:

\begin{figure}[h]
\centering \includegraphics[width=1\linewidth]{Figure3_DNS}
\caption{Same as figure 3 but keeping only proper DNS.}
\end{figure}

First, we observe that only a small number of simulations are not DNS. Secondly,
keeping only DNS does not changes the results presented in the present study. Note that
it is always possible to define an effective Reynolds number based on hyper-viscosity,
and to define the buoyancy Reynolds number as $Re F_h^{-2}$ (where $F_h$ is the Froude
number) \cite{brethouwer_scaling_2007}. \\  We plan to share these datasets, that
could be improved by increasing the resolution (and decreasing the hyper-viscosity) on
non-DNS runs. Yet, as the reviewer's remark suggests, these simulations should be
regarded with caution, in particular when looking at small scales or flows' statistics.
For this reason, we added in the method section the sentence:

\Add{``Simulations with $\kmax \eta \lesssim 1$ remain affected by hyper-viscosity. In
that case, small-scales and flow statistics should be analyzed carefuly since
dissipation is not only due to usual viscosity like in a direct numerical simulation
(typically when $\kmax \eta > 1$). Simulations that do not satisfy the criteria $\kmax
\eta > 1$ are few in our dataset, and we checked that they do not change the results
presented here.''}

\itemit{Shear modes are removed from all simulations, but how is this achieved exactly?
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:

``\Add{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 \ul{by fixing nonlinear transfers to these modes to zero}}.''

\itemit{It is mentionned that the forcing is not delta-correlated in time but that its
correlation time is consistent with the dispersion relation of internal gravity waves
at that particular angle $\theta_k$. Can the authors give more details about the
forcing scheme? In particular, if there is still a stochastic process involved, I guess
it means that not all of the energy is injected in the form of waves. Do the authors
think that using a purely harmonic forcing, thereby injecting all of the energy in the
form of internal gravity waves, would change something to the results presented here?}

TODO: describe better the forcing in the article + add added text here.

- no proper stochastic process involve. Phase and amplitude of the forced wavenumbers
random and changed in time ...

- a bit of Python to compute time spectra of such process.

- this aspect of the forcing will be discussed in a paper describing more globally the
dataset.

TODO: add the question on possible effect of harmonic forcing. We don't know.
Interesting open question that could be investigated in future studies.


\itemit{Figure 7: could you show the spatial structure of the dominant vortical mode by
filtering the coefficients in Fourier space and transform it back to physical space?
More generally, some visualisations of the typical flows considered, with and without
vortical modes, would be helpful in my opinion.}

Here is a 3D figure showing the toroidal velocity of the dominant vortical mode The
figure of the dominant vortical mode is from the simulation at resolution $n_h = 320$.
\begin{figure}[h]
\centering
\includegraphics[width=1\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. In our opinion, the last figure is not
useful enough to be included in the manuscript.

Yet, following reviewer's suggestion, we include the following figure (now figure 4)
showing the buoyancy field for simulations with $(N,\R_i) = (40,20)$, which are studied
in detail in our manuscript:
\begin{figure}[h]
\centering \includegraphics[width=1\linewidth]{../../input/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}
We also add the following text to describe this figure:

\Add{``Figure~4 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~4(a)) so the buoyancy has a different
structure than when vortical modes are absent (Figure~4(b)). Without vortical modes the
dynamics in the horizontal direction is irrotational and the buoyancy field has a
larger amplitude.''}.

\itemit{Figure 9: what is the purpose of the kz=kh line?}

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.

\itemit{Figure 16: this is an important plot in my opinion, it might be useful to
include other studies? In particular, it is rather suprising that [46] is so far from
the dynamical regime considered in this paper (knowing that both studies were looking
for traces of WWT in simulations). Are they no other studies filling the gap between
these two extremes?}

For studies motivated by geophysical flows, the focus is more on large buoyancy
Reynolds $\R$ number, or at least not very small. This is why most of recent
simulations are for $1 \lesssim \R$ (see e.g. \cite{kimura_energy_2012,
bartello_sensitivity_2013, maffioli_vertical_2017}). Some older studies could have fill
the gap between our study and the one of \cite{reun_parametric_2018}, possibly because
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.
\cite{waite_stratified_2004} forced vortical modes so these it is not well suited for
observing internal gravity wave turbulence. We don't see how to extract $\R$ from
\cite{waite_stratified_2004, waite_stratified_2006} tables. In
\cite{lindborg_energy_2006, waite_stratified_2011}, different values of the viscosity
where used on the vertical and horizontal so these studies are more difficult to
compare to real flows.

We are aware of three studies that fill the gap between our study and
\cite{reun_parametric_2018}. The first one is \cite{brethouwer_scaling_2007}, the
second one is \cite{waite_potential_2013}, and the third one to \cite{lam_energy_2021},
so we include them in figure 16 (now figure 17). We also added the following remarks in
the discussion:

\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}. 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 a decreasing function of $\R
< 1$ at a constant $F_h \ll 1$, and if a threshold below which both shear and vortical
modes are stable exists.''}

\itemit{The shear modes are always suppressed while the vortical modes are either left
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
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
for large $\R$.


\itemit{Note that Le Reun et al. [46] did not have to remove neither shear nor vortical
modes (while the same authors had to remove the geostrophic flow in rotating systems to
observe regimes reminiscent of inertial wave turbulence), which seems like an important
information to recall in the present study interested in the role of such slow flows on
the dynamics. Similarly to the rotating case and the spontaneous emergence of
geostrophic flows, there must be some threshold below which both shear and vortical
modes are stable and won't grow if not directly forced?}

This is indeed an important information that must be mentioned. Consequently, we add
the following sentences in the discussion:

\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 are stable and won't grow if not directly forced. Such a situation would
be analogous to rotating flows where a threshold below which geostrophic modes are
stable exists \cite{reun_experimental_2019}.''}

\end{enumerate}

\subsection*{Additional small comments:}

\begin{itemize}

\itemit{Introduction, line 2: typo "aims to provides"}

Corrected.

\itemit{p.2 second paragraph: typo "difficultly"}

Corrected.

\itemit{p.4, section II, first sentence: typo "hypervisity", define buoyancy, what does
"kinematic pressure" mean?}

We add this definition in the method section: \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 +
\frac{\mathrm{d}\bar{\rho}}{\mathrm{d}z}z$.}

The kinematic pressure is the pressure divided by fluid density.

\itemit{Ambiguous notation between real part and Reynolds}

We replace the notation for the real part by $\Re$ to avoid ambiguous notation.

\end{itemize}

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

\end{document}
