Triad near-resonant instability of vertically bounded internal waves in non-uniform stratification

Abstract

Vertically bounded, horizontally propagating internal waves may become unstable through triad resonant instability, in which two sibling waves in background noise draw energy from a parent internal tide. If the background stratification is uniform, then the condition for pure resonance between the parent and sibling wave frequencies and horizontal and vertical wavenumbers can be found semi-analytically from the roots of a polynomial expression. In non-uniform stratification, determining the frequencies and horizontal wavenumbers for which resonance occurs is less straightforward. We develop a theory for near-resonant excitation of a pair of sibling waves from a low-mode internal wave in which the proximity to pure resonance is characterised by the discrepancy between the forced sibling wave frequencies and the natural frequency of these modes. Knowing this discrepancy can be used methodically to determine pure resonance conditions. This inviscid theory is compared with numerical simulations of effectively inviscid waves. For comparison with laboratory experiments, the theory is adapted to include viscous effects both in the bulk of the fluid and at the side walls of the tank. We find that our theoretical predictions for frequencies and wavenumbers of the fastest growing sibling waves are generally consistent between theory, simulations and experiments, though theory overpredicts the growth rate observed in experiments. In all cases, the growth rate of sibling waves decreases with decreasing parent wave frequency, becoming negligibly small in experiments if the parent wave has frequency less than $\approx 0.7$ of the buoyancy frequency at the surface.

1. Introduction

Internal baroclinic tides are generated within the ocean by the movement of the surface barotropic tide over submarine topography. Globally, the power going into the internal tide is approximately 1 TW (Munk & Wunsch 1998; Wunsch & Ferrari 2004). The waves are first generated as horizontally and vertically propagating beams, but transform in the far field to horizontally propagating low vertical mode internal tides (MacKinnon et al. 2017). How energy is transferred from these tides to small dissipative scales is an ongoing area of research. Numerous mechanisms have been proposed, including interaction with continental slopes and shelves, interaction with eddies and mean flows, and resonant wave–wave interactions such as parametric subharmonic instability (Alford et al. 2016).

The case of parametric subharmonic instability is a particular example of triad resonant instability (TRI), which allows for the growth from the background ambient noise of an infinitesimal pair of sibling waves that extract energy from the parent internal wave. These sibling waves are in temporal and spatial resonance with the parent, where their frequencies and wavenumbers either sum or subtract to those of the parent wave. In uniform stratification, TRI has been well studied in the cases of horizontally and vertically propagating plane waves (e.g. Lombard & Riley 1996), wave beams (e.g. Clark & Sutherland 2010; Bourget et al. 2013; Dauxois et al. 2018; Fan & Akylas 2021; Grayson, Dalziel & Lawrie 2022) and vertically confined, horizontally propagating low mode internal waves (e.g. Joubaud et al. 2012; Sutherland & Jefferson 2020). In the case of vertically bounded, internal waves in non-uniform stratification, TRI requires resonance between the frequencies and horizontal wavenumbers of the parent and sibling waves (e.g. Thorpe 1966; Davis & Acrivos 1967; Akylas & Kakoutas 2023). If there is background rotation, then resonance with inertial sibling waves is enhanced if the parent wave frequency is twice the Coriolis frequency (Young, Tsang & Balmforth 2008). With and without background rotation, it is well established that a parent internal mode in non-uniform stratification self-interacts to excite superharmonics (Sutherland 2016; Wunsch 2017; Baker & Sutherland 2020; Sutherland & Dhaliwal 2022). However, the resonant growth of a pair of sibling waves from a parent mode in non-uniform stratification without background rotation has not been as well studied. A general treatment of weakly nonlinear interactions between a parent mode and a superposition of modes having various horizontal wavenumbers and vertical structure was considered by Varma & Mathur (2017), who focused on identifying wavenumber and frequency combinations at which resonance occurs. This work was extended to look at energy transfer between three resonant modes as they propagate horizontally (Varma, Chalamalla & Mathur 2020). Neither of those studies presented explicit analytical solutions for the exponential growth rate of resonant infinitesimal amplitude sibling waves. Such growth was considered for vertically propagating waves in unbounded non-uniformly stratified fluid (Gururaj & Guha 2020), but the assumption of weakly varying stratification limited its extension to vertically bounded low mode internal tides. In a follow-up paper, Gururaj & Guha (2022) examined TRI of vertically bounded internal modes passing over slowly varying topography, with flat bottom topography being a special case. Explicit expressions for the growth rate were found in the case of uniform stratification, with their examination of non-uniform stratification focusing on the effect of topography detuning resonant waves. A rigorous treatment predicting the growth rate of sibling waves from a parent wave in a vertically bounded domain was developed by Akylas & Kakoutas (2023), who adapted the Floquet analysis of Fan & Akylas (2021) to examine resonance between vertical modes in non-uniform stratification. Numerical examination of the results, including accounting for spatial modulation of the waves, was restricted to waves in uniform stratification.

In most of the theoretical studies mentioned above, the influence of viscosity was ignored. However, viscous effects are not necessarily negligible in laboratory experiments (Davis & Acrivos 1967; Bourget et al. 2013; Gururaj & Guha 2020). With this in mind, we develop an approach different to that of Akylas & Kakoutas (2023) to derive an equation for the evolution and growth of a pair of sibling modes in near-resonance with a parent mode from which it is numerically straightforward to locate horizontal wavenumbers and vertical mode numbers of sibling waves in pure resonance and which grow exponentially. This approach is readily adapted to include the effects of viscosity in the bulk of the fluid as well as dissipation due to waves moving against the no-slip side walls of the tank. The theory is tested against numerical simulations that neglect viscosity and against laboratory experiments. The theory for near-resonant sibling excitation is derived in § 2. Experiments and their analyses are described in § 3. The numerical model is described in § 4, with a comparison between theory, experiments and simulations being provided therein. In § 5, we summarise our results, discussing the successes and limitations of theory while suggesting avenues for future research.

2. Theory

Here, we derive a theory for TRI of horizontally periodic, vertically bounded internal waves in non-uniform stratification. Our approach differs from that of Akylas & Kakoutas (2023), but we arrive at the the same result in special circumstances. Using the parent wave amplitude as a perturbation parameter, Akylas & Kakoutas (2023) ultimately derived an equation for the forcing of the vertical structure of one sibling due to the interaction between the parent and the other sibling. The equations reduced to an eigenvalue problem for the vertical structure of sibling waves having horizontal wavenumbers in pure resonance with the parent wave. Using a solvability condition, the growth rate of the sibling waves was found as it depended upon the parent wave amplitude and the degree of detuning of the sibling waves from pure resonance. Here, we derive a time evolution equation for the amplitude of sibling waves in near resonance where their forcing frequencies are shifted from their natural frequencies. We then use the result to predict the growth rate of sibling waves in pure resonance. For the sibling pair considered by Akylas & Kakoutas (2023), in which the difference of their horizontal wavenumbers equals the parent wavenumber, we find the same prediction for the growth rate in pure resonance. However, for sibling waves that are not horizontally long compared with the parent, we find that the largest growth rate occurs when the sum of the sibling horizontal wavenumbers (one being positive, the other being negative) equals the parent wavenumber. We go on to include the influence of viscosity upon the sibling wave growth rates.

The fluid is assumed to be Boussinesq, incompressible and two-dimensional, with motion in the $x$ – $z$ plane. We neglect background rotation, and to begin with, the influence of viscosity is ignored. In a stationary frame of reference, the governing equation of motion can be written as a linear differential operator $L$ acting on the streamfunction $\psi$ , being driven by nonlinear forcing $\mathcal N$ :

(2.1) \begin{equation} \underbrace { \big (\partial _{tt} (\partial _{xx} + \partial _{zz}) + N^2(z)\, \partial _{xx}\big )}_{\equiv \,L} \psi = \underbrace { \boldsymbol{\nabla }\boldsymbol{\cdot }\left [ \partial _t (\boldsymbol{u} \zeta ) - \partial _x (\boldsymbol{u} b) \right ]}_{\equiv {\mathcal \, N}}, \end{equation}

in which $\boldsymbol{\nabla }=(\partial _x,\partial _z)$ , $\boldsymbol{u}=(u,w)=(-\partial _z\psi ,\partial _x\psi )$ is the velocity, $\zeta =\partial _z u - \partial _x w$ is the spanwise vorticity, and $b$ is the buoyancy. Here, $N^2(z)$ is the background buoyancy frequency, given by $N^2(z) = -({g}/{\rho _0})({\rm d}\bar {\rho }/{{\rm d}z})$ , in which $\bar {\rho }(z)$ is the background density profile, $\rho _0$ is the characteristic density, and $g$ is gravity. These equations are solved in an infinite horizontal domain with free-slip and no normal flow upper and lower boundary conditions at $z=0$ and $z=-H$ , respectively.

2.1. Mode structure and dispersion relation

We seek the vertical structure and dispersion relation of horizontally periodic, vertically bounded, small-amplitude internal modes.

