# A Posteriori Error Estimates for a Nonconforming Finite Element Discretization of the Stokes–Biot System

Authors: Koffi Wilfrid Houédanou
Journal: Discrete Dynamics in Nature and Society
 (2022)

---

## Abstract

This paper presents an a posteriori error estimator for a (piecewise linear) nonconforming finite element approximation of the problem defining the interaction between a free fluid and poroelastic structure. The free fluid is governed by the Stokes equations, while the flow in the poroelastic medium is modeled using the Biot poroelasticity system. Equilibrium and kinematic conditions are imposed on the interface. The approach utilizes the same nonconforming Crouzeix–Raviart element discretization on the entire domain. For this discretization, we derive a residual indicator based on the jumps of the normal derivative of the nonconforming approximation. Lower and upper bounds form the main results with minimal assumptions on the mesh.


---

## Body

## 1. 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 [1] 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 [2–10] 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–19]. A numerical study of the problem, using Navier–Stokes equations for the fluid, is presented in [11, 20], 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.Nonphysical pressure oscillations are observed in finite element calculations of Biot’s poroelastic equations in low-permeable media. These pressure oscillations may be understood as a failure of compatibility between the finite element spaces, rather than elastic locking. In [21], Joachim Berdal Haga et al. have presented evidence to support this view by comparing and contrasting the pressure oscillations in low-permeable porous media with those in low-compressible porous media. As a consequence, it is possible to use established families of stable mixed elements as candidates for choosing finite element spaces for Biot’s equations. Through comparison with the displacement-solid pressure mixed formulation of linear elasticity, they identify the spurious pressure modes as a specific consequence of a vanishing Brezzi inf-sup constant. Since the Brezzi inf-sup condition for the poroelastic equations takes on a similar form as in, e.g., the mixed linear elasticity or Stokes problem, this identification opens up the field to a plethora of stable element candidates. These can be used directly for the basic solid displacement-fluid pressure two-field formulation of poroelasticity, or in combinations for the various three- and four-field formulations involving solid pressure and/or fluid velocity [21].Finite element analysis of an arbitrary Lagrangian–Eulerian method for Stokes/parabolic moving interface problem with jump coefficients has been studied in [22]. The authors in [23] study a numerical solution of the coupled system of the time-dependent Stokes and fully dynamic Biot equations. They establish the stability of the scheme and derive 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 [24], Jing Wen and Yinnian He consider a strongly conservative discretization for the rearranged Stokes–Biot model based on interior penalty discontinuous Galerkin method and mixed finite element method. The existence and uniqueness of a solution of the numerical scheme have been presented. Then, the analysis of stability and priori error estimates have been derived. The numerical examples under uniform meshes, which will 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 [25], while Ambartsumyan et al. study in [26] flow and transport in fractured poroelastic media using Stokes flow in the fractures and the Biot model in the porous media. In [27], 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 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 equi-distribute the computational effort and optimize the approximation efficiency. Since the pioneering work of Babuška and Rheinboldt [28–31], adaptive finite element methods based on a posteriori error estimates have been extensively investigated.In the article [32], the author studies a stabilized nonconforming mixed finite element method using the Crouzeix–Raviart element for the Stokes–Biot problem. Considering the mixed formulation of the Darcy problem, the fluid velocity and pressure are treated as functions defined in the entire domain. Existence, uniqueness of the finite element solution of the corresponding discrete problem, and a priori estimates have been shown. The proofs use the standard theory for mixed problems. The approach presented is independent of the normal vectors of the interior edges in both regions, thus making the resulting finite element matrix sparser.We use a nonconforming finite element method that has so many advantages for the velocities and piecewise constant for the pressures in both the Stokes and Biot regions and apply a stabilization term penalizing the jumps over the element edges of the piecewise continuous velocities. Indeed, one can construct finite element methods where the incompressibility condition is exactly satisfied (cf. Fortin [33]), but this leads to the use of complex elements of limited applicability (e.g., oil and gas extraction for conforming case). Thus, in the work [32], Houédanou has constructed and studied the finite element method using simpler elements where the incompressibility condition is only approximatively satisfied (cf. definition of an operator divh (22)).So, in this paper, we have found it very convenient to use nonconforming finite elements which violate the interelement continuity condition of the velocities.In the work [34, 35], 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 ℝd. The approach utilizes a Lagrange multiplier method to impose weakly the interface conditions. To our best knowledge, there is no a posteriori error estimation for the strongly coupled mixed formulation (19) ([32], Section 3) of the coupled Stokes–Biot problem where a nonconforming finite element method is used. Here, we develop such a posteriori error analysis. The technique used to derive the a posteriori error estimator is the residue-based approach developed by Verfürth [36] in the late 80s and 90s. This paper could have several applications in dynamical systems, Wavelet analysis, and generalized harmonic analysis [37–39].This article is mainly composed of six sections. Some preliminaries and notations are given in Section2. In Section 3, the a posteriori error estimators are derived. The efficiency result is derived using the technique of bubble function introduced by Verfürth [36] and used in a similar context by Carstensen [40, 41]. These analytical tools are recalled in Section 4. The main results are given in Section 5. We offer our conclusion and the further works in Section 6.


## 2. Preliminaries and Notation

### 2.1. 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. The two regions are separated by an interface Γfp=∂Ωf∩∂Ωp. Let Γ⋆=∂Ω⋆∖Γfp, ⋆=f,p. Each interface and boundary is assumed to be polygonal d=2 or polyhedral d=3. We denote by nf (resp., np) the unit outward normal vector along ∂Ωf (resp., Ωp). Note that, on the interface Γfp, we have nf=−np. Figure 1 gives a schematic representation of the geometry.Figure 1
Global domainΩ consisting of the fluid region Ωf and the poroelastic media region Ωp separated by the interface Γfp.For any functionv defined in Ω, since its restriction to Ωf or Ωp could play a different mathematical role (for instance, their traces on Γfp), we will set vf=v|Ωf and vp=v|Ωp. In Ω, we denote by u the fluid velocity and by p the pressure and let ηp be the displacement in Ωp. Let μ>0 be the fluid viscosity, let f∈L2Ωd be the body force terms, and let g be the external source or sink terms satisfying the compatibility condition ∫Ωgxdx=0.LetDu and σfu,p denote, respectively, the deformation rate tensor and the stress tensor:(1)Du=12∇u+∇uT,andσfu,p=−pI+2μDu.Equations (1)–(13) consist of the model of the coupled Stokes and Biot flow problem that we will study in this paper.In the free fluid regionΩf, u,p satisfy the Stokes equations:(2)−∇⋅σfu,p=finΩf,(3)∇⋅u=ginΩf(4)u=0onΓf.Letσeηp and σpηp,pp be the elastic and poroelastic stress tensors, respectively:(5)σ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 modified static Biot system [27]:(6)−∇⋅σpηp,pp=finΩp,(7)μK−1u+∇p=0inΩp,(8)α∇⋅ηp+∇⋅u=ginΩp,(9)u⋅nd=0onΓp,(10)ηp=0onΓp.K is the symmetric and uniformly positive definite rock permeability tensor, satisfying, for some constants 0<kmin≤kmax,(11)∀ξ∈ℝd,kminξTξ≤ξTKxξ≤kmaxξTξ,∀x∈Ωp.Following [1], the interface conditions on the fluid-poroelasticity interface Γfp are mass conservation, the balance of stresses, and the Beavers–Joseph–Saffman (BJS) condition [42] modeling slip with friction are as follows:(12)uf⋅nf+up⋅np=0onΓfp,(13)σfnf+σpnp=0onΓfp,(14)−σfnf⋅nf=pp,onΓfp,(15)−σfnf⋅τf,j=μαBJFKj−1uf⋅τf,jonΓfp,where τ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. We note that continuity of flux constraints the normal velocity of the solid skeleton, while the BJS condition accounts for its tangential velocity.


### 2.2. Strongly Coupled Weak Formulation

We begin this section by introducing some useful notations. We first introduce some Sobolev spaces [43] 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 semi-norm on HmΩd are given by(16)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 λ, σ one defines(17)λ,σE≔∫Eλσds.For an open subsetF of the entire domain Ω, i.e., F⊆Ω, we define the space Hdiv;F by(18)Hdiv;F≔v∈L2Fd:divv∈L2F,with a norm(19)vHdiv;F≔vL2Fd2+divvL2F21/2,∀v∈Hdiv;F.To present a variational form of the coupled problem, we define the following three spaces for the velocityu, the structure displacement ηp, and the pressure:(20)H≔v∈Hdiv;Ω:vf∈H1Ωfd,v=0onΓf,v⋅np=0onΓp,equipped with the norm:(21)vH≔v1,f2+vHdiv;Ωp21/2,Xp≔ξp∈H1Ωpd:ξp=0onΓp,with the norm(22)ξpXp≔ξp1,p,(23)M≔02Ω×L02Ωp,equipped with the norm QM≔Q10,Ω2+Q20,Ωp21/2,∀Q=Q1,Q2∈M.Note that the vector-valued functions inH have (weakly) continuous normal components on Γfp (a consequence of Theorem 1.2.5 of [44], p. 27).We setℍ=H×Xp equipped with the product norm(24)Vℍ≔vH+ξpXp,∀V=v,ξp∈ℍ.Let us further introduce two bilinear forms:(25)A:ℍ×ℍ⟶ℝ,U,V↦AU,Vdefineby,AU,V≔2μDu,DvΩf+μK−1u,vΩp+2μpDηp,DξpΩp+λp∇⋅ηp,∇⋅ξpΩp+∑j=1d−1μαBJSKj−1uf⋅τf,j,vf⋅τf,jΓfp,B:ℍ×M⟶ℝ,V,Q↦BV,Qwith,BV,Q≔−Q1,divvΩ−αQ2,divξpΩp,whereQ=Q1,Q2,and two linear forms(26)L:ℍ⟶ℝ,V↦LV≔f,vΩ,(27)G:M⟶ℝ,Q=Q1,Q2↦GQ≔−g,Q1Ω.The weak formulation of the coupled problem (1)–(13) can be stated as follows: find U,P∈ℍ×M with U=u,ηp and P=p,pp such that(28)AU,V+BV,P=LV∀V=v,ξp∈ℍBU,Q=GQ∀Q=Q1,Q2∈M.Note that, iff and g are of mean zero, (28) directly implies that (1)–(11) hold (the differential equations being understood in the distributional sense), while the interface conditions (12) and (13) are imposed in a weak sense.This problem has a unique solution as proved in ([32], Theorem 3.1).Theorem 1.
Iff∈L2Ωd and g∈L02Ω, then there exists a unique solution U,P∈ℍ×M to problem (19).


### 2.3. Discontinuous Galerkin Discretization

