Christine Bernardi
email: bernardi@ann.jussieu.fr

Ajmia Younes Orfi
email: ayounesorfi@yahoo.com

Finite element discretization of the time dependent axisymmetric Darcy problem

Keywords: Darcy's equations, axisymmetric domain, time discretization, finite element discretization, a posteriori analysis

We consider the Darcy problem in an axisymmetric three-dimensional domain with data which are axisymmetric. The solution satisfies a system of equations in the meridian domain. We propose a discretization of this problem in the case of an axisymmetric solution. This discretization relies on a backward Euler's scheme for the time variable and finite elements for the space variables. We prove a priori error estimates and a posteriori error estimates both for the time steps and the meshes and we present some numerical experiments which are in good agreement with the analysis.

Introduction

Let Ω be a bounded three-dimensional domain which is invariant by rotation around an axis. The boundary Γ of this domain is divided into two parts Γp and Γu . We are interested in the following model, suggested by Rajagobal [START_REF]On a hierarchy of approximate models for flows of incompressible fluids through porous solids[END_REF],

                 ∂ t ȗ + αȗ + gradp = f in Ω×]0, T [, divȗ = 0 in Ω×]0, T [, p = pb on Γp ×]0, T [, ȗ.n = g on Γu ×]0, T [, ȗ = ȗ0
in Ω at t = 0.

(

) 1
where the unknowns are the velocity ȗ and the pressure p of the fluid. The data are the quantities f, g, the pressure on the boundary pb and the initial values of the velocity ȗ0 . The parameter α is a positive constant representing the drag coefficient. If the problem is set in a domain which is symmetric by rotation around an axis, it is proved in [START_REF] Bernardi | Spectral Methods for Axisymmetric Domains[END_REF] that, when using a Fourier expansion with respect to the angular variable, a three-dimensional problem is equivalent to a system of two-dimensional problems on the meridian domain, each problem being satisfied by Fourrier coefficient of the solution. Here we are going to present an axisymmetric model, and we propose a discretization of this problem in the case of an axisymmetric solution, i.e. only for the Fourier coefficient of order 0. We recall that the problem considered in [START_REF] Bernardi | Discretisation of an unsteady flow through a porous solid modeled by Darcy's equations[END_REF] which is similar to the present problem is restricted to a boundary condition for the pressure and the case where the domain is a general two-or three-dimensional with a Lipschitz-continuous boundary. In the present study, mixed boundary conditions are considered, and we treat the problem in the simpler case of the meridian domain Ω where the data are axisymmetric, in the sense of [START_REF] Bernardi | Spectral Methods for Axisymmetric Domains[END_REF]Sec. II.3]. So, by using cylindrical coordinates, we can write a variational formulation of this problem in the meridian domain. We prove the well-posedness and some regularity properties of an axisymmetric solution for such a system. Next, we propose a time semidiscrete problem that relies on the backward Euler's scheme. We prove that this problem has a unique solution and derive error estimates. Concerning the space discretization, we consider a conforming finite element method which leads to a well-posed discrete problem for which we prove a priori error estimates. We introduce two families of error indicators, one for the time semi-discretization and another one for the space discretization. We also prove a posteriori error estimates which are optimal according to the standard criteria, see [START_REF] Verfürth | A Review of A Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques[END_REF]. Therefore the error indicators that we propose seem appropriate to perform time and space adaptivity in an efficient way. In a final step, we propose an algorithm for solving the resulting system and present some numerical experiments.

An outline of the paper is as follows: • In Sectin 2, we write a variational formulation of problem [START_REF] Achdou | A priori and a posteriori analysis of finite volume discretizations of Darcy's equations[END_REF] in the case of an axisymmetric solution, and we prove its well-posedness.

• Sectin 3 is devoted to the description and a priori analysis of the discrete problem in the meridian domain Ω.

• In Sectin 4, two families of error indicators are proposed and the a posteriori analysis of the discrete problem is carried out.

• In Sectin 5, we present some numerical experiments.

The two-dimensional formulation

We are interested in modeling a flow through a bounded and symmetric domain Ω with respect to the z axis. We use cylindrical coordinates (r, θ, z). So the axisymmetric domain Ω is the three-dimensional set obtained by rotating the two-dimentional domain Ω, called meridian domain, around the axis r = 0. The domain Ω is defined as:

Ω = {(r, θ, z); (r, z) ∈ Ω ∪ Γ 0 and θ ∈] -π, π]} ,
where Γ 0 is the intersection of Ω with the axis r = 0. For simplicity, we assume that Γ 0 is the union of a finite number of segments with positive measure. The fluid is modeled by the axisymmetric Darcy equations [START_REF] Achdou | A priori and a posteriori analysis of finite volume discretizations of Darcy's equations[END_REF] in the domain Ω, we suppose that the boundary conditions and the external forces are axisymmetric and that their angular component is zero. So, we are interested with problem in the special case of an axisymmetric geometry and for axisymmetric data. The two-dimensional axisymmetric boundary Γ of the physical domain Ω is a Lipschitz-continuous boundary and is divided into two parts Γp and Γu , also with Lipschitz continous boundaries. The part of the boundary Γp has a positive surface measure. Γu = Γ \ Γp is the union of a finite number of surface elements. Setting Γ = ∂Ω \ Γ 0 and rotating Γ around the axis r = 0 gives back Γ, and Γ 0 is a kind of artificial boundary. We also introduce the two parts Γ p and Γ u = Γ \ Γ p of the boundary of Γ. The unit outward normal vector n on Γ is obtained by rotating the unit outward vector n on Γ. An axisymmetric function p on Ω depends only on the radial and axial coordinates, therefore we associate a function p on Ω such that p(r, z) = p(r, 0, z). An axisymmetric vector field ȗ on Ω depends on (r, z). For any vector field ȗ, we denote by ȗr , ȗθ , ȗz its radial, angular and axial components, which are functions of r and z, therefore we associate a vector field u = (u r , u θ , u z ) on Ω such that u r = ȗr , u θ = ȗθ , and u z = ȗz . Relying on the isomorphisms proved in [6, Chap. II], we can write and analyze the variational formulation of the reduced two-dimensional problems. We recall that f r , f θ , f z denote the cylindrical components of f, which are independent of θ. p b and g are also independent of θ.

We introduce the following operators

gradp =   ∂ r p 0 ∂ z p   and divu = ∂ r u r + 1 r u r + ∂ z u z . Then problem (1) is equivalent to                          ∂ t u r + α u r + ∂ r p = f r in Ω×]0, T [, ∂ t u θ + α u θ = f θ in Ω×]0, T [, ∂ t u z + α u z + ∂ z p = f z in Ω×]0, T [, ∂ r u r + r -1 u r + ∂ z u z = 0 in Ω×]0, T [, p = p b on Γ p ×]0, T [, u r n r + u z n z = g on Γ u ×]0, T [, (u r , u θ , u z ) = (u 0r , u 0θ , u 0z ) in Ω at time t = 0. (2)
This problem reduces to a system of two uncoupled problems that we treat separately: the simplest one is a scalar coercive equation for the angular velocity u θ and the other is a saddle-point type problem for (u r , u z , p). ) ,

where u θ belongs to H 1 ( 0, T ; L 2 1 (Ω)

) . So from now on we only consider the reduced problem

                   ∂ t u r + α u r + ∂ r p = f r in Ω×]0, T [, ∂ t u z + α u z + ∂ z p = f z in Ω×]0, T [, ∂ r u r + r -1 u r + ∂ z u z = 0 in Ω×]0, T [, p = p b on Γ p ×]0, T [, u r n r + u z n z = g on Γ u ×]0, T [, (u r , u z ) = (u 0r , u 0z )
in Ω at time t = 0.

(

) 3
In order to write the variational formulation of problem (3) and according to [6, Sec. II.2], we consider the weighted Sobolev spaces

L 2 1 (Ω) = { v : Ω → R measurable; ∫ Ω |v(r, z)| 2 r dr dz < +∞ } , equipped with the norm ∥v∥ L 2 1 (Ω) = (∫ Ω |v(r, z)| 2 r dr dz ) 1 2 ,
and also the space

H 1 1 (Ω) = { v ∈ L 2 1 (Ω); ∂ r v ∈ L 2 1 (Ω) et ∂ z v ∈ L 2 1 (Ω) } ,
which is provided with the seminorm and norm

|v| H 1 1 (Ω) = (∥∂ z v∥ 2 L 2 1 (Ω) + ∥∂ r v∥ 2 L 2 1 (Ω) ) 1 2 , ∥v∥ H 1 1 (Ω) = (∥v∥ 2 L 2 1 (Ω) + |v| 2 H 1 1 (Ω) ) 1 2 .
The trace on Γ u is defined in a nearly standard way see [START_REF] Bernardi | Polynomials in weighted Sobolev spaces: basics and trace liftings[END_REF]Sec. 2]. If H s 1 (Γ u ), s ≥ 0, stands for the scale of Sobolev spaces built from

L 2 1 (Γ u ) = { g : Γ u → R mesurable; ∫ Γu g 2 (τ ) r(τ ) dτ < +∞ } ,
(where r(τ ) denotes the distance of the point with tangential coordinate τ to the axis r = 0), the trace operator:

v -→ v |Γu is continuous from H 1 1 (Ω) onto H 1 2 1 (Γ u ) see [6, Chap. II]. Let also H 1 2 1 (Γ p ) the space of traces of functions in H 1 1 (Ω) on Γ p . The variational space H 1 1⋄ (Ω) = { q ∈ H 1 1 (Ω); q = 0 sur Γ p } ,
is a Hilbert space for the scalar product associated with the previously defined norm.

We denote by (., .) 1 the scalar product on L 2 1 (Ω) 2 . The variational formulation of the saddlepoint type problem is written:

Find (u = (u r , u z ), p) in H 1 (0, T ; L 2 1 (Ω) 2 ) × L 2 (0, T ; H 1 1 (Ω)) such that u(., 0) = u 0 in Ω, (4)
for a.e. t,

0 ≤ t ≤ T, p(., t) = p b on Γ p , ( 5
)
∀v ∈ L 2 1 (Ω) 2 , a(∂ t u, v) + α a(u, v) + b(v, p) = ∫ Ω f(r, z).v(r, z) r dr dz, ∀q ∈ H 1 1⋄ (Ω), b(u, q) = ∫ Γu g(τ ) q(τ ) r(τ ) dτ, ( 6
)
where the bilinear forms a(., .) and b(., .) are defined by:

