Hyperbolicity of the modulation equations for the Serre-Green-Naghdi model

TL;DR

This paper derives modulation equations for the Serre-Green-Naghdi water wave model and demonstrates their strict hyperbolicity, indicating stable wave train behavior, supported by numerical tests.

Contribution

The paper introduces the derivation of modulation equations for the SGN model and proves their strict hyperbolicity for all wave amplitudes, confirming modulational stability.

Findings

Modulation equations for SGN are strictly hyperbolic at all amplitudes.

Numerical tests confirm the modulational stability of wave trains.

SGN equations are a justified approximation of water wave dynamics.

Abstract

Serre-Green-Naghdi equations (SGN equations) is the most simple dispersive model of long water waves having "good" mathematical and physical properties. First, the model is a mathematically justified approximation of the exact water wave problem. Second, the SGN equations are the Euler-Lagrange equations coming from Hamilton's principle of stationary action with a natural approximate Lagrangian. Finally, the equations are Galilean invariant which is necessary for physically relevant mathematical models. We have derived the modulation equations to the SGN model and show that they are strictly hyperbolic for any wave amplitude, i.e., the periodic wave trains are modulationally stable. Numerical tests for the full SGN equations are shown. The results confirm the modulational stability analysis.

Equations

$$h_t+(hu{)}_x=0,$$

$$(hu{)}_t+(h{u}^2+p{)}_x=0,$$

$${\left(he\right)}_t+{\left(hue+pu\right)}_x=0,$$

$$p=\frac{g{h}^2}{2}+\frac{1}{3}{h}^2\frac{D^{2}h}{D{t}^2},e=\frac{u^2}{2}+\frac{gh}{2}+\frac{1}{6}{\left(\frac{Dh}{Dt}\right)}^2,\frac{D}{Dt}=\frac{\partial }{\partial t}+u\frac{\partial }{\partial x}.$$

$$h(\xi ,X,T)$$

$$u(\xi ,X,T)$$

$$p(\xi ,X,T)$$

$$e(\xi ,X,T)$$

$$\displaystyle-Dh_{0\xi}+(h_{0}u_{0})_{\xi}+\varepsilon\big{(}h_{0T}+(h_{0}u_{0})_{X}\big{)}=-\varepsilon\big{[}-Dh_{1\xi}+(h_{0}u_{1}+h_{1}u_{0})_{\xi}\big{]}+O(\varepsilon^{2}),$$

$$\displaystyle-D(h_{0}u_{0})_{\xi}+(h_{0}u_{0}^{2}+p_{0})_{\xi}+\varepsilon\big{(}(h_{0}u_{0})_{T}+(h_{0}u_{0}^{2}+p_{0})_{X}\big{)}=$$

$$\displaystyle\qquad-\varepsilon\big{[}-D(h_{0}u_{1}+h_{1}u_{0})_{\xi}+(2h_{0}u_{0}u_{1}+h_{1}u_{0}^{2}+p_{1})_{\xi}\big{]}+O(\varepsilon^{2}),$$

$$\displaystyle-D(h_{0}e_{0})_{\xi}+\left(h_{0}u_{0}e_{0}+p_{0}u_{0}\right)_{\xi}+\varepsilon\big{(}(h_{0}e_{0})_{T}+\left(h_{0}u_{0}e_{0}+p_{0}u_{0}\right)_{X}\big{)}=$$

$$\displaystyle\qquad\varepsilon\big{[}-D(h_{1}e_{0}+h_{0}e_{1})_{\xi}+(h_{1}u_{0}e_{0}+h_{0}u_{1}e_{0}+h_{0}u_{0}e_{1}+p_{0}u_{1}+p_{1}u_{0})_{\xi}\big{]}+O(\varepsilon^{2}).$$

$$\frac{\partial }{\partial t}=-D\frac{\partial }{\partial \xi }+\epsilon \frac{\partial }{\partial T},\frac{\partial }{\partial x}=\frac{\partial }{\partial \xi }+\epsilon \frac{\partial }{\partial X}.$$

