A Shallow Water Model for Water Flow on Vegetated Hillslope

TL;DR

This paper develops a PDE-based shallow water model incorporating vegetation effects to simulate water flow on vegetated hillslopes, validated through numerical tests against analytical and real-world data.

Contribution

It introduces a novel PDE model that accounts for vegetation porosity in water flow simulation on hillslopes, using finite volume discretization.

Findings

Model accurately predicts water flow dynamics in vegetated terrains.

Numerical results align well with laboratory and real basin data.

The model captures key physical properties of water-vegetation interactions.

Abstract

The hillslope hydrological processes are very important in watershed hydrology research. In this paper we focus on the water flow over the soil surface with vegetation in a hydrographic basin. We introduce a PDE model based on general principles of fluid mechanics where the unknowns are water depth and water velocity. The influence of the plant cover to the water dynamics is given by porosity (a quantity related to the density of the vegetation), which is a function defined over the hydrological basin. Using finite volume method for approximating the spatial derivatives, we build an ODE system which constitutes the base of the discrete model we will work with. We discuss and investigate several physical relevant properties of this model. Finally, we perform different quantitative validation tests by comparing numerical results with exact solutions or with laboratory measured data. We also consider some qualitative validation tests by numerically simulating the flow on a theoretical vegetated soil and on a real hydrographic basin.

Tables

Table 1. Table 1: Data for Riemann problem

	Non-existence	Uniqueness	Two solutions
${(h,u)}_L$	$(0.2,5.0)$
${(h,u)}_R$	$(0.6,0.4)$	$(0.6,1.34)$	$(0.6,0.4)$
${(\theta ,z)}_L$	$(0.3,1.0)$	$(0.8,1.0)$
${(\theta ,z)}_R$	$(1.0,1.2)$

Equations

$$x^3=z(x^1,x^2),(x^1,x^2)\in \Omega$$

$$\begin{array}{rl}{\partial }_t(\theta h)+{\partial }_b(\theta h{v}^b)& =\mathfrak{M},\\ {\partial }_t(\theta h{v}_a)+{\partial }_b(\theta h{v}_a{v}^b)+\theta h{\partial }_{a}w& ={\mathfrak{t}}_a^p+{\mathfrak{t}}_a^s,a=1,2,\end{array}$$

$$w=g\left[z(x^1,x^2)+h\right]$$

$${\mathfrak{t}}_a^p={\alpha }_{p}h\left(1-\theta \right)|\mathbf{v}|v_a,$$

$${\mathfrak{t}}_a^s=\theta {\alpha }_s(h)|\mathbf{v}|v_a,$$

$${\mathfrak{t}}_a^p+{\mathfrak{t}}_a^s:=-\mathcal{K}(h,\theta )|\mathbf{v}|v_a,$$

$$\mathcal{K}(h,\theta )={\alpha }_{p}h\left(1-\theta \right)+\theta {\alpha }_s(h).$$

$$\mathcal{E}:=\frac{1}{2}|\mathbf{v}{|}^2+g\left(x^3+\frac{h}{2}\right),$$

$${\partial }_t(\theta h\mathcal{E})+{\partial }_b\left(\theta h{v}^b\left(\mathcal{E}+g\frac{h}{2}\right)\right)=\mathfrak{M}\left(-\frac{1}{2}|\mathbf{v}{|}^2+w\right)-\mathcal{K}|\mathbf{v}{|}^3.$$

$${\partial }_a(x^3+h)=0,v_a=0,a=1,2,$$

$$\begin{array}{rl}{\partial }_t\underset{{\omega }_i}{\int }\theta h\mathrm{d}x+\underset{\partial {\omega }_i}{\int }\theta h\mathbf{v}\cdot \mathbf{n}\mathrm{d}s=& \underset{{\omega }_i}{\int }\mathfrak{M}\mathrm{d}x,\\ {\partial }_t\underset{{\omega }_i}{\int }\theta h{v}_a\mathrm{d}x+\underset{\partial {\omega }_i}{\int }\theta h{v}_a\mathbf{v}\cdot \mathbf{n}\mathrm{d}s+\underset{{\omega }_i}{\int }\theta h{\partial }_{a}w\mathrm{d}x=& -\underset{{\omega }_i}{\int }\mathcal{K}|\mathbf{v}|v_a\mathrm{d}x,a=1,2.\end{array}$$

$$\underset{{\omega }_i}{\int }\theta h\mathrm{d}x\approx {\sigma }_i{\theta }_i{h}_i,\underset{{\omega }_i}{\int }\theta h{v}_a\mathrm{d}x\approx {\sigma }_i{\theta }_i{h}_i{v}_{ai},\underset{{\omega }_i}{\int }\mathcal{K}|\mathbf{v}|v_a\mathrm{d}x\approx {\sigma }_i{\mathcal{K}}_i|\mathbf{v}{|}_i{v}_{ai},$$

$$\underset{{\omega }_i}{\int }\theta h{\partial }_{a}w\mathrm{d}x=\underset{\partial {\omega }_i}{\int }(w-w_i)\theta h{n}_a\mathrm{d}s+\mathcal{O}(|{\omega }_i|h^2)\approx \underset{\partial {\omega }_i}{\int }(w-w_i)\theta h{n}_a\mathrm{d}s,$$

$$\tilde{\psi }{|}_{\partial \omega (i,j)}:=\underset{\partial \omega (i,j)}{\int }\psi \mathrm{d}s.$$

$$\begin{array}{rl}{\sigma }_i{\partial }_t({\theta }_i{h}_i)+\sum _{j\in \mathcal{N}(i)}\widetilde{\theta h{v}_n}{|}_{\partial \omega (i,j)}& ={\sigma }_i{\mathfrak{M}}_i,\\ {\sigma }_i{\partial }_t({\theta }_i{h}_i{v}_{ai})+\sum _{j\in \mathcal{N}(i)}\widetilde{\theta h{v}_a{v}_n}{|}_{\partial \omega (i,j)}+\sum _{j\in \mathcal{N}(i)}\widetilde{(w-w_i)\theta h}{n}_a{|}_{\partial \omega (i,j)}& =-{\sigma }_i{\mathcal{K}}_i|\mathbf{v}{|}_i{v}_{ai},\end{array}$$

