SOME MULTIPLE FLOW DIRECTION ALGORITHMS FOR OVERLAND FLOW ON GENERAL MESHES

Abstract. After recalling the most classical multiple flow direction algorithms (MFD), we establish their equivalence with a well chosen discretization of Manning–Strickler models for water flow. From this analogy, we derive a new MFD algorithm that remains valid on general, possibly non conforming meshes. We also derive a convergence theory for MFD algorithms based on the Manning–Strickler models. Numerical experiments illustrate the good behavior of the method even on distorted meshes.

(see [18,21,25]). Some attempts have been made to apply MFD algorithms on more general meshes, with an emphasis on triangular ones (see [24,27]). We refer the reader to [9] for a comparison of some of those methods, and to the references in the aforementioned papers for a more complete overview of the tremendous literature on MFD algorithms, on which we will not try to be exhaustive.
All those algorithms proceed sequentially from the higher cells to the lower ones, however as was already noticed in [23], using an MFD algorithm is in fact equivalent to solving a linear system using a specific cell ordering. The key idea underlying the present paper comes from a deeper study of the linear system associated with one of the earliest MFD algorithm. Indeed, we will explain that this linear system is completely equivalent to a two point flux finite volume scheme (see [10]) applied to a stationary Manning-Strickler model for water flow. It was then clear that replacing the TPFA scheme that requires a strong orthogonality hypothesis on meshes to remain valid by more advanced flux approximation schemes would allow us to derive MFD algorithms adapted to general meshes. Moreover, this equivalence will also allow us to derive a theoretical framework within which we will be able to study the convergence properties of MFD algorithms. In our numerical experiments, we have considered two flux reconstructions inspired by two finite volume schemes: the hybrid finite volume (see [11,12]) and the virtual volume method (or conservative first order virtual element method, see [5]). Of course, other choices could have been made (for instance VAG finite volumes [13,14] or discontinuous Galerkin methods [7]), however those two schemes will be sufficient to illustrate our approach.
The paper will be organized as follows: after describing the data and meshes, we recall the most classical multiple flow direction algorithms, and reformulate them in a more algebraic fashion. Next, using this reformulation we explain how they are linked with the TPFA scheme for a family of Manning-Strickler models. We elaborate on this basis to overcome the mesh limitations induced by the TPFA scheme, introducing a new family of MFD algorithms that will remain valid on a huge class of meshes, including those with hanging nodes. We then study the convergence properties of all methods and conclude by some numerical illustrations.

