A priori error estimates for a discretized poro-elastic–elastic system solved by a fixed-stress algorithm

. We consider a poro-elastic region embedded into an elastic non-porous region. The elastic displacement equations are discretized by a continuous Galerkin scheme, while the ﬂow equations for the pressure in the poro-elastic medium are discretized by either a continuous Galerkin scheme or a mixed scheme. Since the over-all system of equations is very large, a ﬁxed-stress algorithm is used at each time step to decouple the displacement from the ﬂow equations in the poro-elastic region. We prove a priori error estimates for the resulting Galerkin scheme as well as the mixed scheme, with the expected order of accuracy, provided the algorithm is sufﬁciently iterated at each time step. These theoretical results are conﬁrmed by a numerical experiment performed with the mixed scheme. A complete analysis including a posteriori estimates for both the Galerkin and the mixed scheme has been done but is too long to appear here.


Introduction
Starting from the pioneering work of von Terzaghi [1] and Biot [2,3], a growing demand for assessing fluid structure interactions has driven the development of numerical models coupling flow and geomechanics.Here important applications include subsidence events, carbon sequestration, ground water remediation, hydrocarbon production and hydraulic fracturing, enhanced geothermal systems, solid waste disposal, and biomedical heart modeling.In this paper we consider the Biot model that consists of a geomechanics equation coupled to a flow model with the displacement, pressure, and flow velocity as unknowns.In contrast to solving the Biot system fully implicitly, we focus on iterative schemes that allow the decoupling of the flow and mechanics equations and offer several attractive features (such as use of existing flow and mechanics codes, use of appropriate pre-conditioners and solvers for the two models, and ease of implementation).The design of iterative schemes however is an important consideration for an efficient, convergent, and robust algorithm.
Four main iterative coupling schemes for coupling flow with mechanics were introduced by Settari and Mourits [4,5] and implemented and studied by Gai et al. [6,7], Gai [8], Dean et al. [9], Girault et al. [10], Almani et al. [11][12][13][14].They are known as the fixed-stress split, fixed-strain split, undrained-split, and drained split schemes (Mikelic ´and Wheeler [15], Mikelic ´et al. [16] and Kim et al. [17,18]).In the fixed-stress and fixed-strain split schemes, the flow problem is solved first, followed by the mechanics problem, and a constant mean total stress or constant strain is assumed during the flow solve respectively.In contrast, in the drained and undrained split coupling schemes, the mechanics problem is solved first, followed by the flow problem, and a constant fluid pressure and a constant fluid mass (i.e., fluid content of the medium) is assumed during the mechanics solve respectively.Both the fixed-stress split and undrained split schemes are shown to be contractive (Mikelic ´and Wheeler [15], Mikelic ´et al. [16]]), whereas the fixed-strain split and drained split schemes are only conditionally stable (Kim et al. [17,18]).In addition, we note here that these iterative coupling schemes can be interpreted as pre-conditioner techniques for solving the fully coupled system.For instance, the work of Gai et al. [7,19] involves interpreting the fixed-stress split scheme as a physics-based preconditioning strategy applied to a Richardson fixed point iteration.The same preconditioning technique was applied to the fully coupled system in the work of Castelletto et al. [20,21].
Extensions of the fixed-stress split and undrained split iterative schemes to the multirate case, in which flow takes multiple fine time steps within one coarse mechanics time step, can be found in the work of Almani et al. [11] and Kumar et al. [22].In addition, extensions to a nonlinear case can be found in the work of White et al. [23], Borregales et al. [24], and da Silva et al. [25], and a multiscale extension of the fixed-stress split scheme is illustrated in the work of Dana et al. [26].Moreover, the work of Bause et al. [27] and Borregales et al. [28] explored space-time methods of the corresponding iterative coupling schemes, and work of Rodrigo et al. [29] considered the stability analysis of the discretization schemes.
Here we restrict our attention to fixed-stress iterative coupling, analyze both Galerkin and mixed finite element for flow and Galerkin for elasticity.We derive a contraction mapping, stability estimates and a priori error estimates for the discretized problem that incorporates convergence of the iteration at each time step.This result appears to be the first that treats the entire spatial-time system.A priori optimal rates of convergence are obtained when convergence of the iteration is achieved at each time step.Due to preference of mass conservation for flow, mixed finite elements are chosen for implementation.Galerkin approximation can be modified by either post processing, see Odsaeter et al. [30], or selecting Enriched Galerkin, see Lee et al. [31], to obtain mass conservation.
While the results presented here apply to the poroelastic system, a novel feature of this work involves the discretized poro-elastic-elastic system.Such systems arise in coupled flow and poromechanics phenomena representing hydrocarbon production or CO 2 sequestration in deep subsurface reservoirs.Here the spatial domain in which fluid flow occurs is generally much smaller than the spatial domain over which significant deformation occurs.The typical approach is to either impose an overburden pressure directly on the reservoir thus treating it as a coupled problem domain or to model flow on a huge domain with zero permeability cells to mimic the no flow boundary condition on the interface of the reservoir and the surrounding rock.The former approach precludes a study of land subsidence or uplift and further does not mimic the true effect of the overburden on stress sensitive reservoirs whereas the latter approach has huge computational costs.A field example that has motivated this work on poro-elastic-elastic sytems is the In Salah Gas Joint Venture in Algeria, a CO 2 storage project operational from 2004 to 2011, see Iding and Ringrose [32].The Krechba field was part of the In Salah gas field development and was the largest onshore CO 2 storage site in the world.At Krechba, CO 2 from several gas fields was removed from the production stream and injected into the deep saline carboniferous formation away from the producers.The large volume of injected CO 2 caused ground upheaval.Surface deflection was recorded using remote sensing sensors on board of an INSAR (Interferometric Synthetic Aperture Radar) satellite.Here applying a compositional flow and geomechanics simulator up to the surface is extremely computationally costly due to the depth of the aquifer.
To address flow and mechanics modeling in deep subsurface reservoirs, a rigorous derivation of the interface conditions between a poro-elastic medium (the pay-zone) and an elastic body (the non-pay zone) was derived in Mikelic ´& Wheeler [33].Here effective equations represent a two-scale system for three displacements and two pressures, coupled through the interface conditions with the Navier equations for the elastic displacement in the nonpay zone.The following appropriate interface conditions at the interface between an elastic and a poro-elastic medium were determined: (i) the displacement continuity, (ii) the effective normal stress continuity and (iii) the normal Darcy velocity zero from the poro-elastic side.This theoretically supports our computational approach for treating the poro-elasticity-elasticity domains in this paper.
This work is organized as follows.In the subsection below we establish notation.In Section 2, a continuous time model involving the decoupling of the model into elastic and poro-elastic domains with interface conditions is formulated both in primal and mixed form.The primal formulation, complete with the fixed-stress splitting algorithm, is fully discretized in Section 3. Convergence of the algorithm, which invokes stability of the scheme, is also established in Section 3. The a priori error equalities, inequalities, and error estimates are derived in Section 4. Section 5 is devoted to the discretization of the mixed formulation and a short survey of its convergence and numerical analysis.Computational results for the mixed scheme that confirm these error bounds are presented in Section 6.

