David Lannes
email: david.lannes@math.u-bordeaux.fr

Fabien Marche
email: fabien.marche@math.univ-montp2.fr

NONLINEAR WAVE-CURRENT INTERACTIONS IN SHALLOW WATER

Keywords: Water waves, shallow water, Green-Naghdi, Boussinesq, nonlinear dispersive equations, vorticity, solitary waves, Finite-Volume discretization

 that the motion of the waves could be described using an extended Green-Naghdi system. In this paper we propose an analysis of these equations, and show that they can be used to get some new insight into wave-current interactions. We show in particular that solitary waves may have a drastically different behavior in the presence of vorticity and show the existence of solitary waves of maximal amplitude with a peak at their crest, whose angle depends on the vorticity. We also show some simple numerical validations. Finally, we give some examples of wave-current interactions with a non trivial vorticity field and topography effects.

1. Introduction 1.1. General setting Several models have been derived for the description of nearshore dynamics. One of the most widely spread is certainly the Nonlinear Shallow Water (NSW) model which is a nonlinear hyperbolic system coupling the time evolution of the surface elevation ζ to the vertically averaged horizontal component of the velocity v. This system is derived from the free surface Euler equations by averaging in the vertical direction and neglecting all the terms of order O(µ), where the shallowness parameter µ is defined as

µ = H 2 0 L 2 =
(typical depth) 2 (horizontal length scale) 2 .