In this section, we will use the nonconforming Crouzeix–Raviart piecewise linear finite element approximation for velocity and piecewise constant approximation for pressure and establish the existence and uniqueness of a finite element solution of the discrete problem.LetTh be a family of triangulations of Ω¯ with nondegenerate elements (i.e. triangles for d=2 and tetrahedrons for d=3). For any K∈Th, we denote by hK the diameter of K and ρK the diameter of the largest ball inscribed into K.We set(29)h=maxK∈ThhK,andσh=maxT∈ThhKρK.We assume that the family of triangulations is regular, in the sense that there existsσ0>0 such that σh⩽σ0, for all h>0. We also assume that the triangulation is conforming with respect to the partition of Ω into Ωf and Ωp, namely, each K∈Th is either in Ωf or in Ωp (see Figures 2–4 for illustration).Figure 2
Isotropic elementK in ℝ2.Figure 3
Example of conforming mesh inℝ2.Figure 4
Example of nonconforming mesh inℝ2.LetThf and Thp be the corresponding induced triangulations of Ωf and Ωp. For any K∈Th, we denote by EK (resp. NK) the set of its edges d=2 or faces d=3 (resp. vertices) and set Εh=∪K∈ThEK, Nh=∪K∈ThNK. For A⊂Ω¯, we define(30)ΕhA=E∈Εh:E⊂A.Notice thatEh can be split up in the form(31)Eh=EhΩf+∪EhΩp∪Eh∂Ωp,where Ωf+=Ωf∪Γf. Note that EhΓfp is included in Eh∂Ωp.With every edgeE∈Eh, we associate a unit vector nE such that nE is orthogonal to E and equals to the unit exterior normal vector to ∂Ω if E⊂∂Ω. For any E∈Eh and any piecewise continuous function φ, we denote by φE its jump across E in the direction of nE:(32)φEx≔limt⟶0+φx+tnE−limt⟶0+φx−tnE,for an interioredge/faceE,−limt⟶0+φx−tnE,for a boundaryedge/faceE.Based on the above notations, we introduce a variant of the nonconforming Crouzeix–Raviart piecewise linear finite element space:(33)Hh≔vh∈L2Ωd:vh|K∈ℙ1Kd∀K∈Th,vhE,1E=0∀E∈ΕhΩf+,vh⋅nEE,1E=00.2cm∀E∈ΕhΩp∪Εh∂Ωp,Xph≔ξph∈L2Ωpd:ξph|K∈ℙ1Kd∀K∈Thp,ξphE,1E=0∀E∈ΕhΩ¯p.ForX⊆Ω, we set(34)EhX≔qh∈L02X:qh|K∈ℙ0K∀K⊂X,K∈Th,and we define(35)Mh≔EhΩ×EhΩp⊂M,ℍh≔Hh×Xph⊈ℍ,where ℙmK is the space of the restrictions to K of all polynomials of degree less than or equal to m.The spaceMh is equipped with the norm ∥⋅∥M while the norm on ℍh will be specified later on. The choice of Hh is more natural since the space Hh approximates only Hdiv;Ωp and not H1Ωpd, while our a posteriori error analysis is only valid in this larger space.Let us introduce the discrete divergence operatordivh∈ℒHh;EhΩ∩ℒH;L02Ω by(36)divhvh|K=divvh|K,∀K∈Th,or divh∈ℒXph;EhΩp∩ℒXp;L02Ωp by(37)divhξph|K=divξph|K,∀K∈Thp.Then, forUh=uh,ηph∈ℍh, Vh=vh,ξph∈ℍh and Qh=Q1h,Q2h∈Mh, we can introduce the following two bilinear forms:(38)AhUh,Vh≔∑K∈Thf2μDuh,DvhK+μK−1uh,vhΩp+∑K∈Thp2μpDηph,DξphK+λpdivhηph,divhξphΩp+∑j=1d−1μαBJSKj−1ufh⋅τf,j,vfh⋅τf,jΓfp,BhVh,Qh≔−Q1h,divhvhΩ−αQ2h,divhξphΩp.Then, we propose the following discrete problem: findUh,Ph∈ℍh×Qh with Uh=uh,ηph∈Hh×Xph and Ph=ph,pph∈EhΩ×EhΩp such that(39)AhUh,Vh+BhVh,Ph+JUh,Vh=,LVh∀Vh∈ℍh,BhUh,Qh=,GQh∀Qh∈Mh.This is the natural discretization of the weak formulation (19) only with the penalizing term JUh,Vh added, where Vh=vh,ξph. We define the bilinear form J⋅,⋅ following the decomposition of Eh:(40)JUh,Vh=JΩf+uh,vh+JΩpuh,vh+J∂Ωpuh,vh+JΩp+ηph,ξph,where(41)JΩf+uh,vh≔1+2μ∑E∈ΕhΩf+hE−1∫EuhE⋅vhEds,JΩpuh,vh≔∑E∈ΕhΩphE−1∫EuhE⋅vhEds,J∂Ωpuh,vh≔∑E∈Εh∂ΩphE−1∫Euh⋅nEEvh⋅nEEds0.2cmandJΩp+ηph,ξph≔∑E∈ΕhΩp+hE−1∫E1+2μpηphE⋅ξphEds.Here,hE is the length d=2 or diameter d=3 of E. Note that each element of Eh only contributes with one jump term in JUh,Vh.We are now able to define the norm onℍh:(42)Vhℍh≔vhHh2+ξphXph2+JVh,Vh1/2,where(43)vhHh≔∑K∈Thfvh1,K2+∑j=1d−1vfh⋅τj,vfh⋅τjΓfp+vhΩp2+divhvhΩp21/2,(44)ξphXph≔∑K∈Thpξph1,K21/2.The following results hold (Theorems 4.1 and 5.1 in32).Theorem 2.
There exists a unique solutionUh,Ph∈ℍh×Mh to a discrete problem (24) and if the solution U,P∈ℍ×M of the continuous problem (19) is smooth enough, then we have(45)U−Uhℍh∪ℍ+P−PhM≲hu2,Ωf+u2,Ωp+ηp2,Ωp+p1,Ωf+p1,Ωp.
Here and below, in order to avoid excessive use of constants, the abbreviationx≲y stands for x≤cy, with c a positive constant independent of x, y and Th.


## 2.1. 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. The two regions are separated by an interface Γfp=∂Ωf∩∂Ωp. Let Γ⋆=∂Ω⋆∖Γfp, ⋆=f,p. Each interface and boundary is assumed to be polygonal d=2 or polyhedral d=3. We denote by nf (resp., np) the unit outward normal vector along ∂Ωf (resp., Ωp). Note that, on the interface Γfp, we have nf=−np. Figure 1 gives a schematic representation of the geometry.Figure 1
Global domainΩ consisting of the fluid region Ωf and the poroelastic media region Ωp separated by the interface Γfp.For any functionv defined in Ω, since its restriction to Ωf or Ωp could play a different mathematical role (for instance, their traces on Γfp), we will set vf=v|Ωf and vp=v|Ωp. In Ω, we denote by u the fluid velocity and by p the pressure and let ηp be the displacement in Ωp. Let μ>0 be the fluid viscosity, let f∈L2Ωd be the body force terms, and let g be the external source or sink terms satisfying the compatibility condition ∫Ωgxdx=0.LetDu and σfu,p denote, respectively, the deformation rate tensor and the stress tensor:(1)Du=12∇u+∇uT,andσfu,p=−pI+2μDu.Equations (1)–(13) consist of the model of the coupled Stokes and Biot flow problem that we will study in this paper.In the free fluid regionΩf, u,p satisfy the Stokes equations:(2)−∇⋅σfu,p=finΩf,(3)∇⋅u=ginΩf(4)u=0onΓf.Letσeηp and σpηp,pp be the elastic and poroelastic stress tensors, respectively:(5)σ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 modified static Biot system [27]:(6)−∇⋅σpηp,pp=finΩp,(7)μK−1u+∇p=0inΩp,(8)α∇⋅ηp+∇⋅u=ginΩp,(9)u⋅nd=0onΓp,(10)ηp=0onΓp.K is the symmetric and uniformly positive definite rock permeability tensor, satisfying, for some constants 0<kmin≤kmax,(11)∀ξ∈ℝd,kminξTξ≤ξTKxξ≤kmaxξTξ,∀x∈Ωp.Following [1], the interface conditions on the fluid-poroelasticity interface Γfp are mass conservation, the balance of stresses, and the Beavers–Joseph–Saffman (BJS) condition [42] modeling slip with friction are as follows:(12)uf⋅nf+up⋅np=0onΓfp,(13)σfnf+σpnp=0onΓfp,(14)−σfnf⋅nf=pp,onΓfp,(15)−σfnf⋅τf,j=μαBJFKj−1uf⋅τf,jonΓfp,where τ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. We note that continuity of flux constraints the normal velocity of the solid skeleton, while the BJS condition accounts for its tangential velocity.


## 2.2. Strongly Coupled Weak Formulation

We begin this section by introducing some useful notations. We first introduce some Sobolev spaces [43] 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 semi-norm on HmΩd are given by(16)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 λ, σ one defines(17)λ,σE≔∫Eλσds.For an open subsetF of the entire domain Ω, i.e., F⊆Ω, we define the space Hdiv;F by(18)Hdiv;F≔v∈L2Fd:divv∈L2F,with a norm(19)vHdiv;F≔vL2Fd2+divvL2F21/2,∀v∈Hdiv;F.To present a variational form of the coupled problem, we define the following three spaces for the velocityu, the structure displacement ηp, and the pressure:(20)H≔v∈Hdiv;Ω:vf∈H1Ωfd,v=0onΓf,v⋅np=0onΓp,equipped with the norm:(21)vH≔v1,f2+vHdiv;Ωp21/2,Xp≔ξp∈H1Ωpd:ξp=0onΓp,with the norm(22)ξpXp≔ξp1,p,(23)M≔02Ω×L02Ωp,equipped with the norm QM≔Q10,Ω2+Q20,Ωp21/2,∀Q=Q1,Q2∈M.Note that the vector-valued functions inH have (weakly) continuous normal components on Γfp (a consequence of Theorem 1.2.5 of [44], p. 27).We setℍ=H×Xp equipped with the product norm(24)Vℍ≔vH+ξpXp,∀V=v,ξp∈ℍ.Let us further introduce two bilinear forms:(25)A:ℍ×ℍ⟶ℝ,U,V↦AU,Vdefineby,AU,V≔2μDu,DvΩf+μK−1u,vΩp+2μpDηp,DξpΩp+λp∇⋅ηp,∇⋅ξpΩp+∑j=1d−1μαBJSKj−1uf⋅τf,j,vf⋅τf,jΓfp,B:ℍ×M⟶ℝ,V,Q↦BV,Qwith,BV,Q≔−Q1,divvΩ−αQ2,divξpΩp,whereQ=Q1,Q2,and two linear forms(26)L:ℍ⟶ℝ,V↦LV≔f,vΩ,(27)G:M⟶ℝ,Q=Q1,Q2↦GQ≔−g,Q1Ω.The weak formulation of the coupled problem (1)–(13) can be stated as follows: find U,P∈ℍ×M with U=u,ηp and P=p,pp such that(28)AU,V+BV,P=LV∀V=v,ξp∈ℍBU,Q=GQ∀Q=Q1,Q2∈M.Note that, iff and g are of mean zero, (28) directly implies that (1)–(11) hold (the differential equations being understood in the distributional sense), while the interface conditions (12) and (13) are imposed in a weak sense.This problem has a unique solution as proved in ([32], Theorem 3.1).Theorem 1.
Iff∈L2Ωd and g∈L02Ω, then there exists a unique solution U,P∈ℍ×M to problem (19).


## 2.3. Discontinuous Galerkin Discretization

In this section, we will use the nonconforming Crouzeix–Raviart piecewise linear finite element approximation for velocity and piecewise constant approximation for pressure and establish the existence and uniqueness of a finite element solution of the discrete problem.LetTh be a family of triangulations of Ω¯ with nondegenerate elements (i.e. triangles for d=2 and tetrahedrons for d=3). For any K∈Th, we denote by hK the diameter of K and ρK the diameter of the largest ball inscribed into K.We set(29)h=maxK∈ThhK,andσh=maxT∈ThhKρK.We assume that the family of triangulations is regular, in the sense that there existsσ0>0 such that σh⩽σ0, for all h>0. We also assume that the triangulation is conforming with respect to the partition of Ω into Ωf and Ωp, namely, each K∈Th is either in Ωf or in Ωp (see Figures 2–4 for illustration).Figure 2
Isotropic elementK in ℝ2.Figure 3
Example of conforming mesh inℝ2.Figure 4
Example of nonconforming mesh inℝ2.LetThf and Thp be the corresponding induced triangulations of Ωf and Ωp. For any K∈Th, we denote by EK (resp. NK) the set of its edges d=2 or faces d=3 (resp. vertices) and set Εh=∪K∈ThEK, Nh=∪K∈ThNK. For A⊂Ω¯, we define(30)ΕhA=E∈Εh:E⊂A.Notice thatEh can be split up in the form(31)Eh=EhΩf+∪EhΩp∪Eh∂Ωp,where Ωf+=Ωf∪Γf. Note that EhΓfp is included in Eh∂Ωp.With every edgeE∈Eh, we associate a unit vector nE such that nE is orthogonal to E and equals to the unit exterior normal vector to ∂Ω if E⊂∂Ω. For any E∈Eh and any piecewise continuous function φ, we denote by φE its jump across E in the direction of nE:(32)φEx≔limt⟶0+φx+tnE−limt⟶0+φx−tnE,for an interioredge/faceE,−limt⟶0+φx−tnE,for a boundaryedge/faceE.Based on the above notations, we introduce a variant of the nonconforming Crouzeix–Raviart piecewise linear finite element space:(33)Hh≔vh∈L2Ωd:vh|K∈ℙ1Kd∀K∈Th,vhE,1E=0∀E∈ΕhΩf+,vh⋅nEE,1E=00.2cm∀E∈ΕhΩp∪Εh∂Ωp,Xph≔ξph∈L2Ωpd:ξph|K∈ℙ1Kd∀K∈Thp,ξphE,1E=0∀E∈ΕhΩ¯p.ForX⊆Ω, we set(34)EhX≔qh∈L02X:qh|K∈ℙ0K∀K⊂X,K∈Th,and we define(35)Mh≔EhΩ×EhΩp⊂M,ℍh≔Hh×Xph⊈ℍ,where ℙmK is the space of the restrictions to K of all polynomials of degree less than or equal to m.The spaceMh is equipped with the norm ∥⋅∥M while the norm on ℍh will be specified later on. The choice of Hh is more natural since the space Hh approximates only Hdiv;Ωp and not H1Ωpd, while our a posteriori error analysis is only valid in this larger space.Let us introduce the discrete divergence operatordivh∈ℒHh;EhΩ∩ℒH;L02Ω by(36)divhvh|K=divvh|K,∀K∈Th,or divh∈ℒXph;EhΩp∩ℒXp;L02Ωp by(37)divhξph|K=divξph|K,∀K∈Thp.Then, forUh=uh,ηph∈ℍh, Vh=vh,ξph∈ℍh and Qh=Q1h,Q2h∈Mh, we can introduce the following two bilinear forms:(38)AhUh,Vh≔∑K∈Thf2μDuh,DvhK+μK−1uh,vhΩp+∑K∈Thp2μpDηph,DξphK+λpdivhηph,divhξphΩp+∑j=1d−1μαBJSKj−1ufh⋅τf,j,vfh⋅τf,jΓfp,BhVh,Qh≔−Q1h,divhvhΩ−αQ2h,divhξphΩp.Then, we propose the following discrete problem: findUh,Ph∈ℍh×Qh with Uh=uh,ηph∈Hh×Xph and Ph=ph,pph∈EhΩ×EhΩp such that(39)AhUh,Vh+BhVh,Ph+JUh,Vh=,LVh∀Vh∈ℍh,BhUh,Qh=,GQh∀Qh∈Mh.This is the natural discretization of the weak formulation (19) only with the penalizing term JUh,Vh added, where Vh=vh,ξph. We define the bilinear form J⋅,⋅ following the decomposition of Eh:(40)JUh,Vh=JΩf+uh,vh+JΩpuh,vh+J∂Ωpuh,vh+JΩp+ηph,ξph,where(41)JΩf+uh,vh≔1+2μ∑E∈ΕhΩf+hE−1∫EuhE⋅vhEds,JΩpuh,vh≔∑E∈ΕhΩphE−1∫EuhE⋅vhEds,J∂Ωpuh,vh≔∑E∈Εh∂ΩphE−1∫Euh⋅nEEvh⋅nEEds0.2cmandJΩp+ηph,ξph≔∑E∈ΕhΩp+hE−1∫E1+2μpηphE⋅ξphEds.Here,hE is the length d=2 or diameter d=3 of E. Note that each element of Eh only contributes with one jump term in JUh,Vh.We are now able to define the norm onℍh:(42)Vhℍh≔vhHh2+ξphXph2+JVh,Vh1/2,where(43)vhHh≔∑K∈Thfvh1,K2+∑j=1d−1vfh⋅τj,vfh⋅τjΓfp+vhΩp2+divhvhΩp21/2,(44)ξphXph≔∑K∈Thpξph1,K21/2.The following results hold (Theorems 4.1 and 5.1 in32).Theorem 2.
There exists a unique solutionUh,Ph∈ℍh×Mh to a discrete problem (24) and if the solution U,P∈ℍ×M of the continuous problem (19) is smooth enough, then we have(45)U−Uhℍh∪ℍ+P−PhM≲hu2,Ωf+u2,Ωp+ηp2,Ωp+p1,Ωf+p1,Ωp.
Here and below, in order to avoid excessive use of constants, the abbreviationx≲y stands for x≤cy, with c a positive constant independent of x, y and Th.


