2.1. The linear dispersion relation

We now study the behaviour and dispersion relation of small-amplitude solutions. When the amplitudes of both the free surface and vorticity interface are zero, an exact solution is known to the boundary-value problem (2.1), given by $\zeta _1=H$ , $\zeta _2=0$ and $\phi _1=\phi _2=\psi _1=\psi _2=0$ . Thus, in the study of small-amplitude solutions we consider

(2.4) \begin{equation} \zeta _1 \sim H+\epsilon \bar {\zeta }_1, \quad \zeta _2 \sim \epsilon \bar {\zeta }_2 , \quad \phi _1 \sim \epsilon \bar {\phi }_1, \quad \phi _2 \sim \epsilon \bar {\phi }_2, \quad \psi _1 \sim \epsilon \bar {\psi }_1, \quad \psi _2 \sim \epsilon \bar {\psi }_2, \end{equation}

which hold as $\epsilon \to 0$ . Further, the Bernoulli constant is expanded as $B_1 \sim H+ \epsilon \bar {B}_1$ . In order to solve for the perturbation quantities in (2.4) above, we first have to obtain the governing equations at $O(\epsilon )$ . Note that the boundary conditions are evaluated at $y=\zeta _1(x,t)$ and $y=\zeta _2(x,t)$ , which are also expanded in $\epsilon$ ; these conditions are resolved by Taylor expanding each term.

The equations obtained at $O(\epsilon )$ are ${\nabla} ^2 \bar {\phi _1}=0$ for $0 \leqslant y \leqslant H$ , and ${\nabla} ^2 \bar {\phi }_2=0$ for $-\infty \lt y \leqslant 0$ . The far-field conditions are $(\bar {\phi }_2,\bar {\psi }_2) \to (0,0)$ as $y \to - \infty$ . The boundary conditions holding at $y=0$ are $\bar {\zeta }_{2t}=\bar {\phi }_{2y}$ , $\bar {\phi }_{2y}=\bar {\phi }_{1y}$ , $\bar {\phi }_{2x}=\bar {\phi }_{1x}+ \bar {\zeta }_2$ and $\bar {\psi }_2=\bar {\psi }_1$ . Those holding at $y=H$ are $\bar {\zeta }_{1t}=\bar {\phi }_{1y}-H\bar {\zeta }_{1x}$ and $\bar {\phi }_{1t}+H\bar {\phi }_{1x}-\bar {\psi }_1+\bar {\zeta }_1=\bar {B}_1$ . The solutions to these can be written as

(2.5a) \begin{align} \{ \bar {\zeta }_1 , \bar {\zeta }_2 \} &= \{ \alpha , \beta \}\mathrm{e}^{ {i}k(x - ct)}+\text{c.c.}, \\[-12pt]\nonumber \end{align}

(2.5b) \begin{align} \{ \bar {\phi }_1 , \bar {\psi }_1 \} &= \big [ \{ \mu , \mathrm{i}\mu \}\mathrm{e}^{k y}+\{ \nu , -\mathrm{i}\nu \}\mathrm{e}^{-k y}\big ]\mathrm{e}^{{i}k(x - ct)} +\text{c.c.}, \\[-12pt]\nonumber \end{align}

(2.5c) \begin{align} \{ \bar {\phi }_2 , \bar {\psi }_2 \} &=\{ \gamma , \mathrm{i}\gamma \}\mathrm{e}^{ {i}k(x - ct)}\mathrm{e}^{k y}+\text{c.c.}, \end{align}

where $\text{c.c.}$ denotes complex conjugate and $k=2 \pi n /L$ is the wavenumber ( $n$ is a positive integer). The solutions in (2.5b ) and (2.5c ) satisfy Laplace’s equation, the Cauchy–Riemann equations applied to the analytic functions $\bar {\phi }_1+\mathrm{i}\bar {\psi }_1$ and $\bar {\phi }_2+\mathrm{i}\bar {\psi }_2$ , as well as the decay conditions $(\bar {\phi }_2,\bar {\psi }_2) \to (0,0)$ . In (2.5), $c \in \mathbb{C}$ is the complex-valued unknown whose real and imaginary parts are the speed and growth rate of the perturbations, respectively.

Substituting solutions (2.5) into the $O(\epsilon )$ equations yields the algebraic equations

(2.6) \begin{align} &\gamma =\mu -\nu , \qquad \mathrm{i}k\gamma =\beta + \mathrm{i}k(\mu +\nu ), \qquad \gamma =-\mathrm{i}c\beta ,\nonumber \\ &\mathrm{i}(H-c)\alpha =\mu \mathrm{e}^{k H}-\nu \mathrm{e}^{-k H} \qquad \textrm {and}\nonumber \\ &\mathrm{i}k(H-c)[\mu \mathrm{e}^{kH} + \nu \mathrm{e}^{-kH}]-\mathrm{i}[\mu \mathrm{e}^{kH} - \nu \mathrm{e}^{-kH}]+\alpha =0. \end{align}

In solving this set of five equations for the five unknown constants, non-trivial solutions exist only when the dispersion relation

(2.7) \begin{equation} (H-c+1)\big (2 kc-1+\mathrm{e}^{-2kH}\big ) =k(H-c)^2 \big (2k c-1-\mathrm{e}^{-2kH}\big ) \end{equation}

is satisfied. This is a cubic equation for the unknown $c$ , in terms of $H$ and $k=2 \pi n/L$ . The real part, $\text{Re}[c]$ , corresponds to the wavespeed of time-dependent perturbations to the flat wave state, and the imaginary part, $\text{Im}[c]$ , to the growth rate of the perturbations.

Figure 2.

Solutions to the cubic linear dispersion relation (2.7) are shown for $H=1$ in $(a)$ , and for $L=5$ in $(b)$ . Values of Re $[c]$ are shown in black and values for Im $[c]$ in grey. In each figure there is a range of parameters for which unstable perturbations exist with Im $[c]\gt 0$ , and the real parts of the wavespeed in these cases are shown with dotted lines. This occurs when two solutions of the cubic have the same real part and are complex conjugates of one another. The points of marginal stability are indicated with circles, and these are studied further across the parameter space in figure 3.

