HIGH-ORDER NONLINEARSCHRÖDINGER EQUATION FOR THE ENVELOPEOF SLOWLY MODULATED GRAVITY WAVESON THE SURFACE OF FINITE-DEPTH FLUIDAND ITS QUASI-SOLITON SOLUTIONS

# High-Order Nonlinear Schrödinger Equation for the Envelope of Slowly Modulated Gravity Waves on the Surface of Finite-Depth Fluid and Its Quasi-Soliton Solutions

## Abstract

We consider the high-order nonlinear Schrödinger equation derived earlier by Sedletsky [Ukr.  J.  Phys.  48(1), 82 (2003)] for the first-harmonic envelope of slowly modulated gravity waves on the surface of finite-depth irrotational, inviscid, and incompressible fluid with flat bottom.  This equation takes into account the third-order dispersion and cubic nonlinear dispersive terms.  We rewrite this equation in dimensionless form featuring only one dimensionless parameter , where is the carrier wavenumber and is the undisturbed fluid depth.  We show that one-soliton solutions of the classical nonlinear Schrödinger equation are transformed into quasi-soliton solutions with slowly varying amplitude when the high-order terms are taken into consideration.  These quasi-soliton solutions represent the secondary modulations of gravity waves.

nonlinear Schrödinger equation, gravity waves, finite depth, slow modulations, wave envelope, quasi-soliton, multiple-scale expansions.
###### pacs:
47.35.Bb
\udk

532.591 \razd\secx \autorcolI.S. Gandzha, Yu.V. Sedletsky, D.S. Dutykh

## I Introduction

The nonlinear Schrödinger equation (NLSE)

 Aτ=−a1Aχ−ia2Aχχ+ia0,0,0A|A|2 (1)

arises in describing nonlinear waves in various physical contexts, such as nonlinear optics Turitsyn_Review_2012 (), plasma physics Infeld_Rowlands_2000 (), nanosized electronics Crutcher_2012 (), ferromagnetics chen1993 (), Bose–Einstein condensates Zakharov_Review_2012 (), and hydrodynamics Zakharov_1971 (); Yuen_Lake_1975 (); Dias_Kharif_Review_1999 (); Onorato_Review_2013 ().  Here, is the direction of wave propagation, is time, is the complex first-harmonic envelope of the carrier wave, and the subscripts next to denote the partial derivatives.  NLSE takes into account the second-order dispersion (term with ) and the phase self-modulation (term with ).  The coefficients , , and take various values depending on the particular physical context under consideration.

In the general context of weakly nonlinear dispersive waves, this equation was first discussed by Benney and Newell Benney_Newell_1967 ().  In the case of gravity waves propagating on the surface of infinite-depth irrotational, inviscid, and incompressible fluid, NLSE was first derived by Zakharov Zakharov_1968 () using the Hamiltonian formalism and then by Yuen and Lake Yuen_Lake_1975 () using the averaged Lagrangian method.  The finite-depth NLSE of form (1) was first derived by Hasimoto and Ono Hasimoto_Ono_1972 () using the multiple-scale method and then by Stiassnie and Shemer Stiassnie_Shemer_1984 () from Zakharov’s integral equations. Noteworthy is also the recent paper by Thomas et al. Thomas_Kharif_2012 () who derived the finite-depth NLSE for water waves on finite depth with constant vorticity.

Under certain relationship between the parameters, when

 a2a0,0,0 < 0, (2)

NLSE admits exact solutions in the form of solitons, which exist due to the balance of dispersion and nonlinearity and propagate without changing their shape and keeping their energy Dodd_solitons ().  In this case, the uniform carrier wave is unstable with respect to long-wave modulations allowing for the formation of envelope solitons.  This type of instability is known as the modulational or Benjamin–Feir instability Zakharov_MI_2009 () (it was discovered for the first time in optics by Bespalov and Talanov Bespalov_1966 ()).  In the case of surface gravity waves, condition (2) holds at , being the carrier wavenumber and being the undisturbed fluid depth.  In addition to theoretical predictions, envelope solitons were observed in numerous experiments performed in water tanks Yuen_Lake_1975 (); Yuen_Lake_1982 (); Su_1982 (); Shemer_2010 (); Chabchoub_2011 (); Onorato_Review_2013 (); Slunyaev_PF_2013 (); Shemer_2013 ().

