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

TL;DR

This paper presents an embedded-hybridized discontinuous Galerkin method for the coupled Stokes-Darcy system, achieving optimal convergence and divergence-conforming velocities, verified through numerical experiments.

Contribution

The paper introduces a novel EDG-HDG discretization that ensures mass conservation and divergence conformity for coupled Stokes-Darcy problems.

Findings

Achieves optimal convergence rates.

Velocity error independent of pressure.

Validated on unstructured grids with various polynomial orders.

Abstract

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. \emph{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.

Tables

Table 1. Table 1: Errors and rates of convergence in Ω s superscript Ω 𝑠 \Omega^{s} (top half) and Ω d superscript Ω 𝑑 \Omega^{d} (bottom half) for the velocity and pressure fields using polynomial degree k = 3 𝑘 3 k=3 for the test case described in section   6.1 .

Stokes solution (${\Omega }^s$)
Cells	${\parallel u_h-u\parallel }_{{\Omega }^s}$	Rate	${\parallel p_h-p\parallel }_{{\Omega }^s}$	Rate	${\parallel \nabla \cdot u_h\parallel }_{{\Omega }^s}$
$\mu =1$, $\kappa =1$
152	1.7e-6	4.8	2.9e-4	3.6	2.3e-14
578	8.4e-8	4.3	3.6e-5	3.0	4.8e-14
2416	4.4e-9	4.3	4.5e-6	3.0	9.3e-14
9580	2.3e-10	4.3	5.6e-7	3.0	1.9e-13
$\mu ={10}^{-6}$, $\kappa =1$
152	9.0e-6	7.6	5.1e-5	3.2	7.5e-12
578	1.2e-7	6.3	6.1e-6	3.1	7.0e-12
2416	4.5e-9	4.7	7.4e-7	3.0	7.0e-12
9580	2.3e-10	4.3	9.0e-8	3.0	6.9e-12
$\mu =1$, $\kappa ={10}^3$
152	1.7e-6	4.8	2.9e-4	3.6	2.3e-14
578	8.4e-8	4.3	3.6e-5	3.0	4.8e-14
2416	4.4e-9	4.3	4.5e-6	3.0	9.3e-14
9580	2.3e-10	4.3	5.6e-7	3.0	1.9e-13
$\mu ={10}^{-6}$, $\kappa ={10}^3$
152	1.6e-6	4.8	5.1e-8	3.2	1.3e-14
578	8.3e-8	4.3	6.1e-9	3.1	1.4e-14
2416	4.4e-9	4.3	7.4e-10	3.0	5.6e-14
9580	2.3e-10	4.3	9.0e-11	3.0	4.9e-14
Darcy solution (${\Omega }^d$)
Cells	${\parallel u_h-u\parallel }_{{\Omega }^d}$	Rate	${\parallel p_h-p\parallel }_{{\Omega }^d}$	Rate	${\parallel \nabla \cdot u_h+{\Pi }_Q{f}^d\parallel }_{{\Omega }^d}$
$\mu =1$, $\kappa =1$
152	3.1e-6	4.7	3.5e-5	3.4	1.2e-12
578	2.0e-7	4.0	4.7e-6	2.9	3.7e-12
2416	1.2e-8	4.1	5.8e-7	3.0	1.5e-11
9580	7.5e-10	4.0	7.3e-8	3.0	6.0e-11
$\mu ={10}^{-6}$, $\kappa =1$
152	3.0e-6	4.7	3.5e-5	3.4	9.4e-13
578	2.0e-7	3.9	4.7e-6	2.9	3.6e-12
2416	1.2e-8	4.1	5.8e-7	3.0	1.5e-11
9580	7.5e-10	4.0	7.3e-8	3.0	5.9e-11
$\mu =1$, $\kappa ={10}^3$
152	3.1e-6	4.7	3.5e-8	3.4	2.1e-12
578	2.0e-7	4.0	4.7e-9	2.9	3.6e-12
2416	1.2e-8	4.1	5.8e-10	3.0	1.0e-10
9580	7.5e-10	4.0	7.3e-11	3.0	6.2e-11
$\mu ={10}^{-6}$, $\kappa ={10}^3$
152	3.1e-6	4.7	3.5e-8	3.4	1.9e-12
578	2.0e-7	3.9	4.7e-9	2.9	3.6e-12
2416	1.2e-8	4.1	5.8e-10	3.0	1.4e-10
9580	7.5e-10	4.0	7.3e-11	3.0	6.3e-11

Equations

$$-\nabla \cdot 2\mu \epsilon (u)+\nabla p$$

$${\kappa }^{-1}u+\nabla p$$

$$-\nabla \cdot u$$

$$u$$

$$u\cdot n$$

$$u^s\cdot n$$

$$p^s-2\mu \epsilon (u^s)n\cdot n$$

$$-2\mu {\left(\epsilon (u^s)n\right)}^t$$

$$\begin{split}V_{h}&\mathrel{\mathop{\mathchar 58\relax}}=\mathinner{\bigl{\{}v_{h}\in\mathinner{[L^{2}(\Omega)]}^{\rm dim}\mathrel{\mathop{\mathchar 58\relax}}\ v_{h}\in\mathinner{[P_{k}(K)]}^{\rm dim}\ \forall\ K\in\mathcal{T}\bigr{\}}},\\ Q_{h}&\mathrel{\mathop{\mathchar 58\relax}}=\mathinner{\bigl{\{}q_{h}\in L^{2}(\Omega)\mathrel{\mathop{\mathchar 58\relax}}\ q_{h}\in P_{k-1}(K)\ \forall\ K\in\mathcal{T}\bigr{\}}}\cap L^{2}_{0}(\Omega),\\ Q_{h}^{j}&\mathrel{\mathop{\mathchar 58\relax}}=\mathinner{\bigl{\{}q_{h}\in L^{2}(\Omega^{j})\mathrel{\mathop{\mathchar 58\relax}}\ q_{h}\in P_{k-1}(K)\ \forall\ K\in\mathcal{T}^{j}\bigr{\}}},\quad j=s,d,\end{split}$$

$$\begin{split}\bar{V}_{h}&\mathrel{\mathop{\mathchar 58\relax}}=\mathinner{\bigl{\{}\bar{v}_{h}\in\mathinner{[L^{2}(\Gamma_{0}^{s})]}^{\rm dim}\mathrel{\mathop{\mathchar 58\relax}}\ \bar{v}_{h}\in\mathinner{[P_{k}(F)]}^{\rm dim}\ \forall\ F\in\mathcal{F}^{s},\ \bar{v}_{h}=0\ \mbox{on}\ \Gamma^{s}\bigr{\}}}\cap\mathinner{\bigl{[}C^{0}(\Gamma_{0}^{s})\bigr{]}}^{\rm dim},\\ \bar{Q}_{h}^{j}&\mathrel{\mathop{\mathchar 58\relax}}=\mathinner{\bigl{\{}\bar{q}_{h}^{j}\in L^{2}(\Gamma_{0}^{j})\mathrel{\mathop{\mathchar 58\relax}}\ \bar{q}_{h}^{j}\in P_{k}(F)\ \forall\ F\in\mathcal{F}^{j}\bigr{\}}},\quad j=s,d.\end{split}$$

$$B_h(({\mathbf{u}}_h,{\mathbf{p}}_h),({\mathbf{v}}_h,{\mathbf{q}}_h))={\int }_{{\Omega }^s}{f}^s\cdot v_{h}dx+{\int }_{{\Omega }^d}{f}^d{q}_{h}dx\forall ({\mathbf{v}}_h,{\mathbf{q}}_h)\in {\mathbf{X}}_h,$$

$$\begin{gathered} B_h(({\mathbf{u}}_h,{\mathbf{p}}_h),({\mathbf{v}}_h,{\mathbf{q}}_h))\text{mathchar}58=a_h({\mathbf{u}}_h,{\mathbf{v}}_h)+\sum _{j=s,d}\left(b_h^j({\mathbf{p}}_h^j,v_h)+b_h^{I,j}({\bar{p}}_h^j,{\bar{v}}_h)\right) \\ +\sum _{j=s,d}\left(b_h^j({\mathbf{q}}_h^j,u_h)+b_h^{I,j}({\bar{q}}_h^j,{\bar{u}}_h)\right). \end{gathered}$$

$$a_h({\mathbf{u}}_h,{\mathbf{v}}_h)\text{mathchar}58=a_h^s({\mathbf{u}}_h,{\mathbf{v}}_h)+a_h^d(u_h,v_h)+a_h^I({\bar{u}}_h,{\bar{v}}_h),$$

$$a_h^s(\mathbf{u},\mathbf{v})$$

$$-\sum _{K\in {\mathcal{T}}^s}{\int }_{\partial K}2\mu \epsilon (u)n^s\cdot (v-\bar{v})ds-\sum _{K\in {\mathcal{T}}^s}{\int }_{\partial K}2\mu \epsilon (v)n^s\cdot (u-\bar{u})ds,$$

$$a_h^d(u,v)$$

$$a_h^I(\bar{u},\bar{v})$$

$$b_h^j({\mathbf{p}}^j,v)$$

$$b_h^{I,j}({\bar{p}}^j,\bar{v})$$

$${\hat{\sigma }}_h^s\text{mathchar}58=2\mu \epsilon (u_h^s)-{\bar{p}}_h^s\mathbb{I}-\frac{2\beta \mu }{h_K}(u_h^s-{\bar{u}}_h^s)\otimes n,$$

$$-\nabla \cdot u_h={\chi }^d{\Pi }_Q{f}^d\forall x\in K,\forall K\in \mathcal{T},$$

$$[[u_h\cdot n]]$$

$$u_h\cdot n$$

$$\begin{split}V&\mathrel{\mathop{\mathchar 58\relax}}=\mathinner{\bigl{\{}v\mathrel{\mathop{\mathchar 58\relax}}v^{s}\in\mathinner{[H^{2}(\Omega^{s})]}^{\rm dim},\ v^{d}\in\mathinner{[H^{1}(\Omega^{d})]}^{\rm dim},\ v=0\ \text{on}\ \Gamma^{s},\ v\cdot n=0\ \text{on}\ \Gamma^{d},\ v^{s}\cdot n=v^{d}\cdot n\ \text{on}\ \Gamma^{I}\bigr{\}}},\\ Q&\mathrel{\mathop{\mathchar 58\relax}}=\mathinner{\bigl{\{}q\in L_{0}^{2}(\Omega)\ \mathrel{\mathop{\mathchar 58\relax}}\ q^{s}\in H^{1}(\Omega^{s}),\ q^{d}\in H^{2}(\Omega^{d})\bigr{\}}},\end{split}$$

$$\mathbf{V}(h)\text{mathchar}58={\mathbf{V}}_h+\mathbf{V},\mathbf{Q}(h)\text{mathchar}58={\mathbf{Q}}_h+\mathbf{Q},\mathbf{X}(h)\text{mathchar}58=\mathbf{V}(h)\times \mathbf{Q}(h).$$

$$\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{v}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}_{v,s}^{2}\mathrel{\mathop{\mathchar 58\relax}}=\sum_{K\in\mathcal{T}^{s}}\mathinner{\bigl{(}\mathinner{\lVert\nabla v\rVert}_{K}^{2}+h_{K}^{-1}\mathinner{\!\left\lVert v-\bar{v}\right\rVert}^{2}_{\partial K}\bigr{)}},$$