For simplicity, we begin by examining the properties of the vertical mode-1 ( $j_0=1$ ) parent wave. This wave has prescribed horizontal wavenumber $k_0$ and streamfunction amplitude $A_0$ . The streamfunction of the parent wave is thus given by

(2.2) \begin{equation} \psi _0 = \dfrac {1}{2} A_0\, \hat {\psi }_0(z)\, {\rm e}^{{\rm i}\phi _0} + \text{c.c.}, \quad \phi _0 = k_0 x - \omega _0 t, \end{equation}

in which c.c. denotes the complex conjugate. Putting this expression in the linear differential equation $L\psi =0$ gives the eigenvalue problem

(2.3) \begin{equation} \hat {\psi }_0^{\prime \prime } = -k_0^2 \left (\frac {N^2(z)}{\omega _0^2}-1\right ) \hat {\psi }_0, \end{equation}

in which primes on the left-hand side denote $z$ -derivatives. The no normal flow boundary conditions at the top and bottom of the domain require $\hat {\psi }_0(-H)=\hat {\psi }_0(0)=0$ . The boundary value problem is solved through a Galerkin method. For given horizontal wavenumber $k_0$ , we extract the frequency $\omega _0$ of the vertical mode-1 wave and its vertical structure $\hat {\psi }_0(z)$ , which is positive everywhere within the domain, having maximum value 1. From this solution, we can compute the velocity field, the spanwise vorticity and the buoyancy, as listed in table 1. In particular, the expression for $\zeta _0$ makes use of (2.3), and the expression for $b_0$ follows from the linearised buoyancy equation $\partial _t b_0 = -N^2(z)\, w_0$ .

Table 1. Polarisation relations for the parent wave (left-hand column) and sibling waves (right-hand column). In these expressions, $c_0\equiv \omega _0/k_0$ , $c_\pm \equiv \omega _\pm /k_\pm$ , $C_\pm \equiv \varOmega _\pm /k_\pm$ , $\phi _0 = k_0\,x -\omega _0\,t$ and $\phi _\pm = k_\pm \,x -\varOmega _\pm \,t$ .

We seek how a pair of sibling waves interact with the parent mode resulting in their possible near-resonant growth. The waves have horizontal wavenumbers $k_+$ and $k_-$ satisfying the resonance condition $k_+ \pm k_- =k_0$ . The vertical mode numbers of the sibling waves are denoted by $j_+$ and $j_-$ , such that $j_\pm -1$ represents the number of nodes of the vertical structure of the streamfunction of the waves within the interior of the domain. Explicitly, the vertical structure is given by $\hat {\psi }_\pm ^{\prime \prime } = -k_\pm ^2 ( ({N^2(z)}/{\omega _\pm ^2})-1 ) \hat {\psi }_\pm$ , with $\hat {\psi }_\pm (-H)=\hat {\psi }_\pm (0)=0$ , in which it is understood that $\hat {\psi }_\pm$ depends upon $j_\pm$ as well as $k_\pm$ , with corresponding dispersion relation $\omega _\pm = \omega _\pm (k_\pm , j_\pm )$ . In uniform stratification, the vertical modes must satisfy the additional resonance condition $j_+ \pm j_- = j_0$ . But this is not necessarily the case in non-uniform stratification (Akylas & Kakoutas 2023).

The modes have natural frequencies $\omega _+$ and $\omega _-$ . However, the interaction of the parent with one of the siblings does not necessarily excite the other sibling at its natural frequency. For this reason, we denote the near-resonant forcing frequencies of the sibling waves by $\varOmega _+$ and $\varOmega _-$ , which satisfy $\varOmega _+ + \varOmega _- = \omega _0$ . Explicit expressions for $\varOmega _\pm$ in terms of $\omega _\pm$ are determined in § 2.3. The streamfunction of the sibling waves is thus given by

(2.4) \begin{equation} \psi _\pm = \dfrac {1}{2} A_\pm\, \hat {\psi }_\pm (z)\, {\rm e}^{{\rm i}\phi _\pm } + \text{c.c.}, \quad \phi _\pm = k_\pm x - \varOmega _\pm t. \end{equation}

Here, the streamfunction amplitudes $A_\pm = A_\pm (T)$ are functions of a slow time variable $T$ (defined in § 2.3). The amplitudes are anticipated to grow exponentially should the sibling waves be in near resonance with the parent wave, meaning that the forcing frequency of each sibling wave is close to the natural frequency of each mode: $\varOmega _\pm \approx \omega _\pm$ . In pure resonance, $\varOmega _\pm = \omega _\pm$ . The corresponding polarisation relations for these waves are listed in table 1.

2.2. Weakly nonlinear forcing

For a sufficiently small-amplitude parent wave, the dominant interaction is between the parent and a pair of sibling waves (Akylas & Kakoutas 2023), which presumably grow out of background noise. Using the ‘pump-wave’ approximation (Craik & Adam 1978; Young et al. 2008; Gururaj & Guha 2020), we assume that the influence of the sibling waves upon the parent itself is negligible, so $A_0$ remains constant. We begin by evaluating the nonlinear forcing $\mathcal N$ defined in (2.1), resulting from the parent interacting with the ‘ $+$ ’ sibling to drive the ‘ $-$ ’ sibling. Resonance between their horizontal wavenumbers is required so that $k_+\pm k_-=k_0$ . We focus first on the $k_+ + k_- = k_0$ resonance, for which the quadratic interaction between the part of the parent with phase $\phi _0$ and the part of the $+$ sibling with complex conjugate phase $-\phi _+$ results in forcing the $-$ sibling with phase $\phi _0-\phi _+ = [(k_0-k_+) x - (\omega _0-\varOmega _+)t] = k_-\, x-\varOmega _-\, t= \phi _-$ . The corresponding four products of fields appearing in $\mathcal N$ are listed in table 2. Hereafter, if not written explicitly, it is understood that $N^2$ is generally a function of $z$ .

Table 2. Products (left-hand column) appearing in the nonlinear forcing between the parent and $+$ sibling, and resulting expressions (right-hand column) that multiply the phase ${\rm e}^{{\rm i}\phi _-}$ . The star superscript denotes the complex conjugate.

Combining these results in the nonlinear forcing gives the following expression for the driving of the part of the $-$ sibling with phase $\phi _-$ :

(2.5) \begin{align} {\mathcal N}_{-}^+ &= \frac {1}{4}A_0 A_+^\star \left [ \left (\varOmega _- k_+ \left (\frac {1}{c_+^2} -\frac {1}{c_0^2}\right ) + k_- k_+ \left (\frac {1}{C_+}-\frac {1}{c_0}\right )\right ) N^2 \hat {\psi }_0^\prime \hat {\psi }_+ \right. \notag\\ & \left. \quad {}+ \left (\varOmega _- k_0 \left (\frac {1}{c_+^2} -\frac {1}{c_0^2}\right ) + k_- k_0 \left (\frac {1}{C_+}-\frac {1}{c_0}\right )\right ) N^2 \hat {\psi }_+^\prime \hat {\psi }_0 \right. \notag\\ & \left.\quad {}+ \left (\varOmega _- \left (\frac {k_0}{c_+^2} -\frac {k_+}{c_0^2}\right ) +k_- \left (\frac {k_0}{C_+} -\frac {k_+}{c_0}\right )\right ) \frac {{\rm d}N^2}{{\rm d}z} \hat {\psi }_0 \hat {\psi }_+ \right] {\rm e}^{{\rm i}\phi _-}. \end{align}

In these expressions, we have defined the horizontal phase speeds $c_0=\omega _0/k_0$ , $c_\pm =\omega _\pm /k_\pm$ and $C_\pm =\varOmega _\pm /k_\pm$ . The nonlinear forcing of the $+$ sibling, ${\mathcal N}_{+}^+$ , is given by the above expression after exchanging $+$ and $-$ signs in the subscripts. In both cases, the superscript $+$ notation on $\mathcal N$ indicates that this forcing results from the interaction involving the sum of sibling horizontal wavenumbers: $k_+ + k_-=k_0$ .

For the resonant interaction with $k_+-k_-=k_0$ , the forcing of the $-$ sibling results from nonlinear products that are the complex conjugates of the values appearing in table 2. For the forcing of the $+$ sibling, the nonlinear products do not involve complex conjugates. It is thus straightforward to derive the nonlinear forcing ${\mathcal N}_{\pm }^-$ by analogy with (2.5).

2.3. Linear response

We now examine the response to the nonlinear forcing, focusing upon deriving a time evolution equation for the streamfunction amplitude of the $-$ sibling, $A_-$ . The time evolution is assumed to be slow compared to the inverse frequency $\varOmega _-^{-1}$ . Thus we define a slow variable $T=\epsilon t$ . Applying the linear operator $L$ to $\psi _-$ , and neglecting terms of $O(\epsilon ^2)$ , gives

