Accurate discretization of poroelasticity without Darcy stability

In this manuscript we focus on the question: what is the correct notion of Stokes–Biot stability? Stokes–Biot stable discretizations have been introduced, independently by several authors, as a means of discretizing Biot’s equations of poroelasticity; such schemes retain their stability and convergence properties, with respect to appropriately defined norms, in the context of a vanishing storage coefficient and a vanishing hydraulic conductivity. The basic premise of a Stokes–Biot stable discretization is: one part Stokes stability and one part mixed Darcy stability. In this manuscript we remark on the observation that the latter condition can be generalized to a wider class of discrete spaces. In particular: a parameter-uniform inf-sup condition for a mixed Darcy sub-problem is not strictly necessary to retain the practical advantages currently enjoyed by the class of Stokes–Biot stable Euler–Galerkin discretization schemes.


Introduction
In this note, we consider a three-field formulation of the time-dependent Biot equations describing flow through an isotropic, porous and linearly elastic medium, reading as: find the elastic displacement u, the Darcy flux z and the (negative) fluid pressure p such that for a given body force f , source s, and given g (typically g = 0) over a domain Ω ⊂ R d (d = 1, 2, 3).The isotropic elastic stress tensor σ(u) = µε(u) + λtrε(u) where tr is the matrix trace.The material parameters are the elastic Lamé parameters µ and λ, the Biot-Willis coefficient α, the storage coefficient c 0 ≥ 0 and the hydraulic conductivity κ = K/µ f > 0, in which K is the material permeability, and µ f is the fluid viscosity.Moreover, ε denotes the (row-wise) symmetric gradient, div is the divergence, ∇ is the gradient, and ∂ t denotes the (continuous) time-derivative.
The three field formulation (1a)-(1c) combines one scalar, time-dependent partial differential equation and two, stationary, vector partial differential equations.This combination of time-dependent and time-independent equations can lead to non-trivial issues when considering discretizations of the time derivative; as a result: several splitting scheme approaches have been proposed [7,11,13,15,24].In this manuscript we will focus on a monolithic approach, namely a straightforward backward Euler scheme, where all unknowns are solved for simultaneously.In the case of monolithic time discretization schemes: robustness with respect to material parameters in spatial discretizations of (1) is a central concern and has been the topic of several recent investigations; c.f. e.g.[13,14,22].A notable difficulty, both practically and theoretically, is that the parameter λ may be very large, while κ and c 0 may be very small.The former corresponds to the (nearly) incompressible regime, while the latter two correspond to the (nearly) impermeable and low-storage regimes, respectively.Special care is required in the formulation and analysis of discretizations of (1) to retain stability and convergence within these parameter ranges.In this context, the novel and important concept of Stokes-Biot stability has emerged [13,16,22] as a guide for the design of wellposed Euler-Galerkin discretizations of (1) for nearly impermeable or low-storage materials; i.e., as κ, c 0 → 0. In this note, we aim to shed further light on this concept: we show that the original Stokes-Biot stability condition can be relaxed and we formalize a new concept of minimal Stokes-Biot stability.
To this end and more formally, we begin by considering a three-field variational formulation of the related system of (time-independent) equations: find u ∈ U , z ∈ W , and p ∈ Q such that τ κ −1 (z, w) + τ (div w, p) = τ (g, w) ∀ w ∈ W, (2b) (div u, q) + (τ div z, q) − (c 0 p, q) = (τ s + div ū − c 0 p, q) , ∀ q ∈ Q, (2c) for given f, g, s, ū, p and with (•, •) denoting the standard L 2 -inner product over the domain Ω.The system (2) is a representation of the equations to be solved at each discrete time t after e.g. an implicit Euler time discretization of (1) with time step τ > 0 and a prescribed set of homogeneous boundary conditions.We let τ = 1 for now without loss of generality, and refer to (2) as (a mixed variational formulation of) the steady Biot equations.
A-priori error estimates, in appropriate parameter-dependent norms, have been established for both non-conforming [13,16] and conforming [22] discretizations of (1) or (2) satisfying these conditions.To exemplify, consider the finite element pairings P d 2 × RT 0 × DG 0 (product space of continuous piecewise quadratic vector fields, lowest order Raviart-Thomas elements and piecewise constants) and P d 1 × RT 0 × DG 0 .The former pairing satisfies conditions (i) and (ii) above (for given κ > 0), and is observed to converge even for κ ≪ 1, see e.g.[22] or Table 1a below.The latter pairing, which violates condition (i), can easily fail to converge when κ is sufficiently small.However, the question remains: are the Stokes-Biot stability conditions (i)-(ii), and in particular condition (ii), necessary to guarantee convergence?This question is key since the Darcy stability condition can easily fail to hold uniformly in κ ≪ 1, as also noted in e.g.[13,Rmk. 5].More precisely, the continuous mixed Darcy problem (4) does not satisfy the Babuska-Brezzi conditions [6] with bounds independent of 0 < κ ≪ 1 in the standard H(div) × L 2 norm.To compensate, permeability-weighted flux and pressure norms, such as e.g.κ −1/2 H(div) × κ 1/2 L 2 , have been suggested as viable alternatives [13].However, the resort to a permeability-weighted pressure space is not entirely satisfactorily as the relation between ( 2) and the Stokes equations (3) in the limit κ → 0 points at p ∈ L 2 rather than p ∈ κ 1/2 L 2 .Moreover, numerical experiments demonstrate convergence of the pressure in the L 2 -norm even for diminishing κ, see e.g.Table 1a for the pairing P d 2 ×RT 0 ×DG 0 .Conversely, consider the pairing P d 2 × P d 1 × DG 0 which violates the Darcy condition (ii) for any κ > 0 and thus does not satisfy the Stokes-Biot stability conditions.However, in numerical experiments, see Table 1b, the pairing appears stable with the displacement and pressure comparable rates as for P d 2 × RT 0 × DG 0 , even for c 0 = 0 and small κ.These observations thus call into question the precise role of the Darcy stability assumption in conforming mixed finite element discretizations of (1) or (2).
In this manuscript we advance the claim that the Darcy stability assumption, of the Stokes-Biot stability conditions, is not necessary; at least not in the case of Euler-Galerkin discretizations of (1) or (2).Specifically, we discuss that it is not needed to prove the stability of such schemes, when 0 ≤ c 0 ≤ 1 or 0 < κ ≪ 1, nor is it needed to establish the arguments of an a-priori error analysis.Instead, what we will observe is that the following two assumptions are key: (I) the displacement-pressure pairing U h ×Q h is a stable pair for the incompressible Stokes equations (3); and that (II) the inclusion div W h ⊆ Q h holds.We return to, and formalize, these minimal Stokes-Biot stability conditions in Section 3.
Moreover, we establish suitable norms with respect to which a uniform-in-κ Darcy stability condition for the pairing (W h , Q h ) need not, as initially considered [13,22], be assumed.An Euler-Galerkin scheme is well posed in these norms and we show that the corresponding a priori error convergence rates hold in the limit as κ → 0 and coincide with canonically expected rates for well known mixed three-field finite element paradigms; e.g.first order for discretizations using linear or Raviart-Thomas type flux approximations, etc. Remark 1.For elasticity systems in primal form, the incompressible case λ ≫ µ is efficiently dealt with by introducing an auxiliary variable, the 'solid pressure' p s = λ∇ • u; this is known as Herrmann's method [12].In the case of poroelasticity, the solid pressure can be replaced by the total pressure, λ∇ • u − p, to achieve robustness with respect to λ.As this topic is well covered in [15,16,18,17,20] we will not address this issue in this note but remark  and flux z − z h div / z div (bottom three rows) for varying κ on a series of uniform meshes T h with mesh size h.The last column 'Rate' denotes the order of convergence using for the last two values in each row.The exact solutions ũ, p, z, defined in Section 6, were represented by continuous piecewise cubic interpolants in the error computations.Similar results were obtained for κ = 10 −2 , 10 −6 , 10 −10 (data not shown).c 0 = 0. (A): that robustness would require the introduction of an auxiliary variable also for the three field formulation (1), see e.g.[15].
The remainder of this manuscript is organized as follows: Section 2 describes basic spaces and notation that will be used throughout; Section 3 overviews the current view of Stokes-Biot stability [13,16,22]; Section 4 introduces a slight relaxation on the Stokes-Biot stable conditions and recalls a well-posedness argument for Euler-Galerkin discrete schemes; Section 5 is aimed at a priori estimates for discretizations satisfying the relaxed conditions; finally, Section 6 is a numerical example demonstrating the retention of convergence behaviour as κ → 0.  [19,23,25].We will consider discretizations of Ω by simplicial complexes of order d.All triangulations, T h of Ω, will be assumed to be shape regular with the maximal element diameter, also referred to as the mesh resolution or mesh size, of T h denoted by h.
We let L 2 (Ω; R d ), H(div, Ω) and H 1 (Ω; R d ) denote the standard Sobolev spaces of squareintegrable fields over Ω, fields with square-integrable divergence, and fields with squareintegrable gradient, respectively, and define the associated standard norms denoting the standard L 2 (Ω)-inner product.We will frequently drop the arguments Ω and R d from the notation when the meaning is clear from the context.The notation H 1 Γ (Ω) represents those functions in H 1 with zero trace on Γ ⊆ ∂Ω.Similarly, H Γ (div, Ω) denotes fields in H(div, Ω) with zero (normal) trace on Γ ⊆ ∂Ω in the appropriate sense [4].We also define the standard space of square-integrable functions with zero average: We will also use parameter-weighted norms.For a Banach space X and real parameter α > 0, the space αX signifies X equipped with the α-weighted norm f αX = α f X .Finally, for a coercive and continuous bilinear form a : V × V → R, we will also write