## 3. Error Estimators

In order to solve the Stokes–Biot coupled problem by efficient adaptive finite element methods, reliable and efficient a posteriori error analysis is important to provide appropriate indicators. In this section, we first define the local and global indicators, and then the lower and upper error bounds are derived (see Sections5.2.1 and 5.2.2).


### 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, the exact element residuals are denoted by(46)R1=f+∇⋅σfuh,phinK∈Thf,R2=g−∇⋅uhinK∈Thf,R3=μK−1uh+∇phinK∈Thp,R4=f+∇⋅σp(ηph,pph)) inK∈Thp,R5=g−α∇⋅ηph−∇⋅uhinK∈Thp,R6j=μαBJFKj−1ufh⋅τf,j+σfuh,phnf⋅τf,jonE∈Εh∂K∩Γ¯fp,R7=pph+σfufh,phnf⋅nfonE∈Εh∂K∩Γ¯fp.As it is common, these exact residuals are replaced by some finite-dimensional approximation called approximate element residualri,K, i∈1,4:(47)ri,K∈ℙmKdonK∈Thl,l∈f,p.This approximation is here achieved by projectingf on the space of piecewise constant functions in Ωl, more precisely for all K∈Th, we take(48)fK,l=1K∫Kfxdx,l∈f,p,∀K∈Thl.Finally, the global functionfh is defined by(49)fh,l=fK,linK,∀K∈Thl.Hence,(50)r1,K≔fK,f+∇⋅σfuh,phinK∈Thf,r4,K≔fK,p+∇⋅σpηph,pphinK∈Thp,with ul,h≔uh|Ωl and pl,h≔ph|Ωl,l=f,p.Next, the gradient jump in a normal direction is introduced by the following:(51)JE,nE≔σfuh,ph⋅nEE,for an interioredge/faceE∈ΕhΩf,0,for a  boundaryedge/faceE∈ΕhΓf,(52)GE,nE≔σpηp,h,pp,h⋅nEE,for an interioredge/faceE∈EhΩp,0,for a boundaryedge/faceE∈EhΓp.Definition 1.
(residual error estimator). The residual error estimator is locally defined by(53)ΥK≔∑i=111Υi,K21/2for eachK∈Th,where(54)Υ1,K2≔hK2r1,KK2ifK∈Thf,hK2∥r4,K∥K2ifK∈Thp,Υ2,K2≔R3K2ifK∈Thp,0ifK∈Thf,Υ3,K2≔curlR3K2ifK∈Thp,0ifK∈Thf,Υ4,K2≔R2K2,ifK∈Thf,R5K2,ifK∈Thp,Υ5,K2≔∑E∈Εh∂K∩Γ¯fphE∑j=1d−1R6jE2,Υ6,K2≔∑E∈Εh∂K∩Γ¯fphER7E2Υ7,K2≔∑E∈Εh∂K∩Ω¯fhEJE,nEE2,ifK∈Thf,∑E∈Εh∂K∩Ω¯phEGE,nEE2+phEE2,ifK∈Thp,Υ8,K2≔∑E∈Εh∂K∩ΩphE−1uhEE2,Υ9,K2≔∑E∈Εh∂K∩∂ΩphE−1uh⋅nEEE2,Υ10,K2≔∑E∈Εh∂K∩Ωf+hE−11+2μuhEE2Υ11,K2≔∑E∈Εh∂K∩Ωp+hE−11+2μpηphEE2.
The global residual error estimator is given by(55)Υ≔∑K∈ThΥK21/2.
Furthermore, the local and global approximation terms are denoted by(56)ΨK≔hKf−fhK,∀K∈Th,,(57)Ψ≔∑K∈ThΨK21/2.


## 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, the exact element residuals are denoted by(46)R1=f+∇⋅σfuh,phinK∈Thf,R2=g−∇⋅uhinK∈Thf,R3=μK−1uh+∇phinK∈Thp,R4=f+∇⋅σp(ηph,pph)) inK∈Thp,R5=g−α∇⋅ηph−∇⋅uhinK∈Thp,R6j=μαBJFKj−1ufh⋅τf,j+σfuh,phnf⋅τf,jonE∈Εh∂K∩Γ¯fp,R7=pph+σfufh,phnf⋅nfonE∈Εh∂K∩Γ¯fp.As it is common, these exact residuals are replaced by some finite-dimensional approximation called approximate element residualri,K, i∈1,4:(47)ri,K∈ℙmKdonK∈Thl,l∈f,p.This approximation is here achieved by projectingf on the space of piecewise constant functions in Ωl, more precisely for all K∈Th, we take(48)fK,l=1K∫Kfxdx,l∈f,p,∀K∈Thl.Finally, the global functionfh is defined by(49)fh,l=fK,linK,∀K∈Thl.Hence,(50)r1,K≔fK,f+∇⋅σfuh,phinK∈Thf,r4,K≔fK,p+∇⋅σpηph,pphinK∈Thp,with ul,h≔uh|Ωl and pl,h≔ph|Ωl,l=f,p.Next, the gradient jump in a normal direction is introduced by the following:(51)JE,nE≔σfuh,ph⋅nEE,for an interioredge/faceE∈ΕhΩf,0,for a  boundaryedge/faceE∈ΕhΓf,(52)GE,nE≔σpηp,h,pp,h⋅nEE,for an interioredge/faceE∈EhΩp,0,for a boundaryedge/faceE∈EhΓp.Definition 1.
(residual error estimator). The residual error estimator is locally defined by(53)ΥK≔∑i=111Υi,K21/2for eachK∈Th,where(54)Υ1,K2≔hK2r1,KK2ifK∈Thf,hK2∥r4,K∥K2ifK∈Thp,Υ2,K2≔R3K2ifK∈Thp,0ifK∈Thf,Υ3,K2≔curlR3K2ifK∈Thp,0ifK∈Thf,Υ4,K2≔R2K2,ifK∈Thf,R5K2,ifK∈Thp,Υ5,K2≔∑E∈Εh∂K∩Γ¯fphE∑j=1d−1R6jE2,Υ6,K2≔∑E∈Εh∂K∩Γ¯fphER7E2Υ7,K2≔∑E∈Εh∂K∩Ω¯fhEJE,nEE2,ifK∈Thf,∑E∈Εh∂K∩Ω¯phEGE,nEE2+phEE2,ifK∈Thp,Υ8,K2≔∑E∈Εh∂K∩ΩphE−1uhEE2,Υ9,K2≔∑E∈Εh∂K∩∂ΩphE−1uh⋅nEEE2,Υ10,K2≔∑E∈Εh∂K∩Ωf+hE−11+2μuhEE2Υ11,K2≔∑E∈Εh∂K∩Ωp+hE−11+2μpηphEE2.
The global residual error estimator is given by(55)Υ≔∑K∈ThΥK21/2.
Furthermore, the local and global approximation terms are denoted by(56)ΨK≔hKf−fhK,∀K∈Th,,(57)Ψ≔∑K∈ThΨK21/2.


## 4. Analytical Tools

### 4.1. Some Technical Results

Our a posteriori analysis requires some analytical results that are recalled.The first one concerns a sort of Helmholtz decomposition of elements ofH. Recall first that if d=3,(58)H0curl,Ωp=ψ∈L2Ωp3:curlψ∈L2Ωp3andψ×n=0on∂Ωp.Theorem 3.
(see reference [2], p. 708). Any v∈H admits the Helmholtz-type decomposition(59)v=v0+v1,where v0,v1∈H but satisfying v0∈H1Ωd,(60)v1=0,in,Ωf,curlβp,in,Ωp,where βp∈H01Ωp if d=2, while βp∈H1Ωp3∩H0curl,Ωp, if d=3, with the estimate(61)v01,Ω+βp1,Ωp≲vH.
The second result that we need is a regularity result for the solutionu∈H of (28).Theorem 4.
(see [2], p. 710). Let U,P∈ℍ×M be the unique solution of (19) with U=u,ηp∈H×Xp. If K∈C0,1Ω¯pd×d, then there exists δ>0 such that(62)u|Ωp∈H12+δΩpd.
Note that the regularity ofu∈H1/2+δΩpd, with δ>0 allows to give a meaning to JΩpu,w+J∂Ωpu,w for all w∈H∪Hh and hence to show that JU,W=0 for all W=w,ξp∈ℍ∪ℍh.
Let us finish this section by an estimation of the nonconformity error (see [2], Theorem 3.3).Theorem 5.
For anyUh=uh,ηph∈ℍh, we have(63)infWh∈ℍh∩ℍUh−Whℍh2≲JUh,Uh.


### 4.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:(64)PcbTh≔v∈C0Ω¯:v|K∈ℙ1K,∀K∈Thandv=0on∂Ω,In addition, we will make use of a vector-valued version ofICl0, that is, ICl0:H01Ωd⋅PcbThd, which is defined componentwise by ICl0. The following lemma establishes the local approximation properties of ICl0 (and hence of ICl0). (for a proof, see [45], Section 3).Lemma 1.
There exist constantsC1,C2>0, independent of h, such that for all v∈H01Ω there hold(65)v−ICl0vK≤C1hKv1,ΔK,∀K∈Th,0.2cmandv−ICl0vE≤C2hE1/2v1,ΔE,∀E∈Εh,where ΔK≔∪K′∈Th:K′∩K≠∅ and ΔE≔∪K′∈Th:K′∩E≠∅.


### 4.3. Inverse Inequalities

In order to derive the lower error bounds, we proceed similarly as in [40, 41], 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, respectively (see (1.5) and (1.6) in [46]). In particular, bK satisfies bK∈ℙ3K, suppbK⊆K, bK=0on∂K, and 0≤bK≤1onK. Similarly, bE∈ℙ2K, suppbE⊆ωE:=K′∈Th:E∈EK′, bE=0on∂K∖E and 0≤bE≤1inω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 componentwise application of L, is denoted by L. Additional properties of bK, bE and L are collected in the following lemma (see [36]).Lemma 2.
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(66)qK≲qbK1/2K≲qK,∀q∈ℙkK,qbK1,K≲hK−1qK,∀q∈ℙkK,pE≲bE1/2pE≲pE,∀p∈ℙkE,LpK+hELp1,K≲hE1/2pE,∀p∈ℙkE.Lemma 3 (continuous trace inequality).
There exists a positive constantβ1>0 depending only on σ0 such that(67)v∂K2≤β1vKv1,K,∀K∈Th,∀v∈H1Kd.


