Instabilities in internal gravity waves

1.   Introduction

Internal gravity waves are propagating disturbances in stably stratified fluids, with buoyancy being the restoring force mechanism ^[1,2]. They can transport momentum and energy to large distances both horizontally and vertically, and are hence a significant consideration in energy budgets of stratified flows. From a fundamental point of view, internal waves possess several interesting properties that distinguish them from several other waves in fluids. For example, the linear internal wave dispersion relation shows that the phase and group velocities of plane internal waves are orthogonal to each other, and the direction of energy propagation is independent of the wave vector magnitude. Classical wave phenomena like scattering ^[3,4], diffraction ^[5] and interference ^[6,7], therefore have interesting manifestations in internal waves. Flow instabilities represent another classical phenomenon, which in addition to transferring energy to different spatial scales, redistributes internal wave energy to different temporal scales too. In this paper, we review the existing knowledge on instabilities in internal gravity waves that occur in continuous stable stratifications, with a focus on oceanic applications ^[8,9].

Internal waves are ubiquitous in the temperature- and salinity- stratified oceans, owing to continuous input of energy into them from tides and winds. While barotropic tides interact with ocean floor topography to generate internal waves at tidal frequencies ^[10], the action of wind on the upper ocean mixed layer generates near-inertial internal waves ^[11]. The details of these internal wave generation mechanisms are reasonably well understood ^[12,13], and from a mathematical viewpoint, perturbation methods, ray tracing and Green function approaches have been insightful to study topographic effects in internal waves ^[14,15,16]. Around 1TW of energy is understood to be irreversibly put into each of internal tides and wind-generated near-inertial waves, respectively ^[17], which are thought to be sufficient to cause the required deep ocean heat addition for maintaining the global circulation ^[18]. As a result, there are several ongoing efforts on understanding the spatio-temporal variability and mechanisms associated with internal wave dissipation. In addition, internal waves are so prevalent that they are a significant part of the ocean state even far from their generation sites ^[19]. Finally, an accurate parameterization of sub-grid scale internal wave processes and associated turbulence & mixing is critical for climate models ^[20]. Instabilities in internal waves, the focus of this review, represent one of the important pathways towards dissipation ^[12,21].

Internal waves, due to the dynamical relevance of both shear and buoyancy in them, are interesting from an instabilities viewpoint. Also, since they do not fit classical one-dimensional descriptions, the relevance of gravitational ^[22] and shear ^[23,24] instabilities is unclear a priori. In addition, time dependence in shear and buoyancy, such as those that occur in internal waves, can nontrivially influence the dynamics and instability characteristics ^[25,26,27]. Wave-wave interactions based on weakly nonlinear expansions ^[9,28,29] have been the primary approach towards understanding instabilities in internal waves though linear stability calculations based on Floquet theory ^[30] have provided significant insights too. Several studies using one of the two aforementioned approaches have focused on plane internal waves of infinite spatial extent ^[31]. In the ocean, however, internal wave energy is often put into a range of wavenumbers, even if at the same frequency. This results in the generation of internal wave beams of finite spatial extent, the instability characteristics of which have received much attention in the last two decades ^[32]. Finite-depth effects are important too, both in the ocean ^[1] and laboratory experiments ^[33,34], and internal wave instabilities in finite-depth domains are also an active area of ongoing research. With recent advances in detailed experimental and numerical modeling, and significant improvements in sophisticated field measurements, mathematical advances on internal wave instabilities under various settings are necessary for better understanding and subsequently driving the field forward.

In this review, we discuss internal wave instabilities in plane waves, wave beams and modes, along with a pedagogical discussion of linear internal waves of different spatial forms and basic ideas of triadic resonance. Section 2 introduces the internal wave dispersion relation, along with a description of different spatial forms of internal waves. Section 3 presents an introductory discussion of triadic resonance in internal waves. In Section 4, the current knowledge of instabilities in internal waves is reviewed, with plane waves, wave beams and modes being discussed separately. An ocean perspective is provided in Section 5, followed by concluding remarks in Section 6.

2.   Internal waves: governing equations and dispersion relation

In this section, we present the governing equations, followed by the internal wave dispersion relation and a description of various fundamental spatial forms in which internal waves can occur.

2.1.   Governing equations

The inviscid governing equations for an incompressible flow within the Boussinesq approximation are ^[35]

where ${\bf{u}}$ is the velocity field of the flow, $t$ is time, $p$ the pressure deviation from the hydrostatic pressure distribution, $\rho_{ref}$ a constant reference density and $b = (\rho_0(z)-\rho)g/\rho_{ref}$ is the buoyancy field. Here, $\rho$ is the total density field, $\rho_0(z)$ the background density profile (assumed stable, i.e., $d\rho_0/dz < 0$, throughout this paper), and $-g{\hat{{\bf{e}}}_{\bf{z}}}$ being the gravity vector. The density differences are assumed to satisfy $|\rho-\rho_{ref}|\ll\rho_{ref}$ and $|\rho_0-\rho_{ref}|\ll\rho_{ref}$ so that the Boussinesq approximation is valid. In Eq (2.3), $w$ is the velocity component along ${\hat{{\bf{e}}}_{\bf{z}}}$, the vertical direction. The Brunt-Väisälä frequency $N$ (also referred to as the buoyancy frequency) is related to the background density profile by $N(z) = \sqrt{(-g/\rho_{ref})(d\rho_0/dz)}$. In summary, Eqs (2.1)–(2.3) are the mass conservation (in the incompressible limit), momentum equation and incompressibility condition, respectively.

Assuming the flow to be two-dimensional, the velocity field can be described in terms of a streamfunction $\psi$ as ${\bf u} = u {\hat{{\bf{e}}}_{\bf{x}}} + w {\hat{{\bf{e}}}_{\bf{z}}} = (-\partial\psi/\partial z, \partial\psi/\partial x)$, where $x$ and $z$ are the horizontal and vertical coordinates, respectively. Equations (2.2)–(2.3), upon eliminating $p$, can now be re-written as ^[36]

The Jacobian operator $J$ is defined as $J(A, B) = (\partial A/\partial x)(\partial B/\partial z) - (\partial B/\partial x)(\partial A/\partial z)$, and $\nabla^2 = \partial^2/\partial x^2 + \partial^2/\partial z^2$ is the Laplacian operator. Equations (2.4)–(2.5) represent a useful form since the linear and nonlinear terms (in $\psi$ and $b$) have been separated as the left and right hand sides, respectively. Equations (2.4)–(2.5) can further be recast as

from which the linear internal wave equation (in terms of $\psi$ alone) is recovered in the limit of neglecting the nonlinear terms on the right hand side.

2.2.   Dispersion relation

In this section, we discuss the internal wave dispersion relation, as obtained from the linearized form of Eq (2.6), in the limit of a uniform stratification ($N(z) = N_0$, where $N_0$ is a constant) in an unbounded domain. To obtain the dispersion relation, the nonlinear terms in the right hand side of Eq (2.6) are ignored, and a monochromatic plane wave solution of the form

is sought, where $\Psi$ is a constant complex amplitude, c.c. denotes complex conjugate and ${\bf{x}} = x{\hat{{\bf{e}}}_{\bf{x}}} + z{\hat{{\bf{e}}}_{\bf{z}}}$ is the position vector. Here, ${\bf{k}} = k_x {\hat{{\bf{e}}}_{\bf{x}}} + k_z{\hat{{\bf{e}}}_{\bf{z}}}$ is the wave vector and $\omega$ (assumed real) is the frequency. The dispersion relation then follows as ^[36]