Notation
In the sequel, we shall use the following functional notation; for the sake of simplicity, we define the spaces in three dimensions.In a region X, C 0 X À Á denotes the space of continuous functions in X.The scalar product of L 2 (X) is denoted by ðÁ; ÁÞ X 8f ; g 2 L 2 ðXÞ; ðf ; gÞ X ¼ and if the domain of integration is clear from the context, we suppress the index X.Let (k 1 , k 2 , k 3 ) denote a triple of non-negative integers, set |k| = k 1 + k 2 + k 3 and define the partial derivative o k by Then, for any non-negative integer m, recall the classical Sobolev space (cf.Adams [34] or Nec ˇas [35]) equipped with the following seminorm and norm (for which it is a Hilbert space) This definition is extended to any real number s ¼ m þ s 0 for an integer m !0 and 0 < s 0 < 1 by defining in dimension d the fractional semi-norm and norm: The reader can refer to Lions and Magenes [36] and Grisvard [37] for properties of these spaces.In particular, we have the following trace property in a domain X with Lipschitz continuous boundary oX: If v belongs to H s (X) for some s 2 1 2 ; 1, then its trace on oX belongs to H sÀ 1  2 ðoXÞ and there exists a constant C such that Finally, if C is a subset of oX with positive measure, |C| > 0, we say that a function g in H where jvj and d(x, C) denotes the distance to C. All the above definitions are directly extended to vector functions, with the same notation.The space H(div, X) is the Hilbert space equipped with the graph norm.The normal trace v Á n of a function v of H(div, X) on oX belongs to H À 1 2 ðoXÞ, the dual space of H 1 2 ðoXÞ, see for instance Girault and Raviart [38].The same result holds when C is a part of oX and is a closed surface.But if C is not a closed surface, then v Á n belongs to the dual of H 1 2 00 ðCÞ.We also recall Korn's inequality, valid for all functions v in H 1 ðXÞ d that vanish on C, for a constant C 1 (X, C) depending only on X and C.
Here e(v) denotes the strain tensor.When combined with Poincare ´'s inequality, also valid for all functions v in H 1 (X) that vanish on C, with another constant C 2 (X, C) depending only on X and C, we recover the full H 1 norm in the left-hand side of (3): We also have the interpolation inequality H 1 ðXÞ : Thus by combining this with (4), (3), and (5), we derive the trace inequality for all functions v in H 1 (X) d that vanish on C As usual, for handling time-dependent problems, it is convenient to consider measurable functions defined on a time interval ]a, b[ with values in a functional space, say X (cf.[36]).More precisely, let ||AE|| X denote the norm of X ; then for any number r, 1 r 1, we define with the usual modification if r = 1.It is a Banach space if X is a Banach space, and for r = 2, it is a Hilbert space if X is a Hilbert space.We denote derivatives with respect to time by o t and we define for instance Similarly, C 0 ð½a; b; X Þ is the space of continuous functions in [a, b] with values in X and 2 Domain and model formulations Let X be a bounded, connected, Lipschitz domain in R d , d = 2, 3, and let X 1 be a connected proper subset of X with a Lipschitz boundary that does not intersect the boundary of X, as in Figure 1.We set We are interested in the situation where a poro-elastic model holds in X 1 (the pay-zone) while an elastic model holds in X 2 (the non-pay zone) which is a non-porous region, see Figure 1.This work extends readily to more general configurations, but for simplicity, we focus on this situation.In reference [33], the model was derived from first principles: 1.In X 1 , the quasi-static Biot's system from consolidation theory.
2. On the boundary of X 1 , the interface conditions between a poro-elastic medium and an elastic medium, i.e., the boundary conditions at a closed outer boundary for the quasi-static Biot system.3.In X 2 , the elasticity equations.
Let C 12 denote the boundary of X 1 and let n 12 be the unit normal on C 12 exterior to X 1 .Let the boundary of X, oX be partitioned into two disjoint open regions, not necessarily connected, but with a finite number of connected components, when d = 3, we suppose that the boundaries of C D and C N are both Lipschitz-continuous.We denote by n X the unit outward normal vector to oX.To simplify, we assume that the measure of C D is positive: |C D | > 0. For the sake of brevity, we present the problem in the threedimensional case.Let f be the body force in X.In the non-pay zone X 2 , the governing equations are those of linear elasticity.Recalling the Cauchy stress tensor for linear elasticity, In the pay-zone X 1 , we use Biot's consolidation model for a linear elastic, homogeneous, isotropic, porous solid saturated with a slightly compressible fluid.The unknowns are the solid's displacement u and the fluid's pressure p.This model is based on a quasi-static assumption, namely it assumes that the material deformation is much slower than the flow rate, and hence the second time derivative of the displacement (i.e., the acceleration) is zero.We refer to [39] for the full derivation of the system of equations, but as an example, let us make precise the fluid's density.The fluid's density q f is related to the pressure by where p r is a reference pressure, q f,r > 0 a constant reference density relative to p r , and c f the fluid compressibility.But the fluid compressibility c f is assumed to be small (e.g. of the order of 10 À8 or 10 À9 ), and c f (p À p r ) is also small.Therefore this term is neglected and we replace q f by the constant q f,r .With this simplification, similar linearization, and arguing as in [39], we obtain the following system of equations a.e. in X 1 • ]0,T [ This system is complemented by an initial condition Here k > 0 and G > 0 are the Lame ´coefficients, which for simplicity are assumed to be constant, a > 0 is the Biot-Willis constant, which is usually around one, l f is the fluid's viscosity assumed to be constant, u 0 is the initial porosity, M > 0 is the Biot modulus satisfying, see Coussy [40], where K b is the drained bulk modulus, g is the gravitation constant, g is the signed distance in the vertical direction, q is a volumetric fluid source term, and K is the permeability tensor, assumed to be symmetric, uniformly bounded, and uniformly positive definite, i.e., each eigenvalue k i of K is real and there exist two constants k min > 0 and k max > 0 such that a:e: x 2 X 1 ; Strictly speaking, the initial condition should be given as However, in practice, the pressure is either measured or computed through a hydrostatic assumption and the initial displacement is computed satisfying (10).When the data are sufficiently smooth, as stated below, initializing the pressure is sufficient to determine the solution.
As the boundary of X 1 is reduced to C 12 , that is entirely contained in X, the only boundary conditions on C 12 are in fact transmission conditions between the two regions.The equations themselves are steady, but the unknowns depend on time through the transmission conditions that are proposed below.On the boundary oX, we prescribe mixed boundary conditions: For the transmission conditions on the interface, we recall the definition of jump through C 12 in the direction of the normal, for any smooth enough function v, Then, we prescribe the following transmission conditions, see [33,41]: there is no flow at the interface Concerning data regularity, we assume that f belongs to , and p 0 belongs to H 1 (X).These are not the most general data assumptions, but they are a convenient simplification for the applications we have in mind.Note that both f and t N are continuous in time and therefore f(0) and t N (0) are well-defined.Finally, the mean stress " r which will be used in the algorithm is defined by It is the hydrostatic stress borne by the composite.More precisely, the product of k and the divergence of the displacement is borne by the skeletal solid, and the product of the Biot-Willis constant and the pore pressure is borne by the pore fluid.The negative sign in front of the pore pressure is due to the sign convention where tensile stresses are considered positive.

