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

Authors: Koffi Wilfrid Houédanou; Jamal Adetola
Journal: Journal of Mathematics
 (2021)

---

## Abstract

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.


---

## Body

## 1. Introduction

The 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. This 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.The 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]. The 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–10]. These multiphysics models exhibit features of coupled Stokes–Darcy flows and fluid-structure interaction (FSI) [11–15].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 [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. The 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]. The authors in [19] studied a numerical solution of the coupled system of the time-dependent Stokes and fully dynamic Biot equations. They 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. The existence and uniqueness of solution of the numerical scheme have been presented. Then, the analysis of stability and priori error estimates has been derived. The 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–Thomas 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. They 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–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)). The 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.The schedule of the paper is as follows. Section2 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. 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 [29] and used in similar context by C. Carstensen [30]. Finally, this paper is summarized with further works in Section 4.


## 2. 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Ω⊂ℝ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. The extension to nonconnected regions is straightforward. Let Γfp=∂Ωf∩∂Ωp (see Figure 1).Figure 1
Global domainΩ consisting of the fluid region Ωf and the poroelastic media region Ωp separated by the interface Γfp.Letu∗,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 Duf and σfuf,pf denote, respectively, the deformation rate tensor and the stress tensor:(1)Duf=12∇uf+∇ufT,σfuf,pf=−pfI+2μDuf.In the free fluid regionΩf, uf,pf satisfies the Stokes equations:(2)−∇⋅σfuf,pf=ff,inΩf×0,T,(3)∇⋅uf=qf,inΩf×0,T,where T>0 is the final time. Let σeηp and σpηp,pp be the elastic and poroelastic stress tensors, respectively:(4)σeηp=λp∇⋅ηpI+2μpDηp,σpηp,pp=σeηp−αppI,where 0<λmin≤λpx≤λmax and 0<μmin≤μpx≤μmax are the Lamé parameters and 0<α≤1 is the Biot–Willis constant. The poroelasticity region Ωp is governed by the quasi-static Biot system [23]:(5)−∇⋅σpηp,pp=fp,μK−1up+∇pp=0,inΩp×0,T,(6)∂∂ts0pp+α∇⋅ηp+∇⋅up=qp,inΩp×0,T,where s0≥0 is a storage coefficient and K is the symmetric and uniformly positive definite rock permeability tensor, satisfying, for some constants 0<kmin≤kmax.(7)∀ξ∈ℝd,kminξTξ≤ξTKxξ≤kmaxξTξ,∀x∈Ωp.Following [1], the interface conditions on the fluid-poroelasticity interface Γfp are mass conservation, balance of stresses, and the Beavers–Joseph–Saffman (BJS) condition [31] modeling slip with friction:(8)uf⋅nf+∂ηp∂t+up⋅np=0,onΓfp×0,T,(9)σfnf+σpnp=0,onΓfp×0,T,(10)−σfnf⋅nf=pp,onΓfp×0,T,(11)−σfnf⋅τf,j=μαBJFKj−1uf−∂ηp∂t⋅τf,j,onΓfp×0,T,where nf and np 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 Γfp, Kj=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.The above system of equations needs to be complemented by a set of boundary and initial conditions. LetΓf=∂Ωf∩∂Ω and Γp=∂Ωp∩∂Ω. Let Γp=ΓpD∪ΓpN. We assume that for simplicity homogeneous boundary conditions,(12)uf=0,onΓf×0,T,ηp=0,onΓp×0,T,pp=0,onΓpD×0,T,up⋅np=0,onΓpN×0,T.To avoid the issue with restricting the mean value of the pressure, we assume thatΓpD>0. We also assume that ΓpD is not adjacent to the interface Γfp, i.e., distΓpD,Γfp≥s>0. Nonhomogeneous 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 nonhomogeneous pressure data would lead to an additional boundary term. We further say that the initial conditions are as follows:(13)ppx,0=pp,0x,ηpx,0=ηp,0x,inΩp.Equations (2)–(10) consist of the model of the coupled Stokes and Biot flows problem that we will study below.


### 2.2. Weak Formulation

In this part, we first introduce some Sobolev spaces [32] and norms. If W is a bounded domain of ℝd and m is a nonnegative integer, the Sobolev space HmW=Wm,2W is defined in the usual way with the usual norm ⋅m,W and seminorm .m,W. In particular, H0W=L2W, and we write ⋅W for ⋅0,W. Similarly, we denote by ⋅,⋅W the L2WL2Wd or L2Wd×d inner product. For shortness, if W is equal to Ω, we will drop the index Ω, while for any m≥0, ⋅m,∗=⋅m,Ω∗, .m,∗=.m,Ω∗, and .,.∗=⋅,⋅Ω∗, for ∗=f,p. The space H0mΩ denotes the closure of C0∞Ω in HmΩ. Let HmΩd be the space of vector valued functions v=v1,…,vd with components vi in HmΩ. The norm and the seminorm on HmΩd are given by(14)vm,Ω≔∑i=0dvim,Ω21/2,vm,Ω≔∑i=0dvim,Ω21/2.For a connected open subset of the boundaryE⊂∂Ωf∪∂Ωp, we write .,.E for the L2E inner product (or duality pairing); that is, for scalar valued functions λ and σ, one defines(15)λ,σE≔∫Eλσds.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(16)Vf=vf∈H1Ωfd:vf=0onΓf,Wf=L2Ωf,Vp=vp∈Hdiv;Ωp:vp⋅np=0onΓpN,Wp=L2Ωp,Xp=ξp∈H1Ωpd:ξp=0onΓp,where Hdiv;Ωp is the space of L2Ωpd-vectors with divergence in L2Ωp with a norm(17)vHdiv;Ωp2≔vΩp2+∇⋅vΩp2.We define the global velocity and pressure spaces as(18)V=v=vf,vp∈Vf×Vp,W=w=wf,wp∈Wf×Wp,with norms(19)vV2≔vfΩf2+vpHdiv;Ωp2,wW2=wfΩf2+wpΩp2.The 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(20)afuf,vf≔2μDuf,DvfΩf,apDup,vp≔μK−1up,vpΩp,apeηp,ξp≔2μpDηp,DξpΩp+λp∇⋅ηp,∇⋅ξpΩp,be the bilinear forms related to Stokes, Darcy, and the elasticity operator, respectively. Let(21)b∗v,w≔−∇⋅v,wΩ∗.Integration by parts in (2) and the two equations in (5) leads to the interface term as follows:(22)IΓfp=−σfnf,vfΓfp−σpnp,ξpΓfp+pp,vp⋅npΓfp.Using the first condition for balance of normal stress in (9), we set(23)λ=−σfnf⋅nf=pp,onΓfp,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(24)IΓfp=aBJSuf,∂tηp;vf,ξp+bΓvf,vp,ξp,λ,where(25)aBJSuf,ηp;vf,ξp=∑j=1d−1μαBJSKj−1uf−ηp⋅τf,j,vf−ξp⋅τf,jΓfp,bΓvf,vp,ξp;μ=vf⋅nf+ξp+vp⋅np,μΓfp.For the well-posedness ofbΓ, we require that λ∈∧=Vp⋅np|Γfp′. According to the normal trace theorem, since vp∈Vp⊂Hdiv;Ωp, then vp⋅np∈H−1/2∂Ωp. Furthermore, since vp⋅np=0 on ΓpN and disΓpD,Γfp≥s>0, then vp⋅np∈H−1/2Γfp, (see, e.g., [33]). Therefore, we take ∧=H1/2Γfp.The Lagrange multiplier variational formulation is fort∈0,T, find uft∈Vf, pft∈Wf, upt∈Vp, ppt∈Wp, ηpt∈Xp, and λt∈∧, such that pp0=pp,0 and ηp0=ηp,0, and for all vf∈Vf, wf∈Wf, vp∈Vp, wp∈Wp, ξp∈Xp, and μ∈∧,(26)afuf,vf+apDup,vp+apeηp,ξp+aBJSuf,∂tηp;vf,ξp+bfvf,pf+bpvp,pp+αbpξp,pp+bΓvf,vp,ξp;λ=ff,vfΩf+fp,ξpΩp,(27)s0∂tpp,wpΩp−αbp∂tηp,wp−bpup,wp−bfuf,wf=qf,wfΩf+qp,wpΩp,(28)bΓuf,up,∂tηp;μ=0,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 af⋅,⋅, apD⋅,⋅, and ape⋅,⋅ are coercive and continuous in the appropriate norms. In particular, there exist positive constants cf, cp, ce, Cf, Cp, and Ce such that(29)cfvfH1Ωf2≤afvf,vf,∀vf∈Vf,(30)afvf,qf≤CfvfH1ΩfqfH1Ωf,∀vf,qf∈Vf,(31)cpvpL2Ωp2≤apDvp,vp,∀vp∈Vp,(32)apDvp,qp≤CpvpL2ΩpqpL2Ωp,∀vp,qp∈Vp,(33)ceξpH1Ωp2≤apeξp,ξp,∀ξp∈Xp,(34)apeξp,ζp≤CeξpH1ΩpζpH1Ωp,∀ξp,ζp∈Xp,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.
There exists a unique solutionuf,pf,up,pp,ηp,λ∈L∞0,T;Vf×L∞0,T;Wf×L∞0,T;Vp×W1,∞0,T;Wp×W1,∞0,T;Xp×L∞0,T;Λ to problems (26)–(28).


