# A posteriori error analysis for a new fully-mixed isotropic discretization to the stationary Stokes-Darcy coupled problem

In this paper we develop an a posteriori error analysis for the stationary Stokes-Darcy coupled problem approximated by conforming finite element method on isotropic meshes in R^d, d∈{2,3}. The approach utilizes a new robust stabilized fully mixed discretization developed by Jiaping Yu et al. (Advances in Difference Equations, SpringerOpen Journal, 2018). 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.

There are no comments yet.

## Authors

• 4 publications
12/03/2019

### Residual-based a posteriori error estimates for a conforming mixed finite element discretization of the Monge-Ampère equation

In this paper we develop a new a posteriori error analysis for the Monge...
04/21/2020

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

In this work we develop an a posteriori error analysis of a conforming m...
04/21/2020

### An a posteriori error analysis for the equations of stationary incompressible magnetohydrodynamics

Magnetohydrodynamics (MHD) is a continuum level model for conducting flu...
06/04/2020

### Explicit a posteriori and a priori error estimation for the finite element solution of Stokes equations

For the Stokes equation over 2D and 3D domains, explicit a posteriori an...
05/05/2021

### Solvability of Discrete Helmholtz Equations

We study the unique solvability of the discretized Helmholtz problem wit...
08/18/2020

### An a posteriori error estimate of the outer normal derivative using dual weights

We derive a residual based a-posteriori error estimate for the outer nor...
06/18/2018

### Quantifying discretization errors for soft-tissue simulation in computer assisted surgery: a preliminary study

Errors in biomechanics simulations arise from modeling and discretizatio...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 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 [24] 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 treatement 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 [8], modified by Saffman in [51], and later mathematically justified in [36, 37, 38, 47].

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 examples [2, 26, 29, 30, 41, 57, 48, 49, 58] for some mathematical analysis. The authors in [57] studied two different mixed formulations: the first one imposes the weak continuity of the normal component of the velocity field on the interface, by using a Lagrange multiplier; while the second one imposes the strong continuity in the functional space. Later on we 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 separately. The works in [50, 26, 29, 30, 57, 44, 35] are based on the weakly coupled formulation. Researches on the strongly coupled formulation have been focused on the development of an unified discretization, that is, 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 the Stokes problems. In [2, 3], a conforming, unified finite element has been proposed for the strongly coupled mixed formulation. Superconvergence analysis of the finite element methods for the Stokes-Darcy system was studied in [15]. Other less restrictive discretizations as the non-conforming unified approach [50, 41] or the discontinuous Galerkin approach have been proposed in [39, 48, 49]. Due to its discontinuous nature, some discretizations for the coupled Darcy-Stokes problem may break the strong coupling in the discrete level [48, 49], as they impose the normal continuity across the interface via interior penalties.

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 Babuska and Rheinboldt [6], adaptive finite element methods based on a posteriori error estimates have been extensively investigated.

An a priori error analysis is performed with some numerical tests confirming the convergence rates.

A posteriori error estimations have been well-established for both the mixed formulation of the Darcy flow [12, 10, 42], and the Stokes flow [1, 7, 14, 22, 25, 34, 46, 52, 55, 56]. However, only few works exist for the coupled Darcy-Stokes problem, see for instance [16, 5, 21, 28, 44]. The paper [16, 44] concern the strongly coupled mixed formulation where a conforming and nonconforming finite element methods have been used and [21, 5] concern the weakly coupled mixed formulation while [21] uses the primal formulation on the Darcy side. The authors in [28] employ a fully-mixed formulation where Raviart-Thomas elements have been used to approximate the velocity in both the Stokes domain and Darcy domain, and constant piecewise for approximate the pressure.

In [58], 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-know MINI elements () to approximate the velocity and pressure in the conduit for Stokes equation. To capture the fully mixed technique in the porous medium region linear Lagrangian elements, have been used for hydraulic (piezometric) head and Brezzi-Douglas-Marini () piecewise constant finite elements have been used for Darcy velocity. To our best knowledge, there is no a posteriori error estimation for the fully-mixed discretization proposed in [58]. Here we develop such a posteriori error analysis. The a posteriori error estimate is based on a suitable evaluation on the residual of the finite element solution. We further prove that our a posteriori error estimator is both reliable and efficient. The difference between our paper and the reference [28] is that our discretization uses MINI elements () to approximate the velocity and pressure in the conduit for Stokes equations, -Lagrange elements to approximate hydraulic (piezometric) head and Brezzi-Douglas-Marini () piecewise constant finite elements have been used for Darcy velocity. As a result, additional term is included in the error estimator that measure the stability of the method. In order to treat appropriately this stability term, we further need a special Helmholtz decomposition [44, Theorem 3.1], a regularity result [44, Theorem 3.2] and an estimate of the stability error [44, Theorem 3.3].

