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.4a–c)\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}$$

where

(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.26a–c)\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 Grimshaw1988,  Reference 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.26a–c) 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(a–c) uses the stratification profile $N^{(1)}$ with the parameters $N_{{max}} = 2N_{b}$, $W_{p} = H/200$ and $z_c = H/80$; figures 3(d–f) 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 3( f). 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(a–c) 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(a–c), 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}$, $\nu _{6x} = 10^{8}\, {m}^{6}\, {s}^{-1}$ and $\nu _{6z} = 81\, {m}^{6}\, {s}^{-1}$. Equations (2.29a) and (2.29b) are solved using the fourth-order Runge Kutta method for time stepping and a second-order-accurate discretization scheme for the term ${\partial a_j}/{\partial x}$, where the scheme is forward or backward depending on the group speed direction of the particular wave. We estimate the validity of two different cases.

For Case-1, self-interaction of a plane wave in the presence of a constant $h$ is simulated. The parameters of the simulation are $H=3000$ m and $N_b=10^{-3}\, {s}^{-1}$. The frequency of the parent wave is taken as $\omega _1/N_b = 0.447$, while $f=0$ is chosen. The stratification profile (4.1) is considered, and the parameters used are:

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

Mode-3 of the parent wave frequency ($\omega _1$) in this scenario self-interacts resonantly with mode-2 of $2\omega _1$ (Varma & Mathur Reference Varma and Mathur2017). The parent wave streamfunction input (initial condition) to the full numerical simulation is given by

(7.2)\begin{equation} \varPsi = A_1 {{\phi}_1} \sin(k_1 x) \exp(-(x-x_c)^2/W_1^2), \end{equation}

where $|A_1| = 0.06$, $W_1 \rightarrow \infty$ is chosen (note that $x_c$ can be any value since $W_1 \rightarrow \infty$), and $\varPsi$ is used to obtain $(u,w)$. The numerical code is also initialized with the corresponding buoyancy frequency for the streamfunction given in (7.2).

For estimating the energy of parent and daughter waves from the numerical simulations, only the potential energy of the waves is considered. This is valid because when $f=0$, energy is partitioned equally between potential energy and kinetic energy. For evaluating the potential energy of the two waves, we take the Fourier transform of $B$ in the $x$-direction. Then by simply isolating the $k_3$ and $k_1$ wavenumbers, the respective fields due to wave-3 and wave-1 can be obtained for all time. The resulting energy evolution of the waves for Case-1 is shown in figure 12(a). At the end of the simulation, the parent wave energy was observed to be $86\,\%$ of the total energy in the Boussinesq equations simulation, while the reduced-order model predicted that $85.5\,\%$ of the total energy will be contained in the parent wave at the specified time interval. Moreover, at the end of the simulation, the daughter wave's energies in the Boussinesq equations simulation and the reduced-order model are $13.0\,\%$ and $14.4\,\%$, respectively.

Figure 12. (a) The energy evolution of waves for Case-1 from reduced-order equations and numerical simulations. The superscript ${(N)}$ denotes the results from numerical simulation of 2-D Boussinesq equations. (b) The energy evolution of waves for Case-2 from reduced-order equations and numerical simulations in the time span $t^*=30$ to $t^*\approx 40$, where $t^* \equiv \omega _1t/2{\rm \pi}$. (c) Fourier transform of $B$ at $\eta = -0.42$ and $t^{*}=40$. Here, W1 and W3 indicate wave-1 and wave-3, respectively.

In Case-2, we consider the self-interaction of a parent wave packet travelling in the presence of a slowly varying bathymetry. The parameters considered are $H=3000$ m and $N_b=2.5\times 10^{-3}\, {s}^{-1}$. The following stratification profile is used:

1. (i) $N^{(14)}$: $N_{{max}} = 5N_b$, $W_p = 3H/100$, $z_c = H/10$.

Here, $\omega _1/N_b = 0.447$ with $f=0$ are chosen. Mode-3 of $\omega _1$ self-interacts resonantly with mode-2 of $2\omega _1$ for $h/H \in [-1,-0.8]$. The bathymetry is given by

