Ali Haidar
email: ali.haidar@umontpellier.fr

Fabien Marche
email: fabien.marche@umontpellier.fr

François Vilar

Francois Vilar
email: francois.vilar@umontpellier.fr

A posteriori Finite-Volume local subcell correction of high-order discontinuous Galerkin schemes for the nonlinear shallow-water equations

Keywords: Non-linear shallow water equations, discontinuous Galerkin scheme, local a posteriori subcell FV correction, well-balancing, positivity preservation, arbitrary high-order

come

A posteriori Finite-Volume local subcell correction of high-order discontinuous Galerkin schemes for the nonlinear shallow-water equations

Introduction

The nonlinear shallow water (NSW) 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] are one of the most widely used set of equations for simulating long wave hydrodynamics, under the assumption that the vertical acceleration of the fluid can be neglected. Given a smooth parametrization of the topography variations b : R → R, and denoting by H the water height, u the horizontal (depth-averaged) velocity and q = Hu the horizontal discharge (see Fig. 1), the NSW equations may be written as follows:

∂ t H + ∂ x q = 0, ∂ t q + ∂ x uq + 1 2 gH 2 = -gH∂ x b. (1a) (1b)
Considering their hydrostatic nature, in comparison to the dispersive nature of more sophisticated models such as the Boussinesq-type models, hyperbolic integral forms of the NSW equations generally provide an accurate representation of steep-fronted flows, such as dam breaks, flood waves or bores propagation in the surf zone. 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 and to forecast coastal inundations.

To allow a proper description of such phenomena, accurate and robust numerical methods have to be considered. Great efforts have been made since the sixties in order to produce accurate approximations of weak solutions of the NSW equations and a large variety of numerical methods have been developed, including Finite-Volumes (FV) (2; 43; 5; 10; 104; 39), Finite-Elements (FE) (77; 96; 83; 9), spectral methods (55; 71; 80) or residual distribution methods (88; 87; 6). Among these numerical strategies, the Godunov-type FV methods are particularly praised, thanks to their low computational cost and their shock-capturing ability, which allows to preserve the discontinuous or steeply varying gradients that may occur in sharp-fronted and trans-critical shallow water flows, see for instance (89; 5; 95; 32; 79; 13; 69) among others and also some references herein. Many of them particularly focus on the issue of balancing the flux gradient and the topography source source term (7; 78; 44; 70; 69; 20; 60; 23; 74). However, 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 for instance (58; 64; 84). In recent years, high-order discontinuous Galerkin (DG) methods have become very popular to approximate the solutions of various linear and nonlinear partial differential equations. Considering the approximation of hyperbolic conservation laws, DG methods combine the background of FE methods, FV methods and Riemann solvers, allowing to take into account the physic of the problem, and they have been successfully validated in many domain of applications. An arbitrary order of accuracy in space can be obtained with the use of high-order polynomials within elements, allowing to keep the stencil compact, along with being able to handle complex geometries through the use of unstructured general meshes and h/p-adaptivity. Moreover, they are highly parallelizable, and exhibit local conservation and strong stability properties. We refer the reader to (28; 29) for a general background. Several DG methods have been designed for the NSW equations since the early 2000s, see for instance (65; 106; 1; 4; 59; 12; 76; 90; 46; 41; 108; 107; 61) and some references hereafter. However, while DG methods may be mature enough to accurately handle some realistic problems in various applications, they still suffer from the lack of nonlinear stability. In particular, highorder DG methods may produce spurious oscillations in the presence of discontinuities or steeply varying gradients, i.e. Gibbs phenomenon. These nonphysical overshoots may lead to nonphysical solutions. Another challenging issue encountered in the approximation of the NSW equations is the preservation of the water height positivity, closely related to the issue of the occurrence and propagation of wet/dry fronts that may occur in dam breaks, flood waves or run-up over coastal shores. Considering the convex set of physical admissible states Θ, defined as follows:

Θ = (H, q) ∈ R 2 ; H ≥ 0 , (2)
a minimal nonlinear stability requirement is to ensure that this set is preserved at the discrete levels. But the use of high-order polynomials is generally not straightforwardly compatible with such a requirement and standard numerical methods may produce negative values for the water height H.

Generally speaking, robustness issues may be among the main remaining challenges for the use of high-order methods in realistic problems for many domains of applications, and in recent years, several approaches have been proposed to stabilize high-order approximations. These techniques mainly rely on two different paradigms that we referred to as a priori and a posteriori. In the so-called a priori framework, the correction procedure is applied before advancing the numerical piecewise polynomial solution further in time. So first, a troubled zone indicator is used to find where a correction is required (see [START_REF] Qiu | A comparison of troubled-cell indicators for Runge-Kutta discontinuous Galerkin methods using weighted essentially nonoscillatory limiters[END_REF] for a review of such troubled elements sensors). Then, sufficient efforts are made on the numerical solution or on the numerical scheme to be sure that one will be able to carry the computation out to the next time step. Among others a priori correction techniques, we could mention artificial viscosity methods (81; 97; 42; 49; 62), where some dissipative mechanism is added in shock regions, borrowing ideas from the streamline upwind Petrov Galerkin (SUPG) and Galerkin least-squares methods. Some other very popular limiting techniques can be gathered and referred to as slope and moment limiters (26; 14; 17; 63; 111; 57; 35; 66). In the former ones, as in (27; 26), the polynomial approximated solution is flattened around its mean value to control the solution jumps at cell interfaces. A smooth extrema detector is then generally used to prevent the limitation technique to spoil the accuracy in regions where no limiting is required. Moments limiters, mainly based on [START_REF] Biswas | Parallel, adaptive finite element methods for conservation laws[END_REF] and further developed in [START_REF] Burbeau | A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods[END_REF], can be seen as the extension of the aforementioned slope limiters to the case of very high orders of accuracy. In those limiting strategies, the different moments of the polynomial solution are successively scaled in a decreasing sequence, from the higher degree to the lower one, allowing the preservation of the solution accuracy, as well as ensuring the solution boundedness near discontinuities. The high-order DG limiter [START_REF] Krivodonova | Limiters for high-order discontinuous Galerkin methods[END_REF], generalized moment limiter [START_REF] Yang | A parameter-free generalized moment limiter for high-order methods on unstrucured grids[END_REF], hierarchical Multi-dimensional Limiting Process (MLP) (57; 56) and vertex-based hierarchical slope limiters (35; 66) all derive from (27; 14; 17), and thus fall into this category. Now, another limiting strategy that deserves to be mentioned is the (H)WENO limiting procedure (86; 8; 115; 67; 116), where the DG polynomial is substituted in troubled regions by a reconstructed (H)WENO polynomial. An alternative way to treat this spurious oscillations issue may be to use a solution filtering method, see for instance (98; 94; 75; 82), which aim at removing high wave-number oscillations. Those filtering procedures are generally done in an ad hoc fashion, filtering being applied "as little as possible, but as much as needed". Last but not least, some original subcell FV shock capturing techniques in the frame of DG schemes (51; 22; 92; 30) have recently gained in popularity. In [START_REF] Huerta | A simple shock-capturing technique for high-order discontinuous galerkin methods[END_REF], the authors use a convex combination between highorder DG schemes and first-order finite volumes on a sub-grid, allowing them to retain the very high accurate resolution of DG in smooth areas and ensuring the scheme robustness in the presence of shocks. Similarly, in (92; 30), after having detected the troubled zones, cells are then subdivided into subcells, and a robust first-order finite volume scheme is performed on the sub-grid in troubled cells.

The a priori paradigm has already and extensively proved in the past its high capability and feasibility, as in the aforementioned articles. Those techniques are a priori in the sense that only the data at time t n are needed to perform the limitation procedure. Then, the limited solution is used to advance the numerical scheme in time to t n+1 . The "worst case scenario" has to be generally considered as a precautionary principle. Furthermore, let us emphasize that most of the a priori correction procedures previously quoted do not ensure a maximum principle or the positivity-preservation of the solution. Generally, additional effort has to been made specifically on that matter, as for example by means of positivity-preserving limiters (112; 114). Specifically concerning the issue of positivity preservation in DG methods for the NSW equations, various a priori strategies have been introduced recently: a free-boundary treatment in mixed Eulerian-Lagrangian elements is introduced in (15) to locate the wet/dry interface, a fixed mesh method with a local conservative slope modification technique based on a redistribution of the fluid and cut-off in discharge is presented in [START_REF] Ern | A well-balanced Runge-Kutta discontinuous Galerkin method for the shallow-water equations with flooding and drying[END_REF], local first moment limitations without mass adding in (16; 60; 59; 99), high-order accuracy a priori polynomial reconstruction and limitation to enforce a strict maximum principle on mean values in (110; 109), in [START_REF] Duran | Recent advances on the discontinuous Galerkin method for shallow water equations with topography source terms[END_REF] for the so-called pre-balanced formulation of the NSW equations, or in (72) for a formulation with implicit time stepping. Let us finally mention [START_REF] Meister | A positivity preserving and well-balanced DG scheme using finite volume subcells in almost dry regions[END_REF] where a Finite-Volume subcell approaches has been adopted. Now, the paradigm of a posteriori correction is different in the way that first an uncorrected candidate solution is computed at the new time step. The candidate solution is then checked according to some criteria (for instance positivity, discrete maximum principle, . . . ). If the solution is considered admissible, we go further in time. Otherwise, we return to the previous time step and correct locally the numerical solution by making use of a more robust scheme. Because the troubled zone detection is performed a posteriori, the correction can be done only where it is absolutely necessary. Furthermore, let us emphasize that in a posteriori correction procedures, the maximum principle preservation or positivity preservation is included without any additional effort. Indeed, the whole procedure will be positivity-preserving as soon as the numerical scheme used as a correction procedure is. Consequently, all the a posteriori techniques that make use of FV scheme as correction method will then be positivity-preserving. Recently, some new a posteriori limitations have arisen.