$$|||\mathbf{v}{|||}_{v^{\prime },s}^2\text{mathchar}58=|||\mathbf{v}{|||}_{v,s}^2+\sum _{K\in {\mathcal{T}}^s}{h}_K^2{\left|v\right|}_{H^2(K)}^2.$$

$$|||\mathbf{v}{|||}_v^2\text{mathchar}58=|||\mathbf{v}{|||}_{v,s}^2+{\Vert v\Vert }_{{\Omega }^d}^2+{\Vert {\bar{v}}^t\Vert }_{{\Gamma }^I}^2,$$

$$|||\mathbf{v}{|||}_{v^{\prime }}^2\text{mathchar}58=|||\mathbf{v}{|||}_v^2+\sum _{K\in {\mathcal{T}}^s}{h}_K^2{\left|v\right|}_{H^2(K)}^2=|||\mathbf{v}{|||}_{v^{\prime },s}^2+{\Vert v\Vert }_{{\Omega }^d}^2+{\Vert {\bar{v}}^t\Vert }_{{\Gamma }^I}^2.$$

$$|||{\mathbf{q}}^j{|||}_{p,j}^2\text{mathchar}58={\Vert q\Vert }_{{\Omega }^j}^2+\sum _{K\in {\mathcal{T}}^j}{h}_K{\Vert {\bar{q}}^j\Vert }_{\partial K}^2,|||\mathbf{q}{|||}_p^2\text{mathchar}58=\sum _{j=s,d}|||{\mathbf{q}}^j{|||}_{p,j}^2.$$

$${\Vert v\Vert }_{\partial K}\le C_{T,1}{h}_K^{-1/2}{\Vert v\Vert }_K\forall v\in P_k(K),$$

$$\mathinner{\!\left\lVert v\right\rVert}_{\partial K}^{2}\leq C_{T,2}\mathinner{\bigl{(}h_{K}^{-1}\mathinner{\!\left\lVert v\right\rVert}_{K}^{2}+h_{K}\mathinner{\!\left\lVert v\right\rVert}_{1,K}^{2}\bigr{)}}\quad\forall v\in H^{1}(K),$$

$$\mathinner{\!\left\lVert v\right\rVert}_{\Gamma^{I}}\leq C_{c,T}\mathinner{\lVert\nabla v\rVert}_{\Omega^{s}}\quad\forall v\in\mathinner{\bigl{\{}v\in H^{1}(\Omega^{s})\ \mathrel{\mathop{\mathchar 58\relax}}\ v=0\ \text{on}\ \Gamma^{s}\bigr{\}}},$$

$$B_h((\mathbf{u},\mathbf{p}),(\mathbf{v},\mathbf{q}))={\int }_{{\Omega }^s}{f}^s\cdot vdx+{\int }_{{\Omega }^d}{f}^{d}qdx\forall (\mathbf{v},\mathbf{q})\in \mathbf{X}(h).$$

$$\begin{array}{rl}{a}_h^s(\mathbf{u},\mathbf{v})& =\sum _{K\in {\mathcal{T}}^s}{\int }_{K}2\mu \epsilon (u)\text{mathchar}58\epsilon (v)dx-\sum _{K\in {\mathcal{T}}^s}{\int }_{\partial K}2\mu \epsilon (u)n^s\cdot (v-\bar{v})ds\\ & =-\sum _{K\in {\mathcal{T}}^s}{\int }_{K}2\mu (\nabla \cdot \epsilon (u))\cdot vdx+\sum _{K\in {\mathcal{T}}^s}{\int }_{\partial K}2\mu \epsilon (u)n^s\cdot \bar{v}ds.\end{array}$$