### 2.3. Finite Element Discretization

LetThf and Thp 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 nonmatching at the interface Γfp. For the discretization of the fluid velocity and pressure, we choose finite element spaces Vf,h⊂Vf and Wf,h⊂Wf, 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 Vp,h⊂Vp and Wp,h⊂Wp 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 as follows:(35)Vh≔vh=vf,h,vp,h∈Vf,h×Vp,h,Wh=wh=wf,h,wp,h∈Wf,h×Wp,h.We employ a conforming Lagrangian finite element spaceXp,h⊂Xp to approximate the structure displacement. Note that the finite element spaces Vf,h, Vp,h, and Xp,h satisfy the prescribed homogeneous boundary conditions on the external boundaries. For the discrete Lagrange multiplier space, we take(36)Λh=Vp,h⋅np|Γfp.The semidiscrete continuous-in-time problem reads that givenpp,h0 and ηp,h0, for t∈0,T, find uf,ht∈Vf,h, pf,ht∈Wf,h, up,ht∈Vp,h, pp,ht∈Wp,h, ηp,ht∈Xp,h, and λht∈Λh such that for all vf,h∈Vf,h, wf,h∈Wf,h, vp,h∈Vp,h, wp,h∈Wp,h, ξp,h∈Xp,h, and μh∈Λh,(37)afuf,h,vf,h+apDup,h,vp,h+apeηp,h,ξp,h+aBJSuf,h,∂tηp,h;vf,h,ξp,h+bfvf,h,pf,h+bpvp,h,pp,h+αbpξp,h,pp,h+bΓvf,h,vp,h,ξp,h;λh=ff,h,vf,hΩf+fp,h,ξp,hΩp,(38)s0∂tpp,h,wp,hΩp−αbp∂tηp,h,wp,h−bpup,h,wp,h−bfuf,h,wf,h=qf,wf,hΩf+qp,wp,hΩp,(39)bΓuf,h,up,h,∂tηp,h;μh=0.We will takepp,h0 and ηp,h0 to be suitable projections of the initial data pp,0 and ηp,0.We introduce the errors for all variables as follows:(40)ef≔uf−uf,h,ep≔up−up,h,es≔ηp−ηp,h,efp≔pf−pf,h,epp≔pp−pp,heλ≔λ−λh.The following results hold cf. [23].Theorem 2 (a priori error estimation).
There exists a unique solutionuf,h,pf,h,up,h,pp,h,ηp,h,λh in L∞0,T;Vf,h×L∞0,T;Wf,h×L∞0,T;Vp,h×W1,∞0,T;Wp,h×W1,∞0,T;Xp,h×L∞0,T;Λh of the discrete problems (37)–(39), and if the solution uf,pf,up,pp,ηp,λ∈Vf×Wf×Vp×Wp×Xp×Λ of the continuous problems (26)–(28) is smooth enough, then we have(41)efL20,T;H1Ωf+epL20,T;L2Ωp+esL∞0,T;H1Ωp+epL∞0,T;L2Ωp+efL20,T;L2Ωf+eλL20,T;Λh≤Cuf,pf,up,pp,ηp,λhr,r>0.Remark 1.
ForWh=vf,h,wf,h,vp,h,wp,h,ξp,h,λh∈L∞0,T;Vf,h×L∞0,T;Wf,h×L∞0,T;Vp,h×W1,∞0,T;Wp,h×W1,∞0,T;Xp,h×L∞0,T;Λh, we can subtract (37)–(39) to (26)–(28) to obtain the Galerkin orthogonality relation for all t∈0,T:(42)afeft,vf,ht+apDest,ξp,ht+aBJSeft,∂test;vf,ht,ξp,ht+bfvf,ht,eft+bpvp,ht,ept+αbpξpt,ept+bΓvft,vpt,ξpt,eλt+s0∂tept,wp,ht−αbp∂ttest,wp,ht−bpept,wp,ht−bfeft,wf,ht=0.


### 2.4. Fully Discrete Formulation

For the time discretization, we employ the backward Euler method. Letτ be the time step, T=Nτ, and let tn=nτ, 0≤n≤N. Let δτun≔un−un−1/τ be the first-order (backward) discrete time derivation, where un≔utn. Then, the fully discrete model reads that given pp,h0=pp,h0 and ηp,h0=ηp,h0, find uf,hn∈Vf,h, pf,hn∈Wf,h, up,h∈Vp,h, pp,hn∈Wp,h, ηp,hn∈Xp,h, and λhn∈Λh, 1≤n≤N, such that for all vf,h∈Vf,h, wf,h∈Wf,h, vp,h∈Vp,h, wp,h∈Wp,h, ξp,h∈Xp,h, and μh∈Λh,(43)afuf,hn,vf,h+apDup,hn,vp,h+apeηp,hn,ξp,h+aBJSuf,hn,δτηp,hn;vf,h,ξp,h+bfvf,h,pf,hn+bpvp,h,pp,hn+αbpξp,h,pp,hn+bΓvf,h,vp,h,ξp,h;λhn=ffn,vf,hΩf+fpn,ξp,hΩp,(44)s0δτpp,h,wp,hΩp−αbpδτηp,hn,wp,h−bpup,hn,wp,h−bfuf,hn,wf,h=qfn,wf,hΩf+qpn,wp,hΩp,(45)bΓuf,hn,up,hn,δτηp,hn;μh=0.We introduce the discrete-in-time norms as follows:(46)ϕl20,T;X2≔τ∑n=0NϕnX21/2,ϕl∞0,T;X2≔max0≤n≤NϕnX,and the errors for all variables-in-time as follows:(47)efn≔uf−uf,hn,epn≔up−up,hn,esn≔ηp−ηp,hn,efpn≔pf−pf,hn,eppn≔pp−pp,hn,eλn≔λ−λhn.Remark 2.
Forn∈0,…,N and Whn=vf,hn,wf,hn,vp,hn,wp,hn,ξp,hn,μhn∈Vf,h×Wf,h×Vp,h×Wp,h×Xp,h×Λh, we can subtract (37)–(39) to (43) and (44) to obtain the Galerkin orthogonality relation:(48)afefn,vf,hn+apDesn,ξp,hn+aBJSefn,δτesn;vf,hn,ξp,hn+bfvf,hn,efn+bpvp,hn,epn+αbpξpn,epn+bΓvfn,vpn,ξpn,eλn+s0δτepn,wp,hn−αbpδτesn,wp,hn−bpepn,wp,hn−bfefn,wf,hn=0.


## 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Ω⊂ℝ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. The extension to nonconnected regions is straightforward. Let Γfp=∂Ωf∩∂Ωp (see Figure 1).Figure 1
Global domainΩ consisting of the fluid region Ωf and the poroelastic media region Ωp separated by the interface Γfp.Letu∗,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 Duf and σfuf,pf denote, respectively, the deformation rate tensor and the stress tensor:(1)Duf=12∇uf+∇ufT,σfuf,pf=−pfI+2μDuf.In the free fluid regionΩf, uf,pf satisfies the Stokes equations:(2)−∇⋅σfuf,pf=ff,inΩf×0,T,(3)∇⋅uf=qf,inΩf×0,T,where T>0 is the final time. Let σeηp and σpηp,pp be the elastic and poroelastic stress tensors, respectively:(4)σeηp=λp∇⋅ηpI+2μpDηp,σpηp,pp=σeηp−αppI,where 0<λmin≤λpx≤λmax and 0<μmin≤μpx≤μmax are the Lamé parameters and 0<α≤1 is the Biot–Willis constant. The poroelasticity region Ωp is governed by the quasi-static Biot system [23]:(5)−∇⋅σpηp,pp=fp,μK−1up+∇pp=0,inΩp×0,T,(6)∂∂ts0pp+α∇⋅ηp+∇⋅up=qp,inΩp×0,T,where s0≥0 is a storage coefficient and K is the symmetric and uniformly positive definite rock permeability tensor, satisfying, for some constants 0<kmin≤kmax.(7)∀ξ∈ℝd,kminξTξ≤ξTKxξ≤kmaxξTξ,∀x∈Ωp.Following [1], the interface conditions on the fluid-poroelasticity interface Γfp are mass conservation, balance of stresses, and the Beavers–Joseph–Saffman (BJS) condition [31] modeling slip with friction:(8)uf⋅nf+∂ηp∂t+up⋅np=0,onΓfp×0,T,(9)σfnf+σpnp=0,onΓfp×0,T,(10)−σfnf⋅nf=pp,onΓfp×0,T,(11)−σfnf⋅τf,j=μαBJFKj−1uf−∂ηp∂t⋅τf,j,onΓfp×0,T,where nf and np 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 Γfp, Kj=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.The above system of equations needs to be complemented by a set of boundary and initial conditions. LetΓf=∂Ωf∩∂Ω and Γp=∂Ωp∩∂Ω. Let Γp=ΓpD∪ΓpN. We assume that for simplicity homogeneous boundary conditions,(12)uf=0,onΓf×0,T,ηp=0,onΓp×0,T,pp=0,onΓpD×0,T,up⋅np=0,onΓpN×0,T.To avoid the issue with restricting the mean value of the pressure, we assume thatΓpD>0. We also assume that ΓpD is not adjacent to the interface Γfp, i.e., distΓpD,Γfp≥s>0. Nonhomogeneous 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 nonhomogeneous pressure data would lead to an additional boundary term. We further say that the initial conditions are as follows:(13)ppx,0=pp,0x,ηpx,0=ηp,0x,inΩp.Equations (2)–(10) consist of the model of the coupled Stokes and Biot flows problem that we will study below.


