A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves

TL;DR

This paper introduces a reformulated Green-Naghdi model tailored for internal wave propagation in the Camassa-Holm regime, enhancing numerical resolution and dispersion properties, validated through a second order splitting scheme and simulations.

Contribution

The paper presents a new reformulation of the Green-Naghdi model that improves numerical stability and dispersion characteristics for internal wave simulations.

Findings

Enhanced frequency dispersion in the new model

Second order splitting scheme effectively solves the system

Numerical simulations validate the improved model and method

Abstract

In this paper we introduce a new reformulation of the Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves over a flat topography derived by Duch\^ene, Israwi and Talhouk. These new Green-Naghdi systems are adapted to improve the frequency dispersion of the original model, they share the same order of precision as the standard one but have an appropriate structure which makes them much more suitable for the numerical resolution. We develop a second order splitting scheme where the hyperbolic part of the system is treated with a high-order finite volume scheme and the dispersive part is treated with a finite difference approach. Numerical simulations are then performed to validate the model and the numerical methods.

Tables

Table 1. Table 1: Propagation of a solitary wave over a flat bottom: relative L 2 superscript 𝐿 2 \displaystyle L^{2} -error table for the conservative variables.

FV1-DF2-Euler

MUSCL-DF2-RK2

WENO5-DF4-RK4

Equations

$$U^{n+1}=S_1(\Delta t/2)S_2(\Delta t)S_1(\Delta t/2)U^n,$$

$${\mathfrak{T}}[\epsilon\zeta]V\ =\ q_{1}(\epsilon\zeta)V\ -\ \mu\nu\partial_{x}\Big{(}q_{2}(\epsilon\zeta)\partial_{x}V\Big{)}.$$

$$\displaystyle\Omega_{1}^{t}\$$

$$\displaystyle\Omega_{2}^{t}\$$

