**2. Governing equations of the problem**

2 Will-be-set-by-IN-TECH

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

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

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

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

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

occur also for sidebands located beyond the Benjamin-Feir instability domain.

numerical simulations which agreed with the theory and experiments of [2].

instability was able to spawn a rogue wave.

restabilizes in the sense of Lyapunov.

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 *φ*. 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

$$
\phi\_{\rm xx} + \phi\_{\rm zz} = 0 \quad \text{for} \quad -\infty < z < \eta(\mathbf{x}, t) \tag{1}
$$

$$
\nabla \phi \to 0 \quad \text{for} \ z \to -\infty \tag{2}
$$

4 Will-be-set-by-IN-TECH 378 Numerical Simulation – From Theory to Industry Numerical Simulations of Water Waves' Modulational Instability Under the Action of Wind and Dissipation <sup>5</sup>

$$
\eta\_t + \phi\_\mathbf{x} \eta\_\mathbf{x} - \phi\_z - 2\nu \eta\_{\mathbf{x}\mathbf{x}} = 0 \quad \text{for} \quad z = \eta(\mathbf{x}, t) \tag{3}
$$

modulational perturbations

This condition can be rewritten as follows

encountering the modulational instability.

*∂η*˜ *<sup>∂</sup>*˜*<sup>t</sup>* <sup>+</sup>

*∂φ*˜ *<sup>∂</sup>*˜*<sup>t</sup>* <sup>+</sup>

**4.1. Mathematical formulation**

modulational instability restabilizes in the sense of Lyapunov.

*<sup>x</sup>*˜ = *<sup>k</sup>*0*x*, *<sup>z</sup>*˜ = *<sup>k</sup>*0*z*, *<sup>η</sup>*˜ = *<sup>k</sup>*0*η*, ˜*<sup>t</sup>* = *gk*0*t*, *<sup>φ</sup>*˜ = *<sup>φ</sup>*/

Hence, the kinematic and dynamic boundary conditions become

<sup>2</sup> <sup>+</sup> *<sup>η</sup>*˜ <sup>+</sup> *<sup>P</sup>*˜

Following [8], we introduce the velocity potential at the free surface *φ*˜

*∂φ*˜ *∂x*˜ *∂η*˜ *<sup>∂</sup>x*˜ <sup>−</sup> *∂φ*˜

<sup>∇</sup>*φ*˜2

**4. Fully nonlinear approach: The High Order Spectral method**