$$\frac{D}{Dt}=\displaystyle\frac{\partial}{\partial t}+u\displaystyle\frac{\partial}{\partial x}=\big{(}u-D\big{)}\displaystyle\frac{\partial}{\partial\xi}+\varepsilon\Big{(}\displaystyle\frac{\partial}{\partial T}+u\displaystyle\frac{\partial}{\partial X}\Big{)}.$$

$$-D{h}_{0\xi }+(h_0{u}_0{)}_{\xi }=0,$$

$$-D(h_0{u}_0{)}_{\xi }+(h_0{u}_0^2+p_0{)}_{\xi }=0,$$

$$-D(h_0{e}_0{)}_{\xi }+{\left(h_0{u}_0{e}_0+p_0{u}_0\right)}_{\xi }=0.$$

$$\displaystyle h_{0T}+(h_{0}u_{0})_{X}=-\big{[}-Dh_{1\xi}+(h_{0}u_{1}+h_{1}u_{0})_{\xi}\big{]},$$

$$\displaystyle(h_{0}u_{0})_{T}+(h_{0}u_{0}^{2}+p_{0})_{X}=-\big{[}-D(h_{0}u_{1}+h_{1}u_{0})_{\xi}+(2h_{0}u_{0}u_{1}+h_{1}u_{0}^{2}+p_{1})_{\xi}\big{]},$$

$$\displaystyle(h_{0}e_{0})_{T}+\left(h_{0}u_{0}e_{0}+p_{0}u_{0}\right)_{X}=\big{[}-D(h_{1}e_{0}+h_{0}e_{1})_{\xi}+(h_{1}u_{0}e_{0}+h_{0}u_{1}e_{0}+h_{0}u_{0}e_{1}+p_{0}u_{1}+p_{1}u_{0})_{\xi}\big{]}.$$

$$(\overline{h_0}{)}_T+(\overline{h_0{u}_0}{)}_X=0,$$

$$(\overline{h_0{u}_0}{)}_T+(\overline{h_0{u}_0^2}+\overline{p_0}{)}_X=0,$$

$$(\overline{h_0{e}_0}{)}_T+{\left(\overline{h_0{u}_0{e}_0}+\overline{p_0{u}_0}\right)}_X=0.$$

$$-D{h}^{\prime }+(hu{)}^{\prime }=0,$$

$$-D(hu{)}^{\prime }+(h{u}^2+p{)}^{\prime }=0,$$

$$-D(he{)}^{\prime }+{\left(hue+pu\right)}^{\prime }=0.$$

$$-D{h}^{\prime }+(hu{)}^{\prime }=0.$$

$$h(u-D)=m=const.$$

$$p=i-\frac{m^2}{h},i=const.$$

$$\frac{D^{2}h}{Dt^{2}}=(u-D)\big{(}(u-D)h^{\prime}\big{)}^{\prime}=\frac{(u-D)h}{h}\left(\frac{(u-D)h}{h}h^{\prime}\right)^{\prime}=\frac{m}{h}\left(\frac{mh^{\prime}}{h}\right)=\frac{m^{2}}{h}\left(\frac{h^{\prime}}{h}\right)^{\prime}.$$

$$\frac{D^{2}h}{D{t}^2}=\frac{m^2}{h}{\left(\frac{h^{\prime }}{h}\right)}^{\prime }.$$

$$p=\frac{g{h}^2}{2}+\frac{1}{3}{m}^{2}h{\left(\frac{h^{\prime }}{h}\right)}^{\prime }.$$

$$\frac{m^2}{h}+\frac{g{h}^2}{2}+\frac{1}{3}{m}^{2}h{\left(\frac{h^{\prime }}{h}\right)}^{\prime }=i.$$

$$\frac{h^{\prime }}{h^3}+\frac{g{h}^{\prime }}{2m^2}+\frac{1}{3}\frac{h^{\prime }}{h}{\left(\frac{h^{\prime }}{h}\right)}^{\prime }=\frac{i{h}^{\prime }}{m^2{h}^2},$$

$$\frac{1}{6}{\left[{\left(\frac{h^{\prime }}{h}\right)}^2\right]}^{\prime }+\frac{h^{\prime }}{h^3}+\frac{g{h}^{\prime }}{2m^2}=\frac{i{h}^{\prime }}{m^2{h}^2}.$$

