A posteriori error analysis for a Lagrange multiplier method for a Stokes/Biot fluid-poroelastic structure interaction model

In this work we develop an a posteriori error analysis of a conforming mixed finite element method for solving the coupled problem arising in the interaction between a free fluid and a fluid in a poroelastic medium on isotropic meshes in $\mathbb{R}^d$, $d\in\{2,3\}$. The approach utilizes the semi-discrete formulation proposed by Ilona Ambartsumyan et al. in (Numerische Mathematik, October 2017). The a posteriori error estimate is based on a suitable evaluation on the residual of the finite element solution. It is proven that the a posteriori error estimate provided in this paper is both reliable and efficient. The proof of reliability makes use of suitable auxiliary problems, diverse continuous inf-sup conditions satisfied by the bilinear forms involved, Helmholtz decomposition, and local approximation properties of the Cl\'ement interpolant. On the other hand, inverse inequalities, and the localization technique based on simplexe-bubble and face-bubble functions are the main tools for proving the efficiency of the estimator. Up to minor modifications, our analysis can be extended to other finite element subspaces yielding a stable Galerkin scheme.


INTRODUCTION
In this paper, we develop an a posteriori error analysis for solving the interaction of a free incompressible viscous Newtonian fluid with a fluid within a poroelastic medium. This is a challenging multiphysics problem with applications to predicting and controlling processes arising in groundwater flow in fractured aquifers, oil and gas extraction, arterial flows, and industrial filters. 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.
We adopt the Stokes equations to model the free fluid and the Biot system [2] for the fluid in the poroelastic media. In the latter, the volumetric deformation of the elastic porous matrix is complemented with the Darcy equation that describes the average velocity of the fluid in the pores. The model features two different kinds of coupling across the interface: Stokes-Darcy coupling [3,4,5,6,7,8,9,10,11] and fluid-structure interaction (FSI) [12,13,14,15,16].
The 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 [17]. A numerical study of the problem, using a Navier-Stokes equations for the fluid, is presented in [12,18], utilizing a variational multiscale approach to stabilized the finite element spaces. The problem is solved using both a monolithic and a partitioned approach, with the latter requiring subiterations between the two problems.
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. They are essential to design adaptive mesh refinement algorithms which aquidistribute the computational effort and optimize the approximation efficiency. Since the pioneering work of Babuška and Rheinboldt [19,20,21,22], adaptive finite element methods based on a posteriori error estimates have been extensively investigated.
In [1], semidiscrete continuous-in-time 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-Thomas and Douglas-Marini elements. An a priori error analysis is performed with some numerical tests confirming the convergence rates. To our best knowledge, there is no a posteriori error estimation for the Stokes/Biot fluid-poroelastic structure interaction model for finite element methods. Here we develop such a posteriori error analysis for the semi-discrete conforming finite element methods. We have got a new family of a local indicator error Θ K (see Definition 3.3, eq. (3.15)) and global Θ (eq. 3.19). We prove that our indicators error are efficiency and reliability, and then, are optimal. The 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 [23] and used in similar context by C. Carstensen [24].
The paper is organized as follows. Some preliminaries and notation are given in Section 2. In Section 3, the a posteriori error estimates are derived. We offer our conclusion and the further works in Section 4.

