Numerical Simulations of Water Waves' Modulational Instability Under the Action of Wind and Dissipation

1. Introduction

The seek of uniform, propagative wave train solutions of the fully nonlinear potential equations has been a major topic for centuries. [5] was the first to propose an expression of such waves, the so called Stokes’ waves. However, pioneer works of [6] emphasized that such waves might be unstable, providing a geometric condition for this stability problem. Later on, [1] showed analytically that Stokes’ waves of moderate amplitude are unstable to long wave perturbations of small amplitude travelling in the same direction. This instability is named the Benjamin-Feir instability (or modulational instability). This result was derived independently by [7] in an averaged Lagrangian approach, and by [8] who used an Hamiltonian formulation of the water wave problem. Using this approach, the latter author derived the nonlinear Schrödinger equation (NLS), and confirmed the previous stability results.

Within the last fifty years, the study of this instability became central for fundamental and applied research. The modulation instability is one of the most important mechanisms for the formation of rogue waves [9]. A complete review on the various phenomena yielding to rogue waves can be found in the book of [10]. In the absence of forcing and damping, Stokes’ waves of specific initial steepness are submitted to this instability, when they encounter perturbations of specific wave numbers [11, 12]. In this case, they encounter a nonlinear quasi-recursive evolution, the so called Fermi-Pasta-Ulam recurrence phenomenon ([13]). This phenomenon corresponds to a series of modulation - demodulation cycles, during which initially uniform wave trains become modulated, leading possibly to the formation of a huge wave. Modulation is due to an energy transfer from the wave carrier to the unstable sidebands. In the wave number space, these unstable sidebands are located in a finite narrow band centered around the carrier wave number. During the demodulation, the energy returns to the fundamental component of the original wave train. Using the Zakharov equation, [14] questions the relevance of the Benjamin-Feir index to indicate the intensity of modulational instability. Indeed, this index is often used to quantify the intensity of interactions between a carrier wave and the finite amplitude sidebands. However, [14] emphasized that nonlinear interactions occur also for sidebands located beyond the Benjamin-Feir instability domain.

A damped nonlinear Schrödinger equation (dNLS) was derived by [15] who revisited the Benjamin-Feir instability in the presence of dissipation. They studied numerically the evolution of narrow bandwidth waves of moderate amplitude. More recently [2] investigated theoretically the modulational instability within the framework of the dNLS equation and demonstrated that any amount of dissipation stabilizes the modulational instability in the sense of Lyapunov. Namely, they showed that the zone of unstable region, in the wavenumber space, shrinks as time increases. As a result, any initially unstable mode of perturbation will finally become stable. [2] have confirmed their theoretical predictions by laboratory experiments for waves of small to moderate amplitude. Later, [3] developed fully nonlinear numerical simulations which agreed with the theory and experiments of [2].

From the latter study we could conclude that dissipation may prevent the development of the Benjamin-Feir instability. This effect questions the occurrence of modulational instability of water wave trains in the field. [16] speculated about the effect of dissipation on the early development of rogue waves and raised the question whether or not the Benjamin-Feir instability was able to spawn a rogue wave.

Nevertheless, these authors did not take the effect of wind into account. When considering the occurrence of modulational instability in the field, the role of wind upon this instability in the presence of dissipation needs to be addressed. Based on this assumption, [4] derived a forced and damped nonlinear Schrödinger equation (fdNLS), and extended the analysis of [2] when wind input is introduced. The influence of wind was introduced through a pressure term acting at the interface, in phase with the wave slope, accordingly to Miles’ theory [17]. This quasi-laminar theory of wind wave amplification is based on the Miles’ shear flow instability. This mechanism of wave amplification is a resonant interaction between water waves and a plane shear flow in air which occurs at the critical height where the wind velocity matches the phase velocity of the surface waves. Stokes waves propagating in the presence of such a forcing, when not submitted to modulational instability, encounter an exponential growth. They demonstrated, within the framework of fdNLS equation, that Stokes’ waves were unstable to modulational instability as soon as the friction velocity is larger than a threshold value. Conversely, for a given friction velocity it was found that only carrier waves presenting frequencies (or wavenumbers) lower than a threshold value are subject to Benjamin-Feir instability. Otherwise, due to dissipation, modulational instability restabilizes in the sense of Lyapunov.

As it was mentioned, this physical result is based on the solution of an approached model, the fdNLS equation. Thus, a proper verification is required. However, the phenomenon at hand is based on the long-time behavior of the modulated wave train when propagating in the presence of wind and dissipation. This remark explains the difficulty to provide an experimental verification of the theory. This physical problem is then especially well adapted for a numerical verification. This verification was performed in a first time by [18], who investigated the development of the modulational instability under wind action and viscous dissipation within the framework of fully nonlinear potential equations. This work is an extension of that of [3] when wind input is considered. Later on, [19] emphasized that the equations empirically introduced by [3] were not completely representative of the dispersion relation in the presence of damping, and corrected the equations to overcome this problem, in accordance with the demonstration of [20] and [21].