$$\begin{array}{rl}\sum _{K\in {\mathcal{T}}^s}{\int }_{\partial K}2\mu \epsilon (u)n^s\cdot \bar{v}ds& ={\int }_{{\Gamma }^I}2\mu \epsilon (u)n^s\cdot \bar{v}ds\\ & =-{\int }_{{\Gamma }^I}\alpha {\kappa }^{-1/2}(u^s{)}^t\cdot {\bar{v}}^{t}ds+{\int }_{{\Gamma }^I}\left(2\mu \left(n^s\cdot \epsilon (u)n^s\right)n^s\right)\cdot \bar{v}ds.\end{array}$$

Full text

††journal: Journal of Computational and Applied Mathematics

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 [1, 2], and were mathematically justified in [3]. We refer to [4] for an overview of the model.

Well-posedness of the weak formulation of the Stokes–Darcy problem can be found in [5] for the primal form, and in [6] 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. [5, 6, 7, 8, 9, 10, 11]. Other devised finite element methods include discontinuous Galerkin (DG) methods [12, 13, 14, 15, 16, 17], hybridizable discontinuous Galerkin (HDG) methods [18, 19, 20, 21], weak Galerkin methods (WG) [22], and weak virtual element methods (WVEM) [23].

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 [14, 24] where they are referred to as ‘strongly conservative’. For the Stokes region [14, 24] 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, [19] 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 [25]).

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 [26], 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 $\Omega\subset\mathbb{R}^{\rm dim}$ be a bounded polygonal domain with ${\rm dim}=2,3$, boundary $\partial\Omega$ and boundary outward unit normal $n$. We assume that $\Omega$ is divided into two non-overlapping regions, $\Omega^{s}$ and $\Omega^{d}$, such that $\Omega=\Omega^{s}\cup\Omega^{d}$ and $\Omega^{s}$ and $\Omega^{d}$ are a union of polygonal subdomains. We denote the polygonal interface between $\Omega^{s}$ and $\Omega^{d}$ by $\Gamma^{I}$. Furthermore, we define the external boundary of $\Omega^{s}$ by $\Gamma^{s}\mathrel{\mathop{\mathchar 58\relax}}=\partial\Omega\cap\partial\Omega^{s}$, and the external boundary of $\Omega^{d}$ by $\Gamma^{d}\mathrel{\mathop{\mathchar 58\relax}}=\partial\Omega\cap\partial\Omega^{d}$. See fig. 1.

Given the kinematic viscosity $\mu\in\mathbb{R}^{+}$, forcing term $f^{s}\mathrel{\mathop{\mathchar 58\relax}}\Omega^{s}\to\mathbb{R}^{\rm dim}$, permeability constant $\kappa\in\mathbb{R}^{+}$ and the source/sink term $f^{d}\mathrel{\mathop{\mathchar 58\relax}}\Omega^{d}\to\mathbb{R}$, the Stokes–Darcy system for the velocity field $u\mathrel{\mathop{\mathchar 58\relax}}\Omega\to\mathbb{R}^{\rm dim}$ and pressure $p\mathrel{\mathop{\mathchar 58\relax}}\Omega\to\mathbb{R}$ is given by

[TABLE]

where $\varepsilon(u)\mathrel{\mathop{\mathchar 58\relax}}=\mathinner{\left(\nabla u+(\nabla u)^{T}\right)}/2$ is the strain rate tensor and $\chi^{d}$ is the characteristic function that has the value 1 in $\Omega^{d}$ and 0 in $\Omega^{s}$. We will also frequently denote the velocity and pressure in $\Omega^{j}$ by $u^{j}$ and $p^{j}$, respectively, for $j=s,d$.

Let $n$ denote the unit normal vector on the interface between the two domains, $\Gamma^{I}$, pointing outwards from $\Omega^{s}$. On the interface $\Gamma^{I}$ we prescribe

[TABLE]

where the tangential component of a vector $w$ is denoted by $(w)^{t}\mathrel{\mathop{\mathchar 58\relax}}=w-(w\cdot n)n$. Equation 2c is the Beavers–Joseph–Saffman law [1, 2], where $\alpha>0$ 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. 1 and 2, and establish some of its key properties.

3.1 Preliminaries

For $j=s,d$, let $\mathcal{T}^{j}\mathrel{\mathop{\mathchar 58\relax}}=\mathinner{\left\{K\right\}}$ be a triangulation of $\Omega^{j}$ into non-overlapping cells $K$. For brevity, we consider the case of matching meshes at the interface $\Gamma^{I}$. Furthermore, let $\mathcal{T}\mathrel{\mathop{\mathchar 58\relax}}=\mathcal{T}^{s}\cup\mathcal{T}^{d}$. The diameter of a cell $K$ is denoted by $h_{K}$ and $h$ denotes the maximum diameter over all cells. The outward unit normal vector on the boundary of a cell, $\partial K$, is denoted by $n$. An interior facet $F$ is shared by adjacent cells, $K^{+}$ and $K^{-}$, while a boundary facet is part of $\partial K$ that lies on $\partial\Omega$. The set and union of all facets are denoted by $\mathcal{F}\mathrel{\mathop{\mathchar 58\relax}}=\mathinner{\left\{F\right\}}$ and $\Gamma_{0}$, respectively. By $\mathcal{F}^{I}$ we denote the set of all facets that lie on $\Gamma^{I}$. For $j=s,d$, we denote by $\mathcal{F}^{j}$ the set of all facets that lie in $\overline{\Omega}^{j}$. Finally, for $j=s,d$, we denote the union of all facets in $\overline{\Omega}^{j}$ by $\Gamma_{0}^{j}$.

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

[TABLE]

where $P_{k}(D)$ denotes the space of polynomials of degree $k$ on domain $D$ and $L^{2}_{0}(\Omega)\mathrel{\mathop{\mathchar 58\relax}}=\{q\in L^{2}(\Omega)\mathrel{\mathop{\mathchar 58\relax}}\int_{\Omega}q\operatorname{d\!}x=0\}$. On $\Gamma_{0}^{s}$ and $\Gamma_{0}^{d}$, we consider the finite element spaces:

[TABLE]

Note that $\bar{q}_{h}^{s}\in\bar{Q}_{h}^{s}$ and $\bar{q}_{h}^{d}\in\bar{Q}_{h}^{d}$ do not necessarily coincide on the interface $\Gamma^{I}$. Note also that functions in $\bar{V}_{h}$ are continuous on $\Gamma_{0}^{s}$, while functions in $\bar{Q}_{h}^{j}$ are discontinuous on $\Gamma_{0}^{j}$, for $j=s,d$.

