{'fraction_non_alphanumeric': 0.04873159604617405, 'fraction_numerical': 0.024192990944909103, 'mean_word_length': 4.145092460881934, 'pattern_counts': {'":': 1, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A priori error estimates of two fully discrete coupled schemes for Biot's consolidation model

This paper concentrates on a priori error estimates of two fully discrete coupled schemes for Biot's consolidation model based on the three-field formulation introduced by Oyarzua et al. (SIAM Journal on Numerical Analysis, 2016). The spatial discretizations are based on the Taylor-Hood finite elements combined with Lagrange elements for the three primary variables. For time discretization, we consider two methods. One uses the backward Euler method, and the other applies a combination of the backward Euler and Crank-Nicolson methods. A priori error estimates show that the two schemes are unconditionally convergent with optimal error orders. Detailed numerical experiments are presented to validate the theoretical analysis.


Introduction
Biot's consolidation model [1,2] describes the interaction between fluid flow and mechanical deformation of a porous medium, which is a solid structure composed of pores. This model has wide applications in many scenarios including biomechanics [14], petroleum engineering [15], etc. Since the exact analytic solution is hard to obtain, a variety of methods are carried out to approximate numerical solutions for the system: finite volume methods [20], virtual element methods [7], and mixed finite element methods [8,16].
In [24], it is pointed out that standard finite element methods for solving the classical two-field Biot's model may suffer from Poisson locking and pressure oscillations. Therefore, various reformulations of Biot's model are proposed to overcome these numerical difficulties. For example, fluid flux arising from the inherent Darcy law is introduced as a new variable to obtain a three-field formulation of Biot's model in [24,13]. A new four-field formulation is proposed in [25]. Another three-field model is established in [10] by introducing two pseudo-pressures. In [22,18], an intermediate variable, called "total pressure", is introduced to derive a three-field reformulation of Biot's model. The advantages of such a three-field reformulation exist in that it avoids to use H(div) space, and the classical inf-sup stable Stokes finite elements combined with Lagrange elements can be applied for the spatial discretization. By taking the advantage of such a three-field formulation, some relevant algorithms and analyses are carried out. For instance, Qi et al. [23] derive optimal-order error estimates. An unconditional convergent algorithm is proposed in [17].
In this work, we study two fully discrete coupled schemes for solving the three-field formulation of Biot's model. In Method 1, the backward Euler method [22,14,23,7] is applied for the time discretization, which is only of first-order convergence in time. Inspired by [17], we come up with Method 2 by using a combination of the backward Euler and Crank-Nicolson methods for a second-order convergence in time. Different from the algorithm in [17], our algorithm applies a unified scheme for all time steps rather than uses different schemes for the first time step and the rest steps. Despite some existing analyses, we note that some works only assume that the Biot-Willis constant α = 1 and the specific storage coefficient c 0 = 0. The only second order in time paper is [17]. However, the analyses and numerical experiments there are not consistent. In this work, we conduct rigorous analysis for both Method 1 and Method 2. Optimal a priori estimates of the two methods are obtained by more careful usage of energy estimates and Grönwall's inequality. Furthermore, our analyses are valid for more general physical parameters. We comment here that most existing analyses only discuss the L 2 error estimates for the fluid pressure. In this work, we also provide the H 1 error analysis for the fluid pressure.
The rest of the paper is organized as follows. In Section 2, we present a three-field formulation of Biot's consolidation model and the corresponding weak formulation. In Section 3, we introduce finite element spaces, projection operators, two fully discrete coupled schemes, and some useful propositions. A priori estimates of Method 1 and Method 2 are given in Section 4. Numerical experiments are carried out to validate the theoretical results in Section 5. Conclusions and outlook are given in Section 6.

Mathematical formulations
Let Ω ⊂ R d (d = 2 or 3) be a bounded polygonal domain with boundary ∂Ω. The classical Sobolev spaces are denoted by H k (Ω) with norm · H k (Ω) . We denote H k 0,Γ (Ω) for the subspace of H k (Ω) with the vanishing trace on Γ ⊂ ∂Ω, and use (·, ·) and ·, · to denote the standard L 2 (Ω) and L 2 (∂Ω) inner products, respectively. In this paper, we use C to denote a generic positive constant independent of mesh sizes, and use x y to denote x ≤ Cy. The governing equations describing the quasi-static Biot system are given as follows where Here, the primary unknowns are the displacement vector of the solid u and the fluid pressure p, the coefficient α > 0 is the Biot-Willis constant which is close to 1, f is the body force, c 0 ≥ 0 is the specific storage coefficient, K represents the hydraulic conductivity, ρ f is the fluid density, g is the gravitational acceleration, Q s is a source or sink term, I is the identity matrix, and Lamé constants λ and µ are computed from the Young's modulus E and the Poisson ratio ν: .
Suitable boundary and initial conditions should be provided to complete the system. Assuming that ∂Ω = Γ d ∪ Γ t = Γ p ∪ Γ f with |Γ d | > 0 and |Γ p | > 0, for the simplicity of presentation, we consider the following boundary conditions: where n is the unit outward normal to the boundary. We comment here that the discussion can be easily extended to nonhomogeneous boundary condition cases. For the initial conditions, we consider Next, we introduce an intermediate variable called "total pressure": Then, (1)-(2) can be rewritten as the following three-field formulation of Biot's consolidation model After the reformulation, we can still apply the boundary conditions (3)-(6) and initial conditions (7) with ξ(0) = αp 0 − λdivu 0 here. For ease of presentation, we assume g = 0 in the rest of the paper.
(Ω) and M = H 1 0,Γp (Ω), and define the bilinear forms Multiplying (8)-(10) by test functions, integrating by parts and applying boundary conditions (3)-(6) yield the following variational problem: for a given t ≥ 0, find (u, ξ, p) ∈ V × W × M such that The well-posedness of problem (11)-(13) is established in [22]. We note that the Korn's inequality [21] holds on V , that is, there exists a constant Furthermore, the following inf-sup condition [5] holds: there exists a constant β > 0 depending only on Ω and Γ d such that

Finite element discretization and numerical schemes
Let T h be a partition of the domain Ω into triangles in R 2 or tetrahedra in R 3 , and h be the maximum diameter over all elements in the mesh. We define finite element spaces on T h where k ≥ 2 and l ≥ 1 are two integers. In this work, the Taylor-Hood element, which consists of the pair (V h , W h ), and Lagrange finite element are adopted to the pair (u, ξ) and p, respectively. Based on the three-field formulation (8)-(10) and discrete spaces V h , W h , and M h , we define two projection operators [19,23]. First, we introduce the Stokes projection operator Second, we define the elliptic projection operator R p : M → M h as follows If u ∈ H k+1 0,Γ d (Ω), ξ ∈ H k (Ω), and p ∈ H l+1 0,Γp (Ω), then the following error estimates hold true for the Stokes projection operator and the elliptic projection operator [6].
As we use a stable Stokes element pair, the corresponding finite element spaces satisfy the following discrete inf-sup condition, i.e., there exists a positive constant β independent of h such that An equidistant partition 0 = t 0 < t 1 < · · · < t N +1 = T with a step size ∆t is considered for the time discretization. For simplicity, we introduce the notations u n = u(t n ), ξ n = ξ(t n ) and p n = p(t n ). Suitable approximation of initial conditions u 0 h = R u u 0 , ξ 0 h = R ξ ξ 0 , and p 0 h = R p p 0 is considered here. Following [22,14,23,7], we present the first coupled scheme using the backward Euler method for the time discretization as follows.
The second coupled scheme considers a combination of the backward Euler and Crank-Nicolson methods for time discretization, which is inspired from [17]. In the entire time interval, we solve the following problem.
We comment here that such a scheme does not complicate the calculations. Next, we state the following basic propositions.
Then, the following estimate for the L 2 -norm of the last term holds true.
Proof. The first part of the conclusion comes from the Taylor expansion theorem.
Using the Cauchy-Schwarz inequality, we have which implies (28) directly.

Proposition 3.2. Let B be a symmetric bilinear form, there holds
which immediately implies the following inequality

Main results
In this section, we present our main results describing a priori error estimates of the proposed schemes. Here, we decompose error terms as We also define Then, we present a priori error estimates for both schemes. The proof of each method consists of three parts. In the first part, we focus on the H 1 norm of e h,N +1 u , L 2 norm of e h,N +1 ξ , and L 2 norm of e h,N +1 p . In the second part, we complete the H 1 norm estimate of e h,N +1 p . In the third part, we draw our main conclusions describing the error bound.
By using the assumptions of the projection operators (15), (16), and (17), the above equations can be rewritten as Differentiating (12) with respect to t at the (n + 1)-th time step, then multiplying the resulted equation by ∆t, we derive that For (34), we write out the schemes for t = t n+1 and t = t n , then take a difference between the two resulted equations, we see that Using the definitions of D n+1 ξ and D n+1 p , we can reformulate (37) as follows After using the (n + 1)-th and n-th time steps of (12) Taking the summation of (40), (41), and (42) over the index n from 0 to N yields Using the definitions of a 2 (·, ·), a 3 (·, ·), and c(·, ·), we can simplify LHS 1 by the identity Applying (31) and (44), we obtain the following lower bound estimate for LHS 1 Next, we bound the terms E i for i = 1, 2, · · · , 5. We use the Cauchy-Schwarz inequality, the Young's inequality, (18), (20) and (28) to estimate E 1 , E 2 , and E 3 with an ǫ 1 > 0 as follows.
We note that (54)  and ǫ 2 D n+1 Applying (54) and (55) to the above equation, we claim that the proof is completed.

Benchmark tests
In this section, we present numerical experiments in two dimensions to validate the theoretical predictions described in Section 4. All computations are implemented by using the open-source software FreeFEM++ [12].
Example 1. Let the domain Ω = [0, 1] 2 and the final time is T = 1.0. We choose the body force f , the source/sink term Q s , initial conditions and Dirichlet boundary data on ∂Ω = Γ d = Γ p such that the exact solution is as follows: Following [9], the physical parameters are given by: We apply a small mesh size h = 1 64 and take polynomial orders k = 3, l = 2 for the spatial discretization so that the spatial error is not dominant. To check the orders of convergence in time, we only refine the time step size ∆t. In Table 1 and Table 2, we present the results of errors and convergences rates for Method 1 and Method 2, respectively. We observe that the orders of H 1 error of u, L 2 error of ξ, L 2 and H 1 errors of p are all around 1 in Table 1, and are all around 2 in Table 2    Other physical parameters will vary to test the robustness of our numerical schemes. Numerical results for this example are summarized in Tables 3-10. Among them, we examine Method 1 in Tables 3, 4, 7, 8, and Method 2 in Tables 5, 6, 9, 10. Since we have already verified the error orders in time for both schemes in Example 1, our focus here is the verification of the spatial error orders. To verify the spatial error orders as analyzed in Theorem 4.3, we take ∆t of an order O(h 2 ) ( Table 3 and Table  7) or O(h 3 ) ( Table 4 and Table 8) for Method 1. Similarly, to verify the spatial error orders as analyzed in Theorem 4.6, we take ∆t of an order O(h) ( Table 5 and Table  9) or O(h 2 ) ( Table 6 and Table 10) for Method 2. Firstly, we fix ν = 0.3 and K = 1.0. The numerical results for errors and convergence orders using k = 2, l = 1 and k = 3, l = 2 are presented in Table 3,  5 and Table 4, 6, respectively. We refine the mesh size (and vary the corresponding time step size) to present the numerical results. From Table 3, it is clearly shown that the convergence e N +1      Secondly, we fix ν = 0.49999 and K = 10 −6 to test the robustness of the proposed schemes with respect to the key physical parameters. The numerical results for errors and convergence orders using k = 2, l = 1 and k = 3, l = 2 are presented in Table 7, 9 and Table 8, 10, respectively. By checking the error results and convergence rates one table by one table, one can verify the theoretical analysis provided in Section 4. From these tables, it is shown clearly that all energy norm errors decrease with the optimal convergence orders. By comparing the results in Tables 3-6 with the corresponding results in Tables 7-10, we conclude that our schemes are robust with respect to the Poisson ratio ν and the hydraulic conductivity K. Table 7: Errors and convergence rates given by Method 1 for Example 2 using k = 2 and l = 1 with ν = 0.49999 and K = 10 −6 .  Table 8: Errors and convergence rates given by Method 1 for Example 2 using k = 3 and l = 2 with ν = 0.49999 and K = 10 −6 .  Table 9: Errors and convergence rates given by Method 2 for Example 2 using k = 2 and l = 1 with ν = 0.49999 and K = 10 −6 .

Conclusions and outlook
In this paper, we present a priori estimates of the two fully discrete coupled schemes for the three-field formulation of Biot's consolidation model. The theoretical results show that both schemes are unconditionally convergent with optimal error orders. We comment here that Method 2 achieves a second-order convergence in time without complicating the calculations. The detailed numerical experiments are carried out to verify the predictions of error estimates. In future work, we plan to develop some decoupled algorithms [14,11] and the corresponding analysis based on the theory studied in this work.