$$\frac{1}{6}{\left(\frac{h^{\prime }}{h}\right)}^2-\frac{1}{2h^2}+\frac{gh}{2m^2}=-\frac{i}{m^{2}h}+ϵ,ϵ=const,$$

$${\left(h^{\prime }\right)}^2=3-\frac{6i}{m^2}h+6ϵh^2-\frac{3g}{m^2}{h}^3.$$

$$F_3(h)=3-\frac{6i}{m^2}h+6ϵh^2-\frac{3g}{m^2}{h}^3,$$

Full text

We have derived the modulation equations to the SGN model and show that they are strictly hyperbolic for any wave amplitude, i.e., the periodic wave trains are modulationally stable. Numerical tests for the full SGN equations are shown. The results confirm the modulational stability analysis.

1 Introduction

One usually uses two methods to obtain the modulation equations for a reversible dispersive system: the Whitham averaged Lagrangian method [44] or averaging of the corresponding conservation laws (cf. [5, 28]). Both of them give the same system for slowly varying wave train characteristics (wave length and amplitude, for example). If the corresponding system of modulation equations is hyperbolic (elliptic), one says that the corresponding wave trains are modulationally stable (unstable). The relation between the modulational instability and, for example, the classical spectral instability is not completely understood. An important result for a class of Hamiltonian systems has been obtained in [4]: the hyperbolicity of modulation equations is necessary for the spectral stability of periodic traveling waves.

The modulational instability (wavetrain instability) is often generically called “Benjamin–Feir instability” even if, formally, this last instability concerns only about the surface gravity water waves (Benjamin and Feir [2] and Zakharov [45] in the case of deep-water waves, and Benjamin [3] in the case of finite-depth waves). We refer interested readers to the article [46] where the history of the modulational instability theory is presented. In the case of small amplitudes, the hyperbolicity (ellipticity) condition can easily be formulated in terms of the non-linear amplitude dependent dispersion relation (see [44], chapter $\displaystyle 15$). A recent application of such an approach can be found in [36] for a “conduit” equation. However, the study of hyperbolicity of the modulation equations in the case of large amplitude solutions is a more difficult problem. For integrable systems, the hyperbolicity of modulation equations and existence of the Riemann invariants were established, for example, for Korteweg-de Vries equation (KdV equation) [44], for nonlinear Schrödinger equation (NLS equation) [39], for sine-Gordon equation [14, 25], for Benjamin-Ono equation [6] (see also the book [28] for further references). For non-integrable systems, one can cite [26] where the Whitham equation was studied and modulational instability for short enough waves was shown, or [33] and [27] where the the regions of modulational and spectral stability for roll waves to the Saint-Venant equations were determined. In general, a “right choice” of unknowns in which the modulation equations are written is necessary to have explicit (or almost explicit) expressions of the corresponding characteristic values. Such a choice is not at all obvious.

The aim of this work is to study the modulational stability of periodic waves to Serre-Green-Naghdi equations (SGN equations). One-dimensional SGN equations can be written in the Eulerian coordinates in the form [42, 43, 21, 22]:

[TABLE]

with

[TABLE]

Here $\displaystyle h$ is the fluid depth, $\displaystyle u$ is the averaged over the fluid depth velocity, $\displaystyle p$ is the integrated over the fluid depth pressure. If $\displaystyle L_{0}$ is a characteristic wave length, and $\displaystyle H_{0}$ is the characteristic water depth, we define the dimensionless small parameter $\displaystyle\beta=H_{0}^{2}/L_{0}^{2}$. The SGN equations are obtained by depth-averaging the Euler system and keeping in the resulting set of equations only first order terms in $\displaystyle\beta$ without making any assumptions on the amplitude of the waves. The third equation in (1) (the energy equation) is a consequence of the mass and momentum equations (the first two equations). The fourth conservation law (generalized Bernoulli conservation law) can also be written here (cf. [16, 18]).

