Error Estimation of Euler Method for the Instationary Stokes–Biot Coupled Problem

In this paper, we study a finite element computational model for solving the interaction between a fluid and a poroelastic structure that couples the Stokes equations with the Biot system. Equilibrium and kinematic conditions are imposed on the interface. A mixed Darcy formulation is employed, resulting in continuity of flux condition of essential type. A Lagrange multiplier method is used to impose weakly this condition. With the obtained finite element solutions, the error estimators are performed for the fully discrete formulations.


Introduction
e paper presents a computation friendly finite element formulation of coupled fluid motions at a free-poroelastic interface, where Navier-Stokes equations are employed for free fluid and Darcy's law is used for the permeable material. A reliable and efficient a posteriori error estimator is analysed. is is a challenging multiphysics problem that has a wide range of applications, including processes arising in gas and oil extraction from naturally or hydraulically fractured reservoirs, designing industrial filters, and blood-vessel interaction. In these applications, it is important to model properly the interaction between the free fluid with the fluid within the porous medium and to take into account the effect of the deformation of the medium. For example, geomechanical effects play an important role in hydraulic fracturing, as well as in modeling phenomena such as subsidence and compaction. e free fluid region can be modeled by the Stokes or Navier-Stokes equations, while the flow through the deformable porous medium is modeled by the quasi-static Biot systems of poroelasticity [1]. e two regions are coupled via dynamic and kinematic interface conditions, including balance of forces, continuity of normal velocity, and a no slip or slip with friction tangential velocity condition [2][3][4][5][6][7][8][9][10].
ese multiphysics models exhibit features of coupled Stokes-Darcy flows and fluid-structure interaction (FSI) [11][12][13][14][15]. e well-posedness of the mathematical model based on the Stokes-Biot system for the coupling between a fluid and a poroelastic structure is studied in [16]. A numerical study of the problem, using a Navier-Stokes equations for the fluid, is presented in [11,17], utilizing a variational multiscale approach to stabilize the finite element spaces. e problem is solved using both a monolithic and a partitioned approach, with the latter requiring subiterations between the two problems.
Finite element analysis of an arbitrary Lagrangian-Eulerian method for Stokes/parabolic moving interface problem with jump coefficients has been studied in [18]. e authors in [19] studied a numerical solution of the coupled system of the time-dependent Stokes and fully dynamic Biot equations. ey established stability of the scheme and derived error estimates for the fully discrete coupled scheme. Numerical errors and convergence rates for smooth problems as well as tests on realistic material parameters have been presented. In [20], Jing Wen and Yinnian He considered a strongly conservative discretization for the rearranged Stokes-Biot model based on the interior penalty discontinuous Galerkin method and mixed finite element method. e existence and uniqueness of solution of the numerical scheme have been presented. en, the analysis of stability and priori error estimates has been derived. e numerical examples under uniform meshes which well validate the analysis of convergence and the strong mass conservation are presented. A staggered finite element procedure for the coupled Stokes-Biot system with fluid entry resistance has been studied by Bergkamp et al. in [21] while Ambartsumyan et al. in [22] studied flow and transport in fractured poroelastic media using Stokes flow in the fractures and the Biot model in the porous media. In Section 6 in [23], fully discrete continuous approximation has been proposed for the weak coupled mixed formulation. For the discretization of the fluid velocity and pressure, the authors have used the finite elements which include the MINI-elements, the Taylor-Hood elements, and the conforming Crouzeix-Raviart elements. For the discretization of the porous medium problem, they choose the spaces that include Raviart-omas and Brezzi-Douglas-Marini elements. An a priori error analysis is performed with some numerical tests confirming the convergence rates.
A posteriori error estimators are computable quantities, expressed in terms of the discrete solution and of the data that measure the actual discrete errors without the knowledge of the exact solution.
ey are essential to design adaptive mesh refinement algorithms which aqui-distribute the computational effort and optimize the approximation efficiency. Since the pioneering work of Babuška and Rheinboldt [24][25][26][27], adaptive finite element methods based on a posteriori error estimates have been extensively investigated.
In [28], we have proposed a family error indicator for semidiscrete approximation for the Stokes-Biot system. To the best of our knowledge, there is no a posteriori error estimation for the Stokes-Biot fluid-poroelastic structure interaction model for fully discrete finite element methods. Here, we develop such a posteriori error analysis for the fully discrete conforming finite element methods. We have got a new family of a local indicator error Θ K (see equation (59) in Definition 3) and global Θ (equation (60)). e difference between this paper and our previous work [28] is that our discretization is fully discrete formulation. As an advantage, the error indicators in this work are more accessible to computation. e schedule of the paper is as follows. Section 2 is devoted to notations and basic results that are used throughout the document. Our main results regarding a posteriori error analysis are stated in Section 3. We prove that our indicator errors are efficient, reliable, and then optimal.
e global inf-sup condition is the main tool yielding the reliability. In turn, the local efficiency result is derived using the technique of bubble function introduced by R. Verfürth [29] and used in similar context by C. Carstensen [30]. Finally, this paper is summarized with further works in Section 4.