Mixed variational formulation
We retain the displacement equation but use a mixed formulation for the flow because it leads to locally conservative discrete schemes.Although p is continuous, the mixed formulation relaxes its regularity and allows us to take p in L 1 ð0; T ; L 2 ðX 1 ÞÞ.We associate with the pressure in X 1 an auxiliary Darcy velocity z defined by and the space for the reservoir velocity is normed by With the same data regularity, the mixed variational formulation reads: Find u 2 L 1 ð0; T ; V Þ, p 2 L 1 ð0; T ; L 2 ðX 1 ÞÞ, and z 2 L 2 ð0; T ; ZÞ, such that subject to the initial condition (11):

Continuous Galerkin approximation and algorithm
In this and the next sections, the theory is developed for the primal formulation.To simplify the presentation, we describe here the case d ¼ 3, but everything carries over to d ¼ 2. Let P k denote the space of polynomials of total degree less than or equal to k in d variables.From now on, we assume that X is a Lipschitz polyhedron and that the boundaries of C D and C N are polygonal.Let T h be a regular family of conforming triangulations of X consisting of tetrahedra, E, of maximum diameter h, and such that no element adjacent to oX intersects both C D and C N and the interior of no element intersects C 12 .We denote by T h;1 the restriction of T h to X 1 .The restriction to simplicial elements is only a matter of convenience; the analysis below extends readily to other shapes such as hexahedra or prisms.As usual, we assume that these triangulations are regular in the sense of Ciarlet [43]: There exists a constant j > 0, independent of h, such that where h E h denotes the diameter of E and q E denotes the diameter of the ball inscribed in E. On this triangulation, we introduce the standard finite element spaces and we set As the exact solution may possibly be not very smooth, it is approximated by Scott & Zhang interpolants (see [44]), say Considering the degree of the polynomial functions in X h and M h , these interpolants have the following quasi-local approximation errors: with constants C independent of E and h, where Á E is a small patch of elements including E containing the values used in computing the approximation.
Regarding approximation in time, we divide the interval ½0; T into N equal subintervals with length Át and set t n ¼ nÁt.The choice of equal time steps is also a simplification; the material below extends readily to variable time steps.For the data, we set With these spaces, the fully discrete split problem is: and march in time, set n :¼ n + 1 and return to step (1).
Considering the uniform positive definiteness of the permeability tensor K and the unique solvability of ( 20) and (21) for given p, it is clear that this algorithm generates a unique sequence.

