Semi-explicit integration of second order for weakly coupled poroelasticity

We introduce a semi-explicit time-stepping scheme of second order for linear poroelasticity satisfying a weak coupling condition. Here, semi-explicit means that the system, which needs to be solved in each step, decouples and hence improves the computational efficiency. The construction and the convergence proof are based on the connection to a differential equation with two time delays, namely one and two times the step size. Numerical experiments confirm the theoretical results and indicate the applicability to higher-order schemes.


Introduction
This paper is devoted to the construction and analysis of a semi-explicit time discretization scheme of second order for linear poroelasticity [DC93,Sho00].The poroelastic equations can be characterized as a coupled system consisting of an elliptic and a parabolic equation and appear, e.g., in the field of geomechanics [Bio41,Zob10].In many applications, this coupling is rather weak, which is also a central assumption in this paper to guarantee convergence.For the temporal discretization of elliptic-parabolic problems such as poroelasticity, one mainly considers implicit schemes such as the implicit Euler method [EM09] or higher-order schemes [Fu19].This is primarily due to the fact that a semi-discretization in space yields a differential-algebraic equation for which explicit time-stepping schemes cannot be used [KM06].
Then again, one is interested in a decoupled approach, in which the elliptic and parabolic equations can be solved sequentially.Such a decoupling does not only replace the solution of a large system by two smaller subsystems to be solved but also enables the application of standard preconditioners [LMW17].Moreover, the decoupling of the systems favors a codesign paradigm, allowing the usage of highly optimized software packages for the porous media flow (the parabolic equation) and the mechanical problem (the elliptic equation) separately, and in addition, includes a linearization step if the permeability depends on the displacement, cf.[AM22].One attempt in this direction are iterative decoupling methods such as the fixed-stress, fixed-strain, or drained splitting schemes; see, e.g., [AS92,WG07,KTJ11,MW13].These schemes come along with an additional inner iteration in each time step that is required to guarantee convergence [SBK + 19] and, additionally, require a careful selection of tuning parameters.In [CR18], an alternative method based on an additional stabilization term rather than an inner iteration is proposed.To combine the advantages of the monolithic and iterative coupling methods, a semiexplicit time-stepping scheme was introduced in [AMU21], which decouples the equations and does not require an additional inner iteration or stabilization parameters.The semiexplicit scheme equals the implicit Euler discretization up to a term with a time shift in one of the equations.The perception of this scheme in terms of delay equations allows proving convergence of the method if a weak coupling condition is satisfied.This condition is independent of the step size and can be quantified explicitly.
In this paper, we extend these ideas to construct a novel higher-order decoupling time integrator for coupled elliptic-parabolic problems.To this end, we first construct a delay equation that differs from the coupled system only by a term of higher order and, afterwards, use an implicit scheme of the same order to discretize the delay equation in time.The details are presented in Section 3. As for the semi-explicit scheme introduced in [AMU21], the equations are still decoupled and no inner iteration or additional tuning parameters are required.To simplify the notation and the analysis, we focus on constructing a second-order scheme but give pointers throughout the manuscript on how this idea can be extended to higher orders.For the time discretization of the related delay equation, we use the second-order backward differentiation formula (BDF).As in the first-order case, our method depends on a weak coupling condition, which we explicitly quantify via the theory of delay differential-algebraic equations.
To prove second-order convergence of the proposed method, we consider in Section 4 a related problem, namely the convergence of the well-known BDF-2 scheme applied to an operator equation with two fixed delays.These delays equal τ and 2τ , with τ denoting the time step size of the method.Since we focus on time discretization, the convergence analysis is given on operator level, i.e., without a spatial discretization.Corresponding results for the fully discrete scheme can be obtained by the introduction of appropriate Ritz projections, cf.[ACM + 20, AMU21].We conclude our presentation with three numerical examples in Section 5.
Notation.We write a b to indicate the existence of a generic constant C, independent of spatial and temporal discretization parameters, such that a ≤ Cb.

