# A Posteriori Error Analysis for a New Fully Mixed Isotropic Discretization of the Stationary Stokes-Darcy Coupled Problem

Authors: Houédanou Koffi Wilfrid; Adetola Jamal; Allaoui Mohamed
Journal: Abstract and Applied Analysis
 (2020)

---

## Abstract

In this paper, we develop an a posteriori error analysis for the stationary Stokes-Darcy coupled problem approximated by conforming the finite element method on isotropic meshes inℝd, d∈2,3. The approach utilizes a new robust stabilized fully mixed discretization. The a posteriori error estimate is based on a suitable evaluation on the residual of the finite element solution plus the stabilization terms. It is proven that the a posteriori error estimate provided in this paper is both reliable and efficient.


---

## Body

## 1. Introduction

There are many serious problems currently facing the world in which the coupling between groundwater and surface water is important. These include questions such as predicting how pollution discharges into streams, lakes, and rivers making its way into the water supply. This coupling is also important in technological applications involving filtration. We refer to the nice overview [1] and the references therein for its physical background, modeling, and standard numerical methods. One important issue in the modeling of the coupled Darcy-Stokes flow is the treatment of the interface condition, where the Stokes fluid meets the porous medium. In this paper, we only consider the so-called Beavers-Joseph-Saffman condition, which was experimentally derived by Beavers and Joseph in [2], modified by Saffman in [3], and later mathematically justified in [4–7].There are three popular formulations of the coupled Darcy-Stokes flow, namely, the primal formulation, the mixed formulation in the Darcy region, or the fully mixed formulation, see, for example, the works [8–16] for some mathematical analysis. Two different mixed formulations have been studied by Layton et al. in [15]. The first one enforces the weak continuity of the normal component of the velocity field on the interface, by employing a Lagrange multiplier, while the second one imposes the strong continuity in the functional space. We can call these two mixed formulations the weakly coupled formulation and the strongly coupled formulation, respectively. The weakly coupled formulation gives more freedom in the choice of the discretization in the Stokes side and the Darcy side singly. The authors in [9–11, 15, 17–19] employ the weakly coupled formulation. The works on the strongly coupled formulation have been based on the development of a unified discretization, namely, the Stokes side and the Darcy side are discretized using the same finite element. This approach simplifies the numerical implementation, only if the unified discretization is not significantly more complicated than the commonly used discretizations for the Darcy and Stokes problems. The authors in [8, 20] have proposed a conforming, unified finite element for the strongly coupled mixed formulation. Superconvergence analysis of the finite element methods for the Stokes-Darcy system was studied by Chen et al. in [21]. Other less restrictive discretizations as the nonconforming unified approach [12, 19] or the discontinuous Galerkin (DG) approach have been analyzed in [13, 14, 22]. Due to its discontinuous nature, some (DG) discretizations for the coupled Darcy-Stokes problem may break the strong coupling in the discrete level [13, 14], as they impose the normal continuity across the interface via interior penalties.The adaptive techniques have become indispensable tools and unavoidable in the field of study behavior of the error committed during solving partial differential equations (PDE). 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 equidistribute the computational effort and optimize the approximation efficiency. Since the pioneering work of Babuška and Rheinboldt [23], adaptive finite element methods based on a posteriori error estimates have been extensively investigated.A posteriori error estimations have been widely studied for both the mixed formulations of the Darcy flow [24–26] and the Stokes flow [27–36]. However, only a few works exist for the coupled Darcy-Stokes problem, see for instance [18, 37–40]. The works in [18, 38] concern the strongly coupled mixed formulation where a Hdiv conforming and nonconforming finite element methods have been employed. The papers [37, 39] concern the weakly coupled mixed formulation while [39] uses the primal formulation on the Darcy side. The authors in [40] employ a fully mixed formulation where Raviart-Thomas elements have been used to approximate the velocity in both the Stokes domain and the Darcy domain, and constant piecewise for approximating the pressure.In [16], a stabilized finite element method for the stationary-mixed Stokes-Darcy problem has been proposed for the fully mixed formulation. The authors have used the well-known MINI elements (P1b−P1) to approximate the velocity and pressure in the conduit for the Stokes equation. To capture the fully mixed technique in the porous medium region linear Lagrangian elements, P1 has been used for hydraulic (piezometric) head and Brezzi-Douglas-Marini (BDM1) piecewise constant finite elements have been used for Darcy velocity. An a priori error analysis is performed with some numerical tests confirming the convergence rates. However, to our best knowledge, they did not talk about the adaptative method for the fully mixed discretization proposed in [16]. In this case, our main objective is to perform an a posteriori error analysis by constructing reliable and efficient indicator errors. The a posteriori error estimate is based on a suitable evaluation on the residual of the finite element solution. We prove that our indicator errors are efficiency and reliability and then are optimal. The difference between our paper and the reference [40] is that our discretization uses MINI elements (P1b−P1) to approximate the velocity and pressure in the conduit for Stokes equations; P1-Lagrange elements to approximate hydraulic (piezometric) head and Brezzi-Douglas-Marini (BDM1) piecewise constant finite elements have been used for Darcy velocity. As a result, an additional term is included in the error estimator that measures the stability of the method. In order to treat appropriately this stability term, we further need a special Helmholtz decomposition ([18], Theorem 3), a regularity result ([18], Theorem 4), and an estimate of the stability error ([18], Theorem 5).The paper is organized as follows. Some preliminaries and notation are given in Section2. The efficiency result is derived using the technique of bubble function introduced by Verfürth [41] and used in a similar context by Carstensen [25, 42]. In Section 3, the a posteriori error estimates are derived. We offer our conclusion and the further works in Section 4.


## 2. Preliminaries and Notations

### 2.1. Model Problem

We consider the model of a flow in a bounded domainΩ⊂ℝdd=2or3, consisting of a porous medium domain Ωp, where the flow is a Darcy flow, and an open region Ωf, where the flow is governed by the Stokes equations. The two regions are separated by an interface Γ=∂Ωp∩∂Ωf. Let Γl=∂Ωl\Γ, l=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 Γ, we have nf=−np. Figure 1 shows a sketch of the problem domain, its boundaries, and some other notations.Figure 1
Global domainΩ consisting of the fluid region Ωf and the porous media region Ωp separated by the interface Γ.The fluid velocity and pressureufx and px are governed by the Stokes equations in Ωf:
(1)−2ν∇·Duf+∇p=ffinΩf,∇·uf=0inΩf,where T=−pI+2νDuf denotes the stress tensor, and Duf=1/2∇uf+∇ufT represents the deformation tensor. The porous media flow is governed by the following Darcy equations on Ωp through the fluid velocity upx and the piezometric head ϕx:
(2)up=−K∇ϕinΩp,∇·up=fpinΩp,We impose impermeable boundary conditions,up·np=0 on Γp, on the exterior boundary of the porous media region, and no-slip conditions, uf=0 on Γf, in the Stokes region. Both selections of boundary conditions can be modified. On Γ, the interface coupling conditions are conservation of mass, balance of forces, and a tangential condition on the fluid region’s velocity on the interface. The correct tangential condition is not completely understood (possibly due to matching a pointwise velocity in the fluid region with an averaged or homogenized velocity in the porous region). In this paper, we take the Beavers-Joseph-Saffman (-Jones), see [2–7], interfacial coupling:
(3)uf·nf+up·np=0onΓ,(4)−nf·T·nf=p−2νnf·Duf·nf=ρgϕonΓ,(5)−nf·T·τj=−2nf·Duf·τj=ατj·Kτjuf·τj,1≤j≤d−1onΓ.This is a simplification of the original and more physically realistic Beavers-Joseph conditions (in whichuf·τj in (2.8) is replaced by uf−up·τj; see [2]). Here, we denote ff and fp are the body forces in the fluid region and source in the porous region, K is the symmetric positive define (SPD) hydraulic conductivity tensor, and α is the constant parameter.We shall also assume that all material and fluid parameters defined above are uniformly positive and bounded, i.e.,(6)0≤kmin≤λK≤kmax<∞.


### 2.2. Notations and the Weak Formulation