The NSW equations are however not fully satisfactory since they neglect all the dispersive effects that play a very important role in many situations, and in particular during the shoaling phase. These dispersive terms are of order O(µ) and are therefore neglected by the NSW equations. Keeping them in the equations, and neglecting only the O(µ 2 ) terms, one obtains a more accurate -but mathematically and numerically more complicated -set of equations known as the Serre [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], or Green-Naghdi [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF][START_REF] Kim | A derivation of the Green-Naghdi equations for irrotational flows[END_REF], or fully nonlinear Boussinesq [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. I. Highly nonlinear unsteady waves[END_REF] equations. We shall refer to these models here as the Green-Naghdi (GN) equations. Contrary to the weakly nonlinear Boussinesq models that go back to Boussinesq himself, no smallness assumption is made on the size of the surface perturbations. We refer to [START_REF] Lannes | The water waves problem: mathematical analysis and asymptotics[END_REF] for a rigorous derivation and a mathematical justification (in the sense that their solutions remain close to the exact solution of the free surface Euler equations) of all these models. If the numerical approximation of various Boussinesq-type equations has attracted a lot of attention for the last 20 years (see for instance among the recent studies [START_REF] Bjorkavag | Wave breaking in Boussinesq models for undular bores[END_REF][START_REF] Erduran | Further application of hybrid solution to another form of Boussinesq equations and comparisons[END_REF][START_REF] Eskilsson | Spectral/hp discontinuous Galerkin methods for modelling 2D Boussinesq equations[END_REF][START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF][START_REF] Engsig-Karup | DG-FEM solution for nonlinear wave-structure interaction using Boussinesq-type equations[END_REF][START_REF] Ricchiuto | Upwind Residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF][START_REF] Roeber | Boussinesq-type model for energetic breaking waves in fringing reef environments[END_REF][START_REF] Shiach | A temporally second-order accurate Godunov-type scheme for solving the extended Boussinesq equations[END_REF][START_REF] Soares-Frazao | A second-order semi-implicit hybrid scheme for one-dimensional Boussinesq-type waves in rectangular channels[END_REF]), it is mostly recently that discrete formulations for the GN equations have been proposed. Denoting by d the horizontal dimension, we can refer for instance, in the case d = 1, to [START_REF] Bonneton | Recent advances in Serre-Green-Naghdi modelling for wave transformation, breaking and runup processes[END_REF][START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF] for hybrid Finite-Volume (FV) and Finite-Difference (FD) discretizations, [START_REF] Duran | Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations[END_REF][START_REF] Li | High order well-balanced CDGFE methods for shallow water waves by a Green-Naghdi model[END_REF][START_REF] Panda | Discontinuous Galerkin methods for solving Boussinesq-Green-Naghdi equations in resolving non-linear and dispersive surface water waves[END_REF] for discontinuous-Galerkin (dG) formulations, [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part II: Boundary conditions and validations[END_REF] for a compact FV approach, [START_REF] Mitsotakis | On the Galerkin/Finite-Element Method for the Serre Equations[END_REF] for a Finite-Element (FE) approach on flat bottom or [START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation[END_REF] for an hybrid FV-FE formulation. There is even less studies in the case d = 2, see [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF][START_REF] Metayer | A numerical scheme for the GreenNaghdi model[END_REF][START_REF] Shi | A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation[END_REF]. These equations have also been adapted to handle wave breaking by adding an artificial viscous term to the momentum equation (see for instance [START_REF] Chen | Boussinesq modeling of wave transformation, breaking, and run-up. II: 2d[END_REF][START_REF] Cienfuegos | A wave-breaking model for Boussinesq-type equations including mass-induced effects[END_REF][START_REF] Kennedy | Boussinesq modeling of wave transformation, breaking, and run-up. I: 1[END_REF][START_REF] Schäffer | A Boussinesq model for waves breaking in shallow water[END_REF]) or by locally switching to the NSW equations in the vicinity of broken waves and using shock capturing schemes [START_REF] Kazolea | Numerical treatment of wave breaking on unstructured finite volume approximations for extended Boussinesq-type equations[END_REF][START_REF] Tissier | A new approach to handle wave breaking in fully non-linear Boussinesq models[END_REF][START_REF] Tonelli | Finite volume scheme for the solution of 2D extended Boussinesq equations in the surf zone[END_REF]. We refer to [START_REF] Bonneton | Recent advances in Serre-Green-Naghdi modelling for wave transformation, breaking and runup processes[END_REF] for a recent review on these aspects.

The Green-Naghdi and more generally most of the Boussinesq-type models rely on the assumption that the flow is irrotational or almost irrotational. Such an assumption is satisfied in most configurations but may fail in the surf zone where wave breaking can create vorticity currents (rip currents) or in present of some underlying current. The difficulty to describe wave motion in the presence of vorticity is that the dynamics of the flow is in general genuinely (d + 1)-dimensional while in the irrotational framework the dynamics is only d-dimensional (vertical averaging has been used to remove the vertical variable).

It is shown in [START_REF] Chen | Boussinesq modeling of longshore currents[END_REF] that the GN equations can describe rotational flows with purely vertical vorticity. In horizontal dimension d = 1, it is shown in [START_REF] Musumeci | The flow in the surf zone: a fully nonlinear Boussinesqtype of approach[END_REF] that vorticity is responsible for the presence of an additional term in the momentum equation, which is coupled to the standard vorticity equation (see also [START_REF] Kim | A depth-integrated model for weakly dispersive, turbulent, and rotational fluid flows[END_REF][START_REF] Kim | A σ-coordinate transport model coupled with rotational Boussinesqtype equations[END_REF] for a related approach).

Following the original approach by Green and Naghdi [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF], several authors [START_REF] Ertekin | Waves caused by a moving disturbance in a shallow channel of finite width[END_REF][START_REF] Shields | On direct methods in water-wave theory[END_REF][START_REF] Zhao | Application of higher-level GN theory to some wave transformation problems[END_REF] assumed a polynomial structure of the velocity profile and solved the mass and momentum equations projected on such a basis of functions. This approach is compatible with the presence of vorticity. An interesting and recent refinement for a better treatment of the surf zone consists in coupling this approach with shallow water asymptotics [START_REF] Zhang | Boussinesq-Green-Naghdi rotational water wave theory[END_REF].

We follow in this paper another approach initially developed in [START_REF] Castro | Fully nonlinear long-waves models in the presence of vorticity[END_REF] where it is shown that additional terms are necessary in the momentum equation in the presence of vorticity. Contrary to other approaches, these additional terms are determined through the resolution of d-dimensional evolution equations and do not require the resolution of the (d + 1)-dimensional vorticity equation. The procedure is reminiscent of turbulence theory with the difference that no artificial closure is needed here: based on the controls on the solutions of the full Euler equations established in [START_REF] Castro | Well-posedness and shallow water stability for a new hamiltonian formulation of the water waves equations with vorticity[END_REF], one can show that the cascade of equations is finite at the precision of the model (we also point out the related work [START_REF] Richard | Gavrilyuk Modelling turbulence generation in solitary waves on shear shallow water flows[END_REF] where a physical modeling of the closure is used instead to handle turbulent bores). The resulting equations are an extended GN system with additional advection-like equations for the vorticity related terms. This paper is concerned with the horizontal one dimensional case d = 1. Its first goal is to show that this approach can be used to get some new insight into wavecurrent interactions; we show for instance that the behavior of solitary waves can be drastically different in the presence of vorticity, leading to extremal peaked solitary waves with an angle at the crest that depends on the vorticity. The second goal of this paper is to propose a simple numerical scheme to numerically solve these extended GN equations and to highlight that despite the fact that these equations of motion are purely d-dimensional, they can be used to reconstruct the internal velocity field, even in the presence of non trivial vorticity and topography.

The models

In the case where the horizontal dimension d is equal to one and following [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF], the irrotational GN equations can be formulated as

∂ t h + ∂ x • (hv) = 0, (I + T)[∂ t (hv) + ∂ x (hv 2 )] + gh∂ x ζ + hQ 1 (v) = 0, (1.1)
where we recall that ζ is the elevation of the wave with respect to the rest level and that v is the vertically averaged horizontal velocity, while g stands for the acceleration of gravity and h is the total water height

h = H 0 + ζ -b,
where {z = H 0 -b(x)} is a parametrization of the bottom (see Figure 1.1). Finally, the linear operator

T = T[h,b] and the quadratic form Q 1 (•) = Q 1 [h,b](•) are defined by T[h,b]w = hT [h,b]( 1 h w), (1.2)
T [h,b]w = - 1 3h ∂ x (h 3 ∂ x w) + 1 2h ∂ x (h 2 w∂ x b) -h 2 ∂ x w∂ x b + w(∂ x b) 2 , (1.3) Q 1 [h,b](v) = 2 3h ∂ x (h 3 (∂ x v) 2 ) + h(∂ x v) 2 ∂ x b + 1 2h ∂ x (h 2 v 2 ∂ 2 x b) + v 2 ∂ 2 x b∂ x b. (1.4)
This formulation does not require the computation of any third-order derivative, allowing for more robust numerical computations, especially when the waves become steeper.

Note also that if one removes the operator T and the nonlinearity Q 1 from the second equation in (1.1), the model reduces to the standard NSW equations; these two terms accounts therefore for the O(µ) dispersive and nonlinear terms specific to the GN equations.

One of the main features of the Green-Naghdi model is that it allows the description of d + 1 dimensional waves (d being the horizontal dimension) by a set of d-dimensional equations (independent on the vertical variable z), hereby leading to considerable gains in mathematical simplicity and computational time. The d-dimensional nature of the flow is due to the fact that the flow is assumed to be irrotational; indeed, the velocity field U in the fluid domain then derives from a scalar velocity potential Φ (i.e. U = ∇ X,z Φ) and as remarked by Zakharov [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF] and Craig-Sulem [START_REF] Craig | Numerical simulation of gravity waves[END_REF] the free surface (d + 1)dimensional Euler equations can then be reduced to an Hamiltonian system coupling the surface elevation ζ to ψ, the trace at the surface of the velocity potential. Both ζ and ψ depend only on time and on the (d-dimensional) horizontal variable X. The Green-Naghdi equation being obtained by an asymptotic expansion in terms of the shallowness parameter µ of the free surface Euler equations (see [START_REF] Alvarez-Samaniego | Large time existence for 3d water-waves and asymptotics[END_REF][START_REF] Lannes | The water waves problem: mathematical analysis and asymptotics[END_REF] for a full mathematical justification of this approximation), it is no surprise that they are also d dimensional.

In presence of vorticity, the situation is drastically different since the dynamics of the vorticity ω = curl U is in general fully (d + 1)-dimensional. The Zakharov-Craig-Sulem formulation has recently been generalized in [START_REF] Castro | Well-posedness and shallow water stability for a new hamiltonian formulation of the water waves equations with vorticity[END_REF] to the rotational case; this generalization, also formally hamiltonian, couples the evolution of ζ and ψ as in the irrotational case 1 , but this evolution is now also coupled to the evolution of the vorticity field which depends in general on all the space variables. One should therefore expect that generalizations of the Green-Naghdi equations in presence of vorticity have a full (d + 1)dimensional dependence in the space variables, hereby implying a considerable increase of computational time. It has been shown recently in [START_REF] Castro | Fully nonlinear long-waves models in the presence of vorticity[END_REF] that this is not the case. In the case of a constant vorticity, that is, when curl U = (0,ω,0) T with ω(t,x,z) = ω 0 = cst, this is not surprising because there is no z dependence coming from the equation on the vorticity. The vorticity field however induces a shear which, together with the dispersive effects, make the horizontal velocity depart from its vertical average. Because of this effect, the Green-Naghdi equations (1.1) must be replaced 2 by

     ∂ t h + ∂ x (hv) = 0, (1 + T)[∂ t (hv) + ∂ x (hv 2 )] + gh∂ x ζ + hQ 1 (v) + ∂ x ( 1 12 h 3 ω 2 0 ) + hC(ω 0 h,v) = 0, (1.5)
with C(ω 0 h,v) obtained by taking v = ω 0 h in the following expression

C(v ,v) = - 1 6h ∂ x 2h 3 v ∂ 2 x v + ∂ x (h 3 v )∂ x v . (1.6)
For the case of a general vorticity, that is, when (in horizontal dimension d = 1), curl U = (0,ω,0) T with ω(t,x,z) = ∂ z u -∂ x w

(and U = (u,0,w) T ), the vorticity ω satisfies the transport equation

∂ t ω + (u∂ x + w∂ z )ω = 0, (1.7)
in which the z dependence cannot be removed. The fact that one can however derive zindependent Green-Naghdi type models in this framework is therefore more surprising. The one-dimensional GN system in the presence of a general vorticity is then given 3 by

               ∂ t h + ∂ x (hv) = 0, (1 + T)[∂ t (hv) + ∂ x (hv 2 )] + gh∂ x ζ + hQ 1 (v) + ∂ x E + hC(v ,v) = 0, ∂ t v + v∂ x v + v ∂ x v = 0, ∂ t E + v∂ x E + 3E∂ x v + ∂ x F = 0, ∂ t F + v∂ x F + 4F ∂ x v = 0.
(1.8)

1 Note however that in presence of vorticity, the velocity field U does not derive from a scalar potential, and that an alternative definition is needed for ψ. Namely, it is defined such that ∇ψ is the projection onto (horizontal) gradient vector field of the horizontal component of the tangential velocity at the surface. 2 The model (1.5) actually holds under a rather weak smallness assumption on the topography. More precisely, it is assumed that εβµ In full generality, an additional topography term is needed, see (1.27) in [START_REF] Castro | Fully nonlinear long-waves models in the presence of vorticity[END_REF]. 3 The same smallness assumption as for (1.5) is made on the topography.

Let us briefly comment on this model (we refer to [START_REF] Castro | Fully nonlinear long-waves models in the presence of vorticity[END_REF] for more details). Due to the presence of the vorticity, there is a vertical dependence of the velocity field inside the fluid domain that may interact nonlinearly with the vertical dependence coming from the dispersive terms. After vertical averaging, this interaction is responsible for the term hC(v ,v) in the momentum equation; it is given by (1.6). The difference with the case of a constant vorticity is that v is now defined as a second order momentum of the vorticity induced shear velocity,

v = 12 h 3 ζ -H0+b (z + H 0 -b) 2 v * sh with v * sh = - ζ z ω + 1 h ζ -H0+b
ω.

(1.9)

Even though v is defined in terms of ω, it is not necessary to solve the (1 + 1)dimensional vorticity equation (1.7) to compute it; indeed, it is shown in [START_REF] Castro | Fully nonlinear long-waves models in the presence of vorticity[END_REF] that it can be determined from its initial value by solving the third equation of (1.8).

Similarly, the term ∂ x ( 1 12 h 3 ω 2 0 ) that appears in (1.5) is now replaced by ∂ x E, where E is a second order tensor (represented by a function here since we are in dimension 1) describing the self quadratic interaction of the vorticity induced shear inside the fluid domain,

E = ζ -H0+b (v * sh ) 2 .
(1.10)

Here again, one wants to be able to compute E without appealing to the vorticity equation (1.7). The strategy adopted in [START_REF] Castro | Fully nonlinear long-waves models in the presence of vorticity[END_REF] is inspired by an analogy with turbulence theory and recent works on roll waves and hydraulic jumps [START_REF] Richard | A new model of roll waves: comparison with Brocks experiments[END_REF][START_REF] Richard | The classical hydraulic jump in a model of shear shallow-water flows[END_REF]. The tensor E is viewed as a "Reynolds" tensor where the "averaging" is in the present case the vertical integration. Looking for an equation on E one obtains a cascade of equations involving tensors of increasing order; but unlike turbulence theory, there is no need for an artificial closure of the cascade. Indeed it can be proved that the contribution of the fourth order and higher tensors are below the overall O(µ 2 ) precision of the model and can therefore be neglected. The last two equations in (1.8) furnish this finite cascade of equations on the second order tensor E and the third order tensor F (here again represented by a function in dimension 1) defined as

F = ζ -H0+b (v * sh ) 3 . (1.11)
The generalization of (1.8) to two-dimensional surfaces is also given in [START_REF] Castro | Fully nonlinear long-waves models in the presence of vorticity[END_REF] but our focus is here on the analysis of some properties of (1.8) as well as the development of a numerical code to compute its solutions.

Remark 1.1. As explained above, the model (1.8) is precise up to O(µ 2 ) terms; lowering the precision to O(µ 3/2 ), one can work with the simpler model

     ∂ t h + ∂ x (hv) = 0, (1 + T)[∂ t (hv) + ∂ x (hv 2 )] + gh∂ x ζ + hQ 1 (v) + ∂ x E = 0, ∂ t E + v∂ x E + 3E∂ x v = 0.
(1.12)

Remark 1.2. Keeping the precision O(µ 2 ), a simplified model can also be obtained in the situation where F is initially almost equal to zero (this is the case when the vorticity is constant or in the situation considered in §5.5 for instance). Removing F from (1.8) one then obtains the reduced model

         ∂ t h + ∂ x (hv) = 0, (1 + T)[∂ t (hv) + ∂ x (hv 2 )] + gh∂ x ζ + hQ 1 (v) + ∂ x E + hC(v,v ) = 0, ∂ t v + v∂ x v + v ∂ x v = 0, ∂ t E + v∂ x E + 3E∂ x v = 0.
(1.13)

Organization of the paper

In Section 2, we study the existence of solitary waves for the GN system with vorticity (1.8). We show in §2.1 that the existence of smooth solitary waves can be reformulated as an ODE problem. The existence of solutions is then established in §2.2 where we also comment on the qualitative differences with the irrotational case. For instance, while there are solitary waves of arbitrary amplitude for the standard GN equations, there are configurations with non trivial vorticities for which solitary waves cannot exceed a critical amplitude. Solitary waves of critical amplitude are then studied in §2.3 where we show that these extremal solitary waves have a peak at their crest, whose angle depends on the vorticity.

We then present in Section 3 the numerical scheme we propose to solve (1.8). After a simple renormalisation of the system using the mass conservation equation, we present in §3.1 a simple splitting scheme inspired by previous works on the standard GN equations. This splitting involves a conservative propagation step and a dispersive correction step. The conservative step is studied in §3.2; as in the irrotational case, it is of hyperbolic type but because of the extra unknowns due to the vorticity, its structure is more complicated. In particular, there are now three wave speeds (instead of two in the irrotational case). A corresponding finite volume scheme is proposed, for which higher order extensions are constructed. The study of the dispersive step being similar to the irrotational case, we just briefly recall the main points in §3.3. Section 4 is then devoted to the numerical validation of this scheme. The different kinds of smooth solitary waves predicted in Section 2 are numerically observed in §4.1 and used to evaluate the convergence rate. We also observe numerically in §4.2 the existence of the extremal peaked solitary waves exhibited in §2.3. We provide in §4.3 a numerical simulation involving a non-flat topography; we use this example to show that the vorticity may play a considerable role on the shoaling of waves.

Finally, we detail in Section 5 how the system of equations (1.8) can be used to describe the dynamics of the (d + 1) velocity field U = (u,w) T at any time, up to a O(µ 3 2 ) accuracy, and show how the previous discrete formulation of Section 3 may be easily modified to perform the corresponding velocity reconstructions. This process is illustrated by two prospective examples of wave-current interactions, involving a non trivial vorticity field and topography effects.

Solitary waves

We investigate here the existence of solitary waves for the Green-Naghdi equations with vorticity (1.8). We show in §2.1 that smooth solitary waves must satisfy a second order ODE. This ODE is solved in §2.1, while extremal peaked solutions are studied in §2.3.

Derivation of the ODE for the shape of the solitary waves

Our purpose here is to show that smooth solitary waves, if they exist, must satisfy a second order ODE. We consider here flat bottoms (i.e. b = 0), and the system (1.8) can therefore be written

                     ζ t + (hv) x = 0, v t + g∂ x ζ + v∂ x v + 1 h E x - 1 6h 2h 3 v v xx + (h 3 v ) x v x x = 1 3 1 h h 3 (v xt + vv xx -v 2 x ) x v t + (vv ) x = 0, E h 3 t + v E h 3 x + 1 h 3 F x = 0, F h 4 t + v F h 4 x = 0,
(2.1) with h = H 0 + ζ. We look for solitary waves solutions to (2.1), i.e. solutions of the form

(ζ,v,v ,E,F )(t,x) = (ζ,v,v ,E,F )(x -ct),
for some constant c ∈ R, and with ζ and v vanishing at infinity, over a current that might not vanish at infinity, that is, we assume that lim

±∞ (ζ,v) = 0 and lim ±∞ (v ,E,F ) = (v ∞ ,E ∞ ,F ∞ ) for some constants E ∞ , F ∞ , v ∞ .
Such solutions should satisfy (for the sake of clarity, we do not underline the functions in the expressions below)

                       [-(c -v)h] x = 0, -(c -v)v x + g∂ x ζ + E x h - 1 6h 2h 3 v v xx + (h 3 v ) x v x x = - 1 3h h 3 ((c -v)v xx + v 2 x ) x (c -v)v x = 0, -(c -v) E h 3 x + 1 h 3 F x = 0, (c -v) F h 4 x = 0. (2.
2) Integrating the first equation, and using the fact that ζ and v vanish at infinity, one readily deduces

(c -v)h = cH 0 . (2.3)
Multiplying the second equation by h and integrating in x, we therefore get

-cH 0 v + g 2 (h 2 -H 2 0 ) + (E -E ∞ )- 1 6 2h 3 v v xx + (h 3 v ) x v x (2.4) = - 1 3 h 2 cH 0 v xx + h 3 v 2 x ,
and we need to determine v and E. Let us proceed first with v . From the third equation in (2.2), we have

(c -v)v = cv ∞ , which, together with (2.3), yields v = h H 0 v ∞ . (2.5)
We now turn to derive an expression for E. From the last equation, we get that

F = h 4 H 4 0 F ∞ . (2.6)
Together with (2.3), this allows one to rewrite the fourth equation as

-cH 0 E h 3 x + 4 F ∞ H 4 0 hh x = 0,
and therefore

E = E ∞ + ( h 3 H 3 0 -1)E ∞ + 2 F ∞ c (h 2 -H 2 0 )h 3 H 5 0 .
(2.7)

Plugging (2.5) and (2.7) into (2.4) we obtain

-cH 0 v + g 2 (h 2 -H 2 0 ) + ( h 3 H 3 0 -1)E ∞ + 2 F ∞ c (h 2 -H 2 0 )h 3 H 5 0 - h 2 3H 0 v ∞ [h 2 v x ] x = - 1 3 h 2 cH 0 v xx + h 3 v 2 x . Since (2.3) implies that h 2 v x = cH 0 h x and v = c h-H0 h , we deduce further that -c 2 H 0 h -H 0 h + g 2 (h 2 -H 2 0 ) + ( h 3 H 3 0 -1)E ∞ + 2 F ∞ c (h 2 -H 2 0 )h 3 H 5 0 - 1 3 cv ∞ h 2 h xx = - c 2 3 H 2 0 h 1 h h x x .
This leads us to the following definition of a (smooth) solitary wave:

Definition 2.1. A solitary wave of speed c for (2.1) is a mapping

(t,x) ∈ R 2 → (ζ,v,v ,E,F )(x -ct) such that there exists h ∈ C 2 (R) and E ∞ > 0, v ∞ ∈ R and F ∞ ∈ R such that ζ = h -H 0 , v = c h -H 0 h , v = h H 0 v ∞ , (2.8)
E = h 3 H 3 0 E ∞ + 2 F ∞ c (h 2 -H 2 0 )h 3 H 5 0 , F = h 4 H 4 0 F ∞ (2.9)
and h solves the ODE

1 3 c(cH 2 0 -v ∞ h 2 )h xx = h -H 0 2h (2c 2 H 0 -gh(h + H 0 )) -( h 3 H 3 0 -1)E ∞ (2.10) -2 F ∞ c (h 2 -H 2 0 )h 3 H 5 0 + c 2 3 H 2 0 h 2 x h
on R and satisfies lim ±∞ h = H 0 . The function h is called the profile of the solitary wave.

Existence of smooth solitary waves

We first consider here the case where F ∞ = 0 and prove the existence of solitary waves in the sense of Definition 2.1.

Proposition 2.2. Let E ∞ > 0, v ∞ ∈ R and F ∞ = 0. Let also h max > H 0 .
i. Up to translations, there can be at most two solitary waves of maximal height h max for (2.1); if they exist, they have opposite speed ±c, with

c = gh max + h max (h max + 2H 0 ) H 3 0 E ∞ 1/2
.

ii. The solitary wave of speed c (resp. -c) exists if and only if the following condition holds

c(cH 2 0 -v ∞ h 2 max ) > 0 (resp. -c(-cH 2 0 -v ∞ h 2 max ) > 0).
The profile of the solitary wave then attains its maximal value at a unique point x max and it is symmetric with respect to the axis x = x max and decaying on the half-line x > x max .

Proof.

Step 1. We derive here an expression for h 2 x in terms of h. For later investigations, we deal with the general case F ∞ ∈ R here. Multiplying the differential equation (2.10) by h x and dividing by h 2 , we get

- c 2 3 H 2 0 h 3 x h 3 + 1 3 c 2 H 2 0 1 h 2 h xx h x - 1 3 cv ∞ h xx h x = h -H 0 2h 3 (2c 2 H 0 -gh(h + H 0 ))h x -( h 3 H 3 0 -1) h x h 2 E ∞ -2 F ∞ c (h 2 -H 2 0 ) H 5 0 hh x .
After integrating in x, this yields

1 6 c cH 2 0 -v ∞ h 2 ) h 2 x h 2 = c 2 2 1 h 2 (h -H 0 ) 2 - 1 2 g 1 h (h -H 0 ) 2 - 1 2H 3 0 1 h (h -H 0 ) 2 (h + 2H 0 )E ∞ - F ∞ 2c (h 2 -H 2 0 ) 2 H 5 0
(the integration constant has been chosen in order to respect the constraint that ζ and its derivatives vanish at infinity), or equivalently

c 3 cH 2 0 -v ∞ h 2 h 2 x = (h -H 0 ) 2 c 2 -gh - h(h + 2H 0 ) H 3 0 E ∞ - h 2 (h + H 0 ) 2 H 5 0 F ∞ c . (2.11)
Step 2. Expressions for the velocity and qualitative analysis. By definition, if h is the profile of a solitary wave then it is a C 2 -function and its derivative must vanish at its maximum. The formula (2.11) then provides directly the only two (recall that F ∞ is assumed to be zero here) possible values for the speed c. Since the function h → gh + h(h+2H0) H 3 0 is strictly increasing on R + , we also deduce from (2.11) that h x cannot vanish at another point, and therefore that the maximum of h is attained at a unique point x max , and further, that h has to be monotonous on both sides of x max . The fact that it is decaying on x > x max follows from the condition that h → H 0 < h max at infinity. Finally, the fact that h is symmetric with respect to x max follows from the simple observation that if h solves (2.10) for x ≥ x max with boundary conditions h(x max ) = h max and h x (x max ) = 0 then x → h(2x max -x) furnishes a solution for x ≤ x max .

Step 3. Existence of a solitary wave of speed c = c or c = -c. If cH 2 0 -v ∞ h 2 max = 0, then the Cauchy-Lipschitz theorem furnishes a local solution with boundary condition h(x max ) = h max and h x (x max ) = 0. If moreover c(cH 2 0 -v ∞ h 2 max ) > 0, then it is easy to deduce from (2.10) that this local solution satisfies h (x max ) < 0 and therefore that the solution attains a local maximum at x max . Proceeding as in Step 2, one gets that this local solution is symmetric with respect to x max and decaying on x > x max . Moreover, one always has h > H 0 ; indeed, if one had h(x 0 ) = H 0 for some x 0 ∈ R, then one would have h x (x 0 ) = 0 by (2.11), and by uniqueness, one would have h ≡ H 0 , which is absurd. Therefore h decays to some limit as x → ∞, and this limit is necessarily H 0 by (2.11). The identity (2.11) also shows that h x remains bounded, so that no blow up of h nor h x can occur and the solution of the ODE (2.10) is global.

Step 4. Non existence of a solitary wave of speed

c = c or c = -c. If c(cH 2 0 - v ∞ h 2
max ) < 0, then it is easy to deduce from (2.11) that no solitary wave can exist. The only case left to investigate is therefore the critical case cH 2 0 -v ∞ h 2 max = 0. In this case, one gets from (2.10) that

c 2 3 H 2 0 h (x max ) h(x max ) = h max -H 0 2h max (gh max (h max + H 0 ) -2c 2 H 0 ) + ( h 3 max H 3 0 -1)E ∞ > 0
which contradicts the assumption that h is a C 1 -function attaining its maximum at x max .

Remark 2.1. If in addition to the assumption F ∞ = 0 we take E ∞ = v ∞ = 0 in the statement of Proposition 2.2, then one has v = E = F = 0 in Definition 2.1 so that the solitary waves are the same as in the irrotational setting for which it is well known that explicit solitary waves exist. More precisely, for any maximal amplitude h max > H 0 , there exists two (up to translations) solitary waves of speed c = ± √ gh max and with the same profile h + = h -given by the resolution of (2.10) and which can in this particular case be computed explicitly,

h(x) = H 0 + εH 0 (sech( x λ )) 2 with λ = 2 √ 3 1 + ε ε H 0 , (2.12)
and where we denoted

h max = H 0 (1 + ε).
Remark 2.2. If in addition to the assumption F ∞ = 0 we take v ∞ = 0 but consider the case E ∞ > 0, the situation is qualitatively the same as in Remark 2.1: for any maximal amplitude h max = H 0 (1 + ε) > H 0 , there exists two solitary waves of same shape and of opposite speed c = ±c. The only difference is that the speed c is larger than in the irrotational case,

c = gh max + h max (h max + 2H 0 ) H 3 0 E ∞ 1/2 , (2.13)
and that the solitary waves becomes narrower as E ∞ increases (this follows easily from the comparison principle for ODEs); see Figure 2.1. Remark 2.3. If we assume that F ∞ = 0 but E ∞ > 0 and v ∞ > 0 (the case v ∞ < 0 can be treated in a similar way) in the statement of Proposition 2.2 then there are two major qualitative changes with respect to the situation considered in Remark 2.2. The first one is that right going solitary waves do not exist for any maximal amplitude h max > H 0 . Indeed, the criterion given in the second point of the proposition is always satisfied for the left-going solitary wave, but requires for the right-going one that

-6 -4 -2 0 2 4 6 1 1.5 2 x (m) (m) E =1 m 3 .s -2 E =5 m 3 .s -2 Irrotational Green-Naghdi
cH 2 0 -v ∞ h 2 max > 0.
or equivalently, using the explicit expression of c given in Proposition 2.2,

gH 0 + ( h + 2) E ∞ H 0 > (v ∞ ) 2 h3 with h = h max H 0 .
In the case where v ∞ > 0, the criterion given in the statement of Proposition 2.2 for the existence of solitary waves can therefore be restated as: left-going solitary waves always exist, but right-going solitary waves exist if and only if h max < h crit where the critical height h crit is given by h crit = H 0 hcrit with hcrit the only positive root of the polynomial

P (X) = (v ∞ ) 2 X 3 -gH 0 + (X + 2) E ∞ H 0 .
The second qualitative change with respect to the situation previously considered is that the shape of the two solitary waves of speed c = ±c are not the same since the ODE in Definition 2.1 does no longer depend on c through c 2 only. We refer to Figure 2.2 for an illustration of this behavior, in which we show the right and left-going waves profiles for increasing values of h max in the vicinity of h crit , for v ∞ = 1m.s -1 . We also highlight the influence of increasing values of v ∞ on the left-going waves profiles for a given value of h max .

We recall that we assumed in Proposition 2.2 that F ∞ = 0. Let us now give a brief discussion about the general case E ∞ > 0, v ∞ > 0 and F ∞ > 0 (one can treat the case v ∞ < 0 and/or F ∞ < 0 in a similar way). The presence of F ∞ = 0 implies that the possible speeds for solitary waves of maximal amplitude h max are found by solving the third order polynomial

X 3 + pX + q with p = -gh max + h max (h max + 2H 0 ) E ∞ H 3 0 (2.14) q = -h 2 max (h max + H 0 ) 2 F ∞ H 5 0 , (2.15)
(this is a simple consequence of (2.11)). Defining as previously ε by the discriminent ∆ = -(4p 3 + 27q 2 ) of this polynomial is always positive provided that the following smallness condition holds for F ∞ F 2 ∞ < 27 4

h max H 0 = 1 + ε -4 -3 -2 -1 0 1 2 3 4 1 1.5
H 5 0 (1 + ε)(2 + ε) 4 g + (3 + ε) E ∞ H 2 0 3 (2.16)
(this condition is satisfied for all realistic configurations). The polynomial (2.14) has then three distinct roots. Since the coefficient of X 2 is equal to zero, the sum of the three roots is necessarily equal to zero; moreover, their product has the sign of -q, and therefore the sign of F ∞ . If F ∞ > 0, then one has one positive root 0 < c + and two negative roots -c -,2 < -c -,1 < 0. There are therefore possibly two left going solitary waves, and a right going one. The right-going wave is subject to the same constraint h max < h crit as in the case E ∞ > 0, v ∞ > 0 and F ∞ = 0. In addition, (2.11) shows that the function

ϕ c : h → c 2 -h -ε 2 µh(h + 2)E ∞ -ε 3 µ 3/2 h 2 (h + 1) 2 F ∞ c (2.17)
must be positive for all H 0 ≤ h < h max . We can now state the following proposition where for the sake of simplicity, we considered only the case v ∞ > 0 and F ∞ > 0. The cases where these quantities are negative can be treated similarly.

Proposition 2.3. Let E ∞ > 0, v ∞ > 0 and F ∞ > 0.
Let also h max = H 0 (1 + ε) with ε > 0 and assume that (2.16) is satisfied. i. Up to translations, there can be at most three solitary waves of maximal height h max for (2.1) and corresponding speeds -c -,2 < -c -,1 < 0 < c + given by the roots of (2.14).

ii. The solitary wave of speed c + exists if and only if the following conditions hold

c + H 2 0 -v ∞ h 2 max > 0 and ∀h ∈ [H 0 ,h max ), ϕ c+ (h) > 0, (2.18)
where ϕ c is as defined in (2.17).

iii. The solitary wave of speed c -,j (j = 1,2) exists if and only ϕ -c-,j > 0 on [H 0 ,h max ).

Remark 2.4. Though there could be in principle a third solitary wave arising in the case F ∞ > 0, we could not exhibit any configuration where this is the case because the condition ϕ -c-,1 > 0 on [H 0 ,h max ) is never fulfilled. In practice, there are as in the case F ∞ = 0 one left going and one right going solitary wave, of different shape and of respective speed -c -,2 and c + . The profiles of the corresponding solitary waves are shown on Figure 2.3. Note that smaller values of F ∞ have to be taken to obtain the profiles of the left-going waves (Figure 2.3 (b)) in order to fulfill the condition ϕ -c-,2 > 0 on [H 0 ,h max ). Note also that additional solitary waves profiles can be observed, for both the constant vorticity model (1.5) and the general model (1.8), in §4.

Existence of peaked solitary waves

As seen in the previous section, when E ∞ > 0, v ∞ > 0 and F ∞ = 0, solitary waves have speed ±c, with

c = c(h max ) = gh max + h max (h max + 2H 0 ) H 3 0 E ∞ 1/2
and the maximal amplitude h max of the right-going solitary wave cannot exceed a critical value h crit corresponding to the only positive root of the equation

c(h crit )H 2 0 -v ∞ h 2 crit = 0.
Figure 2.2 suggests that the shape of the solitary waves tend to form an angle at their crest as their amplitude become close to the maximal amplitude. A byproduct of the analysis of the previous section is that there cannot exist any smooth solitary wave of maximal amplitude h crit . However, we show here that it is possible to obtain a peaked solitary wave of maximal amplitude in the following sense:

Definition 2.4. A peaked solitary wave of speed c, centered at x 0 ∈ R, for (2.1), is a mapping

(t,x) ∈ R 2 → (ζ,v,v ,E,F )(x -ct -x 0 ) such that there exists h ∈ C(R), with h | R + ∈ C 2 ([0,∞)), h | R -∈ C 2 ((-∞,0]), and
E ∞ > 0, v ∞ ∈ R and F ∞ ∈ R such that ζ = h -H 0 , v = c h -H 0 h , v = h H 0 v ∞ , (2.19)
E = h 3 H 3 0 E ∞ + 2 F ∞ c (h 2 -H 2 0 )h 3 H 5 0 , F = h 4 H 4 0 F ∞ , (2.20)
and such that h solves the ODE (2.10) on R + and R -, and satisfies lim ±∞ h = H 0 .

The proposition below proves the existence of peaked solitary waves in the case F ∞ = 0. Such a property could also be established for F ∞ = 0 (according to Proposition 2.3 there is also a critical maximal wave in some cases when F ∞ = 0), but the proof would be more technical and since no new phenomena arises in this case, we decide not to treat it. Some examples of peaked solitary waves are plotted on Figure 2.4, on which we highlight the influence of the value of v ∞ > 0 on the critical amplitude h crit . Proposition 2.5. Let E ∞ > 0, v ∞ > 0 and F ∞ = 0. For all x 0 ∈ R, there exists a unique peaked solitary wave centered at x 0 of critical maximal amplitude h crit and speed c = c(h crit ). It is even, decaying on both sides of the crest, and its angle at the crest is 2θ, with tanθ = 3 2

E ∞ cH 3 0 v ∞ (h crit -H 0 ) 2 1 + (v ∞ ) 2 E ∞ h 2 crit H 0 -1/2 .
Proof. We focus here on the case x ≥ 0; the case of negative values of x can be treated similarly. Without loss of generality, we also assume that x 0 = 0. By definition of h crit , one can write, for all h, and

cH 2 0 -v ∞ h 2 = -v ∞ (h -h crit )(h + h crit ) -4 -3 -2 -1 0 1 2 3 4 1 1.5
c 2 -gh - h(h + 2H 0 ) H 3 0 E ∞ = - E ∞ H 3 0 (h -h crit )(h + (v ∞ ) 2 E ∞ h 3 crit H 0 ).
On can therefore rewrite (2.11) under the form

c 3 v ∞ (h -h crit )(h + h crit )h 2 x = E ∞ H 3 0 (h -H 0 ) 2 (h -h crit )(h + (v ∞ ) 2 E ∞ h 3 crit H 0 ) or equivalently c 3 v ∞ (h + h crit )h 2 x = E ∞ H 3 0 (h -H 0 ) 2 (h + (v ∞ ) 2 E ∞ h 3 crit H 0 ).
Since h crit is by definition the maximal value of h and since h cannot reach the value H 0 (otherwise it would be identically equal to H 0 ), this ODE is equivalent to

h x = -3 E ∞ cH 3 0 v ∞ (h -H 0 ) 2 h + (v ∞ ) 2 E∞ h 3 crit H0 h + h crit 1/2 .
Existence of a local solution is therefore given by the standard Cauchy-Lipschitz theorem; the fact that the solution is global and tends to H 0 at infinity is then easily established as in Step 3 of the proof of Proposition 2.2. When evaluated at the origin (i.e replacing h by h crit ) in the above formula

Numerical method

We introduce now a simple numerical method to approximate the solutions of system (1.8) and illustrate the propagation of some of the various wave profiles exhibited in §2. This approach is inspired by some of our previous works [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF][START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF], and we give some details about the new ingredients introduced to account for the specificities of (1.8).

Splitting

To build our numerical method, denoting E = h 2 Ẽ et F = h 3 F and using the mass conservation equation, we rewrite system (1.8) in the equivalent form:

               h t + (hv) x = 0, (1 + T) (hv) t + (hv 2 ) x + ghζ x + (h 2 Ẽ) x + hQ 1 (v) + hC(v,v ) = 0, v t + (vv ) x = 0, Ẽt + (v Ẽ) x + 3 F h x + h Fx = 0, Ft + (v F ) x = 0. (3.1)
Then, in the spirit of [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF][START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2d simulations[END_REF], we decompose the solution operator S(•) associated to the formulation (3.1), at each time step by the second order splitting scheme

S(δ t ) = S 1 (δ t /2)S 2 (δ t )S 1 (δ t /2), (3.2)
where S 1 and S 2 are respectively associated to the transport part and dispersive perturbation of the Green-Naghdi equations (3.1). More precisely:

• S 1 (t) is the solution operator associated to the conservative propagation step

                   h t + (hv) x = 0, (hv) t + (hv 2 ) x + ghζ x + (h 2 Ẽ) x = 0, v t + (vv ) x = 0, Ẽt + (v Ẽ) x = 0, Ft + (v F ) x = 0. (3.3)
• S 2 (t) is the solution operator associated to the dispersive correction,

                   h t = 0, (hv) t -ghζ x -(h 2 Ẽ) x + (1 + T) -1 ghζ x + (h 2 Ẽ) x + hQ 1 (v) + hC(v,v ) = 0, v t = 0, Ẽt + 3 F h x + h Fx = 0, Ft = 0.
(3.4) As detailed in the next subsection, S 1 (t) is discretized using a finite-volume approach, and we use a finite-difference approach for S 2 (t).

Remark 3.1. The proposed splitting scheme is classically of second order accuracy in time, and it is shown in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF] that the corresponding semi-discretized dispersion relation approaches the exact dispersion relation of the Green-Naghdi at order 2 in time. Of course, other approaches can be used, like the unsplitted approaches proposed in [START_REF] Duran | Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations[END_REF] or [START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation[END_REF][START_REF] Popinet | A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations[END_REF] for instance.

We use the following notations in the following:

• the numerical one-dimensional domain Ω is uniformly divided into

N x cells (C i ) 1≤i≤Nx such that C i = [x i-1 2 ,x i+ 1 2 ]
, where (x i+ 1 2 ) 0≤i≤Nx are the N x + 1 nodes of the regular grid. We denote by x i the center of C i , • we denote by δ x the cell size (constant in this work) and by δ t the chosen time step (to be specified according to a relevant CFL-like condition), • we denote by wi n the averaged value of an arbitrary quantity w on the i th cell C i at time t n = nδ t .

Discretization of the conservative step

We focus on the discretization of system (3.3) which can be written in compact form as follows

∂ t W + ∂ x F(W) = S(W,b), (3.5)
with W = (h,hv,v , Ẽ, F ) and

F(W) = (hv, hv 2 + p(h, Ẽ), vv , v Ẽ, v F ), S(W,b) = (0, -ghb x , 0, 0, 0), (3.6)
with p(h, Ẽ) = g 2 h 2 + h 2 Ẽ. Neglecting the bottom variations, the study of the associated algebra shows that the system is hyperbolic with the following eigenvalues:

λ 1 = v -gh + 3h Ẽ, λ 2 = v + gh + 3h Ẽ, λ 3 = λ 4 = λ 5 = v.
(3.7)

1 st -order FV discretization of the homogeneous system

We first study a first order conservative spatial discretization of the homogeneous system associated with (3.5):

Wn+1 i -Wn i + δt δ x F Wn i , Wn i+1 -F Wn i-1 , Wn i = 0 (3.8)
where (u,v) → F(u,v) is a numerical flux function consistent with the physical flux w → F(w). For the numerical validations shown in §4, and considering the 3 waves structure of the algebra, we have implemented a HLLC-type Riemann solver, see for instance [START_REF] Batten | On the choice of wavespeeds for the HLLC Riemann solver[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF].

Discretization of the topography and high-order extension

The discretization of the topography source term occurring in (3.5) is done following the well-balanced approach for the Saint-Venant equations described in [START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF], allowing to preserve the motionless steady states corresponding to

ζ = 0, v = 0, v = 0, E = 0, F = 0.
(3.9)

Note that one of the main properties of this approach is that whenever the initial solver satisfies some classical stability properties, it yields a simple and fast well-balanced scheme that preserves the positivity of the water height.

As shown in previous studies [START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF][START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[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], the use of high-order schemes is mandatory for the study of dispersive water waves, to avoid as much as possible to pollute the dispersive properties of the model with some dispersive truncation errors associated with 2 nd order schemes. Based on discrete finite-volume cell averaging Wn i at time t n = nδt we use in this work 3 rd and 5 th -order accuracy WENO reconstructions, following [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF], together with the weight splitting method [START_REF] Shi | A technique of treating negative weights in WENO schemes[END_REF].

Spatial discretization of S 2 (•), time discretization and boundary conditions

Following the approach developed in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF], system (3.4) is discretized using 4 th finitedifferences. The resulting matrix, for the discretization of the linear operator T is the same as in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF]. As far as time discretization is concerned, we choose to use explicit methods. The systems corresponding to S 1 and S 2 are integrated in time using third or fourth-order SSP-Runge-Kutta scheme [START_REF] Gottlieb | Strong stability preserving high order time discretization methods[END_REF]. For the sake of simplicity, we only use periodic and Neumann boundary conditions here, adapting the ghost-cells methods detailed in [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF].

Remark 3.2. The whole numerical strategy can be straightforwardly applied to the constant vorticity model (1.5), the medium amplitude equations (1.12) and the reduced model with F = 0 (1.13), adapting the approximate Riemann solvers to the corresponding hyperbolic part.

Numerical validation

In this section, we use the analysis of solitary waves performed in Section 2 to validate our numerical scheme. Various kinds of smooth solitary waves exhibited in Section 2 are numerically observed in §4.1 and used to evaluate the convergence rate. As shown in §4.2, our code is accurate enough to capture also the extremal peaked solitary waves exhibited in §2.3. We then treat in §4.3 an example with a non trivial topography which allows us to show that vorticity may have a considerable influence on the shoaling phase. The specified convergence rates are obtained using a discrete L 2 -error.

Solitary waves propagation

In the next test cases of this subsection, we consider H 0 = 1 and we compute the propagation of several solitary waves in a computational domain of 200m long. We consider several set of values for the wave amplitude εH 0 and the triplet (E ∞ ,v ∞ ,F ∞ ), allowing to cover the various configurations detailed in §2. Unless stated otherwise, we use WENO3 reconstructions, a SSP-RK3 scheme and we set the CFL number to 0.8.

Constant vorticity model

We consider here the Green-Naghdi model (1.5) with constant vorticity. We consequently have curl U = (0,ω,0) T with ω(t,x,z) = ω 0 = cst. The analysis of solitary waves solutions is provided in Proposition 2.2. Indeed, the Green-Naghdi model with constant vorticity can be obtained from the general equations (1.8), setting E = 1 12 h 3 ω 2 0 , v = hω 0 , F = 0, and consequently neglecting the equations on Ẽ and v which are equivalent to the continuity equation. Solitary wave solutions for (1.5) belong to the situation depicted in Remark 2.3. In the present test case, we set ω 0 = 0.3s -1 and we study the propagation of a solitary wave of relative amplitude ε = 0.2, initially centred at x 0 = 25m. The initial water height h 0 (x) = h(0,x) is computed as a solution of equation (2.11) and the corresponding velocity is initialized as

v 0 (x) = c(1 - H 0 h 0 (x) ) with
where c is obtained from (2.13):

c 2 = gh max + h max (h max + 2H 0 ) 12 ω 2 0 . (4.1)
Remark 4.1. From Proposition 2.2, we observe that with the choice ω 0 > 0, we have existence of right-going solitary waves (propagating at speed c > 0) for all amplitudes a. On the contrary, the existence of left-going waves is ruled by the additional compatibility condition

cH 0 -ω 0 h 2 max > 0.
We only focus here on the right-going wave, and we show on Figure 4.1 the corresponding profile h 0 , together with the usual profile for the classical Green-Naghdi equations.

We also compute the propagation of this solitary wave on the time interval ]0,T ], with T = 50s and show the corresponding profiles at several locations along the channel. For these pictures, we set δx = 0.125m. Performing a numerical convergence analysis, we compute the L 2 errors at T = 3s for ζ and hv on a sequence of refined meshes and we observe mean orders of convergence, obtained with linear regressions, of 3.26 and 3.3 respectively for ζ and hv.