Solutions of the cubic dispersion relation (2.7) are shown in figure 2 for $n=1$ , firstly for $H=1$ in figure 2 $(a)$ and for $L=5$ in figure 2 $(b)$ . In each case, we see that there is a range of the remaining parameter ( $L$ or $H$ ) for which the time-dependent perturbation is unstable. For instance, in figure 2 $(a)$ , which has $H=1$ , there are three marginally stable modes (solutions with Im $[c]=0$ ) for $0\leqslant L \leqslant 2.1589$ and $L \geqslant 2.5814$ . When $2.1589\lt L\lt 2.5814$ , there is one marginally stable mode and two complex-conjugate modes with $\text{Im}[c] \neq 0$ ; these correspond to stable and unstable modes with the same value of the wavespeed, $\text{Re}[c]$ .

Figure 3.

Values of $H$ and $L$ for which small-amplitude perturbations to the base state, $\zeta _1(x,t)=H$ and $\zeta _2(x,t)=0$ , are marginally stable are shown. This corresponds to two repeated real-valued roots to the cubic dispersion relation (2.7) for $c$ . We have taken $n=1$ which results in one wave period within the periodic domain. Black curves show the edge of the marginally stable parameter region, within which the linear solutions are unstable. Black circles correspond to the four locations of marginal stability seen in figure 2. This figure is comparable to figure 1 from Caponp et al. (Reference Caponp, Yuen, Milinazzo and Saffman1991).

The values of $H$ and $L$ that correspond to the edge of the region of instability may be determined by further studying the cubic dispersion relation (2.7). At the edge of the instability region, the complex-conjugate roots have $\text{Im}[c]=0$ and therefore are repeated roots. In writing this dispersion relation in the form $c^3 + a_2 c^2 + a_1 c +a_0=0$ , where $a_0$ , $a_1$ and $a_2$ are known functions of $L$ and $H$ , we identify a repeated root by factorising the cubic as $(c-\alpha )^2(c-\beta )=0$ , where $c=\alpha$ is the repeated root. Equating these two representations of the cubic dispersion relation yields three constraints for the two unknowns $\alpha$ and $\beta$ . Two of these constraints determine $\alpha =-a (a_2 + \sqrt {a_2^2-3a_1} )/3$ and $\beta =-a_2-2\alpha$ , and the remaining constraint requires $\alpha ^2 (a_2 + 2\alpha )-a_0=0$ , which forms an implicit relation for $L$ and $H$ . We solve this constraint numerically by Newton iteration, the results of which are shown in figure 3 with $n=1$ . Two branches of parameter values exist as shown in this figure. For every value of $H$ , there exists a range of values of $L$ for which the time-dependent perturbation of the base state is unstable.

We can also study the three solutions of dispersion relation (2.7) under the limits of $k \to 0$ and $k \to \infty$ . Firstly, for the long-wave limit of $k \to 0$ we obtain

(2.8a) \begin{equation} c \sim H -H^2 k \quad \textrm {and}\quad c \sim \pm k^{-1/2} + \frac {H^2}{2}k, \end{equation}

and for short waves with $k \to \infty$ we obtain

(2.8b) \begin{equation} c \sim \frac {1}{2}k^{-1} \quad \textrm {and}\quad c \sim H \pm k^{-1/2}. \end{equation}

In each case, the two latter modes correspond to gravity waves on infinite depth in the presence of a mean current.

2.2. Travelling-wave solutions of small amplitude

When the wavespeed $c$ is real valued, the small-amplitude disturbances studied in § 2.1 are travelling waves. The corresponding solutions for each interface are then given by

(2.9a) \begin{align} \zeta _1(x,t) &\sim H + A \cos {\big (k (x-ct)\big )}, \\[-12pt]\nonumber \end{align}

(2.9b) \begin{align} \zeta _2(x,t) &\sim -\frac {2k(H-c)\mathrm{e}^{-kH}}{(2kc-1)+\mathrm{e}^{-2kH}} A\cos {\big ( k (x-ct)\big )}, \end{align}

where $0 \leqslant A \ll 1$ , and $k=2 \pi n/L$ . We have considered the case for which $A$ is non-negative without any loss of generality – the case with negative $A$ is captured automatically due to translational invariance. Recall that $\zeta _1(x,t)$ is the free surface that bounds the rotational fluid with a region of uniform pressure, and $\zeta _2(x,t)$ is the vorticity interface between the two fluids. Two interesting cases emerge in the linear solutions (2.9). Firstly, if $c =c_1=(1- \mathrm{e}^{-2kH})/(2k)$ , then the solution for $\zeta _1(x,t)$ in (2.9a ) is constant, $\zeta _1=H$ . This corresponds to a travelling wave on the vorticity interface when the top free surface is replaced by a rigid boundary, which was studied by Rayleigh (Reference Rayleigh1894). Therefore, $c=c_1$ represents the linear speed of the shear wave mode. Secondly, if $c=c_2=H$ , then the solution for $\zeta _2(x,t)$ in (2.9b ) is identically zero. This corresponds to a linear surface wave and a zero-amplitude vorticity interface.

Figure 4.

Branches of travelling-wave solutions are shown with solid lines in the $(H,c)$ plane for fixed $L=5$ , $n=1$ and $A=10^{-5}$ from the amplitude condition (2.10). The condition $c_1=(1-\mathrm{e}^{-2kH})/(2k)$ which results in $\zeta _1(x)=0$ is shown dotted, and $c_2=H$ which results in $\zeta _2(x)=0$ is shown dashed. Panels $(a)$ to $(i)$ show the nine labelled solution profiles. In these, the free surface $\zeta _1(x)$ is shown in black, and the vorticity interface $\zeta _2(x)$ is shown dashed. Note that we have plotted $\zeta _1-H$ in order to see the surface profile more clearly. The upper branch with the largest value of $c$ corresponds to the solutions calculated by Milinazzo & Saffman (Reference Milinazzo and Saffman1990).

Figure 4. Long description

The image contains a combination of line graphs and solution profiles. Panel A shows a line graph with the horizontal axis labeled H and the vertical axis labeled c. The graph includes solid, dotted, and dashed lines representing different conditions for the wavespeed. Key points labeled a through i mark specific solution profiles. Panels B to I display the nine labeled solution profiles, each showing the free surface in black and the vorticity interface in dashed lines. The horizontal axis in these panels is labeled x over L, and the vertical axis is labeled y over L. Each profile illustrates the surface and vorticity interface displacements for different parameter values. The upper branch with the largest value of c corresponds to the solutions calculated by Milinazzo & Saffman (1990).

