A combined finite volume - finite element scheme for a dispersive shallow water system

. We propose a variational framework for the resolution of a non-hydrostatic Saint-Venant type model with bottom topography. This model is a shallow water type approximation of the incompressible Euler system with free surface and slightly diﬀers from the Green-Nagdhi model, see [12] for more details about the model derivation. The numerical approximation relies on a prediction-correction type scheme initially introduced by Chorin-Temam [16] to treat the incompressibility in the Navier-Stokes equations. The hyperbolic part of the system is approximated using a kinetic ﬁnite volume solver and the correction step implies to solve a mixed problem where the velocity and the pressure are deﬁned in compatible ﬁnite element spaces. The resolution of the incompressibility constraint leads to an elliptic problem involving the non-hydrostatic part of the pressure. This step uses a variational formulation of a shallow water version of the incompressibility condition. Several numerical experiments are performed to conﬁrm the relevance of our approach.

Such an assumption produces important consequences over the structure and complexity of the model.Indeed, Eq. ( 1) implies that the pressure p is no longer the Lagrange multiplier of the incompressibility constraint and p can be expressed, for free surface flows, as a function of the water depth of the fluid.Therefore, the hydrostatic assumption implies that the resulting model, even though it describes an incompressible fluid, has common features with models arising in compressible fluid mechanics.
In geophysical problems, the hydrostatic assumption coupled with a shallow water type description of the flow is often used.Unfortunately, these models do not represent phenomena containing dispersive effects for which the non-hydrostatic contribution cannot be neglected.And more complex models have to be considered to take into account this kind of phenomena, together with numerical methods able to discretize the high order derivative terms coming from the dispersive effects.Many shallow water type dispersive models have been proposed such as KdV, Boussinesq, Green-Naghdi, see [21,14,6,30,31,27,19,26,2,3,13].The modeling of the non-hydrostatic effects for shallow water flows does not raise insuperable difficulties but their discretization is more tricky.Numerical techniques for the approximaion of these models have been recently proposed [15,11,28].
The model studied in the present paper has been derived and studied in [12].Its numerical approximation based on a projection-correction strategy [16] is described in [1].In [1], the discretization of the elliptic part arising from the non-hydrostatic terms is carried out in a finite difference framework.It is worth noticing that the numerical scheme given in [1] is endowed with robustness and stability properties such as positivity, well-balancing, discrete entropy and wet/dry interfaces treatment.
The main contents of this paper is the derivation and validation of the correction step in a variational framework.Since the derivation in a 2d context of the model proposed in [12] does not raise difficulty, the results depicted in this paper pave the way for a discretization of the 2d model over an unstructured mesh, and we will often maintain general notations as far as possible.
Notice that the non-hydrostatic model we consider slightly differs from the wellknown Green-Naghdi model [21] but the numerical approximation proposed in this paper can also be used for the numerical approximation of the Green-Naghdi system.
Let Ω ⊂ R, be a 1d domain (an interval) and Γ = Γ in ∪ Γ out its boundary (see figure 1).The non-hydrostatic model derived in [12,1] reads where H is the water depth, z b the topography and p nh the non-hydrostatic part of the pressure.The variables denoted with a bar recall that this model is obtained performing an average along the water depth of the incompressible Euler system with free surface.The velocity field is denoted ū = (ū, w) t with ū (resp.w) the horizontal (resp.vertical) component.
We denote η = H + z b the free surface of the fluid.In addition, we give the following notations with n the unit outward normal vector at Γ (in 1d, n = ± 1), n represents the unit outward normal vector of the domain covered by the fluid, namely Ω × [z b , η].We also consider the gradient operator Figure 1.Notations and domain definition.
The smooth solutions of the system (2)-( 5) satisfy moreover an energy conservation law Note that (5) represents a shallow water version of the divergence free constraint, for which the non hydrostatic pressure pnh plays the role of a Lagrange multiplier.Notice that considering pnh = 0 and neglecting (4), the system (2)-( 3), (5) reduces to the classical Saint-Venant system.
The paper is organized as follows.First we give a rewriting of the model and we present the prediction-correction method, the main part being the variationnal formulation of the correction part.Then in Section 3, we detail the numerical approximation.Finally, in Section 4, numerical simulations validating the proposed discretization techniques are presented.

