<figure><img src="image_1.png" /></figure>

2007 International Conference in Honor of Claude Lobry

# On a 1D-Shallow Water Model: Existence of solution and numerical simulations.

Olivier Besson<sup>a</sup> — Soulèye Kane<sup>b</sup> — Mamadou Sy<sup>c,*</sup>

<sup>a</sup> Univ. Neuchâtel, Inst. Math., 11, R. Emile Argand, CH-2007, Neuchâtel, Suisse.
Olivier.Besson@unine.ch,

<sup>b</sup> UCAD. Dakar, Dept. Math.-Info, Avenue Cheikh Anta Diop, Dakar, Sénégal.
Souleye.Kane@unine.ch,

<sup>c</sup> Univ. Gaston Berger de Saint-Louis, UFR SAT, LANI, BP 234 Saint-Louis, Sénégal.

* Corresponding author
Mamadou.Sy@ugb.sn

**ABSTRACT.** The study of a 1D-shallow water model, obtained in a height-flow formulation, is presented. It takes viscosity into account and can be used for the flood prediction in rivers. For a linearized system, the existence and uniqueness of a global solution is proved. Finally, various numerical results are presented regarding the linear and non linear case.

**RÉSUMÉ.** Nous dérivons les équations de Saint-Venant complètes avec la formulation hauteur-débit. La viscosité est prise en compte dans le modèle. Pour le système linéarisé, l'existence et l'unicité de solution globale sont montrées. Des résultats numériques sont présentés aussi bien dans le cas linéaire que non linéaire.

**KEYWORDS :** Shallow water, river, flood, height-flow, Saint-Venant.

**MOTS-CLÉS :** Equation de Saint-Venant, rivière, crue, hauteur-débit.
# 1. Introduction

The geophysical flows in thin domains are described by the shallow water system, see for example [2], [3], [16], and [22]. The classical shallow water system where the viscosity is neglected is of current use. But the mathematical study of this model seems to be not so satisfactory, see e.g. [9].

Most of the developed models use a section-flow formulation or a height-velocity one. For example in [9], the height-velocity formulation gives a shallow water system including frictions, viscosity and the Coriolis-Boussinesq coefficient. This model is numerically validated and satisfactory results are obtained for the dam break case.

In [5], a global existence result of weak solutions in a $\mathbb{R}^2$ periodic domain is given. The reader interested with the derivation of the shallow water equations can consult [2], [3], and [17]. The model used in our paper is heavily inspired by [2]. The considered flow is one dimensional, in a channel with a parallelepipedic cross-section. The derivation of a shallow water system for channels with other shape or for other kind of channel are presented in [11].

This paper is divided into three parts. In the first part, the derivation of a shallow water system in one dimension, using a height-flow formulation, is presented. Then some existence results of solutions are given in the second one. Finally numerical results for the flow in a channel with different scenarios are presented.

# 2. The model

River flows are governed by the shallow water equations. These equations are obtained using some lateral integration (i.e. in width and depth of the river) of the incompressible Navier-Stokes equations.

The considered Navier-Stokes system has the following form.

$$ \frac{\partial u}{\partial t} + (u \cdot \nabla)u - \operatorname{div}\sigma = f \quad \text{in } (0, T) \times \Omega, \quad (1) $$

$$ \operatorname{div} u = 0 \quad \text{in } (0, T) \times \Omega, \quad (2) $$

with the initial condition

$$ u(0, x) = u_0(x) \quad \text{in } \Omega, \quad (3) $$

and the boundary conditions (see fig. 1)

$$ \sigma n = 0 \quad \text{on } (0, T) \times \Gamma_s, \quad (4) $$

$$ u \cdot n = 0 \quad \text{on } (0, T) \times \Gamma_b, \quad \text{and} \quad (5) $$

$$ \sigma n \cdot \tau_i = f_{b,i} \quad \text{on } (0, T) \times \Gamma_b, \quad i = 1, 2, \quad (6) $$

where $u = (u_1, u_2, u_3)$ is the velocity field, $\sigma = -p + \frac{\mu}{2}(\nabla u + \nabla u^t)$ is the stress tensor, $p$ is the pressure, $f$ represents the external forces, and $f_b$ is some traction vector
on the river bed. Instead of using traction conditions on the bottom of the river, one can also use a homogeneous Dirichlet condition:

