Alexandre Ern
email: ern@cermics.enpc.fr

Sébastien Meunier
email: sebastien.meunier@edf.fr

Mathematical Modelling and Numerical Analysis Will be set by the publisher Modélisation Mathématique et Analyse Numérique

Keywords: 1991 Mathematics Subject Classification. 65M60, 65M15, 74F10 Finite element method, a priori error analysis, a posteriori error analysis, Hydro-Mechanical coupling

We analyze Euler-Galerkin approximations (conforming finite elements in space and implicit Euler in time) to coupled PDE systems in which one dependent variable, say u, is governed by an elliptic equation and the other, say p, by a parabolic-like equation. The targeted application is the poroelasticity system within the quasi-static assumption. An abstract setting is proposed to identify a natural energy norm for the PDE system. Two a posteriori error analyzes are performed, both yielding reliable upper error bounds in the sense that all the constants are specified. The first analysis hinges directly on the stability of the continuous problem and can be used to estimate the dominant term associated with the p-component in the energy norm. The second analysis is an extension of the elliptic reconstruction technique introduced by Makridakis and Nochetto [SIAM J. Numer. Anal., 41(4), [1585][1586][1587][1588][1589][1590][1591][1592][1593][1594] 2003] for linear parabolic problems. It is used here to derive an a posteriori error estimate for the u-component in the energy norm that exhibits an optimal convergence rate with respect to mesh size. Numerical results are presented to illustrate the performance of the various estimators.

Introduction

The main motivation for this work is the performance assessment by numerical simulations of underground storage facilities for nuclear waste. The processes involved in near-field models are extremely complex and can involve multi-phase, multicomponent flows through a porous medium subjected to thermal, hydraulic, chemical, and mechanical couplings. Here, we focus on poroelasticity problems involving hydro-mechanical couplings. We consider a linearly elastic and porous medium Ω saturated by a slightly compressible and viscous fluid within the so-called quasi-static assumption in which inertia effects in the elastic structure are negligible. Given a simulation time T > 0, the problem consists of finding a displacement field u : [0, T ] × Ω → R 3 and a pressure field p : [0, T ] × Ω → R such that

-∇•σ(u) + b∇p = f, in [0, T ] × Ω, (1)
Here, σ(u) = 2λ 1 ε(u) + λ 2 (∇•u)I is the so-called effective stress tensor, ε(u) = 1 2 (∇u + (∇u) t ) the linearized strain tensor, λ 1 and λ 2 the Lamé coefficients, I the identity matrix in R 3 , b the Biot-Willis coefficient, M the Biot modulus, κ the permeability of the medium, while f and g are given data. The system (1)-( 2) is supplemented by initial and boundary conditions discussed below. The poroelasticity system can be traced back to the pioneering work of Terzaghi [START_REF] Terzaghi | Theoretical Soil Mechanics[END_REF] and Biot [START_REF] Biot | General theory of three-dimensional consolidation[END_REF]. Equations ( 1)-( 2) respectively express the balance of momentum and the conservation of mass. The quasi-static assumption means that the term ρ∂ tt u (where ρ denotes the density of the elastic structure) has been neglected in the momentum balance. The Biot modulus combines compressibility and porosity effects; it is often assumed to be very large when dealing with the socalled Biot's consolidation problem, but this assumption will not be made here. For the sake of simplicity, we will assume that the coefficients b and M are given constants. A mathematical analysis of the system (1)- [START_REF] Bergam | A posteriori analysis of the finite element discretization of some parabolic equations[END_REF], including existence and uniqueness of strong and weak solutions based on the theory of linear degenerate evolution equations in Hilbert spaces, has been carried out by Showalter [START_REF] Showalter | Diffusion in deformable media[END_REF][START_REF] Showalter | Diffusion in poro-elastic media[END_REF]. Boundary conditions can be prescribed by considering two partitions of the boundary. The first partition is used for the displacement field (either the displacement itself or a traction force is prescribed), while the other partition is used for the pressure field (either the pressure itself or a flux is prescribed). For the sake of simplicity, we assume here that any portion of the boundary is clamped or drained, i.e. at least a Dirichlet condition is enforced on the displacement or on the pressure everywhere. Furthermore, an initial condition must be enforced on the quantity 1 M p + b∇•u. Although the evolution problem related to (1)-( 2) is essentially of parabolic type under minimal smoothness requirements on the data, we refer to it as a coupled elliptic-parabolic problem to stress the fact that Equ. ( 1) is of elliptic type for the displacement and Equ. ( 2) is of parabolic type for the pressure.

In the present work, we assume that the data (including boundary and initial conditions) are smooth enough for a strong solution to exist, and we shall mainly be concerned with the a posteriori error analysis of Euler-Galerkin approximations to the exact solution obtained by using a backward Euler scheme in time and conforming finite element methods in space. In the sequel, we restrict our attention to the case where the exact solution is smooth up to the initial time, as is customary in the a posteriori error analysis of Euler-Galerkin approximations to parabolic problems such as the heat equation. The a priori analysis of Euler-Galerkin approximations for Biot's consolidation problem is covered in the work of Murad, Loula, and coworkers [START_REF] Murad | Improved accuracy in finite element analysis of Biot's consolidation problem[END_REF][START_REF] Murad | On stability and convergence of finite element approximations of Biot's consolidation problem[END_REF][START_REF] Murad | Asymptotic behavior of semidiscrete finite-element approximations of Biot's consolidation problem[END_REF], including the semi-discrete and fully discrete cases and long-time behavior. The problem under scrutiny here is somewhat different since we do not assume that the Biot modulus takes very large values, i.e. we do not discard the pressure time-derivative in [START_REF] Bergam | A posteriori analysis of the finite element discretization of some parabolic equations[END_REF]. As a result, we shall briefly address below the a priori error analysis of the Euler-Galerkin approximation to the evolution problem (1)- [START_REF] Bergam | A posteriori analysis of the finite element discretization of some parabolic equations[END_REF]. The natural stability norm associated with this problem controls on the one hand the L ∞ t (H 1 x )-norm (L ∞ in time and H 1 in space) of the displacement and the L ∞ t (L 2 x )-norm of the pressure and on the other hand the L 2 t (H 1 x )-norm of the pressure. This key feature implies that error estimates with optimal convergence orders in space require the use of different polynomial degrees in the finite element spaces for the displacement and for the pressure, namely one degree higher for the displacement than for the pressure. Note that the use of different polynomial degrees for approximating the displacement and the pressure is motivated here solely by the derivation of optimal error estimates, as opposed to Biot's consolidation problem where the approximation of the initial data by solving a Stokes problem plays also a role. The technique of Wheeler [START_REF] Wheeler | A priori L 2 error estimates for Galerkin approximations to parabolic partial differential equations[END_REF] originally designed to obtain optimal L ∞ t (L 2 x ) a priori error estimates for the heat equation can be adapted to fit the present framework to derive optimal error estimates for the displacement field. The same technique has already been used in [START_REF] Murad | Improved accuracy in finite element analysis of Biot's consolidation problem[END_REF][START_REF] Murad | On stability and convergence of finite element approximations of Biot's consolidation problem[END_REF][START_REF] Murad | Asymptotic behavior of semidiscrete finite-element approximations of Biot's consolidation problem[END_REF] for Biot's consolidation problem.

