---
abstract: 'Here we show that there exist internal gravity waves that are inherently unstable, that is, they cannot exist in nature for a long time. The instability mechanism is a one-way (irreversible) harmonic-generation resonance that permanently transfers the energy of an internal wave to its higher harmonics. We show that, in fact, there are countably infinite number of such unstable waves. For the harmonic-generation resonance to take place, nonlinear terms in the free surface boundary condition play a pivotal role, and the instability does not obtain for a linearly-stratified fluid if a simplified boundary condition such as rigid lid or linear form is employed. Harmonic-generation resonance presented here also provides a mechanism for the transfer of the energy of the internal waves to the higher-frequency part of the spectrum where internal waves are more prone to breaking, hence losing energy to turbulence and heat and contributing to oceanic mixing.'
author:
- 'Y. Liang, Ahmad Zareei, M.-Reza Alam[^1]'
title: Inherently Unstable Internal Gravity Waves due to Resonant Harmonic Generation
---

Introduction
============

Internal gravity waves, outcome of perpetually agitated density-stratified oceans, are known to play a critical role in the dynamics of our planet’s energy balance: they absorb energy to form, carry energy over long distances as they propagate, and release energy where they break [@Staquet2002]. The latter phenomenon usually gives rise to considerable mixing [c.f. @Ferrari2008] whereby nutrients also get distributed, which is vital for a wide range of marine life [@Boyd2007; @Harris2012].

More than a century long research has shed a lot of light on various features of internal gravity waves. Nevertheless, many aspects of their inception and fate is yet not well understood [e.g. @Alford2015]. Specifically, the precise mechanism that transfers energy from longer waves to the high-frequency part of the spectrum, where internal waves are more prone to breaking, is yet a matter of dispute. Aside from linear processes such as interaction of internal waves with the seabed topography and sloped continental shelves , several nonlinear instability mechanisms have also been put forward. For instance, we now know that internal waves may undergo instability due to triad resonance . All discovered destabilizing mechanisms for an internal wave (few named above), however, have one thing in common that they require some type of perturbations in order to get initiated. These perturbations can come from, for instance, seabed corrugations or presence of other waves forming resonance triads.

Here, we show that there are internal gravity waves in the ocean that are *inherently* unstable, that is, they simply cannot sustain their form. Through the mechanism studied here, specific internal waves *naturally* (without requiring any perturbation) give up their energy *permanently* to their higher harmonics through a one-way irreversible harmonic-generation resonance mechanism.

Governing Equations and the Dispersion Relation
===============================================

Consider the propagation of internal waves in an inviscid, incompressible, adiabatic and stably stratified fluid of density $\rho(x,y,z,t)$, bounded by a free surface on the top and a rigid seafloor at the depth $h$. Let’s define a Cartesian coordinate system with $x,y$-axes on the mean free surface and $z$-axis positive upward. Newton’s second law, conservation of mass, and conservation of energy provide five equations for the evolution of the components of the velocity vector $\b{u} = \{u, v, w\}$, density $\rho$, and the pressure $p$. These governing equations together with three boundary conditions (two kinematic boundary conditions on the free surface and the seabed, and one dynamic boundary condition on the free surface) uniquely determine the five unknowns and the surface elevation $\eta(x,y,t)$ [e.g. @Thorpe1966].

We assume internal waves are small perturbations to a stable background state at equilibrium. Therefore, density can be written as $\rho(x,y,z,t) = \bar{\rho}(z)+ \rho' (x, y, z,t)$ where $\bar{\rho}(z)$ is the background (unperturbed) density. Similarly, we define a pressure perturbation $p'$ via $p=\bar p(z)+p'(x,y,z,t)$ such that $d\bar p(z)/dz=-\bar{\rho}(z) g$. With some standard manipulation, the governing equations can be written in terms of either of the five variables involved in this problem. We choose to write the equation, as is customary, in terms of the vertical component of the velocity, $w$. These equations then read [see e.g. @Thorpe1966 or Appendix] &\^2w+N\^2\^2\_Hw=(,’),&-h&lt;z&lt;,\[101\]\
&-g\^2\_H w = (,p’,), & z=0,\[102\]\
& w=0, &z=-h.\[103\] where $\nabla^2_H=\p^2/\p x^2+\p^2/\p y^2$ is the horizontal Laplacian, $N^2=-{g/\rho_0~ \d \bar{\rho}(z)}/{ \d z}$ is the Brunt-V[ä]{}is[ä]{}l[ä]{} frequency in which $\rho_0=\bar{\rho}(z=0)$ is the density on the free surface, and $\E,\F$ are nonlinear functions of their arguments.

To perform a perturbation analysis, we assume that the solution to can be expressed in terms of a convergent series, i.e. $$\begin{aligned}
\label{110}
w(\x,t)=\ep w^{(1)}(\x,t)+\ep^2 w^{(2)}(\x,t)+\O(\ep^3),\end{aligned}$$ where $\ep\ll1$ is a measure of steepness of the waves involved and $w^{(i)}\sim \O(1)$. Similar expressions hold for $u,v,\rho'$ and $p'$. Substituting into and collecting terms of the same magnitude, then at the leading order $\O(\ep)$ the linearized equations are obtained.

We focus our attention here on the two-dimensional problem with a linear mean density profile, i.e. $\bar{\rho}(z)=\rho_0 (1-a z)$ which gives a constant Brunt-V[ä]{}is[ä]{}l[ä]{} frequency $N=\sqrt{ga}$ [c.f. e.g. @Martin1972]. Looking for a progressive wave solution of the leading order (linearized) equation in the form $w^{(1)}=W(z)\sin({\bf k}\cdot {\bf x} -\omega t)$ the following dispersion relations result:

$$\begin{aligned}
\label{2000}
  \hspace{-.70cm} \D(k,\omega)=
\begin{cases}
    \omega^2-\f{g k}{\sqrt{1-{N^2}/{\omega^2}}}\tanh\lp kh \sqrt{1-{N^2}/{\omega^2}}\rp=0,&  \omega>N\\
    \omega^2-\f{g k}{\sqrt{N^2/\omega^2-1}}\tan\lp kh \sqrt{N^2/\omega^2-1} \rp=0,              & \omega<N
\end{cases}\end{aligned}$$