2.2.
Intersections and sums of Hilbert spaces.Let X ⊂ Z and Y ⊂ Z be two Hilbert spaces with a common ambient Hilbert space Z.The intersection space, denoted X ∩ Y , is a Hilbert space with norm For instance, to illustrate our notation, the norm on the intersection space κ −1/2 L 2 ∩ H(div) is given by and is also a Hilbert space.See e.g.[3, Ch. 2] for a further discussion of sum and intersection spaces.
2.3.Operators.For a given time step size τ , times t m−1 and t m and fields u m ≈ u(t m ) and u m−1 ≈ u(t m−1 ), we will make use of a discrete derivative notation (5) 2.4.Finite element spaces.Now, suppose that Ω ⊂ R d is a polygonal and let C k (Ω) denote the space of k-continuously differentiable functions defined on Ω.Let D ⊆ Ω and let P k (D) ⊂ C ∞ (D) denote the set of polynomials of total degree k defined on D. Let T h be a simplicial triangulation of Ω and let T ∈ T h be any simplex; we denote the restriction of a function f to T ∈ T h by f T .The notation for the Lagrange elements of order k used here is then The discontinuous Galerkin spaces of order k relax the overall continuity requirement of the Lagrange finite element spaces; they are defined by A comprehensive discussion on Lagrange and discontinuous Galerkin elements and their interpolation properties can be found in e.g.[8] and [21] respectively.We will also make use of the Brezzi-Douglas-Marini and Raviart-Thomas finite element spaces [4,Sec. 2.3].Throughout the rest of the manuscript we use the notation P k , DG k , BDM k and RT k in reference to the spaces defined above; that is, we drop the additional mesh domain specification.We will also abuse notation, and reuse the notation P k for displacement or flux vector fields with P k components whenever the context is clear.
We follow [23] and remark: that under light regularity assumptions on the source data, (i.e.f , g, and h), and upon assuming the initial fluid content, η(x, 0), is in L 2 (Ω) then there exists a unique solution to (1) satisfying the boundary conditions above [23,25].
2.6.Material parameters.In applied settings the material parameters of Biot's equations, (1a)-(1c) often vary in space and can even be discontinuous.In the analysis here we will assume that the parameters µ, λ, α, κ, and c 0 are spatially constant.For simplicity and without loss of generality for α ≫ 0, we set α = 1.This view can either be interpreted literally or as having divided (1a)-(1c) through by α to obtain rescaled material parameters.
Following the work of previous authors [13,14,16,22] we will be interested in the case where 0 < µ and 0 < λ are fixed, 0 ≤ c 0 ≤ 1 and 0 < κ ≤ 1 are spatially constant but otherwise arbitrarily selected.This assumption allows salient analytic points, which are the focus of the present work, to be readily made; in particular [14] has defined specific parameter-dependent norms suitable for the analysis of conforming Euler-Galerkin discretizations of (1a)-(1c).The arguments we advance can be extended to the case where the parameters µ, λ are uniformly bounded functions in L ∞ ((0, T ]; C 2 (Ω)), c 0 is a non-negative H 1 function with image in [0, 1] and κ is a symmetric, positive-definite matrix with entries in H 1 .