Stability of the scheme and convergence of the algorithm
From the work of Mikelic ´and Wheeler in [45], see also (Girault et al. [10]), it can be shown that the fixed-stress algorithm ( 32)-( 37) produces a contracting sequence.More precisely, let Then, we have for all n ! 1 and all ' ! 2, and since bk > 1 and is independent of n and h, (40) implies that the sequence of the differences between two consecutive iterates of " r n;' h is contracting, uniformly in n and h.As a consequence, we have which implies convergence for fixed h, Á t, and n.

Initial mean stress difference at each time step
It follows from ( 41) that convergence of the algorithm, uniform in time and space, depends on suitable bounds for " r n;1 h À " r nÀ1 h .We shall see below that this relies on the stability of the scheme, i.e., the uniform boundedness of the sequence of discrete solutions.Recall that Let us start with the second term.The following proposition shows that a bound for the second term reduces to a bound for the first term.Proposition 3.1.We have for all n, 1 n N, ) is the product of the constants of (3) and ( 4), and C is the constant of (6).
Proof.The equation satisfied by where r is defined in (7).By testing this equation with u n;1 h À u nÀ1 h and applying (3), (4), and (6), we derive Then (42) follows by suitable applications of Young's inequality.
The next proposition treats the first term.The result is easily derived by testing (35) when ' ¼ 1 with h h ¼ p n;1 h À p nÀ1 h , adding and subtracting the term ðKr p nÀ1 h ; r h h Þ X 1 and suitably applying Young's inequality.
Proposition 3.2.We have for all n, 1 n N , where b is defined by (39).Hence a bound for " r n;1 h À " r nÀ1 h readily follows by combining Propositions 3.1 and 3.2.Lemma 3.1.We have for all n, 1 n N, with the constants C and C of Proposition 3.1, In view of the assumptions on the data, it follows that jj" r n;1 h À " is also bounded.This raises the issue of stability.