Poroelastic Equations
In this section, we introduce the equations of linear poroelasticity and the corresponding abstract formulation as an elliptic-parabolic problem.We consider a bounded Lipschitz domain Ω ⊆ R d , d ∈ {2, 3}, in which we seek the displacement field u : [0, T ] × Ω → R d and the pore pressure p : [0, T ] × Ω → R. For a given time horizon T > 0, the system equations read Therein, σ denotes the stress tensor with Lamé coefficients λ and µ, the permeability κ, the Biot-Willis fluid-solid coupling coefficient α, the Biot modulus M , and the fluid viscosity ν; see [Bio41,Sho00].The righthand sides f and g are the volumetric load and the fluid source, respectively, modeling an injection or production process.Throughout this paper, we assume homogeneous Dirichlet boundary conditions, i.e., we set u = 0 and p = 0 on (0, T ] × ∂Ω. 2.1.Abstract formulation.For an abstract formulation of (2.1), we introduce the Hilbert spaces which include the assumed Dirichlet boundary conditions.With the respective dual spaces of V and Q denoted by V * and Q * , (V, H V , V * ) as well as (Q, H Q , Q * ) form Gelfand triples with dense embeddings; see [Zei90,Ch. 23.4] for more details.Moreover, we define the bilinear forms with the classical double dot notation and the symmetric gradient ε(u) used in continuum mechanics.With this, the weak formulation of (2.1) can be written as follows: seek for all test functions v ∈ V, q ∈ Q. Correspondingly, we assume that the right-hands sides satisfy f : [0, T ] → V * and g : [0, T ] → Q * and denote with •, • the respective dual pairings.We would like to emphasize that it is sufficient to prescribe initial data for p, since equation (2.2a) defines a consistency condition for p 0 and u 0 (which is uniquely solvable for u 0 ; see the forthcoming discussion on the properties of the bilinear forms).System (2.2) may also be written in operator form in the dual spaces of V and Q.For this, let A, B, C, and D denote the operators corresponding to the bilinear forms a, b, c, and d, respectively.Then, (2.2) is equivalent to It remains to discuss the properties of the bilinear forms.The bilinear form a : V ×V → R is symmetric, elliptic, and bounded, i.e., there exist positive constants c a , C a > 0 such that for all u, v ∈ V. We would like to emphasize that a is well-known from the theory of linear elasticity and that the ellipticity follows from Korn's inequality [Cia88,Th. 6.3.4].
Similarly, b : Q × Q → R is symmetric, elliptic, and bounded in Q, i.e., there exist positive for all p, q ∈ Q.The bilinear form c : H Q × H Q → R simply involves the multiplication by a (positive) constant and, hence, defines an inner product in the pivot space H Q .In more detail, there exist positive constants c c , C c > 0 such that c(p, p) ≥ c c p 2 HQ , c(p, q) ≤ C c p HQ q HQ for all p, q ∈ H Q .The remaining bilinear form d : V × H Q → R models the coupling and is continuous, i.e., there exists a positive constant for all u ∈ V and p ∈ H Q .
Remark 2.1.System (2.2) can also be used to model linear thermoelasticity, which consideres the displacement of a material due to temperature changes [Bio56].More generally, system (2.2) is an elliptic-parabolic system, where the elliptic part (modeled by a) and the parabolic part (modeled by b and c) are coupled through the bilinear form d.
2.2.Spatial discretization.Although this paper is mainly concerned with the temporal discretization, we shortly comment on the finite element discretization of (2.2).For more details, we refer to [EM09].In order to transfer the convergence results of this paper to the fully discrete system, one may consider spatial projection operators corresponding to the elliptic bilinear forms a and b; see [AMU21].
Considering finite-dimensional subspaces V h ⊆ V and Q h ⊆ Q, one seeks approximations u h ≈ u and p h ≈ p.Here, the parameter h represents the mesh size of the triangulation used in the construction of V h and Q h .A direct spatial discretization of (2.2) then leads to the differential-algebraic equation Therein, K a and K b denote the stiffness matrices corresponding to the bilinear forms a and b, respectively.Due to the assumptions discussed above, K a and K b can be assumed to be symmetric and positive definite.Moreover, M c equals the mass matrix corresponding to c, which is thus also symmetric and positive definite, and D is a rectangular matrix corresponding to d.
Using standard P 1 finite elements to define V h and Q h , one obtains the expected convergence rates of order one in the energy norms and order two in the L 2 -norms.For more precise results also on higher-order approximations, we again refer to [EM09].
2.3.Temporal discretization of first order.The standard way to discretize system (2.2) in time is the application of the implicit Euler scheme.This results in a timestepping scheme of order one as shown in [EM09].
As already mentioned in the introduction, the differential-algebraic structure rules out the possibility of a fully explicit discretization in time.In [AMU21], however, a semiexplicit scheme was introduced.Considering an equidistant decomposition of [0, T ] with step size τ , this scheme reads for all v ∈ V, q ∈ Q.Here, u n and p n denote the approximations of u(t n ) and p(t n ), t n = nτ , respectively.Note that, in contrast to the implicit Euler discretization, the first equation contains p n rather than p n+1 .Hence, the two equations decouple and can be solved sequentially.It is shown in [AMU21] that this maintains the first-order convergence as long as the weak coupling condition For the convergence analysis, the connection of (2.2) to a related delay system is used.This idea is also applied in the following sections to construct a semi-explicit scheme of order two, leading to a slightly more restrictive weak coupling condition.