In this part, 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 L2WL2WN or L2Wd×d inner product. For shortness, if W is equal to Ω, we will drop the index Ω, while for any m≥0, ·m,l=·m,Ωl, .m,l=.m,Ωl and .,.l=·,·Ωl, for l=f,s. 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
(7)vm,Ω≔∑i=0Nvim,Ω21/2,vm,Ω≔∑i=0Nvim,Ω21/2.For a connected open subset of the boundaryE⊂∂Ωs∪∂Ωd, we write .,.E for the L2E inner product (or duality pairing), that is, for scalar-valued functions λ, σ, one defines:
(8)λ,σE≔∫Eλσds.By setting the space(9)Hdiv≔Hdiv;Ωp=vp∈L2Ωpd:∇·vp∈L2Ωp,we introduce the following spaces:
(10)Xf≔vf∈H1Ωfd:vf=0onΓf,Qf≔L2Ωf,Xp≔vp∈Hdiv;Ωp:vp·np=0onΓp,Qp≔L2Ωp.For the spacesXf and Xp, we define the following norms:
(11)vf1≔vfΩf2+vf1,Ωf2,withvf1,Ωf=∇vfΩf,∀vf∈Xf,vpdiv≔vpΩp2+∇·vpΩp2,∀vp∈Xp.The variational formulation of the steady-state Stokes-Darcy problem (1)–(5) reads as find uf,p;up,ϕ∈Xf,Qf;Xp,Qp satisfying
(12)afuf,vf−bfvf,p+cΓvf,ϕ=ff,vfΩf,∀vf∈Xf,(13)bfuf,q=0,∀q∈Qf,(14)apup,vp−bpvp,ϕ−cΓvp,ϕ=0,∀vp∈Xp,(15)bpup,ψ=ρgfp,ψΩp,ψ∈Qp,where the bilinear forms are defined as
(16)afuf,vf≔2νDuf,DvfΩf+∑j=1d−1ατj·Kτj<uf·τj,vf·τj>Γ,apup,vp≔ρgK−1up,vpΩp,bfvf,p≔p,∇·vfΩf,bpvp,ϕ≔ρgϕ,∇·vpΩp,cΓvf,ϕ≔ρg<ϕ,vf·nf>Γ.After introducing, forU=uf,p,up,ϕ∈Xf×Qp×Xp×Qp=H and V=vf,q,vp,ψ∈Xf×Qp×Xp×Qp,
(17)LU,V≔afuf,vf−bfvf,p+bfuf,q+apup,vp−bpvp,ϕ+bpup,ψ+cΓvf−vp,ϕ,FV≔ff,vfΩf+ρgfp,ψΩp,the weak formulation (12)–(15) can be equivalently rewritten as follows: find U∈H satisfying
(18)LU,V=FV,∀V∈H.It is easy to verify that this variational formulation is well-posedness [16].To end this section, we recall the following Poincaré, Korn’s, and the trace inequalities, which will be used in the later analysis; there exist constantCp, CK, and Cv, only depending on ωf such that for all vf∈Xf,
(19)vf≤Cpvf1,vf≤CKDvfΩf,vfΓ≤CvvfΩf1/2vf1,Ωf1/2.Besides, there exists a constantC~v that only depends on Ωp such that for all ψ∈Qp,
(20)ψL2Γ≤C~vψΩp1/2ψ1,Ωp1/2.


### 2.3. Fully Mixed Isotropic Discretization

First, we consider the family of triangulationsTh of Ω, consisting of Thf and Thp, which are regular triangulations of Ωf and Ωp, respectively, where h>0 is a positive parameter. We also assume that on the interface Γ the two meshes of Thf and Thp, which form the regular triangulation Th≔Thf∪Thp, coincide.The domain of the uniformly regular triangulationΩ¯f∪Ω¯p is such that Ω¯=∪K:K∈Th and h=maxK∈ThhK. There exist positive constants c1 and c2 satisfying c1h≤hK≤c2ρK. To approximate the diameter hK of the triangle (or tetrahedral) K, ρK is the diameter of the greatest ball included in K. Based on the subdivisions Thf and Thp, we can define finite element spaces Xfh⊂Xh, Qfh⊂Qf, Xph⊂Xp, and Qph⊂Qp. We consider the well-known MINI elements P1b−P1 to approximate the velocity and the pressure in the conduit for Stokes equations [44]. To capture the fully mixed technique in the porous medium region linear Lagrangian elements, P1 are used for hydraulic (piezometric) head and Brezzi-Douglas-Marini (BDM1) piecewise constant finite elements are used for Darcy velocity [45].In the fluid region, we select for the Stokes problem the finite element spacesXfh,Qfh that satisfy the velocity-pressure inf-sup condition: there exists a constant Cf>0, independent of h, such that,
(21)inf0≠qh∈Qfhsup0≠vfh∈Xfhbfvfh,qhvfh1,ΩfqfhΩf≥Cf.In the porous region, we use the finite element spacesXph,Qph that also satisfy a standard inf-sup condition: there exist a constant Cp>0 such that for all ϕh∈Qph,
(22)inf0≠ϕh∈Qphsup0≠vph∈Xphbpvph,ϕhvfhdivϕhΩp≥Cp.Then, the finite element discretization of (18) is to find Uh∈Hh=Xfh×Qfh×Xph×Qph such that
(23)LUh,Vh+JΓUh,Vh=FVh,∀Vh∈Hh.This is the natural discretization of the weak formulation (18) except that the stabilized term JΓUh,Vh is added. This bilinear form JΓ.,. is defined by
(24)JΓUh,Vh≔δhufh−uph·nf,vfh−vph·nfΓ,0<h<1.We are now able to define the norm onHh:
(25)Vh≔vfh1,Ωf2+qfhΩf2+vphdiv2+ψhΩp2+h−1vfh−vphΓ2.We have the following results (see [16], Theorem 2 and Theorem 3):Theorem 1.
There exists a unique solutionUh∈Hh to problem (23) and if the solution U∈H of the continuous problem (18) is smooth enough, then we have
(26)U−Uhh≤CUh.
Below, in order to avoid excessive use of constants, the abbreviationx≲y stand for x≤cy, with c a positive constant independent of x, y, and Th.Remark 2. (Galerkin orthogonality relation).
LetU=uf,p,up,ϕ∈H be the exact solution and Uh=ufh,ph,uph,ϕh∈Hh be the finite element solution. Then, for any Vh=vfh,ph,vph,ψh∈Hh, and using technical regularity result Theorem 4 below, we can subtract (18) to (23) to obtain the Galerkin orthogonality relation:
(27)LhU−Uh,Vh=LhU,Vh−LhUh,Vh=LU,Vh−LhUh,Vh=FVh−FVh=0.
Thus, we have the relation:(28)2νDef,DvfhΩf+∑j=1d−1ατj·K·τjef·τj,vfh·τjΓ−εp,∇·vfhΩf+qh,∇·efΩf+ρgK−1ep,vphΩp−λϕ,∇·vphΩp+ψh,∇·epΩp+λϕ,vΓ=0,where here and below, the errors in the velocity and in the pressure of Stokes equations and errors in the hydraulic and Darcy velocity equations are, respectively, defined by
(29)ef≔uf−ufh,εp≔p−ph,ep≔up−uph,
λ
ϕ
=
ϕ
−
ϕ
h
.


## 2.1. Model Problem

We consider the model of a flow in a bounded domainΩ⊂ℝdd=2or3, consisting of a porous medium domain Ωp, where the flow is a Darcy flow, and an open region Ωf, where the flow is governed by the Stokes equations. The two regions are separated by an interface Γ=∂Ωp∩∂Ωf. Let Γl=∂Ωl\Γ, l=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 Γ, we have nf=−np. Figure 1 shows a sketch of the problem domain, its boundaries, and some other notations.Figure 1
Global domainΩ consisting of the fluid region Ωf and the porous media region Ωp separated by the interface Γ.The fluid velocity and pressureufx and px are governed by the Stokes equations in Ωf:
(1)−2ν∇·Duf+∇p=ffinΩf,∇·uf=0inΩf,where T=−pI+2νDuf denotes the stress tensor, and Duf=1/2∇uf+∇ufT represents the deformation tensor. The porous media flow is governed by the following Darcy equations on Ωp through the fluid velocity upx and the piezometric head ϕx:
(2)up=−K∇ϕinΩp,∇·up=fpinΩp,We impose impermeable boundary conditions,up·np=0 on Γp, on the exterior boundary of the porous media region, and no-slip conditions, uf=0 on Γf, in the Stokes region. Both selections of boundary conditions can be modified. On Γ, the interface coupling conditions are conservation of mass, balance of forces, and a tangential condition on the fluid region’s velocity on the interface. The correct tangential condition is not completely understood (possibly due to matching a pointwise velocity in the fluid region with an averaged or homogenized velocity in the porous region). In this paper, we take the Beavers-Joseph-Saffman (-Jones), see [2–7], interfacial coupling:
(3)uf·nf+up·np=0onΓ,(4)−nf·T·nf=p−2νnf·Duf·nf=ρgϕonΓ,(5)−nf·T·τj=−2nf·Duf·τj=ατj·Kτjuf·τj,1≤j≤d−1onΓ.This is a simplification of the original and more physically realistic Beavers-Joseph conditions (in whichuf·τj in (2.8) is replaced by uf−up·τj; see [2]). Here, we denote ff and fp are the body forces in the fluid region and source in the porous region, K is the symmetric positive define (SPD) hydraulic conductivity tensor, and α is the constant parameter.We shall also assume that all material and fluid parameters defined above are uniformly positive and bounded, i.e.,(6)0≤kmin≤λK≤kmax<∞.


## 2.2. Notations and the Weak Formulation