(7.3)\begin{equation} h ={-}H + 0.1H \left[ \tanh((x-x_{t1})/W_{t1}) + \tanh((x_{t2}-x)/W_{t2})\right], \end{equation}

where $W_{t1} = 2.7H$, $W_{t2} = W_{t1}/1.3$, $x_{t1} = 25H$ and $x_{t2} = 83H$ were considered, where $100H$ is the domain length in the $x$-direction. The bathymetry shape can be visualized in figure 13.

Figure 13. Horizontal velocity plot ($u$) from Case-2 simulation at (a) $t^{*}=0$ and (b) $t^{*}=30$. The stratification profile shape used is also shown for visual purposes. A faint but clear mode-2 trail left by the mode-3 parent wave can be seen in (b). The shape of the stratification profile used in Case-2 is given by green curves in (b).

The parent wave streamfunction of the form (7.2) is used with $|A_1| = 0.022$, $W_1 = 3.57H$, and $x_c = 10.2H$ is chosen. Here, $k_1$ is evaluated at $h=-H$ for the initial conditions. The $x_c$ value is chosen such that the wave packet is just at the bottom of the ‘plateau-like’ topography at $t=0$ (as shown in figure 13a). The bathymetry is considered to be slowly varying so that the wave packet scattering by the bathymetry is negligible. Here, the same procedure of energy evaluation as for Case-1 is followed, except that the energy is evaluated from $t^{*} = 30$, when the entire energy of both wave packets is almost confined to the top of the plateau region (where $h = -0.8H$). This makes the energy evaluation straightforward. We consider $B$ only in the range $x/H \in [33, 68]$ (where most of the energy is contained), and then again perform a Fourier transform to separate the energy of the daughter and parent wave packets.

In Case-2, since wave packets are considered, the Fourier transform of $B$ would not have a sharp peak at $k_1$ and $k_3$. Instead, a smoother peak in $k$-space would be produced as shown in figure 12(c). We define a non-dimensional wavenumber $\tilde {k}$ as $\tilde {k} \equiv k/k_1$. For evaluating the energy of both wave packets, amplitude ($|a_{\tilde {k}}|$) in a finite range of $\tilde {k}$ is considered. The energy contained in $\tilde {k}\in [0.6,1.4]$ is considered as the energy of the parent wave packet, while the energy in $\tilde {k}\in [1.6,2.4]$ is considered as the daughter wave's energy. For example, the $\tilde {k}$ ranges considered for the parent and daughter waves are highlighted in figure 12(c) using coloured dashed lines for a specific $\eta$ and $t^{*}$. For the parent wave, $|a_{\tilde {k}}|$ between the blue dashed lines is considered. Similarly, for evaluating the energy of the daughter wave packet, we consider the amplitude $|a_{\tilde {k}}|$ between the red dashed lines. The energy evolution of the wave packets is shown in figure 12(b). The parent wave packet energies in the Boussinesq equations simulation and the reduced-order model were observed to be $88.1\,\%$ and $88.9\,\%$, respectively, at the end of the simulations. At the same time, the daughter wave's energies in the Boussinesq equations simulation and the reduced-order model are $8.95\,\%$ and $10.6\,\%$, respectively.

To obtain the numerical results in figure 11, (7.1a)–(7.1d) along with the equation for meridional velocity ($v$) were solved. A vertical hyperviscous term was not used, while a different horizontal hyperviscous term was used: $\nu _{12x}\,\partial ^{12}/\partial x^{12}$, where $\nu _{12x} = 1.4\times 10^{25}\, {m}^{12}\, {s}^{-1}$. The kinematic viscosity was chosen to be $\nu = 10^{-3} \, {m}^{2}\, {s}^{-1}$. The primary wave was forced by using a forcing function in the $u$-momentum equation, which sends a constant-amplitude mode-1 wave train onto the small-amplitude topography.

8. Summary and conclusion