where $\theta = \tan^{-1}{(k_x/k_z)}$. In the wave propagation regime ($\omega \le N_0$), $\theta$ is real and represents the angle that the wave vector ${\bf{k}}$ makes with the vertical $z-$axis. As $\omega$ varies from 0 to $N_0$, $\theta$ varies from 0 to $\pi/2$. For $\omega > N_0$, referred to as the evanescent regime, $k_x/k_z$ is imaginary, hence corresponding to exponential decay of the wave in space.

The internal wave dispersion relation in Eq (2.8) has several interesting properties. Firstly, for a given $N_0$, the frequency $\omega$ specifies (via the dispersion relation) only the wave vector orientation $\theta$, and not the wave vector magnitude $|{\bf{k}}|$. In addition, for given values of $\omega$ and $N_0$, four possible solutions exist for $\theta$, with each of the corresponding wave vectors being in a different quadrant on the $(k_x, k_z)$ plane. Interestingly, the phase velocity ${\bf{c}}_{\bf{p}} = (\omega/|{\bf{k}}|^2){\bf{k}}$ and the group velocity ${\bf{c}}_{\bf{g}} = (\partial\omega/\partial k_x, \partial\omega/\partial k_z)$ are orthogonal to each other and have opposite signs in their $z-$components for plane internal waves. In other words, $\theta$, in addition to representing the wave vector orientation with respect to the $z-$axis, also represents the energy propagation direction with respect to the horizontal $x-$axis. Physically, small ($\gtrsim 0$) and large ($\lesssim N_0$) values of $\omega$ correspond to shallow and steep internal waves, respectively, with respect to the energy propagation direction.

2.3.   Internal wave forms

In the ocean or the laboratory, the forcing mechanism for internal waves is often at one (or sometimes a few) specific frequency. For example, the principal lunar semidiurnal (so called M2) frequency contains the most energy for internal tides in the ocean; more details are given in Section 5. As seen in Section 2.2, fixed values of $N_0$ and $\omega$ determine the internal wave orientation via the dispersion relation, whereas the specific wavenumbers that are excited depend on spatial scales that are imposed by the forcing mechanism or the domain boundaries. While energy is put into a continuous range of wavenumbers in oceanic settings, idealized theoretical and experimental studies on individual wavenumbers have provided significant insights. Furthermore, the finite depth nature of the ocean or laboratory set-ups reduce the possible wavenumbers to a discrete set. In this subsection, we describe four different spatial internal wave forms, which have either been considered in idealized theoretical and experimental studies or observed in realistic ocean settings. Sections 3 & 4 will then focus on various instabilities that could occur in these internal wave forms.

2.3.1.   Plane wave

A plane internal wave is described by Eq (2.7), characterized by unique values of frequency $\omega$, wavenumbers ($k_x, k_z$) and amplitude $\Psi$. As described in Section 2.2, the phase velocity of a plane wave is along $k_x{\hat{{\bf{e}}}_{\bf{x}}} + k_z{\hat{{\bf{e}}}_{\bf{z}}}$, whereas its group velocity is along ${\rm{sgn}}(k_xk_z)(k_z {\hat{{\bf{e}}}_{\bf{x}}} - k_x {\hat{{\bf{e}}}_{\bf{z}})}$. Figure 1(a) shows an experimental realization of a plane wave, though a continuous range of wavenumbers around a desired values is excited ^[37]. Specifically, Gostiaux et al. ^[37] used a novel internal wave generator in a uniform stratification to excite a single frequency, with the dominant wavenumber being fixed by the wavelength in the wave generator (see the wave form shown in blue in Figure 1(a)). While the energy propagation direction (along ${\bf{c}}_{\bf{g}}$) is fixed by $\omega$ and $N_0$, only one of the two possible solutions for $\theta$ (note that causality allows only for left-to-right propagating waves) in Eq (2.8) is predominant as a result of the upward phase velocity associated with the moving wave form imparted by the wave generator. Owing to the finite-width nature of the forcing from the wave generator, the spatial Fourier spectrum of the generated wave field shows a spread around (instead of a sharp peak at) the dominant wavenumber. More details of plane wave generation by the novel internal wave generator, including the spectrum in the wavenumber space, can be found in ^[38]. Theoretical/numerical studies on exact plane waves, and experimental studies on quasi-plane-waves such as the one shown in Figure 1(a), are discussed in section 4.1 from the perspective of their instability characteristics.

2.3.2.   Wave beam

Unlike a quasi-plane-wave discussed in Section 2.3.1, the internal wave field generated by typical finite-size forcing mechanisms contain a non-negligible amount of energy in a relatively wider range of wavenumbers. For example, a cylinder oscillating at a fixed frequency in a uniform stratification generates four internal wave beams (Figure 1(b)), each of which has an orientation $\theta$ that is consistent with the dispersion relation in Eq (2.8). Each wave beam is a combination of plane waves over a continuous range of wavenumbers, resulting in the energy being present only over a finite width (along the common wave vector direction) in space. The cylinder diameter imposes a length scale on the finite width of the wave beams. Several studies, following the pioneering experiments of ^[39], have investigated the internal wave field around an oscillating cylinder ^{[40,41,42,43]}. In general, the stream function associated with a unidirectional internal wave beam, whose energy propagation direction (along ${\bf{c}}_{\bf{g}}$ in figure 1b) makes an angle $\theta$ with the horizontal, can be written as

where ${\bf{k}} = |{\bf{k}}|(\sin\theta$ ${\hat{{\bf{e}}}_{\bf{x}}} - \cos\theta$ ${\hat{{\bf{e}}}_{\bf{z}})}$, and the amplitude distribution $\Psi(|{\bf{k}}|)$ depends on the forcing characteristics. A discussion of instabilities in internal wave beams, and how they may differ from instabilities in plane waves, is presented in section 4.2.

2.3.3.   Finite-depth effects

Introducing horizontal boundaries, say at $z = 0$ and $z = H$, discretizes the possible wavenumbers owing to the boundary conditions. Specifically, for no-normal-flow boundaries at $z = 0$ and $z = H$ in a uniform stratification, the allowed horizontal wavenumbers are $k_x^n = n\pi/(H\cot\theta)$, with the corresponding vertical mode shape for the stream function being $\phi_n(z) = \sin(n\pi z/H)$, where $n$ is the mode number. An individual mode can then be written as

where $\Psi_n$ is the modal amplitude. The expression in Eq (2.10) can also be constructed as a sum of upward and downward propagating plane waves with the same horizontal wavenumber, suggesting that an individual mode can behave quite differently from an individual plane wave. Figure 1(c) shows an experimental realization of a left-to-right propagating mode $1$ internal wave, again excited using the novel internal wave generator ^[38].

Any internal wave field in a finite-depth stratified fluid satisfying the linearised governing equations can be represented as a summation over one or more individual vertical modes, each having its own amplitude. Such a superposition often results in beam-like features, an example of which is shown in Figure 1(d). In Figure 1(d), the internal wave field generated by barotropic tidal forcing over bottom topography (shown in black) is shown. Internal wave energy is clearly seen to be significant only over finite-width wave beam-like features. In addition to the horizontal boundaries, the introduction of features like a sloping wall ^[45] or a double-ridge bottom topography ^[46] can raise the possibility of internal wave attractors, which are closed loops formed by the internal wave ray paths. The internal wave attractor has also been analyzed using tools from spectral theory and microlocal analysis ^[47]. Instability characteristics of such internal wave attractors, while being of interest in recent studies ^[48,49,50], will not be discussed further in this review.