At the bifurcation point (), when the modulational instability changes to stability, NLSE of form (1) is not sufficient to describe the wavetrain evolution since the leading nonlinear term vanishes.  In this case, high-order nonlinear and nonlinear-dispersive terms should be taken into account.  In the case of infinite depth, such a high-order NLSE (HONLSE) was first derived by Dysthe Dysthe_1979 ().  It includes the third-order dispersion () and cubic nonlinear dispersive terms (, , asterisk denotes the complex conjugate) as well as an additional nonlinear dispersive term describing the input of the wave-induced mean flow (some of these terms were introduced earlier by Roskes Roskes_1977 () without taking into consideration the induced mean flow).  This equation is usually referred to as the fourth-order HONLSE to emphasize the contrast with the third-order NLSE.  Janssen Janssen_1983 () re-derived Dysthe’s equation and corrected the sign at one of the nonlinear dispersive terms. Hogan Hogan_1985 () followed the earlier work by Stiassnie Stiassnie_1984 () to derive the similar equation for deep-water gravity-capillary waves with surface tension taken into account.  Selezov et al.  Selezov_2003 () extended the HONLSE derived by Hogan to the case of nonlinear wavetrain propagation on the interface of two semi-infinite fluids without taking into account the induced mean flow.  Worthy of mention is also the paper by Lukomsky Lukom_1995 () who derived Dysthe’s equation in a different way.  Later, Trulsen and Dysthe Trulsen_Dysthe_1996 () extended the equation derived by Dysthe to broader bandwidth by including the fourth- and fifth-order linear dispersion.  Debsarma and Das Debsarma_Das_2005 () derived a yet more general HONLSE that is one order higher than the equation derived by Trulsen and Dysthe.  Gramstad and Trulsen Gramstad_Trulsen_PF_2011 () derived a set of two coupled fourth-order HONLSEs capable of describing two interacting wave systems separated in wavelengths or directions of propagation.  Zakharov and Dyachenko Zakharov_Dyachenko_2010 (); Zakharov_Dyachenko_2011 (); Zakharov_Dyachenko_2012 () made a conformal mapping of the fluid domain to the lower half-plane to derive a counterpart of Dysthe’s equation in new canonical variables (the so-called compact Dyachenko–Zakharov equation Fedele_Dutykh_2012_JETP (); Fedele_Dutykh_2012 ()).

Original Dysthe’s equation was written for the first-harmonic envelope of velocity potential rather than of surface profile.  In the case of standard NLSE, this difference in not essential because in that order the first-harmonic amplitudes of the velocity potential and surface displacement differ by a dimensional factor only, which is not true anymore in the HONLSE case, as discussed by Hogan Hogan_1986 ().  Keeping this in mind, Trulsen et al. Trulsen_Dysthe_2000 () rewrote Dysthe’s equation in terms of the first-harmonic envelope of surface profile while taking into account the linear dispersion to an arbitrary order.

In the case of finite depth, the effect of induced mean flow manifests itself in the third order, so that the NLSE is generally coupled to the equation for the induced mean flow Benney_Roskes_1969 ().  However, Davey and Stewartson Davey_Stewartson_1974 () showed that these coupled equations are equivalent to the single NLSE derived by Hasimoto and Ono Hasimoto_Ono_1972 ().  On the other hand, such an equivalence is not preserved for high-order equations.  The first attempt to derive a HONLSE in the case of finite depth was made by Johnson Johnson_1977 (), but only for , when the cubic NLSE term vanishes.  The similar attempt was made by Kakutani and Michihiro Kakutani_1983 () (see also a more formal derivation made later by Parkes Parkes_1987 ()).  A general fourth-order HONLSE for the first-harmonic envelope of surface profile was derived by Brinch-Nielsen and Jonsson Brinch-Nielsen_1986 () in coupling with the integral equation for the wave-induced mean flow. Gramstad and Trulsen Gramstad_Trulsen_JFM_2011 (); Gramstad_2014 () derived a fourth-order HONLSE in terms of canonical variables that preserves the Hamiltonian structure of the surface wave problem.