In this part, 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 L2WL2WN or L2Wd×d inner product. For shortness, if W is equal to Ω, we will drop the index Ω, while for any m≥0, ·m,l=·m,Ωl, .m,l=.m,Ωl and .,.l=·,·Ωl, for l=f,s. 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
(7)vm,Ω≔∑i=0Nvim,Ω21/2,vm,Ω≔∑i=0Nvim,Ω21/2.For a connected open subset of the boundaryE⊂∂Ωs∪∂Ωd, we write .,.E for the L2E inner product (or duality pairing), that is, for scalar-valued functions λ, σ, one defines:
(8)λ,σE≔∫Eλσds.By setting the space(9)Hdiv≔Hdiv;Ωp=vp∈L2Ωpd:∇·vp∈L2Ωp,we introduce the following spaces:
(10)Xf≔vf∈H1Ωfd:vf=0onΓf,Qf≔L2Ωf,Xp≔vp∈Hdiv;Ωp:vp·np=0onΓp,Qp≔L2Ωp.For the spacesXf and Xp, we define the following norms:
(11)vf1≔vfΩf2+vf1,Ωf2,withvf1,Ωf=∇vfΩf,∀vf∈Xf,vpdiv≔vpΩp2+∇·vpΩp2,∀vp∈Xp.The variational formulation of the steady-state Stokes-Darcy problem (1)–(5) reads as find uf,p;up,ϕ∈Xf,Qf;Xp,Qp satisfying
(12)afuf,vf−bfvf,p+cΓvf,ϕ=ff,vfΩf,∀vf∈Xf,(13)bfuf,q=0,∀q∈Qf,(14)apup,vp−bpvp,ϕ−cΓvp,ϕ=0,∀vp∈Xp,(15)bpup,ψ=ρgfp,ψΩp,ψ∈Qp,where the bilinear forms are defined as
(16)afuf,vf≔2νDuf,DvfΩf+∑j=1d−1ατj·Kτj<uf·τj,vf·τj>Γ,apup,vp≔ρgK−1up,vpΩp,bfvf,p≔p,∇·vfΩf,bpvp,ϕ≔ρgϕ,∇·vpΩp,cΓvf,ϕ≔ρg<ϕ,vf·nf>Γ.After introducing, forU=uf,p,up,ϕ∈Xf×Qp×Xp×Qp=H and V=vf,q,vp,ψ∈Xf×Qp×Xp×Qp,
(17)LU,V≔afuf,vf−bfvf,p+bfuf,q+apup,vp−bpvp,ϕ+bpup,ψ+cΓvf−vp,ϕ,FV≔ff,vfΩf+ρgfp,ψΩp,the weak formulation (12)–(15) can be equivalently rewritten as follows: find U∈H satisfying
(18)LU,V=FV,∀V∈H.It is easy to verify that this variational formulation is well-posedness [16].To end this section, we recall the following Poincaré, Korn’s, and the trace inequalities, which will be used in the later analysis; there exist constantCp, CK, and Cv, only depending on ωf such that for all vf∈Xf,
(19)vf≤Cpvf1,vf≤CKDvfΩf,vfΓ≤CvvfΩf1/2vf1,Ωf1/2.Besides, there exists a constantC~v that only depends on Ωp such that for all ψ∈Qp,
(20)ψL2Γ≤C~vψΩp1/2ψ1,Ωp1/2.


## 2.3. Fully Mixed Isotropic Discretization

First, we consider the family of triangulationsTh of Ω, consisting of Thf and Thp, which are regular triangulations of Ωf and Ωp, respectively, where h>0 is a positive parameter. We also assume that on the interface Γ the two meshes of Thf and Thp, which form the regular triangulation Th≔Thf∪Thp, coincide.The domain of the uniformly regular triangulationΩ¯f∪Ω¯p is such that Ω¯=∪K:K∈Th and h=maxK∈ThhK. There exist positive constants c1 and c2 satisfying c1h≤hK≤c2ρK. To approximate the diameter hK of the triangle (or tetrahedral) K, ρK is the diameter of the greatest ball included in K. Based on the subdivisions Thf and Thp, we can define finite element spaces Xfh⊂Xh, Qfh⊂Qf, Xph⊂Xp, and Qph⊂Qp. We consider the well-known MINI elements P1b−P1 to approximate the velocity and the pressure in the conduit for Stokes equations [44]. To capture the fully mixed technique in the porous medium region linear Lagrangian elements, P1 are used for hydraulic (piezometric) head and Brezzi-Douglas-Marini (BDM1) piecewise constant finite elements are used for Darcy velocity [45].In the fluid region, we select for the Stokes problem the finite element spacesXfh,Qfh that satisfy the velocity-pressure inf-sup condition: there exists a constant Cf>0, independent of h, such that,
(21)inf0≠qh∈Qfhsup0≠vfh∈Xfhbfvfh,qhvfh1,ΩfqfhΩf≥Cf.In the porous region, we use the finite element spacesXph,Qph that also satisfy a standard inf-sup condition: there exist a constant Cp>0 such that for all ϕh∈Qph,
(22)inf0≠ϕh∈Qphsup0≠vph∈Xphbpvph,ϕhvfhdivϕhΩp≥Cp.Then, the finite element discretization of (18) is to find Uh∈Hh=Xfh×Qfh×Xph×Qph such that
(23)LUh,Vh+JΓUh,Vh=FVh,∀Vh∈Hh.This is the natural discretization of the weak formulation (18) except that the stabilized term JΓUh,Vh is added. This bilinear form JΓ.,. is defined by
(24)JΓUh,Vh≔δhufh−uph·nf,vfh−vph·nfΓ,0<h<1.We are now able to define the norm onHh:
(25)Vh≔vfh1,Ωf2+qfhΩf2+vphdiv2+ψhΩp2+h−1vfh−vphΓ2.We have the following results (see [16], Theorem 2 and Theorem 3):Theorem 1.
There exists a unique solutionUh∈Hh to problem (23) and if the solution U∈H of the continuous problem (18) is smooth enough, then we have
(26)U−Uhh≤CUh.
Below, in order to avoid excessive use of constants, the abbreviationx≲y stand for x≤cy, with c a positive constant independent of x, y, and Th.Remark 2. (Galerkin orthogonality relation).
LetU=uf,p,up,ϕ∈H be the exact solution and Uh=ufh,ph,uph,ϕh∈Hh be the finite element solution. Then, for any Vh=vfh,ph,vph,ψh∈Hh, and using technical regularity result Theorem 4 below, we can subtract (18) to (23) to obtain the Galerkin orthogonality relation:
(27)LhU−Uh,Vh=LhU,Vh−LhUh,Vh=LU,Vh−LhUh,Vh=FVh−FVh=0.
Thus, we have the relation:(28)2νDef,DvfhΩf+∑j=1d−1ατj·K·τjef·τj,vfh·τjΓ−εp,∇·vfhΩf+qh,∇·efΩf+ρgK−1ep,vphΩp−λϕ,∇·vphΩp+ψh,∇·epΩp+λϕ,vΓ=0,where here and below, the errors in the velocity and in the pressure of Stokes equations and errors in the hydraulic and Darcy velocity equations are, respectively, defined by
(29)ef≔uf−ufh,εp≔p−ph,ep≔up−uph,
λ
ϕ
=
ϕ
−
ϕ
h
.


## 3. A Posteriori Error Analysis

### 3.1. Some Technical Results

Our a posteriori analysis requires some analytical results that are recalled. We define the space(30)H=v∈Hdiv,Ω:v∣Ωf∈Xfandv∣Ωp∈Xpwith the norm
(31)vH≔vf1,Ωf2+vpΩp2+∇·vpΩp.The first one concerns a sort of Helmholtz decomposition of elements ofH. Recall first that if d=3,
(32)H0curl,Ωp=ψ∈L2Ωp3:curlψ∈L2Ωp3andψ×n=0on∂Ωp.Theorem 3. (see Ref. [18], page 708).
Anyv∈H admits the Helmholtz-type decomposition
(33)v=v0+v1,where v0,v1∈H but satisfying v0∈H1Ωd,
(34)v1=0inΩf,curlβpinΩp,where βp∈H01Ωp if d=2, while βp∈H1Ωp3∩H0curl,Ωp if d=3, with the estimate
(35)v01,Ω+βp1,Ωp≲vH.The second result that we need is a regularity result for the solutionU=uf,p,up,ϕ∈H of (18) is the following theorem:Theorem 4. (see [18], page 710).
LetU∈H be the unique solution of (18). If K∈C0,1Ω¯pd×d, then there exists ε>0 such that
(36)u∣Ωp∈H1/2+εΩpd.Let us finish this section by an estimation of the stability error (see [18], Theorem 5):Theorem 5.
For anyUh=ufh,ph,uph,ϕh∈Hh, we have
(37)infWh∈Hh∩HUh−Whh2≲JΓUh,Uh.


### 3.2. Error Estimator

In order to solve the Stokes-Darcy coupled 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 and then the lower and upper error bounds are derived in Section3.3.


#### 3.2.1. Error Equations

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. Thus, we define the error equations: letU=uf,p,up,ϕ∈H be the exact solution and Uh=ufh,ph,uph,ϕh∈Hh be the finite element solution. Then, for any Vh=vfh,qh,vph,ψh∈Hh and V=vf,p,vp,ψ∈H, using the Helmholtz decomposition (Theorem 3), we have
(38)LhU−Uh,V=LhU−Uh,V−Vh(39)=∑K∈ThfRKfUh,V−VhK+∑K∈ThpRKpUh,V−VhK,where
(40)RKfUh,V−VhK=ff+2ν∇·Dufh−∇ph,vf−vfhK−q−qh,∇·ufhK−∑E∈Eh∂K∩Γ∑j=1d−12νnf·Dufh·τj+ατj·K·τjufh·τj,vf−vfh·τjE+∑E∈Eh∂K∩Γph−2νnf·Dufh·nf−ρgϕh,vf−vfh·nfE,(41)RKpUh,V−VhK=curlρgK−1uph+∇ϕh,βp−βphK+ρgfp−∇·uph,ψ−ψhK−∑E∈Eh∂K∩ΩpρgK−1uph+∇ϕh×nE,ψ−ψhE+∑E∈Eh∂K∩ΩpρgϕhnEE,βp−βphE−∑E∈Eh∂K∩ΓρgϕhnEE,v−vhEE.


