Adaptive poromechanics computations based on a posteriori error estimates for fully mixed formulations of Biot’s consolidation model

This paper is concerned with the analysis of coupled mixed finite element methods applied to the Biot's consolidation model. 
We consider two mixed formulations that use the stress tensor and Darcy velocity as primary variables as well as the displacement and pressure. The first formulation is with a symmetric stress tensor while the other enforces the symmetry of the stress weakly through the introduction of a Lagrange multiplier. The well-posedness of the two formulations is shown through Galerkin's method and suitable a priori estimates. The two formulations are then discretized with the backward Euler scheme in time and with two mixed finite elements in space. We present next a general and unified a posteriori error analysis which is applicable for any flux- and stress-conforming discretization. Our estimates are based on $H^{1}(\Omega)$-conforming reconstruction of the pressure and a suitable $\left[H^{1}(\Omega)\right]^{d}$-conforming reconstruction of the displacement; both are continuous and piecewise affine in time. These reconstructions are used to infer a guaranteed and fully computable upper bound on the energy-type error measuring the differences between the exact and the approximate pressure and displacement. The error components resulting from the spatial and the temporal discretization are distinguished. They are then used to design an adaptive space--time algorithm. Numerical experiments illustrate the efficiency of our estimates and the performance of the adaptive algorithm.

stabilized discretization for the classical P 1 -RT 0 -P 0 approach was recently proposed in [51]. One can also see [13,25,51,46,47] and the references therein for more discretization schemes. More recently, the mixed finite element methods as a numerical approximation of the four-field formulation have received increasing attention. In this formulation, two pairs of mixed finite elements spaces, one the for linear elasticity and one for mixed Poisson problem, are used. This formulation uses the stress, displacement, flux, and pressure as primary unknowns. The extra unknowns can be seen as a disadvantage against the two-or three-field formulations, regarding the computational cost, but there are reasons to prefer this approach. First, these methods are widely used in computational fluid mechanics because they keep conservation law of mass and balance of forces and produce continuous normal fluxes, regardless of mesh quality. Secondly, the stress and the Darcy velocity as variables are of primary interest in many applications. This eliminates the need to use postprocessing techniques to produce the flux and stress from the numerical pressure and displacement. The availability of these variables is necessary whenever the consolidation is coupled with transport e.g., applications involving multiphase flow [16] or thermal convection [20,21,22,24]. This formulation allows also the use of domain decomposition techniques, where in particular, transmission conditions across physical interfaces, involving the stress and flux are imposed (cf. [2,36,54]).
In this work we consider coupled mixed methods for Biot's consolidation model, and we shall mainly be concerned with the a posteriori error analysis for this model problem. Today, rigorous a posteriori error estimates are well developed for a large class of simple linear parabolic model problems [31,56,42]. These developments often serve as starting point for extensions to diverse complex applications, including nonlinear effects, time-dependent loads or multi-physics phenomena [23,28,29,48]. The challenges for applying the framework of a posteriori error estimates to the Biot problem are numerous. First, the Biot system is a time-dependent and strongly coupled problem. Second, this system does not have an "easy"monotone structure or maximum principle. To the best of our knowledge, contrarily to the case of a priori error estimates [38,59], no results are available for a posteriori error estimates for the Biot's consolidation model in the context of discretizations with conforming flux and stress fields, such as mixed finite elements methods. Residual-based error estimates for the two-field formulation can be found, e.g., in [30]. Recently, a posteriori error analysis of this later formulation has been derived in [48], where the discrete solution is achieved using Taylor-Hood H 1 -conforming finite elements in space (using piecewise polynomials of order k ≥ 1 for the pressure and of order (k + 1) for the displacement) and a backward Euler scheme in time. The estimates are based on flux and stress reconstructions to compute the error indicators that bound the dual norm of the residual. Therein, the spatial and the temporal discretization errors are distinguished, then used in an algorithm adapting the mesh and the time so as to equilibrate these error sources.
The purpose of this article is twofold. Firstly, to fill the gap of lacking rigorous well-posedness results for the mixed formulation, we analyze in the multi-dimensional case the existence and uniqueness of a weak solution to two mixed formulations presented in [38,59] and show their equivalence in the continuous case. To this aim, we use Galerkin's method together with the theory of DAEs (Differential Algebraic Equations) and suitable a priori energy estimates. Secondly, we propose a general and unified a posteriori error estimation for the Biot's consolidation model discretized with fully mixed finite elements method (MFE) in-space and a backward Euler scheme in-time. The error estimation is applicable, in a larger sense, for any locally conservative method, such as cell-centered finite volume scheme, multipoint mixed finite element, mimetic finite difference and hybrid high-order discontinuous Galerkin [17,35,45,50]. The present theory readily extends to conforming methods using equilibrated flux and stress reconstructions ( cf. [48]). For the present mixed setting, the error estimation requires H 1 (Ω)-conforming reconstruction of the pressure and a H 1 (Ω) d -conforming reconstruction of the displacement, both are continuous and piecewise affine in-time.
These reconstructions are used to infer a guaranteed upper-bound on the error between the exact and the approximate solution. Precisely, we will show that an energy-type-norm difference between the exact and the approximate pressure and displacement can be bounded by the dual norm of the residuals [23,27,28,31]. The estimation is then established using this error measure by deriving a guaranteed and fully computable upper bound on the dual norm of the residuals, without unknown constant (cf. [1,3,34,48]). The bound uses easily, fully and locally calculable estimators. We will also show how a posteriori error estimators can distinguish between the space and time errors. We show in particular how to distinguish the pressure and displacement contribution in the overall errors. This is essential for the development of space-time adaptive marching algorithm and particularly for the development of adaptive iterative coupling schemes based on time-splitting techniques, e.g., see the fixed-stress scheme in [13,19].
The paper is organized as follows. In Section 2, we present the physical model and the assumptions on the effective parameters. Section 3 recalls the four-and five-field formulations of the Biot problem as well as the relevant functions spaces. We prove the well-posedness of the aforementioned formulations and show in particular their equivalence. In Section 4, we present the discrete formulations, by combining in-space two mixed finite elements, and the backward Euler scheme in-time. In Section 5, we give the main result of the paper, which states that the energy-type-norm of the difference between the exact weak and the approximate solutions can be bounded by a fully computable a posteriori error estimate. Section 6 is devoted to the proof of the a posteriori error estimate. In Section 7, the a posteriori error estimate is elaborated by distinguishing its space and time error components. This is efficiently used to propose a space-time adaptive marching algorithm. Finally, illustrative numerical results are shown in Section 8.