Let us mention the so-called MOOD technique, (25; 33; 34). Through this procedure, the order of approximation of the numerical scheme is locally reduced in an a posteriori sequence until the solution becomes admissible. In (37; 36; 54), a subcell FV technique similar to the one presented in [START_REF] Sonntag | Shock capturing for discontinuous Galerkin methods using finite volume subcells[END_REF] has been applied to the a posteriori paradigm. Practically, if the numerical solution in a cell is detected as bad, the cell is then subdivided into subcells and a first-order FV, or alternatively other robust scheme (second-order TVD FV scheme, WENO scheme, . . . ), is applied on each subcell. Then, through these new subcell mean values, a high-order polynomial is reconstructed on the primal cell. Related strategies applied to dispersive and turbulent shallow water flows have been introduced in (18; 19).

Nonetheless, in all the aforementioned limitation techniques, a priori and a posteriori, in the troubled cells the high-order DG polynomial is either globally modified in the cell, or even discard as it is in the (H)WENO limiter or any a posteriori correction procedure. One of the main advantage of high-order scheme is to be able to use coarse grids while still being very precise. But even in the case where the troubled zone, as the vicinity a shock for instance, is very small regarding the characteristic length of a cell, the DG polynomial will be globally modified. In [START_REF] Vilar | A posteriori correction of high-order discontinuous galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF], we have introduced a new conservative technique to overcome this issue, by modifying the DG numerical solution only locally at the subcell scale. This correction procedure has been designed first to avoid the occurrence of non-admissible solution, to be maximum principle preserving, or in the context of systems positivity-preserving, and to prevent the code from crashing (for instance avoiding NaN in the code). The corrected scheme is also conservative, at the subcell level. Secondly, the correction permits to essentially avoid the appearance of spurious oscillations. Thirdly, it allows us to retain as much as possible the high accuracy and subcell resolution of DG schemes, by minimizing the number of subcells in which the solution has to be recomputed. Practically, the correction procedure presented in [START_REF] Vilar | A posteriori correction of high-order discontinuous galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF] only modifies the DG solution in troubled subcell regions without impacting the solution elsewhere in the cell. It is also worth mentioning that the whole procedure is totally parameter free, and behaves properly from 2nd order to any order of accuracy.

Making use of the a posteriori local subcell correction introduced [START_REF] Vilar | A posteriori correction of high-order discontinuous galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF], the main objective of this paper is to develop a novel shock-capturing, positivity-preserving and well-balanced DG method for the NSW equations equations with topography source term by using specific local flux correction at the subcell level, with a posteriori numerical admissibility detectors. The well-balanced property for NSW equations, first introduced in ( 11), has been widely studied in recent years. Following the ideas of ( 38), we use the so-called pre-balanced formulation of the NSW equations. Indeed, the alternative formulation of the NSW equations in deviatoric form, obtained by subtracting an equilibrium solution, and introduced in ( 89) is interesting as it leads to a balanced set of hyperbolic equations that does not require specific numerical algorithms to obtain a well-balanced property. However, such a formulation is given in terms of free surface elevation above the still water level, and is therefore not suitable to model cases involving dry areas (the still water depth is undefined in dry areas). In (68; 69), a new formulation given in terms of the total free surface elevation η = H +b (see Fig. 1) allows to alleviate this drawback. Indeed, observing that

1 2 g∂ x H 2 + gH∂ x b = 1 2 g∂ x (η 2 -2ηb) + gη∂ x b,
we obtain the so-called pre-balanced form of the NSW equations, given in a compact form:

∂ t v + ∂ x F(v, b) = B(v, ∂ x b), (3)
where v : R × R + → Θ is the vector of conservative variables, F : Θ × R → R 2 is the flux function and B : Θ × R → R 2 is the topography source term, defined as follows:

v = η q , F(v, b) = q uq + 1 2 g(η 2 -2 η b) , B(v, ∂ x b) = 0 -g η ∂ x b . (4)
In the next section, we introduce a new discrete formulation with an arbitrary order of accuracy for equations (4), following a classical DG formalism. Then, exploiting the fact that such DG formulation can be regarded as a FV-like scheme on a sub-partition with particular (reconstructed) high-order interface fluxes, we show that it is possible to slightly and locally modify these interface fluxes in order to enforce some nonlinear stability property through the use of some a posteriori admissibility sensors. We also show that such an approach is fully compatible with the preservation of motionless steady states at the subcell level, provided that local mean-values of the solution on the sub-partition are carefully reconstructed, adapting the ideas of (69; 38). In the third section, we assess the ability of this new method, called a posteriori Local Subcell Correction (LSC) in the following, to remove nonphysical instabilities in the vicinity of discontinuities, and ensure that the mean values of the water height on the sub-partition remain positive. We highlight, on several test cases, the particularly interesting subcell resolution ability of the resulting scheme. h ω , where h ω is the length of element ω. For a given mesh element ω i ∈ T h , we also note

ω i = x i-1 2 , x i+ 1 2
and by x i its barycenter.

Given an integer polynomial degree k ≥ 1, we consider the broken polynomial space

P k (T h ) := v ∈ L 2 (Ω), v |ω ∈ P k (ω), ∀ ω ∈ T h ,
where P k (ω) denotes the space of polynomials in ω of total degree at most k, with dim(P k (ω)) = k+1. Piecewise polynomial functions belonging to P k (T h ) are denoted with a subscript h in the following, and for any ω ∈ T h and v h ∈ P k (T h ), we may use the convenient shortcut: v ω := v h|ω when no confusion is possible.

For any mesh element ω ∈ T h and any integer k ≥ 0, we fix a basis for P k (ω) denoted by

Ψ ω = ψ ω j j∈ 1, k+1 .
A basis for the global space P k (T h ) is obtained by taking the Cartesian product of the basis for the local polynomial spaces:

Ψ h = × ω∈T h Ψ ω = ψ ω j j∈ 1, k+1 ω∈T h .
Note that we have:

supp(ψ ω j ) ⊂ ω, ∀ω ∈ T h , ∀j ∈ 1, k + 1 .
We introduce the following shortcut notations for smooth enough scalar-valued functions v, w:

v, w T h := ω∈T h v, w ω , v, w ω := ω v(x)w(x)dx, ∀ω ∈ T h , v ∂ω i := v(x i+ 1 2 ) -v(x i-1 2
), ∀ω i ∈ T h .

For ω ∈ T h , we denote p k ω the L 2 -orthogonal projector onto P k (ω) and p k T h the L 2 -orthogonal projector onto P k (T h ). Similarly, we denote I k ω the element nodal interpolator into P k (ω). The corresponding nodal distributions in elements are chosen to be the approximate optimal nodes introduced in [START_REF] Chen | Approximate optimal points for polynomial interpolation of real functions in an interval and in a triangle[END_REF], which have better approximation properties than equidistant distributions, and include, for each element, the elements boundaries into the interpolation nodes. The global I k T h interpolator into P k (T h ) is obtained by gathering the local interpolating polynomials defined on each elements.

We also define the broken gradient operator ∇ k h :

P k (T h ) → P k (T h ) such that, for all v h ∈ P k (T h ), (∇ k h v h ) |ω := ∂ x (v h|ω ), ∀ ω ∈ T h .
For any mesh element ω i ∈ T h , we introduce a sub-partition T ω i into k + 1 open disjoint subcells:

ω i = k+1 m=1 S ω i m ,
where the subcell

S ω i m = x ω i m-1 2 , x ω i m+ 1 2 is of size S ω i m = x ω i m+ 1 2 -x ω i m-1 2
, with the convention

x ω i 1 2 = x i-1 2 and x ω i k+ 3 2 = x i+ 1 2
, see Fig. 2. When considering a sequence of neighboring mesh elements ω i-1 , ω i , ω i+1 , the convenient convention S ω i 0 := S ω i-1

k+1 and S ω i k+2 := S ω i+1 1 may be used.

x i-1 2 x i+ 1 2 x ω i 1 2
x ω i To define the sub-resolution basis functions, required in § 2.3 and initially introduced in (101), we introduce for a given mesh element ω ∈ T h the following set of subcell indicator functions

k+ 3 2 x ω i 3 2 x ω i m-1 2 x ω i m+ 1 2 S ω i m
1 ω m , m ∈ 1, k + 1 , with: 1 ω m (x) = 1 if x ∈ S ω m , 0 if x ∈ S ω m , ∀m ∈ 1, k + 1 .
Then, the set of sub-resolution basis functions φ ω m ∈ P k (ω), m ∈ 1, k + 1 are defined as follows:

φ ω m = p k ω (1 ω m ), ∀m ∈ 1, k + 1 . (5)
Finally, for all ω ∈ T h we also introduce the set of piecewise constant functions on the sub-grid:

P 0 (T ω ) := v ∈ L 2 (ω), v |S ω m ∈ P 0 (S ω m ), ∀ S ω m ∈ T ω .
Concerning time discretization, for a given final computational time t max > 0, we consider a partition (t n ) 0≤n≤N of the time interval [0, t max ] with t 0 = 0, t N = t max and t n+1t n =: ∆t n . More details on the computation of the time step ∆t n and on the time marching algorithms are given in § 2.4.

For any sufficiently regular scalar-valued function of time w, we let w n := w(t n ).

Remark 1. The degrees of freedom are classically chosen to be the functionals that map a given discrete unknown belonging to P k (T h ) to the coefficients of its expansion in the selected basis. Specifically, the degrees of freedom applied to a given function v h ∈ P k (T h ) return the real numbers

v ω j with j ∈ 1, k + 1 and ω ∈ T h , (6)
such that

v ω = v h|ω = k+1 j=1 v ω j ψ ω j , ∀ω ∈ T h .
With a little abuse in terminology, we refer hereafter to the real numbers [START_REF] Arpia | Well balanced residual distribution for the ALE spherical shallow water equations on moving adaptive meshes[END_REF] as the degrees of freedom associated with v h and we note v ω ∈ R k+1 the vector gathering the degrees of freedom associated with v ω .