Notice that the solutions with $c_1=(1- \mathrm{e}^{-2kH})/(2k)$ and $c_2 =H$ are obtained as possible solutions from the linear dispersion (2.7) relation only in the limits of $H \to \infty$ and $H \to 0$ with $k \to \infty$ such that $kH$ is finite. These values of $c$ are shown in the $(c,H)$ plane in figure 4 for $n=1$ and $L=5$ , where the solid lines are solutions of the dispersion relation with $\text{Im}[c]=0$ . Here, the values of $c_1$ are indicated by the dotted line and those of $c_2$ by the dashed line. By expanding the dispersion relation as $H \to \infty$ , the two solution branches parallel to $c_2=H$ as $H \to \infty$ in figure 4 may be obtained as $c \sim H -1/(2k) \pm \sqrt {(1+4k)/(4k^2)} +O(H^{-1})$ , and therefore represent two gravity modes on an infinite-depth fluid with constant vorticity.

As a measure of the solution amplitude, we employ the sum of individual amplitudes defined by

(2.10) \begin{equation} A=\frac {\text{max}(\zeta _1)-\text{min}(\zeta _1) + \text{max}(\zeta _2)-\text{min}(\zeta _2)}{L}, \end{equation}

which is positive definite and tends to zero in the small-amplitude limit of $\epsilon \to 0$ from (2.4). Note that there are many ways to define an amplitude condition for these solutions. For instance, in the study by Milinazzo & Saffman (Reference Milinazzo and Saffman1990), the condition $\text{max}(\zeta _1)-\text{min}(\zeta _1)$ was used, but this has the issue of permitting highly nonlinear vorticity interfaces with $\zeta _2 =O(1)$ in the limit of small amplitude. Alternative choices are the momentum and energy derived in Appendix A.2, which are conserved quantities of the time-dependent system. However, the value of these is finite and non-zero in the small-amplitude limit; moreover when shifted by a constant to ensure that the momentum and energy tend to zero in the small-amplitude limit, they fail to be positive definite.

Figure 5.

Branches of linear solutions are shown in the $(H,c)$ plane for $L=5$ and $A=10^{-5}$ . The three branches with $n=1$ were shown in figure 4. Here, $n$ represents the number of wave periods within the periodic domain. Also shown are the solution branches with $n=2$ (dashed) and $n=3$ (dotted). Solution branches with different wavenumbers intersect with one another, which leads to the bifurcations of nonlinear solutions studied in § 4.

Small-amplitude solutions are shown across each of the three solution branches in figure 4 for $L=5$ . These solutions were calculated with an amplitude of $A=10^{-5}$ from (2.10). Along the solution branch with the largest value of the wavespeed $c$ corresponding to the gravity-wave mode travelling in the same direction as the shear flow, the amplitude of the free surface $\zeta _1(x)$ is large and the amplitude of the vorticity interface $\zeta _2(x)$ is small. This is seen in profiles figure 4(a–c) here, as $H \to \infty$ along this branch we have $\zeta _2 \to 0$ . The remaining two solution branches consist of solutions with a larger-amplitude vorticity interface. The branch starting from $c=0$ corresponds to the shear wave mode while the remaining branch represents the gravity-wave mode travelling in the opposite direction to the shear flow. These two modes collide as the value of $H$ increases and subsequently lose their stability. Along the solution branch containing $d$ , $e$ and $f$ , or before the collision occurs, the free surface and vorticity interface have the same polarity, whereas along branch $g$ , $h$ , $i$ this polarity is opposite to one another. Along these two branches, we have that the $\zeta _1(x) \to 0$ as $H \to 0$ starting from solution $d$ , and also as $H \to \infty$ starting from solution $i$ . Further, on the third solution branch we have $\zeta _2(x) \to 0$ as $H \to \infty$ when starting from solution $g$ .

All solutions shown in figure 4 have $n=1$ (which enters into the wavenumber $k=2 \pi n/L$ ), and so have one wave period within the periodic domain. The linear dispersion relation (2.7) involves $k$ , and so branches computed with different wavenumbers have different values of the wavespeed $c$ . The linear branches of solutions with $n=1$ , $n=2$ and $n=3$ are shown in figure 5 for $L=5$ . We see that solution branches with different wavenumbers intersect with one another at multiple places. This means that at certain parameter values, there exist linear modes with different wavenumbers that travel at the same speed. In general, if it exists, this so-called resonance between gravity and shear wave modes with wavenumbers $k_n=2\pi n/L$ and $k_m=2\pi m/L$ can occur at two wavespeeds determined from (2.7) by

(2.11) \begin{equation} \frac {\big (2 k_nc-1+\mathrm{e}^{-2k_nH}\big )}{n \big (2k_n c-1-\mathrm{e}^{-2k_nH}\big )}= \frac {\big (2 k_mc-1+\mathrm{e}^{-2k_mH}\big )}{m \big (2k_m c-1-\mathrm{e}^{-2k_mH}\big )}\,. \end{equation}

Analogous behaviour arises for gravity-capillary waves on infinite depth as first shown by Wilton (Reference Wilton1915). A weakly-nonlinear theory was developed by Chen & Saffman (Reference Chen and Saffman1979) to resolve this issue for gravity-capillary waves, and it was shown that when linear dispersion curves cross there emerge bifurcations between these, which was later confirmed in a numerical work by Chen & Saffman (Reference Chen and Saffman1980). The same phenomenon emerges in our model formulation due to the crossing of linear solution branches. A numerical method to study this phenomenon for fully nonlinear solutions is now developed.

3.1. The time-dependent conformal mapping

Figure 6.

The time-dependent mapping from the physical $(x,y)$ plane to the conformal $(\xi ,\eta )$ plane is shown. Each fluid region, $-\infty \lt y \leqslant \zeta _2$ and $\zeta _2 \leqslant y \leqslant \zeta _1$ , involves a distinct conformal map. The function $\gamma (\xi ,t)$ is introduced in § 3.2 to enforce continuity of physical quantities across each mapped region at $\eta =0$ .

Figure 6. Long description

A diagram illustrating the time-dependent mapping from the physical plane to the conformal plane for two fluid regions. Panel A: The left side of the diagram shows the physical plane with two distinct fluid regions, labeled as y = ξ1(x, t) and y = ξ2(x, t). The regions are depicted with curved lines representing the fluid surfaces at different times. The x and y axes are labeled, indicating the spatial dimensions. Panel B: The right side of the diagram shows the conformal plane with a grid representing the mapped regions. The η and ξ axes are labeled, indicating the conformal dimensions. The top boundary is labeled as η = h(t), and the bottom boundary is labeled as η = 0. Arrows indicate the mapping functions x(ξ1, η1, t), y(ξ1, η1, t), x(ξ2, η2, t), and y(ξ2, η2, t), showing the transformation from the physical plane to the conformal plane for each fluid region.

