Hostname: page-component-797576ffbb-k7d4m Total loading time: 0 Render date: 2023-12-03T21:48:46.901Z Has data issue: false Feature Flags: { "corePageComponentGetUserInfoFromSharedSession": true, "coreDisableEcommerce": false, "useRatesEcommerce": true } hasContentIssue false

Resonant and near-resonant internal wave triads for non-uniform stratifications. Part 2. Vertically bounded domain with mild-slope bathymetry

Published online by Cambridge University Press:  20 June 2022

Saranraj Gururaj*
School of Science and Engineering, University of Dundee, Dundee DD1 4HN, UK
Anirban Guha
School of Science and Engineering, University of Dundee, Dundee DD1 4HN, UK
Email address for correspondence:


Weakly nonlinear internal wave–wave interaction is a key mechanism that cascades energy from large to small scales, leading to ocean turbulence and mixing. Oceans typically have a non-uniform density stratification profile; moreover, submarine topography leads to a spatially varying bathymetry ($h$). Under these conditions and assuming mild-slope bathymetry, we employ multiple-scale analysis to derive the wave amplitude equations for weakly nonlinear wave–wave interactions. The waves are assumed to have a slowly (rapidly) varying amplitude (phase) in space and time. For uniform stratifications, the horizontal wavenumber ($k$) condition for waves (1, 2, 3), given by ${k}_{(1,a)}+{k}_{(2,b)}\!+\!{k}_{(3,c)}\!=\!0$, is unaffected as $h$ is varied, where $(a,b,c)$ denote the mode number. Moreover, the nonlinear coupling coefficients (NLC) are proportional to $1/h^2$, implying that triadic waves grow faster while travelling up a seamount. For non-uniform stratifications, triads that do not satisfy the condition $a=b=c$ may not satisfy the horizontal wavenumber condition as $h$ is varied, and unlike uniform stratification, the NLC may not decrease (increase) monotonically with increasing (decreasing) $h$. NLC, and hence wave growth rates for weakly nonlinear wave–wave interactions, can also vary rapidly with $h$. The most unstable daughter wave combination of a triad with a mode-1 parent wave can also change for relatively small changes in $h$. We also investigate higher-order self-interactions in the presence of a monochromatic, small-amplitude topography; here, the topography behaves as a zero-frequency wave. We derive the amplitude evolution equations and show that higher-order self-interactions might be a viable mechanism of energy cascade.

JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (, which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
© The Author(s), 2022. Published by Cambridge University Press

1. Introduction

Low-mode long-wavelength internal gravity waves in oceans can travel thousands of kilometres from their generation site without dissipation (Zhao et al. Reference Zhao, Alford, Girton, Rainville and Simmons2016). The energy in these long waves can cascade to small scales through a variety of mechanisms, such as nonlinear interactions among the waves (MacKinnon & Winters Reference MacKinnon and Winters2005; MacKinnon et al. Reference MacKinnon, Alford, Sun, Pinkel, Zhao and Klymak2013), scattering through interaction with the seafloor topography (Legg & Adcroft Reference Legg and Adcroft2003), and scattering through interaction with turbulent quasi-geostrophic flows (Kafiabad, Savva & Vanneste Reference Kafiabad, Savva and Vanneste2019). This transfer of energy to small scales will lead eventually to turbulence and mixing, which is essential for maintaining the meridional overturning circulation (Munk Reference Munk1966). In weakly nonlinear wave–wave interactions, an internal gravity wave can become unstable via resonant triad interactions if it has the largest frequency in the triad (Hasselmann Reference Hasselmann1967); through this mechanism, energy is transferred irreversibly from a high-frequency and low-wavenumber primary wave to lower-frequency and higher-wavenumber secondary waves. In a resonant internal wave triad, a wave of angular frequency $\omega _3$ and wave vector $\boldsymbol {k}_{3}$ can transfer its energy resonantly to two ‘daughter’ waves when both the conditions $\boldsymbol {k}_{3} = \boldsymbol {k}_{1} + \boldsymbol {k}_{2}$ and $\omega _3 = \omega _1 + \omega _2$ are met (Davis & Acrivos Reference Davis and Acrivos1967; Hasselmann Reference Hasselmann1967; Phillips Reference Phillips1967). The former condition is a consequence of the quadratic nonlinearity of the Navier–Stokes equations.

Since ocean's density stratification is non-uniform, recent efforts have been directed towards understanding energy transfer in non-uniformly stratified fluids. In Varma & Mathur (Reference Varma and Mathur2017), the conditions for the existence of resonant, weakly nonlinear wave–wave interactions in a non-uniform stratification were studied. They proved that resonant triads and self-interactions can exist if (i) they satisfy the horizontal wavenumber condition, and (ii) each wave's functional form in the $z$-direction is non-orthogonal to the nonlinear forcing terms. Wunsch (Reference Wunsch2017) studied self-interaction of an internal wave mode in the presence of a non-uniform stratification; the latter was simplified using a three-layer model. It was shown that the amplitude of the superharmonic wave, which is forced by the self-interaction of a parent wave, can be highly sensitive to changes in the stratification profile characteristics such as pycnocline depth and strength. Moreover, Liang, Zareei & Alam (Reference Liang, Zareei and Alam2017) showed that self-interaction also occurs in the presence of uniform stratification, provided that the nonlinear terms in the free surface boundary condition are taken into account. Self-interaction was also studied numerically by Sutherland (Reference Sutherland2016), and it was observed that in the presence of non-uniform stratification, self-interaction of an internal wave mode was more dominant than triadic interactions for low Coriolis frequency. Furthermore, Baker & Sutherland (Reference Baker and Sutherland2020) studied self-interaction of a mode under angular frequency mismatch, and found that the daughter wave (superharmonic wave) can return its energy to the parent wave.

Apart from the weakly nonlinear wave–wave interactions, wave–topography interactions (where the topography is of small amplitude) have also been studied extensively. In Buhler & Holmes-Cerfon (Reference Buhler and Holmes-Cerfon2011), the decay of a mode-$1$ internal tide due to its interaction with a small-amplitude sea floor topography was studied using ray-tracing. It was shown that if the bottom bathymetry is ‘resonant’ (see § 6 for more detail), then the internal mode-$1$ wave interacts with the bathymetry and resonantly gives its energy to the higher modes. The topography in this case acts like a stationary wave with zero angular frequency. In Couston, Liang & Alam (Reference Couston, Liang and Alam2017), this scattering process was explored in a three-dimensional setting, where the mode-$1$ internal wave was incident obliquely on a small-amplitude bottom topography. Buhler & Holmes-Cerfon (Reference Buhler and Holmes-Cerfon2011) and Li & Mei (Reference Li and Mei2014) also focused on scattering of an internal gravity wave by a small-amplitude, stationary, zero-mean random topography. Li & Mei (Reference Li and Mei2014) consider topographies that vary in zonal and meridional directions. Both studies, under realistic parameters, estimate a decay length scale of about 500–1000 km for the mode-1 wave.

In Mathur, Carter & Peacock (Reference Mathur, Carter and Peacock2014), a Green's function approach along with numerical simulations was used to study internal gravity wave scattering under the linear inviscid limit in the presence of constant and non-constant buoyancy frequency in a two-dimensional setting for large-amplitude topographies. Height of the topography and criticality ($C_r$) were the two main factors that influence internal gravity wave scattering. In general, subcritical ($C_r < 1$) topographies were found to scatter the incoming wave lesser than supercritical topographies ($C_r > 1$). Critical topographies ($C_r \approx 1$) were the most proficient in scattering the incoming wave. Scattering of large-amplitude waves – with breaking and the ensuing kinetic energy dissipation – is a very important quantity to study since it provides an estimate for local diffusivity. Internal wave breaking due to different types of topographies was focused on in Legg (Reference Legg2014). In particular, a condition for internal wave breaking was given using the incoming wave's Froude number ($Fr$). The Froude number for a mode-1 wave is defined as:

(1.1)\begin{equation} Fr = \frac{U{\rm \pi}\alpha}{H\omega}, \end{equation}

where $U$ and $\omega$ are respectively the peak horizontal velocity and frequency of the wave, $H$ is the depth of the domain, and $\alpha$ is the slope of the wave. As the mode-1 wave shoals up a large-amplitude topography, its $Fr$ increases. It is empirically determined that if a wave's $Fr$ reaches a range 0.3–1 due to shoaling, then the wave is prone to breaking. A wave's local Froude number can also increase due to reflection from a topography. Highly nonlinear features such as bores were observed in regions of the topography where the local Froude number was greater than $1$ (Legg & Adcroft Reference Legg and Adcroft2003). Scattering and dissipation due to large-amplitude highly supercritical topographies were focused on in Klymak et al. (Reference Klymak, Buijsman, Legg and Pinkel2013). Interestingly, it was observed that a mode-1 wave loses a maximum ${\sim }20\,\%$ of its energy at an isolated tall supercritical topography.

This paper is ‘Part 2’ of Gururaj & Guha (Reference Gururaj and Guha2020), in which the effect of non-uniform (albeit slowly varying) stratification on internal wave triads in an unbounded domain was studied theoretically and numerically. It was shown that the variation in stratification profile may affect significantly the nonlinear coupling coefficients, (vertical wavenumber) detuning, and group speed of the wave packets constituting a triad, and hence the ensuing energy transfer. Different triads were also observed to undergo different amounts of detuning for the same change in the background stratification. The present paper extends the paradigm explored in Gururaj & Guha (Reference Gururaj and Guha2020) to a vertically bounded domain with a mild slope bathymetry, while there is no more restriction for the non-uniform stratification to be slowly varying. To the best of our knowledge, this is the first work that considers the effect of the variation of the ocean depth on internal wave triads. A simplified schematic of the set-up is given in figure 1. The motivation behind this study stems from the simple fact that ocean depth varies spatially (hence consideration of constant depth might be an over-simplification), hence waves can move from one depth to another while they are interacting in a medium of varying stratification. While the effect of change in fluid depth on resonant and near-resonant interactions between three distinct waves is the primary focus of this paper, we have also studied higher-order self-interactions among the internal waves in the presence of a small-amplitude topography. As an analogy, such higher-order interactions have been studied in Alam, Liu & Yue (Reference Alam, Liu and Yue2009) for surface waves. Wave–wave interactions in non-uniformly stratified vertically bounded domains have been considered previously in various studies (Young, Tsang & Balmforth Reference Young, Tsang and Balmforth2008; Baker & Sutherland Reference Baker and Sutherland2020; Varma, Chalamalla & Mathur Reference Varma, Chalamalla and Mathur2020). In this paper, reduced-order models for wave–wave interactions in a region of varying $h$ are derived, hence spatially varying nonlinear coupling coefficients, group speed and detuning are all involved. Moreover, the equations can model wave–wave interactions of wave trains or finite width wave packets in a region of varying $h$. In § 6, we have also derived (and validated numerically) reduced-order equations that model higher-order self-interactions in the presence of a small-amplitude topography. Equations in § 6 can also be used to model standard resonant self-interactions in the presence of slowly varying large-amplitude topographies.

Figure 1. (a) A general schematic of the problem to be studied, showing the streamfunction field of three wave packets interacting in the presence of a varying bathymetry $h(x)$. Here, $H$ is the mean depth (equivalent to the depth in flat bathymetry situation), while $H_b(x)$ denotes the submarine topography shape. (b) The stratification profile used in constructing the modes. The same non-uniform stratification model is used throughout the paper.

The paper is organized as follows. In § 2, we derive the amplitude evolution equations of the constituent waves of a triad in the presence of a slowly varying bathymetry using the Boussinesq Navier–Stokes equations in the $f$-plane. To derive these equations, the streamfunction, buoyancy perturbation and meridional velocity due to each wave are assumed to be a product of a slowly varying amplitude and a rapidly varying phase that are functions of space and time. In §§ 3 and 4, the effects on the horizontal wavenumber condition when waves interact in a region of varying ocean depth in the presence of uniform and non-uniform stratification are studied, respectively. In § 5, we study the effect of ocean depth variation on the rate of energy transfer in triadic interactions and self-interactions in the presence of a non-uniform stratification. In § 6, we analyse higher-order self-interactions of a wave in the presence of small-amplitude monochromatic topography. In § 7, the reduced-order equations derived in this paper are validated by solving the full Boussinesq equations using an open source code Dedalus (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020). The paper has been summarized in § 8.

2. Derivation of the governing equations in terrain-following coordinates

The incompressible, inviscid, two-dimensional (2-D) (in the $x$$z$ plane) Navier–Stokes equations on the $f$-plane under the Boussinesq approximation – hereafter referred to as the Boussinesq equations – can be expressed in terms of the perturbation streamfunction $\psi$, meridional velocity $v$ (along the $y$-direction), and perturbation buoyancy $b$ as follows:

(2.1a)$$\begin{gather} \frac{\partial}{\partial t}(\nabla^{2}\psi) + \frac{\partial b}{\partial x} - f\, \frac{\partial v}{\partial z} ={-}\{\nabla^{2}\psi,\psi\}, \end{gather}$$
(2.1b)$$\begin{gather}\frac{\partial v}{\partial t} + f\,\frac{\partial \psi}{\partial z} ={-}\{v,\psi\}, \end{gather}$$
(2.1c)$$\begin{gather}\frac{\partial b}{\partial t} - N^{2}\,\frac{\partial \psi}{\partial x} ={-}\{b,\psi\}. \end{gather}$$

Here, $N^{2}(z) \equiv -(g/\rho ^{*})({\rm d} \bar {\rho }/{\rm d} z)$ is the squared buoyancy frequency, $\bar {\rho }$ is the base density profile, $\rho ^{*}$ is the reference density, and $g$ is the acceleration due to gravity (directed along $-z$). The perturbation buoyancy is defined as $b \equiv -g\rho /\rho ^{*}$, where $\rho$ is the perturbation density.

The operator $\{G_1,G_2\} \equiv (\partial G_1/\partial x)(\partial G_2/\partial z) - (\partial G_1/\partial z)(\partial G_2/\partial x)$ denotes the Poisson bracket, and $f$ is the Coriolis frequency. Viscous effects have been neglected owing to the fact that we consider waves with long wavelengths.

The fluid domain is bounded at the top ($z=0$) by a rigid lid (i.e. zero vertical velocity, leading to the boundary condition $\psi (x,0) = 0$). The bottom boundary at $z = h(x)$ satisfies the impenetrable boundary condition $\psi (x,h(x)) = 0$.

Instead of solving the fully nonlinear equations (2.1a)(2.1c) numerically, we combine (2.1a)–(2.1c) into a single equation and employ a multiple-scale analysis. To this end, we perform $\partial (2.1a)/ \partial t + f\,\partial (2.1b)/\partial z -\partial (2.1c)/\partial x$, which results in

(2.2)\begin{equation} \frac{\partial^{2}}{\partial t^{2}}(\nabla^{2}\psi) + N^{2}\, \frac{\partial^{2} \psi}{\partial x^{2}} + f^{2}\,\frac{\partial^{2} \psi}{\partial z^{2}} ={-}\frac{\partial}{\partial t} (\{\nabla^{2}\psi,\psi\}) + \frac{\partial }{\partial x}(\{b,\psi\}) - f\,\frac{\partial }{\partial z}(\{v,\psi\}) . \end{equation}

Following the approach of Maugé & Gerkema (Reference Maugé and Gerkema2008), we now change the governing equations to terrain-following coordinates, where a new variable ($\eta$) is defined as

(2.3)\begin{equation} \eta \equiv{-} \frac{z}{h(x)}. \end{equation}

According to the definition (2.3), the bottom boundary condition at $z=h(x)$ would now be enforced at $\eta = -1$, while the surface boundary condition at $z=0$ remains unaltered, except that it is now at $\eta =0$. The governing equations, which are in the $x$$z$ coordinates, need to be transformed into the $x$$\eta$ coordinates. The correspondences between the variables in the $x$$z$ and $x$$\eta$ coordinate systems are as follows:

(2.4ac)\begin{equation} \psi(x,z,t) \Rightarrow \varPsi(x,\eta,t),\quad b(x,z,t) \Rightarrow B(x,\eta,t),\quad v(x,z,t) \Rightarrow \mathcal{V}(x,\eta,t). \end{equation}

On transforming the differential operators from $x$$z$ coordinates to $x$$\eta$ coordinates and substituting the transformed variables in (2.2), we arrive at

(2.5)\begin{align} \left[\frac{\partial^{2} }{\partial t^{2}}({L}_{xx}+{L}_{\eta\eta}) + N^{2}({-}h(x)\,\eta){L}_{xx} + f^2{L}_{\eta\eta}\right]\varPsi &={-}\frac{\partial}{\partial t} \left[\mathcal{J}\{({L}_{xx}+{L}_{\eta\eta})\varPsi,\varPsi\}\right] \nonumber\\ &\quad + {L}_{x}(\mathcal{J}\{B,\varPsi\}) - f\,{L}_{\eta}(\mathcal{J}\{\mathcal{V},\varPsi\}), \end{align}

where the operators $L_{x}, L_{\eta }, L_{xx}, L_{\eta \eta }$ and $\mathcal {J}\{G_1,G_2\}$ have the following definitions:

(2.6a)$$\begin{gather} L_{x} \equiv \frac{\partial}{\partial x} + \frac{\partial \eta}{\partial x}\,\frac{\partial}{\partial \eta},\quad L_{\eta} \equiv{-}\frac{1}{h}\, \frac{\partial}{\partial \eta},\quad L_{\eta\eta} \equiv \frac{1}{h^{2}}\,\frac{\partial^{2}}{\partial \eta^{2}}, \end{gather}$$
(2.6b)$$\begin{gather}L_{xx} \equiv \frac{\partial^{2}}{\partial x^{2}} + \frac{\eta^{2}}{h^{2}}\left(\frac{\partial h}{\partial x}\right)^{2}\frac{\partial^{2}}{\partial \eta^{2}} - 2\,\frac{\eta}{h}\,\left(\frac{\partial h}{\partial x}\right)\frac{\partial^{2}}{\partial \eta\,\partial x} + \frac{\eta}{h}\left[\frac{2}{h}\left(\frac{\partial h}{\partial x}\right)^{2} - \frac{\partial^{2} h}{\partial x^{2}}\right] \frac{\partial}{\partial \eta}, \end{gather}$$
(2.6c)$$\begin{gather}\mathcal{J}\{G_1,G_2\} \equiv {L}_{x} (G_1)\,{L}_{\eta} (G_2) - {L}_{\eta} (G_1)\,{L}_{x} (G_2). \end{gather}$$

For performing multiple-scale analysis, we assume wave-like perturbations, and the streamfunction due to the $j$th wave ($j = 1,2,3$) is given according to the following ansatz:

(2.7)\begin{equation} \varPsi_{j} = a_{j}(\epsilon_{x} x,\epsilon_{t} t)\,\varXi_{j}(x,\eta,t) + \mathrm{c.c.}, \end{equation}

where ‘c.c.’ denotes the complex conjugate, $a_{j}$ is the slowly varying complex amplitude, and $\varXi _j(x,\eta,t)$ is the rapidly varying phase part of the $j$th wave. The small parameters $\epsilon _{t}$ and $\epsilon _{x}$ are respectively used to denote the weak variation of the amplitude function with time and the streamwise ($x$) direction. The amplitude is assumed to be an ${O}(\epsilon _{a})$ quantity, where $\epsilon _{a}$ is a small parameter. The bathymetry ($h$), which is simply the negative of the fluid depth, is assumed to be of the form

(2.8)\begin{equation} h ={-}H + \epsilon_h\,H_b(k_b x), \end{equation}

where $H$ represents the mean depth of the fluid domain, $H_b$ denotes the submarine topography shape, $\epsilon _h$ is its amplitude, and $k_b^{-1}$ represents the length scale of the bathymetry. We always assume the bathymetry to have a ‘mild slope’; for this we use an analogue condition of that used for surface gravity waves (Meyer Reference Meyer1979; Kirby Reference Kirby1986):

(2.9)\begin{equation} {\frac{1}{\mathcal{K}_{j}}}\,\frac{\partial h}{\partial x} = {O}(\epsilon_h \epsilon_k) \ll {O}(1), \end{equation}

where $\mathcal {K}_j\equiv k_j h$ is the non-dimensional horizontal wavenumber ($k_j$ being the horizontal wavenumber) of the $j$th internal wave. Moreover, the relation $k_{j}^{-1} = \epsilon _k k_b^{-1}$ is used in (2.9), which implies that either of the parameters, $\epsilon _h$ or $\epsilon _k$, could be a small quantity, while the other could potentially be an ${O}(1)$ quantity. We note in passing that the mild slope condition in our case can still lead to internal gravity wave scattering. Internal wave scattering is largely dependent on the slope of the wave, which is almost constant; wave slope is dependent on $N$, which is nearly constant away from the pycnocline – even for higher modes whose horizontal wavenumber is much larger. Scaling analysis to find the relations between these small parameters is given in Appendix B.

2.1. Leading-order analysis

Next we substitute (2.7) in (2.5). At the leading order (${O}(\epsilon _{a})$), the governing equation (2.5) reduces to

(2.10)\begin{equation} \left(\frac{\partial^2}{\partial x^2} + \frac{1}{ h^{2}}\,\frac{\partial^{2}}{\partial \eta^{2}}\right) \frac{\partial^2 \varXi_j}{\partial t^{2}} + N^2(- h(x)\,\eta)\,\frac{\partial^2 \varXi_j}{\partial x^2} + \frac{f^2}{ h^{2}}\,\frac{\partial^2 \varXi_j}{\partial \eta^{2}} = 0. \end{equation}

Hereafter, we drop the argument of $N^2$, assuming that it is implied. Furthermore, assuming $\varXi _j = \widehat {\varXi }_j(x,\eta )\,{\rm e}^{-{\rm i}\omega _j t}$, where $\omega _j\in \mathbb {R}^+$ is the angular frequency of the $j$th internal wave, (2.10) simplifies to

(2.11)\begin{equation} \left[(N^2-\omega_j^{2})\,\frac{\partial^2}{\partial x^2} - \frac{\omega_j^{2}-f^2}{h^{2}}\,\frac{\partial^2}{\partial \eta^2}\right]\widehat{\varXi}_j = 0. \end{equation}

For a mild slope bathymetry (see Appendix B for details), we can use variable separation to solve (2.11) at the leading order. To this end, we assume $\widehat {\varXi }_j = {\phi }_j(\eta ;x)\,{P}_j(x)$, which leads to

(2.12)\begin{equation} \frac{h^{2}}{{P}_j}\,\frac{\partial^2 {P}_j}{\partial x^2} = \frac{\omega_j^{2}-f^2}{N^2-\omega^{2}_j}\,\frac{1}{{\phi}_j}\, \frac{\partial^2 {\phi}_j}{\partial \eta^2} ={-}\mathcal{K}_{j}^{2}, \end{equation}

where ${\phi }_j$ depends parametrically on $x$ via $h$. We emphasize that in the $x$$\eta$ coordinates, the presence of bathymetry makes $N$ also a function of $x$; see figure 2 for clarity.

Figure 2. The effective change in the stratification profile $N(z)$ when the coordinates are changed from (a) $x$$z$ to (b) $x$$\eta$, in the presence of bathymetry. For the latter case, if $N$ is a function of $z$ in $x$$z$, then it becomes a function of both $\eta$ and $x$ in $x$$\eta$. The $N$ profiles corresponding to the top of a seamount and an abyssal plain region have been denoted, respectively, by blue and green lines.

Two separate equations, one for ${P}_j$ and the other for ${\phi }_j$, can be formed from (2.12):

(2.13a)$$\begin{gather} \left[\frac{\partial^2}{\partial x^2} + \frac{\mathcal{K}_{j}^{2}}{h^{2}}\right]{P}_j = 0, \end{gather}$$
(2.13b)$$\begin{gather}\mathcal{L}_j {\phi}_j \equiv \left[\frac{\partial^2}{\partial \eta^2} + \mathcal{K}_{j}^{2} \chi_j^2 \right] {\phi}_j = 0, \end{gather}$$

where $\chi _j \equiv \sqrt {({N^2-\omega ^{2}_j)}/{(\omega _j^{2}-f^2})}$ is defined for convenience. The boundary conditions for (2.13b) are ${\phi }_j = 0$ at $\eta = 0,-1$. The non-dimensional horizontal wavenumber of the $j$th wave, i.e. $\mathcal {K}_j$, is the set of eigenvalues obtained from (2.13b), which can vary in $x$ when $N$ is a function of $x$ in $x$$\eta$ coordinates. An important point to note is the convention used in our study. While positive (negative) $k_j$ implies waves propagating along $+x$ ($-x$), owing to the fact that $h$ is negative, $\mathcal {K}_{j}$ follows the exact opposite convention. This means that a negative (positive) $\mathcal {K}_{j}$ implies that the wave is travelling along the $+x$ ($-x$) direction. Moreover we notice that (2.13b) does not depend explicitly on $h$. The only way (2.13b) can be influenced by $h$ is through $N$ when the latter varies in the $z$-direction (in $x$$z$ coordinates). However, for a uniform stratification, i.e. $N=\textrm {constant}$, eigenvalues of (2.13b) are independent of $h$. In this case, the eigenvalues are given by

(2.14)\begin{equation} \mathcal{M}_j \equiv \mathcal{K}_j \chi_j = n {\rm \pi}, \end{equation}

where $n\in \mathbb {Z}^+$. We also observe that the quantity $\mathcal {M}_j$ behaves like the vertical wavenumber of the wave that is non-dimensionalized by the local bathymetry $h$.

Meanwhile, ${P}_j$ at the leading order of the Wentzel–Kramers–Brillouin (WKB) approximation given by

(2.15)\begin{equation} {P}_j = \exp\left\{{\rm i} \int_0^{x} \frac{\mathcal{K}_j(x')}{h(x')}\,{{\rm d}\kern0.06em x}'\right\}. \end{equation}

We introduce a function $\beta _j (\epsilon _k x)$ such that ${P}_j$ is corrected to $P_j/\beta _j$. This slow varying function $\beta _j (\epsilon _k x)$, which acts as a correction to the first-order WKB solution (2.15), is given in (2.22). We note in passing that $P_j/\beta _j$ is still a solution of (2.13a) in the leading order even after the above-mentioned correction. To normalize the eigenfunction of the waves obtained from (2.13b), every wave's ${\phi }_j$ is constrained to satisfy

(2.16)\begin{equation} \frac{1}{2} \int_{{-}1}^{0} \frac{1}{h^2} \left[\mathcal{K}_j^2 {\phi}_j^2 + \left(\frac{\partial {\phi}_j}{\partial \eta} \right)^2 \right] \partial \eta = 1. \end{equation}

After this normalization, waves having the same amplitude ($a_j$) will also have the same energy density at a given $h$, provided that $\beta _j=1$.

The meridional velocity and the buoyancy perturbation at the leading order can be obtained by respectively converting (2.1b) and (2.1c) into the $x$$\eta$ coordinates and then substituting the streamfunction ansatz (2.7):

(2.17)$$\begin{gather} \mathcal{V}_{j} = {\rm i}\,\frac{ f}{h\omega_{j}}\,\frac{a_{j}}{\beta_j}\, \frac{\partial \phi_j}{\partial \eta}\,P_{j}\,{\rm e}^{-{\rm i}\omega_{j}t} + \mathrm{c.c.}, \end{gather}$$
(2.18)$$\begin{gather}B_{j} = {\rm i}\,\frac{N^{2}}{\omega_{j}}\,\frac{a_{j}}{\beta_j}\,\frac{\partial P_{j}}{\partial x}\,\phi_j\,{\rm e}^{-{\rm i}\omega_{j}t} + \mathrm{c.c}. \end{gather}$$

2.2. Second-order analysis

2.2.1. Amplitude evolution equations for a resonant triad in non-uniform stratification

Triad interaction between three internal waves occurs at ${O}(\epsilon ^2)$. Below, we describe the detailed derivation that leads finally to the amplitude evolution equations (2.24a)–(2.24c) of the waves constituting a triad.

After substituting the streamfunction (2.7), meridional velocity (2.17), and buoyancy perturbation (2.18) in (2.5), the equation for the $j$th wave can be written as

(2.19)\begin{equation} a_{j}\,\frac{\mathfrak{P}_{j}}{\beta_j}\,\mathcal{L}_j \phi_j ={-}\mathcal{F}_j, \end{equation}

where $\mathfrak {P}_{j} \equiv P_j\,{\rm e}^{-{\rm i} \omega _j t}$, $\mathcal {L}_j$ has been defined in (2.13b), and

(2.20)\begin{align} & \mathcal{F}_j \equiv \underbrace{{\rm i}\,\frac{\partial a_j}{\partial t}\left(\phi_j \mathcal{K}_{j}^{2} - \frac{\partial^{2} \phi_j}{\partial \eta^{2}} \right)\frac{2 \omega_j}{h^{2}} \left(\frac{\mathfrak{P}_{j}}{\beta_j}\right) + 2{\rm i}(N^{2}-\omega_j^2) \left(\frac{\mathcal{K}_{j}}{h}\,\phi_j\,\frac{\partial a_j}{\partial x}\right) \left(\frac{\mathfrak{P}_{j}}{\beta_j}\right)}_{Linear\ term 1} \nonumber\\ &\quad +\underbrace{ {\rm i}(N^{2}-\omega_j^2)\,\frac{\mathcal{K}_{j}}{h} \left[ 2\,\frac{\partial \phi_j}{\partial x} + \frac{\phi_jh}{\mathcal{K}_{j}}\, \frac{\partial }{\partial x}\left(\frac{\mathcal{K}_{j}}{h} \right) -\frac{2\eta}{h}\,\frac{\partial h}{\partial x}\,\frac{\partial \phi_j}{\partial \eta}-2\, \frac{\phi_j}{\beta_j}\,\frac{{\rm d} \beta_j}{{\rm d}\kern0.06em x}\right] \left( a_{j}\,\frac{\mathfrak{P}_{j}}{\beta_j}\right)}_{Linear\ term 2} - {NL}_j. \end{align}

Here, $\mathcal {F}_j$ is the collection of all the linear and nonlinear (${NL}_j$) terms at ${O}(\epsilon ^2)$ that have the phase of the $j$th wave. Equation (2.19) can have a non-trivial solution when $\mathcal {F}_j$ is orthogonal to the adjoint solutions of the linear operator $\mathcal {L}_j$; this procedure is outlined in Craik (Reference Craik1971). The complete mathematical proof for using such a condition is given in detail in Ince (Reference Ince1956, § 9.34). Following Craik (Reference Craik1971), $\mathcal {F}_j$ is multiplied by $\phi _j$ (since $\mathcal {L}_j$ is a self-adjoint operator, $\phi _j$ is also the solution of the adjoint of $\mathcal {L}_j$) and then integrated in the $\eta$ direction inside the boundary limits. This results in

(2.21)\begin{align} & 2\left[{\rm i} \omega_j\,\frac{\partial a_j}{\partial t} \left(\gamma^{(1)}_j \mathcal{K}_{j}^{2} - \gamma^{(2)}_j \right) \frac{1}{h^{2}} + {\rm i}\gamma^{(3)}_j\left( \frac{\mathcal{K}_{j}}{h}\, \frac{\partial a_j}{\partial x}\right)\right]\frac{\mathfrak{P}_{j}}{\beta_j} \nonumber\\ &\quad +{\rm i}\,\frac{\mathcal{K}_{j}}{h}\left[ 2\gamma^{(4)}_j + \frac{h\gamma^{(3)}_j}{\mathcal{K}_{j}}\, \frac{\partial }{\partial x}\left(\frac{\mathcal{K}_{j}}{h} \right) -\gamma^{(5)}_j\,\frac{2}{h}\, \frac{\partial h}{\partial x} - \frac{2\gamma^{(3)}_j}{\beta_j}\, \frac{{\rm d} \beta_j}{{\rm d}\kern0.06em x} \right]a_j\,\frac{\mathfrak{P}_{j}}{\beta_j} = \int_{{-}1}^{0} {NL}_j \phi_j\,{\rm d}\eta, \end{align}

where $\gamma ^{(n)}_j$ are functions that vary in the $x$-direction and are obtained after integration in the $\eta$ direction; $\gamma ^{(n)}_j$ are provided in Appendix A. Up to this point, $\beta _j$ is an arbitrary function, and for convenience, we define $\beta _j$ such that the second square-bracketed term in the left-hand side of (2.21) vanishes identically. It also implies that ‘Linear term 2’ in (2.20) also vanishes identically. In mathematical terms, this means that

(2.22)\begin{equation} \beta_j = \exp \left\{ \int^{x}_0 \frac{h}{2\mathcal{K}_{j}}\, \frac{1}{\gamma^{(3)}_j} \left[ \left(2\gamma^{(4)}_j - \gamma^{(5)}_j\, \frac{2}{h}\,\frac{\partial h}{\partial x} \right) \left(\frac{\mathcal{K}_{j}}{h} \right) + \gamma^{(3)}_j\, \frac{\partial }{\partial x}\left(\frac{\mathcal{K}_{j}}{h} \right) \right]{{\rm d}\kern0.06em x} \right\}. \end{equation}

For constant $N$, $\beta _j$ can be simplified analytically to $\beta _j = h(x)/h(0) = -h(x)/H$, where it is assumed that $h(0)=-H$. We note in passing that the equivalent of the $\beta _j$ functions was derived in Lahaye & Llewellyn Smith (Reference Lahaye and Llewellyn Smith2020) using a different approach. For this particular choice of $\beta _j$, if the amplitudes $a_j$ are $x$-invariant, then the energy flux will be $x$-invariant as well, regardless of the modal shape or depth. More importantly, a wave packet's maximum amplitude does not change when $\beta _j$, given by (2.22), is used in (2.21). This invariance of the maximum value of the $a_j$ with varying $h$ is very useful in estimating wave growth rates in our study, in which a major focus is on wave interactions in a region of varying $h$.

Next, we outline the procedure to obtain $\int _{-1}^{0} {{NL}}_{j} \phi _j \,{\rm d} \eta$ in (2.21) to complete the amplitude evolution equations. The streamfunction, meridional velocity and buoyancy frequency ansatz are substituted in the nonlinear terms of (2.5). The resultant resonant nonlinear terms, after omitting non-resonant terms, can be written in a compact form as

(2.23a)$$\begin{gather} \int_{{-}1}^{0} {NL}_1 \phi_1 \,{\rm d} \eta = \left[\int_{{-}1}^{0}\left(\widehat{{NL}}_{(\varPsi,1)} + \widehat{{NL}}_{(B,1)} + \widehat{{NL}}_{(\mathcal{V},1)}\right)\phi_1 \,{\rm d} \eta\right] \frac{a_3\bar{a}_2}{\beta_2\beta_3}\,\mathfrak{P}_3\bar{\mathfrak{P}}_2, \end{gather}$$
(2.23b)$$\begin{gather}\int_{{-}1}^{0} {NL}_2 \phi_2 \,{\rm d} \eta = \left[\int_{{-}1}^{0}\left(\widehat{{NL}}_{(\varPsi,2)} + \widehat{{NL}}_{(B,2)} + \widehat{{NL}}_{(\mathcal{V},2)}\right)\phi_2 \,{\rm d} \eta\right] \frac{a_3\bar{a}_1}{\beta_1\beta_3}\,\mathfrak{P}_3\bar{\mathfrak{P}}_1, \end{gather}$$
(2.23c)$$\begin{gather}\int_{{-}1}^{0} {NL}_3 \phi_3 \,{\rm d} \eta = \left[\int_{{-}1}^{0}\left(\widehat{{NL}}_{(\varPsi,3)} +\widehat{{NL}}_{(B,3)} + \widehat{{NL}}_{(\mathcal{V},3)}\right)\phi_3 \,{\rm d} \eta\right]\frac{a_1a_2}{\beta_1\beta_2}\,\mathfrak{P}_1\mathfrak{P}_2. \end{gather}$$

We define ${{NL}}_{(*,j)} \equiv \int _{-1}^{0}\widehat {{NL}}_{(*,j)} \phi _j \,{\rm d} \eta$ for convenience; their expressions are provided in Appendix A. Note that ${{NL}}_{(*,j)}$ is used directly in the amplitude evolution equations given below in (2.25b).

The amplitude evolution equations for the three internal gravity waves are obtained finally after equating the left-hand side of (2.21) with its right-hand side, where the latter has been expressed in terms of (A2)(A4):

(2.24a)$$\begin{gather}\frac{\partial a_{1}}{\partial t} + c^{(g)}_{(x,1)}\,\frac{\partial a_{1}}{\partial x} =\mathfrak{N}_1{a}_{3}\bar{a}_{2}\exp\left\{\int_{0}^x{\rm i}(\mathcal{K}_3-\mathcal{K}_1-\mathcal{K}_2)/h \,{{\rm d}\kern0.06em x}' + {\rm i}\,{\rm \Delta} \omega\,t\right\}, \end{gather}$$
(2.24b)$$\begin{gather}\frac{\partial a_{2}}{\partial t} + c^{(g)}_{(x,2)}\, \frac{\partial a_{2}}{\partial x} = \mathfrak{N}_2{a}_{3}\bar{a}_{1}\exp\left\{\int_{0}^x{\rm i}(\mathcal{K}_3-\mathcal{K}_1-\mathcal{K}_2)/h \,{{\rm d}\kern0.06em x}' + {\rm i}\,{\rm \Delta} \omega\,t\right\},\end{gather}$$
(2.24c)$$\begin{gather}\frac{\partial a_{3}}{\partial t} + c^{(g)}_{(x,3)}\,\frac{\partial a_{3}}{\partial x} = \mathfrak{N}_3{a}_{1}{a}_{2}\exp\left\{\int_{0}^x{\rm i}(\mathcal{K}_1+\mathcal{K}_2-\mathcal{K}_3)/h \,{{\rm d}\kern0.06em x}' - {\rm i}\,{\rm \Delta} \omega\,t\right\},\end{gather}$$


(2.25a)$$\begin{gather} c^{(g)}_{(x,j)} = \left[\frac{2{\rm i}\mathcal{K}_j\gamma^{(3)}_j }{h\mathfrak{D}_j}\right],\quad \text{in which}\ \mathfrak{D}_j = 2{{\rm i}\omega_j} \left(\gamma^{(1)}_j \mathcal{K}_{j}^{2} - \gamma^{(2)}_j \right)/h^2, \end{gather}$$
(2.25b)$$\begin{gather}\mathfrak{N}_{j} = \frac{1}{\mathcal{D}_j}\left[{NL}_{(\mathcal{V},j)} + {NL}_{(B,j)} + {NL}_{(\varPsi,j)} \right]. \end{gather}$$

In the above equation,

(2.26ac)\begin{equation} \mathcal{D}_1 = \mathfrak{D}_1\,\frac{\beta_2\beta_3}{\beta_1},\quad \mathcal{D}_2 = \mathfrak{D}_2\,\frac{\beta_1\beta_3}{\beta_2},\quad \mathcal{D}_3 = \mathfrak{D}_3\,\frac{\beta_1\beta_2}{\beta_3}. \end{equation}

The coefficient $c^{(g)}_{(x,j)}$ denotes the (weakly varying) horizontal group speed, and $\mathfrak {N}_j$ denotes the nonlinear coupling coefficient of the $j$th wave; $\mathfrak {N}_j$ determines the rate of energy transfer between the waves. Also, ${\rm \Delta} \omega \equiv \omega _1+\omega _2 - \omega _3$ denotes the detuning in the frequency. The arguments of the exponential terms in (2.24a)–(2.24c) denote both the detuning in the horizontal wavenumber condition and the frequency condition. For a pure resonant triad, $\mathcal {K}_3-\mathcal {K}_1-\mathcal {K}_2 = 0$ and ${\rm \Delta} \omega = 0$. When $\mathcal {K}_3-\mathcal {K}_1-\mathcal {K}_2 \neq 0$ or ${\rm \Delta} \omega \neq 0$, the triad is said to be detuned. The equations are valid only when both ${\rm \Delta} \omega /\omega _j \ll 1$ and ${\rm \Delta} \mathcal {K}/\mathcal {K}_j \ll 1$ are satisfied; that is, the equations are valid only in the vicinity of resonance. Analytical methods have also been developed for studying wave–wave interactions in non-resonant regimes in the presence of a slowly varying background shear flow (e.g. Grimshaw Reference Grimshaw1988Reference Grimshaw1994; Voelker, Akylas & Achatz Reference Voelker, Akylas and Achatz2021), where the wave train can pass through non-resonant regimes and resonant regimes. However, this is not within the scope of this paper. To summarize, amplitude ($a_j$) in the wave amplitude equations can vary because of the group speed term, or the nonlinear term. The group speed term is responsible for the advection of a wave packet, while the nonlinear term is responsible for energy transfer among the waves. Moreover, a wave's energy density changes because of its motion through a region of varying $h$. Both $\phi _j$ and $\beta _j$ are heavily involved in the change in energy density that occurs in a wave due to its motion through a region of varying $h$. Note that the evolution of $a_j$ does not provide complete information of the changes in a wave's quantities. This is because $\varPsi _j = a_j\phi _j/\beta _j P_j\,{\rm e}^{-{\rm i}\omega _j t}$, where $\phi _j$ and $\beta _j$ themselves are functions of $x$.

For a triad, the parent wave is always wave-3, while the daughter waves (subharmonic waves) are wave-1 and wave-2. For self-interactions, we use a different convention; see § 2.2.2. To determine how fast the daughter waves grow, a growth rate parameter $\sigma$ is defined as

(2.27)\begin{equation} \sigma \equiv \sqrt{\mathfrak{N}_1\mathfrak{N}_2 A_3^2}, \end{equation}

where $A_3$ is the parent wave's amplitude, which is held constant. To obtain this expression, the pump wave approximation of Craik, Adam & Stewartson (Reference Craik, Adam and Stewartson1978) is used. Pump wave approximation is a strong assumption that is valid only at initial times where the parent wave has much more energy than the daughter waves. Equation (2.27) reveals that the growth rate is dependent directly on the nonlinear coupling coefficients. If we ignore the nonlinear terms, then (2.24a)–(2.24c) model the movements of internal wavepackets over a mild-slope bathymetry. We emphasize here that wave scattering is not included in these equations. The amplitude variation of internal waves was analysed recently by Lahaye & Llewellyn Smith (Reference Lahaye and Llewellyn Smith2020) (the authors focused on internal wave scattering, which is essentially a linear mechanism). While we have restricted our study to mild-slope conditions, we have extended the previous works by including the physics of (i) finite width wave packets, (ii) nonlinearity, and (iii) detuning in the horizontal wavenumber condition, and hence investigation of both resonant (zero detuning) and near-resonant conditions. In this paper, we focus mainly on the variation of detuning, and growth rates (using pump wave approximation) with $h/H$ for wave–wave interactions. Even though (2.24a)(2.24c) allow finite width wave packets, we do not discuss them significantly since they have been studied in Gururaj & Guha (Reference Gururaj and Guha2020). The combined effect of nonlinear coupling coefficients, group speed and detuning has been discussed in Gururaj & Guha (Reference Gururaj and Guha2020).

The main results of scaling analysis, detailed in Appendix B, are summarized here. The relation between the small parameters is given by

(2.28)\begin{equation} \epsilon_{t} \sim \frac{{\mathfrak{N}}}{\omega }\,\epsilon_a - {\widehat{c}_g} \epsilon_x, \end{equation}

where ${\widehat {c}_g}$ is a non-dimensional term that gives a scale of the group speed. Equation (2.28) provides the scaling for ‘Linear term 1’, and ${NL}_j$ given in (2.20). These are also the final terms that are present in the wave amplitude equations (2.24a)(2.24c). The wave amplitude ($a_j$) can evolve due to the group speed term or the nonlinear term. Note that if $\mathfrak {N}$ or $\epsilon _a$ is reduced (implying that nonlinear coupling coefficients or amplitude is reduced), then we can expect nonlinear effects to decrease. However, if $\epsilon _x$ is reduced (which means that packet width is increased), then the effect of group speed, which advects the packets, is reduced.

2.2.2. Amplitude evolution equations for self-interaction in non-uniform stratification

Self-interactions can be considered as a special case of triad interactions. During resonant self-interactions, an internal wave gives its energy spontaneously to another internal wave that has twice its frequency and horizontal wavenumber (Wunsch Reference Wunsch2017). In non-uniform stratification, a resonant self-interaction occurs when both $(\omega, k)$ and $(2\omega, 2k)$ satisfy the dispersion relation. The evolution equations for the self-interaction of a mode can be obtained from (2.24a)(2.24c) after some straightforward modifications. The complete set of governing equations for the self-interaction of a mode in the presence of a mild-slope bathymetry $h$ is

(2.29a) $$\begin{gather} \frac{\partial a_{(3,s)}}{\partial t} + c^{(g)}_{(x,3)}\, \frac{\partial a_{(3,s)}}{\partial x} = \mathcal{N}_3\,{a}_{(1,s)}^2 \exp\left\{{\int_{0}^x{\rm i}(2\mathcal{K}_1-\mathcal{K}_3)/h \,{{\rm d}\kern0.06em x}' - {\rm i}\,{\rm \Delta} \omega_s\,t}\right\},\end{gather}$$
(2.29b)$$\begin{gather}\frac{\partial a_{(1,s)}}{\partial t} + c^{(g)}_{(x,1)}\,\frac{\partial a_{(1,s)}}{\partial x} = \mathcal{N}_1\,{a}_{(3,s)}\,\bar{a}_{(1,s)} \exp\left\{{\int_{0}^x{\rm i}(\mathcal{K}_3-2\mathcal{K}_1)/h\,{{\rm d}\kern0.06em x}' + {\rm i}\, {\rm \Delta} \omega_s\,t}\right\},\end{gather}$$

where the subscript ‘$s$’ denotes self-interaction. Moreover, ${\rm \Delta} \omega _s = 2\omega _1 - \omega _3$. Unlike the triad case, the parent wave for self-interaction is wave-1, while the daughter (superharmonic) wave is wave-3. The notation throughout this paper follows the convention that wave-3 always has the highest frequency (hence for triads, wave-3 becomes the parent wave). The functions $c^{(g)}_{(x,j)}$ are the same as the expressions given in (2.25a). The functions $\mathcal {N}_j$, which are the nonlinear coupling coefficients for the self-interaction process, are given by

(2.30a)$$\begin{gather} \mathcal{N}_{1} = \mathfrak{N}_2, \end{gather}$$
(2.30b)$$\begin{gather}\mathcal{N}_{3} = \frac{2\mathcal{K}_{1}^3}{h^4\mathcal{D}_3}\left(\frac{\varGamma^{(4)}}{\omega_1}\right) -\frac{2f^2}{h^4\mathcal{D}_3}\left(\frac{\varGamma^{(3)}_1\mathcal{K}_1}{\omega_1}\right) +\frac{\mathcal{K}_{1}\omega_{3}}{h^4\mathcal{D}_3}\left(\zeta_1\omega_1^2\varGamma^{(1)}_1 - \zeta_1\varGamma^{(2)}_1-\varGamma^{(3)}_1\right), \end{gather}$$

where $\zeta _j \equiv \mathcal {K}_j^2/{(\omega _j^2-f^2)}$ is defined for convenience. Here, all $\varGamma, \mathcal {D}_j$ terms in (2.30a) and (2.30b) are evaluated using (A5) and (2.26ac) by simply considering all ‘2’ subscripts as ‘1’ – for example, substituting $\beta _1$ for $\beta _2$ in $\mathcal {D}_j$, and similarly substituting $\phi _1$ for $\phi _2$ in $\varGamma$ expressions. This is because in self-interaction, wave-2 is the same as wave-1.

Equations (2.29a)(2.29b) can predict the growth of the daughter wave and the consequent decay of the parent wave. For obtaining the growth rate of the daughter waves, we use the pump wave approximation and hence treat the parent wave's amplitude $a_{(1,s)}$ as constant. This yields (assuming plane waves in the $x$-direction)

(2.31)\begin{equation} a_{(3,s)} = \left[\mathcal{N}_3\,{a}_{(1,s)}^2 \right] t, \end{equation}

where the term in square brackets denotes the growth rate. From this equation it is evident that $\mathcal {N}_3$ acts as a proxy to the growth rate.

2.3. Energy evaluation

The time-averaged energy density for an internal gravity wave over its time period is given by

(2.32)\begin{equation} \langle {TE}_{j}\rangle = \frac{\omega_{j}}{2{\rm \pi}}\int_{0}^{{2{\rm \pi}}/{\omega_{j}}} \frac{\rho_{0}}{2}\left[\left( \frac{\partial \psi_{j}}{\partial z}\right)^{2} + \left(\frac{\partial \psi_{j}}{\partial x}\right)^{2} + v_j^2 + \frac{b_j^2}{N^2} \right] {\rm d} t. \end{equation}

The domain-integrated total energy is given by

(2.33)\begin{equation} \widehat{{TE}}_{j} =\int_{0}^{D} \int_{h}^{0}\langle {TE}_{j}\rangle\,{\rm d} z\,{{\rm d}\kern0.06em x} = \int_{0}^{D} \int_{{-}1}^{0} \langle TE_{j}\rangle({-}h(x))\,{\rm d}\eta\,{{\rm d}\kern0.06em x}. \end{equation}

After some simplification, we arrive at

(2.34)\begin{equation} \widehat{{TE}}_{j} = \int_{0}^{D} \int_{{-}1}^{0} -\frac{2}{h} \left[ \mathcal{K}_{j}^2 \phi_j^{2} + \left(\frac{\partial \phi_j}{\partial \eta}\right)^2 \right] |a_j|^2\, \frac{\rho_0}{\beta_j^{2}}\,{\rm d}\eta\,{{\rm d}\kern0.06em x}, \end{equation}

where $D$ is the length of the domain in the $x$-direction. We non-dimensionalize $\widehat {{TE}}_{j}$ with the initial energy of the parent wave (abbreviated as ‘$Pw$’): ${E}_{j} = {\widehat {{TE}}_{j}}/{\widehat {{TE}}_{Pw}|_{t = 0}}$. Note that $Pw=3$ (i.e. wave-3) for triads, and $Pw=1$ (i.e. wave-1) for self-interactions.

3. Triad interactions in a uniform stratification in the presence of a mild-slope bathymetry

In this section, we consider resonant and near-resonant triads in a uniform background stratification in the presence of a mild-slope bathymetry. Here, we will consider briefly the horizontal wavenumber triad condition in uniform stratification as $h$ is varied. Without any loss of generality, the triad condition for the horizontal wavenumber is

(3.1)\begin{equation} \mathcal{K}_{3} = \mathcal{K}_{1} + \mathcal{K}_{2}. \end{equation}

However, using (2.13b), it can be seen that the $\mathcal {K}_j$ are constants for uniform stratification. As a result, triad conditions are satisfied everywhere in the domain, provided that the conditions are satisfied perfectly for any given domain height.

3.1. Effect of bathymetry on the nonlinear coupling coefficients of resonant triads

Here, we focus on the nonlinear coupling coefficients in resonant triad interactions (i.e. no detuning) in the presence of a uniform stratification and a weakly varying bathymetry. Equation (2.27) revealed that the growth rate of the daughter waves is dependent on the nonlinear coupling coefficients. For constant $N$, the nonlinear coupling coefficients ($\mathfrak {N}_j$) in (2.25b) can be further simplified as

(3.2a)\begin{align} \mathfrak{N}_{1} &= \frac{{\rm i} H}{2h^2\omega_1\kappa_1\kappa_2\kappa_3} \left[{N^{2}(\mathcal{K}_{3}-\mathcal{K}_{2})}\left\{\left(\frac{\mathcal{K}_{3}}{\omega_{3}}- \frac{\mathcal{K}_{2}}{\omega_{2}}\right) \left(\mathcal{K}_2\mathcal{M}_3-\mathcal{K}_3\mathcal{M}_2\right)\right\} \right. \nonumber\\ &\quad +\omega_1\left\{(\mathcal{K}_{2}\mathcal{M}_{3}-\mathcal{K}_{3}\mathcal{M}_{2}) \left({\mathcal{M}^{2}_{2}}+\mathcal{K}^{2}_{2}-\mathcal{K}^{2}_{3}-{\mathcal{M}^{2}_{3}}\right) \right\}, \nonumber\\ &\quad \left.{}+{f^2}\left(\mathcal{M}_3-\mathcal{M}_2\right) \left\{\left( \frac{\mathcal{K}_3}{\omega_3}+\frac{\mathcal{K}_2}{\omega_2} \right) \left(\mathcal{M}_3\mathcal{M}_2\right)-\left( \frac{\mathcal{K}_2}{\omega_3}+ \frac{\mathcal{K}_3}{\omega_2} \right) \left(\mathcal{M}^2_2+\mathcal{M}^2_3\right)\right\}\right], \end{align}
(3.2b)\begin{align} \mathfrak{N}_{2} &= \frac{{\rm i} H}{2h^2\omega_2\kappa_1\kappa_2\kappa_3} \left[{N^{2}(\mathcal{K}_{3}-\mathcal{K}_{1})}\left\{\left(\frac{\mathcal{K}_{3}}{\omega_{3}}- \frac{\mathcal{K}_{1}}{\omega_{1}}\right)\left(\mathcal{K}_1\mathcal{M}_3-\mathcal{K}_3\mathcal{M}_1\right)\right\} \right. \nonumber\\ &\quad +\omega_2\left\{(\mathcal{K}_{1}\mathcal{M}_{3}-\mathcal{K}_{3}\mathcal{M}_{1}) \left({\mathcal{M}^{2}_{1}}+\mathcal{K}^{2}_{1}-\mathcal{K}^{2}_{3}-{\mathcal{M}^{2}_{3}}\right) \right\}, \nonumber\\ &\quad \left.{}+{f^2}\left(\mathcal{M}_3-\mathcal{M}_1\right) \left\{\left( \frac{\mathcal{K}_1}{\omega_1}+\frac{\mathcal{K}_3}{\omega_3} \right) \left(\mathcal{M}_1\mathcal{M}_3\right)-\left( \frac{\mathcal{K}_1}{\omega_3}+ \frac{\mathcal{K}_3}{\omega_1} \right) \left(\mathcal{M}^2_1+\mathcal{M}^2_3\right)\right\}\right], \end{align}
(3.2c)\begin{align} \mathfrak{N}_{3} &= \frac{{\rm i} H}{2h^2\omega_3\kappa_1\kappa_2\kappa_3} \left[{N^{2}(\mathcal{K}_{1}+\mathcal{K}_{2})}\left\{\left(\frac{\mathcal{K}_{1}}{\omega_{1}}- \frac{\mathcal{K}_{2}}{\omega_{2}}\right)\left(\mathcal{K}_2\mathcal{M}_1-\mathcal{K}_1\mathcal{M}_2\right)\right\} \right. \nonumber\\ &\quad +\omega_3\left\{(\mathcal{K}_{2}\mathcal{M}_{1}-\mathcal{K}_{1}\mathcal{M}_{2}) \left({\mathcal{M}^{2}_{1}}+\mathcal{K}^{2}_{1}-\mathcal{K}^{2}_{2}-{\mathcal{M}^{2}_{2}}\right)\right\} \nonumber\\ &\quad \left.{}+{f^2}\left(\mathcal{M}_1+\mathcal{M}_2\right) \left\{\left( \frac{\mathcal{K}_1}{\omega_2}+\frac{\mathcal{K}_2}{\omega_1}\right) \left(\mathcal{M}^2_1+\mathcal{M}^2_2\right) - \left( \frac{\mathcal{K}_1}{\omega_1}+ \frac{\mathcal{K}_2}{\omega_2} \right) \left(\mathcal{M}_1\mathcal{M}_2\right)\right\}\right], \end{align}

where $\kappa _j = \sqrt { \mathcal {M}^2_j+\mathcal {K}^2_j}$. Note that the above expressions are obtained only when the vertical wavenumber condition is satisfied. The terms inside the square brackets are constant and hence do not vary with the bathymetry $h$. The fact that $\mathcal {K}_j$ and $\mathcal {M}_j$ are constants for a constant $N$ is given in (2.14). For constant $N$, $\beta _j = -h(x)/H$, which has been used in (3.2a)(3.2c), and this results finally in $\mathfrak {N}_j \propto 1/h^2$. Hence for waves travelling from a given fluid depth to a lesser depth (i.e. as the waves climb up a seamount), the nonlinear coupling coefficients, and hence the growth rates, increase following the inverse square rule.

In summary, for a uniform stratification, if three modes satisfy the resonant triad condition at a particular domain height, then they would satisfy the resonant triad condition for any domain height. Moreover, we also showed that the nonlinear coupling coefficients increase (decrease) as the fluid depth decreases (increases) following an inverse square law.

4. Triad interactions and self-interactions in a non-uniform stratification in the presence of a mild-slope bathymetry: detuning effects

In § 3, it was shown that in the presence of a uniform stratification, if the triad condition is satisfied between three modes at a particular $h$, then it is satisfied for all $h$. However, in non-uniform stratification, such a simple outcome is not possible. In certain types of triads, there can be a heavy mismatch in the horizontal wavenumber condition as the waves involved in the triad interact in a region of varying domain height. This may affect the energy transfer between the waves.

In this section, we study the factors that decide the detuning (or mismatch) between the horizontal wavenumber of the waves as $h$ is varied in the presence of non-uniform stratification. Here, as well as in the rest of this paper, we will consider a Gaussian function to represent the buoyancy frequency

(4.1)\begin{equation} N(z) = N_{b} + N_{{max}}\exp\left[-\left\{(z-z_{c})/W_{p}\right\}^{2}\right], \end{equation}

where the parameters $N_b, N_{{max}}, W_{p}$ and $z_{c}$ are varied. This kind of profile (see figure 1b) is a simplified representation of oceanic stratification and is used widely in the literature; see Grisouard, Staquet & Gerkema (Reference Grisouard, Staquet and Gerkema2011), Mathur et al. (Reference Mathur, Carter and Peacock2014) and Varma & Mathur (Reference Varma and Mathur2017). We choose stratification profiles such that the pycnocline is above the topography. If the topography cuts the pycnocline, then internal wave scattering may be significant, as shown in Hall, Huthnance & Williams (Reference Hall, Huthnance and Williams2013).

4.1. Effect of varying $h$ on the horizontal wavenumber condition for waves satisfying $f \ll \omega _j \ll N_b$

First, we study the class of triads for which the angular frequencies of the constituent waves obey the condition $f \ll \omega _j \ll N_b$. It is assumed that the parent wave (angular frequency $\omega _3$) gives its energy to two subharmonic daughter waves of angular frequencies $\omega _1$ and $\omega _2$, respectively; that is, the conditions $\omega _1<\omega _3$ and ${\omega _2}<{\omega _3}$ are always assumed. A parameter $\alpha \in (0,1)$ is defined such that $\omega _1 = \alpha \omega _3$ and $\omega _2 = (1-\alpha )\omega _3$. Two different types of interactions, Class-1 and Class-2, are defined for which a parent wave can form a triad with the subharmonic daughter waves.

4.1.1. Class-1 interactions

We consider three waves with angular frequencies $\omega _3, \omega _1, \omega _2$ such that $\omega _3 = \omega _1 + \omega _2$. Furthermore, we assume that at a particular $h$, the horizontal wavenumber condition is satisfied between mode $i$ of wave-1, mode $j$ of wave-2, mode $k$ of wave-3, i.e.

(4.2)\begin{equation} \mathcal{K}_{3(k)} = \mathcal{K}_{1(i)} + \mathcal{K}_{2(j)}, \end{equation}

where $i$, $j$ and $k$ are not all equal. This constitutes a Class-1 interaction. Now, if the stratification profile changes (the stratification profile will change in $x$$\eta$ coordinates provided that $h$ is varying), then the wavenumbers ${\mathcal {K}}_{1(i)}, {\mathcal {K}}_{2(j)}, {\mathcal {K}}_{3(k)}$ will also change. However, for a given change in $h$, all the wavenumbers may not change in a way such that condition (4.2) is satisfied. For example, if $\mathcal {K}_{1(i)} = \mathrm {func}(h)$, then it is possible that $\mathcal {K}_{2(j)} \neq c\,\mathrm {func}(h)$, where $c$ and $\mathrm {func}(\,)$ denote an arbitrary constant and function, respectively. Therefore, even though the triad condition may be satisfied at a particular $h$, it may not be satisfied for all $h$. Hence Class-1 triads might get detuned as they interact in a region of varying $h$.

To measure the detuning (or mismatch) in the horizontal wavenumber, we define a new variable ${\rm \Delta} \mathcal {K}$:

(4.3)\begin{equation} {\rm \Delta} \mathcal{K} \equiv \frac{\mathcal{K}_{3(k)}-\mathcal{K}_{1(i)}- \mathcal{K}_{2(j)}}{\mathcal{K}_{{min}}},\end{equation}

where $\mathcal {K}_{{min}}$ is the minimum wavenumber of the three wavenumbers at a particular $x$-coordinate. This ${\rm \Delta} \mathcal {K}$ acts basically as a non-dimensional measure of the detuning between the waves, and for a resonant triad, ${\rm \Delta} \mathcal {K} = 0$.

We now study how different (non-dimensional) wavenumbers ${\mathcal {K}}_{3(n)}$ of frequency $\omega _3$ change as $h$ is varied in the presence of a non-uniform stratification. To obtain ${\mathcal {K}}_{3(n)}$ for a given stratification profile, we solve (2.13b) for $h/H \in [-1,-0.2]$. The functional form of $h$, as long as it is mildly varying, does not influence the wavenumbers or detuning at a particular $h$. The non-uniform stratification profile given by (4.1) is used throughout this paper. The stratification profiles are chosen such that $N_{{max}}=(2 N_b, 4 N_b,\ldots,12N_b)$, $W_p=(H/200, 2H/200,\ldots,5H/200)$ and $z_c=(H/80, H/40, H/20, H/10)$; and we consider all possible ($120$) combinations. Moreover, $\omega _3 = 0.1N_b$ and $f=0$ are used consistently for all combinations. Figure 3 shows the variation of $\widehat {{\mathcal {K}}}_{3(n)} \equiv {\mathcal {K}}_{3(n)}(h)/{\mathcal {K}}_{3(n)}(H)$ with $h/H$ for different modes $n$. Figures 3(ac) uses the stratification profile $N^{(1)}$ with the parameters $N_{{max}} = 2N_{b}$, $W_{p} = H/200$ and $z_c = H/80$; figures 3(df) use the profile $N^{(2)}$ given by $N_{{max}} = 10N_{b}$, $W_{p} = H/50$ and $z_c = H/10$. Note that $N^{(1)}$ has a sharp pycnocline, while $N^{(2)}$ has a larger $W_p$ resulting in a wider pycnocline. For profiles where all three parameters are low (e.g. $N^{(1)}$), $\widehat {{\mathcal {K}}}_{3(1)}$ is nearly constant for some range of $h/H$ and then starts decreasing. This can be seen in figure 3(a), where the first five modes exhibit this behaviour. Moreover, for profiles where all $z_c, W_p, N_{{max}}$ are high (e.g. $N^{(2)}$), $\widehat {{\mathcal {K}}}_{3(1)}$ decreases almost linearly with $h/H$, as can be seen clearly in figure 3(d). For any profile, $\widehat {{\mathcal {K}}}_{3(1)}$ always monotonically decreases as the fluid depth is reduced for $h/H \in [-1, -0.2]$. However, this behaviour does not hold for any mode other than mode-$1$. For example, for $z_c = H/10$ (regardless of $W_p, N_{{max}}$), $\widehat {{\mathcal {K}}}_{3(2)}$ increases for some $h/H$ as fluid depth is reduced; see figure 3(d) (blue curve). Similar behaviour is observed for modes $3$, $4$ and $5$ when $W_p$ is low. In summary, the variation of $\widehat {{\mathcal {K}}}_{3(1)}$ with $h/H$ can be different from that of the higher mode's wavenumber, which can result in detuning.

Figure 3. Variation of $\widehat {{\mathcal {K}}}_{3(n)}$ with $h/H$ for two different stratification profiles. For stratification profile $N^{(1)}$: (a) modes 1–5, (b) modes 6–10, and (c) modes 11–15. For stratification profile $N^{(2)}$: (d) modes 1–5, (e) modes 6–10, and ( f) modes 11–15.

For profiles with high $W_p$, the $\widehat {{\mathcal {K}}}_{3(n)}$ for $n>10$ start to collapse on each other; see figure 3f). In such scenarios, since $\widehat {{\mathcal {K}}}_{3(n)}$ remains nearly the same, ${\rm \Delta} \mathcal {K}$ will not be induced by the difference in higher modes’ $\widehat {{\mathcal {K}}}_{3(n)}$. In general, it was observed that as $W_p$ is reduced, $n$ has to be higher for the modes to collapse on each other.

The interaction of a mode-1 internal wave (wave-3) with different modes in the presence of two different non-uniform stratification profiles is considered next. These profiles are a part of the $120$ profiles that we mentioned already. Sample results are shown in figure 4 in which the frequencies and stratification profile parameters are as follows:

  1. (i) $N^{(3)}$: $\omega _3 = 0.1N_b$, $f=0$,$N_{{max}} = 10N_b$, $W_p = H/100$ and $z_c = H/10$;

  2. (ii) $N^{(4)}$: $\omega _3 = 0.1N_b$, $f=0$,$N_{{max}} = 10N_b$, $W_p = H/50$ and $z_c = H/20$.

Figure 4. Variation of detuning (${\rm \Delta} \mathcal {K}$) with $h/H$ for various modal interactions pertaining to the stratification profiles (a) $N^{(3)}$, and (b) $N^{(4)}$. The legends indicate what daughter waves were involved in the triad interactions, where $\alpha \equiv \omega _1/\omega _3$.

For each profile, we have shown five different modal interactions. Figure 4 reveals clearly that the detuning can be quite sensitive to the changes in the domain height.

4.1.2. Class-2 interactions: a special case of triad interactions

The interaction in Class-2 is between the $n$th modes (where $n\in \mathbb {Z}^+$) of different waves constituting a triad. For example, if mode-$1$ with frequency $\omega _1$, mode-$1$ with frequency $\omega _2$, and mode-$1$ with frequency $\omega _3$ form a triad, then it is classified as a Class-2 interaction. This kind of triad is possible when $f \ll \omega _j \ll N_b$. To show how this interaction is possible, we consider the eigenproblem concerning the $n$th mode of the $j$th wave:

(4.4)\begin{equation} \frac{\partial^2 \phi_{j(n)} }{\partial \eta^2} + \mathcal{K}_{j(n)}^{2} \chi_j^2\phi_{j(n)} \approx \frac{\partial^2 \phi_{j(n)} }{\partial \eta^2} + \left(\frac{\mathcal{K}_{j(n)}}{\omega_j}\right)^2N^2\phi_{j(n)}= 0, \end{equation}

where we used $\chi _j \approx N/\omega _j$ (under the approximation $f \ll \omega _j \ll N_b$), and the system is solved using the boundary conditions $\phi _{j(n)} = 0$ at $\eta = 0$ and $\eta = -1$. However, by Sturm–Liouville theory, for a given operator (here $\partial ^2/\partial \eta ^2$) and weight function (here $N(z)^2$), the $n$th eigenvalue (here $\mathcal {K}_{j(n)} /\omega _j$) is unique, i.e. $\mathcal {K}_{j(n)} /\omega _j=\mathrm {constant}$ for all $j$. Therefore, if the triad condition for frequency $\omega _3 = \omega _2 + \omega _1$ is valid, then this implies automatically the validity of the wavenumber condition $\mathcal {K}_{3(n)} = \mathcal {K}_{1(n)} + \mathcal {K}_{2(n)}$.

The situation mentioned above is true for all stratification profiles satisfying $f \ll \omega _j \ll N(z)$ (at all $z$-locations). This is especially important because in the presence of a bathymetry, the stratification profile changes in the $x$-direction in $x$$\eta$ coordinates. However, for Class-2 interaction, all three non-dimensional wavenumbers (eigenvalues) are the same functions of $h$ since they are the same eigenvalues divided by their frequency. Thus the resonant triad condition will be still be satisfied even if $h$ is varied significantly (i.e. variation of $h$ will not cause detuning). Note that in the parameter regime $f \ll \omega _j \ll N(z)$, only Class-2 self-interactions were observed in numerical experiments of Sutherland (Reference Sutherland2016), hence Class-2 triads may always be dominated by self-interactions, resulting in the parent wave's energy transfer to the superharmonics instead of subharmonics. As a result, Class-2 triads may not be as relevant practically as Class-2 self-interactions.

4.2. Effect of bathymetry on horizontal wavenumber condition for Class-1 self-interaction

Detuning can also be introduced during a self-interaction process as $h$ is varied. Following the same terminology as before, we classify self-interactions as Class-1 and Class-2. As shown by Wunsch (Reference Wunsch2017), Class-2 self-interactions will always be slightly detuned, where the detuning increases as $f$ increases. This is due to the fact that in a non-uniform stratification, if the $n$th mode of $\omega$ satisfies the dispersion relation, then the $n$th mode of $2\omega$ will be able to satisfy it only approximately.

Following (4.3), the detuning for a self-interaction process is defined as

(4.5)\begin{equation} {\rm \Delta} \mathcal{K}_s = \frac{\mathcal{K}_{3(k)}-2\mathcal{K}_{1(i)}}{\mathcal{K}_{1(i)}}, \end{equation}

where wave-3 is the superharmonic (daughter) wave, while wave-1 is the parent wave, i.e. $\omega _1=\omega _3/2$ (following the convention used throughout this paper that wave-3 has the highest frequency). The amplitude evolution equations for a self-interaction process are discussed in § 2.2.2.

For the range $f \ll \omega _j \ll N(z)$, the Class-2 self-interaction process will follow principles similar to those outlined in § 4.1.2. As mentioned in § 4.1.2, significant variations in $h$ for this frequency range will not introduce detuning in Class-2 self-interactions. We now study the other end of the parameter space, where $(N^{2}-\omega _j^{2})\approx N^{2}$, which was the basic approximation used in much of our analysis in § 4.1, is no longer valid. Hence out of Class-1 and Class-2 self-interactions, only the latter are possible. This would mean that as the domain height changes, the detuning introduced could be significant. Interestingly, though, if the wavenumbers involved in the self-interaction change with $h/H$ in a similar way, then the detuning is insignificant; see figure 5. The frequencies and the stratification profile parameters used here are:

Figure 5. (a) Variations of $\mathcal {K}_{3(2)}$ and $\mathcal {K}_{1(3)}$ with $h/H$ for the parameter set $N^{(6)}$. (b) Variation of detuning with $h/H$ for the same case. (c) The detuning for three different self-interaction combinations for the parameter set $N^{(7)}$. Here, the notation (Pa)(Db) implies ‘Parent wave’ with mode-a, and ‘Daughter wave’ with mode-b.

  1. (i) $N^{(6)}$: $N_{{max}} = 10N_b$, $W_p = H/100$ and $z_c = H/10$;

  2. (ii) $N^{(7)}$: $N_{{max}} = 10N_b$, $W_p = H/100$ and $z_c = H/20$;

and $f=0$ always. Figure 5(a) uses the set $N^{(6)}$, and shows the variation of the horizontal wavenumber of mode-2 ($\omega _3 = 0.89N_b$) and mode-3 ($\omega _1 = \omega _3/2$). These modes satisfy the condition for resonant self-interaction. We observe that these two wavenumbers behave quite similarly for a wide range of $h/H$, hence the detuning, shown in figure 5(b), is small (and constant for an appreciable range), in spite of the fact that it is a Class-1 interaction. The same phenomenon is also shown for several other self-interaction combinations in figure 5(c), where the parameter set $N^{(7)}$ is used.

The detuning for all the combinations shown stays constant for a certain range of $h/H$. We note in passing that Class-1 triad interactions may also give rise to a small detuning for a range of $h/H$, provided that all the modes involved behave in a similar way. However, this is a more stringent condition than a self-interaction process, where only two waves are involved. Even though equations derived in § 2 are valid only when ${\rm \Delta} \mathcal {K} \ll 1$ (or ${\rm \Delta} \mathcal {K}_s \ll 1$), there are a significant number of interactions where ${\rm \Delta} \mathcal {K}$, or ${\rm \Delta} \mathcal {K}_s$, is a small quantity even for ${O}(1)$ changes in depth and the wavenumber. For example, for interactions shown in this section, and for Class-2 interactions, detuning can stay as a small quantity even for ${O}(1)$ changes in depth. However, we do note that in several triad interactions, detuning can be sensitive to $h$, and in those cases, ${O}(1)$ changes in depth cannot be modelled accurately by the wave amplitude equations.

To summarize, in the presence of a non-uniform stratification, we divide triad and self-interactions into two classes: Class-1 and Class-2. Class-1 interactions contain waves whose mode numbers are not all the same, while Class-2 interactions contain waves that are the $n$th modes of their respective frequencies. Class-1 interactions may undergo detuning with the variation in $h$, irrespective of the frequency. However, interestingly, certain Class-1 self-interactions do not undergo detuning as $h$ is varied inside a certain range. For both triads and self-interactions, Class-2 interactions can exist only for $f\ll \omega _j \ll N(z)$, and do not get detuned as $h$ is varied.

5. Variation of growth rates and nonlinear coupling coefficients with depth for non-uniform stratification

In this section, we focus on the effects of domain height variation on the growth rate ($\sigma$) of triads, and the nonlinear coupling coefficient $\mathcal {N}_3$, which provides a measure of the growth of the daughter wave in a self-interaction. The non-uniform stratification profile (4.1) will be used in this section.

5.1. Variation of growth rates with domain height for triads

Triad interactions are important for the decay of internal waves near the $28.9^{\circ }$ latitude (MacKinnon & Winters Reference MacKinnon and Winters2005; MacKinnon et al. Reference MacKinnon, Alford, Sun, Pinkel, Zhao and Klymak2013), specifically the mode-1 wave, which is the most energy containing mode (Vic et al. Reference Vic, Garabato, Green, Waterhouse, Zhao, Mélet, De Lavergne, Buijsman and Stephenson2019). Here, we study this phenomenon in the high-latitude region ($f/\omega _3\geq 0.3$) for varying $N(z)$ and $h$. The mode-1 wave (which, being the parent wave, is wave-3) can decay forming various triad combinations; we restrict the subharmonic daughter waves (wave-1 and wave-2) up to mode-50. Moreover, for studying growth rates in the presence of varying $h$, the triads are identified separately at different $h/H$. This is because a triad combination at a particular $h/H$ value may not satisfy the horizontal wavenumber condition at a different $h/H$ (as explained in § 4). Three main branches of triads are considered here for the mode-1 internal wave:

(5.1)\begin{equation} \underbrace{|\mathcal{K}_3| = |\mathcal{K}_2| - |\mathcal{K}_1|}_{{Branch\text{-}1}} \quad {or}\quad \underbrace{|\mathcal{K}_3| = |\mathcal{K}_1| - |\mathcal{K}_2|}_{Branch\text{-}2} \quad {or}\quad \underbrace{|\mathcal{K}_3| \approx |\mathcal{K}_1| + |\mathcal{K}_2|}_{Branch\text{-}3}. \end{equation}

For Branch-1(2) triads, the wavenumber of wave-2(1) is larger in magnitude than that of wave-1(2). The only possible Branch-3 interaction is a Class-2 interaction, where both the daughter waves are also mode-1 of their respective frequencies. However, this interaction, like the Class-2 self-interaction, also undergoes heavy detuning for high $f$ values. Therefore, Branch-3 being an inefficient energy transfer pathway, we restrict our focus to Branch-1 and Branch-2. Triads are studied for $f/\omega _3 = (0.3, 0.4, 0.45)$ in the presence of various stratification profiles. The triads are computed for $\alpha \in [0.31,0.5]$, $\alpha \in [0.41,0.5]$ and $\alpha \in [0.455,0.5]$ for $f/\omega _3 = 0.3$, $0.4$ and $0.45$, respectively (see § 4.1 for the definition of $\alpha$).

Figure 6 shows the non-dimensionalized growth rate contour for a mode-1 wave. All growth rates $\sigma$ are non-dimensionalized with a reference growth rate value $\sigma _{{ref}}$, where the latter denotes the maximum growth rate for all Branch-1 triads at $h=-H$ (hence the value of $A_3$ does not impact the results shown). The frequency of the mode-1 wave is $\omega _3/N_b = 0.2$, while $f/\omega _3 = 0.4$ is taken. The stratification profile is given by:

  1. (i) $N^{(8)}$: $N_{{max}} = 10N_b$, $W_p = H/50$, $z_c = H/20$.

Figure 6. Contours of non-dimensional growth rate ($\sigma /\sigma _{{ref}}$) of triads formed between mode-$n$, mode-$m$ and mode-1 (i.e. wave-3). Blue and red, respectively, represent Branch-1 and Branch-2 triads, and both colours represent positive values. For Branch-1 triads, mode-$m(n)$ is wave-1(2), while for Branch-2 triads, mode-$n(m)$ is wave-1(2).

Branch-1(2) triads have the higher (lower) frequency daughter wave propagating in the same direction as the parent wave. Figure 6 reveals that from both branches, the highest growth rates are centred around $n \approx m$. However, the majority of the white region contains resonant triads, but their growth rates are significantly lower in comparison to that clustered around $n \approx m$. Note that the central region is asymmetric between Branch-1 and Branch-2 triads, and this is purely a consequence of the internal wave's dispersion relation. When the lower frequency daughter wave (wave-1) travels in the same direction as the parent wave (i.e. Branch-2), wave-1's mode number ($n$) should always be higher than wave-2's mode number ($m$) for the triad condition to be satisfied. However, for Branch-1, where wave-2 travels in the same direction as the parent wave, the mode number of wave-2 ($n$) need not be higher than wave-1's mode number ($m$).

The clustering around $n \approx m$ is observed consistently for any setting or stratification profile considered in our study. As a result, instead of focusing on all possible triads, we choose a specific line of interaction near the central region and plot the growth rate along that line of interaction. For example, the interaction lines $(n,n)$, $(n+1,n)$ and $(n,n+1)$ are plotted for $n\in (1,50)$ in figures 7(ac) for Branch-1 triads, and $(n+4,n)$ in figure 7(d) for Branch-2 triads. The notation $(a,b)$ means wave-1(2) is mode-$a(b)$. The notation is same for both branches. The dominant nature of the interaction lines $(n,n)$, $(n+1,n)$ and $(n,n+1)$ has also been observed in Young et al. (Reference Young, Tsang and Balmforth2008) while studying the stability of mode-1 internal waves in the presence of near inertial daughter waves (with frequency $f$). Furthermore, figure 7 also reveals that the different lines are sensitive to $h$. For completeness, we explore another stratification profile given by:

  1. (i) $N^{(9)}$: $N_{{max}} = 10N_b$, $W_p = H/50$, $z_c = H/80$;

and the corresponding plots are in figure 8. Both figures 7 and 8 show that the growth rates along different lines of interaction have a significant oscillatory nature with $n$. In general, line $(n,n)$ has the largest amplitude of oscillations. More importantly, the growth rate of a modal combination can change significantly as $h$ changes. For example, figure 7(a) shows that the most unstable modal combination at $h=-H$ is $(5,5)$. However, for $h=-0.8H$, the most unstable triad is the modal combination $(4,4)$. Moreover, the combination $(5,5)$ has approximately $0.25$ times the $(4,4)$ growth rate at $h=-0.8H$. This behaviour can be seen for the line $(n,n)$ in both figures 7 and 8. This means effectively that the growth rate of certain daughter wave combinations can be sensitive to changes in $h$ (especially the combinations that involve lower modes). Such combinations may not be effective in a region of varying $h$ because of the significant drop in the growth rates. However, sensitivity to $h$ is reduced slowly as the mode number is increased for both the branches. Even though Branch-2 triads have considerably lower growth rates for the profile $N^{(8)}$ and $N^{(9)}$, for different profiles (not displayed here), Branch-2 can have $\sigma$ comparable to that of the Branch-1 triads.

Figure 7. Plots of non-dimensional growth rates of mode-1 triads for profile $N^{(8)}$, with $f/\omega _3 = 0.4$ and $\omega _3/N_b = 0.2$: (a) line $(n,n)$ of Branch-1; (b) line $(n+1,n)$ of Branch-1; (c) line $(n,n+1)$ of Branch-1; and (d) line $(n+4,n)$ of Branch-2.

Figure 8. Plots of non-dimensional growth rates of mode-1 triads for profile $N^{(9)}$, with $f/\omega _3 = 0.4$ and $\omega _3/N_b = 0.2$: (a) line $(n,n)$ of Branch-1; (b) line $(n+1,n)$ of Branch-1; (c) line $(n,n+1)$ of Branch-1; and (d) line $(n+6,n)$ of Branch-2.

5.1.1. Effect of variation of $f/\omega _3$ and $\omega _3/N_b$ on different branches

For both stratification profiles used in § 5.1, $f/\omega _3 = (0.3, 0.45)$ for $\omega _3/N_b = (0.2,0.7)$ is explored (hence a total of four different combinations). For $N^{(8)}$, in all four cases, the qualitative behaviour of all Branch-1 lines is similar to figure 7. However, the maximum growth rate has a significant increase from $h=-H$ to $h=-0.8H$ for $\omega _3/N_b = 0.7$. Moreover, Branch-2 triads’ maximum growth rate significantly increases in two cases of $f/\omega _3 = 0.3$ in comparison to $f/\omega _3 = 0.4$. For $N^{(9)}$, the qualitative behaviour of line $(n,n)$ is similar to what was observed in $f/\omega _3 = 0.4$. In general, the maximum growth rate is the $(n,n)$ modal combination. Interestingly, it is found that the maximum growth rate among all triads increased nearly twice from $h=-H$ to $h=-0.8H$ for $f/\omega _3 = 0.45$, $\omega _3/N_b = 0.7$. A significant increase in maximum growth rate is also found for $\omega _3/N_b = 0.2$ for the same $f$. For $f/\omega _3 = 0.3$, $\omega _3/N_b = 0.7$, the behaviour of line $(n+1,n)$ has a significant oscillation with $n$, similar to line $(n,n)$. Note that this is different from the line $(n+1,n)$ shown in figure 8. In general, it is also observed that reducing the fluid depth increases the maximum growth rate of all possible triads even without considering the $\beta$ term of the parent wave amplitude.

5.2. Variation of nonlinear coupling coefficient with domain height for self-interaction process

Here, we restrict to self-interaction of internal gravity waves that do not experience significant detuning ${\rm \Delta} \mathcal {K}$ with changes in $h$. In this subsection, we focus mainly on the superharmonic wave's nonlinear coupling coefficient $\mathcal {N}_3$ given in (2.30b).

5.2.1. Class-1 interactions

As mentioned previously in § 4.2, some Class-1 self-interactions can have negligible detuning even for a finite range of $h/H$. We study the variation of $\tilde {\mathcal {N}}_3$ under such circumstances; the different interactions considered (denoted by $\mathcal {I}_p$) are as follows:

\[\begin{equation} \mathcal{I}_1 -[P3,D2] \quad \mathcal{I}_2 -[P4,D2] \quad \mathcal{I}_3 -[P4,D3] \quad \mathcal{I}_4 -[P5,D3] \quad \mathcal{I}_5 -[P5,D4] \nonumber\\ \mathcal{I}_6 -[P6,D3] \quad \mathcal{I}_7 -[P6,D4] \quad \mathcal{I}_8 -[P6,D5] \quad \mathcal{I}_9 -[P7,D4] \quad \mathcal{I}_{10} -[P7,D5] \end{equation}\]

Here, the notation [P$m$,D$n$] denotes that the parent wave is the $m$th mode, and the daughter (superharmonic) wave is the $n$th mode.

The stratification profiles are chosen such that $N_{max}=(2 N_b, 5 N_b, 10N_b)$, $W_p=(H/200, H/100, H/50)$ and $z_c=(H/40, H/20, H/10)$. For the profiles considered, we study variations of $\tilde {\mathcal {N}}_3\equiv |\mathcal {N}_3|/\max (|\mathcal {N}_3|)$ for interactions that strictly satisfy $|{\rm \Delta} \mathcal {K}_s| < 0.01$ for $h/H \in [-1,-0.8]$. Figure 9 shows variations of $\tilde {\mathcal {N}}_3$ for two Class-1 self-interactions. Figure 9 reveals that interactions can have a non-monotonic variation of $\tilde {\mathcal {N}}_3$ with $h/H$. Moreover, the figure reveals that even relatively small changes in $h/H$ can lead to significant variations in $\tilde {\mathcal {N}}_3$. For all interactions $\mathcal {I}_p$ considered in the list above, small changes in $h/H$ can cause significant change in $\tilde {\mathcal {N}}_3$. Some generic features are summarized below. Sensitivity of $\tilde {\mathcal {N}}_3$ to $h/H$ increases as $z_c$ is increased for a given $W_p, N_{{max}}$. For every interaction, there are specific combinations of $W_p, N_{{max}}$ where this behaviour is not exhibited. Moreover, increasing $W_p$ also increases the sensitivity of $\tilde {\mathcal {N}}_3$ to changes in $h/H$. Increasing $N_{{max}}$ for given $(z_c, W_p)$ also increases the sensitivity of $\tilde {\mathcal {N}}_3$ to $h/H$.

Figure 9. Variation of $\tilde {\mathcal {N}}_3$ with $h$ for Class-1 self-interactions: (a) $\mathcal {I}_1$ and (b) $\mathcal {I}_{10}$. Each panel shows plots for two different stratification profiles (for details, see legend).

5.2.2. Class-2 interactions

Initially, we study the variation of $\mathcal {N}_3$ with $h$ for Class-2 self-interactions. To this end, we consider $\mathcal {N}_3$ of the first five modes for $27$ different stratification profiles. Similar to the case of Class-1 self-interactions in § 5.2.1, the stratification profiles are chosen such that $N_{{max}}=(2 N_b, 5 N_b, 10N_b)$, $W_p=(H/200, H/100, H/50)$ and $z_c=(H/40, H/20, H/10)$, where we consider all possible ($3^3$) combinations. Out of the $3^3=27$ combinations, nine profiles are chosen for plotting in figure 10 and thereby elucidating the effect of each individual parameter in the stratification profile. For all cases, $\omega _3 = 0.1N_b$ and $f=0$. For some higher modes, the nonlinear coupling coefficient has a band-like structure; there exists some range of $h/H$ where $\tilde {\mathcal {N}}_3$ is significantly higher in magnitude than that corresponding to other values of $h/H$. For example, the mode $n=5$ corresponding to $N_{{max}}=5N_b$ in figure 10(c) reveals a large increase in $\tilde {\mathcal {N}}_3$ near $h/H \approx -0.75$, while it is much lower at either end. Wunsch (Reference Wunsch2017) also observed such a banded structure in the self-interaction of different modes as the stratification profile was changed. The reason behind the direct analogy between our observations and that of Wunsch (Reference Wunsch2017) is straightforward – when an internal wave travels to a different domain height, essentially it travels to a different stratification profile.

Figure 10. The variation of non-dimensionalized nonlinear coupling coefficient ($\tilde {\mathcal {N}}_3$) of wave-3 (the superharmonic wave) as $h$ is varied. Altogether, there are nine blocks, each consisting of five units, and each unit represents the mode number (given by $n$). For each block, $(N_{{max}},W_p,z_c)$ is fixed. The figure is further subdivided into three horizontal panels (each panel consisting of three blocks): (a) $N_{{max}}=2N_b$, $W_p = H/200$, and $z_c$ is varied; (b) $N_{{max}}=5N_b$, $z_c = H/20$, and $W_p$ is varied; and (c) $z_c = H/10$, $W_p=H/50$, and $N_{{max}}$ is varied.

For mode-1, when $z_c, W_p, N_{{max}}$ are all on the lower side, we observe that $\tilde {\mathcal {N}}_3 \propto 1/h^4$. For higher $N_{{max}}$, even lower values of $z_c$ and $W_p$ do not have the property $\tilde {\mathcal {N}}_3 \propto 1/h^4$. In general for modes $>1$, proportionality to $1/h^4$ is lost faster as $z_c, W_p, N_{{max}}$ are increased. In several profiles, $\tilde {\mathcal {N}}_3$ of higher modes is also more sensitive to changes in $h$ than $\tilde {\mathcal {N}}_3 \propto 1/h^4$.

6. Higher-order self-interactions in the presence of a small-amplitude monochromatic topography

The focus of this section is on higher-order self-interactions between a parent wave of frequency $\omega _1$ and a superharmonic daughter wave of frequency $\omega _3=2\omega _1$ in the presence of a small-amplitude monochromatic topography. In such scenarios, the topography can act as a ‘zero frequency wave’, which can lead to resonant higher-order interactions, and is similar to the Class-2 studied in Alam et al. (Reference Alam, Liu and Yue2009) for surface gravity waves. Such higher-order self-interactions might be important for mode-1 internal waves propagating in regions where $f>\omega _d/2$ (where $\omega _d$ is the semidiurnal frequency), since triad interactions involving two (subharmonic) daughter waves is not possible. Moreover, resonant self-interaction for a mode-1 internal wave of frequency $\omega _d$ is also not possible when $f>\omega _d/2$ for any stratification profile, as a consequence of its dispersion relation (Wunsch Reference Wunsch2017). This arises from the fact that a mode-1 parent wave with frequency $\omega _d$ and a superharmonic daughter wave with frequency $2\omega _d$ fail to satisfy the horizontal wavenumber condition for a self-interaction process. Note that higher-order interactions are different from Bragg resonance focused on in Buhler & Holmes-Cerfon (Reference Buhler and Holmes-Cerfon2011), which is also a mechanism via which a parent mode-$1$ wave can decay by transferring its energy to the higher modes. In a standard Bragg resonance, resonant wave–topography interaction occurs if the bottom topography has a wavenumber $k_b$ such that $(\omega _1,k_b \pm k_1)$ satisfies the dispersion relation.

To study higher-order self-interactions, we follow the streamfunction ansatz used in Couston et al. (Reference Couston, Liang and Alam2017) for studying internal wave Bragg resonance, and in Lahaye & Llewellyn Smith (Reference Lahaye and Llewellyn Smith2020) for studying internal wave scattering due to interaction with a large-amplitude topography. This ansatz for the streamfunction of the $j$th wave is

(6.1)\begin{equation} \varPsi_{j} = \mathcal{A}_{j}(x)\,\phi_j(\eta;x)\,{\rm e}^{-{\rm i}\omega_j t} + \mathrm{c.c}. \end{equation}

The corresponding buoyancy frequency and meridional velocity are given by

(6.2)$$\begin{gather} B_{j} = \frac{{\rm i} N^{2}}{\omega_{j}}\,\frac{\partial \mathcal{A}_j}{\partial x}\,\phi_j\,{\rm e}^{-{\rm i}\omega_{j}t} + \mathrm{c.c}., \end{gather}$$
(6.3)$$\begin{gather}\mathcal{V}_{j} = \frac{{\rm i} f}{\omega_{j}}\, \frac{\mathcal{A}_j}{h}\,\frac{\partial \phi_j}{\partial \eta}\, {\rm e}^{-{\rm i}\omega_{j}t} + \mathrm{c.c}. \end{gather}$$

The above-mentioned ansatz can also be used to study systems where the detuning is ${\rm \Delta} \mathcal {K}_s \sim {O}(1)$ in the presence of a flat or slowly varying bathymetry. The functions $\phi _j$ are same as the functions used in § 2 and are obtained by solving (2.13b). Here we consider only a small-amplitude topography whose wavenumber is comparable to the parent wave, i.e. $\epsilon _h \ll {O}(1)$ and $\epsilon _k \sim {O}(1)$.

To study higher-order interactions of a parent wave propagating in the presence of a small-amplitude topography, we also consider the linear scattering of the parent wave. Note that the linear scattering of the parent wave on its own is not resonant (we do not consider topography wavenumbers that allow resonant Bragg scattering) and hence over a long distance has a negligible effect on the parent wave's amplitude. However, even the non-resonant linear interaction of the parent wave with the topography, which leads to higher modes with $\omega _1$ frequency, can impact significantly the growth of the superharmonic wave. To derive the linear scattering of the parent wave as it moves through a topography, we assume the streamfunction of the waves to be

(6.4)\begin{equation} \varPsi_{1} = \sum_{n=1}^{n=M_n} \mathcal{A}_{(1,n)}(x)\,\phi_{(1,n)}(\eta;x)\,{\rm e}^{-{\rm i}\omega_1 t} + \mathrm{c.c}., \end{equation}

where $M_n$ is the maximum mode number after which the series is truncated, and $\phi _{(1,n)}$ is the $n$th eigenfunction of the $\omega _1$ frequency. The streamfunction ansatz (6.4) is substituted in (2.5), and similar to § 2, the linear terms of (2.5) are multiplied by $\phi _{(1,n)}$ and integrated in the $\eta$ direction. This leads to $M_n$ ordinary differential equations, where the $n$th differential equation is given by

(6.5)\begin{align} & \gamma_{n}^{(3)}\left[\frac{\partial^2 \mathcal{A}_{(1,n)}}{\partial x^2} + \mathcal{K}^2_{(1,n)}\,\frac{\mathcal{A}_{(1,n)}}{h^2} \right] \nonumber\\ &\quad ={-}\sum_{m=1}^{m=M_n}\left[\frac{2\gamma_{(m,n)}^{(5)}+\gamma_{(m,n)}^{(6)}}{h^2} \left(\frac{\partial h}{\partial x}\right)^2 - \frac{2\gamma_{(m,n)}^{(7)}}{h}\, \frac{\partial h}{\partial x}\right] \mathcal{A}_{(1,m)} \nonumber\\ &\qquad -\sum_{m=1}^{m=M_n} \left[ \gamma_{(m,n)}^{(8)} - \frac{\gamma_{(m,n)}^{(5)} }{h}\left(\frac{\partial^2 h}{\partial x^2}\right)\right] \mathcal{A}_{(1,m)} \nonumber\\ &\qquad -\sum_{m=1}^{m=M_n} 2\left[\gamma_{(m,n)}^{(4)} - \frac{\gamma_{(m,n)}^{(5)}}{h}\,\frac{\partial h}{\partial x}\right] \frac{\partial \mathcal{A}_{(1,m)}}{\partial x}, \end{align}

where $\mathcal {K}_{(1,n)}$ is the corresponding eigenvalue of $\phi _{(1,n)}$. Moreover, $\gamma ^{(*)}_{(m,n)}$ are evaluated using the expressions given in Appendix A. The above equations are similar to the equations derived in Lahaye & Llewellyn Smith (Reference Lahaye and Llewellyn Smith2020), except that we do not consider waves that travel in the direction opposite to the parent wave since they are assumed to be negligible. Now that we have the full wave spectrum with $\omega _1$ frequency by solving (6.5), we model the evolution of the superharmonic wave. For simplicity, the feedback to the parent wave is neglected, which is analogous to the pump wave approximation used in § 2.

The streamfunction of the superharmonic wave $\varPsi _3$ is substituted in (2.5), and the linear terms are multiplied by $\phi _{3}$ and integrated in the $\eta$ direction. This leads to

(6.6)\begin{align} {LIN}_3 &\equiv \left[ \left(\gamma_3^{(3)}\,\frac{\partial^2 \mathcal{A}_3}{\partial x^2} \right) + \mathcal{K}^2_3\left(\frac{\mathcal{A}_3}{h^2}\,\gamma_3^{(3)} \right) \right] {\rm e}^{-{\rm i}\omega_3 t} \nonumber\\ &\quad +2\left[\frac{\gamma_3^{(5)}}{h^2}\left(\frac{\partial h}{\partial x}\right)^2 - \frac{\gamma_3^{(7)}}{h}\,\frac{\partial h}{\partial x}\right] \mathcal{A}_3 \,{\rm e}^{-{\rm i}\omega_3 t} \nonumber\\ &\quad +\left[\frac{\gamma_3^{(6)}}{h^2}\left(\frac{\partial h}{\partial x}\right)^2 - \frac{\gamma_3^{(5)} }{h}\left(\frac{\partial^2 h}{\partial x^2} \right)+ \gamma_3^{(8)} \right]\mathcal{A}_3\,{\rm e}^{-{\rm i}\omega_3 t} \nonumber\\ &\quad + 2\left(\gamma_3^{(4)} -\frac{\gamma_3^{(5)}}{h}\,\frac{\partial h}{\partial x}\right) \frac{\partial \mathcal{A}_3}{\partial x}\,{\rm e}^{-{\rm i}\omega_3 t}, \end{align}

where ${LIN}_3$ contains only linear terms, and models the propagation of superharmonic wave in the presence of a topography. Note that a superharmonic wave cannot exchange energy with higher modes of $\omega _3$. Now we move on to deriving the nonlinear terms that force the superharmonic wave. Since we are focusing on higher-order interactions, all nonlinear terms (including terms containing $x$-direction derivatives of $\phi$ and $h$), which have the same angular frequency as the superharmonic wave, are retained. In the terrain-following coordinates, this would, however, lead to a large number of terms that need to be evaluated. This issue can be circumvented following the procedure outlined below. The nonlinear terms in the terrain-following coordinates are given by the right-hand side of (2.5).

We assume that the superharmonic wave is forced nonlinearly by the $\omega _1$ spectrum. To model this, we substitute $\varPsi _{1}, B_{1}, \mathcal {V}_{1}$ into the nonlinear terms of (2.5). Note that we can obtain $B_{1}$ and $\mathcal {V}_{1}$ from (6.2) and (6.3), respectively. After the substitution, similar to the linear terms of wave-3, the nonlinear terms are multiplied by $\phi _3$ and integrated in the $\eta$-direction within the domain limits. The resultant expression obtained is

(6.7)\begin{equation} \langle {NL}_3 \rangle = \int_{{-}1}^{0}\phi_3 \left[{\rm i} \omega_3\, \mathcal{J}\{({L}_{xx}+{L}_{\eta\eta})\varPsi_1,\varPsi_1\} + {L}_{x}(\mathcal{J}\{B_1,\varPsi_1\}) - f\,{L}_{\eta} (\mathcal{J}\{\mathcal{V}_1,\varPsi_1\}) \right] {\rm d}\eta. \end{equation}

Therefore the final superharmonic wave equation can be written in a compact form:

(6.8)\begin{equation} {LIN}_{3} = \langle {NL}_3 \rangle. \end{equation}

In (6.6) and (6.8), instead of splitting $\mathcal {A}_j$ into a product of slowly varying amplitude and rapidly varying phase part, we simply solve the equations numerically by retaining $\mathcal {A}_j$ as it is. This is mainly because, as mentioned above, the number of nonlinear terms would be significantly high in terrain-following coordinates. For high ratios of $f/\omega _1$ (for example, north of the critical latitude), the parent wave cannot self-interact resonantly with the superharmonic wave in the presence of a flat bottom. However, a resonant higher-order self-interaction can occur provided that the topography has a wavenumber $k_b$ such that

(6.9)\begin{equation} k_b = k_3 - 2k_1, \end{equation}

where $k_3$ is the wavenumber of the superharmonic wave, and the $k_1$ is the wavenumber of the parent wave. In such scenarios, the daughter wave's amplitude will grow consistently. However, this being a higher-order interaction, the growth rate of the daughter wave (consequently, the decay of the parent wave) can be expected to be slower than a resonant self-interaction.

To elucidate and validate the higher-order self-interaction process, we perform numerical simulations by solving the complete 2-D Boussinesq equations and comparing the output with the results of the reduced-order model derived in this section. We run three simulations where the parent and daughter waves’ frequencies are held fixed. They are denoted by Case-1, Case-2 and Case-3. For all the simulations, the parent wave frequency is $\omega _1/N_b = 0.2$, where $\omega _1$ is the semi-diurnal frequency, i.e. $\omega _1 = 1.4 \times 10^{-4}\, {s}^{-1}$. Both the parent and daughter waves are mode-1 of their respective frequencies ($\omega _1$ and $2\omega _1$). These parameters would result in a significant detuning between the two waves at high $f$ values. The bathymetry profile is given by

(6.10)\begin{equation} h ={-}H + \left[\epsilon_h H \sin{(k_b( x-x_c))}\right]\times\frac{1}{1+(x-x_c)^{32}/W_T^{32}}, \end{equation}

where $\epsilon _h = 0.01$, $H=3000$ m, and domain length $L=540H$ and $x_c=L/2$ are held fixed across all three simulations. The stratification profile parameters, $f$ value, incoming maximum velocity of the parent mode ($u_{{max}}$), and $W_T$ for the three simulations are given in table 1.

Table 1. The stratification profile parameters, Coriolis frequency $W_T$, and velocity amplitude of the three waves for Case-1, Case-2 and Case-3.

The results after solving the reduced-order model and 2-D Boussinesq equations for the above-mentioned parameters are shown in figure 11. The 2-D Boussinesq equations are solved using Dedalus (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020). More details on the simulations are given at the end of § 7. In figures 11(ac), the amplitude of the daughter wave is observed to be slowly increasing due to the higher-order self-interaction. Moreover, the daughter wave's amplitude also rapidly oscillates because of the non-resonant standard self-interaction process between the parent wave and the daughter wave. In the absence of a varying bathymetry, only the rapid non-resonant interaction would be present without any consistent growth in the daughter wave's amplitude. Therefore, we have shown that for scenarios where Bragg resonances are not resonant, higher-order interaction might be a possible mechanism that can scatter the energy of the mode-1 internal wave.

Figure 11. Higher-order self-interaction of mode-1 internal waves in the presence of a monochromatic bathymetry for three different cases: (a) Case-1, (b) Case-2, and (c) Case-3. The topography profile (not to scale) is shown for all three cases.

7. Numerical validation

In this section, we provide numerical validations for the reduced-order equations (2.29a) and (2.29b) derived through multiple-scale analysis for two different cases. This is done by solving the 2-D Boussinesq equations in terrain-following coordinates using the open-source, pseudo-spectral code Dedalus (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020). The above-mentioned equations in primitive variables along with viscous and hyperviscous terms (the latter terms damping much smaller scales than the former) are

(7.1a)$$\begin{gather} \frac{\partial u}{\partial t} + L_{x}(P) + u\,L_{x}(u) + w\,L_{\eta}(u) = \nu\,L_{\eta\eta}(u) + \left(\frac{\nu_{6z}}{H^6}\,\frac{\partial^6 u}{\partial \eta^6}+{\nu_{6x}}\, \frac{\partial^6 u}{\partial x^6}\right), \end{gather}$$
(7.1b)$$\begin{gather}\frac{\partial w}{\partial t} + L_{\eta}(P) + u\,L_{x}(w) + w\,L_{\eta}(w) = \nu\,L_{\eta\eta}(w) + B, \end{gather}$$
(7.1c)$$\begin{gather}\frac{\partial B}{\partial t} + N^2w + u\,L_{x}(B) + w\,L_{\eta}(B) = \nu\,L_{\eta\eta}(B) + \left(\frac{\nu_{6z}}{H^6}\,\frac{\partial^6 B}{\partial \eta^6}+ {\nu_{6x}}\, \frac{\partial ^6 B}{\partial x^6}\right), \end{gather}$$
(7.1d)$$\begin{gather}L_{x}(u) + L_{\eta}(w) = 0. \end{gather}$$

Here, $(u,w)$ = $(L_{\eta }(\varPsi ),-L_{x}(\varPsi ))$, meaning that here the velocity field is defined in $x$$\eta$ instead of $x$$z$ coordinates. In all our simulations, $\nu = 10^{-5} \, {m}^{2}\, {s}^{-1}$,