$$u = 0 \quad \text{on } (0, T) \times \Gamma_b. \qquad (7)$$

<figure><img src="image_2.png" /><figcaption>**Figure 1. River section**</figcaption></figure>

The physical nature of rivers allows to simplify these equations. Indeed they are characterized by the predominance of their length $L$, compared to the width $l$ and the depth $H$. This predominance induces the flow direction and other phenomena described by the following different assumptions.

## 2.1. Assumptions

The following assumptions are used in our model.

1) The depth and the width are very small compared to the river length. So the flow is essentially one dimensional and is parallel to the domain walls and bottom. Thus when the stream line curvature is small, the vertical and lateral acceleration are negligible compared to the longitudinal one and the pressure distribution is hydrostatic.

2) The variation of the water height is very small.

3) The geometry of the domain is fixed, thus deposit and leaching effects of the sediment are neglected.

4) The slope of the domain bank is very small.

5) The friction effects on the shore will be estimated by some heuristic formula like the Manning or the Chezy one.

6) The domain is assumed to be rectilinear. This allows to consider a parallelepipedic domain with length $L$ and width $l$. The water surface for each section of the river is assumed to be horizontal. The extra height effects in a sloping domain are not taken account in the analysis and are assumed to have very small influence on the result.

7) The momentum and the energy fluxes all along the river section resulting from the non uniformity of the velocity distribution will be estimated using a mean velocity and some correction coefficients which are function of the position along the water course and of the height.

8) The gravity force is the only one taken into account. So the influence of the Coriolis force is neglected.

These assumptions allow some simplifications in the system (1-7).
## 2.2. Simplified equations

Let $U_2$ (resp. $U_1, U_3$) represents the characteristic velocity in the flow direction, (resp. in the other directions). If $T$ is a characteristic time, then $U_1 = l/T$, $U_2 = L/T$, and $U_3 = H/T$.

Let $t = Tt'$, $u_1 = U_1v_1$, $u_2 = U_2v_2$, $u_3 = U_3v_3$, $x = lx'$, $dx = ldx'$, $y = Ly'$, $z = Hz'$, $dz = Hdz'$, $p = p'U_2^2$. The variables $t'$, $x'$, $y'$, $z'$, $p'$, $v_1$, $v_2$, $v_3$ are without dimension. Taking into account the previous assumptions, using variables without dimension, redefining $x' = x$, $y' = y$, $z' = z$, $t' = t$, $p' = p$ and assuming that the mainstream is along the $y$-axis, we get for the equation (1):

$$ \frac{\partial v_2}{\partial t} + v_2 \frac{\partial v_2}{\partial y} - \left( \frac{L}{lU_2^2} \frac{\partial \sigma_{12}}{\partial x} + \frac{1}{U_2^2} \frac{\partial \sigma_{22}}{\partial y} + \frac{L}{lU_1^2} \frac{\partial \sigma_{32}}{\partial z} \right) = f, \quad (8) $$

$$ \frac{\partial p}{\partial y} = 0, \quad (9) $$

$$ \frac{U_2^2}{\rho l} \frac{\partial p}{\partial z} = -g, \quad (10) $$

$$ \operatorname{div} v = 0. \quad (11) $$

We have

$$ \sigma_{22} = -p + \frac{\mu U_2}{L} \frac{\partial v_2}{\partial y}, $$

$$ \sigma_{12} = \mu \left( \frac{U_2}{l} \frac{\partial v_2}{\partial x} + \frac{U_1}{L} \frac{\partial v_1}{\partial y} \right), \quad \sigma_{32} = \mu \left( \frac{U_2}{H} \frac{\partial v_2}{\partial z} + \frac{U_3}{L} \frac{\partial v_3}{\partial y} \right). $$

The reader can consult [22] for more informations. Let $\Omega(y, t)$ be the river section:

$$ \Omega(y,t) = \{(x, y, z) \in \mathbb{R}^3 \text{ s.t. } 0 \le x \le l, \ 0 \le z \le h(y, t)\}. $$

The one dimensional shallow water model is obtained from equations (8)-(11) by integration over $\Omega(y, t)$.

The equation (11) give

$$ v_3(x, y, h) = - \int_0^{h(y,t)} \partial_x v_1 \, dz - \int_0^{h(y,t)} \partial_y v_2 \, dz + v_3(x, y, 0). $$

