    # On the structure preserving high-order approximation of quasistatic poroelasticity

We consider the systematic numerical approximation of Biot's quasistatic model for the consolidation of a poroelastic medium. Various discretization schemes have been analysed for this problem and inf-sup stable finite elements have been found suitable to avoid spurios pressure oscillations in the initial phase of the evolution. In this paper, we first clarify the role of the inf-sup condition for the well-posedness of the continuous problem and discuss the choice of appropriate initial conditions. We then develop an abstract error analysis that allows us to analyse some approximation schemes discussed in the literature in a unified manner. In addition, we propose and analyse the high-order time discretization by a scheme that can be interpreted as a variant of continuous-Galerkin or particular Runge-Kutta methods applied to a modified system. The scheme is designed to preserve both, the underlying differential-algebraic structure and energy-dissipation property of the problem. In summary, we obtain high-order Galerkin approximations with respect to space and time and derive order-optimal convergence rates. The numerical analysis is carried out in detail for the discretization of the two-field formulation by Taylor-Hood elements and a variant of a Runge-Kutta time discretization. Our arguments can however be extended to three- and four field formulations and other time discretization strategies.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

In linear quasistatic theory  the consolidation of a poroelastic solid, which is fully saturated by an incompressible fluid, is usually described by Biot’s equations

 −div(2μϵ(u)+λdiv(u)I)+α∇p =f, (1) αdiv(\lx@overaccentset% \textbulletu)−div(κ∇p) =g, (2)

together with appropriate initial and boundary conditions. Here is the solid displacement of the porous medium and is the pressure of the residing fluid whereas denotes the density of external forces and is the fluid source density. The Biot parameter is usually close to one and the hydraulic conductivity is assumed strictly positive.

Due to its many applications, e.g., in geosciences or biomathematics, the theoretical and numerical analysis of (1)–(2) has attracted significant interest in the mathematical literature. The existence of unique solutions for the Biot system has been established by Zenisek [24, 25]

under some regularity conditions on the data via discretization with finite elements and the implicit Euler method using a-priori estimates and compactness arguments. Showalter

 established well-posedness for different formulations of the Biot model by the method of semi-groups. Murad and Loula [15, 16] investigated the Galerkin approximation by stable and unstable finite element pairs and established decay estimates for the discrete error via improved energy estimates. Mixed finite element approximations of the three-field formulation in which the seapage velocity is introduced as a new variable, were considered by Phillips and Wheeler [18, 19, 20]. Yi  investigated the discretization of a four-field formulation in which also the elastic stress field is introduced additionally. In  Kanschat and Riviere considered the approximation of the three-field formulation with a non-conforming approximation of the elastic deformation by finite elements. While most of the previous papers only utilized low order approximations in time, Bause et al. [2, 13] considered the efficient implementation of high order time approximations for poroelasticity, but no convergence analysis was conducted. Various available results concern the discretization of the static Biot systems that arise after time discretization, see e.g., [11, 17] in which stability with respect to model and discretization parameters has been studied. It is well-accepted by now [16, 20, 28], that inf-sup stable finite element pairs should be used to avoid spurious pressure oscillations that might appear in the initial phase of the simulations.

In this paper, we consider the systematic construction of high-order approximations for the quasistatic system (1)–(2) by Galerkin methods in space and time. We give a short and concise proof of well-posedness of the continuous problem which motivates our functional analytic setting and provides guidelines for the choice of initial conditions. The regularity requirements for the data are based on the physically relevant energy-dissipation structure of the evolution problem and compatibility conditions for the initial values are derived from the differential-algebraic structure. We then consider the Galerkin approximation in space and discuss the properties of the resulting differential-algebraic equations. In particular, we show that the index of the differential-algebraic system [5, 14] is one, independent of the approximation spaces, while a discrete inf-sup condition, i.e., the surjectivity of the discrete divergence operator, is required to guarantee the well-posedness of the semi-discretization under natural compatibility conditions for the initial conditions. Otherwise, additional non-physical constraints for the initial conditions arise. Afterwards, we show that the discrete error, i.e., the difference between the Galerkin semi-discretization and an appropriate elliptic projection, only depends on the approximation error in the displacement which implies a certain super-convergence for the discrete error. We then consider the time-discretization by a strategy which is capable of preserving both, the particular differential-algebraic structure and the energy-dissipation property of the underlying problem. The resulting scheme can be interpreted as a continuous Galerkin method or a variant of certain Runge-Kutta methods applied to a reformulation of the problem in which the algebraic equation is differentiated. The latter interpretation also allows us for an efficient implementation.