#### 3.2.2. Residual Error Estimators

Definition 6 (a posteriori error indicators).
The residual error estimator is locally defined by(42)ΘK=ΘK,f2+ΘK,p21/2foreachK∈Th,where
(43)ΘK,f2=hK2ffh+2ν∇·Dufh−∇phK2+∇·ufhK2+∑E∈Eh∂K∩ΓhE∑j=1d−12νnf·Dufh·τj+ατj·K·τjufh·τjE2+∑E∈Eh∂K∩ΓhEph−2νnf·Dufh·nf−ρgϕhE2,(44)ΘK,p2=hK2curlρgK−1uph+∇ϕhK2+ρgfp−∇·uphK2+∑E∈Eh∂K∩ΩpρgK−1uph+∇ϕh×npEE2+∑E∈Eh∂K∩Ω¯phEρgϕhnpEE2+∑E∈Eh∂K∩ΓδhEhufh−uph·nfEE2.
The global residual error estimator is given by(45)Θ≔∑K∈ThΘK21/2.
Furthermore, denote the local and global approximation terms by(46)ζK≔hKff−ffhK,∀K∈Thf,ρgfp−fphK,∀K∈Thp,ζ≔∑K∈ThζK21/2,where the global function ffh:Ωf→ℝd is defined by
(47)ffh∣K=fK=1K∫Kffxdx,∀K∈Thf,while in Ωp, we take fh∣K=fK for all K∈Thp, as the unique element of ℙ1K such that
(48)∫KfKxqxdx=∫Kfxqxdx∀q∈ℙ1K.Remark 7.
The residual character of each term on the right-hand sides of (43) and (44) is quite clear since if ufh,ph,uph,ϕh∈Hh would be the exact solution of (18), then they would vanish.


#### 3.2.3. Analytical Tools

(1) Inverse Inequalities. In order to derive the lower error bounds, we proceed similarly as in [25, 42] (see also [46]), 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 [41]). 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∂KE, and 0≤bE≤1inωE. We also recall from [47] 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 [47]).Lemma 8.
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
(49)qK≲qbK1/2K≲qK,∀q∈ℙkK,(50)qbK1,K≲hK−1qK,∀q∈ℙkK,(51)pE≲bE1/2pE≲pE,∀p∈ℙkE,(52)LpK+hELp1,K≲hE1/2pE,∀p∈ℙkE.(2) Continuous Trace Inequality.Lemma 9 (continuous trace inequality).
There exists a positive constantβ1>0 depending only on σ0 such that
(53)v∂K2≤β1vKv1,K,∀K∈Th,∀v∈H1Kd.(3) Clément Interpolation Operator. In order to derive the upper error bounds, we introduce the Clément interpolation operator ICl0:H01Ω→PcbTh that approximates optimally nonsmooth functions by continuous piecewise linear functions:
(54)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 [48], Section 3.Lemma 10.
There exist constantsC1,C2>0, independent of h, such that for all v∈H01Ω, there hold
(55)v−ICl0vK≤C1hKv1,ΔK,∀K∈Th,v−ICl0vE≤C2hE12v1,ΔE,∀E∈Eh,where ΔK≔∪K′∈Th:K′∩K≠∅ and ΔE≔∪K′∈Th:K′∩E≠∅.


### 3.3. Optimality ofΘKK∈Th

#### 3.3.1. Reliability Result

Theorem 11 (reliability ofΘ).
LetU=uf,p,up,ϕ∈H be the exact solution of (18) and Uh=ufh,ph,uph,ϕh∈Hh be the finite element solution of (23). There exist a constant Crel>0 such that the following estimate holds:
(56)U−Uhh≤CrelΘ+ζ.Proof.
We takeψh=0=qh in error equation (38)–(41) and we obtain
(57)LhU−Uh,V=LhU−Uh,V−Vh=∑K∈ThfRKfUh,V−VhK+∑K∈ThpRKpUh,V−VhK,where
(58)RKfUh,V−VhK=ffh+2ν∇·Dufh−∇ph,vf−vfhK−q,∇·ufhK+ff−ffh,vf−vfhK−∑E∈Eh∂K∩Γ∑j=1d−12νnf·Dufh·τj+ατj·K·τjufh·τj,vf−vfh·τjE+∑E∈Eh∂K∩Γph−2νnf·Dufh·nf−ρgϕh,vf−vfh·nfE,(59)RKpUh,V−VhK=curlρgK−1uph+∇ϕh,βp−βphK+ρgfph−∇·uph,ψK+ρgfp−fph,ψK−∑E∈Eh∂K∩ΩpρgK−1uph+∇ϕh×nE,ψE+∑E∈Eh∂K∩ΩpρgϕhnEE,βp−βphE−∑E∈Eh∂K∩ΓρgϕhnEE,v−vhEE.
The inf-sup condition ofLh leads to
(60)∥U−Uh∥h≤CsupV∈H\0∣LhU−Uh,V∣∥V∥h.
Now, using the error equation (57)–(59), Cauchy-Schwarz inequality, and the Clément operator of Lemma 10, we deduce the estimate (56). The proof is complete.


#### 3.3.2. Efficiency Result