We know that $v_3(x, y, h) = \frac{dh}{dt}$.

Thus

$$ \int_0^l \frac{\partial h}{\partial t} dx = \int_0^h v_1(x=0) dz - \int_0^h v_1(x=l) dz - \int_0^l \int_0^h \partial_y v_2 dz dx + \int_0^l v_3(x, y, 0) dx. $$

We define $Q(t, y) = \int_0^l \int_0^{h(t,y)} v_2 dx dz$, the flow through $\Omega(y, t)$. Then we have

$$ l \frac{\partial h}{\partial t} + \frac{\partial Q}{\partial y} = h v_1(y=0) - h v_1(y=l) + l v_3(x, y, 0). $$
We assume for this model that $v_1$ is the same to the left bank and right bank; and the infiltration is null.

The first shallow water equation will be

$$l \frac{\partial h}{\partial t} + \frac{\partial Q}{\partial y} = 0.$$

The equation (8) gives

$$\int_{\Omega} \partial_t v_2 \, dx dz + \int_{\Omega} v_2 \partial_y v_2 \, dx dz - \frac{L}{U_2^2 l} \int_{\Omega} \partial_x \sigma_{12} \, dx dz - \frac{1}{U_2^2} \int_{\Omega} \partial_y \sigma_{22} \, dx dz - \frac{L}{U_2^2 l} \int_{\Omega} \partial_z \sigma_{32} \, dx dz = \int_{\Omega} f \, dx dz.$$

If we denote by $I_i$, $i = 1, 2, 3, 4, 5$ the left integrate terms of the previous equality, we have

$$I_1 = \int_0^l \int_0^{h(y,t)} \partial_t v_2 \, dx dz = \int_t^l \int_0^{h(y,t)} v_2 \, dx dz - \int_0^l v_2(z=h) \partial_t h \, dx.$$

Then

$$I_1 = \partial_t Q + v_2(z = h) \partial_y Q.$$

$$I_2 = \int_0^l \int_0^h v_2 \partial_y v_2 \, dx dz = \frac{1}{2} \partial_y \int_0^l \int_0^h v_2^2 \, dx dz - \frac{1}{2} v_2^2(z = h) \partial_y h.$$

$$I_3 = \int_0^l \int_0^h \partial_x \sigma_{12} \, dx dz = \int_0^h (\sigma_{12}(x=l) - \sigma_{12}(x=0)) \, dz.$$

$$I_5 = \int_0^l \int_0^h \partial_z \sigma_{32} \, dx dz = \int_0^l (\sigma_{32}(z=h) - \sigma_{32}(z=0)) \, dx.$$

$$
\begin{align*}
I_4 &= \int_0^l \int_0^h \partial_y \sigma_{22} \, dx dz \\
&= \partial_y \int_0^l \int_0^h \sigma_{22} \, dx dz - \int_0^l \sigma_{22}(z=h) \partial_y h \, dx \\
&= \partial_y \int_0^l \int_0^h -p + \frac{\mu U_2}{L} \partial_y v_2 \, dx dz - \int_0^l \sigma_{22}(z=h) \partial_y h \, dx \\
&= -\partial_y \int_0^l \int_0^h p \, dx dz + \frac{\mu U_2}{L} \partial_y^2 \int_0^l \int_0^{h(y,t)} v_2 \, dx dz - \frac{\mu U_2}{L} \partial_y \int_0^l v_2(z=h) \partial_y h - \int_0^l \sigma_{22}(z=h) \partial_y h \, dx \\
&= -\partial_y \int_0^l \int_0^h p \, dx dz + \frac{\mu U_2}{L} \partial_y^2 Q - \frac{\mu U_2 l}{L} \partial_y (v_2(z=h) \partial_y h) - \int_0^l \sigma_{22}(z=h) \partial_y h \, dx.
\end{align*}
$$
Let $\mathcal{T}$ the tangent of the surface. We have $\mathcal{T} = \begin{pmatrix} 0 \\ 1 \\ \partial_y h(y, t) \end{pmatrix}$ and the normal of $\mathcal{T}$

is denoted by

$$n = \begin{pmatrix} 0 \\ -\partial_y h(y, t) \\ 1 \end{pmatrix}$. The adimensional normal is recall $n = \begin{pmatrix} 0 \\ \frac{-l}{L} \partial_y h(y, t) \\ 1 \end{pmatrix}$.