![Plot of the dimensionless frequency $\omega/N$ as a function of dimensionless wavenumber $kh$ of free internal waves (i.e. $\D(k,\omega)=0$) in a fluid of linearly stratified density $\rho(z) =\rho_0(1 - az)$, with $ah$ = 0.05. Associated with each wavenumber there is one surface wave and an infinite number of internal wave modes (blue solid-line branches). Frequency of internal waves cannot exceed the Brunt-V[ä]{}is[ä]{}l[ä]{} frequency $N$, and all branches of the dispersion relation curve are capped at $\omega/N$=1. We also plot contours of $\D(2k,2\omega)$=0 (red dash-dotted lines) whose intersections with blue lines (shown by red circles) mark waves whose second harmonics are also solutions to the dispersion relation. These second harmonics are at the intersections of contours of $\D(k/2,\omega/2)$=0 (green dashed lines) and $\D(k,\omega)$=0 and are marked by black squares. The second harmonic of the wave at “***a***" (mode 2), is the wave “***A***" (mode 1) and so on. Note that second harmonic waves are at least one mode lower than the original waves.[]{data-label="fig1"}](disp94.eps){width="8cm"}

Solutions to the above dispersion relation identify permissible frequency and wavenumber of free propagating waves. Contours of $\D(k,\omega)=0$ are shown in figure \[fig1\] in which we plot the dimensionless frequency $\omega/N$ as a function of dimensionless wavenumber $kh$ (blue solid curves). For $\omega>N$ only one solution exists in the first quadrant (with its mirrors in the other quadrants). This solution corresponds to a wave whose associated fluid particle motion is maximum near the free surface and decreases as the depth increases. Therefore this is basically a classical *surface* wave which is a little perturbed because of stratification. For $\omega<N$ there is an infinite number of solutions to . The first member of this set, is the continuation of the surface wave branch (the left-most branch in figure \[fig1\]), but the rest identify waves with associated fluid particle activities that are minimum near the free surface and the seabed, but gain one (or more) maximum/maxima somewhere inside the fluid domain. Therefore these branches show *internal* waves. The number of maxima in the amplitude of velocity along the vertical water line determines the mode number of the branch (the first three are marked on the right-side of figure \[fig1\] with arrows). A cut-off frequency $\omega/N$=1 sets an upper frequency limit for internal waves. It is to be noted that the dispersion relation , although has a different form, is in fact graphically very close to the one under rigid lid assumption $\omega=N/\sqrt{1+(n\pi/kh)^2}$. But clearly the inclusion of the effect of the free surface in the former has changed its form.

Harmonic Generation
===================

With the linear solution to and its properties at hand, we move to the second order equation by collecting $\O(\ep^2)$ terms obtained from the substitution of into . The second order equation has the exact same form as of the leading order equation on its left-hand side, but with nonlinear terms, arising from nonlinear functions $\E,\F$, on its right-hand side. These nonlinear terms are multiplication of the linear solution (and its derivatives) and, it turns out that, they constitute forcing terms with wavenumber and frequency $(2k,2\om)$. Now if $\D(2k,2\om)$=0, then it means that these forcing terms have a harmonic which is the same as the natural harmonic of the linear system. This is a resonance scenario through which a new second order solution may emerge and may grow large enough to the extent that it becomes comparable to the leading order solution (beyond which the naive expansion is not valid anymore). We would like to note a subtle point here that at the second order the right-side of is identically zero [see e.g. @Tabaei2005]. Therefore a potential harmonic-generation resonance is only possible through the nonlinear terms in the free surface boundary condition . If a rigid-lid assumption or a linearized form of the free surface boundary condition is employed (which is usually the case in the investigation of internal waves) the resonance harmonic generation will simply not obtain unless a non-uniform stratification or Boussinesq terms are considered.

To see whether it is possible to satisfy the resonance condition required for the second harmonic of an internal wave to exist as a free propagating wave, we also plot in figure \[fig1\] the contours of $\D(2k,2\om)$=0 (red dash-dotted curves) for first four internal wave modes. Intersections of these contours with the contours of $\D(k,\om)$=0 (solid blue curves) identify waves whose second harmonics are also solution to the dispersion relation. Some of these intersections are marked in figure \[fig1\] by red circles and are identified by lowercase characters. It is easy to find the corresponding second harmonic by multiplying a designated (red circle) frequency and wavenumber by a factor of two, or alternatively by finding intersections of contours of $\D(k/2,\om/2)$=0 (green dashed curves) and $\D(k,\omega)$ (denoted by black squares and uppercase characters). Specifically, “***A***" is the second harmonic of $``\textit{\textbf{a}}"$, $``\textit{\textbf{B}}"$ is the second harmonic of $``\textit{\textbf{b}}"$ and so on. It is to be noted that the second harmonic of an internal wave always belongs to a lower mode than the mode of the original wave. For example, second harmonic of wave “***a***" (mode 2) is the wave “***A***" (mode 1), and second harmonic of wave “***c***" (mode 3) is the wave “***C***" (mode 1).

A comprehensive collection of all internal waves (for $kh<$12) whose second harmonics are also free propagating waves is shown in figure \[fig2\]. We are basically plotting in this figure all intersections points of contours of $\D(k,\om)=0$ and contours of $\D(2k,2\om)=0$. Clearly there will be an infinite (but countable) number of such solutions. Through the approach described above, a similar exercise can be performed for waves with their third harmonic being free waves. These solutions are marked by blue circles in figure \[fig2\], and waves with their fourth harmonics lie on the dispersion relation curves are shown by green triangles, and this search can continue indefinitely.

![There are countably infinite number of internal waves that are unstable to their second harmonic. Physically this means that specific incident internal waves of wavenumber and frequency ($k,\omega$) will give up their energy permanently (in a one-way irreversible process, c.f. equation ) to their second harmonic ($2k,2\om$). The necessary condition for this to happen is $\D(k,\om)=\D(2k,2\om)$=0. These waves, for the parameters of figure \[fig1\], are shown here by red squares. Similar story holds for another set of waves that are unstable to their third harmonic (blue circles, necessary codition $\D(k,\om)=\D(3k,3\om)$=0), and fourth harmonic (green triangles, $\D(k,\om)=\D(4k,4\om)$=0) and so on. Instability to higher harmonics are clearly much weaker when compared with the instability to the second harmonic.[]{data-label="fig2"}](2nd3rd4th_v3.eps){width="8cm"}

Results and Discussions
=======================