4.1.2. General vorticity model -Case E ∞ > 0, v ∞ > 0, F ∞ = 0
In this second case, we consider the propagation of a solitary wave non-constant vorticity in the case F ∞ = 0. The third order term F remains uniformly equal to 0 during the propagation, and we can therefore consider the reduced model (1.13). The existence of solitary waves solutions is still ruled by Proposition 2.2. There are 2 solitary waves of opposite velocities given by (2.13), but with different profiles. We choose here E ∞ = 0.2m 3 .s -2 , v ∞ = 1.6m.s -1 and again, we set successively ε = 0.2 and ε = 0.3. Such choices ensure the existence of both left and right-going waves, leading to the profiles shown on the left of For these computations, we have set δx = 0.07m. We observe the preservation of the initial profiles, together with a very low numerical dissipation for the considered time interval. Again, we perform a numerical convergence analysis on a sequence of refined meshes for both cases. In the case ε = 0.2, we observe some convergence rates of 2.9 for hv and 2.88 for ζ, Ẽ and v . Similar results are obtained in the case ε = 0.3, with slightly decreased rates. We also obtain a very similar behavior for the left going wave

4.1.3. General vorticity model -Case E ∞ > 0, v ∞ > 0, F ∞ > 0
We now focus on the case F ∞ > 0 and consider the general vorticity model (1.8). The existence of solitary wave solutions is now ruled by Proposition 2.3 and we have existence of 2 solitary waves, provided that conditions (2.18) are fulfilled. We also recall that the waves speeds are now obtained as the minimum and maximum roots of the 3 rd order polynomial (2.14). We consider the following set of parameters .s -3 and ε = 0.5, ensuring the existence of both left and right going waves and leading to the following wave speeds: c -≈ -3.852m.s -1 and c + ≈ 3.93m.s -1 . The corresponding profiles for ζ are shown on Figure 4.4 and Figure 4.5, together with some snapshots of the corresponding time evolutions. For these computations, we set δx = 0.05m. The corresponding convergence analysis is performed and we obtain convergence rates of 2.72 for both v , F , Ẽ and ζ and 2.78 for hv. We obtain the same behavior for the left going wave. Remark 4.2. As observed in the previous test cases, convergence rates classically seem to decrease with respect to the wave amplitude. Also, note that when F ∞ = 0, we observe a slight increase of the L 2 -errors for Ẽ when compared with the values obtained for v , F and ζ, and the corresponding convergence rates are slightly deteriorated. This phenomenon may be related to the naive finite-difference approximation of the nonlinear and non-conservative terms occurring in the corresponding source terms.

