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
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
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.
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
The norm of is chosen such that and (A3) yields
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.
Let Assumption 1 hold. Then for any , , and , there exists a unique weak solution
If, in addition, and , then
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.
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.
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
Inserting this into equation (10) yields the Schur complement problem
From the properties of the operators in Lemma 3, we deduce the following assertions.
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
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.
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
and by integrating with respect to time, we further get
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
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);
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.
The following property of the dynamical system will serve as the basic tool for the stability analysis of approximation schemes in later sections.
Formal differentiation in time yields
which after rearrangement of the terms already yields the assertion of the lemma. ∎
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
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.
with some constant ,
If (A3h) is not valid, additional compatibility conditions for and are required.
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
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
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
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 solutionfollows for any choice of initial conditions and . By integration of the first equation, one obtains
with . Hence solves the original system (21) if, and only if, , i.e., when the compatibility condition
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
with being surjective and , both being positive definite. Differentiation of the second equation then leads to the system
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
has to be enforced, which is only caused by the inappropriate numerical approximation.
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
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.,
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.
With similar arguments as on the continuous level, we further obtain the following discrete energy–dissipation estimate for the discrete error.
Using , the embedding estimate , the definition of the norm , and Young’s inequality allow us to estimate
Then integrating with respect to time and using yield
The first estimate then follows by noting that and . From the discrete error equation (28) and continuity of the bilinear form , we deduce that
The discrete inf-sup condition (A3h) then leads to the second estimate of the lemma. ∎
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.
Let Assumption 1 hold. Then
If also the discrete inf-sup stability condition (A3h) holds, then additionally
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.
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
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