A different conformal map must be used for each of the two fluid layers. For the irrotational fluid confined within $-\infty \lt y \leqslant \zeta _2(x,t)$ , we consider the map that takes $(x,y)$ to $(\xi _2,\eta _2)$ , such that the fluid is located in the range $-\infty \lt \eta _2 \leqslant 0$ . For the rotational fluid located in the range $\zeta _2(x,t) \leqslant y \leqslant \zeta _1(x,t)$ , we consider the map that takes $(x,y)$ to $(\xi _1,\eta _1)$ , such that the fluid is located in the range $0 \leqslant \eta _1 \leqslant h(t)$ . This mapping is depicted in figure 6. The conformal method developed in this section was originally derived for a single bounded fluid layer (Choi & Camassa Reference Choi and Camassa1999) and was extended to the case of two unbounded fluid layers by Murashige & Choi (Reference Murashige and Choi2022). Our implementation requires the combination of an unbounded lower layer for the irrotational fluid, and a bounded upper region for the vorticity layer.

The solutions in the lower fluid layer are then given by $x_2(\xi _2,\eta _2,t)$ , $y_2(\xi _2,\eta _2,t)$ , $\phi _2(x_2(\xi _2,\eta _2,t),y_2(\xi _2,\eta _2,t),t)$ and $\psi _2(x_2(\xi _2,\eta _2,t),y_2(\xi _2,\eta _2,t),t)$ . Each of these satisfies Laplace’s equation within the conformal domain (i.e. ${\nabla} ^2 y_2 = \partial ^2_{\xi _2} y_2 + \partial ^2_{\eta _2} y_2=0$ ). Further, we have the decay conditions $x_2\sim \xi _2$ , $y_2 \sim \eta _2$ , $\phi _2 \to 0$ and $\psi _2 \to 0$ as $\eta _2 \to -\infty$ , the first two of which are obtained from requirements of the map, and the latter two from (2.1g ). Since the functions $x_2+\mathrm{i}y_2$ and $\phi _2 + \mathrm{i}\psi _2$ are harmonic, it suffices to solve only for the imaginary part of each and then use the Cauchy–Riemann equations to obtain the solution for the real part. The solutions of these conformal boundary-value problems are then given by

(3.1a) \begin{align} \{ x_2,y_2 \} &=\{\xi _2 , \eta _2\}+\sum _{\substack {m=-\infty \\ m \neq 0}}^{\infty } \big \{ \! - \mathrm{i}\boldsymbol{\cdot }\text{sgn}(m)a_m(t), a_m(t)\big \}\mathrm{e}^{ |m|k \eta _2}\mathrm{e}^{{i} m k \xi _2}, \\[-12pt]\nonumber \end{align}

(3.1b) \begin{align} \{ \phi _2,\psi _2 \} &=\sum _{\substack {m=-\infty \\ m \neq 0}}^{\infty } \big \{ \! - \mathrm{i}\boldsymbol{\cdot }\text{sgn}(m)b_m(t),b_m(t)\big \}\mathrm{e}^{|m|k \eta _2}\mathrm{e}^{{i} m k \xi _2}, \end{align}

where $\text{sgn}(m)$ is the signum function and $k=2 \pi n/L$ .

In the upper fluid layer, we denote the solutions by $x_1(\xi _1,\eta _1,t)$ , $y_1(\xi _1,\eta _1,t)$ , $\phi _1(x_1(\xi _1,\eta _1,t),y_1(\xi _1,\eta _1,t),t)$ and $\psi _1(x_1(\xi _1,\eta _1,t),y_1(\xi _1,\eta _1,t),t)$ . These satisfy Laplace’s equation for $0 \leqslant \eta _1 \leqslant h(t)$ , as well as matching constraints with the solutions (3.1) of the lower fluid layer at $\eta _1=0$ . A difficulty that arises with this matching is that the horizontal conformal coordinates are not equal at the vorticity interface $\eta _2=0$ and $\eta _1=0$ , and require the introduction of an intermediary variable (which itself is an unknown of the problem). This process is detailed further in § 3.2. The solutions in the upper fluid layer may be obtained as

(3.2a) \begin{align} y_1 &=\eta _1+c_0(t) +\sum _{\substack {m=-\infty \\ m \neq 0}}^{\infty } \big [c_m(t)\mathrm{e}^{m k \eta _1}+d_m(t)\mathrm{e}^{-mk \eta _1}\big ]\mathrm{e}^{{i} m k \xi _1}, \\[-12pt]\nonumber \end{align}

(3.2b) \begin{align} \phi _1 &= e_0(t) \xi _1 +f_0(t) \eta _1+g_0(t)+\sum _{\substack {m=-\infty \\ m \neq 0}}^{\infty }\big [e_m(t)\mathrm{e}^{m k \eta _1}+f_m(t)\mathrm{e}^{-m k \eta _1}\big ]\mathrm{e}^{{i} m k \xi _1}. \end{align}

Solutions for $x_1$ and $\psi _1$ can be obtained by integrating the Cauchy–Riemann equations $x_{1\xi _1}=y_{1\eta _1}$ and $\psi _{1 \xi _1}=-\phi _{1 \eta _1}$ , with the constant of integration for $x_1$ determined by the condition $x_1(0,0,t)=0$ .

In this formulation, we restrict solutions (3.1) in the lower fluid layer to $\eta _2=0$ , and (3.2) solutions in the upper fluid layer to both $\eta _1=0$ and $\eta _1=h(t)$ . The solutions on the boundaries are represented with capital letters, for instance those for the vertical coordinates $y_1$ and $y_2$ are

(3.3a) \begin{align} Y_{1U}(\xi _1,t)=y_1(\xi _1,h(t),t), \quad Y_{1L}(\xi _1,t)=y_1(\xi _1,0,t), \quad Y_{2U}(\xi _2,t)=y_2(\xi _2,0,t). \end{align}

The subscript number refers to the original fluid layer, while the subscript letter denotes whether the solution has been evaluated at the upper or lower surface that bounds the respective layer. For instance, $Y_{1U}$ denotes the solution $y_1(\xi _1,\eta _1,t)$ evaluated at the upper interface $\eta _1=h(t)$ of layer one. Equivalent definitions are also introduced for $X_{1U}$ , $X_{1L}$ , $X_{2U}$ , $\varPhi _{1U}$ , $\varPhi _{1L}$ , $\varPhi _{2U}$ , $\varPsi _{1U}$ , $\varPsi _{1L}$ and $\varPsi _{2U}$ .