For notational purposes, we introduce the spaces $\boldsymbol{V}_{h}\mathrel{\mathop{\mathchar 58\relax}}=V_{h}\times\bar{V}_{h}$, $\boldsymbol{Q}_{h}\mathrel{\mathop{\mathchar 58\relax}}=Q_{h}\times\bar{Q}_{h}^{s}\times\bar{Q}_{h}^{d}$ and $\boldsymbol{Q}_{h}^{j}\mathrel{\mathop{\mathchar 58\relax}}=Q_{h}^{j}\times\bar{Q}_{h}^{j}$ for $j=s,d$. Function pairs in $\boldsymbol{V}_{h}$, $\boldsymbol{Q}_{h}$ and $\boldsymbol{Q}_{h}^{j}$, for $j=s,d$, will be denoted by $\boldsymbol{v}_{h}\mathrel{\mathop{\mathchar 58\relax}}=(v_{h},\bar{v}_{h})\in\boldsymbol{V}_{h}$, $\boldsymbol{q}_{h}\mathrel{\mathop{\mathchar 58\relax}}=(q_{h},\bar{q}_{h}^{s},\bar{q}_{h}^{d})\in\boldsymbol{Q}_{h}$ and $\boldsymbol{q}_{h}^{j}\mathrel{\mathop{\mathchar 58\relax}}=(q_{h},\bar{q}_{h}^{j})\in\boldsymbol{Q}_{h}^{j}$. Furthermore, we set $\boldsymbol{X}_{h}\mathrel{\mathop{\mathchar 58\relax}}=\boldsymbol{V}_{h}\times\boldsymbol{Q}_{h}$.

3.2 Method

The discrete form for the Stokes–Darcy system in eqs. 1 and 2 reads: given the forcing term $f^{s}\in\mathinner{[L^{2}(\Omega^{s})]}^{\rm dim}$, the source/sink term $f^{d}\in L^{2}(\Omega^{d})$, the kinematic viscosity $\mu\in\mathbb{R}^{+}$ and the permeability $\kappa\in\mathbb{R}^{+}$, find $(\boldsymbol{u}_{h},\boldsymbol{p}_{h})\in\boldsymbol{X}_{h}$ such that

[TABLE]

where

[TABLE]

The bilinear form $a_{h}(\cdot,\cdot)$ is defined as

[TABLE]

where

[TABLE]

and where $\beta>0$ is a penalty parameter. The bilinear forms $b_{h}^{j}(\cdot,\cdot)$ and $b_{h}^{I,j}(\cdot,\cdot)$ are defined as:

[TABLE]

3.3 Properties of the numerical scheme

Setting $\bar{v}_{h}=0$ and $\boldsymbol{q}_{h}=0$ in eq. 5 demonstrates cell-wise momentum balance eq. 1a subject to weak satisfaction of the boundary condition provided by $\bar{u}_{h}$, and a cell-wise statement of Darcy’s law eq. 1b subject to weak satisfaction of the boundary condition provided by $\bar{p}_{h}^{d}$ and a Neumann boundary condition on $\Gamma^{I}$. Setting $v_{h}=0$ and $\boldsymbol{q}_{h}=0$ in eq. 5 shows that the formulation imposes normal continuity weakly across facets of the ‘numerical’ Stokes stress tensor:

[TABLE]

where $\mathbb{I}$ is the identity tensor. Setting $\boldsymbol{v}_{h}=0$ and $\bar{q}_{h}^{s}=\bar{q}_{h}^{d}=0$ in eq. 5 and noting that $\nabla\cdot u_{h}\in P_{k-1}(K)$, the numerical scheme imposes pointwise mass conservation, i.e.,

[TABLE]

where $\Pi_{Q}$ is the standard $L^{2}$-projection operator onto $Q_{h}$. Finally, setting $\boldsymbol{v}_{h}=0$, $q_{h}=0$ and $\bar{q}_{h}^{l}=0$ with $l\neq j$ in eq. 5 and noting that $\llbracket u_{h}\cdot n\rrbracket,u_{h}\cdot n,(u_{h}-\bar{u}_{h})\cdot n\in P_{k}(F)$ on each $F\in\mathcal{F}$, we find that $u_{h}$ is $H(\text{div})$-conforming, i.e.,

[TABLE]

where $\llbracket\cdot\rrbracket$ is the usual jump operator and $n$ the unit normal vector on $F$.

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

[TABLE]

and $X\mathrel{\mathop{\mathchar 58\relax}}=V\times Q$. We let $\bar{V}$ be the trace space of $V$ restricted to $\Gamma_{0}^{s}$ and $\bar{Q}$ be the trace space of $Q$ restricted to $\Gamma_{0}$. We introduce the trace operator $\gamma_{V}\mathrel{\mathop{\mathchar 58\relax}}V\to\bar{V}$ to restrict functions in $V$ to $\Gamma_{0}^{s}$, and the trace operators $\gamma_{Q^{j}}\mathrel{\mathop{\mathchar 58\relax}}Q^{j}\to\bar{Q}^{j}$ to restrict functions in $Q^{j}$ to $\Gamma_{0}^{j}$. We remark that $\gamma_{Q^{s}}(p^{s})\neq\gamma_{Q^{d}}(p^{d})$ on the interface $\Gamma^{I}$. Where no ambiguity arises we omit the subscript when using the trace operator. For notational purposes we also introduce $\boldsymbol{V}\mathrel{\mathop{\mathchar 58\relax}}=V\times\bar{V}$, $\boldsymbol{Q}\mathrel{\mathop{\mathchar 58\relax}}=Q\times\bar{Q}$ and

[TABLE]

For $j=s,d$ we denote by $\boldsymbol{V}^{j}(h)$ and $\boldsymbol{Q}^{j}(h)$ the restriction of, respectively, $\boldsymbol{V}(h)$ and $\boldsymbol{Q}(h)$ to $\Omega^{j}$.

We use various norms throughout, which are defined now. On $\boldsymbol{V}^{s}(h)$ we define

[TABLE]

where $\mathinner{\!\left\lVert\cdot\right\rVert}_{D}$ denotes the standard $L^{2}$-norm on domain $D$, and

[TABLE]

On $\boldsymbol{V}(h)$ we introduce

[TABLE]

and

[TABLE]

Note that the norms $\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\cdot\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}_{v}$ and $\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\cdot\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}_{v^{\prime}}$ are equivalent on $\boldsymbol{V}_{h}$, see [27, eq. (5.5)].

Finally, for $\boldsymbol{q}^{j}\in\boldsymbol{Q}^{j}(h)$ with $j=s,d$ and $\boldsymbol{q}\in\boldsymbol{Q}(h)$, we define

[TABLE]

We will make use of various standard estimates. In particular, use will be made of the trace inequalities for $K\in\mathcal{T}$, [28, Lemma 1.46, Remark 1.47]

[TABLE]

and the following straightforward extensions of the continuous trace inequality [29, Theorem 1.6.6]

[TABLE]

and

[TABLE]

where $C_{T,1},C_{T,2},C_{c,T}>0$ are independent of $h_{K}$.

4.2 Consistency

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

Lemma 1** (Consistency)**

If $(u,p)\in X$ solves the Stokes–Darcy system eqs. 1 and 2, then letting $\boldsymbol{u}=(u,\gamma(u))$ and $\boldsymbol{p}=(p,\gamma(p^{s}),\gamma(p^{d}))$,

[TABLE]

Proof 1

We consider each form in the definition of $B_{h}$ separately. Since $u=\gamma(u)$ on cell boundaries in $\overline{\Omega}^{s}$, and integration by parts,

[TABLE]

Using smoothness of $u$, single-valuedness of $\bar{v}$, and eq. 2c, we note that

[TABLE]

Combining with eq. 19,

[TABLE]