This work aims to emphasize how numerical simulations can provide useful information to validate long term results based on weakly nonlinear theory. Furthermore, the numerical approach presented here constitutes an extension of the results of [4] to higher orders of nonlinearity and larger band spectra, too. The long time evolution of modulated wave trains can be investigated in a way not allowed by fdNLS equation. The numerical simulations enable to produce results concerning the long time behavior of the modulated wave train. Especially, the phenomenon of permanent frequency downshift will be investigated.

In section 2, the governing equations of the problem are presented. Section 3 presents the weakly nonlinear model obtained by [4], and summarizes their results. The numerical model used to investigate the long time behavior of the modulated wave train is developed in section 4. The initial conditions used to support the numerical strategy for validating the theory introduced by [4] is presented in section 5. Finally, the results obtained are described in section 6.

2. Governing equations of the problem

The approach used in this study is based on the potential flow theory. The fluid is assumed to be incompressible, inviscid, and animated by an irrotational motion. Thus, the fluid velocity derives from a potential$\varphi$. However, non-potential effects due to wind and viscosity can be taken into account through a modification of the boundary conditions at the surface.

The wind has already been introduced in the dynamic boundary condition through a pressure term acting at the free surface in several numerical potential models. Among them, one may cite [22], [23] and [24] who introduced and discussed this approach for BIEM methods and [25], [26], and [27] who extended it to HOS methods. The pressure term used here is based on the Miles’ theory [17], accordingly to the approach of [4]. The viscosity was introduced heuristically by [3] who used the HOS method to address the question raised in [2] on the restabilisation of the Benjamin-Feir instability of a Stokes wave train in the presence of dissipation. The introduction was made through the addition of a damping term in the dynamic boundary condition. However, a proper derivation of the kinematic and dynamic boundary condition in the presence of viscosity was made by [20], and later on by [21]. A modification of both kinematic and dynamic condition was found, resulting in a slight difference in the dispersion relation, as it was discussed by [19]. Finally, the system of equations corresponding to the potential theory, in the presence of wind and viscous damping reads

$${\varphi }_{xx}+{\varphi }_{zz}=0\mathrm{ }\mathrm{ }\mathrm{ }for\mathrm{ }\mathrm{ }\mathrm{ }-\mathrm{\infty }<z<\eta (x,t)$$E1

$$\nabla \varphi \to 0\mathrm{ }\mathrm{ }\mathrm{ }for\mathrm{ }\mathrm{ }\mathrm{ }z\to -\mathrm{\infty }$$E2

$${\eta }_t+{\varphi }_x{\eta }_x-{\varphi }_z-2\nu {\eta }_{xx}=0\mathrm{ }\mathrm{ }\mathrm{ }for\mathrm{ }\mathrm{ }\mathrm{ }z=\eta (x,t)$$E3

$${\varphi }_t+\frac{1}{2}\left[({\varphi }_x{)}^2+({\varphi }_z{)}^2\right]+g\eta =-\frac{P_a}{\rho }-2\nu {\varphi }_{zz}\mathrm{ }\mathrm{ }\mathrm{ }for\mathrm{ }\mathrm{ }\mathrm{ }z=\eta (x,t),$$E4

where $\varphi (x,z,t)$refers to the velocity potential, $\eta (x,t)$is the free surface elevation, $P_a(x,t)$is the atmospheric pressure due to the wind action, applied at the free surface, and where$g$, $\rho$, and $\nu$are respectively the gravity, the water density and the water kinematic viscosity. In this system of equations, the influence of wind has to be specified. In the absence of wind, the term $P_a/\rho$is equal to zero. Otherwise, the wind action is modeled through the term initially introduced by [17], which reads

$$P_a(x,t)=\frac{{\rho }_{air}\beta u_{\mathrm{*}}^2}{{\kappa }^2}\frac{\partial \eta }{\partial x}(x,t),$$E5

where ${\rho }_{air}$ is the air density, $u_{\mathrm{*}}$the friction velocity, $\kappa$is the von Karman constant, and $\beta$a parameter depending on the friction velocity $u_{\mathrm{*}}$ and the wave carrier velocity$c_0$.

3. Weakly nonlinear approach: The nonlinear Schrödinger equation