## 4.1. Some Technical Results

Our a posteriori analysis requires some analytical results that are recalled.The first one concerns a sort of Helmholtz decomposition of elements ofH. Recall first that if d=3,(58)H0curl,Ωp=ψ∈L2Ωp3:curlψ∈L2Ωp3andψ×n=0on∂Ωp.Theorem 3.
(see reference [2], p. 708). Any v∈H admits the Helmholtz-type decomposition(59)v=v0+v1,where v0,v1∈H but satisfying v0∈H1Ωd,(60)v1=0,in,Ωf,curlβp,in,Ωp,where βp∈H01Ωp if d=2, while βp∈H1Ωp3∩H0curl,Ωp, if d=3, with the estimate(61)v01,Ω+βp1,Ωp≲vH.
The second result that we need is a regularity result for the solutionu∈H of (28).Theorem 4.
(see [2], p. 710). Let U,P∈ℍ×M be the unique solution of (19) with U=u,ηp∈H×Xp. If K∈C0,1Ω¯pd×d, then there exists δ>0 such that(62)u|Ωp∈H12+δΩpd.
Note that the regularity ofu∈H1/2+δΩpd, with δ>0 allows to give a meaning to JΩpu,w+J∂Ωpu,w for all w∈H∪Hh and hence to show that JU,W=0 for all W=w,ξp∈ℍ∪ℍh.
Let us finish this section by an estimation of the nonconformity error (see [2], Theorem 3.3).Theorem 5.
For anyUh=uh,ηph∈ℍh, we have(63)infWh∈ℍh∩ℍUh−Whℍh2≲JUh,Uh.


## 4.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:(64)PcbTh≔v∈C0Ω¯:v|K∈ℙ1K,∀K∈Thandv=0on∂Ω,In addition, we will make use of a vector-valued version ofICl0, that is, ICl0:H01Ωd⋅PcbThd, which is defined componentwise by ICl0. The following lemma establishes the local approximation properties of ICl0 (and hence of ICl0). (for a proof, see [45], Section 3).Lemma 1.
There exist constantsC1,C2>0, independent of h, such that for all v∈H01Ω there hold(65)v−ICl0vK≤C1hKv1,ΔK,∀K∈Th,0.2cmandv−ICl0vE≤C2hE1/2v1,ΔE,∀E∈Εh,where ΔK≔∪K′∈Th:K′∩K≠∅ and ΔE≔∪K′∈Th:K′∩E≠∅.


## 4.3. Inverse Inequalities

In order to derive the lower error bounds, we proceed similarly as in [40, 41], 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, respectively (see (1.5) and (1.6) in [46]). In particular, bK satisfies bK∈ℙ3K, suppbK⊆K, bK=0on∂K, and 0≤bK≤1onK. Similarly, bE∈ℙ2K, suppbE⊆ωE:=K′∈Th:E∈EK′, bE=0on∂K∖E and 0≤bE≤1inω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 componentwise application of L, is denoted by L. Additional properties of bK, bE and L are collected in the following lemma (see [36]).Lemma 2.
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(66)qK≲qbK1/2K≲qK,∀q∈ℙkK,qbK1,K≲hK−1qK,∀q∈ℙkK,pE≲bE1/2pE≲pE,∀p∈ℙkE,LpK+hELp1,K≲hE1/2pE,∀p∈ℙkE.Lemma 3 (continuous trace inequality).
There exists a positive constantβ1>0 depending only on σ0 such that(67)v∂K2≤β1vKv1,K,∀K∈Th,∀v∈H1Kd.


## 5. Main Results

We setX:=ℍ×M and Xh:=ℍh×Mh and define on X, the continuous bilinear form B by(68)BU,W≔AU,V+BV,P+BU,QforU=U,PandforW=V,Q.We also define on the discrete spaceℍh, the form,(69)BhUh,Wh≔AhUh,Vh+BhVh,Ph+BhUh,Qh+JUh,VhforUh=Uh,Ph,Wh=Vh,Qh.The spacesX and Xh are equipped with the product norms:(70)U,P=Uℍ+PMandUh,Phh=Uhℍh∪ℍ+PhMrespectively.To prove local efficiency forω⊂Ω, let us denote by(71)W,Qh,ω2=∑K⊂ω¯∩Ω¯fv1,K2+∑K⊂ω¯∩Ω¯pξp1,K2+∑K⊂ω¯∩Ω¯pvK2+divhvK2+vf×nΓfp∩ω¯2+∑K⊂ω¯JKW,W+Qω2,∀W,Q=v,ξp,Q∈Hh∪H×Xp∪Xph×M,where(72)JKW,W=1+2μ∑E∈ΕhΩf+∩ΕKhE−1vEE2+∑E∈ΕhΩp∩ΕKhE−1vEE2+∑E∈Εh∂Ωp∩ΕKhE−1v⋅nEEE2+∑E∈Εh∂Ωp∩ΕKhE−11+μpξpEE2.


### 5.1. Optimality Result

The main result of this paper can be stated as follows:(1)
Reliability ofΥKK∈Th: the a posteriori error estimator Υ is consider reliable if it satisfies(73)U−Uh,P−Phh≲Υ+Ψ.(2)
Efficiency ofΥKK∈Th:Under the assumptions of Theorem4, the following lower error bound holds(74)ΥK≲U−Uh,P−Phh,w˜K+∑K′∈w˜KΨK′,where w˜K is a finite union of neighboring elements of K.


### 5.2. Proof of the A Posteriori Error Estimate

#### 5.2.1. Proof of the Reliability Estimate

In this section, we shall prove estimate (60).Let us start with the following result.Lemma 4.
Let the assumptions of Theorem4 be satisfied. Then, for all W=V,Q∈ℍ×M, we have the estimate(75)BhU−Uh,W≲Υ˜+ΨWh,where the estimator Υ˜ is defined by(76)Υ˜≔∑K∈Th∑i=17Υi,K212.Proof.
LetW=V,Q∈ℍ×M with V=v,ξp and Q=Q1,Q2. By Theorem 3, v admits the decomposition (49) with v0,v1∈H and satisfies the properties stated in Theorem 3. Then we take Wh=Vh,0∈ℍh×Mh where Vh=vh,ξph with vh=v0,h+v1,h, where v0,h=ICl0v0 and(77)v1,h=0,in,Ωf,curlICl0ψ,in,Ωp,andξph=ICl0ξp.
In2D, ICl0ψ is the standard Clément interpolant of ψ, while in 3D, we take the vectorial Clément interpolant from [2] that satisfies the same estimate as the standard one (see [2]). Note that v0,h belongs to Hh∩H01Ωd while v1,h simply belongs to H∩Hh (ICl0ψ being in H01Ωp if d=2 and ICl0ψ∈H1Ωp3∩H0curl,Ωp if d=3), its curl belongs to H0div,Ωp, hence v1,h, its extension by zero in Ωf, stays in H0div,Ω. With these definitions and noticing that divv−vh=divv0−v0,h and that JUh,Vh=0, we may write(78)BhU−Uh,W=BhU−Uh,W−Wh=AhU−Uh,V−Vh+BhV−Vh,P−Ph+BhU−Uh,Q=AhU,V−Vh+BhV−Vh,P+BhU,Q−AhUh,V−Vh+BhV−Vh,Ph+BhUh,Q=LV−Vh+GQ−AhUh,V−Vh+BhV−Vh,Ph+BhUh,Q=f,v−vhΩ−g,Q1Ω−AhUh,V−Vh+BhV−Vh,Ph+BhUh,Q=∑K∈Thf,v−vhK−g,Q1K−AhUh,V−Vh−BhV−Vh,Ph−BhUh,Q=∑K∈Thf,v−vhK−∑K∈Thfg,Q1K−∑K∈Thf2μDuh,Dv−vhK−∑K∈ThμK−1uh,v−vhK−∑K∈Thp2μpDηph,Dξ−ξphK−∑K∈Thpλpdivhηph,divξp−ξphK−∑j=1d−1μαBJSKj−1ufh⋅τf,j,vf−vfh⋅τf,jΓfp+∑K∈Thph,divv0−v0hK+∑K∈Thpαpph,divξp−ξphK+∑K∈ThQ,divuhK+∑K∈ThpαQ2,divηphK.
Integrate by parts element by element and add boundary (resp. internal) terms that appear on the same edge (resp. at the same element) and reminding thatv=v0 and vh=v0,h in Ωf, we obtain(79)BhU−Uh,W=∑K∈ThfR1,v−vhK−∑K∈ThpR3,v−vhK+∑K∈ThpR4,ξp−ξphK−∑K∈ThfR2,Q1K−∑K∈ThpR5,Q2K+∑E∈ΕhΓfpR7,v0−v0h⋅nfE−∑E∈ΕhΓfp∑j=1d−1(R6j,v0−v0h⋅τj)E−∑E∈ΕhΩf+JE,nE,v0−v0hE−∑E∈ΕhΩp+GE,nE,ξp−ξphE+∑E∈ΕhΩppphE,v0−v0,h⋅nEE.
We now introduce the approximationfh of f for appropriated terms and we have(80)BhU−Uh,W=∑K∈Thfr1,v−vhK−∑K∈ThpR3,v0−v0,hK+∑K∈Thpr4,ξp−ξphK−∑K∈ThfR2,Q1K−∑K∈ThpR5,Q2K+∑E∈ΕhΓfpR7,v0−v0h⋅nfE−∑E∈ΕhΓfp∑j=1d−1(R6j,v0−v0h⋅τj)E−∑E∈ΕhΩf+JE,nE,v0−v0hE−∑E∈ΕhΩp+GE,nE,ξp−ξphE+∑E∈ΕhΩppphE,v0−v0,h⋅nEE+∑K∈Thf−fh,v−vhK−∑K∈ThpR3,v1−v1,hK.
Now for a triangleK∈Thp, we recall that(81)v1−v1,h=curlψ−ICl0ψinK,and use Green’s formula to get(82)∑K∈ThpR3,v1−v1,hK=∑K∈ThpR3,curlψ−ICl0ψK−∑K∈Thp[curlR3,ψ−ICl0ψK+curlR3×n,ψ−ICl0ψ∂K].
We deduce the error equation,(83)BhU−Uh,W=∑K∈Thfr1,v−vhK−∑K∈ThpR3,v0−v0,hK+∑K∈Thpr4,ξp−ξphK−∑K∈ThfR2,Q1K−∑K∈ThpR5,Q2K+∑E∈ΕhΓfpR7,v0−v0h⋅nfE−∑E∈ΕhΓfp∑j=1d−1R6j,v0−v0h⋅τjE−∑E∈Εh(Ωf+)JE,nE,v0−v0hE−∑E∈ΕhΩp+GE,nE,ξp−ξphE+∑E∈ΕhΩppphE,v0−v0,h⋅nEE+∑K∈Thf−fh,v−vhK+∑K∈Thp[curlR3,ψ−ICl0ψK−curlR3×n,ψ−ICl0ψ∂K].
Cauchy–Schwarz inequality and the approximation properties of Lemma1 imply the required estimate and finish the proof.
The second result of this section is given by the following lemma.Lemma 5.
Under the assumptions of Theorem4, the following estimation holds:(84)U−Uh,P−Phh≲Υ˜+Ψ+infWh∈ℍ∩ℍh×MhUh−Whh,where Uh=Uh,Ph and Υ˜ is given by (76).Proof.
For an arbitraryWh∈ℍ∩ℍh×Mh, the inf-sup condition of B on ℍ×M leads to(85)U−Whh≲supW∈ℍ×MBU−Wh,WWh,hence(86)U−Whh≲supWh∈ℍ×MBhU−Uh,W+BhUh−Wh,WWh.
Combining the estimates (62) and (67), it comes(87)U−Whh≲Υ˜+Ψ+supW∈ℍ×MBhUh−Wh,WWh.
The continuity of the operatorBh implies that(88)U−Whh≲Υ˜+Ψ+Uh−Whh.
Thus, by the triangular inequality we deduce that(89)U−Uhh≲Υ˜+Ψ+Uh−Whh,∀Wh∈ℍ∩ℍh×Mh,or equivalently,(90)U−Uh,P−Phh≲Υ˜+Ψ+infWh∈ℍ∩ℍh×MhUh−Whh.
Thus, this lemma holds.
Combining Theorem5 and estimate (65), we have the main result in this section.Theorem 6.
Under the assumptions of Theorem4, the a posteriori error estimator Y satisfies (60).


#### 5.2.2. Proof of the Efficiency Estimate