a(u, v) = (u, v) 1 = ∫ Ω (u r (r, z).v r (r, z) + u z (r, z).v z (r, z)) r dr dz, b(v, q) = (v, gradq) 1 = ∫ Ω (v r (r, z).∂ r q(r, z) + v z (r, z).∂ z q(r, z)) r dr dz.
It is readily checked that the forms a(., .) and b(., .) are continuous on

L 2 1 (Ω) 2 × L 2 1 (Ω) 2 and L 2 1 (Ω) 2 × H 1 1 (Ω) respectively.
The kernel of the bilinear form b(., .) is

V(Ω) = { v ∈ L 2 1 (Ω) 2 ; ∀q ∈ H 1 1⋄ (Ω), b(v, q) = 0 } , is characterized by V(Ω) = { v ∈ L 2 1 (Ω) 2 ; div v = 0 and v.n = 0 on Γ u } ,
and its orthogonal in L 2 1 (Ω) 2 is defined by

V(Ω) ⊥ = { v ∈ L 2 1 (Ω) 2 ; ∀w ∈ V(Ω), ∫ Ω v . w r dr dz = 0 } .
Proving the well-posedness of problem (4-5-6) relies on the ellipticity of a(., .) and on an inf-sup condition of Babuška and Brezzi type on the form b(., .). We begin with this condition. Lemma 1. There exists a constant β > 0 such that the following inf-sup condition holds

∀q ∈ H 1 1⋄ (Ω), sup v=(vr,vz)∈L 2 1 (Ω) 2 b(v, q) ∥v∥ L 2 1 (Ω) 2 ≥ β |q| H 1 1 (Ω) . ( 7
)
Proof. Let q be any function in H 1 1⋄ (Ω). The idea is to choose v equal to gradq, so that

b(v, q) = ∥∂ r q∥ 2 L 2 1 (Ω) + ∥∂ z q∥ 2 L 2 1 (Ω) = |q| 2 H 1 1 (Ω) , and ∥v∥ L 2 1 (Ω) 2 = |q| H 1 1 (Ω)
. This gives the desired inf-sup condition. To make this condition complete, the weighted Poincaré-Friedrichs condition ensures the equivalence between the norms

|.| H 1 1 (Ω) and ∥.∥ H 1 1 (Ω) on H 1 1⋄ (Ω), see [2, Lem. 3.1].

Lemma 2. The following ellipticity property holds

∀v = (v r , v z ) ∈ L 2 1 (Ω) 2 , a(v, v) ≥ ∥v∥ 2 L 2 1 (Ω) 2 . ( 8
)
Proof. We have

a(v, v) = ∫ Ω |v(r, z)| 2 r dr dz = ∥v∥ 2 L 2 1 (Ω) 2 ,
which implies the desired inequality.

Thanks to Lemmas 1 and 2, we easily derive the next theorem, see [START_REF] Bernardi | Spectral Methods for Axisymmetric Domains[END_REF].

Theorem 1. For any data

(f, p b , g) ∈ L 2 (0, T ; L 2 1 (Ω) 2 ) × L 2 (0, T ; H 1 2 1 (Γ p )) × L 2 (0, T ; L 2 1 (Γ u ))
and u 0 ∈ L 2 1 (Ω) 2 , problem (4-5-6) has a unique solution

(u = (u r , u z ), p) ∈ H 1 (0, T ; L 2 1 (Ω) 2 ) × L 2 (0, T ; H 1 1 (Ω)).
Moreover the 4-tuple (u r , u θ , u z , p) is equal to (R -θ ȗ, p), where (ȗ, p) is the solution of problem [START_REF] Achdou | A priori and a posteriori analysis of finite volume discretizations of Darcy's equations[END_REF] with axisymmetic data and R θ the rotation with angle θ with respect to the axis r = 0 in R 3 .

Remark 1. The third equation in problem (3) can equivalently be written

∂ r (r u r ) + ∂ z (r u z ) = 0.
Hence, the function (r u r , r u z ) is divergence-free in the standard sense, so it is the curl of a function φ. Setting: φ = r ψ, we derive the existence of a scalar potential ψ such that

u r = ∂ z ψ and u z = - 1 r ∂ r (r ψ) on Ω.
Now we will find a new simpler variationel formulation which is equivalent to problem (4-5-6). For this let L denote a lifting operator, which is continuous from

H 1 2 1 (Γ p ) into H 1 1 (Ω)
, the existence of such operateur is established in [START_REF] Jerisson | The inhomogeneous Dirichlet problem in Lipschitz domains[END_REF], for instance. Since p b ∈ L 2 (0, T ; H 1 2 1 (Γ p )) we denote by pb the function defined for a.e. t, 0 ≤ t ≤ T,

pb (t) = L (p b (t)). (9)
The function pb belongs to L 2 (0, T ; H 1 1 (Ω)) and satisfies

∥p b ∥ L 2 (0,T ;H 1 1 (Ω)) ≤ c 0 ∥p b ∥ L 2 (0,T ;H 1 2 1 (Γp)) . ( 10
)
The last equation in problem (4-5-6) can be viewed as a non-homogeneous constraint; let us show that we can find a function of L 2 1 (Ω) 2 that satisfies this constraint. For this, define the linear operator B :

L 2 1 (Ω) 2 → H 1 1⋄ (Ω) ′ by ⟨Bv, q⟩ H 1 1⋄ (Ω) ′ ×H 1 1⋄ (Ω) = b(v, q). ( 11
)
It follows from inf-sup condition [START_REF] Bernardi | Discretisation of an unsteady flow through a porous solid modeled by Darcy's equations[END_REF] and [10, Chap. I, Lem. 4.1], that this operator is an isomorphism from

V(Ω) ⊥ onto H 1 1⋄ (Ω) ′ and ∀v ∈ V(Ω) ⊥ , ∥Bv∥ H 1 1⋄ (Ω) ′ ≥ β ∥v∥ L 2 1 (Ω) 2 . Hence, for g ∈ H 1 1⋄ (Ω) ′ there exists a unique u b ∈ V(Ω) ⊥ such that ∀q ∈ H 1 1⋄ (Ω), b(u b , q) = ⟨g, q⟩ Γu ,
and

β ∥u b (., t)∥ L 2 1 (Ω) 2 ≤ ∥g(., t)∥ L 2 1 (Γu) . ( 12
)
When setting u ⋄ = uu b and p * = ppb , we observe that

u ⋄ ∈ H 1 (0, T ; V(Ω)) and p * ∈ L 2 (0, T ; H 1 1⋄ (Ω)). Moreover, if (u, p) is a solution of (4-5-6) then u ⋄ ∈ H 1 (0, T ; V(Ω))
is the unique solution of the simpler variational problem:

Find u ⋄ = (u r⋄ , u z⋄ ) in H 1 (0, T ; V(Ω)) such that, for a.e. t, 0 ≤ t ≤ T, { u ⋄ (., 0) = u 0 -u b0 = u ⋄0 in Ω, ∀v ∈ V(Ω); (∂ t u ⋄ , v) 1 + α (u ⋄ , v) 1 = (f, v) 1 -(v, gradp b ) 1 , (13)
where

u b0 = u b (., 0) in Ω.
Conversely, let us prove that, if u ⋄ is a solution of ( 13), then there exists a unique pair (4-5-6). For this, we integrate the second equation in ( 13) between 0 and t, we define the functional for all v ∈ L 2 1 (Ω) 2 :

(u = u ⋄ + u b , p = p * + pb ) in H 1 (0, T ; L 2 1 (Ω) 2 ) × L 2 (0, T ; H 1 1 (Ω)) solution of problem
L t (v) = ∫ t 0 ((f(., s), v) 1 -α(u ⋄ (., s), v) 1 -b(v, pb (s)))ds-(u ⋄ (., t), v) 1 +(u ⋄0 , v) 1 .
For all t ∈ [0, T ], L t is a continuous linear functional on L 2 1 (Ω) 2 and, according to [START_REF] Johnson | An a posteriori error estimate and adaptive timestep control for a backward Euler discretization of a parabolic problem[END_REF], L t (v) = 0 ∀v ∈ V(Ω). Hence, see [9, Chap. V, Thm. 1.3], for each t ∈ [0, T ] , there exists a unique function P (t) in H 1 1⋄ (Ω) such that: [START_REF]On a hierarchy of approximate models for flows of incompressible fluids through porous solids[END_REF] with respect to t, and setting p * (t) = d dt P (t), we obtain (6) with p = p * + pb . This gives immediately [START_REF] Bernardi | Polynomials in weighted Sobolev spaces: basics and trace liftings[END_REF]. Then we can conclude the following corollary.

∀v ∈ L 2 1 (Ω) 2 , L t (v) = b(v, P (t)), (14)
|P (t)| H 1 1 (Ω) ≤ sup v∈L 2 1 (Ω) 2 L t (v) ∥v∥ L 2 1 (Ω) 2 . Now, differentiating
Corollary 1. The variational problems (4-5-6) and ( 13) are equivalent.

Lemma 3. For any data

(f, p b , g) ∈ L 2 (0, T ; L 2 1 (Ω) 2 ) × L 2 (0, T ; H 1 2 1 (Γ p )) × L 2 (0, T ; L 2 1 (Γ u ))
and u 0 ∈ L 2 1 (Ω) 2 , we have the following a priori estimates for the velocity u solution of problem (4-5-6), for a.e.

t ∈ [0, T ] , ∥u∥ L ∞ (0,t;L 2 1 (Ω) 2 ) ≤ c(∥u 0 ∥ L 2 1 (Ω) 2 + ∥g(., 0)∥ L 2 1 (Γu) + ∥f∥ L 2 (0,t;L 2 1 (Ω) 2 ) (15)
+ ∥p b ∥ L 2 (0,t;H 1 2 1 (Γp)) + ∥g∥ L ∞ (0,t;L 2 1 (Γu))
), with a constant c that only depends on Ω and T .

Proof. Taking v = u ⋄ in (13) gives (∂ t u ⋄ , u ⋄ ) 1 + α(u ⋄ , u ⋄ ) 1 = (f, u ⋄ ) 1 -(u ⋄ , gradp b ) 1 .
By using Cauchy-Schwarz inequality yields

1 2 d dt ∥u ⋄ ∥ 2 L 2 1 (Ω) 2 + α ∥u ⋄ ∥ 2 L 2 1 (Ω) 2 ≤ ∥f∥ L 2 1 (Ω) 2 . ∥u ⋄ ∥ L 2 1 (Ω) 2 + |p b | H 1 1 (Ω) . ∥u ⋄ ∥ L 2 1 (Ω) 2 . Using Young's inequality: a b ≤ a 2 2α + α 2 b 2
, we obtain:

d dt ∥u ⋄ ∥ 2 L 2 1 (Ω) 2 ≤ 1 α (∥f∥ 2 L 2 1 (Ω) 2 + |p b | 2 H 1 1 (Ω)
). Integrating this inequality between 0 and t, using the fact that u

⋄ = u -u b in Ω and u ⋄ (., 0) = u 0 -u b0 yields ∥u(., t) -u b (., t)∥ 2 L 2 1 (Ω) 2 ≤ ∥u 0 -u b0 ∥ 2 L 2 1 (Ω) 2 + 1 α (∥f(., t)∥ 2 L 2 (0,t;L 2 1 (Ω) 2 ) + |p b (., t)| 2 L 2 (0,t;H 1 1 (Ω))
). By triangle inequality, estimates [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations. Theory and Algorithms[END_REF] and ( 12) we obtain

∥u(., t)∥ 2 L 2 1 (Ω) 2 ≤ c(∥u 0 ∥ L 2 1 (Ω) 2 + ∥g(., 0)∥ L 2 1 (Γu) + ∥f∥ L 2 (0,t;L 2 1 (Ω) 2 ) (16)
+ ∥p b ∥ L 2 (0,t;H 1 2 1 (Γp)) + ∥g(., t)∥ L 2 1 (Γu)
). This gives the desired estimate [START_REF] Scott | Finite element interpolation of non-smooth functions satisfying boundary conditions[END_REF].

We refer to [7, Thm. 2.4], for the detailed proof of the next Theorem.

Theorem 2. For any data

(f, p b , g) ∈ L 2 (0, T ; L 2 1 (Ω) 2 ) × L 2 (0, T ; H 1 2 1 (Γ p )) × H 1 (0, T ; L 2 1 (Γ u )) and u 0 ∈ L 2 1 (Ω) 2
, the unique solution

(u = (u r , u z ), p) ∈ H 1 (0, T ; L 2 1 (Ω) 2 ) × L 2 (0, T ; H 1 1 (Ω))
of problem (4-5-6), satisfies the a priori estimate for a.e. t ∈ [0, T ] ,

∥u∥ H 1 (0,T ;L 2 1 (Ω) 2 ) + ∥p∥ L 2 (0,T ;H 1 1 (Ω)) (17) ≤ c(∥u 0 ∥ L 2 1 (Ω) 2 + ∥f∥ L 2 (0,T ;L 2 1 (Ω) 2 ) + ∥p b ∥ L 2 (0,t;H 1 2 1 (Γp)) + ∥g∥ H 1 (0,t;L 2 1 (Γu)) ).

The discrete problem and its a priori analysis

We split the discretization into two steps: First a semi-discretization in time, and next the full discretization. At each step, we prove a priori error estimates.

The time semi-discrete problem

We introduce a partition of the interval 2 , satisfaying divu 0 = 0 in Ω, we consider the following scheme:

[0, T ] into subintervals [t n-1 , t n ], 1 ≤ n ≤ N , such that 0 = t 0 < t 1 < ... < t N = T.
) ∈ C 0 ( 0, T ; L 2 1 (Ω) 2 ) × C 0 ( 0, T ; H 1 2 1 (Γ p ) ) , g ∈ C 0 (0, T ; L 2 1 (Γ u )) and u 0 ∈ L 2 1 (Ω)
Find (u n ) 0≤n≤N ∈ ( L 2 1 (Ω) 2 ) N +1 and (p n ) 1≤n≤N ∈ (H 1 1 (Ω)) N such that u 0 = u 0 in Ω, (18) ∀n, 1 ≤ n ≤ N p n = p n b on Γ p , ( 19
) ∀v ∈ L 2 1 (Ω) 2 , (u n , v) 1 + α τ n (u n , v) 1 = (u n-1 , v) 1 -τ n (v, gradp n ) 1 + τ n (f n , v) 1 , ∀q ∈ H 1 1⋄ (Ω), (u n , gradq) 1 = ⟨g n , q⟩ Γu , ( 20
)
where

f n = f(., t n ), g n = g(., t n ) and p n b = p b (., t n ).
Now we will find a new simpler variational formulation which is equivalent to problem (18-19-20). For this we use the lifting operator L introduced in (9), verifying

∥p n b ∥ H 1 1 (Ω) ≤ c 0 ∥p n b ∥ H 1 2 1 (ΓP ) , ( 21
)
and the linear operator B introduced in (11) which is an isomorphism from

V(Ω) ⊥ into H 1 1⋄ (Ω) ′ . Hence, there exists a unique u n b ∈ V(Ω) ⊥ such that ∀q ∈ H 1 1⋄ (Ω), b(u n b , q) = ⟨g n , q⟩ Γu ,
and

β ∥u n b (., t)∥ L 2 1 (Ω) 2 ≤ ∥g n (., t)∥ L 2 1 (Γu) . ( 22
)
When setting u n ⋄ = u nu n b , then u n ⋄ belongs to V(Ω), and the sequence

(u n ⋄ ) 0≤n≤N ∈ (V(Ω)) N +1
is a solution of the simpler variational problem:

Find (u n ⋄ ) 0≤n≤N in (V(Ω)) N +1 such that for all n, 1 ≤ n ≤ N,    u 0 ⋄ = u 0 -u 0 b in Ω, ∀v ∈ V(Ω), (u n ⋄ , v) 1 + ατ n (u n ⋄ , v) 1 = (u n-1 ⋄ , v) 1 + τ n (f n , v) 1 -τ n b(v, pn b ). (23)
Conversely, if (u n ⋄ ) 0≤n≤N is a solution of (23), we define the linear mapping

L n (v) = ∫ Ω (f n .v)(r, z)rdrdz -b(v, pn b ) - 1 τ n (u n ⋄ -u n-1 ⋄ , v) 1 -α(u n ⋄ , v) 1 .
The mapping: 2 , and ∀v ∈ V(Ω) L n (v) = 0, thus the inf-sup condition [START_REF] Bernardi | Discretisation of an unsteady flow through a porous solid modeled by Darcy's equations[END_REF] implies that there exists a unique element p n * (t) in H 1 1⋄ (Ω) such that:

v → L n (v) is a continuous linear functional on L 2 1 (Ω)
∀v ∈ L 2 1 (Ω) 2 , b(v, p n * ) = L n (v), |p n * | H 1 1 (Ω) ≤ sup v∈L 2 1 (Ω) 2 L n (v) ∥v∥ L 2 1 (Ω) 2 . Therefore, (u n = u n ⋄ + u n b , p n = p n * +
) ∈ C 0 ( 0, T ; L 2 1 (Ω) 2 ) × C 0 ( 0, T ; H 1 2 1 (Γ p ) ) , g ∈ L 2 (0, T ; L 2 1 (Γ u )) and u 0 ∈ L 2 1 (Ω) 2 , satisfaying divu 0 = 0 in Ω, problem (18-19-20) has a unique solution (u n , p n ) such that: ∀ n, 0 ≤ n ≤ N, u n ∈ L 2 1 (Ω) 2 and ∀ n, 1 ≤ n ≤ N, p n ∈ H 1 1 (Ω)
. Moreover the sequence of velocities (u n ) 0≤n≤N satisfies:

∥u n ∥ L 2 1 (Ω) 2 ≤ ∥u 0 ∥ L 2 1 (Ω) 2 + √ 2 α ( n ∑ m=1 τ m (∥f m ∥ 2 L 2 1 (Ω) 2 + c 2 0 ∥p m b ∥ 2 H 1 2 1 (Γp)
))

1 2

(24)

+ ∥g n ∥ L 2 1 (Γu) + ∥g(., 0)∥ L 2 1 (Γu) , ( n ∑ m=1 τ m u m -u m-1 τ m 2 L 2 1 (Ω) 2 ) 1 2 ≤ √ 2α(∥u 0 ∥ L 2 1 (Ω) 2 + ∥g(., 0)∥ L 2 1 (Γu) ) (25) +2 ( n ∑ m=1 τ m (∥f m ∥ 2 L 2 1 (Ω) 2 + c 2 0 ∥p m b ∥ 2 H 1 2 1 (Γp)
)

) 1 2 + √ 2 β ( n ∑ m=1 1 τ m g m -g m-1 2 L 2 1 (Γu) ) 1 2 .
And the sequence of pressures (p n ) 1≤n≤N satisfies:

( n ∑ m=1 τ m |p m | 2 H 1 1 (Ω) ) 1 2 ≤ c(∥u 0 ∥ 2 L 2 1 (Ω) 2 + n ∑ m=1 τ m (∥f m ∥ 2 L 2 1 (Ω) 2 + ∥p m b ∥ 2 H 1 2 1 (Γp) ) + ∥g(., 0)∥ L 2 1 (Γu) ) + n ∑ m=1 1 τ m g m -g m-1 2 L 2 1 (Γu) ) 1 2 . ( 26
)
Proof. Clearly, problem (23) has a unique solution:

u n ⋄ = u n -u n b in V(Ω)
, which yields the existence and uniqueness of the solution

(u n = u n ⋄ + u n b , p n = p n * + pn b ) of problem (18-19-20). 1) To prove estimate (24), we choose v = u n ⋄ in (23), this gives ∥u n ⋄ ∥ 2 L 2 1 (Ω) 2 -(u n-1 ⋄ , u n ⋄ ) 1 + ατ n ∥u n ⋄ ∥ 2 L 2 1 (Ω) 2 = τ n (f n , u n ⋄ ) 1 -τ n (u n ⋄ , gradp n b ) 1 . Using the fact that (a 2 -ab) = 1 2 (a 2 -b 2 + (a -b) 2
) and the Cauchy-Schwarz inequality we obtain

1 2 ( ∥u n ⋄ ∥ 2 L 2 1 (Ω) 2 -u n-1 ⋄ 2 L 2 1 (Ω) 2 + u n ⋄ -u n-1 ⋄ 2 L 2 1 (Ω) 2 ) + ατ n ∥u n ⋄ ∥ 2 L 2 1 (Ω) 2 ≤ τ n ( ∥f n ∥ L 2 1 (Ω) 2 + |p n b | H 1 1 (Ω) ) ∥u n ⋄ ∥ L 2 1 (Ω) 2 . Therefore, Young's inequality gives ∥u n ⋄ ∥ 2 L 2 1 (Ω) 2 -u n-1 ⋄ 2 L 2 1 (Ω) 2 + u n ⋄ -u n-1 ⋄ 2 L 2 1 (Ω) 2 + ατ n ∥u n ⋄ ∥ 2 L 2 1 (Ω) 2 ≤ 2τn α ( ∥f n ∥ 2 L 2 1 (Ω) 2 + |p n b | 2 H 1 1 (Ω)
) . Summing this inequality over n, 1 ≤ n ≤ N we obtain