Stokes-Biot Model Problem.
We consider a multiphysics model problem for free fluid's interaction with a flow in a deformable porous media, where the simulation domain Ω ⊂ R d , d � 2, 3, is a union of nonoverlapping regions Ω f and Ω p . Here, Ω f is a free fluid region with flow governed by the Stokes equations and Ω p is a poroelastic material governed by the Biot system. For simplicity of notation, we assume that each region is connected. e extension to nonconnected regions is straightforward. Let Γ fp � zΩ f ∩ zΩ p (see Figure 1).
Let (u * , p * ) be the velocity-pressure pair in Ω * , * � f, p, and let η p be the displacement in Ω p . Let μ > 0 be the fluid viscosity, let f * be the body force terms, and let q * be external source or sink terms. Let D(u f ) and σ f (u f , p f ) denote, respectively, the deformation rate tensor and the stress tensor: In the free fluid region Ω f , (u f , p f ) satisfies the Stokes equations: where T > 0 is the final time. Let σ e (η p ) and σ p (η p , p p ) be the elastic and poroelastic stress tensors, respectively: where 0 < λ min ≤ λ p (x) ≤ λ max and 0 < μ min ≤ μ p (x) ≤ μ max are the Lamé parameters and 0 < α ≤ 1 is the Biot-Willis constant. e poroelasticity region Ω p is governed by the quasistatic Biot system [23]: where s 0 ≥ 0 is a storage coefficient and K is the symmetric and uniformly positive definite rock permeability tensor, satisfying, for some constants 0 < k min ≤ k max .
∀ξ ∈ R d , Following [1], the interface conditions on the fluidporoelasticity interface Γ fp are mass conservation, balance of stresses, and the Beavers-Joseph-Saffman (BJS) condition [31] modeling slip with friction: Journal of Mathematics where n f and n p are the outward unit normal vectors to zΩ f and zΩ p , respectively, τ f,j , 1 ≤ j ≤ d − 1, is an orthogonal system of unit tangent vectors on Γ fp , K j � (Kτ f,j ) · τ f,j , and α BJS ≥ 0 is an experimentally determined friction coefficient.
Here, equations (8) and (9) represent mass conservation, equation (10) is the balance of normal forces, and equation (11) is the Beavers-Joseph-Saffman conditions. We note that continuity of flux constraints the normal velocity of the solid skeleton, while the BJS condition accounts for its tangential velocity.
e above system of equations needs to be complemented by a set of boundary and initial conditions. Let We assume that for simplicity homogeneous boundary conditions, To avoid the issue with restricting the mean value of the pressure, we assume that |Γ D p | > 0. We also assume that Γ D p is not adjacent to the interface Γ fp , i.e., dist(Γ D p , Γ fp ) ≥ s > 0. Nonhomogeneous displacement and velocity conditions can be handled in a standard way by adding suitable extensions of the boundary data. e pressure boundary condition is natural in the mixed Darcy formulation, so nonhomogeneous pressure data would lead to an additional boundary term. We further say that the initial conditions are as follows: Equations (2)-(10) consist of the model of the coupled Stokes and Biot flows problem that we will study below.

Weak Formulation.
In this part, we first introduce some Sobolev spaces [32] and norms. If W is a bounded domain of R d and m is a nonnegative integer, the Sobolev space H m (W) � W m,2 (W) is defined in the usual way with the usual norm ‖ · ‖ m,W and seminorm |.| m,W . In particular, For a connected open subset of the boundary E ⊂ zΩ f ∪ zΩ p , we write 〈., .〉 E for the L 2 (E) inner product (or duality pairing); that is, for scalar valued functions λ and σ, one defines In the following, we derive a Lagrange multiplier type weak formulation of the system, which will be the basis for our finite element approximation. Let where We define the global velocity and pressure spaces as Figure 1: Global domain Ω consisting of the fluid region Ω f and the poroelastic media region Ω p separated by the interface Γ fp .
with norms e weak formulation is obtained by multiplying the equations in each region by suitable test functions, integrating by parts the second-order terms in space, and utilizing the interface and boundary conditions. Let be the bilinear forms related to Stokes, Darcy, and the elasticity operator, respectively. Let Integration by parts in (2) and the two equations in (5) leads to the interface term as follows: Using the first condition for balance of normal stress in (9), we set which will be used as a Lagrange multiplier to impose the mass conservation interface condition (8). Utilizing the BJS condition (10) and the second condition for balance of stresses in (9), we obtain where For the well-posedness of b Γ , we require that λ ∈∧ � (V p · n p|Γ fp ) ′ . According to the normal trace theo- , (see, e.g., [33]). erefore, we take ∧ � H 1/2 (Γ fp ).
e Lagrange multiplier variational formulation is for and λ(t) ∈∧, such that p p (0) � p p,0 and η p (0) � η p,0 , and for all v f ∈ V f , w f ∈ W f , v p ∈ V p , w p ∈ W p , ξ p ∈ X p , and μ ∈∧, where we used the notation z t � (z/zt). e assumptions on the fluid viscosity μ and the material coefficients K, λ p , and μ p imply that the bilinear forms a f (·, ·), a D p (·, ·), and a e p (·, ·) are coercive and continuous in the appropriate norms. In particular, there exist positive constants c f , c p , c e , C f , C p , and C e such that where (29), (30), (33), and (34) hold true thanks to Poincaré inequality and (33) and (34) also relies on Korn's inequality. In summary, from Corollary 3.1 in [23] (p. 7), the following result holds.