The quantities for the lower fluid layer are related by the Hilbert transform (see below), giving

(3.4) \begin{equation} X_{2U}=\xi _2 - \mathscr{H}\,\,[Y_{2U}], \quad \varPhi _{2U}=-\mathscr{H}\,\,[\varPsi _{2U}]. \end{equation}

The harmonic relations (3.4) relate the real and imaginary parts of an analytic function that decays to zero as $\eta _2 \to -\infty$ . Note that the relation $\varPsi _{2U}=\mathscr{H}\,\,[\varPhi _{2U}]$ is also required for the time-dependent mapping, in which an evolution equation is first obtained for $\varPhi _{2U}$ . The operator $\mathscr{H}\,$ appearing in (3.4) is the Hilbert transform, defined by

(3.5) \begin{equation} \mathscr{H}\,\,[Y] =\frac {1}{\pi }\,{\int \!\!\!\!\!\!{-}}_{-\infty }^{\infty }\frac {Y(\tilde {\xi })}{\tilde {\xi }-\xi } \mathrm{d}{\tilde {\xi }}. \end{equation}

Analogous harmonic relations are obtained for the quantities for the upper fluid layer, which are given by

(3.6) \begin{equation} X_{1L}=\xi _1 + \mathscr{T}\,\,[Y_{1L}] -\mathscr{S}\,\,[Y_{1U}], \quad X_{1U}=\xi _1+\mathscr{S}\,\,[Y_{1L}]-\mathscr{T}\,\,[Y_{1U}], \end{equation}

(3.7) \begin{equation} \varPhi _{1L}^{\prime }=A_0 + \mathscr{T}\,\,[\varPsi _{1L}^{\prime }] -\mathscr{S}\,\,[\varPsi _{1U}^{\prime }], \quad \varPhi _{1U}^{\prime }=A_0+\mathscr{S}\,\,[\varPsi _{1L}^{\prime }]-\mathscr{T}\,\,[\varPsi _{1U}^{\prime }], \end{equation}

where $A_0=(L h)^{-1}\int _{-L/2}^{L/2} [\varPsi _{1U}-\varPsi _{1L}] \mathrm{d}\xi _2$ . The operators $\mathscr{T}$ and $\mathscr{S}$ are principal-valued integrals defined by

(3.8) \begin{equation} \mathscr{T}\,\,[Y]={\int \!\!\!\!\!\!{-}}_{-\infty }^{\infty }\frac {Y(\tilde {\xi })}{2h} \coth {\bigg [\frac {\pi (\tilde {\xi }-\xi )}{2h}\bigg ]} \mathrm{d}{\tilde {\xi }}, \quad \mathscr{S}\,\,[Y]={\int \!\!\!\!\!\!{-}}_{-\infty }^{\infty }\frac {Y(\tilde {\xi })}{2h} \tanh {\bigg [ \frac {\pi (\tilde {\xi }-\xi )}{2h}\bigg ]} \mathrm{d}{\tilde {\xi }}, \end{equation}

where $h=L^{-1}\int _{-L/2}^{L/2}(Y_{1U}-Y_{1L})\mathrm{d}\xi _1$ . In practise, we calculate each of the three principal-valued integrals spectrally with the Fourier multipliers detailed in § 3.5.

By using the chain rule, expressions can be obtained for partial derivatives with respect to $x$ and $y$ in terms of the conformal coordinates $\xi _1$ and $\eta _1$ , or $\xi _2$ and $\eta _2$ . The Cauchy–Riemann equations are then used to eliminate derivatives with respect to $\eta _1$ and $\eta _2$ . This results in expressions which only involve differentiation with respect to $\xi _1$ or $\xi _2$ , which are used to convert the two-dimensional boundary-value problem (2.1) into a one-dimensional boundary-integral formulation. We use primes $(\prime )$ to denote differentiation with respect to either $\xi _1$ or $\xi _2$ , depending on the dependence of the function. For instance, $Y_{1U}^{\prime }$ denotes $\partial _{\xi _1}Y_{1U}$ . We then obtain $\zeta _{1x}=Y_{1U}^{\prime }/X_{1U}^{\prime }$ , $\zeta _{2x}=Y_{1L}^{\prime }/X_{1L}^{\prime }$ and $\zeta _{2x}=Y_{2U}^{\prime }/X_{2U}^{\prime }$ for the derivatives of each interface. Derivatives of the velocity potential are found to be $\phi _{1x}=(Y_{1U}^{\prime } \varPsi _{1U}^{\prime } + X_{1U}^{\prime } \varPhi _{1U}^{\prime })/J_{1U}$ and $\phi _{1y}=(Y_{1U}^{\prime } \varPhi _{1U}^{\prime }- X_{1U}^{\prime } \varPsi _{1U}^{\prime })/J_{1U}$ , with equivalent results holding for $\phi _{1x}$ , $\phi _{1y}$ , $\phi _{2x}$ and $\phi _{2y}$ evaluated at the lower interface, $\eta _1=\eta _2=0$ . Here, we have defined $J_{1U}(\xi _1,t)=(X_{1U}^{\prime })^2+(Y_{1U}^{\prime })^2$ and $J_{2L}(\xi _2,t)=(X_{2L}^{\prime })^2+(Y^{\prime }_{2L})^2$ , which are the Jacobians of each conformal map.

3.2. Matching conditions at the vorticity interface

At the vorticity interface $y=\zeta _2(x,t)$ , which is specified in conformal coordinates by both $\eta _1=0$ and $\eta _2=0$ , there are four matching conditions that must hold. These are $x_2=x_1$ , $y_2=y_1$ and continuity of horizontal and vertical velocities across the interface. These latter two conditions are given by $\phi _{2x}=\phi _{1x}+y_1$ and $\phi _{2y}=\phi _{1y}$ from (2.1f ). Only three of these have to be enforced on solutions of the system, from which the other conditions are automatically satisfied. These three conditions are

(3.9a) \begin{align} X_{2U}(\hat {\xi }+\gamma (\hat {\xi },t),t) = X_{1L}(\hat {\xi }-\gamma (\hat {\xi },t),t), \\[-12pt]\nonumber \end{align}

(3.9b) \begin{align} Y_{2U}(\hat {\xi }+\gamma (\hat {\xi },t),t) = Y_{1L}(\hat {\xi }-\gamma (\hat {\xi },t),t), \end{align}

as well as continuity of the horizontal velocity, which when written in terms of the conformal solutions using the chain rule yields