To determine the strength of the resonance (i.e. the rate of growth of the resonant wave), and the dynamics of the energy interplay between a wave and its second harmonic, here we perform a multiple scale perturbation analysis. The basic assumption is that the amplitude of the original wave and its second harmonic are both functions of spatial variables and time (i.e. $\x,t$), *as well as* a slow spatial variable in the direction of propagation $x_1=\ep x$. Physically speaking, we allow the amplitude of both waves to *slowly* vary as waves propagate. Mathematically speaking this is written as $$\begin{aligned}
\label{210}
w(x,x_1,z,t)=\ep w_1(x,x_1,z,t)+\ep^2 w_2(x,x_1,z,t)+\O(\ep^3)\end{aligned}$$ where $w_i\sim\O(1)$. Expressions of the same form are assumed to hold for other variables $u,\rho',p'$ and $\eta$. Similar to regular perturbation methodology, described earlier in the paper to gain insight into the problem, by substituting into the governing equation and collecting terms of the same order in $\ep$. We assume the original wave with wavenumber and frequency $k,\omega$ has the amplitude $\A_1(x_1)$ and the amplitude of resonant second harmonic wave with wavenumber and frequency $2k,2\omega$ is $\B_2(x_1)$. At the second order, applying a compatibility condition (to avoid unbounded solutions) the following two equations governing spatial evolution of $\A_1(x_1)$ and $\B_2(x_1)$ emerge (see Appendix for details of derivation) &=\_1\^2(x\_1),\[201\]\
&=\_1(x\_1)\_2(x\_1),\[202\] where $$\begin{aligned}
&\alpha=-\f{6m_1^2\om\cos m_2h}{ g k(2m_2h+\sin 2 m_2 h) },\nn\\
&\beta=-\f{k[\sin m_2 h(4m_1\cos^2 m_1h -3m_1) -2m_2\cos m_2h\sin 2m_1h]}{2\omega(2m_1h+\sin 2m_1h)}. \nn\end{aligned}$$ in which $m_1=k\sqrt{N^2-\omega^2}/\omega$ and $m_2=k\sqrt{N^2-4\omega^2}/\omega$.

Spatial evolution of the normalized amplitude of the original wave $\A_1/\A_{10}$ (where $\A_{10}=\A_1(x=0)$) and its second harmonic $\B_2/\A_{10}$ as a function of spatial distance of propagation $x/\lambda_i$ (where $\lambda_i=2\pi/k$= wavelength of the original wave) is shown in figure \[fig3\] respectively by blue-dashed line and solid-red line. In accordance with the case presented in figure \[fig1\], we choose $ah$=0.05 and if we consider that waves are propagating in a water of depth $h$=1 km, then $N$=0.02 rad/s. For this case, figure \[fig3\]a corresponds to the point “***a***" ($kh$=3.144772,$\omega/N$=0.4472136) with its second harmonic at “***A***" in figure \[fig1\], and figure \[fig3\]b corresponds to point “***c***" ($kh$=5.132426,$\omega/N$=0.4780914) with its second harmonic at “***C***". In the former case the interaction is (spatially) faster by a factor of about two, while in the latter case the relative amplitude of the resonant wave is higher. The most striking aspect of the solution is that the interaction is one-way. This can be seen from in which, because $\A_1^2$ is always positive, then $\d \B_2/d x_1$ can never change sign. As a result the magnitude of $\B_2$ can only increase (the sign of $\alpha$ only contributes to 0 or $\pi$ radian phase shift to the wave), and that’s why the direction of energy can never change. This is in contrast to typical triad resonance interactions [e.g. @Alam2010; @Alam2012c] and harmonic generation in shallow water waves [e.g. @Alam2007] where energy initially flows from original waves to resonant waves, but then when the amplitude of resonant waves is large enough the flow of energy reverses. Here, energy only goes from the original wave to the second harmonic and stays there permanently. The original wave will be gone forever.

Our multiple scales results conserve energy, as expected. In figure \[fig3\]c (which corresponds to the case presented in figure \[fig3\]b) we plot the energy flux ($E_f=E\times C_g$ where $E$ is the energy per unit area and $C_g$ is the group velocity) normalized by the energy flux of the original wave at the beginning ($E_{f,\A_{10}}$). Plotted are the normalized energy flux of the original wave $E_{f,\A_1}$ (blue-dashed line), the second harmonic $E_{f,\B_2}$, and the summation of the fluxes $E_{f,total}=E_{f,\A_1}+E_{f,\B_2}$. As expected from the conservation of energy the latter is constant and equal to unity.

![Spatial evolution of the amplitude of the original wave (blue dashed line) and its resonant second harmonic (red solid line) correspond to the point “***a***" in figure \[fig1\] ($kh$=3.144772, $\omega/N$=0.4472136, $aA_{10}/N$=0.0023, fig. a), and the point “***c***" in figure \[fig1\] ($kh$=5.132426, $\omega/N$=0.4780914, $aA_{10}/N$=0.0023, fig. b). In figure (a) energy goes from mode 2 to mode 1, whereas in figure (b) energy goes from mode 3 to mode 1. Figure (c) shows the energy flux of each wave as well as the sum of the energy fluxes (green dash-dotted line). As expected from energy conservation, the overall energy flux is unchanged in the domain of interaction. []{data-label="fig3"}](amp_hg05_h1000_pointa_v7.eps "fig:"){width="7cm"} (-205,82)[(a)]{}\
![Spatial evolution of the amplitude of the original wave (blue dashed line) and its resonant second harmonic (red solid line) correspond to the point “***a***" in figure \[fig1\] ($kh$=3.144772, $\omega/N$=0.4472136, $aA_{10}/N$=0.0023, fig. a), and the point “***c***" in figure \[fig1\] ($kh$=5.132426, $\omega/N$=0.4780914, $aA_{10}/N$=0.0023, fig. b). In figure (a) energy goes from mode 2 to mode 1, whereas in figure (b) energy goes from mode 3 to mode 1. Figure (c) shows the energy flux of each wave as well as the sum of the energy fluxes (green dash-dotted line). As expected from energy conservation, the overall energy flux is unchanged in the domain of interaction. []{data-label="fig3"}](amp_hg05_h1000_pointc_v7.eps "fig:"){width="7cm"} (-205,82)[(b)]{}\
![Spatial evolution of the amplitude of the original wave (blue dashed line) and its resonant second harmonic (red solid line) correspond to the point “***a***" in figure \[fig1\] ($kh$=3.144772, $\omega/N$=0.4472136, $aA_{10}/N$=0.0023, fig. a), and the point “***c***" in figure \[fig1\] ($kh$=5.132426, $\omega/N$=0.4780914, $aA_{10}/N$=0.0023, fig. b). In figure (a) energy goes from mode 2 to mode 1, whereas in figure (b) energy goes from mode 3 to mode 1. Figure (c) shows the energy flux of each wave as well as the sum of the energy fluxes (green dash-dotted line). As expected from energy conservation, the overall energy flux is unchanged in the domain of interaction. []{data-label="fig3"}](energy_hg05_h1000_pointc_v7.eps "fig:"){width="7cm"} (-200,73)[(c)]{}