E ∞ = 1/12m 3 .s -2 , v ∞ = 1m.s -1 , F ∞ = 1/12m

Peaked solitary waves

Let us now briefly illustrate the propagation of peaked solitary waves, exhibited in Proposition 2.5. The numerical simulation of such waves with singularities is a difficult problem as the solution is not smooth. Indeed, the damping introduced by the numerical methods rapidly smooth out the peak and thereby deform and delay the wave. These issues may be avoided by locally refining the mesh and/or increasing the scheme's accuracy in the vicinity of the waves crest. The development of a more involved discrete formulation allowing for h/p-adaptivity, based on the recently introduced dG method for the GN equations [START_REF] Duran | Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations[END_REF], is left for future work.

To achieve this, we consider the simplified case F ∞ = 0, neglecting the third order tensor F , and still set H 0 = 1m. Choosing the following values v ∞ = 3m.s -1 and E ∞ = 9/12m 3 .s -2 , we set the relative amplitude in order to be very close to the critical case. In this particular configuration, we have existence of a right-propagating wave only for h max < h crit = 1.1045m. Consequently, we consider here a wave of amplitude 0.1044m, initially centered at x = 11m. We show the corresponding initial profile on the left side of Figure 4.6. The corresponding propagation is shown on the right side. For this simulation, in order to minimize the numerical diffusion in the vicinity of the singularity we increase the numerical resolution, setting δx = 2.510 -3 m and use WENO5 reconstructions with a SSP-RK4 time marching scheme.