The a posteriori error analysis of evolution problems related to poroelasticity is a much less explored field. The first aim is to derive error estimates that are reliable in the sense that they yield an upper bound for the approximation error (the difference between the exact solution and the discrete solution) that is fully computable in terms of known quantities, i.e. the problem data, universal constants from polynomial interpolation theory, and the discrete solution itself. A further important feature of the error estimates is their so-called optimality in the sense that the estimates are bounded from above by the a priori error estimates. A third desirable feature of the estimates is that they provide an efficient tool to adapt the discretization parameters (mesh size and time step). In the present work, we derive fully reliable and partially optimal a posteriori error estimates of residual type for the poroelasticity system (1)-( 2). The analysis hinges on the stability properties of the continuous problem, so that the natural norms with which to control the approximation error are the L ∞ t (H 1 x )-norm for the displacement and the L ∞ t (L 2 x )-and L 2 t (H 1 x )-norms for the pressure. The direct approach yields error estimates that have a form analogous to those derived for the heat equation by Picasso [START_REF] Picasso | Adaptive finite elements for a linear parabolic problem[END_REF], Chen and Feng [START_REF] Chen | An adaptive finite element algorithm with reliable and efficient error control for linear parabolic problems[END_REF] and Bergam, Bernardi and Mghazli [START_REF] Bergam | A posteriori analysis of the finite element discretization of some parabolic equations[END_REF]. The error estimate consists of three terms, a time error indicator (evaluated from the pressure differences at two consecutive time steps), a space error indicator (evaluated from the finite element residuals for the displacement and the pressure) and a data oscillation term. Despite these similarities with the heat equation, there are however important differences. Firstly, since stability is achieved by testing the momentum balance against the time-derivative of the displacement, the space error indicator contains a term involving the finite element residuals of the time-derivative of the displacement. This also makes the analysis of time-dependent meshes much more cumbersome than for the heat equation. Therefore, we have chosen to restrict ourselves to fixed meshes and to postpone the study of discretizations with time-dependent meshes to future work. The second important difference is that the estimates derived within the direct approach are not optimal for the approximation error on the displacement when the discretization of the latter involves higher-degree polynomials than for the pressure. To tackle this difficulty, we draw inspiration from the elliptic reconstruction technique introduced for linear parabolic problems by Makridakis and Nochetto [START_REF] Ch | Ellitpic reconstruction and a posteriori error estimates for elliptic problems[END_REF] and further analyzed by Lakkis and Makridakis [START_REF] Lakkis | Elliptic reconstruction and a posteriori error estimates for fully discrete linear parabolic problems[END_REF]. This technique, which can be regarded as the counterpart of the elliptic projection method introduced by Wheeler for the a priori error analysis, is designed to obtain optimal a posteriori error estimates in the L ∞ t (L 2 x )-norm (and other higher-order norms) for linear parabolic problems. In the present work, we extend this technique to coupled elliptic-parabolic problems such as the poroelasticity equations to derive a posteriori error estimates for the displacement that exhibit optimal convergence behavior. Other approaches to derive L ∞ t (L 2 x ) a posteriori error estimates for parabolic equations can be found, among others, in the work of Eriksson and Johnson [6,[START_REF] Eriksson | Adaptive finite element methods for parabolic problems. II. Optimal error estimates in l∞l 2 and l∞l∞[END_REF] and of Thomée [START_REF] Thomée | Galerkin Finite Element Methods for Parabolic Problems[END_REF] based on duality techniques, the work of Babuška, Feistauer, and Šolín [1] using a double integration in time, and the work of Verfürth [START_REF] Verfürth | A posteriori error estimates for nonlinear problems. L r (0, T ; L ρ (ω))-error estimates for finite element discretizations of parabolic equations[END_REF] for quasi-linear parabolic equations.

This paper is organized as follows. §2 presents the setting under scrutiny in an abstract framework. This setting allows us to pinpoint the mathematical structure of the poroelasticity system that plays a relevant role in the subsequent analysis. §3 is devoted to the a priori error analysis. The main result is Theorem 3.1 whose proof relies on the technique of elliptic projection introduced by Wheeler. The arguments are similar to those considered in [START_REF] Murad | Improved accuracy in finite element analysis of Biot's consolidation problem[END_REF][START_REF] Murad | On stability and convergence of finite element approximations of Biot's consolidation problem[END_REF][START_REF] Murad | Asymptotic behavior of semidiscrete finite-element approximations of Biot's consolidation problem[END_REF] for Biot's consolidation problem, and are briefly presented here for completeness. §4 deals with the a posteriori error analysis. The first result, Theorem 4.1, relies on the stability of the continuous problem and is established by proceeding in a way similar to the heat equation. It yields a reliable error estimate for the displacement in the L ∞ t (H 1 x )-norm and for the pressure in the L ∞ t (L 2 x )-and L 2 t (H 1 x )-norms. The error upper bound consists of data, time, and space error indicators. The time error indicator is optimal while the space error indicator comprises three terms, two of which yield local lower bounds of the approximation error. The last one exhibits an optimal convergence behavior with respect to mesh size when compared to the pressure error in the L 2 t (H 1 x )-norm, but its optimality with respect to the time step cannot be established because it depends on the time-derivative of the discrete displacement field. The second result, Theorem 4.2, yields reliable error estimates for the displacement error in the L ∞ t (H 1 x )-norm (and also for the pressure in the L ∞ t (L 2 x )-norm) and these estimates exhibit the optimal convergence order with respect to mesh size. The proof relies on an adaptation of the technique of elliptic reconstruction introduced recently by Makridakis and Nochetto. §5 contains numerical results illustrating the performance of the various a posteriori error estimators. Finally, §6 draws some conclusions.

The setting

The continuous problem

Let V a and V d be two Hilbert spaces respectively equipped with symmetric, continuous and coercive bilinear forms a and d. The norms induced by these forms are denoted by • a and • d respectively. Let V a (resp., V d ) be the dual space of V a (resp., V d ) with duality product denoted by •, • a (resp., •, • d ) and norm

• a = sup 0 =v∈Va | •, v a |/ v a (resp., • d = sup 0 =q∈V d | •, q d |/ q d ).
Let L a (resp., L d ) be a Hilbert space equipped with a scalar product (•, •) La (resp., (•, •) L d ) with dense and continuous injection V a → L a (resp.,

V d → L d ).
Identifying L a (resp., L d ) with its dual space, the following injections hold:

V a → L a ≡ (L a ) → V a (resp., V d → L d ≡ (L d ) → V d ).
Let c be a symmetric, continuous and coercive bilinear form defined over L d × L d inducing a norm • c such that for all q ∈ V d , q c ≤ γ q d and q L d ≤ γ q d . Finally, let b be a continuous bilinear form defined over V a × L d with continuity constant β, i.e., for all (v, q)

∈ V a × L d , |b(v, q)| ≤ β v a q c .
The elements of the spaces defined above are functions of the space variable x. In the sequel, we shall deal with functions of time and space. The time variable varies over the interval [0, T ] for a fixed T > 0. Henceforth, L 2 t (Z) denotes the vector space of functions f in space and time such that for a.e.

(V a ) × H 1 t (V d ) such that for all t ∈ [0, T ], a(u, v) -b(v, p) = f, v a , ∀v ∈ V a , (3)
c(∂ t p, q) + b(∂ t u, q) + d(p, q) = g, q d , ∀q ∈ V d , (4)
completed with the initial condition p(0) := p 0 . Note that in the present setting, Equ. (3) holds up to t = 0, thus uniquely determining the initial value of u in terms of p 0 and f (0). Letting u 0 := u(0) ∈ V a , the a priori bound u 0 a ≤ β p 0 c + f (0) a is readily inferred by taking v := u 0 in (3). Our first result is an a priori estimate for the strong solution. This will allow us to identify the natural stability norm for problem (3)-(4). Proposition 2.1. In the above setting, the following holds for a.e. t ∈ [0, T ],

1 2 u(t) 2 a + 1 2 p(t) 2 c + 1 2 t 0 p(s) 2 d ds ≤ 2 2 sup s∈[0,T ] f (s) a + T 0 ∂ t f (s) a ds 2 + T 0 g(s) 2 d ds + u 0 2 a + p 0 2 c . (5)
Proof. Since (u, p) is a strong solution, for a.e.

t ∈ [0, T ], v = ∂ t u is in V a and q = p is in V d . Using these test functions in (3)-(4) yields 1 2 d t u 2 a + 1 2 d t p 2 c + p 2 d = f, ∂ t u a + g, p d . Hence, 1 2 d t u 2 a + 1 2 d t p 2 c + 1 2 p 2 d ≤ f, ∂ t u a + 1 2 g 2 d .
Let t ∈ (0, T ). Integrating the above inequality over (0, t) and integrating by parts in time the term f, ∂ t u a yields

1 2 u(t) 2 a + 1 2 p(t) 2 c + 1 2 t 0 p(s) 2 d ds ≤ f (t), u(t) a -f (0), u(0) a - t 0 ∂ t f (s), u(s) a ds + 1 2 t 0 g(s) 2 d ds + 1 2 u 0 2 a + 1 2 p 0 2 c .
As a result, it is inferred that for a.e. t ∈ [0, T ],

1 2 u(t) 2 a + 1 2 p(t) 2 c + 1 2 t 0 p(s) 2 d ds ≤ Aσ T (u) + B 2 , with σ T (u) = sup s∈[0,T ] u(s) a , A = 2 sup s∈[0,T ] f (s) a + T 0 ∂ t f (s) a ds and B 2 = 1 2 T 0 g(s) 2 d ds + 1 2 u 0 2 a + 1 2 p 0 2 c . Hence, 1 2 σ T (u) 2 ≤ Aσ T (u) + B 2 so that σ T (u) 2 ≤ 4(A 2 + B 2
). This yields

1 2 u(t) 2 a + 1 2 p(t) 2 c + 1 2 t 0 p(s) 2 d ds ≤ A 2 + 1 4 σ T (u) 2 + B 2 ≤ 2(A 2 + B 2 ), yielding (5)
.

An important consequence of Proposition 2.1 is the uniqueness of the strong solution of (3)-( 4). In the sequel, we assume the existence of the strong solution.

Remark 2.1. Proposition 2.1 holds in the slightly more general setting where

g ∈ L 2 t (V d ), p 0 ∈ L d , and p ∈ L 2 t (V d ) ∩ H 1 t (L d ).
Application to poroelasticity. The evolution problem (1)-( 2) fits the present setting. For the sake of simplicity, we consider homogeneous Dirichlet boundary conditions both for the displacement and the pressure. Then, letting

V a = [H 1 0 (Ω)] 3 , L a = [L 2 (Ω)] 3 , V d = H 1 0 (Ω) and L d = L 2 (Ω), we define the bilinear forms a(u, v) = Ω σ(u):ε(v), b(v, p) = Ω bp∇•v, (6)
c(p, q) = Ω 1 M pq, d(p, q) = Ω κ∇p•∇q. (7)
The coercivity of a on V a × V a (resp., d on V d × V d ) results from Korn's First Inequality (resp., Poincaré's Inequality). The bilinear form b is clearly continuous on V a × L d with β = b(M/E) 1/2 where the Young modulus E is associated with the coercivity of the bilinear form a on V a . Proposition 2.1 means that the natural stability norm controls the displacement in the L ∞ t (H 1 x )-norm and the pressure in the L ∞ t (L 2 x )-and L 2 t (H 1 x )-norms. This stability result hinges on the fact that the momentum balance can be tested by the time-derivative of the displacement.