(2.6) \begin{equation} L\psi _- \approx {\rm i}\frac {\varOmega _-}{c_-^2} \left [\epsilon \frac {{\rm d}A_-}{{\rm d}T} -\frac {1}{2}{\rm i}\omega _-\frac {\omega _-}{\varOmega _-} \left (\frac {\varOmega _-^2}{\omega _-^2}-1\right ) A_- \right ]\,N^2 \hat {\psi }_-\, {\rm e}^{{\rm i}\phi _-} + \text{c.c.}. \end{equation}

The corresponding expression for the $+$ sibling has the $-$ subscripts replaced with $+$ subscripts. At near resonance, we expect $\varOmega _-\approx \omega _-$ and $\varOmega _+\approx \omega _+$ . This inspires us to define the small parameter $\epsilon$ such that

(2.7) \begin{equation} 2\epsilon = \frac {\varOmega _-^2}{\omega _-^2}-1 = \frac {\varOmega _+^2}{\omega _+^2}-1. \end{equation}

Given $\varOmega _+ + \varOmega _- = \omega _0$ , we arrive at the following explicit expressions for the near-resonant frequencies and $\epsilon$ :

(2.8) \begin{equation} \varOmega _+ = \frac {\omega _0}{1 + \omega _-/\omega _+},\quad \varOmega _- = \frac {\omega _0}{1 + \omega _+/\omega _-},\quad \epsilon = \frac {1}{2}\left (\frac {\omega _0^2}{(\omega _+ + \omega _-)^2} -1\right)\!. \end{equation}

Should the dispersion relation allow pure resonance to occur ( $\varOmega _+=\omega _+$ , $\varOmega _-=\omega _-$ , $\omega _+ + \omega _-=\omega _0$ ) for some $k_\pm$ and $j_\pm$ , then $\epsilon = 0$ .

2.4. Evolution equations

For given horizontal wavenumber, the vertical modes $\hat {\psi }_\pm (z)$ are orthogonal with respect to weight $N^2$ . Thus the ordinary differential equation for the time evolution of the sibling waves can be found explicitly. For example, the evolution of the $-$ sibling in the $k_++k_-=k_0$ interaction can be found by multiplying both sides of $L\psi _- = {\mathcal N}_{-}^+$ by $\hat {\psi }_-$ and integrating in $z$ over the vertical domain.

Following this procedure for the $+$ and $-$ sibling evolution equations leads to the equations governing the time evolution of the sibling wave amplitudes. For the triad interaction $k_+ + k_-=k_0$ , the evolution equations are

(2.9) \begin{equation} \frac {{\rm d}A_\pm }{{\rm d}T} - {\rm i} \omega _\pm \frac {\omega _\pm }{\varOmega _\pm }A_\pm = -{\rm i} \frac {1}{\epsilon } M_{\pm }^+ A_0 A_\mp ^\star . \end{equation}

For the triad interaction $k_+ - k_- = k_0$ , the evolution equations are

(2.10) \begin{equation} \frac {{\rm d}A_+}{{\rm d}T} - {\rm i}\omega _+ \frac {\omega _+}{\varOmega _+}A_+ = - {\rm i}\frac {1}{\epsilon } M_{+}^- A_0 A_-, \quad \frac {{\rm d}A_-}{{\rm d}T} -{\rm i} \omega _- \frac {\omega _-}{\varOmega _-}A_- = - {\rm i}\frac {1}{\epsilon } M_{-}^- A_0^\star A_+. \end{equation}

In these equations, the star denotes the complex conjugate, and

(2.11) \begin{align} M_{\pm }^s &\equiv \frac {1}{4} \frac {c_\pm ^2}{\varOmega _\pm } \left[ \left (\varOmega _\pm k_\mp \left (\frac {1}{c_\mp ^2} -\frac {1}{c_0^2}\right ) + k_- k_+ \left (\frac {1}{C_\mp }-\frac {1}{c_0}\right )\right ) I_1 \right. \notag\\ & \left. \quad {} + S_\pm (s)\,k_0 \left (\varOmega _\pm \left (\frac {1}{c_\mp ^2} -\frac {1}{c_0^2}\right ) + k_\pm \left (\frac {1}{C_\mp }-\frac {1}{c_0}\right )\right ) I_{2\pm } \right. \notag\\ & \left.\quad {} + \left (\varOmega _\pm \left (\frac {S_\pm (s)\, k_0}{c_\mp ^2} -\frac {k_\mp }{c_0^2}\right ) +k_\pm \left (\frac {S_\pm (s) k_0}{C_\mp } -\frac {k_\mp }{c_0}\right )\right ) I_3 \right ]/I_{0\pm }, \end{align}

in which $s=\pm$ , corresponding to the triad interaction $k_+ + s k_- = k_0$ . The sign $S_\pm (s)$ is defined so that $S_\pm (+)=1$ and $S_\pm (-)=\pm 1$ : the sign in front of $k_0$ is positive in $M_\pm ^+$ and $M_+^-$ , but is negative in $M_-^-$ . The integrals over the domain depth are

(2.12) \begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad I_{0\pm } = \int _{-H}^0 N^2\, \hat {\psi }_\pm ^2 \,{\rm d}z, \end{align}

(2.13) \begin{align}& I_1 = \int _{-H}^0 N^2\, \hat {\psi }_0^\prime \hat {\psi }_+ \hat {\psi }_-\,{\rm d}z, \!\!\quad I_{2\pm } = \int _{-H}^0 N^2\, \hat {\psi }_0 \hat {\psi }_\mp ^\prime \hat {\psi }_\pm \,{\rm d}z, \!\!\quad I_3 = \int _{-H}^0 \frac {dN^2}{{\rm d}z}\, \hat {\psi }_0 \hat {\psi }_+ \hat {\psi }_-\,{\rm d}z. \end{align}

2.5. Instability

For the $k_+ + k_- = k_0$ triad interaction, we combine the equations in (2.9) to derive a single equation in one of the amplitudes $A_+$ or $A_-$ . In particular, eliminating $A_+^\star$ from the equations for $A_-$ and $A_+^\star$ gives

(2.14) \begin{equation} \frac {{\rm d}^2A_-}{{\rm d}T^2} + {\rm i} \left (\omega _+\frac {\omega _+}{\varOmega _+} - \omega _-\frac {\omega _-}{\varOmega _-} \right ) \frac {{\rm d}A_-}{{\rm d}T} + \left (\omega _+\omega _-\frac {\omega _+\omega _-}{\varOmega _+\varOmega _-} - \frac {1}{\epsilon ^2} M_+^+ M_-^+\, |A_0|^2 \right ) A_- = 0. \end{equation}

From this, and its counterpart for $A_+$ , we can derive an expression for the growth rate (and frequency) of the sibling waves by assuming $A_\pm \propto {\rm e}^{\varSigma _\pm T} = {\rm e}^{\sigma _\pm t}$ , with $\sigma = \epsilon \varSigma$ , requiring $|\epsilon |\ll 1$ . Making use of (2.8), we find

(2.15) \begin{equation} \sigma _\pm ^+ = \frac {1}{2}\omega _0\left [ \mp {\rm i} \frac {\epsilon }{\sqrt {1+2\epsilon }}\frac {\omega _+ -\omega _-}{\omega _0} +\left (4 \frac {M_+^+M_-^+}{\omega _0^2}\, |A_0|^2 - \frac {\epsilon ^2}{(1+2\epsilon )^2}\right )^{1/2}\right]\!, \end{equation}

in which we have ignored the negative square root solution. The $+$ superscript on $\sigma$ indicates that this prediction arises from the $k_++k_-=k_0$ interactions. The real part of this expression (which is the same for both sibling waves) gives their exponential growth rate, provided that the expression in the discriminant is positive. In particular, in pure resonance ( $\epsilon =0$ ), $\sigma ^+=\sigma _\pm ^+$ is real and the growth rate is largest:

(2.16) \begin{equation} \sigma ^+ \approx \sqrt {M_+^+M_-^+}\, |A_0|. \end{equation}

For the $k_+ - k_- = k_0$ triad interaction, the sibling wave amplitude evolution equations are

(2.17) \begin{equation} \frac {{\rm d}^2A_\pm }{{\rm d}T^2} - {\rm i}\left (\omega _+\frac {\omega _+}{\varOmega _+} + \omega _-\frac {\omega _-}{\varOmega _-} \right ) \frac {{\rm d}A_\pm }{{\rm d}T} - \left (\frac {1}{\epsilon ^2} M_+^- M_-^-\, |A_0|^2 - \omega _+\omega _-\frac {\omega _+\omega _-}{\varOmega _+\varOmega _-} \right ) A_\pm = 0. \end{equation}

As before, seeking solutions with $A_\pm \propto {\rm e}^{\epsilon \sigma _\pm t}$ gives

(2.18) \begin{equation} \sigma _\pm ^- = \frac {1}{2}\omega _0\left [ \mp {\rm i} \frac {\epsilon }{1+2\epsilon} +\left (-4 \frac {M_+^-M_-^-}{\omega _0^2}\, |A_0|^2\ -\ \frac {\epsilon ^2}{1+2\epsilon } \left (\frac {\omega _+-\omega _-}{\omega _0}\right )^2\right )^{1/2}\right ]\!. \end{equation}