∥u n ⋄ ∥ 2 L 2 1 (Ω) 2 + n ∑ m=1 u m ⋄ -u m-1 ⋄ 2 L 2 1 (Ω) 2 + α n ∑ m=1 τ m ∥u m ⋄ ∥ 2 L 2 1 (Ω) 2 ≤ u 0 ⋄ 2 L 2 1 (Ω) 2 + 2 α n ∑ m=1 τ m ( ∥f m ∥ 2 L 2 1 (Ω) 2 + |p m b | 2 H 1 1 (Ω)
) , then,

∥u n ⋄ ∥ L 2 1 (Ω) 2 ≤ u 0 ⋄ L 2 1 (Ω) 2 + √ 2 α ( n ∑ m=1 τ m ( ∥f m ∥ 2 L 2 1 (Ω) 2 + ∥p m b ∥ 2 H 1 1 (Ω) ) ) 1
2 .

Therefore by using the triangle inequality, the fact that u n ⋄ = u n -u n b , the initial condition of problem (23), the estimates (21) and ( 22), we obtain the desired estimate (24).

2)

In order to prove estimate (25), we take

v = (u n ⋄ -u n-1 ⋄ ) ∈ V(Ω) in problem (23), yields u n ⋄ -u n-1 ⋄ 2 L 2 1 (Ω) 2 + ατ n ( ∥u n ⋄ ∥ 2 L 2 1 (Ω) 2 -(u n ⋄ , u n-1 ⋄ ) 1 ) = τ n (f n , u n ⋄ -u n-1 ⋄ ) 1 -τ n (u n ⋄ -u n-1 ⋄ , gradp n b ) 1 .
By using the fact that (a

2 -ab) = 1 2 (a 2 -b 2 + (a -b) 2
) and Cauchy-Schwarz inequality we obtain

1 2 u n ⋄ -u n-1 ⋄ 2 L 2 1 (Ω) 2 + ατn 2 ( ∥u n ⋄ ∥ 2 L 2 1 (Ω) 2 -u n-1 ⋄ 2 L 2 1 (Ω) 2 + u n ⋄ -u n-1 ⋄ 2 L 2 1 (Ω) 2 ) ≤ τ 2 n ( ∥f n ∥ 2 L 2 1 (Ω) 2 + |p n b | 2 H 1 1 (Ω)
) .

Multiplying by 2 τn and summing over n, with

1 ≤ n ≤ N yields n ∑ m=1 τ m u m ⋄ -u m-1 ⋄ τ m 2 L 2 1 (Ω) 2 + α n ∑ m=1 u m ⋄ -u m-1 ⋄ 2 L 2 1 (Ω) 2 -α u 0 ⋄ 2 L 2 1 (Ω) 2 ≤ 2 n ∑ m=1 τ m ( ∥f m ∥ 2 L 2 1 (Ω) 2 + |p m b | 2 H 1 1 (Ω) ) , whence, n ∑ m=1 τ m u m ⋄ -u m-1 ⋄ τ m 2 L 2 1 (Ω) 2 ≤ α u 0 ⋄ 2 L 2 1 (Ω) 2 (27) +2 n ∑ m=1 τ m ( ∥f m ∥ 2 L 2 1 (Ω) 2 + ∥p m b ∥ 2 H 1 1 (Ω)
) .

On the other hand by using the triangle inequatity, the fact that u n ⋄ = u nu n b for any n, 0 ≤ n ≤ N , and estimate (22), we obtain

u n -u n-1 L 2 1 (Ω) 2 ≤ u n ⋄ -u n-1 ⋄ L 2 1 (Ω) 2 + 1 β g n -g n-1 L 2 1 (Γu)
. Dividing the square of this inequality by τ n , summing over n with 1 ≤ n ≤ N , using estimates (27), (21) and the initial condition of problem (23) we obtain (25).

3) To prove estimate (26), we choose v = gradp n in the first equation of (20) and apply Cauchy-Schwarz inequality we obtain

|p n | H 1 1 (Ω) ≤ u n -u n-1 τn L 2 1 (Ω) 2 + α ∥u n ∥ L 2 1 (Ω) 2 + ∥f n ∥ L 2 1 (Ω) 2 .
Multiplying the square of this inequality by τ n , summing over n, using estimates (24) and (25), we obtain (26).

Remark 2. Let Π τ denote the operator which associates with any continuous function

v ∈ [0, T ] the constant function Π τ v equal to v(t n ) on each interval ]t n-1 , t n ], 1 ≤ n ≤ N. Then, estimate (24) is equivalent to the following sup 0≤m≤n ∥u m ∥ L 2 1 (Ω) 2 ≤ c(∥u 0 ∥ L 2 1 (Ω) 2 + ∥Π τ f∥ L 2 (0,tn;L 2 1 (Ω) 2 )
(28)

+ ∥g(., 0)∥ L 2 1 (Γu) + ∥Π τ p b ∥ L 2 (0,tn;H 1 2 1 (Γp))
+ ∥Π τ g∥ L 2 (0,tn;L 2 1 (Γu)) ).

In order to state the a priori error estimate, we observe that the family (e n ) 0≤n≤N , with e n = u(., t n ) -u n satisfies e 0 = 0 and also, by integrating ∂ t u between t n-1 and t n and subtracting (20) from [START_REF] Bernardi | Spectral Methods for Axisymmetric Domains[END_REF]

at time t n ,    ∀v ∈ L 2 1 (Ω) 2 , (e n , v) 1 + ατ n (e n , v) 1 = (e n-1 , v) 1 + τ n (ϵ n , v) 1 -τ n b(v, p(., t n ) -p n ), ∀q ∈ H 1 1⋄ (Ω), b(e n , q) = 0. ( 29
)
Where the consistency error ϵ n is given by

ϵ n = 1 τ n ∫ tn tn-1 (∂ t u)(s) ds -(∂ t u)(t n ).
We assume that the velocity u of problem (4-5-6) belongs to H 2 ( 0, T ; L 2 1 (Ω) 2 ) , then we can conclude this section, by recalling the main results concerning the a priori estimates, which are proven in [7, Prop. 3.2 and Cor.

3.1], for n, 1 ≤ n ≤ N : (i) ∥e n ∥ L 2 1 (Ω) 2 ≤ 1 √ 3α |τ | ∥u∥ H 2 (0,tn;L 2 1 (Ω) 2 ) .
(ii)

( n ∑ m=1 τ m e m -e m-1 τ m 2 L 2 1 (Ω) 2 ) 1 2 ≤ 1 √ 3 |τ | ∥u∥ H 2 (0,tn;L 2 1 (Ω) 2 ) , ( n ∑ m=1 τ m |p(., t m ) -p m )| 2 H 1 1 (Ω) ) 1 2 ≤ 1 √ 3 |τ | ∥u∥ H 2 (0,tn;L 2 1 (Ω) 2 ) .

The time and space discrete problem

We now describe the space discretization of problem (18-19-20). For each n, 0 ≤ n ≤ N , let (T nh ) h be a regular family of triangulations of Ω by closed triangles, in the usual sense that • for each h, Ω is the union of all elements of T nh , • Γ p is the union of whole edges of elements of T nh ,

• there exists a constant σ > 0 independant of h, n and T such that, for all T in T nh , hT ρT ≤ σ, where h T is the diameter of T, and ρ T the diameter of its inscribed circle, • h n the maximum of the diameters of the elements of T nh , • E nh is the set of all edges e of elements T of T nh , • E 0 nh is the subset of E nh which elements are not contained in ∂Ω,

• V nh : is the set of vertices of the elements of T nh , • V 0 nh : is the subset of V nh which elements are inside Ω,

• V b nh = V nh \ V 0 nh : is the subset of V nh of boundary vertices.
For each triangle T and nonnegative integer k, we denote by P k (T ) the space of restrictions to T of polynomials with degree ≤ k. At each time step, the discrete space of velocities is:

X nh (Ω) = { v h ∈ L 2 1 (Ω) 2 / ∀ T ∈ T nh , v h T ∈ P 0 (T ) 2 } ,
its interpolation operator is the orthogonal projection operator Π nh : L 2 1 (Ω) 2 → X nh associated with the scalar product of L 2 1 (Ω) 2 and verify, for every 0 ≤ s ≤ 1

∀v ∈ H s 1 (Ω), ∥v -Π nh v∥ L 2 1 (Ω) 2 ≤ Ch s n ∥v∥ H s 1 (Ω) 2 . ( 30
)
We assume that the pressure is continuous whence the choice of discrete space as proposed in [START_REF] Achdou | A priori and a posteriori analysis of finite volume discretizations of Darcy's equations[END_REF]:

M nh (Ω) = { q h ∈ H 1 1 (Ω)/ ∀ T ∈ T nh , q h | T ∈ P 1 (T )} ,
its degrees of freedom are defined at the nodes of V nh and its interpolation operator i nh :

H 1 1 (Ω) → M nh (Ω)
is the standard Lagrange interpolation operator at the nodes of V nh with values in M nh and verify, for every

1 2 < s ≤ 1 ∀q ∈ H s+1 1 (Ω), |q -i nh q| H 1 1 (Ω) ≤ Ch s n ∥q∥ H s+1 1 (Ω) , (31)
finally to approximate functions with zero trace on Γ p , we set

M 0 nh (Ω) = {q h ∈ M nh (Ω)/ q h = 0 on Γ p } .

Variational formulation of the discrete problem

For every data (f, p b ) which belongs to C 0 (0, T ;

L 2 1 (Ω) 2 ) × C 0 (0, T ; H s+ 1 2 1 (Γ p )), s > 1 2
, g belongs to C 0 (0, T ; L 2 1 (Γ u )) and u 0 belongs to L 2 1 (Ω) 2 satisfies divu 0 = 0 in Ω, the discret problem constructed by the Galerkin method from (18)-( 19)-(20) reads:

Find (u n h ) 0≤n≤N ∈ (X nh (Ω)) N +1 and (p n h ) 1≤n≤N ∈ (M nh (Ω)) N such that u 0 h = Π 0h u 0 in Ω, (32) ∀n, 1 ≤ n ≤ N , p n h = i nh p b on Γ p , ( 33
)
∀v h ∈ X nh (Ω), (u n h , v h ) 1 + α τ n (u n h , v h ) 1 + τ n b(v h , p n h ) = (u n-1 h , v h ) 1 + τ n ∫ Ω (f n .v h )(r, z) rdrdz, ∀q h ∈ M 0 nh (Ω), b(u n h , q h ) = ⟨g n , q h ⟩ Γu . ( 34
)
The function v h = (gradq h ) | T belongs to X nh and we have the inf-sup condition

∀q h ∈ M nh (Ω), sup v h ∈X nh b(v h , q h ) ∥v h ∥ L 2 1 (Ω) 2 = |q h | H 1 1 (Ω) . ( 35
)
Let the discrete kernel

V nh (Ω) = { v h ∈ X nh (Ω); ∀q h ∈ M 0 nh (Ω), b(v h , q h ) = 0 } .
The choice of the lifting operator L, which is introduced in (9) yields that,

if p n b belongs to H s+ 1 2 1 (Γ p ), 1 2 < s ≤ 1, then pn b belongs to H s+1 1
(Ω), therefore its interpolate by Lagrange interpolation operator i nh (p n b ) is well-defined. Define the linear operator B : ,q). It follows from inf-sup condition (35) and [10, Chap. I, Lem. 4.1], that this operator is an isomorphism from

X nh → M 0 nh (Ω) ′ by ⟨B nh v, q⟩ M 0 nh (Ω) ′ ×M 0 nh (Ω) = b(v
V nh (Ω) ⊥ onto M 0 nh (Ω) ′ . Hence, for g n ∈ M 0 nh (Ω) ′ there exists a unique u n bh ∈ V nh (Ω) ⊥ such that ∀q ∈ M 0 nh (Ω), b(u n bh , q) = ⟨g n , q⟩ Γu and β ∥u n bh ∥ L 2 1 (Ω) 2 ≤ ∥g n ∥ L 2 1 (Γu) . ( 36
)
When setting

u n ⋄h = u n h -u n bh , thus, problem (32-33-34) is equivalent to find (u n ⋄h ) 0≤n≤N ∈ (V nh (Ω)) N +1 such that u 0 ⋄h = u 0 h -u 0 bh = Π 0h u 0 -u 0 bh , (37)
and, for all n, ∀1 ≤ n ≤ N,

∀v h ∈ V nh (Ω), (u n ⋄h , v h ) 1 -(u n-1 ⋄h , v h ) 1 + ατ n (u n ⋄h , v h ) 1 = (38) τ n (f n , v h ) 1 -τ n b(v h , i nh (p n b )). Theorem 4. For every data (f, p b ) belongs to C 0 (0, T ; L 2 1 (Ω) 2 ) × C 0 (0, T ; H s+ 1 2 1 (Γ p )), s > 1 2 , g belongs to C 0 (0, T ; L 2 1 (Γ u )) and u 0 belongs to L 2 1 (Ω) 2 satisfies divu 0 = 0 in Ω. Then problem (32-33-34) has a unique solution (u n h , p n h ) such that ∀ n, 0 ≤ n ≤ N, u n h ∈ X nh (Ω), and ∀ n, 1 ≤ n ≤ N, p n h ∈ M nh (Ω). Moreover, (u n h ) 0≤n≤N satisfies ∥u n h ∥ L 2 1 (Ω) 2 ≤ ∥u 0 ∥ L 2 1 (Ω) 2 + √ 2 α ( n ∑ m=1 τ m ( ∥f m ∥ 2 L 2 1 (Ω) 2 + c 2 0 ∥p m b ∥ 2 H s+ 1 2 1 (Γp) ) ) 1 2 + ∥g(., 0)∥ L 2 1 (Γu) + ∥g n ∥ L 2 1 (Γu)
. Proof. Applying estimate (24) to problem (32-33-34) and using the fact that Π 0h u 0 2 , we obtain the desired estimate.

L 2 1 (Ω) 2 ≤ u 0 L 2 1 (Ω)

A priori error estimates

To establish error estimates, we insert in the error equation an arbitrary element q n h ∈ M nh (Ω) and we obtain

∀v h ∈ X nh ; (Π nh u n -u n h , v h ) 1 + ατ n (Π nh u n -u n h , v h ) 1 + τ n b(v h , q n h -p n h ) = (u n-1 -u n-1 h , v h ) 1 -τ n b(v h , p n -q n h ), (39)
with the starting value Π 0h u 0 -u 0 h = 0 in Ω.

Proposition 1. We suppose that u 0 ∈ H s 1 (Ω) 2 and the solution

(u n , p n ) ∈ H s 1 (Ω) 2 × H s+1 1 (Ω) 2 , 1 2 < s ≤ 1. Then for all n, 1 ≤ n ≤ N ∥u n -u n h ∥ L 2 1 (Ω) 2 ≤ c(( n ∑ m=1 τ m (h m ) 2s ∥p m ∥ 2 H s+1 1 (Ω) ) 1 2
(40)

+ n ∑ m=0 (h m ) s ∥u m ∥ H s 1 (Ω) 2 ).
Proof. In the error equation (39), we choose

q n h = i nh (p n ) then, p n h -q n h belongs to M 0 nh (Ω), and for any v h ∈ V nh (Ω), (Π nh u n -u n h , v h ) 1 + ατ n (Π nh u n -u n h , v h ) 1 = (u n-1 -u n-1 h , v h ) 1 -τ n b(v h , p n -q n h ), we substract (Π n-1,h u n-1 , v h ) to this equation,we obtain for any v h ∈ V nh (Ω), ((Π nh u n -u n h ) -(Π n-1,h u n-1 -u n-1 h ), v h ) 1 + ατ n (Π nh u n -u n h , v h ) 1 (41) = -τ n b(v h , p n -q n h ) + (u n-1 -Π n-1,h u n-1 , v h ) 1 .
Recall that Π nh is the orthogonal projection operator from L 2 1 (Ω) 2 onto X nh , then Π nh u n | T ∈ P 0 (T ) 2 and ∀p ∈ P 0 (T ) 2 ∫ T (Π nh u nu n )p rdrdz = 0. For an arbitray q h ∈ M 0 nh (Ω), we choose p = gradq h ∈ P 0 (T ) 2 , then

∀q h ∈ M 0 nh (Ω), ∫ T (Π nh u n -u n )gradq h rdrdz = 0,
therefore by using equations (20) with q = q h ∈ M 0 nh (Ω) ⊂ H 1 1⋄ (Ω) and the second equation of (34) we obtain 41), and apply Cauchy-Schwarz inequality we obtain

∀q h ∈ M 0 nh (Ω), ∑ T ∈T nh ∫ T Π nh u n gradq h rdrdz = ∑ T ∈T nh ∫ T u n gradq h rdrdz, = ∑ T ∈T nh ∫ T u n h gradq h rdrdz, = ⟨g n , q h ⟩ Γu , then, Π nh u n -u n h ∈ V nh (Ω). Taking v h = Π nh u n -u n h in (
∥Π nh u n -u n h ∥ L 2 1 (Ω) 2 ≤ Π n-1,h u n-1 -u n-1 h L 2 1 (Ω) 2 + τ n |p n -q n h | H 1 1 (Ω) + u n-1 -Π n-1,h u n-1 h L 2 1 (Ω) 2 . Summing over n yields ∥Π nh u n -u n h ∥ L 2 1 (Ω) 2 ≤ n ∑ m=1 τ m |p m -q m h | H 1 1 (Ω) + n ∑ m=1 u m-1 -Π m-1,h u m-1 L 2 1 (Ω) 2 .
Therefore by triangle inequality and the fact that

q m h = i mh (p m ) we obtain ∥u n -u n h ∥ L 2 1 (Ω) 2 ≤ ∥u n -Π nh u n ∥ L 2 1 (Ω) 2 + ∥Π nh u n -u n h ∥ L 2 1 (Ω) 2 , ≤ n ∑ m=1 τ m |p m -i mh p m | H 1 1 (Ω) + n ∑ m=0 ∥u m -Π mh u m ∥ L 2 1 (Ω) 2 ,
Finaly, we apply the discrete Cauchy-Schwarz inequality, and use estimates ( 30) and ( 31) we obtain the desired a priori estimate.

Proposition 2. If the assumptions of Proposition 1 are satisfied, the following a priori error estimate holds for n, 1 ≤ n ≤ N,

1 τ n Π nh ((u n -u n h ) -(u n-1 -u n-1 h )) + grad(p n -p n h ) L 2 1 (Ω) 2 (42)
≤ c((

n ∑ m=1 τ m (h m ) 2s ∥p m ∥ 2 H s+1 1 (Ω) ) 1 2 + n ∑ m=0 τ m (h m ) s ∥u m ∥ H s 1 (Ω) 2 ) +c(h n ) s ∥p n ∥ H s+1 1 (Ω) .
Proof. We can rewrite the error equation (39) as follows

( 1 τn Π nh ((u n -u n h ) -(u n-1 -u n-1 h )) + grad(q n h -p n h ), v h ) 1 = -α(Π nh (u n -u n h ), v h ) 1 + b(v h , q n h -p n ).
We choose

v h = 1 τn Π nh ((u n -u n h ) -(u n-1 -u n-1 h )) + grad(q n h -p n h ), then 1 τn Π nh ((u n -u n h ) -(u n-1 -u n-1 h )) + grad(q n h -p n h ) L 2 1 (Ω) 2 ≤ c(∥u n -u n h ∥ L 2 1 (Ω) 2 + |q n h -p n | H 1 1 (Ω) ). (43)
On the other hand, by triangle inequality we have

1 τn Π nh ((u n -u n h ) -(u n-1 -u n-1 h )) + grad(p n -p n h ) L 2 1 (Ω) 2 ≤ 1 τn Π nh ((u n -u n h ) -(u n-1 -u n-1 h )) + grad(q n h -p n h ) L 2 1 (Ω) 2 +|q n h -p n | H 1 1 (Ω)
. Finaly by using estimates (43), (40) and the approximate properties of M nh (Ω), we obtain the a priori estimate (42).

A posteriori analysis

For the time discretization and the space discretization, we describe a family of error indicators and prove upper and lower bounds for the error.

The time discretisation

For each n, 1 ≤ n ≤ N, we define the time error indicator, see [START_REF] Bergam | A posteriori analysis of the finite element discretization of some parabolic equations[END_REF] and [START_REF] Johnson | An a posteriori error estimate and adaptive timestep control for a backward Euler discretization of a parabolic problem[END_REF]

η n = ( τ n 3 ) 1 2 u n h -u n-1 h L 2 1 (Ω) 2 . ( 44
)
Let u τ denote the function which is continous, affine on each interval

[t n-1 , t n ], 1 ≤ n ≤ N, such that ∀n, 0 ≤ n ≤ N, u τ (t n ) = u n ,
and p τ denote the piecewise constant function such that

∀n, 1 ≤ n ≤ N, ∀t ∈]t n-1 , t n ], p τ (t) = p(t n ).
Then for all t in ]t n-1 , t n ], the residual equation in variational form reads