Semi-explicit Integration Scheme of Second Order
This section is devoted to the extension of the semi-explicit scheme (2.3) to second order.Following the idea in [AMU21], we first construct a related delay equation and then discretize the delay equation with an implicit scheme of second order.For this, we need to replace the pressure in the first poroelastic equation by a time-delayed term which is second-order accurate.We first consider a Taylor expansion before we introduce discrete derivatives, leading to a system with multiple delays.
In the following, we consider a uniform partition of the time interval [0, T ] with step size τ > 0 such that N := T /τ ∈ N. Hence, we consider time points t n = nτ for n = 0, . . ., N .The approximation of a function y at time t n is then denoted by y n .
3.1.Related delay systems by Taylor expansion.In a first step, we aim to decouple the elliptic and parabolic equation in (2.2) by replacing p in the first equation by a Taylor expansion.For a first-order approximation, one can simply replace p(t) by the delay term ∆ τ p(t) := p(t − τ ), cf.Section 2.3.In the general case, we replace p(t) by the Taylor expansion of order k at time t − τ .This then leads to the delay system for test functions v ∈ V and q ∈ Q.As initial condition, we set p(0) = p(0) = p 0 .In contrast to the original system, however, we need additional information on p (and its derivatives) on the interval [−τ, 0].For this, we introduce a so-called history function The conditions on the derivatives of Φ ensure the consistency of the initial data, i.e., ū(0) = u(0) = u 0 .Since system (3.1) is constructed by the help of a Taylor expansion, it is no surprise that the solutions (u, p) and (ū, p) only differ by a term of order τ k as long as the solution of the delay systems stays stable.We refer to Appendix A for further details.Nevertheless, system (3.1) is not well-suited for the construction of a numerical scheme.This is due to the appearance of temporal derivatives, which turn out to be critical as we illustrate in the sequel.Already in the case of interest, namely k = 2, the resulting delay system (3.1) is of advanced type [BC63] and, hence, only well-posed in a distributional setting [TU19].This can be seen as follows.In operator form, the delay equation reads