Positive growth near pure resonance occurs in this case if $M_+^-M_-^-\lt 0$ , with largest growth at resonance, for which

(2.19) \begin{equation} \sigma ^- \approx \sqrt {-M_+^-M_-^-}\, |A_0|. \end{equation}

This final result was also found by Akylas & Kakoutas (2023).

Our focus is upon evaluating sibling horizontal wavenumbers and vertical mode numbers at which pure resonance occurs, and then finding the corresponding growth rates, $\sigma ^\pm$ , given by (2.16) and (2.19). However, we mention here that (2.15) and (2.18) also place bounds on the degree of frequency detuning from pure resonance that leads to sibling wave growth. For the $k_++k_-=k_0$ interaction, we require $|\epsilon |\lesssim 2\sigma ^+/\omega _0$ ; for the $k_+-k_-=k_0$ interaction, we require $|\epsilon |\lesssim 2\sigma ^-/|\omega _+-\omega _-|$ . These ranges are smaller if the parent wave amplitude, and hence maximum growth rate, is smaller. Akylas & Kakoutas (2023) likewise found amplitude-dependent bounds on the instability of sibling waves, but cast in terms of detuning by $\delta _k$ of the horizontal wavenumbers from resonance: $k_+-k_-=k_0+\delta _k$ .

2.6. Inclusion of viscosity

As in Davis & Acrivos (1967), the results above are readily adapted to include the influence of viscosity upon the growth of sibling waves that are in resonance with the parent wave in the absence of viscosity. We assume that diffusivity is so small compared to kinematic viscosity that its influence can be ignored. We further assume that the Reynolds number based on the parent wave properties is sufficiently large that attenuation of the parent wave is negligible on the time scale for growth of the sibling waves.

Under the above assumptions, we may assume that the vertical structure of sibling waves is negligibly affected by viscosity. Hence the linear expression for the growth of the sibling wave with horizontal wavenumber $k_-$ , given by (2.6), is modified by the addition of the viscous term $(1/2)i\nu \varOmega _- (N^2/c_-^2)^2 A_- + \text{c.c.}$ The $-$ subscripts are replaced with $+$ subscripts for the $k_+$ sibling. Consequently, the expressions on the left-hand sides of (2.9) and (2.10) are modified by the addition of the term $(\omega _0 \delta _\pm /\epsilon ) A_\pm$ , in which

(2.20) \begin{equation} \delta _\pm \equiv \frac {\nu }{2\omega _0}\frac {1}{c_\pm ^2} \frac {I_{4\pm }}{I_{0\pm }}, \quad \mbox{with}\quad I_{4\pm }\equiv \int _{-H}^{0} N^4\hat {\psi }_{\pm }^2\,{\rm d}z.\end{equation}

At resonance ( $\epsilon =0$ ), the sibling wave growth rates are thus given by

(2.21) \begin{equation} \sigma ^\pm \approx -\omega _0 (\delta _+\pm \delta _-)+ \sqrt {\pm M_+^\pm M_-^\pm \,|A_0|^2 + \omega _0^2(\delta _+\mp \delta _-)^2},\quad \mbox{for}\ k_+\pm k_-=k_0, \end{equation}

provided that the discriminant in these expressions is positive. For the $k_+ + k_-=k_0$ interaction, the predicted growth rate $\sigma ^+$ reduces to the prediction of Joubaud et al. (2012) in the case of uniform stratification.

For application to laboratory experiments, we also consider the influence of dissipation occurring at the side walls of the tank as a consequence of the no-slip surfaces dampening oscillatory motion. Damping from flow over the bottom boundary is considered negligible upon waves in the bulk of the domain due to the stratification. Comparing the rate of energy loss by the side walls with the energy of a wave with frequency $\omega$ , spanning a tank of width $W$ , gives the damping time scale $W/(\omega d)$ , in which $d\equiv (2\nu /\omega )^{1/2}$ is the well-known boundary layer thickness associated with motion over an oscillating flat plate. Hence the rate for amplitude decay due to dissipation at both side walls is given by $\omega _0 \unicode{x1D6E5} _\pm$ , in which

(2.22) \begin{equation} \unicode{x1D6E5} _\pm \equiv \frac {\sqrt {\nu \omega _\pm /2}}{W\omega _0}. \end{equation}

A similar expression for this damping rate was derived in the context of deep-water waves (Sauret et al. 2015). The growth rate of resonant waves influenced by bulk and side-wall viscous effects is given by (2.21) with $\delta _\pm$ replaced by $\delta _\pm +\unicode{x1D6E5} _\pm$ .

Although we have computed the influence of viscous damping rates upon the growth of sibling waves, the same procedure can be applied to assess the damping of the parent mode itself. The damping rate due to bulk viscosity is $\omega _0 \delta _0$ , in which $\delta _0 \equiv [\nu /(2 \omega _0 c_0^2)] (\int _{-H}^{0} N^4\hat {\psi }_{0}^2\,{\rm d}z)/(\int _{-H}^{0} N^2\hat {\psi }_{0}^2\,{\rm d}z)$ , and the damping rate due to side-wall dissipation is $\omega _0\unicode{x1D6E5} _0$ , in which $\unicode{x1D6E5} _0 \equiv \sqrt {\nu /(2\omega _0)}/W$ . For the parameters of our laboratory experiments (see § 3), we typically find that $\unicode{x1D6E5} _0$ is an order of magnitude larger than $\delta _0$ , indicating that side-wall dissipation is the most significant viscous effect. The value of $\unicode{x1D6E5} _0$ is less than $1$ %, helping to justify our assumption in theory that the parent wave amplitude can be treated as constant, though we will show that viscous attenuation of the parent wave in experiments non-negligibly influences the sibling wave growth rates.

2.7. Theory results

For a given parent wave, we perform a scan over horizontal wavenumbers of the sibling waves, $k_\pm$ , satisfying $k_+ + k_- = k_0$ and $k_+ - k_-=k_0$ . For each $k_\pm$ , we scan over a range of vertical mode numbers $j_\pm$ of both sibling waves, computing the growth rate $\sigma$ in near-resonant cases satisfying $|\epsilon |\leqslant 0.001$ . If the growth rate is positive, then a convergence routine is used to find nearby values of $k_\pm$ for which $|\epsilon |\lt 10^{-6}$ . This combination of $k_\pm$ and corresponding $j_\pm$ is taken to correspond to pure resonant instabilities.

As an example, here we present the results of these scans for a parent wave produced in one of the laboratory experiments (discussed in § 3). The background stratification has $N$ decreasing linearly with depth below the surface according to

(2.23) \begin{equation} N = N_0 + \unicode{x1D6E5} _N z/H,\quad -H\leqslant z\leqslant 0, \end{equation}

in which $N_0=1.04\,\textrm {s}^{-1}$ , $\unicode{x1D6E5} _N=0.59\,\textrm {s}^{-1}$ and $H=32.5$ cm.

Figure 1. For $w_0=2.0\,\mbox{cm}\, \rm {s^{-1}}$ , $\omega _0=0.72\,\textrm {s}^{-1}$ ( $k_0=0.20\,\textrm {cm}^{-1}$ ), $N_0=1.04\,\textrm {s}^{-1}$ , $\unicode{x1D6E5} _N=0.59\,\textrm {s}^{-1}$ : (a) predicted growth rates $\sigma$ for all resonant mode-number pairs computed in inviscid fluid for $0.5\leqslant k_+/k_0 \leqslant 4.5$ ( $k_- = k_0-k_+$ as closed red circles; $k_- = k_+-k_0$ as crosses), and (b) the corresponding resonant vertical mode numbers for $k_- = k_0-k_+$ , where red (blue) squares give the mode number of $j_+$ ( $j_-$ ). Growth rates predicted in theory (c) including bulk viscosity ( $k_- = k_0-k_+$ as open blue circles; $k_- = k_+-k_0$ as crosses), and (d) also including side-wall damping ( $k_- = k_0-k_+$ as closed blue circles; $k_- = k_+-k_0$ as crosses).

