M Francius

C Kharif

Two-dimensional stability of finite amplitude gravity waves on water of finite depth with constant vorticity

 to fully nonlinear waves. In particular, the restabilization of the Benjamin-Feir modulational instability, whatever the depth, for an opposite shear current is confirmed. For these side-band instabilities, the numerical results show some deviations with the weakly nonlinear theory as the wave steepness of the basic wave and vorticity are increased. Besides the modulational instabilities new instability bands corresponding to quartet and quintet instabilities, which are not side-band disturbances, are discovered. The present numerical results show that with opposite shear currents increasing the shear reduces the growth rate of the most unstable side-band instabilities, but enhances the growth rate of these quartet instabilities, which eventually dominate the Benjamin-Feir modulational instabilities.

Introduction

The huge literature on free surface water waves is mainly based on the assumption of irrotational flows. 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 interaction of waves with vertically sheared currents occurring in nature, one can cite [START_REF] Peregrine | Interaction of water waves and currents[END_REF], [START_REF] Jonsson | Wave-current interactions[END_REF] and [START_REF] Thomas | Wave-current interactions in the near shore region[END_REF].

Periodic gravity waves propagating steadily on a rotational current were studied † Email address for correspondence: marc.francius@mio.osupytheas.fr by many authors, either mathematically or analytically and numerically. Since the pioneering work of [START_REF] Constantin | Exact steady periodic water waves with vorticity[END_REF], 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 the recent rigorous results, the reader can refer to [START_REF] Constantin | Steady periodic water waves with constant vorticity: Regularity and local bifurcation[END_REF] and [START_REF] Kozlov | Dispersion equation for water waves with vorticity and stokes waves on flows with counter-currents[END_REF]. Among the authors using asymptotic methods or purely numerical methods, on can cite [START_REF] Tsao | Behaviour of surface waves on a linearly varying current[END_REF], [START_REF] Dalrymple | A finite amplitude wave on a linear shear current[END_REF], [START_REF] Brevik | Higher-order waves propagating on constant vorticity currents in deep water[END_REF], [START_REF] Simmen | Steady deep-water waves on a linear shear current[END_REF], Teles da Silva & [START_REF] Teles Da Silva | Steep, steady surface waves on water of finite depth with constant vorticity[END_REF] , [START_REF] Kishida | Stokes theory for waves on linear shear current[END_REF], [START_REF] Vanden-Broeck | Periodic waves with constant vorticity in water of infinite depth[END_REF], [START_REF] Swan | A simple analytical model for surface water waves on a depthvarying current[END_REF], [START_REF] Ko | Large-amplitude steady rotational water waves[END_REF], [START_REF] Pak | Free surface waves on shear currents with non-uniform vorticity: Third order solutions[END_REF], [START_REF] Cheng | On the interaction of deep water waves and exponential shear currents[END_REF] , [START_REF] Moreira | Vorticity effects on nonlinear wave-current interactions in deep water[END_REF], [START_REF] Hsu | Gravity-capillary waves in finite depth on flows of constant vorticity[END_REF], [START_REF] Ribeiro-Jr | Flow strurcture beneath rotational water waves with stagnation points[END_REF].

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 side-band 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 as compared to the carrier frequency).

