# Characteristics Time Discretization and Mixed Finite Element Approximation for Shallow Water Equations: A priori Error Estimates *

Fadi El Dabaghi, Abdellah El Kacimi, Bassam Nakhle

► To cite this version:

Fadi El Dabaghi, Abdellah El Kacimi, Bassam Nakhle. Characteristics Time Discretization and Mixed Finite Element Approximation for Shallow Water Equations: A priori Error Estimates *. 2025. hal-05009235

HAL Id: hal-05009235

https://inria.hal.science/hal-05009235v1

Preprint submitted on 27 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Characteristics Time Discretization and Mixed Finite Element Approximation for Shallow Water Equations: A priori Error Estimates*

F. El Dabaghi<fnref n="1" />, A. El Kacimi<fnref n="2" />, and B. Nakhle<fnref n="3" />

preprint under submission 2023-2024

This paper is devoted to the error analysis of a numerical velocity-depth formulation for shallow water equations based on the characteristics method for the time discretization combined to mixed finite element for spatial approximation. It is proved by induction that the proposed numerical mixed formulation has a unique convergent solution, under some regularity assumptions on the exact solution. A priori error estimates are established, via $L^2$ projection. A numerical validation of the method is illustrated on a rectangular domain flow.

**KEY WORDS:** shallow water equations, characteristics, mixed finite element, a priori error estimates

## 1. Introduction

Significant effort has been focused, these last years, on the development of systems, for understanding and modeling water related environmental problems. These tools handle computational aspects through mathematical models, describing various physical phenomena, as for example watershed, flood in rivers, tidal fluctuation, breaking waves on shallow beaches,

<footnote>*Supported by WADI (EC) and CMIFM A.I. n° MA/01/03</footnote>

<footnote n="1">¹INRIA - UR Paris, 2 rue Simone Iff, CS 42112 75589 Paris Cedex 12, France - Contact fdabaghi@yahoo.fr</footnote>

<footnote n="2">²EMI-Ecole Mohamadia d'Ingenieurs, Rabat, Maroc.</footnote>

<footnote n="3">³Universit Libanaise, Beyrouth, Liban.</footnote>
bay and estuary flows etc. Particularly, the planning and mangement of
specific watersheds, that are charaterized with the water scarcity and lack
of groundwater, with the objective of improving the water supply in terms
of quantity (optimal dam, water reuse, aquifers etc) and quality (sediments,
pollution, etc) remains a great challenge in the next millennium.

Various numerical methods serving these different purposes, have been
implemented with the finite element method. This allows using high de-
gree of grid flexibility in the complex fluids domains, arising from irregular
land boundaries or other coastal geometries. We refer, for example, to the
ADCIRC system developed by Luettich et al. [19], or the TELEMAC 2D
system developed by Galland et al. [16], which are based on the finite
element approach.

The shallow water equations (SWE) can be seen as advection-diffusion
problem, taking into account diffusive effects, such as eddy viscosity. How-
ever, convection effects dominates diffusion effects, as for example in chan-
nels and narrow inlets. Moreover, special attention should be paid when
dealing with the continuity equation, which is of hyperbolic type, and does
not contain dissipative terms.

The characteristics finite element method is a natural procedure to han-
dle advection flow. It consists in combining a finite element approach with
a time discretization of Lagrangian derivative along the characteristics.
This provides an upwind scheme which has the advantage of stabilizing
convection, allowing therefore the use of large and accurate time-steps.

For convection-diffusion problems, this approach is formulated and an-
alyzed by Bercovier et al. [4], Ewing et al. [15], Russell [22], and Douglas
et al. [10]. For the incompressible Navier-Stokes problem, similar schemes
have been considered by Pironneau [20], Süli [23], El Dabaghi [14] and Ab-
delwahed et al. [1]. In the framework of SWE problems, this approach
combined with the fractional step method, is adopted by Goutal [18], for
the numerical simulation of transcritical flows, and implemented afterwards
in the TELEMAC system. A similar approach, based on the characteristics
finite element formulation with Chorin-type projection and the fractional
step method, has been considered by Zienkiewc et al. [24]. In [9], Dawson et
al. established, a priori error estimates for the characteristics finite element
formulation. In this work, following Goutal [18], we define the character-
istics mixed finite element formulation, for shallow water equations, and
derive a priori error estimates. Some preliminary numerical results, based on
this approach, are obtained in El Dabaghi et al. ([12],[13]). Our analysis
is inspired from the authors mentioned above. Nevertheless, some difficul-
ties mainly related to a supplementary nonlinearity of the advection term,
appearing in the continuity equation, require a special treatment.

The remainder of this paper is outlined as follows. In section 2, we
introduce some notations and definitions. In section 3, we review the gov-
erning equations for 2D shallow water model. Section 4 describes the temporal discretization by the characteristics method. Section 5 is devoted to the variational formulation of the obtained scheme and some existence and uniqueness results. In section 6, we define the characteristics mixed finite element formulation and give an existence and uniqueness result. In section 7, we set some regularity assumptions and boundedness induction hypotheses, and derive using $L^2$ projection a priori error estimates. Finally in section 8 we present some numerical test cases validating the approach methodology as well as the error estimation results.

## 2. Preliminaries

In this section we recall some standard notations and definitions of the classical Sobolev spaces, that will be used in the sequel. Let $\Omega$ be a bounded polygonal domain of $\mathbb{R}^2$, whose generic point is denoted $\mathbf{x} = (x, y)$. The euclidean inner product in $\mathbb{R}^2$ is denoted by $\mathbf{x} \cdot \mathbf{y}$, with $\mathbf{x}, \mathbf{y} \in \mathbb{R}^2$, and the associated norm by $|\mathbf{x}| = \sqrt{\mathbf{x} \cdot \mathbf{x}}$. The dyadic product '$\otimes$' is defined by

$$ (\mathbf{x} \otimes \mathbf{y})_{i,j} = x_i y_j. $$

The scalar product (or double dot product) of two tensors is defined as

$$ A \cdot B = \sum_{i,j} A_{i,j} B_{i,j}. $$

The $L^2$ inner product is denoted by

$$ (\mathbf{u}, \mathbf{v}) = \int_{\Omega} \mathbf{u} \cdot \mathbf{v} \, d\mathbf{x} $$

where the notation '$\cdot$' refers to either multiplication, dot product or double dot product.

For any nonnegative integer $k$, the sobolev space

$$ H^k(\Omega) = \{w \in L^2(\Omega); D^\alpha w \in L^2(\Omega), \text{ for } |\alpha| \le k\}, $$

is endowed with the usual norm and semi-norm :

$$ \|w\|_k = \left( \sum_{|\alpha| \le k} \int_{\Omega} |D^\alpha w|^2 d\mathbf{x} \right)^{\frac{1}{2}}, \quad |w|_k = \left( \sum_{|\alpha|=k} \int_{\Omega} |D^\alpha w|^2 d\mathbf{x} \right)^{\frac{1}{2}} $$

where $\alpha = (\alpha_1, \alpha_2)$ is a multi-index in $\mathbb{N}^2$, $|\alpha| = |\alpha_1| + |\alpha_2|$, and $D^\alpha$ denote a partial derivative : $D^\alpha = \frac{\partial^{|\alpha|}}{\partial x^{\alpha_1} \partial y^{\alpha_2}}$. The convention $H^0(\Omega) = L^2(\Omega)$ is adopted. Moreover, we shall use the Sobolev space

$$ W^{1,\infty}(\Omega) = \{w \in L^\infty(\Omega); \nabla w \in L^\infty(\Omega)\}, $$
endowed with the norm and semi-norm:

$$ \|w\|_{1,\infty} = \max(\|w\|_{\infty}, \|\nabla w\|_{\infty}), \quad |w|_{1,\infty} = \|\nabla w\|_{\infty}. $$

We shall introduce the space $C^{0,1}(\Omega)$ of Lipshitz continuous functions on $\bar{\Omega}$. Let $\Delta t$ be the time step and $t^n = n\Delta t$, with $n = 0, 1, \dots, N$, ($N = \frac{T}{\Delta t}$). For a sequence of functions $\{w^n\}$ defined at the time levels $n = 1, \dots, N$, we use the notations

$$ \|w\|_{\infty} = \max_{n=1,N} \|w^n\|_{\infty}, \quad \|w\|_{0} = \left( \sum_{n=1,N} \Delta t \|w^n\|_{0}^{2} \right)^{\frac{1}{2}}, $$

$$ \|w\|_{k,\infty} = \max_{n=1,N} |w^n|_k, \quad \|w\|_{k,0} = \left( \sum_{n=1,N} \Delta t |w^n|_k^2 \right)^{\frac{1}{2}}. $$

If $w(x, t)$ is defined on the entire time interval $[0, T]$, we use the same notations, as for example,

$$ \|w\|_{k,\infty} = \max_{t \in [0,T]} |w(\cdot, t)|_k, \quad \|w\|_{k,0} = \left( \int_0^T |w(\cdot, t)|_k^2 dt \right)^{\frac{1}{2}}. $$

We often denote multi-valued functions in boldface type. We adopt analogous notations for norms and inner products, for multi-valued functions or scalar functions, for example,

$$ (\mathbf{u}, \mathbf{v}) = (u_1, v_1) + (u_2, v_2), \quad \|\mathbf{u}\|_0 = (\mathbf{u}, \mathbf{u})^{\frac{1}{2}}. $$

### 3. Problem Setting

The two-dimensional SWE are obtained by vertical averaging of the con- tinuum mass, and momentum balances, given by the three dimensional incompressible Navier-Stokes equations. It provides a model describing flow in fluid domains whose bathymetric depth is much smaller than the characteristic length scale in the horizontal direction. Neglecting the influ- ence of the vertical component of the flow acceleration (hydrostatic pressure assumption), the Saint-Venant SWE, written in the conservative form [3], become

$$ \frac{\partial H}{\partial t} + \nabla \cdot (H \mathbf{u}) = 0, \qquad (3.1) $$

$$ \frac{\partial(H\mathbf{u})}{\partial t} + \nabla \cdot (H\mathbf{u} \otimes \mathbf{u}) - \nabla \cdot (\mu H \nabla \mathbf{u}) + gH \nabla z = H\mathbf{F}, \qquad (3.2) $$
with

$$ \mathbf{F} = -\frac{g}{C^2 H} |\mathbf{u}| \mathbf{u} - \frac{\nabla p_a}{\rho} + \frac{\bar{\tau}}{\rho H} + \mathbf{l} \times \mathbf{u}, $$

where $\mathbf{u} = (u, v)^{\perp}$ is the depth-averaged horizontal velocities, $z$ is the free surface elevation under a reference level, $H = z - z_b$ is the water depth, $z_b$ being the elevation of the bottom under the same reference (see figure 1), $\mu$ is the viscosity coefficient assumed constant, $g$ is the gravity acceleration, $C$ is the Chezy coefficient, $\bar{\tau} = (\bar{\tau}_x, \bar{\tau}_y)^{\perp}$ is the wind stress, $\mathbf{l} \times \mathbf{u} = (-lv, lu)^{\perp}$, $l$ is the Coriolis parameter, $p_a$ is the atmospheric pressure, and $\rho$ is the density assumed constant.

<caption>Figure 1: Water elevation and bathymetry</caption>

This form is one of the most closely related with the derivation of shallow water equations from the 3D Navier-Stokes problem. Taking into account the continuity equation (3.1), we obtain by splitting the advection terms in (3.2), another popular form for the equation of motion:

$$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} - \frac{\mu}{H} \nabla \cdot (H \nabla \mathbf{u}) + g \nabla z = \mathbf{F}. \quad (3.3) $$

If we suppose that $|\frac{\nabla H}{H} \cdot (\mu \nabla \mathbf{u})| \ll |\mathbf{u} \cdot \nabla \mathbf{u}|$, see ([2], [21]), the divergence operator can be written as follows

$$ \frac{\mu}{H} \nabla \cdot (H \nabla \mathbf{u}) \approx \mu \Delta \mathbf{u}. $$

This yields the non conservative form of SWE:

$$ \frac{\partial H}{\partial t} + \nabla \cdot (H \mathbf{u}) = 0, \quad (3.4) $$

$$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} - \mu \Delta \mathbf{u} + g \nabla z = \mathbf{F}. $$

This system of equations is completed with the following initial conditions

$$ H(., 0) = H_0, \quad \mathbf{u}(., 0) = \mathbf{u}_0. $$

## 4. Temporal discretization

We start by writing the system of equations (3.4), as follows

$$ \frac{dH}{dt} + H \nabla \cdot \mathbf{u} = 0, \quad \frac{du}{dt} - \mu \Delta \mathbf{u} + g \nabla z = \mathbf{F}, $$
where $\frac{d}{dt}$ is the Lagrangian derivative defined for $H$ and $\mathbf{u}$ by

$$ \frac{dH}{dt} = \frac{\partial H}{\partial t} + \mathbf{u} \cdot \nabla H, \quad \frac{d\mathbf{u}}{dt} = \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u}. $$

The hyperbolic character of the above equations allows us the use of the characteristics method, to perform the time discretization of the advection terms. For a given $H^n$, $\mathbf{u}^n$, the discretized equations are

$$ \frac{H^{n+1} - H^n \circ \mathbf{X}^n}{\Delta t} + H^n \nabla \cdot \mathbf{u}^{n+1} = 0, \quad (4.1) $$

$$ \frac{\mathbf{u}^{n+1} - \mathbf{u}^n \circ \mathbf{X}^n}{\Delta t} - \mu \Delta \mathbf{u}^{n+1} + g \nabla z^{n+1} = \mathbf{F}^n, \quad (4.2) $$

where $\mathbf{F}^n$ is a semi-implicit discretization of the source term $\mathbf{F}$, specified in the sequel, and $\mathbf{X}^n = \mathbf{X}^n(\mathbf{x}, t^{n+1}; t^n)$, $\mathbf{X}^n(\mathbf{x}, t^{n+1}; t)$ being the trajectories (characteristics) associated with the transport field, defined by