Mesh and data description
Let Ω be a bounded polyhedral connected domain of R 2 , whose boundary is denoted BΩ " ΩzΩ. We recall the usual notations describing a mesh ℳ " p , ℱq of Ω.
is a finite family of connected open disjoint polygonal subsets of Ω (the cells of the mesh), such that Ω " Y P . For any P , we denote by | | the measure of | |, by B " z the boundary of , by ℎ its diameter and by its barycenter. ℱ is a finite family of disjoint subsets of hyperplanes of R 2 included in Ω (the faces of the mesh) such that, for all P ℱ, its measure is denoted | |, its diameter ℎ and its barycenter . For any P , there exists a subset ℱ of ℱ such that B " Y Pℱ . Then, for any P ℱ, we denote by " t P | P ℱ u. Next, for all P and all P ℱ , we denote by , the unit normal vector to outward to , and , " ||´||. The set of boundary faces is denoted ℱ ext , while interior faces are denoted ℱ int . Finally for any P ℱ int , we denote " ||´|| where " t , u. We assume that there exists a subset ℱ ext,in of ℱ ext such that: where BΩ in " t P BΩ | ∇¨ą 0u and we denote of course ℱ ext,out " ℱ ext zℱ ext,in . As usual, ℎ " max P ℎ will denote the mesh size. In what follows, we will assume that our mesh satisfies: (A1) There exists a real number ą 0 and a matching simplicial submesh of ℳ such that for any P , ℎ ď where is the inradius of , and for any P and any P such that Ă , ℎ ď ℎ .
(2.1) There also exists ą 0 independent on ℎ such that for any P 1 p q: || || 2 p q ď´ℎ´1|| || 2 2 p q`ℎ ||∇ || 2 Finally, assumption (A1) implies that for any integer , there exists poly, ą 0 such that for any P and any P p q with P t1, . . . ,`1u, there exists P P p q such that |´| p q`ℎ 1{2 |´| pB q ď poly, ℎ´| | p q for P t0, . . . ,´1u . (2.3) In the estimates that will follow, the constant ą 0 will always denote by convention a quantity independent on the mesh size ℎ, whose value can change from line to line. Also by convention for any P ℱ int we will denote " t , u. In the same way, when considering a cell and one of its interior faces P ℱ X ℱ int , by convention the cell will denote the other element of .
In order to be able to establish error estimates, we assume that the topography, often called the basement in the geological community and thus usually denoted , satisfies P 2,8 pΩq and that there exists ą 0 such that´∆ ě for almost every P Ω. However, from the practical point of view is only known through some pointwise values, very often given as a digitized surface elevation, on which some interpolation and upscaling or downscaling processes have been applied (see Fig. 1). Thus, its Laplacian cannot be expected to be practically computable, and the above assumption should really be considered as an abstract technical requirement for establishing convergence estimates. Thus, our data will more realistically consists in some discrete mesh-based p q P representation of , where each can represent several pointwise values of , complemented by a water source term P 8 pΩq, possibly also only known through a mesh-based representation.

The classical multiple flow direction algorithm
In the geological literature, multiple flow direction algorithms are considered as purely algorithmic ways of distributing water from one cell to another. Thus, they are generally described in a purely algorithmic fashion. Consequently, in this section to introduce the most classical MFD algorithms we will follow the presentation generally found in the geological community, that is to say a purely algorithmic point of view. One of our first tasks will precisely consists in abstracting ourselves from this algorithmic setting. For any P , let be a value for the topography associated with cell . To fix notations, consider for the moment that but other approximations can be considered, for instance " p q for any P if is regular enough. Multiple flow direction algorithms are based on formulae to distribute the water flow from a cell to its neighbours that are topographically lower. The most classical distribution formula consists simply in distributing the flow proportionally to the ratio { of the discrete slope between the high cell and the low cell regarding the total positive slope of the high cell , where the discrete slope is given by " | | p´q and the total positive slope of cell is given by Notice that in many cases of the literature, as the MFD algorithm is applied on a uniform cartesian mesh with the same space step in every direction, the face measure | | is simply omitted (see for instance [16,22]), with no impact on the ratio { . To give a detailed description of the most classical multiple flow direction algorithms, we need to introduce some notations. Let 0 " max t | P u, 0 " t P | " 0 u and´1 " H. The set 0 thus denotes the set of cells with maximum topographic height. Then, for any P N we define the setˆof elements of by setting: ) .
By construction, as is a finite set there exists ą 0 such that´1 ‰ H andˆ´1 " . Thus, for any ě , we have " H andˆ"ˆ´1. To model water sources, we define a family p q P by setting: The classical MFD algorithm (reformulated from [16,22]) then reads as follows (see Fig. 2 for a visual illustration, where the width of the arrows is roughly following the amount of distributed water): (i) For any P , the water influx r is initialized at 0.
(ii) For any P 0 , the water influx r is given by r " | | .
(iii-1) For any Pˆ´1, we distribute the entire water influx r to the neighbours that belong to proportionally to the slope between the cell and its neighbour, i.e.: r ÐÝ r`ÿ Pℱ Xℱint, ă , P^´1 | |r p´q for all P Figure 2. Basic principle of MFD algorithms: water is distributed to lower neighbouring cells proportionally to the slope.
where is the total positive slope in cell , and ÐÝ denotes the action of updating r .
(iii-2) For any P , the water influx is complemented by the local sources by setting The MFD algorithm admits a reformulation as a linear system that will play a key role in the remaining of the paper. To our knowledge, only [23] mentioned this reformulation, although without exhibiting an explicit formula. This link seemed to have been most of the time simply overlooked by the geological community: Theorem 2.1. The MFD algorithm is equivalent to solving the following linear system for the unknown where " ÿ Pℱ Xℱint, ě | | p´q using an ordering for the cells of based on decreasing topography and a lower triangular solver.
Proof. First, notice that as cells P are updated only at step , their expression is complete past this step and given by: However, by construction of the setsˆ, any P such that ą for P , will satisfy Pˆ´1. Thus, the above expression can be simplified in: where the setˆ´1 no longer appears. As for any P , there exists 0 ď ď´1 such that P , and noticing that ą 0 as soon as there exists P ℱ X ℱ int , " t , u such that ą , (2.4) follows by induction. Finally, starting from (2.4), if one chooses an ordering for the cells of based on decreasing topography , it is clear that the above system becomes a lower triangular one for r . It should be obvious at this point that the associated lower triangular solver then coincides exactly with the classical MFD algorithm.
A great advantage of considering the above system rather than its algorithmic counterpart is that it makes clear that triangular solvers are not the only possible linear solvers, which was indeed the main point of [23]. A wide range of linear solvers can be used to solve (2.4), and most importantly parallel solvers than can considerably speed up the solving process on meshes with a huge cell number.
Remark 2.2. Many variants of this classical MFD algorithm exist, which at least to the authors knowledge mainly consist in modifying the way the influx is distributed from an upper cell to its lower neighbours: powers of the slope instead of the slope, probabilistic repartitions, repartitions using a specified number of neighbours for a given mesh structure, etc... With the exception of distribution formulae that use powers of the slope instead of the slope itself, proceeding as we have done for the most classical MFD algorithm it is not difficult to show that all those variants can be rewritten under the form: The coefficient associated with each face will take different values depending on the way one want to modify the influx repartition. Notice that we have chosen to not consider distribution formulae based on powers of the slope as they only had some technical difficulties with no fundamental differences in the analysis.