Influence of vorticity on wave shoaling

In the following test case, we assess the topography terms discretization, and aim at giving a brief insight into the important study of the impact of the vorticity on wave shoaling. We still consider a 200m channel with H 0 = 1m but now, we the topography is varying and defined as follows:

b(x) = H 10 (1 + tanh x -x 1 λ ),
with λ = 20m, and x 1 = 100m. We follow successively the propagation of 3 solitary waves over this uneven bottom. The first one is a classical solitary wave (E ∞ = v ∞ = F ∞ = 0) associated with the irrotational Green-Naghdi equations. In this case, the model (3.1) and the associated numerical approach reduces to the extensively validated framework of [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF][START_REF] Bonneton | Recent advances in Serre-Green-Naghdi modelling for wave transformation, breaking and runup processes[END_REF]. This wave may therefore be used as a reference solution to highlight the influence of the additional vorticity terms.

The second and third ones are solitary wave solutions of (3.1) defined with F ∞ = 0 and (E ∞ ,v ∞ ) respectively set to (8.3310 -2 ,-1) and (0.18,-1.5). The corresponding celerities are given by c 0 ≈ 3.43m.s -1 for the irrotational wave and respectively c 0 ≈ 3.47m.s -1 and c 0 ≈ 3.53m.s -1 for the second and third vorticity waves. The corresponding initial profiles are shown on Figure 4.7-left. We also show on Figure 4.8 some snapshots of the waves evolution along the channel, and particularly as the waves prop- agates over the smooth step located in the vicinity of x 1 . The computation is performed with δx = 0.1m. We can observe the increasing influence of the vorticity on the waves shoaling processes, as the waves amplitude increase with respect to the vorticity magnitude, as well as on the stretching area at the rear of the waves and on the mean level at the front. A zoom on the profiles at t = 35s is shown on Figure 4.8-right. Of course, extensive studies are still needed to systematically analyze and accurately this influence.