To prove local efficiency forw⊂Ω, let us denote by
(61)vf,vph,w2≔∑K⊂w¯∩Ω¯fvf1,K2+∑K⊂w¯∩Ω¯pvpK2+divhvpK2+vf×nΓ∩w¯2+∑K⊂Ω¯JKvf,vf,where
(62)JKvf,vf≔∑E∈EhΩ¯f∩EhKδhE−1vfEE2;ε,λw≔εw+λw.The main result of this subsection can be stated as follows.Theorem 12 (efficiency ofΘ).
Under the assumptions of Theorem4, the following lower error bound holds:
(63)ΘK≲ef,eph,w~K+εp,λϕw~K+∑K′⊂w~KζK′,where w~K is a finite union of neighboring elements of K.Proof.
We begin by bounding each term of the residuals separately.(i)
Element residual inΩf: to estimate hK2ffh+2ν∇·Dufh−∇phK2, we choose in error equation (57)–(59) for each K∈Thf, V=vsK,0,vpK,0 and Vh=0,0,0,0 with vpK=0 on Ωp,(64)vfK=ffh+2ν∇·Dufh−∇phbKonK∈Thf,0onΩf\K,for obtained, LhU−Uh,V=ffh+2ν∇·Dufh−∇phbK1/2K2+ff−ffh,vfKK. Noted that ffh+2ν∇·Dufh−∇phbK∈ℙkK and vanish on ∂K. Because LhU−Uh,V=Def,vfKK in this case, then we have Def,vfKK=ffh+2ν∇·Dufh−∇phbK1/2K2+ff−ffh,vfKK. The first inverse inequality (49) and the Cauchy-Schwarz inequality lead to
(65)ffh+2ν∇·Dufh−∇phK2~ffh+2ν∇·Dufh−∇phbK1/2K2=∫KDefffh+2ν∇·Dufh−∇phbK−ff−ffh,vfKK≤efKffh+2ν∇·Dufh−∇phbK1,K+ff−ffhKvfKK.
Using inverse inequality (50), we deduce the estimate:
(66)hKffh+2ν∇·Dufh−∇phK≲ef,eph,w~K+ζK.(ii)
Divergence element residual inΩf (estimation of ∇·ufhK2): for each K∈Thf, we have,(67)∇ufhK=∇uf−ufhK≲uf−ufh1,K.(iii)
Element residual inΩp: we have for each K∈Thp,(68)ρgfph−∇·uphK=ρgfph−∇·uph−ρgfp−∇·upK=ρgdivhup−uph−ρgfp−fphK≲ef,eph,w~K+ζK.(iv)
Curl element residual inΩp: for K∈Thp, we set CK=curlρgK−1uph+∇ϕh and WK=CKbK. Hence, we notice that curlWK belongs to H and is divergence free; therefore, by Equations (57)–(59), we obtain with Vh=0 and βp=WK, LhU−Uh,WK=RKpUh,WKK=∥curlρgK−1uph+∇ϕhbK1/2∥K2. The first inverse inequality (49) and the Cauchy-Schwarz inequality lead to(69)curlρgK−1uph+∇ϕhK2~curlρgK−1uph+∇ϕhbK1/2K2=∫KρgK−1up−uph+∇ϕ−ϕh·curlWK≤ρgK−1up−uphK+∇ϕ−ϕhK·curlWKK≲ρgK−1epK+∇λϕK·WKK.
Again, the inverse inequality (49) allows to get
(70)curlρgK−1uph+∇ϕhK≲ef,eph,w~K+εp,λϕw~K.(v)
Interface elements onΓ: we fix an edge E included in Γ and for a constant rE fixed later on a unit vector N, we consider(71)WE=rEbEN,that clearly belongs to H. Hence, by residual equation (57)–(59), we obtain with Vh=0,
(72)LhU−Uh,WE=RKffUh,WEKf+RKppUh,WEKp,
whereKf (resp., Kp) is the unique triangle/tetrahedron included in Ω¯f (resp., Ω¯p) having E as edge/face, and
(73)RKffUh,WEKf=ff+2ν∇·Dufh−∇ph,WEKf−q,∇·ufhKf−∑E∈Eh∂Kf∩Γ∑j=1d−12νnf·Dufh·τj+ατj·K·τjufh·τj,WE·τjE+∑E∈Eh∂Kf∩Γph−2νnf·Dufh·nf−ρgϕh,WE·nfE,RKpUh,WEKp=curlρgK−1uph+∇ϕh,βpKp+ρgfp−∇·uph,ψKp−∑E∈Eh∂Kp∩ΩpρgK−1uph+∇ϕh×nE,ψE+∑E∈Eh∂K∩ΩpρgϕhnEE,βpE−∑E∈Eh∂K∩ΓρgϕhnEE,WEE.
TakenWE=0 in Kp, q=0 in Ω and for each j=1,⋯,d−1, rE=2νnf·Dufh·τj+α/√τj·K·τjufh·τj with N=τj. We have
(74)LhU−Uh,WE=ff+2ν∇·Dufh−∇ph,WEKf−2νnf·Dufh·τj+ατj·K·τjufh·τj,WE·τjE=ff+2ν∇·Dufh−∇ph,WEKf−∥rEbE1/2∥E2.
Hence,(75)rEE2~rEbE1/2E2=ff+2ν∇·Dufh−∇ph,WEKf−LhU−Uh,WE=ffh+2ν∇·Dufh−∇ph,WEKf−∫Kf2νDef:DWE+∫KfεpdivWE+ff−ffh,WEKf−∑j=1d−1ατj·K·τjef·τj,WE·τjE.
Inverse inequalities (50) and (51) and Cauchy-Schwarz inequality lead to
(76)hE1/22νnf·Dufh·τj+ατj·K·τjufh·τjE≲ef,eph,wE+εp,λϕwE+∑K⊂wEζK,with wE=Kf∪Kp.
TakenWE=0 in Kf, q=0 in Ω and for each j=1,⋯,d−1, rE=2νnf·Dufh·τj+α/√τj·K·τjufh·τj with N=nf. As before the identities (72)–(74) and the inverse inequalities (50) and (51) lead to
(77)hE1/2∥ph−2νnf·Dufh·nf−ρgϕh∥E≲∥ef,ep∥h,wE+∥εp,λϕ∥wE+∑K⊂wEζK.(vi)
Piezometric head jump inΩp: for each edge/face E∈EhΩp, we consider wE=T1∪T2. As ρgϕhnpE∈ℙ1Ed, we set(78)WE≔ρgϕhnpEbE∈H01wEd.
Using the residual equation (57)–(59), we obtain with Vh=0 and V=0,0,βp,0 where βp=WE:
(79)LhU−Uh,V=ρgϕhnEE,WEE.
Therefore, Cauchy-Schwarz inequality and inverse inequalities (51) and (52) lead to
(80)ρgϕhnEEE2~ρgϕhnEEbE1/2E2=∫EρgϕhnEE·WE=∑i=12ρgϕ−ϕh,∇·WEKi≤∑i=12ρgλϕKi∇·WEKi≲∑i=12λϕKihE−1ρgϕhnEEE.
In additionally, since by regularity Theorem4, the jump of u is zero through all the edges of Ω; hence, we clearly have JΓU−Uh,U−Uh=JΓUh,Uh. Thus,
(81)hE1/2ρgϕhnpEE+δhE1/2hufh−uph·nfEE≲ef,eph,wE+ep,λϕwE.(vii)
Finally, forE∈Eh∂K∩Ωp, we have(82)K−1uph+∇ϕh=K−1uph+∇ϕh−K−1up−∇ϕ=−K−1up−uph+∇ϕ−ϕh.
Thus,(83)K−1uph+∇ϕh×npEE≲ef,eph,K+εp,λϕK.
The estimates (66), (67), (70), (76), (77), (81), and (83) provide the desired local lower error bound.


## 3.1. Some Technical Results

Our a posteriori analysis requires some analytical results that are recalled. We define the space(30)H=v∈Hdiv,Ω:v∣Ωf∈Xfandv∣Ωp∈Xpwith the norm
(31)vH≔vf1,Ωf2+vpΩp2+∇·vpΩp.The first one concerns a sort of Helmholtz decomposition of elements ofH. Recall first that if d=3,
(32)H0curl,Ωp=ψ∈L2Ωp3:curlψ∈L2Ωp3andψ×n=0on∂Ωp.Theorem 3. (see Ref. [18], page 708).
Anyv∈H admits the Helmholtz-type decomposition
(33)v=v0+v1,where v0,v1∈H but satisfying v0∈H1Ωd,
(34)v1=0inΩf,curlβpinΩp,where βp∈H01Ωp if d=2, while βp∈H1Ωp3∩H0curl,Ωp if d=3, with the estimate
(35)v01,Ω+βp1,Ωp≲vH.The second result that we need is a regularity result for the solutionU=uf,p,up,ϕ∈H of (18) is the following theorem:Theorem 4. (see [18], page 710).
LetU∈H be the unique solution of (18). If K∈C0,1Ω¯pd×d, then there exists ε>0 such that
(36)u∣Ωp∈H1/2+εΩpd.Let us finish this section by an estimation of the stability error (see [18], Theorem 5):Theorem 5.
For anyUh=ufh,ph,uph,ϕh∈Hh, we have
(37)infWh∈Hh∩HUh−Whh2≲JΓUh,Uh.


## 3.2. Error Estimator

In order to solve the Stokes-Darcy coupled 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 and then the lower and upper error bounds are derived in Section3.3.


### 3.2.1. Error Equations

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. Thus, we define the error equations: letU=uf,p,up,ϕ∈H be the exact solution and Uh=ufh,ph,uph,ϕh∈Hh be the finite element solution. Then, for any Vh=vfh,qh,vph,ψh∈Hh and V=vf,p,vp,ψ∈H, using the Helmholtz decomposition (Theorem 3), we have
(38)LhU−Uh,V=LhU−Uh,V−Vh(39)=∑K∈ThfRKfUh,V−VhK+∑K∈ThpRKpUh,V−VhK,where
(40)RKfUh,V−VhK=ff+2ν∇·Dufh−∇ph,vf−vfhK−q−qh,∇·ufhK−∑E∈Eh∂K∩Γ∑j=1d−12νnf·Dufh·τj+ατj·K·τjufh·τj,vf−vfh·τjE+∑E∈Eh∂K∩Γph−2νnf·Dufh·nf−ρgϕh,vf−vfh·nfE,(41)RKpUh,V−VhK=curlρgK−1uph+∇ϕh,βp−βphK+ρgfp−∇·uph,ψ−ψhK−∑E∈Eh∂K∩ΩpρgK−1uph+∇ϕh×nE,ψ−ψhE+∑E∈Eh∂K∩ΩpρgϕhnEE,βp−βphE−∑E∈Eh∂K∩ΓρgϕhnEE,v−vhEE.


### 3.2.2. Residual Error Estimators

Definition 6 (a posteriori error indicators).
The residual error estimator is locally defined by(42)ΘK=ΘK,f2+ΘK,p21/2foreachK∈Th,where
(43)ΘK,f2=hK2ffh+2ν∇·Dufh−∇phK2+∇·ufhK2+∑E∈Eh∂K∩ΓhE∑j=1d−12νnf·Dufh·τj+ατj·K·τjufh·τjE2+∑E∈Eh∂K∩ΓhEph−2νnf·Dufh·nf−ρgϕhE2,(44)ΘK,p2=hK2curlρgK−1uph+∇ϕhK2+ρgfp−∇·uphK2+∑E∈Eh∂K∩ΩpρgK−1uph+∇ϕh×npEE2+∑E∈Eh∂K∩Ω¯phEρgϕhnpEE2+∑E∈Eh∂K∩ΓδhEhufh−uph·nfEE2.
The global residual error estimator is given by(45)Θ≔∑K∈ThΘK21/2.
Furthermore, denote the local and global approximation terms by(46)ζK≔hKff−ffhK,∀K∈Thf,ρgfp−fphK,∀K∈Thp,ζ≔∑K∈ThζK21/2,where the global function ffh:Ωf→ℝd is defined by
(47)ffh∣K=fK=1K∫Kffxdx,∀K∈Thf,while in Ωp, we take fh∣K=fK for all K∈Thp, as the unique element of ℙ1K such that
(48)∫KfKxqxdx=∫Kfxqxdx∀q∈ℙ1K.Remark 7.
The residual character of each term on the right-hand sides of (43) and (44) is quite clear since if ufh,ph,uph,ϕh∈Hh would be the exact solution of (18), then they would vanish.


### 3.2.3. Analytical Tools