The Stokes-Biot stability conditions for conforming Euler-Galerkin schemes
Combining the nature of ( 1) with the boundary conditions ( 8), we define the spaces If Dirichlet conditions are imposed for the displacement on the entire boundary and thus the pressure is only determined up to a constant (i.e. if Γ c = ∂Ω) we instead let Q = L 2 0 .We consider the following variational formulation of (1) over the time interval (0, T ]: for a.e.t ∈ (0, T ], find the displacement u, flux z and pressure p such that u(t) ∈ U , z(t) ∈ Z and The bilinear forms in (10) are given by: As noted in [22]: the existence and uniqueness of a solution (u, z, p) to (10), with continuous dependence on f , g and s, has been established by previous authors [19,23,25].
3.1.An Euler-Galerkin discrete scheme.Following [22] we consider Euler-Galerkin discretizations; i.e., conforming finite element spaces in space and an implicit Euler in time, of (10).Let 0 = t 0 < t 1 < . . .< t N = T be a uniform partition of the time interval [0, T ].The constant time step is then τ = τ m = t m − t m−1 .For the function f (t, x), evaluation at t m is denoted by f m = f (t m , x), and similarly for g and s.We define conforming discrete spaces The Euler-Galerkin discrete scheme of Biot's equations then reads as follows: for all v ∈ U h , w ∈ W h and q ∈ Q h , and where we have made use of the discrete derivative notation (5).

3.2.
The Stokes-Biot stability conditions.The Stokes-Biot stability conditions were introduced independently, in slightly different contexts, by several authors [13,16,22] and guide the selection of discrete spaces, U h × W h × Q h , for (13).We recall a succinct statement of the (conforming) Stokes-Biot stability conditions, used in analogous forms by all original authors [13,16,22], here for posterity: Definition 1 (c.f.[22,Defn. 3.1]).The discrete spaces U h ⊂ U , W h ⊂ W and Q h ⊂ Q are called a Stokes-Biot stable discretization if and only if the following conditions are satisfied: (i) The bilinear form a, as defined by (11), is bounded and coercive on We remark that [13,16] were not conforming.More precisely, the Stokes and Darcy stability assumptions of Definition 1 entail that the relevant discrete spaces are stable in the (discrete) Babuška-Brezzi sense [4,6] for the discrete Stokes and Darcy problems, respectively.We will now examine the Darcy stability condition more closely.