Stability of the scheme
Here we investigate the uniform boundedness of the sequences ðu n h ; p n h Þ generated by the algorithm in (11), i.e., when the threshold is met.Thus the displacement equation ( 36) is written (without the superscript '), while the flow equation ( 35) written at level ' ¼ ' n can be reformulated as The following proposition is written under the assumption that the time step satisfies Át 1 2 , but the choice 1  2 is purely a matter of convenience.Indeed, as the asymptotic analysis must be written for Dt < 1, otherwise the order of the scheme is meaningless, then 1  2 can be replaced by any number strictly smaller than one.
where I 0 h and D n h are respectively contributions of the initial terms and data, with the constants C and C of Proposition 3.1.
Proof.Since q f ;r gg is a polynomial of degree one (independent of time), " p m h belongs to M h and ( 46) can be tested at time t m with h h ¼ " p m h .Then by testing (45) with Þ, adding the resulting equations and multiplying by 2Dt, we derive, To apply Gronwall's Lemma, it is convenient to write at level n, " Similarly, , the term involving p n h À p nÀ1 h is balanced by the same term in the left-hand side of (23).With another c > 0, we obtain At level m = n À 1, this splitting is not necessary and we simply write for any c 0 > 0, , the total contribution of " p nÀ1 h from levels n and The levels m n À 2 are treated as level n À 1; a more favorable value of c 0 can be chosen, but for the sake of simplicity, it is not used here.These inequalities are summed over m from 1 to n, but the three terms in the last line of (50) are summed by parts because we have no control on 1  Á t ðu m h À u mÀ1 h Þ.Thus, we write in view of (3), (4), and considering (6), The desired result follows by applying Young's inequality to the above terms.
In view of ( 32) and (33), the initial terms of (48) are easily bounded independently of h.If, in addition to the assumptions of Proposition 3.3, f and t N are H 1 in time then the terms of (49) are bounded independently of n, h, and Á t.Therefore, it remains to control the term " r m h À " r m;'mÀ1 h .This is done by combining the contraction property (41) with Lemma 3.1 and it yields All terms above are easily bounded except for the one involving K 1 2 rðp mÀ1 h À q f ;r ggÞ because it lacks a factor Dt in order to be controlled by the left-hand side of (47).To recover a factor Dt, we use the contraction and assume that the number of iterations is sufficiently large.This implies stability of the scheme.
then there exists a constant C independent of n, h, and Dt, such that Proof.
The assumption (55) implies that the factor of K Therefore this term is controlled by the same term in the left-hand side of (47).With the above regularity of the data, an application of Gronwall's Lemma yields a uniform bound for p n h and eðu n h Þ.In turn, this bound permits to estimate r Á ðu n h Þ and , uniformly with respect to h, n, and Át.Finally, the regularity of qgg implies a similar bound for p n h .Remark 3.1.The exponential factor in the right-hand side of (56) can be avoided by using L 1 À L 1 bounds in time such as This is sound, as long as it is done with data, but the same strategy applied to leads to a relation between ' m and Dt that is less favorable than that given by (55).

A priori error estimates
As expected, the above stability proof will be extended to derive a priori error estimates.As the main ingredient is the contribution of the algorithmic error, we begin by revisiting this error.Throughout this section, we use the notation (31) and for convenience we denote by d the finite difference in time for any function f,

The algorithmic error
Let us start with an arbitrary value of '.Recall that the contraction property yields (41), where We have seen in Proposition 3.1 that with the constants C and C of (42).The next proposition sharpens the bound for the pressure in Proposition 3.2.Below, we use the notation for any function q, " q m ¼ q m À q f ;r gg; ð60Þ P h denotes an arbitrary space approximation operator that takes its values in M h and for any function q in L 1 (0, T), m(q) is the average of q in time, Proposition 4.1.Assuming that the solution is sufficiently smooth in time, we have for all n, 1 n N, where b is defined by (39).
V. Girault  Proof.By proceeding as in Proposition 3.2, subtracting Kr P h ð" p nÀ1 Þ, and Kr mð" pÞ, and recalling that q f,r gg is independent of time, we first write By integrating the flow equation ( 21) over]t n À 1 , t n [, the end of the last line above can be written as Then (62) follows by suitable applications of Young's inequality.

Error equation
In addition to P h , we shall apply to the displacement an arbitrary approximation operator R h that takes its values in V h .Let us integrate the flow equation ( 21) over ]t n À 1 , t n [, tested with h h , and divide both sides by Dt; with the notation (58), (60), and (61), this gives By subtracting (35) at level ' ¼ ' n from this equation (see ( 46)), we obtain for all h h 2 M h , The final error equation for flow is derived by adding and subtracting approximations of p and u, and multiplying all terms by Dt.This yields Now, we turn to the displacement and write (20) tested with v h at time t n .We obtain By observing that p n À p n h ¼ " p n À " p n h , the difference between this equation and (36) again with ' ¼ ' n reads It remains to add and subtract approximations of p and u, to derive the final error displacement equation for all ðqðsÞ À q n ; " e n p Þ X 1 ds !:

Error inequality
To simplify, set By comparing (67) at level m with (50), we see that their left-hand sides have much the same structure.Therefore, to deduce an error inequality from (67), it suffices to suitably group the terms in its right-hand side so as to associate them with corresponding terms in the right-hand side of (50).First, the term Next, we can write Thus, this term can be treated as Hence, proceeding as in the proof of Proposition 3.3, assuming again that Át 1 2 , and summing over m, from 1 to n, the part of the right-hand side of (67) corresponding to these two terms is bounded by The corresponding first group of terms in the left-hand side is Next, we consider the terms involving eðdðe n u ÞÞ and r Á dðe n u Þ that must be summed by parts.For the first term, we write Young's inequality gives The contribution of this term to the right-hand side is and the corresponding group of terms that remain in the left-hand side is To simplify the treatment of the terms involving r Á dðe n u Þ, it is convenient to extend a by zero outside X 1 .With this convention, we write Then we deduce by the above argument The contribution of this term to the right-hand side is and the corresponding group of terms that remain in the left-hand side is There remain the two terms involving the gradient of " e n p ; with the notation (68), they can be written as ðKrð" p n À " pðsÞÞ; r " e n p Þ X 1 ds By Young's inequality, when summed over m, this term has the bound for any c > 0, Thus the contribution of this term to the right-hand side is and the corresponding term that remain in the left-hand side is By collecting the terms L i in the left-hand side and R i in the right-hand side, for 1 i 4, we deduce the following error inequality (compare with Proposition 3.3): Proposition 4.2.Let q belong to C 0 ð0; T ; L 2 ðX 1 ÞÞ, u to C 0 ð0; T ; V Þ, and p to C 0 ð0; T ; H 1 ðX 1 ÞÞ.Assume that D t 1 2 and that a is extended by zero outside X 1 .With the notation (60), (66), and (68), the following bound holds for all n, 1 n N, and any c 2 ]0, 2[, jjdð" e m p Þjj where E 0 h and A n h are respectively contributions of the initial errors and interpolation errors in space and time,

Error estimate
With the notation (66) and ( 68), (63) reads The factor of the term involving K and this factor must be controlled in the left-hand side of (77) by 1 with a parameter c 2 ]0,2[ to be chosen.If we choose for instance c ¼ 1 2 , a sufficient condition for this control is which resembles formula (57) that guarantees stability.However, in contrast to the situation of Theorem 3.1, this assumption does not guarantee a satisfactory error bound.Indeed, with this assumption, the first term in the righthand side of (80) is bounded by Considering that this assumption implies which is not sufficient to guarantee a first order in time that requires a factor (Dt) 2 .The second and last terms of (80) are plagued by the same shortcoming.Note that they reflect the inconsistency of the first iteration of the algorithm at each time step.These considerations suggest to replace the above variant of (57) by the stronger condition, With this assumption, the following error estimate holds: Theorem 4.1.In addition to the assumptions and notation of Proposition 4.2, we suppose that the data satisfy , and o t q 2 L 2 (X 1 • ]0, T[ ), that the solution satisfies ot p 2 L 2 (0, T; H 1 (X 1 )) and otu 2 L 2 (0, T; H 1 (X ) d ), and that the assumption (81) on the number of iterations holds.Then, we have the following error estimate for all n, 1 n N, where Proof.By applying the second part of (81) to the first, second, and last terms of the right-hand side of (80) and the first part to the remaining terms, we derive Now, we substitute (86) into (77) with the choice c ¼ 1 2 .The first sum in the second line of the above right-hand side cancels with the corresponding sum in the left-hand side of (77) with the exception of the first and last terms.More precisely, there remains a term (that is neglected for the moment) in the left-hand side 3 2l f Á tjjK and a term in the right-hand side that is included in the initial error, All other remaining terms, apart from the approximation errors, contain time errors and interpolation errors.The latter are left as such since they depend on the degree of the finite element functions and the regularity of the solution.
For the time errors, we observe that for any function f in L 1 (0, T), Hence, grouping terms and applying (87) for the time errors we deduce the following intermediate result: Then (82) follows from Gronwall's lemma.
Note that (82) does not address the error on the gradient of the pressure; it has been eliminated to decrease the restriction on the number of iterations.But this error is an easy consequence of the error equality (67) and the bounds resulting from (82).Thus we obtain with a constant D that is independent of n, Dt and h, There remains to take into account the interpolation error in space and the initial errors.Suppose that o t p belongs to L 2 ð0; T ; H s 1 ðX 1 ÞÞ and otu belongs to L 2 ð0; T ; is the degree of the finite element polynomials of the displacement, respectively the pressure.By taking into account the approximation properties (30) Ch 2ðminðs 2 ;kÞÀ1Þ sup t20;T ½ juj 2 H s 2 ðXÞ ; Finally, we turn to the initial errors.By (32), the initial pressure error vanishes, The initial displacement errors follow from (65) at n = 0, By testing this equation with v h ¼ R h ðu 0 Þ À u 0 h , standard applications of Young's inequality yield Hence the initial error has the same order of convergence as the other errors in space.The next theorem collects these results.
Theorem 4.2.In addition to the assumptions and notation of Proposition 4.2, we suppose that the data satisfy @ t f 2 L 2 ðXÂ0; T ½Þ d , @ t t N 2 L 2 ðC N Â0; T ½Þ d , and @ t q 2 L 2 ðX 1 Â0; T ½Þ, that the solution satisfies p 2 H 1 ð0; T ; H s 1 ðX 1 ÞÞ and u 2 H 1 ð0; T ; H s 2 ðXÞ d Þ, and that the assumption (81) on the number of iterations holds.Then, the solution of the fully discrete scheme (32)- (37) satisfies the error bounds for all n, 1 n N, with a constant C independent of h and Dt.

Mixed approximation and algorithm
The material of Sections 3 and 4 is easily adapted to the mixed formulation ( 25)-( 27), with initial condition (11), and to avoid repetitions, we sketch here the main ideas.
The displacement is discretized in the same spaces X h and V h , but the velocity and pressure spaces are discretized by a standard mixed pair, (Z h , Q h ), where Z h is a finite element subspace of Z with interpolation operator P h and Q h a finite element subspace of L 2 (X 1 ), with interpolation operator q h , (Z h , Q h ) being compatible in the sense that they satisfy a uniform discrete inf-sup condition.In particular, we make the assumption, commonly used in mixed methods, As previously, we denote " p m h ¼ p m h À q f ;r gg.Starting from the fixed-stress splitting algorithm computes u 0 h 2 V h by solving (33), and " r 0 h by (34), " r 0 h ¼ kr Á u 0 h À ap 0 h : (c) " r n;' h by (37) " r n;' h ¼ kr Á u n;' h À ap n;' h ; (d) compute the corrector to difference in fluid content d c u by and march in time, set n :¼ n + 1 and return to step (1).
It is easy to check that this algorithm generates a unique sequence of discrete functions.

Contraction of the algorithm and stability of the mixed scheme
The contraction proof for the mixed scheme is very similar to that for (32)- (38).With the notation (39) for b and dv n;' ¼ v n;' À v n;'À1 ; we derive by the same argument and again we must estimate jjd" r n;1 h jj L 2 ðX 1 Þ .As the discrete displacement equation ( 36) is unchanged, formula (42) is also valid here and we must derive the analogue of (43).This is readily done by testing (95) at level ' ¼ 1 with h and using the velocity-pressure equation to eliminate the divergence in both terms of the above right-hand side.This gives Thus, we have the analogue of Proposition 3.2, with the same notation.
Proposition 5.1.We have for all n, 1 n N, Then a combination of ( 42) and (99) yields the analogue of (44), which again relies on the stability of the scheme.Now, by using again the velocity-pressure equation, so that ðr Á z n it is easy to check that all steps of the proof of the stability Proposition 3.3 carry over to the mixed scheme without modification, and we have the following proposition.
Proposition 5.2.Under the assumptions and notation of Proposition 3.3, the following bound holds for all n, 1 n N, With this result, by reproducing the argument at the end of Section 3.1, we derive stability of the mixed scheme.
Theorem 5.1.Under the assumptions on the data and the number of iterations (55) of Theorem 3.1, there exists a constant C independent of n, h, and Dt, such that For the first term, since z nÀ1 h À P h ðz nÀ1 Þ 2 Z h , we apply the discrete velocity-pressure relation, whereas the other two terms are left unchanged in order to be combined with q n .Thus, we obtain the next proposition.Proposition 5.3.Assuming that the solution is sufficiently smooth in space and time, we have for all n, 1 n N, À P h ðz nÀ1 ÞÞ in the righthand side of (104) is smaller than its counter part in (36) by a factor of 1   2   .In contrast, the factors of the other terms in the right-hand side involving the interpolation errors are slightly larger, but this does not change their order of convergence.
Next, we turn to the error equation.Its derivation follows the lines of Section 4.2, up to the error of the velocity-pressure relation in which we make use of the assumption (92).Thus, we write and hence the difference between the exact and discrete velocity-pressure equations reads This is an important simplification because it will eliminate the factor r Á f h from the error equation.In addition to the notation e n u and a n u of (66), we set e n p ¼ q h p n ð Þ À p n h ; e n z ¼ P h z n ð Þ À z n h ; a n p ¼ q h p n ð Þ À p n ; a n Then (67) is replaced by 1 M ðjj" e n p jj 2 L 2 ðX 1 Þ À jj" e nÀ1 p jj and slightly different interpolation error Ãn h , of the same order, and provided the solution and data are sufficiently smooth.For the sake of brevity the details are skipped, and (77) is replaced by 1 M jj" e n p jj 2 jjdð" e m p Þjj 2 From here, by comparing with the material of Section 4.4, we easily verify that the error of the mixed scheme has the same order as that of the continuous Galerkin scheme, provided the solution, including the Darcy velocity, is smooth enough.The sufficient condition on the number of iterations (81) is slightly improved.Owing to the gain of the

A numerical result
The problem considered here is that of flow to a rate specified injection well in a confined compressible aquifer of thickness H as shown in Figure 2. The analytical solution for the pore pressure in the pay-zone (or aquifer or flow domain) is given as (Verruijt [46]) where k is the permeability, Q is the injection rate, H is the aquifer thickness, r is the radial coordinate measured from the center of the well, t is the time, E 1 (x) is the exponential integral given by and c is the diffusivity coefficient given by where K b ¼ E 3ð1À2mÞ is the drained bulk modulus and G ¼ E 2ð1þmÞ is the shear modulus.The underlying assumptions in the development of (110) are that there are no horizontal deformations in the aquifer and that the total vertical stress remains constant during the development of the hydrological process.We employ the parameters given in Table 1.The flow domain (aquifer) is at a depth of 244 m with the mechanics domain extending all the way to the traction free surface.The aquifer thickness is 61 m.The lateral extents of both the flow and mechanics domains are 915 m.An underburden with a thickness of 61 m is also provided.No flow boundary conditions are imposed on the pay-zone.Gravity is ignored.Denoting the top surface of the mechanics domain as C T , the boundary conditions for the mechanics subproblem are t ¼ 0 on C T ; u Á n ¼ 0 on oX=C T : The Biot parameter is non-zero only inside the pay-zone and zero outside the pay-zone and formally given as & Since for this computational example the term q represents a well and is not in L 1 ðL 2 ðXÞÞ, we compute estimates in L 1 ðL 2 ðX Ã ÞÞ where X Ã ¼ X À B R .Here B R is a cylinder of radius 92 m (which is 1=10 times the areal size of the domain) and height of the pay-zone.The error norm for the pore pressure in the pay-zone is computed using the midpoint quadrature rule:  where m e is the center of mass of element E, s is time, p(s,m e ) is the analytical solution for pore pressure at time s and position m e , p h (s, m e ) is the numerical solution for pore pressure at m e and T is the total time.As shown in Tables 2 and 3, we observe first order convergence for the pore pressure.The simulation was performed using the in-house reservoir simulator IPARS (http://csm.ices.utexas.edu/IparsWeb/index.html) on the bevo3 supercomputer at the Institute for Computational Engineering and Sciences at the University of Texas at Austin.The flow system is solved using a multipoint flux mixed finite element method (Ingram et al. [47]) and the poromechanics system is solved using a conforming Galerkin method.The solver for both the flow and poromechanics systems is the parallel high performance preconditioners library called Hypre (Falgout and Yang [48]).The convergence tolerance is 1 • 10 À8 .The code takes 1-3 coupling iterations to convergence at every time step.
Table 2. Order of convergence of pore pressure solution with time step of 0.01 day.
T h Nature Graded mesh

Fig. 2 .
Fig. 2. Schematic for single well in an infinite confined aquifer problem.
V.Girault et al.:Oil & Gas Science and Technology -Rev.IFP Energies nouvelles 74,24 (2019) of P h and (29) of R h and the regularity (28) of the triangulation, we deduce, with various constants C independent of n, Dt, and h, jjo t a p jj 2 L 2 ðX 1 Â0;T ½Þ Ch 2 minðs 1 ;mÞ jjo t pjj 2 L 2 ð0;T ;H s 1 ðX 1 ÞÞ ; jjeðo t a u Þjj 2 L 2 ðXÂ0;T ½Þ þ jjr Á ðo t a u Þjj .2 A priori error estimates for the mixed scheme Let us start with the algorithmic error.We retain the notation of Section 4.1, but we replace the approximation operator P h by the operator P h that takes its values in Z h .Proposition 4.1 is replaced by a slightly more complex result.Indeed, by proceeding as in Proposition 3.2, we intro-

Table 3 .
Order of convergence of pore pressure solution with time step of 0.005 day.

Table 1 .
Parameters for single well in an infinite confined acquifer problem.