## 2.2. Weak Formulation

In this part, we first introduce some Sobolev spaces [32] and norms. If W is a bounded domain of ℝd and m is a nonnegative integer, the Sobolev space HmW=Wm,2W is defined in the usual way with the usual norm ⋅m,W and seminorm .m,W. In particular, H0W=L2W, and we write ⋅W for ⋅0,W. Similarly, we denote by ⋅,⋅W the L2WL2Wd or L2Wd×d inner product. For shortness, if W is equal to Ω, we will drop the index Ω, while for any m≥0, ⋅m,∗=⋅m,Ω∗, .m,∗=.m,Ω∗, and .,.∗=⋅,⋅Ω∗, for ∗=f,p. The space H0mΩ denotes the closure of C0∞Ω in HmΩ. Let HmΩd be the space of vector valued functions v=v1,…,vd with components vi in HmΩ. The norm and the seminorm on HmΩd are given by(14)vm,Ω≔∑i=0dvim,Ω21/2,vm,Ω≔∑i=0dvim,Ω21/2.For a connected open subset of the boundaryE⊂∂Ωf∪∂Ωp, we write .,.E for the L2E inner product (or duality pairing); that is, for scalar valued functions λ and σ, one defines(15)λ,σE≔∫Eλσds.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(16)Vf=vf∈H1Ωfd:vf=0onΓf,Wf=L2Ωf,Vp=vp∈Hdiv;Ωp:vp⋅np=0onΓpN,Wp=L2Ωp,Xp=ξp∈H1Ωpd:ξp=0onΓp,where Hdiv;Ωp is the space of L2Ωpd-vectors with divergence in L2Ωp with a norm(17)vHdiv;Ωp2≔vΩp2+∇⋅vΩp2.We define the global velocity and pressure spaces as(18)V=v=vf,vp∈Vf×Vp,W=w=wf,wp∈Wf×Wp,with norms(19)vV2≔vfΩf2+vpHdiv;Ωp2,wW2=wfΩf2+wpΩp2.The 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(20)afuf,vf≔2μDuf,DvfΩf,apDup,vp≔μK−1up,vpΩp,apeηp,ξp≔2μpDηp,DξpΩp+λp∇⋅ηp,∇⋅ξpΩp,be the bilinear forms related to Stokes, Darcy, and the elasticity operator, respectively. Let(21)b∗v,w≔−∇⋅v,wΩ∗.Integration by parts in (2) and the two equations in (5) leads to the interface term as follows:(22)IΓfp=−σfnf,vfΓfp−σpnp,ξpΓfp+pp,vp⋅npΓfp.Using the first condition for balance of normal stress in (9), we set(23)λ=−σfnf⋅nf=pp,onΓfp,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(24)IΓfp=aBJSuf,∂tηp;vf,ξp+bΓvf,vp,ξp,λ,where(25)aBJSuf,ηp;vf,ξp=∑j=1d−1μαBJSKj−1uf−ηp⋅τf,j,vf−ξp⋅τf,jΓfp,bΓvf,vp,ξp;μ=vf⋅nf+ξp+vp⋅np,μΓfp.For the well-posedness ofbΓ, we require that λ∈∧=Vp⋅np|Γfp′. According to the normal trace theorem, since vp∈Vp⊂Hdiv;Ωp, then vp⋅np∈H−1/2∂Ωp. Furthermore, since vp⋅np=0 on ΓpN and disΓpD,Γfp≥s>0, then vp⋅np∈H−1/2Γfp, (see, e.g., [33]). Therefore, we take ∧=H1/2Γfp.The Lagrange multiplier variational formulation is fort∈0,T, find uft∈Vf, pft∈Wf, upt∈Vp, ppt∈Wp, ηpt∈Xp, and λt∈∧, such that pp0=pp,0 and ηp0=ηp,0, and for all vf∈Vf, wf∈Wf, vp∈Vp, wp∈Wp, ξp∈Xp, and μ∈∧,(26)afuf,vf+apDup,vp+apeηp,ξp+aBJSuf,∂tηp;vf,ξp+bfvf,pf+bpvp,pp+αbpξp,pp+bΓvf,vp,ξp;λ=ff,vfΩf+fp,ξpΩp,(27)s0∂tpp,wpΩp−αbp∂tηp,wp−bpup,wp−bfuf,wf=qf,wfΩf+qp,wpΩp,(28)bΓuf,up,∂tηp;μ=0,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 af⋅,⋅, apD⋅,⋅, and ape⋅,⋅ are coercive and continuous in the appropriate norms. In particular, there exist positive constants cf, cp, ce, Cf, Cp, and Ce such that(29)cfvfH1Ωf2≤afvf,vf,∀vf∈Vf,(30)afvf,qf≤CfvfH1ΩfqfH1Ωf,∀vf,qf∈Vf,(31)cpvpL2Ωp2≤apDvp,vp,∀vp∈Vp,(32)apDvp,qp≤CpvpL2ΩpqpL2Ωp,∀vp,qp∈Vp,(33)ceξpH1Ωp2≤apeξp,ξp,∀ξp∈Xp,(34)apeξp,ζp≤CeξpH1ΩpζpH1Ωp,∀ξp,ζp∈Xp,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.
There exists a unique solutionuf,pf,up,pp,ηp,λ∈L∞0,T;Vf×L∞0,T;Wf×L∞0,T;Vp×W1,∞0,T;Wp×W1,∞0,T;Xp×L∞0,T;Λ to problems (26)–(28).


## 2.3. Finite Element Discretization

LetThf and Thp 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 nonmatching at the interface Γfp. For the discretization of the fluid velocity and pressure, we choose finite element spaces Vf,h⊂Vf and Wf,h⊂Wf, 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 Vp,h⊂Vp and Wp,h⊂Wp 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 as follows:(35)Vh≔vh=vf,h,vp,h∈Vf,h×Vp,h,Wh=wh=wf,h,wp,h∈Wf,h×Wp,h.We employ a conforming Lagrangian finite element spaceXp,h⊂Xp to approximate the structure displacement. Note that the finite element spaces Vf,h, Vp,h, and Xp,h satisfy the prescribed homogeneous boundary conditions on the external boundaries. For the discrete Lagrange multiplier space, we take(36)Λh=Vp,h⋅np|Γfp.The semidiscrete continuous-in-time problem reads that givenpp,h0 and ηp,h0, for t∈0,T, find uf,ht∈Vf,h, pf,ht∈Wf,h, up,ht∈Vp,h, pp,ht∈Wp,h, ηp,ht∈Xp,h, and λht∈Λh such that for all vf,h∈Vf,h, wf,h∈Wf,h, vp,h∈Vp,h, wp,h∈Wp,h, ξp,h∈Xp,h, and μh∈Λh,(37)afuf,h,vf,h+apDup,h,vp,h+apeηp,h,ξp,h+aBJSuf,h,∂tηp,h;vf,h,ξp,h+bfvf,h,pf,h+bpvp,h,pp,h+αbpξp,h,pp,h+bΓvf,h,vp,h,ξp,h;λh=ff,h,vf,hΩf+fp,h,ξp,hΩp,(38)s0∂tpp,h,wp,hΩp−αbp∂tηp,h,wp,h−bpup,h,wp,h−bfuf,h,wf,h=qf,wf,hΩf+qp,wp,hΩp,(39)bΓuf,h,up,h,∂tηp,h;μh=0.We will takepp,h0 and ηp,h0 to be suitable projections of the initial data pp,0 and ηp,0.We introduce the errors for all variables as follows:(40)ef≔uf−uf,h,ep≔up−up,h,es≔ηp−ηp,h,efp≔pf−pf,h,epp≔pp−pp,heλ≔λ−λh.The following results hold cf. [23].Theorem 2 (a priori error estimation).
There exists a unique solutionuf,h,pf,h,up,h,pp,h,ηp,h,λh in L∞0,T;Vf,h×L∞0,T;Wf,h×L∞0,T;Vp,h×W1,∞0,T;Wp,h×W1,∞0,T;Xp,h×L∞0,T;Λh of the discrete problems (37)–(39), and if the solution uf,pf,up,pp,ηp,λ∈Vf×Wf×Vp×Wp×Xp×Λ of the continuous problems (26)–(28) is smooth enough, then we have(41)efL20,T;H1Ωf+epL20,T;L2Ωp+esL∞0,T;H1Ωp+epL∞0,T;L2Ωp+efL20,T;L2Ωf+eλL20,T;Λh≤Cuf,pf,up,pp,ηp,λhr,r>0.Remark 1.
ForWh=vf,h,wf,h,vp,h,wp,h,ξp,h,λh∈L∞0,T;Vf,h×L∞0,T;Wf,h×L∞0,T;Vp,h×W1,∞0,T;Wp,h×W1,∞0,T;Xp,h×L∞0,T;Λh, we can subtract (37)–(39) to (26)–(28) to obtain the Galerkin orthogonality relation for all t∈0,T:(42)afeft,vf,ht+apDest,ξp,ht+aBJSeft,∂test;vf,ht,ξp,ht+bfvf,ht,eft+bpvp,ht,ept+αbpξpt,ept+bΓvft,vpt,ξpt,eλt+s0∂tept,wp,ht−αbp∂ttest,wp,ht−bpept,wp,ht−bfeft,wf,ht=0.