3.3.
The Darcy stability condition.The discrete Darcy problem reads as: find with b and c as defined by (11).It is also assumed that b and c are continuous over W × Q and W × W with respect to the relevant norms; i.e. there exist constants C b > 0 and C c > 0, independent of h, such that The assumption of discrete Darcy stability, and thus the existence of solutions to the discrete Darcy problem, has been used to define Galerkin projectors for use in the a-priori analysis of the Biot equations (13) (c.f. for instance [22,Sec. 4.2]).Given z(t) ∈ W and p(t) ∈ Q solving the continuous Biot equations (10), these projectors Π W h z(t) and Π Q h p(t) solve the discrete Darcy problem (4) for all w ∈ W h , q ∈ Q h with right-hand sides given by (g, w) = c(z(t), w) + b(w, p(t)) and (s, q) = b(z(t), q).For an a-priori analysis based on such a Galerkin-projection approach to be optimal, including in the limit as κ → 0, the continuity constants C b , C c and the Babuška-Brezzi stability constants α, β must be independent of 0 < κ ≤ 1.
Attaining κ-independent continuity and stability constants is non-trivial for the Darcy problem, and the norms that are selected for W and Q play a vital role.For instance, the standard pairing H(div) × L 2 with the natural norms is not appropriate as e.g.c is not continuous with respect to the H(div) norm: the continuity bound C c depends on κ.However, the following pairings for W × Q are all meaningful for (4) or its dual In particular, the inf-sup condition (15) holds with inf-sup constant β independent of κ for each of these pairings.We remark that p L 2 +κ 1/2 H 1 ≤ p and p κ 1/2 L 2 ≤ p for κ ≤ 1.The κ-independent inf-sup condition for (A) was recently shown in [2], the inf-sup condition for (B) follows directly by a scaling of the flux by κ −1/2 and the pressure by κ 1/2 .Finally, the inf-sup condition of (C) follows directly from Poincare's inequality with a similar scaling as in (B).The boundedness of b(z, p) can be established for each of the pairings above.The pairing of (C) corresponds to the case of the L 2 ×H 1 formulation of the mixed Darcy problem, i.e. b(z, p) = (z, ∇p) with z ∈ W = L 2 and p ∈ Q = H 1 , but boundedness is proved in the same manner as for (B).In the case of (B): applying Cauchy-Schwarz and the weighted norm definitions immediately gives The case of (A) is complicated by the definition of the sum norm on the pressure space Q, and a one-line argument is not possible without additional context; see [2] for details.
Options (A) and (B) above fit naturally with the variational formulation of ( 13) and spaces (9).In the following, we suggest that a natural norm for the Darcy flux is which is equivalent to the norm of the flux in (A) above for the relevant range of κ when τ > 0. However, both options (A) and (B) have disadvantages.For (B), the pressure norm (on Q) becomes progressively weaker as κ nears 0 while the norm of the flux divergence (on W ) is unnecessarily large compared with e.g.(17).The primary drawback to using (A) is that the pressure norm is implicitly defined.This fact means that an a-priori analysis based on the method of projections is more complex to carry out in practice; it is not clear that standard analytic projection techniques [22] could be used directly when the norm of L 2 + κ 1/2 H 1 is chosen for the pressure space.We will argue instead that an a-priori analysis of (13) based on the use of a Galerkin projection of the form (4) is not necessary; thus alleviating the need for an explicit uniformin-κ Darcy stability condition on (W h , Q h ).Neither (15) nor the saddle-point stability of (4) in general play a role in the well-posedness of (13).Condition (iii) of Definition 1 will thus be replaced by a less restrictive condition.An important consequence of relaxing the uniform-in-κ Darcy stability hypothesis is that the standard L 2 -norm on Q can, and will, be used.

Minimal Stokes-Biot conditions yield well-posed schemes
The primary contributions of this section is our definition of minimal Stokes-Biot stability conditions (Definition 2) and the associated Proposition 1 which shows that Euler-Galerkin discretizations of Biot's equations are well-posed under these conditions.In particular, the minimal Stokes-Biot conditions do not assume that W h × Q h is a stable Darcy pairing.
Throughout this section we assume that U , W and Q are defined by (9).The norm on U is taken to be the usual H 1 (Ω)-norm • 1 , the norm on Q is the standard L 2 -norm • , while the norm |||•||| on W is the weighted norm defined by (17).The norm (17) was first introduced in [14, Sec.3.1].The bilinear forms a, b, c, d are as defined by (11).4.1.Minimal Stokes-Biot conditions.We now introduce our set of minimal Stokes-Biot stability conditions.For clarity and completeness (rather than e.g.brevity), we include the precise stability conditions in the definition here.In essence, between Definitions 1 and 2, only condition (iii) changes.