Under the weakly nonlinear assumption, [START_REF] Johnson | On the modulation of water waves on shear flows[END_REF] 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 correspond to the Stokes wave) was proposed. Using a similar approach and extending the analysis to three-dimensional disturbances, [START_REF] Oikawa | The propagation of nonlinear wave packets in a shear flow with a free surface[END_REF] 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 [START_REF] Oikawa | The propagation of nonlinear wave packets in a shear flow with a free surface[END_REF] 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 to discriminate between the effects of a uniform current (equals 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 [START_REF] Miroshnikov | The Boussinesq-Rayleigh approximation for rotational solitary waves on shallow water with uniform vorticity[END_REF]. These approaches are therefore limited to some regions of the parameter space, defined in the present work by the wave steepness ε, 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, [START_REF] Li | Effects of velocity shear on the stability of surface deep water wave trains[END_REF] studied the two-dimensional side-band 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 er-roneous as noted by [START_REF] Baumstein | Modulation of gravity waves with shear in water[END_REF]. In fact, the latter author investigated the effect of piecewise-linear velocity profiles on the side-band instability of a finite-amplitude gravity wave in deep water. Extension of these studies to the case of finite depth was carried out by [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF], 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 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, [START_REF] Hur | Modulational instability in the whitham surface with surface tension and vorticity[END_REF] have recently studied the modulational stability and instability of periodic traveling waves on a linear shear current in finite depth, using the so-called vorticity-modified Whitham equation that combines 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 traveling gravity waves (with sufficiently small amplitudes) are spectrally unstable to long wavelengths 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 that respect, it was shown in [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF] 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 [START_REF] Hur | Modulational instability in the whitham surface with surface tension and vorticity[END_REF] 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 (1978a) and Longuet-Higgins (1978b), [START_REF] Okamura | The linear stability of finite amplitude surface waves on a linear shearing flow[END_REF] found three kinds of instabilities: two of them related to four-and five-wave resonant interactions and a third type called rotational instability, essentially three-dimensional, due to the effects of a strong nonlinear coupling between the mean-flow response and the fundamental wave, as discovered first by [START_REF] Chow | Instabilities of waves on a free surface[END_REF]. 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 fourwave and five-wave resonances. Though, 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 poor convergence property of their numerical method. Like 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 analyze the nonlinear stability of finite amplitude periodic waves, initially weakly modulated due to the presence of side-band disturbances. [START_REF] Nwogu | Interaction of finite-amplitude waves with vertically sheared current fields[END_REF] 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 narrow band sea states. For flows with constant vorticity, [START_REF] Choi | Nonlinear surface waves interacting with a linear shear current[END_REF] 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, [START_REF] Choi | Nonlinear surface waves interacting with a linear shear current[END_REF] 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 [START_REF] Nwogu | Interaction of finite-amplitude waves with vertically sheared current fields[END_REF] and [START_REF] Choi | Nonlinear surface waves interacting with a linear shear current[END_REF] have 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 it simplifies 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 section 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 section 3, we present the numerical methods used in this study and the results of their validation against available numerical and analytical results. In section 4, the results for subharmonic perturbations are presented, first, for the side-band instability (longwavelength disturbances) and, then, for the shorter disturbances involved in quartet and quintet resonances with the basic wave. We summarize our conclusions in section 5.

Mathematical formulation

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 (-Ω 0 ), which varies linearly in the vertical direction and vanishes at the undisturbed free surface or mean water level (z = 0).

For periodic 2D waves of permanent form on a linear shear current it is convenient to consider a co-ordinate 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 stream function ψ(x, z), unique up to an additive constant, such that the total relative fluid velocities are given by

u = U B (z) -c + ũ = ∂ψ ∂z , w = w = - ∂ψ ∂x
where the tilde quantities are perturbations due to the waves. We note that by defining U B (z) = Ω 0 z, the focus is put on the influence of vorticity on the combined wave-current system rather than those of 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, ω = (∂ x w -∂ 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

∂ 2 ψ ∂x 2 + ∂ 2 ψ ∂z 2 = Ω 0 (2.1)
as the governing equation to be satisfied throughout the fluid, namely -h < z < η(x) where h is the mean depth and z = η(x) represents the unknown position of the free surface relatively 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

ψ = 0 on z = -h (2.2) ψ = -Q on z = η (2.3)
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 give the dynamic free surface condition,

gη + 1 2 u 2 + w 2 = R on z = η (2.4)
where 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 ψ(x, z) and η(x) that are periodic in x with finite wave length λ = 2π/k, where k is the wave number. To guarantee a unique solution, when it exists, the wave height H is specified as

η(0) -η(λ/2) = H (2.5)
where η(0) and η(λ/2) are the crest and trough elevation respectively, and we impose the further condition of zero mean level condition,

λ 0 ηdx = 0 (2.6)
Solutions for steep waves on a linear shear current are not currently obtainable analytically, and numerical approaches are necessary to solve the above nonlinear boundaryvalue problem. In this work we use a simple extension of the Fourier approximation method, originated with [START_REF] Rienecker | A Fourier approximation method for steady water waves[END_REF] for irrotational waves, where a Fourier series is assumed only for the stream function as

ψ(x, z) = Ω 0 2 z 2 -h 2 -c (z + h) + α ∞ j=1 B j sinh jk (h + z) cosh jkh cos jkx (2.7)
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 α = 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 -Ω 0 , the coefficients B j , the wave speed c and the unknown position of the free surface η(x) can be determined numerically, after truncation of the Fourier series and a collocation method as described in section 3.

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 the wave speed c, it follows that we can introduce a generalized velocity potential φ(x, z, t), which satisfies the relations

φ x = ψ z -Ω 0 (z + h) , φ z = -ψ x (2.8)
to guarantees that φ is the harmonic conjugate of the function

ψ h = ψ(x, z, t) -1 2 Ω 0 (z + h) 2 .
The unsteady water wave problem can then be reformulated in the moving frame and, in terms of ψ and φ, the governing equations are given by

ψ xx + ψ zz = Ω 0 , -h < z < η
(2.9)

ψ x = 0 on z = -h (2.10) η t + η x ψ z + ψ x = 0 on z = η(x, t) (2.11) φ t + gη -Ω 0 ψ + 1 2 ψ 2 z + ψ 2 x = f (t) + 1 2 c 2 -Ω 0 ch on z = η(x, t) (2.12)
where f (t) is an arbitrary function of time, which can be set to zero without loss of generality or absorbed in the definition of either φ or ψ 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.10) and (2.11) both express that ψ 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

η = η(x) + η(x, t), ψ = ψ(x) + ψ(x, t) (2.13)
where the quantities with an overbar correspond to the periodic wave solution and the quantities with the tilde denote small perturbations, i.e. |η| ≪ |η| and | ψ| ≪ | ψ|. After linearization of the governing equations (2.9-2.12) we obtain

ψxx + ψzz = 0, -h < z < η (2.14) ψx = 0 on z = -h (2.15) ηt + ψz ηx + ηx ψz + ψzz ηx + ψxz η + ψx = 0 on z = η(x) (2.16) φt + g + ψx ψxz + ψz ψzz -Ω 0 ψz η -Ω 0 ψ + ψz ψz + ψx ψx = 0 on z = η(x) (2.17)
By using the Cauchy-Riemann relations (2.8) it is easy to check that without vorticity (Ω 0 =0) the above equations reduce to those derived by [START_REF] Mclean | Instabilities of finite-amplitude gravity waves on water of finite depth[END_REF].

Since the linearized equations (2.14-2.17) 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

η(x, t) = e -iγt e ipx ∞ j=-∞ a j e ijx (2.18) ψ(x, z, t) = e -iγt e ipx ∞ j=-∞ b j e ijx sinh [k j (z + h)] cosh(k j h) (2.19) φ(x, z, t) = e -iγt e ipx ∞ j=-∞ c j e ijx cosh [k j (z + h)] cosh(k j h) (2.20)
where k j = |p + j|, p an arbitrary real number and γ is an unknown, possibly complex, eigenvalue to be determined by the requirement that the linearized surface equations (2.16) and (2.17) have a non trivial solution. Here we have taken the gravitational acceleration g = 1 and the unperturbed wavelength λ = 2π, which is equivalent to non-dimensionalizing the problem with the spatial and temporal scale l ref = λ/(2π) and t ref = √ gk (see section 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 (λ = 2π) and the disturbance is called superharmonic. When p is not an integer, the disturbance contains components with wavelength greater than 2π and it is called a subharmonic perturbation. If γ is real, then the disturbance is said to be stable, whereas if γ 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 φx = ψz and φz = -ψx , we obtain the following relations between the sets of coefficients b j , c j c j = sign (p + j) b j for |j| ∞, (2.21)

Using this relation and substituting the normal mode decomposition of the general solution into (2.16) and (2.17), we obtain the following system of equations

G(x) ∞ j=-∞ a j e ijx + ∞ j=-∞ H j (x)b j e ijx = iγ ∞ j=-∞ F j (x)b j e ijx (2.22) ∞ j=-∞ K j (x)a j e ijx + ∞ j=-∞ L j (x)b j e ijx = iγ ∞ j=-∞ a j e ijx (2.23)
for 0 x 2π, where

G(x) = g + ψx ψxz + ψz ψzz -Ω 0 ψz (2.24) H j (x) = i (p + j) ψx -Ω 0 sinh [k j (h + η)] + k j ψz cosh [k j (h + η)] cosh (k j h) (2.25) F j (x) = - i sign (p + j) cosh [k j (h + η)] cosh (k j h) (2.26) K j (x) = i (p + j) ψz + ψzz ηx + ψxz (2.27) L j (x) = i (p + j) sinh [k j (h + η)] + k j ηx cosh [k j (h + η)] cosh (k j h) (2.28)
The system of equations (2.22)-(2.23) can be interpreted as a generalized eigenvalue problem for the complex eigenvalue γ and 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 relabeled 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 p < 1, which does not mean that our analysis is confined to this range.

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

γ ± j (p) = sign(p + j) -k j c + ω ± j (2.29) ω ± j (p) = - Ω 0 tanh (k j h) 2 ± gk j tanh (k j h) + Ω 0 tanh (k j h) 2 2 (2.30)
where j = 0, ±1, ±2, • • • is a mode number. Since the flat surface is regarded here as the limit of steady waves on a linear shear current with the waveheight approaching zero, the eigenvalues contain the wave speed c = ω/k. The corresponding eigenvectors are such that only those components with suffices p + j are nonzero. All the eigenvalues (2.29) 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 sign designates copropagating and contrapropagating disturbances relative to the undisturbed flow. Like in the irrotational case, it is expected that for small values of the wave steepness instability bands may happen for modes with p values near the points p where two eigenvalues for zero amplitude are equal,

γ ± j1 (p) = γ ± j2 (p) (2.31)
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.

• 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

γ + m (p) = γ - -m (p) m 1 (2.32)
• 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

γ + m (p) = γ - -m-1 (p) m 0 (2.33)
It is noteworthy that there are no collisions with m 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

ω + m + ω - -m = 2mω (2.34) k m + k -m = 2mk (2.35)
and for the Class II as

ω + m + ω - -m-1 = (2m + 1)ω (2.36) k m + k -m-1 = (2m + 1)k (2.37)
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.32) and (2.33) change as illustrated in Figure 1 for m = 1, 2 and kh = 10. For large positive values of Ω, 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.32) and (2.33) 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 ω = ω ± j and wave number κ = p + j, the excess total energy is given, up to second order, by