The Nonlinear Schrödinger equation can be obtained from the fully nonlinear potential theory by using the multi-scale method. The equations are expanded in Taylor series, around a small parameter, $\epsilon$, the wave steepness. In the presence of forcing and dissipation, this work was performed initially by [15], who obtained a forced and damped version of this equation. The equation obtained is an approximation of the system of equations (1 - 4), correct to the third order in$\epsilon$. Recently, [4] used the forced and damped nonlinear Schrödinger equation (fdNLS),

$$i({\psi }_t+c_g{\psi }_x)-\frac{{\omega }_0}{8k_0^2}{\psi }_{xx}-2{\omega }_0{k}_0^2|\psi {|}^2\psi =i\frac{W{\omega }_0{k}_0}{2g\rho }\psi -2i\nu k_0^2\psi$$E6

to investigate both damping and amplification effects on the Benjamin-Feir instability. Herein, $W={\rho }_{air}\beta u_{\mathrm{*}}^2/{\kappa }^2$represents the wind effect, $c_g={\omega }_0/2k_0$is the group velocity of the carrier wave, and where all the parameters$\nu$, $\rho$, ${\rho }_{air}$, $g$, $u_{\mathrm{*}}$, and $\kappa$are the parameters defined in previous section. Equation (6) describes the spatial and temporal evolution of the envelope of the surface elevation, $\psi$, for weakly nonlinear and dispersive gravity waves on deep water when dissipation, due to viscosity, and amplification, due to wind, are considered. If considering the right hand side of this equation, it can be rewritten as

$$i\left(\frac{W{\omega }_0{k}_0}{2g\rho }-2\nu k_0^2\right)\psi =i\mathcal{K}\psi .$$E7

[4] found that the stability of the envelope depends on the sign of the constant$\mathcal{K}$. For values of$\mathcal{K}<0$, solutions are found to be stable, while for values of$\mathcal{K}\ge 0$, solutions are unstable. Physically, they interpreted this result in terms of frequency of the carrier wave ${\omega }_0$ and friction velocity $u_{\mathrm{*}}$ of the wind over the waves. They plotted the critical curve separating stable envelopes from unstable envelopes. Namely, they showed that for a given friction velocity$u_{\mathrm{*}}$, only carrier wave of frequency ${\omega }_0$ which satisfies the following condition are unstable to modulational perturbations

$$\frac{4\nu {\kappa }^2{\omega }_0}{\beta (\kappa c_0/u_{\mathrm{*}})s{u}_{\mathrm{*}}^2}<1$$E8

This condition can be rewritten as follows

$$\frac{{\mathcal{A}}^2\mathrm{\Omega }}{\beta \left(\mathcal{A}\right)}<1$$E9

where $\mathcal{A}=\kappa c_0/u_{\mathrm{*}}$is associated to wind whereas $\mathrm{\Omega }={\omega }_0/(s{g}^2/4\nu {)}^{1/3}$is associated to dissipation. Note that ${\mathcal{A}}^2\mathrm{\Omega }/\beta \left(\mathcal{A}\right)$ is equal to the ratio between the rate of damping and the rate of amplification and illustrates the competition between dissipative effects and wind input. The non dimensional numbers $\mathcal{A}$and $\mathrm{\Omega }$correspond to the wave age and non dimensional carrier wave frequency. The modulational instability was found to be sustained as soon as the friction velocity is larger than a threshold value. Conversely, for a given friction velocity, it was found that only carrier waves presenting frequencies (or wavenumbers) lower than a threshold value are subject to Benjamin-Feir instability. Otherwise, due to dissipation, modulational instability restabilizes in the sense of Lyapunov.

4. Fully nonlinear approach: The High Order Spectral method

Within the framework of two-dimensional flows, a High-Order Spectral Method is used to solve numerically the basic partial differential equations corresponding to equations (1 - 4). The lateral conditions correspond here to space-periodic conditions. The horizontal bottom condition corresponds to infinite depth. The velocity potential is expanded in a series of eigenfunctions fulfilling both these lateral and bottom conditions. A spectral treatment is well adapted to investigate numerically the long time behavior of periodic water waves encountering the modulational instability.

4.1. Mathematical formulation

We first introduce the following dimensionless variables into equations (1), (2), (3) and (4):

$$\tilde{x}=k_{0}x,\mathrm{ }\mathrm{ }\tilde{z}=k_{0}z,\mathrm{ }\mathrm{ }\tilde{\eta }=k_0\eta ,\mathrm{ }\mathrm{ }\tilde{t}=\sqrt{g{k}_0}t,\mathrm{ }\mathrm{ }\tilde{\varphi }=\varphi /\sqrt{g/k_0^3}\mathrm{ }\mathrm{ }\mathrm{ }and\mathrm{ }\mathrm{ }\mathrm{ }\tilde{p}=p/\left(\rho g/k_0\right),$$E10