Model problem
Let Ω be a connected polygonal domain of R d , d = 2, 3. We denote by ∂Ω its boundary (supposed to be Lipschitz-continuous) and by n the unit normal to ∂Ω, outward to Ω. Let a time interval (0, T ) be given with T > 0. We consider Biot's consolidation equations, modeling flow in deformable porous media. The mathematical form of this problem as it is presented in [15,52] reads: find u and p such that there hold, −∇·θ(u) + α∇p = f , in Ω × (0, T ), (2.1a) ∂ t (c 0 p + α∇·u) − ∇·(K∇p) = g, in Ω × (0, T ), (2.1b) u(·, 0) = u 0 , p(·, 0) = p 0 , in Ω, (2.1c) u = 0, p = 0, on ∂Ω × (0, T ), (2.1d) where f is the body force, µ and λ are the Lamé parameters, α is the Biot-Willis constant, K is the permeability tensor divided by fluid viscosity, c 0 is the constrained-specific storage coefficient, g is the volumetric source term, u 0 is the initial displacement and p 0 is the initial pressure. The function θ denotes the effective stress tensor, i.e., θ(u) := 2µ (u) + λtr( (u))I, where (u) is the linearized strain tensor given by (u) := (∇u+∇ T u)/2 and where tr denotes the trace of matrices. The total stress tensor σ(p, u) is given by σ(p, u) := θ(u) − αpI, presenting the internal forces on surface elements. The momentum balance for the fluid is interpreted as the Darcy law for the volumetric fluid flux w := −K∇p. For the sake of simplicity, we consider only homogeneous boundary conditions, but the results of this paper can be extended to more general case (cf. [48]). In what follows, we assume the following assumptions on the effective parameters: 1. Let c 0 , α, µ, and λ be strictly positive constants.
2. Let K, be a symmetric and uniformly positive definite tensor which satisfies the following assumption: there exist positive constants c K and C K such that for or a.e. x ∈ Ω and for all ζ ∈ R d 0 < c K |ζ| 2 ≤ ζ T K −1 (x)ζ, and |K(x)ζ| ≤ C K |ζ|, (2.2) and whose terms are for simplicity supposed piecewise constant on the meshes defined below and constant in time.
3. We assume that the initial pressure p 0 and the initial displacement u 0 lie in H 1 0 (Ω) and respectively.

respectively.
Results on well-posedness of the problem (2.1) were established by Showalter [52]. In the next section, we present and analyze two mixed formulations of (2.1). The first formulation extends to the multidimensional case the four-field formulation introduced in [59] that uses symmetric space for the stress unknown, while the second formulation generalizes the five-field formulation introduced in [38] that uses the weak stress symmetry formulation of the elasticity system. This latter formulation is advantageous compared to the four-field formulation, since only mixed finite elements with weakly symmetric stress will be used. This is advantageous because it requires less computational costs and as the hybridization techniques of Fraeijs de Veubeke (cf. [9]) lead usually a system with reduced sizes.

Mixed variational formulations
We now introduce the function spaces and their norms that will be used throughout. Firstly, we will use the convention that if V is a space of functions, then we designate by V a space of vector functions having each component in V , and we designate by V the space of tensor functions having each component in V . The space L 2 (Ω) is the space of square-integrable functions endowed with its natural inner product written (·, ·) L 2 (Ω) with associated norm denoted by || · ||. We designate by H 1 (Ω) the usual Sobolev space and by H 1 0 (Ω) for its zero-trace subspace. The corresponding norm and semi-norm are written || · || H 1 0 (Ω) and | · | H 1 0 (Ω) , respectively. In particular, H −1 (Ω) is the dual of H 1 0 (Ω). Further, let H(div, Ω) be the space of vector-valued functions from L 2 (Ω) that admits a weak divergence in L 2 (Ω). Its natural norm is ||v|| div,Ω := ||v|| 2 + ||∇·v|| 2 1 2 .

First variational formulation
Following [59] and the references therein, we derive the first mixed formulation for the problem (2.1). Let where A is the fourth-order compliance tensor given by and known to be bounded and symmetric definite uniformly with respect to x ∈ Ω, and where we used in (3.1b) the following relationship derived by taking the trace-operator on both sides of (3.1a). The mixed variational formulation of the fourfield formulation (3.1), (2.1d)-(2.1c) referred from now to as Problem A reads as follows: for a.e. t ∈ (0, T ), together with the initial condition (2.1c). Note that in (3.4b), we used that (tr(ξ), q) = (ξ, qI) for any matrix ξ and a scalar q. In what follows the constrained specific storage coefficient c 0 is assumed to be strictly positive and satisfies: The condition (3.5) is necessary to show the well-posedness of Problem A. For parameter ranges of practical problems, it is typical that 0 < α ≤ 1, so for d = 3 and α > 1 . We omit any further discussion on the justification for these constraints, other than they are necessary to prove the energy estimates. We also refer to [12,39,21] for a more detailed discussion of scaling of the Biot system.
The proof of this result uses Galerkin's method together with the theory of DAEs (Differential Algebraic Equations) and a priori energy estimates (cf. [18,40,21]). It is detailed in Appendix A. In turn, the well-posedness of problem (3.1), (2.1d)-(2.1c) stems from the Lemma given next.

Second variational formulation
To introduce the second mixed formulation of the problem (3.1), (2.1d)-(2.1c), we introduce a new unknown ζ (=rotation) as the skew-symmetric part of ∇u and extend the definition of A from the symmetric tensor space to all tensors. We denote this extension by A as well, since it will also be given by formula (3.2) (cf. [38] for the 2D case). With these definitions, equation (3.1a) can now be replaced by and then enforcing the symmetry of σ weakly in the sense that its anti-symmetric part is zero tested against a skew-symmetric tensor, i.e., (σ, γ) = 0, ∀γ ∈ Q sk , where Q sk = [L 2 (Ω)] d×d sk denotes the subspace of [L 2 (Ω)] d×d composed of skew symmetric-valued tensors. The mixed formulation for the resulting problem referred from now as Problem B reads as follows: for a.e. t ∈ (0, T ), find (σ(t), u(t), w(t), (c 0 + c r )(∂ t p, q) + c r dα (∂ t σ, qI) + (∇·w, q) = (g, q), ∀q ∈ Q, (3.9b) Furthermore, we can give the following equivalence result. Proof. Let (σ, u, w, p, ζ) be a solution of Problem B, then σ is symmetric, i.e., σ(t) ∈ H s (div, Ω), and therefore (σ(t), On the other hand, if (σ, u, w, p) solves Problem A, then its stems from Lemma 3.2 that u(t) ∈ H 1 0 (Ω) and, if we set ζ(t) to the skew-symmetric part of ∇u(t), then (σ, u, w, p, ζ) as given in (3.10) solves Problem B.

A fully discrete scheme based on MFE in space and the backward Euler scheme in time
In this section, after introducing some notations, we present the fully discrete formulations of Problem A and Problem B. We employ for both formulations a backward Euler scheme in-time, and in-space, two mixed finite elements for the linear elasticity and flow problems.

The temporal and spatial meshes
For integer values N ≥ 0, let (τ n ) 0≤n≤N denote a sequence of positive real numbers corresponding to the discrete time steps such that T = N n=1 τ n . Let t 0 = 0, and t n = n j=1 τ j , 1 ≤ n ≤ N be the discrete times. Let I n = (t n−1 , t n ], 1 ≤ n ≤ N . For every time step 0 ≤ n ≤ N , we let v n h := v hτ (·, t n ) for any sufficiently smooth function v hτ . For the spatial meshes, we consider a discretization T n h of the domain Ω at time t n , consisting of simplicial elements K. We assume that the meshes T n h , 0 ≤ n ≤ N , are conforming, i.e., for K, L ∈ T n h with K = L, then K ∩ L is either an empty set or a common vertex or edge or face. We suppose that the initial mesh T 0 h is introduced to approximate the initial condition. The meshes are then refined or coarsened as time progresses. Typically, the mesh T n h is obtained from T n−1 h by refining some elements and coarsening some other ones. For all 1 ≤ n ≤ N , we denote by T n−1,n h , a common refinement of T n−1 h and T n h .

The discrete functional spaces
Let 0 ≤ n ≤ N be fixed. We denote by P l (K) the space of polynomials on K ∈ T n h of degree less than or equal to l. We use the notation · K for the norm in L 2 (K), K ∈ T n h . The corresponding inner product is (·, ·) K . Let |K| be the Lebesgue measure of K ∈ T n h . We define the broken Sobolev space H 1 (T n h ) as the space of all functions v ∈ L 2 (Ω) such that v| K ∈ H 1 (K), for all K ∈ T n h . The (broken) energy-norm on H 1 (T n h ) is given by where the sign ∇ denotes the element-wise gradient, i.e., the gradient of a function restricted to a mesh element K ∈ T n h . The (broken) energy-norm in L 2 (Ω) is given by Let E be a space of functions defined on Ω. We denote P 1 τ (E) the vector space of functions continuous in time and with values in E. We also denote by P 0 τ (E) the space of functions piecewise constant in time and with values in E. We have then if v hτ ∈ P 1 τ (E), then ∂ t v hτ ∈ P 0 τ (E) is such that for all 1 ≤ n ≤ N , In what follows, we denote respectively by c K,K and C K,K the smallest and the largest eigenvalue of the tensor K in K ∈ T n h . For all 0 ≤ n ≤ N , let Q n h × W n h ⊂ L 2 (Ω) × H(div, Ω) be the Raviart-Thomas-Nédélec mixed finite element spaces of order zero on the mesh T n h (cf. [9]): where RTN 0 (K) denotes the lowest-order Raviart-Thomas-Nédélec finite-dimensional subspace associated with the element K ∈ T n h ; any v h ∈ RTN 0 (K) takes the form [P 0 (K)] d + P 0 (K)x for the example of simplices. The degrees of freedom of v h ∈ RTN 0 (K) correspond to the values of the flux of v h across the faces e ⊂ ∂K.
For all 0 ≤ n ≤ N , let Q n h × W n s,h ⊂ L 2 (Ω) × H s (div, Ω) be the Arnold-Winther mixed finite elements for the lowest-order stresses on the mesh T n h (cf. [10]): where AW 1 (K) denotes the Arnold-Winther stress space of order 1, i.e., any tensor be the Arnold-Falk-Winther mixed finite elements for the lowest-order stresses on the mesh T n h (cf. [11]):  For simplicity of presentation, we assume that f ∈ L 2 (0, T ; L 2 (Ω)) and g ∈ L 2 (0, T ; L 2 (Ω)), respectively. We further assume that both f and g are piecewise-constant-in-time with respect to the temporal meshes introduced in Section 4.1. Otherwise, we refer the reader to Remark 5.5 for a more general setting for the source terms.

The MFE scheme
The discrete form of Problem A reads as follows (cf. [59]): given σ 0 That of Problem B reads (cf. [38]): given σ 0 Let us remark that at the discrete level the two schemes are not equivalent since the discrete stress tensor in the five-field formulation (4.4) will not necessarily inherit the symmetry property from the four-field formulation (4.3). However, the weakening the symmetry constraint in (4.4) permits; i) the use of simpler element spaces for the stress; ii) the use of a local static condensation permitting to reduce the MFE system to a symmetric and positive definite one for the triplet pressure, displacement and rotation with the same way as in [5,6]. This system is smaller and easier to solve than the original saddle point problem, but no further reduction is possible. Although the analysis of this reduction has not been conducted, the static condensation for the coupled problem is expected to follow along the same lines as in [45]. In what follows, we denote by (σ hτ , u hτ , w hτ , p hτ , ζ hτ ) the discrete space-time functions piecewise constant in time given

Postprocessing
We present here an improved approximation of the couple (p hτ , u hτ ) by local postprocessing. This step is also necessary for the energy a posteriori error estimate for the above mixed finite element schemes. We first introduce a postprocessing of the scalar variable p hτ as described in [31]. We define the improved We also let p 0 h = Π V 0 h p 0 denote an approximation to the initial datum, where a typical choice being the L 2 −orthogonal projection onto V 0 h . Such a postprocessing is local and its cost is negligible and valid for the two mixed finite schemes described above (cf. [9,57] for more details). In contrast to the postprocessing of the pressure, we should distinguish that of the displacement inherited from the two schemes. The postprocessing of the displacement from the first scheme (4.3) is defined as follows (cf. [8,34,41] where Π W n s,h and Π Q n h are respectively the L 2 -orthogonal projections onto W n s,h and Q n h . Typically, the spaces where e i ∈ R d denotes the i-th Euclidean unit vector. For both cases, we set u 0 h = Π M 0 h u 0 as approximation of the initial displacement datum. Finally, we define the continuous, piecewise affine in-time functions p hτ and u hτ by p hτ (·, t n ) = p n h , u hτ (·, t n ) = u n h , 0 ≤ n ≤ N. (4.8) The advantage of the above postprocessing is threefold: i) the original solution is improved and ∇ p hτ = 0 as well as ∇ u hτ = cst; ii) we have these estimates: ||K∇ p n h + w n h || K = 0 from (4.5) and ||Aσ n h − ∇ u n h + α 2µ + dλ p n h I + ζ n h || K = 0 from (4.7) or that ||Aσ n h − ∇ u n h + α 2µ + dλ p n h I|| K is negligible from (4.4); iii) the important relationship (3.3) between the stress tensor, displacement and pressure is restored (see Lemma 5.2 for more details).