Sedletsky SedletskyUJP2003 (); SedletskyJETP2003 () used the multiple-scale technique to derive a single fourth-order HONLSE for the first-harmonic envelope of surface profile by introducing an additional power expansion of the induced mean flow.  This equation is the direct counterpart of Dysthe’s equation written in terms of the first-harmonic envelope of surface profile Trulsen_Dysthe_2000 () but for the case of finite depth. Slunyaev Slunyaev_2005 () confirmed the results obtained in SedletskyJETP2003 () and extended them to the fifth order. Grimshaw and Annenkov Grimshaw_2011 () considered a HONLSE for water wave packets over variable depth.

The deep-water HONLSE in the form of Dysthe’s equation was extensively used in numerical simulations of wave evolution Lo_Mei_1985 (); Akylas_1989 (); Ablowitz_2000 (); Ablowitz_2001 (); Clamond_2006 (); Slunyaev_2009 (); Fedele_Dutykh_2011 (); FD_2012 (); Slunyaev_PRE_2013 (). However, no such modeling has been performed in the case of finite depth because of the complexity of equations as compared to the deep-water limit.  The equation derived in SedletskyUJP2003 (); SedletskyJETP2003 () can be used as a good starting point for the simulations of wave envelope evolution on finite depth.  The aim of this paper is (i) to rewrite this equation in dimensionless form suitable for numerical integration and (ii) to observe the evolution of NLSE solitons taken as initial waveforms in the case when the HONLSE terms are taken into consideration for several values of intermediate depth.

This paper is organized as follows.  In Section II, we write down the fully nonlinear equations of hydrodynamics used as the starting point in this study.  In Section III, we formulate the constraints at which the fully nonlinear equations can be reduced to HONLSE. Then we briefly outline the multiple-scale technique used to derive this equation, which is presented in Section IV.  Next, we introduce dimensionless coordinate, time, and amplitude to go over to the dimensionless HONLSE.  As a result, only one dimensionless parameter appears in the equation.  The final step is to pass to the reference frame moving with the group speed of the carrier wave.  In Section V, we present the results of numerical simulations and compare the NLSE and HONLSE solutions. Conclusions are made in Section VI.

## Ii Problem Formulation

We consider the dynamics of potential two-dimensional waves on the surface of irrotational, inviscid, and incompressible fluid under the influence of gravity.  Waves are assumed to propagate along the horizontal -axis, and the direction of the vertical -axis is selected opposite to the gravity force.  The fluid is assumed to be bounded by a solid flat bed at the bottom and a free surface at the top (Fig. 1).  The atmospheric pressure is assumed to be constant on the free surface.  Then the evolution of waves and associated fluid flows is governed by the following system of equations Stoker_1992 (); UJP2013 ():

 Φxx+Φyy=0,\leavevmode\nobreak \leavevmode\nobreak −∞
 Φt+12(Φ2x+Φ2y)+gη=0,\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak y=η(x,t); (4)
 ηt−Φy+ηxΦx=0,\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak y=η(x,t); (5)
 Φy=0,\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak y=−h; (6)

where is the velocity potential (the velocity is equal to ), is the acceleration due to gravity, is time.  Here, (3) is the Laplace equation in the fluid domain, (4) is the dynamical boundary condition (the so-called Bernoulli or Cauchy–Lagrange integral), (5) and (6) are the kinematic boundary conditions (no fluid crosses the free surface and the bottom), the indices , , and designate the partial derivatives over the corresponding variables.  The position of the zero level is selected such that the Bernoulli constant (the right-hand side of Eq. (4)) is equal to zero.

Consider a modulated wavetrain with carrier frequency and wavenumber .  In this case, a solution to Eqs.  (3)–(6) can be looked for in the form of Fourier series with variable coefficients:

 (7)

where stands for complex conjugate (here, we assume the carrier wave to be symmetric), the functions and are assumed to be real by definition. Substituting (7) in (3)–(6) and equating the coefficients at the like powers of , one can obtain a system of nonlinear partial differential equations for the functions and .  Linearization of these equations at gives the dispersion relation for gravity waves:

 ω2=gktanh(kh). (8)