$$ \frac{d\mathbf{X}^n}{dt}(\mathbf{x}, t^{n+1}; t) = \mathbf{u}^n(\mathbf{X}^n(\mathbf{x}, t^{n+1}; t)), \quad \text{for } t^n \le t \le t^{n+1}, $$

$$ \mathbf{X}^n(\mathbf{x}, t^{n+1}; t^{n+1}) = \mathbf{x}, $$

The solution of this backward differential equation is calculated by back-tracing along the trajectory $\mathbf{X}^n(\mathbf{x}, t^{n+1}; t)$, until $t = t^n$, to find the 'foot' of the characteristics.

## 5. Variational formulation

For the sake of simplicity, we consider boundary Dirichlet conditions for the velocity $\mathbf{u}$. Let us set

$$ V = (H_0^1(\Omega))^d, \quad Q = L^2(\Omega), \quad d \in \{1, 2\}, $$

and consider $\mathbf{v}$ in $V$, and $q$ in $Q$. Multiplying each equation in (4.2) by the corresponding component of $\mathbf{v}$ and (4.1) by $q$ and integrating by parts on $\Omega$, we get

$$ \left( \frac{H^{n+1} - H^n \circ \mathbf{X}^n}{\Delta t}, q \right) + (H^n \nabla \cdot \mathbf{u}^{n+1}, q) = 0 $$

$$ \left( \frac{\mathbf{u}^{n+1} - \mathbf{u}^n \circ \mathbf{X}^n}{\Delta t}, \mathbf{v} \right) + \mu (\nabla \mathbf{u}^{n+1}, \nabla \mathbf{v}) - g (H^{n+1}, \nabla \cdot \mathbf{v}) = (\mathbf{F}^n, \mathbf{v}) + g(z_b, \nabla \cdot \mathbf{v}). $$
For simplicity of notation, we omit in this section, frequently the superscript n. Let us consider an element $\bar{u}$ of $H^1(\Omega)$ such that $u = \bar{u}$, a.e., on $\partial\Omega$ and put $\bar{V} = \bar{u}+V$. If we treat the bottom friction term semi-implicitly, with respect to $u$, i.e.,

$$-\frac{g}{C^2H} |\mathbf{u}|_{t=t^{n+1}} = -\frac{\lambda}{H} |\mathbf{u}|_{t=t^{n+1}} \approx -\frac{\lambda}{H^n} |\mathbf{u}^n| \mathbf{u}^{n+1},$$

the corresponding linear system is given by