The a posteriori error estimates
We present in this section our a posteriori error analysis of the two mixed schemes given in the previous section. The idea is to estimate the error between the (exact) weak solution (p, u) of (2.1) and its approximate solution ( p hτ , u hτ ) in an energy norm.

Assumptions on the pressure and displacement reconstructions
The main ingredient to derive our error estimates is to construct the following functions.
Definition 5.1 (Reconstructions). We will call pressure and displacement reconstructions any couple of functions (p hτ ,û hτ ) reconstructed from ( p hτ , u hτ ) such that and satisfying for all 0 ≤ n ≤ N , This definition implies that the pair (p hτ ,û hτ ) is defined by the (N + 1) approximations (p n h ,û n h ) ∈ H 1 0 (Ω) × H 1 0 (Ω), associated with the discrete times {t n } 0≤n≤N . We will prescribe a concrete choice of reconstructions in Section 5.4. Let us define by ϕ(p, u) := c 0 p + α∇·u the fluid content. Then, we will employ the abridged notationφ n h := ϕ(p n h ,û n h ) and ϕ n h := ϕ( p n h , u n h ), for all 0 ≤ n ≤ N . An important result of the above reconstructions is given in the following Lemma: Lemma 5.2 (Properties of (p hτ ,û hτ )). Letp hτ andû hτ be the reconstructed functions as given in Definition 5.1. Then, for all 1 ≤ n ≤ N , there hold Proof. For all 1 ≤ n ≤ N , and for all K ∈ T n h , we have from the definition of ϕ and the notation (4.2) together with the divergence theorem: . Thus, for the first term, we infer that (∂ tp n h , 1) K = (∂ t p n h , 1) K for all K ∈ T n h . Similarly, using (5.3) together with the mean condition (5.1b), it is inferred that Thus, (5.2a) holds true. The local conservation property of the mixed schemes (4.3) and (4.4) implies at each time step n, 1 ≤ n ≤ N , Taking the trace-operator on both sides of (4.7a) or the first equation of (4.6), using the fact that ∇· u n h = tr( ( u n h )) = tr(∇ u n h ), then replacing the result in (5.2a) gives at each time step n, 1 ≤ n ≤ N , Replacing the above result in (5.4), we obtain (5.2b).
The condition (5.2a) guarantees that the mean values of the time derivative of ( p hτ , u hτ ) are preserved by (p hτ ,û hτ ). The condition (5.2b) implies that the postprocessed pair ( p hτ , u hτ ) together with condition (5.2a) preserves the local conservation property of the mixed schemes. These conditions are of crucial importance for the a posteriori error estimate.

The error estimators
Let a time step 1 ≤ n ≤ N , and a mesh element K ∈ T n h be fixed. We define the local residual estimators and the local flux estimators Herein h K denotes the diameter of K, and C P,K is the constant from Poincaré inequality: where π 0 q is the mean value of q on the element K and C P,K = 1/π whenever the element K is convex. Furthermore, to capture the nonconformity from the numerical schemes, we define the local nonconformity estimators Clearly all of the above estimators are local-in-space and in-time. Finally, we define the initial data estimators
The second point for our error estimation is that of the error measure. For this purpose, we introduce the following energy-type error measure: The above norms are well-defined owing to the properties of the weak solution (p, u) and the reconstructed functions (p hτ ,û hτ ), i.e., (p −p hτ , u −û hτ ) ∈ E T . For the quantity (p − p hτ , u − u hτ ) which is not in E T , we extend the definition (5.12), where the gradient and divergence are understood in the broken sense.
We are now ready to state the main result of this Section. Theorem 5.3 (A posteriori error estimate). Let (p, u) be the weak solution of the problem (2.1) given by (5.11). Let ( p hτ , u hτ ) be the approximate pressure and displacement as obtained in Section 4.4. Let (p hτ ,û hτ ) be the reconstructed pressure and displacement as given in Definition 5.1. Then there holds where and where we have set L P = 1 and L U = 1 µ , and for 1 ≤ n, l ≤ N , Clearly, the above error estimate is without unknown constant. Its proof is given in Section 6. We now present several remarks.
Remark 5.4 (Conforming methods). The postprocessing of the pressure and displacement as well as their reconstructions are not needed for conforming methods (e.g., finite elements and vertex-centered finite volumes), but equilibrated reconstructions of the Darcy velocity and of the total stress tensor are required to derive the above a posteriori error estimate (cf. [48]). In that case, we have (p n h , u n h ) ∈ H 1 0 (Ω) × H 1 0 (Ω), and we just set (p hτ ,û hτ ) = (p hτ , u hτ ) and consequently η NC = 0.
Remark 5.5 (Time oscillation of the source terms). For the sake of simplicity, we assumed the source terms f and g to be piecewise constant in time. When this is not true, the following data time-oscillation estimators should be added to the estimate (5.14), η osc,g := g − g X T , and η osc, where g ∈ P 0 τ (L 2 (Ω)) is such that g| In = g n where we have set g n = 1 τ n In g(·, t) dt, and similarly for f . In that case, the residual estimators (5.6) are now given by Remark 5.6 (The initial data estimator). The proof of Theorem 5.3 uses Grönwall's Lemma which implies the appearance of the term e T , and this is only within the approximation of the initial condition. In practice, this does not affect the overall estimate as the initial conditions are known and thus the initial data estimators (they are data oscillation errors) η IC,P and η IC,U can be made small enough or neglected (see the numerical results or [1,23,28] for a similar results). Furthermore, (p hτ ,û hτ )(·, 0) = (p 0 , u 0 ), if the pair p 0 and u 0 are piecewise polynomials, hence η IC,P = η IC,U = 0.
An attractive feature of the estimate (5.14) is the ability to distinguish the contribution of the pressure and displacement in the total error as in [48]. To this purpose, each of the nonconformity estimators η n NC2,K and η N NCF,K is split into two parts using the triangle inequality, i.e.,