Using (4), we have

$$\sigma_{22}(z=h)\partial_y h = \frac{L}{l}\sigma_{23}(z=h). \qquad (12)$$

The normal on the bottom is $n = \begin{pmatrix} 0 \\ 0 \\ -1 \end{pmatrix}$. Thus $\tau_{b,2} = \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}$.

The assumption (4) gives

$$f_{b,2} = -k(v_2)v_2.$$

Then, we get with the equation (6) $\sigma_{23}(z=0) = kU_2v_2$.

If $y$ is fixed, the stress tensor is the same in $x=l$ and $x=0$. Thus $I_3 = 0$.

The equation (10) gives

$$p = -\frac{l}{U_2^2}g(z-h) + \bar{p}_a$$

with $\bar{p}_a$ is the adimensional atmospheric pressure. From where, the first term of $I_4$

$$\partial_y \int_0^l \int_0^h p dxdz = \frac{l^2}{2U_2^2}g \partial_y h^2 + \bar{p}_a l \partial_y h.$$

Using (12), the last term of $I_4$ will be

$$\int_0^l \sigma_{22}(z=h) \partial_y h \, dx = \frac{L}{l} \int_0^l \sigma_{23}(z=h) \, dx.$$

Let us gather the terms of momentum equation, we will have the second full equation of shallow water

$$\partial_t Q - \frac{\mu}{U_2 L} \partial_y^2 Q + \partial_y (l v_2^2 h) + \bar{p}_a \frac{l}{U_2^2} \partial_y h + \frac{l^2 g}{2U_2^2} \partial_y h^2 + \frac{\mu l}{U_2 L} \partial_y (v_2 \partial_y h) + \frac{Lk}{U_2 l h} Q = \int_{\Omega} f. \qquad (13)$$

Or

$$\partial_t Q - \frac{\mu l}{U_2 L} \partial_y (h \partial_y v_2) + \partial_y (l v_2^2 h) + \bar{p}_a \frac{l}{U_2^2} \partial_y h + \frac{l^2 g}{2U_2^2} \partial_y h^2 + \frac{Lk}{U_2 l h} Q = \int_{\Omega} f. \quad (14)$$

The assumptions 1. and 2. show that the term $\frac{\mu l}{U_2 L} \partial_y (v_2 \partial_y h)$ of equation (13) is negligible. If we linearize the term $\partial_y (l v_2^2 h)$ by $lU_2^2 \partial_y (h)$, we obtain for the momentum equation

$$\partial_t Q - \nu \partial_y^2 Q + \beta(h) \partial_y h + K(Q) Q = F_{ext} \qquad (15)$$

where the coefficient $K$ is a friction coefficient, $\nu = \frac{\mu}{U_2 L}$ is the renormalized viscosity,

$$\beta(h) = a + 2bh,$$
with

$$a = lU_2^2 + \frac{\bar{p}_a l}{U_2^2}, \quad b = \frac{gl^2}{2U_2^2}.$$

Let $W = (0, T) \times I$ with $I = (0, L)$.

Assume that $K$ is constant, the shallow water system that we want to study follows

$$\begin{align*}
& l \frac{\partial h}{\partial t} + \frac{\partial Q}{\partial x} = 0 \quad \text{in } W, \\
& \frac{\partial Q}{\partial t} - \nu \frac{\partial^2 Q}{\partial x^2} + \beta(h) \frac{\partial h}{\partial x} + KQ = 0 \quad \text{in } W, \\
& Q(t, 0) = Q_e(t), \quad Q(t, L) = Q_s(t), \\
& h(0, x) = h_0(x), \quad Q(0, x) = Q_0(x).
\end{align*}$$

If $Q = q + Q_s + \frac{1}{L}(L-x)(Q_e - Q_s)$, then our 1-D shallow water model is

$$l \frac{\partial h}{\partial t} + \frac{\partial q}{\partial x} = f_1 \quad \text{in } W, \qquad (16)$$

$$\frac{\partial q}{\partial t} - \nu \frac{\partial^2 q}{\partial x^2} + \beta(h) \frac{\partial h}{\partial x} = f_2 \quad \text{in } W, \qquad (17)$$

$$q(t, x) = 0 \quad \text{on } \partial I, \qquad (18)$$

$$h(0, x) = h_0(x), \qquad (19)$$