To cross validate our results with direct simulation, we use the adaptive Navier-Stokes solver code SUNTANS (Stanford Unstructured Nonhydrostatic Terrain-following Adaptive Navier-Stokes Simulator [@Fringer2006]). As a nonhydrostatic parallel ocean model with the capability to implement nonlinear free surface, SUNTANS has been validated and widely used in studying internal waves (and associated turbulence and mixing) in stratified waters [@Zhang2011; @Kang2010; @Wang2011; @Kang2012a; @Walter2012; @Zhang2011; @Wang2011].

Here we consider propagation of waves in a stratified water of density gradient $a=1\times10^{-3}$ m$^{-1}$, depth 100m and in a domain of horizontal extent 20km. We choose $5\times10^3$ grid points in the $x$ direction and 100 layers in the $z$ direction. We specify the velocity of fluid particle at the left boundary to match that of desired incoming wave. Therefore the left boundary that serves as the incoming wave boundary (wave maker). We set the right-side vertical boundary of the domain as the no-penetration and slip-free, and therefore it will act as a rigid vertical wall. We consider the case of incident internal wave of mode 2 ($kh$=3.147946, $\omega/N$=0.4472136, vertical velocity amplitude $w$=0.1m/s) with its second harmonic being a mode 1 wave. The time step $\delta t/T$=1/35000 where $T$ is the period of the incident wave. We run the simulation until the incident wave arrives at the wall, at which point the amplitude of both waves (incident and its second harmonics) have reached a steady state in the domain of interest $0<x/\lambda_i<20$. As is seen from figure \[fig4\], the curve of $\B_2/\A_{10}$ vs $x/\lambda_1$ gives an initial slope of 0.0107 which is the same as our theoretical prediction .

![Comparison of analytical results with the direct simulation. The physical parameters used are: density gradient $a$=${1\e-3}$ $m^{-1}$, water depth h=100m, dimensionless wave number $kh$=3.147946, dimensionless frequency $\omega/N$=0.4472136, the amplitude of vertical velocity of the incident wave $A_{10}$=0.1 m/s. []{data-label="fig4"}](suntans_v4){width="6.5cm"}

To see how fast this instability evolves temporally, we present in figure \[fig42\] results of the temporal evolution of the parent wave and its second harmonic, with the same physical parameters as in the case in figure \[fig3\]a. The governing equation takes a similar form to except the derivatives that are with respect to the time $t$ and clearly expressions for $\alpha$ and $\beta$ are different. The qualitative trend is as the spatial case, and the figure suggests that we need time $\sim\O(1000)$ times the period of the initial wave to see the majority of the energy transmitted to the second harmonic. With the chosen Brunt-V[ä]{}is[ä]{}l[ä]{} frequency of $N=0.022$s$^{-1}$, the time scale is about 10 days, which is of the same order as that of parametric subharmonic instability of the M2 internal tide ($\sim$ 2-5 days) [c.f. e.g. @Gerkema2006c; @MacKinnon2005a]. This is somewhat expected as both mechanisms are *second-order* resonance instabilities.

We would like to comment here that Parametric Subharmonic Instability is dependent upon perturbation waves in the domain in order to get started [@Davis1967; @muller1986]. In direct simulations, this is achieved by adding random noise to the simulation domain. In experimental studies the required noise already exists in the domain due to unavoidable imperfections. Hence, usually in experiments Parametric Subharmonic Instability is automatically obtained (similar to inevitable Benjamin-Feir instability) and destabilizes “single" internal waves to a number of other subharmonic waves [e.g. @Martin1972]. We would like to emphasize that the underlying mechanism of resonant harmonic generation studied here is different than that of Parametric Subharmonic Instability in that the former mechanism, among other characteristics, is an “inherent" instability that does not require ambient perturbations to get started.

Three snapshots of the overall vertical velocity $w$ are shown in figure \[fig42\]b that highlights the horizontal and vertical structures of the wave field at different times. Energy goes from the original internal wave (b1) slowly to its higher harmonics, first leading to modulation of the original wave (b2), but eventually the entire energy is at the second harmonic (b3) whose frequency and horizontal wavenumber are double the horizontal frequency and wavenumber of the initial wave. The interaction is only between the initial wave and its second harmonic, and no cascading to different waves is ensued.

![Time evolution of internal gravity waves that satisfy the conditions of harmonic generation. The physical parameters used here correspond to figure \[fig3\]a. (a) Time evolution of the amplitudes of the two waves. $T_0$ is the period of the parent wave. (b) Snapshots of the wave field at $t=0$ (b1), $t=500T_0$ (b2), and $1500T_0$(b3). The presented mechanism initially results in an modulation of the primary wave, and eventually all the energy goes to the second harmonic.[]{data-label="fig42"}](HT.eps "fig:"){width="8cm"} (-245,82)[(a)]{}\
![Time evolution of internal gravity waves that satisfy the conditions of harmonic generation. The physical parameters used here correspond to figure \[fig3\]a. (a) Time evolution of the amplitudes of the two waves. $T_0$ is the period of the parent wave. (b) Snapshots of the wave field at $t=0$ (b1), $t=500T_0$ (b2), and $1500T_0$(b3). The presented mechanism initially results in an modulation of the primary wave, and eventually all the energy goes to the second harmonic.[]{data-label="fig42"}](ss.eps "fig:"){width="8cm"} (-245,82)[(b)]{}\