Evolution of the velocity field inside the fluid domain

In the previous numerical computations we computed the surface elevation ζ, the average velocity v as well as v , E and F from their initial value through the resolution of the one-dimensional equations (1.8). We show here how this system of equations can be used to describe the dynamics of the (1 + 1)-dimensional velocity field U = (v,w) and of the surface elevation in terms of their initial value U 0 and ζ 0 .

The initial condition

Denoting by v 0 the vertical average of the horizontal velocity at t = 0 and by v * ,0 its fluctuation, one can write

v 0 = v 0 + v ,0 , w 0 = -∂ x (H 0 + z -b)v 0 -∂ x z -H0+b v ,0 ;
(5.1) the expression for w 0 is deduced from the expression for v 0 thanks to the incompressibility condition and the non penetration condition at the bottom. We also recall that the initial vorticity ω 0 is defined as

ω 0 = ∂ z u 0 -∂ x w 0 .
We shall consider here initial conditions that correspond to shallow water configurations (i.e. µ = H 2 0 /L 2 1). Denoting by v 0 sh the shear velocity induced by the vorticity field ω 0 and by v ,0 sh its fluctuation around its vertical mean value,

v ,0 sh = - ζ 0 z ω 0 + 1 h 0 ζ 0 -H0+b ζ 0 z ω 0 (h 0 = H 0 + ζ 0 -b), (5.2)
it is shown in [START_REF] Castro | Fully nonlinear long-waves models in the presence of vorticity[END_REF], Eq. (2.30), that up to O(µ 3/2 ) terms, one has

