3.1. Results of the measurements

The apparently sinusoidal temporal variation of the interface elevation above the wavemaker is plotted in figure 5 together with the records of the wavemaker motion at different forcing frequencies. An effective resonance is observed at $\alpha$ about 0.93; at lower frequencies, the interface oscillations at this lateral location are approximately in antiphase relative to the wavemaker motion. The phase difference seems to decrease notably as the forcing frequency exceeds the effective resonance; in figure 5(d), the oscillations are nearly in phase.

Figure 6.

(a) Instantaneous interface shape for instants with maximum and minimum values of $\eta$ at $x=0$ . Dashed and solid lines are for the frequencies from figure 4(a) and 4(b), respectively. (b) Root-mean-square values of interface elevation at different forcing frequencies.

Several instantaneous interface shapes are plotted in figure 6(a). While those shapes indeed seem to be dominated by the resonant standing wave with $m=3$ , they clearly contain contributions from additional shorter waves. The appearance of shorter waves is more pronounced for frequencies below the resonant frequency (dashed curves in figure 6 a). The root-mean-square $\langle \eta ^2\rangle ^{1/2}$ values of the interface elevation were calculated separately for each horizontal location $x$ (figure 6 b) and averaged over 220 periods corresponding to the quasi-steady state. For a monochromatic standing wave, the expected values of $\langle \eta ^2 \rangle ^{1/2}$ are zero at nodes spaced by half a wavelength, and reach their maximum at antinodes. Figure 6(b) demonstrates that this pattern describes the actual behaviour only approximately; the values of $\langle \eta ^2 \rangle ^{1/2}$ do not vanish at nodes, and the shapes around the antinodes are quite flat. At frequencies well below the effective resonance, the measured distributions of $\langle \eta ^2 \rangle ^{1/2}$ are nearly flat everywhere. The complexity of the interface wavy shape manifests itself in the wavenumber spectra of the complex amplitude of the interface shape band-pass-filtered at the forcing frequency (figure 7). Along with the global maxima at $m=3$ , the additional weaker local maxima at $m\approx 9$ and $m\approx 15$ are seen. The former can be attributed to the cubic nonlinearity and the interaction of the dominant wavemode with itself, while the latter has no analogy in the surface waves.

Figure 7.

Wavenumber amplitude spectra of the temporal amplitudes of the interface shape band-pass-filtered at the forcing frequency, for (a,b) the cases shown in figure 4(a) and 4(b), respectively.

The frequency spectrum of the interface elevation at the antinode is dominated by a harmonic at the wavemaker driving frequency (figure 8 a); this harmonic, while much weaker, is also notable at the node of the standing wave (figure 8 b). At both the node and the antinode, there is a significant peak at double the wavemaker frequency, which is more pronounced when the forcing frequency is below the effective resonance (figure 8 a).

Figure 8.

Normalised amplitude frequency spectra of interface elevation for wavemaker amplitude $F =5.6$ mm: (a) at the antinode ( $x/\lambda _3=0.5$ ); (b) at the node ( $x/\lambda _3=0.25$ ).

Figure 9.

(a) Interface elevation amplitude in the antinode ( $x/\lambda _3=0.5$ ) for different forcing amplitudes, the coloured symbols represent the pycnocline thickness $\delta$ . (b) Phase difference between the surface elevation and the wavemaker displacement for the cases as in (a).

The variation of the amplitude of the interface elevation with forcing frequency is plotted in figure 9(a) for three values of the wavemaker amplitude obtained in a single series of experimental runs with increasing forcing frequency. The effective resonance downshifts as the forcing amplitude increases. The pycnocline thickness $d$ grows so that $\delta$ changes from about $0.17$ to about $0.25$ as shown in figure 9 by colours of the markers. The decrease of the phase shift between the interface elevation above the wavemaker and the wavemaker from $\pi$ towards zero is plotted for those runs in figure 9(b).

The effect of the pycnocline thickness widening in the course of the experiment is studied in an additional series of experiments carried out for $n=4$ . Those runs started from the lowest frequency with $\alpha =0.92$ , attained the value of forcing frequency $\alpha =0.978$ , after which the values of $\alpha$ were gradually decreased. Figure 10 shows the downshift of the resonant frequency with an increase of $\delta$ .

Figure 10.

Interface elevation amplitude at the antinode ( $n=4$ , $x/\lambda _4=0.25$ , $F=10$ mm). The coloured symbols represent the pycnocline thickness $\delta$ .

3.2. Interim conclusions

Shadowgraph experiments reveal the basic features of internal waves in a closed basin containing two liquids with different densities separated by a relatively thin pycnocline. The complex structure of the wave, excited by a wavemaker operating in the vicinity of the resonant frequency corresponding to the selected two-layer mode, is exposed in the temporal records of the visible interface elevation. The temporal spectra at any longitudinal location are dominated by discrete harmonics corresponding to the wavemaker forcing frequency $\omega$ and its second harmonic 2 $\omega$ . The resonant wave in the vicinity of the third-mode resonant frequency clearly deviates from the sinusoidal shape: the wavenumber spectra of records band-pass-filtered at the forcing frequency (figure 7) contain multiple components with $m$ in the vicinity of the resonant mode $n=3$ , components with $m\approx 3n$ that may be attributed to cubic nonlinearity and $m\approx 15$ that, as will be shown in the following, are related to the finite pycnocline thickness. The effective resonance frequency in the experiments is apparently downshifted from the two-layer model prediction $\omega _n$ ; this downshift increases with an increase of the wavemaker amplitude $F$ .

The thickness of the pycnocline $d$ increases in the course experimental sessions. The increase in $d$ affects notably the effective resonance of the system; it also contributes to the downshift of the effective resonant frequency. It should be emphasised that while the density distribution can be monitored between the runs using the PCP, there is no way to control it.

The shadowgraph experiments have the advantage of wide spatial coverage, but the sequence of recorded images only provides information on the movement of one surface representing the effective interface. This technique is thus incapable of revealing the full internal wave structure inside a pycnocline of finite thickness, suggesting that a different experimental approach is needed. The PIV technique can provide the necessary two-dimensional information; however, carrying out these experiments for a wide range of parameters requires extensive resources. An efficient experimental programme can benefit from guidelines supplied by a theoretical model that gives prior estimates of the wavefield behaviour.

4.1. Two-layer model

Consider two-dimensional flow in the $(x,z)$ plane of two effectively infinite deep layers of homogeneous liquids separated by interface $z=\eta (x,t)$ . The variables are rendered dimensionless as

(4.1) \begin{align} (\tilde {x},{\tilde {z}})=k_n(x,z),\quad \tilde {t}=\omega t .\end{align}

The velocity is scaled by $\omega F$ ; the interface location by $F$ :

(4.2) \begin{align} \tilde {\eta }=\frac {\eta }{F }. \end{align}

The dimensionless wall locations are $\tilde {x}=-\pi$ and $\tilde {x}=(n-1)\pi$ . The mean position of the interface is $\tilde {\eta }=0$ . The dimensionless parameters $r$ and $h$ ,

(4.3) \begin{align} r=k_nR,\quad h=k_nH, \end{align}