## 2.4. Fully Discrete Formulation

For the time discretization, we employ the backward Euler method. Letτ be the time step, T=Nτ, and let tn=nτ, 0≤n≤N. Let δτun≔un−un−1/τ be the first-order (backward) discrete time derivation, where un≔utn. Then, the fully discrete model reads that given pp,h0=pp,h0 and ηp,h0=ηp,h0, find uf,hn∈Vf,h, pf,hn∈Wf,h, up,h∈Vp,h, pp,hn∈Wp,h, ηp,hn∈Xp,h, and λhn∈Λh, 1≤n≤N, such that for all vf,h∈Vf,h, wf,h∈Wf,h, vp,h∈Vp,h, wp,h∈Wp,h, ξp,h∈Xp,h, and μh∈Λh,(43)afuf,hn,vf,h+apDup,hn,vp,h+apeηp,hn,ξp,h+aBJSuf,hn,δτηp,hn;vf,h,ξp,h+bfvf,h,pf,hn+bpvp,h,pp,hn+αbpξp,h,pp,hn+bΓvf,h,vp,h,ξp,h;λhn=ffn,vf,hΩf+fpn,ξp,hΩp,(44)s0δτpp,h,wp,hΩp−αbpδτηp,hn,wp,h−bpup,hn,wp,h−bfuf,hn,wf,h=qfn,wf,hΩf+qpn,wp,hΩp,(45)bΓuf,hn,up,hn,δτηp,hn;μh=0.We introduce the discrete-in-time norms as follows:(46)ϕl20,T;X2≔τ∑n=0NϕnX21/2,ϕl∞0,T;X2≔max0≤n≤NϕnX,and the errors for all variables-in-time as follows:(47)efn≔uf−uf,hn,epn≔up−up,hn,esn≔ηp−ηp,hn,efpn≔pf−pf,hn,eppn≔pp−pp,hn,eλn≔λ−λhn.Remark 2.
Forn∈0,…,N and Whn=vf,hn,wf,hn,vp,hn,wp,hn,ξp,hn,μhn∈Vf,h×Wf,h×Vp,h×Wp,h×Xp,h×Λh, we can subtract (37)–(39) to (43) and (44) to obtain the Galerkin orthogonality relation:(48)afefn,vf,hn+apDesn,ξp,hn+aBJSefn,δτesn;vf,hn,ξp,hn+bfvf,hn,efn+bpvp,hn,epn+αbpξpn,epn+bΓvfn,vpn,ξpn,eλn+s0δτepn,wp,hn−αbpδτesn,wp,hn−bpepn,wp,hn−bfefn,wf,hn=0.


## 3. 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 and global indicators (Section3.1) and then the lower and upper error bounds are derived (Sections 3.4 and 3.5).


### 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.Definition 1.
(exact element residuals). Letn∈0,…,N and Whn=vf,hn,pf,hn,vp,hn,wp,hn,ξp,hn,μhn∈Vf,h×Wf,h×Vp,h×Wp,h×Xp,h×Λh be an arbitrary finite element function. The exact element residuals over a triangle or tetrahedra K∈Th and over E∈EhΓfp are defined by(49)Rf,KWhn=ffn+∇⋅σfvf,hn,wf,hn,Rp,K,1Whn=fpn+∇⋅σpξp,hn,wp,hn,Rp,K,2Whn=μK−1vp,hn+∇wp,hn,Rf,KWhn=qfn−∇⋅vf,hn,Rp,KWhn=qpn−δτs0wp,hn+α∇⋅ξp,hn+∇⋅vp,hn,RE,pf,1Whn=vf,hn⋅nf,E+δτξpn+vp,hn⋅np,E,RE,pf,2Whn=wp,h+σfvf,hn,wf,hnnf,E⋅nf,E,RE,pf,3Whn=σfvf,hn,wf,hnnf,E+σpξp,hn,wp,hnnp,E,RE,pf,4Whnj=σfvf,hn,wf,hnτf,E,j+μαBJSKj−1vf,hn−δτξp,hn⋅τf,E,j.
As it is common, these exact residuals are replaced by some finite-dimensional approximation called approximate element residualr∗,K, r∗,K, ∗∈f,p, RE,pf,l, l∈1,2,3,4. This approximation is here achieved by projecting ff and qf on the space of piecewise constant functions in Ωf and piecewise ℙ1 functions in Ωp; more precisely for all K∈Thf, we take(50)ff,Kn=pKffn,qf,Kn=pKqfn,withPK:L1K⟶ℝsuch thatpKϕ=1K∫Kϕxdx.
While for allK∈Thp, we take fp,Kn and qpn as the unique element of ℙ1Kd, respectively, ℙ1K such that(51)∫Kfp,Knx⋅qxdx=∫Kfpnxqxdx,q∈ℙ1Kd,respectively,(52)∫Kqp,Knxqxdx=∫Kqpnxqxdx,∀q∈ℙ1K.
Thereby, we define the approximate element residuals.Definition 2.
(approximate element residuals). Letn∈0…,N and Whn=vf,hn,pf,hn,vp,hn,wp,hn,ξp,hn,μhn∈Vf,h×Wf,h×Vp,h×Wp,h×Xp,h×Λh be an arbitrary finite element function. Then, the approximate element residuals are defined by(53)rf,KWhn=ff,Kn+∇⋅σfvf,hn,wf,hn,rp,K,1Whn=fp,Kn+∇⋅σpξp,hn,wp,hn,rp,K,2Whn=μK−1vp,hn+∇wp,hn,rf,KWhn=qf,Kn−∇⋅vf,hn,rp,KWhn=qp,Kn−δτs0wp,hn+α∇⋅ξp,hn+∇⋅vp,hn.
Next, introduce the gradient jump in normal direction by(54)JE,nE,fUhn≔2μDuf,hn−pf,hnI⋅nEE,ifE∈EhΩf,JE,nE,pUhn≔2μDηp,hn−pp,hnI⋅nEE,ifE∈EhΩp,where I is the identity matrix of ℝd×d.Definition 3.
(residual error estimators). Letn∈0,…,N and Uhn=uf,hn,pf,hn,up,hn,pp,hn,ηp,hn,λhn be the finite element solution of the fully discrete problems (43)–(45) in Vf,h×Wf,h×Vp,h×Wp,h×Xp,h×Λh.
We define(55)ΘKn≔ΘK,f2Uhn+ΘK,p2Uhn+ΘK,pf2Uhn1/2,where(56)ΘK,f2Uhn≔hK2rf,KUhnL2K2+rf,KUhnL2K2+∑E∈EhΩfhEJE,nE,fUhnL2E2,(57)ΘK,p2Uhn≔hK2rp,K,1UhnL2K2+rp,K,2UhnL2K2+hK2curlrp,K,2UhnL2K2+rp,KUhnL2K2+∑E∈EhΩphEJE,nE,pUhnL2E2,and(58)ΘK,pf2Uhn≔∑E∈EhΓpfhERE,pf,1UhnL2E2+∑E∈EhΓpfhERE,pf,2UhnL2E2+∑E∈EhΓpfhERE,pf,3UhnL2E2+∑E∈EhΓpfhE∑j=1d−1RE,pf,4UhnjL2E2.
Then, the residual error estimator is locally defined by(59)ΘK≔min0≤n≤NΘKn,∀K∈Th.
The global residual error estimator is given by(60)Θ≔∑K∈ThΘK21/2.
Furthermore, denote the local approximation terms by(61)ζKn≔ζK,fUhn2+ζK,pUhn21/2,where(62)ζK,fUhn2≔hK2Rf,KUhn−rf,KUhnL2K2+hK2Rf,KUhn−rf,KUhnL2K2,ifK∈Thf,ζK,pUhn2≔hK2Rp,K,1Uhn−rp,K,1UhnL2K2+hK2Rp,KUhn−rp,KUhnL2K2,ifK∈Thp.
We set(63)ζK≔min0≤n≤NζKn,ζ≔∑K∈ThζK21/2.Remark 3.
The residual character of each term on the right-hand sides of (55)–(58) is quite clear since if Uh would be the exact solution of (2)–(10), then they would vanish.