Definition 2. A family of conforming discrete spaces {U
The bilinear form a is continuous and coercive on U h × U h ; i.e. there exists constants C a > 0 and γ a > 0 independent of h such that [5,6]; i.e. in particular there exists an inf-sup constant β S > 0 independent of h such that A more pragmatic comparison, providing an example of lower computational cost and potential for applied interest, is the discretization P 2 × RT 0 × DG 0 .This discretization is both Stokes-Biot stable and minimally Stokes-Biot stable; of note is that P 2 × P 1 × DG 0 is not Stokes-Biot stable but is minimally Stokes-Biot stable.The P 2 ×RT 0 ×DG 0 discretization was a motivating prototype for discretization studied in [22], which was a minimal displacement enrichment of a P 1 × RT 0 × DG 0 mixed element.The comparison between P 2 × RT 0 × DG 0 and P 2 × P 1 × DG 0 motivates Definition 2, and will be studied in Section 6.
Remark 2. It is interesting to note that the Stokes stability and divergence containment conditions, Defn. 2 (ii) and (iii) above, were shown [14, Thm.1] to be the two essential ingredients for establishing stability, for a non-conforming discretization of (1a)-(1c), using the weighted flux norm (17).Though [14] used inf-sup stable (RT 0 × DG 0 ) elements for the fluxpressure pairing, with κ uniformly bounded above and below, it is clear that the important role of Stokes stability alongside (iii) has been recognized by others working on similar problems.

4.2.
A well posedness result for Euler-Galerkin schemes.The variational problem (13) can be shown to satisfy the Brezzi conditions with respect to the weighted norm (17) for the flux; this is the perspective of [14].The Brezzi stability result in [14], though it does involve some delicate manipulations, holds independently of c 0 ≥ 0.
For completeness: we state a straight-forward existence and uniqueness result, independent of κ > 0 and c 0 ≥ 0, for (13).The result has a simple proof and highlights a critical point: the role of the Stokes stability assumption, Definition 2(ii), in the context of a vanishing storage coefficient; i.e. for c 0 = 0. Proposition 1. Suppose that the discrete spaces U h × W h × Q h are minimally Stokes-Biot stable (Definition 2).Suppose that c 0 ∈ [0, 1] and κ ∈ (0, 1] are arbitrary but fixed and that a choice of initial iterates iterates u 0 h , z 0 h and p 0 h has been uniquely determined.Then, for each discrete time t m , for a positive integer m, there exists a unique (u m h , z m h , p m h ) solving (13).Proof.We proceed inductively and establish the existence and uniqueness at each discrete time t m where 1 ≤ m ≤ N , t 0 = 0 and t N = T .For m > 0, the existence and uniqueness of (u m h , z m h , p m h ) will follow from a trivial nullspace of the square and finite dimensional linear operator defined by (13).Let = 0 and assume that the right-hand side terms of (13) are zero.
Case I: First consider the case (13), multiply (13c) by τ , sum the resulting equations, and use the definitions (11) of the forms c and d, to obtain Now, assume that condition (i) of Definition 2 holds; and in particular that a is coercive over U h in the H 1 -norm with coercivity constant γ a (independent of h).Then  (13), multiply (13c) by τ , sum the resulting equations, and use (20) together with the definitions (11) of the forms c and d to obtain The coercivity and continuity assumptions on a (Definition 2(i)), together with Cauchy's inequality with epsilon, (21) and the gathering of like terms yield for any ǫ > 0. Now pick e.g.
and again thus u m h = z m h = p m h = 0.
Remark 3. The result of Proposition 1 relies only on Definition 2(i)-(ii) and does not depend on Definition 2(iii) (nor Definition 1(iii)).

A priori error estimates for minimally Stokes-Biot stable schemes
In this section, we derive a-priori error estimates for the Euler-Galerkin discrete Biot equations (13).The final result is summarized in Proposition 3 of Section 5.4.We will assume the point of view of minimal Stokes-Biot stability as defined by Definition 2 and that U h contains the continuous nodal Lagrange elements P r for some r ≥ 1.We begin by establishing basic assumptions on the spaces U h , W h and Q h , and define projection operators in Section 5.1.

Projections and approximability.
As in the previous section, let U , W , Q be given by ( 9) with norms • 1 , |||•||| cf.(17), and • , respectively.Assume that the discrete spaces U h × W h × Q h satisfy the assumptions of Definition 2. We denote the (continuous) solutions to (10) at time t m by (u m , z m , p m ) for m = 1, 2, . . ., N while (u m h , z m h , p m h ) represent the solutions of the discrete problem (13).For use in the subsequent error analysis, we make basic assumptions on the spaces U h , W h and Q h , and define projection operators W h : Assume that W h contains (at least) piecewise polynomial (vector) fields of order k = k W ≥ 0. We assume the existence of a generic discrete interpolant Π W h : W → W h satisfying either for w ∈ H k+1 (Ω).The estimates (26) are characteristic of a Raviart-Thomas type, RT k (k = 0, 1, 2, . . .), interpolant whereas (27) could correspond to a continuous Lagrange interpolant of order k ≥ 1 [8].U h : Following [22], we define Π U h : U → U h as a modified elliptic projection satisfying for u ∈ U : where q ∈ Q is given and will, in practice, be selected as the exact pressure solution to (10) at given times.Assume that U h contains (at least) continuous piecewise polynomial (vector) fields of order k U ≥ 1.There then exists an interpolant, for all u ∈ H k U +1 , c.f. e.g [8].Then for u ∈ U we have Using assumption (i) of Definition 2 and (28 Combining the above with assumption (25) gives where q ∈ Q is the fixed function defining the elliptic projection (28).

Interpolation notation and identities.
Following standard notation [13,16,22], the error at time t m > 0 can be decomposed into interpolation errors ρ and approximation errors e: The interpolation errors satisfy the following identities.Since div W h ⊆ Q h and by the definition of the

Similarly, by the definition of Π
where we recall the discrete derivative notation (5).Finally, (28) directly gives Taking the difference between the continuous equations (10) and discrete scheme (13), after multiplying (10b) by τ , combined with the cancellations (31)-(33), yield the following error equations at t m : (e m u , e m z , e m p ) satisfies by way of the general identity (36) and similarly for p.

Discrete approximation error estimates.
In this section we estimate the discrete errors described by (34) in their respective norms; that is, e m u 1 , |||e m z ||| and e m p .In contrast to e.g.[22], we do not make use of the restrictive uniform-in-κ Darcy stability assumption.In turn, the error equations require a more technical analysis and we have adapted related methods originally used to study κ fixed [16] and vanishing (c 0 ) storage coefficient.Despite the more technical approach, the resulting estimates presented in Proposition 3 is directly comparable to related results in the literature; c.f. [13,Lem. 3], [16,Thm. 4.1] and [22,Thm 4.6].We conclude that the concept of minimal Stokes-Biot stability provides analogous error estimates for a more general set of conforming discrete spaces than the original Stokes-Biot stability concept.
During the course of the analysis will make use of the following useful inequality Lemma 1. [16, Lemma 3.2] Suppose that A, B, C > 0 and D ≥ 0 satisfy For any (continuous) symmetric bilinear form a with induced norm • a we have the inequality [9] (39 Using the above, and the symmetry of both a(•, •) and d(•, •), it follows that the left-hand side of ( 38) is bounded below by with inequality constant depending on Estimate of e m u a : Following a technique from [16], let J be the integer index where e m u a (for m = 1, . . ., N ) obtains its maximal value.Summing (44) from m = 1 to m = J, using the maximality assumption, and re-arranging terms yields We can apply Lemma 1 to (45) by taking A = e J u a , B = e J p d and dropping the additional left-hand side term; then we choose Provided sufficient regularity, we have that Bound of τ R m : We now develop a bound for the terms τ R m ; c.f. (35).From the fundamental theorem of calculus and integration by parts we have the general result for any m = 1, . . ., N , assuming sufficient temporal regularity of the field f .We therefore, again under the assumption of sufficient spatial and temporal regularity, have the inequalities which control the first and fourth terms of R m .
For the second term of R m we have div Rearranging the terms of ∂ τ ρ m u , applying the fundamental theorem of calculus and using the commutation of the time derivative with the elliptic projection (28) yields ( 49) For the third term of R m we have, again up to sufficient regularity, that Summarizing, (47)-(50) thus yield And so, the estimate (46) becomes Clearly, by Definition 2(i), this also gives a bound for e m u 1 (depending on γ −1 a ) for m = 1, . . ., N .
Estimate of e m p : The norm e J p d in e.g.(52) vanishes in the limit as c 0 → 0. An alternative bound for e m p can be derived from the Stokes stability assumption, Definition 2(ii).In particular, using (42) and (52) it follows that for each 1 ≤ m ≤ N : We begin by considering the first component and again argue based on maximality.Take the difference of the error equation (34a) at time levels m, m − 1 and dividing by τ to get Similarly taking the difference of (34b) at time levels m and m − 1, and divide by τ 2 to get w = e m z in the above as well as q = −∂ τ e m p in (34c); summing these three equations, using Cauchy-Schwarz on the right-hand side, and coercivity on the left-hand side gives From Definition 2(ii) and (54) we have that

S
∂ τ e m u 1 by the analogue of (42).Using this on the right-most term of the above, alongside Cauchy's inequality with epsilon and choosing epsilon appropriately, gives with inequality constant depending on C a β S γ −1 a .Dropping the positive displacement and pressure left-hand side terms, multiplying both sides by τ , and using the symmetry of c together with the inequality (39) give Let M be the index where e m z 2 c achieves its maximum for 1 ≤ m ≤ N .Summing the above from m = 1 to m = M , using the maximality of e M z c , multiplying both sides by τ and re-arranging yields ( 55) By the fundamental theorem of calculus, we have Applying Hölder's inequality on the right-most term, above, gives Inserting (56) and ( 51) into (55), using Young's inequality on the second term on the righthand side and rearranging yields Estimate of τ 2 div e m z : Now we estimate the second, and final, term in the flux norm (17).Let K denote the index where div e K z is maximal.Using Definition 2(iii), and selecting q = τ div e K z in the error equation (34c Thus, re-arranging terms, using Cauchy-Schwarz and the triangle inequalities, and dividing by div e K z gives where the last inequality follows from the majorization of the terms e K u , e K p , e K−1 u , e K−1 p by the maximum e J u = max j=1,2,...,N e j u and inequality (53).Noting that and employing (52), (51) and taking I = max {J, K} then gives Finally, to establish (37), combine the definition of the weighted flux norm ( 17), ( 52), (53) (57), and (58) and use the fact that the integral from 0 to T majorizes all of the time-integral right-hand sides of the summed expressions.

Convergence estimates.
To specialize the general results of Proposition 2 we will first suppose the exact solutions to (10a)-suitable regularity assumptions.Moreover, we assume the interpolants, discussed in 5.1, satisfy approximation inequalities of a certain order.Towards that end let U h × W h × Q h satisfy the assumptions of Definition 2.
For a reflexive Banach space X, a time interval (a, b) ⊆ R and a measurable f : (a, b) → X we define the canonical space-time norm [10] (59) As in the case of spatial derivatives, the usual Sobolev notation f ∈ H r (a, b; X) means that f ∈ L 2 (a, b; X) and that ∂ t f , ∂ 2 t f , . . ., ∂ r t f are also in L 2 (a, b; X).In the sections that follow we will sometimes use the abbreviations f L 2 X or f H r X to signify (59).Proposition 3. Suppose the assumptions of Proposition 2 hold.Let k ≥ 0 be the greatest integer such that the orthogonal projection, Π Q h : Q → Q h , satisfies (25).Suppose r ≥ 1 is the maximal integer such that P r , the space of continuous Lagrange polynomials of order r, is contained in U h ; suppose an interpolation, from W to W h , satisfying either (26) or (27) exists and let s > 0 be the maximal integer satisfying the respective inequality.Suppose that the exact solutions to (10a)-(10c) satisfy the regularity assumptions and that the initial iterates, (u 0 h , z 0 h , p 0 h ), satisfy the estimates consistent with the projections of section 5.1.Then for c = min{k, r, s} we have where M 1 and M 2 are given by Proof.First, note that since Π Q h satisfies (25) and since P r ⊂ U h then, according to the argument directly preceding (29), the inequality (29) holds.Using the triangle inequality, , along with (60) and the projection estimates of section 5.1, applied to the last three terms above, gives (62) e 0 u 1 + τ 1/2 e 0 z c + e 0 p d h r u(0) Then (61) follows from the triangle inequality, with respect to the error decompositions (30), along with: the discrete error estimates (37); discrete initial iterate error estimates (62); and interpolation estimates ( 25), ( 26)-( 27) and (28).
Remark 4. Further assumptions on the discrete spaces, beyond the minimal Stokes-Biot stability of Defn.2, can lead to slightly different versions of Proposition 3.For instance, if (W h , Q h ) are such that the usual Raviart-Thomas type projection commutation relation holds for each z ∈ W then (div ρ z , q h ) = 0 so that, for instance, the contribution z L 1 H s+1 vanishes from M 1 ; this term arises from div ρ z in (37).This observation is used in [22] where W h = RT 0 is fixed.
6.1.Convergence of a Stokes-Biot stable pairing.We first consider the convergence properties for the pairing U h × W h × Q h = P 2 (T h ) × RT 0 (T h ) × DG 0 (T h ).We report on the relative approximation errors for the displacement ũ(T ) − u h (T ) 1 / ũ(T ) 1 , pressure p(T ) − p h (T ) / p(T ) and flux |||z(T ) − z h (T )|||/|||z(T )||| for a series of uniform meshes T h with mesh size h.The exact solutions ũ, p, z were represented by continuous piecewise cubic interpolants in the error computations.
6.1.1.Vanishing storage c 0 = 0, varying conductivity 0 < κ ≤ 1. (see: Table 2) We observe that the displacement error converges at the expected and optimal rate (2) for κ ranging from 1 down to 10 −12 .Overall, the displacement errors remain essentially unchanged as c 0 and κ vary.(We therefore do not report or discuss these further here.)The behaviour for the flux and pressure errors is less regular.The flux and pressure approximation errors increase as κ decreases, but seem to stabilize i.e. not increase substantially further from κ = 10 −4 to 10 −8 and to 10 −12 .Moreover, for each κ, the pressure and flux errors decrease with decreasing mesh size.Indeed, for h = 1/128, the pressure errors are of similar magnitude for the range of κs tested.6.1.2.Fixed storage c 0 = 1, varying conductivity 0 < κ ≤ 1. (see: Table 3) For this case, we again observe that the flux and pressure approximation errors increase as κ decrease, but seem to stabilize and not increase substantially further from κ = 10 −4 to 10 −8 and 10 −12 .Again, for each κ, the pressure and flux errors decrease with decreasing mesh size and for h = 1/128, the pressure errors are nearly identical for the range of κs tested.4) For this case, we observe nearly uniform behaviour as c 0 decreases.The pressure and flux errors are similar for the range of c 0 s considered, and converge at the optimal and expected rate (1).6.2.Convergence of a minimally Stokes-Biot stable pairing.We now turn to consider the convergence properties for the pairing and again report on the relative approximation errors for the displacement, pressure and flux.This pairing does not satisfy the original Stokes-Biot stability criteria, but do satisfy the minimally Stokes-Biot stable criterion.Numerical results for this minimally Stokes-Biot stable discretization, for the three paradigms considered in Section 6.1, are presented in 6.2.1-6.2.3 alongside specific comparisons to the Standard Stokes-Biot stable case.The results of this comparison will supply computational evidence that Definition 1(iii) can be replaced by Definition 2(iii) while retaining the convergence properties first observed in [13,22].6.2.1.Vanishing storage c 0 = 0, varying conductivity 0 < κ ≤ 1. (see: Table 5) Comparing Table 5 with Table 2, we observe that the performance of the two element pairings is almost surprisingly similar.Again, the displacement converges at the optimal and expected rate (2), the pressure and flux errors increase with decreasing κ, but stabilize, and converge with decreasing mesh size.We further observe that the relative errors for the flux for this element pairing is smaller than for the P 2 × RT 0 × DG 0 case (bottom rows).
Relative approximation errors for the displacement ũ(T ) − u h (T ) 1 / ũ(T ) 1 (top rows), pressure p(T ) − p h (T ) / p(T ) (middle rows) and flux |||z(T ) − z h (T )|||/|||z(T )||| (bottom rows) for varying κ on a series of uniform meshes T h with mesh size h.The last column 'Rate' denotes the order of convergence using for the last two values in each row.c 0 = 0. Timedependent test case as presented given in Section 6.The displacement errors for κ = 10 −4 , 10 −8 were identical to the data presented (κ = 1, κ = 10 Relative pressure approximation errors (top rows), and flux errors (bottom rows), as specified in caption of Table 2.
In this manuscript we have shown that the Stokes-Biot perspective can be relaxed; in particular we differ from previous authors [13,22] by departing from a uniform-in-κ Darcy stability assumption in our analysis.In fact, an analysis based on a uniform-in-κ Darcy stability assumption should take into account pressure-space norms of the form L 2 + kH 1 as described in [2].It is not entirely clear that the L 2 + kH 1 norm can be treated with the usual methods presented here and in other related work [13,14,16,22].By removing the uniform-in-κ assumption: Proposition 3 solidifies, and generalizes, previous convergence estimates [22] and broadens the original view of Stokes-Biot stability to include alternative spaces that may not be Darcy stable; even for a fixed choice of κ.

2. Notation and preliminaries 2 . 1 .
Sobolev spaces and norms.Let Ω ⊂ R d for d = 1, 2, 3 be an open and bounded domain with piecewise C 2 boundary The classical flux-pressure pairings, e.g.RT k × DG k or BDM k+1 × DG k for k = 0, 1, 2, . . ., satisfying Definition 1(iii) also satisfy the conditions of minimal Stokes-Biot stability; in particular Definition 2(iii).However, the minimal Stokes-Biot condition also includes discretizations which are not encompassed by Definition 1.For instance: flux-pressure pairings where the flux is taken from the space of continuous Lagrange polynomials can satisfy Definition 2 while not satisfying Definition 1.An illustration of this can be found in the family of discretizations where the displacement-pressure pairing are of Scott-Vogelius type; these either have the formP k × RT m × DG k−1 or P k × P m × DG k−1 where k ≥ 4 and 0 ≤ m ≤ k − 1.The flux-pressure pairings R m × DG k−1 , for m < k − 1, and P m × DG k−1 , for m ≤ k − 1, are not Darcy stable but do satisfy the minimal Stokes-Biot stability containment condition of Definition 2(iii).

and thus u m h = z m h = p m h = 0 .
Case II: If c 0 = 0, assume that also the Stokes stability assumption (Definition 2(ii)) holds.Then, for each p m h ∈ Q h , there exists (see e.g.[5, p. 136]) a y ∈ U h such that (div y, p m h )

||| e 0 u 1 + e 0 p d + τ 1 2 c
minimally Stokes-Biot stable (by satisfying the assumptions of Definition 2).Then, the discrete approximation errors (e m u , e m z , e m p ) described by (34) satisfy the inequality: (37) e m u 1 + e m p + |||e m z ds + C T τ , with inequality constant depending on C a , γ −1 a and where C m τ ≡ tm 0 div ρ z + ρ ∂tu 1 + τ (c 0 ∂ tt p + ∂ tt u 1 ) ds. Proof.In an analogous fashion as for Proposition 1, multiplying (34c) by τ , selecting v = e m u − e m−1 u , w = e m z , and q = −e m p in (34) and summing gives a(e m u − e m−1 u , e m u ) + τ c(e m z , e m z ) + d(e m p − e m−1 p , e m p ) = τ c(ρ m z , e m z ) − τ R m , e m p .(38)

1 ≤ C a e m u 1 . 2 a 2 a − e m−1 u 2 a + τ e m z 2 c + e m p 2 d − e m−1 p 2 d τ ρ m z 2 c
From the Stokes stability assumption (19) and (34a) we have the estimate (42) β S e m p ≤ sup v∈U h b(v, e m p ) v 1 = sup v∈U h −a(e m u , v) v Then, Cauchy-Schwarz, (42) and the coercivity of a gives (43) τ | R m , e m p | ≤ C a β −1 S γ −1/+ R m e m u a .

2 c 2 c and τ 2
and where J is the index where e J u 1 is maximal.Thus e m p can be bounded by the right hand side of (52), independently of c 0 .Estimate of τ em z : In order to estimate the flux error in the norm defined by (17), i.e. |||e m z |||, it will be advantageous to consider the constituents separately; e.g.τ e m z div e m z 2 .