A Saint-Venant shallow water model for overland flows with precipitation and recharge

We propose a one-dimensional Saint-Venant (also known as open channel or shallow water) equation model for overland flows including a water input--output source term. We derive the model from the two-dimensional Navier--Stokes equations under the shallow water assumption, with boundary conditions including recharge via ground infiltration and runoff. We show that the energy-consistency of the resulting Saint-Venant model is strictly dependent upon the assumed level of rain- or recharge-induced friction. The proposed model extends most extant models by adding more scope depending on friction terms that depend on the rate of water entering or exiting the flow via recharge and infiltration. The obtained entropy relation for our model validate it both mathematically and physically. We compare both models computationally based on a kinetic finite volume scheme; in particular, we provide numerical evidence that the two models may show drastically different results, where the model conditioned on the flow velocity provides what may not be the physically relevant solution because of the lack of the appropriate friction terms depending on a master parameter. We also look at a comparison with real-life data.


Introduction
The quantitative modelling of hydrology of catchment basins and rivers holds a central place in environmental sciences, particularly in connection with water availability, flood risks, and urban sewer systems. This is especially important today in understanding and forecasting the impact of climate variability on the human and natural environment. A watercourse's recharge via precipitation, runoff or springs, as well as its loss via ground infiltration, drainage systems, vegetation, and evaporation, play a essential role in quantifying its dynamics. The most important component of the hydrologic water recharge and loss are the infiltration and precipitation processes. The forecast to predict the motion of water is a difficult task to which substantial effort has been devoted [Grace and Eagleson, 1966, Woolhiser and Liggett, 1967, Zhang and Cundy, 1989, Esteves et al., 2000, Weill et al., 2009, Rousseau et al., 2012. One of the most widely used models to describe the overland motion of watercourses is the classical one-dimensional (1D) Saint-Venant shallow water equation (SWE) (also known as open channel equations). The classical SWE, which is a hyperbolic system of 2 or 3 scalar conservation laws, is most commonly used in its massconservative form, where the addition or subtraction of water occurs exclusively through the boundary conditions. It is important, however, to understand also the addition of water from the "inside" of the domain. In the three-dimensional reality or two-dimensional configurations, the additional water comes from underground water (aquifers, springs, water-table, etc.) or sinks and runoff phenomena (torrential tributaries, surface flow, quick-flow) and direct rain; however, for d = 2 or 3, the reduction of the d dimensional Navier-Stokes to the d − 1 dimensional SWE makes some of these boundary terms, such as the underground springs or the soil's absorption, into "internal source terms" S [Gerbeau and Perthame, 2001, Marche, 2007, Ersoy, 2015. This source term quantifies the amount of water that is added to (S > 0) or subtracted from (S < 0) the flow; in practice this may occur with a variety of mechanisms. Thus, according to Sochala [2008] and Delestre et al. [2012] the source term is added to the SWE resulting in where the unknowns h(t, x) and u(t, x) model, respectively, the height of the water and the velocity of the water (column) at a time-space point (t, x), g the gravitational acceleration (considered a constant g ≈ 9.81 m/s 2 ), Z(x) the topography of the river bed with slope ∂ x Z(x), and k 0 an empirical fluid-wall friction. The reader interested in such questions could consult, e.g., Ponce and Simons [1977], Akan and Yen [1981], Moussa and Bocquillon [2000], Singh [2001] as further references. Our main goal in this paper is to derive, starting from the Navier-Stokes equations with a permeable Navier boundary condition to account for the infiltration and a kinematic one to consider the precipitation, a model akin to (1.1) via vertical averaging under the shallow water assumption. The averaged model that we obtain extends model (1.1), in that it has an additional momentum source term of the form Su − (f R +f I )u, where the total recharge is S := R − I (1.2) with R 0 denoting the recharge on the free surface (briefly called rain term but accounting also for run-off and minor torrential tributaries, ultimately coming from rainfall) and −I the recharge rate due to infiltration from water to ground (I < 0) or recharge from the ground into the water (I > 0); in §2 we give a more detailed discussion about R and I. The terms f R u andf I u that we subtract from the momentum's rate model the friction caused by the addition of water (with velocity 0) attaching to and being advected by the flow. For simplicity, in this paper we assume the most basic constitutive relations for this friction: linear in R for f R and piecewise linear in I forf I . Explicitly, the new model we propose is which is a hyperbolic system of balance laws. We show that system (1.3) possesses a mathematical entropy given by which satisfies the following entropy relation for smooth solutions: We remark that the functionÊ is different from its analogue when S ≡ 0: indeed, if this is the case, replacingÊ(h, u) withĚ(h, u, Z) :=Ê(h, u)+g hZ (the "standard" entropy) would lead to base height invariant balance law, whereas in the case of the balance law, this is not possible withĚ.
In equation (1.3) the friction coefficients due to rainfall recharge, infiltration and kinematics respectively, are f R := αR,f I := α max(0, −I), and k 0 (u) (1.6) with α a coefficient which one should be able to capture from empirical measurements. We also introduce the total head function which will be used to derive an entropy condition. Two important remarks are worth making: (1) The friction terms f R andf I are necessary to avoid paradoxical outcomes such as perpetual motion, their physical interpretation being, roughly speaking, that incoming water, assumed to have horizontal velocity zero, must "stick" to the flow and be transported at velocity u.
(2) Our model (1.3) generalises model (1.1), which is the special case when the friction terms match the extra term S = f R +f I , meaning α = 1 and I < 0 (water enters, but does not exit, the flow from the ground).
(3) Our model (1.3) could be further generalised by assuming friction relations more general than ours, for example, by having two separate friction coefficients instead of the single α, possibly accounting for outgoing fluid. Also the linearity of f R in R andf I in I < 0 could be replaced by more precise constitutive relations, possibly obtained from empirical data.
Connected to these remarks, in our main result, discussed §3.4, we deduce that the energetic-consistency of the model is strictly dependent upon the level of rain-or recharge-induced friction, denoted by α; that is, Ru 2 − min(0, I)u 2 (1.8) which is pertinent with the underlying physics. If the terms Su, f R , andf I are dropped we recover the equations (1.1), for which we show that the above properties are conditional on the water height and may exhibit non-physical solutions. We outline the rest of the article as follows: in §2, as our starting point we present the Navier-Stokes equations and the boundary conditions including recharge, infiltration and corresponding friction terms. In §3 we derive the consequent Saint-Venant equations with recharge and infiltration, including the friction terms. In §4 we adapt a finite volume kinetic scheme of Audusse et al. [2000] and Perthame and Simeoni [2001] to our model and provide extensive numerical testing in §5 of the resulting code. A C and C++ implementation of this code, written by Matthieu Besson, Omar Lakkis and Philip Townsend, is freely available on request (an older version is given by Besson and Lakkis [2013]).

Navier-Stokes equations with infiltration and recharge
Our aim is to construct a mathematical model for overland flows that is consistent with the physical phenomena that can affect the motion of such water. To this purpose, we propose a model reduction of the two-dimensional Navier-Stokes equations leading to an extension of the standard Saint-Venant system. By considering suitably chosen boundary conditions, we take into account the addition and removal of water, either by rainfall (e.g. from runoff onto the top of the water course) or by ground-water infiltration of exfiltration processes (e.g via a porous soil). We start in §2.1 by reviewing the Navier-Stokes equations in the special geometric setting, describing the physics with a wet boundary on the bottom of the water course and a free surface on the top. We then introduce the boundary conditions for each surface in §2.2 and §2.3, respectively.
2.1. Geometric set-up and the two-dimensional Navier-Stokes equations. From a modelling point of view there is no upper bound on the time, but with numerical and practical applications in mind, we will work with an arbitrary final time T > 0. With reference to Fig. 1, we consider an incompressible fluid moving in the time-space box where L > 0 is the horizontal length of the domain. (2.1.1) The absolute height of the surface of the water course and the topography of the bed are modelled, respectively, by the functions whose values measure with respect to a reference horizontal height 0. We define the local height of the water by The wet region is defined as the region in which the fluid resides at each time t ∈ [0, T ] Ω(t) := (x, z) ∈ R 2 : x ∈ (0, L) , Z(x) < z < H(t, x) (2.1.4) water column water's free surface: Figure 1. Schematic setting of the model, color coded with a water color, ground color, and flow color, where by "flow" we mean any source of variation of the water quantity, including the boundary fluxes (indicated with arrows, but in fact scalar quantities).
with its global counterpart (2.1.5) We assume that the viscous flow u satisfies, on the space-time domain Ω, the two-dimensional incompressible Navier-Stokes equation where u = (u, v) is the velocity field, ρ 0 is the density of the fluid (taken to be constant), F = (0, − g) is the external force of gravity with constant g, and σ [u] is the total stress tensor whose matrix given by where p is the pressure and µ > 0 the dynamic viscosity. The (algebraic) tensor product of two vectors a ⊗ b is defined as ab (all vectors are displayed as columns) and the div of a covector/tensor is taken as the row-wise divergence of the associated matrix; in coordinates this means To work with the wet region, we introduce its indicator function x)] for all t, x, z ∈ R.
(2.1.9) with the Iverson notation The function Φ is advected by the flow so its material derivative, with respect to the flow u, must therefore be zero. Moreover, thanks to the incompressibility condition, Φ satisfies the following indicator transport equation (2.1.11) 2.2. The wet boundary. Crucial to our model derivation is the particular situation at the wet boundary, where the effect of infiltration plays a central role. The wet boundary is the set of points (t, x, Z(x)), for t > 0 and 0 < x < L and for which H(t, x) − Z(x) > 0. Since the topography is assumed to be rough it produces friction, and due to its porosity it may absorb water from the bulk by infiltration or inject into the bulk through recharge. Given a set G ∈ R 2 and a point x ∈ ∂G, we denote by t G (x) the unique normalized tangential vector and by n G (x) its outward boundary normal (see Fig. 1 for G = Ω). We take friction into account by considering the following Navier boundary condition on the bottom B := {(t, x, Z(x)) : t > 0 and 0 < x < L}: whilst the recharge-infiltration mechanism is modelled with the following permeable boundary condition: The scalar k(u), which models a general kinematic friction law, is defined by k(ξ) := (C lam + C tur |ξ|), for all ξ ∈ R 2 and some C lam , C tur 0. (2.2. 3) The friction coefficients C lam and C tur correspond, respectively, to the laminar and turbulent friction factors [Wylie and Streeter, 1978, Streeter et al., 1998, Gerbeau and Perthame, 2001, Levermore and Sammartino, 2001, Marche, 2007. Furthermore, we consider the following infiltration friction law which models the friction caused by the water's recharging through the ground with average microscopic velocity-rate zero in the horizontal direction; for simplicity we assume a piecewise linear function of I with coefficient α accounting for the magnitude of the frictional effect. The infiltration function I models the amount of water that leaves (I > 0) or enters (I < 0) the flow per elementary boundary element. Notice that, although I should in principle be thought as a function of h, u, and possibly their derivatives, particularly σ[u], as in the recognised Beavers-Joseph-Saffman model described, for instance, by Beavers and Joseph [1967], Jäger and Mikelić [2000], Saffman [1971], Badea et al. [2010], we ignore this in this paper and consider the function I as a given function of space-time. We also note that the recharge-induced frictionf I is active only when water is entering the flow (I < 0); it is thus zero when water infiltrates the ground (I > 0). We define Ω's tangential and outward unit normal vectors on B by respectively, following the convention that the outward normal is the tangential vector rotated by π/2 counterclockwise. It thus follows that (2.2.1) and (2.2.2) on B can be rewritten, respectively, as 2.3. The free boundary. On the free boundary, i.e z = H(t, x), we neglect all other meteorological phenomena except for precipitation in the form of additional water through direct rainfall and runoff. Assuming a kinematic boundary condition, we set where R(t, x) is the recharge rate due to rainfall. The unit tangential and normal vectors t Ω and n Ω to the free surface can be explicitly computed in terms of H as 2) which leads to the following explicit form of the kinematic boundary condition: We also assume a stress condition on the free surface, given by where f R = αR models the friction effect of the rain droplets on the free surface, with α again representing the magnitude of this effect. Using the tangential and normal vectors as above, this condition becomes