The discrete problem

Problem (3)-( 4) is approximated by an Euler-Galerkin scheme, namely conforming finite elements in space and an implicit Euler scheme in time. Let {V ah } h>0 and {V dh } h>0 be two families of finite-dimensional subspaces of V a and V d respectively. The parameter h refers to the size of an underlying mesh family denoted by {T h } h>0 . Let 0 = t 0 < t 1 < . . . < t N = T be a sequence of discrete times and for all n ∈ {1, . . . , N }, set τ n = t n -t n-1 and I n = (t n-1 , t n ). Henceforth, a superscript n indicates the value taken by any function of space and time at the discrete time t n . For instance,

u n := u(t n ) ∈ V a .
The discrete problem consists of seeking

{u n h } N n=1 ∈ [V ah ] N and {p n h } N n=1 ∈ [V dh ] N such that for all n ∈ {1, . . . , N }, a(u n h , v h ) -b(v h , p n h ) = (f n h , v h ) La , ∀v h ∈ V ah , (8)
c(δ t p n h , q h ) + b(δ t u n h , q h ) + d(p n h , q h ) = (g n h , q h ) L d , ∀q h ∈ V dh , (9)
where

δ t p n h = τ -1 n (p n h -p n-1 h ) and δ t u n h = τ -1 n (u n h -u n-1 h ). Given a pair (u 0h , p 0h ) ∈ V ah × V dh , the initial condition is (u 0 h , p 0 h ) := (u 0h , p 0h ). The data {f n h } N n=1 ∈ [V ah ] N ⊂ [V a ] N and {g n h } N n=1 ∈ [V dh ] N ⊂ [V d ] N are approximations of {f n } N n=1 and {g n } N n=1 respectively.
Lemma 2.1. The discrete problem is well-posed.

Proof. For all n ∈ {1, . . . , N }, Equs. ( 8)-( 9) yield a square linear system for the components of (u n h , p n h ) once bases of V ah and V dh are chosen, so it suffices to prove the uniqueness of the discrete solution. Testing with

v h = u n h -u n-1 h
and q h = τ n p n h and using the fact that a(x, x -y) = 1 2 a(x, x) + 1 2 a(x -y, x -y) -1 2 a(y, y) owing to the symmetry of the bilinear form a (along with a similar property for the bilinear form c) yields

1 2 u n h 2 a + 1 2 u n h -u n-1 h 2 a + 1 2 p n h 2 c + 1 2 p n h -p n-1 h 2 c + τ n p n h 2 d = 1 2 u n-1 h 2 a + (f n h , u n h -u n-1 h ) La + 1 2 p n-1 h 2 c + τ n (g n h , p n h ) L d .
Hence,

1 2 u n h 2 a + 1 2 p n h 2 c + 1 2 τ n p n h 2 d ≤ 1 2 u n-1 h 2 a + 1 2 f n h 2 a + 1 2 p n-1 h 2 c + 1 2 τ n g n h 2 d .
This shows the uniqueness of the solution of the square linear system.

The proof of Lemma 2.1 hints at how the discrete scheme ( 8)-( 9) could be modified if time-dependent meshes were used. In this case, we work with two families {V n ah } N n=0 and {V n dh } N n=0 of finite-dimensional subspaces such that for all n ∈ {0, . . . , N }, V n ah ⊂ V a and V n dh ⊂ V d . The discrete scheme takes the general form ( 8)-( 9) with test functions v h ∈ V n ah and q h ∈ V n dh . However, if the expression for the time-derivative of the displacement is kept unchanged, the argument deployed in the above proof breaks down because it is no longer possible to use

v h = u n h -u n-1 h as a test function since in general v h ∈ V n ah (unless the restrictive assumption V n-1 ah ⊂ V n ah is made for all n ∈ {1, . . . , N }). To circumvent this difficulty, let R n * ah : V a → V n ah be the Riesz projection operator defined such that for all v ∈ V a , a(v -R n * ah (v), v h ) = 0, ∀v h ∈ V n ah , (10)
and use

δ t u n h = τ -1 n (u n h -R n * ah (u n-1 h
)) in [START_REF] Lakkis | Elliptic reconstruction and a posteriori error estimates for fully discrete linear parabolic problems[END_REF]. Then, proceeding as in the proof of Lemma 2.1 with the test functions

v h = u n h -R n * ah (u n-1 h ) ∈ V n ah and q h = τ n p n h ∈ V n dh and observing that a(u n h , v h ) = a(u n h , u n h -u n-1 h ) since u n h ∈ V n
ah , the same stability estimate is recovered.

Continuous and discrete differential operators

To formulate the a posteriori error estimates in the usual form, it is convenient to associate differential operators (in space) with the bilinear forms a, b, c, and d. To this purpose, we define the following continuous operators:

A ∈ L(V a ; V a ) s.t. Av, w a = -a(v, w), B ∈ L(V a ; L d ) s.t. (Bv, q) L d = b(v, q) with adjoint B * ∈ L(L d ; V a ) s.t. B * q, v a = b(v, q), C ∈ L(L d ; L d ) s.t. (Cq, r) L d = c(q, r), and D ∈ L(V d ; V d ) s.t. Dq, r d = -d(q, r).
Observe that C is actually a zero-order operator. With this notation, problem (3)-( 4) can be rewritten in the form

-Au -B * p = f, (11)
C∂ t p + B∂ t u -Dp = g, (12)
these equalities holding for a.e. t ∈ [0, T ] in V a and V d respectively. For the poroelasticity system, Av = ∇•σ(v), Bv = b∇•v, B * q = -b∇q, Cq = 1 M q, and Dq = ∇•(κ∇p). The discrete version of the above operators will also be used:

A h ∈ L(V ah ; V ah ) s.t. (A h v h , w h ) La = -a(v h , w h ), B h ∈ L(V ah ; V dh ) s.t. (B h v h , q h ) L d = b(v h , q h ) with adjoint B * h ∈ L(V dh ; V ah ) s.t. (B * h q h , v h ) La = b(v h , q h ), and D h ∈ L(V dh ; V dh ) s.t. (D h q h , r h ) L d = -d(q h ,
r h ). Observe that duality products have been replaced by L a -and L d -scalar products. The discrete problem ( 8)-( 9) can be rewritten in the form

-A h u n h -B * h p n h = f n h , (13)
Cδ t p n h + B h δ t u n h -D h p n h = g n h , (14)
these equalities holding for all n ∈ {1, . . . , N } in V ah and V dh respectively. For later use, we let f 0 h := -A h u 0h -B * h p 0h , so that (13) also holds for n = 0.

The steady problem

Both the a priori and a posteriori error analysis of the approximation of the steady version of ( 3)-( 4) using the subspaces {V ah } h>0 and {V dh } h>0 will play a role in the error analysis for the time-dependent case. The steady version of ( 3)-( 4) consists of seeking u ∈ V a and p ∈ V d such that

a(u, v) -b(v, p) = f , v a , ∀v ∈ V a , (15)
d(p, q) = g, q d , ∀q ∈ V d , (16)
with data f ∈ V a and g ∈ V d . It is straightforward to verify that the problem ( 15)-( 16) is well-posed owing to its upper triangular structure and the coercivity of the bilinear forms a and d.

The discrete problem consists of seeking u h ∈ V ah and p h ∈ V dh such that

a(u h , v h ) -b(v h , p h ) = f , v h a , ∀v h ∈ V ah , (17)
d(p h , q h ) = g, q h d , ∀q h ∈ V dh . (18)
Here, we do not consider an approximation to the data f and g. The discrete problem is conveniently reformulated using a Riesz projection operator R h :

V a × V d → V ah × V dh such that for all (v, q) ∈ V a × V d , R h (v, q) := (R ah (v, q), R dh (q)) is defined by a(v -R ah (v, q), v h ) -b(v h , q -R dh (q)) = 0, ∀v h ∈ V ah , (19)
d(q -R dh (q), q h ) = 0, ∀q h ∈ V dh . (20)
It is clear that (u h , p h ) solves ( 17)-( 18) if and only if u h = R ah (u, p) and p h = R dh (p). The approximation properties of the operator R h can be found in [START_REF] Murad | Improved accuracy in finite element analysis of Biot's consolidation problem[END_REF]. The result is restated here for completeness.

Lemma 2.2. The following holds for all

(v, q) ∈ V a × V d , v -R ah (v, q) a ≤ inf v h ∈V ah v -v h a + β q -R dh (q) c , (21)
q -R dh (q) d = inf

q h ∈V dh q -q h d . (22)
Proof. Property ( 22) is classical. To establish [START_REF] Verfürth | A posteriori error estimates for nonlinear problems. L r (0, T ; L ρ (ω))-error estimates for finite element discretizations of parabolic equations[END_REF], consider the operator R * ah defined by (10) (the upper index n is dropped since meshes are kept fixed in time). Then, observe that since both R * ah (v) and R ah (v, q) are in