$$q(0, x) = q_0(x), \qquad (20)$$

with

$$f_1(t, x) = \frac{1}{L}(Q_e(t) - Q_s(t)), \qquad (21)$$

$$f_2(t, x) = -Q'_s(t) - \frac{1}{L}(L-x)(Q'_e(t) - Q'_s(t)), \qquad (22)$$

and

$$q_0(x) = Q_0(x) - Q_s(0) - \frac{1}{L}(L-x)(Q_e(0) - Q_s(0)). \qquad (23)$$

### 3. Existence of solution for the linearized problem ($\beta$ is constant)

In this section, we give some existence and uniqueness for solutions of equations (16)-(20) assuming $\beta$ is constant. For a given $q$, the hight $h$ is obtained from equation (16) using some time integration:

$$h(t, x) = h^0(x) - \frac{1}{l} \int_0^t \left( \frac{\partial q}{\partial x} - f_1 \right) ds. \qquad (24)$$

Moreover, if $h$ is given, a weak formulation of equation (17) is
$$ \int_I \frac{\partial q}{\partial t} v dx + \nu \int_I \frac{\partial q}{\partial x} \frac{\partial v}{\partial x} dx - \int_I (ah + \frac{1}{2}bh^2) \frac{\partial v}{\partial x} dx = \int_I f_2 v \, dx. \quad (25) $$

Equation (24) makes sense when $q \in L^1(0, t; W^{1,1}(I))$. In the same way, for equation (25), we need $v \in H_0^1(I)$, and $h \in L^2(0, T; L^4(I))$. Let us define the following space

$$ V = L^2(0, T; H_0^1(I)) \cap C(0, T; H^{-1}(I)), $$

and let $F: V \to V$ be the application defined in the following way:

1) For $q$ in $V$, let $h = h(q)$ be the solution of (24).

2) Then for $h$ given by 1), let $\bar{q} = F(q)$ where $\bar{q}$ is the solution of (25).

We show that the problem (24)-(25) has a unique solution in $V \times H^1(0, T, L^2(I))$ when $\beta$ is constant.

**Theorem 3.1** Let $0 < T < \frac{\nu}{\beta}$, there exists a constant $0 < C(T) < 1$ such that

$$ \forall q_1, q_2 \in V, \quad \|F(q_1) - F(q_2)\|_V \le C(T) \|q_1 - q_2\|_V, $$

i.e. $F$ is a contraction.

**proof**

Let $v \in H_0^1(I)$, $\bar{q}_i = F(q_i)$, and $h_i$, $i = 1, 2$, where $h_i$ is the solution of (24) associated to $q_i$. We have

$$ \int_I \frac{\partial \bar{q}_i}{\partial t} v dx + \nu \int_I \frac{\partial \bar{q}_i}{\partial x} \frac{\partial v}{\partial x} dx - \beta \int_I h_i \frac{\partial v}{\partial x} dx = \int_I f_2 v \, dx, \quad i = 1, 2. $$

So

$$ \int_0^t \int_I v \frac{\partial}{\partial t} (\bar{q}_1 - \bar{q}_2) dx ds + \nu \int_0^t \int_I \frac{\partial}{\partial x} (\bar{q}_1 - \bar{q}_2) \frac{\partial v}{\partial x} dx ds - \beta \int_0^t \int_I (h_1 - h_2) \frac{\partial v}{\partial x} dx ds = 0, $$

but

$$ \beta \int_0^t \int_I (h_1 - h_2) \frac{\partial v}{\partial x} dx ds = \frac{\beta}{l} \int_0^t \int_I \frac{\partial v}{\partial x} \int_0^s \left( \frac{\partial q_1}{\partial x} - \frac{\partial q_2}{\partial x} \right) d\tau dx ds. $$

For $v = \bar{q}_1 - \bar{q}_2$, we get

$$ \frac{1}{2} \int_0^t \frac{d}{ds} \int_I (\bar{q}_1 - \bar{q}_2)^2 dx ds + \nu \int_0^t \int_I \left( \frac{\partial}{\partial x} (\bar{q}_1 - \bar{q}_2) \right)^2 dx ds \\ = \frac{\beta}{l} \int_0^t \int_I \left( \frac{\partial}{\partial x} (\bar{q}_1 - \bar{q}_2) \left( \int_0^s \frac{\partial}{\partial x} (q_1 - q_2) d\tau \right) \right) dx ds, $$