Mathematical justification of this model and some related systems can be found in [34, 31, 41, 29, 8]. A variational formulation of the SGN equations is given in [38, 30, 16]. The linear stability of solitary waves of small amplitude to the SGN equations was established in [30]. Also, it has been mentioned there that numerically, the solitary waves are stable for any wave amplitude. Recent years have seen increased activity in both the study of qualitative properties of the solutions to the SGN system and in the development of numerical discretization techniques [17, 32, 37, 13, 7, 35, 20].

In [10] the Riemann problem for the SGN equations was examined. Earlier, the Riemann problem was mainly studied for integrable systems in [23] for the KdV equation, and in [24, 9] for the NLS equation. Recently, this problem has received much attention for non-integrable systems of equations mainly because of the dispersive shocks commonly present in physics [11]. In [10] the wave number, amplitude, average fluid depth and average velocity have been chosen as primary variables to study the dispersive shocks of the SGN equations. This choice is quite natural, because, for example, the leading edge of the dispersive shock corresponds to the limit of small wave numbers, while the trailing edge is the limit of small amplitudes. So, the asymptotic study in the limit of small wave numbers or small amplitudes is important to predict the solution behaviour.

To capture better the case of moderate and large amplitude waves one can try to use other variables. Even if a priori they may be not necessarily physically tractable, they could be useful to parametrize globally the generic solution to the modulation equations. In particular, it could help to determine the regions of modulational stability and instability for waves of arbitrary amplitude.

The structure of the article is as follows. The system of four modulation equations is derived in Sections $\displaystyle 2$, $\displaystyle 3$, $\displaystyle 4$. In Section $\displaystyle 5$ the averaged quantities are expressed as functions of the roots of the third order polynomial determining the fluid depth behaviour, and the phase velocity. The non-conservative form of the modulation equations and their hyperbolicity analysis are given in Sections $\displaystyle 6$, $\displaystyle 7$. Numerical tests showing the wavetrain stability for the full SGN equations are presented in Section $\displaystyle 8$. Technical details are described in Appendix.

2 Averaging of the conservation laws of the SGN equations

A formal derivation of the modulation equations to the SGN equations (even in a more general formulation which contained, in particular, equations of bubbly fluids) can be found in [15, 20]. However, the analysis of the hyperbolicity for such a general formulation was not performed there. Here we will concentrate on SGN equations and will use the approach based on the averaging of conservation laws.

Suppose that the unknowns $\displaystyle h$, $\displaystyle u$ (and also $\displaystyle p$ and $\displaystyle e$ which are functions of these variables and their derivatives) depend on the rapid travelling coordinate $\displaystyle\xi=x-Dt$ and slow variables $\displaystyle X=\varepsilon x$, $\displaystyle T=\varepsilon t$ (see Appendix A for the details). Here $\displaystyle D$ is the travelling wave velocity. Let us introduce the following $\displaystyle\varepsilon$–expansion ansatz for $\displaystyle h$, $\displaystyle u$, $\displaystyle p$ and $\displaystyle e$:

[TABLE]

Here all the terms are supposed to be $\displaystyle L$-periodic with respect to $\displaystyle\xi$, where $\displaystyle L$ is also a slowly varying function of $\displaystyle X,T$. The substitution of these expansions into (1) yields:

[TABLE]

Here we took into account the following transformations of the partial derivatives with respect to time and space:

[TABLE]

In particular, the material derivative $\displaystyle\frac{D}{Dt}$ reads:

[TABLE]

Let us consider only the zero and first order approximation with respect to $\displaystyle\varepsilon$. The zero-order system reads:

[TABLE]

The first-order system reads:

[TABLE]

Since all the functions $\displaystyle h_{i}$, $\displaystyle u_{i}$, $\displaystyle e_{i}$, and $\displaystyle p_{i}$, $\displaystyle i=0,1$, are $\displaystyle L$-periodic with respect to $\displaystyle\xi$, after averaging the first order equations over the period $\displaystyle L$ one gets the following system:

[TABLE]

Notice that here we used the fact that the averaging procedure and the derivation with respect to slow variables commute (cf. [44, 5, 28]).

3 Stationary periodic solution

We will now show that the equations of zero order approximation (4) admit periodic solutions. We rewrite (4) in the form:

[TABLE]

Here and further the zero index is omitted and primes stand for $\displaystyle\displaystyle\frac{\partial}{\partial\xi}$. The first equation reads:

[TABLE]

The integration gives:

[TABLE]

When we write here and further, $\displaystyle m=const$, for example, we mean that $\displaystyle m$ does not depend on rapid variable $\displaystyle\xi$ : it is a function of just two variables $\displaystyle T$ and $\displaystyle X$. The second equation can be integrated as :

[TABLE]

The second-derivative term in (2) can be transformed. Keeping only the zero powers of $\displaystyle\varepsilon$, we rewrite the second derivative of $\displaystyle h$ as:

[TABLE]

Thus, we have

[TABLE]

Hence, the pressure expression in (2) reads:

[TABLE]

Replacing the pressure expression into (7) one obtains:

[TABLE]

Multiplying both side of the equation by $\displaystyle h^{{}^{\prime}}/m^{2}h^{2}$, we have

[TABLE]

which can be rewritten in the form:

[TABLE]

Integrating the equation once leads to

[TABLE]

Thus, the equation for $\displaystyle h$ is given as:

[TABLE]

Denote the polynomial in the right-hand side as $\displaystyle F_{3}(h)$ :

[TABLE]

or, equivalently:

[TABLE]

where $\displaystyle h_{0}<h_{1}<h_{2}$ are the roots of $\displaystyle F_{3}(h)$. Using Vieta’s formulas one can write $\displaystyle F_{3}(h)$ in the following way:

[TABLE]

where

[TABLE]

Identifying the coefficients of $\displaystyle F_{3}(h)$, one obtains :

[TABLE]

We are searching for the periodic solutions of (8), oscillating between two real positive roots $\displaystyle h_{1}$ and $\displaystyle h_{2}$. Since they are real, the third root $\displaystyle h_{0}$ is real too. Moreover, the last formula in (10) implies that $\displaystyle h_{0}$ is necessarily positive.

The periodic solution that oscillates between $\displaystyle h_{1}$ and $\displaystyle h_{2}$ is given by the formula:

[TABLE]

Here the Jacobi elliptic function $\displaystyle cn(u;k)$ is defined as:

[TABLE]

where $\displaystyle\varphi(u,k)$ is obtained implicitly from the relation

[TABLE]

The wavelength $\displaystyle L$ can explicitly be given as:

[TABLE]

where the interval $\displaystyle[\xi_{1},\xi_{2}]$ has the length $\displaystyle L$ and

[TABLE]

The wavelength is thus completely defined by the roots $\displaystyle h_{0}$, $\displaystyle h_{1}$ and $\displaystyle h_{2}$. The averaging of any arbitrary function of $\displaystyle f(h)$ reads:

[TABLE]

4 Averaged equations

Consider the first-order part of (5) (“zero” index is omitted):

[TABLE]

In the following, we will express all averaged quantities in (14) in terms of four unknowns: $\displaystyle h_{0}$, $\displaystyle h_{1}$, $\displaystyle h_{2}$, and $\displaystyle D$.

The flux in the first equation of (14) is:

[TABLE]

We introduced here the depth averaged velocity $\displaystyle U$. In terms of this velocity the mass equation can be rewritten in standard form:

[TABLE]

Since

[TABLE]

and

[TABLE]

the flux in the second equation of (14) is :

[TABLE]

The last two terms represent a combination of the pressure (defined up to multiplicative constant which is the fluid density) first integrated over the water depth and then averaged over the wave period, and the corresponding quadratic velocity correlation. Together, the two terms form an “effective pressure”. One can prove that such an effective pressure is always positive:

[TABLE]

In the third equation, we need to calculate $\displaystyle\overline{he}$, $\displaystyle\overline{hue}$, and $\displaystyle\overline{pu}$. Let us remark that for the travelling wave solutions one has

[TABLE]

Also, it follows from (8) that

[TABLE]

and

[TABLE]

The averaged energy is:

[TABLE]

Thus, the averaged energy reads:

[TABLE]