We begin by considering a parent wave with horizontal wavenumber $k_0=0.20\,\textrm {cm}^{-1}$ , which has corresponding frequency $\omega _0=0.72\,\textrm {s}^{-1}$ . The maximum vertical velocity of the wave is taken to be $w_0=2.0\,\mbox{cm}\, \rm {s^{-1}}$ , consistent with the wave amplitudes produced in many of the experiments. The corresponding streamfunction amplitude is $A_0=10.0\,\mbox{cm}^2\, \rm {s^{-1}}$ . Neglecting viscous effects, the resonant growth rates associated with different pairs of sibling waves interacting with the parent are plotted in figure 1(a). The largest growth rate across all $k_+$ occurs for $k_+=1.1k_0$ and $k_-=0.1k_0$ . However, there is insignificant growth of resonant instabilities associated with the $k_+ - k_- = k_0$ interaction if $k_+\gtrsim 1.3k_0$ ( $k_-\gtrsim 0.3k_0$ ). It is generally the case for all parent wave frequencies and stratifications examined that the growth rates associated with the $k_+-k_-=k_0$ interactions is negligibly small compared with the $k_++k_-=k_0$ interactions for $k_+\gtrsim 1.5k_0$ . In laboratory experiments, we do not see growth of sibling waves with horizontal wavelengths larger than two wavelengths of the parent wave ( $|k_\pm |\lt 0.5k_0$ ), possibly because the parent wave amplitude attenuates with distance from the wavemaker (see § 3.2). For this reason, in most of what follows, we focus our scan for resonant sibling waves interacting through $k_+ + k_-=k_0$ over the range $1.5\leqslant k_+/k_0\leqslant 4.5$ ( $-0.5\geqslant k_-/k_0\geqslant -3.5$ ). The upper bound of this range is chosen because the associated vertical mode numbers of resonant sibling waves exceeds 16, which is much larger than what is observed in experiments. For example, these mode numbers for each unstable resonant sibling pair are plotted in figure 1(b). Typically, as $|k_\pm |$ increases, so do $j_\pm$ , with the largest growth rates occurring if $|j_+ -j_-|=1$ . By neglecting viscous effects, the maximum resonant growth rate is nearly constant as $|k_\pm |$ becomes large, though the largest value of $\sigma =0.049$ occurs at $k_+=2.56k_0$ . This result is consistent with the predictions for doubly periodic inviscid waves in uniform stratification, for which the growth rate is nearly constant as the sibling wave horizontal wavenumber becomes larger (Bourget et al. 2013; Sutherland & Jefferson 2020). As with periodic waves (Bourget et al. 2013), including the effect of Laplacian viscosity within the bulk of the fluid decreases the maximum growth rate of resonant waves according to (2.21). The largest growth rate for the $k_++k_-=k_0$ interaction occurs at smaller $|k_\pm |$ . For example, as shown in figure 1(c), the fastest growth rate for resonant sibling waves affected by viscosity in the bulk of the fluid is $\sigma =0.048\,\textrm {s}^{-1}$ , occurring for $k_+=1.53k_0$ . The growth rate is further reduced to $\sigma =0.043\,\textrm {s}^{-1}$ if the effects of side-wall dissipation are included (figure 1 d).

We go on to examine the maximum growth rate $\sigma$ of sibling waves across all pure resonant sibling (horizontal) wavenumbers $k_\pm$ and vertical mode numbers $j_\pm$ for a parent wave with frequency in a range between $\omega _0=0.4\,\textrm {s}^{-1}$ and $\omega_0=0.9\,\textrm {s}^{-1}$ . The result is plotted in figure 2. This shows that sibling waves grow faster if the frequency of the parent wave is larger. For the choice of parent wave amplitude having $w_0=1\,\mbox{cm}\, \rm {s^{-1}}$ , bulk viscosity moderately decreases the maximum growth rate if the parent wave frequency is large. Also including side-wall dissipation reduces the growth rate over a wide range of parent wave frequencies such that there is no instability predicted for parent waves with frequency below $0.5\,\textrm {s}^{-1}$ . The influence of viscosity is more pronounced if the parent wave amplitude is smaller, with the largest growth rate dropping for $\omega _0\gtrsim 0.85\,\mbox{s}^{-1}$ if $w_0=0.25\,\mbox{cm}\, \rm {s^{-1}}$ .

Figure 2. Maximum growth rate of sibling waves across all horizontal and vertical wavenumbers as it depends on the frequency of the parent wave, with $N_0 = 1.04\, \textrm {s}^{-1}$ and $\unicode{x1D6E5} _N = 0.59\, \textrm {s}^{-1}$ . The predictions are shown neglecting viscosity (red circles), including viscosity in the bulk (black squares), and also including side-wall dissipation (blue triangles). Solid (open) symbols are plotted for cases where the parent wave has fixed vertical velocity amplitude $2\,\mbox{cm}\, \rm {s^{-1}}$ ( $0.5\,\mbox{cm}\, \rm {s^{-1}}$ ).

3. Laboratory experiments

3.1. Experiment set-up and analyses

Experiments were performed in a $4$ m long, $0.17$ m wide, $0.4$ m deep acrylic tank, composed by joining five components. A schematic front view of the tank can be seen in figure 3(a). The tank was filled with a non-uniform salt stratified fluid following the procedure described in Varma et al. (2024). This procedure is an adaptation of the standard double bucket method (Oster 1965; Hill 2002; Economidou & Hunt 2009). The desired density profile was obtained by imposing time-dependent flow rates using two computer-controlled pumps. A typical profile, measured using a conductivity-temperature probe mounted on a vertical linear traverse, is plotted in figure 3(b). Except for a shallow surface mixed region, this profile is well described by a linearly decreasing buoyancy frequency profile with depth, given by (2.23), with quadratic change in $N^2$ (figure 3 c). Experiments were performed just after filling the tank so that the mixed layers at the top and the bottom of the fluid did not have enough time to be well developed. A horizontally propagating vertical mode parent wave was generated at the left-hand end of the tank using a wavemaker, previously used by Dossmann et al. (2017), Husseini et al. (2020) and Varma et al. (2024). This wavemaker is made of 50 independently controlled plates, driven by 50 different motors. This wavemaker enabled us to control the frequency and amplitude of a vertical mode-1 internal wave. For each experiment, the vertical profile of the wavemaker was chosen to correspond to the predicted vertical structure $\hat {\psi }_0(z)$ corresponding to the measured profile $N(z)$ (see (2.3)). The forcing amplitude was set so that in most experiments, the maximum horizontal displacement of the top plate was $1$ cm, though in some experiments this was reduced to $0.75$ cm and $0.5$ cm. The tank was chosen to be long enough so that experiments could be completed before the return of the parent wave that reflected at the right-hand end of the tank entered the observation window: in typical experiments, the horizontal group velocity of the parent mode was of the order of $c_g\approx 1\,\mbox{cm}\, \rm {s^{-1}}$ , so that the return time of the parent mode was of the order of $10$ minutes, which is much longer than the observed time for onset of TRI when it occurred.

Figure 3. (a) Sketch of the experimental set-up. The depth of non-uniformly stratified fluid is $32.5$ cm. The pattern of dots shown in the background is used to visualise and measure of internal waves using synthetic schlieren. The black rectangle represents a joint between two sections of the tank, and is a region where visualisation is not possible. The internal wave modes are forced by imposing the horizontal displacement at the left boundary. (b) Background density profile measure for Exp. II (in $\textrm {g}\, \rm {cm}^{-3}$ ) measured using a conductivity probe. The red line is a fit of the data points using (2.23). (c) Corresponding profiles of the buoyancy frequency $N(z)$ (red line, $\textrm {s}^{-1}$ ) and the squared buoyancy frequency $N^2(z)$ (blue dashed line, $\textrm {s}^{-2}$ ).

The internal wave field was visualised by measuring instantaneous density gradient perturbations using the synthetic schlieren technique (Dalziel et al. 2000, 2007). This technique has been shown to be very effective to measure internal waves generated by TRI in experiments (Joubaud et al. 2012; Bourget et al. 2013; Grayson et al. 2022). A random pattern of dots behind the back wall of the tank was recorded using a CCD AVT Pike F-505 Camera of resolution 2452 $\times$ 2054 pixels at frame rate 6 or 8 Hz located at a large distance between the tank and the camera. The apparent displacement of the dots caused by the index of refraction perturbations in the flow was estimated using pattern-matching image velocimetry implemented in the PIVLAB software (Stamhuis & Thielicke 2014). From these image displacements, the gradient of the perturbation density field associated with spanwise-uniform internal waves could be measured, and from this, we determined the gradient of the perturbation buoyancy field $(\partial _x b, \partial _z b)$ .

3.2. Experiment results

The experiments examined wave evolution in three different non-uniform stratifications, these sets being denoted by Exp. I, II and III. For each stratification, sets of experiments with different frequencies and/or amplitudes of the parent mode-1 wave were performed. The different parameters of all the experiments are given in table 3. In Exp. I, the forcing amplitudes were generally small, and TRI was observed in only two experiments. Exp. II had the largest change in stratification from top to bottom, with $N^2(0)/N^2(-H) \approx 5.3$ . In Exp. III, the stratification change was $N^2(0)/N^2(-H) \simeq 4.2$ .

Table 3. Fitting parameters $N_0$ and $\unicode{x1D6E5} _N$ for the background buoyancy profile, for frequency $\omega _0$ , vertical velocity amplitude $w_0$ , and measured horizontal wavenumber $k_0$ of the parent wave. Note that the measured value of $k_0$ is consistent with the dispersion relation. We give frequencies $\varOmega _{\pm }$ and estimations of the horizontal $k_{\pm }$ and vertical $m_{\pm }$ wavenumbers and growth rate $\sigma _{\pm }$ of the two secondary waves. The errors for $k_+$ and $k_-$ are approximately $10$ %. The units $\textrm {cm}^{-1}$ and $\textrm {s}^{-1}$ indicate radians per centimetre and radians per second, respectively. The vertical wavenumbers can be compared to $m_0=\unicode{x03C0} /H\approx 0.097\ \textrm {cm}^{-1}$ to give an estimate $j_\pm \approx m_\pm /m_0$ .

