Page 1

INF-SUP STABLE DISCRETIZATION OF THE

QUASI-STATIC BIOT’S EQUATIONS IN POROELASTICITY

CHRISTIAN KREUZER AND PIETRO ZANOTTI

Abstract. We propose a new full discretization of the Biot’s equations in

poroelasticity. The construction is driven by the inf-sup theory, which we

recently developed. It builds upon the four-field formulation of the equations

obtained by introducing the total pressure and the total fluid content. We

discretize in space with Lagrange finite elements and in time with backward

Euler. We establish inf-sup stability and quasi-optimality of the proposed

discretization, with robust constants with respect to all material parameters.

We further construct an interpolant showing how the error decays for smooth

solutions.

1. Introduction

This paper is the second one in a series initiated by [24], regarding the analysis

and the discretization of the quasi-static Biot’s equations in poroelasticity. (See

(2.1) below for the statement of the problem). The series centers around the use of

the inf-sup theory for the stability and the error analysis, with the aim of highlight-

ing the possible advantages stemming from the proposed approach, which appears

to be new in this context.

The inf-sup theory is a framework for the analysis of general linear variational

problems. The main result therein is the so-called Banach-Ne˘cas theorem (see, e.g.,

[14, Section 25.3]), that characterizes the well-posedness of such problems. Suc-

cessful applications of the Banach-Ne˘cas theorem additionally provide a two-sided

stability estimate, entailing that the space of all possible solutions is isomorphic

to that of all possible data, cf. Theorem 2.4 below. Such an estimate is of special

interest, as it is the starting point for the derivation of sharp a posteriori error

estimates, since it establishes the equivalence of error and residual, cf. [43, Sec-

tion 4.1.4]. Moreover, if the same approach applies also to the discretization at

hand, then the resulting stability estimate is the starting point for the derivation

of sharp a priori error estimates, see [2] for conforming discretizations and [4] for

nonconforming ones. The inf-sup theory is useful also in the design of robust pre-

conditioners [20, 30] and in the convergence analysis of some adaptive procedures

[16].

The inf-sup theory is well-established for the analysis and the discretization of

stationary linear equations. We refer to [14] for an overview with several examples.

The situation is substantially different for evolution equations, that are usually an-

alyzed by other techniques. While the application of the inf-sup theory to the heat

equation has been recently considered by various authors (see, e.g., [15, Chapter 71]

2010 Mathematics Subject Classification. Primary 65M15, 65M60; Secondary 74F10, 76S05.

Key words and phrases. Inf-sup stability, quasi-static Biot’s equations, poroelasticity, quasi-

optimality, robustness, a priori analysis.

1

arXiv:2407.02939v1 [math.NA] 3 Jul 2024

2

C. KREUZER AND P. ZANOTTI

and the references therein), we are aware of only a few results for other prototypical

problems, like the wave and the Schrödinger equation [39].

The quasi-static Biot’s equations in poroelasticity are not an exception in this

respect. In fact, on the one hand, some authors have used the inf-sup theory in the

analysis and the discretization of the stationary problem resulting from the semi-

discretization in time, see e.g. [21, 23, 27]. On the other hand, in the quasi-static

case, we are not aware of any result obtained by this approach. For the a posteriori

analysis, Li and Zikatanov [25] used, in a sense, an equivalent argument. For the a

priori analysis, all papers we are aware of build upon a different argument, namely

an energy estimate that seems to go back to a seminal contribution of Zenıšek [46].

A far by exhaustive list includes [3, 8, 9, 22, 26, 33, 34, 44].

To our best knowledge, our series is the first contribution making a systematic

use of the inf-sup theory in the design and in the error analysis of a discretization

of the quasi-static Biot’s equations. Our first paper [24] is devoted to the analysis

of the equations. It establishes the well-posedness as well as a two-sided stability

estimate, which is robust with respect to all material parameters. The latter con-

tributions and the functional setting distinguish our results from previous ones by

Zenıšek [46] and Showalter [38]. In particular, we consider an equivalent four-field

formulation of the equations, that is obtained from the original one by introducing

the so-called total pressure and total fluid content. In [24] we prove also that, in

certain circumstances, additional regularity in space of the data imply correspond-

ing additional regularity in space of the solution. Establishing a similar result in

time is more challenging, due to possible singularities at the initial time; compare

e.g. with [31, 38] and the discussion in Remark (4.10) below. Further previous

contributions to the regularity theory are in [6, 45].

In this paper we propose a backward Euler discretization in time and an abstract

discretization in space of the Biot’s equations. We establish the well-posedness and

a two-sided stability estimate via the inf-sup theory, by mimicking the argument in

[24]. Then we consider a simple realization of the abstract discretization in space,

making use of Lagrange finite elements for all variables. We prove a quasi-optimal

a priori error estimate, meaning that the error is equivalent to (i.e. bounded from

above and below by) the best error. To our best knowledge, it is the first time that

such a result is established for a discretization of the Biot’s equations.

The error notion we consider is motivated by the inf-sup theory and it is relatively

involved, as all variables are coupled in a nontrivial way. Therefore, we further

elaborate on our error estimate, by showing that the best error can be bounded

by a sum of decoupled best errors in standard norms, that are much easier to

investigate. All constants in our results are robust with respect to all material

parameters and do not require any additional regularity beyond the one guaranteed

by the well-posedness of the equations. Finally, we establish first-order convergence,

with respect to the space and time meshsize, for sufficiently smooth solutions.

Organization. In section 2 we state the equations and recall the main results

from [24]. Section 3 establishes the stability of an abstract discretization. Section 4

is devoted to the a priori error analysis for an exemplary concrete discretization,

which is also tested numerically in section 5.

Notation. Throughout the paper, we denote by ∥·∥X the norm of a Hilbert space

X. The dual space X

∗ acts on X through the pairing ⟨·, ·⟩X

. We denote by L2(X),

INF-SUP STABLE DISCRETIZATION OF THE BIOT’S EQUATIONS

3

H1(X) and C0(X), respectively, the Bochner spaces of all L2, H1 and C0 functions

mapping the time interval [0,T] into X. For a measurable set Ω ⊆ Rd, we adopt

the simplified notation (·, ·)Ω and ∥·∥Ω for the scalar product and the norm in

L2(Ω). We write a ≲ b and a ≂ b when there are constants 0 < c ≤ c, possibly

different at different occurrences, such that a ≤ cb and cb ≤ a ≤ cb, respectively.

The hidden constants are independent of the material parameters involved in the

equations. The dependence on other relevant quantities is addressed case by case,

see e.g. Remark 4.1.

2. Inf-sup theory for the Biot’s equations

This section introduces the Biot’s equations and summarizes some results from

[24], that are useful for the construction and the analysis of the discretization in

the next sections.

2.1. Initial-boundary value problem. Let Ω ⊆ Rd, 1 ≤ d ≤ 3, be a bounded

domain with polyhedral and Lipschitz continuous boundary. The flow of a Newto-

nian fluid inside a linear elastic porous medium located in Ω, in the time interval

(0,T) with T > 0, is modeled by the quasi-static Biot’s equations as follows

(2.1a)

−div(2µ∇Su + (λdivu − αp)I) = fu

in Ω × (0,T)

∂t(αdivu + σp) − div(κ∇p) = fp

in Ω × (0,T).

The first equation states the momentum balance for the elastic porous medium,

whereas the second one is the mass balance for the fluid.

The unknown functions in the equations are the displacement u : Ω → Rd of the

elastic porous medium and the pressure p : Ω → R of the fluid. The differential

operator ∇S is the symmetric part of the gradient and I is d×d identity tensor. The

material parameters, denoted by Greek letters, are the Lamé constants µ,λ > 0,

the Biot-Willis constant α > 0, the constrained specific storage coefficient σ ≥ 0

and the hydraulic conductivity κ > 0. Consistently with [24], we assume that all

parameters are constant in Ω × [0,T] for simplicity. We refer to [24, Remark 2.1]

for a discussion on possible extensions.

We are interested in the initial-boundary value problem obtained by prescribing

also the initial condition

(2.1b)

(αdivu + σp)|t=0 = ℓ0

in Ω

as well as the boundary conditions

(2.1c)

u = 0

on Γu,E × (0,T)

(2µ∇Su + (λdivu − αp)I)n = gu

on Γu,N × (0,T)

p = 0

on Γp,E × (0,T)

κ∇p · n = gp

on Γp,N × (0,T)

where Γu,E ∪ Γu,N = ∂Ω=Γp,E ∪ Γp,N and Γu,E ∩ Γu,N = ∅ = Γp,E ∩ Γp,N . The

letter n denotes the outward unit normal vector on ∂Ω. The subscripts ‘E’ and

‘N’ are intended to assist the reader in distinguishing the portions of the boundary

with essential and natural conditions.

We point out that different statements of the initial and of the boundary con-

ditions are sometimes met in the literature. We refer to [24, Remark 2.2-2.3] for a

more extensive discussion on this point.

4

C. KREUZER AND P. ZANOTTI

2.2. Weak Formulation and well-posedness. For convenience, we introduce a

compact notation for the differential operators in the Biot’s equations (2.1), namely

(2.2)

E := −div(2µ∇S)

D := div

L := −div(κ∇).

Notice that E and L act on u and p, respectively, in (2.1a). Comparing also with

the boundary conditions, it looks reasonable that the regularity of u in space in a

weak formulation can be described via

(2.3)

U :=