where$x$, $z$, $\eta$, $t$, $\varphi$and $p$are dimensional variables, and where $k_0$ is a reference wave number. Hence, the kinematic and dynamic boundary conditions become

$$\frac{\partial \tilde{\eta }}{\partial \tilde{t}}+\frac{\partial \tilde{\varphi }}{\partial \tilde{x}}\frac{\partial \tilde{\eta }}{\partial \tilde{x}}-\frac{\partial \tilde{\varphi }}{\partial \tilde{z}}-2\frac{\nu k_0}{c_0}\frac{{\partial }^2\tilde{\eta }}{\partial {\tilde{x}}^2}=0\mathrm{ }\mathrm{ }\mathrm{ }on\mathrm{ }\mathrm{ }\mathrm{ }\tilde{z}=\tilde{\eta }(\tilde{x},\tilde{t}),$$E11

$$\frac{\partial \tilde{\varphi }}{\partial \tilde{t}}+\frac{\nabla {\tilde{\varphi }}^2}{2}+\tilde{\eta }+{\tilde{P}}_a+2\frac{\nu k_0}{c_0}\frac{{\partial }^2\tilde{\varphi }}{\partial {\tilde{z}}^2}=0\mathrm{ }\mathrm{ }\mathrm{ }on\mathrm{ }\mathrm{ }\mathrm{ }\tilde{z}=\tilde{\eta }(\tilde{x},\tilde{t}).$$E12

Following [8], we introduce the velocity potential at the free surface ${\tilde{\varphi }}_s(\tilde{x},\tilde{t})=\tilde{\varphi }(\tilde{x},\tilde{z}=\tilde{\eta }(\tilde{x},\tilde{t}),\tilde{t})$ into these equations, and it comes, after dropping the tilde for sake of readability,

$$\frac{\partial \eta }{\partial t}=-{\varphi }_{s_x}{\eta }_x+w\left(1+{\eta }_x^2\right)+2\frac{\nu k_0}{c_0}{\eta }_{xx}\mathrm{ }on\mathrm{ }z=\eta (x,t)$$E13

$$\frac{\partial {\varphi }^s}{\partial t}=-\eta -\frac{1}{2}{\varphi }_{s_x}^2+\frac{1}{2}{w}^2\left(1+{\eta }_x^2\right)+2\frac{\nu k_0}{c_0}\frac{{\varphi }_{s_{xx}}-w{\eta }_{xx}}{1+{\eta }_x^2}-P_a\mathrm{ }on\mathrm{ }z=\eta (x,t)$$E14

with

$$w=\frac{\partial \varphi }{\partial z}(x,z=\eta (x,t),t)$$E15

The main difficulty in this approach is the computation of the vertical velocity at the free surface,$w$. Following [28], the potential $\varphi (x,z,t)$is written in a finite perturbation series up to a given order$M$,

$$\varphi (x,z,t)=\sum _{m=1}^M\mathrm{‍}{\varphi }^{\left(m\right)}(x,z,t).$$E16

The term ${\varphi }^{\left(m\right)}$ is of order$\mathcal{O}\left({\epsilon }^m\right)$, where$\epsilon$, a small parameter, is a measure of the wave steepness. Then expanding each ${\varphi }^{\left(m\right)}$ evaluated on $z=\eta$in a Taylor series about$z=0$, we obtain

$${\varphi }_s(x,t)=\sum _{m=1}^M\mathrm{‍}\sum _{l=0}^{M-m}\mathrm{‍}\frac{{\eta }^l}{l!}\frac{{\partial }^{\left(l\right)}}{\partial z^l}\left({\varphi }^{\left(m\right)}(x,z=0,t)\right).$$E17

At a given instant of time, ${\varphi }_s$and $\eta$are known, and we can estimate ${\varphi }^{\left(m\right)}$ at each order${\epsilon }^m$:

$${\varphi }^{\left(1\right)}(x,z=0,t)={\varphi }_s(x,t),m=1,$$E18

$$⋮$$E19

$${\varphi }^{\left(m\right)}(x,z=0,t)=-\sum _{l=1}^{m-1}\mathrm{‍}\frac{{\eta }^l}{l!}\frac{{\partial }^l}{\partial z^l}{\varphi }^{(m-l)}(x,z=0,t),m\ge 2.$$E20

The boundary conditions, together with the Laplace equations ${\nabla }^2{\varphi }^{\left(m\right)}=0$ define a series of Dirichlet problems for${\varphi }^{\left(m\right)}$. For $2\pi -$periodic conditions in$x$, say, ${\varphi }^{\left(m\right)}$can be written as follows in deep water

$${\varphi }^{\left(m\right)}(x,z,t)=\sum _{j=1}^{\mathrm{\infty }}\mathrm{‍}{\varphi }_j^{\left(m\right)}{e}^{-jz}{e}^{ijx}.$$E21

