A Duran
email: arnaud.duran@math.univ-montp2.fr

A robust and Well Balanced scheme for the 2D Saint-Venant system on unstructured meshes with friction source term

In the following lines we propose a numerical scheme for the shallow water system supplemented by topography and friction source terms, in a 2d unstructured context. This work proposes an improved version of the well-balanced and robust numerical model recently introduced in [J. Comp. Phys., 235, 565-586, 2013] for the pre-balanced Shallow Water Equations, accounting for varying topography. The present work aims at relaxing the robustness condition and include a friction term. To this purpose, the scheme is modified using a recent method, entirely based on a modified Riemann solver. This approach preserves the robustness and well-balanced properties of the original scheme, and prevents from unstable computations in the presence of low water depths. A series of numerical experiments is devoted to highlight the performances of the resulting scheme. Simulations involving dry areas, complex geometry and topography are proposed to validate the stability of the numerical model in the neighbourhood of wet/dry transitions.

Introduction

The description of hydrodynamic processes is present in a wide variety of scientific fields such as oceanography, river modelling or hydraulic engineering. From a practical point of view, the simulation of tsunami waves, dam break problems or pollutant transport in high density population regions for example has become a major concern. This motivates the introduction of numerical tools providing accurate descriptions of these phenomena, which are notably able to handle correctly dry domains, and are computationally efficient in the perspective of large scale applications. A relevant modelization of such environmental flows can be obtained under hydrostatic and shallow water assumptions in the water waves equations, leading to the well known system of Shallow Water Equations (SWE). Under its conservative form, the SWE consists in a set of PDE's involving the total water height h and the discharge q = t (q x , q y ) as vector variable. In a 2d context, denoting z the bottom height and u = t (u, v) the velocity vector, these equations are commonly written as follows:

U t + ∇ • G(U ) = B(U, z) + E(U ) , (1)
with

U =   h q x q y   , G(U ) =   q x q y uq x + 1 2 gh 2 vq x uq y vq y + 1 2 gh 2   , B(U, z) =    0 -ghz x -ghz y    . (2)
A sketch of a 1d configuration is proposed in Fig. 1, where we also define the new variable η as the total free surface (h + z). In what follows, Θ := {(h, q) , h ≥ 0 , q ∈ R 2 }

(3) will denote the convex set of admissible states.

h(x, t) z(x) η(x, t) = h(x, t) + z(x) u(x, t) As regards the friction source term, the proposed strategy allows some latitude concerning the choice of the law. Thus, the present work will involve the three classical formulations:

Manning : E(U ) = -      0 n 2 q h 10/3 q x n 2 q h 10/3 q y      , (4a)
Darcy :

E(U ) = -      0 d 8 q h 2 q x d 8 q h 2 q y      , (4b)
Linear :

E(U ) = -   0 κq x κq y   , (4c)
where n, d and κ are positive parameters.

Note that in equations [START_REF] Anastasiou | Solution of the 2d shallow water equations using the finite volume method on unstructured triangular meshes[END_REF], some minor physical effects have been neglected, as the source terms related to the wind or rainfall contribution for instance. Numerical schemes taking such phenomena into account can be found in [START_REF] Cea | Unstructured finite volume discretization of bed friction and convective flux in solute transport models linked to the shallow water equations[END_REF][START_REF] Chippada | A Godunov-type finite volume method for the system of shallow water equations[END_REF][START_REF] Delestre | Simulation du ruissellement d'eau de pluie sur des surfaces agricoles[END_REF], or even [START_REF] Delestre | A numerical scheme for a viscous shallow water model with friction[END_REF], where a new shallow water model is introduced, including viscosity. In this work, we also choose to set the Coriolis force to zero; numerical investigations shown that their impact was insignificant within the considered scales. The inclusion of such long scale physical effects is performed in [START_REF] Bokhove | Flooding and drying in finite-element discretizations of shallow-water equations. part 2: Two dimensions[END_REF][START_REF] Bouchut | Frontal geostrophic adjustment and nonlinear wave phenomena in one-dimensional rotating shallow water. part 2. high-resolution numerical simulations[END_REF][START_REF] Hassan | Numerical solution of the rotating shallow water flows with topography using the fractional steps method[END_REF][START_REF] Lukacova-Medvidova | Well-balanced finite volume evolution Galerkin methods for the shallow water equations[END_REF] for instance.