## Iii Slowly Modulated Quasi-Harmonic Wavetrains and Multiple-Scale Expansions

Generally the system of equations for and is by no means more simple than original equations.  It can be simplified when solutions are looked for in a class of functions with narrow spectrum, .  In this case, the problem has a formal small parameter (quasi-monochromaticity condition), with and being slow functions of and .  Accordingly, the wave motion can be classified into slow one and fast one by introducing different time scales and different spatial scales:

 Tn≡μnt,Xn≡μnx. (9)

The derivatives with respect to time and coordinate are expanded into the following series:

 ∂∂t=∞∑n=0μn∂∂Tn,∂∂x=∞∑n=0μn∂∂Xn, (10)

the times and coordinates being assumed to be independent variables.

When there are no resonances between higher harmonics, the amplitudes of Fourier coefficients decrease with increasing number (quasi-harmonicity condition):

 ηn∼εnA,n⩾1,η0∼ε2A,ε<1, (11)

where

 η1≡12εA(x,t). (12)

The parameter can be regarded as a formal small parameter related to the smallness of wave amplitude as compared to the carrier wavelength .  In this case, the unknown functions and can be expanded into power series in the formal parameter :

 (Φn(x,y,t)ηn(x,t))=∞∑m=1εm(Φ(m)n(x,y,t)η(m)n(x,t)). (13)

Multiple-scale expansions (10) and (13) allow the functions and to be expressed in terms of the first harmonic envelope , as described in detail in SedletskyUJP2003 ().  Note that in the procedure described in SedletskyUJP2003 () it is essential to set .

In practice, the quasi-harmonicity condition can be written as

 |kη1|≪1, (14)

and the condition of slow modulation (quasi-monochromaticity) can be formalized as

 ∣∣∣AxkA∣∣∣≪1, (15)

which follows from differentiating the function with respect to .  With these conditions satisfied, the original system of equations (3)–(6) can be reduced to one evolution equation for the first harmonic envelope with the use of small-amplitude expansions (10) and (13).

## Iv High-Order Nonlinear Schrödinger Equation

### iv.1 Equation derived by Sedletsky

Sedletsky SedletskyUJP2003 (); SedletskyJETP2003 () used the above-described multiple-scale procedure to derive the following HONLSE for the first-harmonic envelope (Eq.  (68) in SedletskyUJP2003 ()):

 i(∂A∂(εt)+Vg∂A∂(εx))+ +ε(12ω′′∂2A∂(εx)2+ωk2q3|A|2A)+ +iε2(−16ω′′′∂3A∂(εx)3+ωkQ41|A|2∂A∂(εx)+ +ωkQ42A2∂A∗∂(εx))=0[m/s]. (16)

As compared to the standard NLSE, this equation takes into account additional nonlinear and dispersive terms of order .  Equation (16) was later re-derived by Slunyaev Slunyaev_2005 (), who confirmed the symbolic computations presented in SedletskyUJP2003 (); SedletskyJETP2003 () and extended them to the order.  Here, we restrict our attention to the original equation (16).  The parameters of this equation are given by

 ω=(gkσ)1/2,σ≡tanh(kh), (17a) ω′=∂ω∂k≡Vg=ω2k(1+2khsinh(2kh))= =ω2k(1+1−σ2σkh), (17b) ω′′=∂2ω∂k2=ω4k2σ2((σ2−1)(3σ2+1)k2h2− −2σ(σ2−1)kh−σ2), (17c) ω′′′=∂3ω∂k3=−ω8k3σ3((σ2−1)× ×(15σ4−2σ2+3)k3h3−3σ(σ2−1)(3σ2+1)k2h2− −3σ2(σ2−1)kh−3σ3), (17d) q3=−116σ4ν((σ2−1)2(9σ4−10σ2+9)k2h2+ +2σ(3σ6−23σ4+13σ2−9)kh− −σ2(7σ4−38σ2−9)), (17e) Q41=132σ5ν2((σ2−1)5× ×(3σ6−20σ4−21σ2+54)k5h5− −σ(σ2−1)3(11σ8−99σ6−61σ4+7σ2+270)k4h4+ +2σ2(σ2−1)(7σ10−58σ8+38σ6+52σ4− −181σ2+270)k3h3−2σ3(3σ10+18σ8−146σ6− −172σ4+183σ2−270)k2h2−σ4(σ8−109σ6+517σ4+ +217σ2+270)kh+σ5(σ6−40σ4+193σ2+54))+Δ, (17f) Q42=132σ5ν2(−(σ2−1)5× ×(3σ6+7σ4−11σ2+9)k5h5+ +σ(σ2−1)3(11σ8−48σ6+66σ4+8σ2+27)k4h4− −2σ2(σ2−1)(7σ10−79σ8+282σ6− −154σ4−σ2+9)k3h3+2σ3(3σ10−63σ8+314σ6− −218σ4+19σ2+9)k2h2+σ4(σ8+20σ6−158σ4− −28σ2−27)kh−σ5(σ6−7σ4+7σ2−9))−Δ, (17g) ν=(σ2−1)2k2h2−2σ(σ2+1)kh+σ2. (17h)