Applying integration by parts, noting that $\gamma(p^{j})=p^{j}$ on cell boundaries in $\overline{\Omega}^{j}$, with $j=s,d$, and $n^{d}=-n^{s}$ on $\Gamma^{I}$, results in

[TABLE]

Next,

[TABLE]

where for the second equality we used eq. 1c-eq. 1e, and that $\bar{q}^{s}$, $\bar{q}^{d}$, $u^{s}\cdot n^{s}$ and $u^{d}\cdot n^{d}$ are single-valued on interior facets. Furthermore, note that

[TABLE]

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

[TABLE]

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

4.3 Coercivity and continuity of $a_{h}^{s}(\cdot,\cdot)$ and $a_{h}(\cdot,\cdot)$

In this section we show that $a_{h}(\cdot,\cdot)$ is coercive on $\boldsymbol{V}_{h}$ for sufficiently large penalty parameter $\beta$. We furthermore prove continuity of $a_{h}^{s}(\cdot,\cdot)$ and $a_{h}(\cdot,\cdot)$.

Lemma 2** (Coercivity)**

There exists a constant $C>0$, independent of $h$, and a constant $\beta_{0}>0$ such that for $\beta>\beta_{0}$ and for all $\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}$,

[TABLE]

Proof 2

Using identical steps as in [30, Lemma 4.2] and applying Korn’s inequality [31] it can be shown that

[TABLE]

The result follows by definition of $a_{h}$. \qed

Lemma 3** (Continuity)**

There exists a generic constant $C>0$, independent of $h$, such that for all $\boldsymbol{u},\boldsymbol{v}\in\boldsymbol{V}(h)$,

[TABLE]

Proof 3

Equation 28a* follows by definition of $a_{h}^{s}$ eq. 8a, the Cauchy-Schwarz inequality, the trace inequality eq. 16, Hölder’s inequality for sums and since $\mathinner{\lvert\varepsilon(u)\rvert}_{H^{\ell}(K)}\leq C\mathinner{\!\left\lvert u\right\rvert}_{H^{\ell+1}(K)}$, $\ell=0,1$. Equation 28b follows by definition of $a_{h}$ eq. 7, using the Cauchy–Schwarz inequality, eq. 28a and Hölder’s inequality for sums. \qed*

4.4 The inf-sup condition

To present the inf-sup condition it will be convenient to split the velocity-pressure coupling term $b_{h}^{j}(\boldsymbol{q}_{h}^{j},v_{h})$ in eq. 9a as follows

[TABLE]

where

[TABLE]

for $j=s,d$.

The main result of this section is the following theorem.

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

There exists a constant $c_{{\rm inf}}^{\star}>0$, independent of $h$, such that for any $\boldsymbol{q}_{h}\in\boldsymbol{Q}_{h}$,

[TABLE]

To prove this result we require the definition of two interpolation operators. For the velocity we require the following BDM interpolation operator [32, Lemma 7].

Lemma 4** (BDM interpolation operator)**

If the mesh consists of triangles (${\rm dim}=2$) or tetrahedra (${\rm dim}=3$) there is an interpolation operator $\Pi_{V}\mathrel{\mathop{\mathchar 58\relax}}\mathinner{[H^{1}(\Omega^{j})]}^{\rm dim}\rightarrow V_{h}^{j}$, where $V_{h}^{j}$ is the space of functions in $V_{h}$ restricted to cells in $\mathcal{T}^{j}$, such that for all $u\in\mathinner{[H^{k+1}(K)]}^{\rm dim}$,

1. i.

$\displaystyle\int_{K}q(\nabla\cdot(u-\Pi_{V}u))\operatorname{d\!}x=0$* for all $q\in P_{k-1}(K)$.* 2. ii.

$\displaystyle\int_{F}\bar{q}(u-\Pi_{V}u)\cdot n\operatorname{d\!}s=0$* for all $\bar{q}\in P_{k}(F)$, where $F$ is an edge (${\rm dim}=2$) or face (${\rm dim}=3$) of $\partial K$.* 3. iii.

$\llbracket n\cdot\Pi_{V}u\rrbracket=0$, where $\llbracket\cdot\rrbracket$ is the usual jump operator. 4. iv.

$\mathinner{\!\left\lVert u-\Pi_{V}u\right\rVert}_{m,K}\leq Ch_{K}^{\ell-m}\mathinner{\!\left\lVert u\right\rVert}_{H^{\ell}(K)}$* with $m=0,1,2$ and $m\leq\ell\leq k+1$.*

We also require an interpolation operator $\bar{\Pi}_{V}\mathrel{\mathop{\mathchar 58\relax}}\mathinner{[H^{1}(\Omega)]}^{\rm dim}\to V_{h}\cap\mathinner{[C^{0}(\bar{\Omega})]}^{\rm dim}$, for example the Scott–Zhang interpolant [29, Theorem 4.8.12], with the property that for $v\in\mathinner{[H^{\ell}(K)]}^{\rm dim}$, $1\leq\ell\leq k+1$,

[TABLE]

where $C$ is a generic constant independent of $h$.

Additionally, we require the following two auxiliary results.

Lemma 5

There exists a constant $c_{{\rm inf}}>0$, independent of $h$, such that for any $q_{h}\in Q_{h}$,

[TABLE]

Proof 4

Let $q_{h}\in Q_{h}$. It is well known, e.g., [28, Theorem 6.5], that since $q_{h}\in L^{2}_{0}(\Omega)$ there exists a $v\in[H_{0}^{1}(\Omega)]^{\rm dim}$ such that

[TABLE]

Then, by lemma 4 (i.), since $q_{h}\in P_{k-1}(K)$ for all $K\in\mathcal{T}$,

[TABLE]

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

[TABLE]

On the other hand, since $v\in\mathinner{[H^{1}(\Omega^{d})]}^{\rm dim}$ such that $v=0$ on $\Gamma^{d}$, using lemma 4 (iv.),

[TABLE]

Next, by eqs. 31a and 17,

[TABLE]

Combining eqs. 35, 36 and 37, we obtain

[TABLE]

Therefore, by eqs. 34 and 38,

[TABLE]

where we used eq. 33 for the last inequality. \qed

To prove the next auxiliary result, we introduce an operator to lift $\bar{q}_{h}^{s}\in\bar{Q}_{h}^{s}$ and $\bar{q}_{h}^{d}\in\bar{Q}_{h}^{d}$ to $\Omega^{s}$ and $\Omega^{d}$, respectively. Let $R_{k}(\partial K)\mathrel{\mathop{\mathchar 58\relax}}=\mathinner{\{\bar{q}\in L^{2}(\partial K),\ \bar{q}\in P_{k}(F),\ \forall F\in\partial K\}}$ and let $L\mathrel{\mathop{\mathchar 58\relax}}R_{k}(\partial K)\rightarrow\mathinner{[P_{k}(K)]}^{\rm dim}$ be the BDM local lifting operator which satisfies for all $\bar{q}_{h}\in R_{k}(\partial K)$ the following:

[TABLE]

Furthermore, it holds that

[TABLE]

See [33, Example 2.5.1].

Lemma 6

There exists a constant $\bar{c}_{{\rm inf}}>0$, independent of $h$, such that for any $(\bar{q}_{h}^{s},\bar{q}_{h}^{d})\in\bar{Q}_{h}^{s}\times\bar{Q}_{h}^{d}$,

[TABLE]

Proof 5

