{'fraction_non_alphanumeric': 0.03174288756282962, 'fraction_numerical': 0.011241240993388961, 'mean_word_length': 4.438728790735254, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A well balanced fvc scheme for 2d shallow water flows on unstructured triangular meshes

We consider in this work the numerical resolution of a 2D shallow water system with a Coriolis effect and bottom friction stresses on unstructured meshes by a new Finite Volume Characteristics (FVC) scheme, which has been introduced in the preliminary works that will be cited below. Our main goal is to extend this approach to 2D unstructured formalism while preserving the physical and mathematical properties of the system, including the C-property. First, we present our extension by preserving the advantages of the finite volume discretization such as conservation property and the method of characteristics such as elimination of Riemann solvers. Afterward, an approach was applied to the topography source term that leads to a well-balanced scheme satisfying the steady-state condition of still water. A semi-implicit treatment will also be presented in this study to avoid stability problems for the other source terms. Finally, the proposed finite volume method is verified on several benchmark tests and shows good agreement with analytical solutions and experimental results; moreover, it gives a noticeable accuracy and rapidity improvement compared to the original approaches.


Introduction
Water is a crucial issue for poverty reduction, sustainable development, and achieving the Millennium Development Goals. However, until now, some 2.1 billion people, or 30% of the world's population, still do not have access to a safe water supply (www.who.int/news/18-06-2019).Climate change and uncontrollable human activities cause flood accidents to become more frequent in recent decades. As a result, integrated water resources management is necessary and even indispensable for preventing floods and droughts. The preservation of the environment, the prevention, and control of the impact of natural risks are at the heart of major socio-economic issues.
In water resources management, numerical modelling is still an essential tool, and we can cite numerous numerical modelling applications of free surface flows to the management of water resources, environmental and ecosystem protection: simulation of flows due to a dam break, diversions of floods from a river to a water retention area, the change process of a river bed simulation, sediment or pollutant transport simulation in estuary and coastal environments, (see e.g. [3,4,13] etc. ). In rivers, estuaries and coastal areas, flows are characterized by: great topographical and morphological complexity, strong affection, or even pure advection in the case of a dam break on a flat and slippery bottom (without friction), a variable space scale (starting from a dozen to a few thousand meters) and in time scale (starting from a few minutes to several days).
Consequently, when developing a numerical approach to solve a free surface flow or other models, one encounters major difficulties which result from the physical complexity of the area and numerical calculations. We consider in this paper a derivative as a formal first-order approximation of the three-dimensional free surface incompressible Navier-Stokes equations, using the so-called shallow-water assumption [17]. The difficulty in defining accurate numerical schemes for such hyperbolic systems is related to their non-linear behaviour or, more generally, mathematical structure and the physical phenomena they generate. In particular, the presence of a shock front essentially causes numerical oscillations or artificial scattering, which are due to the treatment of the advection terms in the equations governing the transport of the water mass by a standard method of approximation. Another fundamental point is to get schemes that satisfy the preservation of steady states, such as still water equilibrium in the context of the shallow water system. Different approaches to satisfy the well-balanced property have been proposed (see, e.g. [5,9,10,26,34]), and recent extensions to other types of homogeneous solvers can be found in [12,28]. The extension of ENO and WENO schemes to shallow water equations has been studied in [41]. Unfortunately, most ENO and WENO schemes that solve real flows correctly are still very computationally expensive. On the other hand, numerical methods based on kinetic reconstructions have been studied in [31], but the complexity of these methods is relevant. However, most of the above mentioned works, even if they are unstructured two-dimensional methods they lead to rather complex and time consuming algorithms. Other approaches are more efficient, but to our knowledge, they are limited to one-dimensional problems or to two-dimensional Cartesian meshes.
Our main objective in the present study is to develop a class of Eulerian-Lagrangian methods and to accurately solve shallow water equations in 2D-unstructured meshes without relying on Riemann solvers. The proposed FVC scheme belongs to the class of methods that employ only physical fluxes and averaged states in their formulations. It can be interpreted as a predictor-corrector scheme. In the corrector stage, the considered equations are integrated over an Eulerian time-space control volume, whereas, in the predictor stage, the conservation laws are rewritten in an advective form and integrated along the characteristics defined by the advection velocity field. This approach has shown its effectiveness through several test cases represented in the works, [8,7,6], but the authors of these three papers remained limited in the Cartesian mesh formalism but, as I have mentioned before, the real problems are characterised by a great topographical and geometrical complexity, hence, the limit of this formalism. A more thorough study of the accuracy of this finite volume discretization method on unstructured meshes was always the objective in these works (see the conclusions of [7,6]). Therefore, the unstructured finite volume method doesn't only ensure the conservation of mass, which is an important property in the computation of fluid flows but also allows the complex geometry of the computational domain to be perfectly taken into account. For these reasons, we propose an extension of this scheme in an unstructured mesh.
Another strong point of this discretization method is that the Jacobian matrix of the system that caused the slowness of many approximation schemes is not involved in the calculation. Note also that many approximation schemes in the framework of conservation laws require a solver for the Riemann problem at each time step to reconstruct the numerical flow, which is completely avoided in our FVC scheme. Simply this approach is based on a combination of the characteristics method and a finite volume method. Unlike traditional finite volume methods, this technique integrates our equations along the characteristics curves so that numerical flux are easily calculated. We noticed that this approach has several advantages over others conservation law solution techniques. Indeed, the main features of such a finite volume Eulerian-Lagrangian scheme are on the one hand, the capability to satisfy the conservation property resulting in numerical solutions free from spurious oscillations, and on the other hand, the achievement of strong stability and high accuracy for numerical solutions containing shocks or discontinuities.
We first present a projection of the shallow water system in order to find a form of a transport equation on which we will apply the techniques of the characteristics method for the purpose of evaluating our unknowns at the interfaces of our unstructured mesh. We also present original developments to take into account boundary conditions, to reduce the diffusion of the scheme or to increase the efficiency of the method.
Second, we announce a new reformulation of the FVC scheme adapted to non-uniform triangular mesh "unstructured mesh" as well as the discretization of the flow gradients and source terms while keeping the equilibrium and the C-property. It will also be seen that the proposed scheme has the ability to handle calculations of slowly varying flows as well as rapidly varying flows over continuous and discontinuous bottom beds.
This article is organized as follows: Section 2 will be devoted to the presentation of the mathematical model and its projected velocity system, as well as its mathematical and physical properties. In Section 3, numerical methods are formulated for the reconstruction of the finite volume characteristics scheme on an unstructured meshes and the approximation of the source terms keeping the scheme Well-balanced. Section 4 is devoted to the numerical results of several test examples and some results interpretations. It is shown that our new approach achieves the expected accuracy and robustness. Section 5 contains concluding remarks and an outlook.

Governing equation
The 2D shallow water system for the free-surface flow with the Coriolis effect and bottom friction stress is formulated as where the unknowns are always the water height h(t, x, y) 0 and the horizontal speed mean u(t, x, y) = (u, v) T (t, x, y) ∈ R 2 . The parameter f c is linked to the angular speed of the earth's rotation, g is the gravitational acceleration, r(h, u) has various expressions, for example, the asymptotic derivation mentioned in [21] leads the authors to consider, at first order, a linear friction term. The quadratic form in the Manning-Strickler velocity is nevertheless the most widely used in river flow applications [15,24], so in this study we use the latter approximation such that the bottom's friction term r(h, u) is given by, r(h, u) = (r fx , r fy ) := η 2 gh −1/3 |u|u, such as η is the Manning roughness. The function Z(x, y) represents the bottom profile, see the Fig.1. In order to give the reader a global view on the shallow water system we propose to add other aspects in the right hand side of the second equation of system (1). For example we can add the wind's effect on the free surface, i.e. in the case where the wind is moving with high speed there is the friction term τ (h,ȗ) which is not negligible. The viscosity or diffusion term can also be added if we want to solve turbulence problems in free surface flow (see e.g [42]). In this case, the second equation of system (1) becomes where, τ (h,ȗ) = (τ sx , τ sy ) := 1 2C f |ȗ|ȗ, in whichȗ = (ȗ,v) T represents the wind speed andC f is the coefficient of wind friction with water. ν is the diffusion coefficient associated to the term, − → hu := ( (hu), (hv)) T .
For simplicity, in this study we will not deal with these two terms (τ (h,ȗ) and ν − → hu ), we rewrite the system (1) in a conservative vectorial form where, F(W ) = (F (W ), G(W )) T . Such that, we note A T is the transpose of a matrix A, Note that the equation (3) has to be solved in a bounded spatial domain Ω, with given boundary and initial conditions. In practice, these conditions depend on the phenomenon studied ( see the Section 4 where numerical examples are discussed).

Properties of the system
A simple calculation shows that the system is still written in a quasi-linear form where J F and J G are the Jacobian matrices of the fluxes Following the usual techniques [22,38], we define J F ((α 1 , α 2 )) := α 1 J F +α 2 J G . For any (α 1 , α 2 ) ∈ R 2 , the matrix J F ((α 1 , α 2 )) has three eigenvalues defined by The shallow water system is a first-order hyperbolic system of balance laws, and it is also strictly hyperbolic for h > 0 with real and distinct eigenvalues.

•) Equilibrium
An important property is related to the source terms, and the most studied balance family is related to the presence of topography's source term: the shallow water system admits non-trivial steady-states. They are characterized by For flows in complex geometry, it seems very difficult to numerically preserve all two-dimensional balances, except those that correspond to an area at rest and whose characterization is independent of the dimension considered.
This particular stationary state, known as the resting lake state, is important because many flows in lakes or coastal bays are perturbations around this balance. Therefore it is essential to prevent numerical anomalies from disrupting the approached solution. However, the preservation of the stationary states at the numerical level is not obvious to be achieved, and even the simplest one is not an exception. In fact, (6) correspond to a balance between flow terms and source terms, whose discretization are not correlated.
There are other categories of stationary states resulting from an equilibrium between the pressure term and the Coriolis term, i.e.
For example, at large scales, the atmospheric and oceanic flows bear most of the time the perturbations of this stationary state [30], therefore it is also very important to be represented in the approached solution. This balance presents a complexity added to the balance of the lake at rest because it involves non-zero speeds.
•) Entropy inequality The physical system we are dealing with here is a system of conservation laws, so the energy aspect is very important in this type of system, thus we propose to say a few words on this point in order to provide insight into the treatment of this notion when constructing approximate solutions under the constraint of preserving certain properties, for example, the decay of energy in the presence of friction source terms.
We call entropy solution to the shallow water system (see [20]), a weak solution which satisfies the following entropy inequality where E is a mathematical entropy (which is the mechanical energy see [27]), defined as E(h, u, Z) = h |u| 2 2 + gh 2 2 + ghZ.
This inequality becomes equality for regular solutions, in the absence of energy loss terms, notably friction, and remains inequality, for admissible discontinuous solutions, resulting from classical calculations; we invite you to read subsection 1.1 in [19]. The mechanical energy, which is easily verified as convex with respect to the conservative variables, thus acts here as a mathematical entropy. In the case of the system without source terms and in the 1D problem, the mechanical energy is only one of the entropies that must be associated with the system for the problem to be properly posed. In 2D, or when the system contains the source terms, there is no longer a complete family of mathematical entropies. Therefore, even if the inequality (7) alone is not sufficient for a rigorous mathematical study, it nevertheless ensures the presence of an additional bound on a certain positive function of the system unknowns and can provide information on the choice of a physical solution. The finite volume scheme presented in section 3 verifies innately the conservation properties. The energy decay property remains more difficult to satisfy.

The initial and the boundary conditions
To obtain a well-posed problem, we add to this system some initial condition and boundary conditions. In this paper, we consider only two types of boundaries: solid walls on which we prescribe a slip condition and fluid limits on which we prescribe one or two conditions depending on the type of flow (subcritical, supercritical). In some cases we use the injection and Neumann homogeneous boundary conditions. slip wall no-slip wall subcritical inflow subcritical outflow supercritical inflow supercritical outflow Here the subscripts r and l denote the right and left states respectively at a boundary cell interface, such as the local values of the Froude number is used to determinate whether the flow is subcritical or supercritical at a given time. The subscript b above denotes prescribed physical boundary values.

Finite volume discretization of the model in 2D unstructured formalism
In this section we present the FVC scheme for the discretization of the shallow water system (3). The method consists of two steps and can be interpreted as a predictor-corrector approach. The first step deals with the classical finite volume method, whereas in the second step, the reconstruction of the numerical flux by constructing the intermediate state using the characteristics method [45].

Discretization
The integral form of the system (3) can be written as where Ω is the domain of interest, ∂Ω is the boundary surrounding, n is the normal vector to ∂Ω in the outward direction, dV and dσ are respectively the surface element and the length element. The problem domain is first discretized into a set of triangular cells forming an unstructured computational mesh see Fig.2. The average of conserved variables is stored at the centre of each cell, and the edges of each cell define the faces of a cell control what is called "control volume". We use the following notations: • γ ij , boundary edge between the cells T i and T j , • n ij , unit normal to γ ij , outward to T i such as, n ji = −n ij .
For each triangular control volume, the system (3) is written as where W i is the average quantity of cell T i stored at the cell centre. The flux vector over each edge of the triangular cell and the discrete form of the integral is is the numerical flux computed at the interface between the cells T i and T j .
As explained in the Fig.2, γ ij is the edge surrounding the cell T i and N (i) is the neighbouring triangles of the cell T i . The intermediate solution W ij is reconstructed using the FVC scheme. The time discretization of (9) is performed by a first order explicit Euler scheme. The time domain is divided into N sub-intervals [t n , t n+1 ] with time step ∆t = t n+1 − t n for n = 0, 1, ...., N . W n is the value of a generic function W at time t n . The fully-discrete formulation of the system (3) is given by

The FVC scheme
In this subsection, we present a generalization of this scheme that was introduced in the preliminary works [7,45], the first was done in the Cartesian mesh and the second was done without taking into account the source term related to bathymetry and the balance produced by this term. For the corrector stage, we will use 2D finite volume formalism described in Section 3.1. Finally, the predictor stage and the final reformulation of the FVC scheme will be presented in this subsection.

Construction of the projected speed model
Let us discretize the spatial domain Ω with cells with, ∂T i is the border of the cell T i and N ele is the total number of element.
Integrating the equations (1) over the cell T i , the basic equations of the finite volume method obtained using the divergence theorem are given by where η = (n x , n y ) T the unit outward normal to the surface |T i | of the cell T i and tangential τ = (−n y , n x ) T where the normal velocity u η = un x + vn y and the tangential velocity u τ = vn x − un y (see Fig.3). In order to simplify the system (11), we do the following operations which can be rewritten the system (12) in a differential form as The system (13) can also be reformulated in the transport equation form as where, The aim of using this technique on the local coordinates of the two-dimensional shallow water system (1) in the control volume T i , is to reduce the dimension to an one-dimensional system (14) on each surface |T i | of this control volume.

Flux construction
we reconstruct the numerical flux Φ(W n ij , n ij ) using the method of characteristics. The fundamental idea of this method is to impose a regular grid at the new time level and to backtrack the flow trajectories to the previous time level, for more details see [36,37]. At the previous time level, the quantities that are needed are evaluated by interpolation from their known values on a irregular grid, we'll see about that later.

Method of characteristics
The characteristic curves associated with the equation (14) are solutions of the following Cauchy problem Note that X c (s) is the departure point at time s of a particle that will arrive at the interface γ ij in time t n + α∆t, see the Fig.4. The method of characteristics does not follow the flow particles forward in time, as the Lagrangian schemes do, instead it traces backwards the position at time t n of particles that will reach the points of a fixed mesh at time t n + α∆t. By doing so, the method avoids the grid distortion difficulties that the conventional Lagrangian schemes have. Hence, the solution of (15) can be expressed in an integral form as this integral can be calculated using an integral approximation method. In our approach, we used a Runge-Kutta 3 method to approximate the integral in (16) which is accurate enough. In order to complete the reformulation of the algorithm used, the departure points must be calculated once the characteristic curves are known. Therefore, the solution of the transport equation (14) is given by where U(t n , X c (t n )) is the solutions at the characteristic feet computed by the local least squares interpolation method. In other cases, the integral of the equation (17) can be calculated using a first-order approximation based on the rectangle method, the vector U n is reconstructed at the interfaces using such as, λ = (λ x , λ y ) and R = (R x , R y ). The weights parameters are given by formulas (see the Appendix). To approximate F(U, Z, f c ), (i.e. ∂ η (u η ), ∂ η (h + Z), ... ) we will need to approximate these derivatives in the interfaces, for that we use the diamond cell Fig.6. This cell is constructed by connection of centres of gravity (L, R) of cells T i , T j which share the interface γ ij and its endpoints S, N . We obtain the co-volume SRN L by this construction. One can assume that the gradient is constant on the co-volume SRN L. According to Green-Gauss theorem the approximation leads to where u N , u S , u R , and u L represent respectively the values of the quantity u in the point N , S, R and L. n LR is a unit normal vector of the co-volume face γ LR and |γ LR | is its length. The others co-volume interfaces and their normal vectors are labeled analogically. µ SRN L is the area of the co-volume SRN L.
After the discretization of the source terms (see subsection 3.2.4), the district equation system (18) leads to the following predictor step Once these projected states are calculated, the quantity W ij will be calculated using the following transformations hu n ij = (hu η ) n ij n x − (hu τ ) n ij n y , and hv n ij = (hu τ ) n ij n x + (hu η ) n ij n y .

Well-balanced FVC scheme: the discretization of the bathymetry source term
In order to be able to calculate realistic flows we now consider the case ∇Z = 0 R 2 and introduce a numerical discretization of the source terms. As discussed in paragraph 2.1.1, the treatment of source terms related to bathymetry in the shallow water system poses a challenge in many numerical methods. In our scheme, the approximation of the source term S n i is reconstructed in such a way that the C-property [10] is satisfied, i.e. to maintain a discrete local balance of the continuous stationary state in still water.
•) The hydrostatic balance we prove from the hydrostatic balance that the model of the projected speed preserves the stationary state of the lake at rest •) The projected speed model This result ensures the equilibrium property corresponding to the lake at rest, and therefore it is consistent with the continuous form of the system's equilibrium with bathymetry source term.
which is equivalent to To approximate the source terms, we proceed as follows. First, we decompose the triangle T i into three sub-triangles, as depicted in Fig.7. where N xij = (n ij ) x |γ ij |, and N yij = (n ij ) y |γ ij |. Then, the source term is approximated as with h 1 , h 2 and h 3 are the average values of h over T 1 , T 2 and T 3 respectively.
The same applies to the y-direction. Again the stationary flow condition h 1 +Z 1 = h j +Z j = H = cst, ∀j ∈ N (1) ⇒ h 1 +h j +Z 1 +Z j = 2H and H − . Thus, (29) gives Finally, A similar procedure leads to the following approximations of the other terms in (28) T2 Notice that h l , h k and h j are the average values of h, respectively, on the triangle T l , T k and T j . Summing up, using the fact that (N xij = −N xji ) so, the discretization (28) gives For this reconstruction, the source terms in (27) result in If you have noticed, we will need h 1 , h 2 and h 3 to be able to calculate the values of the integrals in equation (30) but the system. (31) has two equations for the three unknowns. To complete the system, we add the natural conservation equation, h 1 + h 2 + h 3 = 3h i . The following system gives us the values we need Analogously, the bottom values Z j , j = 1, 2, 3 are reconstructed in each sub-triangle of T i as Finally the source terms in (30) are approximated as with a similar equation for the other source terms in the y-direction.

•) Computation of the solution
Finally, we write the formally well-balanced FVC scheme after calculation of the interface values (23) and the bathymetry source term approximation (33) as

Semi implied treatment of friction term source
To avoid stability problems related to the bottom friction source term, a fractional semi implied treatment for this term is proposed. The idea is to evaluate the momentum in the system (3) by decomposing it into two equations where Φ hu (W ) represents the convection terms corresponding to the equations of the momentum. In a first step, a linearized semi implied method is used to integrate the first equation of the system (35) In the second step, the value (hu) i is taken as the initial condition for solving the second equation of (35).

Numerical tests
It is clear from (23) that the scheme is conservative and can compute the numerical flux corresponding to the physical solutions of water flow without relying on Riemann problem solvers. The CFL condition for the explicit scheme (22) can be written A fixed CFL = 0.9 is used and α = 2 if not specified in the test. The used computer is an Intel Core i7-8565U CPU @ 1.80GHz × 8, with 15 GB RAM.
In order to validate our FVC scheme on unstructured meshes to simulate shallow water flows, therefore we present some test cases that are proposed by several authors to validate their model and their numerical approach. The accuracy is demonstrated by comparing numerical solutions produced by the FVC scheme with analytical solutions, especially in tests 4.1, 4.2 and 4.3. To reproduce the calculation results reported in the literature, the source term of the bed is always taken into account. The C-property produced by this term has also been treated in tests 4.4 and. The Coriolis effect was taken into account in test 4.5 and this test's results are in good agreement with those presented in the literature.

Accuracy test example
We test our approach on a problem where the exact solution is known [18]. It can be readily checked that.
gives a smooth solution of the shallow water system (1) without the source terms ( i.e. ∇Z = 0, and f c = 0 ) for any choice of constants, M, a, b, x 0 , y 0 and θ. Initial and boundary condition are set according to the exact solution. We let M = 1 2 , g = 1, a = 0.04, b = 0.02, and (x 0 , y 0 ) = (−20, −10). To test the scheme ability to resolve flows that are not aligned with the computational mesh, we let θ = π 6 . We compute for (x, y) ∈ Ω = [−50, 50] × [−50, 50] up to time t = 100 s to compare the accuracy of our FVC scheme to the SRNH scheme originally introduced in [29,35]. We also propose to see the effect of the choice of the parameter α. In [8] where the authors defined for the first time the FVC scheme in its one-dimensional formulation, were able to show that the parameter α controls the accuracy of the FVC scheme, for more rigorous details we invite you to see the Lemma 3.2 and its proof. Our study reaches the same conclusion on alpha as in Fig.8 we show the convergence order for three different choices of α.  The results in Table 1 show that increasing the number of the cells in the computational domain lead to a decrease in the L 1 error for the water height h, and the discharges hu and hv in all schemes. Faster decay of the error is observed in the FVC scheme than in the SRNH scheme which is only natural a simple inspection of Table 1 also reveals that for meshes with a low number of the cells, the measured computation time is comparable for the SRNH scheme and FVC scheme. However, for meshes with a somewhat large number of cells, the FVC method is the most efficient. For example, a mesh of 161316 cells, the FVC scheme is about two faster than the SRNH scheme this is due to our FVC approach is not based on the calculation of the Jacobian matrix of the system, this matrix intervenes in many Q-scheme type approximation schemes, it is responsible for the slowness of this kind of scheme. Note that the SRNH scheme requires a solver for the Riemann problem at each time step to reconstruct the numerical flux, which is completely avoided in our FVC scheme. We can see that both schemes could reach the designed order of accuracy. In Fig.8 we have plotted the log of L 1 error calculated in Table 1 and the other values calculated for α = 0.5 and then α = 2 against the log of the maximum value of the mesh edges. We find that the L 1 errors of the FVC scheme lie on a slope line 1.4, indicating that the accuracy order of the scheme is about 1.4 for α = 0.5.

Dam-break problem
Flood flows produced by the dam break, segments of dykes, or other structures are torrential in nature with the presence of a discontinuous front propagating downstream and a rarefaction wave propagating upstream. The characteristics of these flows such as velocity, water level and time of arrival of floods must be determined in advance in order to manage floods and reduce their impact on the environment and economic infrastructure. In order to test our approach for problems related to dam break, we carried out a series of test cases proposed in the literature ( see, e.g. [8,14,39] ) etc. The proposed approach is not based on a Riemann solver technique which is very appropriate for a hyperbolic problem whose solution is often represented by a discontinuous front. Thus it will be interesting to examine the accuracy of our approach by simulating torrential flows with the presence of a discontinuity in the velocity profile and the free surface.

•) Description of the problem
We consider a rectangular channel with a flat bottom, Z(x, y) = 0 and no friction, i.e. there is no source terms, the problem is purely hyperbolic. The channel is 1.6 m long and 0.1 m wide (we assume a nondimensionalization problem), the initial conditions are given by A numerical instability is likely to occur for small values of the ratio h r /h l . At t = 0, it is assumed that the dam is abruptly removed causing a shock wave with the presence of a discontinuous front of the water surface propagating downstream. The channel is assumed to be closed on all four sides and the "slip" conditions is imposed on all walls. The computational domain is discretised by a mesh of 41776 triangles with an average size of 0.003m 2 . We will compare the water height and velocity obtained by our approach with the analytical solution which is calculated using the Stoker method [40].
The evolution of the water surface profile is used to examine the behaviour of the FVC scheme in capturing the discontinuous shock front. Henderson [23] notes that when the ratio h r /h l is greater than 0.138, the flow is sub-critical in the whole of the channel. When the ratio h r /h l is smaller than 0.138, the flow is supercritical downstream and sub-critical upstream of the dam. For very small values of h r /h l , the upstream flow regime becomes strongly supercritical, and it may be difficult to capture such a shock wave numerically.

•) Results and discussion
The first simulation concerns a river flow with h r /h l = 0.5. The Fig.9 shows the cross-section at y = 0.05 of the evolution of the water depth and the longitudinal velocity. Excellent agreement is obtained between the numerical and analytical results, this is clear from Table 2 where the L 1 error and the accuracy order of h and hu are presented respectively. The comparison shows that, under this condition, our scheme, can accurately predict the shock wave without creating oscillations.
In [16] the authors have proven that for a ratio h r /h l smaller than 0.05, most of the existing numerical models cannot give accurate results especially on the front. The last simulations, with h r /h l = 0.1 and 0.025 (see Fig. 10 and Fig.11), create supercritical flows downstream and sub-critical flows upstream. When h r /h l 0.025, a slight, non-physical oscillation occurs at the shock front. The ratio h r /h l is largely responsible for the problem of numerical instabilities that occur in the simulation of torrential flow due to dam failure. The difficulty of the problem increases with the decrease of the ratio h r /h l . We conclude that the proposed approach is well able to simulate torrential flows with a good capture of the shock front for the ratio h r /h l very small.

Tidal wave flow over an irregular bed
Here we propose to study a tidal flow on a very irregular bed. This test case has been proposed in several works to validate and to test the C-property of their approximation method [43,5]. So this test example is a good illustration of the significance of the source term treatment for practical applications to natural watercourses. It is well known that in the shallow water equations describing a flow over a very irregular background, the source terms become dominant and may cause undesirable numerical instabilities. Therefore this test case allows us to test the reliability and robustness of the proposed model when the bottom variation is fast and irregular. The bed topography is defined in Appendix (Table of bed elevation Z(x) and its illustration ). The initial and boundary conditions are constructed from the asymptotic analytical solution is given by In order to compare the numerical results of our approach with the analytical solution, we choose two results, at t = 10800 s and at t = 32400 s. In Fig.12 we present a comparison between the approach surface level and the analytical solution at t = 10800 s as well as the water height at the same time using a mesh of 200 grid points in x-direction. We also include in Fig.13 a comparison between the water velocity generated by the FVC scheme and the analytical velocity at t = 10800 s then at t = 32400 s. In Table 3 we present a comparison between the exact solution and the solution generated by FVC scheme using the relative L 1 error. An excellent agreement is obtained between the numerical and analytical solutions. This confirms that the proposed scheme is also accurate for tidal flow over an irregular bed. Moreover, these results are qualitatively in good agreement with those published in [43,5]. The accuracy of the proposed scheme in the treatment of the source terms has been identified. The numerical errors produced by the model remain very low despite the fact that the mesh is coarse and the bed is very irregular. The results of this test case confirm the good performance of our approach in the treatment of the source terms, while avoiding undesirable numerical errors due to the rapid variation of the bed.