∀v ∈ L 2 1 (Ω) 2 ; (∂ t (u -u τ ), v) 1 + α(u -u τ , v) 1 + b(v, p * -Π τ p * ) = (f -Π τ f, v) 1 -α(u τ -u n , v) 1 -b(v, pb -Π τ pb ) (45) ∀q ∈ H 1 1⋄ (Ω); b(u -u τ , q) = ⟨g -Π τ g, q⟩ Γu , (46)
where Π τ is introduced in Remark 2 by

∀t ∈]t n-1 , t n ], (Π τ v)(t) = v(t n ),
and we recall that p = p * + pb . Let the regularity parameter

σ τ = max 1≤n≤N τ n τ n-1 ,
where τ 0 = τ 1 .

The reliability of the indicator

Proposition 3. The following a posteriori error estimate holds, for

1 ≤ n ≤ N, ∥u(., t n ) -u n ∥ L 2 1 (Ω) 2 ≤ √ 2α ( n ∑ m=1 η 2 m ) 1 2 (47)
+ √ 1 α ( ∥f -Π τ f∥ L 2 (0,tn;L 2 1 (Ω) 2 ) + c 0 ∥p b -Π τ p b ∥ L 2 (0,tn;H 1 2 1 (Γp))
)

+ √ α(1 + σ 1 2 τ ) ( n ∑ m=0 τ m ∥u m -u m h ∥ 2 L 2 1 (Ω) 2 ) 1 2
Proof. Applying estimate [START_REF] Verfürth | A Review of A Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques[END_REF] to problem (45-46) at time t = t n and using the fact that (u -u τ )(0) = 0, (g -Π τ g)(0) = 0, (g -Π τ g)(., t n ) = 0 and

u n = u τ (t n ) = Π τ u τ (t) we obtain ∥u(., t n ) -u n ∥ 2 L 2 1 (Ω) 2 ≤ α ∥(u τ -Π τ u τ )∥ 2 L 2 (0,tn;L 2 1 (Ω) 2 ) + 1 α ( ∥f -Π τ f∥ 2 L 2 (0,tn;L 2 1 (Ω) 2 ) + ∥p b -Π τ pb ∥ 2 L 2 (0,tn;H 1 1 (Ω) ) . ( 48
)
To estimate the first term in the right-hand side, we observe that on the interval

]t n-1 , t n ], (u τ -Π τ u τ )(t) = -tn-t τn (u n -u n-1
), thus by integrating this equation between t n-1 and t n and using the fact that τ n = t nt n-1 , we obtain

∥u τ -Π τ u τ ∥ 2 L 2 (tn-1,tn;L 2 1 (Ω) 2 ) = u n -u n-1 2 L 2 1 (Ω) 2 ∫ tn tn-1 ( t n -t τ n ) 2 dt, thus ∥u τ -Π τ u τ ∥ L 2 (tn-1,tn;L 2 1 (Ω) 2 ) = ( τ n 3 ) 1 2 u n -u n-1 L 2 1 (Ω) 2 . ( 49
)
On the other hand, triangle inequality yields

u n -u n-1 L 2 1 (Ω) 2 ≤ ∥u n -u n h ∥ L 2 1 (Ω) 2 + u n h -u n-1 h L 2 1 (Ω) 2 + u n-1 h -u n-1 L 2 1 (Ω) 2 . Multipliying by ( τn 3
) 1 2 and using the expression of the error indicator (44) gives

( τn 3 ) 1 2 u n -u n-1 L 2 1 (Ω) 2 ≤ ( τn 3 ) 1 2 ∥u n -u n h ∥ L 2 1 (Ω) 2 + η n + ( τ n 3 ) 1 2 u n-1 h -u n-1 L 2 1 (Ω) 2 . ( 50
)
The expression of σ τ yields

( τn 3 ) ( τ n-1 3 ) ≤ σ τ , thus ( τn 3 ) 1 2 ≤ ( τn-1 3 ) 1 2 .(σ τ ) 1
2 . Therefore by substituting the previous inequality and (50) in (49) we obtain