Remark 2. For any ω ∈ T h , and any v ω ∈ P k (ω), let denote

v ω m with m ∈ 1, k + 1 ,
the low-order piecewise constant components defined as the mean values of v ω on the subcells belonging to the subdivision T ω , called sub-mean values in the following, which may be gathered in a vector v ω ∈ R k+1 . Whenever a sequence of neighboring mesh elements ω i-1 , ω i , ω i+1 and associated neighboring approximations is considered, the following convenient convention may be used:

v ω i 0 := v ω i-1 k+1 and v ω i k+2 := v ω i+1 1
.

We observe that the degrees of freedom v ω m , m ∈ 1, k + 1 are uniquely defined through the sub-mean values v ω m , m ∈ 1, k + 1 , and reversely. Specifically, considering the local transformation matrix Π ω = π ω m,p m,p defined as:

π ω m,p = 1 |S ω m | S ω m ψ ω p dx, ∀ (m, p) ∈ 1, k + 1 2 , (7)
we have the following relation:

Π ω v ω = v ω and Π -1 ω v ω = v ω .
As a consequence, any polynomial function v h ∈ P k (ω) can be expressed equivalently either in terms of the degrees of freedom v ω , or the sub-means values v ω . Finally, let introduce the (one-to-one) following projector onto the piecewise constant sub-grid space:

π k Tω : P k (ω) -→ P 0 (T ω ) v ω -→ π k Tω (v ω ) := v ω . (8) (9)
In practice, once the transformation matrices Π ω are initialized in a preprocessing step, it is straightforward and computationally inexpensive to switch from one representation to another, see Fig. 3.

x ω i 1 2

x ω i 7) are obviously non-singular. This property would be straightforwardly extended to the multi-dimensional case with Cartesian grids.

k+ 3 2 x ω i 3 2 x ω i m-1 2 x ω i m+ 1 2 v ω i x i+ 1 2 v ω i k+1 v ω i 1 v ω i 2 x i-1 2

Discontinuous Galerkin formulation

Let b h = I k T h (b) denote a globally continuous piecewise polynomial approximation of the topography parametrization and let denote ∇b h = ∇ k h b h , for the sake of simplicity. A straightforward semidiscrete in space DG approximation of (3) reads: find

v h = (η h , q h ) ∈ (P k (T h )) 2 such that, for all ϕ ∈ P k (T h ), ∂ t v h , ϕ T h + A h (v h ), ϕ T h = 0, ( 10
)
where the discrete nonlinear operator A h is defined by

A h (v h ), ϕ T h := - ω∈T h F(v h , b h ), ∂ x ϕ ω + ω∈T h ϕ F ∂ω -B(v h , ∇b h ), ϕ T h , ∀ϕ ∈ P k (T h ). (11)
In [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF], F stands for the interface numerical flux function. Denoting by v - as follows:

F i+ 1 2 := F(v - i+ 1 2 , v + i+ 1 2 , b i+ 1 2 ), (12)
where the numerical flux function chosen here is the simple global Lax-Friedrichs flux:

F(v -, v + , b) := 1 2 F(v -, b) + F(v + , b) -σ(v + -v -) , (13)
with σ := max ω∈T h σ ω and

σ ω := max m |u ω m | + gH ω m .
Remark 4. We require that the volume integral and source term in formula [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF] are exactly computed at motionless steady states. This can be achieved, for the pre-balanced formulation (3)-( 4), by using any quadrature rule exact for polynomials of degree up to 2k, thanks to the pre-balanced formulation ( 3)-( 4). On the other hand, for the classical NSW formulation (1), a quadrature rule exact for polynomials of degree up to 3k is needed.

Remark 5. The topography approximation b h is interpolated by a piecewise polynomial but globally continuous function over the mesh. To achieve this, one can simply choose the elements boundaries among the interpolation points with any corresponding interpolation method. To ensure that the scheme is indeed well-balanced, and particularly in wet/dry context, see § 3.4, we initialize the surface elevation η h in dry areas by setting η h = b h . Thus, water height positivity is also assured in dry areas since h h = η hb h = 0, by construction. We emphasize that as long as h = ηb is non-negative, its subcell mean-values are also non-negative. However, nothing ensures that after performing a L 2 projection of the initial water height, the associated submean values of the L 2 projection h ω are positive. This is the reason why, in wet regions, for the initialization, we start by computing the positive h submean values using ( 14) and then reconstruct the associated polynomial using

Π -1 ω . h ω m = 1 |S ω m | S ω m h(x)dx. (14)
In the following, similarly to what has been done in [START_REF] Vilar | A posteriori correction of high-order discontinuous galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF], we demonstrate an equivalence relation between [START_REF] Bermudez | Upwind schemes for the twodimensional shallow water equations with variable depth using unstructured meshes[END_REF] and a FV-like method on a sub-mesh.

DG formulation as a FV-like scheme on a sub-grid

Let introduce the L 2 -projections of the flux function

F h = p k T h (F(v h , b h )) and of the source term B h = p k T h (B(v h , ∇b h ))
, such that:

F(v h , b h ), ϕ T h = F h , ϕ T h , ∀ϕ ∈ P k (T h ), B(v h , ∇b h ), ϕ T h = B h , ϕ T h , ∀ϕ ∈ P k (T h ). ( 15a
) (15b)
Remark 6. As we mentioned, in DG schemes, volume integral and source term contribution are computed using quadrature rule. This quadrature rule should be used to compute the left hand side of the L 2 projections (15a) and (15b).

From [START_REF] Bermudez | Upwind schemes for the twodimensional shallow water equations with variable depth using unstructured meshes[END_REF], we now have:

∂ t v h , ϕ T h - ω∈T h F h , ∂ x ϕ ω + ω∈T h ϕ F ∂ω -B h , ϕ T h = 0, ∀ϕ ∈ P k (T h ),
or equivalently, using an integration by parts:

∂ t v h , ϕ T h + ω∈T h ∂ x F h , ϕ ω - ω∈T h ϕ F h -F ∂ω -B h , ϕ T h = 0, ∀ϕ ∈ P k (T h ). (16)
Substituting φ ω m , defined in [START_REF] Anastasiou | Solution of the 2d shallow water equations using the finite volume method on unstructured triangular meshes[END_REF], into [START_REF] Bunya | A wetting and drying treatment for the runge-kutta discontinuous galerkin solution to the shallow water equations[END_REF] gives the local equations on mesh element ω ∈ T h :

∂ t v ω , φ ω m ω = -∂ x F ω , φ ω m ω + B ω , φ ω m ω + [(F ω -F) φ ω m ] ∂ω , ∀m ∈ 1, k + 1 .
Since ∂ t v ω , ∂ x F ω and B ω belong to (P k (ω)) 2 and considering [START_REF] Anastasiou | Solution of the 2d shallow water equations using the finite volume method on unstructured triangular meshes[END_REF], it follows that

∂ t v ω m = - 1 |S ω m | F ω ∂S ω m -φ ω m (F ω -F) ∂ω + B ω m , ∀m ∈ 1, k + 1 ,
where v ω m and B ω m are respectively the mean values of v ω and B ω on the subcell S ω m . Let introduce the k + 2 subcells interfaces fluxes F ω m+ 1 2 m∈ 0, k+1 such that:

F ω m+ 1 2 -F ω m-1 2 = F ω ∂S ω m -φ ω m (F ω -F) ∂ω , ∀m ∈ 1, k + 1 ,
so that we have

∂ t v ω m = - 1 |S ω m | F ω m+ 1 2 -F ω m-1 2 + B ω m . (17)
Formulation ( 17) can be seen as a FV-like scheme on subcell S ω m . The values F ω m+ 1 2 m∈ 0, k+1 will be thereafter referred to as reconstructed fluxes. Considering the mesh element

ω i = x i-1 2 , x i+ 1 2 ∈
T h , and setting the first and last reconstructed fluxes to the DG numerical flux values at cell boundaries such as:

F ω i 1 2 := F i-1 2 and F ω i k+3/2 := F i+ 1 2
, the m interior reconstructed fluxes expression is then given by:

F ω i m+ 1 2 = F ω i x ω i m+ 1 2 -C i-1 2 m+ 1 2 F ω i x i-1 2 -F i-1 2 -C i+ 1 2 m+ 1 2 F ω i x i+ 1 2 -F i+ 1 2 , ∀m ∈ 1, k + 1 , (18)
where the C

i± 1 2 m+ 1 2
are explicitly computed in [START_REF] Vilar | A posteriori correction of high-order discontinuous galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF]. Their expression is recalled here for the sake of completeness: 2 m∈ 0, k+1 , [START_REF] Busto | On thermodynamically compatible finite volume methods and path-conservative ader discontinuous galerkin schemes for turbulent shallow water flows[END_REF], gives the equivalence with DG schemes, other choices obviously lead to other schemes. For instance, if one set these constants to zero, except for the first and last to be one, one would then recover the spectral volume method, (105; 50). Let us also mention Flux Reconstruction schemes (FR), also referred to as Correction Procedure via Reconstruction (CPR), with which we share this reconstructed fluxes framework, see for instance (52; 103; 3; 45; 53) or the dedicated paragraph in [START_REF] Vilar | A posteriori correction of high-order discontinuous galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF] for more insight on the analogy between the present theory and Flux Reconstruction schemes.

C i-1 2 m+ 1 2 = k+1 p=m+1 φ ω i p x i-1 2 and C i+ 1 2 m+ 1 2 = m p=1 φ ω i p x i+ 1 2 . ( 19
)
Remark 9. The choice of the sub-partition points, x ω m+ 1 2 m∈ 0, k+1 , has already been discussed in [START_REF] Vilar | A posteriori correction of high-order discontinuous galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF]. It appeared that, regarding the reformulation of DG schemes into subcell Finite-Volume methods, the cell decomposition into subcells does not come into account, as any choice would lead to the same piecewise polynomial solution. However, for the correction procedure introduced in § 3, the sub-division does has a slight impact. Indeed, the use of a non-uniform sub-partition, for instance by means of the Gauss-Lobatto points, leads to better results compared to a uniform sub-division. This is more likely the manifestation of the Runge phenomenon in the context of histopolation, as the histopolation basis functions underlying the sub-mean value representation, are more oscillatory for a uniform cell sub-partition. Consequently, in the remainder, we make use of Gauss-Lobatto points to define the sub-partition points x ω m+ 1 2 m∈ 0, k+1 .

