%#! latex

\documentclass{article}
%\documentclass{jfm}

%\usepackage{bm}
\usepackage{graphicx}
\usepackage{tabularx}
\usepackage[sort&compress]{natbib}
\usepackage{amssymb}
%\usepackage{upmath}
%\usepackage{amsmath}
%\usepackage{framed}

\usepackage{fullpage}

%\addtowidth{\leftmargin}{-0.5in}
%\addtowidth{\pagewidth}{1in}

\providecommand\boldsymbol[1]{\mbox{\boldmath $#1$}}
\newcommand\bldp{\boldsymbol{p}}
\newcommand\bldk{\boldsymbol{k}}
\newcommand\bldx{\boldsymbol{x}}
\providecommand\bnabla{\boldsymbol{\nabla}}

\providecommand\upi{\pi}
\providecommand\upartial{\partial}

\title{%
Energy spectra of internal waves
in stratified fluid:
Direct numerical simulations
}
\author{
Yuri V. Lvov and Naoto Yokoyama\\
{\small  Department of Mathematical Sciences,
 Rensselaer Polytechnic Institute} \\{\small 
 Troy, New York 12180 USA}}


\begin{document}

\maketitle

\begin{abstract}
To investigate formation mechanism of energy spectra,
direct numerical simulations of oceanic internal waves
based on a dynamic equation
are performed.
Our simulations demonstrate that the resulting energy spectra in inertial regions
are not universal.
Instead,
energy accumulates in the horizontally longest waves.
We provide evidence that
formation of energy spectra in inertial regions
is dominated by scale-separated interactions with this accumulation at the longest waves.
These nonlocal interactions in the wavenumber space make it questionable
that a universal spectrum might be observed in direct numerical simulations.
Implications of these findings to the real
oceans and universal oceanographic spectrum are discussed.
\end{abstract}




\section{Introduction}

Oceanic internal waves,
whose restoring force is buoyancy in stratified fluid,
exist widely in oceans.
Internal waves are excited by topographies, tides and atmospheric disturbances.
Their energy is transferred by nonlinear interactions among wavenumbers
from large scales to small scales,
and is dissipated
in the small spatial scales of the Navier--Stokes turbulence generated by breaking of internal waves.
Internal waves play a significant role in the general circulation of oceans
and hence the climate of the Earth.

Energy spectra of internal waves are very broad with
horizontal wavelengths varying from $10$ meters to $10^5$ meters,
vertical wavelengths from $10$ meters to $10^3$ meters,
and time periods from $10^3$ seconds to $10^5$ seconds.
Surprisingly and despite all complexities,
energy spectra of internal waves in the oceans appear to be somewhat universal,
and are given by the Garrett--Munk (GM) spectrum.
The GM spectrum was initially observed in the North Atlantic Ocean
\citep{GM72,GM75,GM_ARF}.

In this paper
we address the question of
whether such a universal spectrum can be reproduced in direct numerical simulations.
Similar universal spectra have been observed in direct numerical simulations
for other wave-interacting systems on water surfaces, in particular, for gravity-wave and capillary-wave systems
\citep{kz_PRL_2002,pushkarev,naoto_jfm,lvov2006dai,dyachenko2004wtk}.
Unlike these systems,
it appears from our direct numerical simulations
that energy spectra of internal waves are not universal.
Rather, accumulation of energy happens around the horizontally longest waves.
We demonstrate that
the energy transfer in inertial wavenumbers is dominated by scale-separated interactions with this accumulation.
Such an accumulation around the longest waves has been observed also in the oceans.
Furthermore,
deviations from the Garrett-Munk spectrum were recently categorized~\citep*{lvov-2004-92}
thus mounting additional doubts to the possibility of any universality of the oceanographic spectrum for internal waves.

\section{Hamiltonian formalism for internal waves}

Equations of motion for incompressible stratified fluid
(i.e. conservation of horizontal momentum, hydrostatic balance, mass
conservation and the incompressibility constraint)
under the assumption of zero potential vorticity
in the isopycnal coordinates, $(\bldx, \rho)$,
can be written as a pair of canonical
Hamiltonian equations,
\begin{equation}
 \frac{\partial \Pi}{\partial t}  = - \frac{\delta {\cal H}}{\delta \phi} \, ,
  \qquad
\frac{\partial \phi}{\partial t}  = \frac{\delta {\cal H}}{\delta \Pi} \, ,
  \label{Canonical}
\end{equation}
where $\Pi = \rho dz/d\rho$ is the fluctuation of stratification profile
around the mean, $-g / N_0^2$,
and $\phi$ is the isopycnal velocity potential.
The Hamiltonian is the sum of kinetic and potential energies,
%
\begin{eqnarray}
 {\cal H} = \int \mathrm{d} \bldx \mathrm{d} \rho
&&
 \left(
 - \frac{1}{2} \left(
-\frac{g}{N_0^2}+\Pi(\bldx, \rho) \right) \, 
  \left|\bnabla \phi(\bldx, \rho) - \frac{f N_0^2 }{g} \bnabla^{{\perp}}\Delta^{-1} \Pi(\bldx, \rho)
  \right|^2
\right.
\nonumber\\
&&
\quad
\left.
+
  \frac{g}{2} \left|\int^{\rho} \mathrm{d}\rho^{\prime} \frac{\Pi(\bldx, \rho^{\prime})}{\rho^{\prime}} \right|^2 
  \right) \, ,
  \label{HTL2}
\end{eqnarray}
%
(see \cite{lvov-2004-195} for complete details). 
%
%
The differential operators on the isopycnal surfaces,
$\bnabla=(\upartial/\upartial x, \upartial/\upartial y)$
and
$\bnabla^{\perp}=(-\upartial/\upartial y, \upartial/\upartial x)$,
are the horizontal gradient operator and the rotation operator, respectively.
Also, $g$ is the acceleration of gravity,
$N_0$ is the buoyancy (Brunt--V\"{a}is\"{a}l\"{a}) frequency,
$f$ is the inertial frequency due to the rotation of the Earth,
and $\rho_0$ is the mean density.

We then perform Fourier transformation and canonical transformation to
the field variable, $a(\bldp)$, as
  \begin{eqnarray}
   a(\bldp) = \sqrt{\frac{\omega}{2 g}} \frac{N_0}{|\bldk|} {\widetilde
\Pi}(\bldp)
    - \mathrm{i} \sqrt{\frac{g}{2 \omega}}\frac{|\bldk|}{N_0} {\widetilde
\phi}(\bldp) \, ,
  \end{eqnarray}
%
with linear coupling of the Fourier components of the stratification profile,
${\widetilde \Pi}$,
and the horizontal velocity potential, ${\widetilde \phi}$.
The three-dimensional wavenumber, $\bldp$, consists of
a two-dimensional horizontal wavenumber in the isopycnal surface, $\bldk$,
and a vertical density wavenumber, $m$.
The linear frequency is given by the dispersion relation,
%
  \begin{eqnarray}
   \omega(\bldp) = \sqrt{f^2 + \frac{g^2}{\rho_0^2 N_0^2}
\frac{|\bldk|^2}{m^2}} \, .
    \label{eq:dispersion}
  \end{eqnarray}
%
The usual vertical wavenumber, $k_z$, and the density wavenumber, $m$, are
related as
$m = - g/(\rho_0 N_0^2) k_z$.
The buoyancy frequency, $N_0$,
and the inertial frequency, $f$, are assumed to be constants.

Then, the equations of motion (\ref{Canonical}) are rewritten as a
canonical equation,
%
\begin{eqnarray}
    \mathrm{i} \frac{\upartial a(\bldp)}{\upartial t} = \frac{\delta {\cal
H}}{\delta a^{\ast}(\bldp)}\label{eq:canonical}
\end{eqnarray}
 with the standard Hamiltonian of three-wave interactive system,
\begin{eqnarray}
{\cal H} &=&  \int \mathrm{d}\bldp \: \omega(\bldp) |a(\bldp)|^2
   \nonumber\\
&& + \int \! \mathrm{d}\bldp \mathrm{d}\bldp_1 \mathrm{d}\bldp_2
   \left(
   \left(
   V_{\bldp_1,\bldp_2}^{\bldp} a(\bldp) a^{\ast}(\bldp_1)
a^{\ast}(\bldp_2) + \mathrm{c.c.}\right)
\right.
\nonumber\\
&&
\left.
\qquad\qquad\qquad\qquad
   +
   \left(
   U_{\bldp,\bldp_1,\bldp_2} a(\bldp) a(\bldp_1) a(\bldp_2) +
\mathrm{c.c.}\right)
   \right)
\, .\label{eq:hamiltonian}
\end{eqnarray}
Here, $\delta/\delta a^{\ast}$ is the functional derivative
with respect to $a^{\ast}(\bldp)$ that is the complex conjugate of $a(\bldp)$
and the abbreviation c.c. denotes complex conjugates.
Matrix elements, $V_{\bldp_1,\bldp_2}^{\bldp}$ and
$U_{\bldp,\bldp_1,\bldp_2}$, have exchange symmetries such that
  $V_{\bldp_1,\bldp_2}^{\bldp} = V_{\bldp_2,\bldp_1}^{\bldp}$
and
  $U_{\bldp,\bldp_1,\bldp_2} = U_{\bldp,\bldp_2,\bldp_1} =
U_{\bldp_1,\bldp,\bldp_2}$.

\section{Numerical methods}
\label{section:DNS}

\begin{table}
  \caption{Numerical parameters.}
  \label{table:parameters}
 \begin{center}
\begin{tabular}{cccccc}
  & modes     & $f$ ($\times 10^{-4}$rad/sec) & $L_{\mathrm{v}}$ ($\times 10$kg/m$^3$)& forcing & initial condition\\
\hline
 I & $512^2 \times 256$  & $\sqrt{2}/3$  & 2.7 & none & GM\\
 II & $1024^2 \times 512$ & $\sqrt{2}/3$  & 2.7 & $\omega \sim 3f$ & white noise\\
 III & $256^3$             & $0.25$        & 5   & $|\bldk|^2 \! + m^2 \! \leq 6^2$ & white noise\\
 IV  & $256^3$             & $1$           & 5   & $|\bldk|^2 \! + m^2 \! \leq 6^2$ & white noise\\
 V & $512^2 \times 256$  & $\sqrt{2}/3$  & 2.7 & none & GM w/o long waves
\end{tabular}
 \end{center}
\end{table}

In direct numerical simulations,
we add external forcing and hyper-viscosity
to the canonical equation (\ref{eq:canonical}),
as forcing and damping are needed to achieve statistically steady states.
Thus the resulting dynamic equation is given by
%
\begin{eqnarray}
 \frac{\upartial a(\bldp)}{\upartial t} = - \mathrm{i} \, \omega(\bldp) a(\bldp)
 + {\cal N}(a(\bldp)) + F(\bldp) - D(\bldp) a(\bldp) \, .
\end{eqnarray}
%
In the simulations,
the linear terms,
which are a linear dispersion term and a dissipation term, $-D(\bldp) a(\bldp)$,
are explicitly calculated.
The nonlinear terms, ${\cal N}(a(\bldp))$, derived from the nonlinear parts of the canonical equation (\ref{eq:canonical}) with Hamiltonian (\ref{eq:hamiltonian})
are obtained by a pseudo-spectral method with the phase shift
under the periodic boundary conditions for all three directions.
The external forcing, $F(\bldp)$, is implemented
by fixing   the amplitudes of several small wavenumbers to be constant in time.
The dissipation is also added as
$D(\bldp) = D_{\mathrm{h}} |\bldk|^8 + D_{\mathrm{v}} |m|^4$.
Here, $ D_{\mathrm{h}}$ and $D_{\mathrm{v}}$ are chosen
so that the dissipation is effective for wavenumbers
larger than the half of maximum wavenumbers.
The wavenumbers are discretized as $\bldp = (2\upi/ L_{\mathrm{h}} \bldk, \: 2\upi/ L_{\mathrm{v}} m)$, where $L_{\mathrm{h}}$ and $L_{\mathrm{v}}$ are
horizontal periodic length and vertical period in the isopycnal coordinates,
and $\bldk$ and $m$ are integer-valued wavenumbers afterward.
Time-stepping is implemented with the fourth-order Runge--Kutta method.
In all the simulations,
the buoyancy frequency and horizontal period are fixed at $N_0 = 10^{-2}$rad/sec and $L_{\mathrm{h}}=10^5$m, respectively.

We perform a series of five numerical experiments that are listed in Table~\ref{table:parameters}.
The total energy per unit periodic box of all the numerical experiments except Run II
is around $3 \times 10^3$J/(kg $\cdot$ m$^2$) which is characteristic of real oceans.
The total energy in the Run II is around $1.2 \times 10^3$J/(kg $\cdot$ m$^2$).

Two-dimensional energy spectra are measured in the shell of radius $k$ as
%
\begin{eqnarray}
E(k, |m|) = \int\limits_{k- 1/2 \leq |\bldk^{\prime}| < k+ 1/2} \mathrm{d}\bldk^{\prime} \sum_{s=\pm 1} \omega |a(\bldk^{\prime}, s m)|^2 \, .
\end{eqnarray}
%
Integrated energy spectra are defined as
$\overline{E}_{\mathrm{int}}(k) = \int \mathrm{d}m E(k,|m|)$ and
$\overline{E}_{\mathrm{int}}(|m|) = \int \mathrm{d}k E(k,|m|)$.
Cross-sectional spectra, $E_m(k)$ and $E_k(|m|)$,
are obtained from the two-dimensional energy spectra
as a function of isopycnal wavenumbers $k$
along a certain density wavenumber $m$
and
as a function of density wavenumbers $m$
along a certain isopycnal wavenumber $k$, respectively.
%

In general in real oceans
only integrated spectra,
$\overline{E}_{\mathrm{time}}(\omega)$ obtained from time series and $\overline{E}_{\mathrm{vertical}}(m)$ from vertical series,
can be measured.
Then the two-dimensional spectrum, $E(\omega,m)$, is composed from integrated spectra.
To do it,
the assumption of separability of two-dimensional spectra is made.
Namely, the functional form,
$E(\omega,m) \propto \overline{E}_{\mathrm{time}}(\omega) \overline{E}_{\mathrm{vertical}}(m)$,
is assumed.
As we will see below in our numerical simulations,
our resulting spectra are not always separable.
Such non-separability is also observed in the real oceans \citep{kurtPC}.
Note that the power-law exponents of one-dimensional spectra and those of two-dimensional spectra do not always coincide in our direct numerical simulation,
as explained below.
We do not have any means to see whether this observation will also hold in the real ocean.

\begin{figure}
 \begin{center}
  \includegraphics[scale=0.71]{gm0_psfrag.eps}%
  \includegraphics[scale=0.71]{gm4_psfrag.eps}
 \end{center}
 \caption{Run I. Cross-sectional spectra $E_k(|m|)$
 of GM spectrum as the initial condition (left)
 and
 of energy spectrum after about 35 hours (right).
 Significant differences in the region $16 < |m| < 64$
 indicates that GM spectrum is not statistically steady.
 }
\label{fig:gm_ns}
\end{figure}
%


\section{Results}
%

In the first numerical experiment, Run I of Table~\ref{table:parameters},
we examine statistical stability of the GM spectrum of a freely decaying system,
i.e. without external forcing.
The GM spectrum with the cut-off in large isopycnal and density wavenumbers
is employed as the initial energy spectrum.
The cut-off is introduced to avoid decreasing accuracy of the pseudo-spectral method.
The initial phases of complex amplitude, $a(\bldp)$, are given by uniformly distributed random numbers in $[ 0,2\upi )$.
Figure~\ref{fig:gm_ns} shows the initial spectrum and the energy spectrum after about 35 hours in ocean time.
If the GM spectrum were to be a universal steady state for this numerical run, it would change very little.
Instead,
the density exponents rapidly change from $-1$ to $-2$ only after $1.5$ days.
This spectral change occurs faster than dissipation affects
with timescale estimated to be about 50 days at $|m|=32$.
Therefore it appears that
the GM spectrum is not a stable universal spectrum in the wavenumber region $16 < |m| < 64$
at least for this numerical experiment.
This suggests that
the diffusion approximation of the kinetic equation \citep{mccomas-1977-7,mccomas-1981-11},
which predicts that energy spectra are rapidly relaxed to the GM spectrum,
may not be fully appropriate.

\begin{figure}
 \begin{center}
  \includegraphics[scale=0.8]{spectra2d_psfrag.eps}
 \end{center}
 \caption{ Run II. Two-dimensional energy spectrum $E(k,|m|)$.
 The contours are plotted every powers of ten from $10^{-10}$ to $1$.
 The straight line from the origin shows the wavenumbers $\omega=2f$.
 The inset is the enlargement near the origin.
 }
 \label{fig:energy_spectrum_2D}
\end{figure}

Hoping to find a statistically steady state,
we then performed a forced-damped simulation, Run II (see Table~\ref{table:parameters}).
The external forcing is added to 24 small wavenumbers whose frequencies are close to $3f$
such that $(\bldk, m) = (\pm 1, 0, \pm 2)$, $(0, \pm 1, \pm 2)$, $(\pm 1, \pm 1, \pm 3)$,
$(\pm 2, 0, \pm 4)$ and $(0, \pm 2, \pm 4)$.
The two-dimensional energy spectrum in near statistically steady state
obtained in this simulation is shown in Fig.~\ref{fig:energy_spectrum_2D}.
We observe that there is a significant energy accumulation around $k \sim 1$,
i.e. the horizontally longest waves.
Note that linear frequency given by dispersion relation is exactly $f$ when $k=0$.
Similarly, accumulation of energy around $k \sim 1$ corresponds to accumulation of energy in near-inertial frequencies.
This is precisely what happens in the ocean!
Indeed all versions of the GM spectra
have an integrable singularity in near-inertial frequencies
such as $1/\sqrt{\omega^2 - f^2}$ \citep{GM72,GM75,GM_ARF}, indicating significant energy contents in near-inertial waves.
Such an accumulation can be explained quantitatively as a result of
a parametric subharmonic instability \citep*{furuich-2005} which
arises in small isopycnal wavenumbers.

\begin{figure}
 \begin{center}
  \includegraphics[scale=0.71]{spectra_k_psfrag.eps}%
  \includegraphics[scale=0.71]{spectra_m_psfrag.eps}
 \end{center}
 \caption{Run II.
 Integrated energy spectra $\overline{E}_{\mathrm{int}}(k)$ and $\overline{E}_{\mathrm{int}}(|m|)$,
 and cross-sectional energy spectra $E_m(k)$ and $E_k(|m|)$.
 GM spectrum has $k^{-2} |m|^{-1}$ in large isopycnal and density wavenumbers.
 The cross-sectional spectra are shown every four curves for visibility.
 }
 \label{fig:spectra}
\end{figure}

Figure~\ref{fig:spectra} shows integrated spectra,
$\overline{E}_{\mathrm{int}}(k)$ and $\overline{E}_{\mathrm{int}}(|m|)$,
and cross-sectional spectra, $E_m(k)$ and $E_k(|m|)$
obtained from two-dimensional energy spectrum shown in Fig.~\ref{fig:energy_spectrum_2D}.
The integrated spectra are
$\overline{E}_{\mathrm{int}}(k) \propto k^{-1.98 \pm 0.02}$
and
$\overline{E}_{\mathrm{int}}(|m|) \propto |m|^{-1.33 \pm 0.31}$.
The exponents are not too far from the large-wavenumber self-similar form of GM spectrum,
that is $k^{-2} |m|^{-1}$.
The exponents are obtained with the least-square method
in the intervals of $k \in [8,128]$ for $\overline{E}_{\mathrm{int}}(k)$
and that of $|m| \in [4,32]$ for $\overline{E}_{\mathrm{int}}(|m|)$.

\begin{figure}
 \begin{center}
  \includegraphics[scale=0.71]{exp_k_psfrag.eps}%
  \includegraphics[scale=0.71]{exp_m_psfrag.eps}
 \end{center}
 \caption{Run II.
 Power-law exponents of each cross-sectional spectrum.
 left: $\alpha(|m|)$, which are the exponents of cross-sectional energy spectra $E_{m}(k)$,
 right: $\beta(k)$, which are the exponents of cross-sectional energy spectra $E_{k}(|m|)$.
 The exponents of integrated spectra,
 $\alpha_{\mathrm{int}}=-1.98$ and $\beta_{\mathrm{int}}=-1.33$, are also shown.
 }
 \label{fig:exp_css}
\end{figure}
%
\begin{figure}
 \begin{center}
  \includegraphics[scale=0.71]{local_d_a_psfrag.eps}%
  \includegraphics[scale=0.71]{local_d_b_psfrag.eps}
 \end{center}
 \caption{Run II.
 Local power-law exponents.
 The exponents are derived from $E_{m=16}(k)$, $E_{m=48}(k)$, $E_{k=32}(m)$ and $E_{k=96}(m)$.
 }
 \label{fig:local_exponents}
\end{figure}
Figure~\ref{fig:exp_css} shows power-law exponents of each cross-sectional spectrum,
$\alpha(|m|)$ and $\beta(k)$,
which are obtained by fitting
$E_{m}(k) \propto k^{\alpha(|m|)}$ in $k \in [16,128]$
and $E_{k}(m) \propto |m|^{\beta(k)}$ in $|m| \in [16,64]$.
The exponents of integrated spectra, $\alpha_{\mathrm{int}}=-1.98$ and $\beta_{\mathrm{int}}=-1.33$,
are also shown in Fig.~\ref{fig:exp_css}.
Note that the fitted regions of cross-sectional spectra are different from those of integrated spectra.
The integrated spectra have different exponents from the ones of cross-sectional spectra in inertial region
since the integrated spectra are greatly affected
by the accumulation in small isopycnal or density wavenumbers.
Especially, the discrepancy is clear in the density exponent $\beta$
because of the strong accumulation around $k \sim 1$.
Apparently,
the cross-sectional exponents are not constant.
Therefore,
the energy spectrum of this numerical simulation cannot be accurately fitted by double-power functions.
It also appears that this spectrum is non-separable.
Moreover,
in contrast to the integrated spectra,
both cross-sectional spectra have much steeper exponents, $\alpha, \beta \sim -3.5$,
than the GM spectrum has
in large isopycnal and density wavenumbers considered as inertial region.
Local exponents of cross-sectional spectra are defined as
\begin{eqnarray}
 \alpha_m^{\mathrm{local}}(k) = \frac{\mathrm{d} \log E_m(k)}{\mathrm{d} \log k}, \qquad
 \beta_k^{\mathrm{local}}(m) = \frac{\mathrm{d} \log E_k(m)}{\mathrm{d} \log m}.
\end{eqnarray}
We observe that
the local exponents are also close to $-3.5$ in large isopycnal and density wavenumbers
shown in Fig.~\ref{fig:local_exponents},
although the local exponents, $\alpha_m^{\mathrm{local}}(k)$ in particular, are jagged.
Poor redistribution of energy from
the accumulation of energy in small isopycnal wavenumbers
and that in small density wavenumbers
can make steeper spectra as well.

\begin{figure}
 \begin{center}
  \includegraphics[scale=0.71]{spectra_k_s1_psfrag.eps}%
  \includegraphics[scale=0.71]{spectra_m_s1_psfrag.eps}
 \end{center}
 \caption{Run III.
 Integrated and cross-sectional spectra when the accumulation of energy around the horizontally longest waves is weak.
 }
 \label{fig:energy_spectra_weak}
\end{figure}
\begin{figure}
 \begin{center}
  \includegraphics[scale=0.71]{spectra_k_s2_psfrag.eps}%
  \includegraphics[scale=0.71]{spectra_m_s2_psfrag.eps}
 \end{center}
 \caption{Run IV.
 Integrated and cross-sectional spectra when the accumulation of energy around the horizontally longest waves is
 moderate.}
 \label{fig:energy_spectra_mod}
\end{figure}
\begin{figure}
 \begin{center}
  \includegraphics[scale=0.71]{exp_k_s1_psfrag.eps}%
  \includegraphics[scale=0.71]{exp_m_s1_psfrag.eps}
\end{center}
  \caption{Run III.
 Power-law exponents when the accumulation of energy around the horizontally longest waves is weak.
 }
  \label{fig:exp_weak}
\end{figure}
\begin{figure}
 \begin{center}
  \includegraphics[scale=0.71]{exp_k_s2_psfrag.eps}%
  \includegraphics[scale=0.71]{exp_m_s2_psfrag.eps}
 \end{center}
  \caption{Run IV.
 Power-law exponents when the accumulation of energy around the horizontally longest waves is moderate.
 }
  \label{fig:exp_mod}
\end{figure}

Next,
we investigate the influence of this accumulation of energy on inertial wavenumbers.
To achieve this goal,
we perform two numerical runs III and IV in Table~\ref{table:parameters}.
These two runs are different by value of inertial frequencies $f$ (see Table~\ref{table:parameters}).
The difference of values of inertial frequencies influences the rate of formation of the accumulation around the horizontally longest waves.
After about $10^3$ days from the initial time
when all the wavenumbers have extremely small energy,
the system is still transient
and the accumulation of energy around the horizontally longest waves has not developed strongly.
The timescales of the formation of the accumulation,
which are determined by nonlinear interactions among the long waves,
are much longer than those of the inertial wavenumbers.
The energy spectrum of the case with $f=0.25 \times 10^{-4}$rad/sec has only weak accumulation
since all the forced wavenumbers have frequencies greater than $4f$,
far from the linear frequencies of the wavenumbers that constitute the accumulation.
The energy spectrum of the case with $f=10^{-4}$rad/sec has moderate accumulation,
which is stronger than that of the case with $f=0.25 \times 10^{-4}$rad/sec
and is weaker than that of Run II.

The integrated and cross-sectional spectra and their power-law exponents
after $10^3$ days are shown
in Figs.~\ref{fig:energy_spectra_weak}, \ref{fig:energy_spectra_mod}, \ref{fig:exp_weak} and \ref{fig:exp_mod}.
The least-square fitting to obtain the power-law exponents is made in the intervals of $k,m \in [8,40]$.
When the accumulation of energy around the horizontally longest waves is weak enough,
the energy spectrum is close to the double-power $E(k,m) \propto k^{-2} |m|^{-2.5}$
in large isopycnal and density wavenumbers.
On the other hand,
the energy spectrum does not have double-power laws
in the region when the accumulation is moderate.
It is caused by {\itshape contamination} of the inertial region due to the accumulation.
Another double-power law $E(k,m) \propto k^{-2} |m|^{-2}$ appears
in small isopycnal and density wavenumbers.
Again, we point out that
the exponents of cross-sectional spectra are also sufficiently different from those of the integrated spectra and the GM spectrum.
Therefore, in combination with Run II,
it appears that
the accumulation of energy around the horizontally longest waves strongly affects
statistical properties in the inertial region.

\begin{figure}
 \begin{center}
  \includegraphics[scale=0.8]{rd1_psfrag.eps}
 \end{center}
 \caption{Relative difference between energy spectra
 after about $1.5$ days
 developed from GM spectrum with and without energy in small isopycnal wavenumbers
 $E_{\mathrm{d}}(k,|m|)$.
 }
 \label{fig:gm_nl}
\end{figure}

To further illustrate this point, and to investigate the importance of nonlocal interactions in the wavenumber space,
we perform Run V of Table~\ref{table:parameters}.
There we choose the initial condition
same as the Run I but with no energy in $k<3$ and $|m| < 16$.
The energy spectrum after about $1.5$ days is denoted by $E_{\mathrm{nl}}(k,|m|)$.
The energy spectrum developed from the GM spectrum in Run I that is shown in Fig.~\ref{fig:gm_ns}(right)
is denoted by $E_{\mathrm{GM}}(k,|m|)$.
Figure~\ref{fig:gm_nl} shows the relative difference defined as
\begin{eqnarray}
%E_{\mathrm{d}}(k,|m|) = (E_{\mathrm{nl}}(k,|m|) - E_{\mathrm{GM}}(k,|m|))/ E_{\mathrm{GM}}(k,|m|).
E_{\mathrm{d}}(k,|m|) = \frac{E_{\mathrm{nl}}(k,|m|) - E_{\mathrm{GM}}(k,|m|)}{E_{\mathrm{GM}}(k,|m|)} \, .
\end{eqnarray}
%
The wavenumbers different in initial conditions appear as a black rectangle in bottom left in Fig.~\ref{fig:gm_nl}.
%
If the nonlocal interactions with the accumulation of energy around the horizontally longest waves were not dominant,
the relative difference in Fig.~\ref{fig:gm_nl} would be small
or slightly negative in inertial wavenumbers
since the nonlinear interactions try to compensate the defect in the small wavenumbers.
%%%%%%%%%%%%%%%%%%%%%%%%%%
Instead, $E_{\mathrm{nl}}(k,|m|)$
has about 50\% more energy in $10 \lesssim |m| \lesssim 40$ drawn in white
and less energy transfer to the dissipation region in $|m| \gtrsim 80$.
This suggests that
small wavenumbers in the accumulation of energy
make energy transfer
from small density wavenumbers to large density wavenumbers
in the inertial region.
Energy transfer to large density wavenumbers is qualitatively consistent to
the induced diffusion \citep{mccomas-1977-7,mccomas-1981-11}.

\section{Concluding remarks}
\label{sec:summary}

We tried to reproduce the {\itshape universal} energy spectrum of oceanic internal waves by direct numerical simulations.
Instead of expected universal behaviour in large wavenumbers,
the simulations show a tendency
for energy to concentrate around the horizontally longest waves that are determined by periodic boundary condition.
It also appears that the power-law exponents of
integrated spectra are different from those of cross-sectional spectra.
This is apparent in density-wavenumber spectra,
whose integrated spectrum is strongly affected
by the accumulation of energy in small isopycnal wavenumbers.
It is shown that the accumulation strongly depends on the local value of inertial frequencies.
%
We also find that
the nonlocal interactions in the wavenumber space seem to be dominant in large isopycnal and density wavenumbers.
%
Energy spectra in these wavenumbers are greatly affected by small isopycnal wavenumbers,
since nonlocal interactions in the wavenumber space are dominant in energy transfer in large isopycnal and density wavenumbers.

On the other hand, energy spectra in the small wavenumbers cannot be universal
since they are
box-size dependent in direct numerical simulations with periodic boundary conditions.
In the real oceans, $\beta$-effect changes $f$
and is important for Rossby waves, which are the longest waves in oceans.
The topographies of seabeds and the alignment of continents
can also alter the longest waves.
The temperature gives seasonal variability
and the atmospheric disturbance also gives fluctuations in large scales.
Indeed, it is known that the large-scale motions cannot be universal in the real oceans.

Therefore,
the small-scale properties cannot be universal
if the nonlocal interactions are dominant
in the energy transfer in large wavenumbers.
It is in contrast to the surface-wave systems,
which have the universal energy spectra in small scales as a result of {\em cascades\/}
even if large-scale behavior is not universal.

\vspace{\baselineskip}
{\bf Acknowledgmets}
%\begin{acknowledgments}
This research is supported by NSF CMG grant 0417724.  Y.~L. was also
supported by NSF CAREER DMS 0134955.  We are gratefull to Kurt Polzin
and Esteban Tabak for helpfull discussions.  We are also thank YITP in
Kyoto University (Japan) for allowing us to use SX8, where numerical
simulations were performed.

%\end{acknowledgments}

\bibliographystyle{jfm}
\begin{thebibliography}{14}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi

%\bibitem[Dyachenko {\em et~al.\/}(2004)Dyachenko, Korotkevich \& Zakharov]{dyachenko2004wtk}
\bibitem[Dyachenko, Korotkevich \& Zakharov(2004)]{dyachenko2004wtk}
{\sc Dyachenko, A.~I., Korotkevich, A.~O. \& Zakharov, V.~E.} 2004 Weak
  turbulent {K}olmogorov spectrum for surface gravity waves. {\em Phys. Rev.
  Lett.\/} {\bf 92}, 134501.

\bibitem[Furuich {\em et~al.\/}(2005)Furuich, Hibiya \& Niwa]{furuich-2005}
{\sc Furuich, N., Hibiya, T. \& Niwa, Y.} 2005 Bispectral
  analysis of energy transfer within the two-dimensional oceanic internal wave
  field. {\em J. Phys. Oceanogr.\/} {\bf 35}, 2104--2109.

\bibitem[Garrett \& Munk(1972)]{GM72}
{\sc Garrett, C. J.~R. \& Munk, W.~H.} 1972 Space-time scales of internal
  waves. {\em Geophys. Fluid Dyn.\/} {\bf 3}, 225--264.

\bibitem[Garrett \& Munk(1975)]{GM75}
{\sc Garrett, C. J.~R. \& Munk, W.~H.} 1975 Space-time scales of internal
  waves: a progress report. {\em J. Geophys. Res.\/} {\bf 80}, 281--297.

\bibitem[Garrett \& Munk(1979)]{GM_ARF}
{\sc Garrett, C. J.~R. \& Munk, W.~H.} 1979 Internal waves in the ocean. {\em
  Annu. Rev. Fluid Mech.\/} {\bf 11}, 339--369.

%\bibitem[Lvov {\em et~al.\/}(2006)Lvov, Nazarenko \& Pokorni]{lvov2006dai}
\bibitem[Lvov, Nazarenko \& Pokorni(2006)]{lvov2006dai}
{\sc Lvov, Y.~V., Nazarenko, S. \& Pokorni, B.} 2006 Discreteness and
  its effect on water-wave turbulence. {\em Physica D\/} {\bf 218}, 24--35.

\bibitem[Lvov {\em et~al.\/}(2004)Lvov, Polzin \& Tabak]{lvov-2004-92}
{\sc Lvov, Y.~V., Polzin, K.~L. \& Tabak, E.~G.} 2004 Energy spectra
  of the ocean's internal wave field: theory and observations. {\em Phys. Rev.
  Lett.\/} {\bf 92}, 128501.

\bibitem[Lvov \& Tabak(2004)]{lvov-2004-195}
{\sc Lvov, Y.~V. \& Tabak, E.~G.} 2004 A {H}amiltonian formulation for
  long internal waves. {\em Physica D\/} {\bf 195}, 106--122.

\bibitem[McComas(1977)]{mccomas-1977-7}
{\sc McComas, C.~H.} 1977 Equilibrium mechanisms within the oceanic internal
  wave field. {\em J. Phys. Oceanogr.\/} {\bf 7}, 836--845.

\bibitem[McComas \& M\"{u}ller(1981)]{mccomas-1981-11}
{\sc McComas, C.~H. \& M\"{u}ller, P.} 1981 The dynamic balance of
  internal waves. {\em J. Phys. Oceanogr.\/} {\bf 11}, 970--986.

\bibitem[Onorato {\em et~al.\/}(2002)Onorato, Osborne, Serio, Resio, Pushkarev,
  Zakharov \& Brandini]{kz_PRL_2002}
{\sc Onorato, M., Osborne, A.~R., Serio, M., Resio, D., Pushkarev, A.,
  Zakharov, V.~E. \& Brandini, C.} 2002 Freely decaying weak turbulence for sea
  surface gravity waves. {\em Phys. Rev. Lett.\/} {\bf 89}, 144501.

\bibitem[Polzin(2007)]{kurtPC}
{\sc Polzin, K.~L.} 2007 Private communication.

\bibitem[Pushkarev \& Zakharov(2000)]{pushkarev}
{\sc Pushkarev, A.~N. \& Zakharov, V.~E.} 2000 Turbulence of capillary waves --
  theory and numerical simulation. {\em Physica D\/} {\bf 135}, 98--116.

\bibitem[Yokoyama(2004)]{naoto_jfm}
{\sc Yokoyama, N.} 2004 Statistics of gravity waves obtained by direct
  numerical simulation. {\em J. Fluid Mech.\/} {\bf 501}, 169--178.

\end{thebibliography}

\end{document}