E = ρgã 2 2κtanhκh ω ω + 1 2 Ω 0 tanhκh (2.38)
where ã denotes the perturbation amplitude. Then, using the invariance of wave action A = E/ω through Galilean transform we can deduce that, in a frame moving uniformly at a velocity c, the excess total energy is

E ′ = ρgã 2 2tanhκh ω + 1 2 Ω 0 tanhκh ω κ -c (2.39)
Equation (2.39) shows that the excess energy is proportional both to the pseudo frequency ω + 1 2 Ω 0 tanhκh and the apparent frequency ω -cκ. 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 (ω/k > c > 0) and for counter-propagating disturbances (ω/k < 0), the excess energy is always positive.

Numerical methods

The numerical calculations consist in two parts, determination of the unperturbed basic flow η, ψ 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.

Computation of steady waves

To calculate the unperturbed wave, the kinematic and dynamic surface conditions (2.3-2.4) 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 η i = η(x i ) where x i = (i -1)dx, i = 1, ..., N + 1 and dx = λ/2N , the above equations yield 2N + 2 equations for the 2N + 4 variables [η 1 , ..., η N +1 , B 1 , ..., B N , c, R, Q] for a given value h of the mean depth. The required two additional equations are (2.5) and (2.6), 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,

l ref = 1 k , t ref = 1 √ gk , c ref = g k
the dimensional variables can be scaled as