Shallow water equation with recharge via vertical averaging
We now proceed to write the Navier-Stokes equations with adapted boundary conditions in non-dimensional form. Next, under an assumption on the shallowness of the ratio of the water height to the horizontal domain (represented by small parameter ε), we formally make an asymptotic expansion of the Navier-Stokes system to the hydrostatic approximation at first order. Finally, we derive the Saint-Venant system through an integration on the water height. Our development follows an approach established by Gerbeau and Perthame [2001], also found in Ersoy [2015].
3.1. Dimensionless Navier-Stokes equations. To derive the shallow water model, we assume that the water height is small with respect to the horizontal length of the domain and that vertical variations in velocity are small compared to the horizontal ones. This is achieved by postulating a small parameter ratio where D, L, V , and U are the scales of, respectively, water height, domain length, vertical fluid velocity, and horizontal fluid velocity. As a consequence the time scale T is such that We also choose the pressure scale to be The rationale for the choice (3.1.3) is that we are focusing on the effect of the horizontal forces as mass per horizontal acceleration which has a force scale of and these forces are applied to vertical boundary scale to give the pressure scale It is convenient to define L, U , and thus T , as finite constants with respect to ε → 0, while D = εL and V = εU . This allows us to introduce the dimensionless quantities of timet, space (x,z), pressurep and velocity field (ũ,ṽ) via the following scaling relationst (3.1.6) We also rescale the laminar and turbulent friction factors, and the infiltration and rainfall rates: Note that in the assumed asymptotic setting, C lam,0 and C tur,0 are constants with respect to ε and this implies that C lam and C tur vanish linearly with ε → 0. Finally, we define the following non-dimensional numbers: Froude's number, Fro := U/ g D, Reynolds's number with respect to µ, Rey := ρ 0 U L/µ.
(3.1.9) and consider the following asymptotic setting where µ 0 is the viscosity.
Using these dimensionless variables in the Navier-Stokes equations (2.1.6) and (2.1.7), and reordering the terms with respect to powers of ε, the dimensionless incompressible Navier-Stokes system reads as follows: 1 div u = 0 (3.1.11) Rey (3.1.14) and 3.1.15,ε, u := From the asymptotic setting (3.1.10), and assuming u has bounded second derivatives, definitions (3.1.14) and (3.1.15) formally lead to 3.1.14,ε, u , 3.1.15,ε, u = O(ε). (3.1.16) On the wet boundary B, considering, along with the scaling relations (3.1.6), the dimensionless Navier boundary condition (2.2.6) implies after noting that Rey = O(ε −1 ) and introducing the asymptotic friction laws on the wet boundary. The permeable boundary condition (2.2.7) reads On the free boundary F, the dimensionless free surface boundary condition (2.3.5) becomes with free surface asymptotic friction law Finally, the kinematic boundary condition (2.3.3) is unchanged.
3.2. First order approximation of the dimensionless Navier-Stokes equations. Dropping all terms of O(ε) and above in equations (3.1.11)-(3.1.21), we deduce the hydrostatic approximation of the dimensionless Navier-Stokes system with the following boundary conditions: (3.2.6) Assuming that the pressure exerted by the rain on the free surface p ε (t, x, H) = p c for some constant p c ∈ R, this becomes (3.2.7) Moreover, identifying terms at order 1 ε in (3.2.2), (3.2.4) and (3.2.5), we obtain the motion by slices decomposition (3.2.10) Noting u ε (t, x) as the mean speed of the fluid over the section [Z(x), H(t, x)], (3.2.11) we are able to use the following approximations and drop the first and higher order terms in ε: where q is the discharge defined by In view of the penetration condition (3.2.4) and the kinematic boundary condition (3.2.5), we deduce the following equation: where the source term S := R − I (3.3.4) measures the gain or loss of water through the rainfall and infiltration rates.
Keeping equations (3.2.7), (3.2.8), and (3.2.12) in mind and thanks to the penetration condition (3.2.4) and the kinematic boundary condition (3.2.5), integrating the left hand side of (3.2.2) for z ∈ [Z(x), H(t, x)], we get (3.3.5) Now, integrating the right hand side of (3.2.2) for z ∈ [Z(x), H(t, x)] using the wet boundary condition (3.2.4) and the free surface boundary condition (3.2.5), we obtain: where the friction factors f R ,f I , and k 0 are defined by formulas (3.1.22) and (3.1.19), respectively. Finally, multiplying both sides of each of (3.3.3), (3.3.5), and (3.3.6) by ρ 0 U 2 /D , we obtain the following Saint-Venant system with recharge:  (d) The pair of functions E, E + g h 2 /2 u with E := hu 2 /2 + g h 2 /2 forms a mathematical entropy/entropy-flux pair for system (??), in that they satisfy the following entropy relation for smooth (h, u): (3.4.3) Proof. We consider each statement in turn.
(a) The eigenvalues of the convection matrix D, the Jacobian matrix of the flux, are Applying the product rule to the first term of (3.4.5) and substituting in the conservation of mass equation, we get from which we can cancel Su on both sides. Using the product rule again, we have that which can be substituted into (3.4.6) to give enabling us to now cancel ∂ x hu 2 . We note that (3.4.9) Substituting this into (3.4.8) and dividing by h throughout, we get (3.4.10) Making the further substitution and grouping derivatives of x, we have as required.
(c) Setting u = 0 in equation (3.4.1), we have which is strictly constant, thereby yielding the lake-at-rest steady state.
(d) We begin by noting that E + g h 2 2 u = (ψ − g Z) hu (3.4.14) and hence (3.4.15) Rearranging (3.4.1) we have (3.4.17) Next we substitute Recalling the definition of E in (1.4), we note that E = (ψ − g h /2 − g Z) h, and hence (3.4.21) To conclude the proof we need to show that the terms in the first line of (3.4.21) all cancel. Expanding the first derivative, we have (3.4.22) The last two terms can be rewritten as Using the definition of ψ, we have (3.4.25) The result is thus proved.
3.5. Corollary (energy growth and decay). Let (h, u) satisfy the Saint-Venant equation with rain (1.3) and let E be defined by (1.4). The following statements hold true: (a) Assuming S is negative, i.e., a loss of water from the watercourse, we have the corresponding total energy decay: , a gain of water by the watercourse, we have Ru 2 − min(0, I)u 2 (3.5.2) that is, the sign of the entropy balance (and therefore whether we have growth or decay) is dependent upon the choice of friction effect α.
3.6. Entropy relation for the existing model. We emphasise that Saint-Venant model (1.3) generalises the Saint-Venant model (1.1), in which Sq/h and the additional friction terms in the conservation of momentum equation are neglected. In particular, for model (1.1), properties (a) and (c) in Theorem 3.4 (hyperbolicity of the system and existence of a steady state) still hold, but the entropy relation is altered and the energy-consistency of the system becomes conditional not on the assumed level of friction but on the water height. Using the same notations for E and ψ, one can prove that model (1.1) instead satisfies Theorem 3.4 with the following modifications: (b) For smooth solutions, the mean velocity u = q/h for system (1.1) satisfies: (3.6.1) (d) System (1.1) satisfies the following "entropy relation" for smooth solutions: This indicates that admissible weak solutions only satisfy the entropy inequality on the set h u 2 (S+2k0(u)) /2 g(S−u∂xZ) , and thus the model is only conditionally well-posed. Therefore, as demonstrated in the mathematical derivation, the term Su (or Sq/h) has to be included in the momentum equation.

