A Well-balanced Finite Volume Scheme for Shallow Water Equations with Porosity Application to Modelling Flow through Rigid Vegetation

Strong interactions exist between flow dynamics and vegetation in open-channel. Depth-averaged shallow water equations can be used for such a study. However, explicit representation of vegetation can lead to very high resolution of the mesh since the vegetation is often modelled as vertical cylinders. Our work aims to study the ability of a single porosity-based shallow water model for these applications. More attention on flux and source terms discretizations are required in order to archive the well-balancing and shock capturing properties. We present a new Godunov-type finite volume scheme based on a simple-wave approximation and compare it with some other methods in the literature. A first application with experimental data was performed.


Introduction
Vegetation is known to play important role in dynamic of open-channel flow.Depending on the shape, density and spatial distribution of vegetation, water depth and flow direction might be significantly modified because the vegetation roughness is much larger than the roughness of river bed [1].Understanding the influence of vegetation on river flow, or more general in environmental fluid mechanics, have been primary interest for decades.Shallow Water (SW) model can provide an accurate representation of physical processes of flow through vegetation.Nevertheless, such an explicit modelling is not suitable in the field, because it leads to very expensive computational cost; furthermore, the real geometry is generally not available or not accurate enough.
It seems more appropriate to use implicit or macroscopic modelling for practical application.Traditional approaches consist in adding a drag force globally or locally into the momentum equation of SW model to enhance the determination of the local velocities.A more advanced macroscopic model, that we are interested here, is to introduce a porosity term into SW model.The porosity, φ, represents the fraction of the plan view area available to flow.For emergent and rigid vegetation, one can consider an isotropic and depth-independent porosity, see Fig. 1.This approach is called single porosity (SP) model [2] whose the mass and momentum conservation equations write where h represents the water depth and u = (u, v) denotes the depth-averaged horizontal velocity; g is the gravity acceleration, b is the bed elevation.The terms τ b and τ d stand for the friction stress and the depth-integrated drag due to vegetation respectively, writing in which n is Manning's coefficient and C d is drag coefficient.The parameter a = 1−φ πD/4 is often termed as frontal area of vegetation of effective diameter D. Since drag force acts upon the fluid occupying only a fraction φ of the total volume, the total drag is thus divided by φ.
Numerical scheme for SP model has been less studied than SW model.In framework of finite volume method, a first scheme was proposed by Guinot and Soares-Frazão [3] in which the authors modified the numerical HLLC flux to account the porosity source term.This method is found to be efficient for shock capturing but it is inaccurate for steady solutions of the model.Next, several approaches of Roe-type have been proposed, one can cited [4,5].These approaches are known to have the same difficulties in preserving the positivity of water depth or when dealing with critical state (sonic point).We would like to mention here a third approach proposed by Finaud-Guyot et al. [6] which particularly holds our attention.Under assumption that all waves are rarefactions, the solution can entirely be determined using Riemann invariants.This solver, namely PorAS, is shown to be very accurate for regular solution, including the steady ones, but has difficulties for estimating shock waves.
We aim to study a robust scheme which inherits, on one hand, the good properties of HLLC solver, such as positivity preserving, shock capturing and easy to implement; on the other hand, the scheme captures accurately steady solutions as with PorAS method.Therefore, we have considered a suitable simple-wave approximation of solution on which exact Riemann invariants are imposed.
The paper is organized as follows: we first recall the two-dimension finite volume formalism whose numerical fluxes at each cell's interface are obtained by solving a projected one-dimension Riemann problem.Next, we detail and analyse the construction of the simplesolver.An analytical test case for illustrating the attractive behaviours of the method is presented.Finally a real application with experimental data was performed.

Numerical scheme
We propose and analyse here a novel finite volume discretization for system (1).It should be convenient to rewrite it under a vectorial form of conservation law with source terms such as in which we have denoted the conservative variable W = (φh, φhu, φhv) T , the fluxes F(W), G(W), the non-conservative source term S(W, b) due to bathymetry and porosity gradients, also the source term S τ (W) accounting the friction and additional drag.