The reconstructed pressurep hτ and displacementû hτ
We make precise here the reconstructions introduced in Definition 5.1. We first reconstruct from the approximate pressure p n h a functionp n h , H 1 0 (Ω)-conforming and from the postprocessed displacement u n h a functionû n h , H 1 0 (Ω)-conforming, both satisfying the corresponding mean value constraint from (5.1b). To do so, we proceed as in [31]; at each Lagrangian node x situated in the interior of Ω, we set where for all K ∈ T n,n+1 h , b K denotes the standard (time-independent) bubble function supported on K, defined as the product of the barycentric coordinates of K, and scaled so that its maximal value is 1, and where I av : P 2 (T n h ) → P 2 (T n h ) ∩ H 1 0 (Ω) is the interpolation operator given by with T x being the set of all elements of T n h sharing the node x. At the Lagrange nodes x situated on the boundary ∂Ω, we setp n h (x) := 0. In order to guarantee that the mean value constraint (5.1b) (first equation) holds true, we choose a n K = The same procedure applied to the approximate displacement u n h ∈ [P m (T n h )] d , m ≥ 2, leads to a H 1 0 (Ω)conforming vector function as we set when x is a Lagrange node situated in the interior of Ω. If x ∈ ∂Ω, we setû n h (x) = 0. Here, I av is given by (Ω), and where b K denotes the vector bubble function supported on K. Similarly, in order to obtain a H 1 0 (Ω)-conforming vector function while maintaining the mean condition (5.1b) (second equation), we choose As usually, we define the continuous, piecewise affine in-time functionsp hτ andû hτ bŷ 6 Proof of the a posteriori error estimate The proof of Theorem 5.3 will be broken down in three steps (i) we first show that using the reconstructions introduced in Definition 5.1, the energy-type-norm (5.12) stating the difference between the exact and the reconstructed solutions can be bounded by the dual norm of the residuals (ii) we then bound the dual norm of the residuals by fully computable estimators (iii) we end the proof by putting the previous pieces together and incorporating the nonconformity error in the overall error using the triangle inequality.

Definition of the residuals
We introduce here the residuals stemming from the weak formulation in Section 5.3. To this aim, we let || · || H 1 0 (Ω) stand for the energy norm and we equip X T with the energy norm That of X T is left equipped with Recalling the weak formulation (5.11), and that (p hτ ,û hτ ) ∈ E T . We define the residuals R P (p hτ ,û hτ ) ∈ X T and R U (p hτ ,û hτ ) ∈ X T such that The dual norms of the residuals are defined in a standard way and given by 6.2 Bounding the error of the displacement and the pressure by the dual norm of the residuals We now proceed to step (i) of the proof, i.e., we show that the dual norm of the residuals (6.3) bound the energy-type norm (5.12). The results are obtained through duality technique [23,28]. To this aim, we first give this intermediate result.
Now, we are able to state the main result of this Section.
Theorem 6.2 (Upper bound on the error by the residuals). Let (p, u) be the weak solution of the problem (2.1) given by (5.11). Let (p hτ ,û hτ ) ∈ E T . Then, there holds Proof. First, using Grönwall's Lemma, we infer from (6.4) Qs e t−s ds , We integrate the above inequality over the interval (0, T ), to obtain Qs e t−s ds dt. (6.12) in both sides of (6.4) for t = T , then we use (6.12) to bound the last term, we end up with (6.10). Proof. For the existence of a weak solution to (5.11), we can adopt the arguments in [52] for a slightly different weak formulation or adapt the results of Theorem 3.1. For the uniqueness of a weak solution, we suppose that (p hτ ,û hτ ) satisfies (5.11). This implies first that R P (p hτ ,û hτ ) = R U (p hτ ,û hτ ) = 0. Thanks to the estimate (6.10), we obtain that (p −p hτ , u −û hτ ) 2 en ≤ 2e T −1 2 ||ϕ(p −p hτ , u −û hτ )(·, 0)|| 2 H −1 (Ω) , and hence the uniqueness of the weak solution for a given initial data, as well as the continuous dependence with respect to the initial data.

Bounding the dual norm of the residuals by the estimators
For step (ii) of the proof, we shall bound the dual norm of the residuals by fully computable a posteriori error estimators. As usual, to work with the nonconforming approximation ( u hτ , p hτ ), the X T − and X T −norms are again extended to piecewise regular-in-space functions, i,e, Proof. Let 1 ≤ n ≤ N be given and let q ∈ X T with |||q||| X T = 1. Adding and subtracting (w n h , ∇q) to (6.2a) and applying the Green theorem, we obtain R P (p hτ ,û hτ ), q X T ,X T = We apply Lemma 5.2, the Poincaré inequality (5.7) and Cauchy-Schwarz inequality, it is inferred that Applying the Cauchy-Schwarz inequality yields (6.14a). Proceeding in a similar way for (6.14b) with v ∈ X T such that |||v||| X T = 1; first recalling that where we used the fact that σ(p hτ ,û hτ ) is symmetric so we replaced (v) by ∇v in the first term on the right-hand side of (6.16). Now, adding and subtracting (σ n h , ∇v) to (6.16), applying the Poincaré inequality (5.7) together with the Cauchy-Schwarz inequality, and using the local conservation property of the mixed scheme for the mechanics problem, we obtain We obtain (6.14b) by applying again the Cauchy-Schwarz inequality.

