Skip to content
Snippets Groups Projects
rebut1.tex 21.3 KiB
Newer Older
\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 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{\kk}{\boldsymbol{k}}
\newcommand{\bnabla}{\boldsymbol{\nabla}}
\newcommand{\bomega}{\boldsymbol{\omega}}
\newcommand{\ok}{\omega_{\kk}}
\newcommand{\bv}{Brunt-V\"ais\"al\"a }
Pierre Augier's avatar
Pierre Augier committed
\newcommand{\eep}{\boldsymbol{e}_{p\kk}}

\newcommand{\hatf}{\hat{f}}
\newcommand{\hatff}{\hat{\boldsymbol{f}}}

\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{\bf #1}}
\newcommand{\addtoman}[1]{\hspace{1cm} \begin{minipage}{14cm} ``{\it #1}'' \end{minipage}}


\begin{document}

\maketitle

\noindent We thank the referee for his/her critical comments, that helped us to
significantly improve our 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.

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?}

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 $\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 $\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
smaller than 0.94 to kill the double aliases, which are the most annoying ones. A value
of 0.8 seems to be a good compromize for which most of the aliases are 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 aliases. Note also than 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
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:

\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
required by WWT?}

Forcing low frequency waves is mainly motivated by geophysical applications: in the
oceans, tidal flows over topography generate low frequency waves (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}$ 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:

\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$.}}


Yet, the referee raise an important point: forcing waves at higher frequencies would
help 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?}

Pierre Augier's avatar
Pierre Augier committed
As already mentioned, the presentation of our methods was a strong weakness of 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
standard incompressible Navier-Stokes equations (without hyperdiffusivity) and
introduce our control and output parameters and non dimensional numbers.

We then present how we use hyperdiffusivity in some (but not all) simulations:
\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 show here a re-plotting of Figure 3 of our manuscript while keeping only DNS
($\kmax\eta \geq 1$):

\begin{figure}[h]
\centering \includegraphics[width=1\linewidth]{Figure3_DNS}
\caption{Same as Figure 3 but keeping only proper DNS.}
Comparing with Figure 3 of the manuscript first show that only a small number of
simulations are not DNS. Secondly, keeping only DNS does not changes the results
presented in the present study.
Pierre Augier's avatar
Pierre Augier committed
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 $1/\eta$ and in the diffusive range linked to
Pierre Augier's avatar
Pierre Augier committed
standard dissipation.

Regarding your question about the importance of hyperviscosity and the potential ``blow
up without it'', we can mention that, except for very small $\kmax\eta$ values,
Pierre Augier's avatar
Pierre Augier committed
removing hyperdiffusion would not lead to simulation blow up but we would observe a
tendency to thermalisation, with accumulation of energy at the larger wavenumbers.

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

\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
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?}

Pierre Augier's avatar
Pierre Augier committed
\begin{figure}[h]
\centering \includegraphics[width=0.8\linewidth]{../../tmp/forcing_vs_time.png}
\caption{Time dependence of the forcing.} \label{fig:forcing}
\end{figure}
Pierre Augier's avatar
Pierre Augier committed
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 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 can
be 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:
Pierre Augier's avatar
Pierre Augier committed

\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. 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.}}
Pierre Augier's avatar
Pierre Augier committed

and in Appendix~A, which has also been improved.


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

Pierre Augier's avatar
Pierre Augier committed
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=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:

\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
Pierre Augier's avatar
Pierre Augier committed
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 added the following text to describe this figure:
\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.}}
\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 $\R \gtrsim 1$ (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
Pierre Augier's avatar
Pierre Augier committed
\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
Pierre Augier's avatar
Pierre Augier committed
\cite{waite_stratified_2004} and \cite{waite_stratified_2006} tables. In
\cite{lindborg_energy_2006} and \cite{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:

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

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

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

\subsection*{Additional small comments:}

\begin{itemize}

\itemit{Introduction, line 2: typo "aims to provides"}
\itemit{p.2 second paragraph: typo "difficultly"}
\itemit{p.4, section II, first sentence: typo "hypervisity", define buoyancy, what does
"kinematic pressure" mean?}
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$ and for the imaginary part by $\Im$
to avoid ambiguous notation.
\end{itemize}

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

\end{document}