2.
The projection scheme for the non-hydrostatic model.Projection methods have been introduced by A. Chorin and R. Temam [33] in order to compute the pressure for incompressible Navier-Stokes equations.These methods, based on a time splitting scheme, have been widely studied and applied to treat the incompressibility constraint (see [24,35,34]).We develop below an analogue of this method for shallow water flow.In order to describe the fractional time step method we use, we propose a rewritting of the model ( 2)-(5).

A rewritting.
Let us introduce the two operators ∇ sw and div sw defined by with v = (v 1 , v 2 ) t .We assume for a while that f and v are smooth enough.The shallow water form of the divergence operator div sw (resp. of the gradient operator ∇ sw ) corresponds to a depth averaged version of the divergence (resp.gradient) appearing in the incompressible Euler and Navier-Stokes equations.Notice that the two operators ∇ sw , div sw defined by ( 10)-( 11) are H and z b dependent and we assume that H and z b are sufficiently smooth functions.One can check that these operators verify the fundamental duality relation These definitions allow to rewrite the model ( 2)-( 5) as with ∇ 0 defined by (7).
The system ( 13)-( 15) can be written in the compact form where we denote and Let be given time steps ∆t n and note t n = k≤n ∆t k .As detailed in [1], the projection scheme for system (16)- (17) consists in the following time splitting The first two equations of (20) consist in the classical Saint-Venant system with topography and the third equation is an advection equation for the quantity H w. Equations ( 21)- (22) describe the correction step allowing to determine the non hydrostatic part of the pressure p n+1 nh and hence giving the corrected state X n+1 .The numerical resolution of (20) -especially the first two equations -has received an extensive coverage and efficient and robust numerical techniques exist, mainly based on finite volume approach, see [8,5].The derivation of a robust and efficient numerical technique for the resolution of the correction step (21)-( 22) is the key point.A strategy based on a finite difference approach has been proposed, studied and validated in [1].Unfortunately, the finite difference framework does not allow to tackle situations with unstructured meshes in 2 or 3 dimensions.It is the key point of this paper to propose a variational formulation of the correction step coupled with a finite volume discretization of the prediction step.
The time discretization in the numerical scheme described above corresponds to a fractional time step strategy with a first order Euler scheme, explicit for the hyperbolic part and implicit for the elliptic part.For hyperbolic conservation laws, the second-order accuracy in time is usually recovered by the Heun method [7,8] that is a slight modification of the second order Runge-Kutta method.More precisely, for a dynamical system written under the form the Heun scheme consists in defining y n+1 by with The model we have to discretize has the general form B(X) = 0, and we propose the numerical scheme with the two steps defined by and Xn+2 = Xn+1 + ∆t n F ( Xn+1 , pn+2 ), ( 29) where pn+1 and pn+2 are the solutions of the elliptic equations derived from the divergence free constraints ( 28) and ( 30) respectively.Notice that ( 27) is a compact form of the fractional scheme ( 20)- (21) where the intermediate step no more appears.It is easy to prove the scheme ( 26) is second order accurate.Indeed assuming F and B are enough smooth, we have or equivalently using ( 25) Now from ( 27) and ( 29), we get Using ( 26), we see that relations (31) and ( 32) are equivalent up to third order terms.
2.2.The correction step.In this part, we consider we have at our disposal a space discretization of Eq. ( 20) solving the hydrostatic part of the model and we focus on the correction step ( 21)-( 22).
2.2.1.Variational formulation.The correction step ( 21)-( 22) writes, For the sake of clarity, in the following we will drop the notation with a bar and we denote p instead of pnh .Likewise we drop the superscript n+1 for the corrected states.Equations ( 34)-( 35) is a mixed problem in velocity/pressure, its approximation leads to the variational mixed problem and such that Introducing the bilinear forms b(u, q) = 0, ∀q ∈ Q.
(41) 2.2.2.The pressure equation.Formally, we take v in the form v = ∇sw q H with q ∈ Q 0 and Then, with ( 12) and (41), we get If we introduce the shallow water version of the Laplacian operator ∆ sw it is natural to consider the new variational formulation: From (43), we deduce The resolution of the equation (44) allows to update the velocity at the correction step (34).
Notice that the equation ( 44) is equivalent to apply the operator div sw to the equation ( 21) and to use the shallow water free divergence (22) to eliminate u.
Remark 2. Notice that Eq. ( 44) has the form of a Sturm-Liouville type equation.

