    # An embedded--hybridized discontinuous Galerkin method for the coupled Stokes--Darcy system

We introduce an embedded--hybridized discontinuous Galerkin (EDG--HDG) method for the coupled Stokes--Darcy system. This EDG--HDG method is a pointwise mass-conserving discretization resulting in a divergence-conforming velocity field on the whole domain. In the proposed scheme, coupling between the Stokes and Darcy domains is achieved naturally through the EDG--HDG facet variables. A priori error analysis shows optimal convergence rates, and that the velocity error does not depend on the pressure. The error analysis is verified through numerical examples on unstructured grids for different orders of polynomial approximation.

## Authors

##### 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

Modelling adjacent free flow and porous media flow is important for a range of applications, e.g., transport of drugs via blood flow in vessels in biomedical engineering, and transport of pollutants via surface/groundwater flow in environmental engineering. The problem can be stated as a system of partial differential equations, with free flow governed by the Stokes equations and porous media flow governed by Darcy’s equations. The interactions at the boundary between the free flow and porous media flow regions were specified by

Beavers and Joseph (1967); Saffman (1971), and were mathematically justified in Mikelic and Jäger (2000). We refer to Discacciati and Quarteroni (2009) for an overview of the model.

Well-posedness of the weak formulation of the Stokes–Darcy problem can be found in Discacciati et al. (2002) for the primal form, and in Layton et al. (2002) for the primal–mixed form. Many different finite element and mixed finite element methods have been proposed to discretize the Stokes–Darcy problem for both formulations, e.g. Discacciati et al. (2002); Layton et al. (2002); Burman and Hansbo (2005); Cao et al. (2010); Gatica et al. (2009); Márquez et al. (2015). Other devised finite element methods include discontinuous Galerkin (DG) methods Çeşmelioğlu and Rivière (2009); Girault and Rivière (2009); Girault et al. (2014); Rivière (2005); Rivière and Yotov (2005), and hybridizable discontinuous Galerkin (HDG) methods Egger and Waluga (2013); Fu and Lehrenfeld (2018); Gatica and Sequeira (2017); Igreja and Loula (2018).

We develop a numerical scheme for which the velocity field is divergence-conforming on the whole domain and for which mass is conserved pointwise. Finite element methods that satisfy these properties were proposed in Girault et al. (2014); Kanschat and Rivière (2010) where they are referred to as ‘strongly conservative’. For the Stokes region Girault et al. (2014); Kanschat and Rivière (2010)

used a divergence-conforming DG space for the velocity and a standard DG space for the pressure. In the Darcy region they used a mixed finite element method. It is well-known, however, that DG methods can be expensive due to the large number of degrees-of-freedom on a given mesh compared to other methods. To reduce the number of globally coupled degrees-of-freedom,

Fu and Lehrenfeld (2018) proposed an HDG method for the Stokes region using a divergence-conforming finite element space for the velocity. Their method results in less globally coupled degrees-of-freedom compared to standard HDG methods as they only enforce continuity of the tangential direction of the facet velocity. Additionally, to reduce the problem size even further, they applied the ‘projected jumps’ method in which the polynomial degree of the tangential facet velocity is reduced by one compared to the cell velocity approximation (see also Lehrenfeld and Schöberl (2016)).

In this paper we propose an embedded–hybridized discontinuous Galerkin (EDG–HDG) finite element method of the primal–mixed formulation of the Stokes–Darcy problem on the whole domain. The EDG–HDG method uses a continuous trace velocity approximation and a discontinuous trace pressure approximation. Due to the continuous trace velocity approximation, the number of globally coupled degrees-of-freedom of the EDG–HDG method is fewer than for a traditional HDG method. However, the main motivation for an EDG–HDG discretization is not that the problem size is smaller, but that ‘continuous’ discretizations are generally better suited to fast iterative solvers. This was demonstrated for the Stokes problem in Rhebergen and Wells (2018), where CPU time and iteration count to convergence was reduced significantly compared to a hybridized method using only discontinuous facet approximations. We will show that the EDG–HDG method proposed in this paper is pointwise mass-conserving and that the resulting velocity field is divergence-conforming. We present furthermore an analysis of the proposed EDG–HDG method for the Stokes–Darcy problem, proving well-posedness, and optimal a priori