η ′ = kη, ψ ′ = k 3 g ψ, c ′ = k g c, R ′ = k g R, Q ′ = k 3 g Q, Ω ′ = Ω 0 √ gk
and the dimensionless equations to be solved become, omitting the primes Formally, the above system of 2N + 4 nonlinear equations can be written as

Ω 2 η 2 i -µ 2 -c (µ + η i ) + N j=1 B j sinh j (µ + η i ) cosh jµ cos jkx i + Q = 0 (3.1) η i + 1 2   N j=1 jkB j sinh j (µ + η i ) cosh jµ sin jkx i   2 + 1 2   Ωη i -c + N j=1 jkB j cosh j (µ + η i ) cosh jµ cos jkx i   2 -R = 0 (3.2) 1 2 η 1 + 1 2 η N +1 + i 2 η i = 0 (3.3) η 1 -η N +1 -2ε = 0 (3.4)
F m (Z) = 0 m = 1, ..., 2N + 4 where Z = [η 1 , ..., η N +1 , B 1 , ..., 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 ε = Hk/2, dimensionless mean depth µ = kh and dimensionless vorticity parameter Ω, 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 [START_REF] Cokelet | Steep gravity waves in water of uniform arbitrary depth[END_REF], for irrotational steady waves (Ω = 0), and the few (rare) results reported by Teles da Silva & [START_REF] Teles Da Silva | Steep, steady surface waves on water of finite depth with constant vorticity[END_REF] for finite amplitude periodic waves on a linear shear current of finite depth. Tables 1 and2 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.

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 [START_REF] Hsu | Gravity-capillary waves in finite depth on flows of constant vorticity[END_REF] 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 ε = 0.20, Ω = 0.5 and µ = 2. It should be realized that although the relative error on c is small, as shown in figure 2, the errors on η and φ 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 stream function.

In general, the performance of most numerical methods for solving exact nonlinear waves deteriorates as the wave gets steeper, and the present method can not escape this fact. Indeed, the present approach fails owing to divergence of the (truncated) series for the stream function, 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 µ, this occurs either with downstream propagating waves (when Ω > 0) or with upstream propagating waves (when Ω < 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 stream function, because the analytic continuation of the stream function 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 behavior 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 [START_REF] Vanden-Broeck | Periodic waves with constant vorticity in water of infinite depth[END_REF][START_REF] Ashton | A non-local formulation of rotational water waves[END_REF] or on the conformal mapping technique [START_REF] Choi | Nonlinear surface waves interacting with a linear shear current[END_REF][START_REF] Ribeiro-Jr | Flow strurcture beneath rotational water waves with stagnation points[END_REF], 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 to compute waves on a linear shear current when an interior stagnation point and a recirculating region exist, either with or without pressure anomalies first mentioned in Teles da Silva [START_REF] Teles Da Silva | Steep, steady surface waves on water of finite depth with constant vorticity[END_REF]. To illustrate this, we show in figure 4 the streamlines pattern for three different periodic waves with Ω = -2 and kh = 1. For the sake of visibility, the spacing between the streamlines is constant in each region separated by the separatrix streamline (Ψ = 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 it 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 center under the crest. As the steepness increases, the center 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-Jr et al. (2017), our results fall in "region 2" of their paper (see their figure 5).

Computation of eigenvalues

Once the unperturbed wave has been computed, the normal modes decompositions (2.18-2.20) are truncated at N Fourier modes to solve the eigenvalue problem. Like [START_REF] Mclean | Instabilities of finite-amplitude gravity waves on water of finite depth[END_REF] 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.22-2.23) into a discrete system of order 4N + 2, which can be written as

Av = γBv where v = (a -N , . . . , a N , b -N , . . . , b N ) t .
Here the matrices A and B are complex functions of the mode p and the unperturbed wave that is specified itself by ε, kh and Ω. For each disturbance with mode number p, the eigenvalues γ 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, . . . , 4N + 2 and j = -N, . . . , N . Instability corresponds to Im (γ 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 p 0.5, owing to the symmetries of the eigenvalue problem and the degeneracy of the normal mode decomposition with 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. 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 [START_REF] Mclean | Instabilities of finite-amplitude gravity waves on water of finite depth[END_REF] and [START_REF] Francius | Three-dimensional instabilities of periodic gravity waves in shallow water[END_REF] for two-dimensional quartet interactions between irrotational waves on finite depth (Ω = 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 Ω = 0, since we have not found in the literature useful quantitative results for the 2D instabilities due to quartet resonances; secondly, the numerical results may confirm some weakly nonlinear predictions like, for instance, the complete restabilization to two-dimensional long-wavelength perturbations, which occurs for sufficiently large (constant) vorticity and kh 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 ε = 0.05, Ω = 0, ±0.2, ±0.4 and two values of the dimensionless depth µ = 10, 2. In either case an increase of positive shear yields a wider range of unstable wavenumbers and enhanced growth-rates. For ε = 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 µ decreases. As is well-known for irrotational waves, the complete restabilization to two-dimensional long-wavelength perturbations occurs for the critical value µ ≈ 1.363. It should also be noticed that for p ≪ 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.

Results

The stability of the periodic waves was studied in details for two values of the dimensionless mean depth, µ = 10 and µ = 2, for a range of values of the vorticity parameter Ω and several values of the wave steepness ε, namely ε = 0.01, 0.05, 0.10, 0.15, 0.20. Actually, for each value of µ and ǫ considered here, our numerical computations of the unperturbed waves cover a large range of shear values in which Ω varies from -2 to the maximal value allowed by the other parameters. In order to facilitate the presentation of our numerical results and to analyze the effects of vorticity on the characteristics of the two-dimensional instabilities under investigation, we have decided to select only few values of ε.

First, we present the numerical results for the quartet instabilities of the modulational type, namely due to side-band wave disturbances (0 < p < 1), and compare them with the results of the approximate analysis of [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF], 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 wavetrain to modulational disturbances are discussed. Emphasis is put, in particular for deep water waves, on the differences between the results from the weakly nonlinear theories and those obtained by numerical calculations.

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 [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF], 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(ξ, τ )

ia τ + La ξξ = M |a| 2 a (4.1)
where ξ = ǫ(x -c g t) and τ = ǫ 2 t are coordinates to describe the slow modulation in space and time of a unidirectional wavetrain with a carrier wave with wavenumber k and frequency ω, c g is the group velocity of the carrier wave. In this framework, ε represents a small dimensionless amplitude parameter and the surface elevation z = η(x, t) can be written, to the leading order, as

η(x, t) = 1 2 ǫa(ξ, τ )exp[i(kx -ωt)] + c.c + O(ǫ 2 ) (4.2)
where 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 governs 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 a(ξ, τ ) = a 0 e -iM a 2 0 τ (4.3) where a 0 is an arbitrary constant, is always modulationally stable to spatial periodic perturbations whose real and imaginary parts are proportional to e iKξ , where ±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 = 2 |M/L|a 0 . Over this range of unstable wavenumbers K, the maximum growth rate of the quartet instability is

γ max = M 1 ωa 2 0 k 2 (4.4)
when

K max = |M 1 /L 1 |a 0 k 2 (4.5)
where L 1 = k 2 L/ω and M 1 = M/(ωk 2 ). To correct a misprint in the definition of L in [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF], the expressions for the dispersive and nonlinear coefficients are given here,

L = ω k 2 σ(2 + X) µ(1 -σ 2 )[1 -µσ + (1 -ρ)X] -σρ 2 M = -ωk 2 (U + V W ) 8(1 + X)(2 + X)σ 4
respectively, where ρ = c g /c p denotes the ratio of the group velocity to the phase velocity of the carrier wave, µ = kh is the dimensionless water depth, σ = tanh(µ) and X = σΩ 0 /ω where Ω 0 represents the dimensional shear. The lengthy expressions for U , V and W can be found in [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF].

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.22)-(2.23), and are shown in figure 6 as a function of the dimensionless vorticity, Ω, for µ = 10 and ε = 0.01, 0.05, 0.10, 0.15, 0.20. For -2 Ω 1 and the two smallest values of ε, a good quantitative agreement between our numerical results and the analytical predictions is observed. In fact, for ε = 0.01, we found some differences between numerical and analytical growth-rates when 2 Ω 4 (not shown here), the latter overestimating the former. For the steepest waves considered here, ε 0.15, the analytical predictions are found to overestimate the characteristics of the most unstable modes, namely p max and γ max , although the variations of these characteristics with Ω are consistent between the two results.

For µ = 2 the same analysis was carried out for ε = 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 Ω (upstream propagating waves), for both p max and γ max , and this even for the steepest waves considered here. However, the two results show a tendency to diverge with increasing (positive) values of Ω, except for the waves with ε = 0.01. In this case, interestingly, a very good agreement between the two results is observed in the range of vorticity -2 Ω 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 (Ω < 0) and large enough negative shear. For small steepness, say ε = 0.01, the numerical results agree very well with those of the weakly nonlinear theory. For downstream propagating waves (Ω > 0) the analytical results of [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF] also predict restabilization of the modulational instability in finite depth, when Ω increases above a positive threshold value that depends only on µ (in the limit p → 0). As shown in figure 7 our numerical results for ε = 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 Ω 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 [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF], 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 → 0 and ε → 0, this theory predicts for µ > 1.363 that the stability boundary is -1 < X X c1 < 0 for upstream propagating waves and 0 < X c2 X for downstream propagating waves (see figure 3 which is related to Ω by the relation

Ω = σ 1/2 Ω 1 + σΩ (4.6)
The critical values X c1 and X c2 depend only on µ and Ω. They can be determined either analytically or graphically, as we have done, by plotting for a fixed µ the variations of the (scaled) nonlinear coefficient M 1 with Ω and finding the critical values for which it vanishes. As discussed in [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF], in the deep water limit σ = 1, X c1 = -2/3 and X c2 increases without bound. We considered waves with µ = 10 and ε = 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 Ω was incremented from -2 to 1 with 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 * = -0.657 (Ω = -1.12), in close agreement with the theoretical value X c ≈ -0.654 (Ω 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.

A similar restabilization occurs with downstream propagating waves (Ω > 0), but as the shear parameter X is increased above the positive threshold value X c2 . Though we have not addressed this issue in detail, namely by comparing the analytical results found by [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF] with the results of the present strictly numerical method, we have compared the marginal stability curve obtained with the analytic solutions of NLS equation with vorticity from [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF] with that obtained by [START_REF] Hur | Modulational instability in the whitham surface with surface tension and vorticity[END_REF] (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 [START_REF] Hur | Modulational instability in the whitham surface with surface tension and vorticity[END_REF] match perfectly the results of [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF] only for downstream propagating waves (Ω > 0) and, say roughly, Ω 1.5. For smaller and decreasing values of Ω, the differences between the two results increases without bounds. More importantly, both analytical results are in total contradiction for upstream propagating waves (Ω < 0). In this case, the model of [START_REF] Hur | Modulational instability in the whitham surface with surface tension and vorticity[END_REF] does not predict the restabilization of the modulational instability whatever the value of the depth, provided the shear is sufficiently large. Though 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 (Ω, µ) parameter space.

Growth rates of two-dimensional instabilities

In our investigations we have also detected other band 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 µ = 2, ak = 0.05 and several values of the shear parameter Ω = 0, ±0.2, ±0.4, ±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 Ω considered here. As we have already seen, the Class I(m = 1) instability region with p < 1 is increased (decreased) as Ω increases (decreases).

Figures 10 and11 show, for upstream propagating waves, the effects of increasing the shear on the instability bands associated to 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, though initially their bandwiths and the maximum growth rates diminish slightly, these increase for sufficiently large negative values of Ω. For Ω = -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 figures 10(c) and 10(d) reveals that this Class I(m = 1) instability region must include p = 2 for some values -0.8 < Ω < -0.4. When this occur, 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 ranges -0.8 < Ω < 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 and13 show, for downstream propagating waves, the effects of increasing the shear on the instability bands associated to Class I(m = 1) and Class II(m = 1) disturbances with p > 1, respectively. The effects of vorticity are different here and seem to be weaker for the waves considered here. As the shear increases, both instability regions moves to the right and, their bandwiths and the maximum growth rates have slightly changed in the range 0 < Ω < 0.8. In contrast with the results obtained for upstream waves, no superharmonic instabilities are detected. It is noticed, however, that the bandwith 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).

Discussion

The results reported in section 4.1 appear to confirm the analytical results of [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF]. These constitute somewhat numerical evidences 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 ε = 0.01 and µ = 10 the restabilization is observed only with upstream propagating waves, while for µ = 2 it occurs both with upstream and downstream propagating waves, in agreement with the analytical results of [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF].

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 Ω < 0 and wave steepness as high as ε = 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 (Ω > 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 µ = 2 and ε = 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 [START_REF] Li | Effects of velocity shear on the stability of surface deep water wave trains[END_REF] using a linear shear current and Baumstein (1998) 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 [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF] in the limit µ → ∞ (σ → 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 [START_REF] Li | Effects of velocity shear on the stability of surface deep water wave trains[END_REF] and [START_REF] Baumstein | Modulation of gravity waves with shear in water[END_REF] applied the same method of multiple scales to the governing equations, like [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF], 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 wavetrains in deep water the feedback from this wave-induced mean flow is due to higher-order effects and do not appear to the order of approximation of NLS equation.

As shown in [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF], however, special care should be taken in the limit µ → ∞ for the NLS equation presented in this work. In this limit they obtained

ia τ + L ∞ a ξξ = M ∞ |a| 2 a (4.7) with L ∞ = - ω k 2 (1 + Ω) 2 (2 + Ω) 3 , M ∞ = ωk 2 8 4 + 6Ω + 6Ω 2 + Ω 3 -Ω 2 (2 + Ω) 2 (1 + Ω) (4.8)
As emphasized in [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF], even though the wave-induced mean velocity approaches zero in the deep water limit, its interaction with the carrier wave yields an O(ǫ 3 ) term that has a non zero finite limit when µ → ∞. 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 ∞ , which cancels only without background vorticity (as expected). When this term is neglected, the expression for M ∞ agrees with that of [START_REF] Baumstein | Modulation of gravity waves with shear in water[END_REF], after adoption of our notation. As noticed by Baumstein (1998), [START_REF] Li | Effects of velocity shear on the stability of surface deep water wave trains[END_REF] obtained a somewhat different formula for the nonlinear coefficient M ∞ which, in our notation, can be expressed

M Li = ωk 2 8(2 + Ω) (8 + 16Ω + 10Ω 2 -Ω 4 ) (4.9)
This should be compared with (4.8), which can also be expressed as

M ∞ = ωk 2 8(1 + Ω) (4 + 10Ω + 8Ω 2 + 3Ω 3 ) (4.10)
Less surprisingly, in all previous analytical studies, the same expression was obtained for the coefficient L ∞ , since the wave-induced mean flow does not affect dispersion to 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 equation (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 µ = 10 and ε = 0.01 and the corresponding analytical predictions obtained with equation (4.8).

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 [START_REF] Li | Effects of velocity shear on the stability of surface deep water wave trains[END_REF] is the only one that predicts restabilization at large vorticities for downstream propagating waves (Ω > 0). As shown in figure 14 by the numerical results, when µ = 10 and ε = 0.01, this restabilization is not confirmed by our numerical results, which are in turn in close agreement, over the range of Ω values investigated here, with the analytical predictions from [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF]. [START_REF] Baumstein | Modulation of gravity waves with shear in water[END_REF] has claimed that "...the shear strength increase first enhances, but later suppresses, the instability.". Surprisingly the analytical results, obtained with (4.8) and neglecting the term associated with the wave-mean-flow interaction (second term in the bracket of M ∞ ), contradict this statement. In fact, for the background flow, [START_REF] Baumstein | Modulation of gravity waves with shear in water[END_REF] has assumed a piecewise linear velocity profile, and thus a layer of uniform vorticity (Ω) and depth (∆) 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 = Ω∆ where ∆ is the shear depth. Hence, the condition that U 0 remains bounded in the limit ∆ → ∞ implies to consider simultaneously the limit Ω → 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 (Ω < 0), and their predictions of the maximum growth rate differ remarkably from those obtained with the model of [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF], which are in very good agreement with the reference numerical results.

This restabilization of upstream propagating waves is also observed when µ = 2, as shown in figure 8 for ε = 0.01. For this case, it is noticed that no other instabilities have been detected in our computations as Ω was varied in the range -2 Ω Ω c1 . For higher wave steepness, however, several bubbles of instabilities were detected, but not shown here, as Ω 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 1) or Class II(m 1) cross integer values of p 1. The properties of these instabilities are currently under investigation and will be reported in a future publication.

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 (µ = 10) and finite depth (µ = 2). The effects of the vorticity on the boundaries of the instability band and the characteristics of the most unstable mode have been analyzed.

For weakly nonlinear waves and moderate values of the vorticity we have rediscovered the analytical results of [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF] 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 ε > 0.10. In finite depth, deviations between the two approaches occur for ε > 0.10 and Ω > 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 [START_REF] Thomas | A nonlinear Schrödinger equation for water waves on finite depth with constant vorticity[END_REF].

We have also identified instabilities associated with quartet resonances of Class I(m = 1) with p > 1 and higher-order quintet resonances 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 ε = 0.05, µ = 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 Ω is less than a certain negative critical value.

Figure 1 .

 1 Figure 1. Resonances loci p as a function of Ω for µ = 10 and different classes: --, Class I(m = 1); ---, Class II(m = 1); -• -, Class I(m = 2); • • • • • • , Class II(m = 2).

Figure 2 .Figure 3 .

 23 Figure 2. Wave velocity c as a function of ε for different values of Ω and µ = 2. Solid lines represent the numerical solution and the dashed lines represent the third-order solution.

Figure 4 .

 4 Figure 4. Streamlines pattern for periodic waves with Ω = -2 and µ = 1: (a) ε = 0.05 for which c = 1.9207810276, (b) ε = 0.25 for which c = 1.9427455675 and (c) ε = 0.45 for which c = 1.998.

Figure 5 .

 5 Figure 5. Growth-rate Im(γ) as a function of p for waves with ε = 0.05, µ = 10 (left) and kh = 2 (right), and different shear values: -, Ω = 0; --, Ω = ±0.2; -.-,Ω = ±0.4. Dotted lines represent the corresponding weakly nonlinear results. The curves above the solid line are for Ω > 0.

Figure 6 .

 6 Figure 6. (a) Wavenumber pmax of the most unstable mode and (b) its growth-rate γmax as a function of Ω, for µ = 10 and ε = 0.01 (stars), 0.05 (asterisks), 0.10 (diamonds), 0.15 (squares), 0.20 (circles); Solid lines represent the weakly nonlinear results.

Figure 7 .

 7 Figure 7. (a) Wavenumber pmax of the most unstable mode and (b) its growth-rate γmax as a function of Ω, for µ = 2 and ε = 0.01 (stars), 0.05 (diamonds), 0.10 (squares), 0.15 (circles); Solid lines represent the weakly nonlinear results.

Figure 8 .

 8 Figure 8. Restabilization of modulational perturbations for wave propagating upstream. Growth rate versus X = σΩ0/ω for different wavenumbers (p ≪ 1) with ε = 0.01 and µ = 10.

Figure 9 .

 9 Figure 9. Growth-rate of quartet instabilities, Class I(m = 1) with p < 1, for µ = 2, ε = 0.05 and for (a) negative and (b) positive values of the shear: --, Ω = 0; ---, |Ω| = 0.2; -• -, |Ω| = 0.4; • • • • • • , |Ω| = 0.8.

Figure 10 .

 10 Figure 10. Growth-rate of quartet instabilities, Class I(m = 1) with p > 1, for µ = 2 and ε = 0.05; (a) for Ω = 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 .

 11 Figure 11. Growth-rate of quintet instabilities, Class II(m = 1) with p > 1, for µ = 2 and ε = 0.05; (a) for Ω = 0; (b) -0.2; (c) -0.4; (d) -0.8; (dashed line) as in figure 10.

Figure 12 .

 12 Figure 12. Growth-rate of quartet instabilities, Class I(m = 1) with p > 1, for µ = 2 and ε = 0.05; (a) for Ω = 0; (b) 0.2; (c) 0.4; (d) 0.8; (dashed line) as in figure 10.

Figure 13 .

 13 Figure 13. Growth-rate of quintet instabilities, Class II(m = 1) with p > 1, for µ = 2 and ε = 0.05; (a) for Ω = 0; (b) 0.2; (c) 0.4; (d) 0.8; (dashed line) as in figure 10.

Figure 14 .

 14 Figure 14. Maximum growth rate of the modulational instability as a function of the dimensionless shear, X = σΩ0/ω, for deep water waves: ---, Baumstein's model; -• -, Li et al.'s model; --, Thomas et al.'s model. For waves with µ = 10 and ε = 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.

Table 1 .

 1 Comparison of results for c 2 of irrotational waves. Results from Cokelet are the most accurate presented by this author.

		Cokelet (1977)	Present results
	ε	µ	c 2	N=16	N=32
	0.0602076 0.6960522	0.615059 0.615059401990042 0.615059401990042
	0.134191 0.7057777	0.666501 0.666501043068160 0.666501043084253
	0.196607 0.7157711	0.727629 0.727628636155880 0.727628677787233
	0.127189 2.3106884	0.997193 0.997192554095718 0.997192554095718
	0.264080 2.3351905	1.05546 1.05545831121995 1.05545831121994
	0.361984 2.3574470	1.12534 1.12533616336992 1.12533616411269
	Ω ε	Teles da Silva	N=16	N=20	N=32
		Peregrine (1988)
	1 0.06 0.5883	0.588263328985539 0.588263329189877 0.588263329193313
	-1 0.25 1.3580	1.35807681735679 1.35807681736914 1.35807681736904
	0.5 1.4421	1.44217541243789 1.44204269298035 -

Table 2 .

 2 Comparison

of c for different finite amplitude periodic waves on a linear shear current with µ = 1.

Acknowledgments

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 n • . ANR-13-ASTR-0007.