The first equation yields
Inserting this in the second equation, we obtain For a given history function, this is a parabolic problem that can be solved on the interval (0, τ ).We thus can proceed iteratively to construct a solution.Nevertheless, due to the term ∆ τ p, this solution loses regularity over time, and is hence not suitable for numerical methods.We refer to [BZ03, AZ18, Ung18] for further details.
To avoid such advanced delay systems, we consider an alternative approach and replace the derivatives by discrete derivatives.This then yields a system with multiple delays.
3.2.Discrete derivatives.Based on the discrete difference operator Dy n+1 := y n+1 −y n , we can write discrete derivatives of any order in a short way.As an example, the difference quotients of order one and two read leading, e.g., to the well-known BDF schemes.For the upcoming analysis, it is convenient to extend the definition of the difference operator also to continuous functions.More precisely, we define Dy := y − ∆ τ y with the time shift ∆ τ introduced above.The resulting order-k approximations of the derivative of a function are summarized in the following lemma.
3.3.Related delay system with multiple delays.As already mentioned, we want to replace the derivatives in (3.1) by discrete derivatives.Focusing on the case k = 2, which will lead to a scheme of second order, we replace τ ∆ τ ṗ by D∆ τ p = ∆ τ p − ∆ 2τ p.This leads to the system for all test functions v ∈ V and q ∈ Q.Note that this is a system with two delays, namely τ and 2τ .Again, we need to discuss the initial data, which includes p(0) = p(0) = p 0 and an appropriate history function Φ defined on [−2τ, 0].To obtain consistency in u, namely ũ(0 Remark 3.2 (approximations of higher order).For general k ≥ 1, one possibility is to replace the derivatives ∆ τ p(j) in (3.1) by approximations of order τ k−j .This then guarantees that the resulting expression is an approximation of the Taylor expansion k−1 j=0 τ j j! ∆ τ p(j) of order k.Note, however, that this leads to a growing number of delays.For k = 3 this yields three delays, whereas k = 4 already needs five delays.The resulting scheme for k = 3 is presented in Section 5.3 Remark 3.3 (parabolic equation with multiple delays).Considering the operator formulation of (3.3) and eliminating the variable ũ by the first equation, we get Note that this is a parabolic equation (of neutral type) with two delays.Hence, we consider here multiple delays rather than higher derivatives.
Motivated by the approximation properties of the Taylor expansion approach, the following theorem shows that the solutions to (2.2) and (3.3)only differ by a term of order two.
Theorem 3.4.Assume sufficiently smooth right-hand sides f and g and a history function Φ satisfying (3.4).Then, the solutions to (2.2) and (3.3) are equal up to a term of order τ 2 , i.e., for almost all t ∈ [0, T ] we have Here, the hidden constant depends on higher derivatives of the history function Φ as well as of p.
Proof.We define e p := p − p and e u := ũ − u.Due to the assumptions on the history function, we conclude that e p (0) = 0 and e u (0) = 0. Considering the difference of (3.3a) and (2.2a), we obtain where L ∞ (−2τ, T ; H Q ) denotes the Bochner space on the time interval (−2τ, T ) with values in H Q .In the same manner, we obtain by the derivatives of (3.3a) and (2.2a) that L ∞ (−2τ,T ;HQ ) .Now we can proceed as in the proof of Proposition A.1, i.e., we consider the test function v = ėu in combination with the difference of (3.3b) and (2.2b).
Remark 3.5.The hidden constant in Theorem 3.4 may become arbitrarily large depending on the ellipticity and continuity constants.This is discussed in more detail in Section 3.5.System (3.3) yields a good starting point for the construction of higher-order discretization schemes.This is subject of the following subsection.
3.4.Semi-explicit integration scheme.In order to obtain a semi-explicit time-stepping scheme, we now apply the BDF-2 scheme to (3.3).To shorten notation, we introduce By Lemma 3.1, we know that 1 2τ BDF 2 u(t) = u(t) + O(τ 2 ).Since the first equation does not contain any derivatives, the temporal discretization is simply given by a function evaluation at time t n+2 (as for the implicit Euler scheme).This discretization yields the semi-explicit scheme for test functions v ∈ V and q ∈ Q.Note that this is a 2-step scheme, calling for initial data p 0 = p(0) and p 1 .In place of the history function, we set p The first condition corresponds to (3.4) and gives the consistency condition for u 0 .The second equation ensures that p 1 and u 1 are consistent.To be precise, this means that the resulting values u 0 , u 1 ∈ V satisfy The proposed scheme (3.6) is indeed semi-explicit, since the first equation defines u n+2 purely by already computed values, i.e., without the knowledge of p n+2 .Inserting this value in the second equation, we then obtain the approximation p n+2 .In operator form, this scheme reads Using once more the invertibility of the operator A, we can eliminate the u-variables in the second equation, leading to We would like to emphasize that this equals the BDF-2 discretization of the delay equation (3.5).This fact will be used in the following convergence result.
Theorem 3.6 (Second-order convergence of the semi-explicit scheme).Assume sufficiently smooth right-hand sides f and g.Moreover, let the operators satisfy the weak coupling condition (3.9) Then the semi-explicit scheme (3.6) converges with order two.More precisely, given p 0 = p(0) and p 1 as a second-order approximation of p(τ ), we can define consistent u 0 and u 1 in the sense of (3.8) such that for all n ≥ 0.
Proof.Given p 0 and p 1 , we define p −2 and p −1 satisfying (3.7) such that u 0 , u 1 are consistent.Moreover, let Φ be a history function with Φ(−2τ ) = p −2 and Φ(−τ ) = p −1 such that (3.4) is satisfied.We can now apply Theorem 3.4 and conclude that the exact solution and the solution of the delay system (3.3)only differ by a term of order two.Hence, it is sufficient to compare the discrete solution given by (3.6) with (p, ũ).We have seen that the presented semi-explicit scheme corresponds to the BDF-2 method applied to the delay equation (3.5).Since the operator C only contains a multiplicative factor, we may consider a simple rescaling leading to the question of the convergence of the BDF-2 scheme applied to the delay system ṗ Note that these two operators are symmetric, elliptic, and continuous in the respective spaces and that the continuity constant of C equals ω = C 2 d /(c a c c ).In Theorem 4.2 of the following section, we show that this implies an estimate of the form for n ≥ 2. Due to the assumption on p 1 , we have E init = p(τ )−p 1 2 +τ p(τ )−p 1 2 b τ 4 .Hence, the above estimate holds for all n ≥ 0.Moreover, considering the difference of equations (3.3a) and (3.6a), we get by the ellipticity of the bilinear form a that ũ

HQ
for n ≥ 0. Finally, due to the consistency conditions for u 0 and u 1 , we further get u(0) = ũ(0) = u 0 and u(τ HQ .The combination of the previous estimates completes the proof. Remark 3.7 (Initial data).In practice, appropriate initial conditions can be realized as follows: given p 0 , one first computes u 0 consistent to (2.2a).Then, p 1 and u 1 can be obtained by a single step of the implicit Euler discretization applied to (2.2).This then guarantees consistency as well as the needed accuracy for p 1 .Before we discuss the convergence of the semi-explicit scheme, we focus on the weak coupling condition (3.9) and its meaning in terms of delay equations.
3.5.Asymptotic stability of the delay system.We analyze the asymptotic stability of the related delay system constructed in Section 3.3 with multiple delays.To simplify the presentation, we consider here the finite-dimensional case after a semi-discretization in space (cf.Section 2.2), and study the neutral delay differential equation corresponding to (3.5), i.e., we study the neutral delay equation GKC03,Thm. 3.20]) for the delay-independent asymptotic stability of the unforced (i.e., gh = 0) delay equation (3.10) is that the spectral radius of the matrix is strictly less than one, i.e., ρ(N 2 ) < 1. Hereby, I denotes the identity matrix.We thus have to compute the eigenvalues of N 2 .Since M c is symmetric and positive definite, the (principle) square root M 1/2 c exists and is symmetric and positive definite.Thus, the matrix is symmetric and hence diagonalizable, i.e., there exists a diagonal matrix Λ and an orthogonal matrix U such that For any eigenvalue λ of Λ, it thus suffices to compute the spectral radius of the matrix which is given by λ + √ λ 2 + λ.Since this is a monotone expression, we conclude and thus ρ(N 2 ) < 1 if and only if ρ(M −1 c DK −1 a D T ) < 1 3 .Consequently, we cannot expect the delay equation to be a reasonable approximation of the non-delay equation if 3 is also a sufficient condition for delay-independent asymptotic stability.Using we observe that a weak coupling condition as in (3.9) is not only a technical requirement, but indeed necessary for convergence.We discuss the details in the error analysis in the next section.