∥u τ -Π τ u τ ∥ L 2 (tn-1,tn;L 2 1 (Ω) 2 ) ≤ ( τn 3 ) 1 2 ∥u n -u n h ∥ L 2 1 (Ω) 2 + η n + ( τ n-1
3 2 , Summing over n with 1 ≤ n ≤ N , the square of this inequality we obtain

) 1 2 .(σ τ ) 1 2 u n-1 h -u n-1 L 2 1 (Ω)
∥u τ -Π τ u τ ∥ 2 L 2 (0,tn;L 2 1 (Ω) 2 ) ≤ 2 n ∑ m=1 η 2 m + (1 + σ τ ) n ∑ m=0 τ m ∥u m -u m h ∥ 2 L 2 1 (Ω) 2 . ( 51
)
Finaly by substituting (51) in (48) and using estimate [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations. Theory and Algorithms[END_REF] we obtain the desired a posteriori error estimate.

Proposition 4. The following a posteriori error estimate holds, for 1 ≤ n ≤ N,

∥∂ t (u -u τ )∥ L 2 (0,tn;L 2 1 (Ω) 2 ) ≤ c(∥f -Π τ f∥ L 2 (0,tn;L 2 1 (Ω) 2 ) + ( n ∑ m=1 η 2 m ) 1 2
(52)

+c 0 ∥p b -Π τ p b ∥ L 2 (0,tn;H 1 2 1 (Γp)) + (1 + √ σ τ )( n ∑ m=0 τ m ∥u m -u m h ∥ 2 L 2 1 (Ω) 2 ) 1 2 ).
Proof. We take v equal to ∂ t (u -u τ ) in ( 45) and apply Cauchy-Schwarz inequality we obtain

1 2 ∥∂ t (u -u τ )∥ 2 L 2 1 (Ω) 2 + α 2 d dt ∥u -u τ ∥ 2 L 2 1 (Ω) 2 ≤ ∥f -Π τ f∥ 2 L 2 1 (Ω) 2 +α 2 ∥u τ -Π τ u τ ∥ 2 L 2 1 (Ω) 2 + ∥p b -Π τ pb ∥ 2 H 1 1 (Ω)
. Integrating between 0 and t n , using estimate [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations. Theory and Algorithms[END_REF] and the fact that (u -u τ )(0) = 0, we obtain

∥∂ t (u -u τ )∥ 2 L 2 (0,tn;L 2 1 (Ω) 2 ) + α ∥(u -u τ )(t n )∥ 2 L 2 1 (Ω) 2 ≤ 2(∥f -Π τ f∥ 2 L 2 (0,tn;L 2 1 (Ω) 2 ) +α 2 ∥u τ -Π τ u τ ∥ 2 L 2 (0,tn;L 2 1 (Ω) 2 ) + c 2 0 ∥p b -Π τ p b ∥ 2 L 2 (0,tn;H 1 2 1 (Γp))
).

This last estimate and (51) implies (52).

Proposition 5. The following a posteriori error estimate holds, for

1 ≤ n ≤ N ∥p -p τ ∥ L 2 (0,tn;H 1 1 (Ω)) ≤ c(∥f -Π τ f∥ L 2 (0,tn;L 2 1 (Ω) 2 ) +c 0 ∥p b -Π τ p b ∥ L 2 (0,tn;H 1 2 1 (Γp)) +( n ∑ m=1 η 2 m ) 1 2 + (1 + √ σ τ )( n ∑ m=0 τ m ∥u m -u m h ∥ 2 L 2 1 (Ω) 2 ) 1
2 ).

Proof. From equation (45) we have for all

v ∈ L 2 1 (Ω) 2 , b(v, p * -Π τ p * ) = (∂ t (u τ -u), v) 1 + α(u τ -u, v) 1 + b(v, Π τ pb -pb ) +(f -Π τ f, v) 1 + α(u n -u τ , v) 1 .
Cauchy-Schwarz inequality and the inf-sup condition [START_REF] Bernardi | Discretisation of an unsteady flow through a porous solid modeled by Darcy's equations[END_REF] yields

β ∥p * -Π τ p * ∥ H 1 1 (Ω) ≤ ∥∂ t (u τ -u)∥ L 2 1 (Ω) 2 +α ∥u τ -u∥ L 2 1 (Ω) 2 +∥p b -Π τ pb ∥ H 1 1 (Ω) + ∥f -Π τ f∥ L 2 1 (Ω) 2 + α ∥u n -u τ ∥ L 2 1 (Ω) 2 .
Integrating between 0 and t n and the fact that u n = Π τ u τ we obtain

∥p * -Π τ p * ∥ 2 L 2 (0,tn;H 1 1 (Ω)) ≤ c(∥∂ t (u τ -u)∥ 2 L 2 (0,tn;L 2 1 (Ω) 2 ) + c 2 0 ∥p b -Π τ p b ∥ 2 L 2 (0,tn;H 1 2 1 (Γp)) + ∥f -Π τ f∥ 2 L 2 (0,tn;L 2 1 (Ω) 2 ) + α 2 ∥Π τ u τ -u τ ∥ 2 L 2 (0,tn;L 2 1 (Ω) 2 )
). Finaly by substituting (51) and (52) in the last inequality, using the fact that pp τ = (p * -Π τ p * ) + (p b -Π τ pb ) and triangle inequality we obtain the desired a posteriori error estimate.

The space discretisation

For each T ∈ T nh , we associate • E T the set of edges of T ,

• E 0 T = E T ∩ E 0 nh , • h e the diameter of e, • E Γu nh = {e ∈ E nh ; e ⊂ Γ u }
, where E 0 nh and E nh are defined in section 3.2. For each n, 1 ≤ n ≤ N and each T ∈ T nh , we define the following error indicators

η n T = 1 τ n u n-1 h -Π nh u n-1 h L 2 1 (T ) 2 and η n ∂T = ∑ E Γu nh h 1 2 e ∥[u n h .n e ] e ∥ L 2 1 (e) ,
where the jumps [u n h .n e ] e are constant on each e. Moreover, in the context of mesh adaptivity, the term u n-1 h -Π nh u n-1 h only differs from zero in the elements T of T nh that are the union of several elements of T n-1,h . Therefore these indicators can be computed readily and explicitly.

We approximate the boundary data p n b by the Lagrange interpolation operator i nh , with values in M 0 nh (Ω), i.e. for each continuous function q belongs to Γ p , i nh q is a piecewise affine function equal to q on each node of V b nh . In order to prove the a posteriori estimates, we first write the residual equations. We recall that p n = p

∈ L 2 1 (Ω) 2 ; (u n -u n h , v) 1 + ατ n (u n -u n h , v) 1 + τ n b(v, p n * -p n * h ) = (u n-1 -u n-1 h , v) 1 + (u n-1 h -Π nh u n-1 h , v) 1 + τ n (f n -Π nh f n , v) 1 -τ n b(v, L(p n b -i nh p n b )), ( 54
) ∀q ∈ H 1 1⋄ (Ω); b(u n -u n h , q) = ⟨g n , q -q n h ⟩ Γu - 1 2 ∑ T ∈T nh ∑ E Γu nh ∫ e [u n h .n e ] e (τ )(q -q n h )(τ )dτ.
Proof. Taking v h equal to χ T e in (34), where χ T is the characteristic function of T and e runs through the canonical basis of base R 2 we obtain

∀T ∈ T nh , (u n h , χ T e) 1 + ατ n (u n h , χ T e) 1 + τ n b(χ T e, p n h ) = (u n-1 h , χ T e) 1 + τ n (f n , χ T e) 1 . Then, u n h + ατ n u n h + τ n gradp n h = u n-1 h + τn meas(T ) ∫ T f n rdrdz = Π nh u n-1 h + τ n Π nh f n . Multiplying this equation by any v ∈ L 2 1 (Ω) 2
, integrating on each T ∈ T nh , and summing over all elements T ∈ T nh , we obtain

∑ T ∈T nh ∫ T (u n h +ατ n u n h +τ n gradp n h ).v rdrdz = ∑ T ∈T nh ∫ T (Π nh u n-1 h + τ n Π nh f n ).v rdrdz,
Finaly subtracting this equality from the first equation of (20) and using the fact that p n = p n * + pn b where pn

b = L(p n b ) and p n * h = p n h -L(i nh p n b )
, we obtain the first equation of (54). On the other hand, the second equations of (20) and (34) gives ∀q n h ∈ M 0 nh (Ω), b(u nu n h , q n h ) = 0. Then by Green's formula, we obtain

∀q n h ∈ M 0 nh (Ω), b(u n -u n h , q) = b(u n -u n h , q -q n h ), = b(u n , q -q n h ) -b(u n h , q -q n h ), = ⟨g n , q -q n h ⟩ Γu - ∑ T ∈T nh ∫ T u n h grad(q -q n h ) rdrdz, = ⟨g n , q -q n h ⟩ Γu - ∑ T ∈T nh ∫ ∂T ∩(Ω∪Γu) (u n h .n)(τ )(q -q n h )(τ ) dτ,
whence the second line in (54).

The reliability of the indicators Proposition 7. The following a posteriori error estimate holds between the solutions (u

n , p n ) 1≤n≤N of problem (18-19-20) and (u n h , p n h ) 1≤n≤N of problem (32-33-34), for n, 1 ≤ n ≤ N , ∥u n -u n h ∥ 2 L 2 1 (Ω) 2 ≤ ∥u 0 -Π 0h u 0 ∥ 2 L 2 1 (Ω) 2 + c(G 2 n + ∑ T ∈T nh (η n ∂T ) 2 ) (55) + 2 α n ∑ m=1 τ m (J 2 m + ∑ T ∈T nh (η m T ) 2 ),
where

J n = ∥f n -Π nh f n ∥ L 2 1 (Ω) 2 , G n = ∥p n b -i nh p n b ∥ H 1 2 1 (Γp)
.

Proof. To simplify, let

w n = u n -u n h , r n * = p n * -p n * h and F n = f n -Π nh f n -gradL(p n b -i nh p n b ) + 1 τn (u n-1 h -Π nh u n-1 h
). Therefore the residual equations (54) become

∀v ∈ L 2 1 (Ω) 2 , (w n , v) 1 + ατ n (w n , v) 1 + τ n b(v, r n * ) = (w n-1 , v) 1 + τ n (F n , v) 1 , ∀q ∈ H 1 1⋄ (Ω), b(w n , q) = ⟨g n , q -R n h q⟩ Γu (56) - 1 2
∑ T ∈T nh ∑ e∈E Γu nh ∫ e [u n h .n e ] e (τ )(q -R n h q)(τ )dτ,
where R n h denotes a Clément type regularization operator with values in M 0 nh such as the Scott and Zhang operator [START_REF] Scott | Finite element interpolation of non-smooth functions satisfying boundary conditions[END_REF]. This operator preserves the zero boundary trace and satisfies for each T ∈ T nh , and e ∈ E Γu nh see [8, Cor. IX.3.9], [START_REF] Scott | Finite element interpolation of non-smooth functions satisfying boundary conditions[END_REF] and also [START_REF] Belhachmi | Weighted Clement operator and application to the finite element discretization of the axisymmetric Stokes problem[END_REF] for the extension to weighted spaces,

∀q ∈ H 1 1 (Ω), ∥q -R n h q∥ L 2 1 (e) ≤ ch 1 2
e ∥q∥ H 1 1 (∆e) , where ∆ e is an appropriate neighbourhood of e. Then from this inequality, there exists a unique µ n ∈ H 1 1⋄ (Ω) such that ∀q ∈ H 1 1⋄ (Ω), (gradµ n , gradq) 1 = b(w n , q), and

|µ n | H 1 1 (Ω) ≤ c    ∑ T ∈T nh ∑ e∈E Γu nh h e ∥[u n h .n e ] e ∥ 2 L 2 1 (e)    1 2

The efficiency of the indicators

We will prove an upper bound for the error indicators. For each T ∈ T nh , let ω T denote the union of triangles in T nh that share at least an edge with T .

Proposition 9. For each n, 1 ≤ n ≤ N and T ∈ T nh ,

η n T ≤ α ∥u n -u n h ∥ L 2 1 (T ) 2 + (u n -u n h ) -(u n-1 -u n-1 h ) τ n + grad(p n -p n h ) 2 L 2 1 (T ) 2 + J n T ,
where

J n T = ∥f n -Π nh f n ∥ L 2 1 (T ) 2 . Proof. Taking v = (u n-1 h -Π nh u n-1 h
)χ T in the first equation of (54), where χ T is the charateristic function of T and using the fact that

p n * -p n * h = p n -p n h -L(p n b -i nh p n b ) we obtain u n-1 h -Π nh u n-1 h 2 L 2 1 (T ) 2 = ατ n (u n -u n h , u n-1 h -Π nh u n-1 h ) 1,T +τ n ( (u n -u n h ) -(u n-1 -u n-1 h ) τ n +grad(p n -p n h ), u n-1 h -Π nh u n-1 h ) 1,T -τ n ∫ T (f n -Π nh f n )(u n-1 h -Π nh u n-1 h )rdrdz.
Cauchy-Schwarz inequality yields Proof. By means of a fixed lifting operator on the reference element T and by using the affine transformation that maps T onto T , we construct for each e ∈ E T a lifting operator L e,T such that for each polynomial φ on e vanishing on ∂e, L e,T φ is a polynomial on T vanishing on ∂T \ e and equal to φ on e. Let b e denote the bubble function on e, i.e., the product of the barycentric coordinates associated with the vertices of e. For each e ∈ E 0 T , we denote by T ′ the other element of T nh that contains e. In the second equation of (54), we take q n h = 0 and q = q n e , with Then, we obtain b(u nu n h , q n e ) = ⟨g n , q n e ⟩ Γu - This gives gives the bound for the second indicator η n ∂T .

u n-1 h -Π nh u n-1 h L 2 1 (T )
In both Propositions 9 and 10, the estimates are local in space and time, so that it can be thought that the indicators η n T and η n ∂T provide a good tool for adapting the mesh.

  pn b ) is a solution of problem (18-19-20). Then we can conclude the following lemma: Lemma 4. The variational problems (18-19-20) and (23) are equivalent. Theorem 3. For any data (f, p b

  T ([u n h .n e ] e b e ) on T, L e,T ′ ([u n h .n e ] e b e ) on T ′ , 0 elsewhere.

Figure 2 :

 2 Figure 2: Isovalues of the initial radial velocity

  For any solutions (u n , p n ) 1≤n≤N of problem (18-19-20) and (u n

n * + pn b and p n * h = p n h -L(i nh p n b ), where pn b = L(p n b ) and p n * h is no longer a piecewise polynomial function. Lemma 5. h , p n h ) 1≤n≤N of problem (32-33-34), then (u nu n h , p np n h ) 1≤n≤N satisfies the residual equation ∀v

  +τ n ∥f n -Π nh f n ∥ L 2 1 (T ) 2 .Finaly by multipling this inequality by1 τn and from the expression of the a posteriori indicator η n T , we get the desired estimate.

	τ n	h	)	+ grad(p n -p n h )	L 2 1 (T ) 2

2 ≤ ατ n ∥u nu n h ∥ L 2 1 (T ) 2 +τ n (u nu n h ) -(u n-1 -u n-1 Proposition 10. For each n, 1 ≤ n ≤ N and T ∈ T nh , η n ∂T ≤ c ∥u nu n h ∥ L 2 1 (ωT ) 2 .

  L e,T ([u n h .n e ] e b e ).χ T + L e,T ′ ([u n h .n e ] e b e ).χ T ′ = 2[u n h .n e ] e b e . -u n ∥ L 2 1 (T ∪T ′ ) 2 |q n e | H 1 1 (T ∪T ′ ) . (58)Recall the following inverse inequality, for each constant λ, see [16, Lem. 3.3], ∪T ′ ) 2 . |L e,T ([u n h .n e ] e b e )| ∪T ′ ) 2 . ∥[u n h .n e ] e ∥

										1 2	T ∈T nh ∑	  	nh e∈E Γu ∑	∫	e	[u n h .n e ] e (τ )q n e (τ )dτ	   .
	On the other hand,
	q n e (τ ) = Thus,
	[u n h .n e ] e b h ∥λ∥ L 2 1 2 e 2 L 2 1 (e) ≤ ∥u n 1 (e) ≤ c λb 1 2 e L 2 1 (e)	-1 2 e ∥λ∥ L 2 1 (e) .
	The expression of η n ∂T and the first inverse inequality, yields  
							η n ∂T ≤ c	  e∈E Γu ∑ nh	h	1 2	1 2 e	L 2 1 (e)	  .
	Using estimate (58) we obtain
									
	η n ∂T ≤ c	 	∑		h e ∥u n 1 2 h -u n ∥ L 2 1 2 1 (T 1 2 H 1 1 (T ∪T ′ )	  .
						e∈E Γu nh
	By substituting the second inverse inequality into this last inequality we obtain  
	η n ∂T ≤ c	 	∑	h e ∥u n 1 4 h -u n ∥	1 2 L 2 1 (T 1 2 L 2 1 (e)	  ,
				e∈E Γu nh
									
	≤ c	 	1 2	∑ e∈E Γu nh	∥u h -u n ∥ L 2 1 (T ∪T ′ ) 2 +	1 2	nh e∈E Γu ∑	h e ∥[u n 1 2 h .n e ] e ∥ L 2 1 (e)	  .

and |L e,T (λb

e )| H 1 1 (T ) ≤ ch e [u n h .n e ] e b

  Figure 1: Isovalues of the initial axial velocity

	Graph u
	IsoValue 0.147869 1.43725 2.29684 3.15643 4.01602 4.87561 5.7352 6.59478 7.45437 8.31396 9.17355 10.0331 10.8927 11.7523 12.6119 13.4715 14.3311 15.1907 16.0503 18.1992
	Graph u
	IsoValue -0.765254 0.393722 1.16637 1.93902 2.71167 3.48432 4.25697 5.02962 5.80228 6.57493 7.34758 8.12023 8.89288 9.66553 10.4382 11.2108 11.9835 12.7561 13.5288 15.4604

) + (η m ∂T ) 2 ) ).