Time marching algorithm

Supplementing [START_REF] Bermudez | Upwind schemes for the twodimensional shallow water equations with variable depth using unstructured meshes[END_REF] with an initial datum v(0, •) = v 0 = (η 0 , q 0 ) t , the time stepping may be carried out using explicit SSP-RK schemes, (47; 91). For instance, writing the semi-discrete equation [START_REF] Bermudez | Upwind schemes for the twodimensional shallow water equations with variable depth using unstructured meshes[END_REF] in the operator form

∂ t v h + A h (v h ) = 0,
we advance from time level n to (n + 1) with the third-order scheme as follows:

v h n,1 = v h n -∆t n A h (v h n ), v h n,2 = 1 4 (3v h n + v h n,1 ) - 1 4 ∆t n A h (v h n,1 ) , v h n+1 = 1 3 (v h n + 2v h n,2 ) - 2 3 ∆t n A h (v h n,2 ) ,
where v h n,i , 1 ≤ i ≤ 2, are the solutions obtained at intermediate stages, ∆t n is obtained from the CFL condition [START_REF] Canestrelli | Well-balanced high-order centred schemes for non-conservative hyperbolic systems. applications to shallow water equations with fixed and mobile bed[END_REF], and the discrete initial data v h 0 is defined as the L 2 projection of the initial datum (see remark 5 for details). As the correction described in the following section make use of both DG scheme on the primal cells ω ∈ T h and FV scheme on the subcells S ω m ∈ T ω , the time step ∆t n is computed adaptively using the following CFL condition:

∆t n = min ω∈T h h ω 2k + 1 , min S ω m ∈Tω |S ω m | σ , ( 20
)
where σ is the constant previously introduced in the global Lax-Friedrichs numerical flux definition, (13).

A posteriori local subcell correction

In this section, we show how it is possible to modify the reconstructed fluxes F m+ 1 2 in a robust way in subcells where the uncorrected DG scheme [START_REF] Bermudez | Upwind schemes for the twodimensional shallow water equations with variable depth using unstructured meshes[END_REF] has failed, either by obtaining negative value for the water height or by generating nonphysical oscillations due to the Gibbs phenomenon in the vicinity of discontinuities. For sake of conciseness in notations, the superscript ω i may be avoided in the following when no confusion is possible.

As high-order Runge-Kutta SSP time marching algorithms may be regarded as convex combinations of first-order forward Euler schemes, we consider in the following, for sake of simplicity, a fully discrete formulation obtained from [START_REF] Bermudez | Upwind schemes for the twodimensional shallow water equations with variable depth using unstructured meshes[END_REF] and a first-order forward Euler scheme. We assume that at time level n the numerical solution v n h is admissible in a sense to be clarified later. We then compute an updated candidate solution v n+1 h through the uncorrected DG scheme [START_REF] Bermudez | Upwind schemes for the twodimensional shallow water equations with variable depth using unstructured meshes[END_REF]. If the candidate v n+1 h is admissible, no correction is needed. Otherwise, the uncorrected DG scheme has produced an updated solution v n+1 h which is not admissible on at least one particular mesh element cell ω i . Looking at the subcell level, and assuming that v n+1 ω is not admissible in the particular subcell S m ∈ T ω i , which is thus called a troubled subcell in the following, the main idea of our a posteriori Local Subcell Correction (LSC) is to replace the incriminated subcell mean value v n+1 m by a new one, denoted with a as follows v ,n+1 m , which is computed using a subcell first-order FV scheme of the form:

v ,n+1 m = v n m - ∆t n |S m | F l m+ 1 2 -F r m-1 2 + ∆t n B m , (21)
with some new subcell lowest-order corrected numerical fluxes

F l m+ 1 2 , F r m-1 2
which are defined hereafter. Indeed, because the uncorrected DG scheme ( 10) is equivalent to the subcell FV-like scheme [START_REF] Burbeau | A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods[END_REF] with high-order reconstructed fluxes [START_REF] Busto | On high order ader discontinuous galerkin schemes for first order hyperbolic reformulations of nonlinear dispersive systems[END_REF], we propose to substitute, at the boundaries of S m , the high-order reconstructed fluxes with first-order FV numerical fluxes. Finally, new degrees of freedom at discrete time t n+1 are computed from the modified set of sub-mean values, now given as a blend of uncorrected values v n+1 m and corrected values v ,n+1 m . This strategy is illustrated in Fig. 4, where the marked subcell is identified with red color. Additionally, to preserve the local conservation property of the resulting scheme, the left and right neighboring subcells, colored in green in Fig. 4, have to be updated too, even if they are flagged as admissible subcells, since we have substituted the reconstructed fluxes F m-1 2 and F m+ 1 2 with corrected ones. In the particular case depicted in Fig. 4 where S m-2 and S m+2 are also flagged as admissible, the sub-mean values v n+1 m+1 and v n+1 m-1 are thus replaced respectively by v ,n+1 m-1 and v ,n+1 m+1 computed through a high-order reconstructed flux on one end and a first-order FV numerical flux on the other end, as follows:

F k+ 3 2 F m-1 2 F1 2 F m+ 1 2 S ω m F3 2 F m+ 1 2 F m-1 2
v ,n+1 m-1 = v n m-1 - ∆t n |S m-1 | F l m-1 2 -F m-3/2 + ∆t n B m-1 , (22)
v ,n+1 m+1 = v n m+1 - ∆t n |S m+1 | F m+3/2 -F r m+ 1 2 + ∆t n B m+1 . (23)
For all the remaining admissible subcells (left in grey on Fig. 4), because the associated reconstructed fluxes are not corrected, they do not require any further computation, and the corresponding submean values are the values obtained through the uncorrected DG scheme, see Fig. 5.

S ω m v ⋆, n+1 m-1 v n+1 1 v ⋆, n+1 m v ⋆, n+1 m+1 v n+1 k+1
± m := b m± 1 2 -max 0, b m± 1 2 -η m , (24) b
- m+1 := b m+ 1 2 -max 0, b m+ 1 2 -η m , b + m-1 := b m-1 2 -max 0, b m-1 2 -η m . (25)
We introduce subcell's interfaces reconstructions for the water height as follows:

H ± m := max 0, η m -b m± 1 2 ,
and for the surface elevation and discharge:

η ± m := H ± m + b ± m , q ± m := H ± m q m H m , (26)
leading to the new subcell's interfaces values:

v ± m := (η ± m , q ± m ).
Using these reconstructed values, we introduce some new FV numerical fluxes on subcell's S m left and right interfaces, denoted by F r m-1 2 and F l m+ 1

2

, as follows: Remark 10. To compute the velocity in the vicinity of dry areas, we classically set a numerical threshold = 10 -8 to numerically define what "a dry cell" is and set the velocity to 0 if h < .

F l m+ 1 2 := F v + m , v - m+1 , b + m +   0 gη + m b + m -b x m+ 1 2   , (27)
F r m-1 2 := F v + m-1 , v - m , b - m +   0 gη - m b - m -b x m-1 2   , (28)
Remark 11. For the uncorrected DG scheme, nothing changes with respect to numerical fluxes and intermediate state variables. The DG fluxes F i-1/2 and F i+1/2 on cell ω i left and right interfaces x i-1/2 and x i+1/2 will always be in pure DG context

F i+1/2 = F v - i+1/2 , v + i+1/2 , b i+1/2 , and F i-1/2 = F v - i-1/2 , v + i-1/2 , b i-1/2
, as already defined in (12).

Flowchart

We summarize the proposed new a posteriori LSC method of DG schemes through the following flowchart:

1. starting from an admissible piecewise polynomial approximate solution v n h ∈ (P k (T h )) 2 , compute the candidate solution v n+1 h ∈ (P k (T h )) 2 using the uncorrected DG scheme (10), 2. for any mesh element ω ∈ T h , compute the candidate associated sub-mean values:

P 0 (T ω ) v n+1 ω = π Tω (v n+1 ω ),
3. for any mesh element ω ∈ T h , for any subcell S m ∈ T ω , check admissibility of the associated sub-mean values v n+1 m , and identify accordingly the sub-partition Otherwise:

T ω = T f ω ∪ T u ω ,
(a) for all ω ∈ T h such that T f ω = ∅, and all S m ∈ T f ω , substitute the corresponding reconstructed fluxes with some corrected fluxes defined in ( 27)-(A.3), as follows:

     F l m+ 1 2 ←-F l m+ 1 2 and F r m+ 1 2 ←-F r m+ 1 2 if either S m or S m+1 is marked, F l m+ 1 2 ←-F m+ 1 2 and F r m+ 1 2 ←-F m+ 1 2 otherwise, (b)
for all ω ∈ T h such that T f ω = ∅, and all S m ∈ T f ω , compute new sub-mean values for the marked subcells and their first neighboring subcells, respectively denoted

v * n+1 m , v * n+1 m-1 , v * n+1 m+1
, by means of a corrected subcell FV scheme as:

v ,n+1 p = v n p - ∆t n |S p | F l p+ 1 2 -F r p-1 2 + ∆t n B p , (29)
for p ∈ m -1, m + 1 . This subcell corrected method (29) will either falls in one of the previously introduced cases ( 21), ( 22) or ( 23 Step 3 of the flowchart is detailed in the next section.

Admissibility criteria

A large number of sensors or detectors have been introduced in the literature, to identify the marked subcells, where some kind of stabilization is required to avoid a loss of robustness. Following (101), we use two admissibility criteria: one for the Physical Admissibility Detection (PAD), another addressing the occurrence of spurious oscillations, namely the Subcell Numerical Admissibility Detection (SubNAD). This last criterion is supplemented with a relaxation procedure to exclude the smooth extrema from the troubled cells.

Physical Admissibility Detection (PAD)

Here, we define a sensor function that :

-Check if the sub-mean values v n+1 m belongs to Θ, see [START_REF] Alcrudo | A high-resolution godunov-type scheme in finite volumes for the 2d shallow-water equations[END_REF].

-Check if there is any N aN values.

Those are the minimum requirements to enforce the code robustness.

Subcell Numerical Admissibility Detection (SubNAD)

In order to tackle the issue of spurious oscillations near discontinuities, we enforce a local Discrete Maximum Principle (DMP), at the subcell level, on the surface elevation as follows:

-Check if, for m = 1, . . . , k + 1, the following inequalities hold:

min η n m-1 , η n m , η n m+1 ≤ η n+1 m ≤ max η n m-1 , η n m , η n m+1 .
The SubNAD criterion relies on a DMP based on subcell mean values, and not the whole polynomial set of values. Furthermore, as the neighboring subcells set used in the SubNAD is reduced to the first left and right subcells, and not all the subcells contained in the DG cell as well as in the the left and right first neighboring DG cells, see [START_REF] Vilar | A posteriori correction of high-order discontinuous galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF], one has to introduce a relaxation mechanism in order to preserve the scheme accuracy in the vicinity of smooth extrema.

Detection of smooth extrema.

In the relaxation procedure proposed in [START_REF] Vilar | A posteriori correction of high-order discontinuous galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF], it is assumed that the numerical solution exhibits a smooth extremon if at least the following linearized version of the surface elevation spatial derivative:

(∂ x η) lin ω i (x) = ∂ x η n+1 ω i + (x -x i ) ∂ xx η n+1 ω i ,
:= ∂ x η n+1 ω i - hω i 2 ∂ xx η n+1 ω i
to be the left boundary value of (∂ x η) lin ω i on cell ω i , as well as

∂ x η L min \ max = min \ max ∂ x η n+1 ω i-1 , ∂ x η n+1 ω i
respectively the minimum and maximum values of the mean derivative around x i-1 2 . We then define the left detection factor α L as follows:

α L =                min 1, ∂ x η L max -∂ x η n+1 ω i ∂ x η L -∂ x η n+1 ω i , if ∂ x η L > ∂ x η n+1 ω i , 1, if ∂ x η L = ∂ x η n+1 ω i , min 1, ∂ x η L min -∂ x η n+1 ω i ∂ x η L -∂ x η n+1 ω i , if ∂ x η L < ∂ x η n+1 ω i .
Introducing the symmetric values

∂ x η R min \ max = min \ max ∂ x η n+1 ω i , ∂ x η n+1 ω i+1 and ∂ x η R := ∂ x η n+1 ω i + hω i 2 ∂ xx η n+1
ω i , the right detection factor α R is obtained in a similar manner. Finally, introducing α := min (α L , α R ) , we consider that the numerical solution presents a smooth profile on the cell ω i if α = 1. In this particular case, the SubNAD criterion is relaxed, allowing the high-order accuracy preservation of smooth extrema.

Remark 12. One can apply the subcell numerical admissibility detection SubNAD and relaxation method detailed above on the Riemann invariants I ± = u ± 2 √ gH instead of the surface elevation η. Actually the simplest choice that also leads to the best results, is to make the detection on the surface elevation η. The detection on the Riemann invariants variables produces a more diffused solution.

Well-balancing property

This section is now dedicated to the demonstration of the well-balanced property of this a posteriori LSC of DG schemes.

Remark 13. Let us note that under the motionless steady-state assumption η h = η e and q h = 0, the following relation holds:

∂ x F(v ω , b ω ) = B(v ω , ∂ x b ω ), ∀ω ∈ T h .
Moreover, as we have

F(v h , b h ) = 0 1 2 gη 2 h -gη h b h ,
we emphasize that, under the steady-state hypothesis, F(v h , b h ) belongs to P k (T h ) 2 and not only to P 2k (T h ) 2 , since η h = η e . Therefore, at steady state,

F h := p k T h (F(v h , b h )) = F(v h , b h ).
As for B, under the same assumptions we have:

B(v h , ∇b h ) = 0 -gη e ∂ x b h ∈ P k (T h ) × P k-1 (T h ) ⊂ P k (T h ) 2 , thus, B h := p k T h (B(v h , ∇b h )) = B(v h , ∇b h ).
We have then the following result:

Proposition 1. The discrete formulation obtained by gathering [START_REF] Bermudez | Upwind schemes for the twodimensional shallow water equations with variable depth using unstructured meshes[END_REF] and the local corrected FV schemes on subcells ( 21), ( 22) and ( 23), together with a first-order Euler time-marching algorithm, preserves the motionless steady states, providing that the integrals of ( 11) are exactly computed for the motionless steady states. Specifically, for all n ≥ 0 and all η e ∈ R, (η n h = η e and q n h = 0) =⇒ η n+1 h = η e and q n+1 h = 0 .

Proof. We consider the scheme (17) on uncorrected subcells, and schemes ( 21), ( 22) and ( 23) on corrected subcells. We have to distinguish three different situations: (i) uncorrected subcell, (ii) neighbor of a marked subcell, (iii) marked subcell. We show in what follows that in all those situations, the corrected DG scheme preserves the motionless steady-states at the subcell level:

∀ω ∈ T h , ∀m ∈ [1, . . . , k + 1], η ω,n m = η e , q ω,n m = 0 =⇒ η ω,n+1 m = η e , q ω,n+1 m = 0.
1. Uncorrected subcell: S m-1 , S m and S m+1 are not marked.

In this case, we consider the uncorrected DG scheme or it's equivalent FV-like scheme with reconstructed fluxes:

v n+1 m = v n m - ∆t n |S m | F m+ 1 2 -F m-1 2 + ∆t n B m , (30)
with F m+ 1 2 and F m-1 2 defined in [START_REF] Busto | On high order ader discontinuous galerkin schemes for first order hyperbolic reformulations of nonlinear dispersive systems[END_REF]. We have, at steady state:

η + i± 1 2 = η - i± 1 2 = η e , q + i± 1 2 = q - i± 1 2 = 0, and b + i± 1 2 = b - i± 1 2 ,
and therefore

F i± 1 2 = F v h (x i± 1 2 ), b h (x i± 1 2
) , so that

F m± 1 2 = F v h ( x m± 1 2 ), b h ( x m± 1 2
) .

As a consequence, we have

F m+ 1 2 -F m-1 2 = Sm ∂ x F(v h , b h ) dx, (31)
and injecting ( 31) into [START_REF] De La Rosa | Hybrid DG/FV schemes for magnetohydrodynamics and relativistic hydrodynamics[END_REF] gives

v n+1 m = v n m - ∆t n |S m | Sm ∂ x F(v h , b h )dx + ∆t n |S m | Sm B(v h , ∇b h )dx = v n m .
2. Neighbor of a troubled subcell: S m , S m-1 are not marked and S m+1 is marked. The corresponding scheme, in this case, is the following:

v , n+1 m = v n m - ∆t n |S m | F l m+ 1 2 -F m-1 2 + ∆t n B m , with F l m+ 1 2 and F m-1 2
respectively defined in ( 27) and [START_REF] Busto | On high order ader discontinuous galerkin schemes for first order hyperbolic reformulations of nonlinear dispersive systems[END_REF]. At steady state, the reconstruction (26) yields η + m = η - m+1 = η c and q + m = q - m+1 = 0. It leads to:

F v + m , v - m+1 , b + m = 1 2 F(v + m , b + m ) + F(v - m+1 , b + m ) = 1 2 0 g (η e ) 2 -2η e b + m
, and then to:

F l m+ 1 2 = 1 2   0 g (η e ) 2 -2η e b x m+ 1 2   = F v h ( x m+ 1 2 ), b h ( x m+ 1 2 ) . (32)
Moreover, as in the previous case:

F m-1 2 = F v h ( x m-1 2 ), b h ( x m-1 2 ) . (33)
Gathering ( 32) and ( 33), we then have

F l m+ 1 2 -F m-1 2 = Sm ∂ x F(v h , b h ) dx, so that v , n+1 m = v n m .
3. Corrected subcell: S m is marked. In this case, the corresponding scheme reduces to [START_REF] Carrier | Water waves of finite amplitude on a sloping beach[END_REF]. Following the lines of the previous cases, we have:

F l m+ 1 2 = F v h ( x m+ 1 2 ), b h ( x m+ 1 2 ) , F r m-1 2 = F v h ( x m-1 2 ), b h ( x m-1 2 ) ,
and therefore

F l m+ 1 2 -F r m-1 2 = Sm ∂ x F(v h , b h ) dx, so that v , n+1 m = v n m
We have just shown that schemes ( 17)-( 21)-( 22)-( 23) do ensure the well-balanced property in wet subcells for all contexts, wet/wet and wet/dry. As for dry subcells, we can also simply show wellbalancing property. Considering a dry zone at time level n, under the assumptions η n = b n and q n = 0, one can easily show that the dry zone stays a dry zone a the next time level n + 1, i.e. η n+1 = b n and q n+1 = 0, by following a very similar procedure as in the previous proofs.

Remark 14. Note that the use of a non-smooth topography parameterization may be allowed, while still ensuring the well-balancing property, at the price of considering interface reconstructions also for the cells interfaces, for the DG scheme, in the spirit of (38; 69).

Preservation of the water height positivity

After computing the candidate solution v n+1 h through the uncorrected DG scheme [START_REF] Bermudez | Upwind schemes for the twodimensional shallow water equations with variable depth using unstructured meshes[END_REF], if we detect a negative sub-mean value on an arbitrary subcell, this subcell is then marked and a new (corrected) sub-mean value is evaluated by means of the first-order subcell FV scheme [START_REF] Carrier | Water waves of finite amplitude on a sloping beach[END_REF]. As a consequence, scheme [START_REF] Carrier | Water waves of finite amplitude on a sloping beach[END_REF] with reconstruction (26) should preserve positivity.