Flow over a non-flat irregular bed
We consider the example of water flow in a two-dimensional channel including an irregular bed, a similar test has been proposed in [5]. The mathematical formulation consists of solving the shallow water system (1)  . The purpose of this test example is to verify the achievement of the C-property for the FVC scheme applied to shallow water flows over non-flat bed.  The C-property of well-balanced FVC scheme on unstructured meshes is checked in this example. As expected the water freesurface remains constant during the simulation time. The velocities and the error presented in Fig.14 also show that the equilibrium of the lake at rest is verified. All this shows that the proposed FVC scheme perfectly preserves the C-property.

Circular dam-break problem
We consider the benchmark problem proposed in [25] to study cyclone/anticyclone asymmetry in nonlinear geostrophic adjustment. We solve the shallow water system (1)  The bottom profile has the following expression: Z(x, y) = 0.3 1 + tanh 3x 2 .
Let's start by looking at the behaviour of this phenomenon in a domain with a flat bottom. The Fig.15 shows the representation of the water level calculated at different times for this test case with Z(x, y) = 0. As can be seen, a hole has formed and water is flowing out of the deepest region as a rarefaction wave progresses outwards. It is clear from the results presented that the initial elliptical mass imbalance evolves in a non-axisymmetric manner. The two expected shock waves are very well captured by the proposed FVC method. These results are qualitatively in good agreement with those published in [7,25]. In Fig.16 we exhibit the results for the velocity field corresponding to the plots Fig.15. As can be seen the two shock waves originated behind the water elevation are slowly spinning clockwise in the computational domain. The velocity field is well represented by the FVC method and re-circulation regions within the flow domain are well captured.  Let's move on to the case of the non-flat bottom to assess the performance of our FVC scheme on unstructured meshes to solve the circular dam-break problem on a non-flat bottom. The Fig.17 shows the calculated results for the water depth at t = 2 s, 8 s and 16 s using two meshes of 10040 and 40146 cells. The corresponding results for velocity field are presented in Fig.18. From a numerical point of view this test example is more difficult than the previous one as the flow is expected to exhibit complex features due to the interaction between the water surface and the bed. As in the previous test a hole has formed and the water drains from the deepest region as a rarefaction wave progresses outwards. However, a slower propagation is detected for the water free-surface in this test compared to the simulations on flat-bottom.

Conclusion and perspectives
We have presented in the first part of this paper the extension of the finite volume-characteristics scheme for solving twodimensional equations of nonlinear conservation laws on unstructured formalism. In the second part, a conservative approach has been presented to approximate the source terms while preserving the stability properties of the homogeneous solver. We note this work was the future goal of many works [6,7,2,1]. Our approach has several advantages, especially solving steady flows without large numerical errors, thus demonstrating that the proposed scheme preserves the balance related to bathymetry term; it can also compute the numerical flux corresponding to the real state of water flow without relying on Riemann problem solvers.
The method's performance has been evaluated for several test examples; furthermore, in [32], the authors concluded the robustness of the finite volume-characteristics method in other types of equations. In this work, we have developed a code for calculating free-surface flows over an irregular bottom for complex geometry. This code is based on the solution of the shallow water equations using the FVC scheme in unstructured meshes. This method applies to problems that represent large source terms due mainly to rapid variation and high background irregularity. Faced with these problems, the methods that use a Riemann approximation solver, which is well suited to the solution of purely hyperbolic equations, often encounter difficulties due to the instability of the numerical solution. In addition to the test cases presented in this paper, there is about fifteen other test cases have been carried out in the context of free-surface flows; each of these tests aims at verifying precisely one or more properties of our code. It should be noted that some extensions of this approach can be used for modeling realistic applications, notably sediment and pollutant transport. Extending this approximation method to a multi-layer shallow-water model and multi-phase flows [2,11] will also aim for future works. We can also adapt this method for solving the shallow water magnetohydrodynamic equations [33,44].