$$ (P_1) \quad \left\{ \begin{array}{ll} \text{find } (\mathbf{u}, H) \text{ in } \bar{V} \times Q & \\ a_1(\mathbf{u}, \mathbf{v}) = \langle \mathbf{f}_1, \mathbf{v} \rangle_{V^*, V}, & \text{ for all } \mathbf{v} \text{ in } V, \\ (\alpha H, q) + (\underline{H} \nabla \cdot \mathbf{u}, q) = (c, q), & \text{ for all } q \text{ in } Q, \end{array} \right. $$

where

$$a_1(\mathbf{u}, \mathbf{v}) = (\alpha \mathbf{u}, \mathbf{v}) + \mu (\nabla \mathbf{u}, \nabla \mathbf{v}) + \frac{g}{\alpha} (\underline{H} \nabla \cdot \mathbf{u}, \nabla \cdot \mathbf{v}) + \lambda \left( \frac{|\mathbf{u}^n| \mathbf{u}}{\underline{H}}, \mathbf{v} \right),$$

$$ \langle \mathbf{f}_1, \mathbf{v} \rangle_{V^*, V} = \left\langle -\frac{\nabla p_a}{\rho} + \frac{\bar{\tau}}{\rho H^n} + \mathbf{l} \times \mathbf{u}^n + \alpha \mathbf{u}^n \circ \mathbf{X}^n, \mathbf{v} \right\rangle_{V^*, V} \\ + g(z_b, \nabla \cdot \mathbf{v}) + \frac{g}{\alpha} (c, \nabla \cdot \mathbf{v}). $$

The following assumption ensures the existence and the uniqueness of the solution for $(P_1)$:

$$ \underline{H} \in L^\infty(\Omega), \underline{H} \ge \underline{H}_0 > 0, \text{ a.e. in } \Omega, \\ \mathbf{u}^n \in (L^2(\Omega))^2, \mathbf{f}_1 \in (H^{-1}(\Omega))^2, H^n \circ \mathbf{X}^n \in L^2(\Omega). \quad (A_1) $$

**Proposition 1** Under assumption $(A_1)$, Problem $(P_1)$ admits a unique solution $(\mathbf{u}, H)$.

**Proof.** For a proof of this result, we refer the reader to El Dabaghi et al. [11]. ■

**Remark 2** Let us notice that the condition $\mathbf{u}^n \in (L^2(\Omega))^2$ guarantees that the symmetric bilinear form $(\mathbf{v}, \mathbf{w}) \to \lambda \left( \frac{|\mathbf{u}^n| \mathbf{w}}{\underline{H}}, \mathbf{v} \right)$ is continuous from $V^2$ on $\mathbb{R}$. Indeed, we check for all $\mathbf{w}_1, \mathbf{w}_2$ and $\mathbf{v}$ in $V$, using the fact that the injection from $L^4(\Omega)$ onto $H_0^1(\Omega)$ is continuous (for $d \le 2$) and the generalized h\"older inequality, that
$$ \lambda \left| \left( \frac{|\mathbf{u}^n| (\mathbf{w}_1 - \mathbf{w}_2)}{\underline{H}} , \mathbf{v} \right) \right| \le \lambda \frac{\|\mathbf{u}^n\|_0}{\underline{H}_0} \|\mathbf{v}\|_{4,\Omega} \|\mathbf{w}_1 - \mathbf{w}_2\|_{4,\Omega} \\ \le C(\lambda, d, \Omega) \frac{\|\mathbf{u}^n\|_0}{\underline{H}_0} \|\mathbf{v}\|_1 \|\mathbf{w}_1 - \mathbf{w}_2\|_1, $$

where $\|\cdot\|_{4,\Omega}$ denotes the norm in $L^4(\Omega)$.

Let us multiply the second equation of $(\mathcal{P}_1)$ by $-\frac{gq}{\underline{H}}$, set $\alpha_g = g\alpha$, $D = -\frac{gc}{\underline{H}}$ and define the following functionals

$$ a_2(\mathbf{u}, \mathbf{v}) = (\alpha \mathbf{u}, \mathbf{v}) + \mu (\nabla \mathbf{u}, \nabla \mathbf{v}) + \lambda \left( \frac{|\mathbf{u}^n| \mathbf{u}}{\underline{H}}, \mathbf{v} \right), $$

$$ b(\mathbf{v}, H) = -g(\nabla \cdot \mathbf{v}, H), \quad e(H, q) = \left(\frac{H}{\underline{H}}, q\right), $$

$$ \langle \mathbf{f}_2, \mathbf{v} \rangle_{V^*, V} = \left\langle -\frac{\nabla p_a}{\rho} + \frac{\bar{\tau}}{\rho H^n} + \mathbf{1} \times \mathbf{u}^n + \alpha \mathbf{u}^n \circ \mathbf{X}^n, \mathbf{v} \right\rangle_{V^*, V} \\ + g(z_b, \nabla \cdot \mathbf{v}). $$

Then Problem $(\mathcal{P}_1)$ is equivalent to

$$ (\mathcal{P}) \quad \left\{ \begin{array}{ll} \text{find } (\mathbf{u}, H) \text{ in } \bar{V} \times Q & \\ a_2(\mathbf{u}, \mathbf{v}) + b(\mathbf{v}, H) = \langle \mathbf{f}_2, \mathbf{v} \rangle_{V^*, V} & , \text{ for all } \mathbf{v} \text{ in } V, \\ -\alpha_g e(H, q) + b(\mathbf{u}, q) = (D, q), & , \text{ for all } q \text{ in } Q. \end{array} \right. $$

The above formulation, who was initiated by Goutal in her thesis [18], is similar to the penalized quasi-Stokes problem ([6], [17]), where the penalization parameter is equal to $\alpha_g$. Since the bilinear forms $e$ and $a_2$ are symmetric, Problem $(\mathcal{P})$ corresponds to the saddle point problem:

$$ \inf_{\bar{\mathbf{v}} \in \bar{V}} \sup_{q \in Q} \frac{1}{2} a_2(\bar{\mathbf{v}}, \bar{\mathbf{v}}) + b(\bar{\mathbf{v}}, q) - \frac{1}{2} \alpha_g e(q, q) - \langle \mathbf{f}_2, \bar{\mathbf{v}} \rangle - (D, q). $$

# 6. Mixed finite element formulation

Let $T_h$ be a regular triangulation of the domain $\Omega$ into elements $T$, with $h_T = \text{diam}(T)$ and $h = \max_{T \in T_h} h_T$. Let $W_{l,h}$ be a finite dimensional subspace of $H^1(\Omega)$ defined on this triangulation, consisting of piecewise polynomials in $\mathbb{R}^d$ of degree less than or equal to $l$. Let us set $V_h = W_{k,h}^d \cap V$, $\bar{V}_h =$
$W_{k,h}^d \cap \bar{V}$ and $Q_h = W_{m,h} \cap Q$, with $k, m \in \mathbb{N}^*$. From now, we fix our attention to the following approximate mixed finite element problems for $(\mathcal{P})$:

$$ (\mathcal{P}_h) \left\{ \begin{array}{l} \text{find } (\mathbf{u}_h^{n+1}, H_h^{n+1}) \text{ in } \bar{V}_h \times Q_h \\ a_{2,h}^n(\mathbf{u}_h^{n+1}, \mathbf{v}) + b(\mathbf{v}, H_h^{n+1}) = \langle \mathbf{f}_{2,h}^n, \mathbf{v} \rangle_{V^*,V}, \quad (6.1) \\ -\alpha_g e_h^n(H_h^{n+1}, q) + b(\mathbf{u}_h^{n+1}, q) = (D_h^n, q), \quad . \end{array} \right. $$

for all $(\mathbf{v}, q)$ in $V_h \times Q_h$, where $\mathbf{X}_h^n = \mathbf{X}_h^n(\mathbf{x}, t^{n+1}; t^n)$, $\mathbf{X}_h^n(\mathbf{x}, t^{n+1}; t)$ being the discretized characteristics defined by

$$ \frac{d\mathbf{X}_h^n}{dt}(\mathbf{x}, t^{n+1}; t) = \mathbf{u}_h^n(\mathbf{X}_h^n(\mathbf{x}, t^{n+1}; t)), \quad \text{for } t^n \le t \le t^{n+1}, $$

$$ \mathbf{X}_h^n(\mathbf{x}, t^{n+1}; t^{n+1}) = \mathbf{x}, $$

$$ D_h^n = -\frac{gH_h^n \circ \mathbf{X}_h^n}{\Delta t H_h^n}, \quad e_h^n(H, q) = (\frac{H}{H_h^n}, q). $$

and

$$ a_{2,h}^n(\mathbf{u}, \mathbf{v}) = (\alpha \mathbf{u}, \mathbf{v}) + \mu (\nabla \mathbf{u}, \nabla \mathbf{v}) + \lambda \left( \frac{|\mathbf{u}_h^n|}{H_h^n}, \mathbf{v} \right). $$

We now give an existence and uniqueness results for $(\mathcal{P}_h)$. To this end, we make the following induction assumption :

$$ H_h^l \in L^\infty(\Omega), \quad H_h^l \ge \underline{H}_0 > 0, \text{ a.e. in } \Omega, \text{ for all } l = 0, \dots, n. \quad (A_1) $$

**Proposition 3** Under assumption $(A_1)$, Problem $(\mathcal{P}_h)$ admits a unique solution $(\mathbf{u}_h^{n+1}, H_h^{n+1})$.

**Proof.** For a proof of the above proposition see [18]. ■

# 7. Error Analysis

The present section is devoted to a priori error estimates for the numerical solution given by the characteristics mixed finite element method proposed in §6. The well-known inequality

$$ ab \le \varepsilon a^2 + \frac{1}{4\varepsilon} b^2, $$

where $a, b \in \mathbb{R}$, $\varepsilon \in \mathbb{R}^*$, is frequently used. In all that follows, we suppose that the exact solution $(\mathbf{u}, H)$ exists, is unique and satisfies regularity assumptions:

$$ (\mathbf{u}, H) \in C([0, T], (C^{0,1}(\Omega))^d) \times (L^\infty([0, T[\times\Omega))^d, $$
$\underline{H}_0 \le H \le \underline{H}_1$, a.e., $\Omega$, for all $t \in [0, T]$,

$$ \|\mathbf{u}\|_{k+1,0}, \|H\|_{m+1,0}, \|\mathbf{u}\|_{\infty}, \|\nabla \mathbf{u}\|_{\infty}, \|\nabla H\|_{\infty} \le M, $$

$$ \left\| \frac{\partial H}{\partial t} \right\|_{\infty}, \left\| \frac{\partial H}{\partial t} \right\|_{m+1,0}, \left\| \frac{\partial \mathbf{u}}{\partial t} \right\|_{0,\infty}, \left\| \frac{d^2 H}{dt^2} \right\|_{0}, \left\| \frac{d^2 \mathbf{u}}{dt^2} \right\|_{0} \le M, $$

where $\frac{d^2}{dt^2}$ is the second Lagrangian derivative.

$$ \|H\|_{m+1,\infty}, \|\mathbf{u}\|_{k+1,\infty}, \left\| \frac{d^2 H}{dt^2} \right\|_{0,\infty}, \left\| \frac{\partial H}{\partial t} \right\|_{m+1,\infty} \le M. $$

We will also make, for $l = 0, \dots, n$, the following boundedness induction hypotheses:

$$ (\mathbf{H}_0) \qquad \frac{\underline{H}_0}{2} \le H_h^l \le \underline{H}_1 + \frac{\underline{H}_0}{2} \quad \text{a.e., } \Omega. $$

$$ (\mathbf{H}_1) \qquad \|\mathbf{u}_h^l\|_{\infty}, \|\nabla H_h^l\|_{\infty}, \|\nabla \mathbf{u}_h^l\|_{\infty} \le M + 1. $$

$$ (\mathbf{H}_2) \qquad \alpha \|H_h^l - H_h^{l-1}\|_{\infty} \le M + 1. $$

We now recall some useful results about the approximation of the Lagrangian derivative, who are largely inspired from Ewing et al. [15], Russell [22], Douglas et al. [10], and Süli [23].

In order to avoid technical boundary difficulties related to the characteristics method we consider homogenous boundary conditions, i.e., $\bar{\mathbf{u}} = 0$. For simplicity of notations, let us introduce :

$$ \mathbf{X}(\mathbf{x}; \tau) = \mathbf{X}(\mathbf{x}, t; \tau), \quad \mathbf{X}^l(\mathbf{x}; \tau) = \mathbf{X}(\mathbf{x}, t^{l+1}; \tau), $$

$$ \mathbf{X}_h^l(\mathbf{x}; \tau) = \mathbf{X}_h^l(\mathbf{x}, t^{l+1}; \tau), \quad J_h^l(\mathbf{x}; \tau) = J_h(\mathbf{x}, t^{l+1}; \tau). $$

In all that follows, $\mathbf{X}(\mathbf{x}, t; \tau)$ are the trajectories defined by

$$ \frac{d\mathbf{X}}{d\tau}(\mathbf{x}, t; \tau) = \mathbf{u}(\mathbf{X}(\mathbf{x}, t; \tau), \tau), $$

$$ \mathbf{X}(\mathbf{x}, t; t) = \mathbf{x}, $$

and $\mathbf{X}^l = \mathbf{X}(\mathbf{x}, t^{l+1}; t^l)$. So, we supposed that $\mathbf{u} \in C([0, T], (C^{0,1}(\Omega))^d)$ to ensure, thanks to Cauchy-Lipshitz theorem, the existence and uniqueness of the characteristics (see Süli [23] for more details).

Proofs of the two following lemmas, can be found in the last mentioned reference.
**Lemma 4** There exists a constant $\delta_0 = \delta_0(M)$ such that for all $\theta$ in $[0, 1]$ and $|t - \tau| \le \delta_0$, the mapping $\mathbf{x} \to (1 - \theta)\mathbf{x} + \theta\mathbf{X}(\mathbf{x}; \tau)$ is a homeomorphism from $\Omega$ onto itself and its Jacobian satisfies

$$ \frac{1}{2} \le J_{\theta}(\mathbf{x}, t; \tau) \le \frac{3}{2}, \text{ a.e. in } \Omega. $$

**Lemma 5** Under Hypothesis (H₁), there exists a constant $\delta_1 = \delta_1(M)$ such that for all $\theta$ in $[0, 1]$ and $\Delta t \le \delta_1$, the mapping $\mathbf{x} \to (1 - \theta)\mathbf{X}_h^l(\mathbf{x}; \tau) + \theta\mathbf{X}^l(\mathbf{x}; \tau)$ is a homeomorphism from $\Omega$ onto itself, and its Jacobian satisfies

$$ J_{h, \theta}^{l}(\mathbf{x}; \tau) \geq \frac{1}{2}, \text{ a.e. in } \Omega. $$

for all $\tau$ in $[t^l, t^{l+1}]$.

**Lemma 6** Under Hypothesis (H₁), there exists constants $\delta_2 = \delta_2(M)$ and $C_0 = C_0(M)$ such that for $\Delta t \le \delta_2$ and all functions $w$ in $L^2(\Omega)$ and $f$ in $H^1(\Omega)$, we have

$$ |\alpha(w, f \circ \mathbf{X}_h^l - f)| \le C_0 \|\nabla f\|_0 \|w\|_0. $$

**Proof.** For a proof of this lemma we refer to Ewing et al. [15], Russell [22], Douglas et al. [10]. ■

**Proposition 7** Under Hypothesis (H₁), there exists constants $\delta_3 = \delta_3(M)$ and $C_1 = C_1(M)$ such that for $\Delta t \le \delta_3$, the mapping $\mathbf{x} \to \mathbf{X}_h^l(\mathbf{x}; \tau)$ is a homeomorphism from $\Omega$ onto itself, and its Jacobian satisfies

$$ \frac{1}{2} \le J_h^l(\mathbf{x}; \tau) \le \frac{3}{2}, \text{ a.e. in } \Omega, $$

$$ \left| 1 - \frac{1}{J_h^l(\mathbf{x}; \tau)} \right| \le C_1 \Delta t, \text{ a.e. in } \Omega, $$

for all $\tau$ in $[t^l, t^{l+1}]$.

**Proof.** The first part of the proof is immediate from Lemma 5, by choosing $\theta = 0$, under the condition $\Delta t \le \delta_1$. The remainder of the proof is based on Liouville's identity (see Chorin [8]). For more details see El Dabaghi et al. [11]. ■

**Proposition 8** Under Hypothesis (H₁), there exists a constant $C_2 = C_2(M)$ such that for $\Delta t \le \delta_3$, we have

$$ \|\mathbf{X}_h^l - \mathbf{X}^l\|_0^2 \le C_2 \Delta t \left( \|\mathbf{u}_h^l - \mathbf{u}(t^l, \cdot)\|_0^2 + \Delta t \int_{t^l}^{t^{l+1}} \left\| \frac{\partial \mathbf{u}}{\partial t} \right\|_0^2 ds \right), $$

where $\delta_3$ is given by Proposition 7.
**Proof.** For a proof of the above inequality, we refer to El Dabaghi et al. [11]. Proposition 7 and Gronwall's formula are used. ■

Let us define another norm in $L^2(\Omega)$ by

$$ \|w\|_{0,l}^2 = e_h^l(w, w) = \left(\frac{w}{H_h^l}, w\right), \text{ for all } l = 0, \dots, N. $$

The following proposition is useful for the sequel.

**Proposition 9** Let us assume that $(\mathbf{H}_1)$ holds. Then for $\Delta t \le \delta_3$, we have

$$ \alpha \left| \|w \circ \mathbf{X}_h^l\|_0^2 - \|w\|_0^2 \right| \le C_1 \|w\|_0^2. $$

If in addition $(\mathbf{H}_0)$ holds, then there exists a constant $C_3 = C_3(M, \underline{H}_0)$ such that

$$ \alpha \left| \|w \circ \mathbf{X}_h^l\|_{0,l}^2 - \|w\|_{0,l}^2 \right| \le C_3 \|w\|_0^2, $$

where $C_1$ and $\delta_3$ are given by Proposition 7.

**Proof.** The proof is based on the same arguments of Ewing et al. [15], Russell [22], Douglas et al. [10]. We now prove the second inequality. Using the change of variables $\mathbf{y} = \mathbf{X}_h^l(\mathbf{x}; t^{l+1})$ and Proposition 7, we show that

$$ \alpha \|w \circ \mathbf{X}_h^l\|_{0,l}^2 - \|w\|_{0,l}^2 = \alpha \int_{\Omega} \frac{w^2(\mathbf{y})}{H_h^l(\mathbf{x})} \left( \frac{1}{J_h^l(\mathbf{x}; t^{l+1})} - 1 \right) d\mathbf{y} \\ + \alpha \int_{\Omega} \left( \frac{1}{H_h^l(\mathbf{x})} - \frac{1}{H_h^l(\mathbf{y})} \right) w^2(\mathbf{y}) d\mathbf{y} $$

Observing that

$$ \alpha |H_h^l(\mathbf{x}) - H_h^l(\mathbf{y})| \le \alpha \|H_h^l \circ \mathbf{X}_h^l - H_h^l\|_{\infty} \\ \le \alpha \|\mathbf{X}_h^l - \mathbf{x}\|_{\infty} \|\nabla H_h^l\|_{\infty} \le \|\mathbf{u}_h^l\|_{\infty} \|\nabla H_h^l\|_{\infty}, $$

the induction hypotheses $(\mathbf{H}_0) - (\mathbf{H}_1)$ and Proposition 7, imply

$$ \alpha \left| \|w \circ \mathbf{X}_h^l\|_{0,l}^2 - \|w\|_{0,l}^2 \right| \le \left( \frac{2C_1}{\underline{H}_0} + \frac{4(1+M)^2}{\underline{H}_0^2} \right) \|w\|_0^2. $$

So we can choose $C_3 = \left( \frac{2C_1}{\underline{H}_0} + \frac{4(1+M)^2}{\underline{H}_0^2} \right)$. Similarly, we derive the first inequality. But, naturally we do not need to use the fact that $\|\nabla H_h^l\|_{\infty}$ is bounded. ■
In order to derive an error estimate for the numerical solution $(\mathbf{u}_h^l, H_h^l)$, we introduce some notations: Let us denote by $\pi_h q$ the $L^2$ projection of $q$ in $Q_h$ and $\pi_h \mathbf{v}$ the $L^2$ projection of $\mathbf{v}$ in $V_h$ for all $\mathbf{v}$ in $V \cap H^{k+1}(\Omega)$ and $q$ in $Q \cap H^{m+1}(\Omega)$. We recall the following standard approximation properties ([5], [17]):

$$
\begin{align}
\|\pi_h \mathbf{v} - \mathbf{v}\|_0 &\le c_0 h^{k+1} |\mathbf{v}|_{k+1} \quad , \quad |\pi_h \mathbf{v} - \mathbf{v}|_1 \le c_0 h^k |\mathbf{v}|_{k+1}, \nonumber \\
\|\pi_h q - q\|_0 &\le c_0 h^{m+1} |q|_{m+1} \quad , \quad |\pi_h q - q|_1 \le c_0 h^m |q|_{m+1}, \nonumber \\
\|\pi_h \mathbf{v} - \mathbf{v}\|_{\infty} &\le c_0 h^{k-\frac{d}{2}+1} |\mathbf{v}|_{k+1} \quad , \quad |\pi_h \mathbf{v} - \mathbf{v}|_{1,\infty} \le c_0 h^{k-\frac{d}{2}} |\mathbf{u}|_{k+1}, \nonumber \\
\|\pi_h q - q\|_{\infty} &\le c_0 h^{m-\frac{d}{2}+1} |q|_{m+1} \quad , \quad |\pi_h q - q|_{1,\infty} \le c_0 h^{m-\frac{d}{2}} |q|_{m+1}. \tag{7.1}
\end{align}
$$

Moreover, we will use the following inverse estimates

$$
\begin{align}
\|\mathbf{v}\|_{\infty} &\le c h^{-\frac{d}{2}} \|\mathbf{v}\|_{0} \quad , \quad \|q\|_{\infty} \le c h^{-\frac{d}{2}} \|q\|_{0}, \nonumber \\
|\mathbf{v}|_{1} &\le c h^{-1} \|\mathbf{v}\|_{0} \quad , \quad |q|_{1} \le c h^{-1} \|q\|_{0}, \tag{7.2} \\
|\mathbf{v}|_{1,\infty} &\le c h^{-1-\frac{d}{2}} \|\mathbf{v}\|_{0} \quad , \quad |q|_{1,\infty} \le c h^{-1-\frac{d}{2}} \|q\|_{0}, \nonumber
\end{align}
$$

for all $\mathbf{v}$ in $V_h$ and $q$ in $Q_h$.

Let us set $\mathbf{u}^l = \mathbf{u}(t^l)$, $H^l = H(t^l)$ and introduce the following notations

$$
\begin{align*}
\varepsilon_u^l &= \mathbf{u}_h^l - \mathbf{u}^l \quad , \quad \varepsilon_H^l = H_h^l - H^l, \\
\mathbf{e}_u^l &= \pi_h \mathbf{u}^l - \mathbf{u}^l \quad , \quad e_H^l = \pi_h H^l - H^l, \\
\mathbf{E}_u^l &= \mathbf{u}_h^l - \pi_h \mathbf{u}^l \quad , \quad E_H^l = H_h^l - \pi_h H^l,
\end{align*}
$$

The initial conditions $\mathbf{u}_h^0$ and $H_h^0$ are such that $\mathbf{E}_u^0 = 0$ and $E_H^0 = 0$.

Note that $\mathbf{E}_u^l = \varepsilon_u^l - \mathbf{e}_u^l$, and $E_H^l = \varepsilon_H^l - e_H^l$.

It is easy to see that the exact solution $(\mathbf{u}, H)$ satisfies

$$
\begin{align}
& (\alpha (\mathbf{u}^{l+1} - \mathbf{u}^l \circ \mathbf{X}_h^l), \mathbf{v}) + a(H_h^l, \mathbf{u}_h^l, \mathbf{u}^{l+1}, \mathbf{v}) + \nonumber \\
& b(\mathbf{v}H^{l+1}) = r_1(\mathbf{v}) + \langle \mathbf{f}^{l+1}\mathbf{v} \rangle_{V^*,V}, \nonumber \\
& -g(\alpha(H^{l+1} - H^l \circ \mathbf{X}_h^l), q)_{0,l} + b(\mathbf{u}^{l+1}, q) = r_2(q), \tag{7.3}
\end{align}
$$
for all $(\mathbf{v}, q)$ in $V \times Q$. Moreover, (6.1) implies

$$ (\alpha (\mathbf{u}_h^{l+1} - \mathbf{u}_h^l \circ \mathbf{X}_h^l), \mathbf{v}) + a(H_h^l, \mathbf{u}_h^l, \mathbf{u}_h^{l+1}, \mathbf{v}) + b(\mathbf{v}, H_h^{l+1}) = \langle \mathbf{f}_h^l, \mathbf{v} \rangle_{V^*, V}, \quad (7.4) $$

$$ -g(\alpha(H_h^{l+1} - H_h^l \circ \mathbf{X}_h^l), q)_{0,l} + b(\mathbf{u}_h^{l+1}, q) = 0, $$

for all $(\mathbf{v}, q)$ in $V_h \times Q_h$. Where

$$ a(q, \mathbf{w}, \mathbf{u}, \mathbf{v}) = (\alpha \mathbf{u}, \mathbf{v}) + \mu (\nabla \mathbf{u}, \nabla \mathbf{v}) + \lambda \left( \frac{|\mathbf{w}| \mathbf{u}}{q}, \mathbf{v} \right), $$

$$ r_1(\mathbf{v}) = \alpha (\mathbf{u}^{l+1} - \mathbf{u}^l \circ \mathbf{X}_h^l, \mathbf{v}) - \left( \frac{d\mathbf{u}^{l+1}}{dt}, \mathbf{v} \right) + a(H_h^l, \mathbf{u}_h^l, \mathbf{u}^{l+1}, \mathbf{v}) - a(H^{l+1}, \mathbf{u}^{l+1}, \mathbf{u}^{l+1}, \mathbf{v}), $$

$$ r_2(q) = g \left( \frac{\frac{dH}{dt}}{H^{l+1}} - \frac{\alpha(H^{l+1} - H^l \circ \mathbf{X}_h^l)}{H_h^l}, q \right), $$

$$ \langle \mathbf{f}_h^l, \mathbf{v} \rangle_{V^*, V} = \left\langle -g \nabla z_b - \frac{\nabla p_a}{\rho} + \frac{\bar{\tau}}{\rho H_h^l} + \mathbf{l} \times \mathbf{u}_h^l, \mathbf{v} \right\rangle, $$

$$ \langle \mathbf{f}^{l+1}, \mathbf{v} \rangle_{V^*, V} = \left\langle -g \nabla z_b - \frac{\nabla p_a}{\rho} + \frac{\bar{\tau}}{\rho H^{l+1}} + \mathbf{l} \times \mathbf{u}^{l+1}, \mathbf{E}_u^{l+1} \right\rangle. $$

Subtracting (7.4) from (7.3), with $\mathbf{v} = \mathbf{E}_u^{l+1}$, $q = E_H^{l+1}$, we get

$$ (\alpha(\varepsilon_u^{l+1} - \varepsilon_u^l \circ \mathbf{X}_h^l), \mathbf{E}_u^{l+1}) + b(\mathbf{E}_u^{l+1}, \varepsilon_H^{l+1}) + a(H_h^l, \mathbf{u}_h^l, \varepsilon_u^{l+1}, \mathbf{E}_u^{l+1}) = r_1(\mathbf{E}_u^{l+1}) + r(\mathbf{E}_u^{l+1}) $$

$$ -g(\alpha(\varepsilon_H^{l+1} - \varepsilon_H^l \circ \mathbf{X}_h^l), E_H^{l+1})_{0,l} + b(\varepsilon_u^{l+1}, E_H^{l+1}) = r_2(E_H^{l+1}), $$

where $r(\mathbf{E}_u^{l+1}) = \langle \mathbf{f}_h^l - \mathbf{f}^{l+1}, \mathbf{E}_u^{l+1} \rangle_{V^*, V}$. And then we have

$$ (\alpha(\mathbf{E}_u^{l+1} - \mathbf{E}_u^l \circ \mathbf{X}_h^l), \mathbf{E}_u^{l+1}) + a(H_h^l, \mathbf{u}_h^l, \mathbf{E}_u^{l+1}, \mathbf{E}_u^{l+1}) + g(\alpha(E_H^{l+1} - E_H^l \circ \mathbf{X}_h^l), E_H^{l+1})_{0,l} = \bar{r}_1(\mathbf{E}_u^{l+1}) + \bar{r}_2(E_H^{l+1}) + r(\mathbf{E}_u^{l+1}), \quad (7.5) $$
with

$$
\begin{align*}
\bar{r}_1(E_u^{l+1}) &= r_1(\mathbf{E}_u^{l+1}) - a(H_h^l, \mathbf{u}_h^l, \mathbf{e}_u^{l+1}, \mathbf{E}_u^{l+1}) - b(\mathbf{E}_u^{l+1}, \varepsilon_H^{l+1}) \\
&\quad - (\alpha(\mathbf{e}_u^{l+1} - \mathbf{e}_u^l \circ \mathbf{X}_h^l), \mathbf{E}_u^{l+1}).
\end{align*}
$$

$$
\bar{r}_2(E_H^{l+1}) = -r_2(E_H^{l+1}) + b(\varepsilon_u^{l+1}, E_H^{l+1}) - g\alpha(e_H^{l+1} - e_H^l \circ \mathbf{X}_h^l, E_H^{l+1})_{0,l}.
$$

Let us set

$$
T_0 = a(H_h^l, \mathbf{u}_h^l, \mathbf{E}_u^{l+1}, \mathbf{E}_u^{l+1}) \quad , \quad T_1 = (\mathbf{E}_u^{l+1} - \mathbf{E}_u^l \circ \mathbf{X}_h^l, \mathbf{E}_u^{l+1}),
$$

$$
T_2 = (E_H^{l+1} - E_H^l \circ \mathbf{X}_h^l, E_H^{l+1})_{0,l} \quad , \quad T_3 = \|\mathbf{E}_u^l\|_0^2 - \|\mathbf{E}_u^l \circ \mathbf{X}_h^l\|_0^2,
$$

$$
T_4 = \|E_H^l\|_{0,l-1}^2 - \|E_H^l\|_{0,l}^2 \quad , \quad T_5 = \|E_H^l\|_{0,l-1}^2 - \|E_H^l \circ \mathbf{X}_h^l\|_{0,l}^2.
$$

We can now prove some technical lemmas which are very useful in the
sequel.

**Lemma 10** Under the previous notations, we have

$$
\begin{equation}
\begin{aligned}
& \frac{\alpha}{2} \left( \|\mathbf{E}_u^{l+1}\|_0^2 - \|\mathbf{E}_u^l\|_0^2 \right) + \frac{g\alpha}{2} \left( \|E_H^{l+1}\|_{0,l}^2 - \|E_H^l\|_{0,l-1}^2 \right) + \mu \|\mathbf{E}_u^{l+1}\|_1^2 \\
& \leq -\frac{\alpha}{2} T_3 - \frac{g\alpha}{2} T_4 - \frac{g\alpha}{2} T_5 + \bar{r}_1(\mathbf{E}_u^{l+1}) + \bar{r}_2(E_H^{l+1}) + r(\mathbf{E}_u^{l+1}).
\end{aligned}
\tag{7.6}
\end{equation}
$$

**Proof.** We first establish lower bounds for $T_0$, $T_1$ and $T_2$. We check easily that

$$
\begin{align}
T_0 &= \mu (\nabla \mathbf{E}_u^{l+1}, \nabla \mathbf{E}_u^{l+1}) + \lambda \left( \frac{|\mathbf{u}_h^l|}{H_h^l} \mathbf{E}_u^{l+1}, \mathbf{E}_u^{l+1} \right) \tag{7.7} \\
    &\geq \mu \|\nabla \mathbf{E}_u^{l+1}\|_0^2. \nonumber
\end{align}
$$

Using the inequality $a(a-b) \ge \frac{1}{2}(a^2-b^2)$, we obtain

$$
\begin{align}
T_1 &\ge \frac{1}{2} \left( \|\mathbf{E}_u^{l+1}\|_0^2 - \|\mathbf{E}_u^l \circ \mathbf{X}_h^l\|_0^2 \right) \nonumber \\
&\ge \frac{1}{2} \left( \|\mathbf{E}_u^{l+1}\|_0^2 - \|\mathbf{E}_u^l\|_0^2 \right) + \frac{1}{2} T_3, \tag{7.8}
\end{align}
$$

and

$$
\begin{align}
T_2 &\ge \frac{1}{2} \left( \|E_H^{l+1}\|_{0,l}^2 - \|E_H^l \circ \mathbf{X}_h^l\|_{0,l}^2 \right) \tag{7.9} \\
&\ge \frac{1}{2} \left( \|E_H^{l+1}\|_{0,l}^2 - \|E_H^l\|_{0,l-1}^2 \right) + \frac{1}{2} T_4 + \frac{1}{2} T_5. \nonumber
\end{align}
$$
Finally, thanks to (7.5), summing (7.7) – (7.9) we obtain (7.6). ■

Let us set $\delta_4 = \min(\delta_0, \delta_1, \delta_2, \delta_3)$. The following lemma gives bounds for the terms $T_3, T_4, T_5$.

**Lemma 11** For all $\Delta t \le \delta_4$, the following assertions are satisfied:

i) Let us assume that $(\mathbf{H}_1)$ holds. Then we have

$$ \frac{\alpha}{2} T_3 \le \frac{C_1}{2} \|\mathbf{E}_u^l\|_0^2. $$

ii) Let us assume that $(\mathbf{H}_0)$ and $(\mathbf{H}_1)$ hold. Then we have

$$ \frac{g\alpha}{2} T_5 \le \frac{gC_3}{2} \|E_H^l\|_0^2, $$

where the constants $C_1$ and $C_3$ are given respectively by Proposition 7 and Proposition 9.

iii) Let us assume that $(\mathbf{H}_0)$ and $(\mathbf{H}_2)$ hold. Then there exists a constant $C_4 = C_4(M, \underline{H}_0)$ such that

$$ \frac{g\alpha}{2} T_4 \le \frac{gC_4}{2} \|E_H^l\|_0^2. $$

**Proof.** Assertions i), ii) are immediate from Proposition 9, and iii) is a consequence of the induction hypothesis $(\mathbf{H}_2)$. ■