End of the proof of Theorem 5.3
In step (iii) of the proof, we combine the upper bound in (6.10) with the results of the previous section to obtain a posteriori estimate for the displacement and pressure errors. First, let us recall that our approximate solution ( p hτ , u hτ ) is not an element of E T . For this reason, we cannot apply (6.10) to (p− p hτ , u− u hτ ) en . Thus, we decompose (p − p hτ , u − u hτ ) en into two parts using the triangle inequality, i.e., (p − p hτ , u − u hτ ) en ≤ (p −p hτ , u −û hτ ) en + (p hτ − p hτ ,û hτ − u hτ ) en , (6.18) and then apply (6.10) to bound (p −p hτ , u −û hτ ) en as follows We then proceed as in [28,Thm. 5.3] or [1, Thm. 6.1] using the results of Section 6.3 to bound the terms involving ||R P (p hτ ,û hτ )|| X t and ||R U (p hτ ,û hτ )|| X t , t ∈ (0, T ), with the estimators. There remains to give a computable upper bound to (p hτ − p hτ ,û hτ − u hτ ) en . We have from the definition (5.12) Qs e t−s ds dt Σs e t−s ds dt. (6.20) For the first two terms of (6.20), we have Now, since the quantity ( ϕ hτ −φ hτ ) is piecewise affine and continuous in time, and since for all 1 ≤ n ≤ N and on each element K ∈ T n h , ( ϕ n h −φ n h ) has zero mean value by (5.1b): For the fourth term, we have first Let v hτ = ϕ hτ −φ hτ , v hτ is piecewise affine and continuous in time, it is then inferred that Now, for all 1 ≤ n ≤ N and on each element K ∈ T n h , v n h has zero mean value by (5.1b), so one can proceed as for (6.21a), to infer Gathering the above results in (6.18) to obtain the estimate (5.14).

Distinguishing the error components and adaptive algorithm
In this section, we first elaborate the error estimate (5.14), so as to distinguish the error components resulting from the spatial and the temporal discretization. Once these error contributions are identified, they will be used to design an adaptive algorithm where the mesh size and time step are locally adapted as in [26,28,31,48].  where δ a = 0 for a = tm and δ sp = 1. We have this result: Theorem 7.1 (A posteriori error estimate distinguishing error components). Let (p, u) be the weak solution of the problem (2.1) given by (5.11). Let ( p hτ , u hτ ) be the approximate pressure and displacement as obtained in Section 4.4. Let (p hτ ,û hτ ) be the reconstructed pressure and displacement as given in Definition 5.1. Then, there holds Proof. The idea is to split the diffusive flux estimator η n F,P,K (t) as well as the stress estimator η n F,U,K (t) into two contributions using the triangle inequality (see [4,31,48]); for all K ∈ T n h , t ∈ I n , We then integrate in time the last terms on the right-hand side of of above inequalities and use the equality of norms in (7.4a), i.e., |||K∇v||| ,K = ||K − 1 2 (K∇v)|| K = |||v||| K . Inserting the resulting inequalities in (5.14) with the use of definition (7.1)-(7.2) leads to the estimate (7.3).

Balancing and space-time adaptivity
In this section we use the error estimate in Theorem 7.1 in an algorithm with adaptive refinement/coarsening of the space-time mesh, based on balancing the space and time errors. To do so, first the space and time error components should be equilibrated by selecting the time step τ n and adjusting the spatial meshes T n h in such a way that γ tm η n sp ≤ η n tm ≤ Γ tm η n sp , (7.5a) where η n J = η n J,P + η n J,U , for J = sp, tm, and where γ tm , Γ tm are two user-given weights, typically close to 1. The goal of this balancing is to have η n tm and η n sp of comparable size instead of getting η n tm much smaller than η n sp . This balancing can be defined locally (element-wise) as γ tm,K η n sp,K ≤ η n tm,K ≤ Γ tm,K η n sp,K , ∀K ∈ T n h , where γ tm,K and Γ tm,K are the local (element-wise) weights. In the next step, we have to adapt the space mesh; refinement/coarsening so that the local spatial error estimators are distributed equally: for all To account for limited computing resources, we fix refinement thresholds h min > 0 and τ min > 0 for the mesh size and the time step, respectively, i.e., we require that for all 1 ≤ n ≤ N , We also need to fix threshold that the error on the time interval I n should not exceed some fixed error e n , for all 1 ≤ n ≤ N . Based on the above considerations and balancing conditions, we propose the following algorithm: Algorithm 1: Complete solution algorithm with adaptive space-time refinement Data: Enter f , g, p 0 , σ 0 and Ω, T , K, c 0 , µ, λ and α.
Result: The quadruplet (σ hτ , u hτ , w hτ , p hτ ). 1 Choose an initial mesh T 0 h , an initial time step τ 0 , and set t 0 = 0; 2 Choose the weights γ tm and Γ tm , and the refinement thresholds τ min and h min ; 3 Calculate the initial data estimators η IC,P and η IC,U /* Initial mesh adaptation loop */; 4 Adapt the mesh T 0 h such that the local initial condition error estimators η IC,P,K and η IC,U,K are distributed equally; 5 while t n ≤ T do /* Time loop */ until the balance(7.5a) or (7.5b) is satisfied or τ n ≤ τ min ; 20 Refine or coarsen the cells K ∈ T n h until (7.5c) is satisfied or h K ≤ h min ; 21 until η n sp + η n tm ≤ e n ;

22
(σ n h , u n h , w n h , p n h ) ← (σ n K n , u n K n , w n K n , p n K n ) K n ∈T n h ; 23 t n ← t n−1 + τ n ; 24 end

Numerical results
This section discusses the computational results of applying our a posteriori error estimates when solving (2.1) using the two mixed schemes (4.3) and (4.4).

Example 1: analytical solution
To illustrate numerically our theoretical results, we consider in Ω = [0, 1] 2 the analytical solution of the Biot system (2.1) (see [19]): The effective parameters are K = I, c 0 = α = 1, and µ = λ = 0.6. We choose source terms, initial and Dirichlet boundary conditions such that (8.1) is the solution of problem (2.1). First, convergence in space and time are evaluated on a series of space or time uniform refinement. These results are then compared with those obtained with the adaptive algorithm (Algorithm 1). We compare the estimated errors for the pressure and displacement with the analytical errors measured in the following energy-norms Σs e t−s ds dt. (8.2b) derived from (5.12). In practice, the dual norms in the above norms or in the energy norm (5.12) are incalculable and then they are approximated by solving auxiliary problems as detailed in [28].

Standard case
For both schemes, we set τ n = 1/128 in order to compare the convergence rates of the spatial discretization errors ignoring influences of time discretization errors. In Table 1, we present the convergence results of Method 1 with AF 1 -elements. The results of Method 2 with AFW 1 -elements are presented in Table 2. The last column in both tables displays the corresponding effectivity indices, given by the ratio of the estimates over the global error measured in the energy norm given by (5.12)  We observe that all the convergence rates are in agreement with the expected convergence rates; they are, in particular, not influenced by the time discretizations as well as the overall error is dominated by the errors from the mechanics problem. This seems to indicate that the schemes with postprocessed solutions have similar convergence properties. We also regard the effectivity indices corresponding to the two methods   as similar and likewise excellent; Method 1, in particular, gives a slightly better efficiency. For the temporal refinement, the two methods exhibit similar convergence properties with optimal effectivity indices (not presented). Figures 1 and 2 depict the distribution of the estimated error and analytical errors at the final time step. Clearly, the two methods reflect the distribution of the analytical error in the whole domain.

