Guaranteed and computable error bounds for approximations constructed by an iterative decoupling of the Biot problem

The paper is concerned with guaranteed a posteriori error estimates for a class of evolutionary problems related to poroelastic media governed by the quasi-static linear Biot equations. The system is decoupled employing the fixed-stress split scheme, which leads to a semi-discrete system solved iteratively. The error bounds are derived by combining a posteriori estimates for contractive mappings with those of the functional type for elliptic partial differential equations. The estimates are applicable for any approximation in the admissible functional space and are independent of the discretization method. They are fully computable, do not contain mesh dependent constants, and provide reliable global estimates of the error measured in the energy norm. Moreover, they suggest efficient error indicators for the distribution of local errors, which can be used in adaptive procedures.


Introduction
Problems defined in a poroelastic medium contribute to a wide range of application areas including simulation of oil reservoirs, prediction of the environmental changes, soil subsidence and liquefaction in the earthquake engineering, well stability, sand production, waste deposition, hydraulic fracturing, CO 2 sequestration, and understanding of the biological tissues in biomechanics. In recent years, mathematical modeling of poroelastic problems has become a highly important topic because it helps engineers to understand and predict complicated phenomena arising in such media as well as assists in preventing possible future financial calamities. However, numerical schemes designed for any of the existing models provide approximations that contain errors of different nature, which must be controlled. Therefore, a reliable quantitative analysis of poroelasticity problems requires efficient and computable error estimates that could be applied for various approximations and computation methods.
Mathematical modeling of poroelasticity is usually based on the Biot model that consists of the elasticity equation coupled with the equation describing the slow motion of the fluid. Computational errors in one part of the model may seriously affect the accuracy of the other one. Therefore, getting reliable and efficient a posteriori error estimates for coupled problems is, in general, a much more complicated task than for a single equation.
The Biot model is a system describing the flow and the displacement in a porous medium by momentum and mass conservation equations. Initially, it was derived at a macroscopic scale (with inertia effects negligible) in the works Terzaghi [1] and Biot [2]. The settlement of different types of soils was predicted in [1], which later was extended to the generalized theory of consolidation [2,3]. A comprehensive discussion of the theory of poromechanics can be found in [4]. Thus, for modeling of the solid displacement u and the fluid pressure p, we consider the system that governs the coupling of an elastic isotropic porous medium saturated with slightly compressible viscous single-phase fluid −div λ (divu) I + 2 µ ε(u) − α p I = f in Q := Ω × (0, T ), ∂ t β p + αdivu − divK∇p = g in Q, (1.1) where Q denotes a space-time cylinder (with bounded domain Ω ⊂ R d , d = {2, 3} having Lipschitz continuous boundary ∂Ω and given time-interval (0, T ), 0 < T < +∞), f ∈ H 1 (0, T ; [ L2(Ω)] d ) and g ∈ L 2 (0, T ; L 2 (Ω)) are the body force and the volumetric fluid source, respectively 1 . The first equation in (1.1) follows from the balance of linear momentum for the total Cauchy stress tensor σ por := σ(u) − α p I that accounts not only u but also the pressure p scaled by the dimensionless Biot-Willis coefficient α > 0. Linear elastic tensor is governed by Hooke's law σ(u) := 2 µ ε(u) + λ trε(u) I = 2 µ ε(u) + λ (divu) I, where ε(u) := 1 2 ∇ u + (∇ u) T is the tensor of small strains. Here, λ, µ > 0 are the Lamé constants proportional to Young's modulus E and Poisson's ratio ν via relations µ = E 2 (1+ν) and λ = E ν (1+ν) (1−2ν) . The second equation is the fluid mass conservation (continuity) equation in Q. Here, β stands for the storage coefficient and K is the permeability tensor assumed to be symmetric, uniformly bounded, anisotropic, and heterogeneous in space and constant in time, i.e., λ K |τ | 2 ≤ K(x)τ · τ ≤ µ K |τ | 2 , λ K , µ K > 0, for all τ ∈ R d . (1.4) For the fluid content β p + α divu, we prescribe the following initial condition η(x, 0) := β p(x, 0) + α divu(x, 0) = β p • + α divu • , where p • and u • are defined in (1.3). To simplify the exposition, we consider only homogeneous BCs, i.e., p D , z N = 0 and u D , t N = 0 for the time being, even though all results are valid for more general assumptions. The work [5] provides the results on existence, uniqueness, and regularity theory for (1.1)- (1.4) in the Hilbert space setting, whereas [6] extends the recent results to a wider class of diffusion problems in the poroelastic media with more general material deformation models. Corresponding a priori error estimates can be found in [7]. Considered system can be understood as the singular limit of the fully dynamic Biot-Allard problem (see the details in [8]), where the acceleration of the solid in the mechanics part of (1.1) is neglected. As the Biot model is a coupled system of partial differential equations (PDEs), we have iterative as well as monolithic approaches used for solving the problem (see, e.g., [9]). For the first approach, the problem can be reformulated with a contractive operator, which naturally yields iteration methods for its solution (see [8]). On each step in the time, the flow problem is considered first. It is followed by solving the mechanics using already recovered pressure. The procedure is iterated until the desired convergence is reached. Different alteration of iterative cycles in flow and mechanics, i.e., single- [10] and multi-rate schemes [11,12], can be considered. The second approach is fully coupled and considers the system with two unknowns simultaneously.
The question of a posteriori error control for the poroelastic models has been already addressed using different techniques. Application of the residual-based error estimates to the coupled elliptic-parabolic problems can be traced back to works [40,41]. Recently, similar error indicators were used in [42,43,44,45,46] for immiscible incompressible two-or multi-phase flows in porous media to address the questions of adaptive stopping criteria and mesh refinement. In [47], authors suggest an a posteriori error estimator based on the appropriate dual problem in space-time for a coupled consolidation problem involving large deformations. In [48,49,50,51], adaptive space-time algorithms, relying on equilibrated fluxes technique, were applied to the Biot's consolidation model (formulated as a system with four unknowns). In this work, however, we turn to the functional error estimates (majorants and minorants), which are fully computable and provide guaranteed bounds of errors arising in the numerical approximations. The derivation of such estimates is based on functional arguments and variational formulation of the problem in question. Therefore, the method does not use specific properties of approximations (e.g., Galerkin orthogonality) and special properties of the exact solution (e.g., high regularity). The estimates do not contain mesh dependent constants and are valid for any approximation in the natural energy class. Moreover, the majorant also yields an efficient error indicator, that provides mesh adaptation.
Our main goal is to deduce efficient a posteriori error estimates for the approximation of the system (1.1) and verify them in application to the numerical problems. In [52], a posteriori error estimates of the functional type has been derived for the stationary Barenblatt-Biot model of porous media. This paper deals with more complicated Biot problem presented by an elliptic-parabolic system of partial differential equations. Our approach is based on the contraction property of the iterative method [53], which is rather general and not just restricted to fixed-stress scheme and functional type estimates of each equation in the Biot system (see, e.g., [54]). To the best knowledge of the authors, it is the first study targeting such a coupling between the elastic behavior of the medium and the fluid flow in a context of functional error estimates. The main result is presented in Theorem 3. Moreover, they not only serve as reliable estimates of the global error measured in the energy norm but also as an efficient indicator of the local error distribution over the computational domain (confirmed by majorant application to each of the decoupled problems in [55] and references therein). The latter property makes functional majorants advantageous in application to automated adaptive mesh generation algorithms.
The paper has the following structure: Section 2 is dedicated to the generalized formulation of the Biot system and its semi-discrete counterpart derived after applying the explicit Euler scheme in time. In Sections 3, we introduce an incremental approach, namely, fixed-stress split scheme, for discretizing considered coupled system. In particular, it justifies the optimal choice of parameters in the iterative scheme and proves that it is a contraction with an explicitly computable convergence rate. Sections 4 and 5 are dedicated to the derivation of auxiliary Lemmas used in the proof of Theorems 2 and 3 with the general estimates for the approximations generated by the fully decoupled iterative approach. Finally, Section 6 contains a collection of examples, which illustrates the application of derived error estimates to the Biot problem.