error estimates.

The remainder of this paper is organized as follows. In section 2 we briefly introduce the coupled Stokes–Darcy problem. The EDG–HDG method to this problem is presented in section 3. Consistency, continuity and well-posedness are shown in section 4 while the main results of this paper, an a priori error analysis, is presented in section 5. Numerical simulations support our theoretical results in section 6, and conclusions are drawn in section 7.

## 2 The Stokes–Darcy system

Let be a bounded polygonal domain with , boundary and boundary outward unit normal . We assume that is divided into two non-overlapping regions, and , such that and and are a union of polygonal subdomains. We denote the polygonal interface between and by . Furthermore, we define the external boundary of by , and the external boundary of by . See fig. 1.

Given the kinematic viscosity , forcing term , permeability constant and the source/sink term , the Stokes–Darcy system for the velocity field and pressure is given by

 −∇⋅2με(u)+∇p =fs in Ωs, (1a) κ−1u+∇p =0 in Ωd, (1b) −∇⋅u =χdfd in Ω, (1c) u =0 on Γs, (1d) u⋅n =0 on Γd, (1e)

where

is the strain rate tensor and

is the characteristic function that has the value 1 in

and 0 in . We will also frequently denote the velocity and pressure in by and , respectively, for .

Let

denote the unit normal vector on the interface between the two domains,

, pointing outwards from . On the interface we prescribe

 us⋅n =ud⋅n on ΓI, (2a) ps−2με(us)n⋅n =pd on ΓI, (2b) −2μ\delε(us)nt =ακ−1/2(us)t on ΓI, (2c)

where the tangential component of a vector is denoted by . Equation 2c is the Beavers–Joseph–Saffman law Beavers and Joseph (1967); Saffman (1971), where is an experimentally determined parameter.

## 3 The embedded–hybridized discontinuous Galerkin method

We present now an embedded–hybridized discontinuous Galerkin (EDG–HDG) method for the Stokes–Darcy system eqs. 2 and 1, and establish some of its key properties.

### 3.1 Preliminaries

For , let be a triangulation of into non-overlapping cells . For brevity, we consider the case of matching meshes at the interface . Furthermore, let . The diameter of a cell is denoted by and denotes the maximum diameter over all cells. The outward unit normal vector on the boundary of a cell, , is denoted by . An interior facet is shared by adjacent cells, and , while a boundary facet is part of that lies on . The set and union of all facets are denoted by and , respectively. By we denote the set of all facets that lie on . For , we denote by the set of all facets that lie in . Finally, for , we denote the union of all facets in by .

We consider the following discontinuous Galerkin finite element function spaces on ,

 Vh:=\cbrvh∈\sbrL2(Ω)dim: vh∈\sbrPk(K)dim, ∀ K∈T,Qh:=\cbrqh∈L2(Ω): qh∈Pk−1(K), ∀ K∈T∩L20(Ω), (3)

where denotes the space of polynomials of degree on domain . On and , we consider the finite element spaces:

 ¯Vh:=\cbr¯vh∈\sbrL2(Γs0)dim: ¯vh∈\sbrPk(F)dim ∀ F∈Fs, ¯vh=0 on Γs∩\sbrC0(Γs0)dim,¯Qjh:=\cbr¯qjh∈L2(Γj0): ¯qjh∈Pk(F) ∀ F∈Fj,j=s,d. (4)

Note that and do not necessarily coincide on the interface . Note also that functions in are continuous on , while functions in are discontinuous on , for .

For notational purposes, we introduce the spaces , and for . Function pairs in , and , for , will be denoted by , and . Furthermore, we set .