Adaptive case
In the second part of this test case, we compare the previously error estimates obtained on uniform and fixed meshes and time steps with the estimates obtained using Algorithm 1. We start with an initial mesh of 311 elements. The initial time step τ 0 was chosen as 0.001 and was not changed by the criterion η IC ≤ e 0 . The other parameters used in the algorithm are γ tm = 0.75 and Γ sp = 1.3 and e n = 5 · 10 −3 τ n . Once more the two methods produce nearly identical results and therefore only those results pertaining Method 2 are presented. In Figure 3 (left and center), we give the initial and final meshes. In Figure 4, we compare the estimated error at the final time to the exact error for Method 2. Again, the distribution of the estimated error reflects the exact one and it is evenly distributed throughout the domain. In Figure 3 (right), we plot the evolution of the error estimates depending on the total number of spacetime unknowns. As shown for both methods, one can reduce the number of unknowns necessary to attain the prescribed precision using the derived estimators and adaptively refined space-time grids. Finally, one has to mention that each problem solves for Method 2 is inexpensive compared to Method 1 as only fewer degree of freedoms are required with weakly imposed stresses. This is particularly important for threedimensional problems, where Method 2 is feasible and easier in practice than Method 1 that has highly number of degrees of freedom (stress with 162 degrees of freedom per element).

Balancing criteria
In order to illustrate the importance of the space-time balancing criterion (7.5a) in the adaptive algorithm, we show by a series of computations its impact on the overall estimate. To this aim, we apply Algorithm 1 with Method 2 on two computations in which the space and time errors are not equilibrated. In the first case, we set γ tm = 2 and Γ sp = 3 and we let e n = 5 · 10 −3 τ n as previously. This case corresponds to an over-refinement in-space. In the second case, we choose γ tm = 1/3 and Γ sp = 1/2 which will induce an over-refinement in-time. In Figure 5 (left), the ratio of the time error over the space error as a function of Figure 5: The ratio of the time error over the space error as a function of the total number of spacetime unknowns (left); over-refinement in-space (green) or in-time (blue) are violating criterion (7.5a) in comparison with a computation honoring (7.5a) (black). A comparison of the induced overall errors (right). the total number of space-time unknowns is depicted for the two situations as well as for the case where the space and time errors are equilibrated. The effect of these choices on the overall estimate is shown in Figure 5 (right). Clearly, one can conclude that violating the balancing criterion (7.5a) reduces the efficiency of the algorithm in terms of precision as well as the total computational cost is considerably increased, thus, confirms that (7.5a) is an essential ingredient for our adaptive algorithm.

Example 2: flow through saturated levee
Coupled fluid-structure interaction problems appear in many geotechnical applications, e.g., earthen dams, levees, and other geotechnical structures. This example considers the behavior of a saturated levee after a flood event. One of the major reasons of levee failure is erosion due to overtopping due to this flood event.
Here, we consider the flow through saturated levee in advance of any overtopping [20]. The geometry of the levee is given in Figure 6, and consists of a lower and an upper part (lower 5 [m] and upper 10 [m], respectively). For the flow problem, we impose a Dirichlet boundary condition p = 1 · 10 2 × x(15 − y)  We test Method 2 through the adaptive algorithm with various thresholds for the errors, i.e., we set initially e n = 8 · 10 −3 τ n refined then 3 times by a factor of 2. The other parameters used in the algorithm are γ tm = 0.7 and Γ sp = 1.25. We plot in Figure 7 the meshes at the final time with e n = 4 · 10 −3 τ n and e n = 1 · 10 −3 τ n . We can see that reducing the threshold e n allows more mesh refinement but only around the seepage face. In Figure 8, we compare the estimated error obtained with space-time adaptivity (with e n = 4 · 10 −3 τ n ) to the estimated error with standard method with uniform space-time meshes. We clearly observe the performance of the adaptive algorithm to evenly distribute the total error over the domain as well as reducing it around the seepage face.
In the last part of this test case we analyze the performance of the adaptive algorithm in terms of precision against the number of unknowns. To this aim, we plot in Figure 9 (left) the estimated error and the actual error calculated with a reference solution approximated on very fine space-time discretization as a function of the total number of space-time unknowns in the adaptive and standard cases. As expected, these results confirm that the use of the adaptive algorithm reduces significantly the number of spacetime-unknowns for approximately the same value of the error. The right part of Figure 9 displays the corresponding effectivity indices, given by the ratio of the estimates over the actual error calculated with a reference solution approximated on a very fine space-time discretization. The effectivity index confirms a moderate overestimation of the total error and proves the accuracy of the estimate (5.14).