On the link with Manning-Strickler models
Consider the following classical transient Manning-Strickler model:ˇˇˇˇˇB where is the water height, a parameter and denotes the outward normal to Ω. The coefficient P 8 pΩq is the inverse of the Manning-Strickler coefficient expressing soil roughness, and is assumed to satisfy 0 ă´ď ď`ă`8 almost everywhere in Ω. Formally, the steady state associated with the above system is the following stationary Manning-Strickler model for overland flow:ˇˇˇˇ´d Denoting , " , "´∇ P 1,8 pΩq and "´∆ ě ą 0 P 8 pΩq, remark that (3.1) can be rewritten:ˇˇˇˇˇ¨∇ ,`, " in Ω , " 0 on BΩ in .
Stationary transport problems of the form (3.2) have of course received a lot of attention in the existing literature, in particular as they correspond to a popular model for neutron transport. Their well-posedness has been considered for instance in [1,2], or more recently in [8,15,17]. Thus, from the results of [8,15] and the regularity hypotheses on , , and , we know that there exists a unique P 8 pΩq solution of (3.2). As the water height only appears through its -th power , in the remaining of the paper we will with a slight abuse of notations directly use instead of . We are now going to explain that the classical MFD algorithm coincides with a well chosen discretization of the stationary Manning-Strickler model. Assume that the mesh is orthogonal, i.e. there exists a family of centroids p q P such that: and let us denote , the distance of to the hyperplane containing for any P ℱ and any P . Then, one can use a two-point finite volume scheme to discretize the steady-state Manning model. Assume that " p q or at least a second order approximation of it. Denoting for any P the discrete water height unknown, if one further assumes that " for any P ℱ ext and P which is generally what is done in practical applications of the MFD algorithm, for any P we get: ÿ Pℱ Xℱint r up p´q " | | where r up " if ě and r up " if ă and the transmissivity is given for instance by the harmonic mean: Immediately we deduce the following equivalence result, which despite its simplicity is the cornerstone of the present paper:  If we take " 1, immediately we see that (3.4) is exactly the MFD equation (2.4).
Surprisingly, this result seems to be absent of the MFD literature. The main reason is probably that formulations of MFD algorithms as linear systems are equally difficult to find. From this equivalence between the classical MFD and the two-point flux approximation (TPFA) of the classical Manning-Strickler model, some useful observations can be made. Probably the most surprising one is that the MFD unknown pr q P can be used instead of to solve the two-point approximation of the Manning-Strickler model. Moreover, existence and uniqueness of solutions of the MFD problem (3.4) are now an immediate consequence of its equivalence with the MFD algorithm without requiring any hypothesis on the topography. Next, let us mention that when modeling overland flow, the quantity of interest is not the water height but the water discharge, i.e. the norm || ℎ ∇ || of the water flux vector. This is the reason why no water height explicitly appears in MFD algorithms. However, a crucial consequence of our equivalence result is that the usual unknown r of the MFD algorithm, while an excellent choice from the algebraic perspective, is probably not the good quantity to represent the water discharge. Indeed, using r " and the consistency of the two-point formula we see that it approximates: here`" maxp0, q. We cannot expect such a r to be an approximation of the norm of ∇ , as it approximates the accumulated influx in a cell which is a mesh dependent quantity. Thus, no kind of convergence let alone approximation properties can be expected for such a quantity in general. To effectively compute a discrete water discharge for each cell P , we reconstruct cellwise the water flux vector by setting: The consistent water discharge is immediately given by " || ||. We will illustrate in the numerical section the convergence deficiency of r , and how has a much better behavior. Obviously, for any P such that ‰ 0, we can compute an equivalent positive water height by setting " r { . Cells where " 0 are cells where all discrete fluxes are ingoing fluxes, thus it is clear that one should either set "`8 or " 0 depending if water effectively reaches the cell or not. From the definition of , we see that the value of for such cells will have no influence on the water discharge and its asymptotic convergence. It is thus clear that one can always define by setting: Notice that using the harmonic mean is not mandatory if one can use an approximation of directly on the faces, leading to the transmissivity " | | { in which case (3.4) will correspond to one of the many variants of the classical MFD.

Some conservative and consistent flux reconstructions on general meshes
From the equivalence between the MFD algorithms and the TPFA scheme for the stationary Manning-Strickler model, an obvious way of generalizing MFD algorithms to general meshes consists in replacing the TPFA flux reconstruction formula by more advanced flux reconstruction techniques that are still valid on general meshes. We follow the idea of [4]: we consider a gradient reconstruction in each cell that will consist in a consistent part plus a stabilization part. However, contrary to [4] where the stabilization was used to obtain coercivity, here we will use this stabilization to enforce conservativity.

A conservative flux reconstruction formula
For any cell P , let us denote " p q the local values associated with the topography , and " R the set of local values, with the number of local values. The operator : 1 p q Þ ÝÑ will of course depend on the considered reconstruction formula, however to simplify notations we consider that the value at always belongs to the set of local values. Those local values represent all the knowledge we have in practice of the field , and once again its Laplacian cannot be expected to be computable. For typical applications such as stratigraphic modelling it consists in cell values complemented by vertex or face values, thus conditioning the schemes we can effectively use. We assume that we are given a gradient reconstruction operator ∇ : Þ ÝÑ R 2 , and we define a reconstruction formula in each cell through the operator Π : Þ ÝÑ P 1 p q, where P 1 p q is the set of first order polynomials on and: Π p q p q "`∇ p q¨p´q and thus ∇Π p q " ∇ p q. Such gradient reconstructions are the usual basis of modern finite volume methods (see [4,5,12]), and many formulae can be found in the literature. Using those reconstructions, it is tempting to define the flux operator , p q by setting: As consistency of the flux will of course rely on the consistency of the gradient reconstruction operators, the above definition will indeed lead to consistent fluxes, however to construct our generalized MFD algorithms conservativity will play a crucial role. This means that we will require the conservativity of the flux: which cannot be satisfied in general (except by the two-point fluxes) with such a simple formula. This is the reason why, following the ideas of [4], we introduce a stabilized gradient operator ∇ , :ˆR Þ ÝÑ R by setting: where ą 0 is a stabilization parameter. Then, for any P and any P ℱ , we define an intermediate flux operatorˆ, :ˆR Þ ÝÑ R by setting: Obviously, for any P ℱ int , we would like to defineˆas the solution of Fortunately, as is positive almost everywhere we immediately deduce that such aˆis always defined and is given by Thus, for any P ℱ int , it is legitimate to define the conservative flux operator , : Ś P Þ ÝÑ R by setting: whereˆis the unique solution of (4.1), while for any P ℱ ext , we simply set: It is tempting to recover conservativity using the following simpler formula forˆ, p q: The main drawback of this alternative formula can be understood in cases where is discontinuous, where this would obviously lead to a very coarse approximation of the flux along the lines of discontinuity of . Provided those discontinuities are resolved by the mesh, our formula forˆ, p q will instead behave as the usual harmonic mean. It is thus a more robust and versatile choice.

Consistency properties of the flux
Denoting p , q the ball centered at of radius and Ω p , q " p , q X Ω, we recall the usual definition of strong consistency for gradient reconstruction operators: ą 0 independent on ℎ such that for any P 2´Ω p , q¯, every 0 ă ď 2ℎ and every P Ω p , q: As an immediate consequence of the above definition, we have, using Taylor's expansion and the density of 8´Ω p , q¯in 2,8 p Ω p , qq that for any P 2,8 pΩq there exists ą 0 depending on but not on ℎ such that for every ℎ ď ď 2ℎ , any P and almost every P Ω p , q: Another immediate consequence is the consistency of our discrete fluxes: Assume that assumption (A1) is satisfied and that the gradient reconstruction operator is strongly consistent. Then, there exists ą 0 independent on ℎ such that for any P , any P ℱ , any P 2,8 pΩq and any P 1,8 pΩq, we have:ˇˇˇˇˇˇˇ1 Proof. By density of 8`Ω˘i n 1,8 pΩq and 2,8 pΩq, it is clear that it suffices to establish the result for P 8`Ω˘a nd P 8`Ω˘. Assuming such regularity, we start by establishing that there exists ą 0 independent on ℎ such that for any P ℱ int : |ˆ´p q | ď ℎ 2 || || 2,8 pΩq and |ˆ´Π p q p q | ď ℎ 2 || || 2,8 pΩq @ P . As by construction, we have Π p q p q "`∇ p q¨p´q, slightly rewriting (4.2) we get:´p Immediately, using the strong consistency of the gradient operator we get that for all P : |Π p q p q´p q | ď ℎ 2 || || 2,8 p Ω p ,ℎ qq and thus for all P : as under assumption (A1), we know (see [3,7]) that there exists ą 0 such that for any p , q P 2 such that ℱ X ℱ ‰ H, we have ℎ {ℎ ď . Moreover, as the gradient operator is strongly consistent, we know that: Meanwhile, we clearly get using again assumption (A1): hich finally establishes that: |ˆ´p q | ď ℎ 2 || || 2,8 pΩq and (4.6) immediately follows using the triangular inequality and the strong consistency of the gradient operator. Then, for any P Ω p , 2ℎ q and any P ℱ int : 1 | | , p q`p q∇ p q¨, " p∇ p q´∇ p qq¨,`p p q´q ∇ p q¨, ,´ˆ´Π p q p q¯.
Immediately, as the gradient reconstruction operator is strongly consistent we get: while Taylor's expansion immediately gives: Under assumption (A1), we have ℎ { , ď ℎ { ď 1{ which combined with (4.6) gives the desired result for P ℱ int . For P ℱ ext , we directly get: 1 | | , p q`p q∇ p q¨, " p∇ p q´∇ p qq¨,`p p q´q ∇ p q¨, and thus applying the same estimates than in the case of interior faces concludes the proof of (4.5).

The generalized MFD algorithm
Given an approximation p q P of the topography and p q P of the source term as data, and using both upwinding for the water height and our discrete fluxes, the most natural discrete formulation of the Manning-Strickler model consists in finding´p q P , p q Pℱ ext,in¯s uch that:ˇˇˇˇˇˇÿ where the upwinding formula up for any set of cell values p q P is defined by: and where we have defined the influx boundary associated with the discrete fluxes by setting ℱ ext,in " t P ℱ ext | , p q ă 0, " t u u .
By construction of the fluxes, for any P ℱ int we have that , p q "´, p q. Mimicking what we have done for the TPFA scheme, we gather faces by upwinding kind and use the fact that up " 0 for all P ℱ ext,in , leading to:¨ÿ Pℱ , , p qě0 , p q‚´ÿ Pℱ Xℱint, , p qą0 , p q " | | .
Defining: " ÿ Pℱ , , p qě0 , p q and r " and noticing that ą 0 as soon as there exists P ℱ such that , p q ą 0, we see that the above equation can be rewritten: For any P ℱ ext,in , we set: , p q and r " .
Notice that the above system is set for the unknowns´pr q P , pr q Pℱ ext,in¯, thus it has the same number of unknowns than the original system (5.1). We still reconstruct a water flux vector cellwise by setting: the consistent water discharge being given by " || ||. Finally, for any cell P , we set as we have done for the TPFA case: while for any P ℱ ext,in , " 0. Notice that this new system allows to correctly define the numerical method, while the system (5.1) does not uniquely define the water height on cells where " 0. Contrary to system (3.4), system (5.2) admits no obvious renumbering of mesh elements that makes the system triangular. However, as it is anyway necessary to enable parallelism, we say that the generalized MFD algorithm will consists in solving the linear system (5.2) using any linear solver of the literature. Notice that as we are using the unknown r rather than , the discrete system (5.2) takes the form of a perturbation of the identity matrix. Thus, we expect its condition number to be relatively good (depending of course of the coefficient ) and in any case much better than if had tried to solve (5.2) directly.

On existence and uniqueness for the generalized MFD algorithm
Existence and uniqueness of a solution to the above system are much less obvious than in the case of the TPFA scheme. To explain this fact, let us denote " t P | ą 0u .
Using the definition of˚, summing (5.2) over P , after a straightforward manipulation of the last sum we get: Thus, the existence of a solution for any second member p q P requires that: In the case of the TPFA scheme, condition (5.6) is necessarily satisfied. If condition (5.6) is not satisfied, then (5.2) can have a solution only if ÿ P | | " 0.
From this observation, it seems clear that establishing an exhaustive existence and uniqueness theory for system (5.2) is a more complex task than what one could expect. For the continuous problem, well-posedness is directly linked to the properties of the topography and in particular to its Laplacian. Unsurprisingly, well-posedness of (5.2) will be directly linked to the properties of the quantity that plays the role of the topography at the discrete level, that is to say the flux family p , p qq P , Pℱ . Without any further assumption on the mesh, far from the asymptotic regime consistency alone cannot be expected to ensure well-posedness for coarse meshes. This means that well-posedness will in general be controlled by the structure of the discrete flux family rather than by the properties of . As in practice the Laplacian of cannot be computed in general, it is anyway better to have an existence result that uses only computable quantities. In fact, not only the necessary condition (5.6) must be satisfied, but it also requires that there does not exist what we will call a flux cycle. For any subset of , we denote: We say that a set of cells Ă form a flux cycle if and only if for any P , p q ą 0 and ÿ Pℱ Xℱ int , , p qă0 , p q ă 0.
The first condition implies that water can enter the cycle but not leave it, while the second condition implies that water will run through the cycle without stopping anywhere. Whether such configurations can exist on coarse meshes for consistent flux families is a difficult question. However, from the physical point of view flux cycles are completely unrealistic as they represent regions where water will at some point go from a topographically low region to a topographically higher one. Moreover, under the positivity hypothesis on the Laplacian of the topography, we know that for the exact fluxes we have: and thus no flux cycle can exist. Thus, the presence of flux cycles should be considered as anomalous and a sign that a too coarse mesh has been used. A natural sufficient condition of existence and uniqueness for (5.2) should then be one that is satisfied by the continuous fluxes and immediately implies that no flux cycle exist.
Proposition 5.1. We say that there exists a flowing path from P z to if there exists r P such that D P N, ą 0 and p q 0ď ď´1 such that P ℱ int @0 ď ď´1 and " t ,`1u @0 ď ď´1 where " 0 r " and , p q ą 0.
If ‰ H and for any cell P z there exists a flowing path to , then (5.2) is well-posed. Proof. As system (5.2) is linear, it suffices to establish uniqueness to ensure the existence of solutions. Thus, assume that r "´pr q P , pr q Pℱ ext,in¯i s solution of (5.2) with zero right-hand side. Then, taking the modulus of (5.2) and summing over P we get: |r | , p q ď 0 which immediately implies that r " 0 for any P . Next, let us denote 0 " and 0 " , and define for any ě 1 by setting: " z ď 0ď ď´1 where " p z˚q Y˚ẘ ith˚" t P | ą 0u and˚˚" P˚| P ℱ X ℱ ext | , p q ą 0 ( ‰ H ( and for any ě 1: ℱ " t P ℱ | X ‰ Hu and ℱ int " t P ℱ X ℱ int | Ă u and ℱ ext " ℱ zℱ int .
If " H, then there is nothing to prove. Otherwise, as for any P , there exists a flow path to , then˚‰ H. Thus, summing over and proceeding as above we obtain that: Thus, we obtain that r " 0 for any P and ‰ H which implies`1 Ĺ . Then, it is clear that the sequence p q ě0 is strictly decreasing in the sense that it satisfies Ĺ´1 or " H because of the flow path assumption. Thus, there exists ą 0 such that " H and thus r " 0, which concludes the proof.

Convergence analysis
The main idea underlying our convergence analysis is that for all consistent fluxes, if the topography is regular enough numerical fluxes will converge to the continuous ones. The major originality of the present analysis regarding for instance finite volumes or discontinuous Galerkin methods for steady-advection diffusion problems of the form (3.2) is that here the coefficients and are approximated through discrete differential operators applied to the topography . Thus, this additional approximation has to be handled carefully to recover the correct convergence estimates.
To establish precise error estimates we first need to choose a discrete norm. To this end, we define , " 1 2`| , p q |`| | || ∇¨, || 8 p p ,2ℎ qqȃ nd set for " p q P : In the above expression, one would more naturally expect to find , p q instead of , . However, their behavior is difficult to control and in particular, it is complex to estimate the way this fluxes stay away from zero when the continuous fluxes ş ∇¨, cancel creating technical difficulties, which , avoids. Immediately, we get that: | ,´| , p q | | " 1 2ˇˇ| , p q´| | || ∇¨, || 8 p p ,2ℎ qqǎ nd thus: where for any P and any P ℱ we have denoted: , " maxˆˇˇˇˇ, p q`ż ∇¨,ˇˇˇˇ,ˇˇ| , p q |´| | || ∇¨, || 8 p p ,2ℎ qqˇ˙.
Notice that from (4.5), we get that there exists ą 0 such that for any P and any P ℱ , we have: The main objective of this section is to establish the following explicit 2 error estimate: ith , defined by (6.2). Then, there exists ℐ ℎ p q P 2 pΩq such that for any P , ℐ ℎ p q| " ℐ p q is constant on and: |´ℐ p q| 2 2 p q`ℎ 1{2 |´ℐ p q| 2 2 pB q ď poly ℎ | | 1 p q . Moreover if pℎq Ñ 0 when ℎ Ñ 0, then there exists ℎ 0 ą 0 and ą 0 depending on , , and the mesh parameters such that for any ℎ ď ℎ 0 : where is reconstructed by (5.4) from the solution of (5.2).
Notice that the existence of ℐ ℎ p q in the above result is an immediate consequence of assumption (A1) and (2.3). As a corollary, we will then be able to establish an explicit error estimate for the water discharge: Corollary 6.2. Under the assumptions of Proposition 6.1, there exists ℎ 0 ą 0 and ą 0 depending on , , and the mesh parameters such that for any ℎ ď ℎ 0 : The error estimates use the surprising factor pℎq instead of ℎ as one could legitimately expect. This comes from the fact that we use approximate fluxes, which is equivalent in some sense to using a quadrature formula for the coefficients of (3.2). As the natural norm for the discrete problem uses those approximate fluxes, it seems unfortunately unavoidable that this impacts the error estimate. In the case where for any P and any P ℱ that is indeed present in the norm ||¨|| ℎ , the , 's are bounded by below by a constant independent on ℎ, which will happen most of the time in practice, then pℎq will behave as ℎ and we retrieve a convergence at rate ℎ 1{2 . Controlling the lower bound for , is the reason why we have used a supremum on a set containing more than in the definition of , . If this supremum is zero, most gradient reconstruction operators and in particular those described in our numerical experiments will provide a discrete gradient that aligns with the continuous one and the discrete flux will be zero too, avoiding most of the problematic cases for the asymptotic behavior of , . In the case of exact fluxes it is known that for the water height the optimal convergence rate to is ℎ 1{2 , even if superconvergence at rate ℎ is very often observed in practice, as revealed by [6,19] (see also [7]). Thus, as additional flux consistency terms could only decrease the convergence rate, it seems clear that our estimate for the water height cannot be expected to be improved in terms of order of convergence.
The proof of Theorem 6.1 being rather lengthy, we will try to decompose it as much as possible. Let us notice that in the special case where ∇ is constant, the discrete flux are exact. Then our problem corresponds to a piecewise constant version of the discontinuous Galerkin method described in [20], and we can in principle follow the steps of their proof. Denoting || || 2 ℎ " ř P || || 2 ,ℎ with obvious notations, the first step of their proof is to establish a local error identity for ||´|| ,ℎ for any in R, using the exactness of the flux. Then, taking " ℐ p q and summing over P , they obtain a global error estimate where the residual terms only involve´ℐp q, which can be straightforwardly estimated using polynomial approximation results.
The general guideline of our proof is the same, however we have to endure some additional technicalities due to the non exactness of the flux. We first start by establishing the equivalent of the local error identity of [20], but keeping our approximate flux in the result. Thus, an additional residual term standing for the consistency error of the flux will appear on the right hand side of our identity. Then, to match the definition of the ||¨|| ℎ norm, each time our estimates will involve the sum of the flux over the full set ℱ , we will replace it by the Laplacian of , and put the difference in a residual term. For an isolated flux term, we will do the same using this time either , if the term is involved in the ||¨|| ℎ norm or by the continuous flux if it is directly a residual term. Thus, when summing the local error identities over P to obtain the global error estimate, residual terms involving the discrete counterpart of the Laplacian will be handled through a technical result described in the following subsection, while other residual terms will be estimated directly using the strong consistency of the flux. The reason for doing so is to obtain an error estimate involving only the ||¨|| ℎ norm on the left-hand side and residual terms involving either´ℐp q or flux consistency error terms on the right-hand side. Finally, each residual term will be estimated using polynomial approximation results and the consistency of the flux.

Local error identity and approximation results for the Laplacian
Let us now establish an abstract local error identity, following the lines of [20]. Lemma 6.3. For any P , any P R and any " p q Pℱ P 2 pB q constant on each P ℱ , we have: where is the solution of the continuous problem (3.1).
In the above local error identity, one can see that the left hand side strongly looks like the -term of the ||¨|| ℎ norm. However, to match the ||¨|| ℎ norm, we need to replace the sum of the fluxes over ℱ by the Laplacian of , as well as isolated flux by , . If this last operation straightforwardly leads to a residual term because of (6.1), however the following lemma precise what kind of link we can expect between the discrete Laplacian and the true Laplacian: Lemma 6.4. For any p q P , one has: Proof. Let us first notice that: Indeed, using , p q`, p q " 0 for any P ℱ int , " t , u we have , p q up " 0 while for P ℱ ext such that , p q ą 0, " t u, up " 0 immediately implies that: , p q up " 0.
Thus, we have: , p q using the definition of up . In the same way: Combining those expressions gives the expected result.

Global error estimate
From the local error identity, and the link between the discrete and continuous Laplacians, we deduce a global estimate for the discrete norm of the error: Lemma 6.5. Using the hypothesis and notations of Proposition 6.1, we have: , p q p´ℐ up p qq pp up´ℐ up p qq´p´ℐ p qqq , p q p´ℐ up p qq pp´ℐ p qq´p up´ℐ up p qqq p| , p q |´, q pp´ℐ p qq´p up´ℐ up p qqq p| , p q |´, q |´ℐ p q| 2 .
Proof. We start from (6.6) with " ℐ p q and " ℐ up p q. Summing over P the second term of the left hand side of (6.6) leads to: , p q p up´ℐ up p qq , p q p´ℐ p qq 2 as both and ℐ p q cancels for P ℱ ext,in . Thus, summing over P the first two terms of the left hand side of (6.6), we obtain: , p q pℎ ,´ℐ p qq , p q p up´ℐ up p qq , p q pℎ ,´ℐ p qq 2 .
Consequently, summing over P the entire left hand side of (6.6) leads to: , p q pp´ℐ p qq´p up´ℐ up p qqq , p q |´ℐ p q| 2 " 1`2`3 .
The first term of the above expression rewrites: while the second term and third term give: , pp´ℐ p qq´p up´ℐ up p qqq p| , p q |´, q |´ℐ p q| 2 " 2`3`4`5 with obvious notations. Immediately, using the hypothesis on and we have: Then, noticing that the second member of the sum over P of (6.6) can be decomposed into three terms 1`2`3 with obvious notations and combining the above results, we have that ř P ℐ " 1`2`3 is equivalent to ||´ℐ ℎ p q|| 2 ℎ ď 1`2`3´1´4´5 as 0`2`3 exactly gives the square of the ||¨|| ℎ norm. Using the fact that both and ℐ p q cancels for P ℱ ext,in , we obtain: , p q p´ℐ up p qq pp up´ℐ up p qq´p´ℐ p qqq , p q p´ℐ p qq p´ℐ p qq , p q p´ℐ up p qq pp´ℐ p qq´p up´ℐ up p qqq which concludes the proof.

Estimation of the residual terms
Establishing Proposition 6.1 now just consists in estimating each term appearing in the second member of the above global estimate. From Lemma 6.5, we see that ||´ℐ ℎ p q|| 2 ℎ ď with obvious notations. The first three terms account for the polynomial approximation error, while the four other terms account for the consistency error of the flux. We now estimate those two families of residual terms separately.
Proof of Proposition 6.1: terms accounting for polynomial approximation error. Using Cauchy-Schwarz inequality and the fact that p up´ℐ up p qq´p´ℐ p qq is constant on gives: As the gradient operator underlying the fluxes is consistent, we know from Proposition 4.5 that there exists ą 0 such that:ˇˇˇˇ1 As soon as ℎ ď ℎ 0 for some fixed ℎ 0 this leads to the existence of some ą 0 independent on ℎ such that Then, using the definition of ℐ ℎ p q, we get ż |´ℐ up p q| 2 ď ℎ up | | 2 if , p q ě 0 and up " if , p q ă 0, and up " if P ℱ ext . From [7], we know that assumption (A1) implies that card pℱ q is bounded by a constant ℱ depending on the mesh parameters but independent on ℎ, and thus there exists ą 0 such that | 1 | 2 ď 2ℎ || || 1,8 pΩq || || 2,8 | | 2 1 pΩq ||´ℐ ℎ p q|| 2 ℎ . For 2 , obviously we have using Cauchy-Schwarz inequality: and thus, proceeding the same way than for 1 , we obtain | 2 | 2 ď 2ℎ || || 1,8 pΩq || || 2,8 | | 2 1 pΩq ||´ℐ ℎ p q|| 2 ℎ . The situation is more delicate for 3 . Indeed, ℐ up p q " 0 for P ℱ ext,in , while " 0 on P ℱ ext,in , thus ℐ up p q is not directly a good approximation of . For any P ℱ ext,in Xℱ ext,in , ℐ up p q " " 0 and thus the contribution to 3 is zero. Then, it just remains to consider the case where P ℱ ext,in and R ℱ ext,in . In this case we havé ş ∇¨, ě 0 while , p q ă 0. However, by definition of the residualˇˇ, p q`ş ∇¨,ˇˇď , which immediately implies thatˇˇş ∇¨,ˇˇď , and | , p q| ď , . Consequently, we get using Cauchy-Schwarz inequality: .
Using the definition of pℎq and the trace inequality on Ω, this leads to: It finally remains to estimate the terms accounting for the consistency error on the fluxes.
Proof of Proposition 6.1: terms accounting for flux consistency error. For 4 , let us first remark that: | | , p q`∇¨,˙pp´ℐ p qq´p up´ℐ up p qqq Using the fact that both and ℐ p q cancel for P ℱ X ℱ ext and , p q ă 0, this rewrites: Cauchy-Schwarz inequality then leads to: ,´p´ℐ p qq´p up´ℐ up p qq 2¯‚ 1{2 .
Using (6.3), (2.2) and the definition of pℎq, we easily get: In the same way, Cauchy-Schwarz inequality and this time the trace inequality on Ω directly gives: From Lemma 6.4, we deduce that: ,´p´ℐ p qq 2´p up´ℐ up p qq 2ÿ P ÿ Pℱ Xℱext, , p qą0 , p´ℐ p qq 2 .
Using p 2´2 q " p´qp`q and Cauchy-Schwarz inequality, we obtain: .

Error estimate for the water discharge
To establish Corollary 6.2, we will need an auxiliary discrete stability result: Lemma 6.6. The following estimates holds: Proof. To establish (6.7), using the definition of up , we obtain multiplying (5.1) by and summing over P : Remark that up " 1 2´2`2¯´1 2 p´u p q 2 which immediately leads to , p q p´u p q 2 .
Using the fact that: , p q 2 we immediately get (6.7) as: Proof of Corollary 6.2. Let us denote p q 0ď ď1 the canonical basis of R 2 . Assume that ℎ ď ℎ 0 , where ℎ 0 is defined in Proposition 6.1 and consider: We begin by establishing that ď maxpℎ, pℎqq 1{2 || || 1 pΩq . (6.8) Indeed, we have: Poincare-Wirtinger's inequality applied component by component (see [7]) gives: where ą 0 is independent on ℎ under assumption (A1), and (6.8) follows from Proposition 6.1 and the above estimate. To conclude, it suffices to remark that by definition: Thus applying Cauchy-Schwarz inequality, we get: sing the bound on the flux from the proof of Proposition 6.1, and where: , p q p´u p q 2 .

Numerical exploration
To assess the asymptotic behavior of the method, we use a very simple analytical solution on which the positivity of the Laplacian is guaranteed (contrary to typical geological applications). Let us consider the square domain Ω "s0, rˆs0, r associated with the topography p , q " 0´´p´0 q 2`p´0 q 2¯w here 0 " 0 " {2, for which´∆ " 4 . Consider the water height defined by p , q " 0 exp´´´p´0q 2`p´0 q 2¯¯.
In practice we take 0 " 1, " 2, " 1{2, 0 " 2 and " 1. We consider seven types of mesh sequences. The first sequence consists in classical Delaunay meshes (2dDelaunay). The second one (2dDualDelaunay) is obtained by considering the dual meshes of a sequence of Delaunay meshes (Fig. 3). The third sequence (2dVoronoi ) is made of Voronoi meshes, possessing the mesh orthogonality property. The fourth sequence (2dKershawBox ) is a sequence of Kershaw meshes, while the fifth one (2dCheckerBoardBox ) is a sequence of checkerboard meshes. These two sequences have only quadrangular cells which are distorted for the sequence 2dKershawBox, while the sequence 2dCheckerBoardBox allows to test the behavior of the method in presence of non conformities. The sixth sequence, named 2dSquareCart, is simply a uniform cartesian mesh with square cells, while the seventh one named 2dRectCart is a uniform cartesian mesh with rectangular cells. These two last sequences are mainly intended to serve as reference. For the generalized MFD algorithm, we consider two consistent gradient reconstructions coming from finite volumes: the first one is the hybrid finite volume gradient operator of [12], that uses faces values for the topography , while the second one is the vertex based finite volume gradient operator (VVM) of [5].
As it is the quantity of practical interest, we will express our convergence results directly in terms of the water discharge instead of using the water height. Let us begin by some results for sequences 2dSquareCart, 2dRect-Cart and 2dVoronoi. Being the only ones that satisfy the mesh orthogonality requirement (see Figs. 4 and 5), we expect that the TPFA scheme will also converge on those sequences. The corresponding results are displayed on Figures 6 and 7, while the associated approximate convergence orders for each scheme are given in Table 1.  Clearly, on the basic cartesian cases 2dSquareCart and 2dRectCart, all the schemes give identical results. A closer look at our generalized flux formula immediately reveals that this is perfectly normal as the flux , degenerate into the two-point flux formula on those cartesian meshes where the barycenter of each internal face is exactly the intersection point between the face and the segment joining the cell centers on each side of the face (see [12]). Moreover, Table 1 reveals that all the methods are superconvergent on the three orthogonal mesh sequences. In [20] it is indeed established in the case of exact fluxes that the approximation of the water height superconverges at least on rectangular meshes, which explains that this observed superconvergence is probably not anomalous. Next, we turn to the other mesh sequences, for which we cannot expect convergence for the TPFA scheme. The corresponding results are displayed on Figures 8 and 9, while the associated approximate convergence orders for each scheme are given in Table 2.
As expected, the TPFA scheme does not converge anymore, however the generalized MFD algorithm is converging for both the hybrid and virtual volume (VVM) gradients. The results for those two variants are in fact relatively close to each other. In particular, the behavior of the method on sequence 2dCheckerBoardBox confirms its ability to deal with non conformities and thus local mesh refinement. Remark that superconvergence       Figure 10. Comparison between influx and water discharge for the TPFA scheme on orthogonal meshes (exact water discharge is bottom right).
is also achieved again on those test cases. To conclude, on Figure 10, we represent on the bottom-right the correct water discharge, while the top-left figure is the water influx r for sequence 2dSquareCart, the top-right figure is the water influx for sequence 2dRectCart and the bottom-left figure is the water influx for sequence 2dVoronoi. These very basic examples illustrate our point on the water influx r : it is clearly a mesh dependent quantity, that cannot be expected to reproduce the behavior of the discharge correctly, even on those very simple cases. If one does not consider the Manning-Strickler framework, one of the tricky features of r is that on very common cartesian mesh sequences such as 2dSquareCart and 2dRectCart it has a relatively nice qualitative behavior. It even seems to converge, although to a mesh sequence dependent quantity. This explains why it has been used inadvertently in the hydrogeology community, despite of its mesh dependent nature: on the widely used cartesian meshes, it has an acceptable behavior that cannot be easily discarded without a reference continuous model, in particular for complex topographies. However, on more general meshes, the situation is very different, even for the TPFA scheme as the case of sequence 2dVoronoi reveals. We display on Figure 11 the water influx r for the four other mesh sequences, in the case of the hybrid gradient. It is completely obvious on those sequences that r should not be considered as a physically relevant quantity.

Conclusion
We have established the equivalence between a classical family of multiple flow direction algorithms and the TPFA scheme applied to a family of stationary Manning-Strickler models. From this equivalence we have proposed an improvement of the derivation of the discrete water discharge based on a consistent flux reconstruction Figure 11. Water influx for hybrid gradient on meshes 2dDelaunay, 2dDualDelaunay, 2dKer-shawBox and 2dCheckerBoardBox.
rather than directly using the water influx intermediate unknown that is not only non convergent in general, but does not approximate the water discharge in any case. Then, using more advanced flux reconstruction schemes, we have proposed a multiple flow direction algorithm adapted to general meshes. A convergence theory that covers both cases was developed, and numerical experiments illustrate the good behavior of the method. Remark that the extension to the case of Manning-Strickler tensors rather than coefficients is straightforward, as would be the extension of the present analysis to multiple flow direction models that use powers of the topographic slope.