# HG changeset patch
# User vlabarre <vincent.labarre@oca.eu>
# Date 1698091604 -7200
#      Mon Oct 23 22:06:44 2023 +0200
# Node ID 2f41535b33658b55690c6bc6306099870d7caa35
# Parent  b58af7713c85719608adaf3efa2b90c2c0347b7e
Questions of forcing + addtoman + cover letter

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
@@ -484,8 +484,8 @@
 
 \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}),
+\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
@@ -597,10 +597,11 @@
 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
+$\sin \theta_f = 0.3$, meaning that relatively slow internal waves are forced. \Remove{The
 correlation time of the forcing is equal to the period of the forced waves $T_c = 2\pi
-/(N \sin \theta_f)$. The forcing scheme\Add{, which does not involve any stochastic
-process,} is described in Appendix~\ref{appendix:forcing}. \Add{Forcing slow waves is
+/(N \sin \theta_f)$.} The forcing scheme is described in Appendix~\ref{appendix:forcing}. \Add{It is not harmonic, and not given by a stochastic differential equation either. Instead, the phase and
+the amplitude of the forced wave numbers are randomly changed every renewal time $T_c$ equal to the period of the forced waves $2\pi
+/(N \sin \theta_f)$.} \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
@@ -968,11 +969,11 @@
 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
-stacked shear layers (Figure~\ref{fig:spectra_khkz}$\rm (a)$)}. When vortical modes are
+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.
+show the same trends. 
 
 \begin{figure}
 \includegraphics[width=1.0\textwidth]{figure8}
@@ -1469,7 +1470,7 @@
 \section{Forcing scheme}
 \label{appendix:forcing}
 
-The forcing used in our simulations is prescribed by the following algorithm:
+\Add{Our forcing is designed to excite waves. Phase and amplitude of the forced wave numbers are randomly changed in time every $T_c$ such that the input kinetic energy is equal to unity. Namely, we use the following algorithm:}
 
 \begin{algorithm}[H]
 $t=0$: generate two random numbers $f_{0\kk}, f_{1\kk} \sim \mathcal{P}$ $\forall \kk$ in the forcing region; $t_0 = 0$; \\
@@ -1481,7 +1482,7 @@
 		genererate $f_{1\kk} \sim \mathcal{P}$
 	}{}
 	$\hatff = \left\{f_{0\kk} - \dfrac{(f_{1\kk} - f_{0\kk})}{2} \left[\cos\left(\dfrac{\pi(t-t_0)}{T_c} \right)+ 1 \right]\right\}~ \eep$ \\
-	$\hatff$ are normalized to ensure $P_K = 1$
+	$\hatff$ are normalized to ensure \Add{$P_K(t) = \sum\limits_{\kk} ~ \Re \left[ \hatff \cdot \hatvv^* \right] =  1$}
 }
 %\caption{}
 \end{algorithm}
@@ -1489,11 +1490,11 @@
 where $\mathcal{P}$ is the probability law which has for distribution
 \begin{equation}
 p_X(x) = \begin{cases}
-\frac{1}{4} ~~ \text{if} ~ Re(x) \in[-1:1]  ~ \text{and} ~ Im(x) \in[-1:1] \\
+\frac{1}{4} ~~ \text{if} ~ \Add{\Re(x)} \in[-1:1]  ~ \text{and} ~ \Add{\Im(x)} \in[-1:1] \\
 0 ~~ \text{otherwise}
 \end{cases}
 \end{equation}
-and $\Delta t$ is the time increment at each time step.
+\Remove{and} $\Delta t$ is the time increment at each time step \Add{ and $T$ is the final time of the simulation}.
 
 
 \section{List of simulations}
@@ -1502,7 +1503,7 @@
 \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{5cm}m{4cm}m{4cm}}
+\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}
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
@@ -37,6 +37,34 @@
 
 \maketitle
 
