\documentclass[12pt]{article}

\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{bm}
\usepackage{natbib}
\usepackage{fullpage}
\newcommand \rewrite {{\bf{rewrite\ here}}}
\newcommand \finish {{\bf{end\ of\  rewrite}}}

\topskip 2cm


\newcommand{\todo}[1]{\vspace{5 mm}\par \noindent
\marginpar{\textsc{ToDo}} \framebox{\begin{minipage}[c]{0.9
\textwidth} \tt #1 \end{minipage}}\vspace{5 mm}\par}


\title{Wave-wave interactions in stratified fluids:\\ A comparison
of approaches.}
%
\author{Yuri V Lvov$^1$, Kurt Polzin$^2$ and Naoto Yokoyama$^3$\\
{\small $^1$ Department of Mathematical Sciences, Rensselaer
Polytechnic Institute, Troy NY 12180}\\
%
{\small $^2$ Woods Hole Oceanographic Institution, MS\#21, Woods Hole, MA 02543}\\
%
{\small $^3$ Department of Mechanical Engineering, Doshisha University, Kyotanabe, Kyoto 610-0394 JAPAN}
}

\begin{document}

\maketitle

\begin{abstract}
Various approaches have been developed over the last four
decades to characterize the magnitude of nonlinear interactions between
triads of internal waves in stratified oceanic
flows. The present manuscript compares some of these approaches and
their predictions for internal wave nonlinearity parameter and Boltzman rate. We demonstrate that, for {\em resonant\/} triads in the limit of long internal waves in hydrostatic balance and in the absence of rotation, these various approaches
predict {\em equivalent\/} rates of energy transfer between waves. However, 
with the inclusion of background rotation and off-resonant interactions,
these approaches lead to qualitatively different predictions. 
In particular, a noncanonical approach in Lagrangian coordinates leads to higher levels of nonlinearity at high frequencies and large wavenumbers than a canonical approach in isopycnal coordinates.  

\end{abstract}

\section{Introduction}

Wave-wave interactions in stratified oceanic flows have been
a subject of intensive research in the last four  decades. Of particular
importance is the existence of a ``universal'' internal-wave spectrum, the Garrett and Munk spectrum. It is
generally perceived that the existence of a universal spectrum is, at least in part 
and perhaps even primarily, the result 
of nonlinear interactions of waves with different wavenumbers. Due to
the quadratic nonlinearity of the underlying primitive equations and the fact that the linear internal-wave dispersion relation can satisfy a three-wave resonance condition, waves interact
in triads. Therefore the question arises: how strongly do waves with wavenumber $\bm{p}$ interact with wavevectors $\bm{p}_1$ and $\bm{p}_2$, where $\bm{p} = \bm{p}_1 + \bm{p}_2$?  
What are the oceanographic consequences of this interaction?

Various approaches have been developed to
characterize the magnitude of such interactions~\citep{H66,K66,K68,McC75,MO75,O74,O76,PR77,MB77,PMW80,Voronovich,Milder,Zeitlin,LT,LT2}
(see Table~\ref{TABLEOFELEMENTS} for a summary of the major distinctions.)

All these approaches represent various attempts to derive a closed equation representing the slow time evolution of the wave field's wave action spectrum.  Such an equation is called a {\em kinetic equation}~\citep{ZLF}.  
%In this paper we concentrate of four of these approaches, namely \citep{MO75,Voronovich,Zeitlin,LT}.  


\begin{table}[htbp]
\label{TABLEOFELEMENTS}
\caption{A list of various kinetic equations.  Results from \cite{O76, MB77, PMW80} are reviewed in \cite{M86}, who state that \cite{O76}, \cite{MB77} and an unspecified Eulerian representation are consistent on the resonant manifold.  \cite{PMW80} utilizes Langevin techniques to assess nonlinear transports.   \cite{M86} characterizes those Langevin results as being mutually consistent with the direct evaluations of kinetic equations presented in \cite{O76, MB77}.  \cite{K68} states (without detail) that \cite{K66} and \cite{H66} give numerically similar results.  A formulation in terms of discrete modes will typically permit an arbitrary buoyancy profile, but obtaining results requires specification of the profile.  Of the discrete formulations, \cite{PMW80} use an exponential profile and the others assume a constant stratification rate.
The kinetic equations marked by $^{\dag}$  are investigated in \S\ref{ResonantInteractions}, while kinetic equations marked by $^{\ddag}$ 
are investigated further in \S\ref{OffResonant}.
}
\label{default}
\begin{center}
\begin{tabular}{cccccc}
\hline
source & coordinate & vertical & rotation & hydro- & special \\
 & system & structure & & static & \\
\hline
\citet{H66} & Lagrangian & discrete & no & no & \\
\citet{K66, K68} & Eulerian & discrete & no & no & non-Hamiltonian \\
\citet{MO75}$^{\dag \ddag}$ & Lagrangian & cont. & yes & no &\\
\citet{McC75, McComas} & Lagrangian & cont. & yes & yes & \\
\citet{PR77} & Lagrangian & cont. & no & no &  Clebsh\\
\citet{Voronovich}$^{\dag}$ & Eulerian & cont. & no & no & Clebsh \\
\citet{PMW80} & Lagrangian & discrete & yes & no & Langevin \\
\citet{Milder} & Isopycnal & n/a & no & no &  \\
\citet{Zeitlin}$^{\dag}$ & Eulerian & cont. & no & no & non-Hamiltonian \\
\citet{LT}$^{\dag}$ & Isopycnal & cont. & no & yes & canonical \\
\citet{LT2}$^{\ddag}$ & Isopycnal & cont. & yes & yes & canonical \\
\hline
\end{tabular}
\end{center}
\end{table}



In this manuscript we concentrate on four of these different versions of the internal-wave kinetic equation:  
%
\begin{itemize}
\item a noncanonical description using Lagrangian coordinates \citep{O74,O76,MO75},
\item  a canonical Hamiltonian description using Clebsh variables in Eulerian coordinates \citep{Voronovich},
\item a dynamical derivation of a kinetic equation without use of Hamiltonian formalisms in Eulerian coordinates \citep{Zeitlin},
\item a canonical Hamiltonian description in isopycnal coordinates \citep{LT,LT2}.
\end{itemize}
%



Our intent is to compare these approaches, and in particular, compare the predictions for the wavenumber-dependent characteristic nonlinear time scale of the Garrett and Munk wave action spectrum. 
To achieve this goal, 
we give necessary background in Section~\ref{Background}, briefly review
approaches of Table~\ref{TABLEOFELEMENTS} in Section~(\ref{VariousApproaches}), and then 
 we demonstrate in Section~\ref{ResonantInteractions} that, under assumption of hydrostatic balance and under the assumption of {\em resonant\/} wave-wave interactions, the interaction matrices associated with the listed approaches are {\em equivalent\/}.  While one, with sufficient experience, might regard this as an intuitive statement, it is far from trivial.  We will then demonstrate in Section \ref{OffResonant} that, if the assumption of resonant wave-wave interactions is relaxed, both quantitatively and qualitatively different transfer rates are predicted. In particular, we show that  the Boltzman rate 
  ${\cal  \epsilon}_{\bm{p}}$, defined below in  (\ref{NonlinearTime}) is, in fact,  representation dependent for near-resonant interactions.  

We have {\em not\/}, at this time, achieved a detailed mathematical understanding of how these differences arise and consequently do not digress into a detailed discussion of why, for example, the radius of convergence of two consecutive series expansions in one coordinate system differs so dramatically from a single series expansion in a different coordinate system.  We conclude in Section \ref{Conclusion}.



\section{Background\label{Background}}
A kinetic equation is a closed equation for the time evolution of the
wave action spectrum in a system of weakly interacting
waves. It is usually derived as a central result of wave turbulence
theory. The concepts of wave turbulence theory
provide a fairly general framework
for studying the statistical steady states in a large class of weakly
interacting and weakly nonlinear many-body or many-wave systems.  In
its essence, classical weak turbulence theory~\citep{ZLF} is a
perturbation expansion in the amplitude of the nonlinearity, yielding,
at the leading order, linear waves, with amplitudes slowly modulated
at higher orders by resonant nonlinear interactions.  This modulation
leads to a resonant redistribution of the spectral energy density among
space- and time-scales, and is described by a kinetic equation. 

Typical assumptions needed for the derivation of kinetic equations are: 
%
\begin{itemize}
\item Weak nonlinearity, 
\item Gaussian statistics of the interacting wave field in wavenumber space and
\item Resonant wave-wave interactions.  
\end{itemize}
%

The derivation of the kinetic equation for general nonlinear systems is well studied and 
understood, and thus will not be repeated  here. Three wave 
kinetic equations take the form~\citep{ZLF,NoisyNazarenko,LLNZ}:
%
\begin{eqnarray}
\frac{d n_{\bm{p}}}{dt} = 4\pi \int
 |V_{\bm{p}_1,\bm{p}_2}^{\bm{p}}|^2 \, f_{p12} \,
\delta_{{\bm{p} - \bm{p}_1-\bm{p}_2}} \, \delta({\omega_{\bm{p}}
-\omega_{{\bm{p}_1}}-\omega_{{\bm{p}_2}}})
d \bm{p}_{12}
\nonumber \\
-4\pi\int
 \, |V_{\bm{p}_2,\bm{p}}^{\bm{p}_1}|^2\, f_{12p}\, \delta_{{\bm{p}_1 - \bm{p}_2-\bm{p}}} \,
  \delta({{\omega_{\bm{p}_1} -\omega_{\bm{p}_2}-\omega_{\bm{p}}}})
\, d \bm{p}_{12}
\nonumber \\
-4\pi\int
 \, |V_{\bm{p},\bm{p}_1}^{\bm{p}_2}|^2\, f_{2p1}\, \delta_{{\bm{p}_2 - \bm{p}-\bm{p}_1}} \,
  \delta({{\omega_{\bm{p}_2} -\omega_{\bm{p}}-\omega_{\bm{p}_1}}})
\, d \bm{p}_{12}
 \, ,\nonumber\\
{\rm with} ~~ f_{p12} = n_{\bm{p}_1}n_{\bm{p}_2} -
n_{\bm{p}}(n_{\bm{p}_1}+n_{\bm{p}_2}) \, .
\label{KineticEquation}
\end{eqnarray}
%
Here $n_{\bm{p}} = n(\bm{p})$ is a three-dimensional wave action spectrum (spectral energy density divided by frequency) and the interacting wavevectors $\bm{p}$, $\bm{p}_1$ and $\bm{p}_2$ are given by
%
$$\bm{p} = (\bm{k}, m),$$
%
i.e.\ $\bm{k}$ is the horizontal part of $\bm{p}$ and $m$ is its vertical component.  We assume the wavevectors are signed variables and wave frequencies $\omega_{\bm{p}}$ are restricted to be positive.  The magnitude of wave-wave interactions 
$V_{\bm{p},\bm{p}_1}^{\bm{p}_2}$
is a matrix representation of the coupling between triad members.  It serves
as a multiplier in the nonlinear convolution 
term in what is now commonly called 
the Zakharov equation -- equation in the Fourier space for the waves 
field variable. This is also an expression that multiplies the cubic 
convolution term in the three-wave Hamiltonian. 

For internal waves in the ocean such kinetic equation was derived by the approaches 
in Table~\ref{TABLEOFELEMENTS}.
The development of a kinetic equation is facilitated by transforming to canonical coordinates in a Hamiltonian framework, for which one can demonstrate that the symmetries and hence conservation principles of the original equation set in the spatial/temporal domains have been preserved in the spectral domain~\citep[e.g.][]{ZLF}.  Finding canonical coordinates, however, can be highly nontrivial.  Transformations to canonical coordinates have been found using Clebsch variables in Eulerian coordinates \citep{Voronovich} and in isopycnal coordinates \citep{LT, LT2}.  Kinetic equations in Lagrangian coordinates start by averaging the Lagrangian of the stratified fluid (averaging of the variational principle) and then transform the Lagrangian to the Hamiltonian. The Lagrangian coordinate kinetic equations considered here are noncanonical.  We also note that it is possible to obtain a kinetic equation directly from the dynamical equations of motion, without the use of the Hamiltonian structure. Such an approach was executed by~\citet{Zeitlin}.  The conservation properties of non-canonical and non-Hamiltonian representations are not guaranteed unless explicitly demonstrated.  The issue of conservation properties is greatly compounded for non-canonical and non-Hamiltonian representations off the resonant manifold.  

A typical restriction is to exclude interactions with potential vorticity carrying members of the fluid dynamical system, for which we refer the reader to \cite{Lelong} and \cite{Zeitlin}.  Even in these extended analyses a plane wave formulation is assumed that eliminates the potential vorticity associated with a slowly varying wave-packet structure \citep{BM05, P08}.  


Note that the kinetic equation allows us to numerically estimate the life time
of any given spectrum. In particular, we can define a wavenumber dependent nonlinear time scale proportional to the inverse Boltzman rate:
%
\begin{equation}
\tau^{\mathrm{NL}}_{\bm{p}} = \frac{n_{\bm{p}}}{\dot n_{\bm{p}}}~.
\label{NonlinearTime}
\end{equation}
%
This time scale characterizes the net rate at which the spectrum changes and can 
be directly calculated from the kinetic equation.  

One can also define the 
characteristic linear time scale, 
$$\tau^{\mathrm{L}}_{\bm{p}}=2\pi/\omega_{\bm{p}}.$$ 
The non-dimensional ratio of these time scales can 
characterize  the level of nonlinearity in the nonlinear system:
%
\begin{equation}
{\cal  \epsilon}_{\bm{p}} = 
\frac
{\tau_{\bm{p}}^{\mathrm{L}}} 
{\tau_{\bm{p}}^{\mathrm{NL}}}
= 
\frac
{2\pi \dot n_{\bm{p}}}
{n_{\bm{p}} \omega_{\bm{p}}}
\label{NonlinearRatio}
\end{equation}
%
We refer to (\ref{NonlinearRatio}) as the nonlinearity parameter.  

The nonlinear parameter serves as a low order 
consistency check for the various kinetic equation derivations. An $O(1)$ value of ${\cal  \epsilon}_{\bm{p}}$ implies that the derivation of the kinetic equation is internally inconsistent.  The Boltzman rate represents the net rate of transfer for wavenumber $\bm{p}$ and is an appropriate measure of nonlinearity for smooth, isotropic and homogeneous spectra.  The individual rates of transfer into and out of $\bm{p}$ maybe significantly larger for spectral spikes \citep{M86} and potentially for smooth, homogeneous but anisotropic spectra.  Estimates of the Boltzman rate and ${\cal  \epsilon}_{\bm{p}}$  require integration of Eq.~(\ref{KineticEquation}). In this manuscript such integration 
is performed numerically.   



In this paper we concentrate on four approaches, 
namely \citet{MO75,Voronovich,Zeitlin,LT,LT2}.
We show that on the resonant manifold they 
produce {\em equivalent\/} results. 


Resonant interaction approximation is self-consistent 
for small level of nonlinearities. However, as the 
nonlinearity parameter increase, near-resonant interactions start 
to play a role. 

For realistic estimates the effects of rotation must be included, and this restricts our investigations to 
two approaches that allow inclusion of background rotations. Therefore, we concentrate in more details on the \cite{LT2} and \cite{MO75} representations.


We show that for the  near-resonant interactions, these two approaches returns
qualitatively different predictions for transfer rates. 
This is the main physical result of the
present paper. 
 

There is a multitude of reasons for possible differences.  First and foremost, we view the distinction between Lagrangian, isopycnal and Eulerian coordinates as the most dynamically significant difference.  The use of a Lagrangian coordinate system requires an expansion in powers of small fluid parcel displacements in addition to an assumption of weak nonlinearity, whereas formulations in isopycnal or Eulerian coordinates require only an assumption of weak nonlinearity.  An issue with extant Lagrangian coordinate representations is that the small amplitude assumption represents an unconstrained approximation whose domain of validity {\em vis-a-vis} the weak interaction approximation is not well defined, \citep{M86}.  A subsidiary issue is that the use of a Lagrangian coordinate system places the nonlinearity in the incompressibility constraint, and a single plane wave is not an exact solution of the equations of motion, \citep{S85}.  Similarly, a single plane wave also does not constitute a solution to the isopycnal equations of motion.  In Eulerian coordinates the nonlinearity is advective and a single plane wave is an exact solution of the equations of motion.  On the other hand, it is a robust observational fact that Eulerian frequency spectra at high vertical wavenumber are contaminated by vertical Doppler shifting:  near-inertial frequency energy is Doppler shifted to higher frequency at approximately the same vertical wavelength.  Use of an isopycnal coordinate system considerably reduces this artifact, \citep{SandP91}.  Thus differences in the approaches may represent physical effects rather than technical issues such as the proper implementation of a potential vorticity conservation statement \citep{Zeitlin}.



We emphasize that our intent is to estimate transport rates for various approaches within a common framework and to compare those results.  Our goal is a qualitative physical explanation of the possible reasons for the similarities and differences rather than a quantitative analytical explanation of how those differences arise.  


\section{Various Approaches \label{VariousApproaches}}

In this section we list the approaches that we use. We do so for completion and to transfer everything to a uniform notation.  Our attention is restricted to the hydrostatic balance case, for which 
%
\begin{equation}
|\bm{k}| \ll |m| \ .
\label{hydrostatic}
\end{equation} 
A minor detail is that the linear frequency has different algebraic representations in isopycnal and  Cartesian coordinates. The Cartesian  vertical wavenumber, $k_z$, and the density wavenumber, $m$, are related as $m = - g/(\rho_0 N^2) k_z$ where $g$ is gravity, $\rho$ is density with reference value $\rho_0$, $N$ is the buoyancy (Brunt--V\"{a}is\"{a}l\"{a}) frequency and $f$ is the Coriolis frequency.  
In isopycnal coordinates the dispersion relation is given by,
%
\begin{eqnarray}
\omega({\bm p})  = \sqrt{f^2 + \frac{g^2}{\rho_0^2 N^2} \frac{|\bm{k}|^2}{m^2}}.
 \label{eq:dispersionISO}
\end{eqnarray}
In Cartesian coordinates, 
\begin{eqnarray}
\omega(\bm{p}) = \sqrt{f^2 + N^2 \frac{|\bm{k}|^2}{k_z^2}}~.
\label{eq:dispersionCAR}
\end{eqnarray}
%
In the limit of $f=0$ these dispersion relations assume the 
form  
\begin{equation}
\omega_{\bm{p}} \propto  \frac{ |\bm{k} |}{|m|} \propto \frac{ |\bm{k} |}{|k_z|} 
\label{dispersion}
\end{equation}
%


\subsection{Kenyon and Hasselmann}
The first kinetic equations for wave-wave interactions in a continuously stratified ocean appear in
\citet{K66}, \citet{H66} and \citet{K68}.  \citet{K68} states (without detail) that \citet{K66} and \citet{H66} give numerically similar results.  We have found that \citet{K66} differs from the four approaches examined below on one of the resonant manifolds, but have not pursued the question further.  
It is possible this difference results from a typographical error in \citet{K66}.  We have not rederived this non-Hamiltonian representation and thus exclude it from this study.  

\subsection{  M\"uller and Olbers}



Matrix elements derived in  \citet{O74} are given by
$|{V^{\bm{p}}_{\bm{p}_1, \bm{p}_2}}^{\mathrm{MO}}|^2 = T^{+} / (4\pi)$ and
$|{V^{\bm{p}_1}_{\bm{p}_2, \bm{p}}}^{\mathrm{MO}}|^2 = T^{-} / (4\pi)$.  We extracted $T^{\pm}$ from the Appendix of \citet{MO75}.  In our notation, in the hydrostatic balance approximation, their matrix elements are given by
%

\begin{align}
|{V^{\bm{p}}_{\bm{p}_1, \bm{p}_2}}^{\mathrm{MO}}|^2 =\frac{(N_0^2-f^2)^2}{32 \rho_0} \omega \omega_1 \omega_2
\left|
 \frac{|\bm{k}| |\bm{k}_1| |\bm{k}_2|}{\omega \omega_1 \omega_2 |\bm{p}||\bm{p}_1||\bm{p}_2|}
\right.
\nonumber\\
\left(
 - \frac{\left(-m_1 \frac{\bm{k}_1 \cdot \bm{k}_2 - i f \bm{k}_2 \cdot \bm{k}_1^{\perp}/\omega_1}{k_1^2} + m_2\right) \left(-m_2 \frac{\bm{k}_1 \cdot \bm{k}_2  - i f \bm{k}_1 \cdot \bm{k}_2^{\perp}/\omega_2}{k_2^2} + m_1\right)}{m}
\right.
\nonumber\\
 - \frac{\left(-m_2 \frac{\bm{k}_2 \cdot \bm{k} + i f \bm{k}_2 \cdot \bm{k}^{\perp}/\omega_2}{k_2^2} + m\right)
 \left(-m \frac{\bm{k}_2 \cdot \bm{k}  - i f \bm{k} \cdot \bm{k}_2^{\perp}/\omega}{k^2} + m_2\right)}{m_1}
%
\nonumber\\
\left.\left.
 - \frac{\left(-m \frac{\bm{k} \cdot \bm{k}_1 - i f \bm{k} \cdot \bm{k}_1^{\perp}/\omega}{k^2} + m_1\right)
 \left(-m_1 \frac{\bm{k} \cdot \bm{k}_1 + i f \bm{k}_1 \cdot \bm{k}^{\perp}/\omega_1}{k_1^2} + m\right)}{m_2}
\right)
\right|^2
.
\label{VMO}
\end{align}

Taking a $f=0$ limit  we get:

\begin{eqnarray}
 |{V^{\bm{p}}_{\bm{p}_1, \bm{p}_2}}^{\mathrm{MO}}|^2 \propto
\frac{|\bm{k}||\bm{k}_1||\bm{k}_2|}{|m m_1 m_2|}
 \left(
 - \frac{1}{m}
 \left(-\frac{m_2 \bm{k}_1 \cdot \bm{k}_2}{|\bm{k}_2|^2} + m_1 \right)
 \left(-\frac{m_1 \bm{k}_2 \cdot \bm{k}_1}{|\bm{k}_1|^2} + m_2 \right)
\right.
\nonumber
\\
\left.
 + \frac{1}{m_1}
 \left(\frac{m_2 \bm{k} \cdot \bm{k}_2}{|\bm{k}_2|^2} - m \right)
 \left(-\frac{m \bm{k}_2 \cdot \bm{k}}{|\bm{k}|^2} + m_2 \right)
 + \frac{1}{m_2}
 \left(-\frac{m \bm{k}_1 \cdot \bm{k}}{|\bm{k}|^2} + m_1 \right)
 \left(\frac{m_1 \bm{k} \cdot \bm{k}_1}{|\bm{k}_1|^2} - m \right)
 \right)^2
\end{eqnarray}

\subsection{Pelinovsky and Raevsky}
An important  paper on internal waves is \citet{PR77}.  Clebsh variables are used to obtain the
interaction matrix elements for both constant stratification rates, $N=\mathrm{const.}$, 
and arbitrary buoyancy profiles, $N=N(z)$.
Not much details are given, but there are some similarities in appearance with 
\citet{Voronovich}. The most significant result is the identification of a scale invariant (non-rotating, hydrostatic) stationary state.  It is stated in the paper that their matrix elements are
equivalent to those derived in their citation [11], which is ~\citet{B75}.
Because \citet{B75} and \cite{PR77} are in Russian and not generally available, we 
refrain from including them in this comparison.  


\subsection {Voronovich}
Voronovich used Clebsh-like variables to derive the Hamiltonian for incompressible stratified flows in the ocean. It is probably the first canonical Hamiltonian structure derived for such kind of flows.  A detailed explanation of Voronovich's method appears in section 7.1 of 
the textbook \citet{Mir}
It is a straightforward task to write down the kinetic equation associated with this Hamiltonian structure. 

We formulate the matrix elements for Voronovich's Hamiltonian using his formula (A.1).  This formula is derived for general boundary conditions. To compare with other matrix elements of this paper, we assume a constant stratification profile and Fourier basis as the vertical structure function $\phi(z)$. That allows us to solve for the matrix elements defined via Eq.~(11) and above it in his paper.
Then the convolutions of the basis functions give delta-functions in vertical wavenumbers.
Vornovich's equation (A.1) transforms into:
%
 \begin{eqnarray}
  |{V^{\bm{p}}_{\bm{p}_1, \bm{p}_2}}^{\mathrm{V}}|^2 \propto
 \frac{|\bm{k}||\bm{k}_1||\bm{k}_2|}{|m m_1 m_2|}
 \left(
 - m
 \left(
 \frac{1}{|\bm{k}| |m|}
 \left(\frac{\bm{k} \cdot \bm{k}_1 |m_1|}{|\bm{k}_1|} + \frac{\bm{k}
 \cdot \bm{k}_2 |m_2|}{|\bm{k}_2|} \right)
 + \frac{\omega_1 + \omega_2 - \omega}{\omega}
 \right)
 \right.
 \nonumber\\
 \left.
 + m_1
 \left(
 \frac{1}{|\bm{k}_1| |m_1|}
 \left(\frac{\bm{k} \cdot \bm{k}_1 |m|}{|\bm{k}|} + \frac{\bm{k}_1
 \cdot \bm{k}_2 |m_2|}{|\bm{k}_2|} \right)
 - \frac{\omega_1 + \omega_2 - \omega}{\omega_1}
 \right)
 \right.
 \nonumber\\
 \left.
 + m_2
 \left(
 \frac{1}{|\bm{k}_2| |m_2|}
 \left(\frac{\bm{k} \cdot \bm{k}_2 |m|}{|\bm{k}|} + \frac{\bm{k}_2
 \cdot \bm{k}_1 |m_1|}{|\bm{k}_1|} \right)
 - \frac{\omega_1 + \omega_2 - \omega}{\omega_2}
 \right)
 \right)^2 . \nonumber \\
 \label{eq:Voronovich}
 \end{eqnarray}

Note that Eq.~(\ref{eq:Voronovich}) shares structural similarities with the interaction matrix elements in {\em isopycnal\/} coordinates, Eq.~(\ref{Hamiltonian}) below.



\subsection{Milder}
An alternative Hamiltonian description was developed in
\citet{Milder}, in isopycnal coordinates without assuming a hydrostatic balance.
The resulting Hamiltonian is an iterative expansion in powers of a
small parameter, similar to the case of surface gravity waves. In
principle, that approach may also be used to calculate wave-wave
interaction amplitudes.  Since those calculations were not done
in \citet{Milder}, we do not pursue the comparison further.

\subsection{Caillol and Zeitlin}
%
A non-Hamiltonian kinetic equation for internal waves was derived in
\citet{Zeitlin}, Eq.~(61). To make it appear equivalent to more
traditional form of kinetic equation, as in \citet{ZLF}, we make
a change of variables $\bm{l}\to -\bm{l}$ in the second line, and
$\bm{k}\to -\bm{k}$ in the third line of (61) of \citet{Zeitlin}.  If
we further assume that all spectra are symmetric, $n(-\bm{p}) =
n(\bm{p})$, then the kinetic equation assumes traditional form, as in Eq.~(\ref{KineticEquation}), see
\citet{MO75,ZLF,LT,LT2}.

The matrix elements according to \citet{Zeitlin} are shown as
$X_{k,l,p}$ and $Y_{k,l,p}^{\pm}$ in Eqs.~(62) and (63), where
$|{V^{\bm{p}}_{\bm{p}_1, \bm{p}_2}}^{\mathrm{CZ}}|^2 = X_{\bm{p}_1,\bm{p}_2,\bm{p}}$ and
$|{V^{\bm{p}_1}_{\bm{p}_2, \bm{p}}}^{\mathrm{CZ}}|^2 = Y_{\bm{p}_1,-\bm{p}_2,\bm{p}}^{+}$.
%
In our notation it reads
%
\begin{eqnarray}
 |{V^{\bm{p}}_{\bm{p}_1, \bm{p}_2}}^{\mathrm{CZ}}|^2 \propto
(|\bm{k}| \mathrm{sgn}(m) + |\bm{k}_1| \mathrm{sgn}(m_1) + |\bm{k}_2| \mathrm{sgn}(m_2))^2
\frac{(m^2 - m_1 m_2)^2}{|m| |m_1| |m_2| |\bm{k}||\bm{k}_1||\bm{k}_2|}
\nonumber
\\ 
\times\left(
\frac{|\bm{k}|^2 - |\bm{k}_1| \mathrm{sgn}(m_1) |\bm{k}_2| \mathrm{sgn}(m_2)}{m^2 - m_1 m_2} m
- \frac{|\bm{k}_1|^2}{m_1}
- \frac{|\bm{k}_2|^2}{m_2}
\right)^2 \, \nonumber \\.\label{eq:VCZ}
\end{eqnarray}

\subsection{Isopycnal Hamiltonian}

Finally, in \citet{LT2} the following wave-wave interaction matrix
element was derived based on a canonical Hamiltonian formulation in isopycnal coordinates:
%

\begin{align}
 |{V^0_{1,2}}   ^{\mathrm{H}}
|^2 = \frac{N^2}{32 g}
\left(
\left(
\frac{k \bm{k}_1 \cdot \bm{k}_2}{k_1 k_2} \sqrt{\frac{\omega_1 \omega_2}{\omega}}
+ \frac{k_1 \bm{k}_2 \cdot \bm{k}}{k_2 k} \sqrt{\frac{\omega_2 \omega}{\omega_1}}
+ \frac{k_2 \bm{k} \cdot \bm{k}_1}{k k_1} \sqrt{\frac{\omega \omega_1}{\omega_2}}
\right.
\right.
\nonumber\\
\left.
\left.
+ \frac{f^2}{\sqrt{\omega \omega_1 \omega_2}}
\frac{k_1^2 \bm{k}_2 \cdot \bm{k} - k_2^2 \bm{k} \cdot \bm{k}_1 - k^2 \bm{k}_1 \cdot \bm{k}_2}{k k_1 k_2}
\right)^2
\right.
\nonumber\\
\left.
+
\left(
f \frac{\bm{k}_1 \cdot \bm{k}_2^{\perp}}{k k_1 k_2}
 \left(\sqrt{\frac{\omega}{\omega_1 \omega_2}} (k_1^2 - k_2^2)
 - \sqrt{\frac{\omega_1}{\omega_2 \omega}}  (k_2^2-k^2)
 - \sqrt{\frac{\omega_2}{\omega \omega_1}} (k^2-k_1^2)\right)
\right)^2
\right)~\ .\nonumber \\ \label{LTV}
\end{align}
%
\citet{LT} is a rotationless limit of \citet{LT2}.
Taking the $f\to 0$ limit, the \cite{LT2} reduces to \cite{LT}, and (\ref{LTV}) reduces to 
%
\begin{eqnarray}
|{V^{\bm{p}}_{\bm{p}_1, \bm{p}_2}}^{\mathrm{H}}|^2 \propto
\frac{1}{|\bm{k}||\bm{k}_1||\bm{k}_2|} \left(
 |\bm{k}| \bm{k}_1 \cdot \bm{k}_2 \sqrt{\left|\frac{m}{m_1 m_2}\right|}
+ |\bm{k}_1| \bm{k}_2 \cdot \bm{k} \sqrt{\left|\frac{m_1}{m_2 m}\right|}
+ |\bm{k}_2| \bm{k} \cdot \bm{k}_1 \sqrt{\left|\frac{m_2}{m m_1}\right|}
 \right)^2 .\nonumber\\
\label{Hamiltonian}
\end{eqnarray}
%

Observe that in this form, these equations share structural similarities with 
Eq.~(\ref{eq:Voronovich}).  

In this section we gave brief review of the various approaches that were 
developed for describing wave-wave interactions of internal waves in the 
ocean. While this review is necessarily brief, this is the first time 
all these papers are cited together by a single manuscript. 

\section{Resonant wave-wave interactions \label{ResonantInteractions}}

How one can compare the function of two vectors $\bm{p}_1$ and $\bm{p}_2$, and their
sum or difference? First one realizes that out of 6 components of $\bm{p}_1$ and $\bm{p}_2$, only relative angles between wavevectors enter into the equation for matrix elements. That is because the matrix elements depend on the inner products of wavevectors. The overall horizontal orientation of the wavevectors does not matter:  relative angles can be determined from a triangle inequality and the magnitudes of the horizontal wavevectors $\bm{k}$, $\bm{k}_1$ and $\bm{k}_2$. Thus the only needed components are $|\bm{k}|$, $|\bm{k}_1|$, $|\bm{k}_2|$, $m$ and $m_1$ ($m_2$ is computed from $m$ and $m_1$).  Further note that in the $f=0$ and hydrostatic limit, all matrix elements become scale invariant functions. That is to say that it is sufficient to choose an arbitrary scalar value for $|\bm{k}|$, and $m$, since only $|\bm{k}_1|/|\bm{k}|$, $|\bm{k}_2|/|\bm{k}|$ and  $m_1/m$ enter the expressions for matrix elements.  We make the particular
(arbitrary) choice that $|\bm{k}|=m=1$ for the purpose of numerical evaluation, and thus the only independent variables to consider are
$|\bm{k}_1|$, $|\bm{k}_2|$ and $m_1$. Finally, $m_1$ is determined from the resonance conditions, as
explained in the next subsection below. As a result, we are left with a matrix element as a function of only two parameters, $k_1$ and $k_2$.  This allows us to easily compare the values of matrix elements on the resonant manifold. 

\subsection{Reduction to the Resonant Manifold}
When confined to the traditional form of the kinetic equation, wave-wave
interactions (scattering) are constrained to the resonant manifolds defined by
\begin{eqnarray}
a)~~
\begin{cases}
 \bm{p} = \bm{p}_1 + \bm{p}_2 \\
 \omega = \omega_1 + \omega_2
\end{cases}
b)~~
\begin{cases}
 \bm{p}_1 = \bm{p}_2 + \bm{p} \\
 \omega_1 = \omega_2 + \omega
\end{cases}
c)~~
\begin{cases}
 \bm{p}_2 = \bm{p} + \bm{p}_1 \\
 \omega_2 = \omega + \omega_1
\end{cases}.
\label{RESONANCES}
\end{eqnarray}
%
To compare matrix elements on the resonant manifold we are going to
use the above resonant conditions and the internal-wave dispersion relation (\ref{dispersion}).
To determine vertical components $m_1$ and $m_2$ of the interacting
wavevectors, one has to solve the resulting quadratic equations. Without restricting generality we choose  $m>0$.  There are two solutions for $m_1$ and $m_2$ given below for each of the three resonance types described above.

Resonances of type (\ref{RESONANCES}a) give
%	
\begin{subequations}
	 \allowdisplaybreaks
 \begin{align}
&
\begin{cases}
 m_1 = \frac{m}{2 |\bm{k}|} \left(|\bm{k}| + |\bm{k}_1| + |\bm{k}_2| + \sqrt{(|\bm{k}| + |\bm{k}_1| + |\bm{k}_2|)^2 - 4 |\bm{k}| |\bm{k}_1|}\right)
 \\
 m_2 = m - m_1.
\end{cases}
,
\label{eq:sol1}
\\
&
\begin{cases}
m_1 = \frac{m}{2|\bm{k}|} \left(|\bm{k}| - |\bm{k}_1| - |\bm{k}_2| - \sqrt{(|\bm{k}| - |\bm{k}_1| - |\bm{k}_2|)^2 + 4 |\bm{k}| |\bm{k}_1|}\right)
\\
m_2 = m - m_1.
\end{cases}
,
\label{eq:sol2}
 \end{align}
	\end{subequations}
Note that because of the symmetry, (\ref{eq:sol1}) translates to (\ref{eq:sol2}) if wavenumbers $1$ and $2$ are exchanged.
 
Resonances of type (\ref{RESONANCES}b) give
\begin{subequations}
\allowdisplaybreaks
\begin{align}
&
\begin{cases}
m_2 = - \frac{m}{2 |\bm{k}|} \left(|\bm{k}| - |\bm{k}_1| - |\bm{k}_2| + \sqrt{(|\bm{k}| - |\bm{k}_1| - |\bm{k}_2|)^2 + 4 |\bm{k}| |\bm{k}_2|}\right)
\\
m_1 = m + m_2.
\end{cases}
,
\label{eq:sol3}
\\
&
\begin{cases}
m_2 = - \frac{m}{2|\bm{k}|} \left(|\bm{k}| + |\bm{k}_1| - |\bm{k}_2| + \sqrt{(|\bm{k}| + |\bm{k}_1| - |\bm{k}_2|)^2 + 4 |\bm{k}| |\bm{k}_2|}\right)
\\
m_1 = m + m_2.
\end{cases}
,
\label{eq:sol4}
\end{align}
\end{subequations}

Resonances of type (\ref{RESONANCES}c) give
\begin{subequations}
\allowdisplaybreaks
\begin{align}
&
\begin{cases}
m_1 = - \frac{m}{2|\bm{k}|} \left(|\bm{k}| - |\bm{k}_1| - |\bm{k}_2| + \sqrt{(|\bm{k}| - |\bm{k}_1| - |\bm{k}_2|)^2 + 4 |\bm{k}| |\bm{k}_1|}\right)
\\
m_2 = m + m_1.
\end{cases}
,
\label{eq:sol5}
\\
&
\begin{cases}
m_1 = - \frac{m}{2|\bm{k}|} \left(|\bm{k}| - |\bm{k}_1| + |\bm{k}_2| + \sqrt{(|\bm{k}| - |\bm{k}_1| + |\bm{k}_2|)^2 + 4 |\bm{k}| |\bm{k}_1|}\right)
\\
m_2 = m + m_1.
\end{cases}
.
\label{eq:sol6}
\end{align}
\end{subequations}
%
Because of the symmetries of the problem, (\ref{eq:sol3}) is equivalent to
(\ref{eq:sol5}), and (\ref{eq:sol4}) is equivalent to (\ref{eq:sol6})
if wavenumbers $1$ and $2$ are exchanged.


\subsection{Comparison of matrix elements}

As explained above, we assume $f=0$ and hydrostatic balance. Such a
choice makes the matrix elements to be scale-invariant functions
that depend only upon $|\bm{k}_1|$ and $|\bm{k}_2|$.
As a consequence of the triangle inequality we need to consider matrix elements only within a ``kinematic box'' defined by $$||\bm{k}_1| - |\bm{k}_2|| < |\bm{k}| < |\bm{k}_1| + |\bm{k}_2|.$$
The matrix elements will have different values depending on the dimensions so that isopycnal and Eulerian approaches will give different values (\ref{eq:dispersionISO})-(\ref{eq:dispersionCAR}). To address this issue in the simplest possible way, we multiply each matrix element by a dimensional number chosen so that all
 matrix elements are equivalent for some specific wavevector. In particular, we choose the scaling constant so 
that $|V(|\bm{k}_1|=1,|\bm{k}_2|=1)|^2=1$. This allows a transparent comparison without worrying about dimensional 
differences between various formulations.



\subsubsection{Resonances of the ``sum'' type (\ref{RESONANCES}a)}

Figure~\ref{FIGRESONANTa}
presents the values of the matrix element 
$|{V^{\bm{p}}_{\bm{p}_1, \bm{p}_2}}_{\mathrm{(\ref{eq:sol2})}}|^2$ on the resonant sub-manifold given explicitly by
(\ref{eq:sol2}).
All approaches
give equivalent
results. This is confirmed by plotting the relative ratio between
these approaches, and it is given by numerical noise (not shown).
The solution (\ref{eq:sol1}) gives the same matrix elements
but with $|\bm{k}_1|$ and $|\bm{k}_2|$ exchanged
owing to their symmetries.

\subsubsection{Resonances of the ``difference'' type (\ref{RESONANCES}b) and
 (\ref{RESONANCES}c)}

We then turn our attention to resonances of ``difference'' type (\ref{RESONANCES}b) for which 
(\ref{RESONANCES}c) could be obtained by symmetrical exchange of the indices.
All the matrix elements 
$|{V^{\bm{p}_1}_{\bm{p}_2, \bm{p}}}_{\mathrm{(\ref{eq:sol3})}}|^2$
on the resonant sub-manifold (\ref{eq:sol3}), 
are shown in Fig.~\ref{FIGRESONANTb}.
All the matrix elements are equivalent.  The relative differences between different approaches are
given by numerical noise (not shown).
Finally, $|{V^{\bm{p}_1}_{\bm{p}_2, \bm{p}}}_{\mathrm{(\ref{eq:sol4})}}|^2$ on the 
resonant sub-manifold (\ref{eq:sol4}) are shown in Fig.~\ref{FIGRESONANTc}.
Again, all the matrix elements
are equivalent.

The solutions (\ref{eq:sol5}) and (\ref{eq:sol6}) give the same matrix elements but with $|\bm{k}_1|$ and $|\bm{k}_2|$ exchanged
as the solutions (\ref{eq:sol3}) and (\ref{eq:sol4})
owing to their symmetries.



\subsubsection{Special triads}
Three simple interaction mechanisms are identified by \citet{MB77} in the limit of an extreme scale separation.  In this subsection we look in closer detail at these special limiting triads to confirm that all matrix elements are indeed asymptotically consistent.  The limiting cases are:

\begin{itemize}
\item
the vertical backscattering of a high-frequency wave by a low frequency wave of twice the vertical wavenumber into a second high-frequency wave of oppositely signed vertical wavenumber. This type of scattering is called elastic scattering (ES). The solution (\ref{eq:sol1}) in the limit  $|\bm{k}_1| \to 0$ corresponds to this type of special triad.
\item
The scattering of a high-frequency wave by a low-frequency, small-wavenumber wave into a second, nearly identical, high-frequency large-wavenumber wave. This type of scattering is called  induced diffusion (ID). The solution (\ref{eq:sol2}) in the limit that $|\bm{k}_1| \to 0$ corresponds to this type of special triad.
\item
The decay of a low wavenumber wave into two high vertical wavenumber waves of approximately one-half the frequency. This is called  parametric subharmonic instability (PSI).  The solution (\ref{eq:sol3}) in the limit that $|\bm{k}_1| \to 0$ corresponds to this type of triad.
\end{itemize}

To study the detailed behavior of the matrix elements in the special triad cases, we choose to present the matrix elements along a straight line defined by
$$
(|\bm{k}_1|, |\bm{k}_2|) = (\epsilon, \epsilon/3+1)  |\bm{k}|.
$$
This line originates from the corner of the kinematic box in Figs.~\ref{FIGRESONANTa}--\ref{FIGRESONANTc} at $(|\bm{k}_1|, |\bm{k}_2|) = (0,|\bm{k}|) $ and has a slope of 1/3.  The slope of this line is arbitrary.  We could have taken $\epsilon/4$ or $\epsilon/2$.  The matrix elements here are shown as functions of $\epsilon$
in Fig.~\ref{FigureThree}.  We see that
all four approaches are again {\em equivalent}  on the resonant manifold for the 
case of special triads.

In this section we demonstrated that all four approaches we
considered produce {\it equivalent} results on the resonant
manifold in the absence of background rotation. This statement is not 
trivial, given the different assumptions and coordinate systems that
have been used for the various kinetic equation derivations. 


\section{Smearing of the resonance manifold and near-resonant Interactions\label{OffResonant}}


\subsection{Nonlinear frequency renormalization as a result of nonlinear wave-wave interactions}
Above we have compared the values of matrix element on the {\em
resonant manifold}.  The resonant interaction approximation is a
mathematical simplification which reduces the complexity of the
problem. In this subsection we examine transfers including near
resonant interactions.  Our interest in near-resonant interactions has
significant physical motivations. For example a major unresolved issue
is the importance of Doppler shifting \citep{Polzin2004a}.  Of particular
interest here is the variable effects of Doppler shifting in different
coordinate systems.  The resonant interaction approximation assumes,
perforce, an expansion in terms of a non-advected wavefield, with
dispersion relation given by Eq.~(\ref{eq:dispersionISO}) or
Eq.~(\ref{eq:dispersionCAR}).  In the limit of extreme time scale
separation between high frequency waves and a low frequency
background, one is tempted to replace the non-advected frequency by its
Doppler shifted intrinsic frequency counterpart, $\omega \rightarrow
\omega - \bm{k} \cdot \overline{\bm{u}}$, in which $\omega$ and $\bm{k}$
are the frequency and wavevector of the high frequency wave and
$\overline{\bm{u}}$ is the velocity field of the low frequency
wavefield.  This is the genesis of the eikonal approach \citep{M86} to
internal wave-wave interactions.  
Then the resonant  approximation is self-consistent
for small values of nonlinearities. 
Indeed, change in the wave amplitude will be
small, and the Doppler shift cancels from the frequency delta function. 
Yet, as nonlinearity increases, the 
near-resonant interactions become more and more pronounced, consequently
the  issue of Doppler shifting more and more important. 
Furthermore, near-resonant
interactions play a major role in numerical simulations on a discrete
grids~\citep{LvovNazarenkoPokorni}, for time evolution of discrete systems
~\citep{Gersh2007}, in acoustic turbulence ~\citep{LLNZ}, surface
gravity waves~\citep{JansenXXX,yuen_lake}, and internal waves
\citep{2006JFM...568..273V,Shrira}.


To take into account the effects of near-resonant interactions 
self-consistently, the energy conserving delta-functions in Eq.~(\ref{KineticEquation}), $\delta({\omega_{\bm{p}} -\omega_{{\bm{p}_1}} - \omega_{{\bm{p}_2}} })$, need to be ``broadened''.  The physical motivation for this
broadening is the following: when the resonant kinetic equation is
derived, it is assumed that the amplitude of each plane wave is
constant in time, or, in other words, that the lifetime of single
plane wave is infinite. Resulting kinetic equation, nevertheless, predicts
that the amplitude of the wave do change. Consequently the   
wave lifetime is finite.  For small level of nonlinearity this 
distinction is not significant, and resonant kinetic equation constitutes a 
self-consistent description. For larger values of nonliterary this is no longer the 
case, and the wave lifetime is finite and amplitude changes need to be taken into account.  
Consequently interactions may not be strictly resonant.
This statement also follows from the Fourier uncertainty principle. 
In other words, the
waves with varying amplitude can not be represented by a single Fourier component.  
This effect is larger for stronger level of nonlinearity parameter. 


Derivation of the kinetic equation with a broadened delta
function is given in details in \citep{LLNZ}, and is not going to be 
repeated  here. The result is that
%
\begin{eqnarray}
\frac{d n_{\bm{p}}}{dt} = 4 \int
 |V_{\bm{p}_1,\bm{p}_2}^{\bm{p}}|^2 \, f_{p12} \,
\delta_{{\bm{p} - \bm{p}_1-\bm{p}_2}} \, {\cal L}({\omega_{\bm{p}}
-\omega_{{\bm{p}_1}}-\omega_{{\bm{p}_2}}})
d \bm{p}_{12}
\nonumber \\
-4\int
 \, |V_{\bm{p}_2,\bm{p}}^{\bm{p}_1}|^2\, f_{12p}\, \delta_{{\bm{p}_1 - \bm{p}_2-\bm{p}}} \,
  {\cal L} ({{\omega_{\bm{p}_1} -\omega_{\bm{p}_2}-\omega_{\bm{p}}}})
\, d \bm{p}_{12}
\nonumber \\
-4\int
 \, |V_{\bm{p},\bm{p}_1}^{\bm{p}_2}|^2\, f_{2p1}\, \delta_{{\bm{p}_2 - \bm{p}-\bm{p}_1}} \,
  {\cal L}({{\omega_{\bm{p}_2} -\omega_{\bm{p}}-\omega_{\bm{p}_1}}})
\, d \bm{p}_{12}
 \, ,\nonumber\\
\label{KineticEquationBroadened}
\end{eqnarray}
%
Here ${\cal L}$ is defined as
\begin{equation}
{\cal{L}}(\Delta\omega)  = 
\frac{\Gamma_{k12}}{(\Delta\omega)^2 + \Gamma_{k12}^2},
\label{scriptyL}
\end{equation}
%
where $\Gamma_{k12}$ is the total broadening of each particular resonance, and is given below.
If the nonlinear frequency renormalization tends to zero, i.e. $\Gamma_{k12} \to 0$, ${\cal L}$ reduces to the delta function: 
%
$$\lim\limits_{\Gamma_{k12}\to 0} {\cal{L}}(\Delta\omega)  =\pi 
\delta(\Delta\omega).$$
Consequently, in the limit resonant interactions (i.e. no broadening) 
(\ref{KineticEquationBroadened}) reduces to Eq.~(\ref{KineticEquation}) .
%
We have shown in \cite{LLNZ} that the broadening in Eq.~(\ref{scriptyL}) is  given by 
%
\begin{equation}
\Gamma_{k12}=\gamma_{\bm{p}}+\gamma_{\bm{p}_1}+\gamma_{\bm{p}_2}.
\label{Gammak12}
\end{equation}
%
It means that the total resonance broadening is the sum of individual
frequency broadening, and can be thus seen as the ``triad
interaction'' frequency. 

The single frequency renormalization can be
calculated {\em self-consistently} from
%
\begin{eqnarray}
\gamma_{\bm{p}} = 4 \int
 |V_{\bm{p}_1,\bm{p}_2}^{\bm{p}}|^2 \, (n_{{\bm{p}}_1}+
n_{{\bm{p}}_2}) \,
\delta_{{\bm{p} - \bm{p}_1-\bm{p}_2}} \, {\cal L}({\omega_{\bm{p}}
-\omega_{{\bm{p}_1}}-\omega_{{\bm{p}_2}}})
d \bm{p}_{12}
\nonumber \\
-4\int
 \, |V_{\bm{p}_2,\bm{p}}^{\bm{p}_1}|^2\, (n_{{\bm{p}}_2}  - n_{{\bm{p}}_1}) 
\, \delta_{{\bm{p}_1 - \bm{p}_2-\bm{p}}} \,
  {\cal L} ({{\omega_{\bm{p}_1} -\omega_{\bm{p}_2}-\omega_{\bm{p}}}})
\, d \bm{p}_{12}
\nonumber \\
-4\int
 \, |V_{\bm{p},\bm{p}_1}^{\bm{p}_2}|^2\, (n_{{\bm{p}}_1}- n_{{\bm{p}}_2})
 \, \delta_{{\bm{p}_2 - \bm{p}-\bm{p}_1}} \,
  {\cal L}({{\omega_{\bm{p}_2} -\omega_{\bm{p}}-\omega_{\bm{p}_1}}})
\, d \bm{p}_{12}
 \ . \nonumber\\
\label{Gamma}
\end{eqnarray}
%
The interpretation of this formula is the following: nonlinear 
wave-wave interactions lead to the change of wave amplitude, which 
in turn makes the lifetime of the waves to be finite. This, in turn, makes 
the interactions to be near-resonant. 

A self-consistent estimate of $\gamma_{\bm{p}} $ requires the
iterative solution of (\ref{KineticEquationBroadened}) and
(\ref{Gamma}) over the entire field: the width of the resonance
(\ref{Gamma}) depends on the lifetime of an individual wave [from
(\ref{KineticEquationBroadened})], which in turn depends on the width
of the resonance (\ref{Gammak12}).  This numerically intensive
computation is beyond the scope of this manuscript.  Instead, we make
the uncontrolled approximation that:
%
\begin{equation}
\gamma_{\bm{p}} = \delta \omega_{\bm{p}}.
\label{GammaFraction}
\end{equation}

We note that this 
choice is made for illustration purposes only, we certainly do not claim that 
it represents a self consistent choice. Below, we will take $\delta$ to be $10^{-2}$ and $10^{-3}$. These values are rather small, therefore we
remain in the closest proximity to the resonant interactions.


\subsection{Characteristic nonlinear time scale of the Garrett and Munk Spectrum}

Estimates of near-resonant transfers are obtained by assuming horizontal isotropy integrating (\ref{KineticEquationBroadened}) over horizontal azimuth:
\begin{eqnarray}
\frac{\partial n_{\bm{p}}}{\partial t} = 4\pi \int
\frac{k_1 k_2}{S_{p12}} |V_{\bm{p}_1,\bm{p}_2}^{\bm{p}}|^2 \, f_{p12} \,
\delta_{{\bm{p} - \bm{p}_1-\bm{p}_2}} \, {\cal L}({\omega_{\bm{p}}
-\omega_{{\bm{p}_1}}-\omega_{{\bm{p}_2}}})
dk_{12} dm_1
\nonumber \\
-4\pi \int
 \, \frac{k_1 k_2}{S_{12p}} |V_{\bm{p}_2,\bm{p}}^{\bm{p}_1}|^2\,
f_{12p}\, \delta_{{\bm{p}_1 - \bm{p}_2-\bm{p}}} \,
 {\cal L} ({{\omega_{\bm{p}_1} -\omega_{\bm{p}_2}-\omega_{\bm{p}}}})
\, dk_{12} dm_1
\nonumber \\
-4\pi \int
 \, \frac{k_1 k_2}{S_{2p1}} |V_{\bm{p},\bm{p}_1}^{\bm{p}_2}|^2\,
f_{2p1}\, \delta_{{\bm{p}_2 - \bm{p}-\bm{p}_1}} \,
 {\cal L}({{\omega_{\bm{p}_2} -\omega_{\bm{p}}-\omega_{\bm{p}_1}}})
\, dk_{12} dm_1
\, ,
\label{IntKineticEquationBroadened}
\end{eqnarray}
where $S_{p12}$ is the area of the triangle $\bm{k} = \bm{k}_1 + \bm{k}_2$.
We numerically integrated (\ref{IntKineticEquationBroadened})
for $\bm{p}$'s which have frequencies from $f$ to $N$ [specifically (33/32, 17/16, 9/8, 5/4, 3/2, 2, 4, ...)$f$] 
and vertical wavenumbers from $4\pi/(2b)$ to $200\pi/(2b)$ ([2,4,6, ... 98] $\pi/b$) .
The limits of integration are restricted by
horizontal wavenumbers from $2\pi/10^5$ to $2\pi/10$ meters$^{-1}$,
vertical wavenumbers from $2\pi/(2b)$ to $2\pi/10$ meters$^{-1}$,
and frequencies from $f$ to $N$.  
The integrals over $k_1$ and $k_2$
are obtained in the kinematic box in $k_1-k_2$ space.
The grids in the {$k_1-k_2$} domain have $2^{17}$ points
that  are distributed heavily around the corner of the kinematic box.
The integral over $m_1$ is obtained with $2^{13}$ grid points,
which are also distributed heavily
for the small vertical wavenumbers whose absolute values are less than $5m$, where $m$ is the vertical wavenumber. 

Below we calculate the nonlinear time scale (\ref{NonlinearTime}) and nonlinearity parameter (\ref{NonlinearRatio}).  To calculate this parameter, we need to choose a 
form of spectral energy density of internal waves. We 
we utilize the Garrett and Munk spectrum as an agreed-upon representation of the internal waves:  
%
\begin{equation}
E(\omega,m) = 
\frac{4 f}{\pi^2 m_{\ast}} 
E_0 \frac{1}{1+(\frac{m}{m_{\ast}})^2 } 
\frac{1}{\omega \sqrt{\omega^2-f^2} }  \, .
\label{GM}
\end{equation}
%
Here the reference wavenumber is given by 
%
\begin{equation} m_{\ast} = \pi j_{\ast} / b,\label{Jstar}\end{equation} 
%
in which the variable $j$ represents the mode number of an ocean with
an exponential buoyancy frequency profile having a scale height of $b$.  

  
We choose the following set of parameters:
\begin{itemize}
\item $b$ = 1300 m in the GM model
\item The total energy is set as: 
\begin{equation}
E_0 = 30 \times 10^{-4} {\rm ~ m}^2 {\rm ~ s}^{-2} . \nonumber
\end{equation}
\item Inertial frequency is given by $f=10^{-4}$rad/sec, and buoyancy frequency is given by $N_0=5 \times 10^{-3}$rad/sec. 
\item The reference density is taken to be 
$\rho_0=10^{3}$kg/m$^{3}$.
\end{itemize}

We then calculate the nonlinearity parameter (\ref{NonlinearRatio}) and the nonlinear time scale (\ref{NonlinearTime}).
To do so we substitute the Garrett and Munk
spectrum (\ref{GM}) into the kinetic equation with broadening
(\ref{KineticEquationBroadened}). For matrix elements we use \cite{MO75}, Eq. (\ref{VMO}), and 
\cite{LT2}, Eq. (\ref{LTV}).  We also use the dispersion relation of internal
waves, (\ref{eq:dispersionISO}) for the isopycnal Hamiltonian, and
(\ref{eq:dispersionCAR}) for Lagrangian coordinates.  We use two
values of $\delta$ in (\ref{GammaFraction}): $\delta=10^{-2}$ and
$\delta=10^{-3}$. 
We therefore make four calculations:
\begin{itemize}
\item {Run I}~\cite{LT2} with $\delta=10^{-3}$
\item {Run II}~\cite{MO75} with $\delta=10^{-3}$
\item {Run III}~\cite{LT2} with $\delta=10^{-2}$
\item {Run IV}~\cite{MO75} with $\delta=10^{-2}$
\end{itemize}
Results appear in Figs.~\ref{NonlinearityParameter} and \ref{NonlinearTimeFigure}.  

For Run 1 the nonlinearity parameter is uniformly small, smaller than $10^{-1}$.  
Such value of the nonlinearity parameter indicates 
that the kinetic equation is a self-consistent approach for the Garrett and Munk 
Spectrum. Increasing values of the nonlinearity parameter are noted with increasing vertical wavenumbers.  
This is consistent with intuition that we have about such systems.
The nonlinear time scale is of the order of one hundred wave periods at low vertical wavenumber and of order ten wave periods at high vertical wavenumber.  We also define a ``zero curve'' - It is the locus of
wavenumber-frequency where the nonlinearity parameter and time-derivative of waveaction is exactly
zero. 
The zero curve clearly delineates a pattern of energy gain for frequencies $f < \omega < 2f$, energy loss for frequencies $2f < \omega < 5 f$ and energy gain for frequencies $5f < \omega < N$.  This seems to be  a characteristic pattern that appears in our calculations. Note that the zero curves are nearly independent of vertical wavenumber.  

The \cite{MO75}, matrix element~(\ref{VMO}), {Run II} ($\delta=10^{-3}$) results are qualitatively similar to Run I.  Factor of 2-3 faster
decorrelation times and levels of nonlinearity are noted in the high-frequency and high-wavenumber part of the spectrum.  

Therefore we conclude that when near-resonant interactions are
included, the transfer rates are representation dependent. Furthermore, 
Lagrangian approaches predict higher level of nonlinearity. 


To investigate in more details results of near-resonant interactions, we 
perform numerical calculations for $\delta=10^{-2}$. Results for the canonical Hamiltonian formulation in isopycnal coordinates {Run III} are nearly identical to those with $\delta=10^{-3}$.  Results for the Lagrangian coordinate representation are both {\em quantitatively} and {\em qualitatively} different.  
The Lagrangian coordinate formulation~(\ref{VMO}) now predicts $O(1)$ nonlinearity for high frequencies, while the isopycnal coordinate formulation still returns 
nonlinearity parameter and much slower decorrelation time estimates.  The zero curves for the Lagrangian coordinate representation are no longer simple functions of frequency at this higher level of nonlinearity.  The zero curves in the isopycnal coordinate system are relatively independent of $\delta$.  


To investigate the differences between approaches and the sensitivity of our 
results to the value of $\delta$, in more details, we plot in Fig.~\ref{Differences} the differences of the nonlinearity parameter for
these runs. In particular, we calculate the differences between
{Run I} and {Run II},
{Run I} and {Run III}, and finally between
{Run II} and {Run IV}.
Differences associated with increased resonance broadening are
minimal, $10^{-3}$ or smaller, for the isopycnal Hamiltonian.  As the
nonlinearity parameter estimates are representation dependent,
differences between isopycnal coordinate and Lagrangian coordinate
representations are much larger and increase with increasing $\delta$.

We have found that transports for the canonical Hamiltonian
representation are not too sensitive to near-resonant interactions.  We
have also found in Section \ref{ResonantInteractions} that all approaches are
equivalent on the resonant manifold. We therefore conclude that all
approaches will converge to Hamiltonian one as delta decreases.
We have not undertaken such calculations as such small
values of $\delta$ would require significant modifications to our
numerical algorithm.

Note that the Fig. \ref{NonlinearTimeFigure}, especially  Runs I and II,   bear a strong resemblance to 
Fig. 4 of \cite{O76}.  These figures contain two positive and one negative lobe with similar boundaries separating these regions, consistent with the characteristic pattern mentioned above.   \cite{O76} does not make the hydrostatic approximation, used the GM75 model as the basis of his evaluations and is constrained to the resonance manifold.  We have made the hydrostatic approximation, base our evaluations on the GM76 model and have included resonance broadening.  
Similarities are also apparent with Fig. 12 of \citep{MB77} and Fig.s 10 and 11 of McComas and M\"uller (1981).  In those resonant evaluations using the GM76 model, the hydrostatic approximation was invoked and interactions with frequencies greater than $N/3$ were excluded.  The major difference is that the zero line separating the positive and negative lobe at high frequencies has moved to $10f$.  


\section{Discussion \label{Conclusion}}

In this paper we have review different approaches for wave-wave
interactions that have been presented in the literature in the last three
decades. Namely, we have concentrate on the approaches of
\citet{MO75,Voronovich,Zeitlin,LT,LT2}. 
%
In the absence of background rotation, we demonstrate that these four approaches produce {\em equivalent\/} results on the resonant manifold. 

This statement is not trivial given the different assumptions
and coordinate systems that have been used for the derivation of the 
various kinetic equations. It points to an internal consistency on the resonant manifold that we
still do not completely understand and appreciate. 

This result is less surprising for the canonical Hamiltonian approaches \citep{Voronovich,LT}. A canonical Hamiltonian representation is the gold standard of wave turbulence.  It guarantees that the symmetries and hence conservation principles of the original equation set in the spatial/temporal domains have been preserved in the spectral domain,\citep[e.g.][]{ZLF}.  Thus, if Voronovich's Clebsh variable representation in Eulerian coordinates and the Lvov and Tabak isopycnal Hamiltonian describe the same physical system, then there is a canonical transformation that
connects these two Hamiltonians. It is well known that such a
canonical transformation reduces to the identity transformation on the
resonant manifold. To prove this statement one constructs a
near-identical canonical transformation, which is applicable for
weakly nonlinear systems (See Appendix A3 in \cite{ZLF}).  
The Hamiltonian on the resonant manifold is invariant under a canonical near-identity transformation. 

That is why Voronovich's matrix elements
(\ref{eq:Voronovich}) look identical to the interaction matrix element
in {\em isopycnal\/} coordinates (\ref{Hamiltonian}) {\em on the
resonant manifold\/}.

We can argue that, while the other two matrix elements
(\cite{Zeitlin} and \cite{MO75}) are not in a canonical Hamiltonian formulation,
they nevertheless do describe the same physical system.  Consequently,
they also can be approximated by a
certain Hamiltonian structure, at least for small nonlinearities.  This is explicitly the case for the noncanonical Hamiltonian of \citep{MO75}.  It appears to be implicitly true of the \citet{Zeitlin} non-Hamiltonian kinetic equation.  Therefore equivalence of the scattering matrix
element on the resonant manifold is an intuitive, yet not trivial
statement. 

On the other hand, it is also intuitive that there will be coordinate representation dependent differences.  
It is a robust observational fact that Eulerian frequency spectra at high vertical wavenumber are contaminated by vertical Doppler shifting:  near-inertial frequency energy is Doppler shifted to higher frequency at approximately the same vertical wavelength.  Use of an isopycnal coordinate system considerably reduces this artifact \citep{SandP91}.  Further differences are anticipated in a fully Lagrangian coordinate system \citep{Pinkel08}.  Thus differences in the approaches may represent physical effects and what is a stationary state in one coordinate system may not be a stationary state in another.  In particular, differences may represent the effects of resonance broadening associated with Doppler shifting.  

We also demonstrate that the isopycnal and Lagrangian coordinate system approaches predict qualitatively different results with the inclusion of the near-resonant
interactions and background rotation.  The canonical Hamiltonian isopycnal formalism
is insensitive to off-resonant interactions:  Broadening the resonance width by an order of magnitude does not create significant differences in the nonlinearity parameter.  The noncanonical Lagrangian coordinate representation is, in contrast, quite sensitive to these changes.

As explained above, the Hamiltonian on the resonant manifold is
invariant under near-identity canonical transformations. The kinetic
equation describes the spectral transfers associated with the
cubic terms of the Hamiltonian and conserves the energy associated
with quadratic terms of the Hamiltonian. The kinetic equation should therefore return
representation independent results on the resonant manifold. This
statement is no longer true for near-resonant interactions.  

Indeed, since the structure of the Hamiltonian may be altered off the resonant
manifold by a near-identity canonical transformations, one {\em should} anticipate 
representation dependent
differences in spectral energy transfer when near-resonant interactions are included.  
Such differences become more
and more significant as nonlinearity increases and cubic parts of the
Hamiltonian become increasingly large.


We would like to suggest that the differences between \citet{MO75} and \citet{LT2} off the resonant manifold represent physical effects.  However, an issue with extant Lagrangian coordinate representations is that they require a small amplitude assumption that represents an unconstrained approximation whose domain of validity {\em vis-a-vis} the weak interaction approximation is not well defined, \citep{M86}.  On the basis of estimates of how horizontal Doppler shifting contributes to isopycnal spectra, we would anticipate that the Lagrangian coordinate stationary state would have typically steeper spectral slopes in the frequency domain than frequency spectra in isopycnal coordinates.  The results presented here indicate that resonance broadening will quickly whiten the high frequency Lagrangian coordinate spectrum, in direct contradiction to our intuition regarding physical effects.  


In this paper we have shown that while on the resonant manifold
(i.e. for weakly nonlinear interactions) all approaches we considered
do agree, inclusion of the near-resonant interactions (for stronger
nonlinearities) should be done with care. Results with near resonant
interactions are representation dependent.  This observations warrants
further study.

\vspace*{\baselineskip}

We thank V.~E. Zakharov for presenting us with a book \citep{Mir} and
for encouragement.  We also thank E.~N. Pelinovsky for providing us
with \citet{PR77}.  This research is supported by NSF CMG grants
0417724 and 0417466.  Y.~L. was also supported by NSF CAREER DMS
0134955.  We are grateful to YITP in Kyoto University for permitting
use of their facility.


\newpage

\begin{figure}[tp]
 \begin{center} 
  \includegraphics{fig/MO.0_psfrag.eps}%
  \includegraphics{fig/V.0_psfrag.eps}
  \includegraphics{fig/Z.0_psfrag.eps}%
  \includegraphics{fig/H.0_psfrag.eps}
\caption{Matrix elements   $|{V^{\bm{p}}_{\bm{p}_1, \bm{p}_2}}_{\mathrm{(\ref{eq:sol2})}}|^2$ given by the solution (\ref{eq:sol2}).
upper left: $|{V^{\bm{p}}_{\bm{p}_1, \bm{p}_2}}_{\mathrm{(\ref{eq:sol2})}}^{\mathrm{MO}}|^2$ according to \citet{MO75},
upper right: $|{V^{\bm{p}}_{\bm{p}_1, \bm{p}_2}}_{\mathrm{(\ref{eq:sol2})}}^{\mathrm{V}}|^2$ according to \citet{Voronovich},
bottom left: $|{V^{\bm{p}}_{\bm{p}_1, \bm{p}_2}}_{\mathrm{(\ref{eq:sol2})}}^{\mathrm{CZ}}|^2$ according to \citet{Zeitlin},
bottom right: $|{V^{\bm{p}}_{\bm{p}_1, \bm{p}_2}}_{\mathrm{(\ref{eq:sol2})}}^{\mathrm{H}}|^2$ according to \citet{LT}.
}
\label{FIGRESONANTa}
 \end{center}
\end{figure}

\begin{figure}[tp]
\begin{center}
  \includegraphics{fig/MO.1_psfrag.eps}%
  \includegraphics{fig/V.1_psfrag.eps}
  \includegraphics{fig/Z.1_psfrag.eps}%
  \includegraphics{fig/H.1_psfrag.eps}
\caption{Matrix elements $|{V^{\bm{p}_1}_{\bm{p}_2, \bm{p}}}_{\mathrm{(\ref{eq:sol3})}}|^2$ given by the solution (\ref{eq:sol3}).
upper left: $|{V^{\bm{p}_1}_{\bm{p}_2, \bm{p}}}_{\mathrm{(\ref{eq:sol3})}}^{\mathrm{MO}}|^2$ according to \citet{MO75},
upper right: $|{V^{\bm{p}_1}_{\bm{p}_2, \bm{p}}}_{\mathrm{(\ref{eq:sol3})}}^{\mathrm{V}}|^2$ according to \citet{Voronovich},
bottom left: $|{V^{\bm{p}_1}_{\bm{p}_2, \bm{p}}}_{\mathrm{(\ref{eq:sol3})}}^{\mathrm{CZ}}|^2$ according to \citet{Zeitlin},
bottom right: $|{V^{\bm{p}_1}_{\bm{p}_2, \bm{p}}}_{\mathrm{(\ref{eq:sol3})}}^{\mathrm{H}}|^2$ according to \citet{LT}.
}
\label{FIGRESONANTb}
 \end{center}
\end{figure}

\begin{figure}[tp]
\begin{center}
  \includegraphics{fig/MO.6_psfrag.eps}%
  \includegraphics{fig/V.6_psfrag.eps}
  \includegraphics{fig/Z.6_psfrag.eps}%
  \includegraphics{fig/H.6_psfrag.eps}
\caption{Matrix elements 
$|{V^{\bm{p}_1}_{\bm{p}_2, \bm{p}}}_{\mathrm{(\ref{eq:sol4})}}|^2$ given by the solution (\ref{eq:sol4}).
upper left: $|{V^{\bm{p}_1}_{\bm{p}_2, \bm{p}}}_{\mathrm{(\ref{eq:sol4})}}^{\mathrm{MO}}|^2$ according to \citet{MO75},
upper right: $|{V^{\bm{p}_1}_{\bm{p}_2, \bm{p}}}_{\mathrm{(\ref{eq:sol4})}}^{\mathrm{V}}|^2$ according to \citet{Voronovich},
bottom left: $|{V^{\bm{p}_1}_{\bm{p}_2, \bm{p}}}_{\mathrm{(\ref{eq:sol4})}}^{\mathrm{CZ}}|^2$ according to \citet{Zeitlin},
bottom right: $|{V^{\bm{p}_1}_{\bm{p}_2, \bm{p}}}_{\mathrm{(\ref{eq:sol4})}}^{\mathrm{H}}|^2$ according to \citet{LT}.
}
\label{FIGRESONANTc}
 \end{center}
\end{figure}

\begin{figure}[tp]
 \begin{center}
  \includegraphics{fig/ES.0_psfrag.eps}
  \includegraphics{fig/ID.0_psfrag.eps}
  \includegraphics{fig/PSI.1_psfrag.eps}
  \caption{upper: Matrix elements $|{V^{\bm{p}}_{\bm{p}_1, \bm{p}_2}}_{\mathrm{ES}}|^2$ given by the solution (\ref{eq:sol1}).
middle: Matrix elements $|{V^{\bm{p}}_{\bm{p}_1, \bm{p}_2}}_{\mathrm{ID}}|^2$ given by the solution (\ref{eq:sol2}).
  bottom: Matrix elements $|{V^{\bm{p}_1}_{\bm{p}_2, \bm{p}}}_{\mathrm{PSI}}|^2$ given by the solution (\ref{eq:sol3}),
  which gives PSI as $|\bm{k}_1| \to 0$ ($\epsilon \to 0$).
The matrix elements here are shown as functions of $\epsilon$ such that $(|\bm{k}_1|, |\bm{k}_2|) = (\epsilon, \epsilon/3+1) |\bm{k}|$. All four versions of the Matrix elements are plotted here:  the appearance of a single line in each figure panel testifies to the similarity of the elements on the resonant manifold.  
}
\label{FigureThree}
 \end{center}
\end{figure}




\begin{figure}[tp] 
\begin{center}
\includegraphics[scale=0.5]{fig/NewEpsilon.eps}
\caption{Nonlinearity parameter (\ref{NonlinearRatio}) for the Garrett
and Munk spectrum (\ref{GM}) calculated via
(\ref{KineticEquationBroadened}).  
%
The upper figures  represent the
value of nonlinearity parameter calculated using \cite{LT2}, equation
(\ref{LTV}) with $\delta=10^{-3}$, {Run I} (upper left)  and $\delta=10^{-2}$
{Run III} (upper right).  The bottom two pictures represent the
value of nonlinearity parameter calculated via \cite{MO75}, (\ref{VMO})
with $\delta=10^{-3}$ {Run II} (bottom left) and $\delta=10^{-2}$
{Run IV} (bottom  right).}
\label{NonlinearityParameter}
\end{center}
\end{figure}





\begin{figure}[tp]
\begin{center}
\includegraphics[scale=0.5]{fig/NewTau.eps}
\caption{Nonlinear time (\ref{NonlinearTime}) for the Garrett
and Munk spectrum (\ref{GM}) calculated via
(\ref{KineticEquationBroadened}).  
%\todo{Naoto, what is white on bottom-right figure?}
The upper figures  represent the
value of nonlinearity parameter calculated using \cite{LT2}, equation
(\ref{LTV}) with $\delta=10^{-3}$, {Run I} (upper left)  and $\delta=10^{-2}$
{Run III} (upper right).  The bottom two pictures represent the
value of nonlinearity parameter calculated via \cite{MO75}, (\ref{VMO})
with $\delta=10^{-3}$ {Run II} (bottom left) and $\delta=10^{-2}$
{Run IV} (bottom  right). On this bottom right figure white region 
to the left of the $0.1$ contour corresponds to extremely fast time scales, 
faster then $0.1$ of a day. On these figures, 
$\omega$ in cpd, $m$ in cycle/m, and nonlinear time $\tau^{\mathrm{NL}}$ 
is measured in days.}
\label{NonlinearTimeFigure}
\end{center}
\end{figure}




\begin{figure}[tp]
\begin{center}
\includegraphics[scale=.5]{fig/NewDifference.eps}
\caption{Differences between nonlinearity parameter
(\ref{NonlinearRatio}) calculated via \cite{MO75} and \cite{LT2} with
$\delta=10^{-3}$, i.e. between {Run I} and {Run II}, (a),
calculated with \cite{LT2} with $\delta=10^{-2}$ and $\delta=10^{-3}$,
i.e. the difference {Run I} and {Run III} (b), and finally
between \cite{MO75} with $\delta=10^{-2}$ and $\delta=10^{-3}$, i.e.
between {Run II} and {Run IV} (c).}
\label{Differences}
\end{center}
\end{figure}

\clearpage

\bibliographystyle{ametsoc}

\begin{thebibliography}{100}

\bibitem[{Brehovsky(1975)}]{B75}
Brehovski. 1975: On interactions of internal and surface waves in the ocean. {\em Oceanology\/},
{\bf 15} (in Russian). This is citation [11] of \citet{PR77}. 

\bibitem[B\"uhler and McIntyre, 2005]{BM05} B\"uhler, O. and M. E. McIntyre, 2005:  Wave capture and wave-vortex duality. {\em J. Fluid Mech.}, {\bf 534}, 67--95.  
         
\bibitem[{Caillol and Zeitlin(2000)}]{Zeitlin} Caillol, P. and Zeitlin, V., 2000: Kinetic equations and stationary energy spectra of weakly nonlinear internal gravity waves. {\em Dynamics of Atmospheres and Oceans\/}, {\bf 32}, 81-112.

\bibitem[{Dyachenko et~al.(1995)Dyachenko, Lvov, and
  Zakharov}]{dyachenko1995fwi}
Dyachenko, A.~I., Y.~V. Lvov, and V.~Zakharov, 1995: Five-wave interaction on
  the surface of deep fluid. \textit{Physica D}, \textbf{87}, 233--261.

\bibitem[{Gershgorin {\it et al.}(2007)}]{Gersh2007} Gershgorin B., Lvov, Y. V., and Cai, D., 2007:
Interactions of renormalized waves in thermalized Fermi-Pasta-Ulam chains. {\em Phys. Rev. E\/} {\bf 75}, 046603.

\bibitem[{Hasselmann(1966)}]{H66} Hasselmann, K., 1966: Feynmann diagrams and interaction rules of wave-wave scattering processes.  {\em Rev. Geophys.\/}, {\bf 4}, 1--32.

\bibitem[{Janssen(2003)}] {JansenXXX} Janssen P. A. E. M., 2003: Nonlinear four-wave interactions and freak waves. {\em J. Phys. Oceanogr.\/}, {\bf 33}, 863--884.

\bibitem[{Kenyon(1966)}]{K66} Kenyon, K. E., 1966: {\em Wave-Wave Scattering for Gravity Waves and Rossby Waves}.  Ph. D. dissertation, UCSD, San Diego; 93 pp.

\bibitem[{Kenyon(1968)}]{K68} Kenyon, K. E., 1968: Wave-wave interactions of surface and internal waves, {\em J. Mar. Res.\/}, {\bf 26}, 208--231.

\bibitem[{Krasitskii(1994)}]{kras}
Krasitskii, V.~P., 1994: On reduced equation in the {H}amiltonian theory of
  weakly nonlinear surface waves. \textit{J. Fluid Mech.}, \textbf{272}, 1--20.

\bibitem[{Lelong and Riley(1991)}]{Lelong} Lelong, M. P. and J. J. Riley, 1991:  Internal wave - vortical mode interactions in strongly stratified flows.  {\em J. Fluid Mech.\/}, {\bf 232}, 1--19.


\bibitem[{Lvov {\it et al.}(1997)}]{LLNZ}
      Lvov, V.~S., Lvov, Y.~V., Newell, A.~C. and Zakharov,V.~E., 1997: Statistical description of acoustic turbulence,
      {\em Phys. Rev. E\/}, {\bf 56}, 390--405.

\bibitem[{Lvov and Tabak(2001)}]{LT} Lvov, Y.~V., and E. G. Tabak, 2001: Hamiltonian formalism and the Garrett and Munk spectrum of internal waves in the ocean.  {\em Phys. Rev. Lett.\/}, {\bf 87}, 169501.

\bibitem[{Lvov and Tabak(2004)}] {LT2} Lvov, Y.~V., and Tabak E.~G., 2004: A Hamiltonian Formulation for Long Internal Waves. {\em Physica D} {\bf 195} 106--122. Also at http://www.rpi.edu/$\tilde{\ }$lvovy.

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


\bibitem[{Lvov and Nazarenko (2004)}]{NoisyNazarenko} 
Y.V. Lvov and S. Nazarenko, `` Noisy spectra, long
correlations, and intermittency in wave turbulence,'' {\it Physical
Review E}, {\bf 69}, 066608.

\bibitem[{Lvov {\it et al.}(2006)}] {LvovNazarenkoPokorni} 
Y.Lvov, S. Nazarenko, B.Pokorni, "Discreteness and its effect on the water-wave turbulence."  Physica D: {\bf 218}, 24-35(2006).




\bibitem[{McComas(1975)}]{McC75} McComas, C.~H., 1975:  {\em Nonlinear Interactions of Internal Gravity Waves}, Ph. D. Thesis, The Johns Hopkins University, Baltimore, 74 pp.

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

\bibitem[{McComas and Bretherton(1977)}] {MB77} McComas, C.~ H. and Bretherton, F.~P.,1977: Resonant interaction of oceanic internal waves.  {\em J. Geophys. Res.\/}, {\bf 82}, 1397--1412.

\bibitem[{Milder(1982)}]{Milder} Milder M., 1982: Hamiltonian description of internal waves. {\em J. Fluid Mech.\/}, {\bf 119} 269--282.

\bibitem[{Miropolsky(1981)}]{Mir} Miropolsky, Y.~Z. 1981: Dinamika vnutrennih gravitacionnih voln v okeane,
(Dynamics of internal gravity waves in the ocean), Leningrad, Gidrometeroizdat, (in Russian). 

\bibitem[{M\"uller and Olbers(1975)}]{MO75} M\"uller, P., and Olbers, D.~J., 1975: On the dynamics of internal waves in the deep ocean,
{\em J. Geophys. Res.\/}, {\bf 80}, 3848--3860.

\bibitem[{M\"uller {\it et al.}(1986)}]{M86} M\"uller, P., Holloway, G., Henyey, F., Pomphrey, N. 1986:   Nonlinear interactions among internal gravity waves.  {\em Rev.\ Geophys.\/}, {\bf 24}, 493--536.

\bibitem[{Olbers(1974)}]{O74} Olbers D.~J., 1974: On the energy balance of small scale internal waves in the deep sea, {\em Hamburg, Geophys. Einzelschriftern\/} {\bf 27}.

\bibitem[{Olbers(1976)}]{O76} Olbers D.~J., 1976:  Nonlinear energy transfer and the energy balance of the internal wave field in the deep ocean.  {\em J. Fluid Mech.\/}, {\bf 74}, 375--399.

\bibitem[Pinkel, 2008]{Pinkel08} R. Pinkel, 2008:  Advection, phase distortion, and the frequency spectrum of finescale fields in the sea.  {\em J. Phys. Oceanogr.}, {\bf 38}, 291--313.    

\bibitem[{Pelinovsky and Raevsky(1977)}]{PR77} Pelinovsky, E. N., and Raevsky, M. A., 1977.  Weak turbulence of the internal waves of the ocean.  {\em Atm. Ocean Phys.-Izvestija\/}, {\bf 13}, 187--193.

\bibitem[{Pomphrey {\it et al.}(1980)}] {PMW80} Pomphrey, N., Meiss, J. D. and Watson, K. D., 1980:  Description of nonlinear internal wave interactions using Langevin methods.  {\em J. Geophys. Res.\/}, {\bf 85}, 1085--1094.

\bibitem[Polzin (2004)]{Polzin2004a}
{\sc Polzin, K. L.} 2004 A heuristic description of internal wave  dynamics.  
  {\em J. Phys. Oceanogr.\/} {\bf 34}, 214--230.

\bibitem[Polzin, 2008]{P08} Polzin, K. L., 2008:  Mesoscle Eddy - Internal Wave Coupling.    I.  Symmetry, wave capture and results from the Mid-Ocean Dynamics Experiment.  {\em J. Phys. Oceanogr.},  in preparation.  

\bibitem[{Sanderson(1985)}]{S85} Sanderson, B., 1985:  A Lagrangian solution for internal waves.  {\em J. Fluid Mech.\/}, {\bf 152}, 191--202.

\bibitem[{Sherman and Pinkel(1991)}]{SandP91} Sherman, J. T., and Pinkel, R., 1991:  Estimates of the vertical wavenumber-frequency spectra of vertical shear and strain., {\em J. Phys. Ocean.\/}, {\bf 21}, 292--303.

\bibitem[{Voronovich(1979)}]{Voronovich} Voronovich, A.~G., 1979: Hamiltonian formalism for internal waves in the ocean. {\em Izvestiya, Atmospheric and Oceanic Physics\/}, {\bf 16} 52--57.

\bibitem[{Voronovich et~al.(2006)Voronovich, Sazonov, and
  Shrira}]{2006JFM...568..273V}
Voronovich, V.~V., I.~A. Sazonov, and V.~I. Shrira, 2006: On radiating solitons
  in a model of the internal wave shear flow resonance. \textit{J. Fluid
  Mech.}, \textbf{568}, 273--301.

\bibitem[{Annenkov and Shrira, V.I.(2006)}]{Shrira} Annenkov,
S.Yu. and Shrira, V.I, 2006: Role of non-resonant interactions in the
evolution of nonlinear random water wave fields.\textit{J. Fluid
  Mech.}, \textbf{561}, 181-207. 

\bibitem[{Yuen and Lake(1982)}]{yuen_lake}
Yuen, H.~C. and B.~M. Lake, 1982: Nonlinear dynamics of deep-water gravity
  waves. \textit{Adv. Appl. Mech.}, \textbf{22}, 67--229.

\bibitem[{Zakharov {\it et al.}(1992)}]{ZLF} Zakharov, V.~E., Lvov, V.~S. and Falkovich G. 1992: {\em Kolmogorov Spectra of Turbulence\/}.  Springer-Verlag.

\end{thebibliography}



\end{document}