PRELIMINARIES AND NOTATIONS
2.1. 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 non-overlapping 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. The extension to non-connected regions is straightforward. Let Γ f p = ∂Ω f ∩ ∂Ω p (see Fig. 1). Figure 1: Global domain Ω consisting of the fluid region Ω f and the poroelastic media region Ω p separated by the interface Γ f p .
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 ) satisfy 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-Wallis constant. The poroelasticity region Ω p is governed by the quasi-static Biot system [1]: where s 0 ≥ 0 is a storage coefficient and K the symmetric and uniformly positive definite rock permeability tensor, satisfying, for some constants 0 < k min ≤ k max , Following [2], the interface conditions on the fluid-poroelasticity interface Γ f p are mass conservation, balance of stresses, and the Beavers-Joseph-Saffman (BJS) condition [25] modeling slip with friction: where n f and n p are the outward unit normal vectors to ∂Ω f , and ∂Ω p , respectively, τ f,j , 1 ≤ j ≤ d − 1, is an orthogonal system of unit tangent vectors on Γ f p , K j = (Kτ f,j ) · τ f,j , and α BJS ≥ 0 is an experimentally determined friction coefficient. We note that continuity of flux constraints the normal velocity of the solid skeleton, while the BJS condition accounts for its tangential velocity.
The above system of equations needs to be complemented by a set of boundary and initial conditions. Let We assume for simplicity homogeneous boundary conditions: To ovoid 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 Γ f p , i.e., dist(Γ D p , Γ f p ) ≥ s > 0. Non-homogeneous displacement and velocity conditions can be handled in a standard way by adding suitable extensions of the boundary data. The pressure boundary condition is natural in the mixed Darcy formulation, so non-homogeneous pressure data would lead to an additional boundary term. We further sat the initial conditions: 2.2. Weak formulation. In this part, we first introduce some Sobolev spaces [26] and norms. If W is a bounded domain of R d and m is a non negative integer, the Sobolev space H m (W ) = W m,2 (W ) is defined in the usual way with the usual norm · m,W and semi-norm |.| m,W . In particular, H 0 (W ) = L 2 (W ) and we write · W for · 0,W . Similarly we denote by (·, ·) W the L 2 (W ) [L 2 (W )] d or [L 2 (W )] d×d inner product. For shortness if W is equal to Ω, we will drop the index Ω, while for any m ≥ 0, · m, = · m,Ω , |.
For a connected open subset of the boundary E ⊂ ∂Ω f ∪ ∂Ω p , we write ., . E for the L 2 (E) inner product (or duality pairing), that is, for scalar valued functions λ, σ 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 Ωp . We define the global velocity and pressure spaces as Ωp . The weak formulation is obtained by mutiplying 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 define Ωp be the bilinear forms related to Stokes, Darcy and the elasticity operator, respectively. Let

Integration by parts in (2.1) and the two equations in (2.4) lead to the interface term
Using the first condition for balance of normal stress in (2.7) we set which will be used as a Lagrange multiplier to impose the mass conservation interface condition (2.6). Utilizing the BJS condition (2.8) and the second condition for balance of stresses in (2.7), we obtain The Lagrange multiplier variational formulation is: Where we used the notation ∂ t = ∂ ∂t . The 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 , C e such that: h and T p h be shape-regular and quasi-uniform partition of Ω f and Ω p , respectively, both consisting of affine elements with maximal element diameter h. The two partitions may be non-matching at the interface Γ f p . 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-Thomas or the Brezzi-Douglas-Marini spaces. The global spaces are: 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 The semi-discrete continuous-in-time problem reads: 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: e f := u f −u f,h , e p := u p −u p,h , e s := η p −η p,h , e f p := p f −p f,h , e pp := p p −p p,h and e λ := λ−λ h .

A-POSTERIORI ERROR ANALYSIS
A 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 and global indicators (Section 3.1) and then the lower and upper error bounds are derived (Sections (3.3) and (3.4)).

3.1.
Residual error estimators. The 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: be an arbitrary finite element function. The exact element residuals over a triangle or tetrahedra K ∈ T h and over E ∈ E h (Γ f p ) are defined for all t ∈]0, T ] 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}. This 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 p,K and q p as the unique element of [P 1 (K)] d respectively P 1 (K) such that: Thereby, we define the approximate element residuals.

Definition 3.2 (Approximate Element Residuals). Let t ∈]0, T ] and
be an arbitrary finite element function. Then, the approximate element residuals are defined for all t ∈]0, T ] by: Next, introduce the gradient jump in normal direction by where I is the identity matrix of R d×d .
Then, the residual error estimator is locally defined by and .
The global residual error estimator is given by Furthermore denote the local approximation terms by The global approximation term is defined by Remark 3.1. The residual character of each term on the right-hand sides of (3.15)-(3.18) is quite clear since if U h would be the exact solution of (2.1)-2.8, then they would vanish.