We now derive a bound for $\bar{r}_1(\mathbf{E}_u^{l+1})$.

**Lemma 12** For a given $\varepsilon_0, \varepsilon_1$ in $\mathbb{R}_+^*$, there exists non-negative constants $C_5 = C_5(\varepsilon_0, \varepsilon_1, \lambda, \mu, M, \underline{H}_0, g)$ such that for all $\Delta t \le \delta_4$,

$$ \begin{aligned} & |\bar{r}_1(\mathbf{E}_u^{l+1})| \\ & \le (\varepsilon_0 + \varepsilon_1) \left\| \nabla \mathbf{E}_u^{l+1} \right\|_0^2 + C_5 \left( \|\mathbf{E}_u^{l+1}\|_0^2 + \|E_H^{l+1}\|_0^2 \right. \\ & \qquad \left. + \|E_H^l\|_0^2 + \|\mathbf{E}_u^l\|_0^2 + \|e_H^{l+1}\|_0^2 + \|e_H^l\|_0^2 + \|\mathbf{e}_u^{l+1}\|_1^2 + \|\mathbf{e}_u^l\|_1^2 \right. \\ & \qquad \left. + \left( \int_{t^l}^{t^{l+1}} \left\| \frac{\partial H}{\partial t} \right\|_0^2 dt + \int_{t^l}^{t^{l+1}} \left\| \frac{\partial \mathbf{u}}{\partial t} \right\|_0^2 dt + \int_{t^l}^{t^{l+1}} \left\| \frac{d^2 \mathbf{u}}{dt^2} \right\|_0^2 dt \right) \Delta t \right). \end{aligned} \quad (7.10) $$

**Proof.** We first derive a bound for $r_1(\mathbf{E}_u^{l+1})$. We check that

$$ \begin{aligned} & \left| \alpha (\mathbf{u}^{l+1} - \mathbf{u}^l \circ \mathbf{X}_h^l, \mathbf{E}_u^{l+1}) - \left( \frac{d\mathbf{u}^{l+1}}{dt}, \mathbf{E}_u^{l+1} \right) \right| \\ & \le \left| \alpha (\mathbf{u}^l \circ \mathbf{X}_h^l - \mathbf{u}^l \circ \mathbf{X}^l, \mathbf{E}_u^{l+1}) \right| + \left| \left( \alpha(\mathbf{u}^{l+1} - \mathbf{u}^l \circ \mathbf{X}^l) - \frac{d\mathbf{u}^{l+1}}{dt}, \mathbf{E}_u^{l+1} \right) \right|. \end{aligned} $$
Using Lemma 8, we get for $\Delta t \le \delta_3$,

$$ \begin{aligned} & |\alpha(\mathbf{u}^l \circ \mathbf{X}_h^l - \mathbf{u}^l \circ \mathbf{X}^l, \mathbf{E}_u^{l+1})| \le \alpha \|\nabla \mathbf{u}^l\|_{\infty} \|\mathbf{X}_h^l - \mathbf{X}^l\|_0 \|\mathbf{E}_u^{l+1}\|_0 \\ & \le M\sqrt{C_2} \left( 2 (\|\mathbf{E}_u^l\|_0^2 + \|\mathbf{e}_u^l\|_0^2) + \Delta t \int_{t^l}^{t^{l+1}} \left\|\frac{\partial \mathbf{u}}{\partial t}\right\|_0^2 dt \right)^{\frac{1}{2}} \|\mathbf{E}_u^{l+1}\|_0. \end{aligned} $$

By standard backward difference error analysis, and thanks to Lemma 4, we establish that, for $\Delta t \le \delta_0$,

$$ \left\| \alpha (\mathbf{u}^{l+1} - \mathbf{u}^l \circ \mathbf{X}^l) - \frac{d\mathbf{u}^{l+1}}{dt} \right\|_0^2 \le \frac{\Delta t}{2} \int_{t^l}^{t^{l+1}} \left\| \frac{d^2 \mathbf{u}}{dt^2} \right\|_0^2 dt; $$

and then, by Cauchy-Shwarz's inequality, we have

$$ \begin{aligned} & \left| \alpha (\mathbf{u}^{l+1} - \mathbf{u}^l \circ \mathbf{X}_h^l, \mathbf{E}_u^{l+1}) - \left( \frac{d\mathbf{u}^{l+1}}{dt}, \mathbf{E}_u^{l+1} \right) \right| \\ & \le (1 + M\sqrt{C_2}) \left( \left( \int_{t^l}^{t^{l+1}} \left\| \frac{d^2 \mathbf{u}}{dt^2} \right\|_0^2 dt + \int_{t^l}^{t^{l+1}} \left\| \frac{\partial \mathbf{u}}{\partial t} \right\|_0^2 dt \right) \Delta t \right. \\ & \left. + \|\mathbf{E}_u^l\|_0^2 + \|\mathbf{E}_u^{l+1}\|_0^2 + \|\mathbf{e}_u^l\|_0^2 \right). \end{aligned} \quad (7.11) $$