Note that ${\varphi }^{\left(m\right)}(x,z,t)$ satisfies automatically the Laplace equation and the boundary condition $\nabla {\varphi }^{\left(m\right)}\to 0$ when$z\to -\mathrm{\infty }$.

4.2. Computation of the vertical velocity

Substitution of equation (20) into the set of equations (18 - 19) provides an expression of the modes${\varphi }_j^{\left(m\right)}$. The vertical velocity at the free surface is then

$$w(x,t)=\sum _{m=1}^M\mathrm{‍}\sum _{l=0}^{M-m}\mathrm{‍}\frac{{\eta }^l}{l!}\frac{{\partial }^{l+1}}{\partial z^{l+1}}{\varphi }^{\left(m\right)}(x,z=0,t)$$E22

This expression might be substituted into the kinematic and dynamic boundary conditions (13) and (14), yielding to the evolution equations for ${\varphi }_s$ and$\eta$. Another version of HOSM developed by [29] can be used. The difference between both methods lies in the way of computing $w$from${\varphi }_s$. [29] assume a power series for $w$as

$$w(x,t)=\sum _{m=1}^M\mathrm{‍}{w}^{\left(m\right)},$$E23

where

$$w^{\left(m\right)}=\sum _{l=0}^{m-1}\mathrm{‍}\frac{{\eta }^l}{l!}\frac{{\partial }^{l+1}}{\partial z^{l+1}}{\varphi }^{(m-l)}(x,z=0,t).$$E24

In fact, the version of [29] differs from the version of [28] not only in the expression of the approximated vertical velocity at the surface, but also in its subsequent treatment in the free surface equations. According to [29], the surface equations must be truncated at consistent nonlinear order if they are to simulate a conservative Hamiltonian system. This requires to treat carefully all nonlinear terms containing $w$in the prognostic equations. In contrast to the series used by [28], those used by [29] are naturally ordered with respect to the nonlinear parameter. The [28] formulation is not consistent, after truncation, with the underlying Hamiltonian structure of the canonical pair of free-surface equations. Thus, the formulation of [29] preserves the Hamiltonian structure of the prognostic equations.

5. Initial conditions for the numerical simulations

