Vincent Duch

Christian Klein

NUMERICAL STUDY OF THE SERRE-GREEN-NAGHDI EQUATIONS AND A FULLY DISPERSIVE COUNTERPART

Keywords: November 16, 2021. 2010 Mathematics Subject Classification. Primary: 65M70, 35Q35, 76B15; Secondary: 35B35 Nonlinear dispersive equations, solitary waves, modulated oscillations, numerical investigation, spectral methods

We perform numerical experiments on the Serre-Green-Naghdi (SGN) equations and a fully dispersive "Whitham-Green-Naghdi" (WGN) counterpart in dimension 1. In particular, solitary wave solutions of the WGN equations are constructed and their stability, along with the explicit ones of the SGN equations, is studied. Additionally, the emergence of modulated oscillations and the possibility of a blow-up of solutions in various situations is investigated.

We argue that a simple numerical scheme based on a Fourier spectral method combined with the Krylov subspace iterative technique GMRES to address the elliptic problem and a fourth order explicit Runge-Kutta scheme in time allows to address efficiently even computationally challenging problems.

1. Introduction 1.1. Motivation. The Serre-Green-Naghdi (SGN) model is a popular model for the propagation of surface gravity waves in coastal oceanography. It is expected to provide a reasonable approximation of the response to gravity forces of a layer of homogeneous incompressible fluid with a free surface (hereafter referred to as the water waves problem) in the so-called shallow-water regime, that is for weakly dispersive but possibly strongly nonlinear flows. It has been derived and studied by many authors, including [START_REF] Barthélemy | Nonlinear shallow water theories for coastal waves[END_REF][START_REF] Bazdenkov | Dispersive effects in two-dimensional hydrodynamics[END_REF][START_REF] Camassa | Long-time effects of bottom topography in shallow water[END_REF][START_REF] Clamond | Practical use of variational principles for modeling water waves[END_REF][START_REF] Fedotova | Energy equation for certain approximate models of long-wave hydrodynamics[END_REF][START_REF] Gavrilyuk | A kinematic conservation law in free surface flow[END_REF][START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF][START_REF] Ionescu-Kruse | Variational derivation of the Green-Naghdi shallow-water equations[END_REF][START_REF] Kim | A derivation of the Green-Naghdi equations for irrotational flows[END_REF][START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF][START_REF] Miles | Weakly dispersive nonlinear gravity waves[END_REF][START_REF] Salmon | Hamiltonian fluid mechanics[END_REF][START_REF] Seabra-Santos | Numerical and experimental study of the transformation of a solitary wave over a shelf or isolated obstacle[END_REF][START_REF] Serre | Contribution à l'étude des écoulements permanents et variables dans les canaux[END_REF][START_REF] Su | Korteweg-de Vries equation and generalizations. III. Derivation of the Korteweg-de Vries equation and Burgers equation[END_REF][START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. I. Highly nonlinear unsteady waves[END_REF]. Its rigorous justification as an asymptotic model in the shallow-water limit has been obtained in [START_REF] Alvarez-Samaniego | A Nash-Moser theorem for singular evolution equations. Application to the Serre and Green-Naghdi equations[END_REF][START_REF] Fujiwara | A shallow water approximation for water waves over a moving bottom[END_REF][START_REF] Israwi | Large time existence for 1D Green-Naghdi equations[END_REF][START_REF] Li | A shallow-water approximation to the full water wave problem[END_REF][START_REF] Makarenko | A second long-wave approximation in the cauchy-poisson problem[END_REF]. In addition to its validity as a model for the water waves problem, the SGN equations have attracted interest as they are natural dispersive generalizations of the equations for isentropic compressible flows and as such can be studied through the Lagrange formalism [START_REF] Gavrilyuk | Multiphase flow modeling via Hamilton's principle[END_REF][START_REF] Gavrilyuk | Stationary shock-like transition fronts in dispersive systems[END_REF][START_REF] Gavrilyuk | Generalized vorticity for bubbly liquid and dispersive shallow water equations[END_REF]. In addition, in the irrotational framework they can be obtained through canonical Hamilton's equations [START_REF] Duchêne | Well-posedness of the Green-Naghdi and Boussinesq-Peregrine systems[END_REF][START_REF] Salmon | Hamiltonian fluid mechanics[END_REF], consistently with Zakharov's Hamiltonian formulation of the water waves problem [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF]. Shortly put, the SGN system enjoys strong structural properties.

In this work, we numerically compare the SGN equations with a model introduced by the first author and collaborators in [START_REF] Duchêne | A new class of two-layer Green-Naghdi systems with improved frequency dispersion[END_REF]. The model is obtained using Hamilton's equations with a modified Hamiltonian, and hence preserves at least part of the structure of the SGN equations, while having the additional property that the dispersion relation of the system coincides exactly with the one of the water waves problem. Models with such properties are often said to be fully dispersive, and have been advocated by G. B. Whitham [START_REF] Whitham | Variational methods and applications to water waves[END_REF] as a way to reproduce -at least qualitatively-in a better way some key properties of the water waves problem, such as wavebreaking and non-smooth travelling waves of extreme height. The price to pay is that the equations include non-local pseudodifferential operators (Fourier multipliers). Whitham's prediction turned out to be valid at least for the unidirectional model which bears his name, as shown by [START_REF] Ehrnström | On Whitham's conjecture of a highest cusped wave for a nonlocal dispersive equation[END_REF][START_REF] Hur | Wave breaking in the Whitham equation[END_REF][START_REF] Saut | The wave breaking for Whitham-type equations revisited[END_REF][START_REF] Truong | Global bifurcation of solitary waves for the Whitham equation[END_REF]. This fact triggered renewed activity on bidirectional models (systems), and we refer to the surveys [START_REF] Carter | Bidirectional Whitham equations as models of waves on shallow water[END_REF][START_REF] Dinvay | A comparative study of bi-directional Whitham systems[END_REF][START_REF] Klein | On Whitham and related equations[END_REF] for more information. The aforementioned model refines systems studied therein (sometimes called Whitham-Boussinesq systems) so as to offer improved precision in strongly nonlinear situations. We refer to it in this work as the Whitham-Green-Naghdi (WGN) model. It has been rigorously justified among other fully dispersive models in [START_REF] Duchêne | Many models for water waves[END_REF][START_REF] Emerald | Rigorous derivation from the water waves equations of some full dispersion shallow water models[END_REF]. 1 In this work, we numerically investigate properties of the SGN and WGN equations in extreme situations. More precisely, we will investigate features and stability of solitary waves with large height and large velocity, as well as solutions whose evolution produces steep gradients. It must be emphasized that both the SGN and WGN systems are expected to provide poor approximations to the water waves problem in the above scenarii since we voluntarily depart from the shallow-water regime of validity (δ 1). Our motivation is theoretical as we aim at extracting information on the role of dispersive properties for such fully nonlinear models. We choose the SGN and WGN equations as our subject of investigations in order to step out of the world of unidirectional scalar (nonlinear and dispersive) equations for which similar studies have been realized [START_REF] Abenda | Numerical solution of the small dispersion limit of the Camassa-Holm and Whitham equations and multiscale expansions[END_REF][START_REF] Grava | Numerical study of a multiscale expansion of Korteweg-de Vries and Camassa-Holm equation[END_REF][START_REF] Grava | A numerical study of the small dispersion limit of the Korteweg-de Vries equation and asymptotic solutions[END_REF][START_REF] Klein | Fourth order time-stepping for low dispersion Korteweg-de Vries and nonlinear Schrödinger equations[END_REF][START_REF] Klein | A numerical approach to blow-up issues for dispersive perturbations of Burgers' equation[END_REF], while retaining strong structural properties. In particular, solitary waves can be identified with critical points of functionals which directly derive from the aforementioned Hamiltonian structure [START_REF] Duchêne | Solitary Wave Solutions to a Class of Modified Green-Naghdi Systems[END_REF]. Moreover, the two systems of equations can (and will) be numerically approximated using identical numerical strategies, specifically Fourier pseudospectral methods.

1.2. The equations. In the (horizontal) one-dimensional setting and flat-bottom framework, the SGN equations read (see e.g. [START_REF] Duchêne | Well-posedness of the Green-Naghdi and Boussinesq-Peregrine systems[END_REF])