Convergence Analysis
In this section, we prove the convergence of the BDF-2 method applied to the delay operator equation Here, B : Q → Q * is an operator with the same properties as B in the previous section and Q is an operator with the same properties as C with continuity constant ω.Similar as in Section 3.3, we assume, besides the initial condition z(0) = z 0 , a given history function Φ ∈ C ∞ ([−2τ, 0], Q) with Φ(0) = z 0 .Moreover, the right-hand side r : [0, T ] → Q * is sufficiently smooth.
For the error analysis, we first require the following lemma.
Lemma 4.1.For a symmetric bilinear form a it holds that a , which completes the proof.
After this preparation, we are now able to formulate the main convergence theorem.
Theorem 4.2 (Convergence of BDF-2 for the delay equation (4.1)).Let B : Q → Q * and C : H Q → H * Q be symmetric, elliptic, and continuous in the respective spaces.Moreover, let ω denote the continuity constant of C satisfying ω ≤ 1/5.Then, the BDF-2 scheme applied to (4.1), i.e., the scheme yields an approximation of second order.To be precise, assuming a sufficiently smooth right-hand side r, a step size τ ≤ 1, and initial data Q contains the initial error.Proof.Assuming sufficient regularity of the right-hand side r, we may assume r ≡ 0, since this does not influence the order of the approximation.Inserting the exact solution within the numerical scheme, we obtain the defect equation Note that, due to the assumptions on the history function and the initial data, we have e −2 = e −1 = e 0 = 0.In the following, we write • b for the norm induced by the operator B, which is equivalent to the Q-norm, and • c for the norm induced by C. Note that the latter is equivalent to the , where we use the short notation • := • HQ .
Step 1: In the first step, we derive an auxiliary estimate for differences of e n .If we multiply (4.3) by 2 and apply De n+2 , we get (4.4) 2 BDF 2 e n+2 , De n+2 Reformulating the terms T 1 and T 2 using (4.2) yields and The above computations inserted in (4.4) yield where we use the weighted Young's inequality four times with some positive constants α, β, γ, δ > 0. Rearranging terms leads to We now set α = 7/8, β = 11/2, γ = 10/11, and δ = 10/7.This leads to Assuming ω ≤ 1/5 and τ ≤ 1, we therefore get Building the sum over n, we get with e −2 = e −1 = e 0 = 0 that In particular, we obtain with E init introduced in the statement of the theorem that Step 2: For the desired estimate of the error itself, we go back to (4.3), multiply the equation by 2, and apply e n+2 .This leads to (4.5) 2 BDF 2 e n+2 , e n+2 With Lemma 4.1, we can rewrite T 1 as For the second term, we directly get T 2 = 4τ e n+2 2 b .The third term is simplified using e n+2 = e n+1 + De n+2 and once more Lemma 4.1, leading to Since the first term on the right-hand side can be once again bounded in terms of E init , this is the assertion.
In the upcoming numerical experiments, we observe that the bound ω ≤ 1/5 is close to the (optimal) critical value where the scheme gets instable.Note that a similar statement is also true for the first-order semi-explicit scheme where the theory requires ω ≤ 1, cf. [AMU21].