We would like to note that the harmonic generation mechanism proposed here occurs for non-uniformly stratified ocean as well. As pointed out previously, a necessary condition for this to occur is that for a given wave of frequency $\omega$ and wavenumber $k$, the two conditions $\D(k,\omega)=0$ and $\D(2k,2\omega)=0$ are satisfied. We present two examples in figure \[fig5\] for two different density profiles: parabolic $\rho(z)=\rho_0 (1+a z^2)$, where $a=1\times10^{-5} \text{m}^{-2}$ (figure \[fig5\]a), and exponential $\rho(z)=\rho_0 [1+\delta-\delta \exp(az)]$, where $a=0.1 \text{m}^{-1}, \delta=0.1$ (figure \[fig5\]b). The Brunt-V[ä]{}is[ä]{}l[ä]{} frequencies increase with depth in one case and decrease in the other. It can been seen from figure \[fig5\] that there are, as well, infinitely countable intersections (i.e., solutions) between the two sets of curves. In these cases, contrary to the constant $N$ case, the linear wave solution does not satisfy the nonlinear Boussinesq-Euler equation and hence inclusion of the nonlinear free surface condition is not necessary for the harmonic generation to appear. The analysis for these two cases can be carried out by following the same logic presented in the appendix, although mathematically more involved and closed-form explicit solutions will be tedious, if not impossible, to obtain.

As discussed above, when density profile is not uniform and hence $N$ is not constant, two sources contribute to the harmonic generation: 1- nonlinearities in the free surface boundary condition, i.e. right-hand side of equation , and 2- nonlinear terms in the momentum equation, i.e. right-hand side of equation . As presented before, for a constant $N$ case the latter source is absent. For nonuniform stratifications, it is of interest to evaluate the relative importance of these two contributors. The motivation is to see whether the classical rigid-lid assumption would be a good approximation in estimating energy exchange due to harmonic generation when the density profile is non-uniform (we already presented this is not the case for a uniform stratification).

We comment on this question briefly by considering the exponential density profile as it is representative of actual pycnoclines. For such a density profile, the vertical structure of an internal wave $W(z)$ are obtained as W(z)=- where, ${\mathit J}$ and ${\mathit Y}$ are Bessel function of the first and second kinds respectively, and =,    =,   p=. By substituting this equation into the right-hand side of equations and the relative importance of the right-hand side of the two equations are obtained.

We consider a wave with the wave number and frequency $kh=14.14$ and $\omega/N_\text{max}=0.30$, taken from one solution in figure \[fig5\](b). It turns out that, in this example and for instance (at $z=-20$m), the relative magnitude of nonlinear terms on the right-hand side of the Boussinesq-Euler equation to those of free surface boundary condition (i.e. to ) is about 30%. This means that the contribution of free surface boundary condition to the harmonic generation may be multiple times larger than the nonlinear terms in the momentum equation. We would like to emphasize that the resonant harmonic generation studied here is different from *non-resonant* harmonic generation that obtains in non-uniform stratifications [e.g. @Sutherland2016]. The former occurs for specific internal waves, transfers the entire energy to the second harmonic and is a one-way process, whereas the latter obtains for all internal waves but since it is non-resonant only transfers a portion of energy from the initial wave to the superharmonics.

![Plots of $\D(k,\omega)=0$, $\D(2k,2\omega)=0$ and $\D(k/2,\omega/2)=0$ for parabolic and exponential density profiles. (a) The density is $\rho(z)=\rho_0 (1+a z^2)$, where $a=1\times10^{-5} \text{m}^{-2}$. Water depth is $h=100$m. The Brunt-V[ä]{}is[ä]{}l[ä]{} frequency $N=\sqrt{-2gaz}$. $N_\text{max}=\sqrt{2gah}$. (b) The density is $\rho(z)=\rho_0 [1+\delta-\delta \exp(az)]$, where $a=0.1 \text{m}^{-1}, \delta=0.1$. Water depth is $h=100$m. The Brunt-V[ä]{}is[ä]{}l[ä]{} frequency $N=\sqrt{ga\delta\exp(az)}$. $N_\text{max}=\sqrt{ga\delta}$.[]{data-label="fig5"}](1.eps "fig:"){width="8cm"} (-235,82)[(a)]{}\
![Plots of $\D(k,\omega)=0$, $\D(2k,2\omega)=0$ and $\D(k/2,\omega/2)=0$ for parabolic and exponential density profiles. (a) The density is $\rho(z)=\rho_0 (1+a z^2)$, where $a=1\times10^{-5} \text{m}^{-2}$. Water depth is $h=100$m. The Brunt-V[ä]{}is[ä]{}l[ä]{} frequency $N=\sqrt{-2gaz}$. $N_\text{max}=\sqrt{2gah}$. (b) The density is $\rho(z)=\rho_0 [1+\delta-\delta \exp(az)]$, where $a=0.1 \text{m}^{-1}, \delta=0.1$. Water depth is $h=100$m. The Brunt-V[ä]{}is[ä]{}l[ä]{} frequency $N=\sqrt{ga\delta\exp(az)}$. $N_\text{max}=\sqrt{ga\delta}$.[]{data-label="fig5"}](2.eps "fig:"){width="8cm"} (-235,82)[(b)]{}\

Similarly, since the linear wave solution is not the exact solution to the fully nonlinear problem if taking into the non-Boussinesq terms in the Euler equation, the instability mechanism can also be initiated without a nonlinear free surface.

We also would like to note that non-Boussinesq effects (and corresponding terms in the governing equation) can also lead to the generation of resonant super harmonics. In other words, harmonic generation can be obtained in uniform stratification with the rigid lid, if non-Boussinesq terms are not neglected. However, as will be shown shortly, the relative importance of such routing of energy is orders of magnitude less than that of nonlinear free surface avenue: Following the same procedure for getting , except not making Boussinesq approximation, we obtain,

&w+\_H(w)- (u)- (v)+g\_H\
&+2+ +2(w)+(w)-()\
&-\[()\]-()=0.

Comparing the Boussinesq terms in the above equation and the quadratic terms on the free surface boundary condition, i.e., , we find that their ratio is order of $a/k$. For the two cases presented above, the ratios are less than 1.6%.\
The harmonic generation through a nonlinear free surface involves two energy transfer processes simultaneously, in one process $(k,\omega)$ wave loses energy to $(2k,2\omega)$ wave and generates the second harmonic; in the other the $(2k,2\omega)$ wave forms a triad resonance with $(k,\omega)$ wave and sends energy back in an opposite way. The triad resonance condition is clearly satisfied since $2\omega-\omega=\omega$ and $2k-k=k$. If at the initial moment we only have $(k,\omega)$ wave (i.e. no second harmonic in the domain) then harmonic generation triumphs over triad resonance in transferring energy, resulting in the combined effect that energy only goes in one direction. But it can be shown that there are initial conditions combinations under which the energy at first goes from $(2k,2\omega)$ wave to the $(k,\omega)$ wave, but even in that case eventually the entire energy is transferred to the $(2k,2\omega)$ wave. This can be shown rigorously through Lyapunov stability theorem. Basically, the objective is to prove that $\A_1=0$ and $\B_2$ reaching its maximum is the asymptotically stable solution of the dynamical system described by equation .