### 3.2. Analytical Tools

#### 3.2.1. 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∈Th and E∈EK, we let bK and bE be the usual simplex-bubble and face-bubble functions (see (1.5) and (1.6) in [29]), respectively. In particular, bK satisfies bK∈ℙ3K, suppbK⊆K, bK=0,on∂K, and 0≤bK≤1,onK. Similarly, bE∈ℙ2K, suppbE⊆ωE≔K′∈Th:E∈EK′, bE=0,on∂K∖E and 0≤bE≤1,inωE. We also recall from [36] that, given k∈ℕ, there exists an extension operator L:CE⟶CK that satisfies Lp∈ℙkK and Lp|E=p,∀p∈ℙkE. A corresponding vectorial version of L, that is, the component-wise application of L, is denoted by L. Additional properties of bK, bE, and L are collected in the following lemma (see [36]).Lemma 1.
Givenk∈ℕ∗, there exist positive constants depending only on k and shape-regularity of the triangulations (minimum angle condition), such that for each simplex K and E∈EK, there hold(64)ϕK≲ϕbK1/2K≲ϕK,∀ϕ∈ℙkK,ϕbK1,K≲hK−1ϕK,∀ϕ∈ℙkK,ψE≲bE1/2ψE≲ψE,∀ψ∈ℙkE,LψK+hELψ1,K≲hE1/2ψE,∀ψ∈ℙkE.Lemma 2 (continuous trace inequality).
There exists a positive constantβ1>0 depending only on σ0 such that(65)v∂K2≤β1vKv1,K,∀K∈Th,∀v∈H1Kd.


#### 3.2.2. Clément Interpolation Operator

In order to derive the upper error bounds, we introduce the Clément interpolation operatorICl0:H01Ω⟶PcbTh that approximates optimally nonsmooth functions by continuous piecewise linear functions:(66)PcbTh≔v∈C0Ω¯:v|K∈ℙ1K,∀K∈Thandv=0,on∂Ω.In addition, we will make use of a vector valued version ofICl0, that is, ICl0:H01Ωd⟶PcbThd, which is defined component-wise by ICl0. The following lemma establishes the local approximation properties of ICl0 (and hence of ICl0), for a proof (see Section 3 in [37]).Lemma 3.
There exist constantsC1,C2>0, independent of h, such that for all v∈H01Ω, there hold(67)v−ICl0vK≤C1hKv1,ΔK,∀K∈Th,0.2cm,v−ICl0vE≤C2hE1/2v1,ΔE,∀E∈Eh,where ΔK≔∪K′∈Th:K′∩K≠∅ and ΔE≔∪K′∈Th:K′∩E≠∅.


#### 3.2.3. Helmholtz Decomposition

Lemma 4 (see [7]).
There existsCp>0 such that every vp∈Hdiv;Ωp can be decomposed as vp=w+curlβ, where w∈H1Ωpd, β∈H1Ωp, ∫Ωpβxdx=0, and(68)w1,Ωp+β1,Ωp≤CpvpHdiv;Ωp.


### 3.3. Error Equation

We setH=Vf×Wf×Vp×Wp×Xp×Λ and Hh=Vf,h×Wf,h×Vp,h×Wp,h×Xp,h×Λh.Forn∈0,…,N, let Un=ufn,pfn,upn,ppn,ηpn,λn∈H and W=vfn,wfn,vpn,wpn,ξpn,μn∈H. We define now the operator A by(69)AUn,Wn≔afufn,vfn+apDupn,vpn+bpvpn,ppn+αbpξpn,ppn+bΓvfn,vpn,ξpn,λn+s0δτppn,wpnΩp−αbpδτηpn,wpn−bpupn,wpn−bfufn,wfn+bΓufn,upn,δτηpn;μn,FWn≔ffn,vfnΩf+fpn,ξpn+qfn,wfnΩf+qpn,wpnΩp.Then, the continuous problems (26)–(28) are equivalent to the following: for n∈0,…,N, find Un∈H such that we have(70)AUn,Wn=FWn,∀Wn∈H.We define the discrete version by the same way: findUhn∈Hh such that(71)AUhn,Whn=FWhn,∀Whn∈Hh.Since for alln∈0,…,N and Whn∈Hh, AUn−Uhn,Whn=0, then from (70), we obtain(72)AUn−Uhn,Wn=AUn−Uhn,Wn−Whn=AUn,Wn−Whn−AUhn,Wn−Whn=FWn−Whn−AUhn,Wn−Whn=ffn,vfn−vf,hnΩf+fpn,ξpn−ξp,hnΩp+qfn,wfn−wf,hnΩf+qpn,wpn−wp,hnΩp−AUhn,Wn−Whn=ffn−ff,hn,vfn−vf,hnΩf+fpn−fp,hn,ξpn−ξp,hnΩp+qfn−qf,hn,wfn−wf,hnΩf+qpn−qp,hn,wpn−wp,hnΩp+ff,hn,vfn−vf,hnΩf+fp,hn,ξpn−ξp,hnΩp+qf,hn,wfn−wf,hnΩf+qp,hn,wpn−wp,hnΩp−AUhn,Wn−Whn.Now, we set(73)AUn−Uhn,Wn=∑K∈ThAKUn−Uhn,Wn.Then, from (72), we can write(74)AKUn−Uhn,Wn≔ffn−ff,hn,vfn−vf,hnΩf∩K+fpn−fp,hn,ξpn−ξp,hnΩp∩K+qfn−qf,hn,wfn−wf,hnΩf∩K+qpn−qp,hn,wpn−wp,hnΩp∩K+ff,hn,vfn−vf,hnΩf∩K+fp,hn,ξpn−ξp,hnΩp∩K+qf,hn,wfn−wf,hnΩf∩K+qp,hn,wpn−wp,hnΩp∩K−AKUhn,Wn−Whn=Rf,KUhn−rf,KUhn,vfn−vf,hnΩf∩K+Rf,KUhn−rf,KUhn,wfn−wf,hnΩf∩K+Rp,K,1Uhn−rp,K,1Uhn,ξpn−ξp,hnΩp∩K+Rp,KUhn−rp,KUhn,wpn−wp,hnΩp∩K+BKUhn,Wn−Whn,where(75)BKUhn,Wn−Whn=ff,Kn,vfn−vf,hnΩf∩K+fp,Kn,ξpn−ξp,hnΩp∩K+qf,Kn,wfn−wf,hnΩf∩K+qp,Kn,wpn−wp,hnΩp∩K−AKUhn,Wn−Whn.


### 3.4. Reliability of the A Posteriori Error Estimator

The first main result is given by the following theorem.Theorem 3 (upper error bound).
LetU=Unn=0N such that Un=ufn,pfn,upn,ppn,ηpn,λn∈H be the exact solution and Uh=Uhnn=0N with Uhn=uf,hn,pf,hn,up,hn,pp,hn,ηp,hn,λhn∈Hh be the finite element solution. Then, there exist a positive constant Cupper such that the error is bounded globally from above by(76)U−Uhl∞0,T;Hh≤CupperΘ+ζ.Proof.
Letn∈0,…,N. We consider the residual equation (73), and we set Wn=vfn,wfn,vpn,wpn,ξpn,μn. We take Whn=vf,hn,0,vp,hn,0,ξp,hn,0 with vf,hn=ICl0vfn and ξp,hn=ICl0ξpn. As vpn∈Hdiv;Ωp, then by Lemma 4, vpn admits the decomposition vpn=wpn+curlβn where wpn∈H1Ωpd and βn∈H1Ωp with ∫Ωpβnxdx=0 and wpnH1Ωpd+βnH1Ωp≤CpnvpnVp. We consider vp,hn=wp,hn+curlβp,hn with wp,hn=IClwpn and βp,hn=IClβpn. Thus, vpn−vp,hn=wpn−wp,hn+curlβpn−βp,hn. Now, we recall(77)BKUhn,Wn−Whn=ff,Kn,vfn−vf,hnΩf∩K+fp,Kn,ξpn−ξp,hnΩp∩K+qf,Kn,wfn−wf,hnΩf∩K+qp,Kn,wpn−wp,hnΩp∩K−AKUhn,Wn−Whn.
Therefore, integrating by parts element by element, we may write(78)BKUhn,Wn−Whn=∑i=112ϒi,K,where(79)ϒ1,K≔rf,KUhn,vfn−vf,hnΩf∩K,ϒ2,K≔rp,K,1Uhn,ξpn−ξp,hnΩp∩K,ϒ3,K≔−rp,K,2Uhn,wpn−wp,hnΩp∩K,ϒ4,K≔−curlrp,K,2Uhn,βpn−βp,hnΩp∩K,ϒ5,K≔rf,KUhn,wfnΩf∩K,ϒ6,K≔rp,KUhn,wpnΩp∩K,ϒ7,K≔−∑E∈EhΩf∩KJE,nE,fUhn,vfn−vf,hnE,ϒ8,K≔−∑E∈EhΩp∩KJE,nE,pUhn,ξpn−ξp,hnE,ϒ9,K≔∑E∈EhΓfpRE,pf,1Uhn,vfn−vf,hn⋅nf,EE,ϒ10,K≔∑E∈EhΓfpRE,pf,2Uhn,vfn−vf,hn⋅nf,EE,ϒ11,K≔∑E∈EhΓfpRE,pf,3Uhn,ξpn−ξp,hnE,ϒ12,K≔−∑E∈EhΓfp∑j=1d−1RE,pf,4Uhnj,vfn−vf,hn⋅τjE.
Thus, from equations (74) and (78), we obtain(80)AKUn−Uhn,Wn≔Rf,KUhn−rf,KUhn,vfn−vf,hnΩf∩K+Rf,KUhn−rf,KUhn,wfn−wf,hnΩf∩K+Rp,K,1Uhn−rp,K,1Uhn,ξpn−ξp,hnΩp∩K+Rp,KUhn−rp,KUhn,wpn−wp,hnΩp∩K+∑i=112ηi,K.
Coercivity of operatorA leads to inf-sup condition as follows:(81)Un−UhnHh≤supWn∈HAUn−Uhn,WnWnH.
By consequently, identities (73) and (80), inf-sup condition of operator A (81), Cauchy–Schwarz inequality, estimation of Lemma 4, and the approximation properties of Lemma 3 imply the required estimate and finish the proof. □