V ah , R ah (v, q) -R * ah (v) 2 a = a(R ah (v, q) -R * ah (v), R ah (v, q) -R * ah (v)) = a(R ah (v, q) -v, R ah (v, q) -R * ah (v)) + a(v -R * ah (v), R ah (v, q) -R * ah (v)) = -b(R ah (v, q) -R * ah (v), q -R dh (q)) ≤ β R ah (v, q) -R * ah (v) a q -R dh (q) c .
Hence, R ah (v, q) -R * ah (v) a ≤ β q -R dh (q) c whence it follows by the triangle inequality that

v -R ah (v, q) a ≤ v -R * ah (v) a + R ah (v, q) -R * ah (v) a ≤ v -R * ah (v) a + β q -R dh (q) c ,
readily yielding [START_REF] Verfürth | A posteriori error estimates for nonlinear problems. L r (0, T ; L ρ (ω))-error estimates for finite element discretizations of parabolic equations[END_REF].

To deduce from Lemma 2.2 asymptotic rates of convergence for the approximation error in terms of the parameter h when the exact solution is smooth enough, we introduce the following assumptions.

Hypothesis 2.1. There exist constants c 1 and c 2 , positive real numbers s a and s d , and subspaces W a ⊂ V a and

W d ⊂ V d respectively equipped with norms • Wa and • W d , such that independently of h, ∀v ∈ W a , inf v h ∈V ah v -v h a ≤ c 1 h sa v Wa , (23)
∀q ∈ W d , inf

q h ∈V dh q -q h d ≤ c 2 h s d q W d . ( 24
)
Hypothesis 2.2. There exist a constant c 3 and a positive real number δ such that for all r ∈ L d , the unique solution φ ∈ V d of the dual problem d(q, φ) = c(r, q) for all q ∈ V d , is such that there is

φ h ∈ V dh satisfying φ -φ h d ≤ c 3 h δ r c . (25)
Hypothesis 2.1 is classical in the context of finite element approximations. It will be used in the a priori error analysis. To keep technicalities at a minimum, a version of Hypothesis 2.1 localized to mesh cells is not considered. Hypothesis 2.2 is an elliptic regularity property associated with the bilinear form d on V d . It is stated here in compact form, the usual statement consisting of assuming that the dual solution φ is in a subspace Y d of V d where the interpolation estimate (25) holds in the form φ -φ h d ≤ c 3 h δ φ Y d . Hypothesis 2.2 will serve both in the a priori and the a posteriori error analysis. In the latter case, a sharper statement localized to mesh cells will be introduced in §4.2. For the time being, we will only use the following important consequence of Hypothesis 2.2:

q -R dh (q) c ≤ c 3 h δ q -R dh (q) d .

(26) Indeed, letting φ be the dual solution associated with r := q -R dh (q) and observing that d(r, φ h ) = 0 for

φ h ∈ V dh yields q -R dh (q) 2 c = c(r, r) = d(r, φ) = d(r, φ -φ h ) ≤ r d φ -φ h d , (27)
whence (26) readily follows. An important consequence of ( 21), [START_REF] Wheeler | A priori L 2 error estimates for Galerkin approximations to parabolic partial differential equations[END_REF], and (26) is that for all

(v, q) ∈ W a × W d , v -R ah (v, q) a ≤ c 1 h sa v Wa + βc 2 c 3 h s d +δ q W d . (28)
For the purpose of computational efficiency, it is reasonable to balance both sources of error in v -R ah (v, q) a . This motivates the following hypothesis.

Hypothesis 2.3. s a = s d + δ =: s.

In the framework of Hypotheses 2.1-2.3, Lemma 2.2 yields for all (v, q)

∈ W a × W d , v -R ah (v, q) a ≤ h s (c 1 v Wa + βc 2 c 3 q W d ), (29)
q -R dh (q) c ≤ h s c 2 c 3 q W d , (30)
q -R dh (q) d ≤ h s-δ c 2 q W d . (31)
As a result, whenever the exact solution of the steady problem ( 15)-( 16) is smooth enough, namely (u, p) ∈ W a × W d , the error p -p h d converges asymptotically as h s-δ while the error u -u h a + p -p h c converges asymptotically as h s . Since δ is positive, this means that the error u -u h a + p -p h c converges at a faster rate than p -p h d . This difference in the convergence rates will be accounted for in the subsequent analysis of the time-dependent problem, the goal being to derive a priori and a posteriori error bounds that are optimal for p n -p n h d on the one hand and for u n -u n h a + p n -p n h c on the other hand. Application to poroelasticity. Consider the model problem ( 1)-( 2) with the displacement (resp., the pressure) approximated in space by continuous Lagrange finite elements of degree k ≥ 1 (resp., l ≥ 1). Then, Hypothesis (2.1) holds with

s a := k, W a := [H 1 0 (Ω) ∩ H k+1 (T h )] 3 , s d := l and W d := H 1 0 (Ω) ∩ H l+1 (T h )
, where for m ≥ 0, H m (T h ) denotes the usual broken Sobolev space of order m. Hypothesis 2.2 means that the steady-state version of the pressure equation yields elliptic regularity, namely for all r ∈ L 2 (Ω), the unique solution φ ∈ H 1 0 (Ω) to the dual problem Ω κ∇φ•∇q = Ω rq for all q ∈ H 1 0 (Ω) is in H 2 (Ω). Then, (25) holds with δ := 1. As a result, Hypothesis 2.3 implies

k = l + 1, (32)
i.e. the polynomial interpolation for the displacement is one degree higher than that for the pressure. The most common choice in practice is k := 2 and l := 1, i.e. continuous piecewise quadratics are used to approximate the displacement and continuous piecewise linears are used to approximate the pressure.

A priori error analysis

The a priori error analysis is performed under the assumption that the exact solution is smooth, namely

u ∈ C 1 t (W a ) ∩ C 2 t (V a ), p ∈ C 1 t (W d ) ∩ C 2 t (L d ). ( 33
)
For all n ∈ {1, . . . , N }, define

C n 1 (u, p) = 2γ 2 c 2 2 c 2 3 sup s∈In ∂ t p(s) 2 W d + 2β 2 γ 2 (c 1 sup s∈In ∂ t u(s) Wa + βc 2 c 3 sup s∈In ∂ t p(s) W d ) 2 , (34)
C n 2 (u, p) = 1 2 γ 2 sup s∈In ∂ 2 tt p(s) 2 c + 1 2 β 2 γ 2 sup s∈In ∂ 2 tt u(s) 2 a , (35)
C n (f, g) = 1 2 f n -f n h 2 a + τ n g n -g n h 2 d . (36)
Moreover, it is assumed that the initial data u 0h and p 0h are chosen such that

u 0 -u 0h a ≤ c 4 h s u 0 Wa and p 0 -p 0h c ≤ c 5 h s p 0 W d , (37)
and we define

C(u 0 , p 0 ) = (c 1 + c 4 ) 2 u 0 2 Wa + (β 2 c 2 2 c 2 3 + 1 2 (c 2 c 3 + c 5 ) 2 ) p 0 2 W d . (38
) One possible choice is u 0h = R ah (u 0 , p 0 ) and p 0h = R dh (p 0 ), in which case we can take C(u 0 , p 0 ) = 0. Theorem 3.1. In the above framework, the following holds for all n ∈ {1, . . . , N },

1 4 u n -u n h 2 a + 1 4 p n -p n h 2 c ≤ h 2s C(u 0 , p 0 ) + n m=1 C m (f, g) + n m=1 [τ m h 2s C m 1 (u, p) + τ 3 m C m 2 (u, p)] + h 2s (c 2 1 u n 2 Wa + (β 2 + 1 2 )c 2 2 c 2 3 p n 2 W d ), (39)
and

n m=1 1 8 τ m p m -p m h 2 d ≤ h 2s C(u 0 , p 0 ) + n m=1 C m (f, g) + n m=1 [τ m h 2s C m 1 (u, p) + τ 3 m C m 2 (u, p)] + n m=1 1 4 τ m h 2(s-δ) c 2 2 p m 2 W d . (40)
Proof. (i) For all n ∈ {1, . . . , N }, let us first estimate the quantities

η n ah = R ah (u n , p n ) -u n h and η n dh = R dh (p n ) -p n h .
Observe that

a(η n ah , v h ) -b(v h , η n dh ) = f n -f n h , v h a , ∀v h ∈ V ah , c(η n dh -η n-1 dh , q h ) + b(η n ah -η n-1 ah , q h ) + τ n d(η n dh , q h ) = τ n g n -g n h , q h d + c(θ n dh , q h ) + b(θ n ah , q h ), ∀q h ∈ V dh ,
where

θ n dh = R dh (p n ) -R dh (p n-1 ) -τ n ∂ t p n and θ n ah = R ah (u n , p n ) -R ah (u n-1 , p n-1 ) -τ n ∂ t u n .
Testing with v h := η n ah -η n-1 ah ∈ V ah and q h := η n dh ∈ V dh yields after some straightforward algebra

1 2 η n ah 2 a + 1 2 η n ah -η n-1 ah 2 a + 1 2 η n dh 2 c + 1 2 η n dh -η n-1 dh 2 c + τ n η n dh 2 d = 1 2 η n-1 ah 2 a + 1 2 η n-1 dh 2 c + f n -f n h , η n ah -η n-1 ah a + τ n g n -g n h , η n dh d + c(θ n dh , η n dh ) + b(θ n ah , η n dh ).
Hence,

1 2 η n ah 2 a + 1 2 η n dh 2 c + 1 4 τ n η n dh 2 d ≤ 1 2 η n-1 ah 2 a + 1 2 η n-1 dh 2 c + C n (f, g) + τ -1 n γ 2 θ n dh 2 c + τ -1 n β 2 γ 2 θ n ah 2 a .
(ii) Let us now estimate the quantities θ n dh and θ n ah . Observe that

θ n dh = - In [∂ t p(s) -R dh (∂ t p(s))]ds - In (s -t n-1 )∂ 2 tt p(s)ds.
Hence, owing to the regularity assumptions on the exact solution and estimate (30),

θ n dh c ≤ τ n h s c 2 c 3 sup s∈In ∂ t p(s) W d + 1 2 τ 2 n sup s∈In ∂ 2 tt p(s) c .
Similarly, using (29),

θ n ah a ≤ τ n h s (c 1 sup s∈In ∂ t u(s) Wa + βc 2 c 3 sup s∈In ∂ t p(s) W d ) + 1 2 τ 2 n sup s∈In ∂ 2 tt u(s) a . Therefore, τ -1 n γ 2 θ n dh 2 c + τ -1 n β 2 γ 2 θ n ah 2 a ≤ τ n h 2s C n 1 (u, p) + τ 3 n C n 2 (u, p).

Summing up the above estimates leads to

1 2 η n ah 2 a + 1 2 η n dh 2 c + n m=1 1 4 τ m η m dh 2 d ≤ 1 2 η 0 ah 2 a + 1 2 η 0 dh 2 c + n m=1 [C m (f, g) + τ m h 2s C m 1 (u, p) + τ 3 m C m 2 (u, p)].
(iii) We now estimate the initial errors η 0 ah and η 0 dh . In the case where u 0h = R ah (u 0 , p 0 ) and p 0h = R dh (p 0 ), it is clear that η 0 ah = 0 and η 0 dh = 0. In the general case, use the triangle inequality to infer

η 0 dh c = R dh (p 0 ) -p 0h c ≤ R dh (p 0 ) -p 0 c + p 0 -p 0h c ≤ h s (c 2 c 3 + c 5 ) p 0 W d .
Similarly,

η 0 ah a ≤ R ah (u 0 , p 0 ) -u 0 a + u 0 -u 0h a ≤ h s ((c 1 + c 4 ) u 0 Wa + βc 2 c 3 p 0 W d ).
Hence,

1 2 η n ah 2 a + 1 2 η n dh 2 c + n m=1 1 4 τ m η m dh 2 d ≤ h 2s C(u 0 , p 0 ) + n m=1 [C m (f, g) + τ m h 2s C m 1 (u, p) + τ 3 m C m 2 (u, p)].
(iv) The conclusion readily results from the triangle inequality and estimates (29)-( 30)-(31).

Theorem 3.1 shows that whenever the exact solution is smooth enough and up to data approximation errors that can be made small enough, u n -u n h a + p n -p n h c converges to order s in space and first-order in time, while (

n m=1 τ m p m -p m h 2 d ) 1/2
converges to order (s -δ) in space and first-order in time.

Application to poroelasticity. When continuous piecewise quadratics (resp., linears) are used to approximate the displacement (resp., the pressure), u n -u n h H 1 + p n -p n h L 2 converges to second-order in space and firstorder in time, while (

n m=1 τ m p m -p m h 2
H 1 ) 1/2 converges to first-order in space and in time.

A posteriori error analysis

This section is devoted to the a posteriori error analysis for the discrete scheme ( 8)-( 9). Two estimates are derived. The first directly relies on the stability of the continuous problem and yields an error estimate that is suitable to measure the error in p in the L 2 t (V d )-norm. The second is based on an adaption of the elliptic reconstruction technique of Makridakis and Nochetto and yields an error estimate that is suitable to measure the error in u in the L ∞ t (V a )-norm and the error in p in the L ∞ t (L d )-norm. As is customary in a posteriori error analysis, we assume in this section that the data (f, g) are in L 2 t (L a ) × L 2 t (L d ).

The direct approach

The a posteriori error analysis relies on the stability of the continuous problem. Therefore, we rewrite the discrete scheme as equations holding a.e. in (0, T ) rather than at the discrete times {t n } N n=1 . To this purpose, let u hτ (resp., p hτ ) be the continuous and piecewise affine function in time such that for all n ∈ {0, . . . , N }, u hτ (t n ) = u n h (resp., p hτ (t n ) = p n h ). Observe that ∂ t u hτ and ∂ t p hτ are defined a.e. in (0, T ). Similarly, let f hτ be the continuous and piecewise affine function in time such that for all n ∈ {0, . . . , N }, f hτ (t n ) = f n h . We will also need to consider piecewise constant functions in time, namely π 0 p hτ (resp., π 0 g hτ ) equal to p n h (resp., g n h ) on I n for all n ∈ {1, . . . , N }. With the above notation, the discrete scheme ( 8)-( 9) yields a.e. in (0, T ),

a(u hτ , v h ) -b(v h , p hτ ) = (f hτ , v h ) La , ∀v h ∈ V ah , (41)
c(∂ t p hτ , q h ) + b(∂ t u hτ , q h ) + d(π 0 p hτ , q h ) = (π 0 g hτ , q h ) L d , ∀q h ∈ V dh . (42)
To formulate the a posteriori error estimate, it is convenient to introduce the Galerkin residual G a (resp., G d ) which is a continuous and piecewise affine function in time with values in V a (resp., piecewise constant function in time with values in V d ) such that a.e. in (0, T ),

G a , v a = (f hτ , v) La -a(u hτ , v) + b(v, p hτ ), ∀v ∈ V a , (43)
G d , q d = (π 0 g hτ , q) L d -c(∂ t p hτ , q) -b(∂ t u hτ , q) -d(π 0 p hτ , q), ∀q ∈ V d . ( 44
)
Define the data, space, and time error estimators

E(f, g) = T 0 (g -π 0 g hτ )(s) 2 d ds + 2 sup s∈[0,T ] (f -f hτ )(s) a + T 0 ∂ t (f -f hτ )(s) a ds 2 , ( 45
)
E dat = u 0 -u 0h 2 a + p 0 -p 0h 2 c + 4E(f, g), (46)
E spc = N m=1 4τ m G m d 2 d + 4 2 sup 0≤m≤N G m a a + N m=1 G m a -G m-1 a a 2 , (47)
E tim = N m=1 1 3 τ m p m h -p m-1 h 2 d . ( 48
)
Theorem 4.1. For all n ∈ {1, . . . , N },

1 2 u n -u n h 2 a + 1 2 p n -p n h 2 c + tn 0 1 4 (p -p hτ )(s) 2 d ds + tn 0 1 2 (p -π 0 p hτ )(s) 2 d ds ≤ E dat + E spc + E tim . ( 49
)
Proof. Let ξ = u -u hτ , ζ = p -p hτ , and ζ * = p -π 0 p hτ . Observe that a.e. in (0, T ),

a(ξ, v) -b(v, ζ) = f -f hτ + G a , v a , ∀v ∈ V a , c(∂ t ζ, q) + b(∂ t ξ, q) + d(ζ * , q) = g -π 0 g hτ + G d , q d , ∀q ∈ V d .
Testing for a.e. t ∈ (0, T ) with v := ∂ t ξ and q = ζ yields

1 2 d t ξ 2 a + 1 2 d t ζ 2 c + 1 2 ζ 2 d + 1 2 ζ * 2 d = f -f hτ + G a , ∂ t ξ a + g -π 0 g hτ + G d , ζ d + 1 2 p hτ -π 0 p hτ 2 d ,
where we have used the fact that owing to the symmetry of d,

d(ζ, ζ * ) = 1 2 d(ζ, ζ) + 1 2 d(ζ * , ζ * ) -1 2 d(ζ -ζ * , ζ -ζ * ). Since f -f hτ + G a ∈ H 1 t (V a ) and g -π 0 g hτ + G d ∈ L 2 t (V d ),
we can proceed as in the proof of Proposition 2.1 and infer for all n ∈ {1, . . . , N }, the bound (details are skipped for brevity)

1 2 u n -u n h 2 a + 1 2 p n -p n h 2 c + tn 0 1 4 (p -p hτ )(s) 2 d ds + tn 0 1 2 (p -π 0 p hτ )(s) 2 d ds ≤ E dat + T 0 4 G d (s) 2 d ds + 4 2 sup s∈[0,T ] G a (s) a + T 0 ∂ t G a (s) a 2 ds + T 0 (p hτ -π 0 p hτ )(s) 2 d ds.
The second and third terms in the right-hand side yield E spc and the last term yields E tim . The final bound (49) results from the fact that G a is piecewise affine, G d is piecewise constant, and that for

s ∈ I m , (p hτ -π 0 p hτ )(s) = τ -1 m (s -t m )(p m h -p m-1 h
).

Remark 4.1. The convergence rate of the upper bound in (49) is optimal for

tn 0 1 2 (p -p hτ )(s) 2 d ds, but not for 1 2 u n -u n h 2 a + 1 2 p n -p n h 2
c . The last term in the left-hand side of (49) is kept only to prove an optimality property of the time error estimator; see Proposition 4.2.

To localize the space error estimator E spc , we need to introduce some additional notation and assumptions. We assume that the various bilinear forms in the model problem ( 3)-( 4) can be localized as follows: for all

(v, q) ∈ V a × V d and for all (ξ, ζ) ∈ V a × V d such that either ξ ∈ V ah or Aξ ∈ L a and either ζ ∈ V dh or Dζ ∈ L d , a(ξ, v) = T ∈T h [-(Aξ, v) La(T ) + (J a ξ, v) La(∂T ) ],
(50)

d(ζ, q) = T ∈T h [-(Dζ, q) L d (T ) + (J d ζ, q) L d (∂T ) ]. (51)
Here, for a mesh cell T ∈ T h , L a (T ) and L d (T ) are local versions of L a and L d respectively, (•, •) La(∂T ) and (•, •) L d (∂T ) are scalar products for functions defined on the boundary ∂T of T and J a and J d are suitable (jump) operators such that

J a ξ = 0 if Aξ ∈ L a and J d ζ = 0 if Dζ ∈ L d . In addition, for all v ∈ V a and for all ζ ∈ L d such that either ζ ∈ V dh or B * ζ ∈ L a , b(v, ζ) = T ∈T h [(v, B * ζ) La(T ) + (J b ζ, v) La(∂T ) ], (52)
with

J b ζ = 0 if B * ζ ∈ L a .
Moreover, the bilinear forms b and c are also localized as follows: for all (v, q, r) ∈

V a × L d × L d , b(v, q) = T ∈T h (Bv, q) L d (T ) and c(q, r) = T ∈T h (Cq, r) L d (T ) . (53)
Finally, we assume that there exist two (Clément-type) interpolation operators i ah : V a → V ah and i dh :

V d → V dh such that for all (v, q) ∈ V a × V d , T ∈T h [h -2 T v -i ah (v) 2 La(T ) + h -1 T v -i ah (v) 2 La(∂T ) ] ≤ c 6 v 2 a , (54)
T ∈T h [h -2 T q -i dh (q) 2 L d (T ) + h -1 T q -i dh (q) 2 L d (∂T ) ] ≤ c 7 q 2 d , (55)
where for all T ∈ T h , h T denotes the diameter of T . This assumption is classical in the context of a posteriori error estimation; see, e.g. [START_REF] Ern | Theory and Practice of Finite Elements[END_REF][START_REF] Verfürth | A Review of a Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques[END_REF]. In the context of poroelasticity where V a = [H 1 0 (Ω)] 3 and V d = H 1 0 (Ω), the usual Clément [START_REF] Clément | Approximation by finite element functions using local regularization[END_REF] interpolation operator (modified to account for homogeneous Dirichlet boundary conditions) or the Scott-Zhang [START_REF] Scott | Finite element interpolation of nonsmooth functions satisfying boundary conditions[END_REF] interpolation operator can be used.

We define the following elementwise and jump residuals for all m ∈ {1, . . . , N },

R m uh = f m h + Au m h + B * p m h , J m uh = J a u m h -J b p m h , (56)
R m ph = g m h -Cδ t p m h -Bδ t u m h + Dp m h , J m ph = J d p m h . ( 57
)
For m = 0, R 0 uh , J 0 uh , and J 0 ph are defined similarly, while we set R 0 ph = (D -D h )p 0h for later use in §4.2. We also define for all m ∈ {0, . . . , N }, the space error estimators

E m u = T ∈T h [h 2 T R m uh 2 La(T ) + h T J m uh 2 La(∂T ) ], (58)
E m p,α = T ∈T h h 2α T [h 2 T R m ph 2 L d (T ) + h T J m ph 2 L d (∂T ) ], (59)
for a real parameter α ≥ 0. We will also use the time-incremental version of these estimators, namely for m ∈ {1, . . . , N },

E m u (δ t ) = T ∈T h τ 2 m [h 2 T δ t R m uh 2 La(T ) + h T δ t J m uh 2 La(∂T ) ], (60)
E m p,α (δ t ) = T ∈T h h 2α T τ 2 m [h 2 T δ t R m ph 2 L d (T ) + h T δ t J m ph 2 L d (∂T ) ], (61)
where

δ t R m uh = τ -1 m (R m uh -R m-1 uh ), δ t J m uh = τ -1 m (J m uh -J m-1 uh )
, and so on. Proposition 4.1. In the above framework, for all m ∈ {0, . . . , N },

G m a 2 a ≤ c 6 E m u , (62)
and for all m ∈ {1, . . . , N },

G m a -G m-1 a 2 a ≤ c 6 E m u (δ t ), (63)
G m d 2 d ≤ c 7 E m p,0 . (64)
Hence,

E spc ≤ N m=1 4c 7 τ m E m p,0 + 16c 6 sup 0≤m≤N E m u + 8c 6 N m=1 ( E m u (δ t )) 1/2 2 . (65)
Proof. The proof is only sketched since it uses classical techniques of a posteriori error analysis. Observe that

G m a a = sup 0 =v∈Va G m a , v a v a = sup 0 =v∈Va G m a , v -i ah (v) a v a , since (41)
and (43) imply that G m a , v h a = 0 for all v h ∈ V ah . Then, use definitions (50) and (52), estimate (54) and a Cauchy-Schwarz inequality to infer (62). The proof of (63) and (64) is similar. Finally, (65) results from the definition of E spc .

We now investigate the optimality of the space and time error estimators derived above, namely that these quantities yield lower bounds for the errors estimated in Equ. (49). The optimality of the time error estimator is straightforward. Proposition 4.2. The following holds

E tim ≤ 2 T 0 (p -p hτ )(s) 2 d ds + 2 T 0 (p -π 0 p hτ )(s) 2 d ds. (66)
Proof. Observe that E tim = T 0 (p hτ -π 0 p hτ )(s) 2 d ds and use the triangle inequality. To investigate the space error estimator, we consider the three terms in the right-hand side of (65). The first two can be bounded using the technique of bubble functions introduced by Verfürth [START_REF] Verfürth | A posteriori error estimations and adaptative mesh-refinement techniques[END_REF][START_REF] Verfürth | A Review of a Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques[END_REF]. Since this technique is well-known, the arguments below are only briefly sketched. To alleviate the notation, we denote by x y the inequality x ≤ cy with positive and mesh-independent constant c. We assume that the bilinear forms a, b, c, and d can be localized elementwise as a(v, w) = T ∈T h a T (v, w) and so on; this induces a localization of the norms in the form v 2 a = T ∈T h v 2 a,T with v 2 a,T = a T (v, v) and so on. For all T ∈ T h , the elementwise residuals R m uh and R m ph are in finite-dimensional spaces P a (T ) and P d (T ), respectively, while for a face F in the mesh, the restrictions of the jump residuals J m uh and J m ph to F are in finite-dimensional spaces P a (F ) and P d (F ), respectively. For all T ∈ T h , we assume that there is a bubble function ν a,T ∈ V a (resp., ν d,T ∈ V d ) with support localized in T such that for all v ∈ P a (T ), v 2

La(T ) (ν a,T v, v) La(T ) and ν a,T v La(T ) + h T ν a,T v a,T ≤ v La(T ) (resp., for all q ∈ P d (T ), q 2 L d (T ) (ν d,T q, q) L d (T ) and ν d,T q L d (T ) + h T ν d,T q d,T ≤ q L d (T )
). Moreover, for a face F ⊂ ∂T , we assume that there is a bubble function ν a,F ∈ V a (resp., ν d,F ∈ V d ) with support localized in F and a lifting operator π a,F (resp., π d,F ) from F to ∆ F (the set formed by the one or two elements of T h to which F belongs) such that for all v ∈ P a (F ), v 2

La(F ) (ν a,F v, v) La(F ) and h -1/2 T π a,F (ν a,F v) La(∆F ) + h 1/2 T π a,F (ν a,F v) a,∆F
v La(F ) (resp., for all q ∈ P d (F ),

q 2 L d (F ) (ν d,F q, q) L d (F ) and h -1/2 T π d,F (ν d,F q) L d (∆F ) + h 1/2 T π d,F (ν d,F q) d,∆F q L d (F ) ).
Proposition 4.3. For all T ∈ T h , let ∆ T denote the set of mesh cells that share at least a face with T . Then, the following holds

h 2 T R m uh 2
La(T ) + h T J m uh 2 La(∂T ) T ∈∆T [h 2 T f m -f m h 2 La(T ) + u m -u m h 2 a,T + p m -p m h 2 c,T ], (67)
τ m [h 2 T R m ph 2 L d (T ) + h T J m ph 2 L d (∂T ) ] T ∈∆T [ Im h 2 T (g -π 0 g hτ )(s) 2 L d (T ) ds (68)
+ u m -u m h -u m-1 + u m-1 h 2 a,T + p m -p m h -p m-1 + p m-1 h 2 c,T + Im (p -π 0 p hτ )(s) 2 d,T ds]. Proof. Letting ψ T = ν a,T R m uh yields R m uh 2 La(T ) (R m uh , ψ T ) La(T ) = (f m h -f m , ψ T ) La(T ) + (A(u m h -u m ) + B * (p m h -p m ), ψ T ) La(T ) = (f m h -f m , ψ T ) La(T ) + a T (u m -u m h , ψ T ) -b T (ψ T , p m -p m h ), so that h T R m uh La(T ) h T f m h -f m La(T ) + u m -u m h a,T + p m -p m h c,T .
Similarly, using ψ F = ν a,F J m uh yields