Inf-sup condition.
To ensure that the saddle problem (40)-( 41) is well posed, the Babuska-Brezzi [9,32] condition has to be satisfied.Denoting by B the weak operator defined by ∀v ∈ V 0 , Bv = b(v, q), ∀q ∈ Q, we have We assume H ≥ H 0 > 0. Because of the positivity of H, it is obvious that the bilinear form a is coercive.Choosing v ∈ (L 2 (Ω)) 2 and q such as ∇ sw q ∈ (L 2 (Ω)) 2 , it follows that ∇ sw q = 0, then q = 0 and ker B t = 0. Indeed, in contrast with Navier-Stokes equations for which the pressure is defined up to an additive constant, the non hydrostatic of the shallow water equations is fully defined.Therefore, the mixed problem (40)-( 41) satisfies the inf-sup condition and admits a unique solution.
2.3.Boundary conditions.In this section, we still consider that the hydrostatic part is provided and we study the compatibility of the boundary conditions between the hydrostatic part and the projection part.Therefore, the compatibility between the pressure and velocity at boundary needs to be studied.To this aim, we first provide the conditions required to impose Dirichlet or Neumann pressure at boundary, and then, we couple these conditions with the hydrostatic part.
We consider a more general case taking the space and we introduce the bilinear form with n the unit outward normal vector defined by (6).In one dimension, c(v, p) = (Hpv Therefore instead of ( 40)-(41), we consider the problem: with ∇ 0 defined by (7) and n s (resp.n b ) the (non-unit) normal vector at the surface (resp.at the bottom) Moreover, we have Hence, to satisfy the divergence free condition, the velocity u should verify Dirichlet condition for the pressure.From the variational formulation (40)-(41) of the projection scheme, a natural boundary condition for the pressure is a Dirichlet condition.At Γ i , (i = in, out), p| Γi = p 0 then c(v, p) = Γi p 0 v 1 n ds and we take v ∈ V , q ∈ Q i with Neumann boundary for the pressure.The Neumann boundary condition for the projection scheme is not natural and to enforce such a condition, the elliptic problem (42)-(43) has to be considered.Taking now q ∈ Q sw , with Many studies have been done to choose an appropriate variational formulation for this problem.In [23] J-L.Guermond explores the different variational formulations in order to enforce a Neumann pressure boundary condition, in [25] some equivalent formulations are given to switch between Neumann and Dirichlet boundary conditions.
Taking the normal component at the boundary Γ i of the momentum equation, it follows that We note ∂H ∂n | Γi = β i , i = in, out.• In case β i = 0, a Neumann boundary condition for the pressure is deduced of a Dirichlet condition for u. ∂p ∂n • In the other cases, it gives a mixed boundary condition Then, in the two cases, we have imposed a Dirichlet velocity condition, that leads to take v ∈ V i and q ∈ Q, with for i = in, out Let us now give the coupling boundary conditions between the prediction step and the correction step.Indeed, in the projection part, boundary conditions need to be set in order to be consistent with the hydrostatic part.
Concerning the prediction step, we consider the well known Saint-Venant system and we assume that the Riemann invariant remains constant along the associated characteristic.This approach has been introduced in [10] and distinguishes fluvial and torrential boundaries depending on the Froude number F r = |u| c .Usual boundary conditions consist to impose a flux q 0 at the inflow boundary and a water depth at the outflow boundary.It is also classical to let a free outflow boundary, setting a Neumann boundary condition for the water depth and for the velocity.In both cases, we give the boundary conditions that have to be set in the correction step.We consider the first situation in which we set a flux at the inflow Γ in and a given depth at the outflow Γ out .Assuming a fluvial flow, this case consists in solving a Riemann problem at the interface Γ in where the global flux is given by q 0 = (q 01 , q 02 ) t = (Hu n+1/2 , Hw n+1/2 ) t .That gives the boundary values from the hyperbolic part.This leads to obtain a Dirichlet condition for the pressure at the left boundary of the correction part.Moreover, if H is given for the outflow, we preconize to give a mixed condition for the pressure that corresponds to the boundary condition ( 48) that leads to take u ∈ V out , with the definition (49) and We now consider the second situation in which we still impose a flux q 0 at the inflow and we set a free outflow boundary.In this case, we assume the two Riemann invariants are constant along the outgoing characteristics of the hyperbolic part (see [10]), therefore, we have a Neumann boundary condition for H n+1/2 and u n+1/2 .

∂H ∂n
Preserving these conditions at the correction step, it gives a Neumann boundary condition for the pressure of type (47) For an inflow given, the functional spaces will be defined by 3. Numerical approximation.

3.1.
Discretization.This section is devoted to the numerical approximation and mainly for the correction step.Let us be given a subdivision of Ω with N vertices x 1 < x 2 < ... < x N and we define the space step ∆x i+1/2 = x i+1 − x i .We also note Prediction part.For the prediction step (20) i.e the hydrostatic part of the model, we use a finite volume scheme.We introduce the finite volume cells C i centered at vertices x i such that Ω = ∪ i=1,N C i .Then, the approximate solution X n i at time t n is solution of the numerical scheme where σ n i = ∆t n ∆xi and F (resp.S) is a robust and efficient discretization of the conservative flux F (X) (resp.the source term S(X)).The time step is determined through a classical CFL condition.Many numerical fluxes and discretizations are available in the literature [8,20,29], we choose a kinetic based solver [5] coupled with the hydrostatic reconstruction technique [4].
Correction part.Concerning the correction step ( 21)-( 22), we consider the discrete problem corresponding to the mixed problem (40)-(41).We approach (V 0 ,Q) by the finite dimensional spaces (V 0h ,Q h ) and we note We also denote by (ϕ i ) i=1,N and (φ l ) l=1,M the basis functions of V 0h and Q h respectively.The finite dimensional spaces will be specified later on.We approximate (u, p) Therefore, we consider the discrete problem: Let us introduce the mass matrix M H given by and the two matrices B t , B defined by , and we denote Therefore, the problem (50)-(51) becomes Assuming that M H is invertible and eliminating the velocity U , we obtain the following equation that is a discretization of the elliptic equation (44) of Sturm-Liouville type governing the pressure p.
We now take into consideration the boundary conditions in the more general problem (45)-(46).The velocity u is approximated by u h ∈ V h , and the discrete problem is then written: Considering the matrix ∆t n C = (c (ϕ i , φ l )) 1≤l≤M,1≤i≤N that contains the boundary terms, the equation (52) becomes This approach is suitable for the finite element approximation that is given in the next section.However, it implies to inverse a mass matrix M H that is not diagonal and depends on the water depth H.In practice, we use the mass lumping technique introduced by Gresho ( [22]) to avoid inverting the mass matrix in projection methods for Navier-Stokes incompressible system.
3.2.Finite element P 1 /P 0 .In this part, the problem is solved by the mixed finite element approximation P 1 /P 0 (see [32]) on the domain Ω = ∪ M l=1 K l ( M = N −1 with N the number of nodes), where the velocity is approximated by a continuous linear function and the pressure is approximated by a discontinuous piecewise constant function over each element 1 , ∀l = 1, . . ., N − 1}, and Using the discretization given in 3.1, we denote by K i+1/2 the finite element [x i , x i+1 ], then the pressure is constant on the finite element K i+1/2 .For the sake of clarity, in this situation, let φ j+1/2 1≤j≤M be the basis functions for the pressure p h , and (ϕ i ) 1≤i≤N the basis functions for the velocity u h such that We note ζ = H +2z b and assume ζ is approximated by a piecewise linear function t , then the shallow water gradient operator is written .
Similarly, the shallow water divergence operator writes In one dimension, this approach corresponds to a staggered-grid finite-difference method where the velocity is computed at the nodes and the pressure is computed at the middle nodes.The discretization we obtain corresponds exactly to the finite difference scheme given in [1], and then, the properties established in [1] are conserved.
3.3.Finite element P 1 -iso-P 2 /P 1 .For the one dimensional P 1 -iso-P 2 /P 1 , we consider two meshes K h (the same as before) and K 2h with K h,i+1/2 = [x i , x i+1 ] and K 2h,j = [x 2j−1 , x 2j+1 ] the finite elements defined on the respective meshes K h and with N the total number of vertices of K h and M = (N − 1)/2 (assuming N odd), the number of vertices of K 2h .Therefore, the approximation spaces V h and Q h are defined by Then, the velocity and the pressure are written In figure 2, the dashed lines are the usual elementary basis functions of P 1 on the mesh K h , while the continuous lines are the basis functions on the mesh K 2h .We can define the divergence operator, for all j = 1, We use a linear interpolation for Hϕ j , and consider that ∆x i = ∆x ∀i = 1, . . ., N for the sake of simplicity.We still approximate ζ by ζ h defined before.The discrete shallow water divergence operator is computed for all nodes x j of the mesh K 2h and therefore, denoting i = 2j − 1, it can be written, ∀j = 1, Similarly, the gradient shallow water operator is obtained for all the nodes x i of the mesh K h .However, we distinguish the gradient at the nodes of the elements K 2h from the ones at the interior.In other words, for all the nodes x i of the mesh K 2h , the gradient operator is defined by On the other hand, for all the nodes x i such that i is even With the discretization of the shallow water operators given below, we are able to validate the scheme for the two first order methods.Nevertheless, notice that in the following section, only the first order has been implemented.