We show easily, thanks to Hypothesis (H₀) that,

$$ \begin{aligned} & |a(H^{l+1}, \mathbf{u}^{l+1}, \mathbf{u}^{l+1}, \mathbf{E}_u^{l+1}) - a(H_h^l, \mathbf{u}_h^l, \mathbf{u}^{l+1}, \mathbf{E}_u^{l+1})| \\ &= \left| \lambda \left( \left( \frac{1}{H^{l+1}} - \frac{1}{H_h^l} \right) |\mathbf{u}^{l+1}| |\mathbf{u}^{l+1}, \mathbf{E}_u^{l+1}| \right) + \lambda \left( \frac{(|\mathbf{u}^{l+1}| - |\mathbf{u}_h^l|)}{H_h^l} |\mathbf{u}^{l+1}, \mathbf{E}_u^{l+1}| \right) \right| \\ & \le 6 \left( \frac{\lambda M^2}{\underline{H}_0^2} + \frac{\lambda M}{\overline{H}_0} \right) \left( \|\mathbf{E}_u^{l+1}\|_0^2 + \|\mathbf{E}_u^l\|_0^2 + \|\mathbf{e}_u^l\|_0^2 + \|\mathbf{E}_H^l\|_0^2 + \|\mathbf{e}_H^l\|_0^2 \right. \\ & \left. + \left( \int_{t^l}^{t^{l+1}} \left\| \frac{\partial H}{\partial t} \right\|_0^2 dt + \int_{t^l}^{t^{l+1}} \left\| \frac{\partial \mathbf{u}}{\partial t} \right\|_0^2 dt \right) \Delta t \right). \end{aligned} \quad (7.12) $$

So, taking into account (7.11) and (7.12), we derive the following bound for $r_1(\mathbf{E}_u^{l+1})$, under the condition $\Delta t \le \min(\delta_0, \delta_3)$,

$$ \begin{aligned} r_1(\mathbf{E}_u^{l+1}) \le {}& C_6 (\|\mathbf{E}_u^{l+1}\|_0^2 + \|\mathbf{E}_H^l\|_0^2 + \|\mathbf{E}_u^l\|_0^2 + \|\mathbf{e}_H^l\|_0^2 + \|\mathbf{e}_u^l\|_0^2 \\ & + \Delta t \int_{t^l}^{t^{l+1}} \left( \left\| \frac{\partial H}{\partial t} \right\|_0^2 + \left\| \frac{\partial \mathbf{u}}{\partial t} \right\|_0^2 + \left\| \frac{d^2 \mathbf{u}}{dt^2} \right\|_0^2 \right) dt), \end{aligned} \quad (7.13) $$
with $C_6 = 2 \min \left( 1 + M\sqrt{C_2}, 6 \left( \frac{\lambda M^2}{\underline{H}_0^2} + \frac{\lambda M}{\underline{H}_0} \right) \right)$.

By virtue of Hypotheses $(\mathbf{H}_0) - (\mathbf{H}_1)$, we show that

$$ 
\begin{aligned}
|a(H_h^l, \mathbf{u}_h^l, \mathbf{e}_u^{l+1}, \mathbf{E}_u^{l+1})| &= \mu \left( \nabla \mathbf{e}_u^{l+1}, \nabla \mathbf{E}_u^{l+1} \right) + \lambda \left( \frac{|\mathbf{u}_h^l|}{H_h^l} \mathbf{e}_u^{l+1}, \mathbf{E}_u^{l+1} \right) \\
&\le \mu \left\| \nabla \mathbf{e}_u^{l+1} \right\|_0 \left\| \nabla \mathbf{E}_u^{l+1} \right\|_0 + \frac{2\lambda(1+M)}{\underline{H}_0} \left\| \mathbf{e}_u^{l+1} \right\|_0 \left\| \mathbf{E}_u^{l+1} \right\|_0 \\
&\le \varepsilon_0 \left\| \nabla \mathbf{E}_u^{l+1} \right\|_0^2 + C_7 \left( \left\| \mathbf{e}_u^{l+1} \right\|_1^2 + \left\| \mathbf{E}_u^{l+1} \right\|_0^2 \right),
\end{aligned}
\quad (7.14) $$

where $C_7 = \min \left( \frac{\lambda(1+M)}{\underline{H}_0}, \frac{\mu^2}{4\varepsilon_0} \right)$. We also check that

$$ 
\begin{aligned}
|b(\mathbf{E}_u^{l+1}, \varepsilon_H^{l+1})| &\le g \left( \left| (\nabla \cdot \mathbf{E}_u^{l+1}, e_H^{l+1}) \right| + \left| (\nabla \cdot \mathbf{E}_u^{l+1}, E_H^{l+1}) \right| \right) \\
&\le \varepsilon_1 \left\| \nabla \mathbf{E}_u^{l+1} \right\|_0^2 + \frac{g^2}{2\varepsilon_1} \left( \left\| e_H^{l+1} \right\|_0^2 + \left\| E_H^{l+1} \right\|_0^2 \right).
\end{aligned}
\quad (7.15) $$

Since $\pi_h \mathbf{u}^l$ is the $L^2$ projection of $\mathbf{u}^l$ in $V_h$ and $\mathbf{E}_u^l \in V_h$, we have

$$ (\mathbf{e}_u^l, \mathbf{E}_u^l) = (\pi_h \mathbf{u}^l - \mathbf{u}^l, \mathbf{E}_u^l) = 0, \text{ for all } l = 0, N, $$

and then, thanks to Lemma 6, we get for $\Delta t \le \delta_2$

$$ 
\begin{aligned}
|\alpha (\mathbf{e}_u^{l+1} - \mathbf{e}_u^l \circ X_h^l, \mathbf{E}_u^{l+1})| &= |\alpha (\mathbf{e}_u^l - \mathbf{e}_u^l \circ X_h^l, \mathbf{E}_u^{l+1})| \\
&\le C_0 \left\| \mathbf{E}_u^{l+1} \right\|_0 \left\| \nabla \mathbf{e}_u^l \right\|_0 \\
&\le \frac{C_0}{2} \left( \left\| \mathbf{E}_u^{l+1} \right\|_0^2 + |\mathbf{e}_u^l|_1^2 \right).
\end{aligned}
\quad (7.16) $$

The proof of (7.10) is achieved by summing (7.13) – (7.16) and taking

$$ C_5 = 2 \max \left( \frac{\lambda(1+M)}{\underline{H}_0}, \frac{\mu^2}{4\varepsilon_0}, \frac{g^2}{2\varepsilon_1}, \frac{C_0}{2}, 2(1+M\sqrt{C_2}), 12 \left( \frac{\lambda M^2}{\underline{H}_0^2} + \frac{\lambda M}{\underline{H}_0} \right) \right), $$

with the restriction $\Delta t \le \min(\delta_0, \delta_2, \delta_3)$. ■

The following lemma gives an estimate of $\bar{r}_2(E_u^{l+1})$.
**Lemma 13** For a given $\varepsilon_3$ in $\mathbb{R}_+^*$, there exists non negative constants $C_8 = C_8(\varepsilon_3, M, \underline{H}_0, \underline{H}_1, g)$ such that for all $\Delta t \le \delta_4$,

$$ \begin{aligned} |\bar{r}_2(E_H^{l+1})| &\le \varepsilon_3 \|\nabla \mathbf{E}_u^{l+1}\|_0^2 + C_8 (\|\mathbf{E}_u^l\|_0^2 + \|E_H^{l+1}\|_0^2 + \|E_H^l\|_0^2 \\ &\quad + \|e_H^l\|_1^2 + \|e_u^l\|_0^2 + |e_u^{l+1}|_1^2 + \alpha^2 \|e_H^{l+1} - e_H^l\|_0^2 \\ &\quad + \left( \int_{t^l}^{t^{l+1}} \left\| \frac{d^2 H}{dt^2} \right\|_0^2 dt + \int_{t^l}^{t^{l+1}} \left\| \frac{\partial H}{\partial t} \right\|_0^2 dt \right) \Delta t). \end{aligned} \quad (7.17) $$

**Proof.** Writing

$$ \begin{aligned} |r_2(E_H^{l+1})| &= g \left| \left( \frac{dH^{l+1}}{dt} - \frac{\alpha(H^{l+1} - H^l \circ X_h^l)}{H_h^l}, E_H^{l+1} \right) \right| \\ &\le g \left| \left( \frac{dH^{l+1}}{dt} - \frac{\alpha(H^{l+1} - H^l \circ X_h^l)}{H_h^l}, E_H^{l+1} \right) \right| \\ &\quad + g \left| \left( \left( \frac{1}{H^{l+1}} - \frac{1}{H_h^l} \right) \frac{dH^{l+1}}{dt}, E_H^{l+1} \right) \right|, \end{aligned} $$

the same arguments as before, for obtaining estimate (7.11), yields for the first term of the right-hand side of the above inequality:

$$ g \left| \left( \frac{dH^{l+1}}{dt} - \frac{\alpha(H^{l+1} - H^l \circ X_h^l)}{H_h^l}, E_H^{l+1} \right) \right| \le \frac{2g(1 + M\sqrt{C_2})}{\underline{H}_0} \left( \|E_H^{l+1}\|_0^2 + \|\mathbf{E}_u^l\|_0^2 + \|\mathbf{e}_u^l\|_0^2 + \left( \int_{t^l}^{t^{l+1}} \left\| \frac{d^2 H}{dt^2} \right\|_0^2 dt + \int_{t^l}^{t^{l+1}} \left\| \frac{\partial \mathbf{u}}{\partial t} \right\|_0^2 dt \right) \Delta t \right), $$

with $\Delta t \le \delta_0$. And for the second term we have

$$ \begin{aligned} & \left| \left( \left( \frac{1}{H^{l+1}} - \frac{1}{H_h^l} \right) \frac{\partial H^{l+1}}{\partial t}, E_H^{l+1} \right) \right| \\ &\le \frac{2}{\underline{H}_0^2} \left\| \frac{\partial H^{l+1}}{\partial t} \right\|_{\infty} \|H^{l+1} - H_h^l\|_0 \|E_H^{l+1}\|_0 \\ &\le \frac{2M}{\underline{H}_0^2} \left( \|e_H^l\|_0^2 + \|E_H^l\|_0^2 + \|E_H^{l+1}\|_0^2 + \Delta t \int_{t^l}^{t^{l+1}} \left\| \frac{\partial H}{\partial t} \right\|_0^2 dt \right). \end{aligned} $$
It follows that

$$|r_2(E_H^{l+1})| \le C_9 (\|\mathbf{E}_u^l\|_0^2 + \|E_H^{l+1}\|_0^2 + \|E_H^l\|_0^2 + \|e_H^l\|_0^2 + \|\mathbf{e}_u^l\|_0^2 + \\ \Delta t \int_{t^l}^{t^{l+1}} \left( \left\| \frac{d^2 H}{dt^2} \right\|_0^2 + \left\| \frac{\partial H}{\partial t} \right\|_0^2 + \left\| \frac{\partial \mathbf{u}}{\partial t} \right\|_0^2 \right) dt), \quad (7.18)$$

where $C_9 = 2 \max \left( \frac{2g(1 + M\sqrt{C_2})}{\underline{H}_0}, \frac{2gM}{\underline{H}_0^2} \right)$. In the same manner as before, we show that

$$|b(\varepsilon_u^{l+1}, E_H^{l+1})| \le \varepsilon_3 \left\| \nabla \mathbf{E}_u^{l+1} \right\|_0^2 + \left( \frac{g^2}{4\varepsilon_3} + \frac{g}{2} \right) \left( \|E_H^{l+1}\|_0^2 + \left\| \nabla \mathbf{e}_u^{l+1} \right\|_0^2 \right). \quad (7.19)$$

Thanks to Lemma 6, we have

$$\begin{aligned} & \alpha \left| \left( \frac{e_H^{l+1} - e_H^l \circ X_h^l}{H_h^l}, E_H^{l+1} \right) \right| \\ & \le \alpha \left( \left| \left( \frac{e_H^{l+1} - e_H^l}{H_h^l}, E_H^{l+1} \right) \right| + \left| \left( \frac{e_H^l - e_H^l \circ X_h^l}{H_h^l}, E_H^{l+1} \right) \right| \right) \quad (7.20) \\ & \le \frac{2}{\underline{H}_0} \|E_H^{l+1}\|_0 \left( \alpha \|e_H^{l+1} - e_H^l\|_0 + C_0 \|\nabla e_H^l\|_0 \right) \\ & \le \frac{2(1+C_0)}{\underline{H}_0} \left( \alpha^2 \|e_H^{l+1} - e_H^l\|_0^2 + \|\nabla e_H^l\|_0^2 + \|E_H^{l+1}\|_0^2 \right), \end{aligned}$$

with $\Delta t \le \delta_2$. Summing inequalities (7.18)-(7.20), Assertion (7.17) is proved with the choice

$$C_8 = 3 \max \left( \frac{2g(1 + M\sqrt{C_2})}{\underline{H}_0}, \frac{2gM}{\underline{H}_0^2}, \frac{g^2}{4\varepsilon_3} + \frac{g}{2}, \frac{2(1 + C_0)}{\underline{H}_0} \right),$$

and the restriction $\Delta t \le \min(\delta_0, \delta_2)$. ■

We now give bound for the therm $r(\mathbf{E}_u^{l+1})$.

**Lemma 14** There exists a constant $C_{10} = C_{10}(\|\mathbf{1}\|_{\infty}, \bar{\tau}, \rho, \underline{H}_0)$ such that