The quantity is the wave group speed.  The parameter is the correction introduced by Slunyaev Slunyaev_2005 () to the coefficients derived in SedletskyUJP2003 (); SedletskyJETP2003 ().  This correction is negligible at (see Appendix A), and we ignore it by keeping .

The free-surface displacement is expressed in terms of as

 η=ε2η0+εRe(Aei(kx−ωt))+
 +ε22Re(η2e2i(kx−ωt))+O(ε3), (18)

where stands for the real part of a complex-valued function.  Here, and are defined as

 η0=σ+2(1−σ2)khνk|A|2, (19a) η2=3−σ28σ3kA2. (19b)

The corresponding velocity potential is written as

 Φ=εΦ0+ε2Re(Φ1ei(kx−ωt))+
 +ε22Re(Φ2e2i(kx−ωt))+O(ε3), (20)

where

 −(y+h)∂A∂xsinh(k(y+h))cosh(kh)), (21a) Φ2=3iω(σ4−1)16σ4cosh(2k(y+h))cosh(2kh)A2. (21b)

The term describes the wave-induced mean flow and is expressed implicitly in terms of its derivatives

 ∂Φ0∂x=εωkγ12σν|A|2+iεωγ28σ2ν2(A∂A∗∂x−A∗∂A∂x), (22a) ∂Φ0∂t=−Vg∂Φ0∂x, (22b)

where

 (23a) γ2=(σ2−1)5k4h4+4σ(σ2−1)2(13σ2+3)k3h3− −2σ2(σ2−1)(3σ4+32σ2−3)k2h2+ +4σ3(2σ4−σ2−5)kh−3σ4(σ2−5). (23b)

Functions (18) and (20) define an approximate solution to the original system of equations (3)–(6) in terms of the first-harmonic envelope , which is found from Eq. (16).

### iv.2 Dimensionless form

Introduce the following dimensionless time, coordinate, and amplitude:

 τ=βt,χ=kx,u=α−1εA, (24)

and being the parameters to be determined.  The relationship between the old and new derivatives is

 ∂∂x=k∂∂χ,∂∂t=β∂∂τ.

Then Eq. (16) is transformed to

 iα(βuτ+kVguχ)+α(12ω′′k2uχχ+ωk2q3|α|2|u|2u)+
 +iα(−16ω′′′k3uχχχ+ωk2Q41|α|2uχ|u|2+
 +ωk2Q42|α|2u2u∗χ)=0[m/s].

Here, the indices and designate the partial derivatives with respect to the corresponding variables.  Taking into account that at all (Fig.  2), divide this equation by so that

 i(βω′′k2uτ+Vgω′′kuχ)+(12uχχ+ωω′′|α|2q3|u|2u)+
 +i(−16ω′′′kω′′uχχχ+ωω′′|α|2Q41uχ|u|2+
 +ωω′′|α|2Q42u2u∗χ)=0

and select the values of and as

 |α|2=−ω′′ω>0,β=−ω′′k2>0. (25)