Let $\bar{q}_{h}^{j}\in\bar{Q}_{h}^{j}$ and set $w_{h}^{j}\mathrel{\mathop{\mathchar 58\relax}}=L\bar{q}_{h}^{j}$ for $j=s,d$. Since $w_{h}^{j}\in\mathinner{[P_{k}(K)]}^{\rm dim}$ for all $K\in\mathcal{T}^{j}$, $j=s,d$, then $w_{h}=w_{h}^{s}+w_{h}^{d}\in\mathinner{[P_{k}(K)]}^{\rm dim}$ for all $K\in\mathcal{T}$. Furthermore, by eq. 41,

[TABLE]

Next, using eq. 40,

[TABLE]

where we used $h_{K}<1<h_{K}^{-1}$. Then, combining eqs. 43 and 44,

[TABLE]

Noting that $w_{h}^{j}\cdot n=\bar{q}_{h}^{j}$, $j=s,d$ we find

[TABLE]

The result follows with \bar{c}_{{\rm inf}}=C\mathinner{\Bigl{(}\tfrac{h_{\min}}{h_{\max}}\Bigr{)}}. \qed

Using lemma 5 and lemma 6, the proof to theorem 1 then follows the steps as [34, Lemma 1] and is therefore omitted.

An immediate consequence of lemma 2 and theorem 1 is existence and uniqueness of a solution to eq. 5, as shown by the following proposition.

Proposition 1** (Existence and uniqueness)**

If $\beta>\beta_{0}$, there exists a unique solution $(\boldsymbol{u}_{h},\boldsymbol{p}_{h})\in\boldsymbol{X}_{h}$ to eq. 5.

Proof 6

Since eq. 5 is linear, it is sufficient to show uniqueness. Let $f^{s}=f^{d}=0$. Let us take $\boldsymbol{v}_{h}=\boldsymbol{u}_{h}$ and $\boldsymbol{q}_{h}=-\boldsymbol{p}_{h}$ in eq. 5. Then $a_{h}(\boldsymbol{u}_{h},\boldsymbol{u}_{h})=0$ which, together with the coercivity result eq. 26, implies that $\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\boldsymbol{u}_{h}\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}_{v}=0$, that is, $\boldsymbol{u}_{h}=0$. Inserting this into eq. 5 with $f^{s}=f^{d}=0$, we find

[TABLE]

By the inf-sup condition in theorem 1 this implies that $\boldsymbol{p}_{h}=0$. \qed

5 Error analysis of the Stokes–Darcy system

In this section we present a priori error analysis of the method in eq. 5. For this we require the standard $L^{2}$-projection operators onto $Q_{h}$ and $\bar{Q}_{h}^{j}$, $j=s,d$ which we denote, respectively, by $\Pi_{Q}$ and $\bar{\Pi}_{Q}^{j}$. It can be shown, e.g. [28], that for $k\geq 0$ and $0\leq\ell\leq k$,

[TABLE]

where $C$ is a generic constant independent of $h$.

It will be convenient to split the error into approximation and interpolation errors:

[TABLE]

where

[TABLE]

To be consistent with our notation, we use $\boldsymbol{\xi}_{u}\mathrel{\mathop{\mathchar 58\relax}}=(\xi_{u},\bar{\xi}_{u})$, $\boldsymbol{\xi}_{p}\mathrel{\mathop{\mathchar 58\relax}}=(\xi_{p},\bar{\xi}_{p}^{s},\bar{\xi}_{p}^{d})$ and $\boldsymbol{\xi}_{p}^{j}\mathrel{\mathop{\mathchar 58\relax}}=(\xi_{p},\bar{\xi}_{p}^{j})$, $j=s,d$. Expressions $\boldsymbol{\zeta}_{u}$, $\boldsymbol{\zeta}_{p}$ and $\boldsymbol{\zeta}_{p}^{j}$ are defined similarly. We first present three results that will be useful in following sections.

Lemma 7** ($\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\cdot\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}_{v,s}$ and $\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\cdot\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}_{v}$-norm

interpolation error estimates)**

Suppose that $u$ is such that $u^{s}\in\mathinner{[H^{\ell}(\Omega^{s})]}^{\rm dim}$, and $u^{d}\in\mathinner{[H^{\ell-1}(\Omega^{d})]}^{\rm dim}$, $2\leq\ell\leq k+1$. Then

[TABLE]

Proof 7

Recall that

[TABLE]

Then, since $\Pi_{V}u\in\mathinner{[P_{k}(F)]}^{\rm dim}$ for any $F\in\mathcal{F}^{s}$ we find by eq. 31b

[TABLE]

Equation 48a* follows by combining eq. 49 and eq. 50 and using once again lemma 4 (iv.). Equation 48b follows from this and lemma 4 (iv.). \qed*

Lemma 8** ($\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\cdot\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}_{p,j}$-norm interpolation error)**

For $j=s,d$, let $p^{j}\in H^{\ell}(\Omega^{j})$, $0\leq\ell\leq k$. Then

[TABLE]

Proof 8

By definition of $\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\cdot\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}_{p,j}$ and the properties of the $L^{2}$-projections $\Pi_{Q}$ and $\bar{\Pi}_{Q}^{j}$ given in eq. 47,

[TABLE]

\qed

Lemma 9** (Error equation)**

There holds

[TABLE]

Proof 9

Subtracting eq. 5 from the consistency equation in lemma 1 we obtain

[TABLE]

for all $(\boldsymbol{v}_{h},\boldsymbol{q}_{h})\in\boldsymbol{X}_{h}$. The result follows simply by splitting the errors. \qed

5.1 Energy-norm error estimates

In this section we determine error estimates for the velocity and pressure in the energy-norm.

Theorem 2** (Energy-norm error estimates)**

Let $(u,p)\in X$ be the solution of the Stokes–Darcy system eqs. 1 and 2 such that $u^{s}\in\mathinner{[H^{k+1}(\Omega^{s})]}^{\rm dim}$ and $u^{d}\in\mathinner{[H^{k}(\Omega^{d})]}^{\rm dim}$ for $k\geq 1$ and let $(\boldsymbol{u}_{h},\boldsymbol{p}_{h})\in\boldsymbol{X}_{h}$ solve eq. 5. Then

[TABLE]

Proof 10

We will first prove

[TABLE]

where $c_{1}>0$ is a constant independent of $h$. Setting $(\boldsymbol{v}_{h},\boldsymbol{q}_{h})=(\boldsymbol{\zeta}_{u},-\boldsymbol{\zeta}_{p})$ in lemma 9, and using coercivity of $a_{h}$ eq. 26,

[TABLE]

We proceed by bounding

[TABLE]

Observe first that $I_{2}$ disappears by the properties of $\Pi_{Q},\bar{\Pi}_{Q}^{j}$ and lemma 4 (i.)-(ii.). Consider now $I_{1}$. By lemma 3 and equivalence of the norms $\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\cdot\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}_{v}$ and $\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\cdot\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}_{v^{\prime}}$ on $\boldsymbol{V}_{h}$

[TABLE]

We next consider $I_{3}$ and note that by lemma 4 (i.) and (ii.)

[TABLE]

Since $(\bar{p}_{h}^{j}-\bar{\Pi}^{j}_{Q}p^{j})\in P_{k}(F)$ for $j=s,d$, we have by lemma 4 (ii.) that

[TABLE]

Equation 54* follows by combining eqs. 55 and 58.*