$$r(\mathbf{E}_u^{l+1}) \le C_{10} (\|\mathbf{E}_u^{l+1}\|_0^2 + \|\mathbf{E}_u^l\|_0^2 + \|\mathbf{e}_u^l\|_0^2 + \|E_H^l\|_0^2 + \|e_H^l\|_0^2 \\ + \Delta t \int_{t^l}^{t^{l+1}} \left\| \frac{\partial \mathbf{u}}{\partial t} \right\|_0^2 dt + \Delta t \int_{t^l}^{t^{l+1}} \left\| \frac{\partial H}{\partial t} \right\|_0^2 dt) \quad (7.21)$$
**Proof.** Similarly as before, we show that

$$
\begin{aligned}
r(\mathbf{E}_u^{l+1}) &= \left\langle -g\nabla z_b - \frac{\nabla p_a}{\rho} + \frac{\bar{\tau}}{\rho H_h^l} + \mathbf{l} \times \mathbf{u}_h^l \right. \\
&\quad \left. - \left( -g\nabla z_b - \frac{\nabla p_a}{\rho} + \frac{\bar{\tau}}{\rho H^{l+1}} + \mathbf{l} \times \mathbf{u}^{l+1} \right), \mathbf{E}_u^{l+1} \right\rangle \\
&\le \|\mathbf{E}_u^{l+1}\|_0 \left( \|\mathbf{l}\|_{\infty} \|\mathbf{u}_h^l - \mathbf{u}^{l+1}\|_0 + \frac{2\bar{\tau}}{\rho H_0^2} \|H_h^l - H^{l+1}\|_0 \right) \\
&\le 3C_{11} \left( \|\mathbf{E}_u^{l+1}\|_0^2 + \|\mathbf{E}_u^l\|_0^2 + \|\mathbf{e}_u^l\|_0^2 + \|\mathbf{u}^l - \mathbf{u}^{l+1}\|_0^2 \right. \\
&\quad \left. + \|E_H^l\|_0^2 + \|e_H^l\|_0^2 + \|H^l - H^{l+1}\|_0^2 \right) \\
&\le 3C_{11} \left( \|\mathbf{E}_u^{l+1}\|_0^2 + \|\mathbf{E}_u^l\|_0^2 + \|\mathbf{e}_u^l\|_0^2 + \|E_H^l\|_0^2 + \|e_H^l\|_0^2 \right. \\
&\quad \left. + \Delta t \int_{t^l}^{t^{l+1}} \left\| \frac{\partial \mathbf{u}}{\partial t} \right\|_0^2 dt + \Delta t \int_{t^l}^{t^{l+1}} \left\| \frac{\partial H}{\partial t} \right\|_0^2 dt \right),
\end{aligned}
$$

where $C_{11} = \max(\|\mathbf{l}\|_{\infty}, \frac{2\bar{\tau}}{\rho H_0})$, and then it is sufficient to take $C_{10} = 3C_{11}$. ■

We are now able to prove the main result of this paper.

**Theorem 15** Let $(H_h^{n+1}, \mathbf{u}_h^{n+1})$ be the numerical solution approximating $(H, \mathbf{u})$ at the time $t = t^{n+1}$. Assume that $(H_0) - (H_2)$ hold. Then, there exists a constant $K = K(\Omega, T, k, m, \underline{H}_0, \overline{H}_1, M, \mu, g, \lambda, \|\mathbf{l}\|_{\infty}, \rho)$ such that, for $\Delta t \le \delta$ and $h \le 1$, we have

$$
\|\mathbf{u}_h - \mathbf{u}\|_{0,\infty} + \|H_h - H\|_{0,\infty} + \sqrt{\mu} \|\mathbf{u}_h - \mathbf{u}\|_{1,0} \le K(\Delta t + h^k + h^m). \quad (7.22)
$$

**Proof.** Choosing $\varepsilon_0 = \varepsilon_1 = \varepsilon_2 = \varepsilon_3 = \frac{\mu}{8}$, multiplying (7.6) by $\Delta t$, and summing over $l = 1, \dots, n$, and taking into account the bounds on
$T_3, T_4, T_5, \bar{r}_1(\mathbf{E}_u^{l+1}), \bar{r}_2(\mathbf{E}_u^{l+1}),$ and $r(\mathbf{E}_u^{l+1})$ we get

$$
\begin{align}
& \frac{1}{2} (\|\mathbf{E}_u^{n+1}\|_0^2 - \|\mathbf{E}_u^1\|_0^2) + \frac{g}{2} (\|E_H^{n+1}\|_{0,n}^2 - \|E_H^1\|_{0,0}^2) + \frac{\mu}{2} \sum_{l=1,n} \Delta t |\mathbf{E}_u^{l+1}|_1^2 \nonumber \\
& \le C_{12} \sum_{l=1,n} \left( \int_{t^l}^{t^{l+1}} \left\| \frac{\partial \mathbf{u}}{\partial t} \right\|_0^2 dt + \int_{t^l}^{t^{l+1}} \left\| \frac{\partial H}{\partial t} \right\|_0^2 dt \right. \nonumber \\
& \qquad \left. + \int_{t^l}^{t^{l+1}} \left\| \frac{d^2 \mathbf{u}}{dt^2} \right\|_0^2 dt + \int_{t^l}^{t^{l+1}} \left\| \frac{d^2 H}{dt^2} \right\|_0^2 dt \right) \Delta t^2 \nonumber \\
& + C_{12} \left( \sum_{l=0,n} \Delta t (\|\mathbf{e}_u^{l+1}\|_1^2 + \|e_H^{l+1}\|_1^2) + \sum_{l=1,n} \alpha \|e_H^{l+1} - e_H^l\|_0^2 \right) \nonumber \\
& + C_{12} \sum_{l=0,n} \Delta t (\|\mathbf{E}_u^{l+1}\|_0^2 + \|E_H^{l+1}\|_0^2), \tag{7.23}
\end{align}
$$

where $C_{12} = 6 \max \left( \frac{C_1}{2}, \frac{gC_3}{2}, \frac{gC_4}{2}, C_5, C_8, C_{10} \right)$. As $\mathbf{E}_u^0 = 0$ and $E_H^0 = 0$,
(7.5) implies by using the same reasoning as before, for $l=0$, that

$$
\alpha \|\mathbf{E}_u^1\|_0^2 + g\alpha \|E_H^1\|_{0,0}^2 + \mu \|\mathbf{E}_u^1\|_1^2 \le \bar{r}_1(\mathbf{E}_u^1) + \bar{r}_2(E_H^1) + r(\mathbf{E}_u^1),
$$

and then similarly we get

$$
\begin{align}
& \alpha \|\mathbf{E}_u^1\|_0^2 + g\alpha \|E_H^1\|_{0,0}^2 + \frac{\mu}{2} |\mathbf{E}_u^1\|_1^2 \nonumber \\
& \le C_{13} \left( \|\mathbf{e}_u^0\|_1^2 + \|e_H^0\|_1^2 + \|\mathbf{e}_u^1\|_1^2 + \|e_H^1\|_1^2 + \|\mathbf{E}_u^1\|_0^2 + \|E_H^1\|_0^2 \right. \nonumber \\
& \qquad \left. + \alpha^2 \|e_H^1 - e_H^0\|_0^2 + \left( \int_{t^0}^{t^1} \left\| \frac{\partial \mathbf{u}}{\partial t} \right\|_0^2 dt + \int_{t^0}^{t^1} \left\| \frac{\partial H}{\partial t} \right\|_0^2 dt \right) \Delta t \right. \nonumber \\
& \qquad \left. + \left( \int_{t^0}^{t^1} \left\| \frac{d^2 \mathbf{u}}{dt^2} \right\|_0^2 dt + \int_{t^0}^{t^1} \left\| \frac{d^2 H}{dt^2} \right\|_0^2 dt \right) \Delta t \right), \tag{7.24}
\end{align}
$$

with $C_{13} = 3 \max \left( \frac{C_1}{2}, \frac{gC_3}{2}, \frac{gC_4}{2}, C_5, C_8, C_{10} \right)$ and the same choice of $\epsilon_0$,
$\epsilon_1$, $\epsilon_2$ and $\epsilon_3$.

Using Hypothesis (H₀), multiplying (7.24) by Δt/4 and summing the corre-
sponding inequality with (7.23), we obtain

$$
\begin{align}
& \|\mathbf{E}_u^{n+1}\|_0^2 + \|E_H^{n+1}\|_0^2 + \mu \sum_{l=0,n} \Delta t \|\nabla \mathbf{E}_u^{l+1}\|_0^2 \nonumber \\
& \le C_{14} \sum_{l=0,n} \left( \int_{t^l}^{t^{l+1}} \left\| \frac{\partial \mathbf{u}}{\partial t} \right\|_0^2 dt + \int_{t^l}^{t^{l+1}} \left\| \frac{\partial H}{\partial t} \right\|_0^2 dt \right. \nonumber \\
& \qquad \left. + \int_{t^l}^{t^{l+1}} \left\| \frac{d^2 \mathbf{u}}{dt^2} \right\|_0^2 dt + \int_{t^l}^{t^{l+1}} \left\| \frac{d^2 H}{dt^2} \right\|_0^2 dt \right) \Delta t^2 \nonumber \\
& + C_{14} \left( \sum_{l=0,n+1} \Delta t \left( \|\mathbf{e}_u^l\|_1^2 + \|e_H^l\|_1^2 \right) + \sum_{l=0,n} \alpha \|e_H^{l+1} - e_H^l\|_0^2 \right) \nonumber \\
& + C_{14} \sum_{l=0,n+1} \Delta t \left( \|\mathbf{E}_u^l\|_0^2 + \|E_H^l\|_0^2 \right), \tag{7.25}
\end{align}
$$

where $C_{14} = \frac{3C_{12}}{\min\left(\frac{1}{2}, \frac{g}{\frac{H_0}{2} + H_1}\right)}$. From the definition of the projection

operator $\pi_h$, and thanks to the approximation properties (7.1), we prove
that

$$
\begin{align}
\alpha \|e_H^{l+1} - e_H^l\|_0^2 &= \alpha \|H^{l+1} - H^l - \pi_h (H^{l+1} - H^l)\|_0^2 \nonumber \\
&\le \alpha c_0 h^{2(m+1)} |H^{l+1} - H^l|_{m+1}^2 \tag{7.26} \\
&\le c_0 h^{2(m+1)} \left( \int_{t^l}^{t^{l+1}} \left| \frac{\partial H}{\partial t} \right|_{m+1}^2 dt \right). \nonumber
\end{align}
$$

Using again the approximation properties (7.1) and (7.26), inequality (7.25)
implies, for $h \le 1$:

$$
\begin{align*}
& \|\mathbf{E}_u^{n+1}\|_0^2 + \|E_H^{n+1}\|_0^2 + \mu \sum_{l=1,n+1} \Delta t \|\nabla \mathbf{E}_u^l\|_0^2 \\
& \le C_{14} K_1 \Delta t^2 + C_{14} K_2 (h^{2k} + h^{2m} + h^{2(m+1)}) \\
& \quad + C_{14} \sum_{l=0,n+1} \Delta t (\|\mathbf{E}_u^l\|_0^2 + \|E_H^l\|_0^2) \\
& \le C_{15} \left( (\Delta t^2 + h^{2k} + h^{2m}) + \sum_{l=0,n+1} \Delta t (\|\mathbf{E}_u^l\|_0^2 + \|E_H^l\|_0^2) \right),
\end{align*}
$$
where $C_{15} = \max(1, K_1, 4K_2) C_{14}$,

$$K_1 = \left\| \left\| \frac{\partial \mathbf{u}}{\partial t} \right\| \right\|_0^2 + \left\| \left\| \frac{\partial H}{\partial t} \right\| \right\|_0^2 + \left\| \left\| \frac{d^2 \mathbf{u}}{dt^2} \right\| \right\|_0^2 + \left\| \left\| \frac{d^2 H}{dt^2} \right\| \right\|_0^2,$$

$$K_2 = c_0 \left( \left\| \mathbf{u} \right\|_{k+1,0}^2 + \left\| H \right\|_{m+1,0}^2 + \left\| \frac{\partial H}{\partial t} \right\|_{m+1,0}^2 \right),$$

Applying the discrete Gronwall's formula, we get

$$\| \mathbf{E}_u^{n+1} \|_0^2 + \| \mathbf{E}_H^{n+1} \|_0^2 + \mu \sum_{l=1,n+1} \Delta t \left\| \nabla \mathbf{E}_u^l \right\|_0^2 \le K_3 (\Delta t^2 + h^{2k} + h^{2m})$$

with $\Delta t \le \delta$, $\delta = \min(\delta_4, \frac{1}{C_{15}})$, and $K_3 = C_{15} \exp(C_{15}T)$. Which implies, thanks to the approximation properties (7.1), Assertion (7.22). Where

$$K = 2\sqrt{K_3} + c_0 \left( \sqrt{\mu} \|\mathbf{u}\|_{k+1,0} + \|\mathbf{u}\|_{k+1,\infty} + \|H\|_{m+1,\infty} \right),$$

and then this achieve the proof of Theorem 15. ■

**Theorem 16** Assume that $(\mathbf{H}_0) - (\mathbf{H}_2)$ hold. Then we have the following assertions:

i) There exists constants $c_1$ and $d_1$ such that if $\Delta t \le c_1 h^{\frac{d}{2}}$, $h \le d_1$ and $k, m > \frac{d}{2}$, then the boundedness induction hypotheses on $H_h^l$ and $\mathbf{u}_h^l$ are true for $l = n + 1$.

ii) There exists constants $c_2$ and $d_2$ such that if $\Delta t \le c_2 h^{\frac{d}{2}+1}$, $h \le d_2$ and $k, m > 1 + \frac{d}{2}$, then the boundedness induction hypotheses on $\nabla H_h^l$, $\nabla \mathbf{u}_h^l$ are true for $l = n + 1$.