### 3.5. Efficiency of the A Posteriori Error Estimator

ForK∈Th, we set(82)HhK=Vf,hK×Wf,hK×Vp,hK×Wp,hK×Xp,hK×ΛhK,where(83)Vf,hK≔vf,h|K:vf,h∈Vf,h,Wf,hK≔wf,h|K:wf,h∈Wf,h;Vp,hK≔vp,h|K:vp,h∈Vp,h,Wp,hK≔wp,h|K:wp,h∈Wp,h;Xp,hK≔ξp,h|K:ξp,h∈Xp,h,ΛhK≔λh|K:λh∈Λh.The family estimatorΘKK∈Th is consider efficient if it satisfies the following theorem.Theorem 4 (lower error bound).
LetU=Unn=0N such that Un=ufn,pfn,upn,ppn,ηpn,λn∈H be the exact solution and Uh=Uhnn=0N with Uhn=uf,hn,pf,hn,up,hn,pp,hn,ηp,hn,λhn∈Hh be the finite element solution. Then, there exist a positive constant Clow such that the error is bounded locally from below for all K∈Th by(84)ΘK≤ClowU−Uhl∞0,T;HhK+∑K′⊂w˜KζK′,where w˜K is a finite union of neighbouring elements of K and ⋅HhK is the product norm.Proof.
The lower bound is proved using the standard elementwise integration by parts, namely, error equation of Section3.3 (i.e., identity (73) and equation (80)) and some inverse estimates of Lemma 1 (cf. [38] for details).


## 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.Definition 1.
(exact element residuals). Letn∈0,…,N and Whn=vf,hn,pf,hn,vp,hn,wp,hn,ξp,hn,μhn∈Vf,h×Wf,h×Vp,h×Wp,h×Xp,h×Λh be an arbitrary finite element function. The exact element residuals over a triangle or tetrahedra K∈Th and over E∈EhΓfp are defined by(49)Rf,KWhn=ffn+∇⋅σfvf,hn,wf,hn,Rp,K,1Whn=fpn+∇⋅σpξp,hn,wp,hn,Rp,K,2Whn=μK−1vp,hn+∇wp,hn,Rf,KWhn=qfn−∇⋅vf,hn,Rp,KWhn=qpn−δτs0wp,hn+α∇⋅ξp,hn+∇⋅vp,hn,RE,pf,1Whn=vf,hn⋅nf,E+δτξpn+vp,hn⋅np,E,RE,pf,2Whn=wp,h+σfvf,hn,wf,hnnf,E⋅nf,E,RE,pf,3Whn=σfvf,hn,wf,hnnf,E+σpξp,hn,wp,hnnp,E,RE,pf,4Whnj=σfvf,hn,wf,hnτf,E,j+μαBJSKj−1vf,hn−δτξp,hn⋅τf,E,j.
As it is common, these exact residuals are replaced by some finite-dimensional approximation called approximate element residualr∗,K, r∗,K, ∗∈f,p, RE,pf,l, l∈1,2,3,4. This approximation is here achieved by projecting ff and qf on the space of piecewise constant functions in Ωf and piecewise ℙ1 functions in Ωp; more precisely for all K∈Thf, we take(50)ff,Kn=pKffn,qf,Kn=pKqfn,withPK:L1K⟶ℝsuch thatpKϕ=1K∫Kϕxdx.
While for allK∈Thp, we take fp,Kn and qpn as the unique element of ℙ1Kd, respectively, ℙ1K such that(51)∫Kfp,Knx⋅qxdx=∫Kfpnxqxdx,q∈ℙ1Kd,respectively,(52)∫Kqp,Knxqxdx=∫Kqpnxqxdx,∀q∈ℙ1K.
Thereby, we define the approximate element residuals.Definition 2.
(approximate element residuals). Letn∈0…,N and Whn=vf,hn,pf,hn,vp,hn,wp,hn,ξp,hn,μhn∈Vf,h×Wf,h×Vp,h×Wp,h×Xp,h×Λh be an arbitrary finite element function. Then, the approximate element residuals are defined by(53)rf,KWhn=ff,Kn+∇⋅σfvf,hn,wf,hn,rp,K,1Whn=fp,Kn+∇⋅σpξp,hn,wp,hn,rp,K,2Whn=μK−1vp,hn+∇wp,hn,rf,KWhn=qf,Kn−∇⋅vf,hn,rp,KWhn=qp,Kn−δτs0wp,hn+α∇⋅ξp,hn+∇⋅vp,hn.
Next, introduce the gradient jump in normal direction by(54)JE,nE,fUhn≔2μDuf,hn−pf,hnI⋅nEE,ifE∈EhΩf,JE,nE,pUhn≔2μDηp,hn−pp,hnI⋅nEE,ifE∈EhΩp,where I is the identity matrix of ℝd×d.Definition 3.
(residual error estimators). Letn∈0,…,N and Uhn=uf,hn,pf,hn,up,hn,pp,hn,ηp,hn,λhn be the finite element solution of the fully discrete problems (43)–(45) in Vf,h×Wf,h×Vp,h×Wp,h×Xp,h×Λh.
We define(55)ΘKn≔ΘK,f2Uhn+ΘK,p2Uhn+ΘK,pf2Uhn1/2,where(56)ΘK,f2Uhn≔hK2rf,KUhnL2K2+rf,KUhnL2K2+∑E∈EhΩfhEJE,nE,fUhnL2E2,(57)ΘK,p2Uhn≔hK2rp,K,1UhnL2K2+rp,K,2UhnL2K2+hK2curlrp,K,2UhnL2K2+rp,KUhnL2K2+∑E∈EhΩphEJE,nE,pUhnL2E2,and(58)ΘK,pf2Uhn≔∑E∈EhΓpfhERE,pf,1UhnL2E2+∑E∈EhΓpfhERE,pf,2UhnL2E2+∑E∈EhΓpfhERE,pf,3UhnL2E2+∑E∈EhΓpfhE∑j=1d−1RE,pf,4UhnjL2E2.
Then, the residual error estimator is locally defined by(59)ΘK≔min0≤n≤NΘKn,∀K∈Th.
The global residual error estimator is given by(60)Θ≔∑K∈ThΘK21/2.
Furthermore, denote the local approximation terms by(61)ζKn≔ζK,fUhn2+ζK,pUhn21/2,where(62)ζK,fUhn2≔hK2Rf,KUhn−rf,KUhnL2K2+hK2Rf,KUhn−rf,KUhnL2K2,ifK∈Thf,ζK,pUhn2≔hK2Rp,K,1Uhn−rp,K,1UhnL2K2+hK2Rp,KUhn−rp,KUhnL2K2,ifK∈Thp.
We set(63)ζK≔min0≤n≤NζKn,ζ≔∑K∈ThζK21/2.Remark 3.
The residual character of each term on the right-hand sides of (55)–(58) is quite clear since if Uh would be the exact solution of (2)–(10), then they would vanish.


## 3.2. Analytical Tools