The paper is organized as follows. Some preliminaries and notation are given in section 2. The efficiency result is derived using the technique of bubble function introduced by R. Verfürth [54] and used in similar context by C. Carstensen [12, 13]. In section 3, the a posteriori error estimates are derived.

## 2. Preliminaries and Notations

### 2.1. Model problem

We consider the model of a flow in a bounded domain , consisting of a porous medium domain , where the flow is a Darcy flow, and an open region where the flow is governed by the Stokes equations. The two regions are separated by an interface Let , . Each interface and boundary is assumed to be polygonal or polyhedral . We denote by (resp.

) the unit outward normal vector along

(resp. ). Note that on the interface , we have . The Figure 1 shows a sketch of the problem domain, its boundaries and some other notations.

The fluid velocity and pressure and are governed by the Stikes equations in :

 {−2ν∇⋅D(uf)+∇p=ff in Ωf,∇⋅uf=0 in Ωf, (1)

where

denotes the stress tensor, and

represents the deformation tensor. The porous media flow is governed by the following Darcy equations on through the fluid velocity and the piezometric head :

 {up=−K∇ϕ in Ωp,∇⋅up=fp in Ωp. (2)

We impose impermeable boundary conditions, on , on the exterior boundary of the porous media region, and no slip conditions, on , 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 competely 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 [36, 37, 38, 47, 51, 8], interfacial coupling:

 uf⋅nf+up⋅np = 0 on Γ (3) −nf⋅T⋅nf=p−2νn% f⋅D(uf)⋅nf = ρgϕ on Γ (4) −nf⋅T⋅τj=−2nf⋅D(uf)⋅τj = α√τj⋅Kτjuf⋅τj,1≤j≤(d−1) on Γ. (5)

This is a simplification of the original and more physically relistic Beavers-Joseph conditions (in which in (2.8) is replaced by ; see [8] ). Here we denote
, -body forces in the fluid region and source in the porous region,
K-symmetric positive define (SPD) hydraulic conductivity tensor,
-constant parameter.

We shall also assume that all material and fluid parameters defined above are uniformly positive and bounded, i.e.,

 0≤kmin≤λ(K)≤kmax<∞.

### 2.2. Notations and the weak formulation

In this part, we first introduce some Sobolev spaces [43] and norms. If is a bounded domain of and is a non negative integer, the Sobolev space is defined in the usual way with the usual norm and semi-norm . In particular, and we write for . Similarly we denote by the or inner product. For shortness if is equal to , we will drop the index , while for any , , and , for . The space denotes the closure of in . Let be the space of vector valued functions with components in . The norm and the seminorm on are given by

 ∥v∥m,Ω:=(N∑i=0∥vi∥2m,Ω)1/2 and |v|m,Ω:=(N∑i=0|vi|2m,Ω)1/2. (6)

For a connected open subset of the boundary , we write for the inner product (or duality pairing), that is, for scalar valued functions , one defines:

 ⟨λ,σ⟩Γ:=∫Γλσds (7)

By setting the space

 Hdiv:=H(div;Ωp)={vp∈[L2(Ωp)]d:∇⋅vp∈L2(Ωp)},

we introduce the following spaces:

 Xf := {vf∈[L2(Ωf)]d:vf=0 on Γf}, Qf := L2(Ωf), Xp := {vp∈H(div;Ωp):vp⋅np=0 on Γp}, Qp := L2(Ωp).

For the spaces and , we define the following norms:

 ∥vf∥1 := √∥vf∥2Ωf+|% vf|21,Ωf, with |vf|1,Ωf=∥∇vf∥Ωf∀vf∈Xf, ∥vp∥div := √∥vp∥2Ωp+∥∇⋅vp∥2Ωp,∀vp∈Xp.

The variational formulation of the steady-state Stokes-Darcy problem (1)-(5) reads as: Find satisfying:

 af(uf,vf)−bf(vf,p)+cΓ(vf,ϕ) = (ff,vf)Ωf ∀vf∈Xf, (8) bf(uf,q) = 0 ∀q∈Qf, (9) ap(up,vp)−bp(vp,ϕ)−cΓ(vp,ϕ) = 0 ∀vp∈Xp, (10) bp(up,ψ) = ρg(fp,ψ)Ωp ψ∈Qp, (11)

where the bilinear forms are defined as:

 af(uf,vf) := 2ν(D(uf),D(vf))Ωf+d−1∑j=1α√τj⋅Kτj⟨uf⋅τj,vf⋅τj⟩Γ ap(up,vp) := ρg(K−1up,vp)Ωp, bf(vf,p) := (p,∇⋅vf)Ωf, bp(vp,ϕ) := ρg(ϕ,∇⋅vp)Ωp, cΓ(vf,ϕ) := ρg⟨ϕ,vf⋅nf⟩Γ.

After introducing, for and ,

 L(U,V) := af(uf,vf)−bf(vf,p)+bf(uf,q) + ap(up,vp)−bp(vp,ϕ)+bp(up,ψ)+cΓ(vf−vp,ϕ), F(V) := (ff,vf)Ωf+ρg(fp,ψ)Ωp, (13)

the weak formulation (8)-(11) can be equivalentl rewritten as follows: Find satisfying

 L(U,V)=F(V),∀V∈H. (14)

It is easy to verify that this variational fromulation is well-posedness.

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 constant , , , only depending on such that for all ,

 ∥vf∥≤Cp|vf|1, |vf|≤CK∥D(vf)∥Ωf, ∥vf∥Γ≤Cv∥vf∥1/2Ωf|vf|1/21,Ωf.

Besides, there exists a constant that only depends on such that for all ,

 ∥ψ∥L2(Γ)≤~Cv∥ψ∥1/2Ωp|ψ|1/21,Ωp. (15)

### 2.3. Fully-mixed isotropic discretization

First, we consider the family of triangulations of , consisting of and , which are regular triangulations of and , respectively, where is a positive parameter. We also assume that on the interface the two meshes of and , which form the regular triangulation coincide.

The domain of the uniformly regular triangulation is such that and . There exist positive constants and satisfying . To approximate the diameter of the trangle (or tetrahedral) , is the diameter of the greatest ball included in . Based on the subdivisions and , we can define finite element spaces , , , . We consider the well-known MINI elements to approximate the velocity and the pressure in the conduit for Stokes equations [4]. To capture the fully-mixed technique in the porous medium region linear Lagrangian elements, are used for hydraulic (piezometric) head and Brezzi-Douglas-Marini () piecewise constant finite elements are used for Darcy velocity [11].
In the fluid region, we select for the Stokes problem the finite element spaces that satisfy the velocity-pressure inf-sup condition: There exists a constant , independent of , such that,

 inf0≠qh∈Qfhsup0≠vhf∈Xfhbf(vhf,qh)|vhf|1,Ωf∥qhf∥Ωf≥βf. (16)

In the porous region, we use the finite element spaces that also satisfy a standard inf-sup condition: There exist a constant such that for all ,

 inf0≠ϕh∈Qphsup0≠vhp∈Xphbp(vhp,ϕh)∥vhf∥div∥ϕh∥Ωp≥βp. (17)

Then the finite element discretization of (14) is to find such that

 L(Uh,Vh)+JΓ(Uh,Vh)=F(Vh) ∀Vh∈Hh. (18)

This is the natural discretization of the weak formulation (14) except that the stabilized term is added. This bilinear form is defined by

 JΓ(Uh,Vh):=δh⟨(uhf−uhp)⋅nf,(v% hf−vhp)⋅nf⟩Γ. (19)

We are now able to define the norm on :

 ∥V∥h:=√∥vhf∥21,Ωf+∥qhf∥2Ωf+∥vhp∥2div+∥ψh∥2Ωp+h−1∥(vhf−% vhp)∥2Γ

We have the following results (see [58, Theorem 2 and Theorem 3]):

###### Theorem 2.1.

There exists a unique solution to problem (18) and if the solution of the continuous problem (14) is smooth enough, then we have:

 ∥U−Uh∥h≤C(U)h. (20)

Below, in order to avoid excessive use of constants, the abbreviation stand for , with a positive constant independent of , and .

###### Remark 2.1.

(Galerkin orthogonality relation) Let be the exact solution and be the finite element solution. Then for any , we can subtract (14) to (18) to obtain the Galerkin orthogonality relation:

 Lh(U−Uh,Vh) = Lh(U,Vh)−Lh(Uh,Vh) = L(U,Vh)−Lh(Uh,Vh) = F(Vh)−F(Vh) = 0

thus, we have the relation:

 2ν(D(ef),D(vfh))Ωf + d−1∑j=1α√τj⋅K⋅τj⟨ef⋅τj,vfh⋅τj⟩Γ−(ϵp,∇⋅vfh)Ωf+(qh,∇⋅ef)Ωf + ρg[(K−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 erors in the hydraulic and Darcy velocity equations are respectively defined by:

 ef:=uf−ufh,ϵp:=p−ph,ep:=up−% uph and λϕ=ϕ−ϕ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

 H={v∈H(div,Ω):v|Ωf∈Xf and v|Ωp∈Xp}

with the norm

 ∥v∥H:=√|vf|21,Ωf+∥vp∥2Ωp+∥∇⋅vp∥Ωp.

The first one concerns a sort of Helmholtz decomposition of elements of . Recall first that if ,

 H0(curl,Ωp)={ψ∈L2(Ωp)3:curlψ∈L2(Ωp)3 and ψ×n=0 on ∂Ωp}.
###### Theorem 3.1.

(Ref. [44, Page 708]) Any admits the Helmholtz type decomposition

 v=v0+v1, (21)

where but satisfying ,

 v1={0 in Ωs,curlψ in Ωd, (22)

where if , while if , with the estimate

 ∥v0∥1,Ω+∥ψ∥1,Ωp≲∥v∥H. (23)

The second result that we need is a regularity result for the solution of (14) is the following theorem:

###### Theorem 3.2.

([44, Page 710]) Let be the unique solution of (14). If , then there exists such that

 u|Ωp∈[H12+ϵ(Ωp)]d.

Let us finish this section by an estimation of the stability error (see [44, Theorem 3.3]):

###### Theorem 3.3.

For any we have

 infWh∈Hh∩H∥Uh−% Wh∥2h≲JΓ(Uh,Uh). (24)

### 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 Section 3.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: Let be the exact solution and be the finite element solution. Then for any and , we have:

 Lh(U−Uh,V) = Lh(U−Uh,V−Vh) = ⎡⎢ ⎢⎣∑K∈Tfh(RfK(Uh),V−Vh))K+∑K∈Tph(RpK(Uh),%V−Vh))K⎤⎥ ⎥⎦,

where

 (RfK(Uh),V−Vh)K = (ffh+2ν∇⋅D(ufh)−∇ph,vf−vfh)K−(q−qh,∇⋅ufh)K − ⎡⎢ ⎢⎣∑E∈Eh(∂K∩Γ)d−1∑j=1⎛⎜ ⎜⎝2νnf⋅D(%ufh)⋅τj+α√τj⋅K⋅τjufh⋅τj,(vf−vfh)⋅τj⎞⎟ ⎟⎠E⎤⎥ ⎥⎦ +

and

 (RpK(Uh),V−Vh)K = (curl(ρgK−1uph+∇ϕh),vp−vph)K + (ρg(fp−∇⋅uph),ψ−ψh)K − ∑EEh(∂K∩Ωp)(ρgK−1uph+∇ϕh)×nE,ψ−ψh)E + ∑E∈Eh(∂K∩Ωp)([ρgϕhnE]E,vp−vph)E − ∑EEh(∂K∩Γ)([ρgϕhnE]E,[v−vh]E)E.

#### 3.2.2. Residual Error Estimators

###### Definition 3.1.

(A posteriori error indicators)

1. Error indicators in Stokes domain : For , we set:

 Θ2K,f = h2K∥ffh+2ν∇⋅D(ufh)−∇ph∥2K+∥∇⋅ufh∥2K + ∑E∈Eh(∂K∩Γ)hE⎧⎪ ⎪⎨⎪ ⎪⎩d−1∑j=1∥2νnf⋅D(ufh)⋅τj+α√τj⋅K⋅τjufh⋅τj∥2E⎫⎪ ⎪⎬⎪ ⎪⎭ + ∑E∈Eh(∂K∩Γ)hE∥ph−2νnf⋅D(ufh)⋅nf−ρgϕh∥2E
2. Error indicators in Darcy domain: For , we set:

 Θ2K,p = h2K∥curl(ρgK−1uph+∇ϕh)∥2K+h2K∥ρg(fp−∇⋅uph)∥2K + ∑EEh(∂K∩Ωp)hE∥[(ρgK−1uph+∇ϕh)×%np]E∥2E + ∑EEh(∂K∩Ωp)hE∥[ρgϕhnp]E∥2E+∑EEh(∂K∩Γ)δhEh∥[(ufh−uph)⋅nf]E∥2E
3. The residual error estimator is locally defined by:

 ΘK=√Θ2K,f+Θ2K,p. (30)

The global residual error estimator is given by:

 Θ:=⎛⎝∑K∈ThΘ2K⎞⎠1/2. (31)

Furthermore denote the local and global approximation terms by

 ζK:={hK∥ff−ffh∥K,∀K∈Tfh,0∀K∈Tph,

and

 ζ:=⎛⎝∑K∈Thζ2K⎞⎠1/2. (32)
###### Remark 3.1.

The residual character of each term on the right-hand sides of (1) and (2) is quite clear since if would be the exact solution of (14), 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 [12] and [13] (see also [27]), 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 , and , we let and be the usual simplexe-bubble and face-bubble functions respectively (see (1.5) and (1.6) in [54]). In particular, satisfies , , , and . Similarly, , ,