To prove this, we first note that can be put in the following form \_2\^2 - \_1\^2 = 0. Therefore, since $\f{\alpha}{\beta}$ can be shown to be always negative(see Appendix B), we have \[3457\] \_2\^2 - \_1\^2 = \^2, where $\G$ is a constant that is determined by the initial conditions. Equation shows that $(\A_1,\B_2)$ are always moving on an ellipse. If we define a Lyapunov function $V(\A_1) = \G  - \alpha/|\alpha|[ \int_{0}^{x}\alpha \A_1(x)^2 dx+\B_2(0)]$, we find that $V(\A_1)\geqslant 0$ and $d V/d x<0$. According to Lyapunov asymptotic stability theorem, the solution of the system converges to $\A_1=0$ from any starting points as $x$ goes to $\infty$.

Conclusion
==========

Here we reported that certain internal gravity waves are inherently unstable and are not able to sustain their form. This new instability mechanism, a result of resonance harmonic-generation, draws the energy of an internal wave and hands it over to its second (or generally higher) harmonic. This resonance is distinguished from the classical triad resonance (and associated subharmonic instability) in that 1- a single wave may undergo the instability without requiring any external perturbation, and 2- the transfer of energy is irreversible and the original wave permanently loses its energy to its second harmonic. Extension of the results presented here to the third and higher harmonic generation is straightforward, but the strength of energy exchange in higher harmonics is much weaker.

Derivation of the Interaction Equation
======================================

Consider the propagation of waves in an inviscid, incompressible, adiabatic and stably stratified fluid of density $\rho(x,y,z,t)$, bounded by a free surface on the top and a rigid seafloor with a depth $h$ at the bottom. We consider a Cartesian coordinate system with $x,y$-axes on the mean free surface and $z$-axis positive upward. Equations governing the evolution of the velocity vector $\b{u} = \{u, v, w\}$, density $\rho$, pressure $p$ and surface elevation $\eta$ under Boussinesq approximation read &\_0=-p-g z,  -h&lt;z&lt;\[g101\]\
&=0,  -h&lt;z&lt;\[g102\]\
&=0,   -h&lt;z&lt;\[g103\]\
&\_[t]{}=(z-),  z=\[g104\]\
&=0,  z=\[g105\]\
&w=0,  z=-h,\[g106\] where $g$ is the gravitational acceleration, and $\rho_0=\bar{\rho}(z=0)$ is the mean density on the free surface with $\bar{\rho}(z)$ the background (unperturbed) density such that $\rho = \bar{\rho}(z)+ \rho' (x, y, z,t)$. Equation is the momentum equation (Euler’s equations), comes from conservation of salt, and is continuity equation that together form five equations for five unknown variables of the problem (three components of velocity ${\bf u}$, pressure $p$ and density $ \rho$). Equations - are boundary conditions on the free surface and the bottom.

Similar to the density perturbation $\rho'$, we define a pressure perturbation $p'$ via $p=\bar p(z)+p'(x,y,z,t)$ such that $d\bar p(z)/dz=-\bar{\rho}(z) g$. For the linear terms of the governing equation to be only in terms of $w$ we calculate $\p/\p z[\g \cdot (\text{\ref{g101}})]-\lap (\text{\ref{g101}})_3$, where (\[g101\])$_3$ denotes the $z$ component of (\[g101\])(likewise, 1,2 for $x,y$ will be used later). We obtain \[2011\] w-u-v+\_H\^2(w)+\_H\^2  ’ =0, where $\lap_{H}=\p^2/\p x^2+\p^2/\p y^2$ is the horizontal Laplacian. Taking the time derivative of and substituting $\rho'$ from the expansion of , i.e., \[g202\] +’+w=0. and denoting $N^2=-{g/\rho_0 \d \bar{\rho}(z)}/{ \d z}$ as the Brunt-V[ä]{}is[ä]{}l[ä]{} frequency, we obtain \[401\] w+N\^2 \_[H]{}w= u+v-\_H\^2(w)+ \_[H]{}(’). Expanding the surface dynamic boundary condition and keeping terms up to the second order we obtain \[302\] =w\_0 g-p’-+\_0 g+ gw ,     [at]{}  z=0. where the last term can be equivalently written as $-N^2\bar\rho_0\eta w$. Here $p'$ can be substituted from an expression obtained from taking the $x$ derivative of (\[g101\])$_1$ added to the $y$ derivative of (\[g101\])$_2$: \[301\] -+\_[H]{}+u+v=0. Substituting $\p p'/\p t$ from in we obtain \[402\] -g\_[H]{}w=& u+v\
&+\_[H]{}g+N\^2w-\_[H]{}p’+.

Therefore the governing equation and boundary conditions correct to $\O(\ep^2)$ reduce to , and . For the ease of referring, we rewrite these equations here: w+N\^2 \_[H]{}w=& u+v-w\
&-w+ \_[H]{}’    -h&lt;z&lt;0,\[3001\]\
-g\_[H]{}w=& u+v+\_[H]{}g+N\^2 w\
& - \_[H]{}p’+     z=0,\[3002\]\
w=&0,  z=-h.\[3003\] Equations and are identical to equations (A2) and (A7) of [@Thorpe1966][^2].

To perform a weakly nonlinear analysis, we assume that internal waves in the system described above are small perturbations from the mean state of water at rest, i.e. $\bf u,\rho',p'\sim {\mathcal O}(\epsilon)$, where $\epsilon$ is a measure of the wave steepness. Considering the two-dimensional problem and assuming that the solution to this problem can be expressed in terms of a convergent series we define w(x,x\_1,z,t)=w\_1(x,x\_1,z,t)+\^2 w\_2(x,x\_1,z,t)+Ø(\^3) where $x_1=\epsilon x$ is the slow spatial variable and $\ep\ll1$ is a measure of steepness of the waves and $w_i\sim\O(1)$. Similar expressions exist for other variables, i.e. $u=\ep u_1+\ep^2u_2+\O(\ep^3)$, $\rho'=\ep \rho'_1+\ep^2\rho'_2+\O(\ep^3)$, $p'=\ep p'_1+\ep^2p'_2+\O(\ep^3)$ and $\eta=\ep \eta_1+\ep^2\eta_2+\O(\ep^3)$ with $u_i,\rho'_i,p'_i\sim\O(1)$ being functions of $x,x_1,z,t$, and $\eta_i\sim\O(1)$ being functions of $x,x_1$ and $t$.

Upon substitution into the governing equation, at the leading order $\O(\ep)$ we obtain &w\_1+N\^2 \_[H]{}w\_1=0 &-h&lt;z&lt;0,\
&-g\_[H]{}w\_1=0& z=0,\
&w\_1=0, &z=-h.

At the second order $\O(\ep^2)$ we have

w\_2+N\^2 \_[H]{}w\_2=&-2+N\^2 w\_1+\_1u\_1\
& +\_1v\_1-\_1w\_1-\_1w\_1\
&+ \_[H]{}(\_1’\_1), -h z 0.\[121\]\
-g\_[H]{}w\_2=&2 g + \_1u\_1+\_1v\_1\
&+\_[H]{}g\_1+N\^2 w\_1 \_1\[122\]\
&- \_[H]{}\_1p’\_1 +\_1, z=0.\
w\_2=&0, z=-h.

We now assume that waves with wavenumber and frequency $(k,\omega)$ and $(2k,2\omega)$ satisfy the internal waves dispersion relation (which is obtained from the linear equation ) (k,)=\^2-kh ,

i.e. $\D(k,\omega)=0$ and $\D(2k,2\omega)=0$, and that they both exist in our domain of interest, though potentially with different amplitudes. The propagating wave solution to the linear equation then obtains as w\_1(x,x\_1,z,t)=&A\_1(x\_1)m\_1(z+h)+B\_1(x\_1)m\_1(z+h)\
+&A\_2(x\_1)m\_2(z+h)+B\_2(x\_1)m\_2(z+h) in which $A_1,A_2,B_1,B_2$ are amplitudes of each wave, $m_1^2=k^2(N^2-\omega^2)/\omega^2$ and $m_2^2=k^2(N^2-4\omega^2)/\omega^2$. Other variables $u,\rho',p'$ and $\eta$ can be found respectively via continuity equation , energy equation , kinematic surface boundary condition and the dynamic free surface boundary condition .

The left-hand side of the second order equation is identical in the form to the first order equation , but the right hand side of is clearly non-zero and is a nonlinear function of the leading order solution $\u_1,\rho_1',p'_1,\eta_1$. It turns out, after substitution, that the right hand side contains terms with harmonics which are the same as the harmonics of the leading order equation (secular terms). A compatibility condition then must be enforced to make sure that the solution does not go unbounded, which is clearly unphysical. This compatibility condition determines the spatial behavior of the coefficients $A_i,B_i$.

While the formulation presented here is general, our primary interest is when an initial wave with wavenumber and frequency $(k,\omega)$ *resonates* its second harmonic $(2k,2\omega)$ whose initial amplitude is zero. Therefore in the following we use the adjectives *original* and *resonant* waves to refer to $(k,\omega)$ and $(2k,2\omega)$ waves respectively. We would like to emphasize that the formulation is general and works for any initial condition of the two waves. We will also comment that the presented approach can be easily extended for third and higher-harmonic generation.

Without loss of generality we assume that $B_1=0$, which only has to do with our choice of coordinate system. But we keep $A_2$ and $B_2$, since they determine the phase of the resonant wave $(2k,2\omega)$ with respect to the original wave $(k,\omega)$. The general solution to the second order problem takes the form w\_2(x,x\_1,z,t)=&C\_[11]{}(x\_1,z)   +C\_[12]{}(x\_1,z)\
+&C\_[21]{}(x\_1,z)+C\_[22]{}(x\_1,z) where $C_{ij}(x_1,z)$’s are to be determined from . Substituting and into , and collecting same sine and cosine terms we obtain four ordinary differential equations for $C_{ij}$ ($i,j$=1,2): -\^2&C\_[1i,zz]{}-m\_1\^2\^2C\_[1i]{}=E\_[1i]{}, & -h&lt;z&lt;0,\
-\^2&C\_[1i,z]{}+gk\^2C\_[1i]{}=F\_[1i]{} & z=0,\
&C\_[1i]{}=0, & z=-h, -4\^2&C\_[2i,zz]{}-4m\_2\^2\^2C\_[2i]{}=E\_[2i]{}, & -h&lt;z&lt;0,\[161\]\
-4\^2&C\_[2i,z]{}+4gk\^2C\_[2i]{}=F\_[2i]{}, & z=0,\[162\]\
&C\_[2i]{}=0, & z=-h,\[163\] where $E_{i1},F_{i1}$ are the coefficients of $\sin(ikx-i\om t)$, and $E_{i2},F_{i2}$ are the coefficients of $\cos(ikx-i\om t)$ in the right-hand side of and respectively. Let’s first consider the equation for $C_{22}$, we obtain E\_[22]{}=-m\_2(z+h) ,     F\_[22]{}=4gkm\_2h for which the solution to that satisfies the boundary condition is C\_[22]{}(x\_1,z)=-(z+h)m\_2(z+h) . Upon substitution into we obtain gk =0 therefore, since the coefficient is nonzero then $\d A_2/\d x_1\equiv$0. Physically speaking, this expression says that the amplitude $A_2$ does not change as waves propagate, or in other words, $A_2$ does not take part in the energy exchange. We, therefore, set $A_2$ equal to zero for the rest of the derivation.

We use the same approach as above for $C_{21}$. We have &E\_[21]{}=m\_2(z+h) ,\
&F\_[21]{}=-4gkm\_2h -23m\^2+\^2 m\_1hA\_1\^2(x\_1) for which C\_[21]{}(x\_1,z)=(z+h)m\_2(z+h) , and upon substitution into we obtain =A\_1\^2(x\_1), where =-.

For $C_{11}(x_1,z)$, we obtain $E_{11}=F_{11}=0$ and therefore the equation for $C_{11}(x_1,z)$ does not provide any extra information on $A_1(x_1)$ and $B_2(x_1)$. For $C_{12}$ we obtain E\_[21]{}=&m\_1(z+h)\
&+ I\^[+]{}(m\_1+m\_2)(z+h)+I\^[-]{}(m\_1-m\_2)(z+h)A\_1(x\_1) B\_2(x\_1)\
F\_[21]{}=&2gkm\_1(z+h) +J A\_1(x\_1)B\_2(x\_1) where I\^[+]{}&=(m\_2+2m\_1)(m\_2\^2+3k\^2-m\_1\^2)\
I\^[-]{}&=(2m\_1-m\_2)(m\_1\^2-3k\^2-m\_2\^2)\
J&=. We obtain C\_[12]{}(x\_1,z)=&-(z+h)m\_1(z+h)\
&++A\_1(x\_1) B\_2(x\_1). Substituting into , we obtain =A\_1(x\_1)B\_2(x\_1), where =-.

If we define the actual amplitudes $\A_1=\ep A_1$ and $\B_2=\ep B_2$ (note that $w=\ep w_1+\O(\ep^2)$) then &=\_1\^2\
&=\_1\_2, with the same $\alpha,\beta$ repeated here: &=-,\
&=-.

Proof of the Sign of $\alpha/\beta$
===================================

For the sign of $\alpha/\beta$, we have $$\text{sign}(\frac{\alpha}{\beta})=\frac{\cos m_2 h}{\sin m_2 h (4 m_1 \cos^2 m_1 h-3 m_1)-2 m_2 \cos m_2 h \sin 2 m_1 h}.$$

It is equvalent to the sign of $$\tan m_2 h (4 m_1 \cos^2 m_1 h-3 m_1)-2 m_2  \sin 2 m_1 h.$$ Rearranging terms of the above expression gives, $$\label{10001}
\cos^2 m_1 h (4 m_1\tan m_2 h -4 m_2\tan m_1 h )-3 m_1\tan m_2 h,$$

After substitution of the disperion relation, $$\omega^2=\frac{gk^2}{m_1}\tan m_1h$$ and $$\omega^2=\frac{gk^2}{m_2}\tan m_2h,$$

can be simplified to

$$-3 m_1\tan m_2 h=-3 m_1 m_2\frac{\omega^2}{gk^2}$$ which is less than zero always.

[10]{}

C. Staquet and J. Sommeria. . , 34:559–593, 2002.

Raffaele Ferrari and Carl Wunsch. Ocean circulation kinetic energy: Reservoirs, sources, and sinks. , 41(1):253, 2008.

Philip W Boyd. Biogeochemistry: iron findings. , pages 10–11, 2007.

Graham Harris. . Springer Science & Business Media, 2012.

Matthew H Alford, Thomas Peacock, Jennifer A MacKinnon, Jonathan D Nash, Maarten C Buijsman, Luca R Centuroni, Shenn-Yu Chao, Ming-Huei Chang, David M Farmer, Oliver B Fringer, et al. The formation and fate of internal waves in the south china sea. , 521(7550):65–69, 2015.

HP Zhang, Ben King, and Harry L Swinney. Resonant generation of internal waves on a model continental slope. , 100(24):244504, 2008.

Russ E Davis and Andreas Acrivos. The stability of oscillatory internal waves. , 30(04):723–736, 1967.

K Hasselmann. A criterion for nonlinear wave stability. , 30(04):737–739, 1967.

S. a. Thorpe. . , 263(1145):563–614, December 1968.

C Henry McComas and Francis P Bretherton. Resonant interaction of oceanic internal waves. , 82(9):1397–1412, 1977.

Chung-Hsiang Jiang and Philip S Marcus. , 102(12):1–4, 2009.

H[é]{}l[è]{}ne Scolan, Eugeny Ermanyuk, and Thierry Dauxois. Nonlinear fate of internal wave attractors. , 110(23):234501, 2013.

M.-R. Alam, Y. Liu, and D. K. P. Yue. . , 624:225–253, March 2009.

Mohammad-Reza Alam, Yuming Liu, and Dick KP Yue. Resonant-wave signature of an oscillating and translating disturbance in a two-layer density stratified fluid. , 675:477–494, 2011.

Scott Wunsch. . , 71:91–97, 2015.

S. A. Thorpe. . , 24:737–751, mar 1966.

S. Martin, W. Simmons, and C. Wunsch. . , 53:17–44, March 1972.

Ali Tabaei, T. R. Akylas, and Kevin G. Lamb. . , 526:217–243, 2005.

Mohammad-Reza Alam, Yuming Liu, and Dick KP Yue. Oblique sub-and super-harmonic bragg resonance of surface waves by bottom ripples. , 643:437–447, 2010.

Mohammad-Reza Alam. A new triad resonance between co-propagating surface and interfacial waves. , 691:267–278, 2012.

Mohammad-Reza Alam and Chiang C Mei. Attenuation of long interfacial waves over a randomly rough seabed. , 587:73–96, 2007.

O. B. Fringer, M. Gerritsen, and R. L. Street. . , 14(3-4):139–173, 2006.

Z. Zhang, O. B. Fringer, and S. R. Ramp. . , 116(C05022):1–26, 2011.

D. Kang and O. Fringer. . , 40(11):2539–2545, 2010.

B. Wang, S. N. Giddings, O. B. Fringer, E. S. Gross, D. A. Fong, and S. G. Monismith. . , 116(C02036):1–23, 2011.

D. Kang and O. Fringer. . , 42(2):272–290, 2012.

R. K. Walter, C. B. Woodson, R. S. Arthur, O. B. Fringer, and S. G. Monismith. . , 117(C07017):1–13, 2012.

Theo Gerkema, Chantal Staquet, and Pascale Bouruet-Aubertot. . , 33(8):L08604, 2006.

J. A. MacKinnon and K. B. Winters. . , 32(15):L15605, 2005.

Peter Muller. . , 24(3):493–536, 1986.

Bruce R Sutherland. . , 793:335–352, 2016.

[^1]: reza.alam@berkeley.edu

[^2]: except that we found three typos there: (1) for $S1$ in (A2), the sign of the last term with $g$ should be positive; (2) for $S3$ in (A7), the first term should be $\pmb{\omega}  \cdot \nabla w$; and (3) the term $g \nabla^2_1\frac{\partial}{\partial x}(\eta w_{z}) $ is missing. These are clearly typos as they do not appear in later expressions of [@Thorpe1966].