4.
Validation with an analytical solution.In [1], [12] some analytical solutions of the model (2)-( 5) have been presented and they allow to validate the numerical method.We consider the propagation of a solitary wave without topography.This solution has the form The solitary wave is a particular case where dispersive contributions are counterbalanced by non linear effects so that the shape of the wave remains unchanged during the propagation.The propagation of the solitary wave has been simulated for the parameters a = 0.4 m, H 0 = 1 m, and d = 1 m over a domain of 45 m with 9000 nodes.At time t = 0, the solitary wave is positioned inside the domain.The results presented in figure 3 show the different fields, namely the elevation, the components of velocity and the total pressure at different times, and the comparison with the analytical solution at the last time.In the projection step, the greatest difficulty is to compute the pressure corresponding to the boundary conditions of the hyperbolic part (as seen in 2.3).The solution near the boundary has been confronted to the analytical solution.In the following result, we set a Neumann boundary condition on the non hydrostatic pressure with the parameters given below.As shown in the figure 4, the pressure is well estimated at the outflow boundary and allows the wave to leave the domain with a good behavior.The inflow boundary condition has been tested with this same test case and gives similar results.We are able to let the solitary wave enter in the domain with a good approximation of the elevation.
The numerical simulations for the first order method are compared with the analytical solution and the L 2 -error has been evaluated over different meshes of sizes from 603 nodes to 6495 nodes (see figure 5).With the parameters given above, it gives a convergence rate close to 1 for the two computations, i.e P 1 /P 0 and P 1 -isoP 2 /P 1 .
Notice that the parameters set to validate the method lead to have a significant non hydrostatic pressure (see the figure 4) and then, the results show the ability of the method to preserve the solitary wave over the time.The numerical results have also been obtained for the Thacker's test presented in [1], with the same rate as the P 1 /P 0 method.Figure 5. Convergence rate of the solitary wave solution for P 1 /P 0 and P 1 -iso-P 2 /P 1 scheme.