Figure 4 illustrates the measured $x$ -component of the buoyancy gradient in an experiment with $N_0=1.04\,\textrm {s}^{-1}$ and $\unicode{x1D6E5} _N=0.59 \,\textrm {s}^{-1}$ (Exp. II), and a parent wave with $\omega _0=0.72\,\textrm {s}^{-1}$ and measured horizontal wavenumber $k_0=0.22\,\mbox{cm}^{-1}$ . The predicted horizontal group velocity of this wave is $0.93\,\mbox{cm}\, \rm {s^{-1}}$ , so that it takes of the order of a minute for the wave to propagate two wavelengths from the wavemaker. Viscosity does act moderately to attenuate the parent wave as it propagates away from the wavemaker. Estimates for the e-folding attenuation distance due to bulk viscosity and side-wall dissipation, respectively, are $c_g/(\omega _0 \delta _0) \simeq 32$ m and $c_g/(\omega _0 \unicode{x1D6E5} _0) \simeq 2.6$ m. Due to the dominant side-wall dissipation, the parent wave amplitude is predicted to decay by approximately $20$ % after propagating two wavelengths.

For a particular experiment, the snapshots shown in figure 4 reveal the spontaneous growth of two sibling waves after dimensionless time $N_0t\gtrsim 300$ . Through a temporal spectral analysis (described below), we determined the dominant frequencies of the sibling waves to be $\varOmega _+\simeq 0.43\,\textrm {s}^{-1}$ and $\varOmega _-\simeq 0.30\,\textrm {s}^{-1}$ , satisfying the temporal near-resonant condition $\varOmega _+ + \varOmega _- \simeq \omega _0$ . To visualise the two sibling waves more clearly, the vertical buoyancy gradient field was band-pass time-filtered about frequencies $\varOmega _+ \pm 0.02\,\textrm {s}^{-1}$ and $\varOmega _- \pm 0.02\,\textrm {s}^{-1}$ , as shown respectively for different times in figure 4 (middle and right-hand columns). Using these filtered wave fields, the horizontal wavenumbers for the parent wave and the two sibling waves, respectively $k_0$ and $k_\pm$ , were estimated (see table 3). The relative sibling wavenumbers in this case were $k_+ \approx 3.2 k_0$ and $k_-\approx 2.3k_0$ . Within error bars of $10$ %, the horizontal resonance condition is satisfied.

Figure 4. Left-hand column: snapshots of the horizontal buoyancy gradient field $\partial _x b$ from an experiment with the stratification of Exp. II (as shown in figure 3) and $\omega =0.72\,\textrm {s}^{-1}$ shown at times $N_0t=100$ (top row), $N_0t=300$ (middle row) and $N_0t=500$ (bottom row). The band-pass time-filtered vertical buoyancy gradient field $\partial _z b$ associated with subharmonic disturbances is shown for $\varOmega _+ = 0.43\ (\pm 0.02)\,\textrm {s}^{-1}$ (middle column) and $\varOmega _- = 0.30\ (\pm 0.02)\,\textrm {s}^{-1}$ (right-hand column).

Unlike the assumptions of theory, the sibling waves that develop in experiments are not pure vertical modes; they manifest as vertically propagating quasi-monochromatic waves that reflect from the upper and lower boundaries of the domain. For this reason, it is challenging to characterise their vertical mode number $j_\pm$ as used in theory. Instead, we estimate the local vertical wavenumber $m_\pm$ (measured in radians per metre), about a position located $10$ cm below the surface. Explicitly, we measured the vertical distance between two nodes of the buoyancy gradient field about this depth. From this vertical half-wavelength, $\lambda _{z\pm }/2$ , we define the wavenumber as $m_\pm = 2\unicode{x03C0} /\lambda _{z\pm }$ . For the experiment shown in figure 4, we find $m_+=1.25\, \textrm {cm}^{-1}$ and $m_-=1.13\, \textrm {cm}^{-1}$ . Defining $m_0 = \unicode{x03C0} /H\simeq 0.097\, \textrm {cm}^{-1}$ , which is the lowest vertical mode wavenumber in uniform stratification, we find $m_+-m_-\approx m_0$ .

To estimate the wave amplitude in experiments, a time series of the horizontal buoyancy gradient field was constructed at a horizontal position $5$ cm to the right of the wavemaker, at depth $z_1=-10$ cm. This temporal signal was then band-pass filtered about the parent wave frequency $\omega _0$ and the observed frequency of sibling waves $\varOmega _\pm$ . In particular, table 3 lists values of the parent wave amplitude in terms of its maximum vertical velocity $w_0$ . This is given in terms of the horizontal buoyancy gradient by

(3.1) \begin{equation} w_0 = \frac {\omega _0}{N^2(z_1)}\frac {\|\partial _x b\|}{k_0}. \end{equation}

The parent wave Reynolds number $\textit{Re}=U_0\kappa _0^{-1}/\nu$ is based upon the velocity magnitude, estimated by $U_0\equiv w_0 (N_0/\omega _0)$ , and the wavenumber magnitude $\kappa _0=\mbox{max}(k_0,\unicode{x03C0} /H)$ . For the experiments presented here, this ranges from $200$ to $1400$ , with the largest values occurring for moderate frequency parent waves ( $\omega \sim 0.7\ (\pm 0.2)\,N_0$ ).

Figure 5. For the experiment shown in figure 4, (a) time series of the normalised horizontal buoyancy gradient band-pass filtered for the parent wave (black) and the sibling waves (blue and red), and (b) the normalised spectrum of the full wave field. The growth of the sibling waves is estimated by the increase in sibling wave amplitude between the times indicated by the two dashed lines.

For the experiment shown in figure 4, the corresponding time series of the band-pass filtered horizontal buoyancy gradient field are shown in figure 5(a). Here, each of the time series corresponding to the parent and sibling waves is normalised by its respective maximum value. The spectrum $S$ of the unfiltered time series signal, shown in figure 5(b), exhibits several peaks, of which we identify the peaks satisfying $\varOmega _+ + \varOmega _- \simeq \omega _0$ , to get the sibling wave frequencies. These are used to construct the band-pass filtered signals shown in figure 5(a). Having separated the signal of the sibling waves from the parent, we see that the the parent wave amplitude saturates before significant growth of the sibling waves. By measuring the increase of the sibling wave amplitudes in time, we determine their growth rate $\sigma$ .

These analyses were performed for all the experiments, with results given in table 3. Generally, we see that the growth of sibling waves through TRI occurs only if the parent wave frequency is sufficiently close to $N_0$ , and its amplitude is sufficiently large. When TRI occurred, the temporal and spatial resonances were reasonably well satisfied within errors.

Overall, our experimental results are qualitatively similar to those of Joubaud et al. (2012), who performed experiments examining horizontally propagating internal modes in uniform stratification. As in our experiments, they found that sibling waves appear as vertically propagating disturbances rather than modes. They also observed a low-frequency cut-off for the occurrence of TRI. This cut-off frequency was larger in their experiments presumably because their forcing amplitude is half of that in our experiments.

4. Numerical simulations

We performed fully nonlinear numerical simulations of horizontally periodic, vertically confined internal modes in non-uniform stratification. These simulations solved the Boussinesq equations in two dimensions on the $x$ – $z$ plane, neglecting rotation. The numerical model used to perform these simulations is described in detail in Sutherland (2016).

4.1. Numerical model

The fully nonlinear equations of motion for a two-dimensional non-rotating Boussinesq fluid are given in terms of the spanwise vorticity $\zeta$ and buoyancy $b$ . These are represented in non-dimensional form based on the domain depth $H$ and characteristic buoyancy frequency $N_0$ , which is taken to be the maximum value of $N(z)$ at the top of the domain: $N_0 = N(0)$ . The equations are

(4.1) \begin{align} \zeta _{t} &= -u \zeta _{x} - w \zeta _{z} - b_{x} + \textit{Re}^{-1}\, \mathcal{D} \zeta, \end{align}

(4.2) \begin{align} b_{t} &= -u b_{x} - w b_{z} - N^{2} w + {(\textit{Re}\,\textit{Pr})}^{-1}\, \mathcal{D} b, \end{align}

in which subscripts denote partial derivatives.

The velocity fields $u$ and $w$ are given in terms of the streamfunction as $\boldsymbol{u}=(u,w)=(-\psi _z,\psi _x)$ . From the spanwise vorticity $\zeta$ , we find the streamfunction $\psi$ by inverting the Laplacian operator in $\zeta = u_{z} - w_{x} = -{\nabla} ^{2} \psi$ . The diffusion operator $\mathcal{D}$ is identical to the Laplacian operator, but only acts on Fourier components with a horizontal wavenumber greater than eight times that of the parent wave. For all simulations, the Prandtl number was $\textit{Pr}=1$ and the Reynolds number was $\textit{Re} \equiv N_0 H^2/\nu = 10^{5}$ . These values were chosen to damp numerical noise, ensuring the stability of the simulation while treating the parent and sibling waves as effectively inviscid. The background buoyancy frequency profile is defined by (2.23), with parameters chosen to be the same as those in the laboratory experiments.