Example 3: flow inside an aneurysmal sac
In this study we aim to investigate the adaptive algorithm for the fluid-dynamics of blood in the presence of an aneurysm. This is one of the most dangerous anomalies that may affect the arteries, and its treatment is very difficult. In most cases, the proposed therapy consists in introducing a metallic multi-layered stent (see Figure 10) inside the porous sac. This device slows down the vorticity in the aneurysm by lowering the flow velocity, pressure, displacement and stresses in the aneurysmal sac. Numerical simulation is a powerful approach for the investigation of the relationship between stenting and aneurysm re-growth and rupture, and can provide scientific guidelines for the design of stent structure which may reduce the velocity and the oscillatory shear stress inside the sac [33,43].
For this test case, the domain Ω is as given in the right of the Figure 10 in which four circular struts are aligned in the aperture of the aneurysmal sac. The initial and boundary conditions are extracted from solving a coupled Biot-Stokes equations in the arterial geometry given in the left of Figure 10. That is, the coupled Biot-Stokes equations were first solved on the extended domain (including the artery area in Figure 10), where the Biot system is used inside the porous aneurysmal sac (lower part) and a Stokes flow is used in the remaining part of the artery (upper part), and where Γ 1 is the interface between the two subdomains where natural transmission conditions are used. Note that when the permeability tends to zero in the poroelasticity problem, a Stokes-type problem is obtained (see [51]). The resulting solution inside the aneurysmal sac is used to generate a realistic time-dependent boundary condition on Γ 1 for the present test case [7]. For this setup, the initial pressure and displacement are set accordingly to u(0) = 0 and p(0) = 10 −1 . On the inlet boundary Γ in a parabolic velocity profile was enforced with peak value of 1000 mm/s, multiplied by sin(π(1 − t/40) 3 ) to obtain a pulse like flow pattern. On the outlet boundary Γ out , we set a free-stress boundary condition and on Γ ext a no-flow boundary condition. On the outer boundary of the porous sac Ω as well as on the boundary of the struts Γ o , we impose no-flow and -displacement in normal direction, i.e., w · n| Γ1∪Γo = 0 and u · n| Γ1∪Γo = 0. The total simulation time is T = 30 s and the time step is τ n = 3 · 10 −4 s. The remaining model parameters are the blood viscosity which is set equal to 3.5 · 10 −3 N/(s·m 2 ), K = 10 −6 , α = 0.1, c 0 = 5.0 · 10 −3 , µ = 4.28 · 10 3 and λ = 1.07 · 10 3 , satisfying (3.5).
With the generated data, we compare the performance of the adaptive algorithm with Method 2 to two computations, both are with fixed meshes and time steps. The parameters used in the algorithm are e n = 3 · 10 −2 τ n , γ tm = 0.8 and Γ sp = 1.25. The first computation is set as a reference solution, where the discretization is chosen very fine and the second one (standard) is chosen in a way to have approximately the same number of space-time unknowns as in the adaptive algorithm. The results related to the three computations are depicted in Table 3  We easily remark that the adaptive solution is much better than the standard one, for approximately the same computational cost. We also notice that the error from the flow problem dominates the total error, which can be considerably reduced if a multi-scale space-time discretization is used. In Figure 11, we plot the first component of the normal stress tensor σ h · n and the first component of the velocity w h along the red line tangent to the inserted stent (as indicated in the right Figure 10), at the final time T . We can see that the results with the adaptive algorithm are very satisfactory compared to the standard one. The adaptive algorithm is particularly effective in detecting the oscillations on the solution induced by the presence of the stent. The approximate solution in the adaptive and standard cases is depicted in Figure 12. In Figure 13, we compare the adaptive estimated error to the standard one at the final time T . The adaptive error is not only evenly distributed throughout the domain but also reduced around the corners and the struts of stent.

Conclusions
Two mixed finite element formulations of the Biot's consolidation model are presented and analyzed. Each method utilizes a pair of mixed finite element methods, one for the linear elasticity and one for the mixed Poisson problem. The methods differ by using either symmetric stress space or weakly symmetric stress through the introduction of Lagrange multiplier enforcing the symmetry of the stress in a weak sense. We analyzed the existence and uniqueness of a week solution to the formulations and showed in particular their equivalence at the continuous level. We then propose a unified and guaranteed a posteriori error estimates for the Biot system when discretized with coupled mixed finite elements method in-space and intime with a backward Euler scheme. The a posteriori error estimate requires conforming reconstructions of   Figure 13. the displacement and pressure; both are continuous and piecewise affine in time. Precisely, regardless of the mixed method, we present a guaranteed and fully computable error estimate that bounds the energy-type error between the unknown exact solution and the approximate solution. This estimate is then decomposed into estimators characterizing the space, and time error components. We also separate the pressure error components from those of displacement errors. This is efficiently used in an adaptive space-time algorithm where the individual error component estimators are used to adjust the time step, and to select the mesh elements to be refined or coarsened. Numerical results illustrate the effectiveness of our estimates and the performance of the adaptive algorithm. Finally, we shall point out that the developments in this paper also apply to any locally conservative scheme such as cell-centered finite volume scheme, mimetic finite difference and hybrid high-order discontinuous Galerkin. The present theory can be extended to derive adaptive iterative coupling schemes based on a posteriori error estimates for Biot's consolidation model , which are the subject of ongoing work.
A Appendix: Proof of Theorem 3.1 The proof will be broken down in several steps: we first construct solutions of certain finite-dimensional approximations of (3.1), then we derive suitable energy estimates in Lemma A.6 and prove then the Theorem.
Therefore, let (σ n , u m , w k , p l ) : [0, T ] 4 → W s,n × Q m × W k × Q l be the solution of the following problem: The proof of this Lemma will be broken down in several steps. First, we write a DAEs system equivalent to (A.1). To do so, we introduce the following notations: We also use the following notations; P l to denote the vector of degrees of freedom of p l (t) with respect to the basis {q i } l 0 , S l to denote the vector of degrees of freedom of σ n (t) with respect to the basis {τ i } n 0 , U m to denote the vector of degrees of freedom of u m (t) with respect to the basis {v i } m 0 and W k to denote the vector of degrees of freedom of w k (t) with respect to the basis {z i } k 0 . Now, using the above notations, and for all (n, m, k, l) ∈ N 4 , problem (A.1) may be rewritten as The above equations yield a linear system of DAEs type for the vector X(t) = (P l (t), S n (t), W k (t), U m (t)): where E(t) = (G p (t), 0, 0, −F u (t)) T and So it suffices to see that the matrix pencil sB + D is nonsingular for some s = 0 to prove the existence of the discrete solution satisfying the initial conditions (A.6). We set s = 1 and show that B + D is nonsingular. The results presented next extends those in [59] to any dimensional. Note that we can write B + D as a block 2 × 2 matrix as follows The following lemma will be useful in proving the invertibility of B + D (see Yi 2014 [59] for 2D case). Let U = W s,n × Q l and R = Q m × W k , and define the bilinear forms φ A : U × U → R, φ B : U × R → R and φ C : R × R → R, where φ A ((σ n , p l ), (τ , q)) = (c 0 + c r )(p l , q) + c r dα (Ip l , τ ) + c r dα (σ n , qI) + (Aσ n , τ ), (A.12) and with kernel spaces given by There exist a constant C = C(λ, µ, c 0 , d) > 0 such that Proof. Denoting by τ = {τ i,j } 1≤i,j≤d , with τ i,j = τ j,i , for i = j, and using the definition of the fourth order compliance tensor (3.2), together with the Young and triangle inequalities yields Letting = dα , we obtain φ A ((τ , q), (τ , q)) ≥ c 0 ||q|| 2 + 1 2µ where all the coefficients are strictly positive for any dimension d = 2, 3.
The proof of the following Lemma is straightforward from the properties of K.
Proceeding as in [59, Lemma 3.4& 3.5], we promptly obtain the following result.

A.2 A priori estimates
The next step is to derive a priori estimates with the aim to pass to the limit in a subsequence of the approximated solution.