(J m uh , ψ F ) La(F ) = T ∈∆F [a T (u m -u m h , π a,F ψ F ) -b T (π a,F ψ F , p m -p m h ) -(R m uh , π a,F ψ F ) La(T ) + (f m h -f m , π a,F ψ F ) La(T ) ],
whence (67) is readily deduced. The proof of (68) is similar.

Estimates (67)-(68) show that the first two terms bounding E spc in (65) yield local lower bounds for the approximation error. A similar result cannot be inferred for the last term in (65) because this term involves time-derivatives of the discrete u-component whereas the error in the u-component is only bounded in the L ∞ t (V a )-norm. Hence, this term behaves optimally with respect to the mesh size, but not necessarily with respect to the time step if the error in the u-component is not smooth in time. This observation is directly linked with the fact that the analysis relies on the natural stability norm for the continuous problem which provides a control on the u-component only in the L ∞ t (V a )-norm. Further theoretical work is needed to tackle this issue. In the present work, we will only verify in the numerical experiments presented in §5 that if the exact solution has a smooth behavior in time, this last term also converges optimally.

Estimates using elliptic reconstruction

For all n ∈ {0, . . . , N }, we define the elliptic reconstruction of (u n h , p n h ) ∈ V ah ×V dh as the functions (U n , P

n ) ∈ V a × V d such that a(U n , v) -b(v, P n ) = a(u n h , P ah v) -b(P ah v, p n h ), ∀v ∈ V a , (69)
d(P n , q) = d(p n h , P dh q) -b(δ t u n h , q -P dh q), ∀q ∈ V d , (70)
where P ah (resp., P dh ) denotes the L a -orthogonal projection from V a onto V ah (resp., the L d -orthogonal projection from V d onto V dh ). Henceforth, we use the convention that δ t u 0 h = 0 and δ t p 0 h = 0. Observe that

AU n + B * P n = A h u n h + B * h p n h , (71)
DP n = D h p n h + (B -B h )δ t u n h . (72)
Indeed, for all q ∈ V d , 71) is proved similarly. The key idea to estimate the errors u n -u n h a and p n -p n h c is to consider the decompositions