(SGN)    ∂ t ζ + ∂ x (hu) = 0, ∂ t u -1 3h ∂ x (h 3 ∂ x u) + g∂ x ζ + u∂ x u = ∂ x u 3h ∂ x (h 3 ∂ x u) + 1 2 h 2 (∂ x u) 2 .
1 More precisely, it is proved to be an approximate model to the water waves system of order O( δ 4 ) in the sense of consistency, where is the "nonlinearity" dimensionless parameter defined as the ratio of the typical amplitude of the wave to the reference layer depth, and δ is the "shallowness" dimensionless parameter defined as the ratio of the reference layer depth to the typical horizontal wavelength. The corresponding precision of the SGN equations is O(δ 4 ), and the one of Whitham-Boussinesq equations is O( δ 2 ). The improvement between the WGN predictions and SGN predictions can be witnessed in [26, Figures 3 and 4] in the context of small-amplitude solitary waves and in [27, §I.5] for time-evolving profiles.

Here, d > 0 is the reference layer depth, g > 0 is the gravitation constant,2 and u(t, x) represents the layer-averaged horizontal velocity, ζ(t, x) (or rather its graph) represents the surface deformation and h(t, x) = d + ζ(t, x) represents the water depth at time t and horizontal position x ∈ R. We refer to [START_REF] Duchêne | Well-posedness of the Green-Naghdi and Boussinesq-Peregrine systems[END_REF] for a description of its canonical Hamiltonian structure. Known conserved quantities of (SGN) are

R f i dx (i = 1, . . . , 4) with densities (1) f 1 = ζ, f 2 = hu, f 3 = gζ 2 + hu 2 + 1 3 h 3 (∂ x u) 2 , f 4 = u - 1 3h ∂ x (h 3 ∂ x u)
representing respectively the mass, momentum (or horizontal impulse), total energy and the rescaled tangential fluid velocity at the free interface [START_REF] Gavrilyuk | A kinematic conservation law in free surface flow[END_REF].

The fully dispersive model introduced in [START_REF] Duchêne | A new class of two-layer Green-Naghdi systems with improved frequency dispersion[END_REF] (when restricted to the one-layer case and neglecting surface tension) is 2 , where F is the Fourier multiplier defined by

(WGN)    ∂ t ζ + ∂ x (hu) = 0, ∂ t u -1 3h ∂ x F(h 3 ∂ x Fu) + g∂ x ζ + u∂ x u = ∂ x u 3h ∂ x F(h 3 ∂ x Fu) + 1 2 h 2 (∂ x Fu)
∀ϕ ∈ L 2 (R), Fϕ(ξ) = F (d|ξ|) ϕ(ξ)
where

F (k) = 3 |k| tanh(|k|) - 3 |k| 2 .
Known conserved quantities of (WGN) are R f i dx (i = 1, . . . , 4) with densities

(2)

f 1 = ζ, f 2 = hu, f 3 = gζ 2 +hu 2 + 1 3 h 3 (∂ x Fu) 2 , f 4 = u- 1 3h ∂ x F(h 3 ∂ x Fu).
Our convention for the Fourier transform is the one for which the following identities hold for sufficiently regular and localized functions g:

∀k ∈ R, g(k) := 1 (2π) 1/2 R e -ikx g(x) dx, ∀x ∈ R, g(x) = 1 (2π) 1/2 R e +ikx g(k) dk.
1.3. Outline. Let us now present the structure of this manuscript. In section 2 we numerically construct solitary wave solutions to the WGN equations. In section 3 we present our numerical approach for the numerical solution of the initial-value problem for the SGN and WGN equations, and test its validity on solitary wave solutions. The stability of solitary waves for both SGN and WGN equations is numerically investigated in section 4. In section 5 we study the emergence of zones of rapid modulated oscillations within solutions to both equations starting from unidirectional, long, smooth and localized initial data. Based on the relationship between the SGN equations and the Camassa-Holm equation, we study in section 6 the behaviour of solutions to the former with initial data leading to finite-time blow-up for the latter. In section 7, we study the possibility of a blow-up for initial data near cavitation, that is vanishing depth. We summarize our findings and add some concluding remarks in section 8.

Solitary waves

In this section we study solitary wave solutions to the (fully dispersive) Whitham-Green-Naghdi equations (WGN), that is solutions of the form

(3) ζ(t, x) = ζ c (x -ct) u(t, x) = u c (x -ct), lim |x|→∞ |ζ c |(x) + |u c |(x) = 0
where the constant c ∈ R is the solitary wave velocity.

It is well-known that for any supercritical velocity c > 1 (recall we set g = d = 1), there exists a smooth solitary wave solution to (SGN) with explicit formula 3 (4)

ζ c (x) = (c 2 -1) sech 2 √ 3 2 c 2 -1 c 2 x , u c (x) = cζ c (x) 1 + ζ c (x)
.

The functions ζ c and u c are smooth, even and positive and have a unique nondegenerate maximum at the origin. Such an explicit formula is of course unexpected for the fully dispersive system (WGN). However, the following result has been shown in [START_REF] Duchêne | Solitary Wave Solutions to a Class of Modified Green-Naghdi Systems[END_REF]:

Proposition 2.1. There exists (ζ (q) , u (q) ) q>0 a one-parameter family of smooth square-integrable functions such that for all q > 0, (ζ cq , u cq ) := (ζ (q) , u (q) ) provides a solitary wave solution to (WGN) with velocity c q > 1, and

c q → 1 and (c 2 q -1) -1 ζ (q) ((c 2 q -1) -1/2 •) -sech 2 ( √ 3 2 •) H 1 → 0 (q → 0)
. We also refer to [START_REF] Duchêne | Solitary Wave Solutions to a Class of Modified Green-Naghdi Systems[END_REF] for numerical computations of WGN solitary waves with small velocities 0 < c -1 1. In the following we numerically investigate the existence and behavior of WGN solitary waves for large velocities. Based on these numerical experiments we conjecture the following. Conjecture 2.2. For all c > 1, there exist smooth and rapidly decaying solitary wave solutions to the Whitham-Green-Naghdi system (WGN) with velocity c and such that the following holds.

(1) For all c > 1, the profiles ζ c , u c are positive on the real line.

(2) For all c > 1, the profiles ζ c , u c have a unique critical point corresponding to their maximum, and it is non-degenerate.

(3) For all c > 1, the profiles ζ c , u c are symmetric about their maximum.

(4) max ζ c ∼ c 2 and max u c ∼ c as c → ∞.
This is in sharp contrast to the celebrated result [START_REF] Amick | On solitary water-waves of finite amplitude[END_REF] on the existence of (peaked) solitary waves of extreme height for the water waves problem, and the corresponding result obtained on the Whitham equation [START_REF] Ehrnström | On Whitham's conjecture of a highest cusped wave for a nonlocal dispersive equation[END_REF][START_REF] Truong | Global bifurcation of solitary waves for the Whitham equation[END_REF] (see also [START_REF] Klein | A numerical approach to blow-up issues for dispersive perturbations of Burgers' equation[END_REF] and references therein for a numerical investigation), and invalidates the naive thinking that this feature relies only on the dispersion relation of the equations linearized about the rest state.

2.1. The equations for solitary waves. Plugging (3) into (WGN) yields for the first equation ( 5)

ζ c = h c u c c = u c c -u c , h c = 1 + ζ c = c c -u c ,
and for the second

(6) u c -c 3h c ∂ x F(h 3 c ∂ x Fu c ) + 1 2 h 2 c (∂ x Fu c ) 2 + cu c -ζ c - u 2 c 2 = 0.
3 In this work we consider only smooth solitary waves maintaining positive depth; see e.g. [START_REF] Jiang | Classification of traveling wave solutions to the Green-Naghdi model[END_REF].

From (5) we infer ( 7)

u c = cζ c h c ,
and plugging (7) into (5) yields a single equation for ζ c , namely

(8) -1 3h 2 c ∂ x F(h 3 c ∂ x F ζc hc ) + 1 2 h 2 c (∂ x F ζc hc ) 2 + ζ c h c - ζ c c 2 - ζ 2 c 2h 2 c = 0.
Similarly, we may use [START_REF] Antonopoulos | On the well-posedness of the Galerkin semidiscretization of the periodic initial-value problem of the Serre equations[END_REF] with [START_REF] Amick | On solitary water-waves of finite amplitude[END_REF] to produce a single equation for u c , namely

(9) - (c -u c ) 2 3c ∂ x F ( c c-uc ) 3 ∂ x Fu c + 1 2 ( c c-uc ) 2 (∂ x Fu c ) 2 + cu c - u c c -u c -1 2 u 2 c = 0.
Finally, we find it convenient to solve (9) using the variable η c = uc c = ζc hc :

(10)

-(1 -η c ) 2 3 ∂ x F ∂ x Fη c (1 -η c ) 3 + 1 2(1 -η c ) 2 ∂ x Fη c 2 +η c - η c c 2 (1 -η c ) -1 2 η 2 c = 0.
Remark 2.3. Equation (8) can be written as δ ζc L = 0 with