Thus, Eq.  (16) takes the dimensionless form

 i(uτ−Vgω′′kuχ)−12uχχ+q3|u|2u+
 +i(16ω′′′kω′′uχχχ+Q41uχ|u|2+Q42u2u∗χ)=0

or, equivalently,

 i(uτ+a1uχ)−a2uχχ+a0,0,0|u|2u+i(−a3uχχχ+ +a1,0,0uχ|u|2+a0,0,1u2u∗χ)=0, (26a) which finally yields uτ=−a1uχ−ia2uχχ+ia0,0,0|u|2u+ +(a3uχχχ−a1,0,0uχ|u|2−a0,0,1u2u∗χ), (26b)

where we used the unified notation introduced by Lukomsky and Gandzha ujp09en (). Here, the coefficients

 a1=−Vgω′′k=−2υ(σ2+σ(1−σ2)kh)>0,a2=12,a3≡−16ω′′′kω′′==112συ((σ2−1)(15σ4−2σ2+3)k3h3−−3σ(σ2−1)(3σ2+1)k2h2−−3σ2(σ2−1)kh−3σ3),a0,0,0≡q3,a1,0,0≡Q41,a0,0,1≡Q42,υ=(σ2−1)(3σ2+1)k2h2−2σ(σ2−1)kh−σ2 (27)

are all real and depend on one dimensionless parameter .  Their behavior as functions of is shown in Fig.  3.  It can be seen that Eq.  (16) is valid at , where the coefficients , , and do not diverge.  At smaller depths, the Korteweg–de Vries equation and its generalizations Johnson_2002 (); Kharif_Enpeli_2003 () should be used.  On the other hand, at large , the infinite-depth limit (Dysthe’s equation) should be used.  Indeed, the following asymptotics are easily obtained at :

 a3=14,a0,0,0=−12,a1,0,0=32,a0,0,1=14. (28)

They coincide with the corresponding coefficients of Dysthe’s equation Trulsen_Dysthe_2000 (), except for the term including the wave-induced mean flow, which cannot be explicitly reconstructed from Eq.  (26) because of the additional power expansion of the wave-induced mean flow made to derive Eq.  (16).  However, this term can be reconstructed from the equations generating Eq.  (16), at the stage when the wave-induced mean flow has not been excluded from the equation for yet SedletskyUJP2003 ().  Taking into account these constraints, we will restrict our attention to the following range of intermediate depths:

 1

### iv.3 Moving reference frame

Equation (26) can be rewritten in the form without the term.  To this end, let us proceed to the reference frame moving with speed (dimensionless group speed):

 ξ=χ−a1τ,T=τ. (30)

The relationship between the derivatives in new and old variables is given by the formulas

 ∂∂χ=∂ξ∂χ∂∂ξ+∂T∂χ∂∂T=∂∂ξ,
 ∂∂τ=∂ξ∂τ∂∂ξ+∂T∂τ∂∂T=−a1∂∂ξ+∂∂T,

so that

 uτ=−ia2uξξ+ia0,0,0|u|2u+
 +(a3uξξξ−a1,0,0uξ|u|2−a0,0,1u2u∗ξ). (31)

This is our target equation for numerical simulations.  It possesses the integral of motion

 I0(τ)=∞∫−∞|u(ξ,τ)|2dξ=const, (32)

which expresses the conservation of wave action.  The derivation of this conservation law is given in Appendix C.  It allows one to trace the relative numerical error of simulations:

 Er(I0)=|I0(τ)−I0(0)|I0(0). (33)

Of particular interest is to reveal any relationship of Eq.  (31) to other HONLSEs derived in different contexts.  In Appendix B, we consider one such equation (the Sasa–Satsuma equation) and prove that Eq.  (31) cannot be reduced to it at any .

### iv.4 Dimensionless free surface displacement and velocity potential

The dimensionless free surface displacement is expressed in terms of as follows:

 ζ≡kη=α0|u|2+α1Re(ueiθ)+2α2Re(u2e2iθ), (34)
 α0=σ+2(1−σ2)khcν,α1=1√c,α2=3−σ28cσ3,