-Hello!
+Dear Dr. Eckart Meiburg,
+
+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.
+
+We have corrected the manuscript by improving the presentation of our methods. The
+corresponding section has been reorganized and strengthened. In particular, we provide
+more precision about the forcing mechanism, how we use hyperviscosity and
+hyperdiffusion, and how we varied the vertical resolution with respect to the
+Brunt-Väisälä frequency. We have also clarified the manuscript to highlight the two
+control parameters of this study: the horizontal Froude number and the buoyancy
+Reynolds number.
+
+Following the reviewers' suggestions, we have added visualizations of the flows
+considered here and provided comments on the role of potential vorticity and potential
+enstrophy. We have also found additional studies that help to fill the parameter space
+on the last figure of the manuscript, and we have commented on them in the discussion
+section.
+
+Finally, we have corrected typos and small errors.
+
+We look forward to hearing from you.
+
+Best regards,
+
+Vincent Labarre, Pierre Augier, Giorgio Krstulovic, and Sergey Nazarenko
+
 
 \end{document}
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
@@ -48,7 +48,8 @@
 \newcommand{\Remove}[1]{{\color{red}\st{#1}}}
 \newcommand{\Comment}[1]{{\color{green}#1}}
 
-\newcommand{\itemit}[1]{\item{\it #1}}
+\newcommand{\itemit}[1]{\item{\bf #1}}
+\newcommand{\addtoman}[1]{\hspace{1cm} \begin{minipage}{14cm} ``{\it #1}'' \end{minipage}}
 
 
 \begin{document}
@@ -75,12 +76,12 @@
 Instead of the exact 2/3 rule for cubic truncation, we use a 0.8 spherical truncation.
 We performed with Jason Reneuve an extensive study on dealiasing methods (also
 including phase shifting). It is not yet published but a notable result is that for
-simulations with $\eta \kmax$ slightly smaller than one, we obtain better results to
+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
 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 $\eta \kmax$ close to one, more wavenumbers to represent the dissipative
+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
 \%$ of the grid points are kept!). On the other hand, the coefficient needs to be
@@ -96,9 +97,9 @@
 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:
 
-``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.}''
+\addtoman{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.}}
 
 \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
@@ -113,20 +114,19 @@
 \cite{nikurashin_legg_mechanism_2011}. In this study, the Brünt-Väisälä frequency is
 typically $N=10^{-3}$ rad/s to represent the main thermocline, while the barotropic
 $M_2$ tide is forced at $\omega = 1.4\times 10^{-4}$ rad/s 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 now motivate
-the low frequency forcing in the new version of the manuscript and included:
+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.
 
-``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 correlation time of the forcing is equal to the period of the forced
-waves $T_c = 2\pi /(N \sin \theta_f)$. \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
+We now motivate the low frequency forcing in the new version of the manuscript and
+included:
+
+\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
+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$.}''
+spectra close to $N$.}}
 
 
 Yet, the referee raise an important point: forcing waves at higher frequencies would
@@ -148,32 +148,31 @@
 
 We then present how we use hyperdiffusivity in some (but not all) 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 (8) and (9),
-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~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 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.}''
+\addtoman{\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 (8)
+and (9), 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$ [31,
+64].} \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~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
+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.}}
 
 We can show here a re-plotting of Figure 3 of our manuscript while keeping only DNS
 ($\kmax\eta \geq 1$):
@@ -189,7 +188,7 @@
 
 Moreover, we decided to discard for this article few simulations 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 $\eta$ and in the diffusive range linked to
+under-resolved, with $\kmax$ close to $1/\eta$ and in the diffusive range linked to
 standard dissipation.
 
 Regarding your question about the importance of hyperviscosity and the potential ``blow
@@ -207,10 +206,9 @@
 energy transfer to shear modes is forbidden. We change a sentence in the method section
 to make this point clearer:
 
-``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}.''
-
+\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
+simulations \Add{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
@@ -227,34 +225,38 @@
 
 Figure~\ref{fig:forcing} shows how the forcing depends in time. We plan to include this
 figure and a discussion about this time dependence in a future article studying
-specifically the first dataset without projection.
+specifically the dataset composed of simulations with 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
+randomly changed in time every $T_c$ such that the kinetic energy injection rate is
+equal to one. However, even if we think that your question about the effect of a purely
+harmonic forcing deserves additional studies, we don't think that we need to have a
+purely harmonic forcing to force waves. The forcing injects or removes energy to the
+waves with a positive average. In our opinion it is nicer if the forcing is time
+correlated and do not force wave modes at $\omega > \omega_k$ (which is why we use a
+time correlated forcing), but anyway the system responds with its own dynamics to a
+forcing, and as long as poloidal velocity (or buoyancy) are forced, internal waves are
+forced. Note that for large $N$, the forcing is weak compared to linear terms. As a
+result, we clearly see in the spatio-temporal poloidal spectra that the waves ($\omega
+\simeq \omega_k$) strongly dominate the forced wavenumbers.
 
-As described in Appendix A, there is no stochastic process involved (phase and
-amplitude of the forced wavenumbers are randomly changed in time) but the forcing is
-not harmonic and therefore does not exactly follow the wave behavior. However, even if
-we think that your question about the effect of a purely harmonic forcing deserves
-additional studies, we don't think that we need to have a purely harmonic forcing to
-force waves. The forcing injects or removes energy to the waves with a positive
-average. In our opinion it is nicer if the forcing is time correlated and do not force
-wave modes at $\omega > \omega_k$ (which is why we use a time correlated forcing), but
-anyway the system responds with its own dynamics to a forcing, and as long as poloidal
-velocity (or buoyancy) are forced, internal waves are forced. Note that for large $N$,
-the forcing is weak compared to linear terms. As a result, we clearly see in the
-spatio-temporal poloidal spectra that the waves ($\omega \simeq \omega_k$) strongly
-dominate the forced wavenumbers.
+We give more details about the forcing in the method section:
 
-We now explicitly mention in the text that forcing scheme (described in Appendix~A)
-does not involve any stochastic process:
+\addtoman{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. \Remove{The correlation time of the forcing is equal to the period of the
+forced waves $T_c = 2\pi /(N \sin \theta_f)$.} The forcing scheme is described in
+Appendix~A. \Add{It is not harmonic, and not given by a stochastic differential
+equation either. Instead, the phase and the amplitude of the forced wave numbers are
+randomly changed every renewal time $T_c$ equal to the period of the forced waves $2\pi
+/(N \sin \theta_f)$.}}
 
-``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
-correlation time of the forcing is equal to the period of the forced waves $T_c = 2\pi
-/(N \sin \theta_f)$. The forcing scheme\Add{, which does not involve any stochastic
-process,} is described in Appendix~A.''
+and in Appendix~A.
 
 
 \itemit{Figure 7: could you show the spatial structure of the dominant vortical mode by
@@ -266,13 +268,20 @@
 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}
+\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. In our opinion, the last figure is not
-useful enough to be included in the manuscript.
+\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)$).
+Energy then tends to be accumulated at the smallest horizontal wave vectors, close to
+shear modes.}}
+
+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
@@ -284,16 +293,14 @@
 \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