L := R (ζ c , c ζc 1+ζc ) dx
where

(ζ, u) := 1 2 ζ 2 - 1 2 (1 + ζ)u 2 - 1 6 (1 + ζ) 3 (∂ x Fu) 2
is the Lagrangian density naturally associated with the Hamiltonian formulation of the equations, and physically corresponds to the difference between the potential energy and the kinetic energy. In particular, the Jacobian of δL is, by definition, symmetric for the L 2 -inner product.

2.2. Numerical construction of solitary waves. We seek numerical approximations of solitary waves for (WGN) with fixed velocity c > 1 through zeroes of a finite-dimensional vector-field accounting for the left-hand side of ( 9) or ( 10):

(11) F(u) = 0.

The vector u = (u(x 1 ), . . . , u(x n )) represents values of the function u at collocation points x n = -πL + n2πL/N , n = 1, . . . , N for x ∈ L[-π, π]. Nonlinear operations are naturally computed at collocation points, while ∂ x F is approximated via a discrete Fourier transform computed efficiently with a Fast Fourier transform (FFT) and multiplication in Fourier space, that is

∀v ∈ R N , ∂ x Fv k := ik k F (k k ) v k
where we denote v = ( v -N/2+1 , . . . , v N/2 ) the coefficients of the Fast Fourier transform of v (which we slightly incorrectly refer to as Fourier coefficients in the following), and k k = (k/L), k = -N/2 + 1, . . . , N/2 the discrete Fourier modes. Here L > 0 is a constant chosen such that u and its relevant derivatives decrease to machine precision (roughly 10 -16 in double precision) and N will be chosen such that the u k decreases to machine precision for large |k|. In the results presented here, multiplying N or LN by a factor 2 or 4 neither improves nor deteriorates the accuracy. For more information on Fourier spectral methods, the reader is referred to [START_REF] Trefethen | Spectral Methods in Matlab[END_REF] and the literature cited therein.

The system of N nonlinear equations [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF] will be solved iteratively by a standard Newton iteration, ( 12)

u (m+1) = u (m) -δu (m) , with (13) Jac(F(u))| u=u (m) δu (m) = F(u (m) ),
where u (m) denotes the mth iterate, and where JacF(u) is the Jacobian of F(u) with respect to u. The initial iterate u (0) will be chosen as the SGN solitary wave given by (4) at collocation points. It may also be constructed by extrapolating from previously computed solitary waves and Lagrange polynomials. For the iteration, we apply and compare two strategies: first we test a Newton-GMRES method, i.e., we solve [START_REF] Camassa | An integrable shallow water equation with peaked solitons[END_REF] with the Krylov subspace iterative method GMRES [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] as for instance in [START_REF] Klein | A numerical approach to blow-up issues for dispersive perturbations of Burgers' equation[END_REF]. Alternatively we solve (13) through standard LU factorization. Notice that due to the translation invariance of the problem, the Kernel of the Jacobian of the continuous (infinite-dimensional) vector-field is non-empty when evaluated at a non-trivial solution u, since ∂ x u is an element of its nullspace. The corresponding spectral projection can be inferred from the symmetry property mentioned in Remark 2.3. Correspondingly, Jac(F(u (m) )) has an extremely small eigenvalue as u (m) converges towards the desired solution. However, by symmetry considerations, we can ensure that at each iterate u (m) and therefore F(u (m) ) are even, and as a consequence its spectral projection onto the corresponding eigenspace vanishes (up to machine precision). Yet we find it advisable to add the aforementioned spectral projection to the Jacobian when solving [START_REF] Camassa | An integrable shallow water equation with peaked solitons[END_REF], although in practice we mostly observe a slight phase shift on the numerical approximation if the spectral projection is not added.