(1) Inverse Inequalities. In order to derive the lower error bounds, we proceed similarly as in [25, 42] (see also [46]), 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 [41]). 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∂KE, and 0≤bE≤1inωE. We also recall from [47] 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 [47]).Lemma 8.
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
(49)qK≲qbK1/2K≲qK,∀q∈ℙkK,(50)qbK1,K≲hK−1qK,∀q∈ℙkK,(51)pE≲bE1/2pE≲pE,∀p∈ℙkE,(52)LpK+hELp1,K≲hE1/2pE,∀p∈ℙkE.(2) Continuous Trace Inequality.Lemma 9 (continuous trace inequality).
There exists a positive constantβ1>0 depending only on σ0 such that
(53)v∂K2≤β1vKv1,K,∀K∈Th,∀v∈H1Kd.(3) Clément Interpolation Operator. In order to derive the upper error bounds, we introduce the Clément interpolation operator ICl0:H01Ω→PcbTh that approximates optimally nonsmooth functions by continuous piecewise linear functions:
(54)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 [48], Section 3.Lemma 10.
There exist constantsC1,C2>0, independent of h, such that for all v∈H01Ω, there hold
(55)v−ICl0vK≤C1hKv1,ΔK,∀K∈Th,v−ICl0vE≤C2hE12v1,ΔE,∀E∈Eh,where ΔK≔∪K′∈Th:K′∩K≠∅ and ΔE≔∪K′∈Th:K′∩E≠∅.


## 3.2.1. Error Equations

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. Thus, we define the error equations: letU=uf,p,up,ϕ∈H be the exact solution and Uh=ufh,ph,uph,ϕh∈Hh be the finite element solution. Then, for any Vh=vfh,qh,vph,ψh∈Hh and V=vf,p,vp,ψ∈H, using the Helmholtz decomposition (Theorem 3), we have
(38)LhU−Uh,V=LhU−Uh,V−Vh(39)=∑K∈ThfRKfUh,V−VhK+∑K∈ThpRKpUh,V−VhK,where
(40)RKfUh,V−VhK=ff+2ν∇·Dufh−∇ph,vf−vfhK−q−qh,∇·ufhK−∑E∈Eh∂K∩Γ∑j=1d−12νnf·Dufh·τj+ατj·K·τjufh·τj,vf−vfh·τjE+∑E∈Eh∂K∩Γph−2νnf·Dufh·nf−ρgϕh,vf−vfh·nfE,(41)RKpUh,V−VhK=curlρgK−1uph+∇ϕh,βp−βphK+ρgfp−∇·uph,ψ−ψhK−∑E∈Eh∂K∩ΩpρgK−1uph+∇ϕh×nE,ψ−ψhE+∑E∈Eh∂K∩ΩpρgϕhnEE,βp−βphE−∑E∈Eh∂K∩ΓρgϕhnEE,v−vhEE.


## 3.2.2. Residual Error Estimators

Definition 6 (a posteriori error indicators).
The residual error estimator is locally defined by(42)ΘK=ΘK,f2+ΘK,p21/2foreachK∈Th,where
(43)ΘK,f2=hK2ffh+2ν∇·Dufh−∇phK2+∇·ufhK2+∑E∈Eh∂K∩ΓhE∑j=1d−12νnf·Dufh·τj+ατj·K·τjufh·τjE2+∑E∈Eh∂K∩ΓhEph−2νnf·Dufh·nf−ρgϕhE2,(44)ΘK,p2=hK2curlρgK−1uph+∇ϕhK2+ρgfp−∇·uphK2+∑E∈Eh∂K∩ΩpρgK−1uph+∇ϕh×npEE2+∑E∈Eh∂K∩Ω¯phEρgϕhnpEE2+∑E∈Eh∂K∩ΓδhEhufh−uph·nfEE2.
The global residual error estimator is given by(45)Θ≔∑K∈ThΘK21/2.
Furthermore, denote the local and global approximation terms by(46)ζK≔hKff−ffhK,∀K∈Thf,ρgfp−fphK,∀K∈Thp,ζ≔∑K∈ThζK21/2,where the global function ffh:Ωf→ℝd is defined by
(47)ffh∣K=fK=1K∫Kffxdx,∀K∈Thf,while in Ωp, we take fh∣K=fK for all K∈Thp, as the unique element of ℙ1K such that
(48)∫KfKxqxdx=∫Kfxqxdx∀q∈ℙ1K.Remark 7.
The residual character of each term on the right-hand sides of (43) and (44) is quite clear since if ufh,ph,uph,ϕh∈Hh would be the exact solution of (18), then they would vanish.


## 3.2.3. Analytical Tools

(1) Inverse Inequalities. In order to derive the lower error bounds, we proceed similarly as in [25, 42] (see also [46]), 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 [41]). 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∂KE, and 0≤bE≤1inωE. We also recall from [47] 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 [47]).Lemma 8.
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
(49)qK≲qbK1/2K≲qK,∀q∈ℙkK,(50)qbK1,K≲hK−1qK,∀q∈ℙkK,(51)pE≲bE1/2pE≲pE,∀p∈ℙkE,(52)LpK+hELp1,K≲hE1/2pE,∀p∈ℙkE.(2) Continuous Trace Inequality.Lemma 9 (continuous trace inequality).
There exists a positive constantβ1>0 depending only on σ0 such that
(53)v∂K2≤β1vKv1,K,∀K∈Th,∀v∈H1Kd.(3) Clément Interpolation Operator. In order to derive the upper error bounds, we introduce the Clément interpolation operator ICl0:H01Ω→PcbTh that approximates optimally nonsmooth functions by continuous piecewise linear functions:
(54)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 [48], Section 3.Lemma 10.
There exist constantsC1,C2>0, independent of h, such that for all v∈H01Ω, there hold
(55)v−ICl0vK≤C1hKv1,ΔK,∀K∈Th,v−ICl0vE≤C2hE12v1,ΔE,∀E∈Eh,where ΔK≔∪K′∈Th:K′∩K≠∅ and ΔE≔∪K′∈Th:K′∩E≠∅.


## 3.3. Optimality ofΘKK∈Th

### 3.3.1. Reliability Result

Theorem 11 (reliability ofΘ).
LetU=uf,p,up,ϕ∈H be the exact solution of (18) and Uh=ufh,ph,uph,ϕh∈Hh be the finite element solution of (23). There exist a constant Crel>0 such that the following estimate holds:
(56)U−Uhh≤CrelΘ+ζ.Proof.
We takeψh=0=qh in error equation (38)–(41) and we obtain
(57)LhU−Uh,V=LhU−Uh,V−Vh=∑K∈ThfRKfUh,V−VhK+∑K∈ThpRKpUh,V−VhK,where
(58)RKfUh,V−VhK=ffh+2ν∇·Dufh−∇ph,vf−vfhK−q,∇·ufhK+ff−ffh,vf−vfhK−∑E∈Eh∂K∩Γ∑j=1d−12νnf·Dufh·τj+ατj·K·τjufh·τj,vf−vfh·τjE+∑E∈Eh∂K∩Γph−2νnf·Dufh·nf−ρgϕh,vf−vfh·nfE,(59)RKpUh,V−VhK=curlρgK−1uph+∇ϕh,βp−βphK+ρgfph−∇·uph,ψK+ρgfp−fph,ψK−∑E∈Eh∂K∩ΩpρgK−1uph+∇ϕh×nE,ψE+∑E∈Eh∂K∩ΩpρgϕhnEE,βp−βphE−∑E∈Eh∂K∩ΓρgϕhnEE,v−vhEE.
The inf-sup condition ofLh leads to
(60)∥U−Uh∥h≤CsupV∈H\0∣LhU−Uh,V∣∥V∥h.
Now, using the error equation (57)–(59), Cauchy-Schwarz inequality, and the Clément operator of Lemma 10, we deduce the estimate (56). The proof is complete.


### 3.3.2. Efficiency Result