iii) There exists constants $c_3$ and $d_3$ such that if $\Delta t \le c_3 h^{\frac{d}{2}+1}$, $h \le d_3$ and $k, m > 1 + \frac{d}{2}$, then the boundedness induction hypothesis on $\alpha(H_h^l - H_h^{l-1})$ is true for $l = n + 1$.

**Proof.** Let us set

$$\omega_0 = \left( \frac{\Delta t}{h^{\frac{d}{2}}} + h^{k-\frac{d}{2}} + h^{m-\frac{d}{2}} \right), \omega_1 = \frac{\omega_0}{h}.$$

We first prove that the boundedness induction hypotheses on $H_h^l$ and $\mathbf{u}_h^l$ for $l = n + 1$. From the inverse estimates (7.2) and (7.22) we show that

$$\|E_H^{n+1}\|_{\infty} \le c h^{-\frac{d}{2}} \|E_H^{n+1}\|_0 \le c \omega_0 K,$$
where K is given by Theorem 15. It follows, using the approximation properties (7.1), that

$$ 
\begin{aligned}
& -|E_H^{n+1}| - |e_H^{n+1}| + H^{n+1} \le H_h^{n+1} \le H^{n+1} + |E_H^{n+1}| + |e_H^{n+1}| \\
& -\omega_2 + \underline{H}_0 \le \underline{H}_1 + \omega_2
\end{aligned}
$$

where $\omega_2 = c\omega_0K + c_0Mh^{m+1-\frac{d}{2}}$. Let us set

$$ c_4 = \frac{H_0}{8} \min \left( \frac{1}{cK}, \frac{1}{c_0 M} \right), d_4 = \min \left( 1, \left( \frac{c_4}{\delta} \right)^{\frac{2}{d}}, c_4^{\frac{2}{2m-d}}, c_4^{\frac{2}{2k-d}} \right). $$

Then for $\Delta t \le c_4h^{\frac{d}{2}}$, $h \le d_4$ and $k, m > \frac{d}{2}$, we have $\Delta t \le \delta$ and

$$ \omega_2 \le \frac{H_0}{2}, $$

which implies:

$$ \frac{H_0}{2} \le H_h^{n+1} \le \underline{H}_1 + \frac{H_0}{2}. $$

By the same reasoning we prove the boundedness induction hypothesis on $\mathbf{u}_h^l$ for $l = n + 1$:

$$ \|\mathbf{u}_h^{n+1}\|_{\infty} \le 1 + M, $$

under restrictions: $\Delta t \le c_5h^{\frac{d}{2}}$, $h \le d_5$ and $k, m > \frac{d}{2}$. So Assertion i) is ensured by taking $c_1 = \min(c_4, c_5)$, and $d_1 = \min(d_4, d_5)$.

We now prove the boundedness induction hypotheses on $\nabla \mathbf{u}_h^l$ for $l = n+1$. By the inverse estiamtes (7.2) we have

$$ 
\begin{aligned}
\|\nabla \mathbf{E}_u^{n+1}\|_{\infty} &\le |\mathbf{E}_u^{n+1}|_{1,\infty} \\
&\le c h^{-\frac{d}{2}-1} \|\mathbf{E}_u^{n+1}\|_0 \\
&\le c \omega_1 K;
\end{aligned}
$$

which implies that

$$ 
\begin{aligned}
\|\nabla \mathbf{u}_h^{n+1}\|_{\infty} &\le \|\nabla \mathbf{u}^{n+1}\|_{\infty} + |\mathbf{E}_u^{n+1}|_{1,\infty} + |\mathbf{e}_u^{n+1}|_{1,\infty} \\
&\le M + c\omega_1 K + c_0 M h^{k-\frac{d}{2}};
\end{aligned}
$$

Let us set

$$ c_6 = \frac{1}{4} \min \left( \frac{1}{cK}, \frac{1}{c_0 M} \right), d_6 = \min \left( 1, \left( \frac{c_6}{\delta} \right)^{\frac{2}{d+2}}, c_6^{\frac{2}{2(m-1)-d}}, c_6^{\frac{2}{2(k-1)-d}} \right). $$
Then for $\Delta t \le c_6 h^{\frac{d}{2}+1}$, $h \le d_6$ and $k, m > \frac{d}{2}+1$ we show that $\Delta t \le \delta$ and

$$c\omega_1 K + c_0 M h^{k-\frac{d}{2}} \le 1,$$

and then

$$\|\nabla \mathbf{u}_h^{n+1}\|_{\infty} \le 1 + M.$$

In the same manner, we prove the boundedness induction hypothesis on $\nabla H_h^l$, for $l = n + 1$, under restrictions: $\Delta t \le c_7 h^{1+\frac{d}{2}}$, $h \le d_7$ and $k, m > 1+\frac{d}{2}$. And then Assertion ii) is satisfied by choosing $c_2 = \min(c_6, c_7)$ and $d_2 = \min(d_6, d_7)$.

It remains to prove the boundedness induction hypothesis on $\alpha(H_h^l - H_h^{l-1})$ for $l = n + 1$. Subtracting (7.4) from (7.3), with $\mathbf{v} = 0$, $q = E_H^{l+1} - E_H^{l+1}$, we get in the same manner as before:

$$g\alpha (E_H^{n+1} - E_H^n \circ \mathbf{X}_h^n, E_H^{n+1} - E_H^n)_{0,n} = \bar{r}_2 (E_H^{n+1} - E_H^n)$$

with

$$\bar{r}_2 (E_H^{n+1} - E_H^{n+1}) = -r_2(E_H^{n+1} - E_H^n) + b (\varepsilon_u^{n+1}, E_H^{n+1} - E_H^n) \\ + g\alpha (e_H^{n+1} - e_H^n \circ \mathbf{X}_h^n, E_H^{n+1} - E_H^n)_{0,n}.$$

The same reasoning of Lemma 13, implies

$$\begin{aligned} & |\bar{r}_2 (E_H^{n+1} - E_H^n)| \\ & \le C_{16} \|E_H^{n+1} - E_H^n\|_0 \left( \left( 2 (\|\mathbf{E}_u^n\|_0 + \|\mathbf{e}_u^n\|_0)^2 + \Delta t \int_{t^n}^{t^{n+1}} \left\| \frac{\partial \mathbf{u}}{\partial t} \right\|_0^2 dt \right)^{\frac{1}{2}} \right. \\ & \qquad \left. + \left( \Delta t \int_{t^n}^{t^{n+1}} \left\| \frac{d^2 H}{dt^2} \right\|_0^2 dt \right)^{\frac{1}{2}} + \sqrt{2} (\|\mathbf{E}_u^n\|_1 + \|\mathbf{e}_u^n\|_1 + \|e_H^n\|_1) \right. \\ & \qquad \left. + \left( \Delta t \int_{t^n}^{t^{n+1}} \left\| \frac{\partial H}{\partial t} \right\|_0^2 dt \right)^{\frac{1}{2}} + \left( \Delta t \int_{t^n}^{t^{n+1}} \left\| \frac{\partial \mathbf{u}}{\partial t} \right\|_0^2 dt \right)^{\frac{1}{2}} \right. \\ & \qquad \left. + \|E_H^n\|_0 + \alpha \|e_H^{n+1} - e_H^n\|_0 \right), \end{aligned} \quad (7.27)$$

where $C_{16} = \max \left( g, \frac{2g}{\bar{H}_0}, \frac{2Mg}{\bar{H}_0^2}, \frac{2Mg\sqrt{C_2}}{\bar{H}_0}, \frac{2gC_0}{\bar{H}_0}, \frac{2\lambda M}{\bar{H}_0}, \frac{2\lambda M^2}{\bar{H}_0^2} \right)$. Moreover we show as before using Lemma 6 that

$$g\alpha \left| (E_H^n - E_H^n \circ \mathbf{X}_h^n, E_H^{n+1} - E_H^n)_{0,n} \right| \le \frac{2gC_0}{\bar{H}_0} \|E_H^{n+1} - E_H^n\|_0 \|\nabla E_H^n\|_0. \quad (7.28)$$
Combining (7.27) and (7.28) and thanks to the regularity of the exact solution we get

$$ 
\begin{aligned}
\alpha \|E_H^{n+1} - E_H^n\|_0 &\le \frac{(2\underline{H}_1 + \underline{H}_0) C_{16}}{2g} \left( \left(M^2 \Delta t^2 + (\|\mathbf{E}_u^n\|_0 + \|\mathbf{e}_u^n\|_0)^2\right)^{\frac{1}{2}} \right. \\
&\qquad \left. + 3M\Delta t + \sqrt{2} (\|\mathbf{E}_u^n\|_1 + \|E_H^n\|_1 + \|\mathbf{e}_u^n\|_1 + \|e_H^n\|_1) \right. \\
&\qquad \left. + \alpha \|e_H^{n+1} - e_H^n\|_0 \right)
\end{aligned}
$$

Let us set

$$ C_{17} = \frac{(2\underline{H}_1 + \underline{H}_0) C_{16}}{2g} \left( \sqrt{M^2 + (K + c_0 M)^2} + (3 + (1 + 2\sqrt{2})c_0)M \right. \\ \left. + 2\sqrt{(1 + c^2)K} \right). $$

Then, according to Theorem 15, the approximation properties (7.1) and the inverses estimates (7.2), we establish for $h \le 1$ that

$$ \alpha \|E_H^{n+1} - E_H^n\|_0 \le C_{17} \left( \frac{\Delta t}{h} + h^{k-1} + h^{m-1} \right), $$

Finally in the same manner as before we show, for $h \le 1$, that

$$ 
\begin{aligned}
\alpha |H_h^{n+1} - H_h^n| &\le \alpha |H^{n+1} - H^n| + \alpha |E_H^{n+1} - E_H^n| + \alpha |e_H^{n+1} - e_H^n| \\
&\le M + cC_{17}\omega_1 + c_0Mh^{m+1-\frac{d}{2}}.
\end{aligned}
$$

Let us set

$$ c_3 = \frac{1}{4}\min\left(\frac{1}{cC_{17}}, \frac{1}{c_0M}\right), \quad d_3 = \min\left(1, \left(\frac{c_3}{\delta}\right)^{\frac{2}{d+2}} c_3^{\frac{2(m-1)-d}{2}}, c_3^{\frac{2(k-1)-d}{2}}\right). $$

Then for $\Delta t \le c_3h^{-1+\frac{d}{2}}$, $h \le d_3$ and $k, m > 1 + \frac{d}{2}$ we show that $\Delta t \le \delta$ and:

$$ cC_{17}\omega_1 + c_0Mh^{m+1-\frac{d}{2}} \le 1 $$

and then

$$ \|\nabla \mathbf{u}_h^{n+1}\|_{\infty} \le 1 + M. $$

This achieve the proof of the last assertion of the theorem. ■

## 8. Numerical results

We present in this section some numerical results, for validating a priori error estimates of Theorem 15. We first fix the time step $\Delta t$, $\Delta t = 10^{-2}$,
we study the error behaviour on the exact solution **u**, **H**. We consider a rectangular domain $\Omega = [0, X_L] \times [0, Y_L]$, with $X_L = 100$, $Y_L = 10$, and sequence of uniform triangular meshes denoted $M_1, M_2, M_4$, whose properties are illustrated in the following table:

<table><caption>Table 1: Uniform meshes $M_1, M_2, M_4$
We are interested to solve numerically the following problem:</caption><thead><tr><th>Mesh</th><th>Nodes number</th><th>elements number</th><th>unknowns number</th></tr></thead><tbody><tr><td>$M_1$</td><td>121</td><td>200</td><td>363</td></tr><tr><td>$M_2$</td><td>441</td><td>800</td><td>1323</td></tr><tr><td>$M_4$</td><td>1681</td><td>3200</td><td>5043</td></tr></tbody></table>

$$ \begin{gathered} \frac{\partial H}{\partial t} + \nabla \cdot (H\mathbf{u}) = F_H, \tag{8.1} \\ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} - \Delta\mathbf{u} + \nabla H = \mathbf{F}_u, \end{gathered} $$

with Dirichlet boundary conditions on **u**, using $P_1$ finite element spaces for the approximation of **u** and **H**, and characteristics method for the time discretization.

The analytical solutions considered in these numerical experiments are given by:

$$ \begin{gathered} \mathbf{u} = c(t) \left( \sin\left(\frac{\pi x}{X_L}\right) + \sin\left(\frac{\pi y}{Y_L}\right), 10^{-1}\sin\left(\frac{\pi y}{Y_L}\right) \right), \\ H = 2. + \frac{1}{4} \left( -\cos(\omega t) + \sin\left(\frac{4\pi x}{X_L}\right) + \sin\left(\frac{2\pi y}{Y_L}\right) \right), \end{gathered} $$

where $\omega = \frac{\pi}{2T}$, $c(t) = 1. + \sin(\omega t)$.

The second members $F_H$ and $F_u$ are calculated from (8.1). The resulting system of linear equations approximating (8.1) is solved on the three meshes $M_1, M_2$ and $M_4$ with the Uzawa-bicgstab preconditioned gradient method. The errors on the analytical solutions **u**, and **H** are measured with the norms:

$$ \|\varepsilon_H^n\|_0, \|\varepsilon_u^n\|_1, \text{ with } n = 1, \dots, N, N = \frac{T}{\Delta t}, \text{ and } T = 10. $$

So Theorem 15 implies

$$ \|\varepsilon_H^n\|_0 + \|\varepsilon_u^n\|_1 \le K (\Delta t + h). \tag{8.2} $$

For selected times $t_1 = 1, t_2 = 2, t_3 = 5, t_4 = 9$, and $t_5 = 10$, we validate a priori estimates (8.2). In figures (2-3) we plot respectively the
Logarithmic of the $L^2$ norm error on $H$ and the $H^1$ norm error on $u$ as functions of the mesh parameter $h$ for the corresponding time iterations to $t_i$, $i = 1, \dots, 5$. Tables (2-3) give the numerical values of the $h$ exponent in a priori estimates.