define the geometry of the wavemaker; $-H$ is the mean vertical coordinate of its centre. For $n=3$ , the experiments correspond to $r=0.26$ , $h=0.71$ . The instantaneous location of the wavemaker centreis

(4.4) \begin{align} {\tilde {z}}_{c}(\tilde {t})=-h+\varepsilon \sin \tilde {t},\quad \varepsilon =k_nF . \end{align}

The flow is assumed irrotational in both layers; the dimensionless velocity potentials $\varphi _1$ and $\varphi _2$ in the upper and lower liquids satisfy the Laplace equation:

(4.5a) \begin{align}&\qquad\qquad\qquad \nabla ^2\varphi _1=0,\quad -\pi \lt \tilde {x}\lt (n-1)\pi, \quad {\tilde {z}}\gt \tilde {\eta }(\tilde {x},\tilde {t}), \end{align}

(4.5b) \begin{align}& \nabla ^2\varphi _2=0,\quad -\pi \lt \tilde {x}\lt (n-1)\pi, \quad {\tilde {z}}\lt \tilde {\eta }(\tilde {x},\tilde {t}),\quad \tilde {x}^2+({\tilde {z}}-{\tilde {z}}_{c}(\tilde {t}))^2\gt r^2. \end{align}

The non-penetration conditions at the rigid boundaries (walls and the wavemaker) are

(4.6a) \begin{align}&\qquad\qquad \frac {\partial \varphi _1}{\partial \tilde {x}}=0,\quad \tilde {x}=-\pi, \ (n-1)\pi, \end{align}

(4.6b) \begin{align}&\qquad\qquad\qquad |\nabla \varphi _1|\to 0,\quad {\tilde {z}}\to \infty, \end{align}

(4.6c) \begin{align}&\qquad\qquad \frac {\partial \varphi _2}{\partial \tilde {x}}=0,\quad \tilde {x}=-\pi, \ (n-1)\pi, \end{align}

(4.6d) \begin{align}& \frac {\partial \varphi _2}{\partial n}=\frac {1}{\varepsilon }\frac {{\rm d}{\tilde {z}}_{c}}{{\rm d}\tilde {t}}\frac {{\tilde {z}}-{\tilde {z}}_{c}}{r^2},\quad \tilde {x}^2+({\tilde {z}}-{\tilde {z}}_{c}(\tilde {t}))^2=r^2, \end{align}

(4.6e) \begin{align}&\qquad\quad\qquad |\nabla \varphi _2|\to 0,\quad {\tilde {z}}\to -\infty .\end{align}

The continuity of the vertical velocity and of the pressure at the interface connects the interface location $\tilde {\eta }$ with the potentials:

(4.7) \begin{align} &\frac {\partial \varphi _1}{\partial {\tilde {z}}}=\frac {\partial \varphi _2}{\partial {\tilde {z}}}\notag\\{\tilde {z}}=\tilde {\eta }(\tilde {x},\tilde {t}):\quad &\frac {\partial \varphi _1}{\partial \tilde {t}}+\frac {\varepsilon }{2}\left (\nabla \varphi _1\right )^2=\frac {\rho _2}{\rho _1}\left [\frac {\partial \varphi _2}{\partial \tilde {t}}+\frac {\varepsilon }{2}\left (\nabla \varphi _2\right )^2\right ]+\frac {(\rho _2-\rho _1)gk_n}{\rho _1\omega ^2}\tilde {\eta } \notag\\ &\frac {\partial \tilde {\eta }}{\partial \tilde {t}}+\frac {\partial \varphi _1}{\partial \tilde {x}}\frac {\partial \tilde {\eta }}{\partial \tilde {x}}=\frac {\partial \varphi _1}{\partial {\tilde {z}}}. \end{align}

Neglecting terms of the order of $(\rho _2-\rho _1)/\rho _1$ and $\varepsilon$ and eliminating $\eta$ yields a simplified linear problem for potentials:

(4.8a) \begin{align} &\qquad\qquad\nabla ^2\varphi _1=0,\quad -\pi \lt \tilde {x}\lt (n-1)\pi, \quad {\tilde {z}}\gt 0, \end{align}

(4.8b) \begin{align} &\nabla ^2\varphi _2=0,\quad -\pi \lt \tilde {x}\lt (n-1)\pi, \quad {\tilde {z}}\lt 0,\quad \tilde {x}^2+({\tilde {z}}+h)^2\gt r^2, \end{align}

(4.8c) \begin{align} &\qquad\qquad\qquad\frac {\partial \varphi _1}{\partial \tilde {x}}=0,\quad \tilde {x}=-\pi, \ (n-1)\pi, \end{align}

(4.8d) \begin{align} &\qquad\qquad\qquad\qquad|\nabla \varphi _1|\to 0,\quad {\tilde {z}}\to \infty, \end{align}

(4.8e) \begin{align} &\qquad\qquad\qquad\frac {\partial \varphi _2}{\partial \tilde {x}}=0,\quad \tilde {x}=-\pi, \ (m-1)\pi, \end{align}

(4.8f) \begin{align} &\qquad\qquad\frac {\partial \varphi _2}{\partial n}=\frac {{\tilde {z}}+h}{r^2}\cos t,\quad \tilde {x}^2+({\tilde {z}}+h)^2=r^2, \end{align}

(4.8g) \begin{align} &\qquad\qquad\qquad\quad|\nabla \varphi _2|\to 0,\quad {\tilde {z}}\to -\infty, \end{align}

(4.8h) \begin{align} &\qquad\qquad\qquad\qquad\frac {\partial \varphi _1}{\partial {\tilde {z}}}=\frac {\partial \varphi _2}{\partial {\tilde {z}}},\quad {\tilde {z}}=0, \end{align}

(4.8i) \begin{align} &\qquad\qquad\quad\frac {\partial ^{2} \varphi _1}{\partial \tilde {t}^2}=\frac {\partial ^2 \varphi _2}{\partial \tilde {t}^2}+\frac {2}{\alpha ^2}\frac {\partial \varphi _2}{\partial {\tilde {z}}},\quad {\tilde {z}}=0. \end{align}

The violation of the symmetry between the upper and the lower liquids is due to the presence of the wavemaker under the interface. It leads to the failure of the solution of Lamb (Reference Lamb1932) as it does not satisfy the boundary condition at the wavemaker (4.8f ). The Lamb solution in the upper liquid is retained:

(4.9) \begin{equation} \varphi _1=\sum \limits _{m=1}^\infty a_m\varphi _{1,m},\quad \varphi _{1,m}=\frac {1}{2}\left [\cos \left (\frac {m}{n}(\tilde {x}+\pi )\right )\exp \left (-\frac {m}{n}{\tilde {z}}\right )\exp (-i\tilde {t})+\text{c.c.}\right ]. \end{equation}

In the lower liquid, the solution is represented as a superposition of the forcing potential $\varphi _0$ and wavemodes $\varphi _{2,m}$ :

(4.10) \begin{equation} \varphi _2=\varphi _0+\sum \limits _{m=1}^\infty b_m\varphi _{2,m}. \end{equation}

The potential $\varphi _0$ satisfies equation (4.8b ) with boundary conditions (4.8e –g ), including the inhomogeneous boundary condition (4.8f ). The potentials $\varphi _{2,m}$ satisfy the homogeneous boundary condition at the wavemaker:

(4.11) \begin{align}&\qquad\qquad\qquad\qquad \frac {\partial \varphi _{2,m}}{\partial n}=0,\quad \tilde {x}^2+({\tilde {z}}+h)^2=r^2\nonumber\\& \varphi _{2,m}=\frac {1}{2}\left \{\left [\cos \left (\frac {m}{n}(\tilde {x}+\pi )\right )\exp \left (\frac {m}{n}{\tilde {z}}\right )+\varphi _{2,m}^{\prime}\right ]\exp (-i\tilde {t})+\text{c.c.} \right \} .\end{align}

The symmetry of the functions $\varphi _0$ and $\varphi _{2,m}^{\prime}$ with respect to ${\tilde {z}}=-h$ is also required for uniqueness; they are represented by series of multipoles placed at the wavemaker centre and their successive mirror reflections with respect to the basin walls (Mogilevskiy et al. Reference Mogilevskiy, Kalenko, Zemach and Shemer2024). The leading terms are

(4.12) \begin{align} &\,\,\, \varphi _0=\frac {1}{2}\left [-\frac {r^2(z+h)}{x^2+(z+h)^2}\exp (-i\tilde {t})+\text{c.c.}\right ];\nonumber\\& \varphi _{2,m}^{\prime}=\frac {1}{2}\left [-\frac {m}{n}\exp \left (-\frac {m}{n}h\right )\frac {r^2(z+h)}{x^2+(z+h)^2}\right ]. \end{align}

The forcing potential $\varphi _0$ introduces disturbances of the order of $r^2/h$ to all spatial modes (figure 11); the wavenumber spectrum has peaks at $m=n,\ 2n, \ldots$ , since the wavemaker is located at the antinode of those modes. The response of the system is defined by the boundary conditions (4.8h ) and (4.8i ) representing the continuity of the vertical velocity and the pressure at the interface. The substitution of (4.9), (4.10) yields the system of linear algebraic equations for $a_m$ and $b_m$ :

(4.13) \begin{equation} \begin{array}{c} \displaystyle { \int \limits _{-\pi }^{(n-1)\pi }\left [ \sum \limits _{m=1}^\infty a_m\frac {\partial \varphi _{1,m}(\tilde {x},0)}{\partial {\tilde {z}}}-\sum \limits _{m=1}^\infty b_m\frac {\partial \varphi _{2,m}(\tilde {x},0)}{\partial {\tilde {z}}}\right .}\\ \displaystyle {\left . -\frac {\partial \varphi _0(\tilde {x},0)}{\partial {\tilde {z}}} \right ]\cos \left (\frac {l}{n}(\tilde {x}+\pi )\right ){\rm d}\tilde {x}=0},\\ \displaystyle { \int \limits _{-\pi }^{(n-1)\pi }\left \{ \sum \limits _{m=1}^\infty a_m\varphi _{1,m}(\tilde {x},0)-\sum \limits _{m=1}^\infty b_m\left [\varphi _{2,m}(\tilde {x},0)-\frac {2}{\alpha ^2}\frac {\partial \varphi _{2,m}(\tilde {x},0)}{\partial {\tilde {z}}}\right ]\right .}\\ \displaystyle {-\left .\left [\varphi _{0}(\tilde {x},0)-\frac {2}{\alpha ^2}\frac {\partial \varphi _{0}(\tilde {x},0)}{\partial {\tilde {z}}}\right ] \right \}\cos \left (\frac {l}{n}(\tilde {x}+\pi )\right )d\tilde {x}=0},\\ l=1\ldots \infty .\end{array} \end{equation}

The free terms associated with $\varphi _0$ are of the order of $r^2/h$ (see (4.12)). It can be shown that the eigenvalue of the system (4.13) is $\alpha =m/n-\Delta \alpha _m$ , where $0\lt \Delta \alpha _m\sim r^2/(\pi n h)$ corresponds to the downshift of the resonant frequency due to the finite wavemaker size. For other values of $\alpha$ , (4.13) has a unique solution.

Figure 11.

(a) Normalised forcing potential at ${\tilde {z}}=0$ and (b) its wavenumber amplitude spectrum.

This simplified two-layer model demonstrates the general approach; however, the predicted downshift of the resonant frequency from $\alpha =1$ is about 1 %, an order of magnitude smaller than that measured in experiments (§ 3). The experimental results plotted in figure 10 indicate the significant effect of the thickness of the pycnocline, which is neglected in the two-layer model.

4.2. Eigenfunctions for continuous stratification

In the linear approximation, the complex amplitude function for the vertical velocity $W_m({\tilde {z}})$ for a wave with the wavemode number $m$ satisfies the Taylor–Goldstein equation (Goldstein Reference Goldstein1931; Taylor Reference Taylor1931) for a given dimensionless buoyancy frequency distribution ${\tilde {N}}({\tilde {z}}$ ):

(4.14) \begin{equation} W_m^{\prime\prime}+\left (\frac {2}{\alpha ^2}{\tilde {N}}^2-1\right )\frac {m^2}{n^2}W_m=0 \end{equation}

and the boundary conditions

(4.15) \begin{equation} W_m\to 0,\quad {\tilde {z}}\to \pm \infty, \end{equation}

where primes denote derivatives with respect to $\tilde {z}$ . The analysis of this Sturm–Liouville problem shows that for any given $m$ , there is an infinite number of solutions $W_m^{j}$ and corresponding $\alpha _{m,j}$ that are orthogonal with the weight function ${\tilde {N}}^2$ (Morse & Feshbach Reference Morse and Feshbach1946):

(4.16) \begin{equation} \int \limits _{-\infty }^\infty W_m^{j1}W_m^{j2}{\tilde {N}}^2{\rm d}{\tilde {z}}=0,\quad j_1\neq j_2. \end{equation}

For the piecewise-constant distribution ${\tilde {N}}^2({\tilde {z}})$ :