Let us now present our results, starting with the case of small and slow solitary waves. For c = 1.1, we treat the equation ( 9) for x ∈ 20[-π, π] and use N = 2 9 collocation points. We use a Krasny filter of the order of 10 -14 , which puts to zero Fourier coefficients with modulus smaller than 10 -14 . We also apply a preconditioner of the form M = Diag((1+k 2 /3)) -which is motivated by the linear dispersion of the SGN equations-i.e., instead of solving iteratively with GMRES Ax = b the equation [START_REF] Camassa | An integrable shallow water equation with peaked solitons[END_REF] in Fourier space, we solve M -1 Ax = M -1 b. The Newton-GMRES code converges within 3 iterations with a residual F(u (3) ) ∞ of the order of 10 -13 . The residual of the initial iterate, F(u (0) ) ∞ with u (0) the SGN solitary wave, is of the order of 10 -2 , GMRES converges in 21 iterations with a relative residual of 10 -11 . The resulting Newton iterate, F(u (1) ) ∞ has a residual of 10 -4 , GMRES again converges within 20 iterations to a relative residual of 10 -11 . The residual of the subsequent Newton iteration, F(u (2) ) ∞ is then of the order of 10 -7 . In the next step, GMRES stagnates with a relative residual of the order of 10 -7 , which is explained by the smallness of the previous residual of the Newton iteration. The ensuing residual of the Newton iteration, F(u (3) ) ∞ , is of the order of 10 -13 , and the iteration is stopped. The iteration thus shows the well-known quadratic convergence of a Newton iteration, loosely speaking the number of correct digits doubles in each iteration. The resulting solitary wave can be seen in Fig. 1 on the left. The WGN solitary wave for c = 1.1 is very close to the SGN solitary wave for the same velocity shown in red in the same figure. On the right of the same figure, the modulus of the Fourier coefficients shows that the solution is numerically resolved to the order of the Krasny filter. We also observe that the exponential decay rate of Fourier coefficients with large Fourier modes is slightly smaller for the WGN solitary wave than it is for the SGN solitary wave. The Newton-GMRES code converges in 4 iterations with a residual of the order of 10 -12 . The resulting solitary wave is shown in Fig. 2 on the left. Again the solution is very close to the SGN solitary wave shown in the same figure in red. The solution is well resolved in Fourier space as can be seen on the right of Fig. 2, and we observe once more that the exponential decay rate of Fourier modes is smaller. The Newton-GMRES iteration behaves similarly to what is described before. Note that it does not converge without a preconditioner due to issues for the high Fourier modes. For even larger values of the velocity such as c = 3, the iteration no longer converges because of GMRES problems for the high Fourier modes. Therefore we switch for larger values of the velocity c to a Newton iteration with a direct numerical factorization of the Jacobian. Still with N = 2 10 collocation points on x ∈ 10[-π, π], we consider the case c = 3. The Newton iteration converges with a direct inversion of the Jacobian in 3 iterations to a residual of the order of 10 -9 . After further iterations, the minimal residual reachable with this method appears to be of the order of 10 -13 . For c = 20 and the same parameters as before, iteration converges normally to the order of 10 -7 . The solution and its Fourier coefficients can be seen in Fig. 3. The deviation from the SGN solitary wave in red is now clearly visible. We again observe a lower rate of decay of Fourier coefficients of the WGN solitary wave, however only for higher Fourier modes (the lower decay rate is apparent starting from much smaller Fourier modes when plotting the corresponding figures for the surface deformation, ζ c ). The Fourier coefficients saturate at the order of 10 -10 indicating problems in the conditioning of the Jacobian which we investigate later on. Part of this can be attributed to the presumption motivated by formula ( 4 For even larger values of the velocity as c = 100, these problems become worse and the Newton iteration does not converge. There is however no indication that the solitary wave might become singular or might not exist for such values. Since u c is proportional to c, we look at the rescaled quantity η c = u c /c = ζ c /(1 + ζ c ) and consider the equation [START_REF] Bazdenkov | Dispersive effects in two-dimensional hydrodynamics[END_REF]. Using again the same parameters as before, the optimal residual for the iteration in this case is of the order of 10 -6 . The solution can be seen in Fig. 4 on the left. The Fourier coefficients on the right of the same figure saturate at the order of 10 -8 which partly explains why no lower residual can be achieved. The lower exponential rate of decay of the WGN solitary wave Fourier coefficients is no longer visible, presumably because it occurs for higher Fourier modes than the numerically well-resolved ones. It is apparent when plotting the corresponding figures for the surface deformation, ζ c .

Let us briefly comment on this decay rate. It is well-known that it relates to the maximal width of the strip around the real axis for which the analytic extension of the function is free of singularities (see e.g. [3, Theorem 7.1]). The authors have no clear understanding on how such properties can be related to the balancing of dispersion and nonlinearity of the corresponding equation. The exponential decay rate of the Fourier transform of the surface deformation of the SGN solitary waves can be explicitly inferred from formula (4). Such is not the case for the velocity variable, to the authors' knowledge, due to the fact that the function z → z/(1+z) is not entire (see e.g. [3, Proposition 7.10]), although numerical computations indicate that the decay rate is the same. The decay rate for the surface deformation and velocity variables of the WGN solitary waves also appear to be the same. They are close to the ones of the SGN solitary waves for small values of the velocity c ≈ 1, and quickly deviate for larger values (about c = 2) to approximately half the ones of the SGN solitary waves. This interesting question will certainly need to be studied more in the future. One reason for the problems in the iterations are due to the function ζ c being of the order c 2 . We show these functions for c = 20 and c = 100 in Fig. 5. It can be seen that this function is always more peaked than the corresponding one for the SGN equations with the same velocity, given by the explicit formula (4). To understand the difficulties in the Newton iterations, we look at the Jacobian for the initial iterate (the SGN solitary wave) for c = 100 which is shown in Fig. 6. We denote the discrete Fourier transform of this Jacobian (from the left and from the right) by Ĵ. It can be seen that it is, except for a peak near the center, essentially constant along the diagonal. In addition there is a plateau of the order of 10 -5 seemingly due to rounding errors since the maximum is of the order of 10 10 . Note that the dominant contribution to the Jacobian is due to

Ĵ * := -(1-ηc) 2 3 ∂ x F 1 (1-ηc) 3 ∂
x F which we refer to as the 'non diagonal' part.

Time evolution

In this section we present and validate the numerical scheme we employ for integrating in time the equations (SGN) and (WGN). In principle we adapt here a classic code for the KdV equation (code 27 in [START_REF] Trefethen | Spectral Methods in Matlab[END_REF]) to the equations studied here. The main new aspect is the inversion of the elliptic operator via GMRES, and we concentrate on this aspect whilst refering to [START_REF] Trefethen | Spectral Methods in Matlab[END_REF] for general properties of the numerical approach. Let us first recall the SGN and WGN equations with a slight reformulation: sufficiently regular solutions to (WGN) satisfy (once again, we set The presence of Fourier multipliers in the WGN equations naturally leads to Fourier pseudospectral methods, already employed in [START_REF] Dutykh | Finite volume and pseudo-spectral schemes for the fully nonlinear 1D Serre equations[END_REF] and described thereafter. One of the difficulties when integrating in time the SGN equations (and in a similar way the WGN equations) is that we are led to solve the elliptic problem ( 15) at each timestep. The aforementioned issue is addressed in particular in [START_REF] Dias | On the fully-nonlinear shallow-water generalized Serre equations[END_REF][START_REF] Duran | Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations[END_REF][START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF] (see also [START_REF] Favrie | A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves[END_REF] and references in [START_REF] Duchêne | Many models for water waves[END_REF] for relaxation approaches) where different approximate models are introduced. In this work, we stick with the original equations and simply note that, thanks to the efficiency of pseudospectral methods, it is not too costly -at least in our one-dimensional framework-to solve the elliptic problem at each time step while maintaining high resolution. To this aim, we observe in our experiments that the aforementioned Krylov subspace iterative method GMRES is highly efficient and converges within a few iterations to the desired accuracy, although the choice of preconditioner may turn out to be crucial. Only for extreme situations not studied here and far from the range of applicability of the equations an inversion via standard LU factorization is found necessary.

g = d = 1) (14)    ∂ t ζ + ∂ x (hu) = 0, ∂ t v + ∂ x ζ + uv -1 2 u 2 -1 2 h 2 (∂ x Fu) 2 =
Let us now be more precise. We use the same Fourier pseudospectral approach as outlined in the previous section, i.e., we approximate the solution u, ζ via discrete Fourier transforms. This means we can treat initial data which are smooth and periodic or in the Schwartz class of rapidly decreasing functions (the latter can be treated within the finite numerical precision as periodic on sufficiently large domains) with spectral accuracy, i.e., with a numerical error exponentially decaying with the number N of Fourier modes.

With this spatial discretization, both SGN and WGN are finite dimensional systems of ODEs coupled with a system of equations of the form ( 16)

                 d ζ dt = G 1 ( ζ, û), dv dt = G 2 ( ζ, û, v), M( ζ)û = v
, where ζ(t), û(t), v(t) are N -dimensional vectors, and M( ζ) is an N -by-N matrix. The two ODEs in system ( 16) will be integrated with the standard explicit fourth order Runge-Kutta method (RK4). Note that this is not trivial since the system will be stiff because of the three derivatives with respect to x on the right hand side of (SGN). The stiffness implies that explicit schemes as the ones applied here can become inefficient because of restrictive stability conditions on the time step. We briefly recall the basic concepts of Dahlquist's stability theory, for more details see Chapter 10 of [START_REF] Trefethen | Spectral Methods in Matlab[END_REF] and references therein. The basic idea is to consider a linear model problem y = λy, λ ∈ C, where λ is some characteristic parameter.

If λ ≤ 0, the exact solution of the model problem is bounded for positive times, and a numerical solution is called stable if it also bounded for positive times. For explicit time evolution schemes as the one used here, this condition defines a domain of stability in the complex z := λh plane (with h the time step), see output 25 in Chapter 10 of [START_REF] Trefethen | Spectral Methods in Matlab[END_REF] for RK4. Considering now an ODE system y = f (t, y) (in our case, y would be formed by the vectors of the Fourier coefficients ζ and v), one has to choose a time step h so that for all eigenvalues λ of linearizations, λh is in the stability domain. Since for the SGN equations (and similarly for the WGN equations), the dominant contribution to the linearization is due to orderone operators 4 and since we use a Fourier discretisation where F∂ x = ik with k = -N/2+1, . . . , N/2, the dominant contribution to λ is of the order of N/2. Thus we have to make sure that h = O(1/N ) which is always ensured in our examples.

In our numerical experiments we do not encounter any instability issues which would be marked by exponentially growing modes in time, see again Chapter 10 in [START_REF] Trefethen | Spectral Methods in Matlab[END_REF]. Note that if the stiffness lies in the linear part of the equations as for the Korteweg-de Vries equation, many efficient time integration schemes are known, see for instance [START_REF] Klein | Fourth order time-stepping for low dispersion Korteweg-de Vries and nonlinear Schrödinger equations[END_REF] and references therein and the mentioned code 27 in [START_REF] Trefethen | Spectral Methods in Matlab[END_REF], but here the stiffness is (also) in the nonlinear part.

The system of linear equations in ( 16) is a convolution in the space of Fourier coefficients: the matrix M( ζ) is constructed using (inverse) Fast Fourier Transform and multiplication in physical space. As already mentioned, the inversion will be done with the Krylov approach GMRES [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]. For high accuracy, we use GMRES with up to 100 iterations and a stopping criterion of the iteration of a relative residual of the order of machine precision. As we will show below at examples, GMRES is less of a problem than in the inversion of the Jacobian in the Newton iteration, though we use, unless otherwise stated, the same preconditioner as there.

3.2.

Validation. The accuracy of the code will be validated as discussed in [START_REF] Klein | Fourth order time-stepping for low dispersion Korteweg-de Vries and nonlinear Schrödinger equations[END_REF]: the spatial resolution is controlled via the decay of the Fourier coefficients which is known to be exponential for analytical functions. Thus the order of magnitude of the highest Fourier coefficients gives an indication of the numerical error. The resolution in time will be controlled via conserved quantities of the equations. As mentioned in the introduction, both equations studied here have conserved quantities which will depend during the numerical computation on time because of unavoidable numerical errors. The numerical conservation of these quantities will give an estimation of the numerical error (it generally underestimates this error by 1-2 orders of magnitude, see the discussion in [START_REF] Klein | Fourth order time-stepping for low dispersion Korteweg-de Vries and nonlinear Schrödinger equations[END_REF]). We shall use here the third quantity in, respectively, (1) and (2).

To test the code we consider the solitary waves with velocity c = 2 as an example. First we consider the SGN equations with the initial data (4) for c = 2. We use N = 2 9 Fourier modes for x ∈ 10[-π, π] and N t = 2000 time steps for t ∈ [0, 1]. The relative conservation of the quantities in (1) is of the order of 10 -14 . The Fourier coefficients can be seen on the left of Fig. 7. They decrease to the order of 10 -13 . The difference between the numerical and the exact solution for t = 1 can be seen on the right of Fig. 7. It is of the order of 10 -12 , and the numerical error is thus 4 while the SGN and WGN equations involve operators of order 3, the formulation (14) (thanks to the regularizing effect of inverting the elliptic problem [START_REF] Carter | Bidirectional Whitham equations as models of waves on shallow water[END_REF] and similarly to the Camassa-Holm equation ( 21)) shows that the SGN and WGN equations can be considered as (quasilinear) systems involving operators of order at most one.

just an order of magnitude larger than what is indicated by the Fourier coefficients and the conserved quantities. To test the WGN equations in a similar way, we first numerically construct the solitary wave for c = 2 with N = 2 10 Fourier modes. Then we use this numerical solution as initial data for the WGN equations. This also assesses the accuracy with which the solitary wave is numerically constructed. Again we apply N t = 2000 time steps for t ∈ [0, 1]. The conserved quantities are relatively conserved to the order of 10 -13 . The Fourier coefficients of the solution at the final time can be seen on the left of Fig. 8, the difference with the numerically constructed solitary wave on the right. Obviously we reach the same accuracy as in the SGN case both in terms of resolution in space as indicated by the decay of the Fourier coefficients for large Fourier modes and the difference between numerical and exact solution. Note that in both examples here, a further increase in resolution both in space or time does not lead to higher accuracy.

Stability of the solitary waves

In this section we study the behavior of solutions to both the SGN and the WGN equations set initially as various perturbations of the solitary waves. Let us describe existing results in the literature. In [START_REF] Li | Linear stability of solitary waves of the Green-Naghdi equations[END_REF][START_REF] Li | Hamiltonian structure and linear stability of solitary waves of the Green-Naghdi equations[END_REF], Li investigates the linear stability of the explicit solitary waves of the SGN equations. The first observation is that while solitary waves are critical points of a functional immediately stemming from the Hamiltonian structure of the equation, the second variational derivative of the functional has an infinite-dimensional (essential) negative spectrum, hence standard tools for the nonlinear stability analysis do not apply. This comment also applies to the WGN equations. Then Li proves that solitary waves of the SGN equations with sufficiently small velocity 0 < c -1 1 are (orbitally) linearly stable (see details therein) for infinitely small and exponentially decaying perturbations. The proof is not easily extended to the WGN equations since the differential nature of the operators is used. Let us also mention that Carter and Cienfuegos numerically studied in [START_REF] Carter | The kinematics and stability of solitary and cnoidal wave solutions of the Serre equations[END_REF] the linear stability of cnoidal waves and found that sufficiently large or steep cnoidal waves exhibit linear instability, with relatively small growth rate. By nature, since the unstable modes are periodic with the period being a multiple of the period of the cnoidal wave, the results do not apply to solitary waves; see the discussion in [START_REF] Carter | The kinematics and stability of solitary and cnoidal wave solutions of the Serre equations[END_REF]. Finally, let us also mention the work [START_REF] El | Unsteady undular bores in fully nonlinear shallow-water theory[END_REF] where the modulational stability of small-amplitude bores of the SGN equations is found.

In this section we consider the case of perturbations of solitary waves with moderate up to large amplitudes (we show figures for velocities c = 2, c = 4 and c = 10) and -consistently with the related numerical experiments provided in [START_REF] Mitsotakis | On the nonlinear dynamics of the traveling-wave solutions of the Serre system[END_REF]-find that these solitary waves appear asymptotically orbitally stable.

We consider first perturbations of the solitary waves with velocity c = 2. We work with N = 2 10 Fourier modes for x ∈ 10[-π, π] and N t = 2000 time steps for the time interval t ∈ [0, 10]. The relative conservation of the third quantity in (1) for the SGN equations remains valid up to the order of 10 -10 , and the corresponding one in (2) for the WGN equations up to the order of 10 -9 .

We first study initial data of the form (ζ(x, t = 0), u(x, t = 0) := (ζ c (x), λu c (x)) with λ ∈ R where (ζ c , u c ) is the solitary wave with velocity c. Similar perturbations of the initial data for ζ instead of u lead to similar results, not represented here. The solution to the SGN equations for these initial data (with c = 2) and λ = 0.99 can be seen in in Fig. 9. There is some radiation propagating to the left, but the final state appears to be a solitary wave of slightly different mass than the perturbed solitary wave.

The fact that a solitary wave is approached is even more obvious looking at the L ∞ norm of the solution (computed on collocation points and thus only an approximation of the L ∞ norm since the maximum might not be taken on a grid point. This explains the apparent fluctuations in the subsequent figures.) plotted on the left of Fig. 10. The L ∞ norm is increasing at the beginning and after a few oscillations appears to reach a final value. Since we approximate a situation on the real line by a periodic setting, the radiation cannot escape to infinity here which means that a final state cannot be reached. Extending the computation to a larger time interval, we observe visible oscillations due to the interaction with radiation starting about t = 17. The case λ = 1.01 is very similar in appearance to Fig. 9, therefore we do not show the corresponding figure, just the L ∞ norm on the right of Fig. 10. It can be seen that the norm decreases here before reaching its final value. In Fig. 11 we show the L ∞ norm of the SGN solution with initial data given by (ζ(x, t = 0), u(x, t = 0)) := (ζ 2 (x), u 2 (x) ± 0.01 exp(-x 2 )), on the left for the minus sign in the initial data, on the right for the plus sign. The situation is very similar to the situation of Fig. 10, the final state appears to be a solitary wave of slightly different amplitude than the initial one.

If we consider the same perturbations for the WGN solitary wave (still with c = 2), the resulting figures are very similar as can be seen in Fig. 12. The final state in each example appears to be a solitary wave of slightly different amplitude than the perturbed solitary wave. Obtaining results for c = 4 is already much more computationally demanding. We need to augment the space domain of computation, and hence the number of modes in order to secure a sufficient accuracy, and also the final time of computation before the main wave reaches a final state. Specifically we use N = 2 12 Fourier modes for x ∈ 30[-π, π] and N t = 10 4 time steps for t ∈ [0, 20]. We show the L ∞ norms for the perturbations of the solitary wave of the form , ζ(x, t = 0) = ζ 4 (x) and u(x, t = 0) = λu 4 (x) for λ = 0.99 and λ = 1.01 in Fig. 13. It can be seen that the oscillations are more pronounced in this case than for the case c = 2 in Fig. 10, but that they decrease in amplitude which seems to indicate that the solitary wave is again stable in this case. The study of even larger values of c becomes even more demanding, but can still be handled with the approach based on GMRES, which confirms it as a powerful tool (notice we use the Fourier multiplier 1 + c 4 k 2 /3 as preconditioner). To study the case c = 10, we use N = 2 11 Fourier modes for x ∈ 30[-π, π] and N t = 10 4 time steps for the indicated time intervals. On the left of Fig. 14 we show the L ∞ norm of the solution to the SGN equations for the initial data u(x, t = 0) = 1.01u 10 (x), ζ(x, t = 0) = ζ 10 (x). The solitary wave appears to be again stable. The same is true for a perturbation of the form u(x, 0) = u 10 (x) + 0.01 exp(-x 2 ), although the observed oscillation of the amplitude is about ten times larger than the initial perturbation for u, and with another factor of ten for the oscillation in ζ. If we consider a larger perturbation of the form u(x, 0) = u 10 (x) + 0.1 exp(-x 2 ), the L ∞ norm on the right of Fig. 14 seems to grow beyond what can be seen as a perturbative regime. We expect no blow-up since the norms decrease after some time, but the example makes clear that SGN solitary waves with large velocities c are in applications more easily affected by perturbations than, for instance, Korteweg-de Vries solitons.

This becomes even clearer if we look at the solutions to the SGN equations for these initial data in Fig. 15. A strong growth especially in ζ cannot be interpreted as a perturbation of the initial structure, yet the structures moving to the left do not appear to be solitary waves of smaller amplitude (the bumps in the figures cannot be fitted to a solitary wave).

In all our experiments, considering the solitary waves of the WGN equations instead of the SGN equations yields no qualitative difference. In particular, we observed no coherent state -other than the main wave-emerging from the perturbed solitary wave, which motivates our assertion that solitary waves appear dynamically stable. We also performed a spectral stability analysis (not represented here) essentially amounting to the large-period limit of that in [START_REF] Carter | The kinematics and stability of solitary and cnoidal wave solutions of the Serre equations[END_REF] for cnoidal waves. We found that, while the linearized SGN and WGN equations about cnoidal waves always exhibit unstable modes, the rate of instability decreases as the period grows, and that no unstable mode for solitary waves is captured in this way.

Emergence of modulated oscillations

Both the SGN and WGN equations reduce to the so-called Saint-Venant or shallow-water system in the limit of infinitely long wavelength. This is easily seen when considering the evolution equations of a one-parameter family of initial data varying on a scale of order 1/δ, for times of order 1/δ, where δ 1. Rescaling the coordinates x → δx, t → δt (and once again setting g = d = 1 by similar rescaling) yields for the SGN equations ( 17)

       ∂ t ζ + ∂ x (hu) = 0, ∂ t u -δ 2 3h ∂ x (h 3 ∂ x u) + ∂ x ζ + u∂ x u = δ 2 ∂ x u 3h ∂ x (h 3 ∂ x u) + 1 2 h 2 (∂ x u) 2
, where in an abuse of notation, we have kept the same notation for the functions depending on δ as for the case δ = 1; and for the WGN equations ( 18)

       ∂ t ζ + ∂ x (hu) = 0, ∂ t u -δ 2 3h ∂ x F δ (h 3 ∂ x F δ u) + ∂ x ζ + u∂ x u = δ 2 ∂ x u 3h ∂ x F δ (h 3 ∂ x F δ u) + 1 2 h 2 (∂ x F δ u) 2
, where F δ is the Fourier multiplier defined by

∀ϕ ∈ L 2 (R), F δ ϕ(ξ) = F (δ|ξ|) ϕ(ξ)
where

F (k) = 3 |k| tanh(|k|) - 3 |k| 2 .
Formally setting δ = 0, both [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. I. Model development and analysis[END_REF] and [START_REF] Clamond | Practical use of variational principles for modeling water waves[END_REF] reduce to the aforementioned shallowwater equations ( 19)

   ∂ t ζ + ∂ x (hu) = 0, ∂ t u + ∂ x ζ + u∂ x u = 0.
This formal derivation can be made rigorous for a class of sufficiently regular initial data satisfying the non-cavitation assumption inf R (1 + ζ) > 0, and over a time interval of order 1 (i.e. 1/δ in non-rescaled coordinates); see [START_REF] Lannes | The water waves problem[END_REF]Section 6.1.2]. For longer times, it is well-known that (19) may generate shock singularities, and the large wavelength assumption becomes invalid before the singularity occurs. In this case it is well-documented that solutions to dispersive modified systems may develop zones of rapid modulated oscillations in place of the shocks, which may eventually evolve into fully developed dispersive shock waves or solitary wave resolution, depending on the asymptotic properties of the data as |x| → ∞. The literature on the subject is vast, and we refer to [START_REF] Grava | A numerical study of the small dispersion limit of the Korteweg-de Vries equation and asymptotic solutions[END_REF] for an overview and references in the case of the Korteweg-de Vries (KdV) equation as a perturbation of the inviscid Burgers (iB) equation, where a complete asymptotic description is available, and to [START_REF] El | Dispersive shock waves and modulation theory[END_REF] for an introduction to the modulation theory for more general equations and an extensive list of references. Let us also specifically mention [START_REF] El | Unsteady undular bores in fully nonlinear shallow-water theory[END_REF][START_REF] Tkachenko | Hyperbolicity of the Modulation Equations for the Serre-Green-Naghdi Model[END_REF] for a description of the Whitham modulation theory in the case of the SGN equations, and [START_REF] Dias | On the fully-nonlinear shallow-water generalized Serre equations[END_REF][START_REF] Gavrilyuk | Stationary shock-like transition fronts in dispersive systems[END_REF][START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF][START_REF] Mitsotakis | On the nonlinear dynamics of the traveling-wave solutions of the Serre system[END_REF][START_REF] Mitsotakis | On the Galerkin/finite-element method for the Serre equations[END_REF][START_REF] Pitt | Behaviour of the Serre equations in the presence of steep gradients revisited[END_REF] for some numerical experiments. In contrast to these works, we do not consider here steplike initial data, but study the appearance of modulated oscillations from smooth rapidly decaying initial data, from which solitary wave resolution is the expected large-time asymptotic behavior.

We will observe the emergence of modulated oscillations from initial data set as unidirectional waves for [START_REF] Constantin | Wave breaking for nonlinear nonlocal shallow water equations[END_REF]. Based on Riemann invariants, the solution to [START_REF] Constantin | Wave breaking for nonlinear nonlocal shallow water equations[END_REF] with initial data satisfying u(x, t = 0) = 2 1 + ζ(x, t = 0) -2 can be described as a simple wave with r(x, t) := u(x, t) + 2 1 + ζ(x, t) -2 satisfying the iB equation ( 20)

∂ t r + (1 + 3 4 r)∂ x r = 0.
Any spatially localized smooth solution to (20) develops a shock in finite time. We consider in the following initial data of the form ζ(x, t = 0) = exp(-(x -x 0 ) 2 ), where x 0 is a constant whose role is to avoid the propagation of oscillations beyond the boundary of the computational domain.

We first consider the SGN equations ( 17) with δ = 0.1. We use N = 2 10 Fourier modes on the computational domain 3[-π, π]. We set x 0 = -3, and compute N t = 10 4 time steps for 0 ≤ t ≤ 5. The solution can be seen as a function of time in Fig. 16, and the formation of modulated oscillations in place of the non-dispersive shock is clearly visible.

SGNgausse01water.pdf

SGNgauss1e1zetawater.pdf A cross-section of the plot on the right of Fig. 16, specifically the solution ζ at time t = 5, is shown on the left of Fig. 17. One sees that the rightmost oscillation has already evolved into a localized solitary wave, indicating the onset of solitary wave resolution. The Fourier coefficients on the right of the same figure show that the solution is numerically well resolved.

For smaller values of δ, the modulated oscillations become more localized with stronger gradients. To treat the same initial data as in Fig. 16 for δ = 10 -2 , we use N = 2 12 Fourier modes for x ∈ 2.5[-π, π] and N t = 10 4 time steps for 0 ≤ t ≤ 1.3. The solution to [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. I. Model development and analysis[END_REF] for t = 1.3 can be seen in Fig. 18.

The fact that the solutions from Fig. 18 are more demanding on computational resources is also clear from Fig. 19 where the Fourier coefficients for the solutions in Fig. 18 are shown. Despite higher resolution the Fourier coefficients decrease to the order of 10 -4 (the relative conservation of the third quantity in (1) is of the order of 10 -10 ) which implies a numerical error well below plotting accuracy.

Using the same settings for the WGN equations ( 18) yields qualitatively similar results. We note however that the exponential decay rate for the WGN equations is smaller than the corresponding one for the SGN equation, consistently with the same observation for solitary waves in the preceding section. Therefore to study the case of δ = 0.1 for the WGN equations, we use the same numerical parameters as for the SGN equations except for a higher number of Fourier modes (N = 2 11 ). The solution at t = 5 can be seen in Fig. 20.

The behavior of Fourier coefficients as shown in Fig. 21 is similar to that shown in Fig. 17, although twice the resolution in Fourier space is needed in order to achieve the same decrease. with such a behavior is the Camassa-Holm equation [START_REF] Camassa | An integrable shallow water equation with peaked solitons[END_REF] (

) (1 -5δ 2 12 ∂ 2 x )∂ t w + ∂ x w + 3 2 w∂ x w -δ 2 4 ∂ 3 x w -5δ 2 24 2(∂ x w)(∂ 2 21
x w) + w∂ 3 x w = 0. Here we chose coefficients from [START_REF] Constantin | The hydrodynamical relevance of the Camassa-Holm and Degasperis-Procesi equations[END_REF] (see [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. I. Model development and analysis[END_REF] therein), where it is proved that ( 21) is a higher order (when compared with [START_REF] Constantin | The hydrodynamical relevance of the Camassa-Holm and Degasperis-Procesi equations[END_REF] or the KdV equation) unidirectional model for the SGN equations ( 17) (and hence the WGN equations (18) as well), provided we set

(22) u = w + δ 2 12 ∂ 2 x w + δ 2 6 w∂ 2 x w, ζ = u + 1 4 u 2 + δ 2 6 ∂ t ∂ x u -δ 2 6 u∂ 2 x u -5δ 2 48 (∂ x u) 2 .
This equation generates both dispersive shock waves (see for instance [START_REF] Abenda | Numerical solution of the small dispersion limit of the Camassa-Holm and Whitham equations and multiscale expansions[END_REF][START_REF] Grava | Numerical study of a multiscale expansion of Korteweg-de Vries and Camassa-Holm equation[END_REF]) and finite-time singularities in the form of surging wavebreaking [START_REF] Constantin | Wave breaking for nonlinear nonlocal shallow water equations[END_REF]. Let us finally mention that ( 21) and ( 17) are of the same kind, namely quasilinear nonlocal dispersive equations involving only differential operators.

The set of initial data for which solutions to (21) lead to finite-time wavebreaking contain any smooth and odd function w 0 = w(•, t = 0) such that w 0 (0) < 0 and w 0 (x) < 0 for x > 0 [START_REF] Mckean | Breakdown of the Camassa-Holm equation[END_REF]; we choose here

w 0 (x) = -x exp(-x 2 ).
Inferring initial data for ζ, u through [START_REF] Dinvay | A comparative study of bi-directional Whitham systems[END_REF] would require to know ∂ t ∂ x u(•, t = 0), however we can approximately solve the equations with a harmless O(δ 4 ) approximation by setting [START_REF] Dorodnitsyn | Symmetries, conservation laws, invariant solutions and difference schemes of the one-dimensional green-naghdi equations[END_REF] u(•, t = 0) =:

u 0 = w 0 + δ 2 12 ∂ 2 x w 0 + δ 2 6 w 0 ∂ 2 x w 0 ζ(•, t = 0) = u 0 + 1 4 u 2 0 -δ 2 6 ∂ 2 x (u 0 + 3 4 u 2 0 ) -δ 2 6 u 0 ∂ 2 x u 0 -5δ 2 48 (∂ x u 0
) 2 where we used [START_REF] Dias | On the fully-nonlinear shallow-water generalized Serre equations[END_REF] to infer that

∂ t u = ∂ t w + O(δ 2 ) = -(∂ x w + 3 2 w∂ x w) + O(δ 2 ) = -(∂ x u + 3 2 u∂ x u) + O(δ 2
). We consider the example δ 2 = 0.1. For the computation we use N = 2 11 Fourier modes for x ∈ 5[-π, π] and N t = 10 4 time steps for t ≤ 10. The SGN solution for the initial data (23) can be seen in Fig. 22. Obviously strong gradients appear, but these do not lead to a shock formation.

This example is numerically challenging because of convergence problems for the GMRES algorithm which imply a pollution of the Fourier coefficients at high wave numbers. Therefore we use a dealiasing according to the 2/3-rule which means that the Fourier coefficients corresponding to the highest one third of the wave numbers are put equal to 0. If this is done, the example is well resolved in space as can be seen in Fig. 23 where the Fourier coefficients of the solution in Fig. 22 are shown for t = 10, and in time as inferred by the relative conservation of the third quantity in (1) to the order of 10 -12 . Once again, using the WGN equations (18) instead of [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. I. Model development and analysis[END_REF] does not modify substantially the behavior of the solution, and we do not show the corresponding pictures.

Near-cavitation initial data

In Sections 5 and 6, we used initial data leading to shocks for simplified models, and observed for each scenario that the solutions to both the SGN and WGN equations remained smooth. This leaves open the important question of global-in-time well-posedness of the equations. Let us recall that the local well-posedness of the SGNCHt10uwater01.pdf SGNCHt10zetawater01.pdf SGN equations has been proved in [START_REF] Li | A shallow-water approximation to the full water wave problem[END_REF] 5 -and later on in more general frameworks in [START_REF] Alvarez-Samaniego | A Nash-Moser theorem for singular evolution equations. Application to the Serre and Green-Naghdi equations[END_REF][START_REF] Duchêne | Well-posedness of the Green-Naghdi and Boussinesq-Peregrine systems[END_REF][START_REF] Fujiwara | A shallow water approximation for water waves over a moving bottom[END_REF][START_REF] Israwi | Large time existence for 1D Green-Naghdi equations[END_REF]-but that the technique in these works (essentially energy methods) does not provide any global-in-time result, in particular due to the fact that the functional setting is not controlled by conserved functionals listed in [START_REF] Abenda | Numerical solution of the small dispersion limit of the Camassa-Holm and Whitham equations and multiscale expansions[END_REF]. It should be mentioned that the Boussinesq system obtained when neglecting the nonlinear dispersive terms in (SGN) is known to be globally well-posed provided that the noncavitation assumption inf R (1 + ζ) > 0 is initially satisfied, by [START_REF] Amick | Regularity and uniqueness of solutions to the Boussinesq system of equations[END_REF][START_REF] Schonbek | Existence of solutions for the Boussinesq system of equations[END_REF]; see also [START_REF] Molinet | The classical Boussinesq system revisited[END_REF]. This result does not generalize easily to the SGN or the WGN equations. Very recently, Bae and Granero-Belinchón showed in [START_REF] Bae | Singularity formation for the Serre-Green-Naghdi equations and applications to abcd-Boussinesq systems[END_REF] that if the non-cavitation assumption initially fails to hold at one single point and some symmetry assumptions 5 very recently improved in [START_REF] Inci | On a Lagrangian formulation of the 1D Green-Naghdi system[END_REF] to allow less regular initial data.

are enforced, then solutions to (SGN) (or rather an equivalent reformulation when the non-cavitation assumption holds) preserve these assumptions for positive times and cannot remain smooth globally in time.

Motivated by this result, we consider the initial data [START_REF] Duchêne | Well-posedness of the Green-Naghdi and Boussinesq-Peregrine systems[END_REF] ζ(x, t = 0) = -0.9 exp(-x 2 ), u(x, t = 0) = -x exp(-x 2 ).

The non-cavitation assumption is valid initially and hence will remain valid as long as the solution remains smooth. Indeed we can define, given any x ∈ R and t > 0, h x,t (τ ) = 1 + ζ(X x,t (τ ), τ ) where X x,t (τ ) is the backward characteristic defined by X x,t (τ ) = u(X x,t (τ ), τ ) and X x,t (t) = x. One then observes that, by the conservation of mass equation, one has for any τ ∈ [0, t], h x,t (τ ) = -h x,t (τ )(∂ x u)(X x,t (τ ), τ ), and hence 1 + ζ(x, t) = 1 + ζ(X x,t (0), t = 0) exp t 0 -(∂ x u)(X x,t (τ ), τ ) dτ > 0. As we see in the numerical results below, such initial data produce very steep gradients, and a possible blowup scenario which deserves to be investigated in more details.

To address this question we use N = 2 12 Fourier modes for x ∈ 2.5[-π, π] with dealiasing and N t = 10 4 time steps for t ≤ 3 to solve the SGN equations for the initial data [START_REF] Duchêne | Well-posedness of the Green-Naghdi and Boussinesq-Peregrine systems[END_REF]. The solution can be seen in Fig. 24. To decide whether a blow-up is possible in this case, we show the L ∞ norms of both ζ and ζ x in Fig. 26 in dependence of time. Whereas the L ∞ norm of ζ grows for some time, it reaches a maximum for t ∼ 2.5 and decreases then. Thus there is no L ∞ blow-up, but the strong gradient can be seen in the same figure on the right. But also the L ∞ norm of the gradient appears to reach a finite maximum. This would indicate that one is close to a cusp-like situation, but that the solution stays smooth in this example.

Note that the solution is well resolved: the resolution in time as indicated by the relative conservation of the conserved quantities is of the order of 10 -9 , and the Fourier coefficients of the solution are shown in Fig. 27. The above results do not change within numerical precision if the computation is repeated with N = 2 14 Fourier modes.

Outlook

In this paper we have investigated numerically several aspects of solutions to Serre-Green-Naghdi type equations. First we have obtained a family of supercritical solitary waves of the fully dispersive system with no upper bound on the admissible velocities, alike the explicit family of solitary waves of the original Serre-Green-Naghdi system. Investigating the dynamic stability of these solitary waves, we have found no sign of instability, even for large velocities. We have also set up several experiments for which solutions to either the original or fully dispersive Serre-Green-Naghdi system develop zones of rapid modulated oscillations and/or steep gradients, but we have never monitored finite-time singularity formations.

On the numerical side, we have shown that an approach based on a Fourier spectral method combined with the Krylov subspace iterative method GMRES is very efficient and allows to study with high accuracy computationally demanding problems. From a numerical point of view, there are two main directions of research worth exploring to further improve the code: first better preconditioners for GMRES adapted to the situations to be studied could increase the efficiency (which would be helpful in higher dimensions) and allow to study even more extreme situations which is mainly interesting from a theoretical point of view. Secondly one could improve the time integration by studying stiff integrators for PDEs with stiffness in the nonlinear part, see for instance [START_REF] Caliari | The Leja method revisited: backward error analysis for the matrix exponential[END_REF] and references therein.

Figure 1 .

 1 Figure 1. Left: solitary wave for the WGN equations for c = 1.1 in blue and the SGN equations for the same velocity in red; right: Fourier coefficients for both solitary waves on the left. Now we consider a larger value of the velocity, c = 2 and use N = 2 10 collocation points for x ∈ 20[-π, π].The Newton-GMRES code converges in 4 iterations with a residual of the order of 10 -12 . The resulting solitary wave is shown in Fig.2on the left. Again the solution is very close to the SGN solitary wave shown in the same figure in red. The solution is well resolved in Fourier space as can be seen on the right of Fig.2, and we observe once more that the exponential decay rate of Fourier modes is smaller. The Newton-GMRES iteration behaves similarly to what is described before. Note that it does not converge without a preconditioner due to issues for the high Fourier modes. For even larger values of the velocity such as c = 3, the iteration no longer converges because of GMRES problems for the high Fourier modes.

Figure 2 .

 2 Figure 2. Left: solitary wave for the WGN equations for c = 2 in blue and the SGN equations for the same velocity in red; right: Fourier coefficients for the solitary waves on the left.

  ) that ζ c -and hence h c = 1 + ζ c = c c-uc -scales as c 2 for large c.

Figure 3 .

 3 Figure 3. Left: solitary wave for the WGN equations for c = 20 in blue and the SGN equations for the same velocity in red; right: Fourier coefficients on the left.

Figure 4 .

 4 Figure 4. Left: solitary wave for the WGN equations for c = 100 in blue and the SGN equations for the same velocity in red; right: Fourier coefficients on the left.

Figure 5 .

 5 Figure 5. The function ζ for the solitary wave for the WGN equations in blue and the SGN equations for the same velocity in red: left c = 20, right c = 100.

Figure 6 .

 6 Figure 6. Left: discrete Fourier transform of the Jacobian for the initial iterate (SGN solitary wave) for c = 20; right: non diagonal part Ĵ * on the left. Thus it appears that the problems in the computation of the solitary waves for very high velocities are due to machine precision errors being increased by the multiplication in physical space by the function 1 (1-ηc) 3 taking large values (of order c 6 ), and by the application of ∂ x F which multiplies Fourier modes by ik k F (k k ) ∼ i 3|k k | for large k. But there is no indication of a maximal velocity for the solitary wave solutions to (WGN). The conclusions of the numerical investigation in this section are summarized in the form of Conjecture 2.2.

Figure 7 .

 7 Figure 7. Left: modulus of the Fourier coefficients of the numerical solution to the SGN equations for solitary wave initial data with c = 2 at t=1; right: the difference between this solution to the SGN equations and the exact solution.

Figure 8 .

 8 Figure 8. Left: modulus of the Fourier coefficients of the numerical solution to the WGN equations for solitary wave initial data with c = 2 at t=1; right: the difference between this solution to the WGN equations and the exact solution.

Figure 9 .

 9 Figure 9. Solution to the SGN equations for the initial data u(x, t = 0) = 0.99u c (x) for c = 2 and ζ(x, t = 0) = ζ 2 (x).

Figure 10 .

 10 Figure 10. L ∞ norms of the solutions to the SGN equations for the initial data ζ(x, t = 0) = ζ 2 (x) and u(x, t = 0) = λu 2 (x), on the left for λ = 0.99, on the right for λ = 1.01.

Figure 12 .

 12 Figure 12. L ∞ norms of the solutions to the WGN equations for the initial data ζ(x, t = 0) = ζ 2 (x) and u(x, t = 0) = λu 2 (x) in the upper row, on the left for λ = 0.99, on the right for λ = 1.01; and for u(x, t = 0) = u 2 (x) ± 0.01 exp(-x 2 ) in the lower row, for the minus sign on the left and the plus sign on the right.

Figure 13 .

 13 Figure 13. L ∞ norms of the solutions to the SGN equations for the initial data u(x, t = 0) = λu 4 (x), on the left for λ = 0.99, on the right for λ = 1.01.

Figure 14 .

 14 Figure 14. L ∞ norms of the solutions to the SGN equations for the initial data ζ(x, 0) = ζ 10 (x), and u(x, t = 0) = 1.01u 10 (x) on the left and u(x, t = 0) = u 10 (x) + 0.1 exp(-x 2 ) on the right.

Figure 15 .

 15 Figure 15. Solution to the SGN equations for the initial data u(x, t = 0) = u 10 (x + 40)) + 0.1 exp(-((x + 40) 2 ), on the left u, on the right ζ.

Figure 16 .

 16 Figure 16. Solution to the SGN equations (17) with δ = 0.1 for ζ(x, t = 0) = exp(-(x + 3) 2 ), u(x, t = 0) = 2 1 + ζ(x, t = 0) -2 in dependence of time; on the left u, on the right ζ.

Figure 17 .

 17 Figure 17. Function ζ of Fig. 16 at t = 5 on the left, and the Fourier coefficients on the right.

Figure 18 .

 18 Figure 18. Solution to the SGN equations (17) with δ = 0.01 for ζ(x, t = 0) = exp(-(x + 3) 2 ), u(x, t = 0) = 2 1 + ζ(x, t = 0) -2 at t = 1.3; on the left u, on the right ζ.

6 .Figure 19 .Figure 20 .

 61920 Figure 19. Fourier coefficients of the solutions shown in Fig. 18.

Figure 21 .

 21 Figure 21. Fourier coefficients of the solutions shown in Fig. 20.

Figure 22 .Figure 23 .

 2223 Figure 22. Solution to the SGN equations for the initial data (23) with δ = 0.1.

Figure 24 .

 24 Figure 24. Solution to the SGN equations for the initial data (24), on the left u, on the right ζ.

Figure 25 .

 25 Figure 25. Solution to the SGN equations for the initial data (24) for t = 3, on the left u, on the right ζ.

Figure 26 .

 26 Figure 26. L ∞ norms of the solution ζ to the SGN equations for the initial data (24) on the left and for its gradient on the right.

Figure 27 .

 27 Figure 27. The Fourier coefficients of the solution to the SGN equations for the initial data (24), on the left u, on the right for ζ.

  0 where we recall that h = 1 + ζ and F is the Fourier multiplier operator with symbol F (k) = Sufficiently regular solutions to (SGN) satisfy the above, replacing F by the identity. By standard elliptic theory[START_REF] Lannes | The water waves problem[END_REF] Lemma 5.45], u is uniquely determined by[START_REF] Carter | Bidirectional Whitham equations as models of waves on shallow water[END_REF] from sufficiently regular (v, ζ) with inf R (1 + ζ) > 0, and we can solve[START_REF] Camassa | Long-time effects of bottom topography in shallow water[END_REF] as evolution equations for (ζ, v). Incidentally, notice v represents a rescaled tangential fluid velocity at the free interface, the fourth quantity in (1) and (2).3.1. The numerical scheme. Many numerical schemes have been proposed for solving the SGN equations; see for instance[START_REF] Aïssiouene | A two-dimensional method for a family of dispersive shallow water models[END_REF][START_REF] Antonopoulos | On the well-posedness of the Galerkin semidiscretization of the periodic initial-value problem of the Serre equations[END_REF][START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. I. Model development and analysis[END_REF][START_REF] Dorodnitsyn | Symmetries, conservation laws, invariant solutions and difference schemes of the one-dimensional green-naghdi equations[END_REF][START_REF] Dutykh | Finite volume and pseudo-spectral schemes for the fully nonlinear 1D Serre equations[END_REF][START_REF] Gavrilyuk | Stationary shock-like transition fronts in dispersive systems[END_REF][START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF][START_REF] Li | High order well-balanced CDG-FE methods for shallow water waves by a Green-Naghdi model[END_REF][START_REF] Mitsotakis | On the nonlinear dynamics of the traveling-wave solutions of the Serre system[END_REF][START_REF] Mitsotakis | On the Galerkin/finite-element method for the Serre equations[END_REF][START_REF] Pitt | Behaviour of the Serre equations in the presence of steep gradients revisited[END_REF].

3

|k| tanh(|k|) -3 |k| 2 ; and v and u are related through the elliptic equation

(15) v = u -1 3h ∂ x h 3 ∂ x Fu .

By scaling arguments (specifically settingζ(x, t) = d ζ( 1 d x, √ gd d t), u(x, t) = √ gdũ( 1 d x, √ gd d t)) it is always possible to set g = d = 1, and we shall do so in the following.

This work is partially supported by the ANR-FWF project ANuI -ANR-17-CE40-0035, the isite BFC project NAANoD, the EIPHI Graduate School (contract ANR-17-EURE-0002) and by the European Union Horizon 2020 research and innovation program under the Marie Sklodowska-Curie RISE 2017 grant agreement no. 778010 IPaDEGAN. 1