It can be also expressed as a function of $\displaystyle U$ (instead of $\displaystyle D$) and $\displaystyle h_{0}$, $\displaystyle h_{1}$, $\displaystyle h_{2}$. The volume average energy is:

[TABLE]

Thus,

[TABLE]

and

[TABLE]

Alos, one has :

[TABLE]

Now we are able to write the modulation equations because all quantities we need are given explicitly:

[TABLE]

The last step would be to replace the integration constants $\displaystyle m$, $\displaystyle i$ and $\displaystyle\epsilon$ by their expressions in terms of invariants $\displaystyle I_{i}$, $\displaystyle i=1,2,3$, using (10):

[TABLE]

The negative (positive) sign of $\displaystyle m$ corresponds to the right (left) facing periodic waves. We introduce the synthesis of both notations in order to obtain the simplest form of the equations. Basically, we will describe everything in terms of $\displaystyle m$, $\displaystyle i$ and $\displaystyle I_{1}$. Using (17), one can eliminate the dependence on $\displaystyle\epsilon$ in (16):

[TABLE]

Complemented by equation (28) for the wavelength (see Appendix A), equations (5) will finally be written as:

[TABLE]

We need now to rewrite (19) in quasilinear form in variables $\displaystyle D$, $\displaystyle h_{0}$, $\displaystyle h_{1}$, and $\displaystyle h_{2}$.

5 Expressions for the main averaged variables

The expressions of $\displaystyle\overline{h}$, $\displaystyle\overline{h^{-1}}$ and $\displaystyle L$ in terms of $\displaystyle h_{0}$, $\displaystyle h_{1}$, and $\displaystyle h_{2}$ are (for proof see Appendix B):

[TABLE]

Then, one can write the following differentials :

[TABLE]

Here $\displaystyle\Phi^{i}$, $\displaystyle\Psi^{i}$ and $\displaystyle\Lambda^{i}$ ($\displaystyle i=0,1,2$) read (for proof see Appendix C):

[TABLE]

The formulas for $\displaystyle\Phi^{k}$, $\displaystyle\Psi^{k}$ and $\displaystyle\Lambda^{k}$, $\displaystyle k=0,1,2$ were verified by hand calculations and with Wolfram Mathematica. One must pay attention to the fact that the complete elliptic integrals we use depend on elliptic modulus $\displaystyle k$, while Wolfram Mathematica uses the definition from Abramovitz and Stegun [1] where the complete elliptic integrals depend on parameter $\displaystyle m=k^{2}$ (do not confound the notations $\displaystyle m$ with $\displaystyle m$ coming from the mass conservation equation).

6 Nonconservative modulation equations

Complemented by (26) (or, equivalently, by equation (28)) (see Appendix A) the modulation equations (19) can be written in the following developed form:

[TABLE]

Or, in matrix form:

[TABLE]

where

[TABLE]

The coefficients of $\displaystyle A$ are given by:

[TABLE]

The coefficients of $\displaystyle B$ are given by:

[TABLE]

The characteristic eigenvalues are the roots of the fourth order polynomial

[TABLE]

To simplify the computations let us remark that the introduction of the depth averaged velocity $\displaystyle U=\overline{uh}/\overline{h}$ allows us to rewrite the mass conservation equation in standard form (15). Then we can use the fact that the equations are Galilean invariant if the depth average velocity is used. So, to check the hyperbolicity for any $\displaystyle U$ is equivalent to check the hyperbolicity for $\displaystyle U=0$. Since

[TABLE]

we will put into the coefficients of the matrices $\displaystyle A$ and $\displaystyle B$ the value of $\displaystyle D$ corresponding to $\displaystyle U=0$:

[TABLE]

Let us also remark that for the hyperbolicity study, one can always take $\displaystyle h_{0}=1$ in the coefficients of the polynomial (21) (if $\displaystyle h_{0}=\alpha>0$, all the eigenvalues will only be multiplied by $\displaystyle\sqrt{\alpha}$). The corresponding symmetry relations are the consequences of the fact that the periodic solution is determined in terms of the third degree polynomial.

7 Hyperbolicity region