### 3.2 Method

The discrete form for the Stokes–Darcy system in eqs. 2 and 1 reads: given the forcing term , the source/sink term , the kinematic viscosity and the permeability , find such that

 Bh((uh,ph),(vh,qh))=∫Ωsfs⋅vh\difx+∫Ωdfdqh\difx∀(vh,qh)∈Xh, (5)

where

 Bh((uh,ph),(vh,qh)):=ah(uh,vh)+∑j=s,d\delbjh(pjh,vh)+bI,jh(¯pjh,¯vh)+∑j=s,d\delbjh(qjh,uh)+bI,jh(¯qjh,¯uh). (6)

The bilinear form is defined as

where

 ash(u,v) :=∑K∈Ts∫K2με(u):ε(v)\difx+∑K∈Ts∫∂K2βμhK(u−¯u)⋅(v−¯v)\difs (8a) −∑K∈Ts∫∂K2με(u)ns⋅(v−¯v)\difs−∑K∈Ts∫∂K2με(v)ns⋅(u−¯u)\difs, adh(u,v) :=∫Ωdκ−1u⋅v\difx, (8b) aIh(¯u,¯v) :=∫ΓIακ−1/2¯ut⋅¯vt\difs, (8c)

and where is a penalty parameter. The bilinear forms and are defined as:

 bjh(pj,v) :=−∑K∈Tj∫Kp∇⋅v\difx+∑K∈Tj∫∂K¯pjv⋅nj\difs, (9a) bI,jh(¯pj,¯v) :=−∫ΓI¯pj¯v⋅nj\difs. (9b)

### 3.3 Properties of the numerical scheme

Setting and in eq. 5 demonstrates cell-wise momentum balance eq. 1a subject to weak satisfaction of the boundary condition provided by , and a cell-wise statement of Darcy’s law eq. 1b subject to weak satisfaction of the boundary condition provided by and a Neumann boundary condition on . Setting and in eq. 5 shows that the formulation imposes normal continuity weakly across facets of the ‘numerical’ Stokes stress tensor:

 ^σsh:=2με(ush)−¯pshI−2βμhK(ush−¯ush)⊗n, (10)

where is the identity tensor. Setting and in eq. 5 and noting that , the numerical scheme imposes pointwise mass conservation, i.e.,

 −∇⋅uh=χdΠQfd∀x∈K, ∀K∈T, (11)

where is the standard -projection operator onto . Finally, setting , and with in eq. 5 and noting that on each , we find that is -conforming, i.e.,

 ⟦uh⋅n⟧ =0 ∀x∈F, ∀F∈F∖FI, (12a) uh⋅n =¯uh⋅n ∀x∈F, ∀F∈FI, (12b)

where is the usual jump operator and the unit normal vector on .

## 4 Consistency, continuity and well-posedness for Stokes–Darcy

### 4.1 Preliminaries

To prove consistency, continuity and stability we require extended function spaces and appropriate norms. We introduce

 V:=\cbrv:vs∈\sbrH2(Ωs)dim, vd∈\sbrH1(Ωd)dim, v=0 on Γs, v⋅n=0 on Γd, vs⋅n=vd⋅n on ΓI,Q:=\cbrq∈L20(Ω) : qs∈H1(Ωs), qd∈H2(Ωd), (13)

and . We let be the trace space of restricted to and be the trace space of restricted to . We introduce the trace operator to restrict functions in to , and the trace operators to restrict functions in to . We remark that on the interface . Where no ambiguity arises we omit the subscript when using the trace operator. For notational purposes we also introduce , and

 V(h):=Vh+V,Q(h):=Qh+Q,X(h):=V(h)×Q(h). (14)

For we denote by and the restriction of, respectively, and to .

We use various norms throughout, which are defined now. On we define

 |||v|||2v,s:=∑K∈Ts\del\norm∇v2K+h−1K\normv−¯v2∂K,

