A Duran

F Marche

Recent advances on the discontinuous Galerkin method for shallow water equations with topography source terms

We consider in this work the discontinuous Galerkin discretization of the nonlinear Shallow Water equations on unstructured triangulations. In the recent years, several improvements have been made in the quality of the discontinuous Galerkin approximations for the Shallow Water equations. In this paper, we first perform a review of the recent methods introduced to ensure the preservation of motionless steady states and robust computations. We then suggest an efficient combination of ingredients that leads to a simple high-order robust and well-balanced scheme, based on the alternative formulation of the equations known as the pre-balanced shallow water equations. We show that the preservation of the motionless steady states can be achieved, for an arbitrary order of polynomial expansion. Additionally, the preservation of the positivity of the water height is ensured using the recent method introduced in [J. Comput.Phys, 50, pp 29-62, 2012]. An extensive set of numerical validations is performed to highlight the efficiency of these approaches. Some accuracy, CPU-time and convergence studies are performed, based on comparisons with analytical solutions or validations against experimental data, for several test cases involving steady states and occurrence of dry areas. Some comparisons with a recent finite-volume MUSCL approach are also performed.

Introduction

The non-linear shallow water equations (SWE in the following) model the dynamic of a free surface shallow layer of homogeneous incompressible fluid. They are used to describe vertically averaged flows in terms of horizontal velocity and depth variation. This set of equations, which can be obtained by asymptotic analysis and depthaveraging of the Navier-Stokes equations [START_REF] De Saint-Venant | Théorie du mouvement non-permanent des eaux, avec application aux crues des rivières et à l'introduction des marées dans leur lit[END_REF][START_REF] Gerbeau | Derivation of viscous saint-venant system for laminar shallow water; numerical validation[END_REF][START_REF] Marche | Derivation of a new two-dimensional viscous shallow water model with varying topography, bottom friction and capillary effects[END_REF], is well-suited for the simulation of geophysical phenomena, such as river and oceanic flows, or even avalanches with suitable source terms. This model is also extensively used in coastal engineering, for the study of nearshore flows involving run-up and run-down on sloping beaches or coastal structures. To allow a proper modelization of such a variety of phenomena, accurate and robust numerical methods have to be considered. Assuming a smooth parameterization of the topography z : R 2 → R, the SWE are defined as follows :

U t + ∇ • H(U, z) = S(U, z) , (1)
with

U =   h q x q y   , H(U) =   q x q y uq x + 1 2 gh 2 vq x uq y vq y + 1 2 gh 2   , S(U, z) =   0 -ghz x -ghz y   ,
where h stands for the water height, u = (u, v) for the horizontal velocity and q = (q x , q y ) = (hu, hv) for the horizontal discharge. Denoting by Θ the convex set of admissible states, defined by:

Θ = (h, hu, hv) ∈ R 3 ; h ≥ 0, (u, v) ∈ R 2 , (2)
we denote U : R 2 × R + → Θ the vector of conservative variables, H : Θ → R 3 the flux function and S : Θ × R → R 3 the topography source term.

Nowadays, a large variety of numerical models are able to produce accurate approximations of weak solutions of [START_REF] Aizinger | A discontinuous Galerkin method for two-dimensional flow and transport in shallow water[END_REF]. Finite Volume (FV in the following) methods are known to be very efficient, notably for their low computational cost and their capability in capturing shocks, see for instance [START_REF] Berthon | A positive preserving high order VFRoe scheme for shallow water equations: a class of relaxation schemes[END_REF][START_REF] Delestre | A numerical scheme for a viscous shallow water model with friction[END_REF][START_REF] Dutykh | The volna code for the numerical modelling of tsunami waves: generation, propagation and inundation[END_REF][START_REF] Gallouët | Some approximate godunov schemes to compute shallow-water equations with topography[END_REF][START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF][START_REF] Guo | A characteristic-based finite volume scheme for shallow water equations[END_REF][START_REF] Hubbard | A 2d numerical model of wave run-up and overtopping[END_REF][START_REF] Nikolos | An unstructured node-centered finite volume scheme for shallow water flows with wet/dry fronts over complex topography[END_REF][START_REF] Wang | A 2d shallow flow model for practical dam-break simulations[END_REF] and also some references herein. However, basic FV methods usually offer low accuracy and one generally needs to use some reconstruction methods to offset the low order of convergence and the diffusive losses (see [START_REF] Jiang | Efficient implementation of Weighted ENO schemes[END_REF][START_REF] Van-Leer | Towards the ultimate conservative difference scheme. V -a second-order sequel to Godunov's method[END_REF][START_REF] Perthame | A variant of Van Leer's method for multidimensional systems of conservation laws[END_REF][START_REF] Wolf | High-order ENO and WENO schemes for unstructured grids[END_REF] for instance). Discontinuous Galerkin (DG in the following) methods have raised great interest during the past twenty years. These methods combine the background of the Finite Element methods, FV methods and Riemann solvers, taking into account the physic of the problem. An arbitrary order of accuracy can be obtained with the use of high-order polynomials within elements and they are able to handle complex geometries with the use of unstructured meshes. They are highly parallelizable, and exhibit nice strong stability properties. The reader is referred to the following pioneering works [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V[END_REF][START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF] for the general backgrounds. It is only recently that the DG approach has been applied to the SWE and we can find a growing number of studies, including flows with shocks, such as dam-break and oblique hydraulic jumps ( [START_REF] Li | The discontinuous Galerkin finite element method for the 2d shallow water equations[END_REF][START_REF] Schwanenberg | Discontinuous Galerkin finite-element method for transcritical twodimensional shallow water flows[END_REF][START_REF] Tassi | Space discontinuous Galerkin method for shallow water flowskinetic and hllc flux, and potential vorticity generation[END_REF]). Several approaches involving arbitrary orders methods on unstructured triangulations have been developed for the SWE [START_REF] Eskilsson | A triangular spectral/hp discontinuous Galerkin method for modelling 2d shallow water equations[END_REF][START_REF] Kubatko | hp Discontinuous Galerkin methods for advection dominated problems in shallow water flow[END_REF], possibly with dynamic p-adaptivity [START_REF] Kubatko | Dynamic p-adaptive Runge-Kutta discontinuous Galerkin methods for the shallow water equations[END_REF], adaptive refinement [START_REF] Remacle | An adaptive discretization of shallow-water equations based on discontinuous Galerkin methods[END_REF], discretizations of the viscous SWE relying on a Local Discontinuous Galerkin (LDG) treatment of the second order derivatives [START_REF] Aizinger | A discontinuous Galerkin method for two-dimensional flow and transport in shallow water[END_REF][START_REF] Dawson | A discontinuous Galerkin method for three-dimensional shallow water equations[END_REF], discretizations of the equations on the sphere [START_REF] Bernard | High-order discontinuous Galerkin schemes on general 2D manifolds applied to the shallow water equations[END_REF][START_REF] Giraldo | Nodal high-order discontinuous Galerkin methods for the spherical shallow water equations[END_REF][START_REF] Lauter | A discontinuous Galerkin method for the shallow water equations in spherical triangular coordinates[END_REF][START_REF] Nair | A Discontinuous Galerkin Global Shallow Water Model[END_REF], or even space-time dG methods [START_REF] Ambati | Space-time discontinuous Galerkin finite element method for shallow water flows[END_REF]. The list is of course non-exhaustive.

More recently, several authors have focused on two interesting issues, which are particularly relevant in many applications: the preservation of the motionless steady states, and the preservation of the water height positivity, to properly handle flooding and drying events. The paper is organized as follows: in the next section, we propose a review of some of the existing methods recently introduced to satisfy these two properties. We also give a review of the main limiting technics introduced to handle discontinuities and prevent the generation of spurious oscillations. In §3, we study a new combination of ingredients that lead to an arbitrary order robust and well-balanced nodal discontinuous-Galerkin discretization of the SWE on unstructured meshes, relying on the so called pre-balanced SWE [START_REF] Kurganov | Central-upwind schemes for the saint-venant system[END_REF][START_REF] Rogers | Adaptive Q-tree Godunov-type scheme for shallow water equations[END_REF] (PBSW equations in the following) and the recent method introduced in [START_REF] Xing | Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations[END_REF][START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high order discontinuous Galerkin schemes for conservation laws on triangular meshes[END_REF]. A local limiting process, allowing the possible occurrence of shocks and contact discontinuities, is also described based on [START_REF] Burbeau | A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods[END_REF]. In §3.8, we establish the main well-balancing and robustness properties of this combination. §4 is devoted to extensive numerical validations in the case of second and third order schemes, including convergence and accuracy analysis, CPU time studies, and comparisons with analytical solutions and experimental data for cases involving steady states preservation and occurrence of dry areas. Some comparisons with the MUSCL FV scheme of [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF] are also performed.

A survey of existing methods

Well-balancing

This property is often referred to as C-property, following [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF]. Defining the free surface η = h + z, we say that the motionless equilibrium states are preserved by a given numerical scheme if the following property holds for all n ∈ N :

η n h ≡ η e q n h ≡ 0 ⇒ η n+1 h ≡ η e q n+1 h ≡ 0 , (3)
where η e is a constant, and w h = (η h , q h ) is the discrete solution produced by the numerical scheme. Nowadays, a large number of FV approaches are able to offer such a property, see among them [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF][START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF][START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF][START_REF] Gallouët | Some approximate godunov schemes to compute shallow-water equations with topography[END_REF][START_REF] Gosse | A well-balanced flux-vector splitting scheme designed for hyperbolic systems of conservation laws with source terms[END_REF][START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF][START_REF] Jin | A steady-state capturing method for hyperbolic systems with geometrical source terms[END_REF][START_REF] Katsaounis | Upwinding sources at interfaces in conservation laws[END_REF][START_REF] Kurganov | Central-upwind schemes for the saint-venant system[END_REF][START_REF] Kurganov | A second-order well-balanced positivity preserving central-upwind scheme for the Saint-Venant system[END_REF][START_REF] Leveque | Balancing source terms and flux gradients in high-resolution godunov methods: the quasi-steady wave-propagation algorithm[END_REF][START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF][START_REF] Lukacova | Well-balanced finite volume evolution Galerkin methods for the shallow water equations[END_REF][START_REF] Russo | Central schemes for conservation laws with application to shallow water equations[END_REF][START_REF] Zhou | The surface gradient method for the treatment of source terms in the shallow-water equations[END_REF] for first and second order accuracy well-balanced FVM and [START_REF] Caleffi | Fourth-order balanced source term treatment in central WENO schemes for shallow water equations[END_REF][START_REF] Castro | Well-balanced high order extensions of godunov's method for semilinear balance laws[END_REF][START_REF] Noelle | Well-balanced finite volume schemes of arbitrary order of accuracy for shallow water flows[END_REF][START_REF] Rnjariic | Balanced finite volume WENO and central WENO schemes for the shallow water and the openchannel flow equations[END_REF][START_REF] Xing | High order finite difference WENO schemes with the exact conservation property for the shallow water equations[END_REF][START_REF] Xing | A new approach of high order well-balanced finite volume WENO schemes and discontinuous Galerkin methods for a class of hyperbolic systems with source terms[END_REF] for some higher order schemes. The development of well-balanced DG schemes for the SWE is recent, and there is very few existing works, especially when considering the 2D case on unstructured grids. In [START_REF] Rhebergen | Discontinuous galerkin finite element methods for hyperbolic nonconservative partial differential equations[END_REF], general space and space-time DG formulations are introduced for hyperbolic nonconservative partial differential equations, and applications are performed for the one-dimensional SWE with topography, regarding the topography as an additional variable in the spirit of [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF].

The resulting space method is shown to preserve the C-property.

A well-balanced method is developed for second order accuracy DG schemes in [START_REF] Kesserwani | Well-balancing issues related to the RKDG2 scheme for the shallow water equations[END_REF][START_REF] Kesserwani | A discontinuous Galerkin algorithm for the two-dimensional shallow water equations[END_REF][START_REF] Kesserwani | Locally limited and fully conserved rkdg2 shallow water solutions with wetting and drying[END_REF], for the 1D and 2D case on rectangular meshes, using the PBSW equations and borrowing some ideas coming from the hydrostatic reconstruction [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF] and [START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF].

In [START_REF] Ern | A well-balanced Runge-Kutta discontinuous Galerkin method for the shallow-water equations with flooding and drying[END_REF], the well-balancing is ensured for polynomial expansions of arbitrary orders and on unstructured meshes, using the ideas of the hydrostatic reconstruction. Non-negative reconstructions of the water height are introduced element-wise, for each edge, together with an additional flux modification term directly accounted for in the weak formulation.

Recently, a more general approach has been introduced in [START_REF] Xing | High order well-balanced finite volume WENO schemes and discontinuous Galerkin methods for a class of hyperbolic systems with source terms[END_REF] for a general class of conservation law with separable source terms, leading to a class of high-order DG methods with the well-balanced property. The key ingredient is a suitable decomposition of the integral of the source terms into a sum of several terms, each of which is discretized independently in a consistent way with the discretization of the corresponding flux derivative terms. This ensures the well-balancing and preserves the high-accuracy of the method. However, the overall algorithm is quite complex both to understand and implement.

An easier and less computationally expensive approach is subsequently introduced by the same authors in [START_REF] Xing | A new approach of high order well-balanced finite volume WENO schemes and discontinuous Galerkin methods for a class of hyperbolic systems with source terms[END_REF], based on a generalization of the methods introduced in [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF][START_REF] Noelle | Well-balanced finite volume schemes of arbitrary order of accuracy for shallow water flows[END_REF].

To conclude this section, let us mention the difficult issue of well-balancing for general moving water steady states. Very few high-order accuracy methods are able to maintain such general steady steady states, see for instance the recent works [START_REF] Castro | High order exactly well-balanced numerical methods for shallow water systems[END_REF][START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equation with moving water[END_REF]. However, it is shown in [START_REF] Xing | On the advantage of well-balanced schemes for moving-water equilibria of the shallow water equations[END_REF] that moving-water well-balanced schemes exhibit some advantages when compared to motionless steady states preserving methods. Very recently, a positivity-preserving high-order well-balanced DG scheme for the SWE that preserves moving water steady states is developed in [START_REF] Xing | Exactly well-balanced discontinuous Galerkin methods for the shallow water equations with moving water equilibrium[END_REF] for the 1d case, and this is, to our knowledge, the first dG method that addresses this issue.

Robustness

The second property aims at numerically preserving the convex set of admissible states (2). This property is fundamental to simulate a large range of phenomena involving a possibly moving wet/dry interface, like for dambreaks or propagating waves reaching the shore in coastal engineering. In the context of FV methods, several methods have been introduced since the past years to preserve the positivity of the water height. Such robust schemes usually rely on positive preserving Riemann solvers, like the HLL and HLLC solvers [START_REF] Harten | On upstream differencing and godunov-type schemes for hyperbolic conservation laws[END_REF][START_REF] Toro | Restoration of the contact surface in the HLL riemann solver[END_REF], the kinetic solver [START_REF] Perthame | Kinetic Formulations of Conservation Laws[END_REF] or the relaxation-VFRoe solver [START_REF] Berthon | A positive preserving high order VFRoe scheme for shallow water equations: a class of relaxation schemes[END_REF]. The main issue is to preserve this feature while introducing the wellbalancing discretization of the source terms and the higher-order accuracy, possibly on 2D unstructured meshes, see for instance [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF][START_REF] Audusse | A well-balanced positivity preserving "second-order" scheme for shallow water flows on unstructured meshes[END_REF][START_REF] Berthon | Efficient well balanced hydrostatic upwind schemes for shallow water equations[END_REF][START_REF] Bryson | Well-balanced positivity preserving central-upwind scheme on triangular grids for the Saint-Venant system[END_REF][START_REF] Castro | Well-balanced numerical schemes based on a generalized hydrostatic reconstruction technique[END_REF][START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF][START_REF] Gallardo | On a well-balanced high-order finite volume scheme for shallow water equations with topography and dry areas[END_REF][START_REF] Kurganov | A second-order well-balanced positivity preserving central-upwind scheme for the Saint-Venant system[END_REF][START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF][START_REF] Marche | Evaluation of well-balanced bore-capturing schemes for 2d wetting and drying processes[END_REF][START_REF] Perthame | On positivity preserving finite volume schemes for Euler equations[END_REF]. Note that even for FV schemes, the issue of robustness for higher order schemes on unstructured meshes is far from trivial, see [START_REF] Berthon | Robustness of MUSCL schemes for 2D unstructured meshes[END_REF][START_REF] Perthame | On positivity preserving finite volume schemes for Euler equations[END_REF]. Again, some of the ideas developed for the FV methods have been adapted to the DG approach. A popular approach, adapted for P 1 polynomial expansions basically consists in locally modifying the slope of the linear expansions if negative values of the water height appear. The local gradients are moderated element-wise until such negative values are avoided, and usually until the water height values are above a tolerance. Provided that the whole Runge-Kutta DG scheme preserves the positivity of the average water height, this approach leads to a robust scheme. This strategy can of course be applied to higher order DG schemes, providing that the discrete solution is first locally projected onto linear polynomials, entailing some loss of accuracy. In [START_REF] Kesserwani | A discontinuous Galerkin algorithm for the two-dimensional shallow water equations[END_REF][START_REF] Kesserwani | Locally limited and fully conserved rkdg2 shallow water solutions with wetting and drying[END_REF] such ideas are used in combination with some positive reconstructions of the water height, inspired from the hydrostatic reconstruction [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF], for 1D and 2D case on rectangular meshes, for the second order accuracy DG scheme (linear polynomial expansion of the solution).

In [START_REF] Ern | A well-balanced Runge-Kutta discontinuous Galerkin method for the shallow-water equations with flooding and drying[END_REF], this strategy is adapted to the unstructured meshes framework. Considering a P 1 expansion of the water height, the idea is to locally post-process the gradients of h in elements such that a positive mean value is observed, but with occurrences of negative values at one or two vertices. The gradients are modified in such a way that the water height vanishes at such vertices. However, the resulting scheme can possibly lead to negative mean water height values, which are arbitrarily set to zero. A similar strategy is introduced in [START_REF] Bunya | A wetting and drying treatment for the Runge-Kutta discontinuous Galerkin solution to the shallow water equations[END_REF], but with a posteriori modifications of the water height that preserve the overall accuracy and the mass and momentum conservation properties. Additionally, a sufficient condition ensuring the positivity of the mean water depth in each element is provided. Alternatively, a strategy allowing robust wetting and drying is introduced in [START_REF] Bokhove | Flooding and drying in discontinuous Galerkin finite-element discretizations of shallow-water equations. part 1: One dimension[END_REF] for space-time DG elements, allowing to identify and accurately discriminate the wet and dry areas and moving the mesh accordingly. The reader is also referred to [START_REF] Schwanenberg | Discontinuous Galerkin finite-element method for transcritical twodimensional shallow water flows[END_REF][START_REF] Trahan | Local time-stepping in Runge-Kutta discontinuous Galerkin finite element methods applied to the shallow-water equations[END_REF] for some other recent studies. More generally, the preservation of the water height positivity for an arbitrary order DG method is not a simple problem and a general method was recently introduced in [109]. This method ensures an accuracy preserving maximum-principle property for DG schemes of arbitrary order, in a general scalar conservation laws framework, under a suitable CFL-like condition. This method is extended to triangular meshes in [START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high order discontinuous Galerkin schemes for conservation laws on triangular meshes[END_REF], to the SWE in 1D and 2D on rectangular meshes in [START_REF] Xing | Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations[END_REF], and very recently to the SWE on triangular meshes in [START_REF] Xing | Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water equations on unstructured triangular meshes[END_REF]. The method, based on the positivity of the associated first order scheme and a simple positivity and accuracy preserving limiter, is detailed in §3. Note that the method introduced in [106] also preserves the C-property.

Limiting strategies

In the context of hyperbolic conservation laws, special interest has to be given to the way discontinuities are handled. High-order numerical schemes produce spurious oscillations near discontinuities, which may, indeed, lead to nonphysical solutions (like negative water height), numerical instabilities and unbounded computational solutions.

A very popular approach is to use a limiting procedure which is capable to detect and control the high variations of the approximated solution, leading to local modifications of the approximate solution at each time step. A very popular approach is the TVB-generalized slope-limiter technic introduced in [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V[END_REF], which is used in many studies, see [START_REF] Bunya | A wetting and drying treatment for the Runge-Kutta discontinuous Galerkin solution to the shallow water equations[END_REF][START_REF] Ern | A well-balanced Runge-Kutta discontinuous Galerkin method for the shallow-water equations with flooding and drying[END_REF][START_REF] Kesserwani | A discontinuous Galerkin algorithm for the two-dimensional shallow water equations[END_REF][START_REF] Li | The discontinuous Galerkin finite element method for the 2d shallow water equations[END_REF][START_REF] Rhebergen | Discontinuous galerkin finite element methods for hyperbolic nonconservative partial differential equations[END_REF] for instance. This approach, relying on a modified minmod function, maintains the conservation of mass element by element, does not degrade the accuracy of the method for piecewise linear approximations and does not flatten smooth extrema. The higher order modes are set to zero when limiting is needed. An interesting approach is proposed in [START_REF] Biswas | Parallel, adaptive finite element methods for conservation laws[END_REF] to the case of higher-order slope limiting, which consists in adaptively and successively differentiating the numerical solution. The derivations raise linear terms that can be limited as in the case of a linear approximation. An efficient parameter-free (and consequently problem independent) alternative, relying on the maxmod function is proposed in [START_REF] Burbeau | A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods[END_REF] and is detailed in §3.5. Of course, many slope limiters used within the FV methods can be adapted to the needs of the DG method, like the van Albada type limiting method [START_REF] Tu | A slope limiting procedure in discontinuous galerkin finite element method for gasdynamic applications[END_REF] in which the gradient in an element is limited using the weighted average of face gradients. Let us also mention [START_REF] Chavent | Mathematical Models and Finite Elements for Reservoir Simulation[END_REF][START_REF] Hoteit | New two-dimensional slope limiters for discontinuous Galerkin methods on arbitrary meshes[END_REF] for some extensions of van Leer's slope limiter for two-dimensional DG method and relying on the solution of a least squares problem.

To deal with spurious oscillations around discontinuities, an alternative approach is to introduce a dissipation operator, see for instance [START_REF] Van Der | Space-time discontinuous Galerkin finite element method with dynamic grid motion for inviscid compressible flows[END_REF]. This can be done by adding the operator into the weak formulation, see [START_REF] Ambati | Space-time discontinuous Galerkin finite element method for shallow water flows[END_REF][START_REF] Tassi | Space discontinuous Galerkin method for shallow water flowskinetic and hllc flux, and potential vorticity generation[END_REF]. However, the main drawback of these approaches is that they can decrease the high order of accuracy of the method, when used in smooth regions of the solution. To avoid this, an important issue is to determine the area of application of the limiting procedure (or the smoothing operator) to avoid a possible loss of accuracy. A popular method, widely used in CFD, is to use a discontinuity detector. Using such a detector, it is possible to limit spurious oscillations only near such discontinuities and the high order of accuracy can be preserved away from discontinuities. This has been investigated in [START_REF] Qiu | A comparison of troubled-cell indicators for Runge-Kutta discontinuous Galerkin methods using weighted essentially nonoscillatory limiters[END_REF], where several discontinuity detectors are reviewed and compared. In the framework of DG method, a very popular approach is the one introduced in [START_REF] Krivodonova | Limiters for high-order discontinuous galerkin methods[END_REF][START_REF] Krivodonova | Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws[END_REF]. They developed a discontinuity detector based on a strong super-convergence properties at the outflow boundaries of each element in smooth regions. This method has been proved to be efficient and is used in several studies, see for instance [START_REF] Ambati | Space-time discontinuous Galerkin finite element method for shallow water flows[END_REF][START_REF] Ern | A well-balanced Runge-Kutta discontinuous Galerkin method for the shallow-water equations with flooding and drying[END_REF][START_REF] Kesserwani | Locally limited and fully conserved rkdg2 shallow water solutions with wetting and drying[END_REF][START_REF] Remacle | An adaptive discretization of shallow-water equations based on discontinuous Galerkin methods[END_REF][START_REF] Tassi | Space discontinuous Galerkin method for shallow water flowskinetic and hllc flux, and potential vorticity generation[END_REF][START_REF] Zhu | Runge-Kutta discontinuous Galerkin method using a new type of WENO limiters on unstructured meshes[END_REF].

In [START_REF] Qiu | Runge-Kutta discontinuous Galerkin method using WENO limiters[END_REF][START_REF] Zhu | Runge-Kutta discontinuous Galerkin method using WENO limiters II: Unstructured meshes[END_REF], an alternative approach is proposed, which uses traditional ENO or WENO methodology as limiters for the DG methods, maintaining the high-order of accuracy. The key ideas are first to identify cells which might need some limiting procedure and then to replace the solution polynomials in those cells with reconstructed polynomials, maintaining the original cells mean values and preserving orders of accuracy as before while being less oscillatory. However, this approach relies on a wide stencil, especially for high order of accuracy, even with the use of Hermite-WENO based limiters [START_REF] Luo | A Hermite WENO-based limiter for discontinuous Galerkin method on unstructured grids[END_REF][START_REF] Qiu | Hermite WENO schemes and their application as limiters for Runge-Kutta discontinuous Galerkin method: one dimensional case[END_REF]. To overcome these drawbacks, a new WENO limiting procedure on structured and unstructured meshes is recently introduced in [START_REF] Zhong | A simple weighted essentially nonoscillatory limiter for Runge-Kutta discontinuous Galerkin methods[END_REF][START_REF] Zhu | Runge-Kutta discontinuous Galerkin method using a new type of WENO limiters on unstructured meshes[END_REF]. The main idea is to reconstruct the entire polynomial, instead of reconstructing point values or moments in the classical WENO reconstructions. A smaller stencil is needed, without extensive usage of meshes geometric informations, leading to a simpler implementation and very efficient results.

3 Pre-Balanced RKDG scheme

In this section, we illustrate the efficiency of a particular combination of selected ingredients, leading to a simple well-balanced and robust DG 2D discretization on unstructured triangulations. Compared to [START_REF] Ern | A well-balanced Runge-Kutta discontinuous Galerkin method for the shallow-water equations with flooding and drying[END_REF][START_REF] Xing | Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations[END_REF][START_REF] Xing | Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water equations on unstructured triangular meshes[END_REF], we make a slightly different choice as we use the PBSW equations instead of the classical SWE, see §3.2 for details. This choice is motivated by the authors previous works [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF][START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF]. As shown in §3.8, the use of this formulation can also lead to the preservation of motionless steady states, for polynomial expansions of arbitrary orders, providing that the interface fluxes are slightly modified, adapting the ideas of [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF] and [START_REF] Xing | Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations[END_REF].

As for the limiting strategy, while the TVB-generalized limiter [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V[END_REF] is widely used, we choose to combine the discontinuity detector [START_REF] Krivodonova | Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws[END_REF] and the moment limiting approach of [START_REF] Burbeau | A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods[END_REF]. Some comparison with the TVB-generalized limiter [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V[END_REF] are performed in the next section. The robustness of the resulting scheme is ensured, adapting the recent approach of [START_REF] Xing | Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations[END_REF][START_REF] Xing | Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water equations on unstructured triangular meshes[END_REF] to the pre-balanced formulation.

Settings and notations

Let T h be a partition of the computational domain Ω into N t triangular elements T l , 1 ≤ l ≤ N t . The element T l has a boundary denoted ∂T l , a unit outward normal nl , an area |T l |, a diameter d l (defined as the length of the largest edge), and a perimeter p l . Let x denote the (x, y) coordinate in Ω. Additionally, we define the following notations, relative to the neighborhood of T l :

• σ(k): the index of the element T σ(k) neighboring T l , for each k ∈ {1, . . . , 3},

• Γ lσ : the boundary interface defined by T l and T σ , so that

∂T l = 3 ∪ k=1 Γ lσ(k) ,
• ℓ lσ : the length of Γ lσ , so that

p l = 3 k=1 ℓ lσ(k) ,
• nlσ : the unit normal to Γ lσ , pointing to T σ .

We aim at computing an approximate vector solution on this triangulation, denoted w h = (η h , q h ). Let us define

V h := {v ∈ L 2 (Ω) | ∀ T ∈ T h , v |T ∈ P d (T )}, (4)
where P d (T ) denotes the space of 2-variables polynomials in T of degree at most d. For the sake of simplicity, it is assumed that Ω is a polygonal domain in two space dimensions, so that T h covers Ω exactly.

PBSW equations

This alternate formulation of the SWE is introduced in [START_REF] Rogers | Adaptive Q-tree Godunov-type scheme for shallow water equations[END_REF] with the purpose of naturally balancing flux gradient and topography source term, using the approximate Riemann solver of Roe. A similar set of equations is separately used in [START_REF] Kurganov | Central-upwind schemes for the saint-venant system[END_REF] and [START_REF] Russo | Central schemes for conservation laws with application to shallow water equations[END_REF] as a basis to introduce central-upwind schemes. Such SWE formulated in terms of free-surface are widely used within the FV framework, see for instance [START_REF] Bryson | Well-balanced positivity preserving central-upwind scheme on triangular grids for the Saint-Venant system[END_REF][START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF][START_REF] Kurganov | Central-upwind schemes for the saint-venant system[END_REF][START_REF] Kurganov | A second-order well-balanced positivity preserving central-upwind scheme for the Saint-Venant system[END_REF][START_REF] Rogers | Adaptive Q-tree Godunov-type scheme for shallow water equations[END_REF][START_REF] Rogers | Mathematical balancing of flux gradient and source terms prior to using Roe's approximate Riemann solver[END_REF][START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF][START_REF] Russo | Central schemes for conservation laws with application to shallow water equations[END_REF], but far less within the DG approach (to our knowledge only in [START_REF] Kesserwani | Well-balancing issues related to the RKDG2 scheme for the shallow water equations[END_REF] and subsequent works but for a second order DG method in 1D or 2D on rectangular meshes). These equations read as follows:

W t + ∇ • H(W, z) = S(W, z) , (5)
with :

W =   η q x q y   , S(W, z) =   0 -gηz x -gηz y   , H(W, z) =   q x q y uq x + 1 2 g(η 2 -2ηz) vq x uq y vq y + 1 2 g(η 2 -2ηz)   .
Remark 3.1. In [START_REF] Kurganov | Central-upwind schemes for the saint-venant system[END_REF][START_REF] Russo | Central schemes for conservation laws with application to shallow water equations[END_REF], the modified SWE are obtained from the classical SWE by replacing h by η -z. In [START_REF] Rogers | Adaptive Q-tree Godunov-type scheme for shallow water equations[END_REF], the SWE are written in terms of free surface elevation ζ above the still water depth h s . Then, the PBSW equations are obtained by splitting and redistributing the gh∇ζ term in a novel way, which preserves the hyperbolic nature of the SWE:

gh∇ζ = 1 2 g∇(ζ 2 + 2ζh s ) -gζ∇h s .
This approach has been generalized in [START_REF] Rogers | Mathematical balancing of flux gradient and source terms prior to using Roe's approximate Riemann solver[END_REF], in which a reformulation is introduced for general conservation laws with source terms, exploiting the deviations from the system's unforced equilibrium. However, ζ and h s can be undefined in dry areas. To overcome this, an alternative splitting of the free surface gradient term, relying on the total free surface elevation η = h + z is derived and applied in [START_REF] Liang | Adaptive quadtree simulation of shallow flows with wet-dry fronts over complex topography[END_REF][START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF]:

gh∇ζ = 1 2 g∇(η 2 -2ηz) + gη∇z.
This is the formulation used in the following.

Discrete formulation

A weak formulation of the problem is obtained by multiplying ( 5) by a test function φ h ∈ V h . The result is integrated on a given element T l and the flux term is integrated by part to obtain :

T l ∂ ∂t W(x, t)φ h (x)dx - T l H(W, z).∇φ h (x)dx + ∂T l H(W, z) • nl φ h (s)ds = T l S(W, z)φ h (x)dx. (6)
The local approximated vector solution w h|T l ∈ (V h ) 3 is expressed as a polynomial of order d on each T l :

w h|T l (x, t) = N d i=1 wl i (t)θ l i (x), ∀x ∈ T l , ∀t ∈ [0, t max ], (7)
where {θ l i } N d i=1 is a polynomial expansion basis for P d (T l ), and { wl i (t)} N d i=1 are the local expansion coefficient vectors with wl i (t) = t (η l i , ( qx ) l i , ( qy ) l i ). Many choices are possible for the expansion basis, and we choose in the following to use a nodal approach: {θ l i } N d i=1 will refer to the interpolant Lagrangian expansion basis, with N d = (d + 1)(d + 2)/2. Let also consider a polynomial expansion of the topography parameterization z :

z h|T l (x) = N d i=1 zl i θ l i (x), ∀x ∈ T l . (8)
Thus, a discrete formulation of ( 6) is obtained by replacing the exact solution W(x, t) by the approximation w h (x, t) and the test function φ h by each element of the expansion basis, successively :

T l N d i=1 d dt wl i (t)θ l i (x) θ l j (x)dx - T l H(w h , z h ).∇θ l j (x)dx + ∂T l H(w h , z h ) • nl θ l j (s)ds = T l S(w h , z h )θ l j (x)dx, 1 ≤ j ≤ N d . (9)
Noting that we have

∂T l H(w h , z) • nl θ l j (s)ds = 3 k=1 Γ lσ(k) H lσ(k) θ l j (s)ds ,
where we have set

H lσ(k) = H(w h , z h ) • nlσ(k) ,
we obtain the following semi-discrete formulation, where the space dependency of the basis functions is omitted, for the sake of simplicity :

N d i=1 d dt wl i (t) T l θ l i θ l j dx - T l H(w h , z h ) • ∇θ l j dx + 3 k=1 Γ lσ(k) H lσ(k) θ l j ds = T l S(w h , z h )θ l j dx, 1 ≤ j ≤ N d . ( 10
)
Remark 3.2. As we use the Lagrangian expansion basis, the expansion coefficients in ( 7) and ( 8) can be regarded as the nodal values at the at the corresponding N p nodes. Consequently, the topography parameterization expansion coefficients in ( 8) are obtained by "reading" the value of the topography at these nodes, for each element.

Numerical flux

Classically, since matching conditions are not enforced on the approximated vector solution w h (x, t) at element interfaces, the boundary flux H(w h , z h ) • n is not uniquely defined. We propose in the following a simple choice for the interfaces numerical fluxes H lσ(k) , inspired from the Finite-Volume well-balanced discretization proposed in [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF], that leads to a well-balanced scheme that preserves motionless steady states. This modified flux can also be seen as the adaptation of the ideas of [START_REF] Xing | Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations[END_REF] to the pre-balanced formulation (5). Let us denote, for a given interface Γ lσ(k) , w - k and w + k respectively the restrictions of w h|T l and w h|T σ(k) to Γ lσ(k) (the interior and exterior traces, with respect to the element T l ) . Similarly, z - k and z + k stand for the interior and exterior values of z h on Γ lσ(k) . For each interface Γ lσ(k) , k = 1, .., 3, we set:

z * k = max(z - k , z + k ), žk = z * k -max(0, z * k -η - k ) (11) and ȟ
- k = max(0, η - k -z * k ), ȟ+ k = max(0, η + k -z * k ), (12)
η-

k = ȟ- k + žk , η+ k = ȟ+ k + žk , (13)
leading to the new interior and exterior values:

w- k = t (η - k , ȟ- k η - k -z - k q - k ), w+ k = t (η + k , ȟ+ k η + k -z + k q + k ). ( 14
)
Now we set

H lσ(k) = H( w- k , w+ k , žk , žk , nlσ(k) ) + Ȟlσ(k) , (15)
as the numerical flux function through the interface between T l and T σ(k) , where:

1. the numerical flux function H(u -, u + , z -, z + , n) is classically computed from some 1d scheme (for example HLL or HLLC schemes [START_REF] Harten | On upstream differencing and godunov-type schemes for hyperbolic conservation laws[END_REF][START_REF] Toro | Restoration of the contact surface in the HLL riemann solver[END_REF] or relaxation schemes [START_REF] Berthon | A positive preserving high order VFRoe scheme for shallow water equations: a class of relaxation schemes[END_REF][START_REF] Jin | The relaxation scheme for systems of conservation laws in arbitrary space dimension[END_REF]), assumed to be conservative, Lipschitz continuous, monotone increasing with respect to u -, monotone decreasing with respect to u + and consistent with the exact flux in the following sense:

H(u, u, z, z, n) = H(u, z).n, ∀ (u, z). (16)
2. Ȟlσ(k) is a correction term needed to ensure flux balancing at motionless steady states, defined as follows:

Ȟlσ(k) =   0 0 g η- k (ž k -z - k ) 0 0 g η- k (ž k -z - k )   • nlσ(k) . (17)
Note that the modified interface fluxes [START_REF] Burbeau | A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods[END_REF] induce perturbations of order d + 1 when compared to the traditional interface fluxes. In the following, we choose to use the global Lax-Friedrichs flux :

H(w -,w + , z -, z + , n) = 1 2 ( H(w -, z -) • n + H(w + , z + ) • n -a(w + -w -)), (18)
with a = max 1≤l≤Ne λ l and

λ l = max ∂T l q h|T l η h|T l -z h|T l • nl + g(η h|T l -z h|T l ) . ( 19
)
Remark 3.3. For Γ ⊂ ∂Ω, we give sense to the definitions above by defining w + according to the boundary conditions.

Limiting procedure

Denoting w h ← ΛΠ h w h the limitation operator, acting on the approximated vector solution w h , the main steps of the limiting process are the following :

♯ 1 detect the shocks using the criterion proposed in [START_REF] Krivodonova | Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws[END_REF], and based on the strong superconvergence phenomena exhibited at element's outflow boundaries. More precisely, for a given element T l , let us denote ∂T l in the inflow part of ∂T l and identify K l in the set of superscript σ of the neighboring elements such that Γ lσ is a part of ∂T l in . Then, we use the following criterium :

I l = σ∈K l in Γ lσ (h l -h σ ) ds d (d+1)/2 l σ∈K l in ℓ lσ h l ∞ . (20)
If I l ≥ 1 then we apply a slope limiter on each scalar component w h|T l of w h|T l . This limiting strategy is described in the next step.

♯ 2 restrict the polynomial order to d = 1 in the element T l such that I l ≥ 1 and limit the slope of the resulting P 1 approximation following the simple and efficient procedure introduced in [START_REF] Burbeau | A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods[END_REF]. Furthermore, this limiting procedure is free of problem dependence and there is no parameters needed to be calibrated unlike in the widely used TVB limiter [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V[END_REF]. Let us quickly recall this procedure with some simplified notations. Considering the four elements patch neighborhood described in Fig. 1, let us denote c and (c i ) i=1..3 the respective barycenters of elements T and its three neighbors (T i ) i=1..3 . We also denote by (m i ) i=1..3 , the midpoints of the shared edges between T and T i respectively. Following [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V[END_REF], for each neighboring element T i , we can find T j and two positive parameters α i and β i such that :

cm i = α i cc i + β i cc j .
Then, for any linear function w h we can write :

w h (m i ) -w h (c) = α i (w h (c i ) -w h (c)) + β i (w h (c j ) -w h (c)).
Considering that the value of the function w h at the element barycenter is nothing but the cell average w h|T , and coming back to the general case of piecewise linear functions w h , the limiting procedure consists in replacing the values (w h|T (m i ) -w h|T ) i=1,..,3 , by limited values, denoted (∆ i ) i=1,..,3 , and defined as follows : ∆ i = 0, we replace (∆ i ) i=1,..,3 with ( ∆i ) i=1,..,3 , defined as follows :

∆ i = maxmod ∆ in i , ∆ out i , (21)
∆ in i =minmod w h|T (m i ) -w h|T , ν α i (w h|Ti -w h|T ) + β i (w h|Tj -w h|T ) , ∆ out i = minmod w h|T (m i ) -w h|T , w h|Ti (m i ) -w h|T ,
∆i = min(1, s + s - )max(0, ∆ i ) -min(1, s - s + )max(0, -∆ i )
with

s + = 3 i=1 max(0, ∆ i ), s -= 3 i=1 max(0, -∆ i ). T T 1 T 2 T 3 m 3 m 2 m 1 c 3 c 2 c 1 c

Additional limiting for robustness

We briefly recall the main ideas of the strategy introduced in [START_REF] Xing | Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations[END_REF][START_REF] Xing | Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water equations on unstructured triangular meshes[END_REF][START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high order discontinuous Galerkin schemes for conservation laws on triangular meshes[END_REF] to enforce a strict maximum principle that ensures the robustness property. The ideas are developed for an explicit first order Euler scheme in time for the sake of clarity. As shown in §3.8, these ideas can be extended to the PBSW equations. We first need to compute the new quadrature points for every element T l . Considering the scheme [START_REF] Biswas | Parallel, adaptive finite element methods for conservation laws[END_REF] for the update of η h , with polynomial expansions of order d, we assume that the line integrals are computed using a α-point Gauss quadrature. The Zhang and Shu quadrature rules are obtained by a transformation of the tensor product of a β-point Gauss-Lobatto quadrature, where β is the smallest integer such that 2β -3 ≥ d, and the α-point Gauss quadrature. This special quadrature includes all α-point Gauss quadrature nodes for each edge Γ lσ(k) , k = 1, 2, 3, involves positive weights and its degree is chosen such that it is exact for the integration of η n h|T l over T l . In the following, let us denote S d T l the set of points of this new quadrature rule and ωβ 1 the weight associated with the first node of the β-point Gauss-Lobatto quadrature. We do not give further details and the reader is referred to [START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high order discontinuous Galerkin schemes for conservation laws on triangular meshes[END_REF] for explicit nodes coordinates and weights. In §4, we focus on P 1 and P 2 approximations of the SWE weak solutions, and we show on Figure [START_REF] Ambati | Space-time discontinuous Galerkin finite element method for shallow water flows[END_REF] the quadrature nodes used for these two orders of approximation on a reference element. The robustness property is then enforced in 2 steps : ♯ 1 First, for each element T l , we compute the polynomial expansions for the water height h n h|T l (x) from η n h|T l (x) and z h|T l (x) (this is achieved straightforwardly, as we are working with the nodal basis). We then extract the quantities m T l = min

x∈S d T l h n h|T l (x).
♯ 2 Next, we need to ensure that h n h|T l (x) ≥ 0, ∀x ∈ S d T l , which is a sufficient condition to ensure a robustness property for the DG scheme [START_REF] Biswas | Parallel, adaptive finite element methods for conservation laws[END_REF], under a suitable CFL-like condition, as shown in §3.8. In practice, we ensure that h n h|T l (x) ≥ ǫ, ∀x ∈ S d T l , where ǫ is a small positive real threshold value. This can be enforced using the accuracy preserving limiter introduced in [START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high order discontinuous Galerkin schemes for conservation laws on triangular meshes[END_REF]. Denoting h n T l the mean value of the water height at time t n on the element T l , and assuming that h n T l ≥ 0, we replace h n h|T l by a conservative linear scaling around the cell average :

ĥn h|T l = θ T l (h n h|T l -h n T l ) + h n T l , (22)
where Remark 3.4. Following steps ♯1, 2 at each time step (or eventually substep if a higher order time-discretization is used), and assuming that h 0 T l ≥ ǫ, we ensure that h n T l remains positive at every time step n. Additionally, the scaling [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V[END_REF] ensures that the values of ĥh|T l remain positive at the (d+1)-point Gauss quadrature nodes used to compute the line integrals of (10), the eigenvalues [START_REF] Castro | Well-balanced numerical schemes based on a generalized hydrostatic reconstruction technique[END_REF] used in the Lax-Friedrichs flux [START_REF] Castro | Well-balanced high order extensions of godunov's method for semilinear balance laws[END_REF] and the CFL condition [START_REF] Giraldo | Nodal high-order discontinuous Galerkin methods for the spherical shallow water equations[END_REF]. In practice, the parameter ǫ is a threshold introduced to identify what is numerically called a dry cell and arbitrarily set the velocity to zero.

θ T l = min ǫ -h n T l m T l -h n T l , 1 .
Remark 3.5. The Zhang and Shu quadrature rules are not used to compute the surface integrals, as we prefer to use more efficient cubature rules on triangles. Practically, its only purpose is to allow the computation of the quantities m T l and therefore enables to enforce the robustness sufficient condition introduced in [START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high order discontinuous Galerkin schemes for conservation laws on triangular meshes[END_REF].

Time discretization

The time stepping is carried out using the explicit second and third-order SSP-RK scheme [START_REF] Gottlieb | Strong stability preserving high order time discretization methods[END_REF]. Writing the semidiscrete equations as d dt w h + A h (w h ) = 0, advancing from time level n to n + 1 are computed as follow for the second-order scheme :

w n,1 h = w n h -∆t n Ãh (w n h ) , w n+1 h = 1 2 (w n h + w n,1 h ) -1 2 ∆t n Ãh (w n,1 h ) , (23)
and for the third-order scheme :

       w n,1 h = w n h -∆t n Ãh (w n h ) , w n,2 h = 1 4 (3w n h + w n,1 h ) -1 4 ∆t n Ãh (w n,1 h ) , w n+1 h = 1 3 (w n h + 2w n,2 h ) -2 3 ∆t n Ãh (w n,2 h ) . (24)
with the formal notation Ãh = A h • ΛΠ h , and ∆t n being the time discretization step, obtained from the following CFL condition [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF] :

max 1≤l≤Nt λ l p l |T l | ∆t n ≤ 1 2d + 1 . ( 25
)

Main properties

We focus here on two suitable properties verified by the scheme [START_REF] Biswas | Parallel, adaptive finite element methods for conservation laws[END_REF]. We first show that the well-balancing property is inherited from the particular interface fluxes discretization [START_REF] Burbeau | A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods[END_REF] for polynomial expansions of arbitrary orders. Next, we show that the method of [START_REF] Xing | Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations[END_REF][START_REF] Xing | Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water equations on unstructured triangular meshes[END_REF] can be adapted to the PBSW equations ( 5) to ensure the robustness.

Proposition 1. The scheme [START_REF] Biswas | Parallel, adaptive finite element methods for conservation laws[END_REF], with the interface fluxes discretization [START_REF] Burbeau | A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods[END_REF], preserves the motionless steady states. In other terms, (

Proof. We adapt the ideas of [START_REF] Xing | Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations[END_REF] to the reconstruction introduced in §3.4. We consider T l an element of T h , together with his neighbors (T σ(k) ) k=1,..,3 . Denoting that (3) is equivalent to the following local formulation of the well-balancing property, for all n ∈ N :

η n h|T l ≡ η n h|T σ(k) ≡ η e q n h|T l ≡ q n h|T σ(k) ≡ 0 , k = 1, 2, 3 ⇒ η n+1 h|T l ≡ η e q n+1 h|T l ≡ 0 , (26)
let's assume that the left brace of (26) holds. Let's also assume that a first order Euler time discretization is used to compute w n+1 h from w n h . Considering [START_REF] Biswas | Parallel, adaptive finite element methods for conservation laws[END_REF], we have to show that the residues

R j = - T l H(w n h , z h ) • ∇θ l j dx + 3 k=1 Γ lσ(k) H lσ(k) θ l j ds - T l S(w n h , z h )θ l j dx, (27)
vanishes for 1 ≤ j ≤ N d , when a motionless steady state (3) is reached. We assume that the integrals in [START_REF] Delestre | A numerical scheme for a viscous shallow water model with friction[END_REF] are computed exactly at steady states, and we notice that we have for w n h at steady state:

∇ • H(w n h , z h ) = S(w n h , z h ).
Looking at [START_REF] Burbeau | A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods[END_REF], and highlighting that, for each interface Γ lσ(k) , we have ηk = η+ k = η e , it is easy to check that H lσ(k) = H(w - k , z - k ) • nl , with the notations introduced in §3.4. Consequently, we have

- T l H(w n h , z h ) • ∇θ l j dx + 3 k=1 Γ lσ(k) H lσ(k) θ l j ds = T l ∇ • H(w n h , z h )θ l j dx, (28)
and R j = 0. Now, let us consider the robustness property. At first order in time, the scheme satisfied by the cell averaged values of the free surface in the dG approximation [START_REF] Biswas | Parallel, adaptive finite element methods for conservation laws[END_REF] is:

η n+1 T l = η n T l - ∆t n |T l | 3 k=1 Γ lσ(k) H η lσ(k) ds, (29)
with H η lσ(k) the first component of the numerical flux, defined as

H η lσ(k) = H η ( w- k , w+ k , žk , žk , nlσ(k) ). ( 30
)
To apply the strategy of [START_REF] Xing | Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations[END_REF][START_REF] Xing | Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water equations on unstructured triangular meshes[END_REF], we first have to show that the first order scheme associated with (29) preserves the positivity of the water height. Let us denote w n l and (w n σ(k) ) 1≤k≤3 the piecewise constant values of the solution at the discrete time t n respectively on T l and its 3 neighbors. Similarly, z l and (z σ(k) ) 1≤k≤3 stand for the corresponding piecewise constant values of the topography. Now, for a given interface Γ lσ(k) , k = 1, .., 3, we set:

w - k = w n l , w + k = w n σ(k) , (31)
z - k = z l , z + k = z σ(k) , (32)
as interior and exterior values, with respect to the element T l . Noting that, for a given element T l we have h n l = η n l -z l , we have the following Lemma:

Lemma 3.6. We consider a first order scheme for the free-surface:

η n+1 l = η n l - ∆t n |T l | 3 k=1 Γ lσ(k) H η lσ(k) ds, (33)
with

H η lσ(k) = H η ( w- k , w+ k , žk , žk , nlσ(k) ), (34)
defined from the Lax-Friedrichs flux (18) and wk , w+ k , žk defined according to ( 11)-( 14) and the piecewise constant values ( 31)- [START_REF] Eskilsson | A triangular spectral/hp discontinuous Galerkin method for modelling 2d shallow water equations[END_REF]. If h n l ≥ 0 and h n σ(k) ≥ 0, 1 ≤ k ≤ 3 then h n+1 l ≥ 0 under the condition:

λ l p l |T l | ∆t n ≤ 1. ( 35
)
Proof. We first remark that according to (13), we have η+

k -η- k = ȟ+ k -ȟ- k , for each interface Γ lσ(k) .
Using the first component of the Lax-Friedrichs flux (18)

H η (w -,w + , z -, z + , n) = 1 2 ((q - x + q + x )n x + (q - y + q + y )n y -a(η + -η -)),
and subtracting z l at both side of (33), we can straightforwardly reproduce the proof of [START_REF] Xing | Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water equations on unstructured triangular meshes[END_REF], writing the right hand side of (33) as a linear combination of positive values.

We can now state the robustness result:

Proposition 2. We consider the (d+1)-th order scheme [START_REF] Biswas | Parallel, adaptive finite element methods for conservation laws[END_REF], together with the interface fluxes (15) and a first order Euler time discretization. We assume that h n T l ≥ 0 and that the quantities h n h|T l (x)

x∈S d T l
are modified according to [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V[END_REF], for all T l . Then we have h n+1 T l ≥ 0 under the condition

λ l p l |T l | ∆t n ≤ 2 3 ωβ 1 . ( 36
)
Proof. Noting that the scheme ( 29) can be written as follows, subtracting z T l at both sides:

h n+1 l = h n l - ∆t n |T l | 3 k=1 Γ lσ(k) H η lσ(k) ds, (37)
the proof of [START_REF] Xing | Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water equations on unstructured triangular meshes[END_REF] can be straightforwardly reproduced, using Lemma 3.6.

Remark 3.7. Even if the limiter [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V[END_REF] does not modify the (d+1)-th order accuracy of the polynomial expansion w h , it is necessary to use a (d+1)-th order time SSP-RK discretization to achieve an overall (d+1)-th order accuracy in space and time. Therefore, steps ♯1, 2 should be applied at every intermediate step of the RK algorithms ( 23) and [START_REF] Cox | Experimental and numerical modelling of surf zone hydrodynamics[END_REF]. To fulfill the CFL restriction ( 36) rigorously, we have to provide an accurate estimation of max 1≤l≤Nt λ l for all the stages of the SSP-RK method. As suggested in [START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high order discontinuous Galerkin schemes for conservation laws on triangular meshes[END_REF], this can be achieved with applying a more stringent CFL condition only when a preliminary computation to the next intermediate time step produces negative water height.

Numerical Validations

In this section, we highlight the good behavior of the combination suggested above. The Lagrangian interpolant P 1 and P 2 basis functions are used, together with the global Lax-Friedrichs flux [START_REF] Castro | Well-balanced high order extensions of godunov's method for semilinear balance laws[END_REF]. We will refer to the limiter of Cockburn and Shu [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V[END_REF] as CS limiter and to the limiter of Burbeau et al [START_REF] Burbeau | A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods[END_REF] described in §3.5 as BSB-limiter.

In the following computations, the parameter ν is set to 2 and the threshold value ǫ for the robustness limiter is set to 10 -9 . In the two dimensional case on unstructured meshes, the CFL condition [START_REF] Giraldo | Nodal high-order discontinuous Galerkin methods for the spherical shallow water equations[END_REF] is always more stringent than the usual CFL [START_REF] Dawson | A discontinuous Galerkin method for three-dimensional shallow water equations[END_REF], and therefore only condition (36) is applied. For the P 1 expansion, we use a 3-point Gauss quadrature rule (i.e. α = 3 with the notations of §3.6) to compute the line integrals of the first equation of [START_REF] Biswas | Parallel, adaptive finite element methods for conservation laws[END_REF]. For the P 2 expansion, we use a 4-point Gauss quadrature rule. In both case a 3-point Gauss-Lobatto rule is enough to build the required Zhang and Shu quadratures (see Fig. 2), leading to the value ωd 1 = 1 6 .

Collapsing of a gaussian profile

The following test is commonly used in order to evaluate the diffusive properties of first and second order schemes, see [START_REF] Audusse | A well-balanced positivity preserving "second-order" scheme for shallow water flows on unstructured meshes[END_REF][START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF] for instance. We study the evolution of the flow resulting from a water drop in a center of a square basin with dimensions [0, 20] × [0, 20]. The computational domain is meshed with a regular triangulation (see Fig. 3), with a discretization step fixed to ∆x = ∆y = 0.25. Solid wall conditions are enforced at the boundaries. The bottom is assumed to be flat, and the initial conditions are given by :

h(x, y) = 2.4 1 + e 1 4 ((x-10) 2 +(y-10) 2 ) , q = 0 . ( 38
)
Free surface profiles provided by the P 1 approximation are shown on Fig. 4 at several times during the simulation between t = 0 s and t = 2 s, together with a P 2 highly resolved solution, standing for a reference. Numerical results obtained with the second order vertex-centered unstructured Finite Volume scheme introduced in [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF] are also displayed. Note that DG and FV computations involve the same number of freedom degrees. There is a slight improvement provided by the current DG approach. ∆x ∆y Figure 3: Collapsing of a gaussian profile -Example of regular mesh.

Dam break problems

This second test case aims at validating the shock capturing abilities of the scheme [START_REF] Biswas | Parallel, adaptive finite element methods for conservation laws[END_REF], together with the efficiency of the limiting strategies in the vicinity of a dry area.

Dam break on a dry bed

Here, the computational domain is a [-10, 10] × [0, 4] rectangular basin, regularly meshed with 5200 nodes (see Fig. 3). The initial condition is given by h = 1 and q = 0 on the left part of the domain and h = 0 on the right part. The evolution of the wet/dry interface is accurately computed. We can observe on Fig. 5 some water depth profiles during the evolution: they exhibit an excellent agreement with the analytical solution. For the records, the results are compared with those provided by the CS limiter. The discrepancies are almost indistinguishable at this level of zoom. We show zoom on Fig. 6, focusing on the top part of the free surface profile at t = 0.6 s. Although the discontinuity sensor is no longer active at this time, the limiting process engaged during the first steps to handle the initial discontinuity impacts the evolution of the flow. We observe a slight advantage of the BSB limiter.

Two-dimensional dam-break

We study now the ability of the scheme in dealing with discontinuities in a real 2D context. In this simulation, the computational domain consists in a square basin with dimensions 200 × 200, and the bottom is flat. The dam is located along the y-direction centreline, and the initial flow is supplied by a 75 m breach on the dam wall, centred at y = 125. We consider an upstream water height of 10 m, while the initial downstream water level is fixed to h d = 5. We run the simulation on an unstructured triangulation involving 13000 nodes, until t = 8 s. As shown in Fig. 7 (top), the propagation of the flood wave is properly described, and the stiffness of the water front profile is correctly handled by the limitation procedure. The numerical predictions are in good agreement with other results reported in the literature, notably those provided in [START_REF] Wang | A comparative study of finite volume methods on unstructured meshes for simulation of 2d shallow water wave problems[END_REF]. Lastly, we run again the simulation setting h d = 0. We can observe on Fig 7 (bottom) a 3D view of the free surface at t = 6.5 s and highlight that the wet/dry interface is handled in a very robust way.

Small perturbation of a steady state

We assess here the well-balanced property, together with the ability of the scheme to accurately describe the evolution of a small disturbance of a motionless steady state, with a varying topography. This test, introduced in [START_REF] Leveque | Balancing source terms and flux gradients in high-resolution godunov methods: the quasi-steady wave-propagation algorithm[END_REF], is widely used to exhibit the efficiency of 2D well-balanced methods in the vicinity of steady states, see for instance [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF][START_REF] Xing | A new approach of high order well-balanced finite volume WENO schemes and discontinuous Galerkin methods for a class of hyperbolic systems with source terms[END_REF]. This simulation involves a rectangular channel of 2 m × 1 m, and the topography variations are given by : z(x, y) = 0.8e -5(x-0.9) 2 -50(y-0.5) 2 .

(39)

Preservation of a motionless steady state

L 1 L ∞ eη 0.264E-14 0.133E-14 eq x 0.101E-13 0.710E-14 eq y 0.813E-14 0.712E-14

Table 1: Preservation of a motionless steady state -Numerical errors at t = 10 s.

First, we consider a lake at rest configuration with a 1 m total free surface. The free surface is constant and the discharge set to zero on the whole computational domain. After a simulation time of 10 s, the motionless steady state is preserved up to the machine accuracy. We show the resulting numerical errors in Tab. 1 computed with the P 1 scheme at t = 10 s for η and q in L 1 and L ∞ norms. FV 2nd order P1 Reference Figure 8: Perturbation of a motionless steady state -Vertical middle section at t = 0.12, 0.24, 0.36, 0.48 s. Comparison with a second order Finite Volume scheme and the reference.

Perturbation of a motionless steady state

Here, we enforce an initial water height of 1 + δ in the rectangular band 0.05 < x < 0.15, and follow the motion of the flow until t = 0.5 s. In this test, δ is set to 0.01, and we use a 3321 nodes mesh. We can observe on Fig. 10 a 3D view of the free surface as the flow reaches the top of the bump in the center of the basin. We obtain comparable results to those obtained in the literature with high order schemes, [START_REF] Bryson | Well-balanced positivity preserving central-upwind scheme on triangular grids for the Saint-Venant system[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF][START_REF] Xing | A new approach of high order well-balanced finite volume WENO schemes and discontinuous Galerkin methods for a class of hyperbolic systems with source terms[END_REF], assessing an accurate computation of the propagation. Focusing on the way the initial discontinuity is handled, Fig. 9 shows the influences of both CS and BSB limiters. Like in 4.2.1, there is a small discrepancy arising from the choice of the limiting process used to handle the initial discontinuity. A slight improvement is observable with the use of BSB approach. In a second time, we can observe on Fig. 8 some middle sections of the free surface profiles obtained with the P 1 scheme as the flow is evolving over the bump. These results are compared with a reference solution, computed with the P 2 scheme on a 20000 nodes mesh, and those issued from the second order FV scheme of [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF] (dashed lines). The number of freedom degrees for the Finite Volume test is set to 10000, so that a relevant comparison can be performed between the FV and DG approaches; the efficiency of the P 1 DG approximation is similar to the MUSCL reconstruction method of the FV scheme for this test. The reader is referred to test cases 4.4 and 4.5 for a more detailed comparison of the two methods.

Subcritical flow over a bump

In this test, we study the ability in converging toward a stationary state and we perform some convergence studies. We consider a flow over a non-uniform bathymetry of the form:

z(x, y) = 0.2 -0.05(x -10) if 8 ≤ x ≤ 12, 0 otherwise. (40)
Three main classes of solutions are available, namely the sub-critical case and the trans-critical cases (with or without shock), described in [START_REF] Goutal | Dam-break wave simulation[END_REF]. We consider here the sub-critical case, with a rectangular basin of 20 × 5. We enforce an inflow discharge of q x = 4, 42, q y = 0, and use an outflow boundary condition of η = h = 2.

1.9

1.92 The steady flow is approximatively reached at t = 150 s. Comparing numerical results with the analytical solution, we observe on Fig. 11 an excellent agreement for both free surface and normal discharge. These results are obtained with the P 1 scheme on a mesh involving 800 elements. Considering the well balancing property, this test also highlight the satisfying performance of the numerical scheme for the preservation of steady states involving flows in motion. We perform an accuracy analysis for the P 1 and P 2 schemes, considering regular meshes, with a discretization step ∆x = 20/n (see Fig. 3) and an increasing level of refinement: n = 20, 40, 80, 160. At the end of the simulation, numerical results are compared to the analytical steady state allowing to compute numerical errors for the free surface and the normal discharge. Fig. 12 shows the L 1 -error with respect to the value of ∆x in logarithmic scale. Optimal convergence is reached for the both schemes. To further investigate the comparison between the P 1 DG scheme and the MUSCL FV scheme of [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF], additional investigations are performed involving the total number of freedom degrees (DoF), numerical errors and corresponding CPU time. We stress out that the subject of CPU time evaluation has to be broached carefully, as it highly depends on the implementation choices. The CFL is fixed to 0.5 for both schemes and at t = 150 s, we compute the numerical L 1 -error together with the corresponding CPU time, for an increasing number of DoF. Results are reported in Tab. 2. Then, we plot on Fig. 13 the numerical error with respect to the CPU time and the CPU time with respect to the number of DoF. Results are plotted in a log-log scale, and we focus on the water height, as similar observations stands for the discharge. We point out the slight advantage of the FV approach in terms of CPU time for a given number of DoF. However, the numerical results provided by the P 1 DG scheme are clearly more accurate for comparable CPU time (or comparable number of DoF).

Carrier and Greenspan transient solution

Here, we study the evolution of a flow in a rectangular channel [-20, 6]×[0, 10] with a sloping bottom, parametrized as follows:

z(x, y) = s x with s = 1 50 , (41)
and we focus on the wet/dry transition evolution, exhibiting the ability of the scheme in describing complex flooding and drying phenomena. Computations are run with the P 1 scheme and numerical results are compared with the analytical solution available in [START_REF] Carrier | Water waves of finite amplitude on a sloping beach[END_REF]. In order to perform some convergence analysis, we use a regular triangulation with a space step fixed initially to ∆x = 0.5. Initial conditions are provided by the exact solution w ex at t = 0 s shown on Fig. 15 (top). The flooding phenomena is governed by the offshore boundary condition, given by w ex (t, -20), and we observe an asymptotic convergence towards a steady state with a planar free surface, observed at approximately t = 40 s, see Fig. 15 (bottom). We show on Fig. 14 the evolution of the analytical and predicted free surface profiles along the middle section during the simulation; we can observe an excellent description of the moving interface location.

We finally perform some convergence rate analysis. Several computations are run on a set of regular meshes with space step ∆x = 2/2 n , n = 0, • • • , 4, and the L 1 error is computed at t = 0.1 s. As shown in Fig. 16, we reach rates of 2 and 2.8 respectively in the P 1 and P 2 case, respectively for free surface and normal discharge.

As for the previous test, we further investigate the efficiency of the P 1 DG scheme when compared to the MUSCL FV scheme of [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF]. The same analysis are performed, at t = 0.5 s, and results are reported in Tab. 3 and Fig. 17. It is interesting to notice that even with the occurrence of the wet/dry fronts, the DG scheme provides the expected rates of convergence, whereas the MUSCL FV scheme only provides a convergence rate of 1.6, see [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF]. This clearly highlights the accuracy preserving feature of the strategy of [START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high order discontinuous Galerkin schemes for conservation laws on triangular meshes[END_REF]. We can observe on Tab. 3 and Fig. 17 that the FV approach is still more efficient when considering the needed CPU time for a fixed number of DoF. However, unlike the previous test case, numerical errors are smaller for the MUSCL FV scheme, for comparable CPU time (or comparable number of DoF).

Tsunami wave on a conical island

We now consider a true two-dimensional case and study the propagation of a solitary wave over a conical island. Laboratory experiments were performed [START_REF] Liu | Runup of solitary waves on a circular island[END_REF]. Several authors have used this test as a benchmark to study the run up phenomena [START_REF] Bollermann | Finite volume evolution galerkin methods for the shallow water equations with dry beds[END_REF][START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF][START_REF] Hubbard | A 2d numerical model of wave run-up and overtopping[END_REF][START_REF] Nikolos | An unstructured node-centered finite volume scheme for shallow water flows with wet/dry fronts over complex topography[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF]. The basin dimensions are [0, 25] × [0, 30], and we use an unstructured triangulation involving 23000 nodes. Denoting r the distance from the center of the domain, we consider an ideal island defined as follows:

z(r) = max 0.625, 0.9 -1/ 4 r if r < 3.6, 0 elsewhere .

We enforce an initial motionless configuration with a mean water depth h 0 = 0.32. The tsunami wave generation is supplied by a time dependent offshore boundary condition, as described in [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF]:

     h(t) = h 0 + αh 0 sech 2 √ gh0 L χ t , χ = 3α 4β (1 + α) , β = (h 0 /L) 2 , L = 15 , α = 0.3 . ( 43
)
We follow the evolution of the solitary wave after the impact at the front side of the island, and observe that all the features of the propagation are properly reproduced, according to the data from the experiment. In Fig. 19, we display time series of the free surface at several locations along the basin and the comparison between our numerical results and the data obtained from wave gauges show reasonable agreement, qualitatively similar to the numerical results found in the references above. We can clearly observe on Fig. 18 the run up issuing from the reflection at the front of the island, the separation in two secondary waves surrounding the island, and finally colliding at the rear side.

Figure 18: Tsunami wave over a conical island -Lateral view of the free surface at times t=5, 6, and 7s. Rear view at t=8, 9 and 10s.

-

Cox Experience

In this test we compare our numerical results with experimental data issuing from the study carried out by Cox [START_REF] Cox | Experimental and numerical modelling of surf zone hydrodynamics[END_REF]. Our goal is to highlight the efficiency of the discontinuity sensor [START_REF] Castro | High order exactly well-balanced numerical methods for shallow water systems[END_REF] and the subsequent limitation process described in §3.5. This experiment was carried out in a wave flume of 33 m long. Waves were generated on a horizontal bottom at a depth h 0 = 0.40, shoaled and broke on a 1/35 sloping beach. The wave height at the wave-maker is a = 0.115 and the period is T = 2.2 s. Measurements of the surface elevation and velocity were taken at four locations along the basin.

The following numerical simulation involves a domain with dimensions [0, 8] × [0, 1], regularly meshed with a space step ∆x = 0.01. The first set of data is used as a time evolving inflow boundary condition at the left side of the domain, to simulate the wave-maker. In order to simulate at best physical conditions, we also introduce a friction term in the equations ( 5), following the Darcy law, with a coefficient n calibrated at 0.01 . After 200 s of simulation approximatively, we obtain some periodic free surface profiles. We can observe a very good description of the wave distortion in the vicinity of the shoreline, with the expected sharped wave profiles. In order to demonstrate the good numerical description of the process, free surface profiles at gauges x = 1.2, x = 2.4 and x = 3.6 are available on Fig. 20. We observe a good agreement with the data provided by Cox. In particular, it is worth emphasizing that the shock detector (20) perfectly allows to capture the discontinuities generated during the breaking of the waves, activating the local limitation only in the sharped wave front areas (see Fig. 21).

Solitary wave on a sloping beach

We finally consider a two-dimensional problem of moving shoreline, using the topography suggested by Zelt [START_REF] Zelt | Tsunamis: the response of harbors with sloping boundaries to long wave exitation[END_REF]. We consider the propagation of a solitary wave breaking on a bay, in a basin of dimensions [-10, 15] × [-10, 10]. The topography is defined as follows:

z(x, y) =    d 0
x -x p l cos(πy/l)/π + x p if x ≥ x p , 0 elsewhere .

The constants are defined as follows : l = 10, d 0 = 1.273, x p = 3 l/π. Considering a cross-shore section, we observe that the slope of the topography slowly decreases from the lateral boundaries to the x-direction centreline (see Fig. 22). The generation of the wave is performed as in test case 4.6 at the left offshore boundary, see formula [START_REF] Harten | On upstream differencing and godunov-type schemes for hyperbolic conservation laws[END_REF], with L = 10, α = 0.02, and a mean water depth h 0 = d 0 . It results an initial steady state configuration, shown on Fig. 23 (top) and perturbed by an incoming wave with an amplitude of H = αh 0 . The basin is meshed with an unstructured triangulation of 15000 elements. We mainly focus on the description of the run up and run down phenomena occurring when the wave reach the shore. In Fig. 23, free surface snapshots are exhibited at several times during the flooding process. According to the bottom profile, part of the energy is gradually transferred in the center of the basin; this phenomena is supplied by the flooding and drying process, which evolution is available in Fig. 25. Maximum run-up and minimum run-down computed during the propagation are available on Fig. 24, where a comparison the out FV-MUSCL approach is also provided; these results highlight the robustness of the model and are in agreement with those observed in other works, see [START_REF] Hubbard | A 2d numerical model of wave run-up and overtopping[END_REF][START_REF] Marche | Theoretical and numerical study of shallow water models[END_REF][START_REF] Özkan-Haller | A Fourier-Chebyshev collocation method for the shallow water equations including shoreline runup[END_REF] for instance.

Conclusion

In this work, we performed a review of the existing DG methods for the SWE and we focus on the issues of limiting, well-balancing and robustness. We then introduce a combination of ingredients leading to a robust and well-balanced discretization of the two dimensional SWE with topography source term, on unstructured triangular grids. This well-balancing property is reached adapting the well-balanced FVM of [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF] to the dG framework and the robustness is achieved adapting the ideas of [START_REF] Xing | Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water equations on unstructured triangular meshes[END_REF] to the PBSW equations. An efficient shock detecting and local limitation process are implemented. The resulting scheme is shown to be accurate and stable and very efficient for the simulation of a large variety of flows involving flooding and drying phenomena.

s

  min(|a|, |b|, |c|) if s = sgn(a) = sgn(b) = sgn(c) , 0 otherwise. maxmod(a, b, c) =    s max(|a|, |b|, |c|) if s = sgn(a) = sgn(b) = sgn(c) , 0 otherwise. and ν > 1. Note that this modification may not be mass conserving. Consequently, if 3 i=1

Figure 1 :

 1 Figure 1: Limitation step : geometric configuration.

Figure 2 :

 2 Figure 2: Nodes locations for the Zhang and Shu quadrature -P 1 and P 2 cases.

Figure 4 :

 4 Figure 4: Collapsing of a gaussian profile -Free surface profiles for the FV and P 1 approaches at t = 0.48, 0.96, 1.44 s along the middle section. Comparison with the reference solution.

Figure 5 :

 5 Figure 5: Dam break on a dry bed -Water depth profiles at t=0.1, 0.6, 1.2 and 1.8s.

Figure 6 :

 6 Figure 6: Dam break on a dry bed -Zoom on the top at t=0.6s.

Figure 7 :

 7 Figure 7: Two-dimensional dam break -Free surface profile at t=7.5s for h d = 5m (top), and at t=6.5s for h d = 0m (bottom).

Figure 9 :

 9 Figure 9: Perturbation of a motionless steady state -Zoom on the discontinuity at t = 0.015 s.

Figure 10 :

 10 Figure 10: Perturbation of a motionless steady state -Snapshot of the free surface at t = 0.3 s.

Figure 11 :

 11 Figure 11: Subcritical flow over a bump -Steady state : values of the free surface and normal discharge along the middle section. Analytic vs numeric.

Figure 13 :

 13 Figure 13: Subcritical flow over a bump -Comparison between P 1 DG and MUSCL FV schemes : L 1 -error with respect to CPU time (in s.) and CPU time (in s.) with respect to the number of DoF.

Figure 14 :

 14 Figure 14: Carrier and Greenspan transient solution -Evolution of η along the x direction centreline.

Figure 15 :

 15 Figure 15: Carrier and Greenspan transient solution -Snapshots of the initial condition (top) and free surface profiles at t = 15 s and t=40 s.

Figure 17 :

 17 Figure 17: Carrier and Greenspan transient solution -Comparison between P 1 DG and MUSCL FV schemes : L 1 -error with respect to CPU time (in s.) and CPU time (in s.) with respect to the number of DoF.

Figure 20 :

 20 Figure 20: Cox Experience -Time history of the water level at x=1.2, 2.4 and 3.6 m.

Figure 21 :

 21 Figure 21: Cox Experience -Snapshot of the water depth profile during the breaking process. Areas in need of limiting are specified along the dotted line.

Figure 23 :

 23 Figure 23: Solitary wave on a sloping beach -Free surface profiles at t=0, 11, 14 and 17s.

Figure 24 :

 24 Figure 24: Solitary wave on a sloping beach -Maximum run up and minimum run down during the evolution.

Figure 22 :

 22 Figure 22: Solitary wave on a sloping beach -Topography and configuration.

1 Figure 25 :

 125 Figure 25: Solitary wave on a sloping beach -Time series of the run up along five cross sections.

  Figure 12: Subcritical flow over a bump -Convergence analysis : L 1 numerical error for the free surface and the discharge.

	P1 η P2 η P1 q x P2 q x
						3
				h			hu
	∆x		P 1	P 2			P 1	P 2
	1		0.55e-4	0.96e-5	0.53e-4	0.30e-4
	0.5		0.10e-4	0.10e-5	0.90e-5	0.33e-5
	0.25	0.25e-5	0.13e-6	0.18e-5	0.33e-6
	0.125	0.57e-6	0.28e-7	0.43e-6	0.60e-7
		P 1 DG				MUSCL-FV
	DoF	CP U (s)	h	DoF	CP U (s)	h
	600	11	5.5e-5	750		8	3.3e-3
	1200	36	1.0e-5	1470	33	8.2e-4
	2400	150	2.5e-6	2910	130	2.1e-4
	4800	671	5.7e-7	5790	630	5.3e-5

Table 2 :

 2 Subcritical flow over a bump -Comparison between P 1 DG and MUSCL FV schemes : numerical error and corresponding CPU time for an increasing number of DoF.

  Figure 16: Carrier and Greenspan transient solution -Convergence analysis : L 1 error for the free surface and the normal discharge.

		P1 η P2 η P1 q x P2 q x
			1
		2
							2.8
						1
				h			hu
		∆x	P 1	P 2		P 1	P 2
		2	4.0e-3	1.3e-3	1.3e-1	5.5e-2
		1	1.0e-3	2.4e-4	3.2e-2	4.9e-3
		0.5	3.3e-4	3.9e-5	1.1e-2	7.5e-4
		0.25	9.9e-5	5.4e-6	3.3e-3	1.0e-4
		0.125	2.2e-5	8.0e-7	7.5e-4	1.9e-5
		P 1 DG			MUSCL-FV
	DoF	CP U (s)	h	DoF	CP U (s)	h
	1200	0.2	4.8e-4	1470	0.2	1.6e-4
	2400	0.8	1.7e-4	2910	0.7	4.5e-5
	4800	3.7	3.9e-5	5790	2.7	2.0e-5
	9600	14.5	9.2e-6	11550	11.0	7.5e-6

Table 3 :

 3 Carrier and Greenspan transient solution -Comparison between P 1 DG and MUSCL FV schemes : numerical error and corresponding CPU time for an increasing number of DoF.

  Figure19: Tsunami wave on a conical island -Time history of the free surface at gauges 3,6,9,16 and 22. Numerical (solid line) vs experimental data (dashed line).

			gauge 3 (6.82 , 13.05)			gauge 6 (9.36 , 13.80)
		0.04				0.04
	η(m)	0 0.02			η(m)	0 0.02
		0.02				-0.02
		5	10	15	20	5	10	15	20
			t(s)				t(s)
			gauge 9 (10.36 , 13.80)			gauge 16 (12.96 , 11.22)
		0.06				0.06
		0.04				0.04
	η(m)	0 0.02			η(m)	0 0.02
		-0.02				-0.02
		-0.04				-0.04
		5	10	15	20	5	10	15	20
			t(s)				t(s)
					gauge 22 (15.56 , 13.80)
				0.06
				0.04
			η(m)	0 0.02
				-0.02
				-0.04
				5	10	15	20
					t(s)

Acknowledgments

The authors would also like to acknowledge additional financial and scientific support of the French INSU-CNRS (Institut National des Sciences de l'Univers-Centre National de la Recherche Scientifique) program LEFE-MANU (Méthodes Mathématiques et Numériques), project SOLi. F.M. also acknowledges partial support from the ANR BLANC Bond.