Theorem 1.
ere exists a unique solution (u f , p f , u p ,

Finite Element Discretization. Let T
f h and T p h be shaperegular and quasi-uniform partition of Ω f and Ω p , respectively, both consisting of affine elements with maximal element diameter h. e two partitions may be nonmatching at the interface Γ fp . For the discretization of the fluid velocity and pressure, we choose finite element spaces V f,h ⊂ V f and W f,h ⊂ W f , which are assumed to be inf-sup stable. Examples of such spaces include the mini-elements, the Taylor-Hood elements, and the conforming Crouzeix-Raviart elements. For the discretization of the porous medium problem, we choose V p,h ⊂ V p and W p,h ⊂ W p to be any of well-known inf-sup stable mixed finite element spaces, such as the Raviart-omas or the Brezzi-Douglas-Marini spaces. e global spaces are as follows: We employ a conforming Lagrangian finite element space X p,h ⊂ X p to approximate the structure displacement. Note that the finite element spaces V f,h , V p,h , and X p,h satisfy the prescribed homogeneous boundary conditions on the external boundaries. For the discrete Lagrange multiplier space, we take e semidiscrete continuous-in-time problem reads that given p p,h (0) and η p, We will take p p,h (0) and η p,h (0) to be suitable projections of the initial data p p,0 and η p,0 .
We introduce the errors for all variables as follows: e following results hold cf. [23]. (u f,h , p f,h , u p,h , p p,h , η p,h , λ

Fully Discrete Formulation.
For the time discretization, we employ the backward Euler method. Let τ be the time step, T � Nτ, and let t n � nτ, 0 ≤ n ≤ N. Let δ τ u n ≔ (u n − u n− 1 /τ) be the first-order (backward) discrete time derivation, where u n ≔ u(t n ). en, the fully discrete model reads that given p 0 b Γ u n f,h , u n p,h , δ τ η n p,h ; μ h � 0.
We introduce the discrete-in-time norms as follows: and the errors for all variables-in-time as follows:

Error Estimation
In order to solve the Stokes-Biot model problem by efficient adaptive finite element methods, reliable and efficient a posteriori error analysis is important to provide appropriated indicators. In this section, we first define the local 6 Journal of Mathematics and global indicators (Section 3.1) and then the lower and upper error bounds are derived (Sections 3.4 and 3.5).

Residual Error Estimators.
e general philosophy of residual error estimators is to estimate an appropriate norm of the correct residual by terms that can be evaluated easier and that involve the data at hand. To this end, define the exact element residuals. e exact element residuals over a triangle or tetrahedra K ∈ T h and over E ∈ E h (Γ fp ) are defined by As it is common, these exact residuals are replaced by some finite-dimensional approximation called approximate element residual r * ,K , r * ,K , * ∈ f, p , R E,pf,l , l ∈ 1, 2, 3, 4 { }. is approximation is here achieved by projecting f f and q f on the space of piecewise constant functions in Ω f and piecewise P 1 functions in Ω p ; more precisely for all K ∈ T f h , we take While for all K ∈ T p h , we take f n p,K and q n p as the unique element of [P 1 (K)] d , respectively, P 1 (K) such that respectively, ereby, we define the approximate element residuals.
Next, introduce the gradient jump in normal direction by where I is the identity matrix of R d×d .

Journal of Mathematics
We define and en, the residual error estimator is locally defined by e global residual error estimator is given by Furthermore, denote the local approximation terms by where We set Remark 3. e residual character of each term on the righthand sides of (55)-(58) is quite clear since if U h would be the exact solution of (2)-(10), then they would vanish.

Inverse Inequalities.
In order to derive the lower error bounds, we proceed similarly as in [30,34] (see also [35]), by applying inverse inequalities and the localization technique based on simplex-bubble and face-bubble functions. To this end, we recall some notation and introduce further preliminary results. Given K ∈ T h and E ∈ E(K), we let b K and b E be the usual simplex-bubble and face-bubble functions (see (1.5) and (1.6) in [29]), respectively. In particular, b K satisfies b K ∈ P 3 (K), supp(b K )⊆K, b K � 0, on zK, and 0 ≤ b K ≤ 1, on K.
Similarly, b E ∈ P 2 (K), supp(b E )⊆ω E ≔ K ′ ∈ T h : E ∈ E(K ′ ) , b E � 0, on zK∖E and 0 ≤ b E ≤ 1, in ω E . We also recall from [36] that, given k ∈ N, there exists an extension operator L: C(E) ⟶ C(K) that satisfies L(p) ∈ P k (K) and L(p) |E � p, ∀p ∈ P k (E). A corresponding vectorial version of L, that is, the componentwise application of L, is denoted by L. Additional properties of b K , b E , and L are collected in the following lemma (see [36]).