Numerical results.
5.1.Dam break problem.We next study the dispersive effect on the classical dam break problem, which is usually modeled by a Riemann problem providing a left state (H L , u L ) and a right state (H R , u R ) on each side of the discontinuity x d ( [20]).However, our numerical dispersive model does not allow discontinuous solutions due to the variational spaces required for H (see also [12]), thus we provide an initial data numerically close to the analytical one To evaluate the non hydrostatic effect, the different fields have been compared with the shallow water solution with the initial data: m over a domain of length 600 m with 30000 nodes.In figure 6, the evolution of the state is shown at time t = 10 s and t = 45 s.The oscillations are due to the dispersive effects but the mean velocity does not change.These results are in adequation with the analysis proposed by Gavrilyuk in [28] for the Green-Naghdi model with the same configuration.

5.2.
Wet-dry interfaces.The ability to treat the wet/dry interfaces is crucial in geophysical problems, since geophysicists are interested in studying the behavior of the water-depth near the shorelines.This implies a water depth tending to zero at such boundaries.To treat the problem, we use the method introduced in [1], considering a minimum elevation H . Therefore, we confront the method with a coastal bottom at the right boundary over a domain of 35 m with 3000 nodes.A wave is generated at the left boundary with an amplitude of 0.2 m and an initial water depth H 0 = 1 m.In figure 7, the arrival of the wave at the coast is shown for times t = 7.91 s, 9.92 s and 10.42 s.  [18,17]) that consist in generating a small amplitude wave at the left boundary of a channel with topography as described in figure 8.At the left boundary, a wave is generated with a period T = 2.02 s and an amplitude of 0.02 m.A free outflow condition is set at the right boundary.The initial free surface is set to be η 0 = 0.4 m, and the measurement readings are saved at the following positions 10.5 m, 12.5 m, 13.5 m, 14.5 m ,15, 7 m and 17.3 m, placed at sensors 1 to 6 (fig.8 ).In such a situation, the non hydrostatic effects have a significant impact on the water depth that cannot be represented by a hydrostatic model.These effects result mainly from the slope of the bathymetry, 10% in this case.In the figure 9, the simulation has been run with the hydrostatic model and the elevation has been compared with measures at the sensor 5.As one can see, the non-hydrostatic pressure has to be taken in consideration to estimate the real water depth variation.
The numerical simulation with the non-hydrostatic model has been run with 15000 nodes on a domain of 49 m over 25 s and the comparisons are illustrated for each sensor (fig.10).The goal of this last result is also to highlight the ability of the model to capture dispersive effects for a geophysical flow with a non negligible pressure.
5.4.Remark on iterative method.We recall that this formulation should allow to extend the method on two dimensional unstructured grids.However, it requires to inverse a system at each time iteration, which will become too costly in two  dimensions.To anticipate the two dimensional problem, this method has been tested using different iterative methods like conjugate gradient and Uzawa methods.In figure 11, we show a comparison of the computing time for the implementation P 1 -iso-P 2 /P 1 and Uzawa method.In one dimension, it is not relevant to use one of these methods, while it will be necessary for the two dimension model.

Conclusion.
In this paper, a variational formulation has been established for the one dimensional dispersive model introduced in [12].The main idea is to give a new framework in which it will be possible to extend the scheme to the two Comparison of the computing time (CPU) for the P 1 /P 0 scheme, the P 1 -iso-P 2 /P 1 and Uzawa method with P 1 -iso-P 2 /P 1 solver with a tolerance 10 −5 .dimensional model.To this aim, the finite-element method has been presented with two approximation spaces.First, the P 1 /P 0 approximation has been done and we recover, as expected, the finite difference scheme, together with the good results proved in [1].Then, the P 1 -iso-P 2 /P 1 approximation has been studied to prepare the two dimensional problem.We have validated the method using several numerical tests and studying the dispersive effect on geophysical situations.

1 .
Introduction.Starting from the incompressible Euler or Navier-Stokes system, the hydrostatic assumption consists in neglecting the vertical acceleration of the fluid.More precisely -and with obvious notations -the momentum along the vertical axis of the Euler equation ∂w ∂t

Figure 2 .
Figure 2. Representation of the basis functions.

Figure 3 .
Figure 3. Propagation of the solitary wave at times 1.00008 s, 1.9009 s, 3.9017 s and 5.9025 s.Comparison with analytical solution at time t = 5.9025 s.

Figure 6 .
Figure 6.The dam break problem, elevation H and velocity u at times t = 10 s and t = 45 s.

Figure 7 .
Figure 7. Propagation of a wave at a wet/dry interface.

Figure 10 .Figure 11 .
Figure 10.Comparison between measured and computed elevations on Dingemans test for the first six sensors.