Proposition 2. Under the CFL condition (20), if ∀ω ∈ T h , ∀S m ∈ T ω , v ω, n m ∈ Θ, then ∀ω ∈ T h , ∀S m ∈ T ω , v ω, n+1 m ∈ Θ.
Proof. As the positivity-preserving property of our a posteriori LSC of DG schemes relies on the positivity of the first-order FV scheme used as the correction method, let us prove that if H n m and H n m±1 are non-negative, then scheme [START_REF] Carrier | Water waves of finite amplitude on a sloping beach[END_REF] does produce a water height H n+1 m also non-negative. Let us first recall the equation corresponding to the time evolution of the discrete surface elevation:

η n+1 m = η n m - ∆t n |S m | F 1 v n,+ m , v n,- m+1 , b + m -F 1 v n,+ m-1 , v n,- m , b - m , ( 34
)
where F 1 represents the first component of the numerical flux F and v n,± m , b n,± m are defined in ( 26) and [START_REF] Chen | Approximate optimal points for polynomial interpolation of real functions in an interval and in a triangle[END_REF]. For sake of simplicity, we drop in the following the superscript n. Equation ( 34) rewrites explicitly as:

η n+1 m = η m - ∆t n 2 |S m | H + m q m H m + H - m+1 q m+1 H m+1 -σ η - m+1 -η + m - ∆t n 2 |S m | H - m q m H m + H + m-1 q m-1 H m-1 -σ η - m -η + m-1 . ( 35
)
Noticing that

η - m+1 -η + m = H - m+1 -H + m as well as η - m -η + m-1 = H - m -H + m-1
, and subtracting b m on both sides of this last expression, equation ( 35) can be reformulated as:

H n+1 m = 1 - 1 2 λ (σ -u m ) H - m H m - 1 2 λ (σ + u m ) H + m H m H m + 1 2 λ (σ + u m-1 ) H + m-1 H m-1 H m-1 + 1 2 λ (σ -u m+1 ) H - m+1 H m+1 H m+1 , (36)
with Remark 15. The proposed positivity criteria are based on subcell values, and indeed, our goal is to show that our scheme preserves the positivity at the subcell level. Hence, the chosen strategy, as it is, does not ensure the pointwise positivity of h at some specific nodes: such a property is not needed. If one requires such pointwise positivity, for some specific reasons, we emphasize that an additional "positivity limiter", as the one provided in [START_REF] Xing | Positivity-preserving high order well-balanced discontinuous galerkin methods for the shallow water equations[END_REF] for instance, can be combined with our approach, to ensure the positivity of the polynomial solution at any chosen points.

u m = q m H m and λ = ∆t n |S m | . Therefore, H n+1

Numerical validations

In this numerical results section, we make use of several widely addressed and challenging test cases to demonstrate the performance and robustness of DG schemes provided the a posteriori local subcell correction presented. In all following test cases, if not stated differently, sub-mean values are displayed. It will allow us to fully illustrate the very precise subcell resolution of our scheme.

A new analytical solution for the NSW equations

This first test case aims at numerically evaluating the rates of convergence of the present a posteriori LSC of DG schemes. To do so, following the methodology introduced in (102) in the context of compressible gas dynamics, we make use here of a new manufactured smooth solution of the NSW equations. Details on the design of such solution can be found in Appendix A. This solution has the very interesting features to achieve any arbitrary regularity, i.e. v(., t) ∈ C Ns (Ω), ∀ t < t c (N s ) and any N s ∈ N * , allowing the study of convergence up to any order of accuracy, while involving almost vanishing water depth, together with a loss of regularity and the occurrence of discontinuous profiles for t ≥ t c .

We consider here the computational domain Ω = [-0.5, 2.5], and the particular case of N s = 3. It follows that the critical time reads t c ≈ 0.44 s, see Appendix A for more details. We initialize the problem with the following initial data:

η 0 = u 2 0 4 g and q 0 = u 3 0 4 g ,
with the following C Ns smooth initial velocity

u 0 (x) = 1 if x ≤ 0, e -x Ns+1 elsewhere.
While the uncorrected DG scheme [START_REF] Bermudez | Upwind schemes for the twodimensional shallow water equations with variable depth using unstructured meshes[END_REF] allows to compute the solution without any robustness issue for small enough values of time, nonphysical oscillations may be generated for larger values of time, leading to the activation of the a posteriori LSC method . A comparison between our fourth-order numerical solution computed on a mesh made of 60 cells, and the analytical solution at t = 0.1 s is shown on Fig. 6. One can see in Fig. 6 that only the cell mean values are displayed. We can also observe how the numerical scheme has very accurately captured to exact solution. In Table 1, we gather the global L 2 -errors as well as the rates of convergence for different order of approximation, computed on the surface elevation at t = 0.1s. As expected, the computed rates of convergence scale as O(k + 1). A similar behavior can be observed for the horizontal discharge q.

k In a second time, we consider a larger final computational time t > t c , so that a right-going discontinuity has developed from the initially regular profile, allowing to check the ability of the proposed a posteriori LSC method to stabilize the computation, namely to get rid of the spurious oscillations as well as enforcing the positivity of the water height. We run the previous case until t = 0.55, with k = 3 and n e = 100 mesh elements. Note that the standard DG method crashes in this case, since nonphysical undershoots would be rapidly amplified. In Fig. 7, a comparison between the a posteriori corrected DG solution and a reference solution obtained with a robust first-order FV method and n e = 10000 mesh elements. This is a challenging computation for high-order methods since small values of water height occur and thus small undershoots generally quickly lead to larger undershoots and possibly loss of positivity. In practice, the sensor starts to be activated when the strong gradient appears, slightly before the apparition of the discontinuity. A particular emphasize is put in Fig. 7 on the location of marked subcells, where uncorrected subcells are plotted by green dots while corrected subcells are plotted with blue squares. We observe that the particular combination of admissibility criteria introduced in §3.3 works quite well in practice, as the detection has been able to accurately track the moving front, and doing so removed the spurious oscillations without impacting smooth areas.

1 2 3 h E η L 2 q η L 2 E η L 2 q η L 2 E η L 2 q η L 2 1 15
To conclude this test, we show on Fig. 8 the numerical results obtained with the a posteriori LSC method with a high-order polynomial approximation k = 8, along with a quite coarse mesh made of 20 elements, at time t = 0.55 s. The use of such a coarse mesh permits to highlight the particularly interesting subcell resolution capabilities of our method, allowing to accurately locate the wet-dry interface inside a mesh element.

Dam-break

In this second test case, we focus on two dam-break problems over flat bottoms. The computational domain is set to Ω = [0, 1] and the first set of initial conditions is defined as follows:

η 0 (x) = 1 if x ≤ 0.5, 0.5 elsewhere, , q 0 = 0, b = 0.
The final time is set to t = 0.075 s. In Fig. 9, on a 50 cells mesh, fourth-order uncorrected DG solution is displayed on the left figure, while the corrected solution is plotted on the right one. This illustrates very clearly that even if the correction has been activated on in a very sharp area in the vicinity of the discontinuity, the solution has still been cleansed from its spurious oscillations. Now, we compare our a posteriori LSC method with the limitation process introduced in (110) (refereed to as PL/TVB method in what follows), which combines the positivity-preserving limiter [START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high-order schemes for conservation laws: survey and new developments[END_REF] with a standard TVB limiter [START_REF] Cockburn | The Runge-Kutta Discontinuous Galerkin Method for Conservation Laws V: Multidimensional Systems[END_REF]. Following [START_REF] Xing | Positivity-preserving high order well-balanced discontinuous galerkin methods for the shallow water equations[END_REF], the constant M involved in the TVB limiter is set to M = 0. The results are plotted in Fig. 10 and Fig. 11. In Fig. 10 and11, one can observe that the present correction technique outperforms the positivity-+ TVB limiter, both in the rarefaction and shock resolution. Finally, to demonstrate how this a posteriori LSC method scales going to very high-orders of accuracy and very coarse meshes, we run the same test with k = 9 and a 10 mesh elements. The corresponding numerical result is shown on Fig. 12. Fig. 12 illustrates the high capability of this a posteriori LSC method to retain the precise subcell resolution of discontinuous Galerkin schemes, allowing the use of very coarse meshes, along with being able to avoid the appearance of spurious oscillations or any unfortunate crash of the code. This figure also displays how the present correction affects the solution only at the subcell level, allowing the resolution of the shock in only one mesh element.

In a second time, we modify the initial conditions as follows:

η 0 (x) = 1 if x ≤ x 0 0 elsewhere , q 0 (x) = 0.
We compute the evolution up to t = 0.05 s, with k = 3 and n e = 50, in order to show the ability of the proposed method to compute the propagation of a wet/dry front. A comparison between the numerical results obtained with the a posteriori LSC method and the analytical solution is shown on Fig. 13 (left). Additionally, we compare these results with those obtained with the PL/TVB limitation process at the same times on Fig. 13 (right), together with zoomed profiles on Fig.

Well-balancing property

In this third test, we focus on the preservation of the motionless steady states. The computational domain is Ω = [0, 1]. The topography profile is defined as follows

b(x) =    A sin (x -x 1 ) • π 75 2 if x 1 ≤ x ≤ x 2 , 0 elsewhere, (37)
where A = 4.75, x 1 = 0.125 and x 2 = 0.875. The initial data is defined as

η 0 (x) = max (3, b(x))
and q 0 (x) = 0.

We evolve this initial configuration in time up to 100000 time iterations, with a fourth-order approximation and 120 mesh elements. The numerical results obtained with the a posteriori LSC method method are shown on Fig. 16. x (m) We highlight in Fig. 16 the particular marked cells, in which the correction has been performed. We emphasize that the steady state is effectively preserved up to the machine accuracy, validating numerically the compatibility of the a posteriori LSC method with the well-balancing property. A similar behavior is reported for other values of k and n e .

Next, we slightly modify the initial condition for the water height in order to have the bump totally submerged: η 0 (x) = 10 and q 0 (x) = 0.

We evolve this initial configuration in time up to 100, 000 time iterations, with a fourth-order approximation and 120 mesh elements. The numerical results obtained with the a posteriori LSC method are shown on Fig. In Table 2, we gather the global L 2 -errors obtained for several orders of approximation, for the surface elevation at t = 50s. As expected the steady state is preserved up to double precision accuracy.

k 1 2 3 h E η L 2 E η L 2 E η L 2 1 15
1.35E-15 1.12E-15 6.32E-16

Transcritical flow over a bump

We focus in this test on a classical transcritical flow without shock, see for instance [START_REF]Proceedings of the 2nd workshop on dam-break wave simulation[END_REF] for a complete description. The computational domain is Ω = [0, 25] (m). The topography profile is defined as follows: b(x) = 0.2 -0.05(x -10) 2 if 8 < x < 12, 0 elsewhere.

In this test, the incoming flow is enforced to be fluvial upstream and becomes torrential at the top of the bump. The initial data is defined as: η 0 (x) = 0.66 m and q 0 (x) = 0 m 2 .s -1 , and we prescribe the following boundary conditions: upstream: q = 1.53 m 2 .s -1 , downstream: h = 0.66 m while the flow is subcritical.

We run this test case with k = 3, n e = 100 and t = 200s. We show on Fig. 18

Carrier and Greenspan's transient solution

This test case, introduced in ( 21), describes the physical process in which the water level near the shoreline of a sloping beach is initially depressed, the fluid held motionless and then released at t = 0.

A transient wave is generated which runs up the beach, before returning to equilibrium state in a slow convergence process, reproducing some interesting conditions for assessing the robustness of the a posteriori LSC method in computing long waves run-up. In [START_REF] Carrier | Water waves of finite amplitude on a sloping beach[END_REF], a hodograph transformation is used to solve the NSW equations and obtain an analytical solution. The transformation makes use of two dimensionless variables (in the following, starred variables denote dimensionless quantities) σ * and λ * which are, respectively, a space-like and a time-like coordinate given by σ * = 4c * , λ * = 2 (u * + t * ) .

Let l be the typical length scale of this specific problem and α the beach slope. The scales used to obtain the nondimensionalized variables are:

x * = x/l, η * = η/(αl), u * = u/ gαl, t * = t/ l/αg, (38)
and the non-dimensional phase speed is given by:

c * = η * -x * . ( 39
)
The initial solution is specified by the following initial conditions:

η * 0 (σ * ) = e 1 - 5 2
a 3 (a 2 + σ * 2 ) 3 2 + 3 2 a 5 (a 2 + σ * 2 ) 5 2
, q * 0 (σ * ) = 0 and

x * = - σ * 2 16 + η * 0 , (40)
where a = 3 2 (1 + 0.9e)

1
2 and e is a small parameter which characterizes the surface elevation profile. The analytical solution is then given by  In view of result displayed in Fig. 19, one can see how accurate DG scheme along with our a posteriori LSC method is, as the numerical solution is extremely close to the exact solution and is able to simulate the return to the equilibrium state. This is due to the ability of our correction method to surgically modified the numerical solution only in the very few concerned subcells, as illustrated on Fig. 20. We finally assess the use of a high-order polynomial approximation (k = 8) on a very coarse mesh (n e = 10) to emphasize the very accurate and interesting subcell resolution ability of the proposed approach. The results obtained at t = 7 s are plotted on Fig. 21.

                      η * (σ, λ) = - u * 2 2 + eR e 1 -2 5/4 -iλ [(1 -iλ) 2 + σ 2 ] 3 2 + 3 2 (1 -iλ) 2 [(1 -iλ) 2 + σ 2 ] 5 2 , u * (σ, λ) = 8e a I m 1 [(1 -iλ) 2 + σ 2 ] 3 2 - 3 4 1 -iλ [(1 -iλ) 2 + σ 2 ] 5 2 , t * = 1 2 aλ -u * and x * = η * - a 2 σ 2 16

Carrier and Greenspan's periodic solution

In this test case, a monochromatic wave is let run-up and run-down on a plane beach. This solution represents the motion of a periodic wave of dimensionless amplitude A * and frequency ω * traveling shoreward and being reflected out to sea generating a standing wave on a plane beach. Recalling the dimensionless quantities [START_REF] Duran | Recent advances on the discontinuous Galerkin method for shallow water equations with topography source terms[END_REF] and [START_REF] Erduran | Performance of finite volume solutions to the shallow water equations with shock-capturing schemes[END_REF], the analytical solution is formulated as follows:

             u * = - A * J 1 (σ * ) sin (λ * ) σ * , η * = A * 4 J 0 (σ * ) cos (λ * ) - u * 2 4 , t * = 1 2 λ * -u * and x * = η * - σ * 2 16 ,
where J 0 and J 1 stand for the Bessel functions of zero and first order. We consider the solution obtained for A * = 0.6 and ω * = 1 (non-breaking wave), together with the length scale l = 20 m and a bottom slope α = 1/30. The value of this solution at t = 0 is supplied as initial condition, and similarly to the previous transient case, the analytical variations of the surface elevation at the left boundary is used as an offshore inlet boundary condition, generating the motion. We refer the reader to (21) for a complete description. We set k = 3 and n e = 50 and we compute the time evolution up to t = 1.5 T , where T is the time period of the periodic forcing. We show on Fig. 22 some snapshots of the free surface elevation plotted at several discrete time in the range [1.25T, 1.5T ] with the a posteriori LSC method , showing a very good agreement between the numerical solution and the analytical one. Additionally, we compare these results with those obtained with the PL/TVB method on Fig. 23 with M = 0 (left) and with M = 32 (right). Let us note that in [START_REF] Xing | Positivity-preserving high order well-balanced discontinuous galerkin methods for the shallow water equations[END_REF], the authors make use of M = 0 in every situations, except for the convergence rate analysis where M = 32 is used. As this test case is for the most part smooth (except at the wet/dry transition point), a non-zero value of M can be used in order to improve the quality of the results, as depicted by Fig. 23. However, even for higher value of M , the PL/TVB limiter is outperformed by the present a posteriori LSC method . In order to emphasize the accuracy of the proposed approach for long time integration, we set t = 15T and show on Fig. 24 the free surface elevation obtained at times t = 14.5T (left) and t = 15T (right), for k = 3 and n e = 50. We observe that such a long time-integration has a negligible impact on the accuracy of the predictions of the shoreline location. Such a result can be reproduced with a high-order approximation k = 8 and a very coarse mesh n e = 10, showing again the ability of our approach to provide a high-order accurate subcell description of the motion, see Fig. 26. In Fig. 25, we show time-series of the shoreline elevation "η s " in the range [0, 6T ] (left) and [0, 15T ] (right). We can see that the minimum and maximum water elevations are accurately computed, even after a large number of periods.

-0.14 -0.12 -0.

Run-up of a solitary wave on a plane beach

The last test case is devoted to the computation of the run-up of a solitary wave on a constant slope. Such run-up phenomena are investigated experimentally and numerically in [START_REF] Synolakis | Standards, criteria, and procedures for noaa evaluation of tsunami numerical models[END_REF]. In this test, a solitary wave traveling from the shoreward is let run-up and run-down on a plane beach, before being fully reflected and evacuated from the computational domain. The topography is made of a constant depth area juxtaposed with a plane sloping beach of constant slope α such that cot(α) = 19.85. The right boundary condition is transmissive. The initial condition is defined as follows:

η 0 (x) = A H 0 sech 2 (γ (xx 1 )) and u 0 (x) = g H 0 η 0 (x), where γ = 3A 4H 0 and x 1 = 4H 0 3A arcosh 1 0.05 is nothing but the initial position of the center of the solitary wave. This test is run with A = 0.019 m, H 0 = 1.0 m, k = 8, n e = 20 and t = 40 s. We show on Fig. 27 the free surface obtained with the a posteriori LSC method at several times in the range [1 s, 40 s], showing once more a very good agreement with the reference solution obtained with a robust FV method on a very fine mesh n e = 10000.

Conclusion

In this paper, we have introduced a new well-balanced high-order discontinuous Galerkin discrete formulation with a Finite-Volume subcell correction patch designed for the NSW equations. This formulation, based on (101), combines the very high accuracy of DG schemes along with a robust correction procedure ensuring the water height positivity as well as addressing the issue of spurious oscillations in the vicinity of discontinuities. This robustness is enforced by means of an a posteriori local subcell correction of the conservative variables. This procedure relies on an advantageous reformulation of DG schemes as a FV-like method on a sub-grid, which makes the correction strategy surgical and flexible, as well as conservative at the subcell level. Indeed, only the non-admissible subcells are marked and subject to correction, retaining as much as possible the very accurate subcell resolution of high-order DG formulations. The proposed strategy is investigated through an extensive set of benchmarks, including a brand new smooth solution for the computation of convergence rates, stabilization of flows with discontinuities, the preservation of motionless steady states, or moving shorelines over varying bottoms. We observe in particular that this approach provides a very accurate description of wet/dry interfaces even with the use of very high-order schemes on coarse meshes. Regarding potential advantages of this a posteriori limiting strategy compared to a priori limiters, because the troubled zone detection is performed a posteriori, the correction can be done only where it is absolutely necessary. Furthermore, positivity preservation of the water height is included without any additional effort, while it is generally not the case of a priori limitations of high-order schemes. Let us further emphasize that this a posteriori LSC method scalability to any order of accuracy is also perfectly natural. Finally, it is important to note that this new correction procedure is totally parameter free. In the future, we have the desire to extend this a posteriori correction technique to a general 2D case on unstructured grids for the Shallow-water system with source term. We also plan to investigate the moving mesh case, based on an arbitrary-Lagrangian-Eulerian (ALE) formalism, in the context of a coupling with a floating object.

Appendix A. New particular smooth solution for the NSW equations

This appendix aims at giving further details on the construction of a new smooth solution, of any arbitrary regularity, of the NSW equations. Following the methodology introduced in (100), we consider a smooth solution v in the context of flat bottom (b = 0), so that the NSW equations rewrite as:

∂ t v + A(v)∂ x v = 0,
where the Jacobian matrix writes as:

A(v) = ∇ v F(v) = 0 1 gH -u 2 2u
.

The eigen-analysis of the matrix A(v) leads to the following pair of eigenvalues λ ± = u ± √ gH and eigenvectors:

E ± = 1 u ± √ gH .
By diagonalizing the NSW system of equations, one finally gets the following Riemann invariants α ± = u ± 2 √ gH, governed by the following conservation laws:

∂ t α ± + λ ± ∂ x α ± = 0. (A.1)
In light of the definition of the Riemann invariants, the system eigenvalues can be reformulated in terms of α ± as follows:

λ ± = α + (2 ± 1) + α -(2 ∓ 1) 4 .
To uncouple the two conservation laws (A.1), we consider a particular flow regime corresponding to the trans-critical particular situation where α -= 0, i.e. u = 2 √ gH. The NSW equations then finally reduce to the following very simple Burgers equation:

Figure 1 :

 1 Figure 1: Free surface flow: main notations

Figure 2 :

 2 Figure 2: Partition of a mesh element ωi in k + 1 subcells

Figure 3 :Remark 3 .

 33 Figure 3: Piecewise polynomial function and associated sub-mean values

left and right traces of v h on interface x i+ 1 2 , and by b i+ 1 2 = b - i+ 1 2 = b + i+ 1 2 2 on interface x i+ 1 2

 222122 the trace of b h , we define the numerical flux function F i+ 1

Figure 4 :

 4 Figure 4: Sketch of the correction of the reconstructed fluxes at subcell boundaries

Figure 5 :/r m± 1 2 .m+ 1 2 : 2 :

 5222 Figure 5: Sketch of sub-mean values before and after correction

where b x m± 1 2 2 and x m-1 2 .

 122 are respectively the interpolated polynomial values of b h at x m+ 1

= π - 1

 1 where T f ω and T u ω respectively refer to the set of flagged (non-admissible) subcells and the set of nonflagged (admissible) subcells (note that the sub-mean values v n+1 m may be obtained either from Step 2. (without correction) or from Step 4. (b) (after correction)), 4. if, for all ω ∈ T h , the identity T u ω = T ω holds, then for all ω ∈ T h , v n+1 ω Tω (v n+1 ω ) is admissible, no additional correction is required and we can go further in time: go to Step 1, starting from v n+1 h instead of v n h .

  ), (c) for all ω ∈ T h such that at least one subcell has been corrected, gather the uncorrected sub-mean values v n+1 m and corrected sub-mean values v * n+1 m in a new element of P 0 (T ω ), which is still denoted v n+1 ω for the sake of simplicity, (d) go to step 3,

  m reads as a convex combination of H m-1 , H m and H m+1 . Furthermore, since by construction 0 ≤ H ± p ≤ H p , ∀p ∈ 1, k + 1 and by respect of the CFL condition (20), λα ≤ 1, and then all the coefficients involved in the convex combination (36) are non-negative. It follows that H n+1 m ≥ 0.