Acknowledgements

We thank Prof. Frédéric Hecht for his help in the numerical part of this work. His comments and suggestions were very useful to complete this research work.

Corollary 2. The folowing a posteriori error estimate holds for 1 ≤ n ≤ N, ∥∂ t (u -u τ )∥ L 2 (0,tn;L 2 1 (Ω) 2 ) + ∥pp τ ∥ L 2 (0,tn;H 1 1 (Ω))

≤ c(∥f -Π τ f∥ L 2 (0,tn;L 2 1 (Ω) 2 ) + ∥p b -Π τ p b ∥ L 2 (0,tn;H

2 ).

The efficiency of the indicator

Proposition 6. Each indicator η n , 1 ≤ n ≤ N, defined in (44) satisfies the following bound

Moreover this estimate is local with respect to the time variable.

Proof. By the expression of the indicator η n , triangle inequality, and the fact that

3

In order to evaluate the first term in the right-hand side, we take v = u τ -u n in (45) and using the fact that p * + pb = p, Π τ p = p τ and u n = Π τ u τ , we obtain

). Integrating this inequality between t n-1 and t n and using (49) we obtain

). Finally by substituting the previous inequality in (53) we obtain the desired estimate.

Some numerical experiments

We present some numerical experiments realized with the code FreeFem++, see [START_REF] Hecht | New development in freefem++[END_REF]. The domain Ω is generated by the L-shaped meridian domain Ω defined by

We denote by Γ p the intersection of Γ with the plan z = 0 and Γ u is equal to Γ \ Γ p .

We work with the data g equal to one on

We take p(0, z) = 0 on Γ p , and the initial values of the velocity u in Ω at t = 0 are such that u 0 z = 0 and u 0 r = 1. The data f are such that f r = 1 and f z = 0.

Finally we take α equal to 0.25, dt = 0.05 and T = 1. Figure 1 presents, the curves of isovalues of the initial axial velocity u 0 z . Figure 2 presents, the curves of isovalues of the initial axial velocity u 0 r . Figure 3 presents, the curves of isovalues of the pressure at time T = 1. Figure 4 presents, the curves of isovalues and directions of the velocity at time T = 1.


A Garcia

S Rodriguez

A Le Gall

S Courrech Du Pont

C Narteau

S Le Mouélic

A Lucas

J Radebaugh

K Arnold

J W Barnes

C Sotin

R H Brown

R D Lorenz

E P Turtle

Global mapping and characterization of Titan's dune fields with Cassini: correlation between RADAR and VIMS observations

Introduction: Titan, the largest moon of Saturn, is the only one to possess a dense atmosphere and a great variety of morphological structures on its surface. Thanks to the Cassini-Huygens mission and its microwaves and infrared imaging capabilities (RADAR in SAR mode, VIMS and ISS) large dunes fields were discovered in the equatorial region. These features represent clues for the understanding of the climatic history on Titan [1,2]. Gathering the spectral ability of the VIMS instrument and the locations and morphology of dunes available with the spatial resolution of RADAR/SAR, we are now able to map the geographic distribution of Titan's dunes, to precisely quantify their coverage and to provide some constraints on the nature of the material that compose them; with the final aim to figure out the mystery of the processes that are implied on this satellite and to understand its complex climatology. In this paper, and according to previous works [3-7], we will present updated global maps in infrared and microwaves wavelength and updated estimates of dune coverage. We will also discuss the degree of correlation between the dunes imaged by the RADAR/SAR and some spectral units appearing in VIMS observations and discuss briefly the implication for dune material and compositional properties.

Dune coverage from RADAR/SAR observations:

From October 2004 to June 2011 (Titan flybys TA-T77), the RADAR/SAR from best to coarser resolution imaged 52% of the surface of Titan (see fig. 1). The highest resolution of RADAR/SAR swaths offers the best detection of dunes fields. We find that dunes cover 11.2% of the SAR images of highest resolution (350m/pixel), which cover in total 40% of Titan's surface. Some dunes are also revealed at lower resolution (700m/pixel) accounting for 0.8% of the whole imaged surface. Putting all together, it finally corresponds to dune coverage of 9.3% of the area imaged by RADAR/SAR. This confirms that dunes are one of the most dominant feature on the Titan's surface (equivalent to 2/3 of the Europe area).

The identification and delineation of the dunes at the global scale confirm the previous studies of the particular distribution of this pattern essentially between 30°N-30°S (as proposed by [1] from TA-T55, see fig. 1 b)) with some isolated dunes at 50°N. The figure 1 b. shows in particular that, at the equator, although the RADAR/SAR has a lower coverage on average, dunes cover 60% of the RADAR swaths. Correlation RADAR/VIMS and discussion: It has been suggested that global maps VIMS in composite color of band ratio (1.59/1.27µm, 2.03/1.27µm, 1.27/1.08µm) enhances spectral heterogeneity between terrains with two particular units: 1) Dark blue units which are suspected to be enrich in fine grain water ice [7, 8] and a dark brown unit which would be strongly correlated to the dunes imaged by the RADAR/SAR [5,6]. In our up-to-date map, we outlined the "blue" and "brown" units from the map of [9] (see fig 2. a)). We compared the latitude distribution of the VIMS "brown" unit with the latitudinal diagram distribution given previously (see fig. 2 b)) and fig. 1 b)). We can conclude that they seem to be highly similar, at least in term of latitudinal extend.

Moreover, by overlapping all the sets of data (RADAR/SAR and VIMS) available so far, we observe an excellent geographic matching between the dunes imaged by the RADAR/SAR and the VIMS "brown" unit. Indeed, 77% of the dunes are found to be present in this unit. With another statistic analysis, we also find that 70% of the area covered by the SAR swaths overlapping the "brown" unit, is covered by dunes. Our calculations consolidate the main hypothesis of correlation between the VIMS "brown" unit and the dunes imaged by the RADAR/SAR and demonstrate the high level of homogeneity of the dunes material, at the global scale. This is fully consistent with the uniformity of dune reflectivity and emissivity given by the RADAR [10,11]. However, 1.5% of dunes are present in the VIMS "blue" unit while 20% of the SAR swaths belonging to this unit are covered by the dunes. The precise nature and origin of these bluish dunes are still unclear. They have been interpreted as enriched in water ice or dunes fields, which might present a lack of sediments in the inter-dune areas [4,12].

According this analysis, we can with reasonably good confidence, extrapolate the total extent of dune fields of Titan to the total area of the VIMS "brown" unit, which extends the dune system to 18.5% of the whole surface of Titan (more than 15 million km 2 , almost two times the total area of the USA). Considering the possible composition of Titan's "sand", dunes therefore represent the largest hydrocarbon reservoir on the surface of Titan with an estimated volume of sediments equivalent to 500.000 km 3 . This is consistent with the range of organic sediment inventory given by [13] who established liquid hydrocarbon reservoir (lakes and seas) as the second organic reservoir on Titan. Further works will be carried out to investigate more deeply with RADAR, VIMS and ISS observations several "blue" and "brown" areas to better understand the physics of their setting and the type of "sand" they are made of with the integration of recent cartographic products [14][15][16].

Figure 1 :

 1 Figure 1: a) Global map of all the RADAR/SAR swaths from TA to T77, with dunes outlined in brown -b) Dune latitudinal distribution: The surface area covered by the dunes is normalized by the surface imaged by the RADAR/SAR TA-T77 (latitude bin=2°). The dashed line corresponds to the fractional surface area of surface imaged.

  References: [1] Lorenz, R.D. et al. (2006), Science, Vol. 312 [2] Lorenz, RD. and Radebaugh, J. (2009), GRL, Vol. 36, L03202 [3] Le Corre, L. et al. (2008), Abstract [4] Le Gall, A. et al. (2011), Icarus 608-624 [5] Barnes, J. et al. (2007), Icarus 186, 242-258 [6] Le Corre, L. (2009), thesis [7] Rodriguez, S. et al. (2006), Planetary and Space science 54, 1510,-1523 [8] Le Mouélic, S. et al. (2008), Journal of Geophysical research, Vol. 113 [9] Brown, R.H. et al. (2009), Springer [10] Neish, C.D., et al. (2010), Icarus 192, 77-91 [11] Janssen, M.A. et al. (2011), Icarus [12] Barnes, J. et al. (2009), Icarus 195, 400-414 [13] Lorenz et al. (2008), GRL Vol 35, L02206 [14] Rubin, D.M. and Hesp, P.A. (2009), Nature Geoscience, Vol. 2 [15] Radebaugh, J. (2009), Nature Geoscience, Vol. 2 [16] Le Mouélic et al. (2012), Planet Space Science,73, 178-190.

Figure 2

 2 Figure 2: a) Map of the brown and blue units of VIMS extracted from the global VIMS map [9] -b) Latitudinal distribution of the brown unit.