where

 θ=kx−ωt=χ−cτ=ξ+(a1−c)τ (35)

is the wave phase and

 c=1|α|2k2=−4σ2υ (36)

is the dimensionless phase speed.  Figure 4 shows the ratio of the phase speed to the group speed as a function of .  This ratio is equal to unity at , and it is twice as large at , in full conformity with the classical water wave theory UJP2013 ().  The wave envelope is written as

 [ζ]envelope=α1|u|+(α0+2α2)|u|2. (37)

The corresponding dimensionless velocity potential is expressed as

 φ≡−1ω′′Φ=φ0+2Re(φ1eiθ)+
 +2Re(φ2e2iθ), (38)
 (φ0)τ=−a1(φ0)ξ,
 φ1=√c2σ((−iu+(σ2+1)kh+σ2σuξ)×
 ×cosh(z+kh)cosh(kh)−(z+kh)sinh(z+kh)cosh(kh)uξ),
 φ2=3i(σ4−1)16σ4cosh(2(z+kh))cosh(2kh)u2,

where is the dimensionless vertical coordinate.  The quasi-harmonicity condition is written as

 |u|√c≪1, (39)

and the quasi-monochromaticity condition is

 ∣∣∣uξu∣∣∣≪1. (40)

Finally, the original equations of hydrodynamics can be written in the following dimensionless form:

 φξξ+φzz=0,\leavevmode\nobreak −∞<ξ<∞,\leavevmode\nobreak −kh
 φτ+12(φ2ξ+φ2z)+c2σζ=0,\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak z=ζ(ξ,τ); (42)
 ζτ−φz+ζξφξ=0,\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak z=ζ(ξ,τ); (43)
 φz=0,\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak z=−kh. (44)

## V Numerical Simulations

In this section, we adopt the split-step Fourier (SSF) technique described in Appendix D to compute solutions to HONLSE (31).  To test the accuracy of our numerical scheme, we start from classical NLSE (1) written in terms of the coordinate .  At (), it has an exact one-soliton solution Zakharov_1971 ():

 u(ξ,τ)=u0exp(iκξ−iΩτ)cosh(K(ξ−ξ0−Vτ)), (45)
 Ω=(K2−κ2)a2,V=−2κa2,K=|u0|√−a0,0,02a2,
 u0∈\mathdsC,κ,ξ0∈\mathdsR.

Here, is the soliton speed, is the complex amplitude, and are the soliton’s wavenumber and frequency, and is the soliton’s initial position.  The amplitude and wavenumber should be selected such that the quasi-harmonicity and quasi-monochromaticity conditions (39), (40) hold true.  In practice, these conditions mean that the soliton amplitude and wavenumber should be small:

 |u0|≪1,κ≪1.

In this study, we restrict our attention by the following choice of parameters:

 u0=0.1,κ=−K(⇒Ω=0),ξ0=0. (46)

Figures 5 and 6 demonstrate that constraints (39) and (40) are readily satisfied in this case.  Note that at we have . In this case, solitons move from left to right with speed exceeding the carrier group speed.

Figure 7 shows a soliton computed for using analytical formula (45) for the initial moment and moment .  The same soliton was taken as the initial condition for the simulation with the SSF technique.  The deviation from the exact solution is seen to be negligible.  Indeed, the numerical error estimated with formula (33) is

 S(2)|τ=10000:\leavevmode\nobreak \leavevmode\nobreak Er(I0)=1.3×10−10%,
 Δrms(uexact,ucomp)=1.0×10−4%,
 S(4)|τ=10000:\leavevmode\nobreak \leavevmode\nobreak Er(I0)=2.3×10−10%,
 Δrms(uexact,ucomp)=3.5×10−9%,

where and designate the order of the SSF technique adopted for calculation (see Appendix D) and

 Δrms(u,g)(τ)= ⎷∫∞−∞(|u(ξ,τ)|−|g(ξ,τ)|)2dξI0

is the relative r.m.s. deviation between two functions.  Thus, our numerical scheme reproduces the exact one-soliton solution to NLSE with high accuracy.