The remainder of the manuscript is organized as follows: In Section 2, we introduce an abstract evolution problem which covers the weak formulation of the Biot system (1)–(2) as a special case, and we establish its well-posedness under mild assumptions on the problem structure and the data. A natural compatibility condition for the initial conditions and the underlying energy-dissipation structure are presented together with the resulting a-priori bounds. In Section 3, we discuss the Galerkin approximation of the abstract problem and investigate the effect of the discrete inf-sup stability on the index of the resulting differential-algebraic system. We then establish an abstract convergence result for inf-sup stable approximations. In Section 4, we propose a time discretization scheme which can be interpreted as an inexact continuous Galerkin approximation or a variant of certain Runge-Kutta collocation methods. In the spirit of , this method can be interpreted in a pointwise sense which allows us to show that it preserves the underlying energy-dissipation structure and to conduct the error analysis with similar arguments as on the continuous level. In Section 5, we apply the results to the discretization of the Biot system by Taylor-Hood finite elements and derive order optimal error estimates in space and time. Numerical tests are presented in Section 6 for illustration of our theoretical results, and some directions for possible extensions and further investigations are discussed in the last section.

## 2. An abstract model problem

We now introduce an abstract model problem which covers the Biot system as a special case, and then establish its well-posedness under mild structural assumptions on the data. In addition, we briefly discuss the regularity of the solution as well as the underlying energy-dissipation identity which us later used to establish a-priori estimates.

### 2.1. Notation and summary of results

Let , , be real Hilbert spaces with compact embedding . By identifying with its dual , we obtain a Gelfand-triple ; see [7, Chapter XVIII]. The symbol will be used to denote the duality product on , , and .

In the following, we consider an abstract evolution problem stated in weak form as

 a(u(t),v)−b(v,p(t)) =⟨f(t),v⟩,∀v∈V, t>0, (3) b(\lx@overaccentset\textbulletu(t),q)+k(p(t),q) =⟨g(t),q⟩,∀q∈Q, t>0. (4)

Here , , and are given bilinear forms and , are functions of time with values in and , respectively. To establish the well-posedness for the corresponding Cauchy problem, we make the following assumptions.

###### Assumption 1.
• is bounded, elliptic, and symmetric.

• is bounded, elliptic, and symmetric.

• is bounded and inf-sup stable; see below.

We further use and as scalar products on and , and thus obtain

 a(u,v)≤∥u∥V∥v∥Vanda(v,v)=∥v∥2V, (5) k(p,q)≤∥p∥Q∥q∥Qandk(q,q)=∥q∥2Q. (6)

The norm of is chosen such that and (A3) yields

 b(v,q)≤Cb∥v∥V∥q∥Q0andsup∥v∥V=1b(v,q)≥β∥q∥Q0, (7)

for some . All estimates hold uniformly with respect to their arguments.

The following theorem summarizes the basic results about well-posedness and regularity of the solution for the above system with appropriate initial conditions.

###### Theorem 2.

Let Assumption 1 hold. Then for any , , and , there exists a unique weak solution

 u∈C(0,T;V) with Bu∈H1(0,T;Q∗), p∈L2(0,T;Q)∩H1(0,T;Q∗)∩C(0,T;Q0),

which satisfies (3)–(4) for a.e. and the initial conditions , where is defined by the compatibility conditions

 a(u0,v)−b(v,p0) =⟨f(0),v⟩,∀v∈V. (8)

If, in addition, and , then

 p∈L∞(0,T;Q)∩H1(0,T;Q0)andu∈H1(0,T;V).