v * ,0 = v ,0 sh + T * v 0 (5.3)
where T * v 0 accounts for the fluctuations due to dispersion and is given by

T v 0 = - 1 2 (z + H 0 -b) 2 - 1 3 (h 0 ) 2 ∂ 2 x v + z -ζ 0 + 1 2 h 0 ∂ x b∂ x v 0 + ∂ x (v 0 ∂ x b) .
Since main goal here is to comment on the effects of the vorticity on the propagation of waves, it is convenient to consider an initial vorticity field ω 0 and to construct the corresponding approximate velocity field U 0 given by (5.1), (5.2) and (5.3).

Reconstruction methodology

Let us now briefly recall the procedure described in [START_REF] Castro | Fully nonlinear long-waves models in the presence of vorticity[END_REF] to recover U at all times from ζ 0 and U 0 . It is of course sufficient to compute U on each level line θ ∈ [0,1], that is,

v θ (t,x) = v t,x,-H 0 + b(x) + θh(t,x) , w θ (t,x) = w t,x,-H 0 + b(x) + θh(t,x) .
Remark 5.1. One can also define the vorticity on the level lines by ω θ (t,x) = ω t,x,-H 0 + b(x) + θh(t,x) . It can be computed in terms of v θ and w θ by the formula

ω θ = 1 h ∂ θ v θ -∂ x w θ + ∂ x (b + θh) h ∂ θ w θ .
It is shown in [START_REF] Castro | Fully nonlinear long-waves models in the presence of vorticity[END_REF] that the velocity fields conserves its initial structure (5.1)-(5.3) and that one has

v θ = v + v sh,θ + T θ v (5.4) w θ = -∂ x h(θv + Q θ + θ 0 T θ vdθ ) + ∂ x (-H 0 + b + θh)v θ . (5.5)
with

T θ v = - 1 2 (θ 2 - 1 3 )h 2 ∂ 2 x v + (θ - 1
2 )h ∂ x b∂ x v + ∂ x (∂ x bv) .
The computation of the different quantities involved in these expressions must be performed as follows:

1. From the initial datas ζ 0 , v 0 and v ,0 sh compute the initial values E 0 , F 0 and v ,0 from their definitions (1.9), (1.10) and (1.11). 2. Compute ζ, v, v , E and F on the time interval [0,T ] by solving the Green-Naghdi equations (1.8).

3. Compute the quantities q θ and Q θ on the same time interval by solving

∂ t q θ + ∂ x (vq θ ) = 0, ∂ t Q θ + ∂ x (vQ θ ) = 0
with initial conditions q 0 θ = ∂ θ v ,0 sh,θ and Q 0 θ (x) = θ 0 v ,0 sh,θ dθ , where we use the notation v ,0 sh,θ (x) = v ,0 sh t,x,-H 0 + b(x) + θh 0 (t,x) . 4. Compute v sh,θ by solving on [0,T ] the equation

∂ t v sh,θ + ∂ x v sh,θ (v + 1 2 v sh,θ ) = 1 h ∂ x E + q θ h ∂ x (hQ θ ),
with initial data v ,0 sh,θ . 5. Use (5.4) and (5.5) to get v θ and w θ .

Implementation hints

From a practical point of view, we aim at computing the velocity field U on some given N z level lines along the vertical layer of fluid. We introduce the discrete increments

θ j = δz + (j -1)δz H 0 , 1 ≤ j ≤ N z , with δz = H 0 N z ,
such that, for some given time t and horizontal coordinate x, the j th level line is located at the vertical coordinate z = -H 0 + b(x) + θ j h(t,x). To compute the velocity field on each of these level lines, we supplement the Green-Naghdi model with general vorticity (3.1) with a set of equations describing the time evolution of the auxiliary quantities q θj , Q θj and v θj :

                               ∂ t h + ∂ x (hv) = 0, (1 + T)[∂ t (hv) + ∂ x (hv 2 )] + gh∂ x ζ + hQ 1 (v) + ∂ x (h 2 Ẽ) + hC(v,v ) = 0, ∂ t v + ∂ x (vṽ ) = 0, ∂ t Ẽ + (v Ẽ) x + 3 F ∂ x h + h∂ x F = 0, ∂ t F + ∂ x (v F ) = 0, ∂ t q θj + ∂ x (vq θj ) = 0, 1 ≤ j ≤ N z , ∂ t Q θj + ∂ x (vQ θj ) = 0, 1 ≤ j ≤ N z , ∂ t v θj + ∂ x v θj (v + 1 2 v θj ) = 1 h ∂ x (h 2 Ẽ) + q θj h ∂ x (hQ θj ), 1 ≤ j ≤ N z . (5.6)
Considering the numerical method introduced in §3, the easiest way to account for these additional equations, from a discrete point of view, is to incorporate the auxiliary quantities q θj ,Q θj and v θj into our splitting approach, and into the approximate Riemann solver for the transport part of the equations.

Remark 5.2. As these new equations are decoupled from the equations on mass and momentum, we choose to keep the simple 3 waves structure of the approximate HLLC Riemann solver. Such a choice may appear as surprising for the quantities v θj , considering the particular form of the associated non-linear evolution equations. However, our numerical investigations have shown that considering each scalar equations separately, the corresponding wave speeds associated with the values of v + 1 2 v θj at interfaces are very close to v (which is in accordance with the vorticity scaling studied here, see [START_REF] Castro | Fully nonlinear long-waves models in the presence of vorticity[END_REF]), and always lie between the considered external wave speed estimates.

To simplify the notations, let us consider the case of one level line, located at the vertical coordinate z = -H 0 + b(x) + θh(t,x), with θ ∈]0,1], and denote by q θ ,Q θ and v θ the associated quantities. The splitting scheme (3.3)-(3.4) is modified as follows:

• S 1 (t) is the solution operator associated to the conservative propagation step, which may be written in the following conservative form

∂ t W θ + ∂ x F(W θ ) = S(W θ ,b), (5.7)
with

W θ = T (h, hv, v , Ẽ, F , q θ , Q θ , v θ ), F(W θ ) = T (hv, hv 2 + p(h, Ẽ), vv , v Ẽ, v F , vq θ , vQ θ , v θ (v + 1 2 v θ )),
and S(W θ ,b) = (0, -ghb x , 0, 0, 0, 0, 0, 0). • S 2 (t) is the solution operator associated to the dispersive correction,

                                     ∂ t h = 0, ∂ t (hv) -gh∂ x ζ -∂ x (h 2 Ẽ) +(1 + T) -1 gh∂ x ζ + ∂ x (h 2 Ẽ) + hQ 1 (v) + hC(v,v ) = 0, ∂ t v = 0, ∂ t Ẽ + 3 F ∂ x h + h∂ x F = 0, ∂ t F = 0, ∂ t q θ = 0, ∂ t Q θ = 0, ∂ t v θ - 1 h ∂ x (h 2 Ẽ) - q θ h ∂ x (hQ θ ) = 0. (5.8)
Note that when compared with (3.4), the second step is only modified with the introduction of a source term in the equation on v θ .

A solitary wave propagating in a constant vorticity field

To illustrate the previous reconstructions, we show here the velocity fields associated with some of the new solitary wave solutions exhibited in §2.1. We choose the simplest case of a constant vorticity field, described by model (1.5) and we set curl U = (0,ω,0) T with ω(t,x,z) = ω 0 . Using (5.2), we obtain the initial horizontal shear velocity field v ,0 sh (x,z) = -ω 0 (ζ 0 -z) + ω 0 2 h 0 , and the initial velocity field v 0 θ et w 0 θ are obtained with straightforward computations using (5.4) and (5.5). This initial velocity field is shown on Figure 5.1 for the values ω 0 = 0.3s -1 and ω 0 = 1.5s -1 . For comparison purpose, we also show the velocity field corresponding to the irrotational solitary wave of similar amplitude. We can observe the vertical decreasing of the velocity's horizontal component magnitude for the rotational solitary wave, while this is left unchanged for the irrotational wave.

A solitary wave arriving in a vorticity region