We next prove that

[TABLE]

where $c_{2}>0$ is a constant independent of $h$. Setting $\boldsymbol{q}_{h}=0$ in the error equation lemma 9, using the projection properties of $\Pi_{Q}$ and $\bar{\Pi}_{Q}$, applying eq. 28b and using the equivalence of the norms $\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\cdot\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}_{v}$ and $\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\cdot\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}_{v^{\prime}}$ on $\boldsymbol{V}_{h}$,

[TABLE]

Equation 59* follows by theorem 1.*

We will now combine eqs. 54 and 59. First note, that by Young’s inequality, we may bound eq. 54 as

[TABLE]

for any $\varepsilon>0$ and $\delta>0$. Multiplying eq. 60 by a constant $\beta>0$ and adding to eq. 59,

[TABLE]

Choosing $\varepsilon<2/c_{1}$, setting $\beta>c_{2}/(1-c_{1}\varepsilon/2)$ and choosing $\delta<2/(\beta c_{1})$ we may write

[TABLE]

resulting in

[TABLE]

The energy-norm error estimates eq. 53b now follow by the triangle inequality, the bounds in eq. 61, and the interpolation error estimates from lemmas 7 and 8. \qed

5.2 $L^{2}$-norm error estimate for the velocity

We will now determine $L^{2}$-norm error estimates for the velocity in the Stokes and Darcy regions separately. To obtain these estimates we consider the dual problem where $(U,P)$ is the solution to eqs. 1 and 2 with $f^{s}=\Psi\in\mathinner{[L^{2}(\Omega^{s})]}^{\rm dim}$ and $f^{d}=0$ [14, eq. (13)]. We will assume that this solution to the dual problem satisfies the following regularity estimates

[TABLE]

which will allow us to prove the following result.

Theorem 3** ($L^{2}$-norm error estimate for the velocity)**

Let $(u,p)\in X$ be the solution of eqs. 1 and 2 such that $u^{s}\in\mathinner{[H^{k+1}(\Omega^{s})]}^{\rm dim}$ and $u^{d}\in\mathinner{[H^{k+1}(\Omega^{d})]}^{\rm dim}$ for $k\geq 1$ and let $(\boldsymbol{u}_{h},\boldsymbol{p}_{h})\in\boldsymbol{X}_{h}$ solve eq. 5. Then

[TABLE]

Different arguments will be used to prove the $L^{2}$-norm error estimates for the velocity in the Stokes region eq. 63a and in the Darcy region eq. 63b. We therefore consider the proofs of these two inequalities separately.

Proof 11** (of Stokes error estimate in eq. 63a)**

Let $(U,P)\in X$ be the solution of the dual problem where $\Psi\mathrel{\mathop{\mathchar 58\relax}}=u-u_{h}$. Setting $\boldsymbol{U}\mathrel{\mathop{\mathchar 58\relax}}=(U,\gamma(U^{s}))$, $\boldsymbol{P}\mathrel{\mathop{\mathchar 58\relax}}=(P,\gamma(P^{s}),\gamma(P^{d}))$, lemma 1 implies

[TABLE]

Setting $\boldsymbol{v}=(u-u_{h},\gamma(u^{s})-\bar{u}_{h})$, $\boldsymbol{q}=(p-p_{h},\gamma(p^{s})-\bar{p}_{h}^{s},\gamma(p^{d})-\bar{p}_{h}^{d})$ in this equation,

[TABLE]

Next, setting $\boldsymbol{v}_{h}=(\Pi_{V}U,\bar{\Pi}_{V}U)$, $\boldsymbol{q}_{h}=(\Pi_{Q}P,\bar{\Pi}_{Q}^{s}P,\bar{\Pi}_{Q}^{d}P)$ in eq. 52, we obtain

[TABLE]

Combining eqs. 64 and 65 leads to

[TABLE]

First, observe that $J_{2}$ vanishes due to eqs. 1c and 1d, single-valuedness of $u$ and eqs. 11 and 12. We will bound the remaining terms starting with $J_{1}$. Using eq. 28b,

[TABLE]

By eqs. 48a and 31a and lemma 4 (iv.),

[TABLE]

where in the last step we used the regularity assumption eq. 62. We find

[TABLE]

We next bound $J_{3}$. Recalling eqs. 1c, 11 and 12, smoothness of $u$ and eq. 1e, we have

[TABLE]

where we used lemmas 8 and 62.

For $J_{4}$, observe that since $\nabla\cdot U=0$ and $\nabla\cdot\Pi_{V}U\in P_{k-1}(K)$,

[TABLE]

where we used $\Pi_{Q}p-p_{h}\in P_{k-1}(K)$ and lemma 4 (i.) in the last step. Using eq. 68 and noting that $U$ is smooth, eq. 12a holds, $\bar{\Pi}_{V}U$ is single-valued and $U=0$ (and so $\bar{\Pi}_{V}U=0$) on $\Gamma^{s}$, we have

[TABLE]

where we applied eq. 67.

Finally we bound $J_{5}$. Using eq. 68 and recalling that $U$ is smooth, eq. 12a holds, $\bar{\Pi}_{V}U$ is single-valued and $U\cdot n=0$ (and so $\Pi_{V}U\cdot n=0)$ on $\Gamma^{d}$,

[TABLE]

where we used eq. 67. The result follows by combining the bounds for $J_{1},\ldots,J_{5}$, cancelling a factor of $\mathinner{\!\left\lVert u-u_{h}\right\rVert}_{\Omega^{s}}$, using theorem 2 and the equivalence of $\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\cdot\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}_{v}$ and $\mathopen{|\mkern-1.5mu|\mkern-1.5mu|}\cdot\mathclose{|\mkern-1.5mu|\mkern-1.5mu|}_{v^{\prime}}$ on $\boldsymbol{V}_{h}$. \qed

To prove eq. 63b we will treat the problem in the Darcy region as a separate problem where the interface conditions are treated as boundary conditions. We will first state a result which constructs a vector function in

[TABLE]

such that its normal component on the interface matches the normal component of the error $u-u_{h}$ on the interface.

Lemma 10

Assume that $u^{s}\in\mathinner{[H^{k+1}(\Omega^{s})]}^{\rm dim}$ and $u^{d}\in\mathinner{[H^{k}(\Omega^{d})]}^{\rm dim}$ with $k\geq 1$. Then there exists a function $w\in H^{\rm div}(\Omega^{d})$ that satisfies $\nabla\cdot w=0$ in $\Omega^{d}$, $w\cdot n=0$ on $\Gamma^{d}$ and $w\cdot n=(u^{s}-u_{h}^{s})\cdot n$ on $\Gamma^{I}$ such that

[TABLE]

where $C>0$ is a constant independent of $h$.

We point out that the proof of lemma 10 relies on pointwise mass conservation eq. 11, $H(\text{div})$-conformity eq. 12 and eq. 63a, and can be found in [14, Lemma 3.1, Lemma 3.2]. We now prove eq. 63b.

Proof 12** (of the Darcy error estimate in eq. 63b)**

In eq. 52 set $v_{h}=0$ in $\Omega^{s}$, $\bar{v}_{h}=0$ and $\boldsymbol{q}_{h}=0$. Then

[TABLE]

by the projection properties of $\Pi_{Q}$ and $\bar{\Pi}_{Q}$.