In both cases, the solution can be bounded by the data in the natural norms.

The remainder of this section is devoted to a proof of these assertions.

### 2.2. Well-posedness

As usual, we associate via , , and to any of the bilinear forms a linear operator , , and its adjoint . This allows us to rewrite (3)–(4) as equivalent operator equations

 Au(t)−B∗p(t) =f(t),t>0, (9) B\lx@overaccentset\textbulletu(t)+Kp(t) =g(t),t>0, (10)

which have to be understood as equations in the sense of and , respectively. From the properties of the bilinear forms, we immediately conclude the following result.

###### Lemma 3.

Let Assumption 1 hold. Then

• and are symmetric, bounded, and boundedly invertible
with , and , .

• is bounded, surjective with closed range, and .

• is bounded and injective with .

Using property (i) of the previous lemma, we can eliminate from (9) leading to

 u(t) =A−1(f(t)+B∗p(t)). (11)

Inserting this into equation (10) yields the Schur complement problem

 BA−1B∗\lx@overaccentset\textbulletp(t)+Kp(t) =g(t)−BA−1\lx@overaccentset\textbulletf(t). (12)

From the properties of the operators in Lemma 3, we deduce the following assertions.

###### Lemma 4.

Let Assumption 1 and the conditions on the data in Theorem 2 hold. Then

• The operator is symmetric, bounded, and elliptic.

• .

Note that, according to property (v), the operator induces a symmetric, bounded, and elliptic bilinear form , and the operator equation (12) can hence be written in an equivalent weak form as

 c(\lx@overaccentset\textbulletp(t),q)+k(p(t),q) =⟨h(t),q⟩,∀q∈Q, t>0. (13)

This is an abstract parabolic equation whose well-posedness can be proven via Galerkin approximation; see [9, Chapter 7] or [7, Chapter XVIII]. We thus obtain

###### Lemma 5.

Let Assumption 1 hold. Then for any and , the reduced evolution problem (13) has a unique weak solution

 p∈L2(0,T;Q)∩H1(0,T;Q∗) (14)

with initial value . Moreover, by embedding.

This proves existence and uniqueness of a solution for problem (13) as well as a-priori estimates in the corresponding norms. Inserting into (11) and using the mapping properties of the operators and yields . The compatibility condition (8) then follows from continuity. Furthermore from (10), one can see that . By equivalence of (9)–(10) with the variational form, we obtain existence and uniqueness of a weak solution for the system (3)–(4) with the given initial conditions. Moreover, we have established the a-priori bounds in the first part of Theorem 2.

### 2.3. Regularity

We now turn to the improved a-priori estimates stated in Theorem 2. Additional regularity of the solution to the reduced problem (13) can be obtained with similar arguments as in [9, Chapter 7]. Assume that and with and . Then from the properties of the operators and , we deduce that with and . By formally testing (13) with , we obtain

 c(\lx@overaccentset\textbulletp,\lx@overaccentset\textbulletp)+ddt12k(p,\lx@overaccentset\textbulletp) =⟨h,\lx@overaccentset\textbulletp⟩,

and by integrating with respect to time, we further get

 ∫t0c(\lx@overaccentset\textbulletp,\lx@overaccentset\textbulletp)dt +k(p(t),p(t))≤k(p(0),p(0))+∫t0⟨h1,\lx@overaccentset\textbulletp⟩−⟨\lx@overaccentset\textbulleth2,p⟩dt+⟨h2,p⟩|t0 ≤32k(p(0),p(0))+C1∫t0∥\lx@overaccentset\textbulletf∥2V∗+∥g1∥2Q∗0+∥\lx@overaccentset\textbulletg2∥2Q∗dt +C2(∥g2(0)∥2Q∗+∥g2(t)∥2Q∗)+12∫t0c(\lx@overaccentset\textbulletp,\lx@overaccentset\textbulletp)dt+12k(p(t),p(t)).

Here we used and , which follow directly from the properties of the bilinear forms and norms. The last two terms in the above estimate can be absorbed into the left hand side, which yields the required estimates for . Inserting this into (11) leads to the improved regularity for . This concludes the proof of Theorem 2.