We solved the equations in a domain with free-slip top and bottom boundaries with horizontally periodic boundary conditions. Fields were represented horizontally by Fourier components, and vertically by a uniformly spaced grid in real space. The vertical resolution of the domain and the number of Fourier components were both set to 2048. To initialise the simulation, eight wavelengths of a parent wave with horizontal wavenumber $k_{0}$ (having corresponding frequency and vertical structure given by (2.3)) were superimposed upon the background. For each of the buoyancy frequency profiles of the experiments, simulations were performed over a same range of parent wave frequencies as in the corresponding experiments. In all simulations, we chose an initial vertical displacement amplitude equivalent to $0.37$ cm in a $32.5$ cm domain, giving maximum vertical velocities ranging between $0.15$ and $0.33\,\textrm {cm}\, \rm {s^{-1}}$ depending upon the parent mode frequency. This range is of the order of the amplitudes measured in laboratory experiments. We further superimposed random noise somewhat arbitrarily inspired by the Garrett–Munk spectrum, with amplitude varying with frequency (using the dispersion relation for internal waves in uniform stratification) and vertical wavenumber $1/(\omega k_{z}^{2/5})$ .

4.2. Analysis methods

In simulations for which the growth of subharmonic sibling waves occurred, we performed diagnostics to measure the growth rates $\sigma$ , frequencies $\varOmega _{\pm }$ , and horizontal wavenumbers $k_{\pm }$ , and to estimate the vertical mode numbers $j_{\pm }$ .

We illustrate the application of diagnostics for a particular simulation with snapshots of the horizontal velocity field in figure 6. This shows the rapid growth of sibling waves between non-dimensional times $N_0t \approx 1400$ and $1800$ .

Figure 6. Snapshots at non-dimensional time $N_0t$ values (a) $0$ , (b) $1400$ , (c) $1600$ and (d) $1800$ of the non-dimensional horizontal velocity field $u_0/N_0H$ from a simulation with $k_0=0.168\,\textrm {cm}^{-1}$ , $\omega _0=0.7\,\textrm {s}^{-1}$ , $N_0=1.04\,\textrm {s}^{-1}$ , and buoyancy frequency decreasing with depth as $\unicode{x1D6E5} _N=0.59\,\textrm {s}^{-1}$ , which contains eight horizontal wavelengths of a mode-1 parent wave that interacts resonantly with the background noise, and excites sibling waves.

The horizontal wavenumbers and frequencies of the sibling waves were determined from horizontal time series of the horizontal velocity at a vertical level $z/H = -0.25$ , near the inflection depth of the parent wave horizontal velocity field indicated, for example, by the dashed lines in figure 6. The time series constructed from this simulation, shown in figure 7(a), reveals that after a delay, sibling waves grow rapidly to substantial amplitude between $N_0t\simeq 1700$ and $N_0t\simeq 2200$ .

Figure 7. Horizontal time series constructed at $z/H = -0.25$ of the development of the sibling waves, with parameters as given in figure 6, showing (a) the time series from the start of the simulation, and (b) the time series starting when the growth of sibling waves becomes apparent.

Focusing upon the sibling waves as they begin to grow to large amplitude, we analysed the horizontal time series between non-dimensional times $N_0t=1700$ and $2200$ , as shown in figure 7(b). Time- and space-transforming the signal over this time range produced a power spectrum of horizontal wavenumber and frequency. (Aliasing noise was reduced by first mirroring the signal in time and then Fourier transforming.) The result, shown in figure 8, clearly shows sharp peaks in the spectra at a small number of frequency–wavenumber combinations. The sum of the frequencies of the two strongest peaks is close to the parent wave frequency, as expected for sibling waves produced by TRI. Here, we denote the peak frequencies by $\varOmega _\pm$ , analogous to theoretical forcing frequencies. Likewise, the horizontal wavenumbers of these sibling waves exhibit pure resonance with $k_+ - k_- = k_0$ .

Figure 8. Magnitude of the power spectra of the windowed horizontal time series presented in figure 7(b), showing power associated with horizontal wavenumber and frequency normalised by those of the parent, respectively. For the $N_0$ , $k_0$ , $\unicode{x1D6E5} _N$ and $\omega _0$ provided in figure 6, the energy is strongly peaked at wavenumbers $k_{+} = 2k_{0}$ and $|k_{-}| = k_{0}$ , with corresponding frequencies such that $\varOmega _{+} = 0.439\,\textrm {s}^{-1}$ ( $\varOmega _+/\omega _0 =0.63$ ) and $\varOmega _{-} = 0.256\,\textrm {s}^{-1}$ ( $\varOmega _-/\omega _0 =0.37$ ).

As with the laboratory experiments, the sibling wave vertical mode numbers $j_\pm$ are not clearly defined in snapshots of the simulations. However, they can be estimated given the observed horizontal wavenumbers $k_\pm$ and frequencies $\varOmega _\pm$ of the sibling waves. From the theoretical dispersion relation found by solving (2.3) for a range of $k_\pm$ and $j_\pm$ , we determine the mode number most closely satisfying $\varOmega _\pm = \omega (k_\pm , j_\pm )$ , for given $\varOmega _\pm$ and $k_\pm$ .

To find the growth rate of the sibling waves, we examine power spectra of the kinetic energy constructed at each time step of the simulations. The power spectra are formed from Fourier transforms in the horizontal and vertical to give $P(k_\pm ,m_\pm )$ . Here, the vertical wavenumbers $m_\pm$ are given in terms of the estimated vertical mode numbers using $m_\pm = m_0 j_\pm$ , in which $m_0=\unicode{x03C0} /H$ . For given $k_\pm$ and $m_\pm$ , we extract the corresponding powers at successive times, which are normalised by the initial power $P_0$ associated with the parent wave. Taking the logarithm of the normalised power gives the plots shown in figure 9. We then find best-fit lines through $\ln {(P_\pm /P_{0})}$ versus $N_{0}t$ over a time range when the sibling waves grow exponentially. The slope of the lines gives exponential growth rate of power, with half this value giving the growth rate of the sibling wave amplitudes. Typically, the growth rates found for the $(k_+,m_+)$ and $(k_-,m_-)$ siblings are comparable. For each simulation, we characterised the growth rate $\sigma$ by their average value.

Figure 9. For the simulation shown in figure 6, log plots versus non-dimensional time of the normalised perturbation kinetic energy of sibling waves with (a) $k_+ = 2k_0$ and $m_+ = 5m_0$ , and (b) $k_- = k_0$ and $m_- = 5m_0$ . The slope of the best-fit (dashed) line during the exponential growth phase ( $1200\leqslant N_0t\leqslant 1500$ ) is shown for each plot. The corresponding mean growth rate of the sibling waves is $\sigma = 0.0049 N_0$ .

Although theory predicts exponential growth at the outset, simulations show a delay before sibling waves grow exponentially out of the background noise. We characterised the onset time $t_{c}$ at which the sibling waves begin to grow exponentially by taking the linear fit used to find the growth rate, and extrapolating back in time to when the line crosses the level of the mean power at early times. For the example shown in figure 9, this mean power level is given by $\ln (P(k_\pm ,m_\pm )/P_0) \simeq -16$ , and the mean onset time is $N_0t_c \simeq 697$ .

4.3. Simulation results, and comparison with theory and experiments

Here, we present the quantitative results for a range of numerical simulations, and we compare them with predictions from theory and the results of laboratory experiments. In all simulations, the initial parent wave amplitude was fixed to have initial maximum vertical displacement $0.37$ cm. Because the simulations are effectively inviscid, we found that varying this choice of amplitude did not significantly influence the observed frequency and wavenumbers of the sibling waves, though it did affect their growth rate. Theoretical predictions were computed for inviscid fluid and including the combined effect of viscosity acting in the bulk and at the side walls. For these calculations, the parent wave amplitude was taken to be the same as the maximum vertical velocity $w_0$ measured in corresponding laboratory experiments.

We expect from theory that the growth rate of sibling waves should exhibit an increasing trend with increasing parent wave frequency at fixed amplitude. This is indeed what is found in simulations and experiments, as shown in figure 10. Because growth rate is dependent upon the initial parent wave amplitude (as given by (2.15) and (2.18)), the theoretical and measured growth rates are normalised by $H/A_\xi$ to provide a proper comparison between simulations, experiments and theory. Comparing inviscid theory to simulations, we generally find excellent agreement if the parent wave frequency is not too large; theory underpredicts the growth rate found in simulations if $\omega _0/N_0 \gtrsim 0.7$ . Including viscosity, theory overpredicts the growth rate of sibling waves observed in laboratory experiments except in Exp. III with $\omega _0/N_0=0.92$ (and $w_0=0.06\,\mbox{cm}\, \rm {s^{-1}}$ ), in which case no growth of sibling waves is predicted even though large growth of high vertical wavenumber sibling waves is observed in experiments.