Numerical Experiments
This section is devoted to the numerical illustration of the convergence result presented in Theorem 3.6 and the necessity of a weak coupling condition.Moreover, we present a semi-explicit method of order three based on the above construction.5.1.Poroelastic example.In the first experiment, we investigate the convergence rates of the semi-explicit second-order scheme (3.6) and compare the results with an implicit second-order scheme based on a BDF-2 discretization.We choose Ω = (0, 1) 2 , T = 1, and set the poroelastic parameters as λ = 7.82 • 10 8 , µ = 1.826 • 10 9 , α = 0.6, M = 10 9 , κ/ν = 8 • 10 −10 .
The computations are based on a finite element implementation in FEniCS, leading to a system as described shortly in Section 2.2.We now investigate the convergence behavior of the semi-explicit scheme (3.6) and compare it with a second-order implicit BDF discretization.For the computation of a reference solution, we choose an implicit midpoint scheme with step size τ ref = 2 −11 and a spatial mesh width h ref = 2 −8 .Since we are mainly interested in the temporal discretization errors, we compute the second-order schemes for step sizes τ ∈ {2 −2 , . . ., 2 −9 } with the fixed spatial parameter h = 2 −8 .
The results are presented in Figure 5.1.Therein, we use the notion p(T ) for the reference solution and p N h for the discrete solution at time T = N τ (and accordingly for u).We observe second-order convergence for both the implicit and the semi-explicit scheme.The implicit method, however, achieves slightly better results in u compared to the semiexplicit one, whereas the results are comparable for the pressure variable p.The main advantage of the semi-explicit scheme lies in the fact that the two poroelastic equations can be solved sequentially, which results in a computational speedup.Moreover, standard preconditioners for elliptic and parabolic systems can be used.Note, however, that the semi-explicit method is only stable if an appropriate coupling condition is fulfilled as indicated in Theorem 3.6.This is further investigated in the following subsection.
5.2.Sharpness of the weak coupling condition.We now present a numerical example to investigate the requirement of the weak coupling condition in Theorem 3.6.To this end, we consider the following toy problem of the form (2.2) with The prefactor of A is chosen in such a way that c a , which equals the smallest eigenvalue of A, is exactly 1.Moreover, we have c c = 1 and for the continuity constant of d we get C d = √ ω.Therefore, we have ω = C 2 d /(c a c c ) as in Theorem 3.6.We test our semi-explicit scheme (3.6) with different step sizes τ and different coupling coefficients ω.The relative errors compared to a fine discretization with an implicit midpoint rule with step size τ ref = 2 −14 are computed at the final time T = 1/2.For the forcing functions, we choose f ≡ [ 1 1 1 ] T and g(t) = sin(t).
The corresponding results are presented in Figure 5.2.While the sufficient condition from Theorem 3.6 reads ω ≤ 1/5 and the delay approach from Section 3.5 demands ω < 1/3, we observe that the critical value for stability is roughly 0.38 and therefore slightly relaxed compared to the theoretical considerations.The experiment shows that a coupling condition as in Theorem 3.6 is indeed necessary and -up to a moderate scaling factorrather sharp.
5.3.Semi-explicit scheme of order 3.As an outlook, we go beyond the presented theory and motivate a possible extension to a semi-explicit third-order scheme.This is done by using the BDF-3 scheme for the delay system in (3.1), see the discussion in Remark 3.2.For k = 3, we have p ≈ 3p τ − 3p 2τ + p 3τ , which yields the semi-explicit 3-step scheme To illustrate the behavior in terms of the convergence rate and the weak coupling condition, we consider again the setting presented in Section 5.2.The corresponding results are shown in Figure 5.3.Note that the error decreases roughly by a factor 8 when halving the step size τ , which indicates a third-order convergence rate.As before, we observe that a suitably small coupling of the two equations is necessary in order to ensure stability.The numerically observed critical point for stability is roughly ω ≤ 1/6 and hence smaller than in the second-order case of Figure 5.2.This indicates that the required coupling condition depends on the order k of the corresponding scheme.Performing a similar analysis of the corresponding delay equation as in Section 3.5 yields that delay-independent asymptotic stability is (numerically) guaranteed for ω < 1/7.

