1 Introduction
The huge literature on free surface water waves is mainly based on the assumption of irrotational flow. Nonetheless, there are many circumstances under which one cannot ignore the rotational character of the flow, namely, in the presence of an underlying vertically sheared current. For example, on deep water when a thin shear layer is generated by wind stress, or in coastal zones where the vertical profiles of mean velocity are typically determined by bottom friction and also surface wind stress. For a review on the interaction of waves with vertically sheared currents occurring in nature, one can cite Peregrine (Reference Peregrine1976), Jonsson (Reference Jonsson1990) and Thomas & Klopman (Reference Thomas, Klopman and Hunt1997).
Periodic gravity waves propagating steadily on a rotational current were studied by many authors, either mathematically or analytically and numerically. Since the pioneering work of Constantin & Strauss (Reference Constantin and Strauss2004), there is a huge literature that concerns the formulation and the mathematical properties of steady periodic surface water waves with vorticity (existence, unicity, bifurcation). For a review on recent rigorous results the reader can refer to Constantin & Varvaruca (Reference Constantin and Varvaruca2011) and Kozlov & Kuznetsov (Reference Kozlov and Kuznetsov2014). Among the authors using asymptotic methods or purely numerical methods, on can cite Tsao (Reference Tsao1959), Dalrymple (Reference Dalrymple1974), Brevik (Reference Brevik1979), Simmen & Saffman (Reference Simmen and Saffman1985), Kishida & Sobey (Reference Kishida and Sobey1988), Teles da Silva & Peregrine (Reference Teles da Silva and Peregrine1988), Vanden-Broeck (Reference Vanden-Broeck1996), Swan & James (Reference Swan and James2001), Ko & Strauss (Reference Ko and Strauss2008), Cheng, Cang & Liao (Reference Cheng, Cang and Liao2009), Pak & Chow (Reference Pak and Chow2009), Moreira & Chacaltana (Reference Moreira and Chacaltana2015), Hsu et al. (Reference Hsu, Francius, Montalvo and Kharif2016) and Ribeiro Jr, Milewski & Nachbin (Reference Ribeiro, Milewski and Nachbin2017).
Although the recent important theoretical developments have confirmed that periodic waves can exist over flows with arbitrary vorticity, it appears that their stability to infinitesimal disturbances and their subsequent nonlinear evolution have not been studied extensively so far. In fact, even in the rather simple case of uniform vorticity (linear shear), few papers have been published on the effect of a vertical shear current on the sideband instability of a uniform wave train over finite depth. It is noted here that this instability, which is related to a four-wave resonance (between two small sidebands and two quanta of the strong carrier wave), is often referred in the literature as the modulational instability, or Benjamin–Feir instability (provided the modulation frequency is small compared to the carrier frequency).
Under the weakly nonlinear assumption, Johnson (Reference Johnson1976) studied the slow modulation of a harmonic wave moving on an arbitrary shear flow, in two dimensions only. Using the method of multiple scales he obtained a two-dimensional nonlinear Schrödinger equation (NLS equation), from which a condition of linear stability of its plane wave solution (which corresponds to the Stokes wave) was proposed. Using a similar approach and extending the analysis to three-dimensional disturbances, Oikawa, Chow & Benney (Reference Oikawa, Chow and Benney1987) obtained a set of envelope evolution equations to study the effect of shear on the directional dependence of the modulational instabilities of weakly nonlinear wave packets. They also derived an analytical expression of the linear growth rate of three-dimensional modulational instabilities, and quantitative results were reported only in the simple case of linear shear. Although the numerical results of Oikawa et al. (Reference Oikawa, Chow and Benney1987) suggest that the shear can modify significantly the three-dimensional characteristics of the modulational instabilities of the plane wave solution, it is emphasized here that their parameterization of the undisturbed linear shear flow in terms of one single dimensionless parameter does not enable us to discriminate between the effects of a uniform current (equal to the value of the current at the surface) and those of the vorticity. In both papers the resulting third-order envelope equations possess coefficients that depend in a complicated way on the shear, and are thus not practical to use. Even in the simple case of a linear shear current, some integrals require numerical evaluation. More significantly, these analytical results are obtained with the assumption that no critical layer can occur, a crucial feature which is known to alter the kinematic nature of the wave motion, either it is a periodic wave or a solitary wave as shown by Miroshnikov (Reference Miroshnikov2002). These approaches are therefore limited to certain regions of the parameter space, defined in the present work by the wave steepness
$\unicode[STIX]{x1D700}$
, the dimensionless depth
$kh$
and a dimensionless parameter (to be defined later) characterizing the (constant) vorticity of the background shear flow. Here
$k$
and
$h$
are the carrier wavenumber and depth, respectively.
For flows with constant vorticity, Li, Hui & Donelan (Reference Li, Hui, Donelan, Horikawa and Maruo1987) studied the two-dimensional sideband instability of a Stokes wave train in deep water. They also obtained an NLS equation for the amplitude modulations, but the coefficient of the nonlinear term was erroneous, as noted by Baumstein (Reference Baumstein1998). In fact, the latter author investigated the effect of piecewise-linear velocity profiles on the sideband instability of a finite-amplitude gravity wave in deep water. Extension of these studies to the case of finite depth was carried out by Thomas, Kharif & Manna (Reference Thomas, Kharif and Manna2012), for two dimensions only, by using the same method of multiple scales but stemming from a different formulation of the governing equations. The resulting coefficients of their NLS equation are given explicitly as a function of the vorticity and depth of the shear layer. Explicit analytical expressions are also given for several characteristics of the two-dimensional modulational instabilities of the Stokes wave. In this paper, it is found that vorticity modifies significantly the modulational instability properties of weakly nonlinear plane waves, namely the growth rate and bandwidth. At third order, the paper also shows the importance of the nonlinear coupling between the mean flow induced by the modulation and the vorticity. Furthermore, it is shown that the plane wave solution can be linearly stable to modulational instability for an opposite shear current (positive vorticity in the present work) independently of the dimensionless depth parameter
$kh$
. For
$kh>1.363$
, in practice, this important property means that large vorticities can suppress the modulational instability in both following and opposite shear currents with constant vorticity.
As pointed out by one referee, Hur & Johnson (Reference Hur and Johnson2015) have recently studied the modulational stability and instability of periodic travelling waves on a linear shear current in finite depth, using the so-called vorticity-modified Whitham equation that combines the full linear dispersion relation of water waves and weak nonlinearity stemming from an approximation of the fully nonlinear ‘breaking operator’ of shallow water theory. For any constant vorticity they have shown that periodic travelling gravity waves (with sufficiently small amplitudes) are spectrally unstable to long-wavelength perturbations if
$kh$
is greater than a critical value that depends on the vorticity, and stable otherwise, similar to the zero vorticity setting. For irrotational gravity waves, however, the complete restabilization to the modulational instability is found to occur when
$kh$
becomes less than 1.146, which differs from the known threshold value 1.363. In this respect, it was shown in Thomas et al. (Reference Thomas, Kharif and Manna2012) that the introduction of constant vorticity does not modify this critical value, which is obtained in the limit of zero vorticity. Despite this failure, the work of Hur & Johnson (Reference Hur and Johnson2015) confirms that sufficiently large vorticities can suppress the modulational instability in the case of following shear currents.
To the best of our knowledge, there do not appear to have been many attempts, with extensive use of numerical techniques, on the study of the stability of finite-amplitude gravity waves on currents with vertical shear. Following the general scheme, pioneered by Longuet-Higgins (Reference Longuet-Higgins1978a ,Reference Longuet-Higgins b ), Okamura & Oikawa (Reference Okamura and Oikawa1989) found three kinds of instabilities: two of them related to four- and five-wave resonant interactions and a third type, called rotational instability which is essentially three-dimensional and due to the effects of a strong nonlinear coupling between the mean-flow response and the fundamental wave, as discovered first by Chow & Benney (Reference Chow and Benney1986). These results have been compared with those for the corresponding irrotational wave, to show that increasing the shear yields larger growth rates for each most unstable mode associated with the four-wave and five-wave resonances. Although, for the first time, the rotational instabilities have been captured numerically from the Euler equations for finite-amplitude waves and constant vorticity, the whole stability domain could not be examined owing to the poor convergence property of their numerical method. As in many published works, it is noted here that the same dimensionless parameter is used for both the surface current and the constant shear, which gives rise to some ambiguity, as already mentioned.
Recently, in the same vein, numerical simulations of the Euler equations and of some high-order approximation, for two-dimensional rotational flows, have been used to analyse the nonlinear stability of finite-amplitude periodic waves, initially weakly modulated due to the presence of sideband disturbances. Nwogu (Reference Nwogu2009) considered an exponentially sheared current and reported few results on the long-time evolution of modulational instabilities of periodic gravity waves in deep water. His numerical results demonstrated that the mean-flow vorticity can significantly affect the growth rate of extreme waves in narrowband sea states. For flows with constant vorticity, Choi (Reference Choi2009) used a conformal mapping method to obtain new evolution equations, which are solved numerically using a pseudo-spectral method to study the Benjamin–Feir instability of a modulated wave train in both positive and negative shear currents. For fixed wave steepness and (unstable) disturbances, Choi (Reference Choi2009) found that the envelope of the modulated wave train grows faster in a positive shear current and slower in a negative one, in comparison with the irrotational case. In contrast to the above mentioned studies on the stability problem, it is emphasized here that both Choi (Reference Choi2009) and Nwogu (Reference Nwogu2009) parameterized the background current with two different dimensionless parameter, one for the surface current and one for the vorticity.
In this work we consider the two-dimensional linear stability analysis of finite-amplitude gravity waves propagating steadily on the free surface of a fluid with constant vorticity and finite depth. As we shall see, the approximation of constant vorticity is twofold. Not only does it simplify considerably the mathematical analysis of two dimensional steady-state waves, which are then necessarily irrotational perturbations of the background uniform shear flow (as a consequence of Kelvin’s theorem), but it also allows a straightforward extension of some existing numerical methods, originally dedicated to the study of irrotational finite-amplitude waves on finite depth.
The paper is structured as follows. In § 2 we formulate the systems of equations both for the determination of steady periodic waves on a linear shear current, and for the evolution of superposed small disturbances, considered in a frame moving with the basic state. By introducing a velocity potential function for the (necessary) disturbances and using a normal-mode analysis, we formulate the eigenvalue problem to be solved. In § 3, we present the numerical methods used in this study and the results of their validation against available numerical and analytical results. In § 4, the results for subharmonic perturbations are presented, first, for the sideband instability (long-wavelength disturbances) and then for the shorter disturbances involved in quartet and quintet resonances with the basic wave. We summarize our conclusions in § 5.
2 Mathematical formulation
2.1 Euler equations and periodic wave solutions
The fluid is assumed to be inviscid, the flow is incompressible, the water depth is taken to be finite and the effect of surface tension is neglected. In the absence of waves, we shall assume a basic parallel shear current
$U_{B}(z)$
of constant vorticity (
$-\unicode[STIX]{x1D6FA}_{0}$
), which varies linearly in the vertical direction and vanishes at the undisturbed free surface or mean water level (
$z=0$
).
For periodic two-dimensional (2-D) waves of permanent form on a linear shear current it is convenient to consider a coordinate system
$(x,z)$
moving with the wave speed
$c>0$
, in which the combined wave–current flow is steady. From the condition of incompressibility we can then introduce a streamfunction
$\unicode[STIX]{x1D713}(x,z)$
, unique up to an additive constant, such that the total relative fluid velocities are given by
$$\begin{eqnarray}u=U_{B}(z)-c+\tilde{u} =\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z},\quad w=\tilde{w}=-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}x},\end{eqnarray}$$
where the tilde quantities are perturbations due to the waves. We note that by defining
$U_{B}(z)=\unicode[STIX]{x1D6FA}_{0}z$
, the focus is placed on the influence of vorticity on the combined wave–current system rather than on the uniform current (Doppler effect).
As is well known for two-dimensional incompressible flows, the vorticity vector is perpendicular to the plane of the motion and, in the absence of vortex stretching, the vorticity of a fluid particle,
$\unicode[STIX]{x1D714}=(\unicode[STIX]{x2202}_{x}w-\unicode[STIX]{x2202}_{z}u)$
, is constant along a streamline. Since in the moving frame of reference the combined wave–current flow of constant vorticity is steady, we have
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z^{2}}=\unicode[STIX]{x1D6FA}_{0}\end{eqnarray}$$
as the governing equation to be satisfied throughout the fluid, namely
$-h<z<\unicode[STIX]{x1D702}(x)$
where
$h$
is the mean depth and
$z=\unicode[STIX]{x1D702}(x)$
represents the unknown position of the free surface relative to the mean surface level
$z=0$
.
The kinematic boundary conditions, which express that both the free surface and the rigid plane bed are streamlines, can be written as
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}=0\quad \text{on }z=-h & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}=-Q\quad \text{on }z=\unicode[STIX]{x1D702}, & \displaystyle\end{eqnarray}$$
where
$Q$
represents the total volume flow rate underneath the wave (per unit span). Using the Bernoulli equation at the free surface and setting the surface pressure to zero without loss of generality gives the dynamic free surface condition,
$$\begin{eqnarray}g\unicode[STIX]{x1D702}+{\textstyle \frac{1}{2}}(u^{2}+w^{2})=R\quad \text{on }z=\unicode[STIX]{x1D702},\end{eqnarray}$$
where
$g$
is the gravitational acceleration and
$R$
is a Bernoulli constant of the flow, which corresponds physically to the specific energy at the free surface (in the moving frame).
The values of the constants
$Q$
and
$R$
are unknowns and need to be determined for solutions
$\unicode[STIX]{x1D713}(x,z)$
and
$\unicode[STIX]{x1D702}(x)$
that are periodic in
$x$
with finite wavelength
$\unicode[STIX]{x1D706}=2\unicode[STIX]{x03C0}/k$
, where
$k$
is the wavenumber. To guarantee a unique solution, when it exists, the wave height
$H$
is specified as
$$\begin{eqnarray}\unicode[STIX]{x1D702}(0)-\unicode[STIX]{x1D702}(\unicode[STIX]{x1D706}/2)=H,\end{eqnarray}$$
where
$\unicode[STIX]{x1D702}(0)$
and
$\unicode[STIX]{x1D702}(\unicode[STIX]{x1D706}/2)$
are the crest and trough elevation respectively, and we impose the further condition of zero mean level condition,
$$\begin{eqnarray}\int _{0}^{\unicode[STIX]{x1D706}}\unicode[STIX]{x1D702}\,\text{d}x=0.\end{eqnarray}$$
Solutions for steep waves on a linear shear current are not currently obtainable analytically, and numerical approaches are necessary to solve the above nonlinear boundary-value problem. In this work we use a simple extension of the Fourier approximation method, originated with Rienecker & Fenton (Reference Rienecker and Fenton1981) for irrotational waves, where a Fourier series is assumed only for the streamfunction as
$$\begin{eqnarray}\unicode[STIX]{x1D713}(x,z)=\frac{\unicode[STIX]{x1D6FA}_{0}}{2}(z^{2}-h^{2})-c(z+h)+\unicode[STIX]{x1D6FC}\mathop{\sum }_{j=1}^{\infty }B_{j}\frac{\sinh jk(h+z)}{\cosh jkh}\cos jkx.\end{eqnarray}$$
Solutions with this form satisfy automatically the Poisson equation and the kinematic boundary condition at the bottom, as well as the periodicity condition. In the above expression
$\unicode[STIX]{x1D6FC}=\sqrt{g/k^{3}}$
has been introduced for convenience, as we shall work with dimensionless values, so that the coefficients
$B_{j}$
are dimensionless. When a solution exists for given wave height
$H$
, mean depth
$h$
and constant vorticity
$-\unicode[STIX]{x1D6FA}_{0}$
, the coefficients
$B_{j}$
, the wave speed
$c$
and the unknown position of the free surface
$\unicode[STIX]{x1D702}(x)$
can be determined numerically, after truncation of the Fourier series and a collocation method as described in § 3.
2.2 Linear stability analysis
By confining attention to two-dimensional flows, it should be realized that any perturbations of a flow with constant vorticity (background linear shear current), either small or of finite amplitude, are necessarily irrotational motions as a consequence of Kelvin’s circulation theorem. Hence, in the frame moving with wave speed
$c$
, it follows that we can introduce a generalized velocity potential
$\unicode[STIX]{x1D719}(x,z,t)$
, which satisfies the relations
$$\begin{eqnarray}\unicode[STIX]{x1D719}_{x}=\unicode[STIX]{x1D713}_{z}-\unicode[STIX]{x1D6FA}_{0}(z+h),\quad \unicode[STIX]{x1D719}_{z}=-\unicode[STIX]{x1D713}_{x}\end{eqnarray}$$
to guarantees that
$\unicode[STIX]{x1D719}$
is the harmonic conjugate of the function
$\unicode[STIX]{x1D713}^{h}=\unicode[STIX]{x1D713}(x,z,t)-(\unicode[STIX]{x1D6FA}_{0}(z+h)^{2})/2$
.
The unsteady water wave problem can then be reformulated in the moving frame and, in terms of
$\unicode[STIX]{x1D713}$
and
$\unicode[STIX]{x1D719}$
, the governing equations are given by
$$\begin{eqnarray}\unicode[STIX]{x1D713}_{xx}+\unicode[STIX]{x1D713}_{zz}=\unicode[STIX]{x1D6FA}_{0},\quad -h<z<\unicode[STIX]{x1D702}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}_{x}=0\quad \text{on }z=-h & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}_{t}+\unicode[STIX]{x1D702}_{x}\unicode[STIX]{x1D713}_{z}+\unicode[STIX]{x1D713}_{x}=0\quad \text{on }z=\unicode[STIX]{x1D702}(x,t) & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{t}+g\unicode[STIX]{x1D702}-\unicode[STIX]{x1D6FA}_{0}\unicode[STIX]{x1D713}+{\textstyle \frac{1}{2}}(\unicode[STIX]{x1D713}_{z}^{2}+\unicode[STIX]{x1D713}_{x}^{2})=f(t)+{\textstyle \frac{1}{2}}c^{2}-\unicode[STIX]{x1D6FA}_{0}ch\quad \text{on }z=\unicode[STIX]{x1D702}(x,t), & \displaystyle\end{eqnarray}$$
where
$f(t)$
is an arbitrary function of time, which can be set to zero without loss of generality or absorbed into the definition of either
$\unicode[STIX]{x1D719}$
or
$\unicode[STIX]{x1D713}^{h}$
. If the motion is steady, as in the case of waves of permanent form, the terms with a time derivative disappear and we recover the governing equations of the steady water wave problem, as it has been presented in the previous section. Note that in this case the kinematic conditions (2.11) and (2.12) both express that
$\unicode[STIX]{x1D713}$
is constant at the bottom and at the free surface, respectively.
To consider the stability of these steady waves to infinitesimal two-dimensional disturbances, we let
$$\begin{eqnarray}\unicode[STIX]{x1D702}=\bar{\unicode[STIX]{x1D702}}(x)+\tilde{\unicode[STIX]{x1D702}}(x,t),\quad \unicode[STIX]{x1D713}=\bar{\unicode[STIX]{x1D713}}(x)+\tilde{\unicode[STIX]{x1D713}}(x,t),\end{eqnarray}$$
where the quantities with an overbar correspond to the periodic wave solution and the quantities with the tilde denote small perturbations, i.e.
$|\tilde{\unicode[STIX]{x1D702}}|\ll |\bar{\unicode[STIX]{x1D702}}|$
and
$|\tilde{\unicode[STIX]{x1D713}}|\ll |\bar{\unicode[STIX]{x1D713}}|$
. After linearization of the governing equations (2.10)–(2.13) we obtain
$$\begin{eqnarray}\tilde{\unicode[STIX]{x1D713}}_{xx}+\tilde{\unicode[STIX]{x1D713}}_{zz}=0,\quad -h<z<\bar{\unicode[STIX]{x1D702}}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D713}}_{x}=0\quad \text{on }z=-h & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D702}}_{t}+\bar{\unicode[STIX]{x1D713}}_{z}\tilde{\unicode[STIX]{x1D702}}_{x}+\bar{\unicode[STIX]{x1D702}}_{x}\tilde{\unicode[STIX]{x1D713}}_{z}+(\bar{\unicode[STIX]{x1D713}}_{zz}\bar{\unicode[STIX]{x1D702}}_{x}+\bar{\unicode[STIX]{x1D713}}_{xz})\tilde{\unicode[STIX]{x1D702}}+\tilde{\unicode[STIX]{x1D713}}_{x}=0\quad \text{on }z=\bar{\unicode[STIX]{x1D702}}(x) & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D719}}_{t}+(g+\bar{\unicode[STIX]{x1D713}}_{x}\bar{\unicode[STIX]{x1D713}}_{xz}+\bar{\unicode[STIX]{x1D713}}_{z}\bar{\unicode[STIX]{x1D713}}_{zz}-\unicode[STIX]{x1D6FA}_{0}\bar{\unicode[STIX]{x1D713}}_{z})\tilde{\unicode[STIX]{x1D702}}-\unicode[STIX]{x1D6FA}_{0}\tilde{\unicode[STIX]{x1D713}}+\bar{\unicode[STIX]{x1D713}}_{z}\tilde{\unicode[STIX]{x1D713}}_{z}+\bar{\unicode[STIX]{x1D713}}_{x}\tilde{\unicode[STIX]{x1D713}}_{x}=0\quad \text{on }z=\bar{\unicode[STIX]{x1D702}}(x).\qquad & \displaystyle\end{eqnarray}$$
By using the Cauchy–Riemann relations (2.9) it is easy to check that without vorticity
$(\unicode[STIX]{x1D6FA}_{0}=0)$
the above equations reduce to those derived by McLean (Reference McLean1982).
Since the linearized equations (2.15)–(2.18) have coefficients that are periodic in
$x$
, it follows from Floquet or Bloch wave theory that its general solution can be represented as a linear combination of normal modes in the form
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D702}}(x,t)=\text{e}^{-\text{i}\unicode[STIX]{x1D6FE}t}\text{e}^{\text{i}px}\mathop{\sum }_{j=-\infty }^{\infty }a_{j}\text{e}^{\text{i}jx} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D713}}(x,z,t)=\text{e}^{-\text{i}\unicode[STIX]{x1D6FE}t}\text{e}^{\text{i}px}\mathop{\sum }_{j=-\infty }^{\infty }b_{j}\text{e}^{\text{i}jx}\frac{\text{sinh}[k_{j}(z+h)]}{\text{cosh}(k_{j}h)} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D719}}(x,z,t)=\text{e}^{-\text{i}\unicode[STIX]{x1D6FE}t}\text{e}^{\text{i}px}\mathop{\sum }_{j=-\infty }^{\infty }c_{j}\text{e}^{\text{i}jx}\frac{\text{cosh}[k_{j}(z+h)]}{\text{cosh}(k_{j}h)}, & \displaystyle\end{eqnarray}$$
where
$k_{j}=|p+j|$
,
$p$
an arbitrary real number and
$\unicode[STIX]{x1D6FE}$
is an unknown, possibly complex, eigenvalue to be determined by the requirement that the linearized surface equations (2.17) and (2.18) have a non-trivial solution. Here we have taken the gravitational acceleration
$g=1$
and the unperturbed wavelength
$\unicode[STIX]{x1D706}=2\unicode[STIX]{x03C0}$
, which is equivalent to non-dimensionalizing the problem with the spatial and temporal scales
$l_{ref}=\unicode[STIX]{x1D706}/(2\unicode[STIX]{x03C0})$
and
$t_{ref}=\sqrt{gk}$
(see § 3.1 for more details).
The real physical disturbance, which corresponds to the real part of the above expressions, are not strictly periodic in
$x$
unless
$p$
is a rational number. When
$p$
is an integer (or
$p=0$
without loss of generality), the wavelength of the disturbance is the same as that of the undisturbed wave (
$\unicode[STIX]{x1D706}=2\unicode[STIX]{x03C0}$
) and the disturbance is called superharmonic. When
$p$
is not an integer, the disturbance contains components with wavelength greater than
$2\unicode[STIX]{x03C0}$
and it is called a subharmonic perturbation. If
$\unicode[STIX]{x1D6FE}$
is real, then the disturbance is said to be stable, whereas if
$\unicode[STIX]{x1D6FE}$
is complex, then the disturbance or its complex conjugate is unstable.
With this normal-mode decomposition and the condition of irrotationality for the disturbances, namely
$\tilde{\unicode[STIX]{x1D719}}_{x}=\tilde{\unicode[STIX]{x1D713}}_{z}$
and
$\tilde{\unicode[STIX]{x1D719}}_{z}=-\tilde{\unicode[STIX]{x1D713}}_{x}$
, we obtain the following relations between the sets of coefficients
$b_{j}$
,
$c_{j}$
$$\begin{eqnarray}c_{j}=\text{sign}(p+j)b_{j}\quad \text{for }|j|\leqslant \infty .\end{eqnarray}$$
Using this relation and substituting the normal-mode decomposition of the general solution into (2.17) and (2.18), we obtain the following system of equations
$$\begin{eqnarray}\displaystyle & \displaystyle G(x)\mathop{\sum }_{j=-\infty }^{\infty }a_{j}\text{e}^{\text{i}jx}+\mathop{\sum }_{j=-\infty }^{\infty }H_{j}(x)b_{j}\text{e}^{\text{i}jx}=\text{i}\unicode[STIX]{x1D6FE}\mathop{\sum }_{j=-\infty }^{\infty }F_{j}(x)b_{j}\text{e}^{\text{i}jx} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathop{\sum }_{j=-\infty }^{\infty }K_{j}(x)a_{j}\text{e}^{\text{i}jx}+\mathop{\sum }_{j=-\infty }^{\infty }L_{j}(x)b_{j}\text{e}^{\text{i}jx}=\text{i}\unicode[STIX]{x1D6FE}\mathop{\sum }_{j=-\infty }^{\infty }a_{j}\text{e}^{\text{i}jx} & \displaystyle\end{eqnarray}$$
for
$0\leqslant x\leqslant 2\unicode[STIX]{x03C0}$
, where
$$\begin{eqnarray}\displaystyle & \displaystyle G(x)=g+\bar{\unicode[STIX]{x1D713}}_{x}\bar{\unicode[STIX]{x1D713}}_{xz}+\bar{\unicode[STIX]{x1D713}}_{z}\bar{\unicode[STIX]{x1D713}}_{zz}-\unicode[STIX]{x1D6FA}_{0}\bar{\unicode[STIX]{x1D713}}_{z} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle H_{j}(x)=\frac{[\text{i}(p+j)\bar{\unicode[STIX]{x1D713}}_{x}-\unicode[STIX]{x1D6FA}_{0}]\text{sinh}[k_{j}(h+\bar{\unicode[STIX]{x1D702}})]+k_{j}\bar{\unicode[STIX]{x1D713}}_{z}\text{cosh}[k_{j}(h+\bar{\unicode[STIX]{x1D702}})]}{\text{cosh}(k_{j}h)} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle F_{j}(x)=-\frac{\text{i}\text{sign}(p+j)\text{cosh}[k_{j}(h+\bar{\unicode[STIX]{x1D702}})]}{\text{cosh}(k_{j}h)} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle K_{j}(x)=\text{i}(p+j)\bar{\unicode[STIX]{x1D713}}_{z}+\bar{\unicode[STIX]{x1D713}}_{zz}\bar{\unicode[STIX]{x1D702}}_{x}+\bar{\unicode[STIX]{x1D713}}_{xz} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle L_{j}(x)=\frac{\text{i}(p+j)\text{sinh}[k_{j}(h+\bar{\unicode[STIX]{x1D702}})]+k_{j}\bar{\unicode[STIX]{x1D702}}_{x}\text{cosh}[k_{j}(h+\bar{\unicode[STIX]{x1D702}})]}{\text{cosh}(k_{j}h)}. & \displaystyle\end{eqnarray}$$
The system of (2.23)–(2.24) can be interpreted as a generalized eigenvalue problem for the complex eigenvalue
$\unicode[STIX]{x1D6FE}$
with the coefficients
$a_{j}$
and
$b_{j}$
as the eigenfunctions. It is important to note that there is a degeneracy in the dependence of the eigenfunction on
$p$
, since
$p$
can be changed to
$p+n$
, where
$n$
is an integer, without changing the eigenfunction provided the coefficients
$a_{j}$
and
$b_{j}$
are relabelled and changed to
$a_{j-n}$
and
$b_{j-n}$
, accordingly. Thus, in principle, there would be no loss of generality in confining ourselves to the range
$0\leqslant p<1$
, which does not mean that our analysis is confined to this range.
2.3 Linear resonant conditions
In the limiting case where the undisturbed state has zero wave height, namely a linear shear current with a flat free surface, the analytical solutions of the eigenvalue problem are given by
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FE}_{j}^{\pm }(p)=\text{sign}(p+j)[-k_{j}c+\unicode[STIX]{x1D714}_{j}^{\pm }] & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}_{j}^{\pm }(p)=-\frac{\unicode[STIX]{x1D6FA}_{0}\text{tanh}(k_{j}h)}{2}\pm \sqrt{gk_{j}\text{tanh}(k_{j}h)+\left(\frac{\unicode[STIX]{x1D6FA}_{0}\text{tanh}(k_{j}h)}{2}\right)^{2}}, & \displaystyle\end{eqnarray}$$
where
$j=0,\pm 1,\pm 2,\ldots$
is a mode number. Since the flat surface is regarded here as the limit of steady waves on a linear shear current with the wave height approaching zero, the eigenvalues contain the wave speed
$c=\unicode[STIX]{x1D714}/k$
. The corresponding eigenvectors are such that only those components with suffices
$p+j$
are non-zero. All the eigenvalues (2.30) lie on the real axis and so the flat surface is spectrally stable. Thus the eigensets represent infinitesimal waves on top of a flat surface and the plus and minus signs designate copropagating and contrapropagating disturbances relative to the undisturbed flow.
As in the irrotational case, it is expected that for small values of the wave steepness, instability bands may occur for modes with
$p$
values near the points
$p$
where two eigenvalues for zero amplitude are equal,
$$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{j_{1}}^{\pm }(p)=\unicode[STIX]{x1D6FE}_{j_{2}}^{\pm }(p)\end{eqnarray}$$
for some integers
$j_{1}$
and
$j_{2}$
. Since we can add any integer to
$p$
, only the difference between
$j_{1}$
and
$j_{2}$
matters and the above condition is usually separated into two classes, depending upon whether
$j_{1}-j_{2}$
is even or odd.
-
(i) For the class I,
$j_{1}-j_{2}=2m$
is even with
$j_{1}=m$
and
$j_{2}=-m$
, and it is found that for
$p>0$
the following possibilities exist (2.33)
$$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{m}^{+}(p)=\unicode[STIX]{x1D6FE}_{-m}^{-}(p)\quad m\geqslant 1.\end{eqnarray}$$
-
(ii) For the class II,
$j_{1}-j_{2}=2m+1$
is odd with
$j_{1}=m$
and
$j_{2}=-m-1$
, and it is found that for
$p>0$
the following possibilities exist (2.34)
$$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{m}^{+}(p)=\unicode[STIX]{x1D6FE}_{-m-1}^{-}(p)\quad m\geqslant 0.\end{eqnarray}$$
It is noteworthy that there are no collisions with
$m\geqslant 1$
between modes propagating in the same direction. The collisions of the eigenvalues can also be interpreted as resonances of two infinitesimal waves with the base flow, when considered from the ‘original’ fixed frame of reference. These resonant conditions can be written for the class I as
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}_{m}^{+}+\unicode[STIX]{x1D714}_{-m}^{-}=2m\unicode[STIX]{x1D714} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle k_{m}+k_{-m}=2mk & \displaystyle\end{eqnarray}$$
and for the class II as
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}_{m}^{+}+\unicode[STIX]{x1D714}_{-m-1}^{-}=(2m+1)\unicode[STIX]{x1D714} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle k_{m}+k_{-m-1}=(2m+1)k. & \displaystyle\end{eqnarray}$$
In this study we have considered the linear stability of the so-called forward mode, for which the wave speed always exceeds the free surface speed, hence
$c>0$
here. As the vorticity parameter varies, the resonances loci (2.33) and (2.34) change as illustrated in figure 1 for
$m=1,2$
and
$kh=10$
. For large positive values of
$\unicode[STIX]{x1D6FA}$
, the resonances loci correspond to integer values of
$p$
, which are associated with superharmonic disturbances. Whatever the type of disturbances, there is no coupling between the disturbances and the undisturbed flow at zero wave height. However, it appears that the linear resonant conditions (2.33) and (2.34) correspond to collisions of two eigenvalues with opposite signature or sign of excess total energy, when considered from the moving frame in which the basic state is a steady flow.
To see this we first compute in the rest frame the excess total energy, which is defined as the difference between the total energy of moving fluid with a wave perturbation and without the perturbation. From the analytical linear solutions for sinusoidal waves on a linear shear current, we can easily show that for a wave disturbance with (intrinsic) frequency
$\tilde{\unicode[STIX]{x1D714}}=\unicode[STIX]{x1D714}_{j}^{\pm }$
and wavenumber
$\unicode[STIX]{x1D705}=p+j$
, the excess total energy is given, up to second order, by
$$\begin{eqnarray}E=\frac{\unicode[STIX]{x1D70C}g\tilde{a}^{2}}{2\unicode[STIX]{x1D705}\text{tanh}\unicode[STIX]{x1D705}h}\tilde{\unicode[STIX]{x1D714}}\left(\tilde{\unicode[STIX]{x1D714}}+\frac{1}{2}\unicode[STIX]{x1D6FA}_{0}\text{tanh}\unicode[STIX]{x1D705}h\right),\end{eqnarray}$$
where
$\tilde{a}$
denotes the perturbation amplitude. Then, using the invariance of wave action
$A=E/\unicode[STIX]{x1D714}$
through Galilean transform we can deduce that, in a frame moving uniformly at a velocity
$c$
, the excess total energy is
$$\begin{eqnarray}E^{\prime }=\frac{\unicode[STIX]{x1D70C}g\tilde{a}^{2}}{2\text{tanh}\unicode[STIX]{x1D705}h}\left(\tilde{\unicode[STIX]{x1D714}}+\frac{1}{2}\unicode[STIX]{x1D6FA}_{0}\text{tanh}\unicode[STIX]{x1D705}h\right)\left(\frac{\tilde{\unicode[STIX]{x1D714}}}{\unicode[STIX]{x1D705}}-c\right).\end{eqnarray}$$
Equation (2.40) shows that the excess energy is proportional both to the pseudo frequency
$\tilde{\unicode[STIX]{x1D714}}+(\unicode[STIX]{x1D6FA}_{0}\text{tanh}\unicode[STIX]{x1D705}h)/2$
and the apparent frequency
$\tilde{\unicode[STIX]{x1D714}}-c\unicode[STIX]{x1D705}$
. For a copropagating perturbation, the pseudo frequency is always positive and we see that the excess energy may be negative if the relative frequency is negative, namely if the perturbation is moving slower than the linear basic wave. In contrast for co-propagating perturbations going faster (
$\tilde{\unicode[STIX]{x1D714}}/k>c>0$
) and for counter-propagating disturbances (
$\tilde{\unicode[STIX]{x1D714}}/k<0$
), the excess energy is always positive.
Figure 1. Resonances loci
$p$
as a function of
$\unicode[STIX]{x1D6FA}$
for
$\unicode[STIX]{x1D707}=10$
and different classes: ——, class I (
$m=1$
); – – –, class II (
$m=1$
); — ⋅ —, class I (
$m=2$
);
$\cdots \cdots$
, class II (
$m=2$
).
3 Numerical methods
The numerical calculations consist of two parts, determination of the unperturbed basic flow
$\bar{\unicode[STIX]{x1D702}},\bar{\unicode[STIX]{x1D713}}$
and subsequent solution of the eigenvalue problem as a function of the mode wavenumber
$p$
. In this section, we present the numerical methods which have been used in this study and the results of their validation against available numerical and analytical results.
3.1 Computation of steady waves
To calculate the unperturbed wave, the kinematic and dynamic surface conditions (2.4)–(2.5) are collocated at
$2N$
points equally distributed over one wavelength, though by symmetry only
$N+1$
points from wave crest to wave trough will be considered. Letting
$\unicode[STIX]{x1D702}_{i}=\unicode[STIX]{x1D702}(x_{i})$
where
$x_{i}=(i-1)\,\text{d}x$
,
$i=1,\ldots ,N+1$
and
$\text{d}x=\unicode[STIX]{x1D706}/2N$
, the above equations yield
$2N+2$
equations for the
$2N+4$
variables
$[\unicode[STIX]{x1D702}_{1},\ldots ,\unicode[STIX]{x1D702}_{N+1},B_{1},\ldots ,B_{N},c,R,Q]$
for a given value
$h$
of the mean depth. The required two additional equations are (2.6) and (2.7), which specify the wave parameters.
For numerical purposes, it is convenient to work with dimensionless variables. Hence by choosing the reference scales, as in the (irrotational) deep water case,
$$\begin{eqnarray}l_{ref}=\frac{1}{k},\quad t_{ref}=\frac{1}{\sqrt{gk}},\quad c_{ref}=\sqrt{\frac{g}{k}}\end{eqnarray}$$
the dimensional variables can be scaled as
$$\begin{eqnarray}\unicode[STIX]{x1D702}^{\prime }=k\unicode[STIX]{x1D702},\quad \unicode[STIX]{x1D713}^{\prime }=\sqrt{\frac{k^{3}}{g}}\unicode[STIX]{x1D713},\quad c^{\prime }=\sqrt{\frac{k}{g}}c,\quad R^{\prime }=\frac{k}{g}R,\quad Q^{\prime }=\sqrt{\frac{k^{3}}{g}}Q,\quad \unicode[STIX]{x1D6FA}^{\prime }=\frac{\unicode[STIX]{x1D6FA}_{0}}{\sqrt{gk}}\end{eqnarray}$$
and the dimensionless equations to be solved become, omitting the primes,
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FA}}{2}(\unicode[STIX]{x1D702}_{i}^{2}-\unicode[STIX]{x1D707}^{2})-c(\unicode[STIX]{x1D707}+\unicode[STIX]{x1D702}_{i})+\mathop{\sum }_{j=1}^{N}B_{j}\frac{\sinh j(\unicode[STIX]{x1D707}+\unicode[STIX]{x1D702}_{i})}{\cosh j\unicode[STIX]{x1D707}}\cos jkx_{i}+Q=0\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D702}_{i}+\frac{1}{2}\left(\mathop{\sum }_{j=1}^{N}jkB_{j}\frac{\sinh j(\unicode[STIX]{x1D707}+\unicode[STIX]{x1D702}_{i})}{\cosh j\unicode[STIX]{x1D707}}\sin jkx_{i}\right)^{2}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\frac{1}{2}\left(\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D702}_{i}-c+\mathop{\sum }_{j=1}^{N}jkB_{j}\frac{\cosh j(\unicode[STIX]{x1D707}+\unicode[STIX]{x1D702}_{i})}{\cosh j\unicode[STIX]{x1D707}}\cos jkx_{i}\right)^{2}-R=0\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{2}\unicode[STIX]{x1D702}_{1}+\frac{1}{2}\unicode[STIX]{x1D702}_{N+1}+\mathop{\sum }_{i\geqslant 2}\unicode[STIX]{x1D702}_{i}=0 & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}_{1}-\unicode[STIX]{x1D702}_{N+1}-2\unicode[STIX]{x1D700}=0. & \displaystyle\end{eqnarray}$$
Formally, the above system of
$2N+4$
nonlinear equations can be written as
$$\begin{eqnarray}F_{m}(Z)=0\quad m=1,\ldots ,2N+4,\end{eqnarray}$$
where
$Z=[\unicode[STIX]{x1D702}_{1},\ldots ,\unicode[STIX]{x1D702}_{N+1},B_{1},\ldots ,B_{N},c,R,Q]$
is the vector of arguments of
$F_{m}$
. When an initial approximation of the solution is known for given values of wave steepness
$\unicode[STIX]{x1D700}=Hk/2$
, dimensionless mean depth
$\unicode[STIX]{x1D707}=kh$
and dimensionless vorticity parameter
$\unicode[STIX]{x1D6FA}$
, the above system of equations can be solved iteratively using Newton’s method, which iterates with quadratic convergence towards the required nonlinear solution. The initial approximation to the solution is assumed to be a linear wave on constant vorticity. Alternatively, for large wave steepness or strong vorticity it may be necessary to use initially the results obtained from previous computations, provided those initial data are not too far from the target solution.
This numerical method has been validated by comparison with both the numerical results of Cokelet (Reference Cokelet1977) for irrotational steady waves (
$\unicode[STIX]{x1D6FA}=0$
) and the few (rare) results reported by Teles da Silva & Peregrine (Reference Teles da Silva and Peregrine1988) for finite-amplitude periodic waves on a linear shear current of finite depth. Tables 1 and 2 show the comparison for the squared wave velocity
$c^{2}$
and
$c$
, since these parameters have traditionally been used as the first basis for comparison between wave theories. A very good agreement is obtained between our numerical results and those reported by these authors.
Table 1. Comparison of results for
$c^{2}$
of irrotational waves. Results from Cokelet are the most accurate presented by this author.
Table 2. Comparison of
$c$
for different finite-amplitude periodic waves on a linear shear current with
$\unicode[STIX]{x1D707}=1$
.
Figure 2. Wave velocity
$c$
as a function of
$\unicode[STIX]{x1D700}$
for different values of
$\unicode[STIX]{x1D6FA}$
and
$\unicode[STIX]{x1D707}=2$
. Solid lines represent the numerical solution and the dashed lines represent the third-order solution.
Figure 3. Surface wave profiles for
$\unicode[STIX]{x1D700}=0.20$
,
$\unicode[STIX]{x1D6FA}=0.5$
and
$\unicode[STIX]{x1D707}=2$
. The solid line represents the numerical solution and the dashed line represents the third-order solution.
In order to check the numerical results in some broader range of parameter space, we have also compared our numerical results with the predictions from a third-order solution, derived recently by Hsu et al. (Reference Hsu, Francius, Montalvo and Kharif2016), to investigate analytically the effects of vorticity on the characteristics of Stokes gravity and capillary–gravity waves in finite depth. Figure 2 shows the comparison for the phase velocity, considering only the forward mode, as a function of the wave steepness and for different fixed values of the shear parameter. The weakly nonlinear results are plotted over the largest range of steepness that could be considered with our numerical method. It is interesting to notice that the weakly nonlinear predictions match well with the numerical results far beyond the weakly nonlinear regime, although the analytical and numerical wave profiles are not entirely coincident. This is illustrated in figure 3 for the surface profile of a periodic wave with
$\unicode[STIX]{x1D700}=0.20$
,
$\unicode[STIX]{x1D6FA}=0.5$
and
$\unicode[STIX]{x1D707}=2$
. It should be realized that although the relative error on
$c$
is small, as shown in figure 2, the errors on
$\unicode[STIX]{x1D702}$
and
$\unicode[STIX]{x1D719}$
can be much larger, as well as those on the Fourier coefficients. In figure 3, the numerical wave profiles have been obtained with
$N=32$
harmonics for the streamfunction.
In general, the performance of most numerical methods for solving exact nonlinear waves deteriorates as the wave gets steeper, and the present method cannot escape this fact. Indeed, the present approach fails owing to divergence of the (truncated) series for the streamfunction, when the computed solution is approaching a limiting configuration in the phase space of the governing dimensionless parameters; even though the interior flow is completely free of singularities and the free surface is smooth. For a given value of
$\unicode[STIX]{x1D707}$
, this occurs either with downstream propagating waves (when
$\unicode[STIX]{x1D6FA}>0$
) or with upstream propagating waves (when
$\unicode[STIX]{x1D6FA}<0$
). In the former case, the failure of convergence is reflected in the non-decrease for large
$N$
of the coefficients of the Fourier series of the streamfunction, because the analytic continuation of the streamfunction out of the physical domain develops singularities approaching the wave crest from above. In the latter case, the method also fails before the surfaces become multivalued, i.e. overhanging waves develop, which is reflected in oscillatory behaviour of the Fourier coefficients in which case Fourier expansions are not appropriate. For any sign of the vorticity, the present method is also expected to break down in the limit of very long waves, when a large number of modes has to be taken.
Despite these limitations, it should stressed that when the present method converges, it provides highly accurate results, even with a low number of modes, which may be readily used for computation of the interior flow variables such as, for instance, fluid velocities and pressure which are required in practical problems. This contrasts with many other numerical methods, such as those based on surface integro-differential equations (Vanden-Broeck Reference Vanden-Broeck1996; Ashton & Fokas Reference Ashton and Fokas2011) or on the conformal mapping technique (Choi Reference Choi2009; Ribeiro et al. Reference Ribeiro, Milewski and Nachbin2017), which require the results to be mapped back to the physical domain.
We shall conclude this section by noticing that the present method also enables us to compute waves on a linear shear current when an interior stagnation point and a recirculating region exist, either with or without the pressure anomalies first mentioned in Teles da Silva & Peregrine (Reference Teles da Silva and Peregrine1988). To illustrate this, we show in figure 4 the streamline patterns for three different periodic waves with
$\unicode[STIX]{x1D6FA}=-2$
and
$kh=1$
. For the sake of visibility, the spacing between the streamlines is constant in each region separated by the separatrix streamline (
$\unicode[STIX]{x1D6F9}=0$
), but the flow rate between two consecutive streamlines differs in each region. For each wave the ratio between the flow rates is fixed to plot eight streamlines inside the recirculating region. As can be seen, the flow patterns, reported in figure 4 only in a half-wavelength interval, correspond to waves with three stagnation points: two saddles at the bottom and a centre under the crest. As the steepness increases, the centre of the separated eddy is shifted upwards, its vertical extent increases and its horizontal extent at the bottom is decreased. Note that here pressure anomalies exist only for the highest wave (not shown). According to the classification proposed very recently by Ribeiro et al. (Reference Ribeiro, Milewski and Nachbin2017), our results fall in ‘region 2’ of their paper (see their figure 5).
Figure 4. Streamline patterns for periodic waves with
$\unicode[STIX]{x1D6FA}=-2$
and
$\unicode[STIX]{x1D707}=1$
: (a)
$\unicode[STIX]{x1D700}=0.05$
for which
$c=1.9207810276$
, (b)
$\unicode[STIX]{x1D700}=0.25$
for which
$c=1.9427455675$
and (c)
$\unicode[STIX]{x1D700}=0.45$
for which
$c=1.998$
.
3.2 Computation of eigenvalues
Once the unperturbed wave has been computed, the normal-mode decompositions (2.19)–(2.21) are truncated at
$N$
Fourier modes to solve the eigenvalue problem. Like McLean (Reference McLean1982) we use a collocation method with
$2N+1$
grid points, equally distributed between two adjacent crests of the unperturbed wave, to transform the infinite-dimensional eigenvalue problem (2.23)–(2.24) into a discrete system of order
$4N+2$
, which can be written as
$$\begin{eqnarray}\unicode[STIX]{x1D63C}\boldsymbol{v}=\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D63D}\boldsymbol{v},\end{eqnarray}$$
where
$\boldsymbol{v}=(\boldsymbol{a}_{-\boldsymbol{N}},\ldots ,\boldsymbol{a}_{\boldsymbol{N}},\boldsymbol{b}_{-\boldsymbol{N}},\ldots ,\boldsymbol{b}_{\boldsymbol{N}})^{\boldsymbol{t}}$
. Here the matrices
$\unicode[STIX]{x1D63C}$
and
$\unicode[STIX]{x1D63D}$
are complex functions of the mode
$p$
and the unperturbed wave that is specified itself by
$\unicode[STIX]{x1D700}$
,
$kh$
and
$\unicode[STIX]{x1D6FA}$
. For each disturbance with mode number
$p$
, the eigenvalues
$\unicode[STIX]{x1D6FE}_{n}$
are obtained with an eigenvalue solver based on the QZ algorithm, as well as their associated eigenvector
$(a_{nj},b_{nj})^{t}$
, for
$n=1,\ldots ,4N+2$
and
$j=-N,\ldots ,N$
. Instability corresponds to
$\text{Im}(\unicode[STIX]{x1D6FE}_{n})>0$
for at least one mode
$n$
, given the values of
$p$
and the unperturbed wave. In computations of the linear stability of the unperturbed wave, we choose to set
$0\leqslant p\leqslant 0.5$
, owing to the symmetries of the eigenvalue problem and the degeneracy of the normal-mode decomposition in regard to the choice of
$p$
. In practice, the effects of modal truncation must be monitored by increasing
$N$
until the relevant eigenvalues have converged. For all the cases considered in the present work, our numerical results have been accepted when the convergence of the eigenvalues has been reached with at least six significant digits for the more difficult cases.
Our numerical method for the truncated version of the eigenvalue problem (2.23)–(2.24) has been checked against the analytical results of Thomas et al. (Reference Thomas, Kharif and Manna2012) for the quartet resonant interactions in weakly nonlinear gravity waves in flows with constant vorticity. In fact, preliminary checks of our numerical implementation of the present method has been done with other strictly numerical results of McLean (Reference McLean1982) and Francius & Kharif (Reference Francius and Kharif2006) for two-dimensional quartet interactions between irrotational waves on finite depth (
$\unicode[STIX]{x1D6FA}=0$
). Since the focus in the present study is on the effect of vorticity on two-dimensional perturbations, the details of this comparison is not shown here. It suffices to say that the comparison with the data shows a very good agreement, provided the steady waves computed with the present method are obtained with sufficiently high accuracy.
Instead, we focus on the comparison with the weakly nonlinear theory that serves in fact two purposes. First, it is a mean to check our numerical results when
$\unicode[STIX]{x1D6FA}\neq 0$
, since we have not found in the literature useful quantitative results for the 2-D instabilities due to quartet resonances; secondly, the numerical results may confirm some weakly nonlinear predictions such as, for instance, the complete restabilization to two-dimensional long-wavelength perturbations, which occurs for sufficiently large (constant) vorticity and
$kh\geqslant 1.363$
. This important analytical result will be discussed in the next section.
Figure 5 shows the comparison of the dimensionless growth rates of the two-dimensional quartet instabilities varying with
$p$
for
$\unicode[STIX]{x1D700}=0.05$
,
$\unicode[STIX]{x1D6FA}=0,\pm 0.2,\pm 0.4$
and two values of the dimensionless depth
$\unicode[STIX]{x1D707}=10,2$
. In either case an increase of positive shear yields a wider range of unstable wavenumbers and enhanced growth rates. For
$\unicode[STIX]{x1D700}=0.01$
, the numerical results agree very well with the analytical results and the corresponding curves are indiscernible (not shown here). From figure 5, we see that the increase of wave steepness increases the difference between the analytical and the numerical results, although the variation of this difference with
$p$
depends on the value of the dimensionless depth (for any given steepness). Figure 5 also shows that the linear growth rates of these quartet instabilities decrease when
$\unicode[STIX]{x1D707}$
decreases. As is well known for irrotational waves, the complete restabilization to two-dimensional long-wavelength perturbations occurs for the critical value
$\unicode[STIX]{x1D707}\approx 1.363$
. It should also be noticed that for
$p\ll 1$
, namely for very long-wavelength perturbations, we have always found an excellent agreement between the numerical results and the analytical results of these two-dimensional instabilities.
Figure 5. Growth rate
$\text{Im}(\unicode[STIX]{x1D6FE})$
as a function of
$p$
for waves with
$\unicode[STIX]{x1D700}=0.05$
,
$\unicode[STIX]{x1D707}=10$
(a) and
$\unicode[STIX]{x1D707}=2$
(b), and different shear values: ——
$\unicode[STIX]{x1D6FA}=0$
; — —,
$\unicode[STIX]{x1D6FA}=\pm 0.2$
; — ⋅ —,
$\unicode[STIX]{x1D6FA}=\pm 0.4$
. Dotted lines represent the corresponding weakly nonlinear results. The curves above the solid line are for
$\unicode[STIX]{x1D6FA}>0$
.
4 Results
The stability of the periodic waves was studied in details for two values of the dimensionless mean depth,
$\unicode[STIX]{x1D707}=10$
and
$\unicode[STIX]{x1D707}=2$
, for a range of values of the vorticity parameter
$\unicode[STIX]{x1D6FA}$
and several values of the wave steepness
$\unicode[STIX]{x1D700}$
, namely
$\unicode[STIX]{x1D700}=0.01,0.05,0.10,0.15,0.20$
. Actually, for each value of
$\unicode[STIX]{x1D707}$
and
$\unicode[STIX]{x1D716}$
considered here, our numerical computations of the unperturbed waves cover a large range of shear values in which
$\unicode[STIX]{x1D6FA}$
varies from
$-2$
to the maximal value allowed by the other parameters. In order to facilitate the presentation of our numerical results and to analyse the effects of vorticity on the characteristics of the two-dimensional instabilities under investigation, we have decided to select only a few values of
$\unicode[STIX]{x1D700}$
.
First, we present the numerical results for the quartet instabilities of the modulational type, namely due to sideband wave disturbances (
$0<p<1$
), and compare them with the results of the approximate analysis of Thomas et al. (Reference Thomas, Kharif and Manna2012), which has shown the disappearance of the Benjamin–Feir modulational instability due to the increase of the shear strength. In particular, we focus on the characteristics of the most unstable modes in order to check the validity of the analytical predictions and also to confirm the restabilization of the modulational instability both for upstream and downstream propagating waves. Secondly, we present preliminary numerical results for disturbances that are not of the modulational type, namely for perturbations of class I (
$m=1$
) with
$p>1$
and class II (
$m=1$
) with
$p>1$
, and discuss their importance as the vorticity is varied. Finally, the results showing the restabilization of uniform wave train to modulational disturbances are discussed. Emphasis is placed, in particular for deep water waves, on the differences between the results from the weakly nonlinear theories and those obtained by numerical calculations.
4.1 Maximum growth rates of the modulational instability
Before we show the numerical results for the most unstable disturbances associated with the quartet resonances, we briefly recall here some analytical results of Thomas et al. (Reference Thomas, Kharif and Manna2012) which are used in the present section. As already mentioned in the introduction, these authors have derived recently a NLS equation when there is a background linear shear flow. This equation is given for the complex envelope amplitude
$a(\unicode[STIX]{x1D709},\unicode[STIX]{x1D70F})$
$$\begin{eqnarray}\text{i}a_{\unicode[STIX]{x1D70F}}+La_{\unicode[STIX]{x1D709}\unicode[STIX]{x1D709}}=M|a|^{2}a,\end{eqnarray}$$
where
$\unicode[STIX]{x1D709}=\unicode[STIX]{x1D716}(x-c_{g}t)$
and
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D716}^{2}t$
are coordinates to describe the slow modulation in space and time of a unidirectional wave train with a carrier wave with wavenumber
$k$
and frequency
$\unicode[STIX]{x1D714}$
,
$c_{g}$
is the group velocity of the carrier wave. In this framework,
$\unicode[STIX]{x1D700}$
represents a small dimensionless amplitude parameter and the surface elevation
$z=\unicode[STIX]{x1D702}(x,t)$
can be written, to the leading order, as
$$\begin{eqnarray}\unicode[STIX]{x1D702}(x,t)={\textstyle \frac{1}{2}}\unicode[STIX]{x1D716}a(\unicode[STIX]{x1D709},\unicode[STIX]{x1D70F})\text{exp}[\text{i}(kx-\unicode[STIX]{x1D714}t)]+\text{c.c.}+O(\unicode[STIX]{x1D716}^{2}),\end{eqnarray}$$
where
$\text{c.c.}$
represents the complex conjugate.
As is well known for the NLS equation (4.1), the sign of the product of the dispersive coefficient
$L$
and the nonlinear coefficient
$M$
govern drastically the amplitude modulations of weakly nonlinear wave packets, as well as those of initially uniform wave trains. In the defocussing regime (
$LM>0$
) the plane wave solution of (4.1), which is given by
$$\begin{eqnarray}a(\unicode[STIX]{x1D709},\unicode[STIX]{x1D70F})=a_{0}\text{e}^{-\text{i}Ma_{0}^{2}\unicode[STIX]{x1D70F}},\end{eqnarray}$$
where
$a_{0}$
is an arbitrary constant, is always modulationally stable to spatial periodic perturbations whose real and imaginary parts are proportional to
$\text{e}^{\text{i}K\unicode[STIX]{x1D709}}$
, where
$\pm K$
represents perturbations of the fundamental wavenumber
$k$
. In contrast, in the focusing regime (
$LM<0$
) the plane wave solution (4.3) is modulationally unstable whenever
$K$
lies in a finite band,
$0<K<K_{c}$
with
$K_{c}=\sqrt{2|M/L|}a_{0}$
. Over this range of unstable wavenumbers
$K$
, the maximum growth rate of the quartet instability is
$$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{max}=M_{1}\unicode[STIX]{x1D714}a_{0}^{2}k^{2}\end{eqnarray}$$
when
$$\begin{eqnarray}K_{max}=\sqrt{|M_{1}/L_{1}|}a_{0}k^{2},\end{eqnarray}$$
where
$L_{1}=k^{2}L/\unicode[STIX]{x1D714}$
and
$M_{1}=M/(\unicode[STIX]{x1D714}k^{2})$
. To correct a misprint in the definition of
$L$
in Thomas et al. (Reference Thomas, Kharif and Manna2012), the expressions for the dispersive and nonlinear coefficients are given here,
$$\begin{eqnarray}\displaystyle & \displaystyle L=\frac{\unicode[STIX]{x1D714}}{k^{2}\unicode[STIX]{x1D70E}(2+X)}[\unicode[STIX]{x1D707}(1-\unicode[STIX]{x1D70E}^{2})[1-\unicode[STIX]{x1D707}\unicode[STIX]{x1D70E}+(1-\unicode[STIX]{x1D70C})X]-\unicode[STIX]{x1D70E}\unicode[STIX]{x1D70C}^{2}] & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle M=\frac{-\unicode[STIX]{x1D714}k^{2}(U+VW)}{8(1+X)(2+X)\unicode[STIX]{x1D70E}^{4}} & \displaystyle\end{eqnarray}$$
respectively, where
$\unicode[STIX]{x1D70C}=c_{g}/c_{p}$
denotes the ratio of the group velocity to the phase velocity of the carrier wave,
$\unicode[STIX]{x1D707}=kh$
is the dimensionless water depth,
$\unicode[STIX]{x1D70E}=\text{tanh}(\unicode[STIX]{x1D707})$
and
$X=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FA}_{0}/\unicode[STIX]{x1D714}$
where
$\unicode[STIX]{x1D6FA}_{0}$
represents the dimensional shear. The lengthy expressions for
$U$
,
$V$
and
$W$
can be found in Thomas et al. (Reference Thomas, Kharif and Manna2012).
Figure 6. (a) Wavenumber
$p_{max}$
of the most unstable mode and (b) its growth rate
$\unicode[STIX]{x1D6FE}_{max}$
as a function of
$\unicode[STIX]{x1D6FA}$
, for
$\unicode[STIX]{x1D707}=10$
and
$\unicode[STIX]{x1D700}=0.01$
(stars),
$0.05$
(asterisks),
$0.10$
(diamonds),
$0.15$
(squares),
$0.20$
(circles); solid lines represent the weakly nonlinear results.
The most unstable modes with dimensionless wavenumber
$p_{max}=K_{max}/k$
in our notation, as well as their growth rates, were obtained numerically by solving the truncated version of the eigenvalue problem (2.23)–(2.24), and are shown in figure 6 as a function of the dimensionless vorticity,
$\unicode[STIX]{x1D6FA}$
, for
$\unicode[STIX]{x1D707}=10$
and
$\unicode[STIX]{x1D700}=0.01,0.05,0.10,0.15,0.20$
. For
$-2\leqslant \unicode[STIX]{x1D6FA}\leqslant 1$
and the two smallest values of
$\unicode[STIX]{x1D700}$
, a good quantitative agreement between our numerical results and the analytical predictions is observed. In fact, for
$\unicode[STIX]{x1D700}=0.01$
, we found some differences between numerical and analytical growth rates when
$2\leqslant \unicode[STIX]{x1D6FA}\leqslant 4$
(not shown here), the latter overestimating the former. For the steepest waves considered here,
$\unicode[STIX]{x1D700}\geqslant 0.15$
, the analytical predictions are found to overestimate the characteristics of the most unstable modes, namely
$p_{max}$
and
$\unicode[STIX]{x1D6FE}_{max}$
, although the variations of these characteristics with
$\unicode[STIX]{x1D6FA}$
are consistent between the two results.
Figure 7. (a) Wavenumber
$p_{max}$
of the most unstable mode and (b) its growth rate
$\unicode[STIX]{x1D6FE}_{max}$
as a function of
$\unicode[STIX]{x1D6FA}$
, for
$\unicode[STIX]{x1D707}=2$
and
$\unicode[STIX]{x1D700}=0.01$
(stars),
$0.05$
(diamonds),
$0.10$
(squares),
$0.15$
(circles); solid lines represent the weakly nonlinear results.
For
$\unicode[STIX]{x1D707}=2$
the same analysis was carried out for
$\unicode[STIX]{x1D700}=0.01,0.05,0.10,0.15$
, and the results are shown in figure 7. The agreement between the two results is remarkably very good for negative values of
$\unicode[STIX]{x1D6FA}$
(upstream propagating waves), for both
$p_{max}$
and
$\unicode[STIX]{x1D6FE}_{max}$
, and this even for the steepest waves considered here. However, the two results show a tendency to diverge with increasing (positive) values of
$\unicode[STIX]{x1D6FA}$
, except for the waves with
$\unicode[STIX]{x1D700}=0.01$
. In this case, interestingly, a very good agreement between the two results is observed in the range of vorticity
$-2\leqslant \unicode[STIX]{x1D6FA}\leqslant 2$
.
Both figures 6 and 7 confirm one prominent feature of the effect of vorticity, namely the restabilization of the modulational instability for upstream propagating waves (
$\unicode[STIX]{x1D6FA}<0$
) at large enough negative shear. For small steepness, say
$\unicode[STIX]{x1D700}=0.01$
, the numerical results agree very well with those of the weakly nonlinear theory. For downstream propagating waves (
$\unicode[STIX]{x1D6FA}>0$
) the analytical results of Thomas et al. (Reference Thomas, Kharif and Manna2012) also predict restabilization of the modulational instability in finite depth, when
$\unicode[STIX]{x1D6FA}$
increases above a positive threshold value that depends only on
$\unicode[STIX]{x1D707}$
(in the limit
$p\rightarrow 0$
). As shown in figure 7 our numerical results for
$\unicode[STIX]{x1D700}=0.01$
also confirm this analytical prediction. As already mentioned for downstream propagating waves, the differences between the analytical and the numerical results increase with increasing nonlinearity, and there is poor qualitative agreement. For upstream propagating waves, however, it should be noticed that the negative value
$\unicode[STIX]{x1D6FA}_{c}$
, below which the modulational instability restabilizes, does not seem to be affected by the increase of nonlinearity of the basic wave.
It is emphasized here that the disappearance of the modulational instability for upstream propagating waves has been predicted very recently by Thomas et al. (Reference Thomas, Kharif and Manna2012), but this has not been confirmed with an (alternative) independent method. For this reason we have pursued our analysis by carrying out computations for long-wavelength perturbations, with the purpose of checking further the weakly nonlinear theory. In the limits as
$p\rightarrow 0$
and
$\unicode[STIX]{x1D700}\rightarrow 0$
, this theory predicts for
$\unicode[STIX]{x1D707}>1.363$
that the stability boundary is
$-1<X\leqslant X_{c1}<0$
for upstream propagating waves and
$0<X_{c2}\leqslant X$
for downstream propagating waves (see figure 3 of Thomas et al.
Reference Thomas, Kharif and Manna2012). Here
$X=\unicode[STIX]{x1D70E}\overline{\unicode[STIX]{x1D6FA}}$
with
$\overline{\unicode[STIX]{x1D6FA}}=\unicode[STIX]{x1D6FA}_{0}/\unicode[STIX]{x1D714}$
, which is related to
$\unicode[STIX]{x1D6FA}$
by the relation
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}=\frac{\unicode[STIX]{x1D70E}^{1/2}\overline{\unicode[STIX]{x1D6FA}}}{\sqrt{1+\unicode[STIX]{x1D70E}\overline{\unicode[STIX]{x1D6FA}}}}.\end{eqnarray}$$
The critical values
$X_{c1}$
and
$X_{c2}$
depend only on
$\unicode[STIX]{x1D707}$
and
$\unicode[STIX]{x1D6FA}$
. They can be determined either analytically or graphically, as we have done, by plotting for a fixed
$\unicode[STIX]{x1D707}$
the variations of the (scaled) nonlinear coefficient
$M_{1}$
with
$\unicode[STIX]{x1D6FA}$
and finding the critical values for which it vanishes. As discussed in Thomas et al. (Reference Thomas, Kharif and Manna2012), in the deep water limit
$\unicode[STIX]{x1D70E}=1$
,
$X_{c1}=-2/3$
and
$X_{c2}$
increases without bound.
We considered waves with
$\unicode[STIX]{x1D707}=10$
and
$\unicode[STIX]{x1D700}=0.01$
and the growth rates of long-wavelength perturbations were obtained numerically with a few very small values of
$p$
, and the shear parameter
$\unicode[STIX]{x1D6FA}$
was increased from
$-2$
to
$1$
in steps of
$0.01$
. The numerical results, shown in figure 8, indicate that the restabilization of the modulational instability occurs at the threshold value,
$X_{\ast }=-0.657$
(
$\unicode[STIX]{x1D6FA}=-1.12$
), in close agreement with the theoretical value
$X_{c}\approx -0.654$
(
$\unicode[STIX]{x1D6FA}_{c}=-1.11$
). As indicated by the analytical results the restabilization of the modulational instability occurs when the range of unstable wavenumbers has shrunk to zero and the stability boundaries have merged at
$p=0$
.
Figure 8. Restabilization of modulational perturbations for waves propagating upstream. Growth rate versus
$X=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FA}_{0}/\unicode[STIX]{x1D714}$
for different wavenumbers (
$p\ll 1$
) with
$\unicode[STIX]{x1D700}=0.01$
and
$\unicode[STIX]{x1D707}=10$
.
A similar restabilization occurs with downstream propagating waves (
$\unicode[STIX]{x1D6FA}>0$
), but as the shear parameter
$X$
is increased above the positive threshold value
$X_{c2}$
. Although we have not addressed this issue in detail, namely by comparing the analytical results found by Thomas et al. (Reference Thomas, Kharif and Manna2012) with the results of the present strictly numerical method, we have compared the marginal stability curve obtained with the analytic solutions of the NLS equation with vorticity from Thomas et al. (Reference Thomas, Kharif and Manna2012) with that obtained by Hur & Johnson (Reference Hur and Johnson2015) (see their equation (3.3) in Theorem 3.1), in order to comply with the demand of one referee. As the issue of concern here is that of the linear stability of the forward mode, in particular to long-wavelength perturbations, it is found that the predictions of Hur & Johnson (Reference Hur and Johnson2015) match perfectly the results of Thomas et al. (Reference Thomas, Kharif and Manna2012) only for downstream propagating waves (
$\unicode[STIX]{x1D6FA}>0$
) and, say approximately,
$\unicode[STIX]{x1D6FA}\geqslant 1.5$
. For smaller and decreasing values of
$\unicode[STIX]{x1D6FA}$
, the differences between the two results increases without bounds. More importantly, both analytical results are in total contradiction for upstream propagating waves (
$\unicode[STIX]{x1D6FA}<0$
). In this case, the model of Hur & Johnson (Reference Hur and Johnson2015) does not predict the restabilization of the modulational instability whatever the value of the depth, provided the shear is sufficiently large. Although interesting, this comparison will be reported elsewhere in detail since in the present work we have not studied the marginal stability curve for the modulational instability in the (
$\unicode[STIX]{x1D6FA},\unicode[STIX]{x1D707}$
) parameter space.
4.2 Growth rates of two-dimensional instabilities
In our investigations we have also detected other bands of instabilities, which are associated with the quartet resonances with
$p>1$
, i.e. not of the modulational type, and the quintet resonances which correspond to instabilities of class II (
$m=1$
). These instabilities were identified by inspecting the dominant components of the eigenvectors associated with their complex eigenvalues. In some region of the parameter space we have observed that the maximum growth rates of some of these instabilities are larger than those of quartet instabilities with
$p>1$
. For this reason, we present here some preliminary results which were obtained for
$\unicode[STIX]{x1D707}=2$
,
$\unicode[STIX]{x1D700}=0.05$
and several values of the shear parameter
$\unicode[STIX]{x1D6FA}=0,\pm 0.2,\pm 0.4,\pm 0.8$
.
Figure 9. Growth rate of quartet instabilities, class I (
$m=1$
) with
$p<1$
, for
$\unicode[STIX]{x1D707}=2$
,
$\unicode[STIX]{x1D700}=0.05$
and for (a) negative and (b) positive values of the shear: ——,
$\unicode[STIX]{x1D6FA}=0$
; – – –,
$|\unicode[STIX]{x1D6FA}|=0.2$
; — ⋅ —,
$|\unicode[STIX]{x1D6FA}|=0.4$
;
$\cdots \cdots$
,
$|\unicode[STIX]{x1D6FA}|=0.8$
.
For comparison of these results with those for the quartet instabilities near
$p=0$
, we show in figure 9 the growth rates of these modulational instabilities versus
$p$
for the different values of
$\unicode[STIX]{x1D6FA}$
considered here. As we have already seen, the class I (
$m=1$
) instability region with
$p<1$
is increased (decreased) as
$\unicode[STIX]{x1D6FA}$
increases (decreases).
Figure 10. Growth rate of quartet instabilities, class I (
$m=1$
) with
$p>1$
, for
$\unicode[STIX]{x1D707}=2$
and
$\unicode[STIX]{x1D700}=0.05$
; (a) for
$\unicode[STIX]{x1D6FA}=0$
; (b)
$-0.2$
; (c)
$-0.4$
; (d)
$-0.8$
. The dashed line indicates the corresponding value of the maximum growth rate of the class I (
$m=1$
) instabilities with
$p<1$
.
Figure 11. Growth rate of quintet instabilities, class II (
$m=1$
) with
$p>1$
, for
$\unicode[STIX]{x1D707}=2$
and
$\unicode[STIX]{x1D700}=0.05$
; (a) for
$\unicode[STIX]{x1D6FA}=0$
; (b)
$-0.2$
; (c)
$-0.4$
; (d)
$-0.8$
; (dashed line) as in figure 10.
Figure 12. Growth rate of quartet instabilities, class I (
$m=1$
) with
$p>1$
, for
$\unicode[STIX]{x1D707}=2$
and
$\unicode[STIX]{x1D700}=0.05$
; (a) for
$\unicode[STIX]{x1D6FA}=0$
; (b)
$0.2$
; (c)
$0.4$
; (d)
$0.8$
; (dashed line) as in figure 10.
Figure 13. Growth rate of quintet instabilities, class II (
$m=1$
) with
$p>1$
, for
$\unicode[STIX]{x1D707}=2$
and
$\unicode[STIX]{x1D700}=0.05$
; (a) for
$\unicode[STIX]{x1D6FA}=0$
; (b)
$0.2$
; (c)
$0.4$
; (d)
$0.8$
; (dashed line) as in figure 10.
Figures 10 and 11 show, for upstream propagating waves, the effects of increasing the shear on the instability bands associated with class I (
$m=1$
) and class II (
$m=1$
) disturbances with
$p>1$
, respectively. As the shear increases, both instability regions moves to the right and, although initially their bandwidths and the maximum growth rates diminish slightly, these increase for sufficiently large negative values of
$\unicode[STIX]{x1D6FA}$
. For
$\unicode[STIX]{x1D6FA}=-0.8$
, figures 10(d) and 11(d) reveal that the maximum growth rate of these two-dimensional instabilities is larger than both that of the most unstable mode of the class I (
$m=1$
) with
$p<1$
and that of the most unstable mode of the class II (
$m=1$
) with
$p>1$
. In each figure the dashed line indicates the corresponding value of the maximum growth rate of the quartet instabilities with
$p<1$
, namely of the modulational type.
Comparing figure 10(c,d) reveals that this class I (
$m=1$
) instability region must include
$p=2$
for some values
$-0.8<\unicode[STIX]{x1D6FA}<-0.4$
. When this occurs, the colliding modes have dominant wavenumbers
$p+m=2+1=3$
and
$p-m=2-1=1$
, which correspond to superharmonic instabilities of the basic wave. This result suggests that this two-dimensional superharmonic instability, which appears due to the effect of vorticity, can play an important role (yet to be explored) in the dynamics of weakly modulated wave trains on a linear shear current.
It should be noticed that the results in figure 11 also indicate three occurrences of superharmonic instabilities, as the shear is varied in the range
$-0.8<\unicode[STIX]{x1D6FA}<0$
. The first occurrence is due to collision of the modes
$p+m=3+1=4$
and
$p-m-1=3-2=1$
, the second one with
$p+m=4+1=5$
and
$p-m-1=4-2=2$
and the third with
$p+m=5+1=6$
and
$p-m-1=5-2=3$
. Although these superharmonic instabilities of class II (
$m=1$
) do not generally correspond to the (dominant) maximum growth rate, they appear together with the most unstable (subharmonic) disturbances and it would be interesting to study the effects of increasing nonlinearity on their characteristics. This point could be explored in future studies.
Figures 12 and 13 show, for downstream propagating waves, the effects of increasing the shear on the instability bands associated with class I (
$m=1$
) and class II (
$m=1$
) disturbances with
$p>1$
, respectively. The effects of vorticity are different and seem to be weaker for the waves considered here. As the shear increases, both instability regions moves to the right and, their bandwidths and the maximum growth rates have slightly changed in the range
$0<\unicode[STIX]{x1D6FA}<0.8$
. In contrast with the results obtained for upstream waves, no superharmonic instabilities are detected. It is noticed, however, that the bandwidth of the class I (
$m=1$
) instability region tends to shrink to zero, and also that the whole instability region approaches
$p=1$
, as shown in figure 12(d).
4.3 Discussion
The results reported in § 4.1 appear to confirm the analytical results of Thomas et al. (Reference Thomas, Kharif and Manna2012). These constitute somewhat numerical evidence of this important property of the modulational instability, namely its restabilization at sufficiently large values of the vorticities.
For small wave steepness, it appears that the sign of the vorticity does not matter in this restabilization of the basic wave, although its effects on the growth rates are of primary importance. For
$\unicode[STIX]{x1D700}=0.01$
and
$\unicode[STIX]{x1D707}=10$
the restabilization is observed only with upstream propagating waves, while for
$\unicode[STIX]{x1D707}=2$
it occurs both with upstream and downstream propagating waves, in agreement with the analytical results of Thomas et al. (Reference Thomas, Kharif and Manna2012).
It can be seen, from figures 6 and 7, that for upstream propagating waves on a linear shear current, the effects of vorticity become dominant over those of nonlinearity as the shear is increased. For
$\unicode[STIX]{x1D6FA}<0$
and wave steepness as high as
$\unicode[STIX]{x1D700}=0.20$
, these results show a very good comparison of the analytical results with our numerical solutions, and therefore suggest that the NLS model could be a powerful tool for predicting qualitative features of modulated upstream waves on a strong linear shear current. In particular, correct predictions of the most unstable modes of class I (
$m=1$
) are recognized, as well as the instability ranges. For the case of downstream propagating waves (
$\unicode[STIX]{x1D6FA}>0$
) and large shear, in contrast, the analytical results become quickly invalid as the nonlinearity is increased. A good agreement can be expected, however, provided the steepness is not too large, as shown in figure 7(b) for
$\unicode[STIX]{x1D707}=2$
and
$\unicode[STIX]{x1D700}=0.01$
. More importantly, the analytically predicted restabilization of the modulational instability is confirmed both for upstream and downstream propagating waves by the present numerical results.
Although there has been preliminary analytical evidence for waves in deep water, from Li et al. (Reference Li, Hui, Donelan, Horikawa and Maruo1987) using a linear shear current and Baumstein (Reference Baumstein1998) using a piecewise-linear current, suggesting that weak velocity shear can enhance modulational instability and strong velocity shear can suppress instability, it should be emphasized here that the predictions of these previous analytical studies are not supported by the present numerical results and, thus, by the weakly nonlinear theory of Thomas et al. (Reference Thomas, Kharif and Manna2012) in the limit
$\unicode[STIX]{x1D707}\rightarrow \infty$
(
$\unicode[STIX]{x1D70E}\rightarrow 1$
). As explained below, the source of disagreement can be immediately understood when one realizes that in deep water the wave-induced mean flow is a key feature in the slow spatial modulations of weakly nonlinear waves.
Inspection of previous analytical studies reveals that although Li et al. (Reference Li, Hui, Donelan, Horikawa and Maruo1987) and Baumstein (Reference Baumstein1998) applied the same method of multiple scales to the governing equations, like Thomas et al. (Reference Thomas, Kharif and Manna2012), the wave-induced mean-flow component of the solution was set arbitrarily to zero in their derivation of the modulation equation. At first sight this assumption sounds like a good idea for small velocity shear and wave steepness, if we keep in mind that for irrotational weakly modulated wave trains in deep water the feedback from this wave-induced mean flow is due to higher-order effects and do not appear at the order of approximation of the NLS equation.
As shown in Thomas et al. (Reference Thomas, Kharif and Manna2012), however, special care should be taken in the limit
$\unicode[STIX]{x1D707}\rightarrow \infty$
for the NLS equation presented in this work. In this limit they obtained
$$\begin{eqnarray}\text{i}a_{\unicode[STIX]{x1D70F}}+L_{\infty }a_{\unicode[STIX]{x1D709}\unicode[STIX]{x1D709}}=M_{\infty }|a|^{2}a,\end{eqnarray}$$
with
$$\begin{eqnarray}L_{\infty }=-\frac{\unicode[STIX]{x1D714}}{k^{2}}\frac{(1+\overline{\unicode[STIX]{x1D6FA}})^{2}}{(2+\overline{\unicode[STIX]{x1D6FA}})^{3}},\quad M_{\infty }=\frac{\unicode[STIX]{x1D714}k^{2}}{8}\left[(4+6\overline{\unicode[STIX]{x1D6FA}}+6\overline{\unicode[STIX]{x1D6FA}}^{2}+\overline{\unicode[STIX]{x1D6FA}}^{3})-\overline{\unicode[STIX]{x1D6FA}}^{2}\frac{(2+\overline{\unicode[STIX]{x1D6FA}})^{2}}{(1+\overline{\unicode[STIX]{x1D6FA}})}\right].\end{eqnarray}$$
As emphasized in Thomas et al. (Reference Thomas, Kharif and Manna2012), even though the wave-induced mean velocity approaches zero in the deep water limit, its interaction with the carrier wave yields an
$O(\unicode[STIX]{x1D716}^{3})$
term that has a non-zero finite limit when
$\unicode[STIX]{x1D707}\rightarrow \infty$
. For one-dimensional modulations this additional term with the same periodicity as the fundamental wave is directly proportional to
$|a|^{2}a$
. More precisely, it corresponds to the second term in the brackets defining the nonlinear coefficient
$M_{\infty }$
, which cancels only without background vorticity (as expected). When this term is neglected, the expression for
$M_{\infty }$
agrees with that of Baumstein (Reference Baumstein1998), after adoption of our notation. As noticed by Baumstein (Reference Baumstein1998), Li et al. (Reference Li, Hui, Donelan, Horikawa and Maruo1987) obtained a somewhat different formula for the nonlinear coefficient
$M_{\infty }$
which, in our notation, can be expressed
$$\begin{eqnarray}M_{Li}=\frac{\unicode[STIX]{x1D714}k^{2}}{8(2+\overline{\unicode[STIX]{x1D6FA}})}(8+16\overline{\unicode[STIX]{x1D6FA}}+10\overline{\unicode[STIX]{x1D6FA}}^{2}-\overline{\unicode[STIX]{x1D6FA}}^{4}).\end{eqnarray}$$
This should be compared with (4.10), which can also be expressed as
$$\begin{eqnarray}M_{\infty }=\frac{\unicode[STIX]{x1D714}k^{2}}{8(1+\overline{\unicode[STIX]{x1D6FA}})}(4+10\overline{\unicode[STIX]{x1D6FA}}+8\overline{\unicode[STIX]{x1D6FA}}^{2}+3\overline{\unicode[STIX]{x1D6FA}}^{3}).\end{eqnarray}$$
Less surprisingly, in all previous analytical studies, the same expression was obtained for the coefficient
$L_{\infty }$
, since the wave-induced mean flow does not affect dispersion in the leading-order approximation.
In order to compare these analytical models and also to clarify the debate on the effects of large (positive) vorticities on the modulational instability, we have computed for each NLS model the maximum growth rates with (4.4) as the vorticity is varied. The results of these computations are shown in figure 14, as well as the strictly numerical results for waves with
$\unicode[STIX]{x1D707}=10$
and
$\unicode[STIX]{x1D700}=0.01$
and the corresponding analytical predictions obtained with (4.10).
Irrespective of the sign of the vorticity, the analytical models yield almost identical results when the velocity shear is sufficiently small. For large vorticities, however, the model of Li et al. (Reference Li, Hui, Donelan, Horikawa and Maruo1987) is the only one that predicts restabilization at large vorticities for downstream propagating waves (
$\unicode[STIX]{x1D6FA}>0$
). As shown in figure 14 by the numerical results, when
$\unicode[STIX]{x1D707}=10$
and
$\unicode[STIX]{x1D700}=0.01$
, this restabilization is not confirmed by our numerical results, which are in turn in close agreement, over the range of
$\unicode[STIX]{x1D6FA}$
values investigated here, with the analytical predictions from Thomas et al. (Reference Thomas, Kharif and Manna2012).
Figure 14. Maximum growth rate of the modulational instability as a function of the dimensionless shear,
$X=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FA}_{0}/\unicode[STIX]{x1D714}$
, for deep water waves: – – –, Baumstein’s model; — ⋅ —, Li et al.’s model; ——, Thomas et al.’s model. For waves with
$\unicode[STIX]{x1D707}=10$
and
$\unicode[STIX]{x1D700}=0.01$
, the symbols correspond to the present numerical results and the dotted line represent the weakly nonlinear results from Thomas et al.’s model.
Baumstein (Reference Baumstein1998) has claimed that ‘…the shear strength increase first enhances, but later suppresses, the instability.’. Surprisingly the analytical results, obtained with (4.10) and neglecting the term associated with the wave–mean-flow interaction (second term in the bracket of
$M_{\infty }$
), contradict this statement. In fact, for the background flow, Baumstein (Reference Baumstein1998) has assumed a piecewise-linear velocity profile, and thus a layer of uniform vorticity (
$\unicode[STIX]{x1D6FA}$
) and depth (
$\unicode[STIX]{x1D6E5}$
) overlies an infinitely deep fluid (at rest). More importantly, his conclusion has been reached from numerical computations with finite values of the shear depth. It should be realized, however, that in the limit of infinite shear depth there is some inconsistency with this approach, owing to his definition of the surface current,
$U_{0}=\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D6E5}$
where
$\unicode[STIX]{x1D6E5}$
is the shear depth. Hence, the condition that
$U_{0}$
remains bounded in the limit
$\unicode[STIX]{x1D6E5}\rightarrow \infty$
implies to consider simultaneously the limit
$\unicode[STIX]{x1D6FA}\rightarrow 0$
. This explains perhaps why Baumstein’s model for deep water waves yields consistent results on a rather short range of (small) vorticities, as shown in figure 14.
In figure 14 it is seen that the other analytical models do not predict, for large vorticities, the restabilization in the case of upstream propagating waves (
$\unicode[STIX]{x1D6FA}<0$
), and their predictions of the maximum growth rate differ remarkably from those obtained with the model of Thomas et al. (Reference Thomas, Kharif and Manna2012), which are in very good agreement with the reference numerical results.
This restabilization of upstream propagating waves is also observed when
$\unicode[STIX]{x1D707}=2$
, as shown in figure 8 for
$\unicode[STIX]{x1D700}=0.01$
. For this case, it is noticed that no other instabilities have been detected in our computations as
$\unicode[STIX]{x1D6FA}$
was varied in the range
$-2\leqslant \unicode[STIX]{x1D6FA}\leqslant \unicode[STIX]{x1D6FA}_{c1}$
. For higher wave steepness, however, several bubbles of instabilities were detected, but not shown here, as
$\unicode[STIX]{x1D6FA}$
was varied in this range. These two-dimensional instabilities are not of the modulational type. Although we have not yet identified all of them, we may rest on our experience to presume that these instabilities are observed as the instability bands of either class I (
$m\geqslant 1$
) or class II (
$m\geqslant 1$
) cross integer values of
$p\geqslant 1$
. The properties of these instabilities are currently under investigation and will be reported in a future publication.
5 Conclusions
With the purpose of investigating the two-dimensional linear stability of finite-amplitude steady waves on a linear shear current, we have developed firstly a Fourier approximation method of these solutions that enables very accurate determination of the basic wave motion. The present work concerns only the forward mode. The results have been compared successfully with other numerical solutions, available in the published literature. An excellent agreement have also been observed between the present numerical results and those obtained with a third-order approximation of the solution, provided the steepness is small and, in the case of upstream propagating waves, the shear is sufficiently small.
The stability analysis of these two-dimensional gravity waves propagating steadily on a vertical shear current of constant vorticity has been carried out in deep water (
$\unicode[STIX]{x1D707}=10$
) and finite depth (
$\unicode[STIX]{x1D707}=2$
). The effects of the vorticity on the boundaries of the instability band and the characteristics of the most unstable mode have been analysed.
For weakly nonlinear waves and moderate values of the vorticity we have rediscovered the analytical results of Thomas et al. (Reference Thomas, Kharif and Manna2012) on long-wavelength instabilities associated with quartet resonances, i.e. class I (
$m=1$
) with
$p<1$
. We have extended their results to higher values of the wave steepness of the basic wave and of the vorticity. In deep water, we found some differences between numerical and analytical growth rates for
$\unicode[STIX]{x1D700}>0.10$
. In finite depth, deviations between the two approaches occur for
$\unicode[STIX]{x1D700}>0.10$
and
$\unicode[STIX]{x1D6FA}>0.5$
. Note that disagreement occurs for smaller values of the vorticity when the wave steepness increases. Furthermore, the numerical results confirm the restabilization of the Benjamin–Feir modulational instability as the shear is increased, both with upstream and downstream propagating waves, in agreement with the analytical results of Thomas et al. (Reference Thomas, Kharif and Manna2012).
We have also identified instabilities associated with quartet resonances of class I (
$m=1$
) with
$p>1$
and higher-order quintet resonances of class II (
$m=1$
) with
$p>1$
. Although these instabilities are weak, increasing the vorticity has substantially different effects depending on the direction of propagation of the finite-amplitude waves along the shear flow. For the case
$\unicode[STIX]{x1D700}=0.05$
,
$\unicode[STIX]{x1D707}=2$
, and the considered range of shear values, the present numerical results reveal that the Benjamin–Feir instability always dominates for downstream propagating waves, whereas the most unstable disturbance of the class I (
$m=1$
) with
$p>1$
has the largest the growth rate for upstream propagating waves, provided the shear parameter
$\unicode[STIX]{x1D6FA}$
is less than a certain negative critical value.
Acknowledgements
We are indebted to two anonymous referees for their careful reading of the manuscript and helpful comments and references. This work was supported by the Direction Generale de l’Armement and funded by the ANR project no. ANR-13-ASTR-0007.