### 3.2.1. 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∈Th and E∈EK, we let bK and bE be the usual simplex-bubble and face-bubble functions (see (1.5) and (1.6) in [29]), respectively. In particular, bK satisfies bK∈ℙ3K, suppbK⊆K, bK=0,on∂K, and 0≤bK≤1,onK. Similarly, bE∈ℙ2K, suppbE⊆ωE≔K′∈Th:E∈EK′, bE=0,on∂K∖E and 0≤bE≤1,inωE. We also recall from [36] that, given k∈ℕ, there exists an extension operator L:CE⟶CK that satisfies Lp∈ℙkK and Lp|E=p,∀p∈ℙkE. A corresponding vectorial version of L, that is, the component-wise application of L, is denoted by L. Additional properties of bK, bE, and L are collected in the following lemma (see [36]).Lemma 1.
Givenk∈ℕ∗, there exist positive constants depending only on k and shape-regularity of the triangulations (minimum angle condition), such that for each simplex K and E∈EK, there hold(64)ϕK≲ϕbK1/2K≲ϕK,∀ϕ∈ℙkK,ϕbK1,K≲hK−1ϕK,∀ϕ∈ℙkK,ψE≲bE1/2ψE≲ψE,∀ψ∈ℙkE,LψK+hELψ1,K≲hE1/2ψE,∀ψ∈ℙkE.Lemma 2 (continuous trace inequality).
There exists a positive constantβ1>0 depending only on σ0 such that(65)v∂K2≤β1vKv1,K,∀K∈Th,∀v∈H1Kd.


### 3.2.2. Clément Interpolation Operator

In order to derive the upper error bounds, we introduce the Clément interpolation operatorICl0:H01Ω⟶PcbTh that approximates optimally nonsmooth functions by continuous piecewise linear functions:(66)PcbTh≔v∈C0Ω¯:v|K∈ℙ1K,∀K∈Thandv=0,on∂Ω.In addition, we will make use of a vector valued version ofICl0, that is, ICl0:H01Ωd⟶PcbThd, which is defined component-wise by ICl0. The following lemma establishes the local approximation properties of ICl0 (and hence of ICl0), for a proof (see Section 3 in [37]).Lemma 3.
There exist constantsC1,C2>0, independent of h, such that for all v∈H01Ω, there hold(67)v−ICl0vK≤C1hKv1,ΔK,∀K∈Th,0.2cm,v−ICl0vE≤C2hE1/2v1,ΔE,∀E∈Eh,where ΔK≔∪K′∈Th:K′∩K≠∅ and ΔE≔∪K′∈Th:K′∩E≠∅.


### 3.2.3. Helmholtz Decomposition

Lemma 4 (see [7]).
There existsCp>0 such that every vp∈Hdiv;Ωp can be decomposed as vp=w+curlβ, where w∈H1Ωpd, β∈H1Ωp, ∫Ωpβxdx=0, and(68)w1,Ωp+β1,Ωp≤CpvpHdiv;Ωp.


## 3.2.1. 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∈Th and E∈EK, we let bK and bE be the usual simplex-bubble and face-bubble functions (see (1.5) and (1.6) in [29]), respectively. In particular, bK satisfies bK∈ℙ3K, suppbK⊆K, bK=0,on∂K, and 0≤bK≤1,onK. Similarly, bE∈ℙ2K, suppbE⊆ωE≔K′∈Th:E∈EK′, bE=0,on∂K∖E and 0≤bE≤1,inωE. We also recall from [36] that, given k∈ℕ, there exists an extension operator L:CE⟶CK that satisfies Lp∈ℙkK and Lp|E=p,∀p∈ℙkE. A corresponding vectorial version of L, that is, the component-wise application of L, is denoted by L. Additional properties of bK, bE, and L are collected in the following lemma (see [36]).Lemma 1.
Givenk∈ℕ∗, there exist positive constants depending only on k and shape-regularity of the triangulations (minimum angle condition), such that for each simplex K and E∈EK, there hold(64)ϕK≲ϕbK1/2K≲ϕK,∀ϕ∈ℙkK,ϕbK1,K≲hK−1ϕK,∀ϕ∈ℙkK,ψE≲bE1/2ψE≲ψE,∀ψ∈ℙkE,LψK+hELψ1,K≲hE1/2ψE,∀ψ∈ℙkE.Lemma 2 (continuous trace inequality).
There exists a positive constantβ1>0 depending only on σ0 such that(65)v∂K2≤β1vKv1,K,∀K∈Th,∀v∈H1Kd.


## 3.2.2. Clément Interpolation Operator

In order to derive the upper error bounds, we introduce the Clément interpolation operatorICl0:H01Ω⟶PcbTh that approximates optimally nonsmooth functions by continuous piecewise linear functions:(66)PcbTh≔v∈C0Ω¯:v|K∈ℙ1K,∀K∈Thandv=0,on∂Ω.In addition, we will make use of a vector valued version ofICl0, that is, ICl0:H01Ωd⟶PcbThd, which is defined component-wise by ICl0. The following lemma establishes the local approximation properties of ICl0 (and hence of ICl0), for a proof (see Section 3 in [37]).Lemma 3.
There exist constantsC1,C2>0, independent of h, such that for all v∈H01Ω, there hold(67)v−ICl0vK≤C1hKv1,ΔK,∀K∈Th,0.2cm,v−ICl0vE≤C2hE1/2v1,ΔE,∀E∈Eh,where ΔK≔∪K′∈Th:K′∩K≠∅ and ΔE≔∪K′∈Th:K′∩E≠∅.


## 3.2.3. Helmholtz Decomposition

Lemma 4 (see [7]).
There existsCp>0 such that every vp∈Hdiv;Ωp can be decomposed as vp=w+curlβ, where w∈H1Ωpd, β∈H1Ωp, ∫Ωpβxdx=0, and(68)w1,Ωp+β1,Ωp≤CpvpHdiv;Ωp.


## 3.3. Error Equation

We setH=Vf×Wf×Vp×Wp×Xp×Λ and Hh=Vf,h×Wf,h×Vp,h×Wp,h×Xp,h×Λh.Forn∈0,…,N, let Un=ufn,pfn,upn,ppn,ηpn,λn∈H and W=vfn,wfn,vpn,wpn,ξpn,μn∈H. We define now the operator A by(69)AUn,Wn≔afufn,vfn+apDupn,vpn+bpvpn,ppn+αbpξpn,ppn+bΓvfn,vpn,ξpn,λn+s0δτppn,wpnΩp−αbpδτηpn,wpn−bpupn,wpn−bfufn,wfn+bΓufn,upn,δτηpn;μn,FWn≔ffn,vfnΩf+fpn,ξpn+qfn,wfnΩf+qpn,wpnΩp.Then, the continuous problems (26)–(28) are equivalent to the following: for n∈0,…,N, find Un∈H such that we have(70)AUn,Wn=FWn,∀Wn∈H.We define the discrete version by the same way: findUhn∈Hh such that(71)AUhn,Whn=FWhn,∀Whn∈Hh.Since for alln∈0,…,N and Whn∈Hh, AUn−Uhn,Whn=0, then from (70), we obtain(72)AUn−Uhn,Wn=AUn−Uhn,Wn−Whn=AUn,Wn−Whn−AUhn,Wn−Whn=FWn−Whn−AUhn,Wn−Whn=ffn,vfn−vf,hnΩf+fpn,ξpn−ξp,hnΩp+qfn,wfn−wf,hnΩf+qpn,wpn−wp,hnΩp−AUhn,Wn−Whn=ffn−ff,hn,vfn−vf,hnΩf+fpn−fp,hn,ξpn−ξp,hnΩp+qfn−qf,hn,wfn−wf,hnΩf+qpn−qp,hn,wpn−wp,hnΩp+ff,hn,vfn−vf,hnΩf+fp,hn,ξpn−ξp,hnΩp+qf,hn,wfn−wf,hnΩf+qp,hn,wpn−wp,hnΩp−AUhn,Wn−Whn.Now, we set(73)AUn−Uhn,Wn=∑K∈ThAKUn−Uhn,Wn.Then, from (72), we can write(74)AKUn−Uhn,Wn≔ffn−ff,hn,vfn−vf,hnΩf∩K+fpn−fp,hn,ξpn−ξp,hnΩp∩K+qfn−qf,hn,wfn−wf,hnΩf∩K+qpn−qp,hn,wpn−wp,hnΩp∩K+ff,hn,vfn−vf,hnΩf∩K+fp,hn,ξpn−ξp,hnΩp∩K+qf,hn,wfn−wf,hnΩf∩K+qp,hn,wpn−wp,hnΩp∩K−AKUhn,Wn−Whn=Rf,KUhn−rf,KUhn,vfn−vf,hnΩf∩K+Rf,KUhn−rf,KUhn,wfn−wf,hnΩf∩K+Rp,K,1Uhn−rp,K,1Uhn,ξpn−ξp,hnΩp∩K+Rp,KUhn−rp,KUhn,wpn−wp,hnΩp∩K+BKUhn,Wn−Whn,where(75)BKUhn,Wn−Whn=ff,Kn,vfn−vf,hnΩf∩K+fp,Kn,ξpn−ξp,hnΩp∩K+qf,Kn,wfn−wf,hnΩf∩K+qp,Kn,wpn−wp,hnΩp∩K−AKUhn,Wn−Whn.