Variational formulation and discretization
We study approximations of the system (1.1), where V ≡ H 1 (0, T ; [H 1 (Ω)] d ) denotes the space for u (field of displacements) and W ≡ H 1 (0, T ; H 1 (Ω)) is the space for the variable p (pressure). The generalized setting of (1.1) reads: find a pair (u, p) where The Biot system of type (2.1)-(2.2) was analyzed by several authors to establish existence, uniqueness, and regularity. The first theoretical results on the existence and uniqueness of a (weak) solution are presented in [56] for the case of β = 0. Further work in this direction can be found in [5,57]. The well-posedness of the quasi-static Biot system is ensured under the above-mentioned assumptions. In fact, [58,59] established contractive results in suitable norms for iterative coupling of (2.1)-(2.2). For an overview of the stability of existing iterative algorithms, we refer the reader to [60,61].
The system (2.1)-(2.2) can be viewed as the two-field formulation of the poroelasticity problem. In numerical analysis, there are alternative approaches to treat such a system as three-and four-field formulations. In the case of the three-field model, an additional variable is introduced to represent the flux in the flow equations, whereas the four-field model considers stress as yet another unknown. Three-field formulation is rather flexible (since it allows different combinations of discretizations) and is accepted by the community as providing more physical approximations of the unknowns than in the two-field case. Recently, four-field formulation received increasing attention from the research community, where both equations were treated by mixed methods. The advantages of the latter representation are local conservation of mass and momentum balance and more accurate representation of fluxes and stresses. The choice of the formulation (from the above-mentioned list) is usually motivated by the considered application as well as the restriction on the computational resources. For instance, the mixed formulation of (2.2) does not only provide the flux that satisfies the local mass conservation property but also generates an effective approximation of this function, which is advantageous to functional type error control. It practically minimizes the majorant related to the pressure function (see (4.2)). The same method/principle works for the reconstruction of the stress field.
The system (2.1)-(2.2) is considered in the time-interval [0, T ] divided by N sub-intervals, such that it forms the corresponding set T N = ∪ N n=1 I n , I n = (t n−1 , t n ). Let u n (x) ∈ V 0 and p n (x) ∈ W 0 , where respectively, the spatial parts of the solution at t = t n . Then, the semi-discrete approximation of (2.1)-(2.2) satisfies the system where τ n = t n −t n−1 . This system generates the following problem to be solved on each step of the time-incremental method: find the pair (u, p) n ∈ V 0 × W 0 (K τ n ∇p n , ∇w) Ω + β (p n , w) Ω + α (div u n , w) Ω = ( g n , w) Ω , ∀ w 0 ∈ W 0 , (2.5) where K τ n := τ n K, the right-hand side of (2.5) is defined as and the pair (u, p) n−1 ∈ V 0 ×W 0 is given by the previous time step. The initial values are chosen as (u, Since from now on, we deal only with the semi-discrete counterpart of Biot problem, we omit the subscript Ω in the scalar product. Moreover, we always consider (2.4)-(2.5) on nth time step, which allows us to omit the superscript n for the rest of the paper and consider the system This work aims to derive a fully guaranteed a posteriori estimates of the error between the obtained approximations (ũ,p) ∈ V 0h × W 0h , where V 0h and W 0h are discretization spaces of conforming approximations of functional spaces V 0 and W 0 , respectively, and the pair of the exact solutions (u, p) of the Biot system, which is accumulated from the errors on all N time steps, i.e., e u := u −ũ and e p := p −p.
On each time step, these errors are measured in terms of the combined norm (e u , e p ) := ||| e u ||| 2 u + ||| e p ||| 2 p . (2.9) In turn, each term of the error norm is defined as follows: ||| e u ||| 2 u := ε(e u ) 2 2µ + div(e u ) 2 λ and ||| e p ||| 2 where w 2 λ := Ω λ w 2 dx, ε(w) 2 2µ := Ω 2µ ε(w) : ε(w) dx, and w 2 Kτ := Ω K τ w · w dx are L 2 -norms respectively weighted with 2µ, λ, and tensor K τ for any scalar-and vector-valued functions w and w. The global bound of the errors e u and e p contains incremental contributions from each time-interval, i.e., For the iterative approach, on each time-step I n , the Biot system is decoupled into two sub-problems related to the linear elasticity and a single-phase flow problem. Then, an iterative procedure is applied to obtain the pair (u i , p i ) = (u, p) i . Next, each equation is discretized and solved, such that instead of (u, p) i the pair (u, p) i h , containing the approximation error of the numerical method, is used. In Section 4, we derive computable a posteriori estimates for this pair of approximate solution.
Functional M := M (n) (Theorem 3) is derived by combining the estimates derived for the contractive mapping [53] and the a posteriori error majorants for the elliptic problems (initially introduced [63,64]). The validity of such estimates is based on the contraction property of specifically constructed linear combination of displacement and pressure α γ divu i − L γ p i , L, γ > 0, (the so-called volumetric mean stress). The selection of the parameters L and γ is justified and explained in Section 3.

Remark 1.
For the alternative monolithic approach, one solves (2.7)-(2.8) simultaneously for pressure and displacement, reconstructing the pair of approximations (ũ,p) = (u h , p h ) = (u, p) h . For this case, the corresponding computable bound of the error between (ũ,p) and the exact solution can be derived. Such a functional error bound is a combination of a posteriori error estimates for each of the unknowns in (2.7) and (2.8) (see, e.g., [65,54] and references therein). For the detailed derivation, we refer the reader to [62].
Remark 2. We note that due to the Korn and Friedrichs' inequalities, both e u 2 and ||| e u ||| 2 u are estimated by ε(e u ) 2 . Moreover, the physical bound on the Lamé parameters is given as d λ + 2 µ > 0, in the most general case, allowing for the first parameter λ to be slightly negative for so-called auxetic materials. In this case, we use the fact that ||| e u ||| 2 u := ε(e u ) 2 2µ + div(e u ) 2 λ (6.5) ≤ (2 µ + d λ) ε(e u ) 2 holds, and work with the positively-weighted norm ε(e u ) 2 . However, as auxetic materials are rare, the added complexity associated with allowing for such cases has been avoided in this paper. Consequently, the proofs below are based on the assumption of non-negative Lamé parameters.

Fixed-stress splitting scheme
Formal application of the iteration method to (2.7)-(2.8) yields the problem to be solved on the ith iteration step: where the flow equation (3.1) is solved for p i , using u i−1 , and the elasticity equation (3.2) is used to reconstruct u i using p i recovered on the previous step. To obtain the initial data for the iteration procedure, we first set the pressure equal to the hydrostatic pressure, i.e., it follows from ∇p 0 = ρ f g with ρ f > 0 denoting the fluid phase density. Whereas u 0 is reconstructed by (3.2) using p 0 . The iteration proposed in (3.1)-(3.2) is known as the fixed strain, and is only conditionally stable.
To stabilize the iteration scheme (3.1)-(3.2), we consider the 'fixed-stress splitting approach', which properties were initially studied in [61] and [58], respectively. This scheme operates with a special quantity: volumetric mean total stress where γ and L are certain positive tuning parameters. These parameters are usually kept constant on each half-time step. The optimal choice of γ and L allows us to prove that this iteration scheme is a contraction in the L 2 -norm δη i 2 , where δη i := η i −η i−1 . Moreover, it accelerates the iteration procedure by reducing the number of iterations. By adding L (p i − p i−1 ) into the right-hand side of (3.1), we rewrite the system (3.1)-(3.2) using the definition (3.3) as Theorem 1 establishes a contraction-type inequality for the norm δη i 2 .
Remark 3. The estimates in Theorem 1, satisfying contraction estimate (3.7), also holds for Galerkin approx- Moreover, in Appendix, we show that similar contraction theorem holds for sequence the {δ(η − η h ) i } ∈ W h , where η i ∈ W is generated by the fixed-stress split iterative scheme defined in (3.4)-(3.5) and η i h ∈ W h is discretization of the latter sequence. Generally, it is important to note that all the theorems and lemmas below are formulated for a pair (u, p) i ∈ V 0 × W 0 that forms a contraction w.r.t. to (u, p) i−1 ∈ V 0 × W 0 and its discrete approximation (u, p) i h ∈ V 0h × W 0h that forms a contraction relative to (u, Remark 4. There exist alternative ways to choose the tuning parameter L. In particular, the physically motivated choice L cl = α 2 λ+2µ/d is considered in [61]. Whereas, [59] suggests L opt = α 2 2 (λ+2µ/d) . The recent study [66] suggests the numerical evidence on the iteration counts w.r.t. the full range of the Lamé parameters for heterogeneous media. Numerical investigation of the optimality of these parameters and comparison with physically and mathematically motivated values from the literature was done in [67]. The authors demonstrated that their optimal value is not only dependent on the mechanical material parameters but on the boundary conditions and material parameters associated with the fluid flow problem.
Remark 5. The inequality (3.7) shows that the sequence {δη i } i∈N is generated by a contractive operator. Therefore, due to the Banach theorem, it tends to a certain fixed point. Moreover, since all the terms in the left-hand side of (3.7) are positive, in practice, {δη i } i∈N might converge with even better contraction rate than q = L β+L .
Corollary 1. From Theorem 1, it follows that ∇δp i and ε(δu i ) in (3.6) are also converging sequences and satisfy respectively.
Corollary 1 is used in the derivation of the error estimate for the term ||| e p ||| 2 p . In particular, it yields the following result based on the estimates for the Banach contraction mappings (see [53,54]). Lemma 1 (Estimates for contractive mapping). Let the assumptions of Theorem 1 hold. Then, we have the estimates By taking a limit m → ∞ and noting that in this case (q m + q m−1 + . . . + 1) → 1 1−q , we arrive at (3.8). The inequality (3.9) is proved by similar arguments.
If in Lemma 1, we consider iteration i and i − m, i > m as two subsequent iterations, a more general version of estimates (3.8) and (3.9) can be formulated.
Lemma 2 (General estimates for contractive mapping). Let the assumptions of Theorem 1 hold. Then, we have the estimates Proof: Proof follows along the lines of the proof of Lemma 1 with m = 1, . . . , i.

Remark 6.
Estimates in Lemma 2 improves the value of q m (1−q m ) 2 and q 2m (1−q m ) 2 if q is close to 1. It might look contra intuitive, but the choice m = i is not always most optimal. Simultaneously with the decrease of quotient with q, the term η i − η 0 2 grows. Therefore, the choice of m must be made carefully. Computing majorant on each step of our iterative algorithm might be computationally expensive, therefore, in demanding applications, the choice of m must be done a priori. Alternatively, with several extra iterations made after reaching the desired convergence in p and u, both q (1−q ) 2 with q = q m and η i − η i−m 2 will decrease, impacting the total values of the majorant.

Estimates of errors generated by the discretization
Before deriving estimates of the approximation errors appearing in the contractive iterative scheme, we need to study discretization errors encompassed in (3.4)-(3.5) for the ith iteration. Henceforth, the pair (u, p) i = (u i , p i ) is considered as the exact solution of (3.4)-(3.5), whereas (u, p) i h = (u i h , p i h ) denotes its approximation computed by a certain discretization method. We aim to derive computable and reliable estimates of the error measured in the terms |||e i p ||| 2 p and |||e i u ||| 2 u .
Majorant of the error in the pressure term. For the first equation (3.4), Lemma 3 presents a computable upper bound of the difference e i p := p i − p i h between the exact solution p i ∈ W 0 and its approximation p i h ∈ W 0 , measured in terms of the energy norm ||| e i p ||| 2 p .
Lemma 3. For any p i h ∈ W 0 , any auxiliary vector-valued function

1)
and any parameter ζ ≥ 0, we have the following estimate Here, where g is defined in (2.6), and is defined via the constant in the trace-type inequality and positive parameters of the Biot model β and L.
Next, we set w = e i p and introduce an auxiliary function z i h ∈ H Σ p N (Ω, div) (cf. (4.1)) satisfying the identity Using the Hölder and Young inequalities, the first term on the right-hand side of (4.5) can be estimated as The second term on the right-hand side of (4.5) is bounded analogously, i.e., where C p Ω (cf. (4.3)) is a constant in the inequality defined in (4.4). By summing up the results of (4.6) and (4.7), we obtain At this point, the term r eq (p i h , z i h ) 2 is not fully computable in the usual sense of functional majorants since it is defined using η i−1 . However, it can be estimated by (without the squares) and apply an approach similar to that was used in the proof of Lemma 1 This yields that The combination of (4.9) and (4.11) yields Remark 7. The numerical reconstruction of the majorant involves several steps. They are motivated by the accuracy requirements imposed on the upper bound of the error. To generate guaranteed bounds with the realistic efficiency index h is approximated by the chosen discretization method recovering the exact solution of (3.4)). However, to obtain the sharpest estimate, functional M p must be optimized w.r.t. z i h and ζ iteratively. This generates an auxiliary variational problem w.r.t. vector-valued function z i h . Alternatively, one can consider the mixed formulation of (3.4) and reconstruct the pair (p i h , z i h ) simultaneously using one of the well-developed mixed methods [69,70]. Then, both variables required for the reconstruction of M p are directly computable, and no additional post-processing (computational overhead) is required.
Majorant in Lemma 3 yields an estimate of the e i p measured in terms of L2-norm.
Corollary 2. For any p i h ∈ W 0 , any auxiliary functions and parameters defined in Lemma 3, the estimate is a constant in Friedrichs' inequality (cf. (6.2)), and λ K is the minimum eigenvalue of the permeability tensor (cf. (1.2)).
Majorant of the error in the displacement term. Current section considers estimates for the error between the exact solution u i ∈ V 0 and its respective approximation u i h ∈ V 0 measured in terms of the energy norm ||| · ||| 2 u (cf. (2.10)). Since p i h is, in fact, used instead of p i , the original problem (3.5) is replaced by with a perturbed right-hand side. Therefore, u i h is an approximation ofũ i instead of u i . In other words, e i u is composed of the error arising due to the original problem is replaced by (4.14), i.e., u i −ũ i , and the errorũ i − u i h arising because (4.14) is solved approximately. By means of the triangle inequality, e i u can be estimated by abovedescribed errors as follows: Here, |||ũ i − u i h ||| 2 u can be estimated by the functional majorant for a class of the elasticity problems (see Lemma 4), whereas ||| u i −ũ i ||| 2 u is controlled by the bound following from the difference of model problems (3.5) and (4.14) (see Lemma 5).
Lemma 4. For any u i h ∈ V 0 approximatingũ i in (4.14), any auxiliary tensor-valued function and any parameter ξ ≥ 0, we have the estimate β, α, µ, λ are characteristics of the Biot model, and is defined through the constants C tr Σ u N and C K in trace-type and the Korn first inequalities respectively.
Proof: For the simplicity of exposition, let us assume the following representation of the elasticity tensor L ε(u) := 2 µ ε(u) + λ div(u). (4.20) Then, the derivation of an a posteriori error estimate for the problem follows the lines presented in [54,Section 5.2]. In particular, considering an approximation u i h ∈ V 0 , we subtract bilinear form (Lε(u i h ), ε(v) from the left-and right-hand side of (4.21) and set v =ũ i − u i h to obtain Next, we set v =ũ i − u i h and add the divergence of the tensor-valued function τ i h ∈ [T Div (Ω)] d×d , i.e., into the left-and right-hand side of (4.22), which results in the identity (4.17). By means of the Hölder and Young inequalities, the first term on the right-hand side of (4.22) can be estimated as The second and the third terms are combined and estimated as follows defined through constants C K and C tr Σ u N in the Korn and trace inequalities defined in (6.4) and (6.3), respectively. By choosing parameters α 1 = (ξ + 1), α 2 = (1 + 1 ξ ), where ξ > 0, we arrive at Consider now (4.20) and the representation of tensor L −1 τ i h through Lame parameters, i.e., Then, the first term on the right-hand side of (4.25) can be rewritten as Taking the latter into account, we arrive at the alternative estimate where r d,µ,λ (u i h , τ i h ) in defined in (4.26) and τ i h = τ i h is an auxiliary stress approximating function reconstructed in the correspondence with u i h .
Remark 8. We note the choice of the auxiliary tensor-function providing the optimal values of the error estimate is τ h := Lε(ũ i ) := 2 µε(ũ i ) + λdiv(ũ i )I. In this case, we can show that the equilibration residual r eq (p i h , τ i h ) vanishes, and the dual one provides the exact representation of the error.
Lemma 5 proceeds with estimation of e i u (cf. (4.13)), accounting for the error arising if (3.5) is replaced by (4.14).

Proof:
As it was noted in (4.15), the error is two-folded and composed from u i −ũ i 2 u and ũ i − u i h 2 u , where the second term is controlled by (4.16) in Lemma 4. The estimate of the first term is derived by considering the difference of (3.5) and (4.14), i.e., The latter one can be estimated from above by the Cauchy inequality, which yields By using the Young inequality with χ ≥ 1 2λ , we arrive at According to Lemma 2, the linear combination in (4.29) can be estimated as By using we obtain Combining (4.16) and (4.30), we arrive at In addition to (4.16), we can obtain the estimate for the error measured in terms of div · 2 -norm.
Corollary 3. For any p i h ∈ W 0 , any u i h ∈ V 0 approximatingũ i in (4.14), as well as any parameters and function defined in Lemma 5, we have where M p (p i h , z i h ) and Mũ((u, p) i h , τ i h ) are defined in (4.12) and (4.16) for any z i h ∈ H Σ p N (Ω, div) and τ i h ∈ [T Div (Ω)] d×d , respectively, and µ is characteristic of the Biot model.
Proof: By using inequality (6.5) and substituting it in (4.16), we arrive at

Estimates of errors generated by the iteration method
Next, we consider guaranteed bounds of errors arising in the process of contractive iterations (3.4)-(3.5) applied to the system (2.7)-(2.8).
Error estimates for pressure term. First, we prove the following result for the error in the flow equation, which can be done in two different ways. For Lemma 6, we consider functions p i h , p i−1 h ∈ W 0 as approximations of two consequent pressures associated with the iterations i and i − 1, whereas u i h , u i−1 h ∈ V 0 are approximations ofũ i andũ i−1 ∈ V 0 in (4.14), respectively. From now on, when we refer to both dual variables z i h and τ i h corresponding to the ith iteration step, we refer to them as a pair (τ , z) i h .
Lemma 6. For p i , p i−1 ∈ W 0 approximating p ∈ W 0 in (2.8), the estimate of the error incorporated in the pressure term on the ith iteration step has the following form where M h u,div and M h p,L 2 are defined in Corollaries 2 and 3 with τ i h ∈ [T Div (Ω)] d×d and z i h ∈ H Σ p N (Ω, div), respectively, q = L β+L , and Parameters α, β, λ, µ f , C F Σ p D , λ K , and τ are characteristics of the semi-discrete Biot model (3.2)-(3.1).

Proof:
We begin by noting that for the error p − p i caused by the iterative scheme The estimate of ∇(p − p i ) 2 Kτ follows from (3.8). To proceed forward, we need to estimate the right-hand side of (3.8), namely η i − η i−1 2 . By adding and extracting the discretized approximations η i−1 h and η i h , we obtain Here, the first term η i h − η i−1 h 2 is fully computable, and by means of relation we obtain the estimate for the second and third terms: For simplicity, we exclude parameter γ by substituting γ 2 = 2 L: Therefore, the estimate of ||| p − p i ||| 2 p can be represented as follows Finally, by substituting λ = α 2 2 L in (5.5), we arrive at (5.1).
Corollary 4. Let conditions of Lemma 2 and Lemma 6 hold. Then, for 1 ≤ m ≤ i, the estimate of the error incorporated in the pressure term on the ith iteration step has the alternative form The norm in the left-hand side of (5.3) can be alternatively estimated by the contractive estimate, which results into alternative error bound independent of the functional majorants. This result is presented below.
Lemma 7. For any p i ∈ W 0 approximating p ∈ W 0 in (2.8), the estimate of the error incorporated in the pressure term on the ith iteration step has the following form Proof: We use inequality (5.2), to obtain However, in this proof, we use contraction properties of sequence {δη i }, to obtain where all terms on the right-hand side except η 1 − η 1 where right-hand side contains computable η 0 := α γ divu 0 − L γ p 0 . Combining (5.8), (5.9), and (5.10), we arrive at Error estimates for displacement term. To derive the upper bound for the error u − u i measured in terms of we exploit the idea analogous to the one used to estimate the error in the pressure term.
Lemma 8. For any u i , u i−1 ∈ V 0 approximating u ∈ V 0 in (2.7), the error in the displacement on the ith iteration step has the following form where M h u,div and M h p,L 2 are defined in Corollaries 2 and 3 for Parameters λ, µ, α are characteristics of the Biot model.
where the right-hand side is controlled by the contractive term η − η i 2 , which follows from (3.9). Therefore, we obtain Analogously to the proof of Lemma 6 (cf. (5.4)), the estimate for the second term in (5.14) results into Again, by substituting λ = α 2 2 L in (5.15), we arrive at (5.16).
Corollary 5. Let conditions of Lemma 2 and Lemma 8 hold. Then, for 1 ≤ m ≤ i, the error in the displacement on the ith iteration step has the alternative form Analogously, contraction properties can be exploited to derive an alternative bound of the error u − u i .

Lemma 9.
For any u i ∈ V 0 approximating u ∈ V 0 in (2.7), the error in the displacement on the ith iteration step has the following form  General estimate for the error in the iterative coupling scheme. Derivation of the reliable estimate for the error in the pressure approximation reconstructed on the ith iteration is based on the combination of two different approaches, i.e., estimates for contractive mapping used for iterative methods and functional error estimates. Theorem 2 presents an upper bound of the error in the pressure term, which combines the above-mentioned techniques.
Proof: To decompose the error ||| e p ||| 2 p in two parts, we apply the triangle inequality The first term in the right-hand side of (5.19) is bounded by (5.1) from Lemma 6, whereas the second term is controlled by (4.2) from Lemma 3. Analogously, by using the triangle rule, we obtain The first term in the right-hand side of (5.20) is controlled by (5.16), whereas the estimate of the second term follows from (4.16). where M p and M u are defined in Theorem 2.

Numerical example
The numerical properties of above estimates are explained in the series of the following numerical examples. We take two tests where we begin with a manufactured solution for the pressure and displacement. These tests pertain to different values of the contraction coefficient q in the fixed stress scheme and study the efficiency of the error estimators. For each of the two tests, we consider both parameters that are academic as well as repeat them for more realistic parameters of the material properties (permeability, Lame coefficients etc). We solve the Biot model at each time step using the fixed stress scheme by choosing different mesh sizes and time steps. Moreover, we also vary the finite element pairs used for solving the problem. All these variations, including material properties, discretization parameters, and the finite element pairs allow us to study the performance of the estimators and show their robustness.
Our approach is to compute the approximations of the pressure and displacement for the test cases and study the convergence of this to the exact solution. The errors are computed by employing different norms including the combination of pressure and displacement as well the individual ones in both L 2 and the energy norms. We choose a time and a spatial mesh size and at this discrete time, we fix the number of fixed stress iterations and study the convergence. The majorants M h p and M h u are obtained via minimisations w.r.t. to the auxiliary functions. To characterise the efficiency of all the above mentioned error estimatators, we use the efficiency index defined as I eff := M/|||e||| 2 , where M is a chosen majorant and |||e||| 2 is error measured in the corresponding error norm. We study the individual contributions of the different terms in the majorants, namely the dual term and reliability/equilibration term. The first one allows us to study the efficiency as the error indicators and the small contribution of the second one ensures the reliability of the total error bounds. This is studied for the different majorants, the one for the discretization error and for the iterative ones and for the different choices of error norms.

Example 1. Simplified parameters.
Let Ω := (0, 1) 2 ∈ R 2 , T = 10. The exact solution of (1.1) is defined as The Lame parameters are µ = λ vol = 1, which leads to λ plane = 2 λ vol µ λ vol +2 µ = 2 3 . We set α = β = 1, C F = 1 √ 2 π , and K [mD/cP] is a unit tensor. From the parameters above, it follows that L = α 2 2 (λ+2µ/d) = 0.3 and q = L β+L = 3 13 ≈ 0.23077. Taking into account that the error estimates depend on the term q 2 (1−q 2 ) , which goes to infinity as q goes to 1, the efficiency of the resulting M u and M p drastically depends on the value of q. In this particular example, the ratio q 2 (1−q 2 ) = 0.05625 is relatively small, which prevents M  First, we generate approximations employing 10 time steps (where the total number of steps in time is denoted by N ) with corresponding size of the step τ = 1.0 and spatial mesh-size h = 1 64 using standard P 1 (polynomial/lagrangian first-order) finite elements (see (6.7)). Let I denote the number of iterations to solve (3.4)-(3.5) on each time-step; we consider I = 5 for this discretization. Table 1 is minimized w.r.t. the function z i h and parameter ζ. The results of this optimization procedure are listed in Table  2. We see that only two iterations are enough to achieve good efficiency I eff := M h p /|||e p ||| 2 = 1.0023. In Table 2, we see the decrease of the dual term m 2 as well as the reliability/equilibration term of the majorant m 2 eq := r eq (p In order to achieve the desired efficiency in the error estimate, the reliability term must be several orders of magnitude less than the dual one, which is satisfied in this case. Analogous minimization is performed for the majorant M h u . After minimization, the terms m 2 d,K and m 2 d,λ,µ become very close to the true errors |||e p ||| 2 and |||e u ||| 2 , respectively, and can be used as error indicators for the local error distribution over the computational domain. To confirm that, Figure 1 presents the distribution of the errors and indicators (generated by majorants) w.r.t. to the numbered finite element cells. On the right side of Figure 1, we depict |||e p ||| 2 and m 2 d,K := r d,K (p i h , z i h ) 2 and on the left side the error in u measured in terms of |||e u ||| 2 as well as the indicator m 2 d,µ,λ := r d,µ,λ (u i h , τ i h ) 2 w.r.t. to the numbered finite element cells. To depict the error green marker is used, whereas for the majorant the red color is chosen. We see that the indicator produces quantitatively efficient error indicator. First, let us assume that the error bounds M   Table 3(a). The results are presented corresponding to meshes with two different mesh-sizes h and τ . As the caption highlights, all these values are obtained for the auxiliary functions are reconstructed by the following finite elements z i h ∈ RT 1 and τ i h ∈ [P 2 ] 2×2 , which is an order higher then usual choice of the finite element approximation spaces for the fluxes and stresses in mixed formulations. However, it is important to obtain better accuracy of the error bounds. To emphasize this importance, Table 4 Table 3(c). They illustrate the accumulated values of the errors and majorant over the whole time interval. We see that even with the decreasing τ , which scales the permeability tensor K, the efficiency of total error estimates stays rather robust. We note here, that each increment of M Realistic parameters. Let us assume now more realistic parameters in the above example (similar to those taken from [66] and [59]). Let exact pressure be scaled as follows: p(x, y, t) = 10 8 x (1 − x) y (1 − y) t. The permeability tensor divided by fluid viscosity is taken as K = 100I [mD/cP], fluid compressibility fixed to 4.7 · 10 −7 [psi −1 ], and initial porosity φ 0 is assumed to be 0.2. The Biot and bulk modulus are M = 1.65 · 10 10 [Pa] and E = 0.594 · 10 9 [Pa], respectively. That results into β = 1 M + c f φ 0 = 9.406 · 10 −8 , which, in turn, yields considerable small L = α 2 2 (λ+2µ/d) = 1.35 · 10 −9 . Such a tuning parameter generates instantaneous convergence of iterative scheme with contractive parameter q = L β+L = 6.73 · 10 −12 . The resulting errors and corresponding estimates are summarised in Table 6 (for τ = 0.1 and difference spacial mesh-sizes h). For such q, even one iteration is enough for convergence. However, one can consider two to five iterations to improve the efficiency of the majorant. The efficiency indices obtained here confirm the efficiency of obtained majorants also for parameters closed to those used in engineering applications.
We consider 10 time-steps of the length τ = 1.0 discretizing considered interval and spatial mesh-size h = 1 64 using standard P 1 finite elements (see (6.7)). The number of iterations to solve the problem on each time-step is set to 12. For the time-interval [t 9 , t 10 ] = [9.0, 10.0], the convergence of the errors in u and p is presented by Table 7.
From one side, we can consider iterations I and I − 1 as subsequent ones with a contraction parameter q = 0.9203. However, by using Lemma 2 instead, let m = 6 so that 6th and 12th iterations be treated as two consecutive steps but with q = q 6 = 0.6027. Then, the constants dependent from the ratio q 2 (1−q 2 ) = 2.3012 attain more acceptable values. Table 8 illustrates the improvement of error majorants efficiency indices as the explained-above method is employed. The local distribution of the error in p and u on each cell of the finite-element discretization is presented for mesh sizes h = 1 / 4 and h = 1 / 8 in Figure 2. One can see the resemblance in the local error distribution for p since the exact solutions for both examples are the same. The local values of error and indicator (generated by the majorant) in u has more uniform distribution.
Realistic parameters. Let us consider more realistic parameters similar to those considered in [71]. Again, let exact pressure be p(x, y, t) M + c f φ 0 = 6.8909 · 10 −5 , which, in turn, yields considerable small L = α 2 2 (λ+2µ/d) = 1.1428 · 10 −8 . Such parameter generates contractive parameter q = L β+L = 1.6584 · 10 −4 . The latter means that the ratio q 2 (1−q) 2 is of order 10 −8 , which does not degenerate the efficiency of the majorant. We summarise the convergence results in Table 9. It confirms that even for more realistic parameters (common for engineering applications), the value of the estimates remains rather efficient.

Conclusion
We analyze semi-discrete approximations of the Biot poroelastic problem and deduce guaranteed and fully computable bounds of corresponding errors. The derivation combines estimates for contraction mappings and functional a posteriori error majorants for the elliptic problems. An extended version of the paper that contains an overview of the related works and derivation of the model as well as detailed proofs can be found in [62].       Let us determine the values of parameters c 0 , γ, and L such that the terms in the left-hand side of (6.17) are positive and contraction in δ(η − η h ) i−1 2 is achieved. It follows from δ(η Comparing (6.18) and (6.17), we arrive at the following condition on the free parameters: By using condition γ 2 = 2L, we obtain ε(δ(u − u h ) i ) 2 2µ + q ∇δ(p − p h ) i 2 Kτ + δ(η − η h ) i 2 ≤ q 2 δ(η − η h ) i−1 2 (6. 19) with q = L β+L and L = α 2 2 λ .