and for $t = T$

$$ \frac{1}{2} \|\bar{q}_1(T) - \bar{q}_2(T)\|_{L^2(I)}^2 + \nu \|\bar{q}_1 - \bar{q}_2\|_V^2 \le \frac{\beta}{l} T \|\bar{q}_1 - \bar{q}_2\|_V \|q_1 - q_2\|_V, $$

so

$$ \|\bar{q}_1 - \bar{q}_2\|_V \le C(T) \|q_1 - q_2\|_V $$

with $C(T) = \frac{\beta}{l\nu}T < 1$.
**Theorem 3.2** Let $\beta$ be a positive constant, and let $T > 0$, $f_1 \in L^1(0, T; L^2(I))$, and $f_2 \in H^{-1}(0, T; H_0^1(I)) \cap L^2(0, T; H^{-1}(I))$.

The problem (24)-(25) admits a unique solution $(q, h) \in V \times H^1(0, T; L^2(I))$.

### Proof

For $q$ given in $V$ we have

$$h(t, x) = h^0(x) - \frac{1}{l} \int_0^t \left( \frac{\partial q}{\partial x} - f_1 \right) ds,$$

so $h \in H^1(0, T; L^2(I))$. For this $h$, equation (25) admits a unique solution $\bar{q} = F(q)$. From theorem 3.1, the application $F$ verifies the hypothesis of the Banach fixed point theorem for $T = T^* < \frac{l\nu}{\beta}$. Thus there is a unique $(q, h) \in V \times H^1(0, T; L^2(I))$, solution of (24), (25). The global existence for any $T$ is then obtained by continuity.

*Remark*: The non linear case is a open problem.

In the follow, we give a numerical result in the case $\beta$ is constant and no constant.

# 4. Numerical results

## 4.1. Introduction

In this section we will give some numerical simulation of the following problem:

$$l \frac{\partial h}{\partial t} + \frac{\partial q}{\partial x} = f_1 \quad \text{dans } W, \qquad (26)$$

$$\int_I \frac{\partial q}{\partial t} v \, dx + \nu \int_I \frac{\partial q}{\partial x} \frac{\partial v}{\partial x} \, dx - \int_I (ah + \frac{1}{2}bh^2) \frac{\partial v}{\partial x} \, dx + \int_I \kappa q v \, dx = \int_I f_2 v \, dx. \quad (27)$$

A Lax-Friedrichs scheme is used to solve equation (26). Equation (27) is solved using a finite element method in space and an implicit Euler scheme in time.

The following methodology is used: if the hight $h$ and the flow $q$ is known at a given time step, the hight at the next step is solved using equation (26), then the flow is obtained using equation (27).

Let us introduce the following notations.

- $L$ is the length of the domain,

- $l$ is its width,

- $T$ is the final time,

- $\nu$ is the viscosity,

- $\kappa$ is the friction coefficient,

- $n$ is the number of spatial steps,

- $m$ is the number of time steps,

- $\delta = L/(n + 1)$ is the space step,
- $\tau = T/m$ is the time step,

- $x_i = i\delta$ is the i-th node,

- $t_k = k\tau$ is the k-th time step,

- $h_i^k$ is the approximation of $h(t_k, x_i)$ at time $t_k$ and position $x_i$,

- $q_i^k$ is the approximation of $q(t_k, x_i)$ at time $t_k$ and position $x_i$.

## 4.2. Approximate problem

Let $V_\delta$ be the usual subspace of $H^1(I)$ associated with finite elements of type $\mathbb{P}_1$:

$$V_{\delta} = \{v \in C^0(\bar{I}) ; \forall i \in \{0, \dots, n\}, v|_{[x_i, x_{i+1}]} \in \mathbb{P}_1\},$$

and set

$$V_{\delta,0} = V_{\delta} \cap H_0^1(I).$$

Let also $\phi_i$, $0 \le i \le n+1$ be the usual basis of $V_\delta$; the set $\{\phi_i, 1 \le i \le n\}$ is the canonical basis of $V_{\delta,0}$.

For equation (27), the following Lax-Friderichs scheme is used (see [12]):

$$h_i^{k+1} = \frac{1}{2}(h_{i-1}^k + h_{i+1}^k) - \frac{\tau}{2l\delta}(q_{i+1}^k - q_{i-1}^k). \quad (28)$$