DP n , q d = -d(P n , q) = -d(p n h , P dh q) + b(δ t u n h , q -P dh q) = (D h p n h , P dh q) L d + ((B -B h )δ t u n h , q) L d = (D h p n h + (B -B h )δ t u n h , q) L d , (73) since D h p n h ∈ V dh . Equ. (
u n -u n h = ω n u -ρ n u , ω n u = U n -u n h , ρ n u = U n -u n , (74)
p n -p n h = ω n p -ρ n p , ω n p = P n -p n h , ρ n p = P n -p n . (75)
The quantities ω n u and ω n p can be bounded by a posteriori error estimates for the steady problem, while the quantities ρ n u and ρ n p can be estimated in terms of ω n u and ω n p and other computable quantities. The analysis requires a refinement of Hypothesis 2.2 by localizing the approximation property to mesh cells. In the context of poroelasticity, this assumption means that the steady-state version of the pressure equation yields elliptic regularity and that the finite element space V dh satisfies the usual approximation properties. Hypothesis 4.1. There exist a constant c 8 and a positive real number δ such that for all r ∈ L d , the unique solution φ ∈ V d of the dual problem d(q, φ) = c(r, q) for all q ∈ V d , is such that there is φ h ∈ V dh satisfying

T ∈T h h -2δ T [h -2 T φ -φ h 2 L d (T ) + h -1 T φ -φ h 2 L d (∂T ) ] ≤ c 8 r 2 c . (76)
We first estimate the quantities ω n u and ω n p . Lemma 4.1. In the above framework, the following holds for all n ∈ {0, . . . , N },

ω n u 2 a ≤ 2c 6 E n u + 2β 2 c 8 E n p,δ , (77)
ω n p 2 c ≤ c 8 E n p,δ . (78
ω n p 2 c ≤ c 8 T ∈T h h 2δ T [h 2 T Dω n p 2 L d (T ) + h T J d ω n p 2 L d (∂T ) ].
Using ( 14) and (72) yields for n ≥ 1,

Dω n p = DP n -Dp n h = D h p n h + (B -B h )δ t u n h -Dp n h = -R n ph ,
and this relation also holds for n = 0 by definition of R 0 ph . In addition, for all n ≥ 0,

J d P n = 0 since DP n ∈ V dh ⊂ L d . As a result, ω n p 2 c ≤ c 8 E n p,δ . (ii) Estimate of ω n u a . Observe that ω n u a = sup 0 =v∈Va a(ω n u , v) v a = sup 0 =v∈Va a(ω n u , v) -b(v, ω n p ) v a + β ω n p c = sup 0 =v∈Va a(ω n u , v -i ah (v)) -b(v -i ah (v), ω n p ) v a + β ω n p c .
owing to (69) since i ah (v) ∈ V ah . Using (50), (52), and (54) leads to

ω n u 2 a ≤ 2c 6 T ∈T h [h 2 T Aω n u + B * ω n p 2 La(T ) + h T J a ω n u -J b ω n p 2 La(∂T ) ] + 2β 2 ω n p 2 c .
Owing to ( 13) and (71),

AU n +B * P n = A h u n h +B * h p n h = -f n h so that Aω n u +B * ω n p = AU n +B * P n -Au n h -B * p n h = -R n uh . Moreover, J a ω n u -J b ω n p = -J n uh since AU n + B * P n ∈ L a .
This yields (77). We now turn our attention to the quantities ρ n u and ρ n p . Define the data, space, and time error estimators

E dat = U 0 -u 0 2 a + P 0 -p 0 2 c + 2E(f, g), (79)
E spc = N m=1 τ -1 m [2γ 2 c 8 (1 + 2β 4 ) E m p,δ (δ t ) + 4β 2 γ 2 c 6 E m u (δ t )], (80)
E tim = 2 3 γ2 τ 1 Cδ t p 1 h + Bδ t u 1 h -D h p 0h -g 1 h 2 L d + N m=2 2 3 γ2 τ 3 m δ t (Cδ t p m h + Bδ t u m h -g m h ) 2 L d . (81)
Lemma 4.2. The following holds for all n ∈ {1, . . . , N },

1 2 ρ n u 2 a + 1 2 ρ n p 2 c ≤ E dat + E spc + E tim . (82)
Proof. The estimates for ρ n u and ρ n p rely on the stability properties of the continuous problem. Thus, it is again convenient to handle equations holding a.e. in [0, T ] rather than at the discrete times {t n } N n=0 . Let U τ (resp., P τ ) be the continuous and piecewise affine function in time such that for all n ∈ {0, . . . , N }, U n τ = U n (resp., P n τ = P n ). Let ω uτ and ω pτ be constructed in a similar way from {ω n u } N n=0 and {ω n p } N n=0 . Define ρ uτ = U τ -u and ρ pτ = P τ -p. Observe that for a.e. t ∈ [0, T ] and for all v ∈ V a ,

a(ρ uτ , v) -b(v, ρ pτ ) = a(U τ , v) -b(v, P τ ) -(f, v) La = a(u hτ , P ah v) -b(P ah v, p hτ ) -(f, v) La = (f hτ , P ah v) La -(f, v) La = (f hτ -f, v) La , while for all q ∈ V d , c(∂ t ρ pτ , q) + b(∂ t ρ uτ , q) + d(ρ pτ , q) = c(∂ t P τ , q) + b(∂ t U τ , q) + d(P τ , q) -(g, q) L d = c(∂ t ω pτ , q) + b(∂ t ω uτ , q) + (Dπ 0 P τ , q) L d + d(P τ , q) + (π 0 g hτ -g, q) L d = c(∂ t ω pτ , q) + b(∂ t ω uτ , q) + (D(π 0 P τ -P τ ), q) L d + (π 0 g hτ -g, q) L d ,
where π 0 P τ is the piecewise constant function in time equal to P n on I n for all n ∈ {1, . . . , N }. Indeed, on each time interval

I n , c(∂ t p hτ , q) + b(∂ t u hτ , q) = (Cδ t p n h + Bδ t u n h , q) L d = (g n h + D h p n h -B h δ t u n h + Bδ t u n h , q) L d = (g n h + DP n , q) L d .
Testing the above equations with v := ∂ t ρ uτ and q := ρ pτ yields

1 2 d t ρ uτ 2 a + 1 2 d t ρ pτ 2 c + ρ pτ 2 d = (f hτ -f, ∂ t ρ uτ ) La + c(∂ t ω pτ , ρ pτ ) + b(∂ t ω uτ , ρ pτ ) + (D(π 0 P τ -P τ ), ρ pτ ) L d + (π 0 g hτ -g, ρ pτ ) L d , so that 1 2 d t ρ uτ 2 a + 1 2 d t ρ pτ 2 c ≤ (f hτ -f, ∂ t ρ uτ ) La + γ 2 ∂ t ω pτ 2 c + β 2 γ 2 ∂ t ω uτ 2 a + γ2 D(π 0 P τ -P τ ) 2 L d + π 0 g hτ -g 2 d .
Proceeding as usual yields for all n ∈ {1, . . . , N } the bound (details are skipped for brevity)

1 2 ρ n u 2 a + 1 2 ρ n p 2 c ≤ E dat + T 0 2γ 2 ∂ t ω pτ (s) 2 c ds + T 0 2β 2 γ 2 ∂ t ω uτ (s) 2 a ds + T 0 2γ 2 D(π 0 P τ -P τ )(s) 2 L d ds.
The second and third terms in the upper bound yield the space error estimator E spc owing to the fact that ω pτ and ω uτ are piecewise affine in time so that

T 0 ∂ t ω pτ (s) 2 c ds = N m=1 τ -1 m ω m p -ω m-1 p 2 c .
A similar bound holds for ω uτ . By linearity and proceeding as in the proof of Lemma 4.1 yields

ω m u -ω m-1 u 2 a ≤ 2c 6 E m u (δ t ) + 2β 2 c 8 E m p,δ (δ t ) and ω m p -ω m-1 p 2 c ≤ c 8 E m p,δ (δ t ).
Finally, the last term in the upper bound yields the time error estimator E tim since

T 0 D(π 0 P τ -P τ )(s) 2 L d ds = N m=1 1 3 τ m D(P m -P m-1 ) 2 L d ,
and for all m ≥ 0,

DP m = Dp m h -R m ph .
Theorem 4.2. For all n ∈ {1, . . . , N },

1 4 u n -u n h 2 a + 1 4 p n -p n h 2 c ≤ E dat + E spc + E tim + c 6 E n u + c 8 ( 1 2 + β 2 ) E n p,δ . (83)
Proof. Use Lemma 4.1, Lemma 4.2, and the triangle inequality.

We will not attempt here to prove lower error bounds for the above error estimators; this goes beyond the present scope. We will verify numerically in the following section that these estimators yield the expected, optimal, order of convergence with respect to mesh size.

h -1 0 u -u hτ a p -p hτ c ( T 0 p -p hτ 2 d ) 1/2 ( T 0 p -π 0 p hτ 2

Numerical results

We consider the following analytical solution of ( 1)-( 2 Tables 1 and2 present the convergence results (under space and time refinement respectively) for the approximation errors measured in various norms. All the convergence rates match those predicted by the a priori error analysis. An important observation is that in all cases, the total error is dominated by the L 2 t (H 1 x )-error on the p-component.

To assess the a posteriori error estimate obtained with the direct approach, see (49), we evaluate the quantities , (84) as well as the efficiency indices

η 1 = N m=1 τ m E m p,0 1
I eff = η 1 + η 2 + η 3 + η 4 ( T 0 [ p -p hτ 2 d + p -π 0 p hτ 2 d ]) 1/2 , I * eff = η 1 + η 2 + η 3 + η 4 u N -u N hτ a + p N -p N hτ c . ( 85
)
Recall that η 1 , η 2 , and η 3 are associated with the space error estimator and that η 4 is associated with the time error estimator. For brevity, we concentrate here on these two estimators. Tables 3 and4 present the results obtained under space and time refinement, respectively. All the observed convergence rates match the theoretical predictions. Moreover, the efficiency index I eff takes values between 2 and 3, indicating that the present estimators behave quite satisfactorily to control the pressure error in the L 2 t (H 1 x )-norm. As expected, the situation is quite different if one attempts to control the displacement error in the L ∞ t (H 1 x )-norm. As reflected by the efficiency index I * eff which increases as the mesh is refined, this latter error converges faster to zero than the estimator derived using the direct approach. 5. A posteriori error estimates using elliptic reconstruction and convergence rates under space refinement ; T = 1, τ = 0.1

To assess the a posteriori error estimate using elliptic reconstruction, see (83), we evaluate the quantities

η 5 = N m=1 τ -1 m ( E m p,1 (δ t ) + E m u (δ t )) 1 2
, η 7 = ( E N u )

1 2 , η 8 = ( E N p,1 ) 1/2 , (86)

η 6 = τ 1 Cδ t p 1 h + Bδ t u 1 h -D h p 0h -g 1 h 2 c + N m=2 τ 3 m δ t (Cδ t p m h + Bδ t u m h -g m h ) 2 c 1 2 , ( 87
)
as well as the efficiency index

J eff = η 5 + η 6 + η 7 + η 8 u N -u N hτ a + p N -p N hτ c . ( 88
)
Recall that η 5 is associated with the space error indicator and that η 6 is associated with the time error indicator. Moreover, η 7 and η 8 control the error between the discrete solution and its elliptic reconstruction; see Lemma 4.1. Table 5 presents the results obtained under space refinement. The observed orders of convergence match theoretical predictions, with a slight super-convergence for η 5 and a slight sub-convergence for η 7 on the (very) coarse meshes. The quantity η 6 , which is related to the time error, remains at a fairly constant value under space refinement. Additional tests (not reported here for brevity) indicate that η 6 converges with order close to 1 under time refinement. Finally, we observe that the efficiency index is quite large, especially on the (very) coarse meshes.

Conclusions

We have analyzed Euler-Galerkin approximations to coupled elliptic-parabolic problems with application to poroelasticity. In particular, we have obtained a posteriori error estimates delivering certified upper bounds for the approximation error in the sense that all the constants in the estimates have been specified. This can be an important feature in numerical simulations aiming at performance assessment of underground waste repositories. Furthermore, we have proposed an extension of the elliptic reconstruction technique introduced by Makridakis and Nochetto for linear parabolic problems to the present setting. This technique yields a posteriori error estimates that converge optimally for the u-component, but at the price of less favorable efficiency indices especially on (very) coarse meshes. Finally, we point out that the present work can be extended in many directions, including the use of time-dependent meshes and adaptive simulations driven by the present a posteriori error estimators.

)ω n p 2 c

 2 Proof. (i) Estimate of ω n p c . Let φ be the dual solution associated with the data r := ω n p in Hypothesis 4.1. Then, = c(r, r) = d(r, φ) = d(r, φ -φ h ), since owing to (70), d(r, φ h ) = d(P n -p n h , φ h ) = 0 for φ h ∈ V dh . Using (51) and (76) leads to

= 2π 2 κ b+ 1 M

 1 ) on the domain Ω = (0, 1) × (0, 1), u(t, x, y) = -exp(-At) 2π cos(πx) sin(πy) sin(πx) cos(πy) , p(t, x, y) = exp(-At) sin(πx) sin(πy), with A , κ = 0.05, b = 0.75, 1 M = 3 28 . The Lamé coefficients are λ 1 = 1 2 and λ 2 = 1 8 , yielding a Poisson ratio ν = 0.4 and a Young modulus E = 7 20 . Convergence rates in space are evaluated on a series of uniformly refined structured triangulations based on a boundary mesh step h 0 .

2 ,

 2

Table 1 .

 1 Errors at final time and convergence rates under space refinement ; T = 0.1,

	d ) 1/2

Table 2 .

 2 Errors at final time and convergence rates under time refinement ; T =1, h 0 = 1/128

Table 3 .

 3 .93 2.57e-2 1.89 2.75e-3 1.57 7.67e-4 0.92 2.18 15.33 16 1.71e-2 0.96 6.63e-3 1.96 7.13e-4 1.94 3.89e-4 0.98 1.70 23.81 32 8.62e-3 0.98 1.68e-3 1.98 1.80e-4 1.99 1.96e-4 0.99 1.45 41.19 A posteriori error estimates using the direct approach and convergence rates under space refinement ; T = 0.1, τ = 2.50e-4 .37e-3 1.02e-4 1.04e-4 3.85e-2 0.90 1.51 4.95 0.1 4.45e-3 9.93e-5 8.25e-5 2.01e-2 0.93 1.62 5.58 0.05 4.49e-3 1.02e-4 7.38e-5 1.03e-2 0.96 1.77 6.69

	h -1 0	η 1		η 2		η 3			η 4		I eff	I * eff
	4	6.34e-2	-	9.53e-2	-	8.17e-3	-	1.45e-3	-	3.13 11.43
	8	3.33e-2 0τ	η 1	η 2		η 3	η 4		I eff I * eff
		0.25 4.33e-3 1.07e-4 1.15e-4 4.70e-2	-	1.48 4.80
		0.2 4

Table 4 .

 4 A posteriori error estimates using the direct approach and convergence rates under time refinement ; T = 1, h 0 = 1/128 .74 7.82e-3 5.70e-3 1.44 5.83e-3 1.90 120.26 16 4.36e-2 2.21 7.58e-3 1.79e-3 1.67 1.50e-3 1.96 97.49 32 1.04e-2 2.07 7.55e-3 4.93e-4 1.86 3.77e-4 1.98 19.59 Table

	h -1 0	η 5		η 6	η 7		η 8		J eff
	4	1.34e-0	-	1.14e-2 1.55e-2	-	2.18e-2	-	169.43
	8	2.01e-1 2

Acknowledgment. The authors are grateful to Olivier Boiteau, Clément Chavant and Gérald Nicolas (EdF R&D) for fruitful discussions and instrumental help with the implementation in Code Aster.

* This work was partially supported by EdF R&D and GDR MOMAS (PACEN/CNRS, ANDRA, BRGM, CEA, EdF, IRSN). The second author was supported by a CIFRE PhD fellowship.