(3.9c) \begin{equation} \frac {Y_{2U}^{\prime } \varPsi _{2U}^{\prime } + X_{2U}^{\prime } \varPhi _{2U}^{\prime }}{J_{2U}}=\frac {Y_{1L}^{\prime } \varPsi _{1L}^{\prime } + X_{1L}^{\prime } \varPhi _{1L}^{\prime }}{J_{1L}}+Y_{1L}. \end{equation}

In (3.9c ) above, the left-hand side is evaluated at $\xi _2=\hat {\xi }+\gamma (\hat {\xi },t)$ , and the right-hand side at $\xi _1=\hat {\xi }-\gamma (\hat {\xi },t)$ . Primes $(\prime )$ denote differentiation with respect to the conformal coordinate $\xi _1$ or $\xi _2$ . Further, the unknown function $\gamma (\hat {\xi },t)$ , introduced from taking $\xi _1=\hat {\xi }-\gamma (\hat {\xi },t)$ and $\xi _2=\hat {\xi }+\gamma (\hat {\xi },t)$ , satisfies $\gamma (\pm L/2,t)=0$ . An analogous method was used by Murashige & Choi (Reference Murashige and Choi2022) to study the nonlinear interface problem between two immiscible fluids with different densities, following Saffman & Yuen (Reference Saffman and Yuen1982). Note that in formulations where neither the horizontal or vertical components of the velocity field are continuous across the interface, such as for a density jump, the normal component of the velocity field would be continuous, but the tangential component would be discontinuous.

3.3. Conserved quantities

There are four conserved quantities of our model formulation: mass of the upper rotational fluid layer, mass of the lower irrotational fluid layer, horizontal momentum and energy. Their conservation laws are derived in Appendix A.2 from solutions of the dimensional two-dimensional boundary-value problem, and in the same section are reduced to a one-dimensional form involving only quantities evaluated at each interface. We now express these in terms of the conformal coordinates, $\xi _1$ and $\xi _2$ . Firstly, the mass conditions (2.2) are written in conformal coordinates by using the chain rule at the lower (internal) interface to find $\mathrm{d}x=X_{2U}^{\prime } \mathrm{d}\xi _2$ , and at the upper free surface to find $\mathrm{d}x=X_{1U}^{\prime } \mathrm{d}\xi _1$ , where primes $(\prime )$ denote partial derivatives with respect to the spatial variable. This yields

(3.10a,b) \begin{equation} \frac {1}{L}\int _{-L/2}^{L/2}Y_{1U}X_{1U}^{\prime }\mathrm{d}\xi _1=H, \quad \int _{-L/2}^{L/2}Y_{2U}X_{2U}^{\prime }\mathrm{d}\xi _2=0. \end{equation}

Next, we consider the momentum $P$ and energy $E$ of the system given as boundary-integral expressions involving integration over $x$ from (A8) and (A11). Upon non-dimensionalisation and writing in conformal coordinates, the expressions

(3.11a) \begin{equation} P=\int _{-L/2}^{L/2} \bigg [\! - Y_{1U}\varPsi _{1U}^{\prime } +\frac {1}{2}(Y_{1U}^2X_{1U}^{\prime }+Y_{1L}^2X_{1L}^{\prime })\Bigg ] \mathrm{d}\xi _1, \end{equation}

(3.11b) \begin{equation} \begin{aligned} E&= \int _{-L/2}^{L/2} \big [(\varPsi _{1U}+Y_{1U}^2)\varPhi _{1U}^{\prime } - (\varPsi _{1L}+ Y_{1L}^2)\varPhi _{1L}^{\prime } \big]\mathrm{d}\xi _1+\int _{-L/2}^{L/2} \varPsi _{2U} \varPhi _{2U}^{\prime } \, \mathrm{d}\xi _2\\ &\hphantom {=}+\int _{-L/2}^{L/2}\bigg [\frac {1}{3}(Y_{1U}^3 X_{1U}^{\prime } - Y_{1L}^3 X_{1L}^{\prime }) +\frac {1}{2}Y_{1U}^2 X_{1U}^{\prime }\bigg ] \, \mathrm{d}\xi _1, \end{aligned} \end{equation}

are obtained.

3.4. The steady conformal mapping

The remainder of this paper considers travelling-wave solutions of the time-independent formulation. We now study solutions that are steady in a frame of reference moving to the right with non-dimensional speed $c$ . The Galilean transformation $x \mapsto x+ct$ results in $\phi _{1t} \mapsto \phi _{1t} - c\phi _{1x}$ and $\zeta _{1t}\mapsto \zeta _{1t} -c\zeta _{1x}$ in boundary conditions (2.1c ) and (2.1d ). Then, by setting $\phi _{1t}$ and $\zeta _{1t}$ to zero, we obtain solutions that are steady in the co-moving frame.

First, we convert the steady kinematic boundary conditions (2.1c ) and (2.1e ) at the free surface and vorticity interface into conformal coordinates, which yields

(3.12) \begin{equation} \varPsi _{1U}^{\prime }=(c- Y_{1U})Y_{1U}^{\prime }, \quad \varPsi _{1L}^{\prime }=(c- Y_{1L})Y_{1L}^{\prime },\quad \varPsi _{2U}^{\prime }=cY_{2U}^{\prime }. \end{equation}

These exact solutions may be integrated to obtain

(3.13) \begin{equation} \varPsi _{1U}=C_1+ cY_{1U} -\frac {1}{2}Y_{1U}^{2}, \quad \varPsi _{1L}=C_2+cY_{1L}-\frac {1}{2}Y_{1L}^{2},\quad \varPsi _{2U}=C_3+cY_{2U}, \end{equation}

where $C_1$ , $C_2$ and $C_3$ are constants of integration. We have $C_2=C_3$ from enforcing $\psi _2=\psi _1+{y^2}/{2}$ at the vorticity jump, and the solution $C_2=C_3=-c L^{-1}\int _{-L/2}^{L/2}Y_{2U}\mathrm{d}\xi _2$ is obtained from (3.1b ). The constant $C_1$ is an unknown to be determined. Next, we consider the dynamic boundary condition at the free surface in (2.1d ), which yields

(3.14) \begin{equation} \frac {(\varPsi _{1U}^{\prime })^2 +(\varPhi _{1U}^{\prime })^2}{2J_{1U}}+( Y_{1U}-c)\frac {Y_{1U}^{\prime } \varPsi _{1U}^{\prime } + X_{1U}^{\prime } \varPhi _{1U}^{\prime }}{J_{1U}}- \varPsi _{1U}+\frac {Y_{1U}}{F^2}=B_1. \end{equation}