(4.17) \begin{equation} {\tilde {N}}^2=\left \{ \begin{array}{cr} 0,& {\tilde {z}}\lt -\delta /2\\ \frac {1}{\delta}& |{\tilde {z}}|\lt \delta /2\\ 0,& {\tilde {z}}\gt \delta /2 \end{array} \right . \end{equation}

the analytical solution for $W_{m}^j$ is

(4.18) \begin{equation} \begin{array}{c} W_{m}^{j}=\left \{ \begin{array}{ll} \exp \left (\dfrac {m}{n}{\tilde {z}}\right ),& {\tilde {z}}\lt -\dfrac {\delta }{2}\\[6pt] \cos (\gamma _{m}^{j} {\tilde {z}})\exp \left (-\dfrac {m}{n}\dfrac {\delta }{2}\right )/\cos \left (\gamma _{m}^{j}\dfrac {\delta }{2} \right ),& |{\tilde {z}}|\lt \dfrac {\delta }{2}\\[6pt] \exp \left (-\dfrac {m}{n}{\tilde {z}}\right ),& {\tilde {z}}\gt \dfrac {\delta }{2} \end{array} \right . \quad j \text{ odd},\\ \\[-2pt] W_{m}^{j}=\left \{ \begin{array}{ll} -\exp \left (\dfrac {m}{n}{\tilde {z}}\right ),& {\tilde {z}}\lt -\dfrac {\delta }{2}\\[6pt] \sin (\gamma _{m}^{j} {\tilde {z}})\exp \left (-\dfrac {m}{n}\dfrac {\delta }{2}\right )/\sin \left (\gamma _{m}^{j}\dfrac {\delta }{2} \right ),& |{\tilde {z}}|\lt \dfrac {\delta }{2}\\[6pt] \exp \left (-\dfrac {m}{n}{\tilde {z}}\right ),& {\tilde {z}}\gt \dfrac {\delta }{2} \end{array} \right .\quad j \text{ even},\\ \end{array} \end{equation}

where

(4.19) \begin{equation} \left (\gamma _{m}^{j}\right )^2=\left (\frac {2}{\left (\alpha _{m}^{j}\right )^2\delta }-1\right )\frac {m^2}{n^2} \end{equation}

are roots of one of the equations:

(4.20a) \begin{align} &\,\gamma \tan \frac {\gamma \delta }{2}=\frac {m}{n}, \end{align}

(4.20b) \begin{align} &\gamma \cot \frac {\gamma \delta }{2}=-\frac {m}{n} \end{align}

(Miropol’sky & Shishkina Reference Miropol’sky and Shishkina2001). Equations (4.20a ) and (4.20b ) correspond to the eigenvalues with odd and even superscript indices. The function $W_{m}^j$ has $j-1$ zeros within the pycnocline; it is symmetric for odd $j$ and antisymmetric for even $j$ with respect to ${\tilde {z}}=0$ . The resonant mode in the two-layer model corresponds to $\alpha _*=\alpha _{n}^1$ , in the limit $\delta \ll 1$ ,

(4.21) \begin{equation} \alpha _*=1-\frac {\delta }{6}+O(\delta ^2) \end{equation}

(Thorpe Reference Thorpe1968). Figure 12 presents the eigenfrequencies $\alpha _m^j$ as well as the eigenfunctions $W_m^j$ for $j=1,\ 2$ . Note that $W_m^{1}$ tends to the Lamb solution as $\delta \to 0$ . For finite $\delta$ , the eigenfrequency is downshifted from that of the Lamb solution by 6–10 % in the present experiments.

Figure 12.

(a) Eigenfrequencies of the first symmetric ( $j=1$ , red) and antisymmetric ( $j=2$ , blue) modes for $\delta =0.2$ (open symbols) and $\delta =0.3$ (filled symbols); asterisks represent the two-layer model. The corresponding eigenfunctions $W_n^1$ (b) and $W_n^2$ (c) for $\delta =0.2$ (solid lines) and $\delta =0.3$ (dashed lines). The experimental ranges of the forcing frequency and its second harmonic are shaded in (a), the dotted line in (b) represents the eigenfunction for the two-layer model and the vertical lines in (b,c) show the boundaries of the pycnocline.

4.3. Linear model for finite pycnocline thickness

The solution procedure follows the two-layer model approach. In this subsection, the simplified form of the potential solution (4.10) in the lower liquid ( ${\tilde {z}}\lt -\delta /2$ ) is used, neglecting terms $\varphi _{2,m}^{\prime}$ in (4.11) associated with the finite wavemaker size:

(4.22) \begin{equation} \varphi _2=\varphi _0+\frac {1}{2}\left [\sum \limits _{m=1}^\infty b_m\cos \left (\frac {m}{n}(\tilde {x}+\pi )\right )\exp \left (\frac {m}{n}{\tilde {z}}\right )\exp (-i\tilde {t})+\text{ c.c.}\right ]. \end{equation}

The vertical velocity in the upper liquid, including the pycnocline, is the solution of the Taylor–Goldstein equation (4.14):

(4.23) \begin{equation} w=\frac {1}{2}\left [\sum \limits _{m=1}^\infty a_m W_m({\tilde {z}})\cos \left (\frac {m}{n}(\tilde {x}+\pi )\right )\exp (-i\tilde {t})+\text{ c.c.}\right ]; \quad W_m\to 0,\ {\tilde {z}}\to \infty . \end{equation}

Equation (4.14) is also valid in the lower liquid, so the solution can be smoothly extended to ${\tilde {z}}\lt -\delta /2$ :

(4.24) \begin{align} W_m=\left \{ \begin{array}{ll} C_1^{m}\exp \left(\dfrac {m}{n}{\tilde {z}}\right)+C_2^{m}\exp \left(-\dfrac {m}{n}{\tilde {z}}\right)&{\tilde {z}}\lt -\delta /2\\[7pt] C_3^{m}\cos (\gamma {\tilde {z}})+C_4^{m}\sin (\gamma {\tilde {z}})&|{\tilde {z}}|\lt \delta /2\\[7pt] \exp \left(-\dfrac {m}{n}{\tilde {z}}\right)&{\tilde {z}}\gt \delta /2 \end{array} \right . \quad \gamma ^2=\left (\dfrac {2}{\alpha ^2\delta }-1\right )\dfrac {m^2}{n^2}. \end{align}

The constants $C_1^{m},\ldots, C_4^{m}$ are found from the continuity conditions of $W_m$ and $W_m^{\prime}$ at ${\tilde {z}}=\pm \delta /2$ ; the coefficient at the exponent for ${\tilde {z}}\gt \delta \gt /2$ is chosen as unity for normalisation. For frequencies different from the eigenvalues $\alpha \neq \alpha _{m}^j$ , the solution is unbounded as $z\to -\infty$ ( $C_2^{m}\neq 0$ ). Examples of functions $W_n$ are shown in figure 13(a).

Figure 13.

(a) The vertical velocity profiles $W_n$ in the upper liquid extended to ${\tilde {z}}\lt -\delta /2$ . (b) The sum of $ {{\rm d}\varPhi _n}/{{\rm d}{\tilde {z}}}$ (in black) and $b_n\exp ({\tilde {z}})$ (in blue) matches $a_nW_n$ (in red); the coefficients are defined by (4.26). Thin vertical lines denote the boundaries of the pycnocline; $\delta =0.2$ , $\alpha =\alpha _*+0.1$ .

Following Morse & Feshbach (Reference Morse and Feshbach1946), the contributions of the forcing potential $\varphi _0$ to the $m$ th spatial mode in the region between the wavemaker top and the lower boundary of the pycnocline $-h+r\lt {\tilde {z}}\lt -\delta /2$ can be represented as

(4.25) \begin{align} \varPhi _m=\int \limits _{-\pi }^{(n-1)\pi }\varphi _0(\tilde {x},{\tilde {z}})\cos \left (\frac {m}{n}(\tilde {x}+\pi )\right ){\rm d}\tilde {x}=Q_1^{m}\exp \left (\frac {m}{n}{\tilde {z}}\right )+Q_2^{m}\exp \left (-\frac {m}{n}{\tilde {z}}\right ). \end{align}

The conditions

(4.26) \begin{align} &\qquad\quad\frac {m}{n}Q_2^{m}=\frac {1}{2}n\pi C_2^{m} a_m,\nonumber\\& -\frac {m}{n}Q_1^{m}+\frac {1}{2}n\pi \frac {m}{n}b_m=\frac {1}{2}n\pi C_1^{m} a_m \end{align}

ensure matching of both vertical and horizontal velocities defined by (4.22) and (4.23) at the vertical locations $-h+r\lt {\tilde {z}}\lt -\delta /2$ (figure 13 b). The amplitudes $A_w^m$ and $A_u^m$ of the vertical and horizontal velocity components for the $m$ th mode are

(4.27) \begin{align} &A_w^m=\left |C_3^{m}a_m\right |=\left |\frac {4 m^2\left (n\gamma \cos (\gamma d/2)+m\sin (\gamma d/2)\right )}{\pi n^4\left [(m^2/n^2-\gamma ^2) \sin (\gamma \delta )+2\gamma m/n\cos (\gamma \delta )\right ] }Q_2^{m}\exp \left (\frac {m}{n}\frac {\delta }{2}\right )\right |,\notag\\ & A_u^m=\left |C_4^{m}a_m\right |=\left |\frac {4 m^2\left (m\cos (\gamma d/2)-\gamma \sin (\gamma d/2)\right )}{\pi n^4\left [(m^2/n^2-\gamma ^2) \sin (\gamma \delta )+2\gamma m/n\cos (\gamma \delta )\right ] }Q_2^{m}\exp \left (\frac {m}{n}\frac {\delta }{2}\right )\right |. \end{align}

The singularity leading to infinite velocity amplitude appears if (4.20a ) or (4.20b ) is satisfied. The amplitude of the vertical velocity at the antinode $(\tilde {x},{\tilde {z}})=(\pi, 0)$ is dominated by the contribution of the $n$ th mode, provided $\alpha$ is close to $\alpha _*$ :

(4.28) \begin{align} A_w=A_w^n=\left |\frac {4 \left (\gamma \cos (\gamma \delta /2)+\sin (\gamma \delta /2)\right )}{\pi n\left [(1-\gamma ^2) \sin (\gamma \delta )+2\gamma \cos (\gamma \delta )\right ] }Q_2^{n}\exp \frac {\delta }{2}\right |. \end{align}

The amplitude is inversely proportional to the detuning $\beta$ :

(4.29) \begin{equation} \beta =\alpha -\alpha _*, \end{equation}

since $(A_w)^{-1}=0$ for $\beta =0$ . The Taylor expansion in the vicinity of $\beta =0$ yields

(4.30) \begin{align} (A_w)^{-1}=\left |\frac {{\rm d}(A_w^n)^{-1}(\gamma _*)}{{\rm d}\gamma }\frac {{\rm d}\gamma }{{\rm d}\beta }\right | |\beta |,\quad \gamma _*=\gamma _n^1=\sqrt {\frac {2}{\alpha _*^2\delta }-1}. \end{align}

The slope of the linear dependence of the inverse amplitude on the detuning depends only on $\delta$ , since $\gamma _*$ and the corresponding eigenfrequency $\alpha _*$ are given by (4.20a ) and (4.19), respectively. Calculations show that this slope varies within 5 % for $0\lt \delta \lt 1$ . The limit for $\delta \to 0$ is

(4.31) \begin{equation} (A_w)^{-1}=\frac {n\pi }{Q_2^{n}} |\beta |. \end{equation}

4.4. Effect of dissipation

Accounting for the dissipation due to the friction at the front and back walls of the basin within the Stokes layer model can be performed similarly, since the corresponding generalisation of the Taylor–Goldstein equation (4.14) retains a similar form (see Appendix A for details):

(4.32) \begin{align} & W_m^{\prime\prime}+\left (\frac {2}{\alpha ^2 (1+s)}{\tilde {N}}^2-1\right )\frac {m^2}{n^2}W_m=0,\notag\\ &\qquad s=\sqrt {2}(1+i)\sqrt {\frac {\mu }{\rho _0 B^2}}\sim 10^{-2}, \end{align}

where $\mu$ is the dynamic viscosity. The form (4.14) is obtained from (4.32) defining the effective detuning as

(4.33) \begin{align} \beta '=\beta +\frac {s\alpha _*}{2}. \end{align}

Equation (4.31) shows that the maximum amplitude is inversely proportional to $|s|$ and corresponds to $\beta =-\alpha _*\textrm {Re}(s)/2\lt 0$ . The Stokes layer dissipation model does not account for viscous interaction between the wavemaker and the liquid; thus, for comparison with the experiments, $s$ can be treated as a free parameter of the order of $10^{-2}$ .

The dissipation in bulk is more significant for the shorter waves. The solution in the upper liquid is governed by the fourth-order equation (see Appendix A for details)

(4.34) \begin{align} & W_m^{IV}-2\frac {m^2}{n^2}W_m^{\prime\prime}+\frac {m^4}{n^4}W_m=-i{Re}\left [(1+s)W_m^{\prime\prime}+\left (\frac {2}{\alpha ^2 }{\tilde {N}}^2 -1 -s\right )\frac {m^2}{n^2}W_m\right ],\notag\\&\qquad\qquad\qquad\qquad\qquad\qquad {Re}=-\frac {\rho _0\omega }{k_n^2\mu }\sim 10^4, \end{align}

with boundary conditions

(4.35) \begin{align} &\qquad\quad W_m\to 0, \quad {\tilde {z}}\to \infty, \notag\\& W_{m}^{\prime\prime}(-\delta /2)-\frac {m^2}{n^2}W_m(-\delta /2)=0,\notag\\ & W_{m}^{\prime\prime\prime}(-\delta /2)-\frac {m^2}{n^2}W_{m}^{\prime}(-\delta /2)=0, \end{align}

that ensure vanishing of the velocity amplitude at infinity and matching to the potential solution in the lower liquid. The solution for $W_m$ is found numerically, then smoothly extended to the region ${\tilde {z}}\lt -\delta /2$ . The above-described matching procedure is then applied.

Correction to the eigenfunctions $\varphi _{2,m}$ to account for the finite wavemaker size by imposing boundary conditions (4.11) does not change the procedure significantly. The contribution of the eigenfunctions to the $m$ th spatial mode is a sum of $\exp (m{\tilde {z}}/n)$ and $\exp (-m{\tilde {z}}/n)$ . Matching of the vertical velocity and its $\tilde {z}$ derivative at any point ${\tilde {z}}_c\in (-h+r,-\delta /2)$ ensures the matching of the solutions in the entire region $-h+r\lt {\tilde {z}}\lt -\delta /2$ . The coefficients $a_m$ and $b_m$ in the decompositions (4.23) and (4.10) are found from the system of equations

(4.36) \begin{align} &\,\, \frac {1}{2}n\pi a_lW_l({\tilde {z}}_c)-\int \limits _{-\pi }^{(n-1)\pi }\left [\sum \limits _{m=1}^\infty b_m\frac {\partial \varphi _{2,m}(\tilde {x},{\tilde {z}}_c)}{\partial {\tilde {z}}} +\frac {\partial \varphi _0(\tilde {x},{\tilde {z}}_c)}{\partial {\tilde {z}}} \right ]\cos \left (\frac {l}{n}(\tilde {x}+\pi )\right ){\rm d}\tilde {x}=0,\notag\\& \frac {1}{2}n\pi a_lW_l^{\prime}({\tilde {z}}_c)-\int \limits _{-\pi }^{(n-1)\pi }\left [ \sum \limits _{m=1}^\infty b_m\frac {\partial ^2 \varphi _{2,m}(\tilde {x},{\tilde {z}}_c)}{\partial \tilde {z}^2}+\frac {\partial ^2 \varphi _{0}(\tilde {x},{\tilde {z}}_c)}{\partial \tilde {z}^2} \right ]\cos \left (\frac {l}{n}(\tilde {x}+\pi )\right ){\rm d}\tilde {x}=0,\notag\\ &\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad l=1\ldots \infty . \end{align}

For arbitrary density gradient distribution ${\tilde {N}}({\tilde {z}})$ , a boundary ${\tilde {z}}_*$ is selected so that ${\tilde {N}}({\tilde {z}})\ll 1$ for ${\tilde {z}}\lt {\tilde {z}}_*$ . The boundary conditions (4.35) are set at ${\tilde {z}}_*$ , and the matching conditions (4.36) are also applied there. The analytical inviscid results for the piecewise-constant ${\tilde {N}}({\tilde {z}})$ as well as the numerical computations based on the present model confirm that the results practically do not depend on the particular choice of ${\tilde {z}}_*$ .

4.5. Numerical results

The numerical calculations are performed for the measured Gaussian buoyancy frequency distribution (2.5); the vertical velocity distribution is given by (4.23) and the corresponding horizontal velocity is

(4.37) \begin{align} u=\frac {1}{2}\left [\sum \limits _{m=1}^\infty \frac {n}{m} a_m W_m({\tilde {z}})'\sin \left (\frac {m}{n}(\tilde {x}+\pi )\right )+\text{ c.c.}\right ]. \end{align}

Figure 14 shows an example of wavenumber amplitude spectra of vertical and horizontal velocities at ${\tilde {z}}=0$ accounting for the Stokes layer and both dissipation mechanisms discussed above. The wavenumber domains with dominant contributions to the vertical and the horizontal velocity components are clearly separated, $m=3$ and $m\approx 17$ , respectively, close to the wavenumbers that satisfy the dispersion relation for the symmetric and antisymmetric eigenfunctions. The effect of the bulk dissipation on the longer-mode amplitude is largely insignificant but is observed for the shorter modes. The dependence of the resonant amplitude modes on the detuning $\beta$ agrees with the simplified model predictions (4.31) as confirmed by the comparison of the inverse amplitude values in figure 15.

Figure 14.

Wavenumber spectra of vertical (red) and horizontal (blue) velocity components at $z=0$ at an effective resonance $\alpha =0.94$ for $\delta =0.24$ ; open symbols account for dissipation in Stokes layers only, filled symbols for bulk dissipation as well.

Figure 15.

Inverse amplitude of the vertical velocity as a function of the detuning.

The amplitudes of the antisymmetric short modes with relatively large values of $m$ depend on the detuning $\beta$ in a complex way. Three modes with largest amplitudes are shown in figure 16 for two different pycnocline dimensionless thicknesses $\delta$ ; the range of the corresponding values of $m$ shifts to the left as $\delta$ increases, in agreement with figure 12. For all modes in figure 16, a local amplitude maximum is located at $\beta \approx -\textrm {Re}(s)/2$ . The appearance of this maximum can be attributed to the fact that due to the presence of the finite-size wavemaker eigenfunctions $\varphi _{2,m}$ have complicated wavenumber spectra at any fixed vertical location $\tilde {z}$ . The response curve for the mode with the eigenfrequency closest to the resonance eigenfrequency has a single maximum. The neighbouring modes exhibit additional resonances; the corresponding frequencies increase monotonically with $m$ . The amplitudes of these resonances are proportional to the forcing at the corresponding modes (figure 11 b), and thus the maximum amplitudes of the modes with $m=n,\ 2n, \ldots$ are larger than those of the neighbouring ones.

Figure 16.

Amplitudes of horizontal velocity of dominant antisymmetric modes for (a) $\delta =0.2$ and (b) $\delta =0.3$ .

Experimental validation of the predicted properties of the wavefield structure within a finite-thickness pycnocline cannot be obtained using the shadowgraph technique described in § 3. Therefore, additional series of measurements in which the PIV technique was applied have been performed.

5.1. Procedure

The PIV set-up consists of a continuous 532 nm laser (1.5 W) with light sheet generation optics (see figure 1). To minimise the optical distortions caused by density gradient in the imaged area (Dalziel et al. Reference Dalziel, Hughes and Sutherland2000), the laser sheet is placed at a distance of 10 mm, sufficiently close to the front wall out of the Stokes layer. A 3.2 MPixel USB 3.1 camera (FLIR Blackfly) operating at 5 Hz recorded the flow field. Silver-coated hollow 10 μm glass spheres served as the tracers. LaVision Davis software was used to compute the two-dimensional instantaneous velocity field; the interrogation cell size was 48 $\times$ 48 pixels with 50 $\%$ overlap, corresponding to a spatial resolution of 0.13 mm px⁻¹. The temporal resolution defined by the camera frame rate corresponds to approximately 30 frames per wave period. The imaged area spans almost half a wavelength of the third standing-wave mode and is centred at its antinode. Supplementary movie 2 shows the maps of the vertical and horizontal velocity components normalised by the wavemaker velocity amplitude.

The experiments were carried out for amplitude of the wavemaker displacement $F$ $ =7.5$ mm, $n$ $=3$ , forcing frequency corresponding to 0.905 $\leqslant \alpha \leqslant$ 0.946 and $0.18\lt$ $\delta$ $ {\lt } 0.43$ . As in the shadowgraph experiments, the mean position of the wavemaker centre was 350 mm above the bottom of the tank. Before each series of experiments, the basin was refilled as described in § 2. To decouple the effect of frequency detuning and the pycnocline thickness $d$ , in each series of PIV experiments, the forcing frequency was kept constant while in the course of the successive runs $d$ gradually increased. At every $\alpha$ , the experimental series consisted of 12–20 runs. The first run started from rest, then the wavemaker motion was initiated and lasted for 300 cycles; after that, the wavemaker stopped and paused for 10 min during which the waves fully decayed allowing the initiation of the next run. The density profile was measured by a PCP before the initiation of the series and between the runs at the end of the wavemaker pause.

The consistency between PIV and shadowgraph measurements and the theoretical model is examined by the comparison of normalised vertical velocity amplitude distribution along the horizontal line ${\tilde {z}}=\delta$ in figure 17. The velocity amplitudes are normalised by the corresponding values at the third-mode antinode $x/\lambda _3=0.5$ ; the curves are plotted at the corresponding effective resonant conditions. A good agreement is observed, although the coverage of two experimental datasets is different.

Figure 17.

Vertical velocity amplitude along the upper pycnocline boundary ${\tilde {z}}=\delta$ , normalised by the corresponding value at $x/\lambda _3=0.5$ for PIVmeasurement, shadowgraph measurement and theory at the corresponding effective resonant conditions ( $\alpha =0.917,$ $\delta$ $=0.32$ ; $\alpha =0.931,$ $\delta$ $=0.22$ ; $\alpha =0.925,$ $\delta$ $=0.306$ , respectively).

Figure 18 presents a summary of the experimental conditions in the present experiments, along with the theoretically obtained dispersion relation $\alpha (k_3d)$ and the theoretical estimates of the effect of the finite wavemaker size and the dissipation. Each set of horizontal circles at a constant forcing frequency represents a single series of experimental runs as specified above. It demonstrates the natural widening of the pycnocline during each series. The colour of the circles corresponds to the amplitude of the vertical velocity oscillations at the antinode ( $x/\lambda _3=0.5$ ) after a quasi-steady state is attained for the given operational parameters.

Figure 18.

Diagram of the experimental conditions. The symbols represent the experimental results, with the colour showing the vertical velocity amplitude at the antinode ( $x/\lambda _3=0.5$ ) normalised by the wavemaker velocity amplitude. The lines correspond to the eigenfrequency of the rectangular basin (red), the resonant frequency for the basin with the immersed wavemaker (blue) and the effective resonant frequency accounting for dissipation (black).

5.2. Particle image velocimetry data analysis

5.2.1. Response curves

Numerous experiments carried out under diverse experimental conditions allow decoupling of the effect on the wavefield of the pycnocline thickness and of the forcing frequency. The horizontal and vertical sections of the diagram in figure 18 yield the experimentally measured variation of the system response to the forcing as a function of the dimensionless pycnocline thickness $\delta$ at fixed frequency, and of the dimensionless forcing frequency $\alpha$ at fixed pycnocline thickness. Typical examples of the resulting curves are plotted in figure 19. The system response is represented by the amplitude of the vertical velocity component at the third-mode antinode normalised by the amplitude of the wavemaker velocity $A_w(x/\lambda _3=0.5,z=0)$ . The amplitude variation with the pycnocline thickness for a fixed frequency does not change notably with $\alpha$ . The effective resonance frequency decreases with an increase in $\delta$ , while the maximum amplitude exhibits only a weak change. The normalised amplitude at the effective resonance is significantly smaller than those measured for the similarly excited surface waves (Mogilevskiy et al. Reference Mogilevskiy, Kalenko, Zemach and Shemer2024) that represent a limit case $\delta$ $\to 0$ with negligible dissipation in the pycnocline.

The viscous model yields similar results (figure 20). The downshift of the effective resonance frequency is somewhat underestimated, while the normalised amplitude at the effective resonance is nearly twice larger than in the experiments. As in the experiments (figure 19 a), the maximum of the normalised amplitude practically does not vary with $\alpha$ .

Figure 19.

The amplitude of the vertical velocity oscillations at the antinode ( $x/\lambda _3=0.5$ ) of the third-mode standing wave from PIV measurements: (a) as function of the pycnocline thickness for different forcing frequencies; (b) as function of frequency for constant pycnocline thickness $\delta$ $=0.305$ .

Figure 20.

The results of the theoretical model accounting for the full dissipation for parameters as in figure 19.

Measurements of the vertical velocity amplitude at the antinode are summarised in figure 21 as a function of the detuning $\beta =\alpha -\alpha _*(\delta )$ . Note that the figure contains about 160 data points, each representing a separate experimental run more than 10 minutes long. Both shadowgraph- and PIV-derived results collapse onto close vicinity of a single curve. The measured data agree well with the linear theory (4.31) far below the resonance. The discrepancy of the measured maximum amplitude and of the corresponding detuning cannot be attributed solely to the effect of dissipation; thus, the effect of nonlinearity may also be significant. Indeed, selection of the value of the dissipation coefficient $s$ to fit the maximum amplitude within the linear computational model does not cause a notable shift in the detuning at the effective resonance. To resolve this discrepancy by accounting for nonlinearity, advantage can be taken of close similarity between nonlinear resonant free-surface waves (Mogilevskiy et al. Reference Mogilevskiy, Kalenko, Zemach and Shemer2024) and those at the interface between the two liquids. Adjustments of their model to the present problem result in the estimation of the downshift of the effective resonance frequency $(k_nF r^2/(h-\delta ))^{2/3}\approx 0.02,$ in good agreement with experiments.

Figure 21.

Normalised amplitude of the vertical velocity at the antinode as a function of the detuning. The colour of the experimental points corresponds to the pycnocline thickness $\delta$ ; the value of dissipation coefficient $s_0$ corresponds to experiments.

5.2.2. Symmetric and antisymmetric modes in the velocity field at the forcing frequency

The series of experimental runs carried out at a constant forcing frequency $\alpha =0.917$ (marked by a rectangular box in figure 18) is now considered. To analyse the wavefield at the forcing frequency, the PIV signal was band-pass-filtered at each spatial location $(x,z)$ independently, resulting in periodic in time functions:

(5.1) \begin{align} u^1(t,x,z)=A_u(x,z)\cos (\omega t+\theta _u),\quad w^1(t,x,z)=A_w(x,z)\cos (\omega t+\theta _w). \end{align}

The contour maps of the spatial distribution of the dimensionless amplitudes and phases of the vertical $w$ and the horizontal $u$ velocity components are presented in figure 22. The distributions are dominated by the mode with $m=3$ , although the short-scale features are also present.

Figure 22.

Spatial distributions of the velocity components for $\delta$ $=0.282$ band-pass-filtered at the forcing frequency $\alpha =0.917$ . (a,b) Amplitudes $A_{w,u}$ ; (c,d) phase shift $\theta _{w,u}$ in time (in $\pi$ units) relative to the wavemaker displacement; (a,c) vertical velocity component $w$ ; (b,d) horizontal velocity component $u$ . Supplementary movie 3 shows the variation of maps of both band-pass-filtered velocity components.

The dispersion relation (figure 12) indeed allows coexistence at the forcing frequency of the dominant symmetric mode with $m=3$ and of the higher antisymmetric modes with $m$ ranging from 15to 17 (depending on the dimensionless pycnocline thickness $\delta$ ). To examine this conjecture, the third spatial mode is first identified in the recorded velocity distributions, and the assumption of the dominance of this standing-wave mode in the wavefield is verified. To this end, for each frame corresponding to instant $t$ and each $z$ independently, the distributions of both velocity components normalised by the wavemaker velocity amplitude $\omega F$ were fitted to

(5.2) \begin{align} u^{fit}(t,x,z)=u_3(t,z)\cos (k_3 x+\gamma _u(t,z)),\!\!\!\quad w^{fit}(t,x,z)=w_3(t,z)\cos (k_3 x+\gamma _w(t,z)), \end{align}

defining functions $u^{fit}(t,x,z)$ and $w^{fit}(t,x,z)$ that are periodic in time with period $T_3$ and in the horizontal direction with the period $\lambda _3$ ; $u_3(t,z)$ , $w_3(t,z)$ , $\gamma _u(t,z)$ and $\gamma _w(t,z)$ are the fitting parameters. For each location $(x,z)$ , the time-dependent periodic functions $u^{fit}$ and $w^{fit}$ have amplitudes $A_{u,w}^{fit}(x,z)$ and phase shift $\theta _{u,w}^{fit}(x,z)$ relative to the wavemaker displacement plotted in figure 23. The patterns of the phase maps correspond to the third-mode standing wave with a practically constant $\theta _{w}^{fit}(x,z)$ . The values of $\theta _{u}^{fit}(x,z)$ change sharply by $\pi$ in the vicinity of the corresponding nodes.

Figure 23.

As in figure 22 for the fitted velocity amplitudes $A^{fit}_{u,w}$ and phases $\theta ^{fit}_{u,w}$ . Supplementary movie 4 shows the temporal variation of maps of both fitted velocity components.

The profiles of the normalised vertical velocity component amplitude $A_w(x,z)$ at the antinode $x/\lambda _3=0.5$ of the band-pass-filtered record (figure 24 a) are compared with the corresponding distributions of the normalised fitted amplitudes $A_w^{fit}$ (figure 24 b) at multiple experimental runs with different values of the pycnocline thickness $\delta$ . The shapes of the velocity amplitude profiles $A_{w}(z)$ differ notably from those of the fitted distributions. The profiles obtained for different experimental conditions $A_w^{fit}$ appear to be symmetric around $z=0$ ; when normalised by their maxima, they nearly collapse, with shapes that agree well with the theoretical eigenfunction of the third mode (figure 24 c).

Figure 24.

(a) Profiles of the vertical velocity component amplitude. (b) The fitted amplitude of the third-mode wave for different pycnocline thickness $\delta$ . (c) The fitted distributions normalised by their maxima are compared with theoretical eigenfunction for $m=3$ , $\delta$ $=0.305$ and $\alpha =0.917$ (the black bold dashed line).

The contribution to the wavefield of the antisymmetric modes at the wavemaker forcing frequency is evaluated by computing at each instant the residual defined as

(5.3) \begin{align} u^{res}(t,x,z)=u^1-u^{fit},\quad w^{res}(t,x,z)=w^1-w^{fit}. \end{align}

The maps of amplitudes $A_{w,u}^{res}$ of the vertical and horizontal velocity components $w^{res}$ and $u^{res}$ presented in figure 25(a,b) exhibit a spatially quasi-periodic horizontal structure, with mean wavelength corresponding to $m$ about 16, in agreement with the dispersion relation (figure 12). Note that the amplitudes $A_w^{res}$ of the residual of the vertical velocity component $w^{res}$ practically vanish along the centreline of the pycnocline $z=0$ , where the maxima of the horizontal velocity amplitudes $A_u^{res}$ are located. The corresponding phase maps are presented in figure 25(c) and 25(d), respectively. The phase shift of $\pi$ between the two maxima of the residual of the vertical velocity amplitude confirms that the wave is antisymmetric around $z=0$ .

Figure 25.

As in figure 22 for the residual velocity amplitudes $A^{res}_{u,w}$ and phases $\theta ^{res}_{u,w}$ . Supplementary movie 5 shows the temporal variation of maps of both residual velocity components.

The estimated vertical profiles of the normalised amplitudes of the residuals of the vertical and horizontal velocity components are now compared with the theoretical eigenmodes for $m=16$ for different values of the pycnocline thickness $\delta$ . The profiles are plotted for the antinode locations of the corresponding velocity components: $x/\lambda _3=0.5$ and $x/\lambda _3=0.537$ , respectively (figure 26 a,b). The theoretical curve computed for $m=16$ and $k_3d=0.306$ is compared with the experimental profiles normalised by their maxima and agree reasonably well.

Figure 26.

Profiles of normalised amplitudes of (a) residual vertical velocity component at $x/\lambda _3=0.5$ and (b) residual horizontal velocity component at $x/\lambda _3=0.537$ for different pycnocline thickness. In (a,b), theoretical curves correspond to $m=16$ and $\delta$ $=0.306$ ; colours as in figure 24.

The adopted data processing procedure allows decomposition of the complex wavefield into coexisting different symmetric and antisymmetric spatial modes at the single forcing frequency. The wave is dominated by a single symmetric eigenmode, over which a weaker and shorter antisymmetric wave is superposed. The effective wavemode number $m$ of the antisymmetric mode was obtained from the mean distance between the consecutive maxima of the residual vertical velocity amplitude (figure 25 a). The measured values agree with the theoretically predicted values based on the dispersion relation for the range of $\alpha$ and $\delta$ adopted in the experiments (figures 12 a and 27 a). The variation of the measured normalised amplitudes of the antisymmetric mode with the pycnocline thickness $\delta$ for different constant forcing frequencies is presented in figure 27(b).

Figure 27.

(a) Normalised frequency $\alpha$ as a function of $\delta$ for antisymmetric modes with $15\leqslant m \leqslant 17$ . The markers correspond to the experimentally observed modes. (b) The response curves of the antisymmetric mode as a function of the pycnocline thickness $\delta$ for different values of $\alpha$ .

5.2.3. Second frequency harmonic

The shadowgraph experiments presented in § 3 indicate that the interface has a complex structure not just in space but also in time, with the second harmonic in the frequency that cannot be neglected (figure 8). Its excitation can be attributed to nonlinear interaction of the primary wave with itself, as well as to the temporal variation of the wavemaker immersion depth relative to the pycnocline; see (4.5b ) and boundary condition (4.6d ) that account for the variable wavemaker position.

The PIV-derived temporal variation of the vertical velocity component measured at the third-mode antinode $x/\lambda _3=0.5$ is dominated by the contributions of the two lowest-frequency harmonics $\alpha$ and $2\alpha$ . Figure 28 supports this statement for the forcing frequency $\alpha =0.917$ and different values of $\delta$ .

Figure 28.

Decomposition of the temporal record of vertical velocity into that at forcing frequency $\alpha$ and at double that frequency $2\alpha$ for (a–d) different values of $\delta$ .

The main features of the second harmonic can be predicted based on the linear dispersion relation (figure 12 a). A contour map of the vertical velocity component amplitude band-pass-filtered at double the wavemaker frequency presented in figure 29(a) is characterised by a quasi-periodic wave pattern. The vertical amplitude profile has a single maximum within the pycnocline at $z\approx 0$ which is a fingerprint of a symmetric eigenmode. The estimated mode number $m\approx 20$ is in agreement the linear theory. The two-layer model predicts a second-harmonic wave with mode number $m=12$ ; accounting for the finite thickness of the pycnocline upshifts the resonant mode number considerably. For the range of $\delta$ and $\alpha$ in the present experiments, the slope of the dispersion curve in figure 12(a) at the vicinity of $m$ corresponding to the second frequency harmonic is quite small. This enables coexistence in the wavefield of multiple neighbouring resonant standing modes, resulting in complex measured spatial distributions of the vertical velocity component amplitudes (figure 29 b).

Figure 29.

(a) Contour map of vertical velocity amplitude of the second harmonic $A_w^{(2)}$ (normalised by the wavemaker velocity amplitude) for $\delta$ $=0.342$ and $\alpha =0.903$ . (b) The spatial distribution of the normalised amplitude profiles of the second harmonic of the vertical velocity component at $z=0$ , for different $\delta$ . Supplementary movie 6 shows the temporal variation of the map of the band-pass-filtered vertical component.