Figure 10. Normalised average growth rate of the sibling waves as it depends on the normalised frequency of the parent wave in three different stratification profiles corresponding to (a) Exp. I, (b) Exp. II and (c) Exp. III. We compare simulations (red squares) and experiments where TRI was seen (blue triangles) with theory including viscosity (blue circles) and neglecting viscosity (red crosses).

Figure 11 compares the near-resonant frequencies of the sibling waves found in simulations and experiments to those in theory. In the last case, we find the resonant sibling waves with the largest growth rate across all horizontal wavenumbers and vertical mode numbers, plotting their corresponding frequencies. This calculation is performed with and without viscous effects included. We find that the frequencies of the sibling waves grow proportionally to the parent wave and around a line $\varOmega _\pm \approx 0.5 \omega _0$ , which would mark the special case of parametric subharmonic instability. At lower relative parent wave frequencies, the sibling waves have frequencies close to half that of the parent. If the parent wave has frequency closer to $N_0$ , then simulations and experiments show a greater spread between the values of $\varOmega _+$ and $\varOmega _-$ , though their sum is close to $\omega _0$ . This spread is also evident in theory, with the spread being narrower if viscous effects are neglected. With viscous effects included, theory is in closer alignment with experiments, as expected.

Figure 11. Normalised frequencies of sibling waves, $\varOmega _{\pm }$ , as they depend on the normalised frequency of the parent wave with background stratification corresponding to (a) Exp. I, (b) Exp. II and (c) Exp. III, as given in table 3. We plot values corresponding to the maximum growth rate measured in simulations (squares) and experiments (triangles), and predicted by theory including viscous effects (solid circles) and neglecting viscosity (crosses). Values of the $+$ sibling are plotted in blue; values of the $-$ sibling are plotted in red. The dashed lines indicate where the sibling wave frequency is half that of the parent.

Figure 12 shows the measured and predicted horizontal and vertical wavenumbers of the $+$ sibling waves from simulations, experiments and theory. For experiments, the estimated vertical wavenumber $m_+$ is normalised by $m_0=\unicode{x03C0} /H$ , in order to make a more direct comparison with the vertical mode number $j_+$ . Because $m_+$ is determined from measurements around $z=-10$ cm where the background is relatively strongly stratified, the value of $m_+/m_0$ likely overestimates the actual vertical mode number should the observed waves be manifest as modes. Though the results for the $-$ sibling waves are not plotted, we generally find that the near resonance condition is satisfied such that $k_++k_-\simeq k_0$ , with $k_-\lt 0$ . Our simulation, experiment and theory results show a large variation in the wavenumbers, with no clear dependence upon the relative frequency of the parent wave, though there is an apparent decrease (increase) in the sibling wave horizontal (vertical) wavenumber with increasing parent wave frequency. In comparison with simulations and inviscid theory, experiments and viscous theory generally have smaller horizontal and vertical wavenumbers.

Figure 12. The normalised horizontal and vertical wavenumbers of the $+$ sibling as they depend on the normalised frequency of the parent wave with background stratification corresponding to (a) Exp. I, (b) Exp. II and (c) Exp. III, as given in table 3. Plotted values show measurements from simulations (red squares) and experiments (blue triangles), and predictions from theory including viscosity (blue circles) and without viscosity (red crosses).

Figure 13. Normalised onset times from simulations of the development of the two sibling waves as it depends on the normalised frequency of the parent wave. The legend and the colour of the point indicate the type of stratification profile (given in table 3) corresponding to each point.

Finally, we examine the onset times of TRI measured in simulations, as shown in figure 13. Generally, the onset time is smaller if the parent wave frequency is larger. While theory predicts immediate exponential growth, in simulations and experiments, the sibling waves evolve from background noise. A tentative explanation for the longer onset times arising from lower frequency parent waves in simulations (and experiments) is that the predicted vertical structure of the sibling waves generally has lower vertical wavenumber, hence larger vertical scale. It takes more time to plant the seed of a sibling wave if spatially large coherent scales are to develop from noise.

5. Discussion and conclusions

We have developed a theory for the near-resonance growth of sibling waves in non-uniform stratification from a low vertical mode parent wave with and without viscosity. Rather than the Floquet analysis approach of Akylas & Kakoutas (2023), which derived the growth rate and frequency of sibling waves from the solvability condition of a differential eigenvalue problem, we used orthogonality of modes to derive time evolution equations for the amplitude of sibling waves. Our approach was readily adapted to include viscous effects not just in the bulk of the fluid but also due to side-wall dissipation. In the study of Akylas & Kakoutas (2023), the parent wave amplitude was used as a perturbation parameter that determined the growth rate at resonance as well as its decay for near-resonant triads. Instead, we use the triadic frequency mismatch $\epsilon$ as our perturbation parameter. While rigorously the growth rate is a function of parent wave amplitude and $\epsilon$ , in our analyses we focus on the case of pure resonance ( $\epsilon =0$ ). The advantage of using $\epsilon$ as a parameter is that it provides an efficient means numerically to locate sibling wave horizontal wavenumbers and vertical mode numbers for which pure resonance occurs.

Our prediction for pure resonant sibling wave excitation from the $k_+-k_-=k_0$ interaction agrees with the prediction of Akylas & Kakoutas (2023), though we show that the dominant growth rate occurs for the $k_++k_-=k_0$ interaction if $|k_\pm |\gt 0.5k_0$ . These results, particularly the predicted growth rate of sibling waves, are validated quantitatively by numerical simulations of effectively inviscid waves. In the theory including viscous effects, both in the bulk of the fluid and due to side-wall dissipation, the predicted sibling wave growth rate is smaller than corresponding inviscid theory. However, these predictions are much larger than measurements from most laboratory experiments. This is likely due to the viscous attenuation of the parent wave itself, whose amplitude moderately decreases with distance from the wavemaker. Thus contrary to the assumptions of theory, the parent wave is not horizontally periodic, inhibiting the development of sibling waves with long horizontal wavelength, and the decreasing amplitude of the parent would further lower the growth rate of the sibling waves. The viscous theory does improve the prediction of the sibling wave frequencies, these differing more from the subharmonic frequency $\omega _0/2$ in comparison with the predictions of inviscid theory, particularly if the parent mode frequency is closer to the maximum background buoyancy frequency $N_0$ . The predicted and measured wavenumbers of the sibling waves show more scatter, though there is a general trend for the sibling wave horizontal wavenumbers to decrease, and the vertical mode numbers to increase, with increasing parent wave frequency. The scatter, particularly for inviscid theory, is attributed to the fact that the growth rates of sibling waves are similar across a wide range of sibling wave horizontal (and hence vertical) wavenumbers. Although the parent mode propagated in non-uniform stratification, in neither simulations nor experiments was there evidence for self-interaction leading to significant growth of superharmonics. This is anticipated because in most cases, the parent mode frequency was more than half the maximum background buoyancy frequency, thus precluding the generation of frequency-doubled disturbances.

Perhaps the most significant finding of this study is that the (theoretical, simulated and experimental) growth rate of sibling waves becomes negligibly small if the parent wave frequency is small. Furthermore, in simulations, the onset time for the development of sibling waves becomes longer if the parent wave frequency is smaller. For the background stratifications and amplitudes considered in experiments, negligible growth of sibling waves is predicted if $\omega /N_0\lesssim 0.7$ . By extrapolation of theory and simulations, this work suggests that sibling waves manifest as pure resonant modes may not develop through TRI of the oceanic internal tide for which $\omega _0\ll N_0$ . Not only would their growth rate be small, but the time for them to develop from background noise would be delayed by a significant onset time. Nonetheless, with forcing by the low-frequency internal tide, and with consideration of background rotation, it may be possible for large vertical wavenumber sibling wave packets to grow near resonantly over ranges at mid-depth where the stratification is quasi-uniform (Gururaj & Guha 2020). Such development is evident in the numerical study of Hazewinkel & Winters (2011), who potentially seeded such disturbances as an initial condition. It is also evident in our laboratory experiments for which sibling waves emerge as vertically propagating disturbances rather than vertical modes. Such manifestation of sibling waves helps to explain why strong growth of sibling waves with high vertical wavenumber are observed in Exp. III for a parent wave having frequency $0.92N_0$ even though no growth is predicted by viscous theory: the assumption in theory that sibling waves have a vertical mode structure is over-restrictive. What is similar between our experiments and the simulations of Hazewinkel & Winters (2011) is that the parent waves are forced at one side of the domain, whereas our simulations are initialised with horizontally periodic waves. Separately, we have performed simulations that locally force rightward propagating parent waves, with these giving evidence for sibling waves growing with non-modal structure. This suggests the need to adapt theory to account for the structure of horizontally modulated parent waves. A detailed examination of these simulations is the subject of future work.

Despite what agreement there is between theory and corresponding laboratory experiments and numerical simulations, the discrepancy between the predictions extrapolated to oceanic internal tides motivates future work to explore TRI theory including more realistic stratification, the effects of background rotation, and the horizontal non-uniformity of the low-mode internal tide.