$$\begin{array}[]{l}\widetilde{\theta hv_{n}}|_{\partial\omega(i,j)}\approx l_{(i,j)}\theta h_{(i,j)}(v_{n})_{(i,j)},\\ \widetilde{\theta hv_{a}v_{n}}|_{\partial\omega(i,j)}\approx l_{(i,j)}\theta h_{(i,j)}(v_{a})_{(i,j)}(v_{n})_{(i,j)},\\ \widetilde{(w-w_{i})\theta h}n_{a}|_{\partial\omega(i,j)}\approx l_{(i,j)}(w_{(i,j)}-w_{i})\theta h^{s}_{(i,j)}(n_{a})_{(i,j)},\end{array}$$

$$\begin{array}[]{l}(v_{a})_{(i,j)}=\displaystyle\frac{v_{a\,i}+v_{a\,j}}{2},\quad a=1,2,\\ (v_{n})_{(i,j)}=\boldsymbol{v}_{(i,j)}\cdot\boldsymbol{n}_{(i,j)},\\ w_{(i,j)}=\displaystyle\frac{w_{i}+w_{j}}{2},\end{array}$$

$$\theta h_{(i,j)}=\left\{\begin{array}[]{ll}\theta_{i}h_{i},&{\rm if}\;(v_{n})_{(i,j)}>0,\\ \theta_{j}h_{j},&{\rm if}\;(v_{n})_{(i,j)}<0.\end{array}\right.$$

$$(\theta h)^{s}_{(i,j)}=\left\{\begin{array}[]{ll}\theta h_{(i,j)},&{\rm if}\;(v_{n})_{(i,j)}\neq 0\\ \theta_{i}h_{i},&{\rm if}\;(v_{n})_{(i,j)}=0\;{\rm and}\;w_{i}>w_{j}\\ \theta_{j}h_{j},&{\rm if}\;(v_{n})_{(i,j)}=0\;{\rm and}\;w_{i}\leq w_{j}\end{array},\right.$$

$$(\theta h)^{s}_{(i,j)}=\left\{\begin{array}[]{ll}\theta_{i}h_{i},&{\rm if}\;w_{i}>w_{j}\\ \theta_{j}h_{j},&{\rm if}\;w_{i}\leq w_{j}\end{array}\right.$$

$$(\theta h{)}_{(i,j)}^s=\frac{{\theta }_i{h}_i+{\theta }_j{h}_j}{2}.$$

$$\begin{array}{rl}{\sigma }_i\frac{\mathrm{d}}{\mathrm{d}t}({\theta }_i{h}_i)+\sum _{j\in \mathcal{N}(i)}{l}_{(i,j)}\theta h_{(i,j)}(v_n{)}_{(i,j)}& ={\sigma }_i{\mathfrak{M}}_i,\\ {\sigma }_i\frac{\mathrm{d}}{\mathrm{d}t}({\theta }_i{h}_i{v}_{ai})+\sum _{j\in \mathcal{N}(i)}{l}_{(i,j)}\theta h_{(i,j)}(v_a{)}_{(i,j)}(v_n{)}_{(i,j)}+& \\ +\frac{1}{2}\sum _{j\in \mathcal{N}(i)}{l}_{(i,j)}(w_j-w_i)(\theta h{)}_{(i,j)}^s{n}_a{|}_{(i,j)}& =-{\sigma }_i{\mathcal{K}}_i|\mathbf{v}{|}_i{v}_{ai}.\end{array}$$

$$\begin{array}[]{l}\boldsymbol{v}_{\partial\omega_{i}\cap\Gamma}=\boldsymbol{v}_{i},\\ h_{\partial\omega_{i}\cap\Gamma}=\left\{\begin{array}[]{ll}h_{i},&{\rm if}\;\boldsymbol{v}_{i}\cdot\boldsymbol{n}|_{\partial\omega_{i}\cap\Gamma}>0,\\ 0,&{\rm if}\;\boldsymbol{v}_{i}\cdot\boldsymbol{n}|_{\partial\omega_{i}\cap\Gamma}<0.\end{array}\right.\end{array}$$

$${\mathcal{E}}_i(h_i,{\mathbf{v}}_i)={\theta }_i\left(\frac{1}{2}|\mathbf{v}{|}_i^2{h}_i+\frac{1}{2}g{h}_i^2+g{x}_i^3{h}_i\right).$$

$${\sigma }_i\frac{\mathrm{d}{\mathcal{E}}_i}{\mathrm{d}t}={\sigma }_i\left(\left(w_i-\frac{1}{2}|\mathbf{v}{|}_i^2\right)\frac{\mathrm{d}({\theta }_i{h}_i)}{\mathrm{d}t}+⟨{\mathbf{v}}_i,\frac{\mathrm{d}({\theta }_i{h}_i{\mathbf{v}}_i)}{\mathrm{d}t}⟩\right),$$

$${\sigma }_i\frac{\mathrm{d}}{\mathrm{d}t}{\mathcal{E}}_i+\sum _{j\in \mathcal{N}(i)}{l}_{(i,j)}⟨{\mathcal{H}}_{(i,j)},{\mathbf{n}}_{(i,j)}⟩=-{\sigma }_i{\mathcal{K}}_i|\mathbf{v}{|}_i^3,$$

$${\mathcal{H}}_{(i,j)}=\frac{1}{2}\theta h_{(i,j)}\left(w_i{\mathbf{v}}_i+w_j{\mathbf{v}}_j+⟨{\mathbf{v}}_i,{\mathbf{v}}_j⟩{\mathbf{v}}_{(i,j)}\right).$$

$$\mathcal{H}=\theta h\mathbf{v}\left(\frac{1}{2}|\mathbf{v}{|}^2+w\right)$$

$$\begin{array}[]{rcl}\sigma_{i}\displaystyle\frac{\rm d}{{\rm d}t}{\cal E}_{i}&=&-(w_{i}-\displaystyle\frac{1}{2}|\boldsymbol{v}|_{i}^{2})\sum\limits_{j\in{\cal N}(i)}l_{(i,j)}\theta h_{(i,j)}(v_{n})_{(i,j)}-\\ &&-\left<\boldsymbol{v}_{i},\sum\limits_{j\in{\cal N}(i)}l_{(i,j)}\theta h_{(i,j)}\boldsymbol{v}_{(i,j)}(v_{n})_{(i,j)}\right>-\\ &&-\displaystyle\frac{1}{2}\left<\boldsymbol{v}_{i},\sum\limits_{j\in{\cal N}(i)}l_{(i,j)}(w_{j}-w_{i})(\theta h)^{s}_{(i,j)}\boldsymbol{n}_{(i,j)}\right>-\\ &&-\sigma_{i}{\cal K}_{i}|\boldsymbol{v}|_{i}^{3}.\end{array}$$

$$\begin{array}[]{rcl}w_{i}\displaystyle\sum\limits_{j\in{\cal N}(i)}l_{(i,j)}\theta h_{(i,j)}(v_{n})_{(i,j)}&=&\displaystyle\sum\limits_{j\in{\cal N}(i)}l_{(i,j)}\theta h_{(i,j)}(v_{n})_{(i,j)}\displaystyle\frac{w_{i}+w_{j}}{2}+\\ &&+\displaystyle\sum\limits_{j\in{\cal N}(i)}l_{(i,j)}\theta h_{(i,j)}(v_{n})_{(i,j)}\displaystyle\frac{w_{i}-w_{j}}{2}\end{array}$$

Full text

Keywords: hydrological process, non-homogeneous hyperbolic system, shallow water equations, finite volume method, well-balanced scheme, porosity

2010 MSC: 35L03, 35L60, 35Q35, 65M08, 74F10.

1 Introduction

The mathematical modeling of hydrodynamic processes in hydrographic basins is a challenging subject. There are two main difficulties: the complexity of the phenomena and the multitude of the factors involved in the process. However, there are plenty of performant models dedicated to some specific aspects of the hydrodynamic processes only. To review the existent mathematical models is beyond the purposes of this paper, but we can group them into two large classes: physical base models and regression models. The most known regression models are the unit hydrograph [12] and universal soil loss [43, 50]. From the first class, we mention here a few well known models: SWAT [47], SWAP [46] and KINEROS [51]. Due to the complexity and heterogeneity of the processes (see [38]), models in this class are not purely physical because they need additional empirical relations. The main difference between models here is given by the nature of the empirical relations. For example, in order to model the surface of the water flow, SWAP and KINEROS use a mass balance equation and a closure relation, while SWAT combines the mass balance equation with the momentum balance equation. A very special class of models are cellular automata which combine microscale physical laws with empirical closure relations in a specific way to build up a macroscale model, e.g. CAESAR [8, 24]. In this paper, we introduce a model that extends the shallow water equations and takes into account the presence of vegetation on the soil surface. Such kind of extension is not new; we remind here some applications that use similar models: flow through rigid vegetation [32], flash food propagation in the urban area [18, 19] and tsunami advance in the coastal region [48]. Our model results as an asymptotic limit of a space averaged version of the mass conservation equation and linear momentum balance equations [27, 28]. It takes into account the plant vegetation density, water-plant interaction, soil topography and water-soil interaction. One note that if the cover plant is uniform distributed on the soil, there are not variations in the plant type and plant density, the model reduces to the classical shallow water equation. For numerical applications, we introduce a discrete version of the model. First, we use a finite volume method to discretize the space variable and space derivatives. The resulting ordinary differential system of equations sets up the base of the full discrete model. Then, we introduce a fractional time step scheme to discretize the time variable and time derivatives appearing in the ODE model. Since water flow is the main driving force in more complex phenomena like soil erosion or pollutant spreading, one needs a water flow modulus to be easily incorporated in a large computational package. For this reason, we develop a simple discrete model with low algebraic calculations and reduced memory requirements. For a coupled model water flow

• soil erosion and for the influence of the emergent vegetation on soil erosion, the reader is referred to [22] and [2], respectively.

Acknowledging that there must be a balance between the accuracy of the model and the computational cost of it, we underline that we will focus on a first order numerical method dedicated especially to practical applications on the field. This option is a compromise between the precision of the measured data on the field and the accuracy of the method. High order schemes require higher computing resources and can become prohibitive (due to the enormous increase of the volume of processed data) when they are applied at the natural scale of hydrographic basins.

In Section 2 we introduce the PDE model; the reader is referred to [27, 28] for more details concerning this model. Section 3 is dedicated to the numerical approximation of the solution for this model. Section 3.1 deals with the FVM space discretization of the continuum PDE model. The method is described in some details to facilitate a less familiar reader with the FVM method to understand the discrete version of the continuum model and to allow him to develop a proper approximation scheme. In Section 3.2 we investigate some physical relevant qualitative properties of this ODE system: monotonicity of the energy, positivity of the water depth function $h$ and well-balanced properties of the scheme. In Section 3.3 we obtain the full discrete version of our continuous model. In fact, this last version provides the numerical result of the model.

Section 4 is devoted to the validation of the numerical method and the continuum PDE model. Section 4.1 deals with the internal validation. To analyze the performance of the numerical scheme, we compare the numerical results with the exact solutions of PDE in three cases: flow over a bump, the Tacker solution and the solution of the Riemann problem of the shallow water equations. The solutions of Riemann problem of the shallow water equation with topography and vegetation require a deep discussion. If in the case of the shallow water equation without topography and vegetation there is no doubt concerning the solution of the problem, in the case of shallow water equations with topography and vegetation the definition of solution of the problem require new concepts, measure solution and path connection, [10, 25]. When we compare a numerical solution with a theoretical one, there are two relevant aspects for us: the dependence of the theoretical solution on the definition of the path connection and the non-uniqueness of the theoretical solution.

Four more tests are considered in Section 4.2 for external validation. The first two tests are performed in order to compare the numerical results with measured data in laboratory. The third one is a qualitative test and emphasizes the influence of the plant cover on the asymptotic behavior of the dynamical system defined by the numerical scheme. The last test is a numerical simulation that mimics a real phenomenon, the dumping effect of the vegetation on the downhill water flow.

Last section is dedicated to final remarks and conclusions.

2 Mathematical Model

The model we discuss here is a simplified version a more general model of water flow on a hillslope introduced in [28]. Let $\Omega\in\mathbb{R}^{2}$ be a connected bounded open set. Assume that the soil surface is represented by

[TABLE]

and the first derivatives of the function $z(\cdot,\cdot)$ are small quantities. The unknown variables of the model are the water depth $h(t,x)$ and the two components $v_{a}(t,x)$ of the water velocity $\boldsymbol{v}$. The density of the plant cover is quantified by the porosity function $\theta:\Omega\rightarrow[0,1]$ which equals to $1$ for bare soil and [math] for a complete sealant cover plant. The model reads as

[TABLE]

where

[TABLE]

stands for the free water surface level, and $g$ for the gravitational acceleration. The contribution of rain and infiltration to the water mass balance is taken into account by $\mathfrak{M}$.

The term $\mathfrak{t}^{p}_{a}$ quantifies the resistance opposed by plants to the water flow. We use an experimental relation as in [3, 41]

[TABLE]

where $\alpha_{p}$ is a non-negative constant depending on the type of vegetation.

The term $\mathfrak{t}^{s}_{a}$ quantifies the frictional force exerted by the interaction water-soil. In this case, one can use experimental relations as Manning, Chézy or Darcy-Weisbach formula. The Darcy-Weisbach expression has the advantage to be non-singular if the water depth becomes zero. We analyze the model for

[TABLE]

where $\alpha_{s}(h)$ is a non-negative function characteristic to a given soil surface.

Combining the two relations, one can write

[TABLE]

with the function ${\cal K}(h,\theta)$ given by

[TABLE]

Note that there is an energy function ${\cal E}$ given by

[TABLE]

satisfying the conservative equation

[TABLE]

In the absence of mass source, a steady state of the system is given by

[TABLE]

whose physical meaning is a lake.

When dealing with constant porosity function $\theta<1$, the system (1) is similar to the classical shallow water equations, while for $\theta=1$, we recover the classical shallow water equations.

The model (1) is a hyperbolic system of equations with source term, see [27]. Among different features, we require from our approximation scheme to have two properties with physical significance: $h$-positivity and well-balanced (the numerical solution preserves the lake). We speculate that the lake is an attractive steady state of the numerical scheme if one has a conservative equation for the energy as in the continuous case. Schemes with similar properties but in the absence of vegetation were investigated in [7, 42, 16].

3 Numerical Approximation

3.1 Finite Volume Method Approximation of 2D Model

Let $\Omega$ be the domain of the space variables $x^{1}$, $x^{2}$ and $\Omega=\cup_{i}\omega_{i},i=\overline{1,N}$ an admissible polygonal partition, [36]. To build a spatial discrete approximation of the model (1), one integrates the continuous equations on each finite volume $\omega_{i}$ and then defines an approximation of the integrals.

Let $\omega_{i}$ be an arbitrary element of the partition. Relatively to it, the integral form of (1) reads as

[TABLE]

Now, we build a discrete version of the integral form by introducing some quadrature formulas. With $\psi_{i}$ standing for some approximation of $\psi$ on $\omega_{i}$, we introduce the approximations

[TABLE]

where $\sigma_{i}$ denotes the area of the polygon $\omega_{i}$.

The integrals of the gradient of the free surface are approximated by

[TABLE]

where $w_{i}$ is an approximation of $w$ on $\omega_{i}$.

Note that if $\omega_{i}$ is a regular polygon and $w_{i}$ is the cell-centered value of $w$, then the approximation is of second order accuracy for smooth fields and it preserves the null value in the case of constant fields $w$.

Let $\partial\omega(i,j)$ be the common interface between the cells $\omega_{i}$ and $\omega_{j}$. We introduce the notation

[TABLE]

Using the approximations (10) and (11) and keeping the boundary integrals, one can write

[TABLE]

where ${\cal N}(i)$ denotes the set of indexes of all the neighbors of $\omega_{i}$.

The next step is to define the approximations of the boundary integrals in (13). We approximate an integral $\widetilde{\psi}|_{\partial\omega(i,j)}$ of the form (12) by considering the integrand $\psi$ to be a constant function $\psi_{(i,j)}(\psi_{i},\psi_{j})$, where $\psi_{i}$ and $\psi_{j}$ are some fixed values of $\psi$ on the adjacent cells $\omega_{i}$ and $\omega_{j}$, respectively. Thus,

[TABLE]

where $\boldsymbol{n}_{(i,j)}$ denotes the unitary normal to the common side of $\omega_{i}$ and $\omega_{j}$ pointing towards $\omega_{j}$, and $l_{(i,j)}$ is the length of this common side.

The issue is to define the interface value functions $\psi_{(i,j)}(\psi_{i},\psi_{j})$ such that the resulting scheme has certain desired properties. Among them, we impose the scheme to be well-balanced and to preserve the positivity of the water depth function $h$.

Well-balanced and h-positive scheme. For any internal interface $(i,j)$, we define the following quantities:

[TABLE]

and

[TABLE]

For the term $(\theta h)^{s}_{(i,j)}$, we analyze the alternatives

[TABLE]

[TABLE]

or

[TABLE]

The semidiscrete scheme takes now the form of the following system of ODEs

[TABLE]

Boundary conditions. Free discharge. We need to define the values of $h$ and $\boldsymbol{v}$ on the external sides of $\Omega$. For each side in $\Gamma=\partial\Omega$, we introduce a new cell (“ghost” element) adjacent to the polygon $\omega_{i}$ corresponding to that side. For each “ghost” element, one must somehow define its altitude and then we set zero values to its water depth. We can now define $h$ and $\boldsymbol{v}$ on the external sides of $\Omega$ by

[TABLE]

Now, the solution is sought inside the positive cone $h_{i}>0,\;i=\overline{1,N}$.

3.2 Properties of the Semidiscrete Scheme

The ODE model (20) is the basis of the numerical scheme for obtaining a numerical solution of the PDE model. Therefore, it is worthwhile to analyze the properties of (20). Unfortunately, not all the properties of the solution of the ODE model can be preserved by its numerical approximation. Numerical schemes preserving properties of some particular solutions of the continuum model were and are intensively investigated in the literature [7, 16, 5, 39]. We analyze such properties for the semi-discrete scheme in the present section and dedicate the next section to the properties of the complete discretized scheme. Depending on the scheme used to approximate the terms $(\theta h)^{s}_{(i,j)}$, the ODE model (20) has different properties. If one uses the formula (17), then the ODE model can have discontinuities in r.h.s. and therefore it is possible that the solution in the classical sense of this system might not exist for some initial data. However, the solution in Filipov sense [20] exists for any initial data. There are initial data for which the solution in the classical sense exists only locally in time. As a reward, there is an energy balance equation and the energy is a monotone function with respect to time. On the other hand, if the alternative formula (18) or (19) is used, then the r.h.s. of the ODE model is defined by continuum functions, but we can not anymore find an energetic function.

3.2.1 Energy Balance

Definitions (15) and (17) yields a dissipative equation in conservative form for the cell energy ${\cal E}_{i}$,

[TABLE]

The time derivative of ${\cal E}_{i}$ can be written as

[TABLE]

where $\left<\cdot,\cdot\right>$ denotes the euclidean scalar product in $\mathbb{R}^{2}$.

Proposition 3.1** (Cell energy equation).**

In the absence of mass source, one has

[TABLE]

where

[TABLE]

Remark 3.1**.**

If $(\theta h,v,w)_{j}=(\theta h,v,w)_{i}$ for any $j\in{\cal N}(i)$, then

[TABLE]

is the continuous energy flux in (7).

Proof.

Using the equality (23), we can write

[TABLE]

Now, one has the identities

[TABLE]

and

[TABLE]

Therefore,

[TABLE]

Similarly, one obtains the identity

[TABLE]

∎

Taking out the mass exchange through the boundary, the definitions of the interface values ensure the monotonicity of the energy with respect to time.

3.2.2 h-positivity and Critical Points

Proposition 3.2** (h-positivity).**

If $\mathfrak{M}\geq 0$ when $h=0$, then the ODE system (20) with (15), (17), (16) and (21) preserves the positivity of the water depth function $h$.

Proof.

One can rewrite the mass balance equations as

[TABLE]

Observe that if $h_{i}=0$ for some $i$, then $\sigma_{i}\displaystyle\frac{\rm d}{{\rm d}t}(\theta_{i}h_{i})\geq 0$. ∎

There are two kinds of stationary points for the ODE model: the lake and uniform flow on an infinitely extended plan with constant vegetation density.

Proposition 3.3** (Stationary point. Uniform flow.).**

Consider $\{\omega_{i}\}_{i=\overline{1,N}}$ to be a regular partition of $\Omega$ with $\omega_{i}$ regular polygons. Let $z-z_{0}=\xi_{b}x^{b}$ be a representation of the soil plane surface. Assume that the discretization of the soil surface is given by

[TABLE]

where $\overline{x}^{b}_{i}$ is the mass center of the $\omega_{i}$ and $\theta_{i}=\theta$. Then, given a value $h$, there is $\boldsymbol{v}$ so that the state $(h_{i},\boldsymbol{v}_{i})=(h,\boldsymbol{v})$, $i=\overline{1,N}$ is a stationary point of the ODE (20).

Proof.

For any constant state $h_{i}=h$ and $(v_{a})_{i}=v_{a}$, the ODE (20) reduces to

[TABLE]

Introducing the representation (25), one writes

[TABLE]

Note that for a regular partition, one has the identity

[TABLE]

where $y_{(i,j)}$ is the midpoint of the common side $\overline{\omega}_{i}\cap\overline{\omega}_{j}$. Taking into account

[TABLE]

we obtain that the velocity is a constant field

[TABLE]

∎

A lake is a stationary point characterized by a constant value of the free surface and a null velocity field over connected regions. A lake for which $h_{i}>0$ for any $i\in\{1,2,\ldots,N\}$ will be named regular stationary point and a lake that occupies only a part of a domain flow will be named singular stationary point.

Proposition 3.4** (Stationary point. Lake.).**

Denote by $\bar{w}$ some real number value. In the absence of mass source, the following properties hold:

(a)* Regular stationary point: the state*

[TABLE]

is a stationary point of ODE (20).

(b)* Singular stationary point: if $(\theta h)^{s}_{(i,j)}$ is calculated through (17) or (18), then the state*

[TABLE]

for some ${\cal I}\subset\{1,2,\ldots,N\}$ is a stationary point. ($\complement{\cal I}$ is the complement of ${\cal I}$.)

Proof.

For the sake of simplicity, in the case of the singular stationary point, we consider that $\displaystyle\Omega_{\cal I}=\cup_{i\in{\cal I}}\omega_{i}$ is a connected domain. Since the velocity field is zero, it only remains to verify that

[TABLE]

for any cell $\omega_{i}$. If $i\in\complement{\cal I}$, then the above sum equals zero since $h^{s}_{(i,j)}=0$, for all $j\in{\cal N}(i)$. If $i\in{\cal I}$, then the sum is again zero because either $h^{s}_{(i,j)}=0$, for $j\in\complement{\cal I}$ or $w_{j}=w_{i}$, for $j\in{\cal I}$. ∎

3.3 Time Integration Scheme

In what follows, we introduce a (first order) fractional time integration scheme in order to integrate the ODE (20).

We introduce some notations

[TABLE]

Now, (20) becomes

[TABLE]

Mass source. We assume that the mass source $\mathfrak{M}$ is of the form

[TABLE]

where $r(t)$ quantifies the rate of the rain and $\iota(t,h)$ quantifies the infiltration rate. The infiltration rate is a continuous function and satisfies the following condition

[TABLE]

The basic idea of a fractional time method is to split the initial ODE into two sub-models, integrate them separately, and then combine the two solutions [35, 45].

We split the ODE (20) into

[TABLE]

and

[TABLE]

A first order fractional step time accuracy reads as

[TABLE]

[TABLE]

Steps (33) and (34) lead to

[TABLE]

To advance a time step, one needs to solve a scalar nonlinear equation for $h$ and a 2D nonlinear system of equations for velocity $\boldsymbol{v}$.

In what follows, we investigate some important physical properties of the numerical solution given by (35): $h$-positivity, well-balanced property and monotonicity of the energy.

3.3.1 h-positivity. Stationary Points

Proposition 3.5** ($h$-positivity).**

There is an upper bound $\tau_{n}$ for the time step $\triangle t_{n}$ such that if $\triangle t_{n}<\tau_{n}$ and $h^{n}>0$, then $h^{n+1}\geq 0$.

Proof.

For any cell $i$ one has

[TABLE]

A choice for the upper bound $\tau_{n}$ is given by

[TABLE]

∎

Proposition 3.6** (Well-balanced).**

The lake and the uniform flow are stationary points of scheme (35).

Proof.

One can prove this result similarly as in propositions 3.3 and 3.4. ∎

Unfortunately, the semi-implicit scheme (35) does not preserve the monotonicity of the energy.

3.3.2 Discrete Energy

The variation of the energy between two consecutive time steps can be written as

[TABLE]

If the sequence $(h,\boldsymbol{v})^{n}$ is given by the scheme (35), we obtain

[TABLE]

where $TB$ and $TS$ stand for the contribution of boundary and mass source to the energy production.

3.3.3 Stability

The stability of any numerical scheme ensures that errors in data at a time step are not further amplified along the next steps. To acquire the stability of our scheme, we have investigated several time-bounds $\tau_{n}$ and we introduced an artificial viscosity in the scheme to force the monotony of the discrete energy.

The artificial viscosity is a known method for improving the performance of a numerical scheme, see for example [36, 30]. So, we add a “viscous” term to ${\cal J}$:

[TABLE]

In our numerical experiments for the case of continuum soil and porosity functions, we found out that the following version of $\nu_{a(i,j)}(h,\boldsymbol{v})$ works very well:

[TABLE]

where

[TABLE]

and

[TABLE]

The case with discontinuous data will be discussed later.

The variation of energy is now given by

[TABLE]

The hyperbolic character of the shallow water equations and the necessity to preserve the positivity of the water depth impose an upper bound for the time step. Both requirements are satisfied if this upper bound is chosen to be

[TABLE]

where

[TABLE]

and $CFL$ is a number between [math] and $1$.

Remark 3.2**.**

An upper bound of type (44) for the time-step is well known in the theory of hyperbolic systems as the Courant–Friedrichs–Lewy condition [36, 5].

4 Model Validation

A rough classification of validation methods splits them into two classes: internal and external. For the internal validation, one analyses the numerical results into a theoretical frame: comparison to analytic results, sensibility to the variation of the parameters, robustness, stability with respect to the errors in the input data etc. These methods validate the numerical results with respect to the mathematical model and not with the physical processes; this type of validation is absolutely necessary to ensure the mathematical consistency of the method.

The external validation methods assume a comparison of the numerical data with measured real data. The main advantage of these methods is that a good consistency of data validates both the numerical data and the mathematical model. In the absence of measured data, one can do a qualitative analysis: the evolution given by the numerical model is similar to the observed one, without pretending quantitative estimations.

4.1 Internal Validation

We analyze the performance of the 1D version of our numerical scheme by considering the following test problems:

flow over a bump, 2. 2.

Thacker problem, 3. 3.

Riemann problems.

4.1.1 Flow Over a Bump

We consider two classical tests related to the steady solution and the propagation of a wave over a bump.

Steady solution

The steady solution was analyzed in [15, 11, 52]. In the absence of all frictional terms and mass sources, a steady solution of the shallow water equations is given by the solution of the algebraic equations

[TABLE]

We consider the case of a soil surface that is flat with a bump

[TABLE]

defined over the domain $\Omega=[0,L]$ of the flow with $L=25$m. The only boundary condition we use is

[TABLE]

The initial state is given by

[TABLE]

The goal is to find out whether or not the discrete dynamical system (35) with (17) and (40) reaches a steady state. If a steady state is reached, then we are interested to see the way this state is determined by the initial water level $\eta$. As in the aforementioned papers, we consider the cases corresponding to values of the left boundary flux $q_{0}=1.53$m2/s and $q_{0}=0.18$m2/s.

Figures 1 and 2 illustrate the behavior of the numerical solution when it develops a moving or a steady shock.

Propagation of a wave

The second test problem concerns the propagation of a wave over a bump. This kind of the test was firstly considered by LeVeque [37]. The soil topography is defined by

[TABLE]

over the domain of flow $\Omega=[-1.25,1.25]$. The initial datum is a flat surface perturbed by a small crest

[TABLE]

There are no boundary conditions. The water level on $\Omega$ at a given moment of time calculated with our numerical scheme is pictured in Figure 3 for $x\in[-0.5,0.5]$.

4.1.2 Thacker Problem

The oscillation of the water level in a convex soil surface or the oscillatory motion of the sea water in the shoreline are very natural phenomena. A very difficult problem here is that a singularity appears at the wet/dry contact between water and soil surfaces. Different numerical schemes were proposed to solve this kind of problem and among them we mention [6].

Thacker problem is a useful test to check the ability of a numerical scheme to catch and highlight the propagation of a wet/dry front. Moreover, this problem can be exactly solved [49, 44]. Thus, once can take advantage of its analytic solution to verify the accuracy of the numerical scheme. Thacker solution for 1-D shallow water equations was obtained in the absence of vegetation ($\theta=1$), for $\mathfrak{t}^{p}=\boldsymbol{0}$ and $\mathfrak{t}^{s}=h\tau\boldsymbol{v}$ in (1), where $\tau$ is a proportionality coefficient, for a soil surface $z(x)$ of parabolic shape, and for a water velocity which does not depend on space variable, i.e. $\boldsymbol{v}(t,x)=\boldsymbol{f}(t)$.

For the case of a 1D flow, the general solution is given by

[TABLE]

where

[TABLE]

$\lambda_{1}$ and $\lambda_{2}$ solve

[TABLE]

and

[TABLE]

In (52), $A$ and $B$ are integration constants. Depending on $z(x)$ and $\tau$, the solutions $\lambda_{1}$ and $\lambda_{2}$ can be real or complex numbers. An oscillating solution is obtained if

[TABLE]

In our numerical experiment, we chose

[TABLE]

where the constants $z_{m}$, $L$ and $u_{0}$ are free parameters.

We compare the results from our numerical scheme (using different formulas (17) or (19) for the interface values of $(\theta h)^{s}$) with the analytic Thacker solution: the evolution of the errors between exact and numerical solutions in Figure 4 and the dynamics of the water surface in Figure 5. In all simulations, we use $z_{m}=10{\rm m}$, $L=6000{\rm m}$, $u_{0}=5{\rm m}\cdot{\rm s}^{-1}$ and $\tau=0.001{\rm s}^{-1}$. The artificial viscosity is defined by (40).

4.1.3 Riemann Problem

We remind that Riemann problem for shallow water equations with topography and vegetation requires to solve the equations for piecewise constant initial data

[TABLE]

and piecewise constant soil surface and porosity functions

[TABLE]

The problem was extensively analyzed [33, 34, 4, 21, 17, 1] for the case of constant vegetation, but there are only few results when the porosity is a variable function [19, 9]. Perhaps, this scarcity of results is due to the fact that the problem is less tractable for the case of jumps in both soil and porosity functions.

Usually, the solution of the problem is sought as a bunch of simple waves (shock waves, rarefaction waves and standing waves) ordered on the base of the entropy principle, [31]. We remark that such a solution may not exist or it can exist (being unique or non-unique). We perform a comparative study on these three cases. The data for each case are given in Table 1.

The choice of these data was suggested by some cases analyzed in [34].

The analytic solutions for the Riemann problem are calculated using the path connection introduced in [25] and the algorithms described in [29].

Figure 6 summarizes the results for the Riemann problem in the three cases (non-existence, uniqueness and two solutions) associated to the initial data given in Table 1.

4.2 External Validation

Numerical data incorporate two main errors: one given by the approximation of the physical phenomena by a mathematical model and the other one given by the approximation of the solution of the mathematical model. Consequently, a good fit of the numerical data with experimental data validates both the mathematical model and the numerical scheme. Thus, for real world applications, it is essential to confront the numerical data with the experimental data.

We test the mathematical model and the numerical solution considering the following problems:

dam break flow over a triangular bump, 2. 2.

downslope flow through rigid vegetation, 3. 3.

flow on a convex-concave vegetated soil surface, 4. 4.

simulation on Paul’s Valley (2D).

The first two tests simulate solutions for experiments with measured data, while the last two tests target the qualitative behavior of the solution.

4.2.1 Dam Break Flow Over a Triangular Bump

First, we consider the CADAM test case of the dam break flow propagation in an impermeable channel with a symmetric triangular bump [40, 14, 23]. The layout of this experiment and the initial state of the water at rest are presented in Figure 7.

Several gauges $G_{i}$ are placed in the channel to measure the water depth evolution once the gate retaining the water volume from the left side is suddenly removed.

In this experiment, the cover plant is absent: $\theta=1$. For the frictional term $\mathfrak{t}^{s}_{a}=\alpha_{s}(h)|\boldsymbol{v}|v_{a}$, we consider the Darcy-Weisbach formula

[TABLE]

where $\tau$ is constant. For this setup, we performed several numerical simulations with different values for $\tau$. The experimental data are extracted from the graphics reported in [23], Figure 5. The results are synthesized in Figure 8.

Note that our solution is in agreement to the measured one (for all four gauges) for $\tau\approx 2e-3$.

We repeated the numerical simulations using Stickler-Manning formula

[TABLE]

with $n^{2}\approx 0.003$m1/3 and we have seen that there is no significant difference from the previous case.

4.2.2 Flow Over a Slope with Vegetation

The flow in the presence of vegetation is more complex than the one on bare soil and perhaps this explains a scarcity of the reported data.

Here, we compare our results with the ones reported in [13]. Briefly, the experimental installation consists of an $18$m long and $1$m width laboratory flume with a longitudinal bottom slope $S=1.05{\rm mm}/{\rm m}$. The vegetation is modelled using uniformly distributed emergent circular cylinders of radius $R_{c}=5$mm. The density of vegetation is $N=81$ cylinder/m2. This type of vegetation allows the estimation of the porosity by

[TABLE]

To obtain the parameters $\alpha_{p}$ and $\alpha_{s}(h)$, we use the steady state values of $h$ and $u$ reported [13] for different values of the water discharge $q=hu$. These values can be theoretically estimated from relations

[TABLE]

One can estimate $\alpha_{s}(h)$ from the data for bare soil and then use it to calculate $\alpha_{p}$ using the data for uniformly distributed vegetation. Thus, if Darcy-Weisbach formula is used, one finds that $\tau$ runs in $[6.508084,9.00375]\cdot 10^{-3}$ and $\alpha_{p}$ takes values in $[73.20,80.20]$. For the case of Strickler-Manning formula, one finds $n^{2}\in[2.893622,3.178546]\cdot 10^{-3}$ and $\alpha_{p}\in[72.27832,77.97572]$.

These estimations indicate that Strickler-Manning formula is more appropriate than Darcy-Weisbach formula to be used for this kind of flow.

For testing numerical data, we choose the experiment of water flow on a partial vegetated bottom surface. We use Strickler-Manning formula with the parameters

[TABLE]

These values correspond to the Strickler $K=57.099$ and drag coefficient $C_{D}=1.1738$.

Figure 9 includes the numerical and experimental data for six different steady configurations. The experimental data are extracted from the graphics reported in [13], Figure 8.

4.2.3 Flow Over a Convex-Concave Vegetated Soil

Surface

There are two reasons for considering the flow over a convex-concave vegetated soil surface: to highlight the influence of the cover plant on the water flow dynamics and to emphasize the ability of our numerical scheme not only to preserve the lake as an equilibrium point, but to asymptotically reach it.

We consider a topography with the surface given by a fourth degree polynomial on $[-L,L]$

[TABLE]

where $L=1500$m, $h_{m}=50$m and $h_{b}=5$m. For this configuration, the mathematical equilibrium point may consist of one or two lakes. Our model is run for two different configurations, without vegetation and with vegetation, respectively. We use Stickler-Manning formula (54) with $n^{2}=0.003{\rm m}^{1/3}$. The initial conditions

[TABLE]

are identical for both configurations. For the configuration with vegetation, we add

[TABLE]

Figure 10 illustrates the numerical solution of the water surface for these two configurations.

4.2.4 Simulation on Paul’s Valley

Unfortunately, we do not have data for the water distribution, plant cover density and measured velocity field in a hydrographic basin to compare our numerical results with. However, to be closer to reality, we have used GIS data for the soil surface of Paul’s Valley and accomplished a theoretical experiment: starting with a uniform water depth on the entire basin and using different cover plant densities, we run the 2D version of our scheme (a flow modulus of ASTERIX) on a hexagonal network. To generate an arbitrary hexagonal network from an arbitrary raster GIS data, we use the method introduced in [26]. Figure 11 shows that the numerical results are consistent with direct observations concerning the water time residence in the hydrographic basin.

5 Discussion, Conclusions and Further Research

The effect of vegetation on the water circulation in a continuum soil-plant-atmosphere can relatively be easy described, but it is hardly quantifiable. The canopy of vegetation intercept and store a certain amount of rain, the stem of plant generate a resistance force to the water running on the soil surface and the plant roots modify the water infiltration rate into soil. To quantify these physical processes and put them in a mathematical model is a big challenge and it requires new mathematical concept, experimental results and numerical simulations.

The mathematical model (1) is simplified as possible as can be to retain essential interaction factors and to be mathematical tractable. The mathematical model (1) belongs to the class of hyperbolic systems and it has certain distinctive feature that raises special difficulties when one elaborates a numerical scheme to approximate its solution. Namely,

(a) the system becomes singular for $h=0$;

(b) the hyperbolic system has a non-conservative form and its coefficients depend on the space variable, non-homogeneous hyperbolic system;

(c) there can be a dependence of the system of the derivatives of discontinuous functions, soil surface and/or porosity functions;

(d) the unknown water depth function $h(t,x)$ must be a positive function.

We appreciate that the finite volume method for the space discretization of the system is the most fitted method. It allows one to treat the non-conservative product, see (11), and then perform the approximations (14-3), (17), (18) and (19). By introducing the flux terms (15-1), (15-2), (16), we were able to ensure the positivity of $h$.

To solve the problem of the singularities introduced by $h=0$, we perform a first order fractional time step method. In the presence of frictional force or drag resistance force, the second equation in the scheme (34) is nonsingular even if $h=0$.

When validating using Riemann problem, Section 4.1.3, we should observe that for defining the standing waves, one needs to augment the shallow water equations with an additional principle: either path connection or an energetic principle.

For us, testing the numerical scheme against the solution of Riemann problem, given by a composition of simple waves, is instructive in two ways. In one way, it is interesting to find out how the numerical scheme works in all the possible situations: non-existence, uniqueness and non-uniqueness of the solutions. On the other way, the numerical scheme does not “see” the second principle invoked in the construction of the standing waves since it is based on the approximation principles of the operators and functions appearing in the shallow water equations. Thus, when one compares the numerical solution to the analytic solution, the question of compatibility between the two kind of principles can be raised. However, the case of the Godunov type approximation schemes is different since these schemes are obtained on the basis of the (exact or approximate) solution of the Riemann problem. We observed that it does not matter whether the quantity $(\theta h)^{s}_{(i,j)}$ is calculated using (18) or (19) in the cases of non-existence or uniqueness of the solution. But in the case of multiple solutions, different approximation formulas of the gradient of the free surface lead to different solutions. Also, we draw attention that by non-existence case we mean ”the non-existence of a composite wave as a solution of the Riemann Problem”, which is different from ”the non-existence of a solution of the Riemann Problem”. Consequently, the numerical solution can be an approximation of the true solution or can be a fictitious mathematical object.

Note that the transition from the non-existence to the existence but not uniqueness of the solution for Riemann problem is controlled only by the size of the jump in $\theta$.

The numerical simulations for the steady solution case of a flow over a bump (Section 4.1.1) suggest the following conjecture.

Conjecture**.**

There is a value $q^{M}$ such that for any $q\in(0,q^{M})$ there are two threshold values $0.2\leq\eta_{1}(q)<\eta_{2}(q)$ such that:

the numerical solution asymptotically reaches a state $(h^{\star},u^{\star})$; 2. 2.

the state $(h^{\star},u^{\star})$ casts into one of the following three distinct categories:

1. (a)

for any $\eta\in(0.2,\eta_{1})$, it develops a moving shock which goes out of the domain $\Omega$ after a while, 2. (b)

for any $\eta\in(\eta_{1},\eta_{2})$, it has a steady shock located in the domain occupied by the bump, 3. (c)

for any $\eta>\eta_{2}$, it is smooth (no shock).

Case (a) of the conjecture is illustrated on the first row of pictures in Figures 1 and 2; case (b) is also illustrated on the second row of pictures in the same figures.

A common problem in many countries, Romania being a typical example, is land deforestation with severe flood effects. The test proposed in Section 4.2.3 was considered to emphasize the protective role of vegetative barrier against such natural disasters. Remark that in the configuration with vegetation, the water does not anymore cross the knoll since it loses part of its energy due to the friction with plants, see Figure 10. Note that our numerical method is not only well-balanced, but it is capable to produce solutions that reach different equilibrium states (depending on the initial data and the values of the parameters). It should be also remarked that the porosity $\theta$ plays the role a control function to get different equilibrium points: one or two lakes.

Water flow in hydrographic basins is more complicated due to variations in soil topography and in vegetation diversity. In contrast to the 1D cases where we were able to easily simulate the soil data, we now need 2D land and vegetation data to be closer to reality. Data for soil surface can be obtained from topographic maps, GIS, LiDAR etc. For the purpose of modeling, one needs an interface to transfer data between different scales. In this sense, we used the tool built and described in [26]. Simulation on Paul’s Valley (with no rain and infiltration) was considered in Section 4.2.4 as a first step to get closer to the interesting cases from ecology, biology, hydrology and terrestrial sciences. An important variable for such a hydrographic basin is the amount of water leaving the basin. This variable can be measured using the relative amount $q$ of water in the basin at the moment of time $t$ given by

[TABLE]

Figure 12 shows the results for the water content in Paul’s Valley basin obtained with our models ASTERIX and with CAESAR-Lisflood-OSE [24].

A general issue relates to whether higher cover plant densities can prevent soil erosion and flood or not. Both pictures show that if the cover plant density is increasing then the decreasing rate $\dot{q}$ of $q$ is smaller. We can think at a “characteristic velocity” of the water movement in the basin and this velocity is in a direct relation with $\dot{q}$. We can now speculate that smaller values of $\dot{q}$ imply softer erosion processes.

This valley belongs to Ampoi’s catchment basin. Flood generally appears when the discharge capacity of a river is overdue by the water coming from the river catchment area. Our pictures show that higher cover plant densities imply smaller values of $\dot{q}$ which in turn give Ampoi River the time to evacuate the water amount flowing from the valley.

As conclusions, we have presented a mathematical model able to take into account the effects introduced by plants on the water flow on soil surface and we have introduced a numerical scheme to approximate its solution. The numerical scheme is of first order in time and space accurate and is able to model a large class of physical phenomena, waves propagation, backward wave reflection induced by obstacles or filter band vegetation, oscillating flow, one or multiple lake formation. The scheme is relatively simple to be implemented in a more complex model like soil erosion.

As current research, we work on elaborating a methodology based on inverse method to estimate the hydrodynamics parameters $\alpha_{p}$ and $\alpha_{s}$ from in situ experimental measurements.

In the future, we intend to build up a coupled model – shallow water - soil erosion – to investigate the role of plants to mitigate the transport of the polluted soil particles from contaminated area to around area.

Acknowledgment

This work was partially supported by a grant of the Ministry of Research and Innovation, CCCDI-UEFISCDI, project number PN-III-P1-1.2-PCCDI-2017-0721/34PCCDI/2018, within PNCDI III.