To prove local efficiency forw⊂Ω, let us denote by
(61)vf,vph,w2≔∑K⊂w¯∩Ω¯fvf1,K2+∑K⊂w¯∩Ω¯pvpK2+divhvpK2+vf×nΓ∩w¯2+∑K⊂Ω¯JKvf,vf,where
(62)JKvf,vf≔∑E∈EhΩ¯f∩EhKδhE−1vfEE2;ε,λw≔εw+λw.The main result of this subsection can be stated as follows.Theorem 12 (efficiency ofΘ).
Under the assumptions of Theorem4, the following lower error bound holds:
(63)ΘK≲ef,eph,w~K+εp,λϕw~K+∑K′⊂w~KζK′,where w~K is a finite union of neighboring elements of K.Proof.
We begin by bounding each term of the residuals separately.(i)
Element residual inΩf: to estimate hK2ffh+2ν∇·Dufh−∇phK2, we choose in error equation (57)–(59) for each K∈Thf, V=vsK,0,vpK,0 and Vh=0,0,0,0 with vpK=0 on Ωp,(64)vfK=ffh+2ν∇·Dufh−∇phbKonK∈Thf,0onΩf\K,for obtained, LhU−Uh,V=ffh+2ν∇·Dufh−∇phbK1/2K2+ff−ffh,vfKK. Noted that ffh+2ν∇·Dufh−∇phbK∈ℙkK and vanish on ∂K. Because LhU−Uh,V=Def,vfKK in this case, then we have Def,vfKK=ffh+2ν∇·Dufh−∇phbK1/2K2+ff−ffh,vfKK. The first inverse inequality (49) and the Cauchy-Schwarz inequality lead to
(65)ffh+2ν∇·Dufh−∇phK2~ffh+2ν∇·Dufh−∇phbK1/2K2=∫KDefffh+2ν∇·Dufh−∇phbK−ff−ffh,vfKK≤efKffh+2ν∇·Dufh−∇phbK1,K+ff−ffhKvfKK.
Using inverse inequality (50), we deduce the estimate:
(66)hKffh+2ν∇·Dufh−∇phK≲ef,eph,w~K+ζK.(ii)
Divergence element residual inΩf (estimation of ∇·ufhK2): for each K∈Thf, we have,(67)∇ufhK=∇uf−ufhK≲uf−ufh1,K.(iii)
Element residual inΩp: we have for each K∈Thp,(68)ρgfph−∇·uphK=ρgfph−∇·uph−ρgfp−∇·upK=ρgdivhup−uph−ρgfp−fphK≲ef,eph,w~K+ζK.(iv)
Curl element residual inΩp: for K∈Thp, we set CK=curlρgK−1uph+∇ϕh and WK=CKbK. Hence, we notice that curlWK belongs to H and is divergence free; therefore, by Equations (57)–(59), we obtain with Vh=0 and βp=WK, LhU−Uh,WK=RKpUh,WKK=∥curlρgK−1uph+∇ϕhbK1/2∥K2. The first inverse inequality (49) and the Cauchy-Schwarz inequality lead to(69)curlρgK−1uph+∇ϕhK2~curlρgK−1uph+∇ϕhbK1/2K2=∫KρgK−1up−uph+∇ϕ−ϕh·curlWK≤ρgK−1up−uphK+∇ϕ−ϕhK·curlWKK≲ρgK−1epK+∇λϕK·WKK.
Again, the inverse inequality (49) allows to get
(70)curlρgK−1uph+∇ϕhK≲ef,eph,w~K+εp,λϕw~K.(v)
Interface elements onΓ: we fix an edge E included in Γ and for a constant rE fixed later on a unit vector N, we consider(71)WE=rEbEN,that clearly belongs to H. Hence, by residual equation (57)–(59), we obtain with Vh=0,
(72)LhU−Uh,WE=RKffUh,WEKf+RKppUh,WEKp,
whereKf (resp., Kp) is the unique triangle/tetrahedron included in Ω¯f (resp., Ω¯p) having E as edge/face, and
(73)RKffUh,WEKf=ff+2ν∇·Dufh−∇ph,WEKf−q,∇·ufhKf−∑E∈Eh∂Kf∩Γ∑j=1d−12νnf·Dufh·τj+ατj·K·τjufh·τj,WE·τjE+∑E∈Eh∂Kf∩Γph−2νnf·Dufh·nf−ρgϕh,WE·nfE,RKpUh,WEKp=curlρgK−1uph+∇ϕh,βpKp+ρgfp−∇·uph,ψKp−∑E∈Eh∂Kp∩ΩpρgK−1uph+∇ϕh×nE,ψE+∑E∈Eh∂K∩ΩpρgϕhnEE,βpE−∑E∈Eh∂K∩ΓρgϕhnEE,WEE.
TakenWE=0 in Kp, q=0 in Ω and for each j=1,⋯,d−1, rE=2νnf·Dufh·τj+α/√τj·K·τjufh·τj with N=τj. We have
(74)LhU−Uh,WE=ff+2ν∇·Dufh−∇ph,WEKf−2νnf·Dufh·τj+ατj·K·τjufh·τj,WE·τjE=ff+2ν∇·Dufh−∇ph,WEKf−∥rEbE1/2∥E2.
Hence,(75)rEE2~rEbE1/2E2=ff+2ν∇·Dufh−∇ph,WEKf−LhU−Uh,WE=ffh+2ν∇·Dufh−∇ph,WEKf−∫Kf2νDef:DWE+∫KfεpdivWE+ff−ffh,WEKf−∑j=1d−1ατj·K·τjef·τj,WE·τjE.
Inverse inequalities (50) and (51) and Cauchy-Schwarz inequality lead to
(76)hE1/22νnf·Dufh·τj+ατj·K·τjufh·τjE≲ef,eph,wE+εp,λϕwE+∑K⊂wEζK,with wE=Kf∪Kp.
TakenWE=0 in Kf, q=0 in Ω and for each j=1,⋯,d−1, rE=2νnf·Dufh·τj+α/√τj·K·τjufh·τj with N=nf. As before the identities (72)–(74) and the inverse inequalities (50) and (51) lead to
(77)hE1/2∥ph−2νnf·Dufh·nf−ρgϕh∥E≲∥ef,ep∥h,wE+∥εp,λϕ∥wE+∑K⊂wEζK.(vi)
Piezometric head jump inΩp: for each edge/face E∈EhΩp, we consider wE=T1∪T2. As ρgϕhnpE∈ℙ1Ed, we set(78)WE≔ρgϕhnpEbE∈H01wEd.
Using the residual equation (57)–(59), we obtain with Vh=0 and V=0,0,βp,0 where βp=WE:
(79)LhU−Uh,V=ρgϕhnEE,WEE.
Therefore, Cauchy-Schwarz inequality and inverse inequalities (51) and (52) lead to
(80)ρgϕhnEEE2~ρgϕhnEEbE1/2E2=∫EρgϕhnEE·WE=∑i=12ρgϕ−ϕh,∇·WEKi≤∑i=12ρgλϕKi∇·WEKi≲∑i=12λϕKihE−1ρgϕhnEEE.
In additionally, since by regularity Theorem4, the jump of u is zero through all the edges of Ω; hence, we clearly have JΓU−Uh,U−Uh=JΓUh,Uh. Thus,
(81)hE1/2ρgϕhnpEE+δhE1/2hufh−uph·nfEE≲ef,eph,wE+ep,λϕwE.(vii)
Finally, forE∈Eh∂K∩Ωp, we have(82)K−1uph+∇ϕh=K−1uph+∇ϕh−K−1up−∇ϕ=−K−1up−uph+∇ϕ−ϕh.
Thus,(83)K−1uph+∇ϕh×npEE≲ef,eph,K+εp,λϕK.
The estimates (66), (67), (70), (76), (77), (81), and (83) provide the desired local lower error bound.


## 3.3.1. Reliability Result

Theorem 11 (reliability ofΘ).
LetU=uf,p,up,ϕ∈H be the exact solution of (18) and Uh=ufh,ph,uph,ϕh∈Hh be the finite element solution of (23). There exist a constant Crel>0 such that the following estimate holds:
(56)U−Uhh≤CrelΘ+ζ.Proof.
We takeψh=0=qh in error equation (38)–(41) and we obtain
(57)LhU−Uh,V=LhU−Uh,V−Vh=∑K∈ThfRKfUh,V−VhK+∑K∈ThpRKpUh,V−VhK,where
(58)RKfUh,V−VhK=ffh+2ν∇·Dufh−∇ph,vf−vfhK−q,∇·ufhK+ff−ffh,vf−vfhK−∑E∈Eh∂K∩Γ∑j=1d−12νnf·Dufh·τj+ατj·K·τjufh·τj,vf−vfh·τjE+∑E∈Eh∂K∩Γph−2νnf·Dufh·nf−ρgϕh,vf−vfh·nfE,(59)RKpUh,V−VhK=curlρgK−1uph+∇ϕh,βp−βphK+ρgfph−∇·uph,ψK+ρgfp−fph,ψK−∑E∈Eh∂K∩ΩpρgK−1uph+∇ϕh×nE,ψE+∑E∈Eh∂K∩ΩpρgϕhnEE,βp−βphE−∑E∈Eh∂K∩ΓρgϕhnEE,v−vhEE.
The inf-sup condition ofLh leads to
(60)∥U−Uh∥h≤CsupV∈H\0∣LhU−Uh,V∣∥V∥h.
Now, using the error equation (57)–(59), Cauchy-Schwarz inequality, and the Clément operator of Lemma 10, we deduce the estimate (56). The proof is complete.


## 3.3.2. Efficiency Result