Inverse inequalities.
In order to derive the lower error bounds, we proceed similarly as in [24] and [28] (see also [29]), 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 simplexe-bubble and face-bubble functions respectively (see (1.5) and (1.6) in [23]). In We also recall from [30] 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 [30]): Lemma 3.1. Given k ∈ N * , there exist positive constants depending only on k and shaperegularity of the triangulations (minimum angle condition), such that for each simplexe K and E ∈ E(K) there hold  In addition, we will make use of a vector valued version of I 0 Cl , that is, which is defined componentwise by I 0 Cl . The following lemma establishes the local approximation properties of I 0 Cl (and hence of I 0 Cl ), for a proof see [31,Section 3].
The first main result is given by the following theorem: Then, there exist a positive constant C rel such that the error is bounded globally from above by: Then the continuous problem (2.14)-(2.16) is equivalent to: Find U ∈ H such that for all t ∈]0, T ], we have: We define the discrete version by the same way: Find U h ∈ H h such that for all t ∈]0, T ], Since for all t ∈]0, T ], and W h ∈ H h A(U(t) − U h (t), W h (t)) = 0, then from (3.31) we obtain Hence, with, .
and β p,h (t) = I Cl (β p (t)). Thus v p (t) − v p,h (t) = (w p (t) − w p,h (t)) + curl(β p (t) − β p,h (t)). Therefore, integrate by parts element by element we may write: Coercivity of operator A leads to inf-sup condition: By consequently , the identity (3.33), inf-sup condition of operator A (3.34), Cauchy-Schwarz inequality, estimation of Lemma (3.4) and the approximation properties of Lemma 3.3 imply the required estimate and finish the proof.
3.4. Efficiency of the a posteriori error estimator. For K ∈ T h , we set The error estimator Θ(U h ) is consider efficient if it satisfies the following theorem: Theorem 3.6 (Lower Error Bound). Let U = (u f , p f , u p , p p , η p , λ) ∈ H be the exat solution and U h = (u f,h , p f,h , u p,h , p p,h , η p,h , λ h ) ∈ H h be the finite element solution. Then, there exist a positive constant C eff such that the error is bounded locally from below for all K ∈ T h by: wherew K is a finite union of neighborning elements of K.
Proof. We begin by bounding each the residuals separately.

DISCUSSION
In this paper we have discussed a posteriori error estimates for a finite element approximation of the Stokes-Biot system. A residual type a posteriori error estimator is provided, that is both reliable and efficient. Many issues remain to be addressed in this area, let us mention other types of a posteriori error estimators and nonconforming finite element methods or implementation and convergence analysis of adaptive finite element methods. Further it is well known that an internal layer appears at the interface Γ f p as the permeability tensor degenerates, in that case anisotropic meshes have to be used in this layer (see for instance [11]). Hence we intend to extend our results to such anisotropic meshes.

NOMENCLATURES
• Ω ⊂ R d , d ∈ {2, 3} bounded domain • Ω p : the poroelastic medium domain • Ω f = Ω Ω d • Γ f p = ∂Ω f ∩ ∂Ω p • Γ = ∂Ω Γ f p , = f, p • n f (resp. n p ) the unit outward normal vector along ∂Ω f (resp. ∂Ω p ) • u f : the fluid velocity in Ω f • p f : the fluid pressure in Ω f • u p , η p : the fluid velocities in Ω p • p p : the fluid pressure in Ω p • In 2D, the curl of a scalar function w is given as usual by curl w := ∂w ∂x 2 , − ∂w ∂x 1 • In 3D, the curl of a vector function w = (w 1 , w 2 , w 3 ) is given as usual by curl w := ∇ × w namely, • P k : the space of polynomials of total degree not larger than k • T h : triangulation of Ω • T h : the corresponding induced triangulation of Ω , ∈ {f, p} • For any K ∈ T h , h K is the diameter of K and ρ K = 2r K is the diameter of the largest ball inscribed into K • h := max