2.4.   Internal wave forms as nonlinear states

The various internal wave forms discussed in Sections 2.3.1–2.3.3 are solutions of the linear internal wave equations, which ignore the right hand side terms in Eqs (2.6) and (2.5). Interestingly, individual plane waves (Section 2.3.1) and unidirectional wave beams (Eq 2.9) in a uniform stratification are also solutions of the fully nonlinear Eqs (2.6) and (2.5) ^[51]. As a result, plane waves and unidirectional wave beams represent fully nonlinear internal wave states, whose instability characteristics are discussed in Section 4. Similarly, it can be verified that individual modes (see Section 2.3.3) are also solutions of the fully nonlinear internal wave equations ^[52]. For a superposition of modes, however, the nonlinear terms in the right hand sides of Eqs (2.6) and (2.5) do not vanish. For example, two different modes at frequencies $\omega_1$ and $\omega_2$ can give rise to terms containing the frequencies $\omega_1 + \omega_2$ and $\omega_1-\omega_2$ when the nonlinear terms $J(\psi, \nabla^2\psi)$ and $J(\psi, b)$ are evaluated. Similarly, nonlinear terms at sum and difference frequencies can emerge upon interaction between plane waves as well, an aspect we discuss further in Sections 3 and 4.

3.   Triadic resonance

Triadic resonance is a mechanism by which a triad of waves, under certain conditions, can spontaneously transfer energy to each other even at arbitrarily small but non zero amplitudes ^[29]. It is instructive to discuss the origins of triadic resonance based on the nonlinear terms in the governing equations. Assuming a linear wave field consisting of two plane internal waves: $\Psi_1 \exp[i({\bf k_1}\cdot{\bf{x}}-\omega_1 t)] + c.c.$ and $\Psi_2 \exp[i({\bf k_2}\cdot{\bf{x}}-\omega_2 t)] + c.c.$, the nonlinear advection term ${\bf u.\nabla u}$ in the momentum Eq (2.2) (or the Jacobian terms in the right-hand side of Eqs (2.6) and (2.5)) would give rise to terms of the form $\exp[i(({\bf k_1 \pm k_2})\cdot{\bf{x}}-(\omega_1 \pm \omega_2) t)] + c.c.$, which would then act as forcing terms in the governing equations at second order in the linear wave amplitudes. If the wave vector and frequency combination, i.e., (${\bf k_1} + {\bf k_2}, \omega_1 + \omega_2$) or (${\bf k_1} - {\bf k_2}, \omega_1-\omega_2$), in any of the aforementioned terms satisfies the dispersion relation, the governing equation for the second order wave field would have resonant forcing terms.

In general, three plane waves satisfying the following temporal and spatial resonance conditions

form a resonant triad, as weakly nonlinear interaction between any two of the waves resonantly forces the third. Here, $s_1 = \pm1$ and $s_2 = \pm1$. A geometric depiction of the spatial resonance condition in Eq (3.1) with $s_1 = s_2 = 1$ is shown in Figure 2(a). In such instances of triadic resonance, the weakly nonlinear wave field (resulting from interactions between any wave pair within the triad) with constant amplitudes would diverge, and the wave amplitudes in the triad would spontaneously evolve. The corresponding amplitude evolution equations can be derived using techniques like the method of multiple scales ^[29,53] or the variational method ^[9,54]. In the inviscid limit, they take the form

where $\Psi_j(\tau)$ is the complex wave amplitude of the $j^{th}$ wave, and $\tau$ is a slow time coordinate (slow compare to the wave periods). In Eqs (3.3)–(3.4), $s_1 = s_2 = 1$ has been assumed, and the signs of the coefficients $\alpha_j$ (constants, which depend on the values of wavenumbers and frequencies of waves in the triad) determine the stability characteristics of the triad. Such derivations assume that the wave amplitudes evolve over time scales that are much larger than those corresponding to the wave frequencies.

Pioneering studies on triadic resonance were performed within the context of surface gravity waves by Phillipps (1960) ^[55] and the 2021 Nobel prize-winning Hasselmann (1962) ^[56], and a detailed review of the early developments can be found in ^[28]. While we have thus far discussed triadic resonance as a case of when the resonance conditions (3.1)–(3.2) are exactly satisfied, quasi-resonance that occurs in the near vicinity of exact triadic resonance is also an important consideration. Off-resonance amplitude evolution equations have been derived too, as discussed in ^[29]. While quasi-resonance could be associated with relatively small secondary wave growth compared to triadic resonance, it's likelihood of occurrence is potentially larger. Finally, apart from the leading order triadic resonance described by resonance conditions (3.1)–(3.2), higher order resonances are possible too. For example, a primary wave specified by $(\omega_3, {\bf k_3})$ could succumb to higher order resonances containing secondary daughter waves 1 & 2 such that

where $n > 1$ is an integer, are satisfied ^[57].

With respect to internal waves, it is convenient to first assume the frequencies to be always positive. It is noteworthy that this assumption causes no loss of generality as long as the wave vector components are allowed to be positive or negative. With positive frequencies, $s_1 = s_2 = -1$ is ruled out in Eq (3.2). Now, choosing the wave indices such that the triadic resonance conditions are written as $\omega_1+\omega_2 = \omega_3$ (hence $\omega_3$ is the largest frequency) and ${\bf k_1} + {\bf k_2} = {\bf k_3}$, resonant triads can be identified on the plane of ($k_{x, 1}/k_{x, 3}, k_{z, 1}/k_{z, 3}$) as

where ($\omega_1, {\bf k_1}$), ($\omega_2, {\bf k_2}$) and ($\omega_3, {\bf k_3}$) have been assumed to satisfy the dispersion relation in Eq (2.8). For a given $\omega_3/N_0$, and using the dispersion relation $k_{z, 3}^2/k_{x, 3}^2 = (N^{2}_0 - \omega^{2}_3)/\omega^{2}_3$, Eq (3.7) can be plotted on the ($k_{x, 1}/k_{x, 3}, k_{z, 1}/k_{z, 3}$) plane, as shown in Figure 2(b) for $\omega_3/N_0 = 0.7$. Three different branches (denoted A, B, C in Figure 2(b)) are identified in the upper half plane, with equivalent branches present in the lower half plane too. In branch C, $k_{x, 1}$ and $k_{x, 2}$ have the same sign as $k_{x, 3}$, whereas $k_{z, 1}$ and $k_{z, 2}$ are of opposite signs. More importantly, $|{\bf{k_1}}|$, $|{\bf{k_2}}|$ and $|{\bf{k_3}}|$ all correspond to comparable wavelengths in branch C. In contrast, branches A & B have ${\bf{k_1}}$ and ${\bf{k_2}}$ in diagonally opposite quadrants, and the corresponding wavelengths can either be comparable or much smaller than that of ${\bf{k_3}}$. Two sample resonant triads one each in branch C and branch B, are indicated in Figure 2(b), with the blue, magenta and black vectors denoting ${\bf k_1}$, ${\bf k_2}$ and ${\bf k_3}$, respectively. Figure 2(c) shows the variation of $\omega_1/\omega_3$ as we move along the triadic resonance curves in Figure 2(b). At the four corner points, namely $(k_{x, 1}/k_{x, 3}, k_{z, 1}/k_{z, 3}) =$ $(0, 0)$, $(1, -1)$, $(1, 1)$ and $(0, 2)$, either wave $1$ or $2$ has the same frequency as wave 3, with the corresponding wave vector inclination with the $x$-axis being $\theta_3$ or -$\theta_3$, where $\theta_3$ is the angle that ${\bf{k_3}}$ makes with the $x$-axis. Finally, in the limit of $|k_{x, 1}/k_{x, 3}|\to\infty$ (occurring on branches A & B), waves $1$ and $2$ both approach half the frequency of wave $3$, with the corresponding wave vectors becoming anti-parallel to each other.