Two-dimension finite volume formalism
Let Ω denote the computational domain discretized by a simplex mesh T h .Given an approximation {φ, b} j of the geometric data on the mesh and assuming that a piecewise constant approximation W n j at time t n = n∆t is known, finite volume scheme consists in computing the updated solution W n+1 j at next time level t n+1 = t n + ∆t.Well-balanced scheme, i.e. that preserves at least the steady state at rest (u = 0, h + b = const.),becomes nowadays a prerequisite criteria for modern numerical method.A classical way to design such a discretization is to solve system (3) in two following steps.
Convection.In order to balance the convective terms and the geometrical source terms, i.e. which contain gradient of φ and b, we solve first the system . The resulting scheme writes ) is an approximation of flux and source terms along the common edge Γ jk of two adjacent cells C j and C k .Thanks to rotational invariance property, that is with R n i j being the rotation matrix and n jk = (n 1 , n 2 ) standing for the outward unit normal vector to Γ jk , the numerical flux can thus be derived from the one-dimension system.This later will be detailed in the next section.
Friction and drag.The next step is to account the friction and drag momentum source terms by solving ∂ t W = −S τ with W(x, 0) = W n+1/2 .Regarding empirical law (2), we used a semi-implicit discretization which writes (5)

A one-dimension Godunov-type scheme
Let W L,R = R n jk W n j,k be the corresponding left-and right-states, we are concerned now to solve the self-similar solution W(x/t) of one-dimension Riemann problem in which S x (W, b) = (0, g 2 h 2 ∂ x φ − gφh∂ x b, 0) T is the geometrical source term in direction n jk .As reported in [6], system (6) has three characteristic fields propagating with speeds It is shown that the porosity and bathymetry remain constant along all these wave curves.They may change only across a stationary wave, and along which the following Bernoulli's relation, also called well-balancing property, has to be satisfied Approximation by simple-solver.We consider therefore a simple-solver W R (x/t) composed by the given data W L , W R and three intermediate states They are separated by four discontinuities waves propagating with velocities λ L ≤ λ 0 = 0 ≤ λ R and λ * as illustrated in Fig. 2. The first order Godunov-type scheme based on this simple-solver can be written as A four-wave approximate solution of Riemann problem (6).
where ∆x stands for the one-dimension mesh size, and for each cell's interface, the left-and right-numerical fluxes with λ * ± standing for the positive and negative parts of λ * .One-and two-dimension fluxes can be linked to each other by relation Determination of intermediate states.The simple-solver has to preserve the averaged value of exact solution at interface, saying under a half-CFL condition, i.e. ∆t max{−λ L , λ R } ≤ ∆x/2.Additional relations on intermediate states can be imposed in order to be consistent with well-balancing property and Riemann invariants of exact solution Now, let W HLL denote the usual HLL-state of the homogeneous Riemann problem, i.e.
for which λ L and λ R are estimations of slowest and fastest wave speeds, intermediate states of the simple-solver can be seen a priori as some perturbations of W HLL due to the source term S x (W, b) := (0, S x , 0) T .To see this, we integrate first the conservation law (6) on the rectangular C = [−∆x/2, ∆x/2] × [0, ∆t], see again Fig. 2. We use next the consistency (9) and equation ( 10) to obtain the following relations, after some manipulations, (φh Therefore, by equation ( 14) resulting from mass consistency, (φh) HLL is nothing that a convex combination of φ L h * L and φ R h * R .Equation (15) expresses momentum consistency and allows to compute intermediate discharge q * from that of HLL-state once given an approximation S x of the source term.We discuss later how such an approximation can be made.
Once q * is known, we turn now to solve intermediate water depths by employing the well-balancing condition (11).This equation can be rewritten under the form Combining it with (14), by which φ R h * R can be seen as function of φ L h * L , the well-balancing condition results thus a nonlinear equation in φ L h * L .Let consider further a natural condition saying that intermediate states have the same regime, i.e. they are both sub-critical or supercritical, this nonlinear equation admits thus an unique and positive solution which can be solved numerically by any iterative method.Solving this equation allows to provide very accurate result, in particular for the case with large porosity gradient and/or with steep bottom.
An alternative approach, which is less accurate but faster and preserves as well steady state at rest, is to replace (11) by a hydrostatic approximation, being

Coupling again with (14) results an explicit expression of intermediate water depths, writing
It could be checked that the fluxes F L,R are non-conservative for second component due to the source term, that is F L φhu F R φhu , while they are conservative for first and third components, i.e.
In practical, we have not to compute λ * since straightforward manipulations show that F φh and λ * have the same sign, and furthermore, the flux F φhv can be expressed under upwinding form Source term approximation.Deriving an appropriate approximation S x is a key point of the scheme.This consists in defining a numerical value of water depth and porosity at the interface, and can be done by investigating well-balancing property.Indeed, considering now the case where W L and W R are steady states at rest, that is u L = u R = 0 and h L + b L = h R + b R , these states are preserved by the scheme if q * = 0. Substituting this into (15) leads to The first test case is one-dimension and aims to assess the shock capturing ability.Next, a real application of the scheme for macroscopic modelling of flow with vegetation is found in the second test case, for which numerical results and experimental data are compared.

Dambreak with large porosity discontinuity
We consider a dambreak on flat bottom but with large porosity discontinuity.This test case was proposed in [3] and consists in solving Riemann problem (6) with initial condition We notice that discontinuity on porosity at the dam position leads to a complex structure of solution compared to the case of SW model, i.e. with constant porosity.Precisely, from left-state W L , the solution starts first with a 1-rarefaction wave following by a stationary discontinuity to reach a critical state W * , that is u * = gh * , just after the dam; it continues again with a 1-rarefaction wave and links finally to right-state W R by a 2-shock, see Fig The proposed scheme with Bernoulli resolution is able to capture correctly the exact solution of Riemann problem while hydrostatic approximation can not, as we can observe on Fig. 3 which is the results at T = 3s computed with ∆x = 10 −2 m.Indeed, this can be understood by the fact that hydrostatic approximation leads to only one intermediate depth because the topography is flat, while solving Bernoulli relation allows to account the discontinuity of porosity.Finally we recall that the construction of PorAS scheme is based on rarefaction waves hypothesis, and, therefore it might present some difficulties when estimating shock waves speed, as noticed in [6] by their authors.

Transition from meadow to wood
We consider here a case of longitudinal transition from meadow to wood in an 18m long and 1m width laboratory flume, and with two types of hydraulic roughness: a bed-roughness figuring a highly submerged dense meadow and emergent macro-roughness figuring a forest.The longitudinal bottom slope was S 0 = 1.05mm/m.Wood-type vegetation was modelled using circular cylinders of diameter D = 10mm, uniformly distributed in staggered rows with density N = 81 cylinders/m 2 , see Fig. 4 (left).We refer to [7] for more details.
Explicit simulation with SW model can provide correct result but leads to very expensive computation cost.Indeed, the circular form of cylinder requires a very high mesh resolution whose the element diameters range typically from 1mm (near the cylinders) to 10mm, see Fig. 4 (right).This restriction on mesh resolution can of course be relaxed when using SP model.In our study, we used an uniform resolution of 10mm, that is the cylinder's diameter.Next, porosity of the elements close to a cylinder is the ratio of cylinder's area, πD 2 /4, and the total surface occupied by these elements.This simple setting results values of porosity ranging from 0.8 to 0.9.Otherwise, porosity of the elements in water region is set to unity.Numerical simulations were made for the case of uniform discharge 0.015m 2 /s.
Strickler roughness coefficient of the meadow was set to K s = 60.24m 1/3 /s based on the measures from several uniform flows on meadow channel without wood transition.Drag coefficient was evaluated as C d = 1.2, see [7] for more details.Our first simulation was performed with this reference set of values.On Fig. 5 we visualize the velocity field and the unit discharge in vegetation region given by SW and SP models.One can observe that SP model is able to cope macroscopic behaviours: the flow is accelerated between two longitudinal rows of cylinders (called fast vein) while it is decelerated after each cylinder.Slight deviation of flow from the main direction around the cylinders is also observed with SP model.nearly constant within vegetation region.In fact, bed friction is negligible compared to drag force.The water depth predicted by SW model compared to the measured one confirms again that the SW model can be used for detailed modelling of interactions between vegetation and the flow.Qualitative behaviours of water depth profile is well captured by SP model.Moreover, a good agreement can be found in vegetation region for the result computed with reference value of friction and drag coefficients (blue curve in Fig. 6).Nevertheless, water depth is overestimated about 2mm in meadow region; this can be improved by imposing a smaller drag coefficient, but in that case, the predictive quality in vegetation part will be reduced.A value C d = 1 seems to be the best compromise for two regions.

Conclusion
We have presented in this paper a novel finite volume scheme of Godunov-type for SP model.The solver is based on a four-wave approximation of Riemann problem.This can be seen as an augmented HLLC scheme which is well-balanced, positivity preserving and shock capturing.Details on practical implementation of the scheme were also discussed.Next, a first application to modelling interactions between rigid vegetation with the flow in a laboratory flume was performed.Good agreement was found both with experimental data and the result provided by an explicit simulation with SW model.

Figure 3 .
Figure 3. Dambreak with very large porosity discontinuity.The hydrostatic approximation converges to wrong solution.

Figure 5 .
Figure 5. Visualization of velocity field and unit discharge (represented by background color) in vegetation region: SW model (top) and SP model (bottom).

Fig. 6
Fig. 6 presents water depths measured along the line y = 0.5m.One can observed that the water depth increases upstream of roughness transition, i.e. in meadow region, and becomes

Figure 6 .
Figure 6.Longitudinal profile of water depth along y = 0.5m.