## 3.4. Reliability of the A Posteriori Error Estimator

The first main result is given by the following theorem.Theorem 3 (upper error bound).
LetU=Unn=0N such that Un=ufn,pfn,upn,ppn,ηpn,λn∈H be the exact solution and Uh=Uhnn=0N with Uhn=uf,hn,pf,hn,up,hn,pp,hn,ηp,hn,λhn∈Hh be the finite element solution. Then, there exist a positive constant Cupper such that the error is bounded globally from above by(76)U−Uhl∞0,T;Hh≤CupperΘ+ζ.Proof.
Letn∈0,…,N. We consider the residual equation (73), and we set Wn=vfn,wfn,vpn,wpn,ξpn,μn. We take Whn=vf,hn,0,vp,hn,0,ξp,hn,0 with vf,hn=ICl0vfn and ξp,hn=ICl0ξpn. As vpn∈Hdiv;Ωp, then by Lemma 4, vpn admits the decomposition vpn=wpn+curlβn where wpn∈H1Ωpd and βn∈H1Ωp with ∫Ωpβnxdx=0 and wpnH1Ωpd+βnH1Ωp≤CpnvpnVp. We consider vp,hn=wp,hn+curlβp,hn with wp,hn=IClwpn and βp,hn=IClβpn. Thus, vpn−vp,hn=wpn−wp,hn+curlβpn−βp,hn. Now, we recall(77)BKUhn,Wn−Whn=ff,Kn,vfn−vf,hnΩf∩K+fp,Kn,ξpn−ξp,hnΩp∩K+qf,Kn,wfn−wf,hnΩf∩K+qp,Kn,wpn−wp,hnΩp∩K−AKUhn,Wn−Whn.
Therefore, integrating by parts element by element, we may write(78)BKUhn,Wn−Whn=∑i=112ϒi,K,where(79)ϒ1,K≔rf,KUhn,vfn−vf,hnΩf∩K,ϒ2,K≔rp,K,1Uhn,ξpn−ξp,hnΩp∩K,ϒ3,K≔−rp,K,2Uhn,wpn−wp,hnΩp∩K,ϒ4,K≔−curlrp,K,2Uhn,βpn−βp,hnΩp∩K,ϒ5,K≔rf,KUhn,wfnΩf∩K,ϒ6,K≔rp,KUhn,wpnΩp∩K,ϒ7,K≔−∑E∈EhΩf∩KJE,nE,fUhn,vfn−vf,hnE,ϒ8,K≔−∑E∈EhΩp∩KJE,nE,pUhn,ξpn−ξp,hnE,ϒ9,K≔∑E∈EhΓfpRE,pf,1Uhn,vfn−vf,hn⋅nf,EE,ϒ10,K≔∑E∈EhΓfpRE,pf,2Uhn,vfn−vf,hn⋅nf,EE,ϒ11,K≔∑E∈EhΓfpRE,pf,3Uhn,ξpn−ξp,hnE,ϒ12,K≔−∑E∈EhΓfp∑j=1d−1RE,pf,4Uhnj,vfn−vf,hn⋅τjE.
Thus, from equations (74) and (78), we obtain(80)AKUn−Uhn,Wn≔Rf,KUhn−rf,KUhn,vfn−vf,hnΩf∩K+Rf,KUhn−rf,KUhn,wfn−wf,hnΩf∩K+Rp,K,1Uhn−rp,K,1Uhn,ξpn−ξp,hnΩp∩K+Rp,KUhn−rp,KUhn,wpn−wp,hnΩp∩K+∑i=112ηi,K.
Coercivity of operatorA leads to inf-sup condition as follows:(81)Un−UhnHh≤supWn∈HAUn−Uhn,WnWnH.
By consequently, identities (73) and (80), inf-sup condition of operator A (81), Cauchy–Schwarz inequality, estimation of Lemma 4, and the approximation properties of Lemma 3 imply the required estimate and finish the proof. □


## 3.5. Efficiency of the A Posteriori Error Estimator

ForK∈Th, we set(82)HhK=Vf,hK×Wf,hK×Vp,hK×Wp,hK×Xp,hK×ΛhK,where(83)Vf,hK≔vf,h|K:vf,h∈Vf,h,Wf,hK≔wf,h|K:wf,h∈Wf,h;Vp,hK≔vp,h|K:vp,h∈Vp,h,Wp,hK≔wp,h|K:wp,h∈Wp,h;Xp,hK≔ξp,h|K:ξp,h∈Xp,h,ΛhK≔λh|K:λh∈Λh.The family estimatorΘKK∈Th is consider efficient if it satisfies the following theorem.Theorem 4 (lower error bound).
LetU=Unn=0N such that Un=ufn,pfn,upn,ppn,ηpn,λn∈H be the exact solution and Uh=Uhnn=0N with Uhn=uf,hn,pf,hn,up,hn,pp,hn,ηp,hn,λhn∈Hh be the finite element solution. Then, there exist a positive constant Clow such that the error is bounded locally from below for all K∈Th by(84)ΘK≤ClowU−Uhl∞0,T;HhK+∑K′⊂w˜KζK′,where w˜K is a finite union of neighbouring elements of K and ⋅HhK is the product norm.Proof.
The lower bound is proved using the standard elementwise integration by parts, namely, error equation of Section3.3 (i.e., identity (73) and equation (80)) and some inverse estimates of Lemma 1 (cf. [38] for details).


## 4. Discussion

In this paper, we have discussed a posteriori error estimates for a finite element approximation of the Stokes–Biot system where homogeneous boundary conditions are employed. The approach utilizes a fully discrete conforming finite element method. A residual type a posteriori error estimator is provided, that is, both reliable and efficient.In a future paper, we will study the influence of nonhomogeneous boundary conditions on the a posteriori error indicators presented in this work. Further, it is well known that an internal layer appears at the interfaceΓfp as the permeability tensor degenerates; in that case, anisotropic meshes have to be used in this layer (see, for instance, [10]). Hence, we intend to extend our results to such anisotropic meshes.


## 5. Nomenclatures

(i)
Ω⊂ℝd,d∈2,3: bounded domain(ii)
Ωp: the poroelastic medium domain(iii)
Ωf=Ω∖Ω¯d(iv)
Γfp=∂Ωf∩∂Ωp(v)
Γ∗=∂Ω∗∖Γfp, ∗=f,p(vi)
nf(respectively, np): the unit outward normal vector along ∂Ωf (respectively, ∂Ωp)(vii)
uf: the fluid velocity in Ωf(viii)
pf: the fluid pressure in Ωf(ix)
up,ηp: the fluid velocities in Ωp(x)
pp: the fluid pressure in Ωp(xi)
In2D, the curl of a scalar function w is given as usual by curlw≔∂w/∂x2,−∂w/∂x1⊤(xii)
In3D, the curl of a vector function w=w1,w2,w3 is given as usual by curlw≔∇×w, namely, curlw≔∂w3/∂x2−∂w2/∂x3,∂w1/∂x3−∂w3/∂x1,∂w2/∂x1−∂w1/∂x2(xiii)
ℙk: the space of polynomials of total degree not larger than k(xiv)
Th: triangulation of Ω(xv)
Th∗: the corresponding induced triangulation of Ω∗, ∗∈f,p(xvi)
For anyK∈Th, hK is the diameter of K and ρK=2rK is the diameter of the largest ball inscribed into K(xvii)
h≔maxK∈ThhKand σh≔maxK∈ThhK/ρK(xviii)
Eh: the set of all the edges or faces of the triangulation(xix)
EK: the set of all the edges (N=2) or faces (N=3) of a element K(xx)
Eh≔∪K∈ThEK(xxi)
NK: the set of all the vertices of a element K(xxii)
Nh≔∪K∈ThNK(xxiii)
ForA⊂Ω¯, EhA≔E∈Eh:E⊂A(xxiv)
ForE∈EhΩ∗, we associate a unit vector nE,∗ such that nE,∗ is orthogonal to E and equals to the unit exterior normal vector to ∂Ω∗, ∗∈f,p(xxv)
ForE∈EhΩ∗, ϕE is the jump across E in the direction of nE,∗(xxvi)
In order to avoid excessive use of constants, the abbreviationsx≲yand x∼y stand for x≤cy and c1x≤y≤c2x, respectively, with positive constants independent of x, y or Th(xxvii)
∂t=∂/∂t(xxviii)
ϕL20,T;X≔∫0TϕsX2ds1/2(xxix)
ϕL∞0,T;X≔esssup0≤t≤TϕtX(xxx)
ϕW1,∞0,T;X≔esssup0≤t≤TϕtX,∂tϕtX(xxxi)
ϕl∞0,T;X2≔max0≤n≤NϕnX(xxxii)
ϕl20,T;X≔τ∑n=0NϕnX21/2(xxxiii)
⋅Hh: product norm on Hh(xxxiv)
⋅HhK: local product norm on HhK.



---

*Source: 5982948-2021-11-25.xml*