where denotes the standard -norm on domain , and

 |||v|||2v′,s:=|||v|||2v,s+∑K∈Tsh2K\envertv2H2(K).

On we introduce

 |||v|||2v:=|||v|||2v,s+\normv2Ωd+\norm¯vt2ΓI,

and

 |||v|||2v′:=|||v|||2v+∑K∈Tsh2K\envertv2H2(K)=|||v|||2v′,s+\normv2Ωd+\norm¯vt2ΓI.

Note that the norms and are equivalent on , see (Wells, 2011, eq. (5.5)).

Finally, for with and , we define

We will make use of various standard estimates. In particular, use will be made of the trace inequalities for , (Di Pietro and Ern, 2012, Lemma 1.46, Remark 1.47)

 \normv∂K≤CT,1h−1/2K\normvK∀v∈Pk(K), (15)

and the following straightforward extensions of the continuous trace inequality (Brenner and Scott, 2010, Theorem 1.6.6)

 \normv2∂K≤CT,2\delh−1K\normv2K+hK\normv21,K∀v∈H1(K), (16)

and

 \normvΓI≤Cc,T\norm∇vΩs∀v∈\cbrv∈H1(Ωs) : v=0 on Γs, (17)

where are independent of .

### 4.2 Consistency

We now prove that the scheme in eq. 5 is consistent with the Stokes–Darcy system in eqs. 2 and 1.

###### Lemma 1 (Consistency)

If solves the Stokes–Darcy system eqs. 2 and 1, then letting and ,

 Bh((u,p),(v,q))=∫Ωsfs⋅v\difx+∫Ωdfdq\difx∀(v,q)∈X(h). (18)

We consider each form in the definition of separately. Since on cell boundaries in , and integration by parts,

 ash(u,v)=∑K∈Ts∫K2με(u):ε(v)\difx−∑K∈Ts∫∂K2με(u)ns⋅(v−¯v)\difs=−∑K∈Ts∫K2μ(∇⋅ε(u))⋅v\difx+∑K∈Ts∫∂K2με(u)ns⋅¯v\difs. (19)

Using smoothness of , single-valuedness of , and eq. 2c, we note that

 ∑K∈Ts∫∂K2με(u)ns⋅¯v\difs=∫ΓI2με(u)ns⋅¯v\difs=−∫ΓIακ−1/2(us)t⋅¯vt\difs+∫ΓI\del2μ\delns⋅ε(u)nsns⋅¯v\difs. (20)

Combining with eq. 19,

 ash(u,v)=−∑K∈Ts∫K2μ(∇⋅ε(u))⋅v\difx−∫ΓIακ−1/2(us)t⋅¯vt\difs+∫ΓI\del2μ\delns⋅ε(u)nsns⋅¯v\difs. (21)

Applying integration by parts, noting that on cell boundaries in , with , and on , results in

 ∑j=s,d\delbjh(pj,v)+bI,jh(γ(pj),¯v)=∑K∈Ts∫K∇p⋅v\difx+∑K∈Td∫K∇p⋅v\difx+∫ΓI(pd−ps)¯v⋅ns\difs. (22)

Next,

 ∑j=s,d\delbjh(qj,u)+bI,jh(¯qj,γ(us))=−∑K∈Ts∫Kq∇⋅u\difx+∑K∈Ts∫∂K¯qsu⋅ns\difs−∫ΓI¯qsus⋅ns\difs−∑K∈Td∫Kq∇⋅u\difx+∑K∈Td∫∂K¯qdu⋅nd\difs−∫ΓI¯qdud⋅nd\difs=∫Ωdqfd\difx, (23)

where for the second equality we used eq. 1c-eq. 1e, and that , , and are single-valued on interior facets. Furthermore, note that