Since the roots of polynomial (12) satisfy the inequality : $\displaystyle 1=h_{0}<h_{1}<h_{2}$, one can parametrize $\displaystyle h_{1}$ and $\displaystyle h_{2}$ as : $\displaystyle h_{1}=s$, $\displaystyle h_{2}=s+\tau$, where $\displaystyle s>1$ and $\displaystyle\tau>0$. Such a parametrization allows us to use a standard “Cartesian” frame for the computation of the eigenvalues of (21). The eigenvalues thus are given explicitly as functions of $\displaystyle s$ and $\displaystyle\tau$. We used Wolfram Mathematica for such a computation. The numerical results show that the eigenvalues are all real in a large region $\displaystyle\Omega=\{(s,\tau)|1<s<100,\;0<\tau<100\}$. First, for each pair $\displaystyle(s,\tau)$ from $\displaystyle\Omega$ we computed the numerical values of matrices $\displaystyle A$ and $\displaystyle B$. Then the corresponding eigenvalues were computed as the roots of the fourth order polynomial (21). Moreover, one can find that the resultant of the corresponding polynomial (21) and its derivative (denoted further as $\displaystyle R$) does not change sign (see Fig. 1). Thus, the polynomial (21) has no multiple real roots. Since all the roots are real and different, the system of modulation equations for the SGN equations is thus strictly hyperbolic. The fact that periodic waves of all lengths are modulationally stable corroborates the results of [30] where the spectral stability of solitary waves (the limit $\displaystyle s\rightarrow 1$) has been proven for small amplitudes and numerically confirmed for large amplitudes.

The whole hyperbolicity region $\displaystyle s>1,\;\tau>0$ is divided by a smooth curve corresponding to $\displaystyle\lambda=0$ into two sub-regions, “grey” and “white” (see Fig. 2). If $\displaystyle m<0$, in the grey sub-region one has three positive and one negative eigenvalues, while in the white sub-region one has two positive and two negative eigenvalues. If $\displaystyle m>0$, the signs of the roots will only change, since we have taken $\displaystyle U$ vanishing.

8 Numerical results

To study the modulational stability of periodic waves to SGN model numerically, we use the following set of parameters: $\displaystyle h_{0}=1\;m$, $\displaystyle h_{1}=1.5\;m$, $\displaystyle h_{2}=2\;m$ (i.e. $\displaystyle s=1.5\;m$ and $\displaystyle\tau=0.5\;m$), and $\displaystyle g=10\;m/s^{2}$, as an example, for the solution of the water height described by (11) over a single wave length $\displaystyle L$ which corresponds to $\displaystyle 7.4163\;m$ approximately. To setup the problem, a wave train is formed initially that consists of $\displaystyle N$ aforementioned single stationary wave solutions (see Section 3) in a domain of size $\displaystyle L_{1}=N\times L$, where the travelling wave speed of this wave train is taken be $\displaystyle D=-m/\overline{h}\approx 3.1688\;m/s$ (this allows us to take vanishing the averaged over period the mass weighted velocity $\displaystyle U=\overline{hu}/\overline{h}$). With that, we then introduce perturbations to the height of the wave train $\displaystyle h(x)$ as

[TABLE]

where $\displaystyle a$ is a small parameter. For this problem, periodic boundary conditions are assumed and used on the left and right of the interval $\displaystyle[0,L_{1}]$.

In the numerical simulations of the SGN model performed below, we take $\displaystyle N=50$, and perturbation amplitudes $\displaystyle a=10^{-j}$ for $\displaystyle j=1,2,3$ in the runs. We show the initial condition of the test in Fig. 3, and the computed solutions at four different times $\displaystyle t=200$, $\displaystyle 400$, $\displaystyle 800$, $\displaystyle 1200$ $\displaystyle s$ in Fig. 4, where both the water height and the phase portrait in the $\displaystyle(h,h\dot{h})$-plane are present. The choice of $\displaystyle h\dot{h}$ variable is natural, because on travelling wave solutions $\displaystyle\displaystyle{h\dot{h}=m\frac{dh}{d\xi}}$, so up to a multiplicative constant, $\displaystyle(h,h\dot{h})$-plane is nothing than the classical phase space $\displaystyle\left(h,\frac{dh}{d\xi}\right)$. The periodic wave remains stable when the smaller values $\displaystyle a=10^{-2}$ and $\displaystyle 10^{-3}$ are used. To see the limit of linear stability, we have also taken a large amplitude perturbation ($\displaystyle a=10^{-1}$). The periodic wave train becomes unstable : we are too far from the classical “small perturbation analysis”. The numerical results are obtained using a hyperbolic-elliptic splitting method proposed by the authors [20] with $\displaystyle 400$ meshes for each wave length.