$$\left\{\begin{array}[]{ll}\partial_{x}^{2}\phi_{i}\ +\ \partial_{z}^{2}\phi_{i}\ =\ 0&\mbox{ in }\Omega^{t}_{i},\ i=1,2,\\ \partial_{t}\phi_{i}+\frac{1}{2}|\nabla_{x,z}\phi_{i}|^{2}=-\frac{P_{i}}{\rho_{i}}-gz&\mbox{ in }\Omega^{t}_{i},\ i=1,2,\\ \partial_{z}\phi_{1}\ =\ 0&\mbox{ on }\Gamma_{\rm t}\equiv\{(x,z),z=d_{1}\},\\ \partial_{t}\zeta\ =\ \sqrt{1+|\partial_{x}\zeta|^{2}}\partial_{n}\phi_{1}\ =\ \sqrt{1+|\partial_{x}\zeta|^{2}}\partial_{n}\phi_{2}&\mbox{ on }\Gamma\equiv\{(x,z),z=\zeta(t,x)\},\\ \partial_{z}\phi_{2}\ =\ 0&\mbox{ on }\Gamma_{\rm b}\equiv\{(x,z),z=-d_{2}\},\\ \lim\limits_{\varepsilon\to 0}\Big{(}P(t,x,\zeta(t,x)+\varepsilon)-P(t,x,\zeta(t,x)-\varepsilon)\Big{)}=-\sigma k(\zeta)&\mbox{ on }\Gamma,\end{array}\right.$$

$$c_0=\sqrt{g\frac{({\rho }_2-{\rho }_1)d_1{d}_2}{{\rho }_2{d}_1+{\rho }_1{d}_2}}.$$

$$\tilde{z}\equiv \frac{z}{d_1},\tilde{x}\equiv \frac{x}{\lambda },\tilde{t}\equiv \frac{c_0}{\lambda }t,$$

$$\tilde{\zeta }(\tilde{t},\tilde{x})\equiv \frac{\zeta (t,x)}{a},\widetilde{{\varphi }_i}(\tilde{t},\tilde{x},\tilde{z})\equiv \frac{d_1}{a\lambda c_0}{\varphi }_i(t,x,z)(i=1,2),$$

$$\gamma =\frac{{\rho }_1}{{\rho }_2},ϵ\equiv \frac{a}{d_1},\mu \equiv \frac{d_1^2}{{\lambda }^2},\delta \equiv \frac{d_1}{d_2},\mathrm{bo}=\frac{g({\rho }_2-{\rho }_1)d_1^2}{\sigma }.$$

$$G^{\mu }\psi$$

$$H^{\mu ,\delta }\psi$$

$$\displaystyle\displaystyle\left\{\begin{array}[]{ll}\left(\ \mu\partial_{x}^{2}\ +\ \partial_{z}^{2}\ \right)\phi_{1}=0&\mbox{ in }\Omega_{1}\equiv\{(x,z)\in\mathbb{R}^{2},\ \epsilon{\zeta}(x)<z<1\},\\ \partial_{z}\phi_{1}=0&\mbox{ on }\Gamma_{\rm t}\equiv\{(x,z)\in\mathbb{R}^{2},\ z=1\},\\ \phi_{1}=\psi&\mbox{ on }\Gamma\equiv\{(x,z)\in\mathbb{R}^{2},\ z=\epsilon\zeta\},\end{array}\right.$$

$$\displaystyle\displaystyle\left\{\begin{array}[]{ll}\left(\ \mu\partial_{x}^{2}\ +\ \partial_{z}^{2}\ \right)\phi_{2}=0&\mbox{ in }\Omega_{2}\equiv\{(x,z)\in\mathbb{R}^{2},\ -\frac{1}{\delta}<z<\epsilon\zeta\},\\ \partial_{n}\phi_{2}=\partial_{n}\phi_{1}&\mbox{ on }\Gamma,\\ \partial_{z}\phi_{2}=0&\mbox{ on }\Gamma_{\rm b}\equiv\{(x,z)\in\mathbb{R}^{2},\ z=-\frac{1}{\delta}\}.\end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle\partial_{t}{\zeta}\ -\ \frac{1}{\mu}G^{\mu}\psi\ =\ 0,\\ \\ \displaystyle\partial_{t}\Big{(}H^{\mu,\delta}\psi-\gamma\partial_{x}{\psi}\Big{)}\ +\ (\gamma+\delta)\partial_{x}{\zeta}\ +\ \frac{\epsilon}{2}\partial_{x}\Big{(}|H^{\mu,\delta}\psi|^{2}-\gamma|\partial_{x}{\psi}|^{2}\Big{)}\\ \displaystyle\hskip 142.26378pt=\ \mu\epsilon\partial_{x}\mathcal{N}^{\mu,\delta}-\mu\frac{\gamma+\delta}{\rm bo}\frac{\partial_{x}\big{(}k(\epsilon\sqrt{\mu}\zeta)\big{)}}{{\epsilon\sqrt{\mu}}}\ ,\end{array}\right.$$

$$\displaystyle\mathcal{N}^{\mu,\delta}\ \equiv\ \frac{\big{(}\frac{1}{\mu}G^{\mu}\psi+\epsilon(\partial_{x}{\zeta})H^{\mu,\delta}\psi\big{)}^{2}\ -\ \gamma\big{(}\frac{1}{\mu}G^{\mu}\psi+\epsilon(\partial_{x}{\zeta})(\partial_{x}{\psi})\big{)}^{2}}{2(1+\mu|\epsilon\partial_{x}{\zeta}|^{2})}.$$

$$\left\{\begin{array}[]{l}\partial_{t}\zeta+\partial_{x}\left(\dfrac{h_{1}h_{2}}{h_{1}+\gamma h_{2}}v\right)\ =\ 0,\\ \\ \mathfrak{T}[\epsilon\zeta]\left(\partial_{t}v+\epsilon\varsigma{v}\partial_{x}{v}\right)+(\gamma+\delta)q_{1}(\epsilon\zeta)\partial_{x}\zeta\\ \qquad\qquad+\frac{\epsilon}{2}q_{1}(\epsilon\zeta)\partial_{x}\left(\frac{h_{1}^{2}-\gamma h_{2}^{2}}{(h_{1}+\gamma h_{2})^{2}}|v|^{2}-\varsigma|v|^{2}\right)=-\mu\epsilon\frac{2}{3}\frac{1-\gamma}{(\gamma+\delta)^{2}}\partial_{x}\big{(}(\partial_{x}v)^{2}\big{)},\end{array}\right.$$

$$v\equiv \frac{1}{h_2}{\int }_{-\frac{1}{\delta }}^{ϵ\zeta (t,x)}{\partial }_x{\varphi }_2(t,x,z)dz-\frac{\gamma }{h_1}{\int }_{ϵ\zeta (t,x)}^1{\partial }_x{\varphi }_1(t,x,z)dz.$$

$${\mathfrak{T}}[\epsilon\zeta]V\ =\ q_{1}(\epsilon\zeta)V\ -\ \mu\nu\partial_{x}\Big{(}q_{2}(\epsilon\zeta)\partial_{x}V\Big{)},$$

$$\lambda =\frac{1+\gamma \delta }{3\delta (\gamma +\delta )},\alpha =\frac{1-\gamma }{(\gamma +\delta {)}^2}\text{mbox}and\beta =\frac{(1+\gamma \delta )({\delta }^2-\gamma )}{\delta (\gamma +\delta {)}^3}.$$

$$\nu =\lambda -\frac{1}{\mathrm{bo}}=\frac{1+\gamma \delta }{3\delta (\gamma +\delta )}-\frac{1}{\mathrm{bo}},$$

$$(\lambda -\frac{1}{\mathrm{bo}}){\kappa }_1=\frac{\gamma +\delta }{3}(2\beta -\alpha ),(\lambda -\frac{1}{\mathrm{bo}}){\kappa }_2=(\gamma +\delta )\beta ,$$

$$(\lambda -\frac{1}{\mathrm{bo}})\varsigma =\frac{2\alpha -\beta }{3}-\frac{1}{\mathrm{bo}}\frac{{\delta }^2-\gamma }{(\delta +\gamma {)}^2}.$$

$$f:X\to \frac{(1-X)({\delta }^{-1}+X)}{1-X+\gamma ({\delta }^{-1}+X)},$$

$$f(ϵ\zeta )=\frac{h_1{h}_2}{h_1+\gamma h_2}\text{mbox}and{f}^{\prime }(ϵ\zeta )=\frac{h_1^2-\gamma h_2^2}{(h_1+\gamma h_2{)}^2}.$$

$$\kappa=\frac{2}{3}\frac{1-\gamma}{(\delta+\gamma)^{2}}\quad\mbox{ and }\quad q_{3}(\epsilon\zeta)=\frac{1}{2}\big{(}f^{\prime}(\epsilon\zeta)-\varsigma\big{)}=\frac{1}{2}\Big{(}\frac{h_{1}^{2}-\gamma h_{2}^{2}}{(h_{1}+\gamma h_{2})^{2}}-\varsigma\Big{)}\ .$$

$$\left\{\begin{array}[]{l}\displaystyle\partial_{t}\zeta+\partial_{x}\big{(}f(\epsilon\zeta)v\big{)}\ =\ 0,\\ \\ \displaystyle{\mathfrak{T}}\left(\partial_{t}v+\epsilon\varsigma v\partial_{x}v\right)+(\gamma+\delta)q_{1}(\epsilon\zeta)\partial_{x}\zeta+\epsilon q_{1}(\epsilon\zeta)\partial_{x}(q_{3}(\epsilon\zeta){v}^{2})+\mu\epsilon\kappa\partial_{x}\big{(}(\partial_{x}v)^{2}\big{)}=0.\end{array}\right.$$

$$\mathcal{P}_{\rm SW}\ \equiv\ \Big{\{}(\mu,\epsilon,\delta,\gamma,\rm bo):\ 0\ <\ \mu\ \leq\ \mu_{\max},\ 0\ \leq\ \epsilon\ \leq\ 1,\ \delta\in(\delta_{\min},\delta_{\max}),\Big{.}\\ \Big{.}\ 0\ \leq\ \gamma\ <\ 1,\ \rm bo_{\min}\leq\rm bo\leq\infty\ \Big{\}},$$

$${\mathcal{P}}_{\mathrm{CH}}\equiv \left\{(\mu ,ϵ,\delta ,\gamma ,\mathrm{bo})\in {\mathcal{P}}_{\mathrm{SW}}:ϵ\le \mathrm{M}\sqrt{\mu }\text{mbox}\mathrm{and}\nu \equiv \frac{1+\gamma \delta }{3\delta (\gamma +\delta )}-\frac{1}{\mathrm{bo}}\ge {\nu }_0\right\},$$

$$\forall U=(\zeta ,v{)}^{\top }\in X^s,|U{|}_{X^s}^2\equiv |\zeta {|}_{H^s}^2+|v{|}_{H^s}^2+\mu |{\partial }_{x}v{|}_{H^s}^2.$$

$$\exists h_{01}>0\text{mbox}suchthat,\underset{x\in \mathbb{R}}{\inf }{h}_1\ge h_{01}>0,\underset{x\in \mathbb{R}}{\inf }{h}_2\ge h_{01}>0.$$

$$\exists h_{02}>0\text{mbox}suchthat,\underset{x\in \mathbb{R}}{\inf }\left(1+ϵ{\kappa }_1\zeta \right)\ge h_{02}>0,\underset{x\in \mathbb{R}}{\inf }\left(1+ϵ{\kappa }_2\zeta \right)\ge h_{02}>0.$$

$$\mathfrak{T}[ϵ\zeta ]V=(q_1(ϵ\zeta )I+\mu \nu T[ϵ\zeta ])V\text{mbox}withT[ϵ\zeta ]V=-{\partial }_x(q_2(ϵ\zeta ){\partial }_{x}V).$$

Full text

Key words : Green–Naghdi model, Nonlinear shallow water, Dispersive waves, Nonlinear interactions, Improved dispersion, Splitting method, Hybrid method, Finite volume, Finite Difference, High order scheme, WENO reconstruction

1. 1 Introduction

2. 2 Full Euler system

3. 3 Green-Naghdi model in the Camassa-Holm regime

4. 3.1 Reformulation of the model

5. 3.2 Improved Green-Naghdi equations

6. 3.3 Choice of the parameter $\displaystyle\alpha$

7. 3.3.1 The dispersion relation associated to the improved GN formulation

8. 3.3.2 The dispersion relation associated to the full-Euler system

9. 3.3.3 Phase velocity agreement

10. 3.4 High frequencies instabilities

11. 4 Numerical methods

12. 4.1 The splitting scheme

13. 4.2 Finite volume scheme

14. 4.2.1 First order finite-volume scheme

15. 4.2.2 Second order finite-volume scheme: MUSCL-RK2

16. 4.2.3 Higher order finite-volume scheme: WENO5-RK4

17. 4.3 Finite difference scheme for the dispersive part

18. 4.4 Boundary conditions

19. 5 Numerical validations

20. 5.1 Numerical validations in the one layer case

21. 5.1.1 Propagation of a solitary wave

22. 5.1.2 Head-on collision of counter-propagating waves

23. 5.1.3 Breaking of a Gaussian hump

24. 5.1.4 Dam-break problem in the one layer case

25. 5.2 Numerical validations in the two layers case

26. 5.2.1 Kelvin-Helmholtz instabilities

27. 5.2.2 Dam-break problem in the two layers case

28. 6 Conclusion

1 Introduction

This study deals with the propagation of internal waves in the uni-dimensional setting located at the interface between two layers of fluids of different densities. The fluids are assumed to be incompressible, homogeneous, and immiscible, limited from above by a rigid lid and from below by a flat bottom. This type of fluid dynamics problem is encountered by researchers in oceanography when they study the wave near the shore. Because of the difference in the salinity of the different layers of water near the shore, it is useful to model the flow of salted water by a two layers incompressible fluids flow. The usual way of describing such a flow is to use the 3D-Euler equations for the different layers adding some thermodynamic and dynamic conditions at the interface. This system will be called the full Euler system. This system of partial differential equations is very rich but very difficult to manipulate both mathematically and numerically. This is the reason why reduced models have been derived to characterize the evolution of the solution in some physical and geometrical specific regimes. Many models for the water wave (air-water interface) system have already been derived and studied in the shallow-water regime, where we consider that the wave length of the flow is large compared to the typical depth. We refer the reader to the following papers [2, 9, 30, 6, 29]. Earlier works have also set a very interesting theoretical background for the two-fluid system see [5, 1, 16, 17, 18], for more details.

One of these reduced models is the Green-Naghdi system of partial differential equations (denoted GN in the following). Many numerical studies of the fully nonlinear GN equations for the one layer case have been proposed in the literature. Let us introduce some of these results. The GN model describing dispersive shallow water waves have been numerically studied in [32] after being written in terms of potential variables in a pseudo-conservative form and using a Godunov scheme. Bonneton et al. studied numerically in [7] the fully nonlinear and weakly dispersive GN model for shallow water waves of large amplitude over variable topography. The original model was firstly adjusted in a suitable way for numerical resolution with an improvement of the dispersive properties, then they proposed to use a finite volume-finite difference splitting scheme that decomposes the hyperbolic and dispersive parts of the equations. This similar approach was earlier introduced in [22] for the solution of the Boussinesq equations. In [11], a three-parameters GN system is derived, yielding further improvements of the dispersive properties. Let us mention also [34], where a highly accurate and stable numerical scheme based on the Galerkin / finite-element method was presented for the Serre system. The method is based on smooth periodic splines in space and an explicit fourth-order Runge-Kutta method in time. Recently, Lannes and Marche introduced in [31] a new class of two-dimensional GN equations over varying topography. These fully nonlinear and weakly dispersive equations have a mathematical structure more suitable for 2D simulations. Using the same splitting strategy initiated in [11], they develop a high order, well balanced, and robust numerical code for these new models. Finally, we would like to mention the recent work of Duchêne, Israwi and Talhouk [19], where they derive a new class of GN models with improved frequency dispersion for the propagation of internal waves at the interface between two layers of fluids. They numerically compute their class of GN models using spectral methods [37] for space discretization and the Matlab solver ode45, which is based on the fourth and fifth order Runge-Kutta-Merson method for time evolution. Their numerical simulations show how the different frequency dispersion of the modified GN models may affect the appearance of high frequency Kelvin-Helmholtz instabilities.

In this paper, we present the numerical resolution of the GN model in the Camassa-Holm (medium amplitude) regime obtained and fully justified by Duchêne, Israwi and Talhouk in [18] using the same strategy initiated in [7, 31]. Let us recall that this model describes the propagation of one-dimensional internal waves at the interface between two layers of ideal fluids, limited by a rigid lid from above and a flat bottom from below. This model is first recast under a new formulation more suitable for numerical resolution with the same order of precision as the standard one but with improved frequency dispersion. More precisely, following [7] we derive a family of GN equations depending on a parameter $\alpha$ to be chosen in order to minimize the phase velocity error between the reduced model and the full Euler system. Then we propose a numerical scheme that decomposes the hyperbolic and dispersive parts of the equations. Following the same strategy adopted in [7, 31], we use a second order splitting scheme. The approximation $U^{n+1}=(\zeta^{n+1},v^{n+1})$ is computed at time $t^{n+1}=t^{n}+\Delta t$ where $\zeta$ represents the deformation of the interface and $v$ represent the shear mean velocity defined in Section 3, in terms of the approximation $U^{n}$ at time $t^{n}$ by solving

[TABLE]

where $S_{1}(.)$ is the operator associated to the hyperbolic part and $S_{2}(.)$ the operator associated to the dispersive part of the GN equations. For the numerical computation of $S_{1}(.)$, we use a finite volume method. We begin by the VFRoe method (see [10, 23, 24]), that is an approximate Godunov scheme. Unfortunately, the VFRoe scheme seems to be very diffusive. To this end, we propose a second-order scheme following the classical “MUSCL” approach [38]. Finally, following [28] we implement a fifth-order accuracy WENO reconstruction in order to reach higher order accuracy in smooth regions and a good resolution around discontinuities since the second order schemes are known to degenerate to first order accuracy at smooth extrema. On the other hand, $S_{2}(.)$ is computed using a finite difference scheme discretized using second and fourth order formulas, whereas for time discretization we use classical second and fourth order Runge-Kutta methods according to the order of the space derivative in consideration. The computation of $S_{2}(\Delta t)$ in the above splitting scheme (dispersive part) requires the inversion of the following symmetric second order differential operator introduced in [18]

[TABLE]

with $q_{i}(X)\equiv 1+\kappa_{i}X$ ($i=1,2$), where $\kappa_{i}$ and $\nu$ are constants depending on three parameters: the ratio between the densities of the two layers, the ratio between their depth and the capillary effect of the interface. Moreover, it is known that third order derivatives involved in this model may create high frequencies instabilities, but the presence of the operator ${\mathfrak{T}}[\epsilon\zeta]^{-1}$ in the second equation stabilizes the model with respect to these perturbations, allowing for more robust numerical computations, see Section 3.4. The invertibility of this operator played also an important role in the well-posedness of these equations (see [27] for the one layer case and [18] for the two layers case). However, this operator is time dependent since it depends on the deformation of the interface $\zeta(t,x)$. In order to avoid the inversion of this operator at each time step keeping its stabilizing effects, we derive a new family of physical models that are equivalent to the standard one in the sense of consistency (same order of precision $\mathcal{O}(\mu^{2})$), where we remove the time dependency of the operator $\mathfrak{T}$. The structure of the time-independent model leads to a small improvement in terms of computational time due to its simple one-dimensional structure. In fact, this strategy was originally initiated for numerical simulations of the fully nonlinear and weakly dispersive GN models in the two-dimensional case [31], in order to reduce significantly the computational time.

We organize the paper as follows. In Section 2, we introduce the non-dimensonalized full Euler system. Section 3 is devoted to recall the GN model in the Camassa-Holm regime, where an equivalent time-independent reformulation is given with an improved frequency dispersion based on the choice of a parameter $\alpha$. The stability issue of this new reformulation is discussed in the one and two layers cases. Section 4 is then devoted to the presentation of the numerical scheme. Firstly the hyperbolic/dispersive splitting is introduced. Then we describe the finite-volume spatial discretization of the hyperbolic part. The first-order finite volume scheme is given, an extension to the second-order accuracy is considered following the “MUSCL” approach and finally high-order accuracy around discontinuities is achieved thanks to the implementation of a fifth-order WENO reconstruction. The finite-difference discretization of the dispersive part is then detailed and the boundary conditions are briefly described. Finally, we present in Section 5 several numerical validations of our model. The one layer case is numerically validated: we compare the accuracy of the first, second and fifth orders of accuracy by studying the propagation of a solitary wave. To evaluate the influence of dispersive nonlinear terms, we consider the head-on collision of counter-propagating waves and the breaking of a Gaussian hump. Dealing with discontinuities is numerically validated by studying the dam-break problem in the one layer case. We highlight then the importance of the choice of the parameter $\alpha$ in improving the frequency dispersion of the model and compare our results with numerical experiments. Finally, we validate the ability of our numerical scheme to deal with discontinuities when considering a dam-break problem in the two layers case.

2 Full Euler system

In this section, we briefly recall the derivation of the full Euler system governing the evolution equations of the two-layers flow and refer to [1, 5, 16, 17] for more details. The two-layers flow considered are assumed to be incompressible, homogeneous, immiscible perfect fluids of different densities under the sole influence of gravity.

The study is restricted to the one-dimensional horizontal variable. The deformation of the interface between the two layers is represented by the graph of function a $\zeta(t,x)$ while the bottom and the top surfaces are assumed to be rigid and flat. The domains of the upper and lower fluid at time $t$ (denoted, respectively, $\Omega_{1}^{t}$ and $\Omega_{2}^{t}$), are given by

[TABLE]

In what follows, we assume that the domains remain strictly connected, that is there exists $h_{0}>0$ such that $d_{1}-\zeta(t,x)\geq h_{0}>0$ and $d_{2}+\zeta(t,x)\geq h_{0}>0$. The density and velocity fields of the upper and lower layers are denoted by $(\rho_{i},{\mathbf{v}}_{i})$, $(i=1,2)$ respectively. We assume the fluids to be incompressible, homogeneous and irrotational so that the velocity fields are divergence free and expressed as gradients of a potential denoted $\phi_{i}$. The fluids being ideal, that is with no viscosity, we may assume that each fluid follows the Euler equations. Assuming that the surface, the bottom or the interface are impenetrable one deduces the kinematic boundary conditions. Finally, the set of equations is completed by the continuity of the stress tensor at the interface.

Altogether, the governing equations are given by the following system:

[TABLE]

where $\displaystyle n$ denotes the unit upward normal vector at the interface.

The function \displaystyle k(\zeta)=-\partial_{x}\Big{(}\frac{1}{\sqrt{1+|\partial_{x}\zeta|^{2}}}\partial_{x}\zeta\Big{)} denotes the mean curvature of the interface and $\displaystyle\sigma$ the surface (or interfacial) tension coefficient.

To go further in the study of the flow of the two layers in order to build a numerical scheme for the asymptotic dynamics, we write the system in dimensionless form. To this end, we introduce dimensionless parameters and variables that reduce the setting to the physical regime under consideration. Firstly, let $\displaystyle a$ be the maximum amplitude of the deformation of the interface and $\displaystyle\lambda$ the wavelength of the interface. Then the typical velocity of small propagating internal waves (or wave celerity) is given by

[TABLE]

Consequently, we introduce the dimensionless variables:

[TABLE]

the dimensionless unknowns:

[TABLE]

and finally the dimensionless parameters:

[TABLE]

In what follows we assume that the depth ratio $\displaystyle\delta$ does not approach zero or infinity which means that the two layers of fluid have comparable depth. Therefore, the choice of the reference vertical length is harmless so we decided to choose $\displaystyle d_{1}$. We would like to mention that $\displaystyle\rm Bo=\mu\rm bo$ where $\displaystyle\rm Bo$ is the classical Bond number. The significations of the different dimensionless parameters are the following :

• •

$\displaystyle\gamma$ represents the ratio between the densities of the two layers,

• •

$\displaystyle\epsilon$ represents the amplitude of the interface between the two layers,

• •

$\displaystyle\mu$ represents the nonlinear effect of the shallowness approximation,

• •

$\displaystyle\delta$ represents the ratio between the depth of the two layers,

• •

$\displaystyle\rm bo$ represents the capillary effect of the interface.

Let us now remark that the system can be rewritten as two evolution equations coupling Zakharov’s canonical variables [40, 15], $\displaystyle(\zeta,\psi)$ representing respectively the deformation of the interface and the trace of the dimensionless upper potential at the interface defined by $\displaystyle\psi\ \equiv\ \phi_{1}(t,x,\zeta(t,x)).$

In order to do so, we define the so-called Dirichlet-Neumann operators. The tildes are removed for the sake of readability.

Definition 2.2** **(Dirichlet-Neumann operators)

Let $\displaystyle\zeta\in H^{t_{0}+1}(\mathbb{R})$, $\displaystyle t_{0}>1/2$, such that there exists $\displaystyle h>0$ with $\displaystyle h_{1}\ \equiv\ 1-\epsilon\zeta\geq h>0$ and $\displaystyle h_{2}\ \equiv\ \frac{1}{\delta}+\epsilon\zeta\geq h>0$, and let $\displaystyle\psi\in L^{2}_{\rm loc}(\mathbb{R}),\partial_{x}\psi\in H^{1/2}(\mathbb{R})$. Then we define:

[TABLE]

where, $\displaystyle\phi_{1}$ and $\displaystyle\phi_{2}$ are uniquely deduced from $\displaystyle(\zeta,\psi)$ as solutions of the following Laplace’s problems:

[TABLE]

At this stage of the model, using the above definition and withouth making any assumption on the different dimensionless parameters, one can rewrite the nondimensionalized version of (2.1) as a system of two coupled evolution equations, which writes:

[TABLE]

where we denote:

[TABLE]

We will refer to (2.5) as the full Euler system.

3 Green-Naghdi model in the Camassa-Holm regime

We now recall the new Green-Naghdi model in the Camassa-Holm regime recently derived by Duchêne, Israwi and Talhouk in [18]. This new model is derived after expanding the different operators of the original Green-Naghdi model with respect to $\displaystyle\epsilon\mbox{ and }\mu$. Then several additional transformations are made using the smallness assumption $\displaystyle\epsilon=\mathcal{O}(\sqrt{\mu})$ in order to produce an equivalent precise system of partial differential equations whose unknowns are $\displaystyle\zeta\mbox{ and }v$ :

[TABLE]

where $\displaystyle h_{1}\equiv 1-\epsilon\zeta$ (resp. $\displaystyle h_{2}\equiv\frac{1}{\delta}+\epsilon\zeta$) denotes the depth of the upper (resp. lower) fluid, and $\displaystyle v$ is the shear mean velocity defined by:

[TABLE]

The operator $\displaystyle{\mathfrak{T}}$ is defined as:

[TABLE]

with $\displaystyle q_{i}(X)\equiv 1+\kappa_{i}X$ ($\displaystyle i=1,2$) and $\displaystyle\nu$, $\displaystyle\kappa_{1}$, $\displaystyle\kappa_{2}$, $\displaystyle\varsigma$ are defined below. Let us first introduce the following constants in order to ease the reading:

[TABLE]

Thus

[TABLE]

[TABLE]

[TABLE]

The function $\displaystyle f$ is defined as follows:

[TABLE]

so one has:

[TABLE]

Additionally, let us denote:

[TABLE]

Problem (3.1) writes now in compact form:

[TABLE]

Let us now recall the regime of validity of the system (3.4) as exactly given in [18]. In the first place, we consider the so-called shallow water regime for two layers of comparable depths:

[TABLE]

with given $\displaystyle 0\leq\mu_{\max},\delta_{\min}^{-1},\delta_{\max},\rm bo_{\min}^{-1}<\infty$.

Two additional key restrictions are necessary for the validity of the model:

[TABLE]

with given $\displaystyle 0\leq M,\nu_{0}^{-1}<\infty$.

Duchêne, Israwi and Talhouk proved in [18], that the system (3.4) is well-posed (in the sense of Hadamard) in the energy space $\displaystyle X^{s}\ =\ H^{s}(\mathbb{R})\times H^{s+1}(\mathbb{R})$, endowed with the norm

[TABLE]

This result is obtained for parameters in the Camassa-Holm regime (3.6) and under the following non-zero depth and ellipticity (for the operator $\displaystyle\mathfrak{T}$) conditions:

[TABLE]

[TABLE]

They also prove that this new asymptotic model model is fully justified by a convergence result in the Camassa-Holm regime.

3.1 Reformulation of the model

One can easily check that the operator $\displaystyle\mathfrak{T}$ can be written under the form:

[TABLE]

Therefore, this operator is time dependant through the presence of the term $\displaystyle\zeta$ and at each time, this operator has to be inverted in order to solve equation (3.4) Thus we want to derive a new model equivalent to (3.1) up to a lower order in $\displaystyle\epsilon\mbox{ or }\mu$ but with a structure adapted to construct a numerical scheme for its resolution.

Let us first denote :

[TABLE]

such that one can rewrite problem (3.4) as :

[TABLE]

As we said before since the operator $\displaystyle(q_{1}(\epsilon\zeta)I+\mu\nu T[\epsilon\zeta])$ depends on $\displaystyle\zeta(t,x)$, one has to remove this time dependency to avoid its inversion at each time.

Firstly, let us write the operator $\displaystyle T[\epsilon\zeta]$ under the form:

[TABLE]

with

[TABLE]

Expanding the term \displaystyle\Big{(}q_{1}(\epsilon\zeta)I+\mu\nu T[\epsilon\zeta]\Big{)}\Big{(}\partial_{t}v+\epsilon\varsigma v\partial_{x}v\Big{)} in problem (3.8) leads to :

[TABLE]

But we have: $\displaystyle\mu\epsilon\nu S[\zeta]\left(\partial_{t}v+\epsilon\varsigma v\partial_{x}v\right)=\mu\epsilon\nu S[\zeta]\partial_{t}v+\mathcal{O}(\mu\epsilon^{2})$. By assumption the term $\displaystyle\mathcal{O}(\mu\epsilon^{2})$ is of order $\displaystyle\mathcal{O}(\mu^{2})$; thus Problem (3.9) is equivalent to the following system:

[TABLE]

3.2 Improved Green-Naghdi equations

As done by many authors, see [39, 33, 13], the frequency dispersion of problem (3.10) can be improved by adding some terms of order $\displaystyle\mathcal{O}(\mu^{2})$ to the momentum equation. The accuracy of the model is not affected by this manipulation since this equation is precise up to terms of the same order as the added ones, namely $\displaystyle\mathcal{O}(\mu^{2})$.

Going back to problem (3.8), one notices that:

[TABLE]

We can now divide (3.11) by $\displaystyle q_{1}(\epsilon\zeta)$ to obtain:

[TABLE]

Let us fix $\displaystyle\alpha\in\mathbb{R}$. Multiplying (3.12) by $\displaystyle(1-\alpha)$ leads to:

[TABLE]

Replacing $\displaystyle\partial_{t}v$ by its expression given in (3.13) in the second equation of (3.10) and regrouping the $\displaystyle\mathcal{O}(\mu^{2})$ terms yields to the following equation:

[TABLE]

After straightforward computations,

[TABLE]

Replacing in (3.14), $\displaystyle\epsilon\kappa_{1}\zeta\left(\partial_{t}v+\epsilon\varsigma v\partial_{x}v\right)$ by \displaystyle\epsilon\kappa_{1}\zeta\Big{(}-(\gamma+\delta)\partial_{x}\zeta-\epsilon\partial_{x}(q_{3}(\epsilon\zeta){v}^{2})-\dfrac{\mathcal{O}(\mu)}{q_{1}(\epsilon\zeta)}\Big{)}, we have:

[TABLE]

One deduces easily that,

[TABLE]

Plugging (3.16) in (3.15) one has:

[TABLE]

Remembering that $\displaystyle\mu\nu T[0](\partial_{t}v+\epsilon\varsigma v\partial_{x}v)+\mu\epsilon Q_{1}(v)+\mu\epsilon\nu S[\zeta]\partial_{t}v=\mathcal{O}(\mu)$, and regrouping the $\displaystyle\mathcal{O}(\mu\epsilon^{2})$ terms, (3.17) becomes:

[TABLE]

Expanding in terms of $\displaystyle\epsilon$ and $\displaystyle\mu$, and looking for the term $\displaystyle\partial_{t}v$, we get:

[TABLE]

Replacing $\displaystyle\partial_{t}v$ by its expression obtained in (3.19) in the last two terms of the equation (3.18) the Green–Naghdi equations with improved dispersion can therefore be written as (up to the order $\displaystyle\mathcal{O}(\mu^{2})$) :

[TABLE]

with

[TABLE]

and

[TABLE]

In [11], a three-parameter family of Green–Naghdi equations in the one layer case is derived yielding additional improvements of the dispersive properties. In this paper we stick to the one-parameter family (3.20). This new formulation does not contain any third-order derivative, thus one can expect more stable and robust numerical computations.

3.3 Choice of the parameter $\displaystyle\alpha$

The choice of the parameter $\displaystyle\alpha$ is motivated by the agreement of the dispersion properties of the full Euler system and the improved Green-Naghdi system (3.20) in term of the dispersion relation. Therefore, we have first to find the dispersion relation for the improved Green-Naghdi system, then the dispersion relation for the full Euler system and finally to find an optimal parameter $\displaystyle\alpha_{opt}$ to ensure a good agreement between these two relations.

3.3.1 The dispersion relation associated to the improved GN formulation

The system (3.20) can be written under the following form:

[TABLE]

where the operators $\displaystyle Q_{1}$, $\displaystyle Q_{2}$ and $\displaystyle Q_{3}$ are explicitly given in (3.7), (3.21) and (3.22). Looking at the linearization of (3.23) around the rest state $\displaystyle(\zeta,v)=(0,0)$, one derives the dispersion relation associated to (3.23).

This relation is obtained by looking for plane wave solutions of the form $\displaystyle(\zeta,v)=(\underline{\zeta},\underline{v})e^{i(kx-wt)}$, with $\displaystyle k$ the spatial wave number and $\displaystyle\omega$ the time pulsation, of the linearized equations.

The first equation of (3.23) gives:

[TABLE]

Thus we first obtain:

[TABLE]

From the second equation of (3.23) we have:

[TABLE]

This equation may be written as:

[TABLE]

Replacing $\displaystyle\underline{v}$ in (3.25) by its expression given in (3.24) we obtain:

[TABLE]

After straightforward computations, we obtain the following dispersion relation of the new GN formulation for $\displaystyle\alpha\geq 1$:

[TABLE]

Defining the linear phase velocity associated to (3.27) as:

[TABLE]

we choose $\displaystyle\alpha$ such that the phase velocity stays close to the reference phase velocity $\displaystyle C_{S}^{p}(k)$ coming from Stokes linear theory. Thanks to this approach, the error on the phase velocity is minimized for any discrete value of $\displaystyle\mu|k|$ and the corresponding local optimal value of $\displaystyle\alpha$, denoted by $\displaystyle\alpha_{opt}$ is computed. In order to do so, we have firstly to obtain the dispersion relation of the original full Euler system.

3.3.2 The dispersion relation associated to the full-Euler system

Let us recall the full-Euler system given in (2.5):

[TABLE]

where

[TABLE]

and

[TABLE]

where $\displaystyle\phi_{1}$ and $\displaystyle\phi_{2}$ are uniquely defined (up to a constant for $\displaystyle\phi_{2}$) as the solutions in $\displaystyle H^{2}(\mathbb{R})$ of the following Laplace’s problems.

[TABLE]

To obtain the dispersion relation associated to the full-Euler system, we first have to linearize it around the rest state that is $\displaystyle\zeta=0\,,\,v=0$. We have thus to compute the two operators $\displaystyle G^{\mu}$ and $\displaystyle H^{\mu,\delta}$ at the rest state. Therefore we have firstly to linearize the two previous equations (3.32)-(3.36), then write the linear system in the wave number space by performing a Fourier transform. After performing some algebraic computations, we will obtain the dispersion relation. Let us remark that the two equations (3.32) and (3.36) are very similar except for the boundary conditions, and the space domain $\displaystyle\Omega_{1}$ and $\displaystyle\Omega_{2}$. Therefore we have chosen to perform the complete analysis for the two equations.

$\displaystyle\bullet$ Linearizing (3.32) around the rest state ($\displaystyle\zeta=0$) gives:

[TABLE]

Applying the Fourier transform with respect to $\displaystyle x$, one has:

[TABLE]

One can easily remark that (3.44) is an ordinary linear differential equation of order 2 whose solution has the following form (since $\displaystyle\mu k^{2}\geq 0$):

[TABLE]

Writing the boundary condition at $\displaystyle z=0$ i.e. $\displaystyle\widehat{\phi_{1}}|_{z=0}=\widehat{\psi}$ and at $\displaystyle z=1$ i.e. $\displaystyle(\partial_{z}\widehat{\phi_{1}})|_{z=1}=0$ we obtain:

[TABLE]

Thus we obtain:

[TABLE]

Therefore as:

[TABLE]

we finally have:

[TABLE]

$\displaystyle\bullet$ Linearizing (3.36) around the rest state ($\displaystyle\zeta=0$) gives:

[TABLE]

Applying the Fourier transform with respect to $\displaystyle x$, one has:

[TABLE]

One can easily remark that (3.53) is an ordinary differential equation of order 2 whose solution has the following form (since $\displaystyle\mu k^{2}\geq 0$):

[TABLE]

Writing the boundary condition at $\displaystyle z=0$ i.e. $\displaystyle(\partial_{z}\widehat{\phi_{2}})|_{z=0}=-\sqrt{\mu}|k|\tanh(\sqrt{\mu}|k|)\widehat{\psi(k)}$ and at $\displaystyle z=-1/\delta$ i.e. $\displaystyle(\partial_{z}\widehat{\phi_{2}})|_{z=-1/\delta}=0$ we obtain:

[TABLE]

Replacing $\displaystyle B(k)$ by its expression, we obtain:

[TABLE]

Since $\displaystyle\widehat{\phi_{2}}(k)|_{z=0}=A(k)$ we have thus: \displaystyle\phi_{2}(k)|_{z=0}=-\dfrac{\tanh\Big{(}\sqrt{\mu}|k|\Big{)}}{\tanh\Big{(}\dfrac{\sqrt{\mu}|k|}{\delta}\Big{)}}\psi(k).

Therefore,

[TABLE]

Having obtained the two operators $\displaystyle G^{\mu}$ and $\displaystyle H^{\mu,\delta}$ at the rest state, (3.45) and (3.54), we can now proceed to compute the dispersion relation for the full Euler system.

One derives the dispersion relation associated to (3.28) by looking for plane wave solutions of the form $\displaystyle(\zeta,\psi)=(\underline{\zeta},\underline{\psi})e^{i(kx-wt)}$ to the linearized equations around the rest state $\displaystyle(\zeta,\psi)=(0,0)$.

The linearisation of the full Euler system, (3.28) at the rest state, $\displaystyle\zeta=0\,,\,\psi=0$, writes under the form:

[TABLE]

Replacing $\displaystyle G^{\mu}[0]\psi$ by its expression given by (3.45) in the first equation of the system (3.55), we firstly have:

[TABLE]

Thus,

[TABLE]

Replacing $\displaystyle H^{\mu}[0]\psi$ by its expression given by (3.54) in the second equation of the system (3.55), we obtain:

[TABLE]

Thus,

[TABLE]

After straightforward computations and replacing $\displaystyle\underline{\zeta}$ by its expression given in (3.56) we finally obtain the following dispersion relation:

[TABLE]

3.3.3 Phase velocity agreement

One can easily remark that the Taylor expansions of the two previous dispersion relations (3.27) and (3.57) are equivalent for small wavenumbers,

[TABLE]

Indeed, for small wavenumber one has:

[TABLE]

and

[TABLE]

Therefore for small wavenumbers and in the whole range of regime, the choice of $\displaystyle\alpha$ does not influence the dispersion relation (3.27). In the desired simulations, we are interested in small wave length (i.e large wavenumbers) dispersion characteristics and thus we would like to find an optimal value of $\displaystyle\alpha$ in order to observe the same dispersion properties for the reduced Green-Naghdi model that the full Euler original problem satisfies. Therefore, we are interested in values of $\displaystyle\alpha$ such that $\displaystyle C_{GN}^{p}(k)=C_{F.E}^{p}(k)$ for large value of $\displaystyle k$. We can thus compute $\displaystyle\alpha_{opt}$ from the equality: $\displaystyle w_{\alpha,GN}^{2}=w^{2}_{F.E}$.

Let us denote $\displaystyle X=\sqrt{\mu}|k|$. The previous equality writes:

[TABLE]

After straightforward computations we obtain the follwing expression for the value of $\displaystyle\alpha_{opt}$:

[TABLE]

In Figure 2 (top), $\displaystyle\alpha_{opt}$ is plotted against the spatial wave number $\displaystyle k$, for $\displaystyle k\in[0,4]$ and $\displaystyle\mu=1$, for the one and two layers cases. For the one layer case we set $\displaystyle\gamma=0$, $\displaystyle\delta=1$ and $\displaystyle\rm bo^{-1}=0$ whereas for the two layers case we set $\displaystyle\gamma=0.95$, $\displaystyle\delta=0.5$ and $\displaystyle\rm bo^{-1}=5\times 10^{-5}$. We would like to mention that $\displaystyle\alpha_{opt}\rightarrow 1$ in both cases when considering very large values of $\displaystyle k$, this fact is confirmed by the numerical simulation done in Section 5.2.2, see Figure 14, where we notice similar observations at time $\displaystyle t=20\ s$. Therefore, we will use an optimal value of $\displaystyle\alpha$ different of $\displaystyle 1$ only when looking at the dispersion effects for intermediate wave numbers. This is highlighted, when comparing to the numerical experiments done in [19] (for more details see Section 5.2.1). In Figure 2 (bottom), the ratio $\displaystyle\dfrac{C^{p}_{GN}(k)}{C^{p}_{S}(k)}$ (i.e linear phase velocity error), is plotted against the spatial wave number $\displaystyle k$ for $\displaystyle k\in[0,4]$, in the one and two layer cases.

3.4 High frequencies instabilities

In what follows, a qualitative discussion on the stability of the improved model (3.20) is given. We study the high frequencies instabilities in the one layer case, as well as in the two layers case. As $\displaystyle\mu<<1$, a simple expansion of the operator $\displaystyle(I+\mu\nu\alpha T[0])^{-1}$ shows that we have:

[TABLE]

Thus, we could replace $\displaystyle Q_{2}(\zeta)$ by \displaystyle\widetilde{Q}_{2}(\zeta)=-S[\zeta]\big{(}(\gamma+\delta)\partial_{x}\zeta\big{)}, and $\displaystyle Q_{3}(\zeta)$ by $\displaystyle\widetilde{Q}_{3}(\zeta)=\kappa_{1}\zeta T[0][(\gamma+\delta)\partial_{x}\zeta]$ in the second equation of (3.20), keeping the same order of precision $\displaystyle\mathcal{O}(\mu^{2})$. This simple expansion would avoid the inversion of $\displaystyle(I+\mu\nu\alpha T[0])$ (resolution of an extra linear system) in the computation of $\displaystyle{Q}_{2}$ and $\displaystyle{Q}_{3}$ but it leads to instabilities. Indeed, this is due to the fact that the two terms $\displaystyle\widetilde{Q}_{2}(\zeta)$ and $\displaystyle\widetilde{Q}_{3}(\zeta)$ contain third order derivatives in $\displaystyle\zeta$ that may create high frequencies instabilities.

Remark 3.58

In order to recover the dimensionalized version of system (3.20), one has to set $\displaystyle\mu=\epsilon=1$ and add the acceleration of gravity term $\displaystyle g$ when necessary to obtain the following system:

[TABLE]

with

[TABLE]

and

[TABLE]

Firstly, we discuss the stability of the model (3.59) in the one layer case without surface tension. To this end, we set $\displaystyle\gamma=0$, $\displaystyle\delta=1$ and $\displaystyle\dfrac{1}{\rm bo}=0$. Thus system (3.59) becomes:

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

When linearizing system (3.62) around constant state solution $\displaystyle(\underline{\zeta},\underline{v})$, one obtains the following linear system in $\displaystyle(\tilde{\zeta},\tilde{v})$:

[TABLE]

Looking for plane wave solution of the form $\displaystyle(\tilde{\zeta},\tilde{v})=e^{i(kx-wt)}(\zeta^{0},v^{0})$ as solution of the above system, one obtains the following dispersion relation:

[TABLE]

When we consider the linearization of our new model around the rest state $\displaystyle(\underline{\zeta},\underline{v})=(0,0)$, the dispersion relation (3.66) becomes:

[TABLE]

In this case the perturbations are always stable if $\displaystyle\alpha\geq 1$. We refer to Section 3.3 for the discussion concerning the choice of $\displaystyle\alpha$ in order to improve the dispersive properties of the model.

Now, replacing $\displaystyle Q_{2}(\zeta)$ defined in (3.63) by \displaystyle\widetilde{Q}_{2}(\zeta)=\partial_{x}\Big{(}\zeta\partial_{x}\big{(}g\partial_{x}\zeta\big{)}\Big{)} and $\displaystyle Q_{3}(\zeta)$ defined in (3.64) by $\displaystyle\widetilde{Q}_{3}(\zeta)=\kappa_{1}\zeta T[0][g\partial_{x}\zeta]$ in the second equation of (3.62), modifies the dispersion relation (3.66) and we obtain instead:

[TABLE]

One can easily check that if $\displaystyle\underline{\zeta}>\frac{(\alpha-1)}{2}$, the numerator of the right hand side of the dispersion relation (3.67) will become negative for sufficiently large values of $\displaystyle k^{2}$. Thus the root of $\displaystyle\tilde{w}$ will become complex inducing a high frequency instability of the model, see Figure 3. On the other hand, the numerator of the r.h.s of the dispersion relation (3.66) is positive for sufficiently large values of $\displaystyle k^{2}$, due to the existence of the term $\displaystyle-\dfrac{2k^{2}\underline{\zeta}}{3(1+\frac{\alpha}{3}k^{2})}$, thus $\displaystyle w$ is always real at high frequencies.

Let us now discuss the stability issue of the model (3.59) in the two layers case without surface tension. When linearizing system (3.59) around motionless steady state solution $\displaystyle(\underline{\zeta}=\text{cst},\underline{v}=0)$ and after following the same method as above, one obtains the following dispersion relation:

[TABLE]

Replacing $\displaystyle Q_{2}(\zeta)$ defined in (3.60) by \displaystyle\widetilde{Q}_{2}(\zeta)=\partial_{x}\Big{(}\zeta\partial_{x}\big{(}g(\gamma+\delta)\partial_{x}\zeta\big{)}\Big{)} and $\displaystyle Q_{3}(\zeta)$ defined in (3.61) by $\displaystyle\widetilde{Q}_{3}(\zeta)=\kappa_{1}\zeta T[0][g(\gamma+\delta)\partial_{x}\zeta]$ in the second equation of (3.59), modifies the dispersion relation (3.68) and we obtain instead:

[TABLE]

In this case, there exists a critical ratio for the depth of the two layers. Indeed, when $\displaystyle\delta^{2}<\gamma$, one has $\displaystyle\kappa_{2}<\kappa_{1}$, thus the perturbations are always stable if $\displaystyle\alpha\geq 1$. Whereas, when $\displaystyle\delta^{2}\geq\gamma$ i.e assuming $\displaystyle\gamma=0.5$, $\displaystyle\delta=0.8$ and $\displaystyle\dfrac{1}{\rm bo}=0$, one has $\displaystyle\kappa_{2}>\kappa_{1}$, thus the perturbations are unstable if $\displaystyle\underline{\zeta}>\dfrac{\alpha-1}{\kappa_{2}-\kappa_{1}}$. In fact, under the previous conditions the numerator of the right hand side of the dispersion relation (3.69) will become negative for sufficiently large values of $\displaystyle k^{2}$. Thus the root of $\displaystyle\tilde{w}$ will become complex inducing a high frequency instability of the model, see Figure 4. On the other hand, the numerator of the r.h.s of the dispersion relation (3.68) is positive for sufficiently large values of $\displaystyle k^{2}$, due to the existence of the term $\displaystyle-\nu\dfrac{(\kappa_{2}-\kappa_{1})k^{2}\underline{\zeta}}{(1+\nu\alpha k^{2})}$, thus $\displaystyle w$ is always real at high frequencies. This ensures the numerical stability of the model (3.59) which will be considered in the rest of this work.

4 Numerical methods

This section is devoted to the numerical methods developped to solve the improved Green-Naghdi equations (3.20). As pointed out by many authors [7, 31] the improved dispersion Green-Naghdi equations (3.20) is well-adapted to the implementation of a splitting scheme separating the hyperbolic and the dispersive parts of the equations. We present in Section 4.1 this splitting scheme inspired by [7, 31]. We explain in Sections 4.2 and 4.3 how we treat respectively the hyperbolic and dispersive parts of the equations.

We decided to treat the hyperbolic part by a finite volume method of Roe type. We will construct first order, second order “MUSCL” type method and finally 5th order WENO method. The high order method is suitable to compute correctly the maximum value of the height $\displaystyle\zeta$ and the discontinuities by limiting the diffusive effects. As we will show in the numerical validations, the high order scheme is also suitable to catch correctly the dispersive effects.

The dispersive part of the proposed splitting method is solved using a classical finite difference method.

4.1 The splitting scheme

Let us recall the improved Green-Naghdi system that we consider:

[TABLE]

$\displaystyle q_{3}$ defined by (3.3),$\displaystyle Q_{1}$, $\displaystyle Q_{2}$ and $\displaystyle Q_{3}$ are defined by (3.7), (3.21) and (3.22).

We decompose the solution operator $\displaystyle S(.)$ associated to (4.1) at each time step $\displaystyle\Delta t$ by the following second order splitting scheme:

[TABLE]

where $\displaystyle S_{1}(.)$ is the solution operator associated to the hyperbolic part, and $\displaystyle S_{2}(.)$ the solution operator associated to the dispersive part of the equations (4.1).

$\displaystyle\bullet\ S_{1}(t)$ is the solution operator associated to the hyperbolic part namely the nonlinear shallow water equations, NSWE:

[TABLE]

Using the definition of $\displaystyle q_{3}(\epsilon\zeta)$ given in (3.3), one can easily check that $\displaystyle\displaystyle{\dfrac{\varsigma}{2}+q_{3}(\epsilon\zeta)=\dfrac{f^{\prime}(\epsilon\zeta)}{2}}$. Thus we rewrite the NSWE system (4.2) in the following condensed form:

[TABLE]

We recall that, $\displaystyle f(\epsilon\zeta)=\displaystyle\frac{h_{1}h_{2}}{h_{1}+\gamma h_{2}}$ and $\displaystyle f^{\prime}(\epsilon\zeta)=\displaystyle\frac{h_{1}^{2}-\gamma h_{2}^{2}}{(h_{1}+\gamma h_{2})^{2}}$, with $\displaystyle h_{1}=1-\epsilon\zeta$ and $\displaystyle h_{2}=1/\delta+\epsilon\zeta$.

$\displaystyle\bullet\ S_{2}(t)$ is the solution operator associated to the remaining (dispersive) part of the equations.

[TABLE]

In this study, $\displaystyle S_{1}$ is computed using a finite volume method while $\displaystyle S_{2}$ is computed using a classical finite-difference method.

In order to discretize system (4.1), the numerical domain is an interval of length $\displaystyle L$ denoted $\displaystyle[0,L]$. Let $\displaystyle N\in\mathbb{N}^{*}$, and let us consider the following mesh on $\displaystyle[0,L]$. Cells are denoted for every $\displaystyle i\in[0,N+1]$, by $\displaystyle m_{i}=(x_{i-1/2},x_{i+1/2})$, with $\displaystyle x_{i}=\displaystyle\frac{x_{i-1/2}+x_{i+1/2}}{2}$ and $\displaystyle h_{i}=x_{i+1/2}-x_{i-1/2}$ the space mesh. The “fictitious” cells $\displaystyle m_{0}$ and $\displaystyle m_{N+1}$ denote the boundary cells and the mesh interfaces located at $\displaystyle x_{1/2}=0$ and $\displaystyle x_{N+1/2}=L$ are respectively the upstream and the downstream ends (see Figure 5).

We denote $\displaystyle\displaystyle\Delta x=\min_{i=1,N}h_{i}$.

We also consider a time discretization $\displaystyle t^{n}$ defined by $\displaystyle t^{n+1}=t^{n}+\Delta t$ with $\displaystyle\Delta t$ the time step.

4.2 Finite volume scheme

For the sake of simplicity in the notations, it is convenient to rewrite the hyperbolic system (4.3) in the following form:

[TABLE]

with the following conservative variables and flux function:

[TABLE]

The Jacobian matrix is given by:

[TABLE]

where $\displaystyle f^{\prime\prime}(\epsilon\zeta)=-\dfrac{2\gamma(h_{1}+h_{2})^{2}}{(h_{1}+\gamma h_{2})^{3}}$.

A simple computation shows that the homogeneous system (4.5) is strictly hyperbolic provided that:

[TABLE]

As a matter of fact, these conditions simply consist in assuming that the deformation of the interface is not too large and imposing a smallness assumption on $\displaystyle\epsilon v$. Notice that theses conditions correspond exactly to the hyperbolicity condition for the shallow water system provided in [26].

As a consequence, the solutions may develop shock discontinuities. In order to rule out the unphysical solutions, the system (4.5) must be supplemented by entropy inequalities, (see for instance [8] and references therein for more details). The Cauchy problem associated to (4.5) is the following:

[TABLE]

We are interested in the approximation of (4.9) by the finite volume method.

We denote $\displaystyle\overline{U}_{i}=(\zeta_{i},v_{i})$, the cell-centered approximation of $\displaystyle U$ on the cell $\displaystyle m_{i}$ at time $\displaystyle t$ given by:

[TABLE]

The piecewise constant representation of $\displaystyle U$ is given by, $\displaystyle U(t,x)=\displaystyle\overline{U}_{i}\mathds{1}_{m_{i}}(x)$.

We denote $\displaystyle\overline{U}_{i}^{n}=(\zeta_{i}^{n},v_{i}^{n})$, the cell-centered approximation of $\displaystyle U$ on the cell $\displaystyle m_{i}$ at time $\displaystyle t^{n}$ given by:

[TABLE]

The spatial discretization of the homogeneous system (4.2) can be recast under the following classical semi-discrete finite-volume formalism:

[TABLE]

where $\displaystyle\widetilde{F}$ is a numerical flux function based on a conservative flux consistent with the homogeneous NSWE:

[TABLE]

VFRoe method.

In what follows, we consider the numerical approximation of the hyperbolic system of conservation laws in the form of (4.5). To this end, we adopt the VFRoe method (see [10, 23, 24]) which is an approximate Godunov scheme. It relies on the exact resolution of the following linearized Riemann problem:

[TABLE]

where $\displaystyle\widetilde{A}(\overline{U}_{i}^{n},\overline{U}_{i+1}^{n})=A\left(\dfrac{\overline{U}_{i}^{n}+\overline{U}_{i+1}^{n}}{2}\right)$.

By solving the linearized Riemann problem we obtain $\displaystyle\overline{U}_{i+1/2}^{*}=U(x=x_{i+1/2},t=t_{n})$, the interface value between two neighbouring cells.

In what follows, we will detail the choice of the numerical flux for different order of approximation. The only change is in the computation of the interface values $\displaystyle\overline{U}_{i+1/2}^{*}$ which depends on the right and the left states of the linearized Riemann problem. For the sake of simplicity, we will still denote by $\displaystyle\widetilde{F}$ these numerical fluxes.

CFL condition

It is always necessary to impose what is called a CFL condition (for Courant, Friedrichs, Levy) on the timestep to prevent the blow up of the numerical values. It comes usually under the form

[TABLE]

where $\displaystyle a_{i+1/2}=\displaystyle\max_{i\in[1,N]}(j=1,2,|\lambda_{j}(\widetilde{U}_{i})|)$ and $\displaystyle\lambda_{j}(\widetilde{U}_{i})$ are the eigenvalues of \displaystyle A\big{(}\widetilde{U}_{i}=\dfrac{\overline{U}_{i}^{n}+\overline{U}_{i+1}^{n}}{2}\big{)}.

The restriction (4.13) enables in practice to compute the timestep at each time level $\displaystyle t_{n}$, in order to determine the new time level $\displaystyle t_{n+1}=t_{n}+\Delta t$ (within this view, $\displaystyle\Delta t$ is not constant, it is computed in an adaptive fashion).

Consistency.

The numerical flux $\displaystyle\widetilde{F}(U_{l},U_{r})$ is called consistent with (4.5) if

[TABLE]

4.2.1 First order finite-volume scheme

The semi-discrete equation (4.10) is discretised by an explicit Euler (in time) method to obtain :

[TABLE]

where the numerical flux is defined directly as the value of the exact flux at the interface value, namely:

[TABLE]

Let us remark that by construction the numerical flux given by (4.2.1) ensures the consistency property.

In the sequel, we will suppose that the space discretisation is uniform.

Algorithm.

In the following, we state the algorithm for computing the discrete values $\displaystyle\overline{U}_{i}^{n+1}$ at $\displaystyle t^{n+1}$. Given the initial data and boundary conditions and the number $\displaystyle CFL\leq 1$, we start with the known discrete averaged values $\displaystyle(\overline{U}_{i}^{n})$ for $\displaystyle i=0,...,N+1$ at $\displaystyle t^{n}$. As long as ($\displaystyle t<T$) one has to do:

1. Computation of $\displaystyle\widetilde{A}_{i}$ for $\displaystyle i=0,...,N$ where $\displaystyle\widetilde{A}_{i}=A\left(\dfrac{\overline{U}_{i}^{n}+\overline{U}_{i+1}^{n}}{2}\right)$.

2. Computation of $\displaystyle r_{i}^{1}$,$\displaystyle r_{i}^{2}$ and $\displaystyle\lambda_{i}^{1}$, $\displaystyle\lambda_{i}^{2}$ set respectively as the eigenvectors and eigenvalues of $\displaystyle\widetilde{A}_{i}$.

3. Computation of $\displaystyle\Delta t$, such that $\displaystyle\dfrac{\Delta t}{\Delta x}\leq\dfrac{CFL}{a_{i+1/2}}$.

4. Computation of $\displaystyle{\overline{U}_{i-1/2}^{*}}$ for $\displaystyle i=1,...,N+1$ by solving the linearized Riemann problem.

In fact we have 3 cases:

$\displaystyle\bullet$ if $\displaystyle\lambda_{i}^{1}$,$\displaystyle\lambda_{i}^{2}$¡0 then $\displaystyle{\overline{U}_{i-1/2}^{*}}=\overline{U}_{i}^{n}$.

$\displaystyle\bullet$ if $\displaystyle\lambda_{i}^{1}$,$\displaystyle\lambda_{i}^{2}$¿0 then $\displaystyle{\overline{U}_{i-1/2}^{*}}=\overline{U}_{i-1}^{n}$.

$\displaystyle\bullet$ if $\displaystyle\lambda_{i}^{1}<0$, $\displaystyle\lambda_{i}^{2}>0$ then for:

[TABLE]

with $\displaystyle R=(r_{i}^{1}|r_{i}^{2})$ and $\displaystyle[U]=\overline{U}^{n}_{i}-\overline{U}^{n}_{i-1}$.

5. Computation of $\displaystyle F({\overline{U}_{i-1/2}^{*}})$.

6. Computation of $\displaystyle\overline{U}^{n+1}_{i}=\overline{U}^{n}_{i}-\dfrac{\Delta t}{\Delta x}(F_{i+1/2}^{n}-F_{i-1/2}^{n})$ for $\displaystyle i=1,...,N$.

We repeat this algorithm for the new level of time $\displaystyle(t^{n+1}+\Delta t)$, until we reach the required final time $\displaystyle T$.

4.2.2 Second order finite-volume scheme: MUSCL-RK2

A drawback of the Roe scheme (such as Godunov) is to be very diffusive. A remedy for this situation is through the extension of the scheme to the second order in space, associated to a second order Runge-Kutta scheme in time. In fact, we would like to reduce both numerical dissipation and dispersion within the hyperbolic component $\displaystyle S_{1}(.)$. To this end, high order reconstructed states at each interface have to be considered, following the classical MUSCL approach [38] (Monotonic Upstream Scheme for Conservation Laws). To prevent the spurious oscillations that would occur around discontinuities or shocks, we suggest to use the “minmod” limiter, designed to generate slope limited, reconstructed left and right states for each cell that are used to calculate the flux at the interfaces. The implementation of this scheme is very easy and provides a natural extension of the Roe scheme described above. In fact, the main interests seem to be, after the tests, a gain of precision and stability.

The steps of the second-order reconstruction are as follows:

1. Using the discrete averaged values $\displaystyle\overline{U}_{i}^{n}$, we construct the slopes $\displaystyle S_{i}^{n}$ using the “minmod” limiter as the reconstruction must be non oscillatory in some sense, see [25]. We consider:

[TABLE]

where the function minmod is defined on $\displaystyle\mathbb{R}^{2}$ by

[TABLE]

2. On the cell $\displaystyle m_{i}=]x_{i-1/2},x_{i+1/2}[$, the solution is approached by:

[TABLE]

3. We compute the numerical flux at the interfaces:

[TABLE]

with:

[TABLE]

4. We compute $\displaystyle\widetilde{U}_{i}^{n+1}$ by the application of (4.15), thus the scheme is given as follows:

[TABLE]

As far as time discretization is concerned, we use the second-order explicit Runge–Kutta “RK2” method which is described in the following. Given the ODE $\displaystyle\dfrac{dy}{dt}=f(t,y)$, one has,

[TABLE]

with $\displaystyle\tilde{y}^{n+1}=y^{n}+hf(t^{n},y^{n})$, and $\displaystyle t^{n+1}=t^{n}+h$.

Applying (4.19) to (4.15), we obtain the following modified scheme “MUSCL-RK2”:

[TABLE]

with

[TABLE]

$\displaystyle-\widetilde{F}$: the numerical flux determined as in the first order VFRoe method, given by (4.2.1).

$\displaystyle-\widetilde{U}^{n+1}_{i}$ is computed as follows:

[TABLE]

$\displaystyle-\widetilde{U}^{n+1,+}_{i}=\widetilde{U}^{n+1}_{i}+\dfrac{\Delta x}{2}\widetilde{S}_{i}^{n+1}.$

$\displaystyle-\widetilde{U}^{n+1,-}_{i+1}=\widetilde{U}^{n+1}_{i+1}-\dfrac{\Delta x}{2}\widetilde{S}_{i+1}^{n+1}.$

$\displaystyle-\widetilde{S}_{i}^{n+1}$ and $\displaystyle\widetilde{S}_{i+1}^{n+1}$ are associated respectively to $\displaystyle\widetilde{U}_{i}^{n+1}$ and $\displaystyle\widetilde{U}_{i+1}^{n+1}$ by (4.18).

4.2.3 Higher order finite-volume scheme: WENO5-RK4

The second order schemes are known to degenerate to first order accuracy at smooth extrema. To reach higher order accuracy in smooth regions and a good resolution around discontinuities, we implement fifth-order accuracy WENO reconstruction, following [28, 36], where Jiang and Shu constructed third and fifth order finite difference WENO schemes in multi-space dimensions with a general framework for the design of the smoothness indicators and nonlinear weights. To automatically achieve high order accuracy and non-oscillatory property near discontinuities, WENO schemes use the idea of adaptive stencils in the reconstruction procedure based on the local smoothness of the numerical solution. We would like to mention also the previous studies [11, 7, 31], where it is shown that for the study of dispersive waves, it is necessary to use high-order schemes to prevent the corruption of the dispersive properties of the model by some dispersive truncation errors linked to second-order schemes. Using the same reconstruction proposed in [7], we consider a cell $\displaystyle m_{i}$, and the corresponding constant averaged value $\displaystyle\overline{U}_{i}^{n}=(\overline{\zeta}_{i}^{n},\overline{v}_{i}^{n})$, with a constant space step $\displaystyle\Delta x$. This approach, provides high order reconstructed left and right values $\displaystyle\overline{U}_{i}^{n,-}$ and $\displaystyle\overline{U}_{i}^{n,+}$, built following the five points stencil, and introduced as follows:

[TABLE]

where $\displaystyle\overline{\delta U}_{i}^{n,+}$ and $\displaystyle\overline{\delta U}_{i}^{n,-}$ are defined as follows:

[TABLE]

[TABLE]

and the coefficients $\displaystyle\dfrac{2}{3}$, $\displaystyle\dfrac{1}{3}$, $\displaystyle\dfrac{-1}{10}$ and $\displaystyle\dfrac{-1}{15}$ are set in order to obtain better dissipation and dispersion properties in the truncature error. We consider the following modified scheme:

[TABLE]

To reduce spurious oscillations near discontinuities, we apply the same limitation procedure as in [7], preserving the scheme positivity and the high order accuracy. Thus scheme (4.24) becomes

[TABLE]

We define the limited high-order values as follows:

[TABLE]

Using the following limiter,

[TABLE]

we define the limiting process as,

[TABLE]

with $\displaystyle\delta\overline{U}_{i+1}^{n}=\overline{U}_{i+1}^{n}-\overline{U}_{i}^{n}$ and $\displaystyle\delta\overline{U}_{i}^{n}=\overline{U}_{i}^{n}-\overline{U}_{i-1}^{n}$ are upstream and downstream variations, and $\displaystyle\overline{\delta U}_{i}^{n,+}$ and $\displaystyle\overline{\delta U}_{i}^{n,-}$ taken from (4.22) and (4.23).

The limited high order reconstructions stated above must be performed on both conservative variables $\displaystyle\overline{U}_{i}^{n}=(\overline{\zeta}_{i}^{n},\overline{v}_{i}^{n})$. We would like to mention that the resulting finite volume scheme preserve motionless steady states, $\displaystyle\zeta=\text{cst}$ and $\displaystyle v=0$.

As far as time discretization is concerned, we use the fourth-order explicit Runge–Kutta “RK4” method which is described in the following. Given the ODE $\displaystyle\dfrac{dy}{dt}=f(t,y)$, one has,

[TABLE]

with $\displaystyle t^{n+1}=t^{n}+h$. Applying (4.2.3) to (4.25), one gets the “WENO5-RK4” scheme.

4.3 Finite difference scheme for the dispersive part

First of all, let us detail how to construct nodal values of the unknowns (which are the ones used for a finite difference discretization) from the cell-averaged value computed by a finite volume scheme and vice versa.

We denote by $\displaystyle U_{i}^{n}$ the nodal value of $\displaystyle U$ at the $\displaystyle i$th node $\displaystyle(x_{i+1/2})_{i\in[0,N]}$ and at time $\displaystyle t_{n}$ i.e. $\displaystyle U_{i}^{n}$ is an approximation of $\displaystyle U(x_{i+1/2},t^{n})$. The finite volume-finite difference mix imply to switch between the cell-averaged and nodal values for each unknown and at each time step. To this end, we use the fifth-order accuracy WENO reconstruction, that allows to approximate the nodal values (i.e finite difference unknowns) $\displaystyle(U_{i}^{n})_{i=1,N+1}$ in terms of the cell-averaged values (i.e finite volume unknowns) $\displaystyle(\overline{U}_{i}^{n})_{i=1,N}$. The relation is given by:

[TABLE]

with adaptations at the boundaries following the method presented in Section 4.4. One can easily recover the relation that allows to determine the cell-averaged values $\displaystyle(\overline{U}_{i}^{n})_{i\in[1:N]}$ in terms of the nodal values $\displaystyle(U_{i}^{n})_{i\in[1:N+1]}$.

We can easily check that (4.28) preserve the steady state at rest and that this formula is precise up to order $\displaystyle\mathcal{O}(\Delta x^{5})$ terms, thus preserving the global order of the scheme.

We can now proceed by explaining how we compute the solution operator $\displaystyle S_{2}(.)$ associated to the dispersive part of the equations. Let us recall the system (4.4) corresponding to the operator $\displaystyle S_{2}(.)$, given in Section (4.1).

[TABLE]

where the operators $\displaystyle Q_{1}$, $\displaystyle Q_{2}$ and $\displaystyle Q_{3}$ are explicitly given in (3.7), (3.21) and (3.22). For the sake of simplicity, we detail the numerical resolution of (4.29) using an explicit Euler in time scheme. Standard extensions to second and fourth order Runge-Kutta method has been used according to the order of the space derivative operators as done in the previous section.

The finite discretization of the system (4.29) leads to the following discrete problem:

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

The system (4.30) is solved at each time step using a classical finite-difference technique, where the matrices $\displaystyle D_{1}$ and $\displaystyle D_{2}$ are the classical centered discretizations of the derivatives $\displaystyle\partial_{x}$ and $\displaystyle\partial^{2}_{x}$ given below.

The first formula is the second-order formula called “DF2”, where the spatial derivatives are given as follows:

[TABLE]

[TABLE]

For time discretization, the second-order formula “DF2” is associated to a second-order classical Runge-Kutta “RK2” scheme, and thus one obtains the “DF2-RK2” scheme.

The second formula is the fourth-order formula called “DF4” where the spatial derivatives are given as follows:

[TABLE]

[TABLE]

For time discretization, the fourth-order formula “DF4” is associated to a fourth-order classical Runge-Kutta “RK4” scheme, and thus one obtains the “DF4-RK4” scheme.

4.4 Boundary conditions

In the following section, we show how to treat the boundary conditions for the hyperbolic and dispersive parts of the splitting scheme. Suitable relations are imposed on both cell-averaged and nodal quantities. We only treat either periodic boundary conditions or reflective boundary conditions. We detail now how we have implemented these boundary conditions for the hyperbolic part and the dispersive part of the numerical scheme.

For the hyperbolic part, we have introduced “ghosts cells” respectively at the upstream and downstream boundaries of the domain. The imposed relations on the cell-averaged quantities are the following:

$\displaystyle\bullet$ $\displaystyle\overline{U}_{-k+1}=\overline{U}_{N-k+1}$, and $\displaystyle\overline{U}_{N+k}=\overline{U}_{k}$, $\displaystyle k\geq 1$, for periodic conditions on upstream and downstream boundaries.

$\displaystyle\bullet$ $\displaystyle\overline{\zeta}_{-k+1}=\overline{\zeta}_{-k}$, $\displaystyle\overline{v}_{-k+1}=-\overline{v}_{-k}$ and $\displaystyle\overline{\zeta}_{N+k-1}=\overline{\zeta}_{N+k}$, $\displaystyle\overline{v}_{N+k-1}=-\overline{v}_{N+k}$, $\displaystyle k\geq 1$, for reflective conditions on the left and right boundaries.

For the dispersive part of the splitting, the boundary conditions are simply imposed on the nodal values that are located outside of the domain, in order to maintain centered formula at the boundaries, while keeping a regular structure in the discretized model:

$\displaystyle\bullet$ $\displaystyle U_{-k+1}=U_{N-k+1}$, and $\displaystyle U_{N+k}=U_{k}$, $\displaystyle k\geq 1$, for periodic conditions on upstream and downstream boundaries.

$\displaystyle\bullet$ $\displaystyle\zeta_{-k+1}=\zeta_{-k}$, $\displaystyle v_{-k+1}=-v_{-k}$ and $\displaystyle\zeta_{N+k-1}=\zeta_{N+k}$, $\displaystyle v_{N+k-1}=-v_{N+k}$, $\displaystyle k\geq 1$, for reflective conditions on upstream and downstream boundaries.

5 Numerical validations

In this section, several numerical tests are performed in both one and two layers cases in order to validate the numerical efficiency and accuracy of the improved Green-Naghdi model (3.59). We first consider several numerical tests in the one layer case without any surface tension i.e $\displaystyle\rm bo^{-1}=0$. We begin by studying the propagation of a solitary wave over a flat bottom. We compare our numerical solution with an analytic one (up to an $\displaystyle\mathcal{O}(\mu^{2})$ remainder) at several times and for different orders of discretizations and show that our numerical scheme is very efficient and accurate. To evaluate the influence of the nonlinear and dispersive terms we study the collision between two solitary waves traveling in opposite directions (head-on collision). We then study the breaking of a Gaussian hump into two solitary waves. Finally, we study the dam-break problem supplemented by a comparison between the second and fifth order accuracy in order to show the ability of the higher order numerical scheme in dealing with discontinuities. We used the value $\displaystyle\alpha=1$ in the aforementioned cases. In fact, we obtained very similar results when performing the same simulations with $\displaystyle\alpha=1.159$. Secondly, two numerical simulations are performed in the two layers cases. In the first one, we compare our results with numerical data from [19], where we show that a very good matching is observable if $\displaystyle\alpha$ is carefully chosen as in Section 3.3. We then consider the dam-break problem in the two layers cases, where we test the ability of the splitting scheme to compute dispersive shock waves with high accuracy. We would like to mention that in all the numerical tests, we use the WENO5 reconstruction for the hyperbolic part of the splitting scheme and a fourth order finite difference scheme “DF4” for the dispersive part, both associated to a fourth-order classical Runge-Kutta “RK4” time scheme. In every numerical simulation presented in the following section, we choose to use a CFL number equal to 1 in the algorithm stated in page 4.2.1, in order to obtain a stable numerical scheme.

5.1 Numerical validations in the one layer case

5.1.1 Propagation of a solitary wave

Here, we test the accuracy of our numerical scheme (3.59) with $\displaystyle\alpha=1$, by using the exact solitary wave solutions of the one layer Green-Naghdi equations in the one-dimensional setting over a flat bottom (see [31]), given in variables with dimensions, by

[TABLE]

where we recall that $\displaystyle d_{2}$ is depth of the fluid when considering the one layer case. Such solitary waves are also solution of the improved Green-Naghdi model (3.59) up to an $\displaystyle\mathcal{O}(\mu^{2})$ remainder. We consider the propagation of a solitary wave initially centered at $\displaystyle x_{0}=20\ m$, of relative amplitude $\displaystyle a=0.2\ m$, over a constant water depth $\displaystyle d_{2}=1\ m$. The computational domain is $\displaystyle 200\ m$ in length and discretized with 1280 cells. The single solitary wave propagates from left to right. In this test, since the solitary wave is initially far from boundaries, the boundary conditions do not affect the computation, thus we choose to impose periodic boundary conditions on each boundary for the sake of simplicity. We compare the water surface profile of our numerical solution provided by the model (3.59) (after setting the parameters $\displaystyle\gamma=0$ and $\displaystyle\delta=1$ corresponding to the one layer case), with the exact one given by (5.1) at several times using the first, second and fifth order discretizations (see Figure 6). One can easily remark that the fifth order discretization “WENO5-DF4-RK4” provides very accurate solutions and reduces both numerical dissipation and dispersion, contrarily to the first order discretization “FV1-DF2-Euler” which seems to be very diffusive. Indeed, looking at the amplitude and shape of the solitary wave at $\displaystyle t=50\ s$ in the bottom of Figure 6, we can observe an excellent agreement between numerical and exact solutions, unlike the second order discretization “MUSCL-DF2-RK2” (middle of Figure 6), generating a less accurate numerical solution. The preservation of the amplitude and shape of the solitary wave computed using the fifth order scheme, even after $\displaystyle 200\ m$ indicate that the governing equations have been accurately discretized in space and time.

In what follows, we quantify the numerical accuracy of our numerical scheme by computing the numerical solution for this particular test case for an increasing number of cells $\displaystyle N$, over a duration $\displaystyle T=5\ s$. Starting with $\displaystyle N=80$ number of cells, we successively multiply the number of cells by two. The relative errors $\displaystyle E_{L^{2}}(\zeta)$ and $\displaystyle E_{L^{2}}(v)$ on the water surface deformation and the averaged velocity presented in Table 1 are computed at $\displaystyle t=5\ s$, using the discrete $\displaystyle L^{2}$ norm $\displaystyle\|.\|_{2}$:

[TABLE]

where $\displaystyle(\zeta_{num},v_{num})$ are the numerical solutions and $\displaystyle(\zeta_{sol},v_{sol})$ are the analytical ones coming from (5.1). Very accurate results are thus obtained as an evaluation of the capacity of our numerical method to compute in a stable way the propagation of a solitary wave.

Remark 5.3

We believe that the main reason for not obtaining the predicted order in each space discretization might be due to the fact that the analytic solution given in (5.1) satisfies the model (3.59) up to an $\displaystyle\mathcal{O}(\mu^{2})$ remainder, that is to say it is an approximate solution. A remedy for this situation could be through an “iterative cleaning” technique that acts to damp the high frequency oscillations (i.e the oscillatory dispersive tails) that appears due to the remainder term of size $\displaystyle\mathcal{O}(\mu^{2})$. This technique has been used by several authors, see for instance [3, 4, 35]. In this paper, we do not try to give some optimal convergence result and the filtering technique is left to future work.

5.1.2 Head-on collision of counter-propagating waves

We now investigate the interaction of solitary waves which allows us to evaluate the impact of non linearities and dispersive terms. To this end, we study the head-on collision of two solitary waves traveling in opposite directions. The initial data for the two counter-propagating solitary waves are given in (5.1). Many authors have set different models and numerical methods in order to numerically study this problem (see [14, 20, 34]). Unlike solitary waves in integrable systems, the one for the full Euler equations are often followed by a non zero residual wave after interactions. The resulting residual has the form of dispersive trailing waves of small amplitude. In this test, we study the interaction of two solitary waves of equal amplitudes propagating in an opposite direction, initially located at $\displaystyle x=100\ m$ and $\displaystyle x=300\ m$. The spatial domain is $\displaystyle 400\ m$ in length discretized using $\displaystyle 1200$ cells. Periodic conditions are imposed on each boundary.

The results are shown in Figure 7 at different propagation times. Before the collision, at times $\displaystyle t=20\ s$ and $\displaystyle t=25\ s$ one can observe two dispersive tails of very small amplitude located to the left and right of each solitary wave. The generation of such dispersive tails is due to the $\displaystyle\mathcal{O}(\mu^{2})$ remainder term as mentioned in Remark 5.3. As expected, the waves collide to reach a maximum height larger than the sum of the amplitudes of the two incident solitary waves at the time $\displaystyle t=27\ s$. After the collision, dispersive tails with small amplitudes appear clearly when zooming at $\displaystyle t=70\ s$, illustrating an appropriate characterization of the nonlinear interactions. Capturing this dynamics validate the high precision of our numerical scheme. The head-on collision is also studied in [34, 20], leading to similar observations.

5.1.3 Breaking of a Gaussian hump

In this test, we consider the following initial data representing a heap of water,

[TABLE]

where $\displaystyle a$ represents the amplitude, $\displaystyle\lambda$ represents the width and $\displaystyle L$ represents the length of the domain. Figure 8, shows the overall behavior of the solutions using $\displaystyle a=0.4\ m$ and $\displaystyle\lambda=40\ m$ and $\displaystyle L=400\ m$ discretized with $\displaystyle 2000$ cells using periodic boundary conditions. The initial hump breaks up into two large solitary waves and smaller dispersive tails. These waves and the dispersive tails travel in opposite directions.

The results shown in Figure 8, tends to confirm the ability of our numerical scheme to capture this dynamics accurately. Indeed, this typical behaviour is also studied in [34], leading to very similar observations.

5.1.4 Dam-break problem in the one layer case

We consider now a dam-break problem in the one layer case in order to test the ability of our numerical scheme to deal with discontinuities. In general, discontinuous initial data of this type generates dispersive shock waves due to the dispersive effects [32]. Analytic and computational studies of the dispersive shock waves in fully-nonlinear dispersive shallow water systems were carried out in [21, 32]. We would like to mention also the previous studies [11, 7, 31], where it is shown that for the study of dispersive waves, it is necessary to use high-order schemes to prevent the corruption of the dispersive properties of the model by some dispersive truncation errors linked to second-order schemes. Indeed, this test is supplemented by a comparison between the second and fifth order accuracy exhibiting the ability of higher order schemes to capture the rapid oscillations in dispersive shock waves. We use the following initial data:

[TABLE]

where $\displaystyle a=0.2091\ m$. The computational domain is the interval $\displaystyle x\in(-700,700)$ and discretized using 2800 cells. We choose to impose periodic boundary conditions.

Figure 9 shows the results of the numerical simulations with two different orders of discretization, “MUSCL-DF2-RK2” and “WENO5-DF4-RK4”, generating two dispersive tails propagating in opposite direction, on the left and right side of the “dam”, and two rarefaction waves that travel towards the center. A zoom on the shock and rarefaction waves at $\displaystyle t=30\ s$ and $\displaystyle t=65\ s$, shows clearly the corruption of the dispersive properties by the second order scheme. The dam-break problem was also studied in [21, 32] and seems to fit well with our numerical observations. The numerical model proposed in [34] computes using a finite element method all the nonlinear dispersive terms, in particular the third order ones. These third order derivatives are present in our model but in order to improve the frequency dispersion we have proposed to factorize these terms, making it possible not to compute them. This formulation has the inconvenience of numerical diffusion. This is the reason why the dispersive tails observed in the dam-break problem in [34] have larger amplitude oscillations.

5.2 Numerical validations in the two layers case

5.2.1 Kelvin-Helmholtz instabilities

In this case, we would like to highlight the importance of the choice of the parameter $\displaystyle\alpha$ in order to improve the frequency dispersion of the model (3.20), through the simulation of a sufficiently regular initial wave, following the numerical experiments performed in [19]. In the aforementioned paper they introduce a new class of Green-Naghdi type models for the propagation of internal waves with improved frequency dispersion in order to prevent high-frequency Kelvin-Helmholtz instabilities. These models are obtained by regularizing the original Green-Naghdi one by slightly modifying the dispersion components using a class of Fourier multipliers. They represent three different choices of the Fourier multipliers, each one yields to a specific Green-Naghi model which they denote as follows:

$\displaystyle\bullet$ “original” as the classical Green-Naghi model introduced in [12].

$\displaystyle\bullet$ “regularized” which is a well-posed system for sufficiently small and regular data, even in absence of surface tension. In addition its dispersion relation fit the one of full Euler system at order $\displaystyle\mathcal{O}(\mu^{3})$.

$\displaystyle\bullet$ “improved” whose dispersion relation is the same as the one of the full Euler system.

Using the spectral methods [37] for the space discretization and the Matlab solver ode45, which is based on the fourth and fifth order Runge-Kutta-Merson method for time evolution, they numerically compute several of their Green-Naghdi systems. Several computations are made, with and without surface tension in order to observe how the different frequency dispersion may affect the appearance of Kelvin-Helmholtz instabilities. In order to compare with the numerical experiments done in [19] we choose the initial data $\displaystyle\zeta(0,x)=-e^{-4|x|^{2}}$ and $\displaystyle v(0,x)=0$ (represented by the dashed lines). The dimensionless parameters are set as follows: $\displaystyle\mu=0.1$, $\displaystyle\epsilon=0.5$, $\displaystyle\delta=0.5$, $\displaystyle\gamma=0.95$ and $\displaystyle\rm bo^{-1}=5\times 10^{-5}$. The computational domain is the interval $\displaystyle x\in(-4,4)$ discretized with 512 cells using periodic boundary conditions. In all the following simulations we compute our numerical solution using the fifth order accuracy scheme “WENO5-DF4-RK4”.

Figures 10 and 11 show the comparisons between our numerical solution for $\displaystyle\alpha_{opt}=1$ (left) and $\displaystyle\alpha_{opt}=1.271$ (right) and the Green-Naghdi models solutions obtained in [19], with a small amount of surface tension, at time $\displaystyle t=2$ and $\displaystyle t=3$ respectively. We observe an excellent agreement between our numerical solution computed for $\displaystyle\alpha_{opt}=1.271$ and both “improved” and “regularized” models at $\displaystyle t=2$ and $\displaystyle t=3$. As expected, at $\displaystyle t=3$ the original model induces Kelvin-Helmholtz instabilities. Meanwhile, the flows predicted by the regularized and improved models and by our model (3.20) with $\displaystyle\alpha_{opt}=1.271$ remain smooth and are very similar. Similarly, Figure 12 shows an excellent agreement between the numerical solutions computed for $\displaystyle\alpha_{opt}=1.271$ with the “improved” and “regularized” models, without surface tension at time $\displaystyle t=2$, while the flow of the original model is completely destroyed due to Kelvin-Helmholtz instabilities. However, at a larger time ($\displaystyle t=5$), figure 13 shows that in absence of surface tension ($\displaystyle\rm bo^{-1}=0$) Kelvin-Helmholtz instabilities will destroy completely both “original” and “improved” models while the numerical solution computed for $\displaystyle\alpha_{opt}=1.271$ and for the “regularized” models remain smooth and very similar.

The overall observations show the importance of the choice of the parameter $\displaystyle\alpha$ as in Section 3.3 in improving the frequency dispersion. Indeed, when choosing $\displaystyle\alpha_{opt}=1.271$, we observe an excellent matching between our numerical solutions and those obtained by the “improved” model before the latter is completely destroyed in absence of surface tension due to the Kelvin-Helmholtz instabilities. As well, our numerical solution matches the one computed by the “regularized” model even for a large time and with or without surface tension. This is not the case when choosing $\displaystyle\alpha_{opt}=1$. In fact, the “improved” model has the same dispersion relation as the one of the full Euler system and the dispersion relation of the “regularized” model fit the one of the full Euler system to an $\displaystyle\mathcal{O}(\mu^{3})$ order. This explains the reason behind the matching when choosing an optimal value for $\displaystyle\alpha$.

5.2.2 Dam-break problem in the two layers case

This simulation concerns a test of the ability of our numerical scheme (3.59) to handle discontinuities when considering a dam-break problem in the two layers case. To this end, we use the same initial data as in the one layer case given by (5.4), where $\displaystyle a=0.2091\ m$. The simulations are performed on the interval $\displaystyle x\in(-700,700)$, discretized with 2800 cells imposing periodic conditions on each boundary. The dimensionless parameters representing the ratio between the depth and the ratio between the densities of the two layers are set respectively as follows: $\displaystyle\delta=0.5$, $\displaystyle\gamma=0.95$. Taking into account a small amount of surface tension, we set $\displaystyle\rm bo^{-1}=5\times 10^{-5}$. We would like to mention that the simulations are computed using $\displaystyle\alpha=1$. As expected, the same simulations performed when choosing $\displaystyle\alpha=1.1498$ yielded the same result since $\displaystyle\alpha_{opt}\rightarrow 1$ when considering large wave numbers as explained in Section 3.3. Figure 14 shows the results of the numerical simulation at several times, generating two dispersive tails propagating towards the center and two rarefaction waves that travel in opposite directions, on each side of the “dam”. A zoom on the dispersive tails is proposed at $\displaystyle t=55\ s$ and $\displaystyle t=75\ s$. Indeed, capturing this phenomenon accurately exhibit the high accuracy of our numerical scheme.

6 Conclusion

In this work, we presented a numerical scheme for the Green-Naghdi model in the Camassa-Holm regime. This model is first reformulated under more appropriate structure, where the time dependency of a second order differential operator present in the model is removed, keeping the stabilizing effects of its inversion. Furthermore, we improved the frequency dispersion of the original model keeping the same order of precision thanks to a parameter $\displaystyle\alpha$ to be precisely chosen. Additionally, our model do not contain third order derivatives that may create high frequencies instabilities. We then propose an efficient, precise and stable numerical splitting scheme that decomposes the hyperbolic and dispersive parts of the equations. The hyperbolic part is treated with a finite-volume method, where we consider the following schemes with different orders of accuracy: the first order finite-volume scheme based on a VFRoe method, the second order finite volume scheme following the classical MUSCL approach and finally a fifth-order WENO reconstruction. On the other hand, we treat the dispersive part with a finite-difference scheme, using second and fourth order formulas. Concerning time discretization we use classical second and fourth order Runge-Kutta methods, according to the order of the space derivative in consideration. Finally, we present several numerical validations in the one and two layers cases, showing the numerical efficiency and accuracy of our scheme and exhibiting its ability to reduce numerical dispersion and dissipation and to deal with discontinuities. The next step of this study may concern the extension of this numerical scheme to a more general configuration of variable topography. We believe that this splitting strategy may be applied in the variable bottom case.