We now analyze the $L^2$ norm error on $H$ and $H^1$ norm error on $u$, as functions of $\Delta t$ on the fixed mesh $M_2$. Four values of $\Delta t$ are considered:

$$ \Delta t_1 = \frac{1}{10}, \quad \Delta t_2 = \frac{1}{20}, \quad \Delta t_4 = \frac{1}{40}, \quad \Delta t_8 = \frac{1}{80}. $$

Figures (4-5) show the error behaviour on $H$ and $u$, with respect to the time step $\Delta t$. Tables (4-5) give the numerical values of the $\Delta t$ exponent. For illustration, we plot in figures (8-15) the exact and analytical solutions on the three meshes $M_1$, $M_2$ and $M_4$ at the time $t_3 = 5$. This gives an idea on the convergence of the numerical solution when the parameter mesh $h$ goes to zero. Figures (16-25) show the time evolution of the numerical solutions on the mesh $M_2$. The selected times, for this example, are $t_1 = 1$, $t_2 = 2$, $t_3 = 5$, $t_4 = 9$, and $t_5 = 10$.

From figures (2-3) and tables (2-3), we observe a supra convergence of the numerical solutions with respect to $h$. Figures (4-5) and tables (4-5) show that the $\Delta t$ exponent for the error on $H$ becomes small between $\Delta t_4$ and $\Delta t_8$. So the error on the exact solution is not improved if ones continues to decrease the time step $\Delta t$, because numerical convergence is established. If the mesh is more refined, better convergence occurs by optimizing the time step $\Delta t$, see for example figures (6-7) and tables (6-7) in the case where the fixed mesh is $M_4$.

## 9. Conclusion and remarks

We have established error estimates that prove the convergence of the numerical solution under the condition $\Delta t \le ch^{1+\frac{d}{2}}$. Some consistent restrictions concerning the time step $\Delta t$ are used in the lemmas mentioned in Section 7. The numerical experiments carried out allows us to validate the error estimates for selected time iterations with respect to the mesh parameter $h$. Furthermore, the same analysis with respect to the time step $\Delta t$ leads some numerical difficulties, for ensuring more precisely the error estimates. This is due undoubtedly to the strong nonlinearity of the advectif terms in shallow water equations.

## References

1. M. Abdelwahed and F. El Dabaghi and D.Ouazar (2002). A virtual numerical simulator for aeration effects in lake eutrophication. Inter-
national Journal of Computational Fluid Dynamics. 16 (2), 119-128.

2. V.I. Agoshkov, A. Quarteroni, F., Saleri (1992). Recent developments in the numerical simulation of shallow water equations. I-Boundary conditions, Diparimento di Mathematica, Poletcnico di Milano, n. 80/P.,

3. J.P. Benque, A. Haugel, P.L. Viollet (1982). Numerical methods in environmental fluid mechanics. In Engineering Applications of Com-putational Hydraulics, ed. M.B. Abbott and J.A. Cunge, Vol. II.

4. M. Bercovier, O. Pironneau (1982). Characteristics and the finite element method, in proceeding of the 4th international syposium of Finite Element methods in Flow Problems, T. Kawai, ed., North-Holland, Amsterdam, 67-73.

5. S.C. Brenner, L.R. Ridgway Scott (1994). The Mathematical Theory of Finite element Methods, Texts in Appl. Math., Vol. 15, Springer-Verlag, New York.

6. F. Brezzi, M. Fortin (1991). Mixed and hybrid finite element method. Springer Series in Computational Mathematics, Springer Verlag New York (1991).

7. Cea J. (1971). Optimisation.,Théorie et algorithmes. Dunod.

8. A.J. Chorin, J.E.Marsden (1984). A Mathematical introduction of fluid mechanics. Berlin Heidelberg, New York, Springer.

9. C.N. Dawson, M.L. Martínez-Canales (2000). Characteristic-Galerkin approximation to a system of shallow water equation. Numer. Math. 86, 239-256.

10. J. Douglas, T.F. Russell (1982). Numerical methods for convection dominated diffusion problems based on combining the method of characteristics with finite element methods or finite difference method. Siam J. Numer. Anal. 19, 871-885.

11. F. El Dabaghi, A. El Kacimi and B. Nakhle (2004). Characteristics time discretization and mixed finite element approximation for Saint-Venant equations: a priori error estimates. Rapport de recherche INRIA, http://www.inria.fr/RRRT/RR-???.html.

12. F. El Dabaghi, A. El Kacimi and B. Nakhle (2003). A Priori Error Analysis of the Characteristics Mixed Finite Element Method for Shallow Water Equations. Proceeding 22nd IASTED International Conference, Modelling, Identification, and Control, MIC'2003, 459-464, Innsbruck, Austria.
13. F. El Dabaghi, A. El Kacimi and B. Nakhle (2003). Quelques préconditionneurs par blocs pour des problèmes de point-selle issus d'un modèle numérique d'eau peu profonde. Actes du 1er TAM-TAM (Tendances pour les Applications des Mathématiques, Tunisie-Algérie-Maroc), Rabat, Maroc.

14. F.El Dabaghi (2000). Numerical Aspects of Aeration Process Modelling in Eutrophised Water Basins. J. of Systems Analysis Modelling Simulation. 39, 1-23.

15. R.E. Ewing, T.F. Russell, M.F. Wheeler (1984), Convergence analysis of an approximation of miscible displacement in porous media by mixed finite elements and a modified method of characteristics. Comp. Methods. in Appl. Mech. and Eng. 47, 73-92, North-Holland.

16. J.C. Galland, N. Goutal, J.M. Hervouet (1991). TELEMAC : A new numerical model for solving shallow water equations. Adv. Water resources, 14 (3), 183-148.

17. V. Girault and P.A. Raviart (1986). Finite element methods for Navier-Stokes equations, Springer Ser. Compt. Math 5, Springer-Verlag, New York.

18. N. Goutal (1987). Résolution des équations de Saint-Venant en régime transcritique par une méthode d'éléments finis. Application aux bancs découvrants. Thèse de doctorat de l'université de Paris VI.

19. R.A. Luettich, J.J. Westerink, Scheffner, N.W. (1991). ADCIRC : An advanced three dimensional circulation model for shelves, coasts, and estuaries. Technical report 1, Department of the army, U.S. Army corps of Engineers, Washington, D.C. 20314-1000.

20. O. Pironneau (1982). On the transport-diffusion algorithm and its application to the Navier-Stokes equations. Numer. Math. 83, 309-332.

21. O. Pironneau (1988). Méthodes des éléments finis pour les fluides, Masson.

22. T.F. Russell (1985). Time stepping along characteristics with incomplete iteration for a Galerkin approximation of miscible displacement in porous media. Siam J., Numer. Anal., 22, 167-183.

23. E. Süli (1988). Stability and convergence of the Lagrange-Galerkin method for Navier-Stokes equations, Num. Math., 53 459-483.
24. O.C. Zienkiewicz, P. Ortiz (1995). A split-characteristic based finite element model for the shallow water equations. Int. J. for Numer. Methods in Fluids 20, 1061-1080.
Figure 2: Log of the $L^2$ norm error on $H$ with respect to $h$

Figure 4: Log of the $L^2$ norm error on $H$ with respect to $\Delta t$, for the fixed mesh $M_2$

Figure 6: Log of the $L^2$ norm error on $H$ with respect to $\Delta t$, for the fixed mesh $M_4$

Figure 8: Exact $H$ at $t = t_3$ sur $M_4$

Figure 10: Calculated $H$ at $t = t_3$ on $M_4$

Figure 12: Calculated $H$ at $t = t_3$ on $M_2$

Figure 14: Calculated $H$ at $t = t_3$ on $M_1$

Figure 16: Calculated $H$ at $t = t_1$ on $M_2$

Figure 18: Calculated $H$ at $t = t_2$ on $M_2$

Figure 20: Calculated $H$ at $t = t_3$ on $M_2$

Figure 22: Calculated $H$ at $t = t_4$ on $M_2$

Figure 24: Calculated $H$ at $t = t_5$ on $M_2$

Figure 3: Log of the $H^1$ norm error on $u$ with respect to $h$

Figure 5: Log of the $H^1$ norm error on $u$ with respect to $\Delta t$, for the fixed mesh $M_2$

Figure 7: Log of the $H^1$ norm error on $u$ with respect to $\Delta t$, for the fixed mesh $M_4$

Figure 9: Exact $u$ at $t = t_3$ on $M_4$

Figure 11: Calculated $u$ at $t = t_3$ on $M_4$

Figure 13: Calculated $u$ at $t = t_3$ on $M_2$

Figure 15: Calculated $u$ at $t = t_3$ on $M_1$

Figure 17: Calculated $u$ at $t = t_1$ on $M_2$

Figure 19: Calculated $u$ at $t = t_2$ on $M_2$

Figure 21: Calculated $u$ at $t = t_3$ on $M_2$

Figure 23: Calculated $u$ at $t = t_4$ on $M_2$

Figure 25: Calculated $u$ at $t = t_5$ on $M_2$
<table><caption>Table 2: Exponent of $h$ for the error on $H$ profile</caption><thead><tr><th></th><th>t<sub>1</sub></th><th>t<sub>2</sub></th><th>t<sub>3</sub></th><th>t<sub>4</sub></th><th>t<sub>5</sub></th></tr></thead><tbody><tr><td>M<sub>1</sub> → M<sub>2</sub></td><td>1.905</td><td>1.905</td><td>1.982</td><td>2.112</td><td>2.141</td></tr><tr><td>M<sub>2</sub> → M<sub>4</sub></td><td>1.975</td><td>1.944</td><td>1.866</td><td>1.798</td><td>1.795</td></tr></tbody></table>

<table><caption>Table 3: Exponent of $h$ for the error on $u$ profile</caption><thead><tr><th></th><th>t<sub>1</sub></th><th>t<sub>2</sub></th><th>t<sub>3</sub></th><th>t<sub>4</sub></th><th>t<sub>5</sub></th></tr></thead><tbody><tr><td>M<sub>1</sub> → M<sub>2</sub></td><td>2.148</td><td>2.029</td><td>1.814</td><td>1.760</td><td>1.758</td></tr><tr><td>M<sub>2</sub> → M<sub>4</sub></td><td>1.043</td><td>1.043</td><td>1.100</td><td>1.194</td><td>1.213</td></tr></tbody></table>

<table><caption>Table 4: Exponent of $Δt$ for the error on $H$ profile, for the fixed mesh $M_2$</caption><thead><tr><th></th><th>t<sub>1</sub></th><th>t<sub>2</sub></th><th>t<sub>3</sub></th><th>t<sub>4</sub></th><th>t<sub>5</sub></th></tr></thead><tbody><tr><td>Δt<sub>1</sub> → Δt<sub>2</sub></td><td>0.594</td><td>0.544</td><td>0.482</td><td>0.469</td><td>0.465</td></tr><tr><td>Δt<sub>2</sub> → Δt<sub>4</sub></td><td>0.116</td><td>0.115</td><td>0.160</td><td>0.206</td><td>0.201</td></tr><tr><td>Δt<sub>4</sub> → Δt<sub>8</sub></td><td>-0.041</td><td>-0.026</td><td>0.028</td><td>0.073</td><td>0.071</td></tr></tbody></table>

<table><caption>Table 5: Exponent of $Δt$ for the error on $u$ profile, for the fixed mesh $M_2$</caption><thead><tr><th></th><th>t<sub>1</sub></th><th>t<sub>2</sub></th><th>t<sub>3</sub></th><th>t<sub>4</sub></th><th>t<sub>5</sub></th></tr></thead><tbody><tr><td>Δt<sub>1</sub> → Δt<sub>2</sub></td><td>0.988</td><td>0.960</td><td>0.904</td><td>0.830</td><td>0.805</td></tr><tr><td>Δt<sub>2</sub> → Δt<sub>4</sub></td><td>0.809</td><td>0.777</td><td>0.705</td><td>0.615</td><td>0.591</td></tr><tr><td>Δt<sub>4</sub> → Δt<sub>8</sub></td><td>0.391</td><td>0.398</td><td>0.380</td><td>0.340</td><td>0.333</td></tr></tbody></table>

<table><caption>Table 6: Exponent of $Δt$ for the error on $H$ profile, for the fixed mesh $M_4$</caption><thead><tr><th></th><th>t<sub>1</sub></th><th>t<sub>2</sub></th><th>t<sub>3</sub></th><th>t<sub>4</sub></th><th>t<sub>5</sub></th></tr></thead><tbody><tr><td>Δt<sub>1</sub> → Δt<sub>2</sub></td><td>1.058</td><td>1.044</td><td>0.980</td><td>0.923</td><td>0.936</td></tr><tr><td>Δt<sub>2</sub> → Δt<sub>4</sub></td><td>0.980</td><td>0.939</td><td>0.823</td><td>0.757</td><td>0.777</td></tr><tr><td>Δt<sub>4</sub> → Δt<sub>8</sub></td><td>0.567</td><td>0.507</td><td>0.453</td><td>0.450</td><td>0.459</td></tr></tbody></table>

<table><caption>Table 7: Exponent of $Δt$ for the error on $u$ profile, for the fixed mesh $M_4$</caption><thead><tr><th></th><th>t<sub>1</sub></th><th>t<sub>2</sub></th><th>t<sub>3</sub></th><th>t<sub>4</sub></th><th>t<sub>5</sub></th></tr></thead><tbody><tr><td>Δt<sub>1</sub> → Δt<sub>2</sub></td><td>1.012</td><td>1.005</td><td>0.993</td><td>0.976</td><td>0.968</td></tr><tr><td>Δt<sub>2</sub> → Δt<sub>4</sub></td><td>1.020</td><td>1.005</td><td>0.976</td><td>0.936</td><td>0.920</td></tr><tr><td>Δt<sub>4</sub> → Δt<sub>8</sub></td><td>0.995</td><td>0.964</td><td>0.905</td><td>0.827</td><td>0.804</td></tr></tbody></table>