Figure 8 shows the evolution of the same one-soliton waveform taken as the initial condition in HONLSE (31).  As compared to the NLSE case, the wave amplitude is smaller, the pulse width is larger, and the wave speed is higher. The wave amplitude does not remain constant and exhibits slow oscillations that can be interpreted as the secondary modulation of the carrier wave.  The amplitude of these oscillations decreases with time (Fig.  9).  Such a solution does not fall under the definition of soliton because it does not preserve the constant amplitude and shape during the evolution.  On the other hand, it moves with nearly the constant speed (Fig.  10) and still possesses the unique property of solitons to exist over long periods of time without breaking.  In view of this unique property, we call such solutions quasi-solitons. The term quasi-soliton was introduced earlier by Zakharov and Kuznetsov Zakharov_1998 (), but in somewhat different context; then Karpman et al.  Karpman_2001 () and Slunyaev Slunyaev_2009 () used it in the same context as in the present study.

Such a behavior of NLSE solitons in the HONLSE case was first described by Akylas Akylas_1989 () in the context of asymptotic modeling and numerical simulations of Dysthe’s equation in the infinite-depth limit.  Growth in the soliton speed corresponds to the well-known carrier frequency downshift observed in deep-water experiments by Su Su_1982 () and in simulations of Dysthe’s equation by Lo and Mei Lo_Mei_1985 ().  Dysthe Dysthe_1979 () pointed out that this phenomenon originates due to the wave-induced mean flow, whose component in the direction of propagation of the wave causes a local Doppler shift.  Here, we proved for the first time that this well-known phenomenon can be observed on finite depth as well.  This result is the main practical achievement of our study.  Figures  11, 12, and 13 demonstrate that the same quasi-soliton solution and frequency downshift are observed at a smaller depth, .

Finally, the free surface profile reconstructed with formula (34) is shown in Fig. 14 for . The dimensionless maximum free surface elevation is about . The case corresponds to wavelengths twice as large as depth, .  The typical depth of the shelf near the north-west shore of the Black Sea varies from 10 to 100 m.  Hence, the wavelength corresponding to falls within the range from 20 to 200 m, which is quite typical of water waves observed on the Black Sea.  For  m, we have  m and .  The corresponding maximum free surface elevation of the wave shown in Fig. 14 is about  m, and the significant wavetrain width is about 2 km.  Thus, quasi-soliton solutions obtained in this study can describe swells propagating on the relatively calm background on seas with intermediate depths.  The typical trough-to-crest height of such swells is about 1 m.

## Vi Conclusions

The HONLSE derived earlier by Sedletsky SedletskyUJP2003 () for the first-harmonic envelope of slowly modulated gravity waves on the surface of finite-depth irrotational, inviscid, and incompressible fluid with flat bottom was rewritten in the dimensionless form suitable for numerical simulations.  One-soliton solutions to NLSE are transformed into quasi-soliton solutions with slowly varying amplitude when the HONLSE terms are taken into consideration.  These quasi-solitons represent the secondary modulations of gravity waves.  They propagate with nearly constant speed and possess the unique property of solitons to exist over long periods of time without breaking.  Their speed was found to be higher than the speed of the NLSE solitons taken as initial conditions in computations.  This phenomenon was observed earlier both in experiment and numerical modeling in the case of deep-water limit Su_1982 (); Akylas_1989 ().  It is related to the frequency downshift originating due to the wave-induced mean flow Dysthe_1979 (); Lo_Mei_1985 ().  The quasi-soliton solutions obtained in this study describe swells propagating on the relatively calm background on seas with intermediate water depth.

The authors are grateful to Dr.  S. S. Rozhkov for initial discussions that motivated us to undertake this study. D. Dutykh would like to acknowledge the hospitality of Institut für Analysis, Johannes Kepler Universität Linz, where this work was performed.

## Appendix A On the Correction Introduced by Slunyaev

Slunyaev Slunyaev_2005 () re-derived HONLSE (16) and introduced a correction,

 Δ=−116σ3ν((σ2−1)4(3σ2+1)k3h3−
 −σ(σ2−1)2(5σ4−18σ2−3)k2h2+
 +σ2(σ2−1)2(σ