A resonant triad in any system often manifests in two different ways. The first is triadic resonant instability, where only one of the waves in the triad represents the primary wave field (finite amplitude), and the other two emerge from noise-level amplitudes as the secondary wave field. A triadic resonant instability with both the secondary wave frequencies being smaller than the primary wave frequency is also referred to as subharmonic resonance. The second manifestation of a resonant triad, namely triadic resonance interaction, involves two waves in the triad being present in the primary wave field, and the third being excited as the secondary wave field. When the excited secondary wave field corresponds to a frequency that equals the sum of the primary wave frequencies, it is referred to as superharmonic resonance. A self-interacting primary wave at frequency $\omega$ generating a secondary wave at frequency 2$\omega$ due to triadic resonance is a specific form of superharmonic resonance. It is noteworthy that triadic resonance due to a self-interacting plane internal wave in a uniform stratification is not possible. In fact, it is evident from Figure 2(b) that two plane waves with the same direction for their wave vectors cannot be part of a resonant triad, which is consistent with a unidirectional wave beam being a solution of the fully nonlinear internal wave equation (see Section 2.4).

Subharmonic and superharmonic resonances are instability mechanisms by which internal wave energy at a given frequency gets transferred to other spatial scales and frequencies. In general, which triad becomes dynamically active for a given primary wave field is a non-trivial question, and is likely a function of the stability characteristics of the corresponding amplitude evolution Eqs (3.3)–(3.4). To estimate growth rates, the amplitude evolution equations are often solved within the pump-wave approximation ^[29], wherein it is assumed that the primary wave field remains at constant amplitude during the early evolution of the secondary wave field. Previous theoretical studies ^[30,58] have argued that resonant triad interactions are the basic mechanism behind every instability in a plane internal wave. A more recent theoretical study has also shown that several instabilities in a plane internal wave are related to triadic resonances at various orders ^[59]. In Section 4, we present and discuss the existing knowledge on instabilities in various internal wave forms presented in Section 2.3.

4.   Instabilities in internal waves

In this section, we review existing literature on instabilities in the different spatial forms of internal waves discussed in Section 2.3. While the focus in this section is on idealized settings of single-frequency primary wave fields, a discussion of more realistic oceanic settings is provided in Section 5.

4.1.   Plane waves

Plane waves are internal waves of a specific frequency and wave vector in an infinite medium, as summarized in section 2.3.1. Early theoretical developments on instabilities in plane waves had two different approaches, namely (ⅰ) investigating triadic resonance and (ⅱ) performing a linear stability analysis. The main results from the two approaches are presented in Sections 4.1.1 and 4.1.2.

4.1.1.   Triadic resonance

As discussed in Section 3, all the triadic resonances that a plane wave is part of, can be identified by moving along the curves shown in Figure 2(b). Two different cases are possible for triadic resonant instability of a single primary wave: (ⅰ) the primary plane wave is of the largest frequency in the triad, i.e., the wave denoted by subscript 3, and (ⅱ) one of the secondary waves is of the largest frequency, i.e., wave 1 or 2 is the primary wave. The secondary waves have a sum and difference interaction in the aforementioned cases (ⅰ) and (ⅱ), respectively. Hasselmann (1967) ^[60] has shown that the case (ⅰ) results in an unstable growth of secondary waves, whereas the case (ⅱ) is neutrally stable. As a result, the following discussion is focused on case (ⅰ), where a primary plane wave of a finite initial amplitude resonantly interacts with two secondary waves (infinitesimally small initial amplitudes) of smaller frequencies, i.e., subharmonic resonance.

For a given primary plane wave frequency $\omega_3$, a continuous range of subharmonic frequencies ($\omega_1, \omega_3-\omega_1$) can be found for the secondary wave pair that is in triadic resonance with the primary wave (Figure 2c). Inviscid amplitude evolution equations, along with the corresponding growth rates for the secondary waves, have been derived for all such subharmonic resonances using the method of multiple scales ^[60]. Of all the possible subharmonic resonances, the scenario with $k_{x, 1}/k_{x, 3}\to\infty$, $k_{z, 1}/k_{z, 3}\to\infty$ and $\omega_1/\omega_3\to0.5$, corresponding to the $k_{x, 1}/k_{x, 3}\to\infty$ limit of the branch B of red curves in Figure 2(b) & (c), has the largest growth rate ^[61]. It is noteworthy that such a scenario is equivalent to the $k_{x, 1}/k_{x, 3}\to-\infty$ limit of the green curves in Figure 2(b) & (c). For the limit with $k_{x, 1}/k_{x, 3}\to\infty$, $k_{z, 1}/k_{z, 3}\to\infty$ and $\omega_1/\omega_3\to0.5$, one of the secondary waves is in the same quadrant as the primary wave, with an infinitesimally short wavelength, while the other secondary wave is anti-parallel to the first. This limit has come to be known as the parametric subharmonic instability (PSI), representing a significant physical mechanism by which primary internal wave energy could be transferred to smaller spatial scales and frequencies. As discussed in the later sections, PSI and related phenomena are observed in other internal wave forms and the ocean too.

4.1.2.   Linear stability analysis

The amplitude evolution Eqs (3.3-3.4) discussed in Sections 3 & 4.1.1 consider weakly nonlinear interactions, and hence sufficiently small wave amplitudes, associated with triadic resonance. An alternative approach is to perform a conventional linear stability analysis, which assumes perturbations of sufficiently small amplitude, of a plane internal wave. A plane internal wave, being a solution to the fully nonlinear equations (Section 2.3.1), represents a legitimate base flow irrespective of its amplitude. As a result, linear stability analysis is potentially useful to understand instabilities in finite-amplitude plane internal waves too.

Mied (1976) ^[62] considered linearized inviscid governing equations for two-dimensional perturbations superimposed on a plane internal wave, to derive a pair of coupled partial differential equations with periodic coefficients. Using Floquet theory ^[63], and considering only those perturbations with an integer value for the non-dimensionalized wavenumber along the base flow wave vector, Mied (1976) ^[62] established a connection between internal wave instabilities and the instabilities present in the Mathieu equation. Analyzing the $\theta = 60^0$ plane wave in detail, Mied (1976) ^[62] concluded that plane waves of any arbitrary amplitude are parametrically unstable, and the parametric instability coincides with subharmonic triadic resonance in the limit of infinitesimally small base flow amplitude. Klostermeyer (1982) ^[64], advancing the numerical method proposed by Mied (1976) ^[62] for solving the Floquet system for two-dimensional perturbations, surveyed all possible perturbation wave vector directions to confirm that the most unstable perturbations for small-amplitude plane internal waves correspond to the secondary waves in PSI. Drazin (1977) ^[65] independently reached a similar conclusion that plane internal waves of any amplitude are parametrically unstable, and also showed that the marginal stability boundaries on the space of plane wave amplitude and wavenumbers can be described by catastrophe theory ^[66,67].