To summarise, the steady boundary-value problem involves the four unknown functions $Y_{1U}(\xi _1)$ , $Y_{1L}(\xi _1)$ , $Y_{2U}(\xi _2)$ and $\gamma (\hat {\xi })$ , which are closed by the four equations (3.14), (3.9a ), (3.9b ) and (3.9c ). There are also six unknown constants: the non-dimensional wavespeed $c$ , non-dimensional wavelength $L$ , non-dimensional fluid mass $H$ , amplitude $A$ from (2.10), Bernoulli constant $B_1$ and constant of integration $C_1$ from (3.13). We fix values for $L$ , $H$ and $A$ , for which imposing the two mass conditions (3.10a ) and (3.10b ) and the amplitude constraint (2.10) yields $c$ , $B_1$ and $C_1$ as unknowns of the problem.

3.5. Numerical implementation of the steady model formulation

We use a pseudospectral method to obtain numerical solutions to the steady system of equations described in § 3.4, in which the solutions are found by collocation. For fixed values of $L$ , $H$ and $A$ , we express each solution as a Fourier series of the form

(3.15) \begin{equation} \{Y_{1U} , Y_{1L}, Y_{2U}\}=\sum _{m=0}^{N} \{a_m,b_m,c_m \} \cos {\big ( m k \hat {\xi }\big )}, \quad \gamma (\hat {\xi })=\sum _{m=1}^{N} d_m\sin {\big ( m k \hat {\xi }\big )}, \end{equation}

where $k=2 \pi n/L$ as above. We use $N=80$ and $n=1$ throughout this work unless otherwise specified. Derivatives and integral transforms are computed spectrally by using the Fourier multipliers ${d}/\mathrm{d}\xi [\mathrm{e}^{{i}m k \xi }]=\mathrm{i}m k \mathrm{e}^{{i}m k \xi }$ , $\mathscr{H}\,\,[\mathrm{e}^{{i}mk \xi _2}]=\mathrm{i} \boldsymbol{\cdot }\text{sgn}(m) \mathrm{e}^{{i}mk\xi _2}$ , $\mathscr{T}\,\,[\mathrm{e}^{{i}mk \xi _1}]=\mathrm{i}\boldsymbol{\cdot }\coth {(m k h)} \mathrm{e}^{{i}mk\xi _1}$ and $\mathscr{S}\,\,[\mathrm{e}^{{i}mk \xi _1}]=\mathrm{i}\boldsymbol{\cdot }\text{csch}{(m k h)}\mathrm{e}^{{i}m k\xi _1}$ . To evaluate the solutions in physical space, we use the inverse Fourier transform with the fast-Fourier-transform algorithm.

To evaluate each matching condition (3.9a ), (3.9b ) and (3.9c ), which involve functions represented by Fourier series evaluated at $\gamma (\hat {\xi })$ having the Fourier expansion from (3.15), we sum each component of the series manually. That is, when evaluating $Y_{1L}(\hat {\xi }-\gamma (\hat {\xi }))$ for instance, we calculate

(3.16) \begin{equation} Y_{1L}(\hat {\xi }-\gamma (\hat {\xi }))=\sum _{m=0}^{N} b_m \cos {\big (mk(\hat {\xi }-\gamma (\hat {\xi }))\big )} \end{equation}

term by term by resumming the series expansion. Manually resumming each term involved in the three matching conditions is the slowest part of the overall numerical scheme.

Due to the assumed symmetry of $Y_{1U}$ , $Y_{1L}$ and $Y_{2U}$ , and the asymmetry of $\gamma (\hat {\xi })$ , only half of the periodic domain needs to be considered. We therefore introduce the collocation points $\xi _i=-L/2 + L(i-1)/(2N+1)$ for $1 \leqslant i \leqslant N+1$ , which results in $4N+6$ unknown constants obtained in total that we need to solve for. These are $\{ Y_{1U}(\xi _{i}) \}_{i=1}^{N+1}$ , $\{ Y_{1L}(\xi _i) \}_{i=1}^{N+1}$ , $\{ Y_{2U}(\xi _i) \}_{i=1}^{N+1}$ , $\{ \gamma (\xi _i) \}_{i=2}^{N+1}$ , $c$ , $B_1$ and $C_1$ . Next, we obtain $4N+6$ algebraic equations to close the system. Three of these constraints are given by the mass and amplitude conditions (3.10a ), (3.10b ) and (2.10). The remaining $4N+3$ constraints come from evaluating the remaining four equations at the collocation points: (3.14), (3.9b ) and (3.9c ) each provide $N+1$ constraints due to being even and (3.9a ) provides $N$ constraints due to being odd.

Once each of the $4N+3$ algebraic constraints have been obtained, the solution is found using Newton iteration to minimise these constraints. We continue iterating until the $L^2$ -norm of the constraints falls below a residual of $10^{-12}$ , which takes approximately 3 s when $N=100$ Fourier modes are used in (3.15). The initial guess used is either taken to be a linear solution from § 2.2, or a previously computed numerical solution obtained with specified values for $L$ , $H$ and $A$ . We typically begin with the linear solution as an initial guess for a small value of $A$ , such as $A=10^{-5}$ , and then sequentially use the computed solution as the initial guess for a larger value of $A$ . Once the desired value has been reached, such as $A=0.03$ , we then employ arclength continuation in the $(H,c)$ plane to explore the bifurcation curves that pass through this computed solution.

We find that significant aliasing errors emerge when computing solutions that are dominated by the vorticity interface. These are resolved by using the 36th-order Fourier filter $\exp {(-36 (1.2 m/N)^{36})}$ during the Newton-iteration procedure, and we found that best performance is obtained when only $Y_{1L}$ , $Y_{2U}$ and $\gamma$ are de-aliased once per iteration, while no dealiasing is applied to $Y_{1U}$ . The 36th-order filter was analysed by Hou & Li (Reference Hou and Li2007) for time-dependent problems, and forms a balance between the suppression of aliasing errors and resolving the remaining Fourier coefficients in the spectrum. The filter used is shown below as a dashed line alongside the Fourier spectrum of each solution in figure 14 and, for the solutions in this figure that required dealiasing, we were able to compute solutions of approximately three times larger amplitude than those computed without dealiasing.

4.1. Solutions of limiting amplitude