From a numerical point of view, one part of the initial condition is obtained by considering a Stokes wavetrain $(\stackrel{-}{\eta },\stackrel{-}{\varphi })$ which is computed using the approach first introduced by [30]. A very high-order Stokes wave of amplitude $a_0$ and wavenumber $k_0$ is calculated iteratively. In the absence of wind and dissipation, the infinitesimal perturbation components $(\eta \mathrm{\text{'}},\varphi \text{'})$ calculated through a perturbative approach developed by [31] correspond to a Benjamin-Feir instability of wavenumber$\delta k$. The perturbed Stokes wave is obtained by adding the infinitesimal perturbations at the sidebands $k_0\pm \delta k$of the fundamental and its harmonics. For fixed values of $(\mathcal{A},\mathrm{\Omega })$ two kinds of initial conditions are used when wind and dissipation are considered. The first kind (unseeded case) corresponds to the unperturbed Stokes’ wave$(\eta ,\varphi )=(\stackrel{-}{\eta },\stackrel{-}{\varphi })$, whereas the second kind (seeded case) corresponds to the perturbed Stokes’ wave$(\eta ,\varphi )=(\stackrel{-}{\eta },\stackrel{-}{\varphi })+\epsilon (\eta \mathrm{\text{'}},\varphi \text{'}+{\stackrel{-}{\varphi }}_z\eta \mathrm{\text{'}})$, with$\epsilon =10^{-3}$. In both cases, we consider a Stokes wavetrain such as $a_0{k}_0=0.11$ and$k_0=5$. The wavenumber of the modulational instability is$\delta k=1$. This choice of the perturbation wave number corresponds to the closest approximation of the most unstable wave number that can be fitted in the computational domain. The order of nonlinearity was taken equal to$M=6$. In other words, nonlinear terms have been retained up to sixth-order. The highest wavenumber taken into account in the simulations is$k_{max}=50$, corresponding to the ninth harmonic of the fundamental wavenumber. The number of mesh points was taken equal to$N=750$, satisfying the stability criterion$N>(M+1)\times k_{max}$. In the absence of wind and damping, the unperturbed initial condition leads to the steady evolution of the Stokes’ wavetrain, whereas the perturbed initial condition leads to the well known Fermi-Pasta-Ulam recurrence. We propagate these initial wavetrains under various conditions of wind and dissipation, to analyze the behavior of the modulational instability of the Stokes wavetrain.

6. Results and comparisons

One of the difficulties involved in this study is to define clearly the stability. Indeed, since Stokes’ waves are propagating under the action of wind and viscosity, this flow cannot be considered stationary nor periodic. Discussing of the combined influence of wind forcing and damping on the modulational instability, however, implies to define a reference flow. In order to do so, we first consider the evolution of the unperturbed Stokes’ waves in the presence of forcing and dissipation (unseeded case). It is checked that the instability does not develop spontaneously in the laps of time considered. Afterwards, we consider the evolution of the initially perturbed Stokes’ wave train under the same conditions of wind forcing and damping (seeded case). The nonlinear evolution of the Stokes’ wavetrain perturbed by the modulational instability in the presence of wind and dissipation is then compared to that of the reference flow. In that way, the deviation from the reference flow can be interpreted in terms of modulational instability, and the influence of wind forcing and dissipation can be analyzed. Following our previous works [18, 19], the evolution of the energy of the perturbation is thus obtained.

Figures 1 and 2 present the time evolution of the amplitudes of three components of the water waves’ spectrum. The mode $k=5$is the fundamental mode, while modes $k=4$and $k=6$are sidebands, respectively the subharmonic and the superharmonic. Each of these figures present to two kinds of initial conditions, namely the unseeded and the seeded cases. The two figures correspond to two different conditions of wind forcing and damping.

Figure 1.

Time evolution of the normalized amplitudes of the fundamental mode ($k=5$), subharmonic mode ($k=4$) and superharmonic mode ($k=6$) for$(\mathcal{A},\mathrm{\Omega })=\left(\mathrm{4,0.59}\right)$. Fundamental mode amplitude ( —), subharmonic mode amplitude ( —) and superharmonic mode amplitude ( —) for an initially unperturbed Stokes’ wave (unseeded case). Fundamental mode amplitude ( —), subharmonic mode amplitude ( —) and superharmonic mode amplitude ( —) for an initially perturbed Stokes’ wave (seeded case). $T$is the fundamental wave period.

Figure 2.

Time evolution of the normalized amplitudes of the fundamental mode ($k=5$), subharmonic mode ($k=4$) and superharmonic mode ($k=6$) for$(\mathcal{A},\mathrm{\Omega })=\left(\mathrm{4,0.61}\right)$. Fundamental mode amplitude ( —), subharmonic mode amplitude ( —) and superharmonic mode amplitude ( —) for an initially unperturbed Stokes’ wave (unseeded case). Fundamental mode amplitude ( —), subharmonic mode amplitude ( —) and superharmonic mode amplitude ( —) for an initially perturbed Stokes’ wave (seeded case). $T$is the fundamental wave period.

Figure 1 shows the time evolution of the normalized amplitudes $a\left(t\right)/a_0$of the fundamental mode$k=5$, subharmonic mode $k=4$and superharmonic mode $k=6$with and without perturbations for the modulational instability. For both cases, the simulations correspond to a wind parameter $\mathcal{A}=4$and to a viscosity parameter$\mathrm{\Omega }=0.59$. Within the framework of the NLS equation, [4] showed that the wave train is unstable to modulational instability for these values of $\mathcal{A}$and$\mathrm{\Omega }$. From this figure, it appears that both wavetrains (unseeded and seeded cases) present a similar evolution during the first hundred periods of propagation, $T$being the fundamental wave period. Then, the behavior of the wavetrain is strongly affected by the development of the modulational instability. For the unperturbed case (unseeded case), the fundamental component increases, since no occurrence of the modulational instability is expected. However, due to the accumulation of numerical errors, the spontaneous occurrence of the modulational instability cannot be avoided, but not before$t=900T$. On figure 2 one can observe the persistence of the modulational instability through the evolution of the free surface. Indeed, it is observed that the wave packet is alternatively modulated, leading to the formation of a large wave, and demodulated, corresponding to a state which is closer from the origin. This is an expression of the Fermi-Pasta-Ulam quasi-recurrence.

Figure 3.

Surface wave profiles at different times, obtained while propagating initial condition corresponding to seeded case with$(\mathcal{A},\mathrm{\Omega })=\left(\mathrm{4,0.59}\right)$. From top to bottom$t/T=1,\mathrm{ }291,\mathrm{ }436,\mathrm{ }596,\mathrm{ }793,\mathrm{ }846$.

Figure 4.

Surface wave profiles at different times, obtained while propagating initial condition corresponding to seeded case with$(\mathcal{A},\mathrm{\Omega })=\left(\mathrm{4,0.61}\right)$. From top to bottom$t/T=1,\mathrm{ }410,\mathrm{ }496,\mathrm{ }601,\mathrm{ }676,\mathrm{ }885$.

Figure 2 corresponds to $(\mathcal{A},\mathrm{\Omega })$=$\left(\mathrm{4,0.61}\right)$. Wind condition is similar to the previous numerical simulation, but the dissipative effect considered is stronger. This case correspond to a linearly stable case of the modulational instability, as obtained by [4] in the framework of the NLS equation. From this figure, one can see that wind energy goes to the subharmonic mode whereas dissipation reduces the fundamental and superharmonic components, as previously observed. However, modulation of modes decrease, and they present a monotonic behavior. For unseeded case, as expected, we observe an exponential decay of the fundamental mode. Note that there is no natural occurrence of the subharmonic mode of the modulational instability as it was found in figure 1. For seeded case, the first maximum of modulation that occurs at $t=410T$is followed by partial damped modulation/demodulation cycles. Figure 3 illustrates the disappearance of the modulational instability through the evolution of the free surface. In this case dissipation prevails over amplification due to wind and [2] have obtained linear and nonlinear stability of modulational perturbations within the framework of the dissipative NLS equation. More specifically they showed that dissipation reduces the set of unstable wavenumbers as time increases. Consequently every mode becomes stable. The result of this numerical simulation agrees with that of [2] and [3] who considered only dissipation. In their approach, a solution is said to be stable if every solution that starts close to this solution at $t=0$remains close to it for all$t>0$, otherwise the solution is unstable. To include nonlinear stability analysis they introduced a norm and considered stability in the sense of Lyapunov.

In our previous work [18], we assumed that the dominant mode describes the main behavior of a wave train, and we introduced a norm measuring the distance between the fundamental modes of the unperturbed and perturbed Stokes wave corresponding to unseeded case and seeded case respectively. However, it is more consistent to consider the energy of the perturbation, as it was stated in [19]. Thus, another norm can be introduced as

$$E_N\left(t\right)=\frac{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathrm{‍}\left(a_{k_S}\right(t)-a_{k_{US}}(t){)}^{2}dk}{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathrm{‍}{a}_{k_{US}}^2\left(0\right)dk},$$E25

where $a_{k_{US}}\left(t\right)$ is the amplitude of the component of water elevation $\eta$of wave number$k$, for the initially unperturbed wave train (unseeded case), and $a_{k_S}\left(t\right)$ is the amplitude of the component of water elevation $\eta$of wave number$k$, for the initially perturbed wave train (seeded case). This norm corresponds to the potential energy of the perturbation. Its value characterizes the deviation of the perturbed solution from the unperturbed solution.

Figure 5.

Time evolution of the norm $E_n$ for $(\mathcal{A},\mathrm{\Omega })=\left(\mathrm{4,0.59}\right)$ ( —) and $(\mathcal{A},\mathrm{\Omega })=\left(\mathrm{4,0.61}\right)$ ( —).

Figure 6.

Theoretical ( —) and numerical ( –$\cdot$–) marginal stability contour lines. ($\circ$) correspond to numerical results obtained in the framework of the equations suggested by [3]. The theoretical curve corresponds to the figure 1 of [4].

Figure 4 shows the time evolution of this norm for two sets of parameters $(\mathcal{A},\mathrm{\Omega })=\left(\mathrm{4,0.59}\right)$ and$(\mathcal{A},\mathrm{\Omega })=\left(\mathrm{4,0.61}\right)$. For the two cases we can observe two regimes. The first regime corresponds to the development of the modulational instability and shows that it is the nonlinear interaction between the fundamental mode and its sidebands which dominates with a weak effect of the wind forcing and the dissipation. The second regime corresponding to the oscillatory evolution of the norm is dominated by the competition between wind forcing and dissipation. The nonlinear interaction between the fundamental mode and the sidebands is not the dominant mechanism. The magenta curve exhibits oscillations around an averaged value growing exponentially, whereas the yellow curve exhibits the same oscillations around a constant value. We can claim that the norm, $E_N$, presents globally exponential growth or asymptotical saturation corresponding to instability and stability respectively. Herein, the stability can be interpreted in terms of asymptotic stability. The first case is said to be unstable whereas the second case corresponds to a stable solution. In the latter case we expect that the solution will remain close to the unperturbed solution. In other words, nonlinear interactions are affected by the non conservative effects that are wind and dissipation, leading to a long time disappearance of these interactions.

Many numerical simulations have been run for various values of the parameters $\mathcal{A}$and$\mathrm{\Omega }$. Figure 6 shows a stability diagram which presents comparison between the present numerical results and those of [4] obtained theoretically. The marginal curve corresponding to the fully nonlinear equations is very close to the theoretical marginal curve obtained within the framework the NLS equation. The region above the critical curve corresponds to stable cases, whereas the region beneath corresponds to unstable cases. Bars in figure 6 correspond to uncertainty on stability or instability. Numerical results obtained in our previous work [18] within the framework of equations suggested by [3] are plotted for the sake of reference ($\circ$). The way of introducing damping effect into the kinematic boundary condition has little influence on the results, especially for young waves. The present numerical simulations demonstrate that the results derived by [4] within the framework of the NLS equation are correct in the context of the fully nonlinear equations.

This result provides a validation of the weakly nonlinear theory obtained in the framework of nonlinear Schrödinger equation. However, the numerical approach allows to investigate the long time evolution of the wave train, taking into account the strongly nonlinear behavior of water waves. One phenomenon especially illustrates this nonlinear behavior: the permanent frequency downshift. This phenomenon was discussed by [32] and [33] within the framework of gravity waves. These authors considered that dissipation due to breaking wave was responsible for this permanent downshift. [34] modeled the phenomenon in the presence of wind and eddy viscosity, and latter on [35] in the presence of only molecular viscosity. All these works are based on equations valid up to fourth order in nonlinearity, or higher. Indeed, it is well known that the frequency downshift cannot be observed in the framework of nonlinear Schrödinger equation, which preserves the symmetry between subharmonic and superharmonic components. In fact, [36] concluded that in the absence of wind and dissipation, it was not possible to observe the phenomenon even with higher order equations.

If going back to figure 1, for the initially perturbed case (seeded case), the development of the modulational instability is responsible for the frequency downshift observed at around$t=500T$. Indeed, one can see that the subharmonic component increases continuously whereas the fundamental and superharmonic component decrease. The superharmonic component decreases faster than the fundamental component. Hence, wind energy goes to the subharmonic mode whereas dissipation reduces the fundamental and superharmonic components. During the modulation process, a broadening of the spectrum is observed, even if not presented here for the sake of clarity. Beyond$t=500T$, the subharmonic mode $k=4$is dominant in the spectrum.

Figure 7.

Nondimensional time ${\tau }_{downshift}=t/T$for which the permanent frequency downshift is observed, plotted as a function of the wave age parameter$\mathcal{A}$. $\mathrm{\Omega }$is chosen here to be on the marginal stability curve presented on figure 6.

To investigate the effect of wind and damping, another series of simulations is performed. Namely, the values of $\mathcal{A}$and $\mathrm{\Omega }$where chosen to fulfill the condition $\mathcal{K}=0$(see equation (7) for the definition of$\mathcal{K}$). This might be interpreted in terms of balance between wind forcing and damping, corresponding to an equilibrium state. Furthermore, in these conditions, the forced and damped nonlinear Schrödinger equation reduces to the canonical NLS equation. It becomes obvious that the weakly nonlinear theory does not predict any variation in behavior, since the NLS equation remains unchanged for each values of $\mathcal{A}$and$\mathrm{\Omega }$. However, the results of the numerical simulations are different. Indeed, figure 7 presents the nondimensional time ${\tau }_{downshift}$ after which the permanent frequency downshift is observed, as a function of$\mathcal{A}$. From this figure, it seems obvious that this time is strongly dependant on the wind speed. The youngest the waves are, the fastest is the downshift.

7. Conclusion

In this work, it was evidence how numerical simulations can provide a good demonstration of a weakly nonlinear theory that cannot be achieved by means of experimental demonstration. In this study, an extension of the work of [4] to the fully nonlinear case was suggested. Within the framework of the NLS equation the latter authors considered the modulational instability of Stokes wave trains suffering both effects of wind and dissipation. The results they obtained show that the modulational instability depends on both frequency of the carrier wave and strength of the wind velocity. They plotted the curve corresponding to marginal stability in the $(\mathcal{A},\mathrm{\Omega })$-plane. Here, a numerical verification is performed, by means of a fully nonlinear approach. The long term behavior of a wave group propagating under both actions of wind and dissipation was obtained thanks to this method. To distinguish stable solutions from unstable solutions, a norm based on the potential energy of the perturbations was introduced. A nonlinear stability diagram resulting from the numerical simulations of the fully nonlinear equation has been given in the $(\mathcal{A},\mathrm{\Omega })$-plane which coincides with the linear stability analysis of [4]. In the presence of wind, dissipation and modulational instability it is found that wind energy goes to the subharmonic sideband whereas dissipation lowers the amplitude of the fundamental mode of the wave train yielding to a permanent frequency-downshifting. This permanent frequency downshift is strongly influenced by the wind and dissipation parameter. If the wave group is at equilibrium in energy input and dissipation, the fdNLS equation reduces to the classical NLS equation, and predict no influence of the wind. However, by considering the asymmetry between wave components, induced by strong nonlinearity (higher than fourth order), a strong influence of the wind and dissipation is observed.