+\addtoman{\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 [22, 30, 31, 66]. 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.''}.
-
+larger amplitude.}}
 
 \itemit{Figure 9: what is the purpose of the kz=kh line?}
 
@@ -330,17 +337,16 @@
 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.''}
+\addtoman{\Add{We observe that [74] 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 [42]. Yet, the study of the
+wave energy ratio is not done while keeping $F_h$ or $\R$ constant in [42]. 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
@@ -363,12 +369,12 @@
 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 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}.''}
+\addtoman{\Add{It is worth mentioning that [46] 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
+[49].}}
 
 \end{enumerate}
 
@@ -387,17 +393,19 @@
 \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 + (\mathrm{d}\bar{\rho} /
-\mathrm{d}z) z$.}
+We add this definition in the method section:
+
+\addtoman{\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 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.
+We replace the notation for the real part by $\Re$ and for the imaginary part by $\Im$
+to avoid ambiguous notation.
 
 \end{itemize}
 
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
@@ -52,7 +52,9 @@
 \newcommand{\Remove}[1]{{\color{red}\st{#1}}}
 \newcommand{\Comment}[1]{{\color{green}#1}}
 
-\newcommand{\itemit}[1]{\item{\it #1}}
+
+\newcommand{\itemit}[1]{\item{\bf #1}}
+\newcommand{\addtoman}[1]{\hspace{1cm} \begin{minipage}{14cm} ``{\it #1}'' \end{minipage}}
 
 \begin{document}
 \maketitle
@@ -79,19 +81,18 @@
 
 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.}''
-
+\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
+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 [27]. 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
@@ -109,36 +110,31 @@
 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
+\addtoman{\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$.}''
+\right)$ is a Lagrangian invariant [69]. 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.}''
+\addtoman{\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}''}
+\addtoman{\Add{We observe that [74] 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 [42]}}
 
 \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
@@ -154,7 +150,7 @@
 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
+\addtoman{\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
@@ -164,14 +160,11 @@
 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$.}
+number \Remove{, the Reynolds number,} and the buoyancy Reynolds number [31] 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$),
@@ -198,10 +191,9 @@
 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}.''
-
+\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
+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
@@ -214,10 +206,10 @@
 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.''
+\addtoman{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}
 
@@ -245,18 +237,18 @@
 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}$.}''.
+after equation (3) by adding \addtoman{\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
+example, the direction of the energy transfer 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
+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.
 
@@ -276,7 +268,7 @@
 (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
+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.