Conclusions
Within this paper, we have constructed a semi-explicit second-order time-integration scheme for linear poroelasticity that decouples the problem and hence is suitable in a co-design paradigm where specialized legacy codes for the elliptic and parabolic equation can be used.The method is constructed by first perturbing the elastic equation with time delays, which equal multiples of the time step size, and then applying BDF-2 to this delay equation.We have proven convergence of this scheme under a suitable weak coupling 0.09 0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17  condition.This coupling condition is, as in the first-order case [AMU21], explicitly quantified via an asymptotic stability analysis of the delay equation.While our work focuses on the second-order scheme, we have demonstrated in a numerical example that the same idea can also be used to construct a third-order scheme, which however requires a more restrictive weak coupling condition as well as an alternative convergence proof.for all test functions v ∈ V, q ∈ Q.Moreover, considering the derivatives of (2.2a) and (3.

Figure 5 . 1 .
Figure 5.1.Relative errors in p (left, measured in the b-norm) and u (right, measured in the a-norm) for the poroelastic example in Section 5.1 at the final time T for fixed h = 2 −8 and varying τ .

Figure 5 . 2 .
Figure 5.2.Relative errors of the second-order semi-explicit method at the final time point T = 1/2 for different coupling parameters ω and different time step sizes τ .

Figure 5 . 3 .
Figure 5.3.Relative errors of the third-order semi-explicit method at the final time point T = 1/2 for different coupling parameters ω and different time step sizes τ .