(

H1(Ω)d/RM if Γu,N = ∂Ω

H1

Γu,E

(Ω)d

otherwise

with RM denoting the rigid body motions. Analogously, for the regularity of p in

space, we consider

(2.4)

P :=

⎛

⎢

⎢k

H1(Ω) ∩ L2

0(Ω)

if Γp,N = ∂Ω

H1

Γp,E

(Ω) ∩ L2

0(Ω) if Γp,N ̸= ∂Ω, Γu,E = ∂Ω,σ = 0

H1

Γp,E

(Ω)

otherwise

where L2

0(Ω) = {q ∈ L2(Ω) | /

Ω

q = 0}. We refer to [24, Remark 2.5] for a

motivation of the nonstandard definition of P in the second case.

Remark 2.1 (Notation). We are aware of the fact, that the abstract notation in-

troduced here (and the subsequent one for all related spaces and operators) makes

the reading possibly harder. Still, it has the advantage that all combinations of the

boundary conditions (as well as the critical case σ = 0) can be treated at the same

time. In our view, this is quite important, because our approach to the analysis

and the discretization of the Biot’s equations is mainly the same in all cases, but

each case may require subtle minor modifications.

The action of the divergence on u in (2.1a) indicates that also the space

(2.5)

D := D(U) =

(

L2

0(Ω) if Γu,E = ∂Ω

L2(Ω) otherwise

plays a relevant role in the Biot’s equations. Similarly, since p is involved in (2.1a)

also without the action of any differential operator in space, we repeatedly make

use of

(2.6)

P =

(

L2

0(Ω) if Γp,N = ∂Ω or Γu,E = ∂Ω,σ = 0

L2(Ω) otherwise,

where the closure is taken with respect to the L2(Ω)-norm. An important point for

our analysis is that both D and P are subspaces of L2(Ω), but their mutual relation

depends on the boundary conditions and on the parameter σ. Therefore, it is useful

introducing

PD : L2(Ω) → D

and

P

P

: L2(Ω) → P,

the L2(Ω)-orthogonal projections onto D and P, respectively.

Remark 2.2 (Functional setting). The diagram in Figure 2.1 summarizes the rela-

tion among the spaces and the operators introduced up to this point. In addition,

we denote by i : P → P the inclusion operator and D∗ and i∗ are the adjoint oper-

ators of D and i, respectively. The spaces D and P are identified with their duals

INF-SUP STABLE DISCRETIZATION OF THE BIOT’S EQUATIONS

5

via the L2(Ω)-scalar product. Thus, the square on the right side of the diagram

involves the Hilbert triplet

(2.7)

P ↩→ P ≡ P

∗

↩→ P∗.

U

∗ oo

E

OO

D∗

U

D

L2(Ω)

PD

~

PP

P

L

//

i

P

∗OO

i∗

D

∗

D

P

P

∗

Figure 2.1. Spaces and operators describing the regularity in

space for the weak formulation (2.8) of the Biot’s equations. The

triple lines on the bottom indicate identification via the L2(Ω)-

scalar product.

With this preparation, we are in position to recall the weak formulation of the

initial-boundary value problem (2.1) introduced in [24, section 2.3], namely

(2.8)

Eu + D∗ptot = ℓu

in L2(U∗)

λDu − ptot − αPDp = 0

in L2(D)

αP

P

Du + σp − m = 0

in L2(P)

∂tm + Lp = ℓp

in L2(P∗)

m(0) = ℓ0

in P∗.

The loads ℓu ∈ L2(U∗) and ℓp ∈ L2(P∗) result from the data fu and fp in the

equations (2.1a) as well as gu and gp in the boundary conditions (2.1c).

Remark 2.3 (Auxiliary variables). Compared to (2.1a), the weak formulation (2.8)

involves two additional unknown variables, namely the total pressure ptot and the

total fluid content m defined by the second and third equation, respectively. Intro-

ducing these variables is not strictly necessary for our analysis, but it substantially

simplifies the definition of the space Y1 in (2.9) below. The use of the total pressure

was observed in [27] to help also the design of robust linear solvers.

The main result in [24] states that (2.8) is uniquely solvable within the closure

Y1 of the space

(2.9)

Y1 := L2(U) × L2(D) × L2(P) × (L2(P) ∩ H1(P∗))

with respect to the norm

∥(Cu, Cptot, Cp, Cm)∥2

1 :=

∫ T

0

(

∥Cu∥2

U

+

1

µ

∥Cptot∥2

Ω + ∥∂t

Cm + LCp∥

2

P∗

)

+ ∥ Cm(0)∥2

P∗

+

∫ T

0

( 1

µ + λ

∥λDCu − Cptot − αPD Cp∥

2

Ω + γ∥αP

PDCu + σCp− Cm∥2

Ω

)

.

(2.10)

6

C. KREUZER AND P. ZANOTTI

Here U and P are equipped with the energy norm

(2.11)

∥·∥2

U

= ⟨E·, ·⟩

U

and

∥·∥2

P

= ⟨L·, ·⟩

P

and the parameter γ is given by

γ :=

⎛

⎢⎢⎢⎢⎢⎢

⎢⎢⎢⎢⎢⎢k

min

{µ + λ

α2

,

1

σ

}

if σ > 0 and P ⊆ D

µ + λ

α2

+

1

σ

if σ > 0 and P⊈D

µ + λ

α2

if σ = 0.

(2.12)

Taking the closure is indeed necessary, because Y1 is not closed with respect to

∥·∥1, cf. [24, Proposition 4.4].

Theorem 2.4 (Well-posedness of the weak formulation). For all possible data

(ℓu,ℓp,ℓ0) ∈ L2(U∗)×L2(P∗)×P∗, the weak formulation (2.8) has a unique solution

y1 = (u, ptot, p, m) ∈ Y1, which satisfies the two-sided stability bound

∥y1∥2

1 ≂

∫ T

0

(∥ℓu∥2

U∗ + ∥ℓp∥2

P∗) + ∥ℓ0∥2

P∗ .

Moreover, we have m ∈ C0(P∗) as well as the norm equivalence

∥y1∥2

1 ≂ ∥y1∥2

1 + ∥m∥2

L∞(P∗) +

∫ T

0

(λ∥Du∥2

Ω + γ−1∥p∥2

Ω) .

All hidden constants depend only on Ω and T.

Proof. Combine [24, Theorem 3.5] with [24, Proposition 4.1].

□

Although we omit the proof of Theorem 2.4, it is worth roughly summarizing

how it is obtained. Indeed, we shall make use of a similar argument in order to

verify the well-posedness of the discretization introduced in the next section, cf.

Theorem 3.9 below.

The weak formulation (2.8) can be viewed as a special instance of the following

linear variational problem: given ℓ ∈ Y∗

2, find y1 ∈ Y1 such that

(2.13)

b(y1,y2) = ⟨ℓ, y2⟩

Y2

∀y2 ∈ Y2.

The test space is obtained by collecting all possible test functions for (2.8), namely

Y2 := L2(U) × L2(D) × L2(P) × L2(P) × P

and it is equipped with the norm

∥(v, qtot, q, n, n0)∥2

2 :=

∫ T

0

(

∥v∥2

U

+ ∥n∥2

P

)

+ ∥n0∥2

P

+

∫ T

0

((µ + λ)∥qtot∥2

Ω + γ−1∥q∥2

Ω) .

(2.14)

The bilinear form b : Y1 × Y2 → R is defined according to the left-hand side of

(2.8) by

b(Cy1,y2) :=

∫ T

0

(

⟨ECu + D∗

Cptot,v⟩

U

+ ⟨∂t

Cm + LCp, n⟩P

)

+ ⟨ Cm(0),n0⟩

P

+

∫ T

0

(

(λDCu − Cptot − αPD Cp, qtot)Ω + (αP

PDCu + σCp− Cm, q)Ω

)

(2.15)

INF-SUP STABLE DISCRETIZATION OF THE BIOT’S EQUATIONS

7

for Cy1 = (Cu, Cptot, Cp, Cm) ∈ Y1 and (v, qtot, q, n, n0) ∈ Y2.

The so-called Banach-Ne˘cas theorem characterizes the well-posedness of the lin-

ear variational problem in terms of boundedness, inf-sup stability and nondegen-

eracy of the form b, see e.g. [14, theorem 25.9]. These properties are verified in

[24, section 3]. Their combination implies the well-posedness of (2.8) as well as the

two-sided stability bound in Theorem 2.4.

Remark 2.5 (Trial functions). As in (2.10) and (2.15), we use hereafter the super-

script ‘∼’ to distinguish a general trial function in Y1 from the solution of the weak

formulation (2.8).

Remark 2.6 (Functions and functionals). Owing to Remark 2.2, we often identify

the functions in D and P with their Riesz representatives in D∗ and P

∗

and vice

versa. For instance, the term D∗ptot from the first equation in (2.8) and even the

space L2(P) ∩ H1(P∗) proposed for the total fluid content must be interpreted in

this way. As usual, we omit to explicitly indicate the Riesz isometries, accepting

some ambiguity, so as to alleviate the notation. We apply the same principles to

the discretization in the next sections.

3. Abstract inf-sup stable discretization

In this section we design a discretization of the weak formulation (2.8) of the

initial-boundary value problem for the Biot’s equations (2.1). The space discretiza-

tion is challenging, due to the nontrivial coupling of the variables and the various

differential operators acting on them. Therefore, we initially work with a general

discretization, so as to allow for the maximal flexibility. We discuss a concrete real-

ization in section 4 below. The time discretization seems less problematic, because

(2.8) involves only one time derivative. Hence, we directly make a concrete choice,

namely the backward Euler scheme.

The space discretization in this section builds upon a number of assumptions

that must be verified case by case. The set of our assumptions is identified by

special tags with the dedicated enumeration (H1), (H2), etc. With a small abuse,

we actually include among the assumptions also the definition of some relevant

constants. In those cases, the size (and not just the existence) of the constants is

the property to be verified in each concrete example.

The main result in this section is the well-posedness established in Theorem 3.9.

We do not attempt at analyzing the error at this level of generality. Indeed, we

believe that the estimates we would obtain either require too many assumptions

or are too much abstract to be of interest. Thus, we prefer to discuss the error

analysis for the concrete example in section 4.

Remark 3.1 (Notation for the discretization). In general, we mark all spaces and

operators related to the discretization in space by the subscript ‘s’ and those related

to the discretization in time by the subscript ‘t’. The combination ‘st’ of the

two subscripts identifies the full discretization in space and time. To alleviate the

notation, we use capitol letters (in place of subscripts) to distinguish the functions

involved in the discretization from those related to the original Biot’s equations.

3.1. Abstract discretisation in space. The general concept of our discretiza-

tion in space consists in replacing all spaces and operators in Figure 2.1 by finite-

dimensional counterparts, while preserving the structure of the diagram therein, cf.

8

C. KREUZER AND P. ZANOTTI

Figure 3.1. In order to discretize the displacement and the pressure, we consider

two finite-dimensional linear spaces

(3.1)

Us

and

Ps

i.e. discrete counterparts of the spaces U in (2.3) and P in (2.4). We replace the

operators E and L in (2.2) by positive definite and self-adjoint linear operators

(3.2)

Es : Us → U∗

s

and

Ls : Ps → P∗

s .

In analogy with (2.11), we equip the two spaces with the energy norms

(3.3)

∥·∥2

Us := ⟨Es·, ·⟩Us

and

∥·∥2

Ps := ⟨Ls·, ·⟩Ps .

Then, the dual norms on U∗

s and P∗

s are given by

(3.4)

∥·∥2

U∗

s

:= sup

V ∈Us

⟨·,V ⟩

Us

∥V ∥Us

= ⟨·, E−1

s

·⟩Us

∥·∥2

P∗

s

:= sup

N∈Ps

⟨·,N⟩Ps

∥N∥Ps

= ⟨·, L−1

s ·⟩Ps .

Notice that Us and Ps are not required to be conforming, i.e. subspaces of U and

P, respectively. Therefore, in order to define an error notion on the sums U + Us

and P + Ps, we assume the following.

(H1)

The norms ∥·∥U and ∥·∥P in (2.11) can be extended to U + Us and

P + Ps with ∥·∥U ≂ ∥·∥Us in Us and ∥·∥P ≂ ∥·∥Ps in Ps.

In order to discretize the space D in (2.5), we consider a linear operator

(3.5)

Ds : Us → L2(Ω)

i.e. a discrete counterpart of the divergence D in (2.2). Then, we set

(3.6)

Ds := Ds(Us).

The proof of Theorem 2.4 given in [24] exploits the norm equivalence µ∥D∗ · ∥2

U∗ ≂

∥·∥2

Ω in D, that is nothing else than boundedness and surjectivity of D, cf. [14,

Lemma C40]. (Note that here D is identified with its dual via the L2(Ω)-scalar

product, cf. Remark 2.2.) In order to reproduce this property at the discrete level,

we assume the following.

(H2)

There are constants c = c(Us) and C = C(Us) with 0 < c ≤ C

and such that c∥·∥2

Ω ≤ µ∥D∗

s · ∥2

U∗

s

≤ C∥·∥2

Ω in Ds.

Actually, this prescribes that Us/Ds is a stable pair for the discretization of the

Stokes equations. Note that also Ds is identified here with its dual space.

The proof of Theorem 2.4 given in [24] exploits also the density of P in P, that

gives rise to the Hilbert triplet in (2.7). Since Ps is finite-dimensional, we are led

to discretize P by Ps itself, giving rise to the triplet

(3.7)

Ps = Ps ≡ P

∗

s = P∗

s .

Also in this case, the identification of Ps with its dual space is made via the L2(Ω)-

scalar product. Hence, the pairing ⟨·, ·⟩

Ps

coincides with (·, ·)Ω upon identifying the

functionals in P∗

s with their Riesz representative.

INF-SUP STABLE DISCRETIZATION OF THE BIOT’S EQUATIONS

9

Remark 3.2 (Discretization of the Hilbert triplet). As Ps coincides with its closure,

the above Hilbert triplet is trivial from the algebraic viewpoint. Still, the spaces

in it play different roles and are equipped with different norms, depending on their

position, in analogy with (2.7). To alleviate the notation, we omit hereafter the

symbol of the closure.

We denote by PDs

and PPs

the L2(Ω)-orthogonal projections onto Ds and Ps,

respectively. In particular, the adjoint P∗

Ps

of the second projection maps functionals

on Ps into functionals on P. This observation is important for the definition of the

error notion in section 3.4 below, because the trial norm ∥·∥1 in (2.10) involves the

dual norm ∥·∥P∗ . Thus we assume the following, in order to keep the norm of P∗

Ps

under control.

(H3)

There are constants c = c(Ps) and C = C(Ps) with 0 < c ≤ C and

such that c∥·∥P∗

s

≤ ∥P∗

Ps

· ∥P∗ ≤ C∥·∥P∗

s

in P∗

s .

The upper bound here is equivalent to the P-stability of the projection PPs . The

lower bound can be formulated as an inf-sup condition and it is equivalent to the

existence of a bounded right inverse of PPs , cf. [23, Proposition 3.5].

Finally, the proof of the inf-sup stability in Lemma 3.8 below for vanishing σ

makes use of the following assumption.

(H4)

The inclusion Ps ⊆ Ds holds true if σ = 0.

Notice that the spaces P and D in (2.6) and (2.5), respectively, satisfy the same

inclusion.

Remark 3.3 (Spurious pressure oscillations). The combination of the assumptions

(H2) and (H4) prescribes that Us/Ps is a stable pair for the discretization of the

Stokes equations if σ = 0. This property was observed both numerically [18] and

theoretically [29] to be important to prevent from spurious pressure oscillations

in certain regimes. Indeed, forgetting the time derivative for a moment (or, more

precisely, discretizing it in time) the Biot’s equations (2.1a) are close to the Stokes

equations for vanishing σ and small κ.

Figure 3.1 summarizes the relation among the spaces and the operators intro-

duced in this section. As announced, the structure is exactly as in Figure 2.1. We

refer to Remark 2.2 for the details of the notation.

U

∗

s

oo

Es

OO

D∗

s

Us

D

L2(Ω)

PDs

}}

PPs

!!

Ps

Ls

//

i

P

∗

sOO

i∗

D

∗

s

Ds

Ps

P

∗

s

Figure 3.1. Spaces and operators describing the regularity in

space for the discretization of the Biot’s equations. The triple lines

on the bottom indicate identification via the L2(Ω)-scalar product.

10

C. KREUZER AND P. ZANOTTI

3.2. Discretization in time. As announced, we consider a simple first-order dis-

cretization in time, namely backward Euler. To this end, we introduce a partition

(3.8)

0 = t0 < t1 < ··· < tJ = T

of the time interval [0,T] with J ≥ 1. We denote the local time intervals and their

length by

(3.9)

Ij := [tj−1,tj]

and

|Ij| := tj − tj−1

respectively, for j = 1,...,J.

For X ∈ {Us, Ds, Ps}, we consider the space of piecewise time-constant functions

on the above partition with values in X

S

0

t (X) := {X ∈ L∞(X) | X|Ij =: Xj ∈ X, j = 1,...,J}.

Whenever useful, we identify a function X ∈ S0

t (X) with the sequence of its values

(Xj)J

j=1 ⊆ X.

The space S0

t (X) is suitable for a first-order discretization in time of L2(X), i.e.

of the first three components in the trial (2.9) and of the first four components in

the test (2.2). The last component in the trial space (the fluid content) is different,

as it involves H1-regularity in time. We discretize it in time by

S

0

t (Ps) × Ps

where the second component plays the role of the initial value. Then we introduce

a discrete counterpart dt : S0

t (Ps) × Ps → S0

t (P∗

s ) of the time derivative ∂t, in the

vein of the backward Euler scheme, i.e.

(3.10)

dt(b

M, b

M0)Ij :=

b

Mj − b

Mj−1

|Ij|

for (b

M, b

M0) ∈ S0

t (Ps) × Ps and for all j = 1,...,J.

The proof of Theorem 2.4 given in [24] and, in particular, the control of the

point values of the total fluid content hinge on the following integration by parts

rule 2 /

T

0

(∂t

Cm, L

−1

Cm

〉

P= ∥ Cm(T)∥2

P∗ −∥ Cm(0)∥2

P∗ , which holds true for Cm ∈ H1(P∗),

cf. [15, Lemma 64.40]. The operator dt defined above satisfies a similar relation.

Lemma 3.4 (Time discrete integration by parts rule). We have

2

∫ tj

0

⟨dt(b

M, b

M0), L−1

s

b

M⟩Ps ≥ ∥b

Mj∥2

P∗

s − ∥b

M0∥2

P∗

s

for all (b

M, b

M0) ∈ S0

t (Ps) × Ps and j = 1,...,J.

Proof. We exploit (3.10), rearrange terms and recall the second part of (3.4). It

results

2

∫ tj

0

⟨dt(b

M, b

M0), CL−1

s M⟩Ps = 2

j

∑

k=1

⟨b

Mk − b

Mk−1, L−1

s

b

Mk⟩Ps

= ∥b

Mj∥2

P∗

s

+

j

∑

k=1

∥b

Mk − b

Mk−1∥2

P∗

s − ∥b

M0∥2

P∗

s

cf. [15, eq. (67.9)]. This readily implies the claimed inequality.

□

INF-SUP STABLE DISCRETIZATION OF THE BIOT’S EQUATIONS

11

3.3. Full discretization. Combining the discretizations in space and time from

the previous sections, we are in position to propose an abstract full discretization

of the initial-boundary value problem (2.1) for the Biot’s equations.

We consider the trial space

(3.11)

Y1,st := S0

t (Us) × S0

t (Ds) × S0

t (Ps) × S0

t (Ps) × Ps.

Inspired by (2.10) and taking the assumption (H1) in section 3.1 into account, we

equip Y1,st with the norm

∥(CU, CPtot, CP, b

M, b

M0)∥2

1,st :=

∫ T

0

(

∥CU∥2

U

+

1

µ

∥ CPtot∥2

Ω + ∥dt(b

M, b

M0) + Ls CP∥2

P∗

s

)

+ ∥b

M0∥2

P∗

s

+

∫ T

0

( 1

µ + λ

∥λDs CU − CPtot − αPDs CP∥2

Ω + γ∥αPPs Ds CU + σ CP − b

M∥2

Ω

)

.

(3.12)

Remark 3.5 (Equivalent trial norm). According to the assumption (H3) in sec-

tion 3.1, we could equivalently replace the discrete dual norm ∥·∥P∗

s

in the definition

of ∥·∥1,st by ∥P∗

Ps

· ∥P∗ . All the results stated in section 3.4 hold true also in this

case, with the only difference that the hidden constants additionally depend on the

constants in (H3). This observation is important in view of the definition of the

error notion in section 4.2.

For the test space, we proceed similarly and set

(3.13)

Y2,st := S0

t (Us) × S0

t (Ds) × S0

t (Ps) × S0

t (Ps) × Ps.

Recalling again the assumption (H1) in section 3.1, we equip Y2,st with the norm

∥·∥2 in (2.14).

Let (Lu,Lp,L0) ∈ S0

t (U∗

s ) × S0

t (P∗

s ) × P∗

s be a discretization of the correspond-

ing data (ℓu,ℓp,ℓ0) ∈ L2(U∗) × L2(P∗) × P∗ in the weak formulation (2.8) of

the Biot’s equations. We consider the following full discretization of (2.8): find

(U, Ptot,P,M,M0) ∈ Y1,st such that

(3.14)

EsU + D∗

s Ptot = Lu

in S0

t (U∗

s )

λDsU − Ptot − αPDs P = 0

in S0

t (Ds)

αPPs DsU + σP − M = 0

in S0

t (Ps)

dt(M,M0) + LsP = Lp

in S0

t (P∗

s )

M0 = L0

in P∗

s .

Note that, also in this case, there are some ambiguities between functions and

functionals, that can be clarified in the vein of Remark 2.6.

Remark 3.6 (Two- and four-field formulation). The weak formulation 2.8 involves

four unknown variables but it can be equivalently rewritten as a two-field weak

formulation of the initial-boundary value problem (2.1), by eliminating the total

pressure ptot and the total fluid content m, cf. [24, Remark 2.7]. Analogously,

we could eliminate the discrete total pressure Ptot and the discrete fluid content

M from (3.14). In this way, we would equivalently obtain a discretization of the

two-field weak formulation, with trial and test spaces given by S0

t (Us)×S0

t (Ps)×Ps.

12

C. KREUZER AND P. ZANOTTI

Since we aim at establishing the well-posedness of (3.14) via the inf-sup theory,

it is convenient viewing it as an instance of the following linear variational problem:

given L ∈ Y∗

2,st, find Y1 ∈ Y1,st such that

(3.15)

bst(Y1,Y2) = ⟨L, Y2⟩

Y2,st

∀Y2 ∈ Y2,st.

Of course, this can be seen as a discretization of (2.13). Here, the bilinear form

bst : Y1,st × Y2,st → R is defined by

bst(CY1,Y2) :=

∫ T

0

(

⟨Es CU + D∗

s CPtot,V ⟩Us + ⟨dt(b

M, b

M0) + Ls CP,N⟩Ps

)

+ ⟨b

M0,N0⟩Ps

+

∫ T

0

(

(λDs CU − CPtot − αPDs CP,Qtot)Ω + (αPPs Ds CU + σ CP − b

M,Q)Ω

)

(3.16)

for CY1 = (CU, CPtot, CP, b

M, b

M0) ∈ Y1,st and Y2 = (V,Qtot,Q,N,N0) ∈ Y2,st.

A remarkable difference between (2.13) and (3.15) is that, in the latter one, we

do not have to explicitly take the closure of the trial space. Indeed, Y1,st is certainly

closed, being finite-dimensional.

3.4. Well-posedness. The goal of this section is to establish the well-posedness

of the discretization (3.14) by means of the inf-sup theory. Since the trial space

Y1,st and the test space Y2,st are finite-dimensional with equal dimension, we can

use a simplified version of the Banach-Ne˘cas theorem, which does not require to

verify the nondegeneracy of the form bst, see e.g. [14, Theorem 26.6]. Therefore,

the well-posedness is equivalent to the properties verified in the next two lemmas.

Lemma 3.7 (Boundedness). The bilinear form bst in (3.16) satisfies

sup

Y2∈Y2,st

bst(CY1,Y2)

∥Y2∥2

≲ ∥CY1∥1,st

for all CY1 ∈ Y1,st. The hidden constant depends only on the constants in the as-

sumptions (H1) and (H2) in section 3.1.

Proof. The claimed bound follows from the Cauchy-Schwarz inequality applied to

each term in the definition of bst, in combination with (3.3) and with the norm

equivalences in the assumptions (H1) and (H2).

□

Lemma 3.8 (Inf-sup stability). The bilinear form bst in (3.16) satisfies

sup

Y2∈Y2,st

bst(CY1,Y2)

∥Y2∥2

≳

(1 + T)− 1

2

(

∥CY1∥2

1,st + ∥b

M∥2

L∞(P∗

s ) +

∫ T

0

(

λ∥Ds CU∥2

Ω + γ−1∥ CP∥2

Ω

)

\1

2

.

for all CY1 = (CU, CPtot, CP, b

M, b

M0) ∈ Y1,st. The hidden constant depends only on the

constants in the assumptions (H1) and (H2) in section 3.1.

Proof. See section 3.5.

□

The combination of the two above lemmas implies the main result of this section.

INF-SUP STABLE DISCRETIZATION OF THE BIOT’S EQUATIONS

13

Theorem 3.9 (Well-posedness of the full discretization). For all possible loads

(Lu,Lp,L0) ∈ S0

t (U∗

s ) × S0

t (P∗

s ) × P∗

s , the equations (3.14) have a unique solution

Y1 = (U, Ptot,P,M,M0) ∈ Y1,st, which satisfies the two-sided stability bound

∥Y1∥2

1,st ≂

∫ T

0

(

∥Lu∥2

U∗

s

+ ∥Lp∥2

P∗

s

)

+ ∥L0∥2

P∗

s

.

(3.17)

Moreover, we have the norm equivalence

(3.18)

∥Y1∥2

1,st ≂ ∥Y1∥2

1,st + ∥M∥2

L∞(P∗

s ) +

∫ T

0

(λ∥DsU∥2

Ω + γ−1∥P∥2

Ω) .

The hidden constants depend only on the final time T and on the constants in the

assumptions (H1) and (H2) in section 3.1.

Proof. The equations (3.14) are equivalent to the linear variational problem (3.15)

with load L = (Lu, 0, 0,Lp,L0) ∈ Y∗

2,st. Then, the simplified version of the Banach-

Ne˘cas theorem for finite-dimensional linear variational problems [14, Theorem 26.6]

implies existence and uniqueness of the solution, in combination with Lemmas 3.7

and 3.8. The identity bst(Y1, ·) = L in Y∗

2,st further implies

∥Y1∥2

1,st + ∥M∥2

L∞(P∗

s ) +

∫ T

0

(

λ∥DsU∥2

Ω + γ−1∥P∥2

Ω

)

≲ ∥(Lu, 0, 0,Lp,L0)∥2

2,∗ ≲ ∥Y1∥2

1,st

according to estimates in Lemmas 3.7 and 3.8. The hidden constants depend only on

the ones therein. Here ∥·∥2,∗ denotes the dual of the test norm. This readily implies

the second claimed equivalence. The first one follows by recalling the definition of

the norm ∥·∥2 in (2.14).

□

Remark 3.10 (Jumps of the discrete fluid content). Recall that the component M

of the solution Y1 in Theorem 3.9 represents the discretization of the total fluid

content m in Theorem 2.4. While the latter one is guaranteed to be continuous

in time, the former one is not, (indeed it is piecewise constant). Still, a careful

inspection at the proofs of Lemmas 3.4 and 3.8 reveals that the left-hand side in

the stability bound (3.17) controls the jumps of M in time, namely

J

∑

j=1

∥Mj − Mj−1∥2

P∗

s

≲ ∥Y1∥2

1,st.

This can be viewed as a discrete counterpart of the continuity of m.

3.5. Inf-sup stability. This section is devoted to the proof of Lemma 3.8. The

argument is similar to the one in [24, section 3.1], which establishes the inf-sup

stability of the form b in (2.15). Nevertheless, some subtleties exist with regard

to the discretization and we therefore include a detailed proof so as to keep the

discussion as much complete and self-contained as possible.

Let CY1 = (CU, CPtot, CP, b

M, b

M0) ∈ Y1,st be given and set

CHDs := λDs CU − CPtot − αPDs CP

and

CHPs := αPPs Ds CU + σ CP − b

M

14

C. KREUZER AND P. ZANOTTI

for shortness. Since the projections PDs and PPs map onto Ds and Ps, respectively,

we have CHDs ∈ S0

t (Ds) as well as CHPs ∈ S0

t (Ps). We consider the test function

Y2,j :=

(

(CU + E−1

s

D∗

s CPtot)χj,

4 max{1,C}

µ + λ

CHDs χj,

4γ

min{1,c}

CHPs χj,

L−1

s (2b

M + dt(b

M, b

M0) + Ls CP)χj, L−1

s

b

M0

)

with j = 1,...,J, where χj : [0,T] → R denotes the indicator function of the inter-

val [0,tj]. Here the constants c and C are as in the assumption (H2) in section 3.1.

We refer to [24, Remark 3.9] for a motivation of the test function.

First of all, notice that we indeed have that Y2,j ∈ Y2,st, i.e. it is an admissible

test function. Indeed, recall the definition of the spaces Y1,st and Y2,st in (3.11)

and (3.13). The operator E−1

s

D∗

s maps Ds into Us and Ls is an isometry between Ps

and P∗

s . Moreover, the indicator function χj is piecewise constant on the partition

(3.8) of [0,T]. Hence, the multiplication by χj preserves the inclusion in Y2,st.

The proof consists of two main steps. First, we establish the lower bound

bst(CY1,Y2,j + Y2,J ) ≳∥CY1∥2

1,st + ∥b

M∥2

L∞(P∗

s ) +

∫ T

0

(

λ∥Ds CU∥2

Ω + γ−1∥ CP∥2

Ω

)

(3.19)

for a suitable index 1 ≤ j ≤ J. Then, we estimate the norm of the test function

(3.20)

∥Y2,j∥2

2 ≲ ∥CY1∥2

1,st + (1 + T)∥b

M∥2

L∞(P∗

s )

for some hidden constant independent of j. The combination of these inequalities

readily implies the inf-sup stability claimed in Lemma 3.8.

We start with the proof of (3.19). Using the test function Y2,j in the definition

(3.16) of the form bst yields

(3.21)

bst(CY1,Y2,j) =

∫ tj

0

⟨Es CU + D∗

s CPtot, CU + E−1

s

D∗

s CPtot⟩Us

(=: I1)

+ 2

∫ tj

0

⟨dt(b

M, b

M0) + Ls CP, L−1

s

b

M⟩Ps

(=: I2)

+ ∥b

M0∥2

P∗

s

+

∫ tj

0

∥dt(b

M, b

M0) + Ls CP∥2

P∗

s

+

∫ tj

0

4 max{1,C}

µ + λ

∥ CHDs ∥2

Ω

+

∫ tj

0

4γ

min{1,c}

∥ CHPs ∥2

Ω.

Notice that the dual norms on the third line are obtained thanks to the second part

of (3.4). We are led to analyze the terms I1 and I2 on the right-hand side.

Regarding the first term, we have

I1 =

∫ tj

0

(

⟨Es CU, CU⟩Us + 2⟨D∗

s CPtot, CU⟩Us + ⟨D∗

s CP, E−1

s

D∗

s CPtot⟩Us

)

.

We rewrite the first summand with the help of the first identity in (3.3). Analo-

gously, we estimate the third summand from below by means of the first identity

in (3.4) and the assumption (H2) in section 3.1

⟨Es CU, CU⟩Us = ∥CU∥2

Us

and ⟨D∗

s CPtot, E−1

s

D∗

s CPtot⟩Us ≥

c

µ

∥ CPtot∥2

Ω.

INF-SUP STABLE DISCRETIZATION OF THE BIOT’S EQUATIONS

15

We estimate the remaining term by observing that the assumption (H2) implies

also µ∥Ds · ∥2

Ω ≤ C∥·∥2

Usin Us. This bound and a Young’s inequality imply

⟨D∗

s CPtot, CU⟩Us = −(Ds CU, CHDs )Ω + λ∥Ds CU∥2

Ω − α(Ds CU, PDs CP)Ω

≥

3λ

4

∥Ds CU∥2

Ω −

1

4

∥CU∥2

Us −

max{1,C}

µ + λ

∥ CHDs ∥2

Ω − α(Ds CU, CP)Ω.

Notice that we were able to omit the projection PDs in the last term thanks to the

inclusion Ds CU ∈ S0

t (Ds). Combining this bound with the identities above reveals

J1 ≥

∫ tj

0

(1

2

∥CU∥2

Us +

c

µ

∥ CPtot∥2

Ω +

3λ

2

∥Ds CU∥2

Ω

−

2 max{1,C}

µ + λ

∥ CHDs ∥2

Ω − 2α(Ds CU, CP)Ω

)

.

Regarding the other critical term in (3.21), we have

I2 = 2

∫ tj

0

(

⟨dt(b

M, b

M0), L−1

s

b

M⟩Ps + (b

M, CP)Ω

)

.

The use of the L2(Ω)-scalar product in the second summand is justified by the

discussion after (3.7). We estimate the first summand from below by invoking

Lemma 3.4

2

∫ tj

0

⟨dt(b

M, b

M0), L−1

s

b

M⟩Ps ≥ ∥b

Mj∥2

P∗

s − ∥b

M0∥2

P∗

s

.

To investigate the second summand, assume first σ > 0. A Young’s inequality

reveals

(b

M, CP)Ω = ( CHPs , CP)Ω + α(PPs Ds CU, CP)Ω + σ∥ CP∥2

Ω

≥ α(Ds CU, CP)Ω +

σ

2

∥ CP∥2

Ω −

1

2σ

∥ CHPs ∥2

Ω.

As before, we omit the projection PPs thanks to the inclusion CP ∈ S0

t (Ps). Alterna-

tively, for general σ ≥ 0, it holds that

(b

M, CP)Ω = ( CHPs , CP)Ω + α(PPs Ds CU, CP)Ω + σ∥ CP∥2

Ω

= −

1

α

( CHPs , CHDs )Ω +

λ

α

( CHPs , Ds CU)Ω −

1

α

( CHPs , CPtot)Ω

+ ( CHPs , CP − PDs CP)Ω + α(Ds CU, CP)Ω + σ∥ CP∥2

Ω.

Bounding the term ( CHPs , CP − PDs CP)Ω conveniently is subtle. Whenever Ps ⊆ Ds,

we have PDs CP = CP, hence ( CHPs , CP − PDs CP)Ω = 0. When the above inclusion fails,

we have σ > 0 due to the assumption (H4) in section 3.1. Then, we apply Young’s

inequality to obtain ( CHPs , CP − PDs CP)Ω ≤ ∥ CHPs ∥2

Ω/(2σ) + σ∥ CP∥2

Ω/2. Combining this

observation with other applications of Young’s inequality, the previous lower bound

and the definition (2.12) of γ, we arrive at

(b

M, CP)Ω ≥ −

max{1,C}

2(µ + λ)

∥ CHDs ∥2

Ω −

c

4µ

∥ CPtot∥2

Ω −

λ

2

∥Ds CU∥2

Ω

+

σ

2

∥ CP∥2

Ω −

3γ

2 min{1,c}

∥ CHPs ∥2

Ω + α(Ds CU, CP)Ω.

16

C. KREUZER AND P. ZANOTTI

By combining the last two bounds with the above identity for I2, we obtain

I2 ≥ ∥b

Mj∥2

P∗

s − ∥b

M0∥2

P∗

s

+

∫ tj

0

(

−

max{1,C}

µ + λ

∥ CHDs ∥2

Ω −

c

2µ

∥ CPtot∥2

Ω −

λ

2

∥Ds CU∥2

Ω + σ∥ CP∥2

Ω

−

3γ

min{1,c}

∥ CHPs ∥2

Ω + 2α(Ds CU, CP)Ω

)

.

After this preparation, we are in position to choose a specific test function in

order to establish (3.19). We let j ∈ {1,...,J} be such that ∥b

Mj∥P∗

s

= ∥b

M∥L∞(P∗

s ).

Then, we insert the previous lower bounds of I1 and I2 into (3.21), resulting in

bst(CY1,Y2,j + Y2,J ) ≥

∫ T

0

(1

2

∥CU∥2

Us +

c

2µ

∥ CPtot∥2

Ω + ∥dt(b

M, b

M0) + Ls CP∥2

P∗

s

)

+

∫ T

0

( 1

µ + λ

∥ CHDs ∥2

Ω + γ∥ CHPs ∥2

Ω +

λ

2

∥Ds CU∥2

Ω + σ∥ CP∥2

Ω

)

+ ∥b

M∥2

L∞(P∗

s ).

In combination with the definition (3.12) of the norm ∥·∥1,st and the assumption

(H1) in section 3.1, this almost establishes (3.19). Indeed, we have

bst(CY1,Y2,j + Y2,J ) ≳∥CY1∥2

1,st + ∥b

M∥2

L∞(P∗

s ) +

∫ T

0

(

λ∥Ds CU∥2

Ω + σ∥ CP∥2

Ω

)

.

The proof that we can indeed replace σ by γ−1 in the last summand follows verbatim

the argument in the proof of [24, Proposition 4.1].

The last step of the proof consists in establishing (3.20). According to the

definition (2.14) of the test norm ∥·∥2, we have for all j = 1,...,J, that

∥Y2,j∥2

2 =

∫ tj

0

(

∥CU + E−1

s

D∗

s CPtot∥2

U

+ ∥L−1

s (2b

M + dt(b

M, b

M0) + Ls CP)∥2

P

)

+

∫ tj

0

(16 max{1,C2}

µ + λ

∥ CHDs ∥2

Ω +

16γ

min{1,c2}

∥ CHPs ∥2

Ω

)

+ ∥L−1

s

b

M0∥2

P

.

We exploit the identities (3.3) and (3.4) and the norm equivalences in the as-

sumptions (H1) and (H2) from section 3.1. We extend also the integrals from the

interval [0,tj] to [0,T]. Hence, we obtain

∥Y2,j∥2

2 ≲

∫ T

0

(

∥CU∥2

U

+

1

µ

∥ CPtot∥2

Ω + ∥b

M∥2

P∗

s

+ ∥dt(b

M, b

M0) + Ls CP∥2

P∗

s

)

+

∫ T

0

( 1

µ + λ

∥ CHDs ∥2

Ω + γ∥ CHPs ∥2

Ω

)

+ ∥b

M0∥2

P∗

s

.

Finally, we recall the definition (3.12) of the norm ∥·∥1,st and exploit the upper

bound /

T

0

∥·∥2

P∗

s

≤ T∥·∥2

L∞(P∗

s ). This establishes (3.20) for some hidden constant

independent of j and concludes the proof.

Remark 3.11 (Time-dependent space discretization). In our framework the space

discretization does not change in time and this is important for the proof of the

inf-sup stability in this section. Indeed, a change in the space discretisation would

imply the change of the operator Ls in time and this would have an effect on our

use of the integration by parts formula from Lemma 3.4. For the same reason,

INF-SUP STABLE DISCRETIZATION OF THE BIOT’S EQUATIONS

17

we argued in [24, Remark 2.1] that the parameter κ in the Biot’s equations (2.1a)

could be allowed to vary in space but not in time.

4. Discretisation with Lagrange elements in space

In this section we propose an exemplary concrete space discretization, based on

H1-conforming Lagrange finite elements for all variables. Hence, we first verify the

assumptions (H1)-(H4) from section 3.1, in order to infer the well-posedness. Then,

we discuss the a priori error analysis.

Our notation and assumptions for the finite elements are as follows. We denote

by T a face-to-face simplicial mesh of Ω. The shape constant of T is given by

(4.1)

max

T∈T

hT

rT

where hT is the diameter of a d-simplex T ∈ T and rT is the diameter of the largest

ball inscribed in T.

We denote by F the set of all faces of T and by Fi the interior faces. We assume

that the mesh is compatible with the boundary conditions (2.1c). This means that

each face F ∈ F \ Fi (i.e. each boundary face) satisfies

either F ⊆ Γu,N or F ⊆ Γu,E

and

either F ⊆ Γp,N or F ⊆ Γp,E.

Hence, we have two partitions of the boundary faces Fu,E ∪ Fu,N and Fp,E ∪ Fp,N ,

where each set F∗ is defined as

(4.2)

F∗ := {F ∈ F | F ⊆ Γ∗}.

We let the meshsize h and the normal n be the piecewise constant functions on

T and on the skeleton of T (i.e. the union of all faces) defined as

h|T := hT

and

n|F := nF

for all T ∈ T and F ∈ F, respectively. Here nF is a fixed unit normal vector of F,

pointing outside Ω if F is a boundary face.

The space Pk(S), k ≥ 0, consist of all polynomials of total degree ≤ k on an

n-simplex S ⊆ Rd with 1 ≤ n ≤ d. The corresponding space of (possibly discontin-

uous) piecewise polynomials on the mesh T is denoted by Pk(T).

4.1. Concrete space discretisation. A concrete realization of the abstract space

discretization from section 3.1 requires a specific choice of the spaces Us and Ps and

of the operators Es and Ls in (3.1) and (3.2), respectively. We have to prescribe

also the operator Ds in (3.5) and to characterize its range Ds in (3.6). Finally, we

need to verify the assumptions (H1)-(H4).

For the discretization of the displacement, we consider H1-conforming Lagrange

finite elements of degree k + 1 with k ≥ 1, i.e., we set

(4.3)

Us := Pk+1(T)d ∩ U.

The intersection with the space U from (2.3) enforces the global continuity (hence

the H1-conformity) as well as the boundary conditions. As usual in conforming

finite elements, we discretize the operator E in (2.2) by Es : Us → U∗

s defined as

(4.4)

⟨Es CU,V ⟩Us := ⟨E CU,V ⟩U

for all CU,V ∈ Us.

18

C. KREUZER AND P. ZANOTTI

Recall that the operator Ds should be a discretization of the divergence and that

the pair Us/Ds with Ds = Ds(Us) should enjoy a discrete Stokes inf-sup condition

in order to satisfy assumption (H2). Given the above definition of Us, one option

for Ds is suggested by the so-called Hood-Taylor pair, cf. [14, section 54.4]. Hence,

we consider H1-conforming Lagrange finite elements of degree k

Ds := Pk(T) ∩ H1(Ω) ∩ D

(4.5)

and we define Ds : Us → Ds by

(4.6)

(Ds CU,Qtot)Ω = (D CU,Qtot)Ω

for all CU ∈ Us and Qtot ∈ Ds. According to (2.5), the intersection with D in (4.5)

simply enforces the vanishing mean value when Γu,E = ∂Ω. Note also that Ds CU is

nothing else than the H-orthogonal projection of D CU onto Ds.

Finally, the assumption (H4) suggests that also the space Ps, for the discretiza-

tion of the pressure and of the total fluid content, should consist of H1-conforming

Lagrange finite elements of degree k. Therefore, we set

(4.7)

Ps := Pk(T) ∩ P.

The intersection with the space P from (2.4) enforces global continuity, the bound-

ary conditions and, possibly, the vanishing mean value. In this case, we define

Ls : Ps → P∗

s in terms of L in (2.2) by

(4.8)

⟨Ls CP,Q⟩Ps := ⟨L CP,Q⟩P

for all CP,Q ∈ Ps.

Remark 4.1 (Hidden constants). In order to simplify the statement of the next

results, it is implicitly understood hereafter that all hidden constants in our esti-

mates potentially depend on the final time T, the shape constant (4.1) of T and the

polynomial degree k. In particular, the latter is arbitrary but fixed in our setting.

The possible dependence on other relevant quantities is addressed case by case.

Having introduced all the spaces and the operators required for the space dis-

cretization in section 3.1, we can verify the validity of the assumptions (H1)-(H4).

Proposition 4.2 (Verification of the assumptions). Let the spaces Us, Ds and Ps

and the operators Es, Ds and Ls be defined by (4.3)-(4.8).

(1) The assumption (H1) holds true and the equivalences therein are actually

identities.

(2) The assumption (H2) holds true and the constants therein depend only on

the quantities mentioned in Remark 4.1, provided that each d-simplex in T

has at least one vertex in the interior of Ω.

(3) The assumption (H3) holds true and the constants therein depend only on

the quantities mentioned in Remark 4.1, provided that the grading of T,

defined as in [12], is strictly less than (

√

2k + d +

√

k)/(

√

2k + d −

√

k).

(4) The assumption (H4) holds true.

Proof. (1) Owing to the inclusions Us ⊆ U and Ps ⊆ P, there is no need to extend

the norms ∥·∥U and ∥·∥P. Moreover, the combination of (3.3) with (4.4) and (4.8)

readily implies ∥·∥Us = ∥·∥U in Us as well as ∥·∥Ps = ∥·∥P in Ps.

INF-SUP STABLE DISCRETIZATION OF THE BIOT’S EQUATIONS

19

(2) The upper bound in (H2) follows from the boundedness of Ds, which satisfies

∥Ds ·∥2

Ω ≤∥D·∥2

Ω ≲ µ−1∥·∥2

Uin Us, according to (2.2), (2.11) and (4.6). The lower

bound can be rephrased as

c∥·∥2

Ω ≤ µ

(

sup

V ∈Us

(DV, ·)Ω

∥V ∥U

)2

in Ds. The definitions (2.2) and (2.11) reveal that this is equivalent to the discrete

Stokes inf-sup stability of the pair Us/Ds, i.e. of the Hood-Taylor pair. The latter

condition is known to hold true under the above assumption on the mesh; see [14,

sections 54.3 and 54.4].

(3) We have

∥P∗

Ps · ∥P∗ = sup

w∈P

⟨·, PPs w⟩

Ps

∥w∥P

= sup

w∈P

(·, PPs w)Ω

∥w∥P

in P∗

s . The inclusion Ps ⊆ P implies the lower bound in (H3) with c = 1, because

PPs

is a projection onto Ps. The upper bound is equivalent to the P-stability of

PPs ; cf. [41]. Owing to (2.11) and (4.7), this is the H1(Ω)-stability of the L2(Ω)-

orthogonal projection onto H1-conforming Lagrange finite element spaces. Such a

condition is known to hold under the asserted grading assumption thanks to [11,

Theorem 4.14(ii)] together with [11, Theorem 4.4] and [11, Remark 4.4].

(4) Let σ = 0. The combination of (4.5) and (4.7) with (2.4) and (2.5) implies,

for Γu,E = ∂Ω that

Ps = Pk(T) ∩ H1

Γp,E (Ω) ∩ L2

0(Ω) and Ds = Pk(T) ∩ H1(Ω) ∩ L2

0(Ω).

For Γu,E ̸= ∂Ω we have instead,

Ps ⊆ Pk(T) ∩ H1

Γp,E (Ω) and Ds = Pk(T) ∩ H1(Ω).

In both cases we have Ps ⊆ Ds, i.e., assumption (H4) is verified.

□

Remark 4.3 (Assumptions for (H3)). We refer to [12, Definitions 1.1] for the exact

definition of mesh grading. Clearly the grading of quasi-uniform meshes equals

1 and thus satisfies the grading condition in Proposition 4.2(3). It follows from

[12, Theorem 1.3], that the grading condition is also satisfied for all polynomial

degrees k ≥ 1 and dimensions d ≤ 6 provided T is obtained by adaptive bisection

of a colored initial mesh; see [12, Assumption 3.1] for the notion of colored mesh.

Note that assumption (H3) can be verified under different assumptions. Indeed,

the H1(Ω)-stability of the L2(Ω)-orthogonal projection onto finite element spaces

has been extensively analyzed by various authors; we refer [11] for an overview of

the existing results.

Having verified all assumptions in section 3.1, we deduce the well-posedness of

the discretization (3.14) with the spaces and the operators proposed in this section.

In particular, the inclusions Us ⊆ U and Ps ⊆ P suggest to define the data in the

discretization by restriction of the data in the weak formulation (2.8)

(4.9)

Lu = ℓ

u|S0

t (Us),

Lp = ℓ

p|S0

t (Ps)

and

L0 = ℓ0|Ps .

Theorem 4.4 (Well-posedness with Lagrange elements). Let the spaces Us, Ds and

Ps and the operators Es, Ds and Ls be defined by (4.3)-(4.8). Assume that T is ob-

tained from a colored initial mesh by newest vertex bisection and that each d-simplex

20

C. KREUZER AND P. ZANOTTI

in T has at least one vertex in the interior of Ω. Then, the discretization (3.14)

with the data (4.9) has a unique solution Y1 = (U, Ptot,P,M,M0) ∈ Y1,st with

∥Y1∥2

1,st ≂ ∥Y1∥2

1,st + ∥M∥2

L∞(P∗

s ) +

∫ T

0

(λ∥DsU∥2

Ω + γ−1∥P∥2

Ω

)

≂

∫ T

0

(

∥Lu∥2

U∗

s

+ ∥Lp∥2

P∗

s

)

+ ∥L0∥2

P∗

s

≲

∫ T

0

(∥ℓu∥2

U∗ + ∥ℓp∥2

P∗) + ∥ℓ0∥2

P∗ .

(4.10)

All hidden constants depend only on the quantities mentioned in Remark 4.1.

Proof. The combination of Theorem 3.9 with Proposition 4.2 implies the existence

and the uniqueness of the solution and guarantees that the two equivalences in

(4.10) hold true. Then, the upper bound follows from the discretization (4.9) of the

data and the inclusions Us ⊆ U and Ps ⊆ P.

□

4.2. Quasi-optimality. In order to quantify the accuracy of the proposed dis-

cretization, we need an error notion that is related to both the trial norm ∥·∥1 in

(2.10) and to the discrete trial norm ∥·∥1,st in (3.12), cf. (4.12) below. A major

issue is that the former one involves the dual norm in P, whereas the latter one

involves the dual norm in Ps.

In order to compare functionals defined on the two spaces, we can use the adjoint

P∗ : P∗

s → P∗ of a bounded projection P : P → Ps. The specific choice of P is critical

because of the term ∂t

Cm + LCp in the norm ∥·∥1. Since this term acts in space

via the pairing ⟨·, ·⟩

P

, one may want to employ the L2(Ω)-orthogonal projection,

because the pivot space in the Hilbert triplet (2.7) is equipped with the L2(Ω)-norm.

Still, according to (2.11), the action of L induces the P-norm, thus suggesting to

make use of the P-orthogonal projection. Of course, similar comments apply to the

counterpart of ∂t

Cm + LCp in ∥·∥1,st.

Both options have advantages and disadvantages. The latter one, being related

to the P-norm, suggests the use of piecewise polynomials of degree k + 1 for the

discretization of P, but this would be critical for the assumption (H4), as Ds consists

of piecewise polynomials of degree k, cf. (4.5)-(4.7). The former option does not

have this issue, therefore we go for it. Still, the derivation of higher-order decay

rates in space appears critical in this case, cf. Remark 4.12.

In accordance with the above discussion and with the inclusions Us ⊆ U and

Ps ⊆ P, we consider the error notion Err : Y1 × Y1,st → [0, +∞) defined as

ERR(Cy1, CY1)2 :=

∫ T

0

(

∥u − CU∥2

U

+

1

µ

∥ptot − CPtot∥2

Ω

)

+

∫ T

0

∥∂t

Cm + LCp− P

∗

Ps (dt(b

M, b

M0) + Ls CP)∥2

P∗

+ ∥ Cm(0) − P∗

Ps b

M0∥2

P∗

+

∫ T

0

1

λ + µ

∥λDCu − Cptot − αPD Cp− (λDs CU − CPtot − αPDs CP)∥2

Ω

+

∫ T

0

γ∥αP

PDCu + σCp− Cm − (αPPs Ds CU + σ CP − b

M)∥2

Ω

(4.11)

INF-SUP STABLE DISCRETIZATION OF THE BIOT’S EQUATIONS

21

for Cy1 = (Cu, Cptot, Cp, Cm) ∈ Y1 and CY1 = (CU, CPtot, CP, b

M, b

M0) ∈ Y1,st. Notice that this

is not a norm on the sum Y1 + Y1,st. For instance, in general, we have D CU ̸= Ds CU

for CU ∈ S0

t (Us) ⊆ L2(U). Still, the above error notion measures the accuracy in the

approximation of all functions and functionals involved in the trial norm ∥·∥1 and

it holds that

(4.12)

ERR(Cy1, 0) = ∥Cy1∥1

and

ERR(0, CY1) ≂ ∥CY1∥1,st.

Indeed, the second equivalence follows from (1) and (3) in Proposition 4.2.

The equations (3.14) with the spaces and the operators defined by (4.3)-(4.8) are

not a conforming Petrov-Galerkin discretization of the weak formulation (2.8). In-

deed, the trial space Y1,st in the former problem is not a subspace of its counterpart

Y1 in the latter one. Nevertheless, we have for the test spaces the inclusion

Y2,st ⊆ Y2.

This property and the definition (4.9) of the load in (3.14) by restriction of the one in

(2.8) ensure that we can still guarantee the fundamental quasi-optimality property

of inf-sup stable conforming Petrov-Galerkin methods by a standard argument; cf.

[2].

Theorem 4.5 (Quasi-optimality). Let all assumptions in Theorem 4.4 be verified.

Denote by y1 ∈ Y1 and Y1 ∈ Y1,st, respectively, the solutions of (2.8) and (3.14)

with (4.9). Then, we have

ERR(y1,Y1) ≲ inf

˜Y1∈Y1,st

ERR(y1, CY1).

(4.13)

The hidden constant depends only on the quantities mentioned in Remark 4.1.

Proof. For CY1 ∈ Y1,st, the triangle inequality and (4.12) imply

ERR(y1,Y1) ≤ ERR(y1, CY1) + ERR(0,Y1 − CY1) ≲ ERR(y1, CY1) + ∥Y1 − CY1∥1,st.

According to the inf-sup stability in Lemma 3.8, we have

∥Y1 − CY1∥1,st ≲ sup

Y2∈Y2,h

bst(Y1 − CY1,Y2)

∥Y2∥2

= sup

Y2∈Yh

2

b(y1,Y2) − bst(CY1,Y2)

∥Y2∥2

,

where the identity follows by comparing problems (2.8) and (3.14). We exploit the

definitions (2.15), (3.16), and (4.4)-(4.8), as well as the invariance of the projection

PPs onto Ps. For y1 = (u, ptot, p, m) and CY1 = (CU, CPtot, CP, b

M, b

M0), this results in

b(y1,Y2) − bst(CY1,Y2) :=

∫ T

0

((

E(u − CU),V 〉

U

+ (ptot − CPtot, DV )

Ω

)

+

∫ T

0

(∂tm + L CP − P∗

Ps (dt(b

M, b

M0) + Ls CP),N

〉

P

+ (m − P∗

Ps b

M0,N0

〉

P

+

∫ T

0

(λDu − ptot − αPDp − (λDs CU − CPtot − αPDs CP),Qtot

)

Ω

+

∫ T

0

(αP

P

Du + σp − m − (αPPs Ds CU + σ CP − b

M),Q)

Ω

.

(4.14)

22

C. KREUZER AND P. ZANOTTI

Then Cauchy-Schwarz inequalities, the bound ∥DV ∥2

Ω ≲ µ−1∥V ∥2

U

and the defini-

tion (2.14) of the test norm ∥·∥2 yield

∥Y1 − CY1∥1,st ≲ ERR(y1, CY1).

Inserting this estimate into the first inequality above concludes the proof.

□

Remark 4.6 (Augmented error notion). Our definition of the error notion is in a

sense minimal, because we have included only the terms that arise by applying the

Cauchy-Schwarz inequality in (4.14). Still, the first equivalence in (4.10) clarifies

that that the statement and the proof of Theorem 4.5 will remain unchanged if we

augment ERR(·, ·) by adding any of the following terms

∥ Cm − P∗

Ps b

M∥2

L∞(P∗),

∫ T

0

λ∥DCu − Ds CU∥2

Ω,

and

∫ T

0

γ−1∥Cp− CP∥2

Ω.

4.3. A priori error estimates. The quasi-optimality stated above has two re-

markable properties and (at least) one clear disadvantage. One the one hand, the

estimate (4.13) holds true for any solution of the weak formulation (2.8) (i.e. no ad-

ditional regularity is required) and the hidden constant is robust with respect to all

material parameters. On the other hand, it is not immediate how (and even if) the

error decays to zero, because the definition (4.11) of ERR(·, ·) involves a nontrivial

coupling of the various components of the solution.

In this section, we investigate the latter aspect. Roughly speaking, we aim at

showing that the best error in the right-hand side of (4.13) is equivalent to a sum

of best errors. In order to preserve the two above-mentioned properties, we make

sure also that the constants in the equivalence are independent of the material

parameters and that no additional regularity of the solution of (2.8) is required

beyond the minimal one, namely y1 ∈ Y1. When such a result is available, the

decay of the error to zero can be easily discussed in terms of classical results from

approximation theory.

The precise statement of our main result is as follows. Notice that the last term

on the right-hand side of (4.15) does not appear in our definition of the error notion,

but we could equivalently include it according to Remark 4.6.

Theorem 4.7 (Best error decoupling). Let all assumptions in Theorem 4.4 be

verified and denote by y1 = (u, ptot, p, m) ∈ Y1 the solution of (2.8). Then we have

that

inf

˜Y1∈Y1,st

ERR(y1, CY1)2 ≲

inf

̂U∈S0

t (Us)

∫ T

0

∥u − CU∥2

U

+

inf

̂Ptot∈S0

t (D)

∫ T

0

1

µ

∥ptot − CPtot∥2

Ω

+

inf

̂

W ∈S0

t (Ps)

∫ T

0

∥∂tm + Lp − P∗

Ps ̂

W∥2

P∗ + inf

̂

M0∈Ps

∥m(0) − P∗

Ps ̂

M0∥2

P∗

+ ε−1

s

inf

̂P∈S0

t (Ps)

∫ T

0

1

γ

∥p − CP∥2

Ω.

(4.15)

The hidden constant depends on the quantities mentioned in Remark 4.1 and on

the constant in (4.20) and εs is defined in (4.24) below.

Proof. See sections 4.4-4.6.

□

INF-SUP STABLE DISCRETIZATION OF THE BIOT’S EQUATIONS

23

A first consequence of Theorem 4.7 is plain convergence: the error converges to

zero as the mesh-size converges to zero, both in space and time. This holds true

irrespective of the regularity of the solution. Moreover, first-order convergence can

be established under additional regularity assumptions.

Corollary 4.8 (First-order convergence). Let all assumptions in Theorem 4.4 be

verified. Denote by y1 = (u, ptot, p, m) ∈ Y1 and Y1 ∈ Y1,st, respectively, the

solutions of (2.8) and (3.14) with (4.9). Assume additionally

(4.16)

u ∈ H1(U) ∩ L2(H2(Ω)d)

ptot ∈ H1(D) ∩ L2(H1(Ω))

p ∈ H1(P) ∩ L2(H1(Ω))

(∂tm + Lp) ∈ H1(P∗) ∩ L2(L2(Ω))

m(0) ∈ L2(Ω).

Then, the error can be bounded from above as follows

ERR(y1,Y1)2 ≲

(

max

Ω

h

)2

(

1

κ

∥m(0)∥2

Ω

+

∫ T

0

(

µ∥∇2u∥2

Ω +

1

µ

∥∇ptot∥2

Ω +

1

εsγ

∥∇p∥2

Ω +

1

κ

∥∂tm + Lp∥2

Ω

)}

+

(

max

j=1,...,J

|Ij|

)2 ∫ T

0

(

∥∂tu∥2

U

+

1

µ

∥∂tptot∥2

Ω +

1

εsγ

∥∂tp∥2

Ω + ∥∂t(∂tm + Lp)∥2

P∗

)

.

The hidden constant depends on the quantities mentioned in Remark 4.1 and on

the constant in (4.20) and εs is defined in (4.24) below.

Proof. Combine Theorem 4.5 with Theorem 4.7. Then, use standard Bramble-

Hilbert-like estimates (see, e.g., [40, Chapter 7]) in order to bound each term on

the right-hand side of (4.15).

□

Remark 4.9 (Regularity in space). According to [24, Theorem 5.4], the solution y1

of (2.8) satisfies the space regularity in (4.16) at least in some circumstances. To

be more precise, we have

u ∈ L2(H2(Ω)2)

ptot ∈ L2(H1(Ω))

p ∈ L2(H1(Ω))

(∂tm + Lp) ∈ L2(L2(Ω))

when the data are such that

ℓu ∈ L2(L2(Ω)2),

ℓp ∈ L2(L2(Ω)),

ℓ0 ∈ L2(Ω),

upon additionally assuming that Ω ⊆ R2 is convex, the boundary conditions (2.1c)

are posed on Γu,E = ∂Ω=Γp,N , and λ ≳ µ.

Remark 4.10 (Regularity in time). According to [31], the time regularity in (4.16)

may fail to hold even for smooth data. For a rough explanation, assume the data are

smooth in time, the solution enjoys the time regularity in (4.16) and, in addition,

we have m ∈ C1(P∗), i.e. we have one more time derivative than in Theorem 2.4.

Then, on the one-hand, the fourth equation in (2.8) implies p ∈ C0(P). On the

other hand, the evaluation of the other equations at t = 0 implies that the pair

(u(0),p(0)) solves the Stokes-like problem

(E + λD∗D)u(0) + αD∗PDp(0) = fu(0) in U∗

αP

P

divu(0) + σp(0) = ℓ0

in P.

24

C. KREUZER AND P. ZANOTTI

In general, the solution of this problem is in U×P. Hence the condition p(0) ∈ P can

hold true only for compatible data, otherwise some singularity must be expected

at the initial time. We refer to the Terzaghi test case discussed in section 5 below

for a numerical illustration.

Remark 4.11 (Grading of the time partition). The constant in Theorem 4.7 depends

on the one in (4.20), which measures the grading of the partition employed for the

discretization in time. The latter constants enters into play because Theorem 4.7

does not assume additional regularity in time of the solution. This calls for a Scott-

Zhang-like interpolation in time, see section 4.4 below for the details. When all

components of the solution are continuous in time, a Lagrange-like interpolation

is possible and the constant in (4.20) does not enter into the error estimation, cf.

[40, Theorem 4.5]. Still, according to Remark 4.10, the continuity at t = 0 is not

obvious.

Remark 4.12 (Decay rate). The space discretization considered in this section is

actually of higher-order, because we make use of H1-conforming Lagrange finite

elements of degree k + 1 for the displacement and of degree k for the other com-

ponents of the solution with k ≥ 1, cf. (4.3)-(4.7). Therefore, the question arises

if a higher-order decay rate with respect to h can be obtained under higher regu-

larity assumptions on the solution. On the one hand, the P∗-norm errors on the

right-hand side of (4.15) appear to be critical in this respect, because the space

Ps, possibly incorporating boundary conditions, is used to approximate functional

from P∗, which have no prescribed boundary conditions. On the other hand, it is

unclear to us if the above-mentioned regularity of the solution can be expected in

general. Indeed, the regularity result mentioned in Remark 4.9, assumes Γp,E = ∅.

Due to the subtlety of this matter, we postpone any further investigation on this

point to future work.

The remaining part of this section is devoted to the proof of Theorem 4.7. For

this purpose, we aim at constructing a bounded interpolant I : Y1 → Y1,st. The

operator I cannot be obtained by just approximating each component of the solu-

tion irrespective of the others, because of the nontrivial coupling of the components

in the error notion (4.11). Thus, to make sure that I is bounded and robust with

respect to the material parameters, we must guarantee that it is compatible with

the coupling. Roughly speaking this means that, when the error notion involves

some combination of the components of the solution, it must be possible to ac-

curately approximate it by the corresponding combination of the components of

the interpolant. We achieve this goal with the help of a number of commutative

diagrams, cf. Figures 4.1 and 4.2.

We divide our construction into three parts. In Section 4.4 we first discuss the

interpolation in time. In Section 4.5 we discuss the interpolation in space. Finally,

in Section 4.6, we combine the interpolation in time and space in order to prove

Theorem 4.7.

4.4. Time interpolation. The error notion ERR(·, ·) involves several terms with

L2 regularity in time, so we initially consider the problem of approximating function

from L2(X) into S0

t (X) for some Hilber space X. Since we cannot control the point

values of a function in L2(X), we cannot use Lagrange interpolation; cf. remark 4.11.

Thus, we resort to Scott-Zhang-like [37] interpolation in the vein of [40, Section 4].

INF-SUP STABLE DISCRETIZATION OF THE BIOT’S EQUATIONS

25

Recall (3.8)-(3.9). For j ∈ {1,...,J}, the polynomial ψj ∈ P1(Ij) defined as

(4.17)

ψj(t) := 6

t − tj−1

|Ij|2

−

2

|Ij|

is such that,

(4.18)

∫

Ij

Qψj = Q(tj)

∀Q ∈ P1(Ij).

We define J : L2(X) → S0

t (X) by

(J x)|Ij :=

∫

Ij

xψj

for x ∈ L2(X) and j = 1,...J.

Since the error notion in (4.11) involves the discrete time derivative dt, it is

useful to introduce another interpolant CJ : L2(X) → S0

t (X), defined as

( CJx)|Ij :=

∫

Ij

(∫ t

tj−1

x(s)ds

\

ψj(t)dt +

∫

Ij−1

(∫ tj−1

t

x(s)ds

)

ψj−1(t)dt

|Ij|

for j = 2,...J and

( CJx)|I1 :=

∫

I1

(∫ t

0

x(s)ds

)

ψ1(t)

|I1|

.

The first part of Lemma 4.13 below ensures that both J x and CJx are near best

approximations of x in the L2(X)-norm. The second part additionally clarifies that

the two interpolants are related, through the weak and the discrete time derivative,

as summarized by the commuting diagram in Figure 4.1.

Lemma 4.13 (Time interpolation). The operators J and CJ defined above are such

that

∫ T

0

∥x − J x∥2

X

≤ 4 inf

X∈S0

t (X)

∫ T

0

∥x − X∥2

X

(4.19a)

∫ T

0

∥x − CJx∥2

X ≲

inf

X∈S0

t (X)

∫ T

0

∥x − X∥2

X

(4.19b)

for all x ∈ L2(X). Moreover, for x ∈ H1(X), we have

(4.19c)

dt(J x, x(0)) = CJ∂tx.

The hidden constant in (4.19b) is an increasing function of

(4.20)

max

j=2,...,J

|Ij−1|

|Ij|

.

Proof. According to [40, Remark 4.1], the operator J is invariant on S0

t (X) and it

is bounded with

∫ T

0

∥J x∥2

X

≤ 4

∫ T

0

∥x∥2

X

for all x ∈ L2(X). The combination of the two properties implies (4.19a).

26

C. KREUZER AND P. ZANOTTI

Invariance of CJ on S0

t (X) follows by elementary calculations employing (4.18).

In order to prove stability for CJ, we observe for x ∈ L2(X) and j ≥ 2 that

∥( CJx)|Ij ∥X ≤

1

|Ij|

∫

Ij

∥x∥X

∫

Ij

|ψj| +

1

|Ij|

∫

Ij−1

∥x∥X

∫

Ij−1

|ψj−1|

≤

(∫

Ij

∥x∥2

X

\1

2 (∫

Ij

|ψj|2

\1

2

+

|Ij−1|

|Ij|

(∫

Ij−1

∥x∥2

X

\1

2 (∫

Ij−1

|ψj−1|2

\1

2

.

By recalling (4.17)-(4.18), we arrive at

∫

Ij

∥ CJx∥2

X

≤ 4

(

1 +

|Ij−1|

|Ij|

)∫

Ij−1∪Ij

∥x∥2

X

.

The same argument applies for j = 1. Summing over j = 1,...,J proves (4.19b).

Finally (4.19c) follows from the fundamental theorem of calculus and (4.18). □

H1(X)

( J , (·)(0) )

∂t

// L2(X)

˜

J

S

0

t (X) × X

dt

// S0

t (X)

Figure 4.1. Commutative diagram representing the relation be-

tween the time interpolants. The second component in the opera-

tor ( J , (·)(0) ) on the left is the evaluation at t = 0.

4.5. Space interpolation. Regarding the approximation in space, the error notion

ERR(·, ·) leads to the problem of interpolating functions from U into Us, from D into

Ds as well as from P and P∗ into Ps. Moreover, the spaces are related via the

operators D and L and their discrete counterparts and the error notion involves

various L2(Ω)-orthogonal projections. Therefore, commutative diagrams like the

one in Figure 4.1 are of interest also in this context.

In order to deal with all these tasks, we invoke the existence of an operator

which maps H1-conforming piecewise polynomials into H2-conforming ones and

enjoys several other properties, whose importance is made clear along the proof of

the next results in this section. The operator is defined on the Lagrange space of

degree k without boundary conditions, namely

S1

k := Pk(T) ∩ H1(Ω).

Moreover, we denote the jumps and the averages across faces by the usual symbols

{{·}} and [·], cf. [14, Section 38.2.1].

INF-SUP STABLE DISCRETIZATION OF THE BIOT’S EQUATIONS

27

Lemma 4.14 (Smoothing operator). There is a linear operator S : S1

k → H2(Ω)

which satisfies the following properties

S(Ps) ⊆ P

(4.21a)

∇SQ · n = 0

on Γp,N

(4.21b)

∫

T

(SQ)N =

∫

T

QN

∀N ∈ Pk(T), T ∈ T

(4.21c)

∫

F

(SQ)N =

∫

F

QN

∀N ∈ Pk−1(F), F ∈ Fi ∪ Fp,N

(4.21d)

as well as the stability estimate

∥Dj(Q − SQ)∥2

T ≲

∑

F∈Fi∪Fp,N , F∩T̸=∅

∫

F

{{h}}3−2j|[∇Q] · n|2

(4.21e)

for all Q ∈ S1

k and j ∈ {0, 1, 2}. The hidden constant depends only on the quantities

mentioned in Remark 4.1.

Proof. See Appendix A.

□

Having the operator S at hand, we define IΩ : L2(Ω) → S1

k via the problem

(4.22a)

(IΩ

Cptot,Q)Ω = (Cptot, SQ)Ω

∀Q ∈ S1

k

for Cptot ∈ L2(Ω). Analogously, we define IP∗ : P∗ → Ps via the problem

(4.22b)

(IP∗ Cm, N)Ω = ⟨ Cm, SN⟩P

∀N ∈ Ps

for Cm ∈ P∗. Note that the main difference between IΩ and IP∗ is in the range. We

introduce also IU : U → Us by

(4.22c)

IUCu := CU + Rs(IΩDCu − Ds CU)

for Cu ∈ U, where CU is the U-orthogonal projection of u onto Us and Rs : Ds → Us is

a right inverse of Ds. The existence and the boundedness of the latter operator are

equivalent to Proposition 4.2(2), cf. [14, Lemma C.42]. Notice also that we have

IΩDCu ∈ Ds in view of Lemma 4.15(2) below.

Let us clarify the logic behind the definitions in (4.22). The interpolant IΩ is

intended for the approximation of the total pressure, IP∗ for the total fluid content

and the pressure, and IU for the displacement. It is important for our purpose that

IP∗ is P∗-stable and this requires the use of an operator S in the right-hand side of

(4.22b), mapping the test functions (at least) into P. (We mention, in passing, that

S is actually required to map into even more regular functions, in order to enforce

the L2(Ω)-stability of the interpolant IP introduced below, cf. Lemma 4.17.) In

principle, the use of S is not needed in (4.22a), when only stability in the L2(Ω)-

norm is required. Still, it is important at some point relating IΩ and IP∗

via a

commutative diagram, cf. Figure 4.2 below. Therefore, we use S also in (4.22a).

Finally, the interpolant IU is nothing else than the U-orthogonal projection plus a

correction, that is necessary for another commutative diagram.

Lemma 4.15 (Space interpolation – Part 1). The operators IΩ, IP∗ and IU defined

in (4.22) enjoy the following properties.

(1) IΩ maps D into Ds and, for all Cptot ∈ L2(Ω), we have

∥Cptot − IΩ

Cptot∥Ω ≲ inf

̂Ptot∈S1

k

∥Cptot − CPtot∥Ω.

28

C. KREUZER AND P. ZANOTTI

(2) For all Cm ∈ P∗ and Cp ∈ P, we have, respectively,

∥ Cm − P∗

Ps IP∗ Cm∥P∗ ≲ inf

̂

M∈Ps

∥ Cm − P∗

Ps ̂

M∥P∗

and

∥Cp− IP∗ Cp∥Ω ≲ inf

̂P∈Ps

∥Cp− CP∥Ω.

(3) For all Cu ∈ U, we have DsIUCu = IΩDCu as well as

∥Cu − IUCu∥U ≲ inf

̂U∈Us

∥Cu − CU∥U.

(4) The following identities hold true in L2(Ω)

PDs IΩ = IΩPD

and

PPs IΩ = IP∗ P

P

.

The hidden constants depend only on the quantities mentioned in Remark 4.1.

Proof. (1) We first check that IΩ indeed maps D into Ds. Recall (2.5) and (4.5).

For D = L2(Ω) there is nothing to prove. If D = L2

0(Ω), then (4.22a) implies, for

Cptot ∈ D, that

(IΩ

Cptot, 1)Ω = (Cptot, S1)Ω = (Cptot, 1)Ω = 0.

Note that the identity S1 = 1 follows from (4.21e) with j = 1. This confirms the

inclusion IΩ

Cptot ∈ Ds. A similar argument reveals that IΩ is invariant on S1

k . In

fact, owing to (4.21c), we have, for Cptot ∈ S1

k , that

(IΩ

Cptot,Q)Ω = (Cptot, SQ)Ω = (Cptot,Q)Ω

∀Q ∈ S1

k .

Finally, for general Cptot ∈ H, the boundedness (4.21e) of S with j = 0, combined

with a discrete trace inequality [13, Lemma 12.8] and an inverse estimate [13,

Lemma 12.1], yields

∥IΩ

Cptot∥Ω = sup

Q∈S1

k

(Cptot, SQ)Ω

∥Q∥Ω

≲ ∥Cptot∥Ω.

The claimed estimate follows by the invariance and the boundedness of IΩ.

(2) For Cm ∈ P∗ and ̂

M ∈ Ps, we have

Cm − P

∗

Ps IP∗ Cm = Cm − P

∗

Ps ̂

M + P∗

Ps (̂

M − IP∗ Cm).

We bound the norm of the second summand on the right-hand side with the help of

Proposition 4.2(3), the inclusion Ps ⊆ L2(Ω) and (4.21e) with j = 1 and a discrete

trace inequality

∥P∗

Ps (̂

M − IP∗ Cm)∥P∗ ≂ ∥̂

M − IP∗ Cm∥P∗

s

= sup

N∈Ps

⟨̂

M − IP∗ Cm, N⟩Ps

∥N∥P

= sup

N∈Ps

⟨P∗

Ps ̂

M − Cm, SN⟩P

∥N∥P

≲ ∥P∗

Ps ̂

M − Cm∥P∗ .

The first claimed estimate follows by combining this bound with the above identity.

The other estimate follows by establishing invariance on Ps as well as boundedness

in the L2(Ω)-norm exactly as in (1).

(3) Recall that the operator Rs in (4.22c) is a right inverse of Ds, i.e. DsRs

is the identity on Ds. Then, the first part of the statement directly follows from

the definition of IU. We verify the second part as in (1), i.e. we prove that IU is

INF-SUP STABLE DISCRETIZATION OF THE BIOT’S EQUATIONS

29

invariant on Us and bounded in the U-norm. For Cu ∈ Us, we have CU = Cu in (4.22c).

Then, according to (4.21c) and (4.6), we infer

(IΩDCu, Qtot)Ω = (DCu, SQtot)Ω = (DCu, Qtot)Ω = (Ds

Cu, Qtot)Ω

∀Qtot ∈ Ds.

This implies IΩDCu = Ds

Cu and, in turn, IUCu = Cu. Next, for general Cu ∈ U, we have

∥IUCu∥U ≲ ∥CU∥U + ∥IΩDCu − Ds CU∥Ω.

Indeed, the boundedness of Rs is equivalent to 4.2(3); see [14, Lemma C.42]. We

conclude the proof of (3) with the help of the boundedness of IΩ, D and Ds in the

respective norms and by the definition of CU as the U-orthogonal projection of Cu.

(4) Recall that PD and PDs are the L2(Ω)-orthogonal projections onto the spaces

D and Ds, defined in (2.5) and (4.5), respectively. For Cptot ∈ L2(Ω), we have

IΩPD Cptot ∈ Ds, in view of (1) above. Then, for all Qtot ∈ Ds, it holds that

(IΩPD Cptot,Qtot)Ω = (PD Cptot, SQtot)Ω = (Cptot, SQtot)Ω = (IΩ

Cptot,Qtot)Ω.

In particular, we can remove PD after the second identity, because S maps Ds into

D, thanks to (4.21c). Hence, the first claimed identity is verified. The proof of the

other one is similar. Recall that P

P

and PPs are the L2(Ω)-orthogonal projections

onto the spaces P and Ps, defined in (2.4) and (4.7), respectively. For Cp ∈ L2(Ω)

and N ∈ Ps, we have

(IP∗ P

P Cp, N)Ω = (P

P Cp, SN)Ω = (Cp, SN)Ω = (IΩ

Cp, N)Ω.

This time we could remove P

P

after the second identity thanks to (4.21a). Thus,

also the second claimed identity is verified.

□

U

IU

D

// D

IΩ

Us

Ds

// Ds

P

IP

L

// P∗

IP∗

Ps

Ls

// P∗

s

L2(Ω)

IΩ

PD

// D

IΩ

S1

k

PDs

// Ds

L2(Ω)

IΩ

PP

// P

IP∗

S1

k

PPs

// Ps

Figure 4.2. Commutative diagrams representing the relations

among the space interpolants.

The interpolants defined up to this point are compatible with the divergence, i.e.

with the operators D and Ds, and with the various L2(Ω)-orthogonal projections,

as summarized in Figure 4.2. Still, the compatibility with the Laplacian, i.e. with

30

C. KREUZER AND P. ZANOTTI

the operators L and Ls, is not guaranteed. Since L and Ls are one-to-one, the

diagram on the upper right corner of Figure 4.2 uniquely defines IP : P → Ps as

(4.23)

IP := L−1

s IP∗ L.

The next lemma reveals that also IP has favorable approximation properties. To

this end, we introduce the following broken H2-norm

∥·∥2

H2(T) := κ

⎛

⎝

∑

T∈T

∥D2 · ∥2

L2(T) +

∑

F∈Fi∪Fp,N

∥{{h}}−1/2

[∇·]∥

2

L2(F)





as well as the constant

(4.24)

εs := inf

N∈Ps

∥LsN∥Ω

∥N∥H2(T)

.

Remark 4.16 (Discrete elliptic regularity). The constant εs in (4.24) is certainly

finite, because Ps is finite-dimensional. The size of the constant is related to the

control of a H2-like norm by the L2-norm of Ls, i.e. our discretization of the

Laplacian. Therefore, the constant is a discrete measure of the elliptic regularity.

Note also that εs can be equivalently interpreted as an inf-sup constant

εs = inf

N∈Ps

sup

̂P∈Ps

⟨LN, CP⟩P

∥N∥H2(T)∥ CP∥Ω

,

where the stability of the standard weak formulation of the Laplacian is analyzed

in a nonstandard (i.e. nonsymmetric) setting, cf. (4.8). We are aware of only a few

results regarding the size of εs. For H1-conforming Lagrange finite elements and

convex Ω, Makridakis established a lower bound in terms of the shape parameter

(4.1) of T, provided that the mesh is not too much graded, see [28]. It is unclear

to us whether the latter condition is necessary or not. When Ω is not convex,

the connection with the elliptic regularity suggests that a lower bound of εs only

in terms of shape regularity might be not possible. As for the question raised in

Remark 4.12, we postpone further investigation on this point to future work.

Lemma 4.17 (Space interpolation – Part 2). The operator IP defined above has a

unique bounded extension (denoted by the same symbol) to P which satisfies for all

Cp ∈ P the estimate

∥Cp− IP Cp∥Ω ≲ ε−1

s

inf

̂P∈Ps

∥Cp− CP∥Ω.

The hidden constant depends only on the quantities mentioned in Remark 4.1.

Proof. As for Lemma 4.15(2)-(3), we verify the claimed estimate by showing that

IP is invariant on Ps and that it is bounded in the L2(Ω)-norm. Note that the

latter property implies also the existence of a unique bounded extension of IP to P,

because P is a dense subspace thereof. Assume first Cp ∈ Ps. Owing to (4.23) and

(4.21a), we have

⟨IP Cp, LsN⟩Ps = ⟨LCp, SN⟩P = ⟨Cp, LSN⟩P

N ∈ Ps.

We recall (2.2) and integrate by parts element-wise [1, eq. (3.6)]

⟨IP Cp, LsN⟩Ps = κ

⎛

⎝−

∑

T∈T

∫

T

(∆Cp)SN +

∑

F∈Fi∪Fp,N

∫

F([∇Cp] · n)SN



 .

INF-SUP STABLE DISCRETIZATION OF THE BIOT’S EQUATIONS

31

Then, we make use of (4.21c)-(4.21d) and we integrate by parts back. It results

⟨IP Cp, LsN⟩Ps = ⟨Cp, LN⟩P= ⟨Cp, LsN⟩

Ps

where the last identity is obtained with the help of (4.8). Since Ls is one-to-one,

we conclude IP Cp = Cp.

Next, for general Cp ∈ P, a similar argument as before yields

∥IP Cp∥Ω = sup

N∈Ps

⟨IP Cp, LsN⟩Ps

∥LsN∥Ω

= sup

N∈Ps

⟨Cp, LSN⟩P

∥LsN∥Ω

.

Recall the inclusion SN ∈ H2(Ω) and (4.21a)-(4.21b). We integrate by parts, then

we invoke (4.21e) with j = 2 and (4.24)

∥IP Cp∥H = κ(Cp, −∆SN)Ω ≲ ∥Cp∥Ω∥N∥H2(T) ≤ ε−1

s ∥Cp∥Ω∥LsN∥Ω.

The combination of this bound with the previous identity implies the announced

boundedness of IP.

□

4.6. Space-time interpolation. We are in position to combine the interpolants

from the two previous sections in order to conclude the proof of Theorem 4.7. To

this end, let us preliminary recall that the operators acting in time commute with

those acting in space. Our main device is the space-time interpolant I : Y1 → Y1,st

defined as

ICy1 := (IUJ Cu, IΩJ Cptot, IP CJC

p, IP∗ J Cm, IP∗ Cm(0)

)

for Cy1 = (Cu, Cptot, Cp, Cm) ∈ Y1.

Our strategy to verify Theorem 4.7 is as follows. We first establish (4.15) for

y1 ∈ Y1. Then we extend the result by density, because the right-hand side of

(4.15) is bounded with respect to the norm ∥·∥1, cf. Theorem 2.4.

First of all, we recall the definitions of Y1 and Y1,st in (2.9) and (3.11), respec-

tively, and notice that I is well-defined. In fact, the discussion in the previous

sections shows that the interpolation of each component acts as follows

L2(U)

J

−→ S0

t (U)

IU

−→ S0

t (Us)

for Cu

L2(D)

J

−→ S0

t (D)

IΩ

−→ S0

t (Ds)

for Cptot

L2(P) ˜

J

−→ S0

t (P)

IP

−→ S0

t (Ps)

for Cp

H1(P∗)

J

−→ S0

t (P∗)

IP∗

−→ S0

t (Ps)

for Cm

P

∗ IP∗

−→ Ps

for Cm(0).

In particular, the second line makes use of Lemma 4.15(1) and the last one exploits

the inclusion Cm ∈ H1(P∗) ⊆ C0(P∗).

In order to verify Theorem 4.7, we bound the left-hand side of (4.15) by

(4.25)

inf

˜Y1∈Y1,st

ERR(Cy1, CY1)2 ≤ ERR(Cy1, ICy1)2.

Then, we recall the definition (4.11) of the error notion and we bound the six

terms therein (denoted by (i), (ii), ... for brevity) one by one. Lemma 4.13 and

Lemma 4.15(3) state that J and IU are bounded linear projections. Therefore,

their combination J IU is a L2(U)-bounded linear projection onto S0

t (Us). This

implies

(i) ≲

inf

̂U∈S0

t (Us)

∫ T

0

∥Cu − CU∥2

U

.

32

C. KREUZER AND P. ZANOTTI

The same reasoning with Lemma 4.15(1) implies

(ii) ≲

inf

̂Ptot∈S0

t (Ds)

∫ T

0

1

µ

∥Cptot − CPtot∥2

Ω.

Regarding the third term, we first recall (4.19c)

dt(IP∗ J Cm, IP∗ Cm(0)) = IP∗ dt(J Cm, Cm(0)) = IP∗ CJ∂t

Cm,

then, we make use of (4.23)

LsIP CJC

p = IP∗ CJ LC

p

and, combining the two identities, we arrive at

dt(IP∗ J Cm, IP∗ Cm(0)) + LsIP CJC

p = IP∗ CJ(∂t

Cm + LCp).

The same argument used in the above bounds for (i) and (ii), with Lemma 4.15(2),

implies

(iii) ≲

inf

̂

W ∈S0

t (Ps)

∫ T

0

∥∂t

Cm + LCp− P

∗

Ps ̂

W∥2

P∗ .

For the fourth term, we directly use Lemma 4.15(2) to obtain

(iv) ≲ inf

̂

M0∈Ps

∥ Cm(0) − P∗

Ps ̂

M0∥2

P∗ .

Concerning the fifth term, we invoke Lemma 4.15(3)

λDsIUJ Cu − IΩJ Cptot − αPDs IP CJC

p =

IΩJ (DCu − Cptot − αPDCp) + αIΩJ PD Cp− αPDs IP CJC

p

The reasoning from the bound of (ii) shows that the first summand in the right-

hand side is a near best approximation of (DCu− Cptot −αPDCp) in S

0

t (Ds). The other

two summands can be rewritten as αPDs (IΩJ −IP CJ)C

p, thanks to Lemma 4.15(4).

Arguing as before, we see that both IΩJ Cp and IP CJC

p are near best approximations

of Cp in S0

t (Ps), the latter one in view of Lemma 4.17. These observations, the

triangle inequality and the definition (2.12) of γ reveal

(v) ≲

inf

̂Qtot∈S0

t (Ds)

∫ T

0

1

λ + µ

∥(DCu−Cptot−αPD Cp)− CQtot∥2

Ω+ε−1

s

inf

̂P∈S0

t (Ps)

∫ T

0

1

γ

∥Cp− CP∥2

Ω.

The sixth and last term can be treated similarly. We invoke again Lemma 4.15(3)

and Lemma 4.15(4), so as to obtain

αPPs DsIUJ Cu + σIP CJC

p − IP∗ J Cm

= IP∗ J (αP

PDCu + σCp− Cm) + σ(IP CJ −IP∗ J )Cp.

Again, all operators yield near best approximation in the respective spaces. In

particular, for IP∗ J , we make use of the second part of Lemma 4.15(2). Thus we

arrive at the following estimate

(vi) ≲

inf

̂Q∈S0

t (Ps)

∫ T

0

γ∥(P

PDCu + σCp− Cm) − CQ∥2

Ω + ε−1

s

inf

̂P∈S0

t (Ps)

∫ T

0

1

γ

∥Cp− CP∥2

Ω

with the help of the definition (2.12) of γ. We insert the above bounds of (i)-(vi)

into (4.25). This establishes (4.15) for Cy1 ∈ Y1 and the right-hand side in the

resulting estimate is bounded in terms of the norm ∥·∥1, cf. Theorem 2.4. Thus,

INF-SUP STABLE DISCRETIZATION OF THE BIOT’S EQUATIONS

33

we can extend the estimate by density from Y1 to Y1. We conclude by noticing

that the bounds of (v) and (vi) simplify to

(v)+(vi) ≲ ε−1

s

inf

̂P∈S0

t (Ps)

∫ T

0

1

γ

∥Cp− CP∥2

Ω

if Cy1 = y1 is the solution of (2.8).

5. Numerical results

In this section we test the performance of the concrete space discretization pro-

posed in Section 4 on two well-established benchmarks for the Biot’s equations.

Owing to Corollary 4.8 and Remark 4.12, we restrict ourselves to the lowest order

case k = 1, corresponding to quadratic (resp. linear) H1-conforming Lagrange fi-

nite elements for the displacement (resp. for the other variables). All experiments

have been implemented with the help of the library ALBERTA 3.0, see [19, 35].

5.1. Terzaghi’s problem. The consolidation of a vertical soil column of total

depth H > 0 is a classical one-dimensional test case in poroelasticity. The column

is subject to compression and it is drained on top, whereas it is impermeable and

no displacement occurs at the bottom. Moreover, there is no other force acting in

the interior of the column, the initial fluid content equals zero and there are no

sources nor sinks. Therefore, in this case, the initial-boundary value problem (2.1)

for the Biot’s equations reads

(5.1)

−(2µ + λ)uzz + αpz = 0,

∂t(αuz + σp) − κpzz = 0, in (0,H) × (0,T)

−(2µ + λ)uz + αp = F,

p = 0, on {0} × (0,T)

u = 0,

pz = 0, on {H} × (0,T)

αuz + σp = 0, in (0,H) × {0}

with the subscripts (·)z and (·)zz denoting the partial derivatives with respect to

the depth z ∈ (0,H) and F the modulus of the force acting on top of the column.

Remarkably, the exact solution of (5.1) is explicitly known, cf. [32, Section 4.1.1].

In particular, the pressure equals

(5.2)

p(z, t) = p0

+∞

∑

m=0

4

(2m + 1)π

sin

((2m + 1)πz

2H

)

exp

(

−

(2m + 1)2π2

Cγkt

4H2

)

with the auxiliary parameters

Cγ :=

( α2

2µ + λ

+ σ

)−1

and

p0 :=

αCγF

2µ + λ

.

The displacement can be derived via the relation

uz =

αp − F

2µ + λ

and the other variables are obtained through their definition, cf. Remark 2.3.

In analogy with [32, Section 4.1.1], we set

H = 1,

T = 1,

F = 103

34

C. KREUZER AND P. ZANOTTI

as well as 1

µ = 41667, λ = 27778, α = 1, σ = 0.1, κ = 10−6.

The explicit knowledge of the exact solution makes this test case particularly

suited to observe the error decay rate. Owing to (2.12), (4.11), Remark 4.6 and

Theorem 4.7, we consider the following error notion

(5.3)

∫ T

0

(

∥u − U∥2

U