Figure 6 :

 6 Figure 6: Test 1 -A new analytical solution for the NSW equations -Free surface elevation computed at t = 0.1 s with the a posteriori LSC method for k = 3 and ne = 60.

Figure 7 :

 7 Figure 7: Test 1 -A new analytical solution for the NSW equations -Free surface elevation computed at t = 0.55 s with the a posteriori LSC method (left) for k = 3 and ne = 100, with a zoom on the discontinuity and wet/dry interface (right).

Figure 8 :

 8 Figure 8: Test 1 -A new analytical solution for the NSW equations -Free surface elevation computed at t = 0.55 s with the a posteriori LSC method (left) for k = 8 and ne = 20, with a zoom on the discontinuity and wet/dry interface (right).

Figure 9 :

 9 Figure 9: Test 2 -Dam break on a wet bottom-Free surface elevation computed at t = 0.075 s with the uncorrected DG method (left) and the a posteriori LSC method (right), with k = 3 and ne = 50 mesh elements.

Figure 10 :Figure 11 :

 1011 Figure 10: Test 2 -Dam break on a wet bottom -Free surface elevation computed at t = 0.075 s -Comparison between a posteriori LSC method and PL/TVB method for k = 3 and ne = 50.

Figure 12 :

 12 Figure 12: Test 2 -Dam break on a wet bottom -Free surface elevation computed at t = 0.075 s -Comparison between a posteriori LSC method (right) and PL/TVB method (left) for k = 9 and ne = 10.

Figure 13 :Figure 14 :Figure est 2 -

 13142 Figure 13: Test 2 -Dam break on a dry bottom -Free surface elevation computed at different times between 0.002s and 0.05 s -Comparison between a posteriori LSC method (left) and PL/TVB method (right) for k = 3 and ne = 50.

Figure 16 :

 16 Figure 16: Test 3 -Preservation of a motionless steady state -Free surface elevation at t = 50s (left), with a zoom on wet/dry interface (right).

17

 17

  .

Figure 17 :

 17 Figure 17: Test 3 -Preservation of a motionless steady state -Free surface elevation at t = 50s.

Figure Test 4 -

 4 Figure Test 4 -Transcritical flow without shock -Free surface elevation and discharge computed at several moments, 3.55s, 20.3s and 200s, with the a posteriori LSC method , for k = 3 and ne = 100.

,Figure 19 :

 19 Figure 19: Test 4 -Carrier and Greenspan's transient solution -Free surface elevation η * /β plotted versus the onshore coordinate x * -Free surface elevation for different values of time in the range [0.5 s, 23 s] (left) and at t = 200 s (right) for k = 3 and ne = 50.

Figure 20 :

 20 Figure 20: Test 4 -Carrier and Greenspan's transient solution -Free surface elevation computed at t = 7 s with the a posteriori LSC method for k = 3 and ne = 50 (left): corrected and uncorrected subcells are respectively plotted with blue squares and green dots, with a zoom on the shoreline (right)

Figure 21 :

 21 Figure 21: Test 4 -Carrier and Greenspan's transient solution -Free surface elevation computed at t = 7 s with the a posteriori LSC method for k = 8 and ne = 10.

Figure 22 : 50 -Figure 23 :

 225023 Figure 22: Test 5 Carrier and Greenspan's periodic solution -Free surface elevation computed for different values of time in the range [1.25T, 1.5T ] with the a posteriori LSC method for k = 3 and ne = 50

Figure 24 :Figure 25 :Figure 26 :

 242526 Figure 24: Test 5 -Carrier and Greenspan's periodic solution -Free surface elevation computed at t = 14.5T (left) and t = 15T (right) with the a posteriori LSC method for k = 3 and ne = 50.

Figure 27 :

 27 Figure 27: Test 6 -Run-up of a solitary wave on a plane beach -Free surface elevation computed for different values of time in the range [1 s, t = 40 s] with the a posteriori LSC method obtained for k = 8 and ne = 20.

  R denote an open segment with boundary ∂Ω. We consider a partition T h = ω 1 , . . . , ω ne of Ω in open disjoint segments ω of boundary ∂ω such that Ω = ω∈T h ω. The partition is characterized by the mesh size h := max

	2. Discrete formulations
	2.1. Discrete settings
	Let Ω ⊂ ω∈T h

  Remark 7. One can see that the reconstructed flux is nothing but the polynomial interior flux F ω , plus some correction terms taking into account the jump in fluxes at cell boundary ∂ω. A simple explicit expression of the correction coefficients[START_REF] Busto | On thermodynamically compatible finite volume methods and path-conservative ader discontinuous galerkin schemes for turbulent shallow water flows[END_REF] only depending on k and the subcells interfaces coordinates x ω Let us note that if this particular definition of C

	2 Remark 8. i± 1 m+ 1
	2 m+ 1	is given in (101).

  are respectively the mean values of (∂ x η h ) |ω i and (∂ xx η h ) |ω i on mesh element ω i . In practice, the DMP relaxation used here works as a vertexbased limiter on (∂ x η) lin ω i . Hence, we set ∂ x η L

	has a monotonous profile, where ∂ x η n+1 ω i	and ∂ xx η n+1 ω i

Table 1 :

 1 Test 1 -A new analytical solution for the NSW equations: L 2 -errors between numerical and analytical solutions and convergence rates for η at time t = 0.1s

		5.91E-4 1.96 2.13E-5 3.19	3.20E-6 4.05
	1 30 1 60 1 120	1.52E-4 2.02 2.33E-6 2.85 3.73E-5 2.02 2.99E-7 2.95 9.21E-6 -4.18E-8 -	1.93E-7 4.18 1.06E-8 3.95 6.91E-10 -

Table 2 :

 2 Test 3 -Preservation of a motionless steady state: L 2 -errors between numerical and exact steady state solutions for η at time t = 50s.

Acknowledgments. F. Marche gratefully acknowledges the support provided under Agence Nationale de la Recherche (ANR) project NABUCO (ANR-17-CE40-0025) and CNRS-LEFE-MANU project D-Wave.

To design a C Ns smooth solution, we initialize the problem with the following initial data:

with the following C Ns smooth initial velocity

The method of characteristics provides us with the expression of the analytical solution, for any given t ∈ [0, t c [:

where the characteristic lines read x(X, t) = 3 2 e -X Ns+1 t + X. For practical applications, to assess the position of the characteristic line origin point X given x and t, one may use an iterative rootfinding process, as Newton's method, to solve the non-linear problem g(X) = 0, where for given x and t function g(X) = 3 2 e -X Ns+1 t + Xx.

The analytical solution, (A.2), is defined ∀ t ∈ [0, t c [, where the critical time at which the characteristic lines cross is defined as follows: