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 Reference Akylas and Kakoutas2023).

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 Reference Akylas and Kakoutas2023), which presumably grow out of background noise. Using the ‘pump-wave’ approximation (Craik & Adam Reference Craik and Adam1978; Young et al. Reference Young, Tsang and Balmforth2008; Gururaj & Guha Reference Gururaj and Guha2020), 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 (Reference Akylas and Kakoutas2023).

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 (Reference Akylas and Kakoutas2023) 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 (Reference Davis and Acrivos1967), 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. (Reference Joubaud, Munroe, Odier and Dauxois2012) 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. Reference Sauret, Boulogne, Cappello, Dressaire and Stone2015). 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. Reference Bourget, Dauxois, Joubaud and Odier2013; Sutherland & Jefferson Reference Sutherland and Jefferson2020). As with periodic waves (Bourget et al. Reference Bourget, Dauxois, Joubaud and Odier2013), 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.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. (Reference Varma, Pacary, Dauxois and Joubaud2024). This procedure is an adaptation of the standard double bucket method (Oster Reference Oster1965; Hill Reference Hill2002; Economidou & Hunt Reference Economidou and Hunt2009). 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. (Reference Dossmann, Pollet, Odier and Dauxois2017), Husseini et al. (Reference Husseini, Varma, Dauxois, Joubaud, Odier and Mathur2020) and Varma et al. (Reference Varma, Pacary, Dauxois and Joubaud2024). 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. Reference Dalziel, Hughes and Sutherland2000, Reference Dalziel, Carr, Sveen and Davies2007). This technique has been shown to be very effective to measure internal waves generated by TRI in experiments (Joubaud et al. Reference Joubaud, Munroe, Odier and Dauxois2012; Bourget et al. Reference Bourget, Dauxois, Joubaud and Odier2013; Grayson et al. Reference Grayson, Dalziel and Lawrie2022). 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 Reference Stamhuis and Thielicke2014). 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. (Reference Joubaud, Munroe, Odier and Dauxois2012), 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.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.