We next investigate the behaviour of solution profiles as they transition to waves of limiting form. By taking the six solutions with $A=0.03$ and $n=1$ from figure 7 that have either $H=0.5$ or $H=2.5$ , we increase the amplitude $A$ until the decay of the Fourier coefficients plateaus. These six solutions correspond to those for $A=10^{-5}$ , denoted by ( $a,d,f)$ and $(c,g,i)$ in figure 4, where $(d,i)$ and $(a,c,f,g)$ represented the shear and gravity-wave modes, respectively. The plateau of the Fourier coefficients is indicative of the clustering of a countably infinite number of singularities in the analytic continuation of the solution profiles, which results in the formation of a corner (which is also a stagnation point) at the wave crest. The specific termination criterion used here is that the largest of the four Fourier coefficients $a_N$ , $b_N$ , $c_N$ and $d_N$ from (3.15) exceeds $10^{-10}$ . In figure 12, we show six near-limiting solutions calculated with $N=100$ Fourier modes. The labelling of these matches that used for the linear solutions in figure 4. These solutions display different transitions to solutions of maximal form. For the solutions shown in $(d)$ and $(i)$ , the vorticity interface $\zeta _2(x)$ first becomes limiting, whereas for the solutions shown in $(a)$ , $(c)$ , $(g)$ and $(f)$ , the free surface $\zeta _1(x)$ first becomes limiting. We observe that the amplitudes of limiting free surface and internal waves differ greatly depending on the parameter values, or wave modes. For instance, the amplitudes of $\zeta _1(x)$ in isolation for each of the six profiles is $A_{1}=(0.0525,0.0414,0.0134,0.257,0.181,0.00432)$ in alphabetical order, whereas those for $\zeta _2(x)$ are $A_{2}=(0.0180,0.000318,0.147,0.0381,0.110,0.114)$ . Note that the amplitude of the limiting gravity (Stokes) wave on an irrotational fluid of infinite depth is $0.1412$ (Cokelet Reference Cokelet1977). This limiting solution forms a corner at the wave crest with an internal angle of $120$ degrees. This also occurs for limiting waves with vorticity (Kozlov & Lokharu Reference Kozlov and Lokharu2021).

Table 1.

Wavespeed, momentum and energy values are recorded, to four decimal places, for the six solutions shown in both figures 12 and 14. The wavespeed $c$ forms an unknown of the problem, while the momentum $P$ and energy $E$ are calculated from (3.11b ) and (3.11a ).

Table 1. Long description

A table with two rows of headers and twelve columns. The first row of headers includes Solution, c, P, E, Solution, c, P, and E. The second row of headers repeats Solution, c, P, and E. The table presents values for wavespeed (c), momentum (P), and energy (E) for different solutions labeled 11(a) to 11(i). Row 1: Solution 11(a), c 1.1623, P 0.6685, E 0.9147, Solution 11(c), c 3.1112, P 15.6526, E 41.8264. Row 2: Solution 11(d), c 0.2838, P 0.7757, E 0.9158, Solution 11(g), c 0.8764, P 14.3236, E 38.5486. Row 3: Solution 11(f), c -0.8090, P 0.0600, E 1.4275, Solution 11(i), c 0.4261, P 15.7176, E 41.7431.

Figure 13.

Streamlines and stagnation points are shown for the six solution profiles from figure 12. Streamlines (shown in grey) are calculated as contours of the streamfunctions $\hat {\psi }_1$ and $\hat {\psi }_2$ defined in (B1). Stagnation points (shown with circles) are calculated from (B2) and (B3). Ten streamlines (solid grey) are shown in the upper fluid layer for each profile, with the value on each contour uniformly dividing the range min $(\hat {\psi }_1)$ to max $(\hat {\psi }_1)$ . The same spacing, $[\text{max}(\hat {\psi }_1)-\text{min}(\hat {\psi }_1)]/11$ , between contour values has been used to plot streamlines in the lower fluid layer. Additional streamlines (dashed grey) show those emerging from stagnation points.

Streamlines and stagnation points for each of the six solutions displayed in figure 12 are shown in figure 13. These were calculated using the method described in Appendix B, which finds values for the internal velocity field within each fluid layer from solutions of the boundary-integral scheme. Interestingly, closed streamlines – and therefore stagnation points – are observed within the shear layer for solutions $(d,g,i)$ . For fluids with constant vorticity, stagnation points typically emerge either within the fluid or on the free surface at the lines of symmetry $x=0$ , $x=\pm L/2$ (Ribeiro et al. Reference Ribeiro, Milewski and Nachbin2017). Solution $(d)$ , which has two wave crests near $x=0$ on the free surface $\zeta _1(x)$ , has three stagnation points: one inside the fluid and two on the free surface. Solutions $(g)$ and $(i)$ each have two internal stagnation points, one at $x=0$ and another at $x=\pm L/2$ .

Figure 14.

The Fourier spectrum for each of the six solutions from figure 12 is shown. The absolute value of each Fourier coefficient for the free surface $\zeta _1(x)$ , corresponding to $a_m$ from (3.15), is shown with filled circles and that for the vorticity interface $\zeta _2(x)$ , corresponding to $c_m$ from (3.15), is shown with hollow pentagons. The Fourier filter $\exp {(-36 (1.2 m/N)^{36})}$ , used for solutions with a slow decay of $b_m$ , is shown dashed in $(d)$ and $(i)$ .

Figure 14. Long description

The image contains nine line graphs labeled (a), (c), (d), (f), (g), and (i). Each graph shows the absolute value of Fourier coefficients on the vertical axis and the index m on the horizontal axis. The graphs compare the Fourier spectrum of the free surface and the vorticity interface. Filled circles represent the Fourier coefficients for the free surface, while hollow pentagons represent those for the vorticity interface. A dashed line indicates the Fourier filter used for solutions with a slow decay of Fourier coefficients. In graphs (a), (c), (f), and (i), the Fourier coefficients for the free surface and vorticity interface show a general decay as m increases. Graphs (d), (g), and (i) include a dashed line representing the Fourier filter, which is applied to resolve aliasing errors. The filter helps balance the suppression of aliasing errors and the resolution of remaining Fourier coefficients in the spectrum.

The Fourier spectrum for each of the six solutions plotted in figure 12 is displayed in figure 14. The decay of these Fourier coefficients indicates which of the two interfaces first develops a corner at the wavecrest as the amplitude continues to increase. For solutions $(d)$ and $(i)$ , it is the vorticity interface $\zeta _2(x)$ that has the slowest decay of Fourier coefficients, which indicates that this interface will become limiting first.