In this section, we shall prove the estimate (61). We bound each term of the residual separately. Since by Theorem 4 the jump of U=u,ηp∈H×Xp is zero through all the edges of Ωp; hence, for all i∈8,9,10,11, we clearly have(91)Υi,K2≲JKUh,Uh=JKUh−U,Uh−U≲U−Uh,P−Phh,K.Hence it remains to estimate the local indicators fori≤7.(1)
Residual elementr1,K in Ωf. Let K∈Thf and set wK≔r1,KbK∈H01Kd, and consider(92)∫Kr1,K⋅wK=∫KfK,f+∇⋅σfuh,ph⋅wK.Introducef and use the weak formulation (19) with V=wK,0∈ℍ to obtain(93)∫Kr1,K⋅wK=∫KfK,f−f⋅wK+∫K(2μDu:DwK−pdivwK)+∫K2μdivDuh−∇ph⋅wK.Integrating by parts in this last term, we obtain(94)∫Kr1,K⋅wK=∫KfK,f−f⋅wK+2μ∫KDu−uh:DwK−∫Kp−phdivwK.Cauchy–Schwarz inequality implies that(95)∫Kr1,K⋅wK≲f−fK,fKwKK+2μ∇u−uhK+p−phK∇wKK.The inverse inequalities (55) and (56) and the obvious relation wKK≤r1,KK imply(96)hKr1,KK≲ΨK+U−Uh,P−Phh,K.(2)
Residual Elementr4,K in Ωp. Let K∈Thp and wK≔r4,KbK∈H01Kd and we consider(97)∫Kr4,K⋅wK=∫KfK,p+∇⋅σpηph,pph⋅wK.Introducef and use the weak formulation (19) with V=0,wK∈ℍ to obtain(98)∫Kr4,K⋅wK=∫KfK,p−f⋅wK+∫K2μpDηp:DwK+∫Kλp∇⋅ηp∇⋅wK−α∫Kpp∇⋅wK+α∫K∇pph⋅wK−∫K∇λp∇⋅ηph⋅wK−∫K∇⋅2μpDηph⋅wK.Integrating by parts in these three last terms, we obtain(99)∫Kr4,K⋅wK=∫KfK,p−f⋅wK+∫K2μpDηp−ηph:DwK+∫Kλp∇⋅ηp−ηph∇⋅wK−α∫Kpp−pph∇⋅wK.Cauchy–Schwarz inequality, the conditions0<λmin≤λpx≤λmax and 0<μmin≤μpx≤μmax for all x∈Ωp, and the inverse inequalities (55) and (56) lead to(100)hKr4,KK≲ΨK+U−Uh,P−Phh,K,From (96) and (100), we deduce(101)Υ1,K≲ΨK+U−Uh,P−Phh,K.(3)
Residual ElementR3 in Ωp. Let K∈Thp and, use the relation μK−1u+∇p=0 in Ωp to obtain(102)R3=μK−1uh+∇ph=−μK−1u+∇p+μK−1uh+∇ph=−μK−1u−uh+∇p−ph.As before, Cauchy–Schwarz inequality leads to(103)Υ2,K=R30,K≲U−Uh,P−Phh,K.(4)
Curl Residual ElementcurlR3 in Ωp: for K∈Thp, we use also the relation μK−1u+∇p=0 in Ωp to obtain(104)curlR3=curlμK−1uh+∇ph=−curlμK−1u+∇p+curlμK−1uh+∇ph=−curlμK−1u−uh+∇p−ph.Using Cauchy–Schwarz inequality, we obtain(105)Υ3,K=curlR30,K≲U−Uh,P−Phh,K.(5)
Residual elementR2 in Ωf: we directly see that(106)g−divuh=divu−divuh=divu−uh.Hence by Cauchy–Schwarz inequality, we conclude(107)R2K≲U−Uh,P−Phh,K.(6)
Residual elementR5 in Ωp: as before, we directly, see that(108)g−α∇⋅ηph+∇⋅uh=α∇⋅ηp−ηph+∇⋅u−uh.Hence, by Cauchy–Schwarz inequality, we have(109)R5K≲U−Uh,P−Phh,K.Inequalities (80) and (81) lead to(110)Υ4,K≲U−Uh,P−Phh,K.(7)
Normal jumpJE,nE in Ωf: for each edge/face E∈EhΩf, we consider wE=K1∪K2. As JE,nE∈ℙ0Ed, we set(111)wE≔−JE,nEbE∈H01wEd.First, the weak formulation (19) with V=wE,0∈ℍ yields(112)AU,V+BV,P=LV,that is equivalent to(113)∫wEf⋅wE=∫wE2μDu:DwE−∫wEpdivwE−∫∂wEσfu,pnE⋅wE.By elementwise partial integration, we further have(114)−∫EJE,nE⋅wE=∫wE2μDuh:DwE−∫wEphdivwE−∑i=12∫Ki−2μdivDuh+∇ph⋅wE.Hence, by the previous identity (83), we obtain(115)−∫EJE,nE⋅wE=∑i=12∫Kif+2μdivDuh+∇ph⋅wE−∫wEDu−uh:DwE+∫wEp−phdivwE.We introduce the approximationfh of f and use the Cauchy–Schwarz inequality and the inverse inequalities (57)–(58) to obtain(116)JE,nEE≲hE1/2∑i=12f−fhKi+r1,KiKi+hE−1/2u−uh1,ωE+p−phωE.The previous bound (74) of r1,Ki and the obvious estimate hE⩽hK imply that(117)hE1/2JE,nEE≲u−uh1,ωE+p−phωE+∑K′⊂ωEhK′f−fhK′.(8)
Pressure Jump inΩp: for each edge/face E∈EhΩp, we consider ωE=K1∪K2. As phE∈ℙ0E, we set(118)wE≔phEbEnE∈H01ωEd.First, we notice that asp∈H1ωE, we have by Green formula(119)∫ωE∇p⋅wE+pdivwE=0.Again by elementwise partial integration, we further have(120)∫EphEwE.nE=∑i=12∫Ki∇ph⋅wE+phdivwE.Taking the difference between these two identities, we obtain(121)∫EphEwE.nE=∑i=12∫Ti∇ph−p⋅wE+ph−pdivwE.Recalling that∇p=−μK−1u and introducing the term μK−1uh, we find(122)∫EphEwE.nE=∑i=12∫Ki∇ph+μK−1u⋅wE+ph−pdivwE=∑i=12∫Ki∇ph+μK−1uh⋅wE+ph−pdivwE+∑i=12∫KiμK−1u−uh⋅wE.Cauchy–Schwarz inequality and inverse inequalities lead to(123)phEE≲∑i=12R3KihE12+ph−pKihE−12+hE12∑i=12K−1u−uhKi.SincehE≤1, then by (103), we deduce that(124)hE12phEE≲p−phωE+K−1u−uhωE.(9)
Normal JumpGE,nE in Ωp: for each edge/face E∈EhΩp, we consider wE=K1∪K2. As GE,nE∈ℙ0Ed, we set(125)wE≔−GE,nEbE∈H01wEd.First the weak formulation (19) with V=0,wE∈H×Xp yields(126)AU,V+BV,P=LV,that is equivalent to(127)∫wEf⋅wE=∫wEμK−1u⋅wE+∫wE2μpDηp:DwE+∫wEλp∇⋅ηp∇⋅wE−α∫wEpp∇⋅wE−∫∂wEσpηp,ppnE⋅wE.By elementwise partial integration, we further have(128)−∫wEGE,nE⋅wE=∫wE2μpDηph:DwE+∫wEλp∇⋅ηph∇⋅wE−α∫wEpph∇⋅wE−∑i=12∫Ki−σpηph,pph⋅wE.By the previous identity, we obtain(129)−∫wEGE,nE⋅wE=−∫wE2μpDηp−ηph:DwE−∫wEλp∇⋅ηp−ηph∇⋅wE+α∫wEpp−pph∇⋅wE+∑i=12∫Kifh+σpηph,pph⋅wE+∫wEμK−1u−uh⋅wE+∫wEμK−1uh+∇pph⋅wE+∑i=12∫Kif−fh⋅wE.The previous bounds ofr4,K, R3, and the obvious estimate hE≤hK imply that(130)hE1/2GE,nEE≲U−Uh,P−Phh,wE+ΨK1+ΨK2.From (84), (85), and (86), we deduce the estimation(131)Υ7,K≲U−Uh,P−Phh,w˜K+∑K′⊂w˜KΨK′.(10)
Interface elements onΓfp (Υ5,K and Υ6,K): to estimate Υ5,K and Υ6,K, we fix an edge E included in Γfp and, for a constant rE fixed later on and a unit vector i, we consider(132)wE≔rEbEi,that clearly belongs toH. We take W=wE,0 and the weak formulation (19) yields(133)∫wEf⋅wE=AU,W+BW,P,that is equivalent to(134)∫wEf⋅wE=∫Kf2μDu:DwE−p∇⋅wE+∫KpμK−1u⋅wE−p∇⋅wE+∑j=1d−1∫ΓfpμαBJSKj−1uf⋅τf,jwE⋅τf,j,whereKf (resp. Kp) is the unique triangle/tetrahedron included in Ω¯f (resp. Ω¯p) having E as edge/face. On the other hand, integrating by parts in Kf and in Kp yields(135)∫Kf2μDuh:DwE−phdivwE+∫KpμK−1uh⋅wE−phdivwE+∑j=1d−1∫ΓfpμαBJSKj−1uf,h⋅τf,jwE,f⋅τf,j=−∫Kf2μdivDuh−∇ph⋅wE+∫KpμK−1uh⋅wE+∇ph⋅wE+∑j=1d−1∫ΓfpμαBJSKj−1uf,h⋅τf,jwE,f⋅τf,j−∫E(phEwE⋅nE−2μDuf,hnE⋅wE.Subtracting this identity to (89), we find(136)∫E(phEwE⋅nE−2μDuf,hnE⋅wE−∑j=1d−1∫ΓfpμαBJSKj−1ufh⋅τf,jwE,f⋅τf,j=∫Kf2μD(u−uh):D(wE)−(p−ph)divwE+∫KpμK−1u−uh⋅wE−p−phdivwE+∑j=1d−1∫ΓfpμαBJSKj−1(uf−ufh)⋅τf,jwE,f⋅τf,j−∫Kff+2μdivDuh−∇ph⋅wE−∫Kp−μK−1uh⋅wE−∇ph⋅wE.In that last terms introducing the element residualr1,K and R3, we arrive at(137)∫E(phEwE⋅nE−2μDuf,hnE⋅wE−∑j=1d−1∫ΓfpμαBJSKj−1uf−ufh⋅τf,jwE,f⋅τf,j=∫Kf2μDu−uh:DwE−p−phdivwE+∫KpμK−1u−uh⋅wE−p−phdivwE+∑j=1d−1∫ΓfpμαBJSKj−1uf−ufh⋅τf,jwE,f⋅τf,j−∫Kff−fh+rf,K⋅wE−∫KpR3⋅wE.(a)
To estimateΥ5,K, for each j=1,…,d−1, we take rE=R6j and i=τj. With this choice, the identity (90) and the inverse inequality (57) yield(138)rEE2≲∫Kf2μDu−uh:DwE−p−phdivwE+∫KpμK−1u−uh⋅wE−p−phdivwE+∑j=1d−1∫ΓfpμαBJSKj−1uf−ufh⋅τf,jwE,f⋅τf,j−∫Kff−fh+r1,K⋅wE−∫KpR3⋅wE.Hence, Cauchy–Schwarz inequality, the inverse inequalities (58), and the estimates of bounds r1,Kl and R3 lead to(139)hE12R6jE≲u−uhh,ωE+p−phh,ωE+ΨKf+ΨKp,withωE=Kf∪Kp.Thus,(140)Υ5,K≲U−Uh,P−Phh,w˜K+∑K′⊂w˜KΨK′.(b)
To estimateΥ6,K, we take rE=R7 and i=nf. As before, the identity (90), the inverse inequalities (57) and (58), and the estimates of bounds rK,l and R3 lead to(141)Υ6,K≲U−Uh,P−Phh,w˜K+∑K′⊂w˜KΨK′.Combining estimates (72), (77), (78), (79), (82), (85), (87), and (92), we have the main result of this section.Theorem 7.
Under the assumptions of Theorem4, the family ΥKK∈Th satisfies (61).


## 5.1. Optimality Result

The main result of this paper can be stated as follows:(1)
Reliability ofΥKK∈Th: the a posteriori error estimator Υ is consider reliable if it satisfies(73)U−Uh,P−Phh≲Υ+Ψ.(2)
Efficiency ofΥKK∈Th:Under the assumptions of Theorem4, the following lower error bound holds(74)ΥK≲U−Uh,P−Phh,w˜K+∑K′∈w˜KΨK′,where w˜K is a finite union of neighboring elements of K.


## 5.2. Proof of the A Posteriori Error Estimate

### 5.2.1. Proof of the Reliability Estimate

In this section, we shall prove estimate (60).Let us start with the following result.Lemma 4.
Let the assumptions of Theorem4 be satisfied. Then, for all W=V,Q∈ℍ×M, we have the estimate(75)BhU−Uh,W≲Υ˜+ΨWh,where the estimator Υ˜ is defined by(76)Υ˜≔∑K∈Th∑i=17Υi,K212.Proof.
LetW=V,Q∈ℍ×M with V=v,ξp and Q=Q1,Q2. By Theorem 3, v admits the decomposition (49) with v0,v1∈H and satisfies the properties stated in Theorem 3. Then we take Wh=Vh,0∈ℍh×Mh where Vh=vh,ξph with vh=v0,h+v1,h, where v0,h=ICl0v0 and(77)v1,h=0,in,Ωf,curlICl0ψ,in,Ωp,andξph=ICl0ξp.
In2D, ICl0ψ is the standard Clément interpolant of ψ, while in 3D, we take the vectorial Clément interpolant from [2] that satisfies the same estimate as the standard one (see [2]). Note that v0,h belongs to Hh∩H01Ωd while v1,h simply belongs to H∩Hh (ICl0ψ being in H01Ωp if d=2 and ICl0ψ∈H1Ωp3∩H0curl,Ωp if d=3), its curl belongs to H0div,Ωp, hence v1,h, its extension by zero in Ωf, stays in H0div,Ω. With these definitions and noticing that divv−vh=divv0−v0,h and that JUh,Vh=0, we may write(78)BhU−Uh,W=BhU−Uh,W−Wh=AhU−Uh,V−Vh+BhV−Vh,P−Ph+BhU−Uh,Q=AhU,V−Vh+BhV−Vh,P+BhU,Q−AhUh,V−Vh+BhV−Vh,Ph+BhUh,Q=LV−Vh+GQ−AhUh,V−Vh+BhV−Vh,Ph+BhUh,Q=f,v−vhΩ−g,Q1Ω−AhUh,V−Vh+BhV−Vh,Ph+BhUh,Q=∑K∈Thf,v−vhK−g,Q1K−AhUh,V−Vh−BhV−Vh,Ph−BhUh,Q=∑K∈Thf,v−vhK−∑K∈Thfg,Q1K−∑K∈Thf2μDuh,Dv−vhK−∑K∈ThμK−1uh,v−vhK−∑K∈Thp2μpDηph,Dξ−ξphK−∑K∈Thpλpdivhηph,divξp−ξphK−∑j=1d−1μαBJSKj−1ufh⋅τf,j,vf−vfh⋅τf,jΓfp+∑K∈Thph,divv0−v0hK+∑K∈Thpαpph,divξp−ξphK+∑K∈ThQ,divuhK+∑K∈ThpαQ2,divηphK.
Integrate by parts element by element and add boundary (resp. internal) terms that appear on the same edge (resp. at the same element) and reminding thatv=v0 and vh=v0,h in Ωf, we obtain(79)BhU−Uh,W=∑K∈ThfR1,v−vhK−∑K∈ThpR3,v−vhK+∑K∈ThpR4,ξp−ξphK−∑K∈ThfR2,Q1K−∑K∈ThpR5,Q2K+∑E∈ΕhΓfpR7,v0−v0h⋅nfE−∑E∈ΕhΓfp∑j=1d−1(R6j,v0−v0h⋅τj)E−∑E∈ΕhΩf+JE,nE,v0−v0hE−∑E∈ΕhΩp+GE,nE,ξp−ξphE+∑E∈ΕhΩppphE,v0−v0,h⋅nEE.
We now introduce the approximationfh of f for appropriated terms and we have(80)BhU−Uh,W=∑K∈Thfr1,v−vhK−∑K∈ThpR3,v0−v0,hK+∑K∈Thpr4,ξp−ξphK−∑K∈ThfR2,Q1K−∑K∈ThpR5,Q2K+∑E∈ΕhΓfpR7,v0−v0h⋅nfE−∑E∈ΕhΓfp∑j=1d−1(R6j,v0−v0h⋅τj)E−∑E∈ΕhΩf+JE,nE,v0−v0hE−∑E∈ΕhΩp+GE,nE,ξp−ξphE+∑E∈ΕhΩppphE,v0−v0,h⋅nEE+∑K∈Thf−fh,v−vhK−∑K∈ThpR3,v1−v1,hK.
Now for a triangleK∈Thp, we recall that(81)v1−v1,h=curlψ−ICl0ψinK,and use Green’s formula to get(82)∑K∈ThpR3,v1−v1,hK=∑K∈ThpR3,curlψ−ICl0ψK−∑K∈Thp[curlR3,ψ−ICl0ψK+curlR3×n,ψ−ICl0ψ∂K].
We deduce the error equation,(83)BhU−Uh,W=∑K∈Thfr1,v−vhK−∑K∈ThpR3,v0−v0,hK+∑K∈Thpr4,ξp−ξphK−∑K∈ThfR2,Q1K−∑K∈ThpR5,Q2K+∑E∈ΕhΓfpR7,v0−v0h⋅nfE−∑E∈ΕhΓfp∑j=1d−1R6j,v0−v0h⋅τjE−∑E∈Εh(Ωf+)JE,nE,v0−v0hE−∑E∈ΕhΩp+GE,nE,ξp−ξphE+∑E∈ΕhΩppphE,v0−v0,h⋅nEE+∑K∈Thf−fh,v−vhK+∑K∈Thp[curlR3,ψ−ICl0ψK−curlR3×n,ψ−ICl0ψ∂K].
Cauchy–Schwarz inequality and the approximation properties of Lemma1 imply the required estimate and finish the proof.
The second result of this section is given by the following lemma.Lemma 5.
Under the assumptions of Theorem4, the following estimation holds:(84)U−Uh,P−Phh≲Υ˜+Ψ+infWh∈ℍ∩ℍh×MhUh−Whh,where Uh=Uh,Ph and Υ˜ is given by (76).Proof.
For an arbitraryWh∈ℍ∩ℍh×Mh, the inf-sup condition of B on ℍ×M leads to(85)U−Whh≲supW∈ℍ×MBU−Wh,WWh,hence(86)U−Whh≲supWh∈ℍ×MBhU−Uh,W+BhUh−Wh,WWh.
Combining the estimates (62) and (67), it comes(87)U−Whh≲Υ˜+Ψ+supW∈ℍ×MBhUh−Wh,WWh.
The continuity of the operatorBh implies that(88)U−Whh≲Υ˜+Ψ+Uh−Whh.
Thus, by the triangular inequality we deduce that(89)U−Uhh≲Υ˜+Ψ+Uh−Whh,∀Wh∈ℍ∩ℍh×Mh,or equivalently,(90)U−Uh,P−Phh≲Υ˜+Ψ+infWh∈ℍ∩ℍh×MhUh−Whh.
Thus, this lemma holds.
Combining Theorem5 and estimate (65), we have the main result in this section.Theorem 6.
Under the assumptions of Theorem4, the a posteriori error estimator Y satisfies (60).


### 5.2.2. Proof of the Efficiency Estimate

In this section, we shall prove the estimate (61). We bound each term of the residual separately. Since by Theorem 4 the jump of U=u,ηp∈H×Xp is zero through all the edges of Ωp; hence, for all i∈8,9,10,11, we clearly have(91)Υi,K2≲JKUh,Uh=JKUh−U,Uh−U≲U−Uh,P−Phh,K.Hence it remains to estimate the local indicators fori≤7.(1)
Residual elementr1,K in Ωf. Let K∈Thf and set wK≔r1,KbK∈H01Kd, and consider(92)∫Kr1,K⋅wK=∫KfK,f+∇⋅σfuh,ph⋅wK.Introducef and use the weak formulation (19) with V=wK,0∈ℍ to obtain(93)∫Kr1,K⋅wK=∫KfK,f−f⋅wK+∫K(2μDu:DwK−pdivwK)+∫K2μdivDuh−∇ph⋅wK.Integrating by parts in this last term, we obtain(94)∫Kr1,K⋅wK=∫KfK,f−f⋅wK+2μ∫KDu−uh:DwK−∫Kp−phdivwK.Cauchy–Schwarz inequality implies that(95)∫Kr1,K⋅wK≲f−fK,fKwKK+2μ∇u−uhK+p−phK∇wKK.The inverse inequalities (55) and (56) and the obvious relation wKK≤r1,KK imply(96)hKr1,KK≲ΨK+U−Uh,P−Phh,K.(2)
Residual Elementr4,K in Ωp. Let K∈Thp and wK≔r4,KbK∈H01Kd and we consider(97)∫Kr4,K⋅wK=∫KfK,p+∇⋅σpηph,pph⋅wK.Introducef and use the weak formulation (19) with V=0,wK∈ℍ to obtain(98)∫Kr4,K⋅wK=∫KfK,p−f⋅wK+∫K2μpDηp:DwK+∫Kλp∇⋅ηp∇⋅wK−α∫Kpp∇⋅wK+α∫K∇pph⋅wK−∫K∇λp∇⋅ηph⋅wK−∫K∇⋅2μpDηph⋅wK.Integrating by parts in these three last terms, we obtain(99)∫Kr4,K⋅wK=∫KfK,p−f⋅wK+∫K2μpDηp−ηph:DwK+∫Kλp∇⋅ηp−ηph∇⋅wK−α∫Kpp−pph∇⋅wK.Cauchy–Schwarz inequality, the conditions0<λmin≤λpx≤λmax and 0<μmin≤μpx≤μmax for all x∈Ωp, and the inverse inequalities (55) and (56) lead to(100)hKr4,KK≲ΨK+U−Uh,P−Phh,K,From (96) and (100), we deduce(101)Υ1,K≲ΨK+U−Uh,P−Phh,K.(3)
Residual ElementR3 in Ωp. Let K∈Thp and, use the relation μK−1u+∇p=0 in Ωp to obtain(102)R3=μK−1uh+∇ph=−μK−1u+∇p+μK−1uh+∇ph=−μK−1u−uh+∇p−ph.As before, Cauchy–Schwarz inequality leads to(103)Υ2,K=R30,K≲U−Uh,P−Phh,K.(4)
Curl Residual ElementcurlR3 in Ωp: for K∈Thp, we use also the relation μK−1u+∇p=0 in Ωp to obtain(104)curlR3=curlμK−1uh+∇ph=−curlμK−1u+∇p+curlμK−1uh+∇ph=−curlμK−1u−uh+∇p−ph.Using Cauchy–Schwarz inequality, we obtain(105)Υ3,K=curlR30,K≲U−Uh,P−Phh,K.(5)
Residual elementR2 in Ωf: we directly see that(106)g−divuh=divu−divuh=divu−uh.Hence by Cauchy–Schwarz inequality, we conclude(107)R2K≲U−Uh,P−Phh,K.(6)
Residual elementR5 in Ωp: as before, we directly, see that(108)g−α∇⋅ηph+∇⋅uh=α∇⋅ηp−ηph+∇⋅u−uh.Hence, by Cauchy–Schwarz inequality, we have(109)R5K≲U−Uh,P−Phh,K.Inequalities (80) and (81) lead to(110)Υ4,K≲U−Uh,P−Phh,K.(7)
Normal jumpJE,nE in Ωf: for each edge/face E∈EhΩf, we consider wE=K1∪K2. As JE,nE∈ℙ0Ed, we set(111)wE≔−JE,nEbE∈H01wEd.First, the weak formulation (19) with V=wE,0∈ℍ yields(112)AU,V+BV,P=LV,that is equivalent to(113)∫wEf⋅wE=∫wE2μDu:DwE−∫wEpdivwE−∫∂wEσfu,pnE⋅wE.By elementwise partial integration, we further have(114)−∫EJE,nE⋅wE=∫wE2μDuh:DwE−∫wEphdivwE−∑i=12∫Ki−2μdivDuh+∇ph⋅wE.Hence, by the previous identity (83), we obtain(115)−∫EJE,nE⋅wE=∑i=12∫Kif+2μdivDuh+∇ph⋅wE−∫wEDu−uh:DwE+∫wEp−phdivwE.We introduce the approximationfh of f and use the Cauchy–Schwarz inequality and the inverse inequalities (57)–(58) to obtain(116)JE,nEE≲hE1/2∑i=12f−fhKi+r1,KiKi+hE−1/2u−uh1,ωE+p−phωE.The previous bound (74) of r1,Ki and the obvious estimate hE⩽hK imply that(117)hE1/2JE,nEE≲u−uh1,ωE+p−phωE+∑K′⊂ωEhK′f−fhK′.(8)
Pressure Jump inΩp: for each edge/face E∈EhΩp, we consider ωE=K1∪K2. As phE∈ℙ0E, we set(118)wE≔phEbEnE∈H01ωEd.First, we notice that asp∈H1ωE, we have by Green formula(119)∫ωE∇p⋅wE+pdivwE=0.Again by elementwise partial integration, we further have(120)∫EphEwE.nE=∑i=12∫Ki∇ph⋅wE+phdivwE.Taking the difference between these two identities, we obtain(121)∫EphEwE.nE=∑i=12∫Ti∇ph−p⋅wE+ph−pdivwE.Recalling that∇p=−μK−1u and introducing the term μK−1uh, we find(122)∫EphEwE.nE=∑i=12∫Ki∇ph+μK−1u⋅wE+ph−pdivwE=∑i=12∫Ki∇ph+μK−1uh⋅wE+ph−pdivwE+∑i=12∫KiμK−1u−uh⋅wE.Cauchy–Schwarz inequality and inverse inequalities lead to(123)phEE≲∑i=12R3KihE12+ph−pKihE−12+hE12∑i=12K−1u−uhKi.SincehE≤1, then by (103), we deduce that(124)hE12phEE≲p−phωE+K−1u−uhωE.(9)
Normal JumpGE,nE in Ωp: for each edge/face E∈EhΩp, we consider wE=K1∪K2. As GE,nE∈ℙ0Ed, we set(125)wE≔−GE,nEbE∈H01wEd.First the weak formulation (19) with V=0,wE∈H×Xp yields(126)AU,V+BV,P=LV,that is equivalent to(127)∫wEf⋅wE=∫wEμK−1u⋅wE+∫wE2μpDηp:DwE+∫wEλp∇⋅ηp∇⋅wE−α∫wEpp∇⋅wE−∫∂wEσpηp,ppnE⋅wE.By elementwise partial integration, we further have(128)−∫wEGE,nE⋅wE=∫wE2μpDηph:DwE+∫wEλp∇⋅ηph∇⋅wE−α∫wEpph∇⋅wE−∑i=12∫Ki−σpηph,pph⋅wE.By the previous identity, we obtain(129)−∫wEGE,nE⋅wE=−∫wE2μpDηp−ηph:DwE−∫wEλp∇⋅ηp−ηph∇⋅wE+α∫wEpp−pph∇⋅wE+∑i=12∫Kifh+σpηph,pph⋅wE+∫wEμK−1u−uh⋅wE+∫wEμK−1uh+∇pph⋅wE+∑i=12∫Kif−fh⋅wE.The previous bounds ofr4,K, R3, and the obvious estimate hE≤hK imply that(130)hE1/2GE,nEE≲U−Uh,P−Phh,wE+ΨK1+ΨK2.From (84), (85), and (86), we deduce the estimation(131)Υ7,K≲U−Uh,P−Phh,w˜K+∑K′⊂w˜KΨK′.(10)
Interface elements onΓfp (Υ5,K and Υ6,K): to estimate Υ5,K and Υ6,K, we fix an edge E included in Γfp and, for a constant rE fixed later on and a unit vector i, we consider(132)wE≔rEbEi,that clearly belongs toH. We take W=wE,0 and the weak formulation (19) yields(133)∫wEf⋅wE=AU,W+BW,P,that is equivalent to(134)∫wEf⋅wE=∫Kf2μDu:DwE−p∇⋅wE+∫KpμK−1u⋅wE−p∇⋅wE+∑j=1d−1∫ΓfpμαBJSKj−1uf⋅τf,jwE⋅τf,j,whereKf (resp. Kp) is the unique triangle/tetrahedron included in Ω¯f (resp. Ω¯p) having E as edge/face. On the other hand, integrating by parts in Kf and in Kp yields(135)∫Kf2μDuh:DwE−phdivwE+∫KpμK−1uh⋅wE−phdivwE+∑j=1d−1∫ΓfpμαBJSKj−1uf,h⋅τf,jwE,f⋅τf,j=−∫Kf2μdivDuh−∇ph⋅wE+∫KpμK−1uh⋅wE+∇ph⋅wE+∑j=1d−1∫ΓfpμαBJSKj−1uf,h⋅τf,jwE,f⋅τf,j−∫E(phEwE⋅nE−2μDuf,hnE⋅wE.Subtracting this identity to (89), we find(136)∫E(phEwE⋅nE−2μDuf,hnE⋅wE−∑j=1d−1∫ΓfpμαBJSKj−1ufh⋅τf,jwE,f⋅τf,j=∫Kf2μD(u−uh):D(wE)−(p−ph)divwE+∫KpμK−1u−uh⋅wE−p−phdivwE+∑j=1d−1∫ΓfpμαBJSKj−1(uf−ufh)⋅τf,jwE,f⋅τf,j−∫Kff+2μdivDuh−∇ph⋅wE−∫Kp−μK−1uh⋅wE−∇ph⋅wE.In that last terms introducing the element residualr1,K and R3, we arrive at(137)∫E(phEwE⋅nE−2μDuf,hnE⋅wE−∑j=1d−1∫ΓfpμαBJSKj−1uf−ufh⋅τf,jwE,f⋅τf,j=∫Kf2μDu−uh:DwE−p−phdivwE+∫KpμK−1u−uh⋅wE−p−phdivwE+∑j=1d−1∫ΓfpμαBJSKj−1uf−ufh⋅τf,jwE,f⋅τf,j−∫Kff−fh+rf,K⋅wE−∫KpR3⋅wE.(a)
To estimateΥ5,K, for each j=1,…,d−1, we take rE=R6j and i=τj. With this choice, the identity (90) and the inverse inequality (57) yield(138)rEE2≲∫Kf2μDu−uh:DwE−p−phdivwE+∫KpμK−1u−uh⋅wE−p−phdivwE+∑j=1d−1∫ΓfpμαBJSKj−1uf−ufh⋅τf,jwE,f⋅τf,j−∫Kff−fh+r1,K⋅wE−∫KpR3⋅wE.Hence, Cauchy–Schwarz inequality, the inverse inequalities (58), and the estimates of bounds r1,Kl and R3 lead to(139)hE12R6jE≲u−uhh,ωE+p−phh,ωE+ΨKf+ΨKp,withωE=Kf∪Kp.Thus,(140)Υ5,K≲U−Uh,P−Phh,w˜K+∑K′⊂w˜KΨK′.(b)
To estimateΥ6,K, we take rE=R7 and i=nf. As before, the identity (90), the inverse inequalities (57) and (58), and the estimates of bounds rK,l and R3 lead to(141)Υ6,K≲U−Uh,P−Phh,w˜K+∑K′⊂w˜KΨK′.Combining estimates (72), (77), (78), (79), (82), (85), (87), and (92), we have the main result of this section.Theorem 7.
Under the assumptions of Theorem4, the family ΥKK∈Th satisfies (61).


## 5.2.1. Proof of the Reliability Estimate

In this section, we shall prove estimate (60).Let us start with the following result.Lemma 4.
Let the assumptions of Theorem4 be satisfied. Then, for all W=V,Q∈ℍ×M, we have the estimate(75)BhU−Uh,W≲Υ˜+ΨWh,where the estimator Υ˜ is defined by(76)Υ˜≔∑K∈Th∑i=17Υi,K212.Proof.
LetW=V,Q∈ℍ×M with V=v,ξp and Q=Q1,Q2. By Theorem 3, v admits the decomposition (49) with v0,v1∈H and satisfies the properties stated in Theorem 3. Then we take Wh=Vh,0∈ℍh×Mh where Vh=vh,ξph with vh=v0,h+v1,h, where v0,h=ICl0v0 and(77)v1,h=0,in,Ωf,curlICl0ψ,in,Ωp,andξph=ICl0ξp.
In2D, ICl0ψ is the standard Clément interpolant of ψ, while in 3D, we take the vectorial Clément interpolant from [2] that satisfies the same estimate as the standard one (see [2]). Note that v0,h belongs to Hh∩H01Ωd while v1,h simply belongs to H∩Hh (ICl0ψ being in H01Ωp if d=2 and ICl0ψ∈H1Ωp3∩H0curl,Ωp if d=3), its curl belongs to H0div,Ωp, hence v1,h, its extension by zero in Ωf, stays in H0div,Ω. With these definitions and noticing that divv−vh=divv0−v0,h and that JUh,Vh=0, we may write(78)BhU−Uh,W=BhU−Uh,W−Wh=AhU−Uh,V−Vh+BhV−Vh,P−Ph+BhU−Uh,Q=AhU,V−Vh+BhV−Vh,P+BhU,Q−AhUh,V−Vh+BhV−Vh,Ph+BhUh,Q=LV−Vh+GQ−AhUh,V−Vh+BhV−Vh,Ph+BhUh,Q=f,v−vhΩ−g,Q1Ω−AhUh,V−Vh+BhV−Vh,Ph+BhUh,Q=∑K∈Thf,v−vhK−g,Q1K−AhUh,V−Vh−BhV−Vh,Ph−BhUh,Q=∑K∈Thf,v−vhK−∑K∈Thfg,Q1K−∑K∈Thf2μDuh,Dv−vhK−∑K∈ThμK−1uh,v−vhK−∑K∈Thp2μpDηph,Dξ−ξphK−∑K∈Thpλpdivhηph,divξp−ξphK−∑j=1d−1μαBJSKj−1ufh⋅τf,j,vf−vfh⋅τf,jΓfp+∑K∈Thph,divv0−v0hK+∑K∈Thpαpph,divξp−ξphK+∑K∈ThQ,divuhK+∑K∈ThpαQ2,divηphK.
Integrate by parts element by element and add boundary (resp. internal) terms that appear on the same edge (resp. at the same element) and reminding thatv=v0 and vh=v0,h in Ωf, we obtain(79)BhU−Uh,W=∑K∈ThfR1,v−vhK−∑K∈ThpR3,v−vhK+∑K∈ThpR4,ξp−ξphK−∑K∈ThfR2,Q1K−∑K∈ThpR5,Q2K+∑E∈ΕhΓfpR7,v0−v0h⋅nfE−∑E∈ΕhΓfp∑j=1d−1(R6j,v0−v0h⋅τj)E−∑E∈ΕhΩf+JE,nE,v0−v0hE−∑E∈ΕhΩp+GE,nE,ξp−ξphE+∑E∈ΕhΩppphE,v0−v0,h⋅nEE.
We now introduce the approximationfh of f for appropriated terms and we have(80)BhU−Uh,W=∑K∈Thfr1,v−vhK−∑K∈ThpR3,v0−v0,hK+∑K∈Thpr4,ξp−ξphK−∑K∈ThfR2,Q1K−∑K∈ThpR5,Q2K+∑E∈ΕhΓfpR7,v0−v0h⋅nfE−∑E∈ΕhΓfp∑j=1d−1(R6j,v0−v0h⋅τj)E−∑E∈ΕhΩf+JE,nE,v0−v0hE−∑E∈ΕhΩp+GE,nE,ξp−ξphE+∑E∈ΕhΩppphE,v0−v0,h⋅nEE+∑K∈Thf−fh,v−vhK−∑K∈ThpR3,v1−v1,hK.
Now for a triangleK∈Thp, we recall that(81)v1−v1,h=curlψ−ICl0ψinK,and use Green’s formula to get(82)∑K∈ThpR3,v1−v1,hK=∑K∈ThpR3,curlψ−ICl0ψK−∑K∈Thp[curlR3,ψ−ICl0ψK+curlR3×n,ψ−ICl0ψ∂K].
We deduce the error equation,(83)BhU−Uh,W=∑K∈Thfr1,v−vhK−∑K∈ThpR3,v0−v0,hK+∑K∈Thpr4,ξp−ξphK−∑K∈ThfR2,Q1K−∑K∈ThpR5,Q2K+∑E∈ΕhΓfpR7,v0−v0h⋅nfE−∑E∈ΕhΓfp∑j=1d−1R6j,v0−v0h⋅τjE−∑E∈Εh(Ωf+)JE,nE,v0−v0hE−∑E∈ΕhΩp+GE,nE,ξp−ξphE+∑E∈ΕhΩppphE,v0−v0,h⋅nEE+∑K∈Thf−fh,v−vhK+∑K∈Thp[curlR3,ψ−ICl0ψK−curlR3×n,ψ−ICl0ψ∂K].
Cauchy–Schwarz inequality and the approximation properties of Lemma1 imply the required estimate and finish the proof.
The second result of this section is given by the following lemma.Lemma 5.
Under the assumptions of Theorem4, the following estimation holds:(84)U−Uh,P−Phh≲Υ˜+Ψ+infWh∈ℍ∩ℍh×MhUh−Whh,where Uh=Uh,Ph and Υ˜ is given by (76).Proof.
For an arbitraryWh∈ℍ∩ℍh×Mh, the inf-sup condition of B on ℍ×M leads to(85)U−Whh≲supW∈ℍ×MBU−Wh,WWh,hence(86)U−Whh≲supWh∈ℍ×MBhU−Uh,W+BhUh−Wh,WWh.
Combining the estimates (62) and (67), it comes(87)U−Whh≲Υ˜+Ψ+supW∈ℍ×MBhUh−Wh,WWh.
The continuity of the operatorBh implies that(88)U−Whh≲Υ˜+Ψ+Uh−Whh.
Thus, by the triangular inequality we deduce that(89)U−Uhh≲Υ˜+Ψ+Uh−Whh,∀Wh∈ℍ∩ℍh×Mh,or equivalently,(90)U−Uh,P−Phh≲Υ˜+Ψ+infWh∈ℍ∩ℍh×MhUh−Whh.
Thus, this lemma holds.
Combining Theorem5 and estimate (65), we have the main result in this section.Theorem 6.
Under the assumptions of Theorem4, the a posteriori error estimator Y satisfies (60).


## 5.2.2. Proof of the Efficiency Estimate

In this section, we shall prove the estimate (61). We bound each term of the residual separately. Since by Theorem 4 the jump of U=u,ηp∈H×Xp is zero through all the edges of Ωp; hence, for all i∈8,9,10,11, we clearly have(91)Υi,K2≲JKUh,Uh=JKUh−U,Uh−U≲U−Uh,P−Phh,K.Hence it remains to estimate the local indicators fori≤7.(1)
Residual elementr1,K in Ωf. Let K∈Thf and set wK≔r1,KbK∈H01Kd, and consider(92)∫Kr1,K⋅wK=∫KfK,f+∇⋅σfuh,ph⋅wK.Introducef and use the weak formulation (19) with V=wK,0∈ℍ to obtain(93)∫Kr1,K⋅wK=∫KfK,f−f⋅wK+∫K(2μDu:DwK−pdivwK)+∫K2μdivDuh−∇ph⋅wK.Integrating by parts in this last term, we obtain(94)∫Kr1,K⋅wK=∫KfK,f−f⋅wK+2μ∫KDu−uh:DwK−∫Kp−phdivwK.Cauchy–Schwarz inequality implies that(95)∫Kr1,K⋅wK≲f−fK,fKwKK+2μ∇u−uhK+p−phK∇wKK.The inverse inequalities (55) and (56) and the obvious relation wKK≤r1,KK imply(96)hKr1,KK≲ΨK+U−Uh,P−Phh,K.(2)
Residual Elementr4,K in Ωp. Let K∈Thp and wK≔r4,KbK∈H01Kd and we consider(97)∫Kr4,K⋅wK=∫KfK,p+∇⋅σpηph,pph⋅wK.Introducef and use the weak formulation (19) with V=0,wK∈ℍ to obtain(98)∫Kr4,K⋅wK=∫KfK,p−f⋅wK+∫K2μpDηp:DwK+∫Kλp∇⋅ηp∇⋅wK−α∫Kpp∇⋅wK+α∫K∇pph⋅wK−∫K∇λp∇⋅ηph⋅wK−∫K∇⋅2μpDηph⋅wK.Integrating by parts in these three last terms, we obtain(99)∫Kr4,K⋅wK=∫KfK,p−f⋅wK+∫K2μpDηp−ηph:DwK+∫Kλp∇⋅ηp−ηph∇⋅wK−α∫Kpp−pph∇⋅wK.Cauchy–Schwarz inequality, the conditions0<λmin≤λpx≤λmax and 0<μmin≤μpx≤μmax for all x∈Ωp, and the inverse inequalities (55) and (56) lead to(100)hKr4,KK≲ΨK+U−Uh,P−Phh,K,From (96) and (100), we deduce(101)Υ1,K≲ΨK+U−Uh,P−Phh,K.(3)
Residual ElementR3 in Ωp. Let K∈Thp and, use the relation μK−1u+∇p=0 in Ωp to obtain(102)R3=μK−1uh+∇ph=−μK−1u+∇p+μK−1uh+∇ph=−μK−1u−uh+∇p−ph.As before, Cauchy–Schwarz inequality leads to(103)Υ2,K=R30,K≲U−Uh,P−Phh,K.(4)
Curl Residual ElementcurlR3 in Ωp: for K∈Thp, we use also the relation μK−1u+∇p=0 in Ωp to obtain(104)curlR3=curlμK−1uh+∇ph=−curlμK−1u+∇p+curlμK−1uh+∇ph=−curlμK−1u−uh+∇p−ph.Using Cauchy–Schwarz inequality, we obtain(105)Υ3,K=curlR30,K≲U−Uh,P−Phh,K.(5)
Residual elementR2 in Ωf: we directly see that(106)g−divuh=divu−divuh=divu−uh.Hence by Cauchy–Schwarz inequality, we conclude(107)R2K≲U−Uh,P−Phh,K.(6)
Residual elementR5 in Ωp: as before, we directly, see that(108)g−α∇⋅ηph+∇⋅uh=α∇⋅ηp−ηph+∇⋅u−uh.Hence, by Cauchy–Schwarz inequality, we have(109)R5K≲U−Uh,P−Phh,K.Inequalities (80) and (81) lead to(110)Υ4,K≲U−Uh,P−Phh,K.(7)
Normal jumpJE,nE in Ωf: for each edge/face E∈EhΩf, we consider wE=K1∪K2. As JE,nE∈ℙ0Ed, we set(111)wE≔−JE,nEbE∈H01wEd.First, the weak formulation (19) with V=wE,0∈ℍ yields(112)AU,V+BV,P=LV,that is equivalent to(113)∫wEf⋅wE=∫wE2μDu:DwE−∫wEpdivwE−∫∂wEσfu,pnE⋅wE.By elementwise partial integration, we further have(114)−∫EJE,nE⋅wE=∫wE2μDuh:DwE−∫wEphdivwE−∑i=12∫Ki−2μdivDuh+∇ph⋅wE.Hence, by the previous identity (83), we obtain(115)−∫EJE,nE⋅wE=∑i=12∫Kif+2μdivDuh+∇ph⋅wE−∫wEDu−uh:DwE+∫wEp−phdivwE.We introduce the approximationfh of f and use the Cauchy–Schwarz inequality and the inverse inequalities (57)–(58) to obtain(116)JE,nEE≲hE1/2∑i=12f−fhKi+r1,KiKi+hE−1/2u−uh1,ωE+p−phωE.The previous bound (74) of r1,Ki and the obvious estimate hE⩽hK imply that(117)hE1/2JE,nEE≲u−uh1,ωE+p−phωE+∑K′⊂ωEhK′f−fhK′.(8)
Pressure Jump inΩp: for each edge/face E∈EhΩp, we consider ωE=K1∪K2. As phE∈ℙ0E, we set(118)wE≔phEbEnE∈H01ωEd.First, we notice that asp∈H1ωE, we have by Green formula(119)∫ωE∇p⋅wE+pdivwE=0.Again by elementwise partial integration, we further have(120)∫EphEwE.nE=∑i=12∫Ki∇ph⋅wE+phdivwE.Taking the difference between these two identities, we obtain(121)∫EphEwE.nE=∑i=12∫Ti∇ph−p⋅wE+ph−pdivwE.Recalling that∇p=−μK−1u and introducing the term μK−1uh, we find(122)∫EphEwE.nE=∑i=12∫Ki∇ph+μK−1u⋅wE+ph−pdivwE=∑i=12∫Ki∇ph+μK−1uh⋅wE+ph−pdivwE+∑i=12∫KiμK−1u−uh⋅wE.Cauchy–Schwarz inequality and inverse inequalities lead to(123)phEE≲∑i=12R3KihE12+ph−pKihE−12+hE12∑i=12K−1u−uhKi.SincehE≤1, then by (103), we deduce that(124)hE12phEE≲p−phωE+K−1u−uhωE.(9)
Normal JumpGE,nE in Ωp: for each edge/face E∈EhΩp, we consider wE=K1∪K2. As GE,nE∈ℙ0Ed, we set(125)wE≔−GE,nEbE∈H01wEd.First the weak formulation (19) with V=0,wE∈H×Xp yields(126)AU,V+BV,P=LV,that is equivalent to(127)∫wEf⋅wE=∫wEμK−1u⋅wE+∫wE2μpDηp:DwE+∫wEλp∇⋅ηp∇⋅wE−α∫wEpp∇⋅wE−∫∂wEσpηp,ppnE⋅wE.By elementwise partial integration, we further have(128)−∫wEGE,nE⋅wE=∫wE2μpDηph:DwE+∫wEλp∇⋅ηph∇⋅wE−α∫wEpph∇⋅wE−∑i=12∫Ki−σpηph,pph⋅wE.By the previous identity, we obtain(129)−∫wEGE,nE⋅wE=−∫wE2μpDηp−ηph:DwE−∫wEλp∇⋅ηp−ηph∇⋅wE+α∫wEpp−pph∇⋅wE+∑i=12∫Kifh+σpηph,pph⋅wE+∫wEμK−1u−uh⋅wE+∫wEμK−1uh+∇pph⋅wE+∑i=12∫Kif−fh⋅wE.The previous bounds ofr4,K, R3, and the obvious estimate hE≤hK imply that(130)hE1/2GE,nEE≲U−Uh,P−Phh,wE+ΨK1+ΨK2.From (84), (85), and (86), we deduce the estimation(131)Υ7,K≲U−Uh,P−Phh,w˜K+∑K′⊂w˜KΨK′.(10)
Interface elements onΓfp (Υ5,K and Υ6,K): to estimate Υ5,K and Υ6,K, we fix an edge E included in Γfp and, for a constant rE fixed later on and a unit vector i, we consider(132)wE≔rEbEi,that clearly belongs toH. We take W=wE,0 and the weak formulation (19) yields(133)∫wEf⋅wE=AU,W+BW,P,that is equivalent to(134)∫wEf⋅wE=∫Kf2μDu:DwE−p∇⋅wE+∫KpμK−1u⋅wE−p∇⋅wE+∑j=1d−1∫ΓfpμαBJSKj−1uf⋅τf,jwE⋅τf,j,whereKf (resp. Kp) is the unique triangle/tetrahedron included in Ω¯f (resp. Ω¯p) having E as edge/face. On the other hand, integrating by parts in Kf and in Kp yields(135)∫Kf2μDuh:DwE−phdivwE+∫KpμK−1uh⋅wE−phdivwE+∑j=1d−1∫ΓfpμαBJSKj−1uf,h⋅τf,jwE,f⋅τf,j=−∫Kf2μdivDuh−∇ph⋅wE+∫KpμK−1uh⋅wE+∇ph⋅wE+∑j=1d−1∫ΓfpμαBJSKj−1uf,h⋅τf,jwE,f⋅τf,j−∫E(phEwE⋅nE−2μDuf,hnE⋅wE.Subtracting this identity to (89), we find(136)∫E(phEwE⋅nE−2μDuf,hnE⋅wE−∑j=1d−1∫ΓfpμαBJSKj−1ufh⋅τf,jwE,f⋅τf,j=∫Kf2μD(u−uh):D(wE)−(p−ph)divwE+∫KpμK−1u−uh⋅wE−p−phdivwE+∑j=1d−1∫ΓfpμαBJSKj−1(uf−ufh)⋅τf,jwE,f⋅τf,j−∫Kff+2μdivDuh−∇ph⋅wE−∫Kp−μK−1uh⋅wE−∇ph⋅wE.In that last terms introducing the element residualr1,K and R3, we arrive at(137)∫E(phEwE⋅nE−2μDuf,hnE⋅wE−∑j=1d−1∫ΓfpμαBJSKj−1uf−ufh⋅τf,jwE,f⋅τf,j=∫Kf2μDu−uh:DwE−p−phdivwE+∫KpμK−1u−uh⋅wE−p−phdivwE+∑j=1d−1∫ΓfpμαBJSKj−1uf−ufh⋅τf,jwE,f⋅τf,j−∫Kff−fh+rf,K⋅wE−∫KpR3⋅wE.(a)
To estimateΥ5,K, for each j=1,…,d−1, we take rE=R6j and i=τj. With this choice, the identity (90) and the inverse inequality (57) yield(138)rEE2≲∫Kf2μDu−uh:DwE−p−phdivwE+∫KpμK−1u−uh⋅wE−p−phdivwE+∑j=1d−1∫ΓfpμαBJSKj−1uf−ufh⋅τf,jwE,f⋅τf,j−∫Kff−fh+r1,K⋅wE−∫KpR3⋅wE.Hence, Cauchy–Schwarz inequality, the inverse inequalities (58), and the estimates of bounds r1,Kl and R3 lead to(139)hE12R6jE≲u−uhh,ωE+p−phh,ωE+ΨKf+ΨKp,withωE=Kf∪Kp.Thus,(140)Υ5,K≲U−Uh,P−Phh,w˜K+∑K′⊂w˜KΨK′.(b)
To estimateΥ6,K, we take rE=R7 and i=nf. As before, the identity (90), the inverse inequalities (57) and (58), and the estimates of bounds rK,l and R3 lead to(141)Υ6,K≲U−Uh,P−Phh,w˜K+∑K′⊂w˜KΨK′.Combining estimates (72), (77), (78), (79), (82), (85), (87), and (92), we have the main result of this section.Theorem 7.
Under the assumptions of Theorem4, the family ΥKK∈Th satisfies (61).


## 6. Summary

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 or implementation and convergence analysis of adaptive finite element methods. 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 [8]). Hence, we intend to extend our results to such anisotropic meshes.



---

*Source: 7472965-2022-03-23.xml*