### 2.4. Choice of initial conditions

Before we proceed, let us discuss in more detail the choice of initial conditions. Under Assumptions 1, the variational problem

 a(u0,v)−b(v,p0) =⟨f0,v⟩,∀v∈V, (15) b(u0,q) =⟨ϕ0,q⟩,∀q∈Q, (16)

is uniquely solvable for all and , which is a direct consequence of Brezzi’s splitting lemma . In the above arguments, we simply chose and specified . By assumption (A1), we can determine from (15), and inserting into the equation (16) determines . Alternatively, one could set , choose freely, and then determine , by solving the coupled system (15)–(16). This again provides initial values and satisfying (8), where can be used as initial value for (12). The following choices of initial conditions are therefore equivalently possible:

• choose and determine by (15);

• choose and determine , by (15)–(16).

While the first choice is the natural one for the reduced problem (13), the second choice seems more natural for the coupled system (3)–(4). As indicated above, both choices are possible and actually equivalent in the considered functional analytic setting.

### 2.5. Energy-dissipation

The following property of the dynamical system will serve as the basic tool for the stability analysis of approximation schemes in later sections.

###### Lemma 6.

Let Assumption 1 hold and let denote a regular solution of (9)–(10) in the sense of Theorem 2. Then

 ddt12a(u,u)+k(p,p) =⟨f,\lx@overaccentset\textbulletu⟩+⟨g,p⟩. (17)
###### Proof.

Formal differentiation in time yields

 ddt12a(u,u) =a(u,\lx@overaccentset\textbulletu)=⟨f,\lx@overaccentset\textbulletu⟩+b(\lx@overaccentset\textbulletu,p)=⟨f,\lx@overaccentset\textbulletu⟩+⟨g,p⟩−k(p,p),

which after rearrangement of the terms already yields the assertion of the lemma. ∎

###### Remark 7.

By integration in time, the validity of the stability estimate can be extended to less regular weak solutions, and these estimates provide an alternative route for proving uniqueness and a-priori estimates for weak solutions.

## 3. Galerkin approximation

For discretization of the variational problem (3)–(4), we now consider Galerkin approximations in finite dimensional sub-spaces and , i.e., we search for semi-discrete functions and satisfying

 a(uh(t),vh)−b(vh,ph(t)) =⟨f(t),vh⟩,∀vh∈Vh, t>0, (18) b(\lx@overaccentset\textbulletuh(t),qh)+k(ph(t),qh) =⟨g(t),qh⟩,∀qh∈Qh, t>0, (19)

together with appropriate initial conditions to be specified below. Using similar arguments as for the analysis on the continuous level, we will show the following result.

###### Theorem 8.

Let Assumption 1 hold. Then (18)–(19) is a regular system of differential-algebraic equations of index . If is inf-sup stable on , i.e.,

• with some constant ,

then a unique solution to the system (18)–(19) exists for any choice of initial values and satisfying

 a(uh,0,vh)−b(vh,ph,0) =⟨f(0),vh⟩,∀vh∈Vh. (20)

If (A3h) is not valid, additional compatibility conditions for and are required.

###### Remark 9.

A quick comparison with Theorem 2 shows that (20) are the natural compatibility conditions for the initial values of the problem under consideration, while the additional conditions required when (A3h) is not valid are artificial; see below. Inf-sup stable approximation spaces , are therefore required to guarantee the well-posedness of the semi-discrete problem without artificial conditions on the initial values.

In the following, we give a detailed proof of the above theorem and then turn to the error analysis of the Galerkin approximations defined in the beginning of the section.

### 3.1. Proof of Theorem 8

Choice of a basis for the spaces , allows to convert the semi-discrete problem into an equivalent system of differential-algebraic equations

 (00B0)(\lx@overaccentset% \textbulletu\lx@overaccentset\textbulletp)+(A−B⊤0K)(up)=(fg), (21)

with and

denoting the coordinate vectors of the discrete solutions and data, and