The numerical model
Our aim is to design a numerical method which can suitably model the shallow water system, extended to include rainfall and infiltration effects. For numerical simulation of the standard shallow water system, a range of methods have been developed, including finite difference methods, which are simpler and easier to implement but at the cost of some accuracy, to the more complex finite volumes, which, whilst being harder to implement, capture the original equations more exactly [Kröner, 1997, LeVeque, 1992, 2002, Toro, 2009. These methods compute the flux between discretised grid cells through approximate Riemann solvers (e.g. Roe, Godunov, HLL), but a drawback of such an approach is that they do not provide all the desirable properties that we would want from our scheme.
4.1. Well balanced schemes. For the standard shallow water equations, a desirable property is the preservation of equilibrium states (called lake at rest), given by h + Z = constant and u = 0.
(4.1.1) Since our Saint-Venant system is no longer a conservation law but a balance law, we have the possibility of water being added to or lost from the lake, and thus this particular equilibrium only holds in the case S = 0, i.e. R = I. We adapt this, therefore, and instead desire that our system preserves the filling the lake state: ∂ t h = R and u = 0; (4.1.2) that is, the rate at which the height of water changes is equal to the rate at which water is added through the rainfall term. Failing to do this would mean a change in mass of water higher or lower than the rate at which it is added, thus violating the balance of mass property of our system. If we wish to maintain these properties, we cannot rely on the usual finite difference or finite volume methods, and thus a well-balanced scheme is required. Such an approach can be found by going back to a kinetic interpretation of the system, as detailed in Perthame and Simeoni [2001], Ersoy [2015]. It is this kinetic reformulation that we will use to derive a numerical method with the properties we wish to have; we introduce a real function χ(ω), defined rigorously below, which we use to turn our Saint-Venant system into a kinetic equation. These kinetic solvers can be modified to preserve the filling-the-lake state, while at the same time maintaining their simplicity and stability properties. One of the direct benefits of using such a approach for the Saint-Venant system is the ability of the kinetic solver to deal with dry soil cases (that is, when h = 0), which will be of importance in ensuring our model continues to function if infiltration causes the water level to fall closer to zero, which might otherwise cause some complications.

Kinetic functions.
We begin with an overview of the kinetic formulation proposed by Perthame and Simeoni [2001] and further developed by Bourdarias et al. [2014], Ersoy [2015]. We consider a kinetic averaging weight function χ : R → R and a kinetic density function M satisfying These functions originate in the kinetic theory where M (t, x, ξ) accounts for the density of particles with speed ξ at the space-time point (t, x). As far as a numerical method goes, the goal is for the derivation of the finite-volume scheme fluxes to be based on M , through the following property which links the macroscopic variables with the microscopic ones.
4.3. Proposition (macroscopic-microscopic relations). Let the functions h, u solve the shallow water system (1.3) and M as in (4.2.2). If h(t, x) > 0 at (t, x) then the following macroscopic-microscopic relations hold For the right-hand side of the conservation of momentum equation in the Saint-Venant system (1.3), we note that it can be rewritten as, [Bourdarias et al., 2011, e.g.], To rewrite in a divergence form, we introduce the nonlinear flux integral operator for each x ∈ (0, L) and system (1.3) becomes The kinetic approach allows to write the system into a single scalar equation with an extra variable; more specifically, we set (0, T ) × (0, L) × R (t, x, ξ) → M (t, x, ξ) as the solution of the following semilinear kinetic partial integro-differential equation where we are using the following moment notation for m = 0, 1, . . . 4.4. Remark (how is the kinetic formulation used). In general, it is easier to find a numerical scheme to solve equation (4.3.5) for M that has the properties we desire, such as entropy stability, than to solve the full shallow water system for h and u. However, in finding M , we can calculate h and hu by virtue of the macro-microscopic relations (proposition 4.3), leading to (h, u) satisfying (and thus solving) (4.3.4). In fact, M needs never be calculated, nor approximated, explicitly, but only the function is used to build the fluxes appearing in a finite volume method as we shall explain in §4.5.

4.5.
Discretisation and kinetic fluxes. The kinetic equation we have for our Saint-Venant system differs from that considered by Perthame and Simeoni [2001] by the additional S M M 0 term. Through solving the kinetic equation for the standard shallow water system, i.e. equation (4.3.5) with S = 0, they developed the following kinetic scheme, itself based on the general method for developing finite volume schemes: The definition of the right and left numerical fluxes for cell c i at discrete time n, F n i±1/2 , will be described below in §4.6. We follow the same process for our Saint-Venant system, which naturally becomes where S n i is a discretisation of the combined rain and infiltration terms. For the choice of time-step, ∆t, we use the following condition: where CFL is the Courant stability constant which lies in (0, 1]. 4.6. Construction of the numerical fluxes. The construction of the numerical fluxes F n i±1/2 is based on the Nemitskii-type operator associated withM given in (4.4.1). While details if our construction may be found in Perthame and Simeoni [2001] or Bourdarias et al. [2014], we here give a quick guideline here for completeness: where the intermediate quantities M ∓ i±1/2 (ξ) are realised as upwinded fluxes: (4.6.2) with M n i±1/2 := M n i (−ξ)1 |ξ| 2 2 g ∆W n i±1/2 + M n i±1 ∓ |ξ| 2 − 2 g ∆W n i±1/2 1 |ξ| 2 2 g ∆W n i±1/2 ; (4.6.3) we use the Iverson notation defined in (2.1.10).
The term ∆W n i±1/2 is the upwinded source term, and provides the jump condition necessary for a particle in one cell to overcome the friction and topography to move to an adjacent cell. Numerically, we calculate this term as: where for each given cell c i = [x i−1/2 , x i+1/2 ] and t > 0 (4.6.5) whereŴ is the operator defined by (4.3.3) and, for φ = h or u, φ n is the cell-wise constant functions with φ n (x) = φ n i for x ∈ c i . The semidiscretised kinetic density, (4.6.6) The discretisation we use in our scheme will be based upon the Barrenblatt kinetic weighting function χ(ω) = 1 π g (2 g −ω 2 ) + for ω ∈ R, (4.6.7) where (X) + stands for the positive part of X, 0 ∨ X [Perthame and Simeoni, 2001, eq.(2.13)].

Numerical tests
The kinetic scheme we use for our numerical method was implemented by extending the code of Besson and Lakkis [2013] to account for the additional source term in (4.5.3), and we present here several simple numerical tests to demonstrate the validity and application of our Saint-Venant system and the associated numerical method. Since the infiltration and precipitation term have almost the same mathematical and numerical difficulties, we will consider only the rain term, and thus take I ≡ 0. Numerical tests with a realistic infiltration term will be considered in a forthcoming paper.
5.1. Discussion with respect to the friction coefficient α. We start by studying the influence of the rain-induced friction effect f R = αR, which is included in Equations (1.3) and omitted in Equations (1.1). As we saw in Corollary 3.5, for S > 0 (i.e. a net gain of water into the system), the entropy relation depends upon the value of α as per (3.5.2). Similarly, the presence of the term f R in the conservation of momentum equation indicates a dependence on α for the change in momentum and velocity. We consider a rainfall-runoff process on a river with zero topography (i.e. Z ≡ 0). We prescribe periodic boundary conditions and assume an initial height and discharge of h(0, x) = q(0, x) = 1 ∀x ∈ [0, 10] (5.1.1) The rainfall intensity is applied uniformly on the river as a function of time: The parameters we use are as follows: final time T = 1, rainfall intensity R 0 = 1, Courant number CFL = 0.95, meshpoints N = 1000. (5.1.3) Using these parameter values and assuming that the spatial derivative in both the mass and momentum equation can be neglected, our extended shallow water system simplifies to: We can solve explicitly for h, q, and u (which we note now only depend on t) as follows: h(t) = t + 1, (5.1.5) For the entropy relation, which simplifies to and comprises both the kinetic energy, K = hu 2 /2, and the potential energy due to water height, g h 2 /2. In our case, we are only interested in the change in kinetic energy, and thus our entropy relation becomes Using these equations, we can plot how the momentum, velocity, and change in kinetic energy depend on the friction level α, giving us a total of seven separate regimes: (i) α < 0: the friction effect acts with the flow; momentum, velocity, and kinetic energy all increase. (ii) α = 0: momentum increases as velocity stays the same; kinetic energy increases at a fixed rate. (iii) 0 < α < 1 /2 : momentum increases as velocity decreases; kinetic energy increases but the rate slows over time. (iv) α = 1 /2: momentum increases as velocity decreases; kinetic energy does not change.
(v) 1 /2 < α < 1: momentum increases as velocity decreases; kinetic energy decreases. (vi) α = 1: momentum is conserved as the friction from rain balances the increase in height; velocity and kinetic energy both decrease. (vii) α > 1: rain friction slows the flow faster than the increase in height; momentum, velocity, and kinetic energy all decrease.
It follows that the only physically reasonable cases are the last two, when α 1. Note that model 5.2. Comparison with real-world data. For our second test, we consider how our numerical scheme compares with data taken from a real-world experiment. The experiment in question concerns a slope with a constant (but non-zero) gradient, an initial height h 0 = 0, and initial discharge q 0 = 0. Rain falls onto the slope uniformly at a constant rate within a given time interval, and we measure the discharge at the downstream edge of the slope.  The topography and rain function are given by: The hydrograph for the experiment, together with the computed values, is provided below. We note that both sets of results compare well. The simulated results depict the three characteristic regions of the test: an initial phase in which the discharge is increasing, a second phase where the discharge stabilises following a peak up to the time at which the rain stops falling, and a third final phase where no rain is falling and the discharge decreases gradually over time.

5.3.
Single-level and three-level cascade. For our final test, we consider a higher intensity rainfall-runoff process on a much shallower slope, and compare 4 m 4 m 4 m Figure 5. Topography for the three-level cascade.
how the water flows when the gradient of the slope is constant across the full domain, and when the gradient decreases from the upstream to downstream end. The experiment is run three times, with the rain falling for T R = 10, 20, and 30 seconds at a constant rate across the full domain, and with rain-induced friction level α = 0, 1, 5 for the first case of a constant slope, and α = 0, 1 for the second case of a decreasing slope. We measure the height of the flow across the entire domain when the rainfall stops, and also measure the height and discharge at the downstream up to the final time.
The parameters of the experiment are given as follows: The simulations show that the level of assumed friction α has a notable effect on the motion of the flow, slowing down the flow and causing the graphs of the height and momentum to be extended over time. This effect becomes, as one might expected, more pronounced for longer rain times. From a flood modelling perspective, an increase in rain-induced friction while result in longer lasting floods, as well as greater devastation due to the increase in the height profile. Comparing the results from the single-and three-level cascade, we note that the overall profile is reasonably similar, but the three-level cascade induces multiple waves to be formed over time. These waves become more pronounced as α increases, though the length of rain time doe snot appear to have any significant impact. It can also be noted that the cascade profile does not seem to have a major effect on the time taken for the height profile to decrease to zero.

Conclusions
We propose an open channel Saint-Venant system of equations (1.3)describing the flow of a spatially unidimensional watercourse with recharge, e.g., from rain runoff, and infiltration. We derive the model from the two-dimensional Navier-Stokes equations coupled with appropriate boundary conditions modelling recharge via rain, runoff and tributaries on one side and groudwater infiltration or recharge. The proposed model has an additional momentum source and friction terms in comparison with earlier models, such as (1.1) proposed by ? which become special cases of our model. The friction terms are obtained naturally by the derivation process and their presence is essential to explain how water entering the flow picks up the velocity of the flow itself and intertially oppose the flow. The existence of these additional terms leads to a model whose energetic consistency depends exclusively on the level of assumed rain-induced friction coefficient, denoted by α, and whenever this term is dropped, the model is conditionally consistent with respect to the water height. For certain regimes, the conditionally consistent model may yield non-physical solutions. We have illustrated the effect these additional terms have on the results in several numerical tests, as well as validating the model against existing experimental data.