Regarding the structure of the system (1), setting u n = un x + vn y for a given normal vector n = (n x , n y ), the Jacobian matrix G n := ∂(G(U ). n) ∂U admits u n -c, u n , u n + c as eigenvalues (with c = √ gh), all distinct if h = 0, in which case the system is strictly hyperbolic. As a consequence, our work naturally follows to the general orientations issuing from the analysis and resolution of such class of conservation law problems (see [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF][START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[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] Gosse | A well-balanced scheme using non-conservative products designed for hyperbolic systems of conservation laws with source terms[END_REF][START_REF] Leveque | A class of approximate riemann solvers and their relation to relaxation schemes[END_REF][START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF] for some reference studies). These works generally rely on Finite Volume Methods (FVM), often preferred for their natural ability in capturing shocks, and an interesting latitude regarding the treatment of the fluxes and source terms.

Recently, a large number of theoretical and numerical investigations are run in order to develop suitable approximations for weak solutions of the shallow water system. In addition of the necessity of handling shocks that may arise from their hyperbolic structure, several additional concerns have to be taken into account when seeking approximations of these equations. A main issue stands in the so called "C-property" [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF], meaning that the numerical solution is able to preserve the particular steady state: η = cte and u = 0 , [START_REF] Benkhaldoun | A new finite volume method for flux-gradient and source-term balancing in shallow water equations[END_REF] usually called lake at rest. Numerical schemes satisfying this property are often referred to as well-balanced schemes. Secondly, an important problem commonly broached in the literature concerns the so called robustness property, which signify that the approximate solution stays in the set of admissible states (3) at each time step. The capacity of preserving the positivity of the water height is fundamental to avoid the possible appearance of spurious oscillations or unbounded computations entailed by the production of non physical solutions. Another central issue deserving to be highlighted, and closely related to robustness, is the ability to guarantee stable computations in the vicinity of dry areas, especially if the model pretends to deal with a friction term under the form (4a) or (4b), that involves divisions by h. The combination of these constraints stands for a discriminating and challenging objective in the construction of a numerical model for SWE, especially in the framework of 2d unstructured meshes.

There is a large range of numerical techniques to obtain approximations of weak solutions for SWE. Unfortunately, none of them may directly lead to a numerical scheme that satisfies the conditions mentioned above, and some modifications consistently need to be introduced. For example, the C-property requires an appropriate discretization of the bed slope source term, which exactly counteracts with the contribution of the numerical fluxes. Since the pioneering work of Bermudez and Vazquez [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF], many techniques have been introduced to overcome such difficulties, notably the hydrostatic reconstruction, initially proposed in [2] and extended on unstructured triangulations [START_REF] Audusse | A well balanced positivity preserving "second-order"scheme for shallow water flows on unstructured meshes[END_REF]. This method, which also preserves the water height positivity, inspired many authors in a large variety of different contexts (see [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] for some 1d applications, and [START_REF] Marche | Evaluation of well-balanced borecapturing schemes for 2D wetting and drying processes[END_REF][START_REF] Noelle | Well-balanced finite volume schemes of arbitrary order of accuracy for shallow water flows[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 2d applications on cartesian grids). In an unstructured framework, FVM [START_REF] Dutykh | The VOLNA code for the numerical modelling of tsunami waves: generation, propagation and inundation[END_REF] or even discontinuous Galerkin schemes [START_REF] Djadel | A discontinuous Galerkin method for the shallow water equation with bathymetric source terms and dry areas[END_REF][START_REF] Ern | A well-balanced Runge-Kutta discontinuous Galerkin method for the shallow-water equations with flooding and drying[END_REF] can also employ this technique. Of course, many other works have been proposed in order to satisfy both robustness and well-balanced properties. Concerning the 1d SWE, some interesting ideas can be found in [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] Gallouët | Some approximate godunov schemes to compute shallow-water equations with topography[END_REF][START_REF] Perthame | A kinetic scheme for the Saint-Venant system with a source term[END_REF][START_REF] Xing | High order finite difference WENO schemes with the exact conservation property for the shallow water equations[END_REF]. Considering this numerical challenge in a 2d context, one can mention the approaches available in [START_REF] Bradford | Finite-volume model for shallow-water flooding of arbitrary topography[END_REF][START_REF] Brodtkorb | Efficient shallow water simulations on GPUs : Implementation, visualization, verification and validation[END_REF][START_REF] Hubbard | Flux difference splitting and the balancing of source terms and flux gradients[END_REF][START_REF] Hubbard | A 2d numerical model of wave run-up and overtopping[END_REF][START_REF] Kesserwani | Locally limited and fully conserved RKDG2 shallow water solutions with wetting and drying[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] Lukacova-Medvidova | Well-balanced finite volume evolution Galerkin methods for the shallow water equations[END_REF] for cartesian grids, and [START_REF] Audusse | Kinetic Schemes for Saint-Venant Equations with Source Terms on Unstructured Grids[END_REF][START_REF] Benkhaldoun | A new finite volume method for flux-gradient and source-term balancing in shallow water equations[END_REF][START_REF] Berthon | Efficient well balanced hydrostatic upwind schemes for shallow water equations[END_REF][START_REF] Bokhove | Flooding and drying in finite-element discretizations of shallow-water equations. part 2: Two dimensions[END_REF][START_REF] Bollermann | Finite volume evolution Galerkin methods for the shallow water equations with dry beds[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] Xing | Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water equations on unstructured triangular meshes[END_REF] on unstructured meshes.

An interesting idea to minimize efforts in obtaining well-balanced and robustness properties, not mentioned so far, consists in considering the total free surface η instead of h as a flow variable. This method finds certainly its roots in the works of Zhou [START_REF] Zhou | The surface gradient method for the treatment of source terms in the shallow-water equations[END_REF] with the so called Surface Gradient Methods, in which the free surface is preferred to the water height in the evaluation of the fluxes contributions within a Godunov type second order scheme. The strategy provides a balanced scheme by mean of a simple centred approach of the source. With the aim of preserving quiescent flows, the idea of exploiting the invariance of η during lake at rest situations has been subsequently followed by several authors, notably Xing and Shu [START_REF] Xing | High order finite difference WENO schemes with the exact conservation property for the shallow water equations[END_REF] which suggested a modification of the Lax Friedrichs flux splitting to eliminate the viscosity term, thus reaching the C-property. In a 2d context, this substitution is performed by Bradford [START_REF] Bradford | Finite-volume model for shallow-water flooding of arbitrary topography[END_REF] for the evaluation of surface integrals in some correction terms. Meanwhile, an alternate formulation of NSW equations, using the free surface elevation ξ above the still water depth, was introduced by Rogers et al [START_REF] Rogers | Adaptive Q-tree Godunov-type scheme for shallow water equations[END_REF] with the objective of obtaining telescoping effects between fluxes and source contributions employing a Roe solver. In a similar way, Kurganov and Levy [START_REF] Kurganov | Central-upwind schemes for the saint-venant system[END_REF] and Russo [START_REF] Russo | Central schemes for conservation laws with application to shallow water equations[END_REF] advanced the idea to directly replace h by η -z in the set of primitive equations.

Nowadays, in view of the efficiency of the approach, and based on these pioneering works, many improvements and derivations have been proposed. Kurganov and Petrova [START_REF] Kurganov | A second-order well-balanced positivity preserving centralupwind scheme for the saint-venant system[END_REF] improved the original scheme [START_REF] Kurganov | Central-upwind schemes for the saint-venant system[END_REF] enforcing simultaneously well balancing and robustness, and Bryson et al [START_REF] Bryson | Well-balanced positivity preserving central-upwind scheme on triangular grids for the Saint-Venant system[END_REF] worked on an extension to unstructured triangulations. We also refer to [START_REF] Singh | Two-dimensional numerical modeling of dam-break flows over natural terrain using a central explicit scheme[END_REF] for another using on 2d cartesian grids, and [START_REF] Brodtkorb | Efficient shallow water simulations on GPUs : Implementation, visualization, verification and validation[END_REF] for an implementation on GPU and additional highlight of performances. In this work, we seek a discretization of SWE from such a point of view, working with the so called pre-balanced formulation (see [START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF][START_REF] Liang | Adaptive quadtree simulation of shallow flows with wetdry fronts over complex topography[END_REF]):

V t + ∇ • H(V, z) = S(V, z) + F (V, z) , (6)
where

V =   η q x q y   , H(V, 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)   . (7)
The bed slope source term is given by:

S(V, z) =   0 -gηz x -gηz y   , (8)
while the friction term is still defined as: F (V, z) = E(U ), using the relation h = η -z. This system is still hyperbolic, and has the same eigenstructure as [START_REF] Anastasiou | Solution of the 2d shallow water equations using the finite volume method on unstructured triangular meshes[END_REF]. It should be noticed that the convenient aspects of this formulation have been successfully exploited in several well-balanced and robust FVM approaches [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][START_REF] Hou | A stable 2d unstructured shallow flow model for simulations of wetting and drying over rough terrains[END_REF][START_REF] Wang | A 2d shallow flow model for practical dam-break simulations[END_REF].

As for the complex issue of the numerical treatment of vanishing water heights on triangulations, lots of numerical methods are able to handle efficiently occurrence of dry areas (see for instance [START_REF] Bokhove | Flooding and drying in finite-element discretizations of shallow-water equations. part 2: Two dimensions[END_REF][START_REF] Bollermann | Finite volume evolution Galerkin methods for the shallow water equations with dry beds[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] Nikolos | An unstructured node-centered finite volume scheme for shallow water flows with wet/dry fronts over complex topography[END_REF][START_REF] Özkan-Haller | A Fourier-Chebyshev collocation method for the shallow water equations including shoreline runup[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF] and references therein). This being so, among the quantity of available works, the problem of a suitable approach for the bed friction source term may deserve investigations yet. The main reason stands in the fact that divisions by h may be involved, and possibly generate instabilities in low water depths areas. As a further restriction, the use of unstructured triangulations makes difficult the implementation of direct extensions of 1d approaches, which is generally not the case on cartesian grids.

To circumvent this difficulty, in some appropriate circumstances, linear terms [START_REF] Anastasiou | Solution of the 2d shallow water equations using the finite volume method on unstructured triangular meshes[END_REF][START_REF] Zhou | The surface gradient method for the treatment of source terms in the shallow-water equations[END_REF] or simplified friction laws [START_REF] Bryson | Well-balanced positivity preserving central-upwind scheme on triangular grids for the Saint-Venant system[END_REF][START_REF] Kurganov | A second-order well-balanced positivity preserving centralupwind scheme for the saint-venant system[END_REF] may appear relevant as a first approximation. However, this solution cannot be considered as totally satisfactory since it does not allow to describe a large range of realistic situations. More severe laws may be accounted for through straightforward pointwize methods [START_REF] Hubbard | A 2d numerical model of wave run-up and overtopping[END_REF][START_REF] Vukovic | Eno and weno schemes with the exact conservation property for one-dimensional shallow water equations[END_REF][START_REF] Ying | Improved hll scheme for 1d dam-break flows over complex topography[END_REF], provided again that the resistance parameters are moderate enough. Such an explicit method is suggested in [START_REF] Nikolos | An unstructured node-centered finite volume scheme for shallow water flows with wet/dry fronts over complex topography[END_REF] in a 2d node-centred unstructured FVM. In more general contexts, when attempting to reach better level of stability, semi implicit or fully implicit treatments are often employed [START_REF] Bristeau | Boundary conditions for the shallow water equations solved by kinetic schemes[END_REF][START_REF] Burguete | Efficient construction of high-resolution tvd conservative schemes for equations with source terms : application to shallow water flows[END_REF][START_REF] Delis | Numerical flood simulation by depth averaged free surface flow models[END_REF]. In the 2d unstructured context, an efficient upwind semi implicit method can be found in [START_REF] Cea | Unstructured finite volume discretization of bed friction and convective flux in solute transport models linked to the shallow water equations[END_REF], with a supplementary effort concerning the balance between source and friction terms to exactly restore uniform channel flows equilibriums.

However, calculations issuing from these methods are generally not exempted from divisions by h either, and underlying stability issues still holds. As stated in [START_REF] Bradford | Finite-volume model for shallow-water flooding of arbitrary topography[END_REF], numerical models turn out to be very sensitive to the calibration of the tolerance value used to define the wet/dry transitions when friction terms are included. This threshold has sometimes to be considerably increased to avoid the perturbation mechanisms that threaten to be activated in very shallow waters, which becomes a dilemma if an accurate description of moving interfaces is simultaneously desired.

Although being generally overshadowed by other concerns in many works, several propositions driven by this problematic started to emerge since the last decade. Motivated by the necessity to have some control on these potential instabilities, Liang and Marche [START_REF] Liang | Adaptive quadtree simulation of shallow flows with wetdry fronts over complex topography[END_REF] proposed to add an additional a posteriori limiting stage to prevent an eventual reversing of the flow, which is evidence that the resistance effects have been locally overestimated. This method is successfully reused in [START_REF] Wang | A 2d shallow flow model for practical dam-break simulations[END_REF] on cartesian grids, and more recently in the frame of unstructured grids [START_REF] Hou | A robust well-balanced model on unstructured grids for shallow water flows with wetting and drying over complex topography[END_REF][START_REF] Hou | A stable 2d unstructured shallow flow model for simulations of wetting and drying over rough terrains[END_REF]. Note that similar stability arguments were also advanced by Murillo et al in the context of triangular meshes [START_REF] Murillo | The influence of source terms on stability, accuracy and conservation in two-dimensional shallow flow simulation using triangular finite volumes[END_REF]. The method ensures that the sign of the discharge remains unchanged under a suitable time step condition. Thereafter, some improvements have been proposed by the authors in order to relax this time constraint [START_REF] Murillo | Time step restrictions for well-balanced shallow water solutions in non-zero velocity steady states[END_REF], with the use of an upwind explicit well-balanced method, combined with an implicit discretization of the friction source term, specifically employed in critical areas to prevent from instabilities. Based on these ideas, another strategy consists in directly limiting the friction force to ensure relevant and stable computations [START_REF] Burguete | Friction term discretization and limitation to preserve stability and conservation in the 1d shallow-water model: Application to unsteady irrigation and river flow[END_REF].

Another way to proceed is borrowed from the recent headways made in the context of Asymptotic Preserving schemes [START_REF] Berthon | Asymptotic preserving hll schemes[END_REF]. The key idea consists of accounting for the source through a modified Riemann solver. The strategy is successfully employed in 1d for the NSW equations [START_REF] Berthon | An efficient scheme on wet/dry transitions for shallow water equations with friction[END_REF], leading to a stable, robust and well-balanced second order scheme. It should be stressed that similar ideas can be encountered in [START_REF] Murillo | Augmented versions of the hll and hllc riemann solvers including source terms in one and two dimensions for shallow flow applications[END_REF], where the authors include the bed slope and bottom friction terms within an augmented version of the HLL solver.

The goal of this paper is to propose a FVM which satisfies all the points mentioned above. More precisely, following the general strategy described in [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF], one of the main the objectives is to extend the well-balanced and robustness properties of the approach developed in [START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF]. Still starting from a 1d discretization of the pre-balanced equations, the construction of the 2d scheme is slightly modified, such that the approximate vector variable update can be seen as a convex combination of 1d schemes. Taking advantage of this construction, the inclusion of friction is also based on a 1d approach, inspired by the ideas used by Berthon et al [START_REF] Berthon | Asymptotic preserving hll schemes[END_REF]. The resistance effects are accounted for through a suitable adaptation of the intermediate states involved in the HLL approximate Riemann solver. An application can be found in the FVM proposed in [START_REF] Berthon | An efficient scheme on wet/dry transitions for shallow water equations with friction[END_REF] for the 1d SWE, providing a robust scheme, and ensuring stable computations for low values of the water height. Similarly, our work is based on the consideration of an alternative formulation of the HLL solver, with the objective of extending the benefits of this method to the 2d unstructured case.

The structure of the present work is the following: first, we propose to detail the construction of the 2d pre-balanced approach in the absence of friction. The main settings and results of the original 1d scheme are restated in Part 1.1, and the extension to the unstructured case is subsequently proposed. The second section is devoted to the inclusion of the bed friction source term. At least, an experimental part will involve the description of some appropriate benchmark tests. Amongst others, they are expected to validate the ability in handling appropriately bathymetry and friction source terms, notably in the vicinity of dry cells. Numerical validations including accuracy and convergence rate analysis are also performed, highlighting the benefits brought by the MUSCL reconstruction [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF].

1 Pre Balanced scheme and 2d unstructured derivation

The original approach

We briefly recall the possible approach of the 1d pre-balanced system [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF], recently proposed in [START_REF] Wang | A 2d shallow flow model for practical dam-break simulations[END_REF]. Let's consider a regular partition of the computational domain Ω with a space step ∆x. In the present Godunov-type Finite Volume framework, using a Euler scheme in time for the sake of simplicity, and denoting V i the flow variable approximation on the cell i, the flow update from time t n to time t n+1 = t n + ∆t is governed by the classical relation:

V n+1 i = V n i - ∆t ∆x H i+ 1 2 -H i-1 2 + ∆tS i , (9)
where H i-1 2 are H i+ 1 2 are the numerical fluxes and S i stands for a discretization of the bed slope. The evaluation of the exchanging fluxes H i-1 2 and H i+ 1 2 is based on the resolution of a local Riemann problem involving reconstructed states at each side of the interfaces i -1/2 and i + 1/2. They take the form:

H i-1 2 = H(V - i-1 2 , V + i-1 2 , zi-1 2 , zi-1 2 ) , H i+ 1 2 = H(V - i+ 1 2 , V + i+ 1 2 , zi+ 1 2 , zi+ 1 2 ) , (10)
where a double dependency with z is introduced, to account for the presence of the topography in the convective fluxes. The evaluation of these intermediate states is explained in the following lines.

Borrowing the pioneering ideas of [2], the reconstruction process starts by the introduction of a new interface value for the bed elevation:

zi+ 1 2 = max(z i , z i+1 ) , (11)
and the following positive reconstruction of the water depth component (see [START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF]):

h - i+ 1 2 = max(0, η i -zi+ 1 2 ), h + i+1/2 = max(0, η i+1 -zi+ 1 2 ) . (12)
From this, new interface values for the free surface and discharge are defined by (see Fig. 2a):

η - i+ 1 2 = h - i+ 1 2 + zi+ 1 2 , q - i+1/2 = h - i+1/2 u i . η + i+ 1 2 = h + i+ 1 2 + zi+ 1 2 , q + i+1/2 = h + i+1/2 u i+1 . (13)
Remark 1.1. In the evaluation of the discharge (13), we implicitly assumed that we have h ± > 0 to extract the velocity from the initial discharge. The case h ± = 0 (or numerically h ± < ǫ) is classically handled by setting the corresponding values of q ± (and u ± if needed) to 0. Cells for which h is falling under this tolerance will be considered as dry.

We follow a similar procedure at the interface i -1/2. As stated in [START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF], these new values can be directly injected in the numerical fluxes, leading to a well-balanced scheme that also ensures non-negative water depths, provided the use of an appropriate Riemann solver. However, this strategy requires the introduction of additional source terms in order to ensure the preservation of motionless steady states [START_REF] Benkhaldoun | A new finite volume method for flux-gradient and source-term balancing in shallow water equations[END_REF] in wet/dry contexts. To avoid this, one can consider the following parameter:

∆ i+ 1 2 = max(0, zi+ 1 2 -η i ) , (14)
specifically introduced to account for such particular configurations (see Fig. 2b). Finally, the new intermediate states for the topography and free surface are obtained after subtraction from the original values:

zi+ 1 2 = zi+ 1 2 -∆ i+ 1 2 , η - i+ 1 2 = h - i+ 1 2 + zi+ 1 2 -∆ i+ 1 2 , η + i+ 1 2 = h + i+ 1 2 + zi+ 1 2 -∆ i+ 1 2 . ( 15
)
z i h i η i z i+1 h i+1 η i+1 i i + 1 zi+ 1 2 h - i+ 1 2 η - i+ 1 2 zi+ 1 2 h + i+ 1 2 η + i+ 1 2 i i + 1 (a)
Case of a wet/wet interface.

z i h i η i η i+1 = z i+1 i i + 1 ∆ i+ 1 2 i i + 1 η - i+ 1 2 η + i+ 1 2 = zi+ 1 2 (b)
Case of a wet/dry interface. Obviously, this additional step [START_REF] Bouchut | Frontal geostrophic adjustment and nonlinear wave phenomena in one-dimensional rotating shallow water. part 2. high-resolution numerical simulations[END_REF] does not affect the flow variables in wet bed applications, since in this case ∆ i+ 1 2 = 0 and zi+ 1 2 = zi+ 1 2 . It is shown in [START_REF] Wang | A 2d shallow flow model for practical dam-break simulations[END_REF] that combined with the following centred discretization of the source:

S i =   0 -g ηi zi+ 1 2 -zi-1 2 ∆x   , (16)
where ηi =

η + i-1 2 + η - i+ 1 2 2
, this reconstruction extends the C-property to configurations involving wet/dry interfaces.

In this work, the evaluation of the numerical fluxes will be performed with the HLL Riemann solver. In addition of ensuring the continuity with the investigations of Section 2, this choice appears to be a good compromise between simplicity and efficiency. Let us also mention that the resulting flux function H is Lipschitz continuous and satisfies consistency and conservation properties. We also have the following robustness result:

Proposition 1. We consider the scheme ( 9), together with the reconstruction steps [START_REF] Berthon | Asymptotic preserving hll schemes[END_REF][START_REF] Bokhove | Flooding and drying in finite-element discretizations of shallow-water equations. part 2: Two dimensions[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF][START_REF] Bouchut | Frontal geostrophic adjustment and nonlinear wave phenomena in one-dimensional rotating shallow water. part 2. high-resolution numerical simulations[END_REF]. Assume that the numerical fluxes are computed with the HLL solver, and that h n i-1 , h n i and h n i+1 are positive. Then, under the following CFL condition:

s ∆t ∆x ≤ τ CF L 2 , ( 17
)
where s = max i∈Z (|u i ± c i |) and τ CF L ≤ 1, we have h n+1 i ≥ 0.
Proof. We introduce the following notations for the HLL wave speeds (see [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF]):

s - i+ 1 2 = min(u i -c i , u i+1 -c i+1 ) , s + i+ 1 2 = max(u i + c i , u i+1 + c i+1 ) , ∆s i+ 1 2 = s + i+ 1 2 -s - i+ 1 2
, and set:

s - i+ 1 2 = min(0, s - i+ 1 2 ) , s + i+ 1 2 = max(0, s + i+ 1 2 ) , ∆s i+ 1 2 = s + i+ 1 2 -s - i+ 1 2 .
The HLL flux hence writes:

H η i+ 1 2 = s + i+ 1 2 q - i+ 1 2 -s - i+ 1 2 q + i+ 1 2 + s + i+ 1 2 s - i+ 1 2 (η + i+ 1 2 -η - i+ 1 2 ) ∆s i+ 1 2 .
According to [START_REF] Bouchut | Frontal geostrophic adjustment and nonlinear wave phenomena in one-dimensional rotating shallow water. part 2. high-resolution numerical simulations[END_REF], we have:

η + i+ 1 2 -η - i+ 1 2 = h + i+ 1 2 -h - i+ 1 2
, leading to:

H η i+ 1 2 = s + i+ 1 2 h - i+ 1 2 (u i -s - i+ 1 2 ) ∆s i+ 1 2 + s - i+ 1 2 h + i+ 1 2 (s + i+ 1 2 -u i+1 ) ∆s i+ 1 2 .
With similar considerations at the other interface, and substracting z i on both sides of the free surface issuing from ( 9), we have:

h n+1 i =h n i 1 - ∆t ∆x h - i+ 1 2 h n i s + i+ 1 2 (u i -s - i+ 1 2 ) ∆s i+ 1 2 - h + i-1 2 h n i s - i-1 2 (s + i-1 2 -u i ) ∆s i-1 2 + h n i ∆t ∆x h - i-1 2 h n i s + i-1 2 (u i-1 -s - i-1 2 ) ∆s i-1 2 - h + i+ 1 2 h n i s - i+ 1 2 (s + i+ 1 2 -u i+1 ) ∆s i+ 1 2 . ( 18
)
Finally, arguing that

h - i+ 1 2 h n i , h + i-1 2 h n i ∈ [0, 1],
we easily obtain:

0 ≤ h - i+ 1 2 h n i s + i+ 1 2 (u i -s - i+ 1 2 ) ∆s i+ 1 2 ≤ s + i+ 1 2 ≤ s , 0 ≤ - h + i-1 2 h n i s - i-1 2 (s + i-1 2 -u i ) ∆s i-1 2 ≤ -s - i-1 2 ≤ s .
As a consequence, under [START_REF] Bristeau | Boundary conditions for the shallow water equations solved by kinetic schemes[END_REF], the second member in equality ( 18) comes down to the sum of two positive values. The proof is complete. Note that a relaxed condition can be obtained with the use of Global Lax Friedrichs fluxes (see [START_REF] Xing | Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations[END_REF]), simply taking τ CF L as CFL estimate at the right side of inequality [START_REF] Bristeau | Boundary conditions for the shallow water equations solved by kinetic schemes[END_REF].

Extension to the 2D unstructured case

We start by a brief introduction of some settings and notations. Computations are run on a dual mesh issuing from a triangulation T . We will denote the cells by C i , and their area by

|C i |.
In what follows, K i stands for the set of subscripts j for which C j is adjacent to C i . We also set Γ ij , the boundary interface between two adjacent cells C i and C j , ℓ ij its length, and n ij the unit normal to Γ ij pointing to C j (Fig. 3). Still considering a Finite Volume framework, we suggest to advance in time by the mean of the usual formula:

T ij T k M i M j n ij Γ ij C i C j
V n+1 i = V n i - ∆t |C i | j∈Ki ℓ ij H ij + ∆tS i , (19)
where H ij , j ∈ K i are the exchanging fluxes between the cell C i and its neighbours C j , and S i denotes the bottom source term discretization. Following the 1d procedure, we write :

H ij = H(V * ij , V * ji , zij , zij , n ij ) , (20)
and aim at proposing a non-negative reconstruction for the Riemann states that allow the restoration of the well-balanced property in wet/dry configurations. Therefore, for a given interface Γ ij , we define a single valued term for the topography:

zij = max(z i , z j ) , (21)
and the water depth reconstruction is given by:

h * ij = max(0, η i -zij ), h * ji = max(0, η j -zji ) . (22)
Then, gathering steps [START_REF] Bollermann | Finite volume evolution Galerkin methods for the shallow water equations with dry beds[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF][START_REF] Bouchut | Frontal geostrophic adjustment and nonlinear wave phenomena in one-dimensional rotating shallow water. part 2. high-resolution numerical simulations[END_REF], we naturally set:

∆ ij = max(0, zij -η i ) , zij = zij -∆ ij , η * ij = h * ij + zij -∆ ij , q * ij = h * ij u i , η * ji = h * ji + zij -∆ ij , q * ji = h * ji u j , (23)
with the necessary caution mentioned in Remark 1.1 for the treatment of vanishing water depths.

To complete the numerical scheme, we introduce the following discretization of the bottom source term:

S i = j∈Ki ℓ ij |C i | 0 g ηij (z i -zij ) n ij , (24)
where

ηij = η i + η * ji 2 . ( 25
)
The motivations for this choice of ηij are twofold. The first point is linked to the general problem of consistency, meaning that the contribution S i actually has to provide a correct approximation of the physical source term issuing from the continuous equations. The second, more specific, concerns the structure of the scheme. Indeed, from such a point of view, this particular choice for ηij has important consequences, since the numerical approach can actually be seen as convex combinations of 1d schemes of the form [START_REF] Berthon | A positive preserving high order VFRoe scheme for shallow water equations: a class of relaxation schemes[END_REF]. To this end, let's split each vertex cell C i in a set of triangles (T ij ) j∈Ki defined by the edge Γ ij and the node M i (Fig. 3). In the following lines, we show how this subgrid can be exploited to obtain a proof of the robustness and C-properties in a considerably simple way, with a straightforward stability criteria. More precisely, we have the following result :

Proposition 2. We consider the scheme [START_REF] Brufau | A numerical model for the flooding and drying of irregular domains[END_REF], together with the reconstruction steps [START_REF] Burguete | Efficient construction of high-resolution tvd conservative schemes for equations with source terms : application to shallow water flows[END_REF][START_REF] Burguete | Friction term discretization and limitation to preserve stability and conservation in the 1d shallow-water model: Application to unsteady irrigation and river flow[END_REF][START_REF] Camarri | A low-diffusion muscl scheme for les on unstructured grids[END_REF], and assume that the numerical fluxes are computed with the HLL solver. Then, the C-property is satisfied. Moreover, if we suppose h n i ≥ 0 and h n j ≥ 0 , ∀ j ∈ K i , then, the water height positivity is preserved under the following CFL condition:

s ij ∆t ℓ ij |T ij | ≤ τ CF L 2 , ∀C i , ∀j ∈ K i , (26)
where

s ij = max (|u i . n ij ± c i | , |u j . n ij ± c j |) and τ CF L ≤ 1.
Proof. The key idea consist in rewriting the scheme [START_REF] Brufau | A numerical model for the flooding and drying of irregular domains[END_REF] as a convex combination of 1d schemes.

In the spirit of [START_REF] Berthon | Efficient well balanced hydrostatic upwind schemes for shallow water equations[END_REF], the current scheme is recast under the following form:

V n+1 i = j∈Ki |T ij | |C i | V n+1 ij , (27)
where each convex component is defined as:

V n+1 ij = V n i - ∆t δ ij H ij -H c ij + ∆tS ij , (28)
with

S ij = 1 δ ij 0 -g ηij (z ij -z i ) n ij , (29)
and where we have set

δ ij = |T ij |/ℓ ij and H c ij = H(V n i , z i ). n ij .
Involving the solver's consistency and the steps (21 -22 -23), one can easily verify that H c ij corresponds to the numerical flux in the direction n ij evaluated with V n i as "left" and "right" Riemann states. It follows that relation [START_REF] Cournède | A Positive MUSCL Scheme for Triangulations[END_REF] is nothing but the 1d scheme (9) applied to the states V n i , V n i , V n j (Fig. 4) in the direction n ij , and with a space step δ ij . Consequently, according to [START_REF] Chippada | A Godunov-type finite volume method for the system of shallow water equations[END_REF], the preservation of the water height positivity is directly inherited from the 1d numerical model (Proposition 1), leading to the criteria [START_REF] Chanson | Analytical solution of dam break wave with flow resistance. application to tsunami surges[END_REF]. The C-property can also be straightforwardly deduced from this construction, arguing that formula [START_REF] Chippada | A Godunov-type finite volume method for the system of shallow water equations[END_REF] preserves the motionless steady states as soon as (28) does.

V - i-1 2 V + i-1 2 V - i+ 1 2 V + i+ 1 2 i -1 i i + 1 ∆x V i V i V * ij V * ji C i C i C j δ ij
Remark 1.2. It should be noticed that this scheme is very close to the one proposed in [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF]. Actually, the only modification consists in the choice of η in the source term, as we made the choice ηij = (η * ji +η * ij )/2 in the previous paper. This slight alternative appears to be advantageous, as it allows to relax the time constraint in a simple way, and directly ensures de C-property.

Remark 1.3. In this work, we suggest to broach the problem of boundary conditions through the consideration of ghost cells at each boundary interface. The value of the solution is enforced on the set of fictive cells, according to the problem under consideration. This strategy is chosen for its simplicity and efficiency in handling Neumann and Dirichlet conditions, which are principally involved in the numerical simulations of the last section. Note that solid wall conditions can also be prescribed, combining a Neumann condition on h and a Dirichlet condition on q.

Remark 1.4. As far as the time discretization is concerned, it's worthwhile to recall that more accurate time-marching algorithms can be used for the advance in time. For the numerical validations, a second order explicit Runge Kutta method will be preferred to the Euler scheme currently used to illustrate the scheme's construction.

Friction scheme 2.1 The 1d case

Now, let us take advantage of formulation [START_REF] Chippada | A Godunov-type finite volume method for the system of shallow water equations[END_REF][START_REF] Cournède | A Positive MUSCL Scheme for Triangulations[END_REF] to include the friction source term. To address such an issue, we propose to employ a 1d friction scheme for the evaluation of each convex component [START_REF] Cournède | A Positive MUSCL Scheme for Triangulations[END_REF]. Let's first assume that the friction law has the following form:

E(U ) = -σ |q| h γ 0 q , (30)
with σ, γ two positive constants. Note that Manning and Darcy friction laws (4a,4b) falls within such a formalism, with σ = n 2 , γ = 10/3 and σ = d/8, γ = 2 respectively. Following the general ideas introduced in [START_REF] Berthon | Asymptotic preserving hll schemes[END_REF] or [START_REF] Berthon | A free streaming contact preserving scheme for the m1 model[END_REF] for such class of source terms, the method relies on the use of a modified HLL scheme. Denoting a ± the minimum and maximum velocity waves involved in the approximate Riemann solver, the inclusion of friction is accounted for through the following corrected intermediate states:

ṼR ( x t , V L , V R ) =          V L if x t ≤ a -, αV * + (1 -α)R(V L ) if min(0, a -) < x t < min(0, a + ) , αV * + (1 -α)R(V R ) if max(0, a -) < x t < max(0, a + ) , V R if x t ≥ a + . ( 31
)
where V * is the classical HLL intermediate state:

V * ( x t , V L , V R ) = a + V R -a -V L -H(V R ) -H(V L ) a + -a - , (32)
and we use the simplified notation:

R(V ) = R η (V ) R q (V ) = η * q -|q|q . ( 33
)
Remark 2.1. As concerns the expression of R given just above, the choice of the first component appears as natural, considering that, at the continuous level, the friction term does not impact the mass equation. Thus, if we denote V R the usual HLL Riemann solver, we have V η R = Ṽ η R , and we consequently have to focus on the evolution of the discharge only. Considering now the discharge, we first remark that the corresponding component of the friction source term can be recast under the form:

E q (U ) = ρ(h) R q (V ) -q , with ρ(h) = σ h γ . ( 34
)
Historically, in the pioneering works mentioned above, modified HLL solvers are precisely introduced to handle source terms under the general form [START_REF] Dutykh | The VOLNA code for the numerical modelling of tsunami waves: generation, propagation and inundation[END_REF]. In several physical situations of interest indeed, the presence of a convex combination linking V * to R (appearing here in the intermediate waves [START_REF] Delis | Numerical flood simulation by depth averaged free surface flow models[END_REF]) allows both to recover the classical solver in absence of source terms, and supplies a technical support for the restoration of particular limit regimes at the discrete level (corresponding to the situation where α tends to zero). Obviously, it requires a consistent calibration of α, offering relevant approximations in both limits. In this work, and based on these objectives, we actually show that this modification allows a proper evaluation of the friction source term, able to cover very low resistance values and extremely rough coefficients.

In a 1d context, we can straightforwardly reproduce the analysis performed in [START_REF] Berthon | Asymptotic preserving hll schemes[END_REF] for this particular solver. We skip the details of the computations, to write the following time evolved value for q = hu:

q n+1 i = q n i - ∆t ∆x α i+ 1 2 H q i+ 1 2 -α i-1 2 H q i-1 2 + ∆t ∆x (1 -α i-1 2 )F +,q i-1 2 + (1 -α i+ 1 2 )F -,q i+ 1 2 . ( 35
)
The consideration of the modified solver hence induces the presence of weighted fluxes and generates a numerical source term, designed to approximate the friction forces. In the above formula, the fluxes H i± 1 2 are provided by the usual HLL scheme, applied to the reconstructed states [START_REF] Berthon | Asymptotic preserving hll schemes[END_REF][START_REF] Bokhove | Flooding and drying in finite-element discretizations of shallow-water equations. part 2: Two dimensions[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF][START_REF] Bouchut | Frontal geostrophic adjustment and nonlinear wave phenomena in one-dimensional rotating shallow water. part 2. high-resolution numerical simulations[END_REF] to guarantee robustness and well-balancing properties. The additional friction terms read as follows:

F +,q i-1/2 = max(0, a - i-1/2 ) q i-1 -R q (V i ) + max(0, a + i-1/2 ) R q (V i ) -q i + H q (V i ) , F -,q i+1/2 = min(0, a - i+1/2 ) q i -R q (V i ) + min(0, a + i+1/2 ) R q (V i ) -q i+1 -H q (V i ) . (36)
We come now to the evaluation of the quantities α i± 1 2 , which are the discrete equivalents of the parameter α ∈ [0, 1] introduced in [START_REF] Delis | Numerical flood simulation by depth averaged free surface flow models[END_REF]. As stated in Remark 2.1, this coefficient mainly controls friction and has to be appropriately defined. One of the first requirements is driven by the fact that we ask the modified solver to degenerate toward the usual HLL solver when the resistance effects tends to zero. On the other hand, the case of dominant friction terms will correspond to the situation where α tends to zero. Another fundamental requirement is related to consistency aspects, since we actually have to ensure that the contribution:

F i := 1 -α i-1 2 ∆x F +,q i-1 2 + 1 -α i+ 1 2 ∆x F -,q i+ 1 2 (37)
furnishes a correct approximation of the corresponding term E q (U ) of the continuous equations. Lastly, keeping in mind the willingness to prevent from computational difficulties in situations of vanishing water depths, the expression of α should ideally remain relevant whenever h tends to (or even reaches) zero. To summarize, the criteria that have to be satisfied are:

• α ∈ [0, 1] .
• lim σ→∞ α = 0 .

• α = 1 in absence of friction.

• α is defined on dry cells.

• Consistency with the friction source term of the continuous equations. Among the possible choices satisfying these properties, one can consider (see also [START_REF] Berthon | An efficient scheme on wet/dry transitions for shallow water equations with friction[END_REF]) :

α = h γ a + -a - h γ a + -a -+ σ∆x , where h = h -+ h + 2 , (38)
from which are extracted the coefficients α i± 1 2 involved in [START_REF] Ern | A well-balanced Runge-Kutta discontinuous Galerkin method for the shallow-water equations with flooding and drying[END_REF], using the exterior and interior values at the interfaces i ± 1/2. Remark 2.2. Thus, with this choice, it's worth mentioning that the friction source term discretization F i proposed in [START_REF] Gallouët | Some approximate godunov schemes to compute shallow-water equations with topography[END_REF] makes it possible to avoid divisions by h, which traditionally threaten the stability of numerical models dealing with laws based on [START_REF] Delestre | A numerical scheme for a viscous shallow water model with friction[END_REF]. The consistency of F i with the continuous equations will be established in the 2d case and is not given here. We refer to [START_REF] Berthon | Asymptotic preserving hll schemes[END_REF] for consistency results in the 1d case.

The 2d case

Now, let us take advantage of the formulation [START_REF] Chippada | A Godunov-type finite volume method for the system of shallow water equations[END_REF][START_REF] Cournède | A Positive MUSCL Scheme for Triangulations[END_REF] to include the friction source term. We propose to carry out the update of each convex component by the mean of [START_REF] Ern | A well-balanced Runge-Kutta discontinuous Galerkin method for the shallow-water equations with flooding and drying[END_REF]. Based on the 1d scheme, for any edge Γ ij , we modify each convex component of the 1d-like scheme in [START_REF] Cournède | A Positive MUSCL Scheme for Triangulations[END_REF] as follows:

hu n+1 ij = hu n i -∆t ℓ ij |T ij | α ij H hu ij -α c ij H c,hu ij + ∆t ℓ ij |T ij | (1 -α c ij )F c,hu ij + (1 -α ij )F hu ij , hv n+1 ij = hv n i -∆t ℓ ij |T ij | α ij H hv ij -α c ij H c,hv ij + ∆t ℓ ij |T ij | (1 -α c ij )F c,hv ij + (1 -α ij )F hv ij , (39)
where the friction parameters are straightly derived from those given in the 1d formulation. In the sequel, we denote R q (V ) = t R hu (V ), R hv (V ) = qq q, and employ similar notations to gather the discharge components of the friction terms F q ij , F c,q ij involved in [START_REF] Gosse | A well-balanced flux-vector splitting scheme designed for hyperbolic systems of conservation laws with source terms[END_REF]. Let us also employ the notations a c,± ij , a ± ij for the velocity waves involved in the computation of H c ij and H ij respectively, and set hij = (h i + h j )/2. Thus, in the reference associated with n ij , each contribution relative to friction is evaluated that way:

α c ij = h γ i a c,+ ij -a c,- ij h γ i a c,+ ij -a c,- ij + σδ ij , α ij = hγ ij a + ij -a - ij hγ ij a + ij -a - ij + σδ ij , F c,q k = max(0, a c,- ij ) q i -R q (V i ) + max(0, a c,+ ij ) R q (V i ) -q i + H q (V i , z i ). n ij , F q k = min(0, a - ij ) q i -R q (V i ) + min(0, a + ij ) R q (V i ) -q j -H q (V i , z i ). n ij , (40)
which are nothing but the equivalents of ( 38) and [START_REF] Gallardo | On a well-balanced high-order finite volume scheme for shallow water equations with topography and dry areas[END_REF] in the current 2d context. Now, keeping in mind that the evolution of the free surface is still governed by the original scheme, and including the bed slope source term introduced in (24), formula [START_REF] Chippada | A Godunov-type finite volume method for the system of shallow water equations[END_REF][START_REF] Cournède | A Positive MUSCL Scheme for Triangulations[END_REF] gives, integrating modifications (39) :

η n+1 i = η n i - ∆t |C i | j∈Ki ℓ ij H η ij , hu n+1 i = hu n i - ∆t |C i | j∈Ki ℓ ij α ij H hu ij -α c ij H c,hu ij + ∆t |C i | j∈Ki ℓ ij F hu ij + ∆tS hu i , hv n+1 i = hv n i - ∆t |C i | j∈Ki ℓ ij α ij H hv ij -α c ij H c,hv ij + ∆t |C i | j∈Ki ℓ ij F hv ij + ∆tS hv i , (41)
where

F q ij = t (F hu ij , F hv ij ) = (1 -α c ij )F c,q ij + (1 -α ij )F q ij . (42)
It is useful to note that, invoking the relation H c ij = H(V i , z i ). n ij , the discharge evolution in the previous set of equations can equally be written as :

hu n+1 i = hu n i - ∆t |C i | j∈Ki ℓ ij α ij H hu ij -α ij H c,hu ij + ∆t |C i | j∈Ki ℓ ij E hu ij + ∆tS hu i , hv n+1 i = hv n i - ∆t |C i | j∈Ki ℓ ij α ij H hv ij -α ij H c,hv ij + ∆t |C i | j∈Ki ℓ ij E hv ij + ∆tS hv i , (43)
where

E q ij = t (E hu ij , E hv ij ) = (1 -α c ij )E c,q ij + (1 -α ij )E q ij , E c,q ij = max(0, a c,- ij ) q i -R q (V i ) + max(0, a c,+ ij ) R q (V i ) -q i , E q ij = min(0, a - ij ) q i -R q (V i ) + min(0, a + ij ) R q (V i ) -q j . (44)
Thanks to the use of convex combinations, the consistency results issuing from the 1d analysis ensure a relevant estimation of the source term appearing in the continuous equations. More precisely, we have the following result :

Proposition 3. Assume that the numerical fluxes are computed using the Rusanov fluxes. In other words, we define the maximum and minimum characteristic speeds in the HLL Riemann solver as :

a ± ij = ±s ij , s ij = max |u i . n ij | + gh i , |u j . n ij | + gh j . (45)
Then, the contribution ∆t

|C i | j∈Ki ℓ ij E q ij is
consistent with the friction source term present in the momentum equations of (6).

Proof. We first notice that in these conditions, [START_REF] Hubbard | A 2d numerical model of wave run-up and overtopping[END_REF] writes :

E q ij = s c ij (1 -α c ij )(R q (V i ) -q i ) + s ij (1 -α ij )(R q (V i ) -q i ) .
We then write :

ℓ ij |C i | (1 -α c ij ) = |T ij | |C i | σ 2s c ij h γ i + σδ ij , ℓ ij |C i | (1 -α ij ) = |T ij | |C i | σ 2s ij hγ ij + σδ ij ,
and study the behaviour of the source term j∈Ki ℓ ij |C i | E q ij as d i = diam(C i ) tends to zero. We can easily check that we have :

lim di→0 ℓ ij |C i | s c ij (1 -α c ij )(R q (V i ) -q i ) = lim di→0 ℓ ij |C i | s ij (1 -α ij )(R q (V i ) -q i ) = |T ij | |C i | σ 2h γ i (R q (V i ) -q i ) .
Consequently :

lim di→0 1 |C i | j∈Ki ℓ ij E q ij = j∈Ki |T ij | |C i | σ h γ i (R q (V i ) -q i ) . = E q (U i ) .
Remark 2.3. Recalling that no additional time constraint has been introduced, that the evolution of the total free surface is not modified, and that the friction terms vanish in absence of motion, robustness and C-property results still holds. Thus, one of the advantages of the current method stands in the fact that the CFL [START_REF] Chanson | Analytical solution of dam break wave with flow resistance. application to tsunami surges[END_REF] does not need to be modified in the presence of friction.

In practice, we take τ CF L = 0.9 in our numerical validations. Again, note that the friction terms evaluation [START_REF] Gosse | A well-balanced scheme using non-conservative products designed for hyperbolic systems of conservation laws with source terms[END_REF] is not problematic when the water depth tends to zero. In some sense, this feature hence allows a uniform calibration of the threshold ǫ used to discriminate dry cells, independent of whether friction terms are included or not. This value will be set to 1.e -6 in our numerical validations.

Improving space accuracy

From a general point of view, the diffusive losses entailed by the consideration of first order schemes generally require the implementation of numerical techniques devoted to increase the space accuracy. This problem arises notably in the context of FVM with the use of cell-averaged approximations. In this section, we give an overview of a formal "second order" extension based on the ideas available in [START_REF] Camarri | A low-diffusion muscl scheme for les on unstructured grids[END_REF][START_REF] Cournède | A Positive MUSCL Scheme for Triangulations[END_REF][START_REF] Koobus | Nonlinear V6 schemes for compressible flow[END_REF]. The resulting scheme is very close to the frictionless case [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF], and we refer to this paper for more details. Following the general strategy of MUSCL methods [START_REF] Van-Leer | Towards the ultimate conservative difference scheme. V -A second-order sequel to Godunov's method[END_REF], a gain in space accuracy is reached by the mean of a preliminary reconstruction of the flow variable at each interface. For a given edge of the initial triangulation T , let's define the downstream and upstream triangles T D and T U , having respectively i and j as a vertex and such that the edge ij intersect their opposite edge (Fig. 5). The reconstruction method is performed on the augmented vector variable V = (η, q x , q y , h), and involves the three following gradients:

T D T U C s C r C q C p C i M i Γ ij C j M j
∆ V - ij = ∇ VTD • ij , ∆ V c ij = Vj -Vi , ∆ V ho ij = 1 3 ∆ V - ij + 2 3 ∆ V c ij , (46)
where ∇ VT stands for the P 1 gradient approximation of V on the triangle T . These quantities are then injected in a three entries continuous limiter, defined as follows:

L(a, b, c) = 0 if sgn(a) = sgn(b), sgn(a) min(|2a|, |2b|, |c|) otherwise . ( 47
)
Introducing the notation

L ij ( V ) = L(∆ V - ij , ∆ V c ij , ∆ V ho ij )
, the new interpolated values Vij and Vji at the edge Γ ij are defined by:

Vij = Vi + 1 2 L ij ( V ) , Vji = Vj - 1 2 L ji ( V ) . ( 48
)
At this stage, new face values V ij = (η ij , hu ij , hv ij ) and V ji = (η ji , hu ji , hv ji ) are extracted from [START_REF] Koobus | Nonlinear V6 schemes for compressible flow[END_REF], and the "left" and "right" states for the topography are recovered as follows: We finally apply steps (21-22 -23), and inject the resulting states in the numerical fluxes and source terms in [START_REF] Hassan | Numerical solution of the rotating shallow water flows with topography using the fractional steps method[END_REF]. The resulting method provides a better space accuracy and the general characteristics of the scheme are preserved. Particularly, no additional source term has to be introduced in order to maintain the well-balanced property, as stated in [START_REF] Duran | On the well-balanced numerical discretization of shallow water equations on unstructured meshes[END_REF]. The reader is also referred to this paper for the proof of robustness.

z ij = η ij -h ij , z ji = η ji -h ji . (49)

Numerical Validations

Accuracy validation

This first test case is devoted to an extensive accuracy analysis and convergence rate study. We consider a rectangular channel with dimensions [0, 1] × [0, 0.1] and study the collapsing of an exponential profile:

η(x, y, 0) = 1 + a η exp - x -0.5 p 2 , (50)
over a topography supplied by the following expression:

z(x, y) = 1 2 (1 -a z exp - x -0.5 p 2 . ( 51
)
The constants are taken to a η = 0.05, a z = p = 0.2, and the flow is initially motionless. Closed boundaries are considered, and the Manning friction coefficient is set to n = 0.1.

As no exact solution is available, a reference solution is computed at t = 0.1s with the MUSCL scheme on a regular mesh characterized by ∆x = ∆y = 1 20 × 2 6 , following the splitting indicated on Fig. 6. A series of computation is then run with an increasing level of refinement: ∆x = ∆y = 1 20 × 2 p , p = 0 • • • 5, with the first and formal second order schemes. These solutions are compared to the reference, through the evaluation of the L 1 , L 2 and L ∞ errors. Numerical results are available on Tab. 1. The expected convergence rates are reached for the L 1 , L 2 norms, and this in spite of the presence of a rough friction coefficient. As concerns the L ∞ error, even though more accurate results are obtained, we observe that the benefits induced by the use of the MUSCL scheme are less significant in terms of convergence rates. Those obtained for both water height and discharge using the L 1 norm can be assessed on Fig. 7, where the error is plotted against the mesh step in a log-log scale.

Lastly, let us mention that several other values of roughness coefficient have been considered, from n = 0 to n = 1, and we did not observe any significant variability as regards the behaviour of the numerical error. Different simulation times have also been tested, leading to very similar results. This preliminary test provides relevant assessment data regarding the efficiency of the MUSCL method in the context of smooth solutions, and highlights the stability of the numerical scheme regarding the inclusion of possibly steep friction terms. We refer to Test 3.3 for an additional error analysis in the presence of dry cells.

Dam break with friction

We consider a classical 1d dam break problem with friction on a flat bottom, and focus on the treatment of dominant resistance effects in the neighbourhood of the wet/dry front. In a 2d framework, the test involves a rectangular flume with dimensions [-10, 10] × [0, 4]. The initial water depth is set to 1m at the left of the dam (x ≤ 0), and 0 elsewhere. For this test, we employ a Darcy friction law (4b) with d = 0.05, and use a 8 241 nodes regular grid. This mesh is built on the model given by Fig. 6, with ∆x = ∆y = 0.1m. Solid wall conditions are assumed at each boundaries.

One can build a relevant approximation of the exact solution decomposing the domain in two regions behind the wave front location : in the first area, friction effects are neglected and the analytical solution is provided by an ideal fluid flow model, while in the wave tip region bed the resistance phenomena are considered as dominant compared to acceleration and inertial effects. Physical and mathematical aspects of such construction are detailed in [START_REF] Chanson | Analytical solution of dam break wave with flow resistance. application to tsunami surges[END_REF]. Fig. 8 (top) shows some profiles of the water depth along the x-direction centreline, until t = 2.4s. We can observe a good agreement with the analytical solution, for both first order and MUSCL schemes. If the benefits of the MUSCL reconstruction principally appear during the handling of the shock for the water depth, we also remark a more accurate evaluation of the front speed on Fig. 8 (bottom). The flood wave front location seems to be accurately computed, which tends to validate the ability of the current scheme to deal with friction when vanishing water heights and dry cells are involved. We can also point out the good concordance with the predictions provided by the approach of Céa et al in [START_REF] Cea | Unstructured finite volume discretization of bed friction and convective flux in solute transport models linked to the shallow water equations[END_REF] for the water depth.

Moving boundary over a parabolic bottom

The following test is adapted from a 1d version [START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF] of the experiments initially set up by Thacker [START_REF] Thacker | Some exact solutions to the nonlinear shallow water wave equations[END_REF]. Recently, several authors proposed some enhancements on the family of Thacker's test cases, introducing some friction terms in the exact solutions. Sampson et al [START_REF] Sampson | Moving boundary shallow water flow above parabolic bottom topography[END_REF] proposed a first advancing for the planar case over a parabolic bottom in 1d, considering a linear friction law (4c). Since then, this test is regularly employed [START_REF] Kesserwani | Well-balanced RKDG2 solutions to the shallow water equations over irregular domains with wetting and drying[END_REF][START_REF] Kesserwani | Locally limited and fully conserved RKDG2 shallow water solutions with wetting and drying[END_REF][START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF], and more recently, a 2d extrapolation has been proposed by Wang et al [START_REF] Wang | A 2d shallow flow model for practical dam-break simulations[END_REF]. We refer to this paper for a complete 2d simulation with friction. Here, computations are run on a 8 640m × 500m channel, regularly meshed with discretisation steps ∆x = 24m and ∆y = 25m, until T = 10 000s. The topography is the following :

z(x, y) = h 0 ( x a ) 2 -1 .
We are accordingly working with a bed friction source term of the form (4c), setting κ = 0.001. Enforcing the y-component of the velocity to zero for the 2d derivation, we consider the following exact solution of the NSW equations :

η(x, y, t) = h 0 + a 2 B 2 e -κt 8g 2 h 0 -sκ sin(2st) + ( κ 2 4 -s 2 )cos(2st) - B 2 e -κt 4g - e -κt/2 g Bs cos(st) + κB 2 sin(st) x , u(t) = Be -κt/2 sin(st), (52)
where s = p 2 -κ 2 and p = 8gh 0 /a 2 . The coefficients are a = 3 000, B = 5 and the still water level h 0 is fixed to 10m. Straightforward computations give ∂ x η(0, x, y) = cte and q(0, x, y) = 0, so that the initial condition consists in a motionless planar flow (Fig. 9). Note that the problem of inflow and outflow boundary conditions is anecdotal for this test, since the flow is not supposed to reach the corresponding boundaries. According to the exact solution, the free surface is supposed to remain planar along the channel, and oscillate with a decreasing amplitude, under the effects of the resistance. As t → ∞, the velocity vanishes and the free surface asymptotically settles back to the water level at rest h 0 . Considering the relative poor number of analytical solutions involving friction flows, we dispose of an interesting test case to assess the capacity of dealing simultaneously with complex topography, friction and dry states. We first focus on the evolution of the water height during a half period (Fig. 10), and observe a very good agreement with the exact solution. Having knowledge of the theoretical location of the moving water interface :

x = a 2 e -κt/2 2gh 0 -Bs cos(st) - κB 2 sin(st) ± a ,
comparisons are also run on the time history of the wet/dry front evolution (Fig. 11), and exhibit the capacity in correctly handling flooding and drying. In conformity with theoretical predictions, the numerical approximation is progressively damped by friction up to the apparition of the motionless steady state. As far as accuracy is concerned, we proceed to quantifications of the L 1 error in a series of simulations involving coarser regular triangulations (∆x = 24, 48, 96, 192), based again on Fig. 6. The error is taken at t=200s and evaluated in L 1 norm for both water height and normal discharge. We can observe the significant improvements provided by the MUSCL scheme in terms of convergence rate on Fig. 12. The numerical error is plotted with respect to ∆x, in a log-log scale for both water height and normal discharge hu. Slopes around 1 and 1.55 are respectively observed for the first order and MUSCL schemes. Numerical results are available in Tab. 2, and clearly confirm the efficiency of the MUSCL reconstruction discussed in §2.3, even when combined with the friction approach.

Periodic subcritical flow

An important step in the validation of a numerical model for shallow water flows is to confirm the ability to converge toward steady states. Some classical test cases involving converging frictionless flows can be found in the literature (see for instance [START_REF] Audusse | A well balanced positivity preserving "second-order"scheme for shallow water flows on unstructured meshes[END_REF][START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF]), but far less involving resistance. In [START_REF] Delestre | A numerical scheme for a viscous shallow water model with friction[END_REF], in a one dimensional framework, an approach able to generate new steady states solutions for a Shallow Water model with viscous and friction terms is developed. Returning to a Manning formulation (4a), we choose to derive one of them for a 2d experiment, considering a 5 000m × 500m channel, in which we study a converging subsonic flow. Numerical predictions are compared to the reference solution provided by the steady state for the free surface :

h ref (x) = 9 8 + 1 4 sin( πx 500 ) . (53)
The topography computation is run following the 1d iterative method described in [START_REF] Delestre | A numerical scheme for a viscous shallow water model with friction[END_REF]. The discharge is enforced at the left boundary, prescribing 2m.s -1 , as well as the water height, set to h ref (0). Neumann conditions are set a the outflow boundary. The Manning coefficient set to n = 0.03, and we run computations on a regular mesh with a space step ∆x = ∆y = 100/3. Fig. 13 shows the evolution of the total water height at several times during the transient part, obtained with the first order scheme; the convergence toward the steady state is clear. It should be stressed that similar observations can be made with the MUSCL scheme. The absence of instabilities despite the possible hazards resulting from dry areas and resistance terms confirms the efficiency of the friction approach. The results are very close to the ones provided by the 1d scheme [START_REF] Berthon | An efficient scheme on wet/dry transitions for shallow water equations with friction[END_REF]. Note that when equilibrium is reached, the flow is not static and friction effects are still present. As a consequence, these results are of interest since they actually show that the current scheme may be able to restore a larger class of steady states than those involved by the C-property.

Two dimensional steady flow with friction

In this test, extracted from [START_REF] Murillo | The influence of source terms on stability, accuracy and conservation in two-dimensional shallow flow simulation using triangular finite volumes[END_REF], we study now the convergence toward a truly 2d steady state in the presence of non trivial topography and friction effects. The computational domain is the square [0, 10] × [0, 10], and we use an unstructured mesh made of 3054 vertices. A constant discharge is assumed on the whole domain :

q x = q y = 0.1 ,
and the water height at steady state is defined as follows :

h(x, y) = d + q x x + q y y , with d = 0.5. The analytical function for the associated topography can be written as :

z(x, y) = - 1 2g q 2 x + q 2 y h 2 -h + 3 7 q x n 2 √ 2 h 7/3 .
The total free surface elevation is set to zero at the beginning of the simulation, the final time being t=300s. In a first time the case of a Manning roughness coefficient n = 0.03 is investigated. This configuration leads to the free surface elevation available on Fig. 14, toward which is expected to converge the solution. A comparison between the first order approximation and the analytical solution is available on Fig. 15 at the end of the simulation, exhibiting a good correspondence. Similar results are obtained with the MUSCL extension of the scheme. As done in [START_REF] Murillo | The influence of source terms on stability, accuracy and conservation in two-dimensional shallow flow simulation using triangular finite volumes[END_REF], we propose to focus on the evolution of the L 1 numerical error to evaluate the relative performances of the first order and MUSCL approaches. In accordance with the observations made by the authors, Fig 16 shows that in the presence of slight perturbations induced by friction, the high order scheme fails in reaching a perfect equilibrium, while the first order solution is able to restore the steady state with a better level of precision.

We consider now the case where the flow is strongly dominated by the resistance terms, taking n = 0.3. A 3d view of the corresponding steady state is proposed in Fig. 17. Again, we initially consider a still water configuration, and run the simulation until t=300s. An overview of the free surface extracted from the first order solution is proposed on Fig 18, together with the exact profile. In spite of the steepness of the friction effects, particularly pronounced in this case, the scheme is able to give stable and accurate results. Note here that, probably tempered by friction, the behaviour of the numerical error rapidly tends to stabilize, even for the MUSCL scheme, still slightly less accurate in this case (see Fig. 19).

Tsunami wave on a sloping beach

Now, let us focus on a 2d case, with the study of a wave surge over a beach with complex topography. This test has initially been performed by Zelt to validate a Lagrangian model for shallow flows [START_REF] Zelt | Tsunamis: the response of harbors with sloping boundaries to long wave exitation[END_REF]. The dimensions of the coastal area are set to [-10, 15] × [-10, 10]. The computational domain is meshed with a 8 413 nodes unstructured triangulation, deliberately refined in the neighbourhood of threatened dry areas (Fig. 20). The high order MUSCL reconstruction is chosen here. Denoting L the half-width of the bay, the bed slope source term is expressed as follows :

z(x, y) =    H 1 + x -x p L cos(πy/L)/π + x p if x ≥ x p , 1 elsewhere.
(

) 54
Considering a cross-shore section, we observe that the slope of the topography slowly decreases from the lateral boundaries to the x-direction centreline. The water depth at rest is initially set to H = 1.273m, and we set x p = 3 L/π. At the left side of the domain, an impulse for the generation of the wave is furnished by the following offshore boundary condition : the bottom profile, steeper at the upper and lower boundaries, part of the energy is gradually transferred into the center of the basin. This preferential direction of the propagation is also supplemented by a precocious run-down issuing from the reflection at the lateral emerging boundaries. Our numerical model is able to correctly describe the flooding phenomena involved in this test, as well as the entire development of the reflected wave until it goes down to the open sea. To illustrate this, we propose on Fig. 22 the time series of the run up along five cross sections, until T = 60s. These results can be compared with other observations available in the literature, as [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. We do not notice any significant discrepancies with these works, or even when compared against the model of Zelt, which attests of a good reproduction of the two-dimensional process. Now, we study the impact of a supplemented physical resistance over the evolution of the water front. To do so, we observe the results obtained with a Manning friction law (4a) for increasing values of κ = n 2 . Middle sections of the total free surface along the x-direction are available at several stages of the propagation on Fig. 23, for increasing friction coefficients (κ = 0, 0.1, 0.25, 0.5 and 1), and clearly highlight the loss of amplitude resulting from the consideration of friction. The reduction of the flooded area can be observed on Fig. 24, offering an overview of the rectangular domain and pointing all the wet cells at the level of the beach for the different values of κ involved in the simulations. It is noteworthy that the same CFL has been used for all these computations, with the same tolerance value as in the frictionless case to identify dry areas, precisely ǫ = 1.e -6. In all cases, in accordance with the objectives sought, we did not have to face any instability problem.

h(t) = H + αHsech 2 { √ gH L χ t} , (55)

Toce River dam break

We propose to study a dam break flow in a realistic river, employing the benchmark test of the Toce river valley (Northern Alps, Italy). Topographic data and reference results are provided by a 1:100 scale physical model, allowing comparisons with our numerical observations. In the same way, inflow boundary conditions are calibrated on the discharge hydrograph used to initiate the motion in the laboratory studies. Opened boundary conditions are set downstream. This benchmark test was used within the CADAM project [START_REF] Soares-Frazão | The Toce River test case: numerical results analysis[END_REF], and is also performed in [START_REF] Brufau | A numerical model for the flooding and drying of irregular domains[END_REF][START_REF] Chang | Numerical simulation of shallow-water dam break flows in open channels using smoothed particle hydrodynamics[END_REF][START_REF] Delis | Numerical flood simulation by depth averaged free surface flow models[END_REF][START_REF] Wang | A 2d shallow flow model for practical dam-break simulations[END_REF][START_REF] Ying | Numerical simulation of flood inundation due to dam and levee breach[END_REF][START_REF] Ying | Improved hll scheme for 1d dam-break flows over complex topography[END_REF]. Computations are run with the first order scheme on an unstructured mesh composed of 30 529 elements, and a Manning roughness coefficient equal to n = 0.0162. We stress out the difficulties represented by this test, since it is expected to describe the evolution of an advancing front dragged in a complicated bed river, and submitted to resistance laws. The geometry of the floodplain is indeed not trivial, and induces a flow with a complex structure, subject to multiple variations of regime throughout its propagation. Some 3d views of the free surface are available in Fig. 25, describing the overall propagation of the flood wave. These indications on flooded areas can be put in comparison with some of the reference previously mentioned to corroborate our results. The main steps of the flow evolution and corresponding times seems to be respected, from the inundation of the central reservoir to the arrival at the downstream basin. The relevance of our predictions can be further assessed analysing the arrival time of the flood wave at several gauges disposed along the bed river, for which we have reference data. We report a good concordance with the data of the physical model, in the same order of quality with the studies above mentioned. This simulation has also been performed with the MUSCL scheme and we note here that the high order resolution helps in obtaining a better matching with the experimental values on the last gauges, as highlighted in Fig. 26.

Malpasset dam break

As a final experiment, we perform the classical simulation of the Malpasset dam break. A 3d view of the narrow gorge of the Reyran River Valley (South of France) is available in Fig. 27, locating the dam and its upstream and downstream regions. In December 1 959, the left side of the dam collapsed under the effects of pressure, after a heavy rain fall; more than 55 000m 3 of water started to pour in the valley, and the generated flood wave reached the city of Frejus, located 12km downstream. The context of this benchmark test for dam break models involves complex bathymetry and geometry, and is de facto supplemented by occurrence of dry areas. As we also have to consider friction effects, this real life application is particularly appropriate to assess the performances of the current scheme. Reconstruction of the disaster during laboratory studies furnished reference data for the arrival time of the water front at several gauges located in the floodplain. Concerning the implementation, we set an upstream total water height of 100m, and assume the bed river initially dry. We make use of the first order scheme on a vertex-centred mesh of 13 541 nodes, and the Manning roughness coefficient is fixed to n = 0.03. We follow the evolution of the flood wave on Fig. 28, where snapshots of the free surface are available as the flood wave reaches gauges 8, 10 and 14. Numerical results are compared with some data issuing from the experiment on Fig. 29a and29b, where the arrival time of the flood wave is available at gauges 6 to 14, as well as the maximum elevation. We also propose a comparison on the maximum elevation with the police measurements at some survey points on Fig. 29c, that also exhibits a good agreement. Note that significant discrepancies are observed on Fig. 29b in absence of friction, highlighting the important role played by these terms in this test.

From a qualitative point of view, numerical predictions are in the same order of precision than other studies reported in the literature (see [START_REF] Brodtkorb | Efficient shallow water simulations on GPUs : Implementation, visualization, verification and validation[END_REF][START_REF] Singh | Two-dimensional numerical modeling of dam-break flows over natural terrain using a central explicit scheme[END_REF][START_REF] Wang | A 2d shallow flow model for practical dam-break simulations[END_REF]). Again we can conclude on the good accuracy and stability properties provided by our friction approach in the context of real life applications.

Conclusion

In this paper we have gathered and exploited some recent results to obtain two-dimensional resolutions of the shallow water system on unstructured meshes, in the presence of friction and varying topography. The numerical model is the result of an appropriate extension of the 1d approaches proposed in [START_REF] Liang | Numerical resolution of well-balanced shallow water equations with complex source terms[END_REF][START_REF] Wang | A 2d shallow flow model for practical dam-break simulations[END_REF] for the pre-balanced SWE. The extension to flows with friction is inspired from other 1d works [START_REF] Berthon | Asymptotic preserving hll schemes[END_REF][START_REF] Berthon | An efficient scheme on wet/dry transitions for shallow water equations with friction[END_REF], with a particular focus on the problem of vanishing water depths. It results from these investigations a simple and well-balanced scheme, able to account for general friction laws and to guarantee stable computations in the presence of partly dry domains. These points stand for a notable feature in the context of unstructured triangulations, and are particularly important in the perspective of practical applications to environmental flows.

Figure 1 :

 1 Figure 1: Non conservative variables for SWE equations.

Figure 2 :

 2 Figure 2: Non-negative reconstruction for motionless steady state configurations. Initial state (left ) and reconstructed state (right ). Arrows indicate the impact on the topography.

Figure 3 :

 3 Figure 3: Vertex centered formalism -focus on the interface Γ ij .

Figure 4 :

 4 Figure 4: Left and right Riemann states for the 1d and 2d schemes.

Figure 5 :

 5 Figure 5: Upstream and Downstream triangles T U and T D .

Figure 6 :

 6 Figure 6: Example of regular mesh.

Figure 7 :

 7 Figure 7: Accuracy validation : Convergence rate analysis for the water height (left ) and the discharge (right ).

Figure 8 :

 8 Figure 8: Dam break with friction : Water depth (top) and velocity (bottom) profiles at t=0.6, 1.2, 1.8 and 2.4s. Analytical solution compared with first order and MUSCL approximations.

Figure 9 :Figure 10 :

 910 Figure 9: Moving boundary over a parabolic bottom : Initial condition.

  h

Figure 11 :

 11 Figure 11: Moving boundary over a parabolic bottom : Time history of the wet/dry interface.

Figure 12 :

 12 Figure 12: Moving boundary over a parabolic bottom : Convergence rate analysis for the water height (left ) and the discharge (right ).

Figure 13 :

 13 Figure 13: Periodic subcritical flow : Free surface profile evolution -Convergence toward h ref .

Figure 14 :

 14 Figure 14: Two dimensional steady flow with friction : 3d view of the exact free surface for n = 0.03.

Figure 15 :

 15 Figure 15: Two dimensional steady flow with friction : Predicted (left ) and exact (right ) free surface for n = 0.03.

Figure 16 :

 16 Figure 16: Two dimensional steady flow with friction : evolution of the numerical error for first order and MUSCL schemes for n = 0.03.

Figure 17 :

 17 Figure 17: Two dimensional steady flow with friction : 3d view of the exact free surface for n = 0.3.

Figure 18 :

 18 Figure 18: Two dimensional steady flow with friction : Predicted (left ) and exact (right ) free surface for n = 0.3.

Figure 19 :

 19 Figure 19: Two dimensional steady flow with friction : evolution of the numerical error for first order and MUSCL schemes for n = 0.3.

Figure 20 :

 20 Figure 20: Tsunami wave on a sloping beach : Mesh of the computational domain.

Figure 21 :

 21 Figure 21: Tsunami wave on a sloping beach : 3d surface profiles at t=11, 14 and 17s (frictionless flow).

1 Figure 22 :

 122 Figure 22: Tsunami wave on a sloping beach : Time series of the run up along several cross sections.

Figure 23 :

 23 Figure 23: Tsunami wave on a sloping beach : Water depth profiles at t=2, 6 and 14s for increasing Manning coefficients.

Figure 24 :

 24 Figure 24: Tsunami wave on a sloping beach : Sketch of flooded areas.

Figure 25 :

 25 Figure 25: Toce River dam break. 3d views of the free surface at t=20, 40, 50 and 60s.

Figure 26 :

 26 Figure 26: Toce River dam break. Free surface history at several gauges for the first order (squares) and MUSCL (circles) schemes. Experimental data appear in solid lines.

Figure 27 :

 27 Figure 27: Malpasset dam break : Initial configuration.

Figure 28 :

 28 Figure 28: Malpasset dam break : Mesh and evolution of the flow at gauges 8, 10 and 14.

  Maximum elevation at survey points.

Figure 29 :

 29 Figure 29: Malpasset dam break : Comparison with experimental data.

Table 1 :

 1 Accuracy validation : Convergence analysis -L 1 numerical error quantification and corresponding convergence rates for n = 0.1.

	∆x	1/20	1/40	1/80	1/160 1/320 1/640 order
	L 1			1 st order
	h	4.7e-4 2.5e-4 1.3e-4 5.9e-5 2.6e-5 8.6e-6	1.1
	q	1.5e-1 7.1e-2 3.0e-2 1.4e-2 6.0e-3 1.9e-3	1.3
				MUSCL
	h	4.5e-5 1.2e-5 3.1e-6 6.2e-7 1.9e-7 2.9e-8	2.1
	q	3.3e-2 1.2e-2 2.0e-3 1.1e-3 3.3e-4 2.4e-5	2.1
	L 2			1 st order
	h	6.1e-4 3.2e-4 1.6e-4 7.6e-5 3.3e-5 1.1e-5	1.1
	q	1.6e-1 7.4e-2 3.1e-2 1.5e-2 6.1e-3 1.9e-3	1.3
				MUSCL
	h	8.0e-5 2.4e-5 7.3e-6 1.7e-6 5.2e-7 1.5e-7	1.8
	q	3.2e-2 1.1e-2 2.9e-3 1.2e-3 3.5e-4 6.0e-5	1.8
	L ∞			1 st order
	h	8.0e-4 4.3e-4 2.1e-4 1.0e-4 4.3e-5 1.5e-5	1.1
	q	1.0e-3 4.8e-4 2.0e-4 9.2e-5 3.8e-5 1.2e-5	1.3
				MUSCL
	h	1.7e-4 6.5e-5 2.5e-5 6.2e-6 1.9e-6 1.1e-6	1.5
	q	2.2e-4 7.6e-5 3.8e-5 1.9e-5 6.7e-6 2.8e-6	1.3

Table 2 :

 2 Moving boundary over a parabolic bottom : Convergence analysis -L 1 numerical error quantification.

	hu


Attribution theory meets compensatory control

Correlational study

Experimental study A

Experimental study B

• Yes! Exposing consumers to the notion of order prior to goal failure prevents them from blaming the brand.

• Next step: Differentiating controllability when locus is firm-related: volitional vs. constrained

Conclusion and next step

Research question: Does "forcing" order into the consumers' environment prevent them from blaming the brand(s) involved in goal failure?

  Folkes, V. S.(1984). Consumer reactions to product failure: An attributional approach. Journal of Consumer Research, 10(4), 398-409. Kay, A. C., Gaucher, D., Napier, J. L., Callan, M. J., & Laurin, K. (2008). God and the government: testing a compensatory control mechanism for the support of external systems. Journal of Personality and Social Psychology, 95(1), 18-35. Kay, A. C., Laurin, K., Fitzsimons, G. M., & Landau, M. J. (2014). A functional basis for structureseeking: Exposure to structure promotes willingness to engage in motivated action. Journal of Experimental Psychology: General, 143(2), 486= 426 online participants Design: two-condition (goal attainment vs. goal failure) between-subjects, based on a memory task

Figure 1 :Figure

 1 Figure 1: Structural model


From euroregions to macroregions: a reasonable step

Introduction

Macro-regional strategies of the European Union are conceived as new instruments for territorial cooperation. They must integrate all the available and existing resources and schemes within the same transnational geographical area: structural funds and territorial cooperation programmes; other national, local and transnational planning schemes; public and private initiatives, etc. These strategies are expected to articulate functional and political complexity and to improve policy e ciency through cooperative and flexible approaches, without creating new institution or funding 1 .

This territorial approach evokes the concept of "mega-region", i.e. a "meso-scale" transnational territorial unit, located inside a supranational "region", such as for instance the European Union, North America, Mercosur, East Asia 2 . Moreover, macro-regional scale was already at stake in the transnational areas designed from the mid-1990's in the strand of territorial cooperation programmes, or in cooperative organisations like the Benelux, the Visegrad Group, the Nordic Councils or the Euro-Mediterranean partnership.

While these macro-regional strategies are expected to foster European territorial cooperation, they raise the question of how they can e ciently interact with the organisations that are considered to be the most successful schemes of territorial cooperation so far: the cross-border cooperative organisations commonly known as 'euroregions' 3 . Whereas around thirty euroregions existed at the beginning of the 1990s, there are currently about a hundred organisations of this type. Poland for instance is a country (cross) bordered by around 15 euroregions and in this context an association of Polish euroregions was founded.

In the framework of the Western Mediterranean macro-region, we can wonder for instance whether two existing euroregions in this area -Euroregion Pyrénées-Méditerranée and Euroregion Alpes-Méditerranée -can be milestones in the constitution of a macro-regional strategy for this area. The case of cultural cooperation is an interesting example to address this question. Indeed, among the various policies that are developed within euroregions, culture appears to be a frequent domain of cooperation. A recent study shows that culture and tourism are the favoured themes of cross-border cooperation 4 .

Cultural cooperation: a vector for macro-regionalisation?

Euroregion Pyrénées-Méditerranée was constituted in 2004 between the French regions of Languedoc-Roussillon, Midi-Pyrénées and the Spanish autonomous communities of Aragon, Catalonia and Balearic Islands. Culture has been a mainstream domain of cooperation, through a specific agenda that includes a culture portal online, specific calls for projects, annual cultural forums. Culture was also one of the first attributions of the European Grouping of Territorial Cooperation (

) that was launched in 2010 to reinforce the structure and governance of the euroregion 5 . Euroregion Alpes-Méditerranée was launched in 2007 and gathers the In this context, different elements suggest that the interaction between both euroregions can help to implement a macro-regional strategy in this area, specifically starting with cultural cooperation 6 : they encompass focal key areas of the Western Mediterranean, with the Barcelona-Marseilles-Genoa Arc and the metropolitan regions of Toulouse, Lyon, Torino; thus constituting a potentially homogenous territorial entity like "North of Southern Europe 7 " or "Golden Mediterranean Belt 8 ". There already exist cooperation schemes that can be integrated within a macro-regional strategy: cooperation programme , cooperation Rhône-Alpes/Catalonia in the framework of the association "Four Motors for Europe", Union for the Mediterranean. These two euroregions show potential geo-cultural interactions: transnational historical and identity links, common fields of "creative excellence" that are already mobilised in the euroregional cultural policies, like for instance audiovisual and image industries.

The two euroregions showed first step of cooperation with a joint participation to the Open Days in Brussels within the grouping "Euroregions for territorial integration" in 2009, 2010 and 2011.

2. Actual and potential macro-regional strategy areas

:

-

The realpolitik limits

Even though the macro-regional interaction between euroregions Pyrénées-Méditerranée and Alpes-Méditerranée can prove to be an option for the constitution of a Western Mediterranean macro-regional strategy, several points 9 show a certain gap between the theoreticalideal level and the practical one.

First, some observers question the very relevance of a macroregional strategy for this area: can this really help reinforcing cooperation, does it have a real added value and does it respond to the specificities of the Mediterranean territorial cooperation context 10 ? Secondly, the evolution of both euroregional partnerships shows that they remain uncertain schemes in several aspects. In Euroregion Alpes-Méditerranée the whole cooperation project slowed down following the change of political majority in the Italian regions. The president of Piedmont, who was elected in 2010 with the support of the extremeright party Lega Nord, declared that the euroregional cooperation was not a priority anymore and most of the projects are in stand by. In Pyrénées-Méditerranée, recent financial cuts and changes in the

The interaction between both euroregions can help to implement a macro-regional strategy in this area, specifically starting with cultural cooperation François Alfonsi, The evolution of EU macroregional strategies: present practice and future prospects, especially in the Mediterranean, European Parliament, Committee of Regional Developement, 2012, p. 25.

Luiza Bialasiewicz, et al., "Re-scaling ' 'rope: macro-regional fantasies in the Mediterranean", European Urban and Regional Studies January, vol. 20, n° 1, 2013, p. 59-76. political majorities could start to impact the initial ambition and objectives of the euroregion.

Developing new territorial cooperation instruments at transnational scale can allow a step further in the territorial construction of the European Union. This achievement can be based on the previously developed cooperation schemes, especially at cross-border scale and in the cultural sector. Nevertheless, real political will and commitment are crucial guarantees of e ciency for these instruments. In this way, planning experts and political leaders could work together to find the balance between giving a response to citizens' expectations, while giving an impulse to meta-geographical projections.

-

  French regions Rhône-Alpes and Provence-Alpes-Côte d'Azur with the Italian regions Piemonte, Liguria and Val d'Aoste. The partners established a cultural agenda based on networking in the field of visual arts and on organising cultural events to promote the euroregion.

›

  Thomas Perrin is a lecturer in Planning and European Studies at Université Lille 1, research centre . Previously, he worked for the French Ministry of Culture and cultural institutes abroad (Morocco, Chad) and was a research assistant at Toulouse University. Thomas holds a PhD in Political Science from Institut d'Études Politiques, University of Grenoble, with a thesis that deals with cultural policy and cooperation within euroregions. In 2011, he was a post-doctoral fellow at the Universitat Autònoma of Barcelona. He is expert for the Compendium of Cultural Policies and Trends in Europe and member of : European Association of Cultural Researchers. His research and teaching deal with cultural policy and cultural relations on the one hand, and with the issues and dynamics of regionalism and territorial cooperation on the other hand. His works have been published in several books and journals. In 2011, he was awarded the Mark Schuster Prize on comparative cultural policy and a Merit Prize in the Committee of the Regions' Doctoral Thesis Competition. In 2012, Toulouse City Council gave him the award of "chercheur novélisé" during the Novela, Toulouse festival of science and knowledge.

  4. Euroregion Pyrénées-Méditerranée /: www.medalp.eu3. Euroregion Alpes-Méditerranée /: Marie-Laure Malaval, Université de Toulouse-Jean-Jaurès :

Stefanie Dühr, Mer Baltique, Danube et stratégies macro régionales: un modèle de coopération transnationale dans l'UE ?, "Notre Europe Études & Recherches", °86

, 2011. 2 Richard Florida, Tim Gulden, Charlotta Mellander, The rise of the mega-region, "Cambridge Journal of Regions, Economy and Society", Volume 1 Issue 3, November 2008, p. 459-476 3 Francesc Morata, Euroregions i integració europea, Documents d'Anàlisi Geogràfica.Vol. 56, °1, 2010, p. 41-56

Thomas Perrin, Culture et Eurorégions. La coopération culturelle entre régions européennes, Éditions de l'Université de Bruxelles, 2013 Hugues de Jouvenel, Maria-Àngels Roque (dir.), Catalogne à l'horizon 2010, Barcelone/Paris: Institut catala d'estudis mediterranis/ Economica, 1994, p. VI. Guy Baudelle, Jean Ollivro Catherine Guy, Les scénarios de l'espace européen, in Baudelle & Castagnède (dir.), Le Polycentrisme en Europe. Une vision de l'aménagement du territoire européen, Editions de l'Aube/ , 2002, p. 130.