Klostermeyer (1991) ^[58] extended the linear stability analysis to three-dimensional perturbations to conclude that three-dimensional instabilities are stronger than two-dimensional ones in finite-amplitude plane internal waves. The three-dimensional instabilities, which were traced to higher order wave interactions including vortical modes, were also suggested to be closely associated with the occurrence of statically unstable regions in the base flow; as a result, static instability is not to be thought of as a mechanism that is independent of wave interactions. Higher order resonances in small-amplitude plane internal waves have been studied using the linear stability framework, with such instabilities having the new possibility to align with the base flow shear ^[57]. More extensive linear stability analyses have further strengthened our understanding of connections between internal wave instabilities and resonant wave interactions at small amplitudes, with complementary energy budget analyses revealing the physical mechanisms underlying various instabilities ^[30,68]. Figure 3(a) reveals one of the important outcomes of linear stability analysis, highlighting the different dominant instabilities on the plane of internal wave amplitude and orientation.

An important conclusion from the aforementioned studies is that several instabilities, both in small- and finite-amplitude internal waves, occur at spatial scales that are smaller than the wavelength of the primary plane wave. Indeed, a two-dimensional study of small-scale perturbations superimposed on the horizontal acceleration of a large-scale internal wave field has revealed the basic characteristics of parametric instability in a plane internal wave ^[69]. Recently, Ghaemsaidi & Mathur (2019) ^[59] performed a local stability analysis ^[70], also known as geometric optics or rapid distortion theory, of a plane internal wave. Such an approach focused on short-wavelength perturbations allows computationally efficient investigations over a wide parameter space spanning both the base flow and perturbation characteristics, and has previously been used to study instabilities in geophysically relevant wave flows ^[71,72,73]. As shown in Figure 3(b), Ghaemsaidi & Mathur (2019) ^[59] delineated the internal wave amplitude and orientation space into various regions of dominant instabilities, and concluded that 2D PSI and 3D transverse instabilities are dominant at small and large internal wave amplitudes, respectively. Importantly, they highlighted the relevance of triadic resonances (of various orders) in small-amplitude internal waves in describing the instabilities of finite-amplitude internal waves. They also showed that three-dimensional transverse instabilities, which are closely related to the occurrence of statically unstable regions in finite-amplitude internal waves, are well-described by the Mathieu equation.

4.1.3.   Numerical simulations and laboratory experiments

Direct numerical simulations have investigated the validity of conclusions reached by triadic resonance and linear stability calculations, and further extended the studies to viscous regimes. Lombard & Riley (1996b) ^[74] performed direct numerical simulations of a plane internal wave at different values of amplitude, orientation and Reynolds number (which quantifies the viscous effects) to report that three-dimensional instabilities are important in finite-amplitude internal waves. Furthermore, they showed that the breakdown to turbulence is driven by both shear and convective instabilities, even at relatively large internal wave amplitudes. Studying the energetics in detail, Koudella & Staquet (2006) ^[75] performed two-dimensional direct numerical simulations of small-amplitude plane internal waves to conclude that parametric subharmonic instability represents an optimal mechanism by which perturbations extract energy from an internal wave that contains oscillating shear and density gradient. Figure 4 shows buoyancy-driven instability at early times, followed by shear-driven Kelvin Helmholtz instability at later times in the two-dimensional direct numerical simulations of Koudella & Staquet (2006) ^[75]. Recently, direct numerical simulations have explored the viscous, nonlinear regimes, and the associated turbulent mixing, that follow local instabilities in a plane internal wave ^[76].

Laboratory experiments ^[77] have also revealed the occurrence of subharmonic triadic resonance in plane internal waves, whose generation in an experiment improved with the advent of the novel internal wave generator ^[37,38]. Using a quasi-plane-wave, containing around two to three wavelengths associated with a specific wavenumber, as the primary wave (see Figure 5 (a)), Bourget et al. (2013) ^[77] observed the evolution of secondary waves generated via triadic resonant instability ^[77]. Figure 5(b) shows the net wave field containing both the primary and the secondary waves. In contrast to the inviscid prediction that the most unstable secondary waves occur at half the primary wave frequency and infinitely large wavenumber, the experiments revealed finite-wavenumber secondary waves at frequencies noticeably far from half the primary wave frequency. The observed difference was attributed to viscous effects. Indeed, earlier studies ^[78] have highlighted that viscous effects would completely suppress the inifinite wavenumber limit of branches A & B in Figure 2(b), and secondary waves whose wavelengths are comparable to the primary wavelength would emerge. Bourget et al. (2013) ^[77] also highlighted the effect of the finite spatial extent of the primary wave on the triadic resonant instability, which we discuss further in Section 4.2.

4.1.4.   Interaction with boundary

The interaction between two plane waves in the primary wave field could potentially excite a superharmonic wave field whose wave vector and frequency are the sum of the corresponding values of the primary waves. As discussed in Section 2.4, the right hand sides of Eqs (2.6) and (2.5) are zero for a plane internal wave in a uniform stratification. As a result, a self-interacting plane internal wave does not excite a superharmonic wave field. The presence of two non-collinear plane wave vectors in the primary wave field, say ${\bf k_1}$ and ${\bf k_2}$ (with $|{\bf k_1}|\ne |{\bf k_2}|$) such as what is shown in Figure 2, could, however, lead to resonant superharmonic generation. The interaction of a plane wave with a boundary, such as bottom topography or free surface, typically creates a primary wave field that contains non-collinear plane wave vectors with the same frequency. Such a scenario then enables the generation of a superharmonic wave field.

Figure 6 shows the schematic of a plane wave incident on a sloping boundary. Conservation of frequency and the along-slope wavenumber upon reflection results in the reflected plane wave being of a shorter wavelength than the incident wave. The interaction between the incident and reflected plane waves could hence result in superharmonic wave generation ^[79], which under certain conditions (depending on primary wave and bottom slope orientations -$\beta$ and $\gamma$, respectively in Figure 6) will represent triadic resonance between the primary incident, reflected plane waves and the superharmonic wave ^[80]. A recent theoretical study ^[81] has shown that the superharmonic wave amplitude increases linearly from the slope at resonance. It is also noteworthy that in the limit of the sloping boundary being horizontal, no superharmonic wave field is generated since the wave vector magnitude of the incident and reflected plane waves are equal ^[33,81].

Dauxois & Young (1999) ^[82] theoretically investigated the near critical incidence ($\beta\approx\gamma$ in figure 6) of an internal wave on a slope, to obtain the initial evolution of the superharmonic wave field and conditions under which overturning of the buoyancy field occurs, resulting in rapid transition to turbulence and enhanced mixing near the slope. A more recent study has shown that the exact solution of the weakly nonlinear and weakly viscous equations at critical incidence is well approximated by a sum of the incident wave, a reflected second harmonic and some boundary layer terms ^[83]. The theoretical predictions of Dauxois & Young (1999) ^[82] for the reflected wave field at the critical incident wave frequency were recently validated in laboratory experiments ^[84]. Laboratory experiments of a quasi-plane wave incident on a slope ^[85] have also shown the excitation of superharmonics and trapped evanescent higher harmonics near the slope. Other experiments, however, have shown the generation of dominant wave induced mean flows near the slope, which suppress the excitation of higher harmonics ^[86]. The reflection of finite-width wave beam on a slope is also an important problem in realistic ocean settings and is discussed in section 4.2.

4.2.   Internal wave beams

As discussed in Section 2.3.2, wave beams of a finite spatial width form when a distribution of unidirectional plane waves is superimposed. The characteristics of such a distribution, like its dominant wavenumber and width, depend on the nature of the forcing. As already pointed out, such wave beams are also solutions to the fully nonlinear inviscid equations of motion. While the study by Bourget et al. (2013) ^[77], depicted in Figure 5, focused on PSI in quasi-plane waves, the effects of finite spatial extent of the primary wave were already evident in their experiments. Specifically, though branch C (Figure 2b) has growth rates comparable to branch B in the presence of viscous effects at laboratory scales, the growth of secondary waves whose wavelength is larger than the primary wave is impeded by the finite spatial extent of the primary wave beam ^[77]. It is useful to recall that the subharmonic secondary waves, owing to their smaller frequency and hence shallower propagation angles, can leave the primary wave beam, as depicted in Figure 5(b). Finite-width effects on PSI in an internal wave beam, using energy balance arguments ^[87] and asymptotic analysis ^[88], have now been addressed in detail, with the growth of the secondary waves being shown to strongly depend on how rapidly they leave the finite spatial extent of the primary wave beam. In summary, a narrow wave beam has a smaller interaction zone, and hence a smaller growth rate, for the secondary waves, in comparison to a plane wave of an infinite spatial extent. Furthermore, laboratory experiments and numerical simulations for large-amplitude internal wave beams have revealed how PSI can lead to wave breaking ^[89]. A detailed discussion of triadic resonant instability in internal wave beams can be found in ^[32].

Internal wave beams have been subject to linear stability analyses too. Performing an asymptotic analysis of the evolution of small-amplitude long-wavelength modulations on an internal wave beam, Kataoka & Akylas (2013) ^[90] found a three-dimensional instability resulting from a resonant interaction between the primary wave beam, three-dimensional perturbations at the beam frequency and a mean flow, which can be significant even far from the vicinity of the primary wave beam. This instability occurs above a threshold steepness of the primary wave beam ^[90]. Considering small-amplitude two-dimensional perturbations (with no assumptions on their spatial scale) on a locally confined internal wave beam, Fan & Akylas (2020) ^[91] performed a Floquet analysis of the linearized perturbation equations to discover an instability beyond a threshold wave beam amplitude. This instability involving two subharmonic perturbations with wavepacket-like spatial structure was observed in experimental realizations of wave beams generated by finite-amplitude oscillations of a cylinder, and cannot be recovered in a classical triadic resonance calculation which assumes a small amplitude for the primary wave beam as well ^[91]. Focusing on the limit of fine-scale perturbations, Fan & Akylas (2021) ^[92] used their earlier approach ^[91] to show that PSI in confined internal wave beams involves secondary waves at 3/2 times the primary wave frequency in addition to those at half the primary wave frequency. Furthermore, Fan & Akylas (2021) ^[92] have shown that small-amplitude internal wave beams that are not susceptible to PSI may develop an instability with a broadband frequency spectrum.

Interacting internal wave beams can allow for the constituent non-collinear plane waves to interact, which as discussed in Section 4.1.4, can lead to superharmonic generation. Indeed, numerical studies have investigated higher harmonic generation, and the resulting wave breaking due to trapped evanescent waves, due to wave-wave interactions in colliding wave beams emanating from two different momentum sources ^[93]. Focusing on the scenario where the higher harmonic waves can propagate and hence radiate away from the interaction region, Tabaei et al. (2005) ^[94] used weakly nonlinear expansions to derive the steady-state wave field associated with the secondary beams whose frequencies are equal to the sum and difference of the primary wave beam frequencies. Jiang & Marcus (2009) ^[95] then proposed selection rules that determine which of the possible four secondary wave beams (that are allowed for a given higher harmonic frequency) would actually appear for a given pair of intersecting primary wave beams. An extension to the case of obliquely intersecting internal wave beams can be found in ^[96].

As for plane waves, reflection from an inclined slope paves way for interaction between the incident and reflected wave beams. Trapping of evanescent higher harmonic waves in the interaction region on a slope has been studied numerically by Javam et al. (1999) ^[97]. In the regime where higher harmonics can propagate, laboratory studies ^[98] have revealed the excitation of superharmonic wave beams when an internal wave beam is incident on a slope of a suitable inclination (Figure 7a). As for PSI in a wave beam, viscous effects are important in wave beam reflection off a slope at laboratory scales ^[81,99], and recent theoretical studies have developed weakly nonlinear models to incorporate viscous effects ^[100]. In addition to superharmonic generation, amplitude intensification due to reflection from a slope can also render the reflected wave beam more susceptible to PSI than the incident wave beam ^[101]. Interestingly, in contrast to a plane wave ^[33,81] (see Section 4.1.4), wave beam reflection from a horizontal boundary can generate higher harmonics too, as reported in both theory ^[94] and simulations ^[102,103]. On a horizontal boundary, the nature of the boundary condition can also influence the higher harmonic generation – for example, mean currents, superharmonic and subharmonic secondary wave generation have all been observed in numerical simulations of wave beam reflection from a free-slip surface ^[104]. Finally, colliding internal wave beams can also be a feature in internal wave generation by bottom topography, as depicted in Figure 7(b) based on the numerical simulation of Lamb (2004) ^[105].

4.3.   Vertical modes

In this section, we discuss both subharmonic and superharmonic resonances that occur in internal waves modes, whose linear regime was discussed in Section 2.3.3. Internal wave mode interactions, and associated amplitude evolution equations at triadic resonance, were considered by Thorpe (1966) ^[106], though the focus was on resonant triads containing internal waves and surface gravity waves. In general, the weakly nonlinear wave field associated with an interaction between modes $m$ and $n$ at frequencies $\omega_1$ and $\omega_2$, respectively, in a uniform stratification contains two different terms: one at frequency $\omega_1 + \omega_2$ with a vertical structure similar to that of mode $|m-n|$, and the other at frequency $\omega_1-\omega_2$ with a vertical structure similar to that of mode $m+n$ ^{[34,106,107,108]}. Such an interaction would then be resonant if the horizontal spatial structure of either of the aforementioned terms corresponds to an internal wave mode. It is worth highlighting that a direct equivalent of the triadic resonance conditions in Section 3 is easily identified only for the frequencies and the horizontal wavenumbers of internal wave modes since they have no propagation in the vertical direction.

With respect to subharmonic resonance, Davis & Acrivos (1967) ^[109] performed a combined theoretical and experimental study to investigate subharmonic secondary wave generation from a primary mode $1$ that propagates within a continuously stratified thin layer between two constant density layers. In the presence of viscous effects, they identified a threshold amplitude above which the primary mode $1$ succumbs to a subharmonic instability, with the subharmonic secondary modes being at frequencies that are noticeably different from half the primary wave frequency. For a given primary mode in a uniform stratification, inviscid theory predicts maximum growth rate for secondary subharmonic waves at half the primary wave frequency and infinitely large mode numbers whose difference equals the primary wave mode number ^[34]. Viscous effects, however, prevent this PSI limit to occur under laboratory conditions. Martin et al. (1972) ^[34], using a paddle-type wave generator that spans the entire depth of a uniformly stratified fluid, observed different subharmonic resonant triads excited simultaneously when a mode $1$ or a mode $3$ primary wave was forced. The subharmonic secondary wave frequencies were again different from half the primary wave frequency, and the corresponding mode numbers were consistent with the most unstable set (which is similar to the branch B identified for plane waves in Figure 2b) of triadic resonances identified by the amplitude evolution equations. More recently, Joubaud et al. (2012) ^[110] have used an internal wave generator ^[38] to generate a primary mode $1$ wave and quantitatively measure the growth rate associated with the dominant subharmonic resonant triad that is observed (Figure 8). The subharmonic frequencies were shifted away from half the primary wave frequency due to viscous effects, and the measured growth rates were then found to be reasonably consistent with the growth rates associated with suitably chosen plane wave resonant triads. Interestingly, subharmonic resonance seems to be readily observed only for large primary mode $1$ frequencies ($\omega_0/N \gtrsim 0.9$) in laboratory settings ^[110], though large primary wave amplitudes have not been studied in detail.

Subharmonic resonances have been studied in horizontally finite domains too. The degeneration of a continuously forced standing internal wave mode via a subharmonic resonant triad was investigated theoretically and experimentally in the presence of viscous damping by McEwan (1971) ^[111]. Specifically, he identified a critical amplitude for a primary mode $1$ (forced using a paddle type generator at one end of the tank) above which it irreversibly loses energy to a pair of subharmonic free wave modes belonging to the most unstable resonant triad. Subsequently, McEwan et al. (1972) ^[112] allowed for multiple subharmonic resonant triads that contain the primary mode to evolve simultaneously, and concluded that the steady state is often governed by the triad with the lowest critical amplitude for the primary wave. Detailed two-dimensional numerical simulations in the weakly viscous regime have then revealed that PSI, instead of a cascade of wave interactions, drives a standing mode towards wave breaking and dissipation ^[113]. Benielli & Sommeria (1998) ^[114] excited standing modes by vertical oscillations of a horizontally and vertically confined tank, and again observed PSI to drive wave breaking and turbulence. However, they also highlighted the importance of boundary layer and three-dimensional instabilities in the overall evolution, and also reported intermittency where the primary wave alternately grows and decays over long times.

With respect to superharmonic resonance, it is instructive to recall that a left-to-right propagating mode $m$ of frequency $\omega$ in a uniform stratification can be constructed as the sum of an upward propagating plane wave with $(k_x, k_z) = (|\bf{k}|\sin\theta, |\bf{k}|\cos\theta)$ and a downward propagating plane wave with $(k_x, k_z) = (|\bf{k}|\sin\theta, -|\bf{k}|\cos\theta)$, where $\theta$ is as defined in Eq (2.8). As mentioned for plane wave reflection on a horizontal boundary (Section 4.1.4), an interaction between such a pair of plane waves with the same wave vector magnitude does not result in superharmonic generation. A combination of modes at the same frequency, however, can result in superharmonic wave generation due to the interaction between their constituent plane waves. Thorpe (1966) ^[106] identified the conditions under which modes $m$ and $n$ at frequency $\omega$ are in triadic resonance with mode $|m-n|$ at frequency $2\omega$. Interestingly, in a uniform stratification, the temporal and horizontal spatial resonance conditions in plane waves (Section 3) are necessary but not sufficient for triadic resonance in internal wave modes ^[108]. Amplitude evolution equations for resonant superharmonic generation by modal interactions have been validated for early times in numerical simulations ^[115] and laboratory experiments (^[116], see Figure 9). Relative importance of superharmonic resonance with respect to subharmonic resonance, viscous effects and off-resonant superharmonic generation are some aspects from the aforementioned studies that would be worthwhile to investigate further in detail. Internal wave modes $m$ and $n$ at different frequencies have also been found to resonantly generate an internal wave mode $m+n$ at the difference frequency. Finally, an interesting recent study by Liang et al. (2017) ^[117] has shown that nonlinearity in the free surface boundary condition can make even an isolated vertical mode unstable via superharmonic resonance.

5.   Ocean perspective

The ocean is more complex than what is represented in the idealized models discussed thus far in this review. In this section, we briefly overview the important factors and questions regarding internal waves in the ocean, focusing on their instabilities.

One significant factor influencing internal wave dynamics in the ocean is the nonuniformity in the vertical stratification profile. Specifically, the ocean typically contains a relatively thin unstratified mixed layer near its surface, where turbulence driven by winds, heat and buoyancy fluxes homogenizes the density in the top few to several tens of metres. A strongly stratified pycnocline then lies between the mixed layer and the more or less uniformly stratified deep waters (see Figure 10(a) for a representative stratification profile in the ocean). In the deeper ocean below the pycnocline, even if the stratification is not uniform, it often varies sufficiently slowly with depth so that a Wentzel–Kramers–Brillouin–Jeffreys (WKBJ) approximation, i.e., a ray tracing approach, allows for reasonably accurate modeling of internal wave propagation ^[1,118]. In the upper ocean, however, the stratification is not necessarily slowly varying, resulting in non-trivial internal wave propagation characteristics even in the linear regime. Several studies have pointed out selective transmission and reflection of internal wave energy across nonuniformly stratified regions ^{[6,7,119,120,121,122,123,124,125]}. Owing to the continuous transmission and reflection, nonuniform stratifications enable interactions between upward and downward propagating internal waves. A classical example in the ocean is that of an internal wave beam incident on the pycnocline, which results in subharmonic ^{[126,127,128]} and superharmonic ^{[129,130,131,132]} resonance depending on parameters such as incident wave beam profile and amplitude, and the pycnocline characteristics.

For a finite-depth ocean, where the rigid lid approximation is mostly valid ^[1], the presence of a pycnocline modifies the vertical mode shapes $\phi_n(z)$ from those for a uniform stratification (Section 2.3.3). Specifically, the mode shapes in a nonuniformly stratified ocean are governed by the Sturm-Liouville equation ^[1]

where $N(z)$ is the stratification profile, and the no-normal-flow boundary conditions require $\phi_n(z = 0) = \phi_n(z = H) = 0$. Individual modes in a nonuniform stratification don't have a unique vertical wavelength, resulting in a non-trivial vertical structure for the weakly nonlinear wave field associated with modal interactions. Triadic resonance in finite-depth nonuniform stratifications has received attention in recent years, with new resonant triads being enabled in stratifications that contain a pycnocline ^[108]; the occurrence of these newer resonances can be understood as a consequence of the alternative theorem for linear differential equations of the second order ^[133]. Interestingly, even individual modes in nonuniform stratifications can be unstable due to resonant self-interaction ^[108,134]. Theoretical studies on amplitude evolution at and off resonance, and numerical simulations have indeed confirmed superharmonic generation by self-interacting modes in nonuniform stratifications ^[52,115,135].

In addition to nonuniform stratification profiles, another factor of dynamical importance in the ocean is the Coriolis effect due to the Earth's background rotation. Fundamentally, the internal wave dispersion relation introduced in Section 2.2 gets modified by background rotation as ^[1]

where $f = 2\Omega_{e}\sin\zeta$ ${\hat{{\bf{e}}}_{\bf{z}}}$, with $\Omega_e$ and $\zeta$ being the Earth's angular velocity of rotation and latitude, respectively. As a result, the two inertia-gravity wave propagation regimes are (ⅰ) $f\leq \omega \leq N_0$ and (ⅱ) $N_0 \leq \omega \leq f$. An interesting limit for PSI is the near-inertial limit of $\omega/2 \approx f$, where the subharmonic waves are at the edge of the aforementioned first wave propagation regime. Indeed, northward propagating internal tides in the ocean can undergo rapid deterioration of their energy at latitudes where $f$ is half the tidal frequency ^[136,137]. Theoretical studies have investigated energy transfer rates to fine-scale subharmonic waves due to PSI in the near-inertial limit in plane waves, modes and internal wave beams ^{[138,139,140,141]}. Towards understanding the lack of strong evidence of near-inertial PSI in field observations, a recent linear stability analysis based on Floquet theory has highlighted non-trivial latitudinal dependence of instabilities in finite-amplitude internal wave beams ^[142]. Specifically, Onuki & Tanaka (2019) ^[142] have shown that the location of maximum instability in finite-amplitude inertia-gravity wave beams can occur at values of $f/\omega$ smaller than 0.5, while the largest PSI growth rate occurs at $f/\omega = 0.5$ for small-amplitude inertia-gravity wave beams.

Maurer et al. (2016) ^[143] have extended the studies of Bourget et al. (2013) ^[77] and Bourget et al. (2014) ^[87] to incorporate the effects of background rotation on TRI in plane inertia-gravity waves and finite-width inertia-gravity wave beams. They report that TRI is enhanced in a specific range of $f$, and their experiments also report the generation of sub-inertial secondary waves, which require further investigations. Investigating triadic resonant instability in the presence of both background rotation and stratification in a horizontally periodic domain, Sutherland & Jefferson (2020) ^[144] observed subharmonic secondary waves in numerical simulations though theory predicts that secondary waves with frequencies larger than the primary wave constitute the most unstable resonant triad. In regards to superharmonic resonance, Varma & Mathur (2017) ^[108] have considered steady-state weakly nonlinear modal interactions in the presence of both background rotation and stratification, and highlighted the strong possibility of triadic resonance resulting from high mode interactions at $\omega\approx f$. Such superharmonic resonances are potentially relevant for near-inertial waves excited in the upper ocean by winds ^[145]. With respect to the second wave propagation regime, it is often studied in the limit of $N_0 = 0$, i.e., inertial waves for $0\le\omega\le f$. Though inertial and internal waves are similar in terms of their dispersion relation, their instability characteristics can be different. Subharmonic resonance in inertial waves have been observed in experiments ^[146], and the three-dimensionality of triadic resonant instability in inertial waves has recently been highlighted ^[147]. Interestingly, a recent study has shown that transient growth and relaxation oscillations can be important in linearly stable resonant triads of inertial waves ^[148].

A few important considerations drive studies on internal waves in the ocean. The first concerns the spectrum associated with fluid motions in the ocean. For example, Figure 10(b) shows the frequency spectrum of kinetic energy measured at a single location near the ocean floor at a site of steep topography ^[149]. At high frequencies, strong peaks are observed at the local Coriolis frequency and the semi-diurnal frequency (indicated a $M_2$ in Figure 10(b)) associated with internal tides. The energy content at the Coriolis frequency could be a consequence of either wind-forced near-inertial motions (generated in the upper ocean) or subharmonic instabilities in the $M_2$ internal tides. Interestingly, higher harmonics, including 2$M_2$, are also observed in the spectrum. The low frequency part of the spectrum is dominated by eddy motions, with the relative fraction of energy in the eddy and wave motions being a function of depth as well. It is worth noting that the slow eddy motions and the fast internal wave motions can exchange energy too ^{[150,151,152]}. Moving away from strong topographic sites, the internal wave frequency band between $f$ and $N$ is usually described by a heuristic model that has come to be known as the Garrett-Munk (GM) spectrum ^[19,153]. The GM spectrum is expected to be a result of nonlinear interactions between various frequencies observed in the spectra such as the one shown in Figure 10(b), with the additional complexity of interaction between various spatial scales that are sometimes introduced at the generation stage itself ^[12,154]. Indeed, energy cascades and wave turbulence resulting from various interactions in internal waves are active areas of current research ^{[50,155,156,157,158,159]}.

6.   Conclusions & future directions

In this paper, we have reviewed the existing knowledge on instabilities in internal waves of various spatial forms. At sufficiently small amplitudes, internal wave instabilities occur in the form of different manifestations of triadic resonance. While short-term evolution of isolated triadic resonances is reasonably well-understood in the inviscid limit, several outstanding questions remain even for seemingly simple settings. Examples of such unresolved questions include (ⅰ) which resonant triads would drive instabilities in viscous finite-amplitude internal waves? – it is worth recalling that classical theoretical models of TRI are often in the weakly viscous, small-amplitude regime. (ⅱ) what is the long-term behavior associated with triadic resonances, which is likely related to the dynamics associated with multiple simultaneously occurring resonant triads ^[34]. (ⅲ) what are the instability dynamics in systems operating away from resonance ^{[135,160,161]}? In addition, the occurrence of a mean flow in three-dimensional viscous internal waves has been recognized in several recent studies ^[32], and the effects of mean flows on triadic resonances are currently being investigated ^{[162,163,164]}. An alternative, and possibly more appropriate for the ocean, approach to study energy transfer rates in a continuous internal wave spectrum through wave-wave interactions is the kinetic equation ^[165,166], which has recently been used to investigate the decay of internal tides near the critical latitude ^[167].

Apart from the various considerations discussed in Section 5 for the ocean, other factors could be important too. Three-dimensionality, associated with both the internal wave forcing mechanisms (3D topography, for example) and the instabilities ^[30,59,168], is one such factor. In the lab, axisymmetric configurations ^{[169,170,171,172]} provide useful insight as an intermediary between two-dimensional and fully three-dimensional configurations. Internal wave interaction with background flows, such as vertical shear ^{[173,174,175]} and upper ocean balanced flows ^{[150,152,176]}, is also a relevant factor in the ocean. Another potentially relevant topic of interest is the non-traditional effect, i.e., the role of the component of Earth's background rotation in the tangential plane of Earth's surface ^[177].

While instabilities represent an important first step towards small-scale turbulence and dissipation in internal waves, understanding and quantifying internal wave driven mixing remains a challenge ^{[18,149,178,179]}. Finally, an accurate parameterization of internal wave driven mixing in large scale ocean models is an essential ingredient of climate modeling ^[20,180]. We conclude this review by stating that internal wave dynamics, including instabilities, are an important consideration in atmospheric ^[181] and astrophysical ^[182] applications too.

Acknowledgments

The authors acknowledge support from the Simons Foundation through Grant No. 651475 (DV and TD), the Indian Ministry of Human Resource Development under the SPARC programme (Sanction Order No. SPARC/2018-2019/P1213/SL), the research initiative Geophysical Flows Lab at IIT Madras and Grant ANR-17-CE30-0003 (DisET). The authors thank Sylvain Joubaud and Yohei Onuki for their helpful comments on the manuscript.

Conflict of interest

The authors declare no conflict of interest.