We consider now an initial vorticity field with vertical and horizontal dependencies, defined as follows: (this scaling of the vorticity corresponds to the regime studied in [START_REF] Castro | Well-posedness and shallow water stability for a new hamiltonian formulation of the water waves equations with vorticity[END_REF][START_REF] Castro | Fully nonlinear long-waves models in the presence of vorticity[END_REF] where the Green-Naghdi models with vorticity have been derived and justified). With this definition, we observe that the vorticity vanishes as x → ±∞. Additionally, its strongest variations of amplitude along the horizontal dimension are located near some point of abscissa x 1 . In the following, we aim at qualitatively observing the behavior of an irrotational solitary wave arriving from afar (where the vorticity is therefore negligible). We consider a channel of 50m long, set H 0 = 1m, x 1 = 25m and take as initial surface elevation the profile of an irrotational solitary wave (2.12) centred at x 0 = 12.5m with ε = 0.25. We highlight that x 0 is chosen far enough from x 1 to ensure that the vorticity in initially negligible in the vicinity of the solitary wave. The choice (5.9) for the initial vorticity field leads to the following initial horizontal shear velocity v ,0 sh (x,z) = ωH 0 sin 2π

ζ 0 -z H 0 - H 0 2π 1 -cos(2π h 0 H0 ) h 0 exp - (x -x 1 ) 2 λ 2
so that the full initial velocity field can again be constructed 4 through (5.1) (see Figure 5.2). We now follow the procedure described in the previous sections. From simple 4 The two integrals in the expression for w 0 can be computed explicitly,

z -H 0 v * ,0 sh = ωH 2 0 6π cos 2π ζ -z H 0 -cos 2π h 0 H 0 - z + H 0 h 0 (1 -cos(2π h 0 H 0 )) exp - (x -x 1 ) 2 λ 2 z -H 0 T v 0 = - 1 6 (z + H 0 ) (z + H 0 ) 2 -(h 0 ) 2 ∂ 2 x v 0 .
computations, one gets v ,0 = 2H 2 0 ω π 3 (h 0 ) 3 (h 0 ) 2 π 2 (1 + 2c 2 ) -3H 2 0 s 2 exp(-

(x -x 1 ) 2 λ 2 )
E 0 = 3H 2 0 ω 2 2π 2 h 0 H 0 h 0 π(s 2 -c 2 )sc -2H 2 0 s 4 + π 2 (h 0 ) 2 exp(-2 (x -x 1 ) 2 λ 2 )

F 0 = 3H 4 0 ω 3 2π 3 (h 0 ) 2 s 2 4π 2 (h 0 ) 2 (1 -2c 2 )c 2 -9πH 0 h 0 (1 -2c 2 )cs + 12H 2 0 (1 -c 2 ) 2 -5π 2 (h 0 ) 2 exp(-3 (x -x 1 ) 2 λ 2 )
where we used the notations s = sin(π ζ 0 H0 ) and c = cos(π ζ 0 H0 ). These functions are shown on Figure 5.3. Note in particular that the function F almost identically vanishes (its sup norm does not exceed 2.10 -9 ), and it is therefore possible to work with the reduced model (1.13). For the computation, we set N x = 500, leading to δx = 0.1m, and N z = 30. The corresponding results are shown on Figure 5.2, where we plot the free surfaces at several times during the propagation, together with the reconstructed velocity fields and the corresponding vorticity. We can observe the impact of the current on the initial velocity field and the free surface profile, and that the initial vorticity area is slightly dragged along and stretched as the wave propagates.

Influence of vorticity on shoaling

We consider here the same configuration as in §5.5 but with the vorticity region now located in an area where the topography varies. More precisely, instead of a flat bottom (b = 0) we now consider a bottom parametrized by z = -H 0 + b(x), with b(x) = H 0 6 (1 + tanh x -x 1 λ ).

One could derive as in §5.5 explicit expressions for the initial values of E, F and v , but they can also easily be numerically computed; they are represented in Figure 5.5. Note that in this configuration, F is still small (of order 10 -3 ) but not as much as in the previous configuration with a flat bottom. For this kind of configurations, it is therefore suitable to work with the full system (1.8). We show the corresponding results for the velocity fields on Figure 5.4.

Fig. 1 . 1 .

 11 Fig. 1.1. Main notations.

Fig. 2 . 1 .

 21 Fig. 2.1. Influence of E∞ on the solitary wave profile for H 0 = 1m, hmax = 2m, v ∞ = 0, F∞ = 0.

Fig. 2 . 2 .

 22 Fig. 2.2. Influence of v ∞ > 0 on the profiles for H 0 = 1m, E∞ = 1m 3 .s -2 and F∞ = 0. (a) Rightgoing solitary waves of nearly critical amplitudes for v ∞ = 1m.s -1 and therefore h crit ≈ 2.42m. (b) Left going solitary waves for the same amplitudes. (c) Right-going solitary waves for hmax = 2m and increasing values of v ∞ . (d) Left-going solitary waves for hmax = 2m and increasing values of v ∞ .

- 3 F =0.6m 4 .s - 3 F- 3 F = 1 m 4 .s - 3 F = 5 m 4 .s - 3 Fig. 2 . 3 .

 3433434323 Fig. 2.3. Influence of F∞ on the solitary waves profiles. (a) Shape of the left and right going solitary waves, and comparison with the usual Green-Naghdi solitary wave, for H 0 = 1m, hmax = 2m, E∞ = 1m 3 .s -2 , v∞ = 1m.s -1 , F∞ = 1m 4 .s -3 . (b) Influence of the value of F∞ on the left-going wave's profile. (c) Influence of the value of F∞ on the right-going wave's profile.

Fig. 2 . 4 .

 24 Fig. 2.4. Influence of v ∞ > 0 on the peaked solitary waves profiles for H 0 = 1m, E∞ = 1m 3 .s -2 and F∞ = 0.

Fig. 4 . 1 .

 41 Fig. 4.1. Propagation of a right-going solitary wave solution for the constant vorticity model -ε = 0.2 and ω 0 = 0.3s -1 : initial profile (on the left) and snapshots at different locations along the channel (on the right). The solitary wave solution of the classical Green-Naghdi equations is plotted in black solid line at t = 0s.

Figure 4 . 2 and

 42 Figure 4.3 respectively for the right-going and the left-going waves in the case ε = 0.3. The right-going wave is initially centred at x = 20m and the left-going wave at x = 170m. Again, we also plot the classical solitary wave profile of same amplitude for comparison purpose. The time evolutions of both right and left going waves are shown on the right of Figure 4.2 and Figure 4.3.

Fig. 4 . 2 .

 42 Fig. 4.2. Case E∞ > 0, v∞ > 0, F∞ = 0 -Right-going wave with ε = 0.3: initial profile (on the left) and snapshots at different locations along the channel (on the right).

Fig. 4 . 3 .

 43 Fig. 4.3. Case E∞ > 0, v∞ > 0, F∞ = 0 -Left-going wave with ε = 0.3: initial profile (on the left) and snapshots at different locations along the channel (on the right).

Fig. 4 . 4 .

 44 Fig. 4.4. Case E∞ > 0, v∞ > 0, F∞ > 0 -Right-going wave with ε = 0.5. Initial profile (on the left) and snapshots at different locations along the channel (on the right).

Fig. 4 . 5 .

 45 Fig. 4.5. Case E∞ > 0, v∞ > 0, F∞ > 0 -Left-going wave with ε = 0.5: initial profile (on the left) and snapshots at different locations along the channel (on the right).

Fig. 4 . 6 .

 46 Fig.4.6. Case E∞ > 0, v∞ > 0, F∞ = 0 -Right-going peaked wave with ε = 0.1044. We also show in black solid line the shape of the solitary wave for the irrotational Green-Naghdi equations.

Fig. 4 . 7 .

 47 Fig. 4.7.Influence of vorticity on wave shoaling: initial wave profiles (t = 0s) on the left and propagated profiles at t = 35s on the right. In the legend, 'Wave 1' refers to the irrotational solitary wave, while 'Wave 2' and 'Wave 3' respectively refers to the cases (E∞ = 8.3310 -2 ,v ∞ = -1) and (E∞ = 0.18,v ∞ = -1.5).

Fig. 4 . 8 .

 48 Fig. 4.8. Influence of vorticity on wave shoaling: evolution of the waves profiles along the channel at times t = 20,30 and 35s.

Fig. 5 . 1 .

 51 Fig. 5.1. Velocity fields in solitary waves with constant vorticity: influence of the vorticity magnitude.

  ω 0 (x,z) = -2πω cos 2π ζ -

Fig. 5 . 2 .

 52 Fig. 5.2. A solitary wave arriving on a vorticity region: snapshot of the free surface profile at times t = 0,2,3,4,4.5 and 5s. The velocity fields structures are plotted with black arrows and the corresponding vorticity is shown according to the left colorbar.

FFig. 5 . 3 .

 53 Fig. 5.3. Initial values for ζ, v and E, F , v for the configuration of §5.5 with H 0 = 1m and ε = 0.25.

Fig. 5 . 4 .

 54 Fig. 5.4. Influence of vorticity on shoaling: snapshot of the free surface profile at times t = 0,2,3,4,4.5 and 5s. The velocity fields structures are plotted with black arrows and the corresponding vorticity is shown according to the left colorbar.

FFig. 5 . 5 .

 55 Fig. 5.5. Initial values for ζ, v and E, F , v for the configuration of §5.6 with H 0 = 1m and ε = 0.25.

  3/2 = O(µ 2 ), where

	ε =	a surf H 0	=	Amplitude of the waves typical depth	,	β =	a bott H 0	=	Amplitude of the bottom variations typical depth	.

  4

Acknowledgments. D. Lannes acknowledges support from the ANR-13-BS01-0003-01 DYFICOLTI. F. Marche acknowledges support from the CNRS project LEFE-MANU SOLi. Both authors acknowledge support from the ANR-13-BS01-0009-01 BOND.