To prove local efficiency forw⊂Ω, let us denote by
(61)vf,vph,w2≔∑K⊂w¯∩Ω¯fvf1,K2+∑K⊂w¯∩Ω¯pvpK2+divhvpK2+vf×nΓ∩w¯2+∑K⊂Ω¯JKvf,vf,where
(62)JKvf,vf≔∑E∈EhΩ¯f∩EhKδhE−1vfEE2;ε,λw≔εw+λw.The main result of this subsection can be stated as follows.Theorem 12 (efficiency ofΘ).
Under the assumptions of Theorem4, the following lower error bound holds:
(63)ΘK≲ef,eph,w~K+εp,λϕw~K+∑K′⊂w~KζK′,where w~K is a finite union of neighboring elements of K.Proof.
We begin by bounding each term of the residuals separately.(i)
Element residual inΩf: to estimate hK2ffh+2ν∇·Dufh−∇phK2, we choose in error equation (57)–(59) for each K∈Thf, V=vsK,0,vpK,0 and Vh=0,0,0,0 with vpK=0 on Ωp,(64)vfK=ffh+2ν∇·Dufh−∇phbKonK∈Thf,0onΩf\K,for obtained, LhU−Uh,V=ffh+2ν∇·Dufh−∇phbK1/2K2+ff−ffh,vfKK. Noted that ffh+2ν∇·Dufh−∇phbK∈ℙkK and vanish on ∂K. Because LhU−Uh,V=Def,vfKK in this case, then we have Def,vfKK=ffh+2ν∇·Dufh−∇phbK1/2K2+ff−ffh,vfKK. The first inverse inequality (49) and the Cauchy-Schwarz inequality lead to
(65)ffh+2ν∇·Dufh−∇phK2~ffh+2ν∇·Dufh−∇phbK1/2K2=∫KDefffh+2ν∇·Dufh−∇phbK−ff−ffh,vfKK≤efKffh+2ν∇·Dufh−∇phbK1,K+ff−ffhKvfKK.
Using inverse inequality (50), we deduce the estimate:
(66)hKffh+2ν∇·Dufh−∇phK≲ef,eph,w~K+ζK.(ii)
Divergence element residual inΩf (estimation of ∇·ufhK2): for each K∈Thf, we have,(67)∇ufhK=∇uf−ufhK≲uf−ufh1,K.(iii)
Element residual inΩp: we have for each K∈Thp,(68)ρgfph−∇·uphK=ρgfph−∇·uph−ρgfp−∇·upK=ρgdivhup−uph−ρgfp−fphK≲ef,eph,w~K+ζK.(iv)
Curl element residual inΩp: for K∈Thp, we set CK=curlρgK−1uph+∇ϕh and WK=CKbK. Hence, we notice that curlWK belongs to H and is divergence free; therefore, by Equations (57)–(59), we obtain with Vh=0 and βp=WK, LhU−Uh,WK=RKpUh,WKK=∥curlρgK−1uph+∇ϕhbK1/2∥K2. The first inverse inequality (49) and the Cauchy-Schwarz inequality lead to(69)curlρgK−1uph+∇ϕhK2~curlρgK−1uph+∇ϕhbK1/2K2=∫KρgK−1up−uph+∇ϕ−ϕh·curlWK≤ρgK−1up−uphK+∇ϕ−ϕhK·curlWKK≲ρgK−1epK+∇λϕK·WKK.
Again, the inverse inequality (49) allows to get
(70)curlρgK−1uph+∇ϕhK≲ef,eph,w~K+εp,λϕw~K.(v)
Interface elements onΓ: we fix an edge E included in Γ and for a constant rE fixed later on a unit vector N, we consider(71)WE=rEbEN,that clearly belongs to H. Hence, by residual equation (57)–(59), we obtain with Vh=0,
(72)LhU−Uh,WE=RKffUh,WEKf+RKppUh,WEKp,
whereKf (resp., Kp) is the unique triangle/tetrahedron included in Ω¯f (resp., Ω¯p) having E as edge/face, and
(73)RKffUh,WEKf=ff+2ν∇·Dufh−∇ph,WEKf−q,∇·ufhKf−∑E∈Eh∂Kf∩Γ∑j=1d−12νnf·Dufh·τj+ατj·K·τjufh·τj,WE·τjE+∑E∈Eh∂Kf∩Γph−2νnf·Dufh·nf−ρgϕh,WE·nfE,RKpUh,WEKp=curlρgK−1uph+∇ϕh,βpKp+ρgfp−∇·uph,ψKp−∑E∈Eh∂Kp∩ΩpρgK−1uph+∇ϕh×nE,ψE+∑E∈Eh∂K∩ΩpρgϕhnEE,βpE−∑E∈Eh∂K∩ΓρgϕhnEE,WEE.
TakenWE=0 in Kp, q=0 in Ω and for each j=1,⋯,d−1, rE=2νnf·Dufh·τj+α/√τj·K·τjufh·τj with N=τj. We have
(74)LhU−Uh,WE=ff+2ν∇·Dufh−∇ph,WEKf−2νnf·Dufh·τj+ατj·K·τjufh·τj,WE·τjE=ff+2ν∇·Dufh−∇ph,WEKf−∥rEbE1/2∥E2.
Hence,(75)rEE2~rEbE1/2E2=ff+2ν∇·Dufh−∇ph,WEKf−LhU−Uh,WE=ffh+2ν∇·Dufh−∇ph,WEKf−∫Kf2νDef:DWE+∫KfεpdivWE+ff−ffh,WEKf−∑j=1d−1ατj·K·τjef·τj,WE·τjE.
Inverse inequalities (50) and (51) and Cauchy-Schwarz inequality lead to
(76)hE1/22νnf·Dufh·τj+ατj·K·τjufh·τjE≲ef,eph,wE+εp,λϕwE+∑K⊂wEζK,with wE=Kf∪Kp.
TakenWE=0 in Kf, q=0 in Ω and for each j=1,⋯,d−1, rE=2νnf·Dufh·τj+α/√τj·K·τjufh·τj with N=nf. As before the identities (72)–(74) and the inverse inequalities (50) and (51) lead to
(77)hE1/2∥ph−2νnf·Dufh·nf−ρgϕh∥E≲∥ef,ep∥h,wE+∥εp,λϕ∥wE+∑K⊂wEζK.(vi)
Piezometric head jump inΩp: for each edge/face E∈EhΩp, we consider wE=T1∪T2. As ρgϕhnpE∈ℙ1Ed, we set(78)WE≔ρgϕhnpEbE∈H01wEd.
Using the residual equation (57)–(59), we obtain with Vh=0 and V=0,0,βp,0 where βp=WE:
(79)LhU−Uh,V=ρgϕhnEE,WEE.
Therefore, Cauchy-Schwarz inequality and inverse inequalities (51) and (52) lead to
(80)ρgϕhnEEE2~ρgϕhnEEbE1/2E2=∫EρgϕhnEE·WE=∑i=12ρgϕ−ϕh,∇·WEKi≤∑i=12ρgλϕKi∇·WEKi≲∑i=12λϕKihE−1ρgϕhnEEE.
In additionally, since by regularity Theorem4, the jump of u is zero through all the edges of Ω; hence, we clearly have JΓU−Uh,U−Uh=JΓUh,Uh. Thus,
(81)hE1/2ρgϕhnpEE+δhE1/2hufh−uph·nfEE≲ef,eph,wE+ep,λϕwE.(vii)
Finally, forE∈Eh∂K∩Ωp, we have(82)K−1uph+∇ϕh=K−1uph+∇ϕh−K−1up−∇ϕ=−K−1up−uph+∇ϕ−ϕh.
Thus,(83)K−1uph+∇ϕh×npEE≲ef,eph,K+εp,λϕK.
The estimates (66), (67), (70), (76), (77), (81), and (83) provide the desired local lower error bound.


## 4. Summary

In this paper, we have discussed a posteriori error estimates for a finite element approximation of the Stokes-Darcy 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.


## 5. Nomenclature

(i)
Ω
⊂
ℝ
d
,
d
∈
2
,
3 bounded domain
(ii)
Ω
p
: the porous medium domain
(iii)
Ω
f
: the fluid region
(iv)
Γ
=
∂
Ω
f
∩
∂
Ω
p
(v)
Ω
f
=
Ω
\
Ω
p
∪
Γ
(vi)
Γ
l
=
∂
Ω
l
\
Γ
,
l
=
f
,
p
(vii)
n
f (resp., np) the unit outward normal vector along ∂Ωf (resp., ∂Ωp)
(viii)
u: the fluid velocity
(ix)
p: the fluid pressure
(x)
In 2D, thecurl of a scalar function w is given as usual by
(84)
curl
w
≔
∂
w
/
∂
x
2
,
−
∂
w
/
∂
x
1
Τ
(xi)
In 3D, thecurl of a vector function w=w1,w2,w3 is given as usual by curlw≔∇×w, namely,
(85)
curl
w
≔
∂
w
3
/
∂
x
2
−
∂
w
2
/
∂
x
3
,
∂
w
1
/
∂
x
3
−
∂
w
3
/
∂
x
1
,
∂
w
2
/
∂
x
1
−
∂
w
1
/
∂
x
2
(xii)
ℙ
k: the space of polynomials of total degree not larger than k
(xiii)
T
h: triangulation of Ω
(xiv)
T
h
l: the corresponding induced triangulation of Ωl, l∈f,p
(xv)
For anyK∈Th, hK is the diameter of K and ρK=2rK is the diameter of the largest ball inscribed into K
(xvi)
h
≔
max
K
∈
T
h
h
K and σh≔maxK∈ThhK/ρK
(xvii)
E
h: the set of all the edges or faces of the triangulation
(xviii)
E
K: the set of all the edges (N=2) or faces (N=3) of a element K
(xix)
E
h
≔
⋃
K
∈
T
h
E
K
(xx)
N
K: the set of all the vertices of a element K
(xxi)
N
h
≔
⋃
K
∈
T
h
N
K
(xxii)
ForA⊂Ω¯, EhA≔E∈Eh:E⊂A
(xxiii)
ForE∈Eh, we associate a unit vector nE such that nE is orthogonal to E and equals to the unit exterior normal vector to ∂Ω
(xxiv)
ForE∈Eh, ϕE is the jump across E in the direction of nE
(xxv)
In order to avoid excessive use of constants, the abbreviationsx≲y and x~y stand for x≤cy and c1x≤y≤c2x, respectively, with positive constants independent of x, y, or Th



---

*Source: 8628739-2020-10-12.xml*