We will first determine a $v_{h}$ so that $a_{h}^{d}((u-u_{h},v_{h}))=0$. Setting $q_{h}=0$ in $\Omega^{s}$, $\bar{q}_{h}^{s}=0$ and $\boldsymbol{v}_{h}=0$ in eq. 52, using lemma 4 (i.) and (ii.), we find

[TABLE]

where the second equality is by eqs. 1c and 11. Next, let $w\in H^{\rm div}(\Omega^{d})$ such that $\nabla\cdot w=0$ in $\Omega^{d}$, $w\cdot n=0$ on $\Gamma^{d}$, $w\cdot n^{d}=(\gamma(u^{s})-u_{h}^{s})\cdot n^{d}$ on $\Gamma^{I}$ as defined in lemma 10. Then, by lemma 4 (i.) and (ii.)

[TABLE]

where the second equality is by eqs. 1e and 12a. Combining eq. 71 and eq. 72, and using eqs. 2a and 12b, we find that

[TABLE]

Therefore, setting $v_{h}=\Pi_{V}u-u_{h}-\Pi_{V}w$ in eq. 70, we find

[TABLE]

Next, by definition of $a_{h}^{d}(\cdot,\cdot)$ and using eq. 74, we find that

[TABLE]

Hence, cancelling a factor of $\kappa^{-1}\|u-u_{h}\|_{\Omega^{d}}$ above and using lemma 4 (iv.) and eq. 69,

[TABLE]

The result follows. \qed

We note that the error estimates for the velocity in theorem 2 and theorem 3 do not depend on the approximation error of the pressure. The EDG–HDG method eq. 5 for the Stokes–Darcy system is therefore pressure-robust [35]. We furthermore remark that the analysis in this paper is easily extended to spatially dependent permeability that varies continuously over $\Omega^{d}$ under the condition that $0<\kappa_{\min}\leq\kappa\leq\kappa_{\max}<\infty$.

6 Numerical examples

The examples in this section have been implemented using the NGSolve finite element library [36]. Both examples are posed on the domain $\Omega=[0,1]^{2}$ with $\Omega^{d}=[0,1]\times[0,0.5]$ and $\Omega^{s}=[0,1]\times[0.5,1]$. The domain is discretized by an unstructured simplicial mesh. The mesh is such that $\mathcal{T}^{s}$ is an exact triangulation of $\Omega^{s}$, $\mathcal{T}^{d}$ is an exact triangulation of $\Omega^{d}$, and cell facets match on the interface $\Gamma^{I}$. We set $k=3$ in eq. 3 and eq. 4 and set the penalty parameter to be $\beta=10k^{2}$.

6.1 Rates of convergence

We solve the Stokes–Darcy system eqs. 1 and 2 with the source terms and boundary conditions chosen such that the exact solution is given by

[TABLE]

with $\alpha=\mu\kappa^{1/2}(1+4\pi^{2})/2$. For $\mu=1$ and $\kappa=1$ we recover the test case introduced in [37]. We introduce this extended version to investigate the effect of $\mu$ and $\kappa$ on the solution.

Table 1 presents computed errors in the velocity and the pressure on $\Omega^{s}$ and on $\Omega^{d}$ for different values of $\mu$ and $\kappa$. We observe optimal rates of convergence for all unknowns, independent of the value of the viscosity and permeability. We observe also that, although the error in the pressure changes when viscosity and permeability changes, the error in the velocity remains the same, as predicted by the analysis. Furthermore, pointwise satisfaction of the mass equation is obtained in both regions. Although not shown here, these observations hold also for other values of $k$ in eq. 3 and eq. 4.

6.2 Coupled surface and subsurface flow

In this section we consider a test case similar to that proposed in [38, Example 7.2]. This test case is representative of surface flow coupled to subsurface flow.

Let the boundary of the Stokes region be partitioned as $\Gamma^{s}=\Gamma^{s}_{1}\cup\Gamma_{2}^{s}\cup\Gamma_{3}^{s}$ where $\Gamma^{s}_{1}\mathrel{\mathop{\mathchar 58\relax}}=\mathinner{\left\{x\in\Gamma^{s}\ \mathrel{\mathop{\mathchar 58\relax}}\ x_{1}=0\right\}}$, $\Gamma^{s}_{2}\mathrel{\mathop{\mathchar 58\relax}}=\mathinner{\left\{x\in\Gamma^{s}\ \mathrel{\mathop{\mathchar 58\relax}}\ x_{1}=1\right\}}$ and $\Gamma^{s}_{3}\mathrel{\mathop{\mathchar 58\relax}}=\mathinner{\left\{x\in\Gamma^{s}\ \mathrel{\mathop{\mathchar 58\relax}}\ x_{2}=1\right\}}$. Similarly, let $\Gamma^{d}=\Gamma^{d}_{1}\cup\Gamma_{2}^{d}$ where $\Gamma^{d}_{1}\mathrel{\mathop{\mathchar 58\relax}}=\mathinner{\left\{x\in\Gamma^{d}\ \mathrel{\mathop{\mathchar 58\relax}}\ x_{1}=0\ \text{or}\ x_{1}=1\right\}}$ and $\Gamma^{d}_{2}\mathrel{\mathop{\mathchar 58\relax}}=\mathinner{\left\{x\in\Gamma^{d}\ \mathrel{\mathop{\mathchar 58\relax}}\ x_{2}=0\right\}}$. We impose the following boundary conditions:

[TABLE]

Note that due to the boundary condition on the pressure on $\Gamma_{2}^{d}$ the pressure solution need not be in $L^{2}_{0}(\Omega)$ as imposed in eq. 3. We set the permeability to

[TABLE]

Other parameters in eqs. 1 and 2 are set as $\mu=0.1$, $\alpha=0.5$, $f^{s}=0$ and $f^{d}=0$. We consider the solution on a mesh with 5754 simplicial cells.

Figure 2 shows the permeability, computed velocity field, and computed pressure field. We observe a similar flow field as presented in [38, Example 7.2]. In the Darcy region $\Omega^{d}$ the flow field avoids areas of low permeability, free-flow is present in the Stokes region $\Omega^{s}$, while the tangential velocity is discontinuous along the interface. To plot the pressure we used two colour scales to illustrate the pressure differences in the Darcy and the Stokes region. We note that the pressure is discontinuous across the interface. Furthermore, the pressure in the Darcy region follows a similar pattern as the permeability while in the Stokes region pressure is highest at the inlet.

7 Conclusions

We have formulated a pointwise mass-conserving and divergence-conforming embedded–hybridized discontinuous Galerkin method for the Stokes–Darcy system. The Stokes equations are coupled to the Darcy equations at an interface by the Beavers–Joseph–Saffman condition. This coupling is handled naturally by the facet variables that occur in the EDG–HDG method and act as Lagrange multipliers enforcing divergence-conformity of the scheme. We analyzed the embedded–hybridized method and proved existence and uniqueness of the solution to the discretization. Additionally, we presented a priori error analysis showing optimal rates of convergence and demonstrated that the error in the velocity is independent of the pressure. Numerical examples support the analysis of the scheme.

Acknowledgements

SR gratefully acknowledges support from the Natural Sciences and Engineering Research Council of Canada through the Discovery Grant program (RGPIN-05606-2015) and the Discovery Accelerator Supplement (RGPAS-478018-2015). Part of this research was done while AC was on sabbatical at the University of Waterloo, Department of Applied Mathematics.