Weakly nonlinear wave–wave interactions is one of the mechanisms through which internal gravity waves energy cascade from large length scales (hundreds of kilometres) to small scales (centimetres to metres). At small length scales, internal waves can give rise to convective or shear instabilities (Koudella & Staquet Reference Koudella and Staquet2006) and cause mixing, thus resulting in increased diffusion in oceans. The 2-D Boussinesq equations are written in terrain-following coordinates ($x$–$\eta$). Using multiple-scale analysis, we derive the amplitude evolution equations for internal gravity waves undergoing weakly nonlinear wave–wave interactions in the presence of varying density stratification (resembling that of actual oceanic scenarios) as well as mild-slope bathymetry in a vertically bounded domain. If the stratification varies with $z$ in the $x$–$z$ coordinates, then it becomes a function of both $x$ and $\eta$ in the $x$–$\eta$ coordinates when bathymetry $h$ varies with $x$. In other words, the effective stratification profile varies with the ocean depth. Both triads and self-interactions are studied, and both pure resonant conditions as well as systems with wavenumber detuning are analysed. The main results of this paper are given in a brief format in figure 14.

Figure 14. A summary diagram shows how different factors such as detuning and nonlinear coefficients can vary for different classes of interaction in the presence of non-uniform stratification. It also provides a brief picture of the higher-order self-interaction studied in § 6.

In the presence of uniform stratification, we show that the horizontal wavenumber triad condition, given by $k_{(1,a)}+k_{(2,b)}+k_{(3,c)}=0$, is not violated due to changes in $h$. Here, $a,b,c$ are the mode numbers of waves $1$, $2$ and $3$, respectively. Moreover, in the presence of uniform stratification, the nonlinear coupling coefficients are inversely proportional to the square of the fluid depth ($\propto 1/h^2$).

For non-uniform stratifications, we define two classes of interaction for both triads and self-interactions. Class-1 involves weakly nonlinear interactions of waves that do not have the same mode number. Class-2 is a special situation that involves interactions of waves with the same mode number, i.e. $a=b=c$. Class-2 triad interactions can exist only in the parameter regime $f \ll \omega _j \ll N$. Moreover, in the same parameter regime, near-resonant Class-2 self-interactions can exist with very low detuning even as $h$ is varied. This is because the wavenumbers involved in a self-interaction change in the same way as $h$ changes. For Class-1 interactions, detuning may be induced in triads and self-interactions if the waves interact in a region of varying $h$. This is because in a vertically bounded domain, the horizontal wavenumbers are not only a function of $h$ but also a function of the mode number. Moreover, the functional dependence of the wavenumber on $h$ may change as the mode number changes. Therefore, in a weakly nonlinear interaction where different mode numbers are involved, there is no constraint for the wavenumbers to satisfy the triad condition in a given range of $h$.

The variation of the growth rate of the daughter waves in both triadic- and self-interactions is studied when $h$ is varied. For both Class-1 and Class-2 self-interactions, it is observed that small changes in $h$ may result in large changes in the growth rate of the daughter waves. This characteristic is observed especially for Class-1 self-interactions. Variation of growth rates with $h$ is studied for triads of a mode-$1$ parent wave in the presence of non-uniform stratification. Triads were identified such that the daughter waves can be up to mode-$50$. For relatively small changes in $h$, the growth rates can vary significantly for triads that involve only lower modes. Moreover, the most unstable daughter wave combination for the same parent wave can also change for relatively small changes in $h$. Unlike uniform stratification, in non-uniform stratification, the growth rates do not have a monotonic behaviour with $h$. This was observed for both triadic- and self-interactions.

Reduced-order equations for higher-order self-interactions of an internal wave in the presence of a small-amplitude monochromatic topography are also derived. In the higher-order self-interaction process, the small-amplitude topography behaves as a zero-frequency wave. It is shown that such higher-order interactions can cause resonant growth of the superharmonic wave. Such higher-order interactions can play a crucial role in the decay of the mode-$1$ internal wave at latitudes greater than $28.9^{\circ }$. This is because sum-type triad interactions are not possible (Olbers, Pollmann & Eden Reference Olbers, Pollmann and Eden2020), and a mode-$1$ internal wave cannot self-interact resonantly for high values of $f$ (Wunsch Reference Wunsch2017).