The approximate solution $h$ of (24) is then interpolated as

$$h(x, t_k) = \sum_{i=0}^{n+1} h_i^k \phi_i(x).$$

For the numerical solution of equation (27), let

$$q(x, t) = \sum_{i=1}^{n} q_i(t) \phi_i(x).$$

For $v = \phi_j$, ($j \in \{1, \dots, n\}$), if $q_i^k$ is the approximation of $q(x_i, t_k)$, the implicit Euler scheme for (27) is

$$\sum_{i=1}^{n} q_i^{k+1} \int_I \phi_i \phi_j \, dx + \sum_{i=1}^{n} \nu \tau q_i^{k+1} \int_I \phi_i' \phi_j' \, dx + \sum_{i=1}^{n} q_i^{k+1} \int_I \tau \kappa \phi_i \phi_j \, dx = \\ \sum_{i=1}^{n} q_i^k \int_I \phi_i \phi_j \, dx + \sum_{i=1}^{n} \int_I \tau (ah(x, t_k) + bh^2(x, t_k)) \phi_j' \, dx + \int_I f_2 \phi_j \, dx.$$

Let $q^k = (q_1^k, \dots, q_n^k)^t$, $A = (A_{ij})_{1 \le i,j \le n}$, $B = (B_{ij})_{1 \le i,j \le n}$, and $L^k = (L_j^k)_{1 \le j \le n}$, with

$$A_{ij} = \int_I (1 + \tau \kappa) \phi_i \phi_j \, dx + \nu \tau \int_I \phi_i' \phi_j' \, dx$$
$$B_{ij} = \int_{I} \phi_{i} \phi_{j} \, dx,$$

$$L^k = \int_I \tau ah(x, t_k) + bh^2(x, t_k))\phi_j' \, dx + \int_I f_2 \phi_j \, dx.$$

Then the approximate solution of $q$ is given as usual by the solution of the linear system.

$$Aq^{k+1} = Bq^k + L^k. \qquad (29)$$

### 4.3. Some numerical results

In this section, some numerical results are presented in the following situation. A canal with length $L = 1m$, and width $l = 5cm$ is considered. This canal is subdivided by a dam into two parts with length $20cm$, and respectively $80cm$. The final time is set to $T = 5s$. So $b = l^2 \cdot g \cdot T^2/(2 \cdot L^2)$ and the coefficient $a$ is set to $a = 2.5 \cdot 10^{-4} \cdot l$.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** The case $\beta$ constant, $\mu = 5.0 \cdot 10^{-3}$, and $\kappa = 1.5$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** The case $\beta$ non constant, $\mu = 5.0 \cdot 10^{-3}$, and $\kappa = 1.5$.</figcaption></figure>

At time $t = 0$ the height in the first part equals $0.3m$, and the other part is empty. The aim of the simulation is to determine the evolution of the height and flow when the dam is removed.

All the simulations have been done with a time step $\tau = 10^{-4}s$ and a space step $\delta = 10^{-4}m$. The friction coefficient equals $\kappa = 1.5m^2s^{-2}$ for the first three simulations,
$\kappa = 0.6m^2s^{-2}$ two next situations and finally $\kappa = 0$ for the last one. The following figures show the flow evolution in the canal for various situations.

We can observe that the wave generated by the dam break is twice reflected on the downstream wall at $x = L$. We can also notice the influence of the function $\beta$: the case when $\beta$ is constant gives rise to much more diffusion on the hight $h$, see fig. 2 and 3. Therefore the non linearity seems to provide a more precise description of the hight. When $\beta$ is not constant, the influence of the viscosity has a strong influence on the fluid hight: see figures 3 and 4. Figure 6 shows the situation of two dam brakes, with a non constant function $\beta$, $\mu = 5.010^{-4}$, $\kappa = 0.6$. In this case we observe the formation of two waves moving in opposite directions. Figure 7 shows the case when the friction coefficient equals zero. In this situation, we can observe some instabilities for a too low viscosity coefficient. Therefore a larger one is used.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** The case $\beta$ non constant, $\mu = 1.010^{-3}$, and $\kappa = 1.5$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** The case $\beta$ non constant, $\mu = 5.010^{-3}$, and $\kappa = 0.6$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** *The case $\beta$ non constant, $\mu = 5.0 \times 10^{-4}$, and $\kappa = 0.6$.*</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** *The case $\beta$ non constant, $\mu = 5.0 \times 10^{-2}$, and $\kappa = 0.0$.*</figcaption></figure>