, , being matrices of appropriate size. From (A1)–(A2), one can see that the matrices and are symmetric and positive definite, and therefore, the matrix

 S(λ)=(00B0)+λ(A−B⊤0K),

is positive definite, e.g., for . As a consequence, the matrix pencil is regular and so the system (21) is a regular differential-algebraic equation; see [5, 14] for details. Differentiating the first equation in (21) leads to

 (A−B⊤B0)(\lx@overaccentset% \textbulletu\lx@overaccentset\textbulletp)+(000K)(up)=(\lx@overaccentset% \textbulletfg). (22)

Note that the matrix in front of the time derivatives is regular if, and only if, is surjective. In that case, (22

) is an (implicit) ordinary differential equation, and the existence of a unique solution

follows for any choice of initial conditions and . By integration of the first equation, one obtains

 Au(t)−B⊤p(t)=f(t)+c,

with . Hence solves the original system (21) if, and only if, , i.e., when the compatibility condition

 Au0−B⊤p0=f(0) (23)

is satisfied. If, on the other hand, is not surjective, then (22) is still a differential-algebraic equation. By change of basis, we may transform the system into

 =⎛⎜⎝\lx@overaccentset\textbulletfg1g2⎞⎟⎠,

with being surjective and , both being positive definite. Differentiation of the second equation then leads to the system

 =⎛⎜⎝\lx@overaccentset\textbulletf\lx@overaccentset\textbulletg1g2⎞⎟⎠. (24)

Using the fact that is surjective and and are positive definite, one can deduce that the matrix in front of the time derivatives is regular, hence (24) is an ordinary differential equation. An inspection of the right hand side shows that none of the equations has been differentiated more than once, and hence the index of (21) is also one in this case. To obtain equivalence of (24) with (21), not only the compatibility condition (23) is required, but the additional artificial condition

 K1p1(0)=g1(0) (25)

has to be enforced, which is only caused by the inappropriate numerical approximation.

The assertions of Theorem 8 follow immediately from the above results by equivalence of the differential-algebraic system (21) with the weak formulation (18)–(19).

### 3.2. Abstract error analysis

We next turn to the a-priori analysis of Galerkin discretizations introduced in the beginning of this section. As usual, see e.g. [22, 23, 26], we decompose the error between the continuous and the semi-discrete solution via

 ∥u−uh∥V≤∥u−~uh∥V+∥~uh−uh∥V, ∥p−ph∥Q≤∥p−~ph∥Q+∥~ph−ph∥Q,

into approximation errors and discrete error components. Following [15, 16], we utilize the variational problem corresponding to the stationary system associated with (3)–(4) to define the approximation and , i.e.,

 a(~uh(t)−u(t),vh)−b(vh,~ph(t)−p(t)) =0,∀vh∈Vh, t>0, (26) k(~ph(t)−p(t),qh) =0,∀qh∈Qh, t>0. (27)

Note that and can be computed by solving elliptic variational problems and error estimates for the elliptic projection can therefore be obtained by standard arguments; see Section 5 below. From the definition of the projections and the discrete solution, we immediately obtain the following discrete error equation.

###### Lemma 10.

Let Assumption 1 hold and denote a solution of (18)–(19) with initial values and . Then the discrete error components and satisfy , , and

 a(δuh(t),vh)−b(vh,δph(t)) =0, ∀vh∈Vh, t>0, (28) b(δ\lx@overaccentset\textbulletuh(t),qh)+k(δph(t),qh) =b(\lx@overaccentset\textbulletu(t)−\lx@overaccentset\textbullet~uh(t),qh), ∀qh∈Qh, t>0. (29)

With similar arguments as on the continuous level, we further obtain the following discrete energy–dissipation estimate for the discrete error.

###### Lemma 11.

Let Assumption 1 hold and be a solution of (18)–(19). Then

 ∥δuh(t)∥2V+∫t0∥δph(s)∥Qds≤∫t0∥B(\lx@overaccentset%\textbulletu(s)−\lx@overaccentset\textbullet~uh(s))∥2Q∗0ds,0≤t≤T.