4*νκ*2*ω*<sup>0</sup> *<sup>β</sup>*(*κc*0/*u*∗)*su*<sup>2</sup>

<sup>A</sup>2<sup>Ω</sup>

where <sup>A</sup> <sup>=</sup> *<sup>κ</sup>c*0/*u*<sup>∗</sup> is associated to wind whereas <sup>Ω</sup> <sup>=</sup> *<sup>ω</sup>*0/(*sg*2/4*ν*)1/3 is associated to dissipation. Note that <sup>A</sup>2Ω/*β*(A) 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 A and Ω 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,

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

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

where *x*, *z*, *η*, *t*, *φ* and *p* are dimensional variables, and where *k*<sup>0</sup> is a reference wave number.

*νk*<sup>0</sup> *c*0

*η*˜(*x*˜, ˜*t*), ˜*t*) into these equations, and it comes, after dropping the tilde for sake of readability,

*∂*2*η*˜

*∂*2*φ*˜

*<sup>∂</sup>z*˜ <sup>−</sup> <sup>2</sup>

*<sup>a</sup>* + 2 *νk*<sup>0</sup> *c*0

 *g*/*k*<sup>3</sup>

<sup>0</sup> and *p*˜ = *p*/ (*ρg*/*k*0), (10)

*<sup>s</sup>*(*x*˜, ˜*t*) = *φ*˜(*x*˜, *z*˜ =

*<sup>∂</sup>x*˜<sup>2</sup> <sup>=</sup> 0 on *<sup>z</sup>*˜ <sup>=</sup> *<sup>η</sup>*˜(*x*˜, ˜*t*), (11)

*<sup>∂</sup>z*˜<sup>2</sup> <sup>=</sup> 0 on *<sup>z</sup>*˜ <sup>=</sup> *<sup>η</sup>*˜(*x*˜, ˜*t*). (12)

∗

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

*<* 1 (8)

*<sup>β</sup>*(A) *<sup>&</sup>lt;* <sup>1</sup> (9)

$$\phi\_t + \frac{1}{2}\left[ (\phi\_x)^2 + (\phi\_z)^2 \right] + g\eta = -\frac{P\_d}{\rho} - 2\nu \phi\_{zz} \text{ for } z = \eta(x, t), \tag{4}$$

where *φ*(*x*, *z*, *t*) refers to the velocity potential, *η*(*x*, *t*) is the free surface elevation, *Pa*(*x*, *t*) is the atmospheric pressure due to the wind action, applied at the free surface, and where *g*, *ρ*, and *ν* 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 *Pa*/*ρ* is equal to zero. Otherwise, the wind action is modeled through the term initially introduced by [17], which reads

$$P\_{\mathbf{d}}(\mathbf{x},t) = \frac{\rho\_{air}\beta u\_{\ast}^{2}}{\kappa^{2}} \frac{\partial \eta}{\partial \mathbf{x}}(\mathbf{x},t),\tag{5}$$

where *ρair* is the air density, *u*<sup>∗</sup> the friction velocity, *κ* is the von Karman constant, and *β* a parameter depending on the friction velocity *u*<sup>∗</sup> 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, *ε*, 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 *ε*. Recently, [4] used the forced and damped nonlinear Schrödinger equation (fdNLS),

$$\dot{a}(\psi\_t + c\_{\mathcal{S}}\psi\_x) - \frac{\omega\_0}{8k\_0^2}\psi\_{xx} - 2\omega\_0 k\_0^2|\psi|^2\psi = i\frac{\mathcal{W}\omega\_0 k\_0}{2\text{g}\rho}\psi - 2i\nu k\_0^2\psi\tag{6}$$

to investigate both damping and amplification effects on the Benjamin-Feir instability. Herein, <sup>W</sup> <sup>=</sup> *<sup>ρ</sup>airβu*<sup>2</sup> <sup>∗</sup>/*κ*<sup>2</sup> represents the wind effect, *cg* <sup>=</sup> *<sup>ω</sup>*0/2*k*<sup>0</sup> is the group velocity of the carrier wave, and where all the parameters *ν*, *ρ*, *ρair*, *g*, *u*∗, and *κ* are the parameters defined in previous section. Equation (6) describes the spatial and temporal evolution of the envelope of the surface elevation, *ψ*, 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{\mathcal{W}\omega\_0 k\_0}{2\mathcal{g}\rho} - 2\nu k\_0^2\right)\psi = i\mathcal{K}\psi. \tag{7}$$

[4] found that the stability of the envelope depends on the sign of the constant K. For values of K *<* 0, solutions are found to be stable, while for values of K ≥ 0, solutions are unstable. Physically, they interpreted this result in terms of frequency of the carrier wave *ω*<sup>0</sup> and friction velocity *u*<sup>∗</sup> 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*∗, only carrier wave of frequency *ω*<sup>0</sup> which satisfies the following condition are unstable to 378 Numerical Simulation – From Theory to Industry Numerical Simulations of Water Waves' Modulational Instability Under the Action of Wind and Dissipation <sup>5</sup> Numerical Simulations of Water Waves' Modulational Instability Under the Action of Wind and Dissipation 379

modulational perturbations

4 Will-be-set-by-IN-TECH

<sup>+</sup> *<sup>g</sup><sup>η</sup>* <sup>=</sup> <sup>−</sup> *Pa*

where *φ*(*x*, *z*, *t*) refers to the velocity potential, *η*(*x*, *t*) is the free surface elevation, *Pa*(*x*, *t*) is the atmospheric pressure due to the wind action, applied at the free surface, and where *g*, *ρ*, and *ν* 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 *Pa*/*ρ* is equal to zero. Otherwise, the wind action is modeled through the term initially

> ∗ *κ*2

where *ρair* is the air density, *u*<sup>∗</sup> the friction velocity, *κ* is the von Karman constant, and *β* a

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, *ε*, 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 *ε*. Recently, [4] used the forced and damped nonlinear Schrödinger equation (fdNLS),

<sup>0</sup>|*ψ*|

to investigate both damping and amplification effects on the Benjamin-Feir instability. Herein,

wave, and where all the parameters *ν*, *ρ*, *ρair*, *g*, *u*∗, and *κ* are the parameters defined in previous section. Equation (6) describes the spatial and temporal evolution of the envelope of the surface elevation, *ψ*, for weakly nonlinear and dispersive gravity waves on deep water when dissipation, due to viscosity, and amplification, due to wind, are considered. If

<sup>2</sup>*ψ* = *i*

<sup>∗</sup>/*κ*<sup>2</sup> represents the wind effect, *cg* <sup>=</sup> *<sup>ω</sup>*0/2*k*<sup>0</sup> is the group velocity of the carrier

W*ω*0*k*<sup>0</sup>

<sup>2</sup>*g<sup>ρ</sup> <sup>ψ</sup>* <sup>−</sup> <sup>2</sup>*iνk*<sup>2</sup>

*ψ* = *i*K*ψ*. (7)

<sup>0</sup>*ψ* (6)

*∂η ∂x*

*Pa*(*x*, *<sup>t</sup>*) = *<sup>ρ</sup>airβu*<sup>2</sup>

parameter depending on the friction velocity *u*<sup>∗</sup> and the wave carrier velocity *c*0.

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

*<sup>ψ</sup>xx* <sup>−</sup> <sup>2</sup>*ω*0*k*<sup>2</sup>

*φ<sup>t</sup>* + 1 2 

introduced by [17], which reads

<sup>W</sup> <sup>=</sup> *<sup>ρ</sup>airβu*<sup>2</sup>

(*φx*)<sup>2</sup> + (*φz*)<sup>2</sup>

*<sup>i</sup>*(*ψ<sup>t</sup>* <sup>+</sup> *cgψx*) <sup>−</sup> *<sup>ω</sup>*<sup>0</sup>

8*k*<sup>2</sup> 0

considering the right hand side of this equation, it can be rewritten as

<sup>W</sup>*ω*0*k*<sup>0</sup>

<sup>2</sup>*g<sup>ρ</sup>* <sup>−</sup> <sup>2</sup>*νk*<sup>2</sup>

0 

[4] found that the stability of the envelope depends on the sign of the constant K. For values of K *<* 0, solutions are found to be stable, while for values of K ≥ 0, solutions are unstable. Physically, they interpreted this result in terms of frequency of the carrier wave *ω*<sup>0</sup> and friction velocity *u*<sup>∗</sup> 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*∗, only carrier wave of frequency *ω*<sup>0</sup> which satisfies the following condition are unstable to

*i*

*η<sup>t</sup>* + *φxη<sup>x</sup>* − *φ<sup>z</sup>* − 2*νηxx* = 0 for *z* = *η*(*x*, *t*) (3)

*<sup>ρ</sup>* <sup>−</sup> <sup>2</sup>*νφzz* for *<sup>z</sup>* <sup>=</sup> *<sup>η</sup>*(*x*, *<sup>t</sup>*), (4)

(*x*, *t*), (5)

$$\frac{4\nu\kappa^2\omega\_0}{\beta(\kappa c\_0/u\_\*)su\_\*^2} < 1\tag{8}$$

This condition can be rewritten as follows

$$\frac{\mathcal{A}^2 \Omega}{\beta(\mathcal{A})} < 1 \tag{9}$$

where <sup>A</sup> <sup>=</sup> *<sup>κ</sup>c*0/*u*<sup>∗</sup> is associated to wind whereas <sup>Ω</sup> <sup>=</sup> *<sup>ω</sup>*0/(*sg*2/4*ν*)1/3 is associated to dissipation. Note that <sup>A</sup>2Ω/*β*(A) 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 A and Ω 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):

$$\mathfrak{x} = k\_0 \mathfrak{x}, \ \mathfrak{z} = k\_0 \mathfrak{z}, \ \mathfrak{y} = k\_0 \mathfrak{y}, \ \mathfrak{f} = \sqrt{\mathfrak{g} k\_0} t, \ \breve{\mathfrak{g}} = \mathfrak{g} / \sqrt{\mathfrak{g}/k\_0^3} \text{ and } \ \mathfrak{p} = \mathfrak{p} / (\mathfrak{g} \mathfrak{g}/k\_0), \tag{10}$$

where *x*, *z*, *η*, *t*, *φ* and *p* are dimensional variables, and where *k*<sup>0</sup> is a reference wave number. Hence, the kinematic and dynamic boundary conditions become

$$2\frac{\partial\vec{\eta}}{\partial\vec{t}} + \frac{\partial\vec{\phi}}{\partial\vec{x}}\frac{\partial\vec{\eta}}{\partial\vec{x}} - \frac{\partial\vec{\phi}}{\partial\vec{z}} - 2\frac{\nu k\_0}{c\_0}\frac{\partial^2\vec{\eta}}{\partial\vec{x}^2} = 0 \quad \text{on} \quad \vec{z} = \vec{\eta}(\vec{x}, \vec{t}), \tag{11}$$

$$\frac{\partial \tilde{\phi}}{\partial \tilde{t}} + \frac{\nabla \tilde{\phi}^2}{2} + \tilde{\eta} + \tilde{P}\_a + 2 \frac{\upsilon k\_0}{c\_0} \frac{\partial^2 \tilde{\phi}}{\partial \tilde{z}^2} = 0 \quad \text{on} \quad \tilde{z} = \tilde{\eta}(\tilde{x}, \tilde{t}). \tag{12}$$

Following [8], we introduce the velocity potential at the free surface *φ*˜ *<sup>s</sup>*(*x*˜, ˜*t*) = *φ*˜(*x*˜, *z*˜ = *η*˜(*x*˜, ˜*t*), ˜*t*) into these equations, and it comes, after dropping the tilde for sake of readability,

$$\frac{\partial \eta}{\partial t} = -\phi\_{\rm s\_x} \eta\_x + w \left(1 + \eta\_x^2\right) + 2 \frac{\nu k\_0}{c\_0} \eta\_{xx} \text{ on } z = \eta(x, t) \tag{13}$$

This expression might be substituted into the kinematic and dynamic boundary conditions (13) and (14), yielding to the evolution equations for *φs* and *η*. Another version of HOSM developed by [29] can be used. The difference between both methods lies in the way of

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

*M* ∑ *m*=1

*w*(*m*)

, (22)

(*x*, *z* = 0, *t*). (23)

, *φ*�

), with *ε* = 10−3. In both cases, we consider a Stokes

) calculated

*w*(*x*, *t*) =

*ηl l*!

*∂l*+<sup>1</sup> *<sup>∂</sup>zl*+<sup>1</sup> *<sup>φ</sup>*(*m*−*l*)

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

From a numerical point of view, one part of the initial condition is obtained by considering a Stokes wavetrain (*η*¯, *φ*¯) which is computed using the approach first introduced by [30]. A very high-order Stokes wave of amplitude *a*<sup>0</sup> and wavenumber *k*<sup>0</sup> is calculated iteratively. In the

through a perturbative approach developed by [31] correspond to a Benjamin-Feir instability of wavenumber *δk*. The perturbed Stokes wave is obtained by adding the infinitesimal perturbations at the sidebands *k*<sup>0</sup> ± *δk* of the fundamental and its harmonics. For fixed values of (A, Ω) two kinds of initial conditions are used when wind and dissipation are considered. The first kind (unseeded case) corresponds to the unperturbed Stokes' wave (*η*, *φ*)=(*η*¯, *φ*¯), whereas the second kind (seeded case) corresponds to the perturbed Stokes'

wavetrain such as *a*0*k*<sup>0</sup> = 0.11 and *k*<sup>0</sup> = 5. The wavenumber of the modulational instability is *δ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 *kmax* = 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) × *kmax*. 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

*m*−1 ∑ *l*=0

of [29] preserves the Hamiltonian structure of the prognostic equations.

absence of wind and dissipation, the infinitesimal perturbation components (*η*�

**5. Initial conditions for the numerical simulations**

, *φ*� + *φ*¯

*zη*�

wave (*η*, *φ*)=(*η*¯, *φ*¯) + *ε*(*η*�

computing *w* from *φs*. [29] assume a power series for *w* as

*w*(*m*) =

where

$$\frac{\partial \phi^s}{\partial t} = -\eta - \frac{1}{2}\phi\_{s\_\mathbf{x}}^2 + \frac{1}{2}w^2\left(1 + \eta\_\mathbf{x}^2\right) + 2\frac{\nu k\_0}{c\_0}\frac{\phi\_{s\_{\mathbf{x}\mathbf{x}}} - w\eta\_{\mathbf{x}\mathbf{x}}}{1 + \eta\_\mathbf{x}^2} - P\_\mathbf{d} \text{ on } \mathbf{z} = \eta(\mathbf{x}, t) \tag{14}$$

with

$$w = \frac{\partial \phi}{\partial z}(\mathbf{x}, z = \eta(\mathbf{x}, t), t) \tag{15}$$

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

$$\phi(\mathbf{x}, z, t) = \sum\_{m=1}^{M} \phi^{(m)}(\mathbf{x}, z, t). \tag{16}$$

The term *<sup>φ</sup>*(*m*) is of order <sup>O</sup>(*εm*), where *<sup>ε</sup>*, a small parameter, is a measure of the wave steepness. Then expanding each *φ*(*m*) evaluated on *z* = *η* in a Taylor series about *z* = 0, we obtain

$$\phi\_{\mathbf{s}}(\mathbf{x},t) = \sum\_{m=1}^{M} \sum\_{l=0}^{M-m} \frac{\eta^{l}}{l!} \frac{\partial^{(l)}}{\partial \mathbf{z}^{l}} \left(\phi^{(m)}(\mathbf{x}, \mathbf{z} = \mathbf{0}, t)\right). \tag{17}$$

At a given instant of time, *φ<sup>s</sup>* and *η* are known, and we can estimate *φ*(*m*) at each order *εm*:

$$\phi^{(1)}(\mathbf{x}, z=0, t) = \phi\_{\mathbf{s}}(\mathbf{x}, t), \quad m = 1,\tag{18}$$

$$\mathbf{1}$$

. .

$$\phi^{(m)}(\mathbf{x}, z=0, t) = -\sum\_{l=1}^{m-1} \frac{\eta^l}{l!} \frac{\partial^l}{\partial z^l} \phi^{(m-l)}(\mathbf{x}, z=0, t), \quad m \ge 2. \tag{19}$$

The boundary conditions, together with the Laplace equations <sup>∇</sup>2*φ*(*m*) <sup>=</sup> 0 define a series of Dirichlet problems for *φ*(*m*) . For 2*π*−periodic conditions in *<sup>x</sup>*, say, *<sup>φ</sup>*(*m*) can be written as follows in deep water

$$\phi^{(m)}(\mathbf{x}, z, t) = \sum\_{j=1}^{\infty} \phi\_j^{(m)} e^{-jz} e^{ij\mathbf{x}}.\tag{20}$$

Note that *φ*(*m*)(*x*, *z*, *t*) satisfies automatically the Laplace equation and the boundary condition <sup>∇</sup>*φ*(*m*) <sup>→</sup> 0 when *<sup>z</sup>* → −∞.

#### **4.2. Computation of the vertical velocity**

Substitution of equation (20) into the set of equations (18 - 19) provides an expression of the modes *<sup>φ</sup>*(*m*) *<sup>j</sup>* . The vertical velocity at the free surface is then

$$w(\mathbf{x},t) = \sum\_{m=1}^{M} \sum\_{l=0}^{M-m} \frac{\eta^l}{l!} \frac{\partial^{l+1}}{\partial z^{l+1}} \phi^{(m)}(\mathbf{x}, z=0, t) \tag{21}$$

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

$$w(\mathbf{x}, t) = \sum\_{m=1}^{M} w^{(m)}\,\_{\text{\textquotedblleft}p} \,\_{\text{\textquotedblleft}p} \tag{22}$$

where

6 Will-be-set-by-IN-TECH

The main difficulty in this approach is the computation of the vertical velocity at the free surface, *w*. Following [28], the potential *φ*(*x*, *z*, *t*) is written in a finite perturbation series up

> *M* ∑ *m*=1

The term *<sup>φ</sup>*(*m*) is of order <sup>O</sup>(*εm*), where *<sup>ε</sup>*, a small parameter, is a measure of the wave steepness. Then expanding each *φ*(*m*) evaluated on *z* = *η* in a Taylor series about *z* = 0,

> *ηl l*! *∂*(*l*) *∂z<sup>l</sup> φ*(*m*)

At a given instant of time, *φ<sup>s</sup>* and *η* are known, and we can estimate *φ*(*m*) at each order *εm*:

*φ*(1)

*ηl l*! *∂l ∂z<sup>l</sup>*

The boundary conditions, together with the Laplace equations <sup>∇</sup>2*φ*(*m*) <sup>=</sup> 0 define a series

∞ ∑ *j*=1

Note that *φ*(*m*)(*x*, *z*, *t*) satisfies automatically the Laplace equation and the boundary condition

Substitution of equation (20) into the set of equations (18 - 19) provides an expression of the

*ηl l*!

*∂l*+<sup>1</sup> *<sup>∂</sup>zl*+<sup>1</sup> *<sup>φ</sup>*(*m*)

*φ*(*m*−*l*)

*<sup>φ</sup>*(*m*)

*m*−1 ∑ *l*=1

(*x*, *z*, *t*) =

*φ*(*m*)

 1 + *η*<sup>2</sup> *x* + 2 *νk*<sup>0</sup> *c*0

> *φsxx* − *wηxx* 1 + *η*<sup>2</sup> *x*

> > (*x*, *z* = 0, *t*)

. For 2*π*−periodic conditions in *<sup>x</sup>*, say, *<sup>φ</sup>*(*m*) can be written as

(*x*, *z* = *η*(*x*, *t*), *t*) (15)

(*x*, *z* = 0, *t*) = *φs*(*x*, *t*), *m* = 1, (18) . . .

(*x*, *z* = 0, *t*), *m* ≥ 2. (19)

*<sup>j</sup> <sup>e</sup>*<sup>−</sup>*jzeijx*. (20)

(*x*, *z* = 0, *t*) (21)

(*x*, *z*, *t*). (16)

*ηxx* on *z* = *η*(*x*, *t*) (13)

− *Pa* on *z* = *η*(*x*, *t*) (14)

. (17)

*∂η*

*∂φ<sup>s</sup>*

to a given order *M*,

we obtain

with

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>−</sup>*<sup>η</sup>* <sup>−</sup> <sup>1</sup>

2 *φ*2 *sx* + 1 2 *w*2 1 + *η*<sup>2</sup> *x* + 2 *νk*<sup>0</sup> *c*0

*φs*(*x*, *t*) =

(*x*, *z* = 0, *t*) = −

*φ*(*m*)

*<sup>j</sup>* . The vertical velocity at the free surface is then

*M* ∑ *m*=1 *M*−*m* ∑ *l*=0

*w*(*x*, *t*) =

*φ*(*m*)

**4.2. Computation of the vertical velocity**

of Dirichlet problems for *φ*(*m*)

<sup>∇</sup>*φ*(*m*) <sup>→</sup> 0 when *<sup>z</sup>* → −∞.

modes *<sup>φ</sup>*(*m*)

follows in deep water

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>−</sup>*φsxη<sup>x</sup>* <sup>+</sup> *<sup>w</sup>*

*<sup>w</sup>* <sup>=</sup> *∂φ ∂z*

*φ*(*x*, *z*, *t*) =

*M* ∑ *m*=1 *M*−*m* ∑ *l*=0

$$w^{(m)} = \sum\_{l=0}^{m-1} \frac{\eta^l}{l!} \frac{\partial^{l+1}}{\partial z^{l+1}} \phi^{(m-l)}(x, z=0, t). \tag{23}$$

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 (*η*¯, *φ*¯) which is computed using the approach first introduced by [30]. A very high-order Stokes wave of amplitude *a*<sup>0</sup> and wavenumber *k*<sup>0</sup> is calculated iteratively. In the absence of wind and dissipation, the infinitesimal perturbation components (*η*� , *φ*� ) calculated through a perturbative approach developed by [31] correspond to a Benjamin-Feir instability of wavenumber *δk*. The perturbed Stokes wave is obtained by adding the infinitesimal perturbations at the sidebands *k*<sup>0</sup> ± *δk* of the fundamental and its harmonics. For fixed values of (A, Ω) two kinds of initial conditions are used when wind and dissipation are considered. The first kind (unseeded case) corresponds to the unperturbed Stokes' wave (*η*, *φ*)=(*η*¯, *φ*¯), whereas the second kind (seeded case) corresponds to the perturbed Stokes' wave (*η*, *φ*)=(*η*¯, *φ*¯) + *ε*(*η*� , *φ*� + *φ*¯ *zη*� ), with *ε* = 10−3. In both cases, we consider a Stokes wavetrain such as *a*0*k*<sup>0</sup> = 0.11 and *k*<sup>0</sup> = 5. The wavenumber of the modulational instability is *δ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 *kmax* = 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) × *kmax*. 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.

0 100 200 300 400 500 600 700 800 900 1000

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

t/T

**Figure 2.** Time evolution of the normalized amplitudes of the fundamental mode (*k* = 5), subharmonic mode (*k* = 4) and superharmonic mode (*k* = 6) for (A, Ω)=(4, 0.61). 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

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

Figure 1 shows the time evolution of the normalized amplitudes *a*(*t*)/*a*<sup>0</sup> 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 A = 4 and to a viscosity parameter Ω = 0.59. Within the framework of the NLS equation, [4] showed that the wave train is unstable to modulational instability for these values of A and Ω. 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* = 900*T*. On figure 3 one can observe the persistence of the modulational instability through the evolution of the free

figures correspond to two different conditions of wind forcing and damping.

0

case). *T* is the fundamental wave period.

0.2

0.4

0.6

0.8

a(t)/a0

1

1.2