9 Conclusions and perspectives

We have derived the modulation equations system to the SGN model and show that it is strictly hyperbolic for arbitrary wave amplitudes, i.e. the periodic wave trains are modulationally stable. This corroborates the results [30] where the linear stability of solitary waves (which can be considered as the limit of periodic waves of large length) has been proven. The existence of the Riemann invariants and nature of the characteristic fields (genuinely degenerate or genuinely nonlinear in the sense of Lax) will be the topic for future research.

A Multiscale decomposition

The classical Whitham method [44] consists in decomposing the scales in the following way (for simplicity, we consider just the velocity variable $\displaystyle u$):

[TABLE]

Here $\displaystyle\theta$ is a fast phase variable, $\displaystyle\Theta$ is a slow phase variable, and $\displaystyle\varepsilon$ is a small parameter. The solution is supposed to be $\displaystyle 2\pi$-periodic with respect to $\displaystyle\theta$. The definitions of the local wave frequency $\displaystyle\omega$ and the local wave number $\displaystyle\kappa$ :

[TABLE]

automatically imply the evolution equation for $\displaystyle\kappa$:

[TABLE]

Written in slow variables $\displaystyle X=\varepsilon x$, $\displaystyle T=\varepsilon t$, equations (24) are equivalent to

[TABLE]

and (25) reads as

[TABLE]

One can also define the travelling wave coordinate $\displaystyle\xi=x-Dt$, and the phase velocity $\displaystyle D=\omega/\kappa$. The solution is decomposed as:

[TABLE]

The wavelength $\displaystyle L$ is defined as:

[TABLE]

Thus, since $\displaystyle\theta\in[0,2\pi]$, then $\displaystyle\xi=\frac{\Theta}{\varepsilon\kappa}\in\left[0,\frac{2\pi}{\kappa}\right]=[0,L]$. So, $\displaystyle u$ is an $\displaystyle L$-periodic function with respect to the travelling coordinate $\displaystyle\xi$:

[TABLE]

Finally, (26) reads:

[TABLE]

The approach employing the travelling coordinate is equivalent to the one using the phase variable. The consistency equation can always be written in any of the two forms presented above: (26) or (28).

B Computation of elliptic integrals

Let

[TABLE]

Then one has :

[TABLE]

Here $\displaystyle K(k)$, $\displaystyle E(k)$ and $\displaystyle\Pi(n,k)$ are the complete elliptic integrals of the first, second and third type, respectively :

[TABLE]

The characteristic $\displaystyle n$ and elliptic modulus $\displaystyle k$ are defined as :

[TABLE]

The definition of averaging (13) implies then :

[TABLE]

The derivatives of the complete elliptic integrals read as follows:

[TABLE]

Or, in terms of $\displaystyle h_{i}$:

[TABLE]

C Differentials of $\displaystyle\overline{h}$, $\displaystyle\overline{h^{-1}}$ and $\displaystyle L$

The differentials of $\displaystyle n$ and $\displaystyle k$ are given by:

[TABLE]

Then, the differentials of $\displaystyle\overline{h}$ and $\displaystyle\overline{h^{-1}}$ are :

[TABLE]

The differential of $\displaystyle L$ is :

[TABLE]

Acknowledgments Part of this work was performed during the visit of SG to the Department of Mathematics of National Taiwan University. SG was partially supported by l’Agence Nationale de la Recherche, France (grant numbers ANR-11-LABEX-0092, and ANR-11-IDEX-0001-02). KMS was partially supported by MOST 105-2115-M-002-013-MY2 and 107-2115-M-002-014.