Regular Article
A priori error estimates for a discretized poroelastic–elastic system solved by a fixedstress algorithm
^{1}
Sorbonne Universités, UPMC Univ. Paris 06, CNRS UMR 7598, Laboratoire JacquesLouis Lions, 75005, Paris, France
^{2}
Center for Subsurface Modeling (CSM), Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin, Austin, TX 78712, USA
^{3}
ARAMCO, Saudi Arabia
^{*} Corresponding author: girault@ann.jussieu.fr
Received:
8
March
2018
Accepted:
24
September
2018
We consider a poroelastic region embedded into an elastic nonporous region. The elastic displacement equations are discretized by a continuous Galerkin scheme, while the flow equations for the pressure in the poroelastic medium are discretized by either a continuous Galerkin scheme or a mixed scheme. Since the overall system of equations is very large, a fixedstress algorithm is used at each time step to decouple the displacement from the flow equations in the poroelastic 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 sufficiently iterated at each time step. These theoretical results are confirmed 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.
© V. Girault et al., published by IFP Energies nouvelles, 2019
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 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 preconditioners 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–14]. They are known as the fixedstress split, fixedstrain split, undrainedsplit, and drained split schemes (Mikelić and Wheeler [15], Mikelić et al. [16] and Kim et al. [17, 18]). In the fixedstress and fixedstrain 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 fixedstress split and undrained split schemes are shown to be contractive (Mikelić and Wheeler [15], Mikelić et al. [16]]), whereas the fixedstrain 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 preconditioner techniques for solving the fully coupled system. For instance, the work of Gai et al. [7, 19] involves interpreting the fixedstress split scheme as a physicsbased 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 fixedstress 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 fixedstress 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 spacetime 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 fixedstress 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 spatialtime 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 poroelastic–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 poroelastic–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 poroelastic medium (the payzone) and an elastic body (the nonpay zone) was derived in Mikelić & Wheeler [33]. Here effective equations represent a twoscale 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 poroelastic medium were determined: (i) the displacement continuity, (ii) the effective normal stress continuity and (iii) the normal Darcy velocity zero from the poroelastic side. This theoretically supports our computational approach for treating the poroelasticity–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 poroelastic domains with interface conditions is formulated both in primal and mixed form. The primal formulation, complete with the fixedstress 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.
1.1 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 Ω, denotes the space of continuous functions in . The scalar product of L^{2}(Ω) is denoted by and if the domain of integration is clear from the context, we suppress the index . Let (k _{1}, k _{2}, k _{3}) denote a triple of nonnegative integers, set k = k _{1} + k _{2} + k _{3} and define the partial derivative ∂^{ k } byThen, for any nonnegative integer m, recall the classical Sobolev space (cf. Adams [34] or Neč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' for an integer m ≥ 0 and 0<s'<1 by defining in dimension the fractional seminorm 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 Ω with Lipschitz continuous boundary ∂Ω: If v belongs to H^{ s }(Ω) for some , then its trace on ∂Ω belongs to and there exists a constant C such thatFinally, if Γ is a subset of ∂Ω with positive measure, Γ > 0, we say that a function g in belongs to if its extension by zero to belongs to . It is a proper subspace of , and is normed by(1)whereand d( x , Γ) denotes the distance to Γ.
All the above definitions are directly extended to vector functions, with the same notation. The space H(div, Ω) is the Hilbert space(2)equipped with the graph norm. The normal trace v · n of a function v of H(div, Ω) on ∂Ω belongs to , the dual space of , see for instance Girault and Raviart [38]. The same result holds when Γ is a part of ∂Ω and is a closed surface. But if Γ is not a closed surface, then v · n belongs to the dual of . We also recall Korn’s inequality, valid for all functions v in that vanish on Γ,(3)for a constant C _{1}(Ω, Γ) depending only on Ω and Γ. Here ε ( v ) denotes the strain tensor. When combined with Poincaré’s inequality, also valid for all functions v in H^{1}(Ω) that vanish on Γ, with another constant C _{2}(Ω, Γ) depending only on Ω and Γ,(4)we recover the full H^{1} norm in the lefthand side of (3):(5)
We also have the interpolation inequalityThus by combining this with (4), (3), and (5), we derive the trace inequality for all functions v in H^{1}(Ω)^{ d } that vanish on Γ(6)
As usual, for handling timedependent 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 · _{ X } denote the norm of ; then for any number r, 1 ≤ r ≤ ∞, we defineequipped with the normwith the usual modification if r = ∞. 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 ∂ _{ t } and we define for instanceSimilarly, is the space of continuous functions in [a, b] with values in X and
2 Domain and model formulations
Let Ω be a bounded, connected, Lipschitz domain in , d = 2, 3, and let Ω_{1} be a connected proper subset of Ω with a Lipschitz boundary that does not intersect the boundary of Ω, as in Figure 1. We setWe are interested in the situation where a poroelastic model holds in Ω_{1} (the payzone) while an elastic model holds in Ω_{2} (the nonpay zone) which is a nonporous 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:

In Ω_{1}, the quasistatic Biot’s system from consolidation theory.

On the boundary of Ω_{1}, the interface conditions between a poroelastic medium and an elastic medium, i.e., the boundary conditions at a closed outer boundary for the quasistatic Biot system.

In Ω_{2}, the elasticity equations.
Fig. 1
Payzone with surrounding rock. 
Let Γ_{12} denote the boundary of Ω_{1} and let n _{12} be the unit normal on Γ_{12} exterior to Ω_{1}. Let the boundary of , 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 Γ_{ D } and Γ_{ N } are both Lipschitzcontinuous. We denote by n _{Ω} the unit outward normal vector to ∂Ω. To simplify, we assume that the measure of Γ_{ D } is positive: Γ_{ D } > 0. For the sake of brevity, we present the problem in the threedimensional case. Let f be the body force in Ω. In the nonpay zone Ω_{2}, the governing equations are those of linear elasticity. Recalling the Cauchy stress tensor for linear elasticity,(7)we have a.e. in Ω_{2} × ]0, T[(8)
In the payzone Ω_{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 quasistatic 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 ρ _{ f } is related to the pressure bywhere p _{ r } is a reference pressure, ρ _{ 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 ρ _{ f } by the constant ρ _{ f,r }. With this simplification, similar linearization, and arguing as in [39], we obtain the following system of equations a.e. in Ω_{1} × ]0, T[(9) (10)
This system is complemented by an initial condition(11)Here λ > 0 and G > 0 are the Lamé coefficients, which for simplicity are assumed to be constant, α > 0 is the BiotWillis constant, which is usually around one, μ _{ f } is the fluid’s viscosity assumed to be constant, φ _{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, η 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 λ _{ i } of K is real and there exist two constants λ _{min} > 0 and λ _{max} > 0 such that(12)Strictly speaking, the initial condition should be given asHowever, 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 Ω_{1} is reduced to Γ_{12}, that is entirely contained in Ω, the only boundary conditions on Γ_{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 ∂Ω, we prescribe mixed boundary conditions:(13)
For the transmission conditions on the interface, we recall the definition of jump through Γ_{12} in the direction of the normal, for any smooth enough function v,(14)Then, we prescribe the following transmission conditions, see [33, 41]:(15) (16) (17)Concerning data regularity, we assume that f belongs to H^{1}(0, T; L^{2}(Ω)^{ d }), q belongs to L^{2}(Ω_{1} × ]0, T[), t _{ N } is in H^{1}(0, T; L^{2}(Γ_{ N })^{ d }), and p _{0} belongs to H^{1}(Ω). 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 welldefined.
Finally, the mean stress which will be used in the algorithm is defined by(18)It is the hydrostatic stress borne by the composite. More precisely, the product of λ and the divergence of the displacement is borne by the skeletal solid, and the product of the BiotWillis 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.
2.1 Primal variational formulation
Let us introduce the space(19)As shown in reference [42], problem (9)–(13), (15)–(17) has the equivalent variational formulation: Find u in L^{∞}(0, T; V ) and p in L^{∞}(0,T; L^{2}(Ω_{1})) ∩ L^{2}(0, T; H^{1}(Ω_{1})) solving(20) (21) We observe that the transmission conditions are such that no jump appears in the variational formulation.
2.2 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 is continuous, the mixed formulation relaxes its regularity and allows us to take in . We associate with the pressure in an auxiliary Darcy velocity z defined by(22)and the space for the reservoir velocity is(23)normed by(24)With the same data regularity, the mixed variational formulation reads: Find , , and , such that(25) (26) (27)subject to the initial condition (11):
3 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 , but everything carries over to . Let denote the space of polynomials of total degree less than or equal to in variables. From now on, we assume that is a Lipschitz polyhedron and that the boundaries of and are polygonal. Let be a regular family of conforming triangulations of consisting of tetrahedra, , of maximum diameter , and such that no element adjacent to intersects both and and the interior of no element intersects . We denote by the restriction of to . 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 , independent of , such that(28)where denotes the diameter of and ρ _{ E } denotes the diameter of the ball inscribed in . On this triangulation, we introduce the standard finite element spaces and we setAs the exact solution may possibly be not very smooth, it is approximated by Scott and Zhang interpolants (see [44]), sayConsidering the degree of the polynomial functions in and , these interpolants have the following quasilocal approximation errors:(29) (30)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(31)With these spaces, the fully discrete split problem is:
Initialization. Set(32)Compute and by solving(33)and setting(34)
Time step. n ≥ 1.

Set , , and .

For , compute
(b) compute the predictor of the difference in fluid content by(36)
(d) compute the corrector to difference in fluid content by
If threshold, increment and return to (a);
if threshold, set(38)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 , it is clear that this algorithm generates a unique sequence.
3.1 Stability of the scheme and convergence of the algorithm
From the work of Mikelić and Wheeler in [45], see also (Girault et al. [10]), it can be shown that the fixedstress algorithm (32)–(37) produces a contracting sequence. More precisely, let(39)
Then, we have for all and all ,(40)and since βλ > 1 and is independent of n and h, (40) implies that the sequence of the differences between two consecutive iterates of is contracting, uniformly in n and h. As a consequence, we have(41)which implies convergence for fixed , , and .
3.1.1 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 . 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 thatLet 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,(42) where C = C _{ 1 } (Ω, Γ)C _{ 2 } (Ω, Γ) is the product of the constants of (3) and (4), and is the constant of (6).
Proof. The equation satisfied by iswhere σ is defined in (7). By testing this equation with and applying (3), (4), and (6), we deriveThen (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 with , adding and subtracting the term and suitably applying Young’s inequality.
Proposition 3.2.
We have for all , ,(43) where β is defined by (39).
Hence a bound for 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 of Proposition 3.1 , (44)
In view of the assumptions on the data, it follows that is bounded provided is also bounded. This raises the issue of stability.
3.1.2 Stability of the scheme
Here we investigate the uniform boundedness of the sequences generated by the algorithm in (11), i.e., when the threshold is met. Thus the displacement equation (36) is written (without the superscript ),(45)while the flow equation (35) written at level can be reformulated as(46)
The following proposition is written under the assumption that the time step satisfies , but the choice is purely a matter of convenience. Indeed, as the asymptotic analysis must be written for Δt < 1, otherwise the order of the scheme is meaningless, then can be replaced by any number strictly smaller than one.
Proposition 3.3.
Let belong to , to , and to . Assume that . For any , let . The following bound holds for all , ,(47) where and are respectively contributions of the initial terms and data, (48) (49) with the constants C and of Proposition 3.1 .
Proof. Since ρ _{ f,r } gη is a polynomial of degree one (independent of time), belongs to M _{ h } and (46) can be tested at time t _{ m } with . Then by testing (45) with , adding the resulting equations and multiplying by 2Δt, we derive,(50)To apply Gronwall’s Lemma, it is convenient to write at level , . For any γ > 0, we have(51)Similarly,(52)Choose γ = 1. Since , the term involving is balanced by the same term in the lefthand side of (23). With another γ > 0, we obtain At level m = n − 1, this splitting is not necessary and we simply write for any γ' > 0, By choosing , the total contribution of from levels n and n − 1 isThe levels m ≤ n − 2 are treated as level n − 1; a more favorable value of γ′ 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 . Thus, we write in view of (3), (4), and considering (6),(53) Finally,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 . If, in addition to the assumptions of Proposition 3.3, and are in time then the terms of (49) are bounded independently of , , and Δt. Therefore, it remains to control the term . This is done by combining the contraction property (41) with Lemma 3.1 and it yields(54)All terms above are easily bounded except for the one involving because it lacks a factor Δt in order to be controlled by the lefthand side of (47). To recover a factor Δt, we use the contraction and assume that the number of iterations is sufficiently large. This implies stability of the scheme.
Theorem 3.1.
Let f belong to , to , and to . Assume that . If for each m, 1 ≤ m ≤ N, the number of iterations satisfies (55) then there exists a constant C independent of n, h, and Δt, such that (56)
Proof. The assumption (55) implies that the factor of satisfies(57)Therefore this term is controlled by the same term in the lefthand side of (47). With the above regularity of the data, an application of Gronwall’s Lemma yields a uniform bound for and . In turn, this bound permits to estimate and , uniformly with respect to , , and Δt. Finally, the regularity of ρgη implies a similar bound for .
Remark 3.1. The exponential factor in the righthand side of (56) can be avoided by using L^{∞} − L^{1} bounds in time such asThis is sound, as long as it is done with data, but the same strategy applied toleads to a relation between and Δt that is less favorable than that given by (55).
4 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 δ the finite difference in time for any function f,(58)
4.1 The algorithmic error
Let us start with an arbitrary value of . Recall that the contraction property yields (41),whereWe have seen in Proposition 3.1 that(59)with the constants C and of (42). The next proposition sharpens the bound for the pressure in Proposition 3.2. Below, we use the notation for any function q,(60)Π_{ 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,(61)
Proposition 4.1.
Assuming that the solution is sufficiently smooth in time, we have for all n, 1 ≤ n ≤ N, (62) where β is defined by (39).
Proof. By proceeding as in Proposition 3.2, subtracting , and , and recalling that ρ _{ f,r } gη is independent of time, we first writeBy integrating the flow equation (21) over]t _{ n − 1}, t _{ n }[, the end of the last line above can be written asThen (62) follows by suitable applications of Young’s inequality.
Thus, by combining (59), (62), and (41), we derive an upper bound for ,(63)
4.2 Error equation
In addition to Π_{ 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 }, and divide both sides by Δt; with the notation (58), (60), and (61), this givesBy subtracting (35) at level from this equation (see (46)), we obtain for all θ _{ h } ∈ M _{ h },The final error equation for flow is derived by adding and subtracting approximations of p and u , and multiplying all terms by Δt. This yields(64)
Now, we turn to the displacement and write (20) tested with v _{ h } at time t _{ n }. We obtainBy observing that , the difference between this equation and (36) again with readsIt remains to add and subtract approximations of and u , to derive the final error displacement equation for all ,(65)Let us set(66)The final error equality is obtained by testing (64) with and (65) by , for 1 ≤ n ≤ N, and multiplying both sides by 2(67)
4.3 Error inequality
To simplify, set(68)By comparing (67) at level m with (50), we see that their lefthand sides have much the same structure. Therefore, to deduce an error inequality from (67), it suffices to suitably group the terms in its righthand side so as to associate them with corresponding terms in the righthand side of (50). First, the termcan be treated asNext, we can writeThus, this term can be treated asHence, proceeding as in the proof of Proposition 3.3, assuming again that , and summing over m, from 1 to n, the part of the righthand side of (67) corresponding to these two terms is bounded by(69)The corresponding first group of terms in the lefthand side is(70)Next, we consider the terms involving and that must be summed by parts. For the first term, we writeYoung’s inequality givesThe contribution of this term to the righthand side is(71)and the corresponding group of terms that remain in the lefthand side is(72)To simplify the treatment of the terms involving , it is convenient to extend by zero outside Ω_{1}. With this convention, we writeThen we deduce by the above argumentThe contribution of this term to the righthand side is(73)and the corresponding group of terms that remain in the lefthand side is(74)There remain the two terms involving the gradient of ; with the notation (68), they can be written asBy Young’s inequality, when summed over m, this term has the bound for any γ > 0,Thus the contribution of this term to the righthand side is(75)and the corresponding term that remain in the lefthand side is(76)
By collecting the terms L _{ i } in the lefthand side and R _{ i } in the righthand side, for 1 ≤ i ≤ 4, we deduce the following error inequality (compare with Proposition 3.3):
Proposition 4.2.
Let q belong to , to , and to . Assume that and that is extended by zero outside . With the notation (60), (66), and (68) , the following bound holds for all n, 1 ≤ n ≤ N, and any γ ∈ ]0, 2[, (77) where and are respectively contributions of the initial errors and interpolation errors in space and time, (78) (79)
4.4 Error estimate
With the notation (66) and (68), (63) reads(80)
The factor of the term involving isand this factor must be controlled in the lefthand side of (77) bywith a parameter γ ∈ ]0, 2[to be chosen. If we choose for instance , a sufficient condition for this control iswhich 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 thatthis assumption implieswhich is not sufficient to guarantee a first order in time that requires a factor (Δt)^{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,(81)
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 ∂ _{ t } f ∈ L^{2}(Ω × ]0, T[)^{ d } , ∂ _{ t } t _{ N } ∈ L^{ 2 }(Γ _{N} × ]0, T[)^{ d } , and ∂ _{ t } q ∈ L^{ 2 }(Ω_{1} × ]0, T[), that the solution satisfies ∂ _{ t } p ∈ L^{2}(0, T; H^{1}(Ω_{1})) and ∂ _{ t } u ∈ L^{2}(0, T; H^{1}(Ω)^{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, (82) where (83) (84) (85)
Proof. By applying the second part of (81) to the first, second, and last terms of the righthand side of (80) and the first part to the remaining terms, we derive(86)Now, we substitute (86) into (77) with the choice . The first sum in the second line of the above righthand side cancels with the corresponding sum in the lefthand side of (77) with the exception of the first and last term. More precisely, there remains a term (that is neglected for the moment) in the lefthand sideand a term in the righthand 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),(87)Hence, grouping terms and applying (87) for the time errors we deduce the following intermediate result:(88)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, Δt and h,(89)
There remains to take into account the interpolation error in space and the initial errors. Suppose that ∂_{ t } p belongs to and ∂_{ t } u belongs to . Recall that k ≥ 1, respectively m ≥ 1, is the degree of the finite element polynomials of the displacement, respectively the pressure. By taking into account the approximation properties (30) of Π_{ h } and (29) of R _{ h } and the regularity (28) of the triangulation, we deduce, with various constants C independent of n, Δt, and h,(90)
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 , 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 ∈ L^{2}(Ω × ]0, T[)^{ d } , ∂ _{ t } t _{ N } ∈ L^{2}(Γ_{ N } × ]0, T[)^{ d } , and ∂ _{ t } q ∈ L^{ 2 }(Ω_{1} × ]0, T[), that the solution satisfies and , 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, (91) with a constant C independent of h and Δt.
5 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}(Ω_{1}), with interpolation operator ρ _{ h }, ( Z _{ h }, Q _{ h }) being compatible in the sense that they satisfy a uniform discrete infsup condition. In particular, we make the assumption, commonly used in mixed methods,(92)
As previously, we denote . Starting from(93)the fixedstress splitting algorithm computes by solving (33), by solving(94)and by (34),
Then for ,

Set , , , and .

For , compute
(a) the pair in by solving(95) (96)
(b) compute the predictor of the difference in fluid content by
(c) by (37)
(d) compute the corrector to difference in fluid content by
If threshold, increment and return to (a);
if threshold, set(97)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.
5.1 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 β andwe derive by the same argument
Therefore(98)and again we must estimate . 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 with , splittingand using the velocitypressure equation to eliminate the divergence in both terms of the above righthand side. This givesThus, we have the analogue of Proposition 3.2, with the same notation.
Proposition 5.1.
We have for all n, 1 ≤ n ≤ N, (99)
Then a combination of (42) and (99) yields the analogue of (44),(100)which again relies on the stability of the scheme.
Now, by using again the velocitypressure equation, so thatit 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,(101)
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 Δt, such that (102)
5.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 Π_{ 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 introduce the term that we write as
For the first term, since , we apply the discrete velocitypressure 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, (103)
Then (63) is replaced by(104)
Note that the factor of in the righthand side of (104) is smaller than its counter part in (36) by a factor of . In contrast, the factors of the other terms in the righthand 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 velocitypressure relation in which we make use of the assumption (92). Thus, we writeand hence the difference between the exact and discrete velocitypressure equations reads(105)
This is an important simplification because it will eliminate the factor from the error equation. In addition to the notation and of (66), we set(106)
Then (67) is replaced by(107)
A comparison with (67) immediately shows that the conclusion of Proposition 4.2 holds with the same initial error and slightly different interpolation error , 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(108)
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 factor mentioned above, it becomes(109)
6 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 payzone (or aquifer or flow domain) is given as (Verruijt [46])(110)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(111)and c is the diffusivity coefficient given bywhere is the drained bulk modulus and 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 m. An underburden with a thickness of m is also provided. No flow boundary conditions are imposed on the payzone. Gravity is ignored. Denoting the top surface of the mechanics domain as , the boundary conditions for the mechanics subproblem are
Fig. 2
Schematic for single well in an infinite confined aquifer problem. 
Parameters for single well in an infinite confined acquifer problem.
The Biot parameter is nonzero only inside the payzone and zero outside the payzone and formally given as
Since for this computational example the term represents a well and is not in , we compute estimates in where . Here is a cylinder of radius m (which is times the areal size of the domain) and height of the payzone. The error norm for the pore pressure in the payzone is computed using the midpoint quadrature rule:(112)where m _{ e } is the center of mass of element E, τ is time, p(τ, m _{ e }) is the analytical solution for pore pressure at time τ and position m _{ e }, p _{ h }(τ, 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 inhouse 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.
Order of convergence of pore pressure solution with time step of 0.01 day.
Order of convergence of pore pressure solution with time step of 0.005 day.
References
 von Terzaghi K. (1944) Theoretical soil mechanics, Chapman, London. [Google Scholar]
 Biot M.A. (1941) Consolidation settlement under a rectangular load distribution, J. Appl. Phys. 12, 5, 426–430. [Google Scholar]
 Biot M.A. (1941) General theory of threedimensional consolidation, J. Appl. Phys. 12, 2, 155–164. [Google Scholar]
 Settari A., Mourits F.M. (1994) Coupling of geomechanics and reservoir simulation models, in: H.J. Siriwardane, M.M. Zema (eds), Computer methods and advances in geomechanics, Balkema, Rotterdam, pp. 2151–2158. [Google Scholar]
 Settari A., Mourits F.M. (1998) A coupled reservoir and geomechanical simulation system, SPE J. 3, 219–226. [Google Scholar]
 Gai X., Dean R.H., Wheeler M.F., Liu R. (2003) Coupled geomechanical and reservoir modeling on parallel computers, The SPE Reservoir Simulation Symposium, Houston, Texas, February 3–5. [Google Scholar]
 Gai X., Sun S., Wheeler M.F., Klie H. (2005) A timestepping scheme for coupled reservoir flow and geomechanics on nonmatching grids, SPE Annual Technical Conference and Exhibition, SPE97054. [Google Scholar]
 Gai X. (2004) A coupled geomechanics and reservoir flow model on parallel computers, PhD Thesis, The University of Texas at Austin, Austin, Texas. [Google Scholar]
 Dean R.H., Gai X., Stone C.M., Minkoff S.E. (2006) A comparison of techniques for coupling porous flow and geomechanics, SPE J. 11, 1, 132–140. [CrossRef] [Google Scholar]
 Girault V., Kumar K., Wheeler M.F. (2016) Convergence of iterative coupling of geomechanics with flow in a fractured poroelastic medium, Comput. Geosci. 20, 997–1011. [Google Scholar]
 Almani T., Kumar K., Dogru A., Singh G., Wheeler M.F. (2016) Convergence analysis of multirate fixedstress split iterative schemes for coupling flow with geomechanics, Comput. Methods Appl. Mech. Eng. 311, 180–207. [Google Scholar]
 Almani T., Kumar K., Singh G., Wheeler M.F. (2016) Stability of multirate explicit coupled of geomechanics with flow in a poroelastic medium, ICES Report 16–12, Institute for Computational Engineering and Sciences, The University of Texas at Austin, Austin, Texas. [Google Scholar]
 Almani T., Kumar K., Wheeler M.F. (2017) Convergence and error analysis of fully discrete iterative coupling schemes for coupling flow with geomechanics, Comput. Geosci. 21, 5, 1157–1172. [Google Scholar]
 Almani T., Dogru A.H., Kumar K., Singh G., Wheeler M.F. (2016) Convergence of multirate iterative coupling of geomechanics with flow in a poroelastic medium, Saudi Aramco J. Technol. [Google Scholar]
 Mikelić A., Wheeler M.F. (2013) Convergence of iterative coupling for coupled flow and geomechanics, Comput. Geosci. 17, 455–461. [Google Scholar]
 Mikelić A., Wang B., Wheeler M.F. (2014) Numerical convergence study of iterative coupling for coupled flow and geomechanics, Comput. Geosci. 18, 325–341. [Google Scholar]
 Kim J., Tchelepi H.A., Juanes R. (2009) Stability, accuracy and efficiency of sequential methods for coupled flow and geomechanics, The SPE Reservoir Simulation Symposium, Houston, Texas, SPE119084. [Google Scholar]
 Kim J., Tchelepi H.A., Juanes R. (2011) Stability and convergence of sequential methods for coupled flow and geomechanics: fixedstress and fixedstrain splits, Comput. Methods Appl. Mech. Eng. 200, 13–16, 1591–1606. [Google Scholar]
 Gai X. (2004) A coupled geomechanics and reservoir flow model on parallel computers, PhD Thesis, The University of Texas at Austin, Austin, Texas. [Google Scholar]
 Castelletto N., White J.A., Tchelepi H.A. (2015) Accuracy and convergence properties of the fixedstress iterative solution of twoway coupled poromechanics, Int. J. Numer. Anal. Methods Geomech. 39, 1593–1618. [Google Scholar]
 Castelletto N., White J.A., Tchelepi H.A. (2014) A unified framework for fully implicit and sequentialimplicit schemes for coupled poroelasticity, ECMOR XIV – 14th European Conference on the Mathematics of Oil Recovery, September 8–11. [Google Scholar]
 Kumar K., Almani T., Singh G., Wheeler M.F. (2016) Multirate undrained splitting for coupled flow and geomechanics in porous media, Springer International Publishing, Cham, pp. 431–440. [Google Scholar]
 White D., Ganis B., Liu R., Wheeler M.F. (2017) A near wellbore study with a DruckerPrager plasticity model coupled with a parallel compositional reservoir simulator, SPE Reservoir Simulation Symposium. [Google Scholar]
 Borregales M., Radu F.A., Kumar K., Nordbotten J.M. (2018) Robust iterative schemes for nonlinear poromechanics, Comput. Geosci. 22, 4, 1021–1038. [Google Scholar]
 da Silva R.O.S., Murad M.A.M., Obregon J.A.L.O. (2018) A new fixedstress split scheme in poroplastic media incorporating general plastic porosity constitutive theories, ECMOR XVI – 16th European Conference on the Mathematics of Oil Recovery, September 3–6. [Google Scholar]
 Dana S., Ganis B., Wheeler M.F. (2018) A multiscale fixed stress split iterative scheme for coupled flow and poromechanics in deep subsurface reservoirs, J. Comput. Phys. 352, 1–22. [Google Scholar]
 Bause M., Radu F.A., Köcher U. (2017) Spacetime finite element approximation of the Biot poroelasticity system with iterative coupling, Comput. Methods Appl. Mech. Eng. 320, 745–768. [Google Scholar]
 Borregales M., Kumar K., Radu F.A., Rodrigo C., Gaspar F.J. (2018) A parallelintime fixedstress splitting method for Biot’s consolidation model, arXiv preprint arXiv:1802.00949. [Google Scholar]
 Rodrigo C., Gaspar F.J., Hu X., Zikatanov L.T. (2016) Stability and monotonicity for some discretizations of the Biots consolidation model, Comput. Methods Appl. Mech. Eng. 298, 183–204. [Google Scholar]
 Odsaeter L.H., Wheeler M.F., Kvamsdal T., Larson M.G. (2017) Postprocessing of nonconservative flux for compatibility with transport in heterogeneous media, Comput. Methods Appl. Mech. Eng. 315, 1, 799–830. [Google Scholar]
 Lee S., Lee Y.J., Wheeler M.F. (2016) A locally conservative enriched Galerkin approximation and efficient solver for elliptic and parabolic problems, SIAM J Sci. Comput. 38, 3, A1404–A1429. [Google Scholar]
 Iding M., Ringrose P. (2010) Evaluating the impact of fractures on the performance of the In Salah CO2 storage site, Int. J. Greenhouse Gas Control 4, 2, 242–248. The Ninth International Conference on Greenhouse Gas Control Technologies. [CrossRef] [Google Scholar]
 Mikelić A., Wheeler M.F. (2012) On the interface law between a deformable porous medium containing a viscous fluid and an elastic body, Math. Models Methods Appl. Sci. 22, 11, 1250031. [Google Scholar]
 Adams R.A. (1975) Sobolev spaces, Pure and Applied Mathematics, Vol. 65, Academic Press, New York/London. [Google Scholar]
 Nečas J. (1967) Les méthodes directes en théorie des équations elliptiques, Masson, Paris. [Google Scholar]
 Lions J.L., Magenes E. (1972) Nonhomogeneous boundary value problems and applications, Vol. I, SpringerVerlag, New York. [Google Scholar]
 Grisvard P. (1985) Elliptic problems in nonsmooth domains, volume 24 of Monographs and Studies in Mathematics, Pitman, Boston, MA. [Google Scholar]
 Girault V., Raviart P.A. (1986) Finite element methods for NavierStokes equations: Theory and algorithms, volume 5 of Springer Series in Computational Mathematics, SpringerVerlag, Berlin. [CrossRef] [Google Scholar]
 Girault V., Wheeler M.F., Ganis B., Mear M.E. (2015) A lubrication fracture model in a poroelastic medium, Math. Models Methods Appl. Sci. 25, 4, 587–645. [Google Scholar]
 Coussy O. (2003) Poromechanics, WileyBlackwell. [CrossRef] [Google Scholar]
 Girault V., Pencheva G.V., Wheeler M.F., Wildey T.M. (2011) Domain decomposition for poroelasticity and elasticity with DG jumps and mortars, Math. Models Methods Appl. Sci. 21, 1, 169–213. [Google Scholar]
 Girault V., Pencheva G.V., Wheeler M.F., Wildey T.M. (2009) Domain decomposition for linear elasticity with DG jumps and mortars, Comput. Methods Appl. Mech. Eng. 198, 21–26, 1751–1765. [Google Scholar]
 Ciarlet P.G. (1991) Basic error estimates for elliptic problems, in: Ciarlet P.G., Lions J.L. (eds), Handbook of numerical analysis, Vol. II, Elsevier Sciences, NorthHolland, Amsterdam, pp.17–351. [Google Scholar]
 Scott L.R., Zhang S. (1990) Finite element interpolation of nonsmooth functions satisfying boundary conditions, Math. Comp. 54, 483–493. [CrossRef] [MathSciNet] [Google Scholar]
 Mikelić A., Wheeler M.F. (2013) Convergence of iterative coupling for coupled flow and mechanics, Comput. Geosci. 17, 3, 455–461. [Google Scholar]
 Verruijt A. (2013) Theory and problems of poroelasticity, Delft University of Technology. [Google Scholar]
 Ingram R., Wheeler M.F., Yotov I. (2010) A multipoint flux mixed finite element method on hexahedra, SIAM J. Numer. Anal. 48, 4, 1281–1312. [Google Scholar]
 Falgout R.D., Yang U.M. (2002) hypre: A library of high performance preconditioners, in: Sloot P.M.A., Hoekstra A.G., Tan C.J.K., Dongarra J.J. (eds), Computational Science – ICCS 2002, Springer, Berlin, Heidelberg, pp. 632–641. [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1
Payzone with surrounding rock. 

In the text 
Fig. 2
Schematic for single well in an infinite confined aquifer problem. 

In the text 