hence summing eqs. 24, 23, 22 and 21 results in

 Bh((u,p),(v,q))=∑K∈Ts∫K\del−2μ(∇⋅ε(u))+∇p⋅v\difx+∑K∈Td∫K\delκ−1u+∇p⋅v\difx+∫ΓI\del2μ\delns⋅ε(u)nsns⋅¯v\difs+∫ΓI(pd−ps)ns⋅¯v\difs+∫Ωdqfd\difx. (25)

The result follows after using eqs. 2b, 1b and 1a. ∎

### 4.3 Coercivity and continuity of ash(⋅,⋅) and ah(⋅,⋅)

In this section we show that is coercive on for sufficiently large penalty parameter . We furthermore prove continuity of and .

###### Lemma 2 (Coercivity)

There exists a constant , independent of , and a constant such that for and for all ,

 ah(vh,vh)≥C|||vh|||2v. (26)

Using identical steps as in (Rhebergen and Wells, 2017, Lemma 4.2) and applying Korn’s inequality Brenner (2004) it can be shown that

 ash(vh,vh)≥C∑K∈Ts\del∥ε(vh)∥2K+h−1K∥vh−¯vh∥2∂K≥C|||vh|||2v,s, (27)

The result follows by definition of . ∎

###### Lemma 3 (Continuity)

There exists a generic constant , independent of , such that for all ,

 ash(u,v) (28a) ah(u,v) ≤C|||u|||v′|||v|||v′. (28b)

Equation 28a follows by definition of eq. 8a, the Cauchy-Schwarz inequality, the trace inequality eq. 16, Hölder’s inequality for sums and since , . Equation 28b follows by definition of eq. 7, using the Cauchy–Schwarz inequality, eq. 28a and Hölder’s inequality for sums. ∎

### 4.4 The inf-sup condition

To present the inf-sup condition it will be convenient to split the velocity-pressure coupling term in eq. 9a as follows

 bjh(qjh,vh)=bj,1h(qh,vh)+bj,2h(¯qjh,vh), (29)

where

 bj,1h(qh,vh)=−∑K∈Tj∫Kqh∇⋅vh\difxandbj,2h(¯qjh,vh)=∑K∈Tj∫∂K¯qjhvh⋅nj\difs,

for .

The main result of this section is the following theorem.

###### Theorem 1 (inf-sup condition)

There exists a constant , independent of , such that for any ,

 (30)

To prove this result we require the definition of two interpolation operators. For the velocity we require the following BDM interpolation operator

(Hansbo and Larson, 2002, Lemma 7).

###### Lemma 4 (BDM interpolation operator)

If the mesh consists of triangles () or tetrahedra () there is an interpolation operator , where is the space of functions in restricted to cells in , such that for all ,

1. for all .

2. for all , where is an edge () or face () of .

3. , where is the usual jump operator.

4. with and .

We also require an interpolation operator , for example the Scott–Zhang interpolant (Brenner and Scott, 2010, Theorem 4.8.12), with the property that for , ,

 \normv−¯ΠVv∂K ≤Chℓ−1/2K\normvHℓ(K), (31a) \normΠVv−¯ΠVv∂K ≤Chℓ−1/2K\normvHℓ(K), (31b)

where is a generic constant independent of .

Additionally, we require the following two auxiliary results.

###### Lemma 5

There exists a constant , independent of , such that for any ,

 cinf\normqhΩ≤supvh∈Vhvh≠0∑j=s,dbj,1h(qh,vh)|||vh|||v. (32)

Let . It is well known, e.g., (Di Pietro and Ern, 2012, Theorem 6.5), that since there exists a such that

 −∇⋅v=qhandcinf\normvH1(Ω)≤\normqhΩ. (33)

Then, by lemma 4 (1), since for all ,

 \normqh2Ω=−∫Ωqh∇⋅v\difx=−∫Ωqh∇⋅ΠVv\difx. (34)

By lemma 4 (4) and eq. 31b,

 |||(ΠVv,¯ΠVv)|||2v,s=∑K∈Ts\del\norm∇(ΠVv)2K+h−1K\normΠVv−¯