© V. Girault et al., published by IFP Energies nouvelles, 2019

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 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–14]. They are known as the fixed-stress split, fixed-strain split, undrained-split, and drained split schemes (Mikelić and Wheeler [15], Mikelić 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 (Mikelić and Wheeler [15], Mikelić 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 poro-elastic 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₂ 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₂ 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₂ storage site in the world. At Krechba, CO₂ 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₂ 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 Mikelić & 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 non-pay 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.

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²(Ω) is denoted by and if the domain of integration is clear from the context, we suppress the index . Let (k ₁, k ₂, k ₃) denote a triple of non-negative integers, set |k| = k ₁ + k ₂ + k ₃ and define the partial derivative ∂ ^k byThen, for any non-negative 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 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 Ω 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 ₁(Ω, Γ) depending only on Ω and Γ. Here ε ( v ) denotes the strain tensor. When combined with Poincaré’s inequality, also valid for all functions v in H¹(Ω) that vanish on Γ, with another constant C ₂(Ω, Γ) depending only on Ω and Γ,(4)we recover the full H¹ norm in the left-hand 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¹(Ω) ^d that vanish on Γ(6)

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 ||·|| _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 Ω₁ 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 poro-elastic model holds in Ω₁ (the pay-zone) while an elastic model holds in Ω₂ (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 Ω₁, the quasi-static Biot’s system from consolidation theory.

2. On the boundary of Ω₁, 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 Ω₂, the elasticity equations.

Let Γ₁₂ denote the boundary of Ω₁ and let n ₁₂ be the unit normal on Γ₁₂ exterior to Ω₁. 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 Lipschitz-continuous. 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 three-dimensional case. Let f be the body force in Ω. In the non-pay zone Ω₂, the governing equations are those of linear elasticity. Recalling the Cauchy stress tensor for linear elasticity,(7)we have a.e. in Ω₂ × ]0, T[(8)

In the pay-zone Ω₁, 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 ρ _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⁻⁸ or 10⁻⁹), 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 Ω₁ × ]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 Biot-Willis constant, which is usually around one, μ _f is the fluid’s viscosity assumed to be constant, φ ₀ 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 Ω₁ is reduced to Γ₁₂, that is entirely contained in Ω, the only boundary conditions on Γ₁₂ 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 Γ₁₂ 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¹(0, T; L²(Ω) ^d ), q belongs to L²(Ω₁ × ]0, T[), t _N is in H¹(0, T; L²(Γ _N ) ^d ), and p ₀ belongs to H¹(Ω). 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 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 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.

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²(Ω₁)) ∩ L²(0, T; H¹(Ω₁)) 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 quasi-local 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.

1. Set , , and .

2. For , compute

(a) by solving(35)

(b) compute the predictor of the difference in fluid content by(36)

(c) by(37)

(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 fixed-stress 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 ₁ (Ω, Γ)C ₂ (Ω, Γ) 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 left-hand 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 left-hand 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 left-hand 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 right-hand side of (56) can be avoided by using L^∞ − L¹ 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¹(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 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 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 right-hand side of (67) corresponding to these two terms is bounded by(69)The corresponding first group of terms in the left-hand 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 right-hand side is(71)and the corresponding group of terms that remain in the left-hand side is(72)To simplify the treatment of the terms involving , it is convenient to extend by zero outside Ω₁. With this convention, we writeThen we deduce by the above argumentThe contribution of this term to the right-hand side is(73)and the corresponding group of terms that remain in the left-hand 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 right-hand side is(75)and the corresponding term that remain in the left-hand side is(76)

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 , 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 left-hand 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 right-hand 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)². 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²(Ω × ]0, T[) ^d , ∂ _t t _N  ∈ L ² (Γ _N × ]0, T[) ^d , and ∂ _t q ∈ L ² (Ω₁ × ]0, T[), that the solution satisfies ∂ _t p ∈ L²(0, T; H¹(Ω₁)) and ∂ _t u  ∈ L²(0, T; H¹(Ω)^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 right-hand 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 right-hand side cancels with the corresponding sum in the left-hand 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 left-hand sideand 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¹(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²(Ω × ]0, T[) ^d , ∂ _t t _N  ∈ L²(Γ _N  × ]0, T[) ^d , and ∂ _t q ∈ L ² (Ω₁ × ]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²(Ω₁), with interpolation operator ρ _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,(92)

As previously, we denote . Starting from(93)the fixed-stress splitting algorithm computes by solving (33), by solving(94)and by (34),

Then for ,

1. Set , , , and .

2. 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 velocity-pressure equation to eliminate the divergence in both terms of the above right-hand 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 velocity-pressure 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 velocity-pressure relation, whereas the other two terms are left unchanged in order to be combined with q ⁿ . 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 right-hand 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 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 writeand hence the difference between the exact and discrete velocity-pressure 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 pay-zone (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 ₁(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 pay-zone. Gravity is ignored. Denoting the top surface of the mechanics domain as , the boundary conditions for the mechanics subproblem are

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 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 pay-zone. The error norm for the pore pressure in the pay-zone 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 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⁻⁸. The code takes 1–3 coupling iterations to convergence at every time step.

References

All Tables

All Figures