If also (A3h) holds, then additionally

 β∥δph(t)∥Q0≤∥δuh(t)∥V.
###### Proof.

From the discrete error equations (28)–(29), one can deduce that

 ddt12a(δuh,δuh) =a(δuh,δ\lx@overaccentset% \textbulletuh)=b(δ\lx@overaccentset% \textbulletuh,δph)=b(\lx@overaccentset% \textbulletu−\lx@overaccentset\textbullet~uh,δph)−k(δph,δph).

Using , the embedding estimate , the definition of the norm , and Young’s inequality allow us to estimate

 b(\lx@overaccentset\textbulletu−\lx@overaccentset\textbullet~uh,δph)≤12∥B(\lx@overaccentset\textbulletu−\lx@overaccentset\textbullet~uh)∥2Q∗0+12k(δph,δph).

Then integrating with respect to time and using yield

 a(δuh(t),δuh(t))+∫t0k(δph(s),δph(s))ds≤∫t0∥B(\lx@overaccentset% \textbulletu(s)−\lx@overaccentset\textbullet~uh(s))∥2Q∗0ds.

The first estimate then follows by noting that and . From the discrete error equation (28) and continuity of the bilinear form , we deduce that

 b(vh,δph(t))=a(δuh(t),vh)≤∥δuh(t)∥V∥vh∥V.

The discrete inf-sup condition (A3h) then leads to the second estimate of the lemma. ∎

###### Remark 12.

Let us emphasize that the particular choice of the elliptic projection in the error decomposition allows us to bound the discrete error by the approximation error in the component alone; this leads to improved error estimates and allows the use of post-processing techniques to obtain approximations for the pressure in polynomial spaces of higher order; see [4, 16] for details.

Using the previous bounds for the discrete error components, we now obtain the following estimates for semi-discrete approximation (18)–(19).

###### Theorem 13.

Let Assumption 1 hold. Then

 ∥u(t)−uh(t)∥2V ≤∥u(t)−~uh(t)∥2V+∫t0∥B(\lx@overaccentset\textbulletu(s)−\lx@overaccentset\textbullet~uh(s))∥2Q∗0ds, (30) ∫t0∥p(s)−ph(s)∥2Qds ≤∫t0∥p(s)−~ph(s)∥2Q+∥B(\lx@overaccentset\textbulletu(s)−\lx@overaccentset\textbullet~uh(s))∥2Q∗0ds. (31)

If also the discrete inf-sup stability condition (A3h) holds, then additionally

 ∥p(t)−ph(t)∥2Q0 ≤∥p(t)−~ph(t)∥2Q0+β−2h∫t0∥B(\lx@overaccentset\textbulletu(s)−\lx@overaccentset\textbullet~uh(s))∥2Q∗0ds. (32)

After choosing the approximation spaces and , this estimate allows to derive quantitative error bounds via corresponding estimates for the elliptic projection. Details for a particular discretization will be given in Section 5.

## 4. Time discretization

We now turn to the discretization of the semi-discrete variational problem (18)–(19) in time. The guiding principle will be to preserve the underlying differential-algebraic structure and energy–dissipation identity elaborated in Section 3 as good as possible.

### 4.1. Preliminaries

Let be a partition of the time interval with time steps and . We denote by the space of piecewise polynomial functions of with values in . Following the notation of , we utilize capital letters to denote piecewise polynomial functions of time in the following.

In the spirit of  and to highlight the preservation of the problem structure, we first give a pointwise definition of our method, which requires the following two projection operators in time: We denote by the -orthogonal projection and by , , the -conforming projection, which is defined by the relations

 ∂tΠ1qu =Π0q−1∂tuand (33) (Π1qu)(tn−1) =u(tn−1),1≤n≤N. (34)

We will refer to (33) as the commuting-diagram property. By integration of this relation and use of (34), one can see that and consequently is continuous on ; thus is an -conforming projection in time. We now consider the following time-discretization of the semi-discrete problem.

###### Problem 14 (Fully discrete scheme).

Find and such that , , and

 <