**5. References**

[1] J-P BENQUÉ, A. HAUGEL, P-L VIOLET. "Engineering Applications of Computational hydraulics, Homage to Alexandre Preissmann, vol. II, 1992. "Nonstandard Asymptotic Analysis", Springer Lecture Notes in Math., vol. 1249, 1987.

[2] A. BERMÜDEZ. "Mathematical techniques for some environnemental problems related to water pollution control", in *Environment, Economics and their Mathematical Models*, edited by J. I. Diaz J. L. Lions, Masson 1994.

[3] A. BERMÜDEZ. Mathematical techniques for some environnemental problems related to water pollution control, in Mathematics Climate and environment, edited by J. I. Diaz J. L. Lions, Masson 1993.

[4] D. BRESCH, B. DESJARDINS, C.K. LIN. On some compressible fluid models: Korteweg, lubrication and shallow water systems. To appear in *Comm. Partial Differential Equations* 28, (2003), no. 3-4, 843-868.

[5] D. BRESCH, B. DESJARDINS. Sur un modèle de Saint-Venant visqueux et sa limite quasi-géostrophique. *C.R.Acad.Paris, Ser.I* 335 (2002), 973-978.

[6] H. BREZIS. Analyse Fonctionnelle, Théorie et Applications, Masson, 1987.

[7] A.J.C. DE SAINT-VENANT Théorie du mouvement non permanent des eaux, avec application aux crues de rivières et à l'introduction des marrées dans leur lit, *C. R. Acad. Sc. Paris*, 73, (1871), 147-154.

[8] L. C. EVANS. *Partial Differential Equation*, Publ. A.M.S., 1998.

[9] F. GERBEAU, B. PERTHAME. Derivation of viscous Saint-Venant system for laminar shallow water; Numerical results. *Discrete and Continuous Dynamical Systems-series B. 1*, No. 1, (2001), 89-102.
[10] S. HUANG, J. ZHAO. Structure of Flood Wave with Viscosity. *Math. Meth. Appl. Sci.*, **23**, 949-963 (2001).

[11] P. M. JACOVKIS, E. G. TABAK. A kinematic Wave Model for Rivers with Flood Plains and Other Irregular Geometries *Mathl. Comput. Modelling*, Vol. 24, No. 11, pp 1-21, 1996.

S. KANE Analyse mathématique et simulation numérique de modèles d'écoulement de fluide incompressible en surface libre et milieu poreux déformable *Thèse Univ. Neuchâtel* (2005)

[12] R.J. LEVEQUE, Numerical Methods for Conservation Laws, *Birkhauser*, 1991.

[13] C.D. LEVERMORE, M. SAMMARTINO. A shallow water model with eddy viscosity for basins with varying bottom topography. *Nonlinearity*, **14**, (2001), No. 6, 1493-1515.

[14] P.-L. LIONS. Mathematical topics in fluid dynamics, Vol 1 Incompressible models, *Oxford University Press*, 1998.

[15] P.-L. LIONS. Mathematical topics in fluid dynamics, Vol 2 Compressible models, *Oxford University Press*, 1998.

[16] P. ORENGA. Un théorème d'existence de solutions d'un problème de shallow water. *Arch. Rat. Mech. Anal.*, **130**, No. 2, (1995), 183-204.

[17] A. PAQUIER. Modélisation et simulation de la propagation de l'onde de rupture de barrage, *Thèse Univ. J. Monnet*, 1995.

[18] O. PIRONNEAU. Méthode des éléments finis pour les fluides, *Masson*, 1988.

[19] J. SIMON. Compact Sets in the Space $L^p(0, T; B)$ *Ann. mat. pura e applicata*, 1987, No IV, Vol. CXLVI, pp. 65-96.

[20] L. SUNDBYE. Global existence for the Dirichlet problem for the viscous shallow water equations. *J. Math. Anal. Appl.*, **202**, No. 1, (1996), 236-258.

[21] R. TEMAM. Navier-Stokes Equations, *North-Holland*, 1991.

[22] R. KH. ZEYTOUNIAN. Modélisation asymptotique en mécanique des fluides newtoniens, *Springer-Verlag*, 1994.