A Robin-type domain decomposition method for a novel mixed-type DG method for the coupled Stokes-Darcy problem

Coupling incompressible flow and porous media flow has drawn great attention over the past years, which has been involved in many practical applications, such as ground water contamination and industrial filtration. This coupled phenomenon can be mathematically expressed by the Stokes-Darcy problem, where the free fluid region is governed by the Stokes equations and the porous media region is described by Darcy’s law, and three transmission conditions are prescribed on the interface (cf. [1, 29]).

The devising of numerical schemes for the coupled Stokes-Darcy problem hinges on a suitable choice of stable pairs for both the Stokes equations and Darcy equations. As it is well known, the standard mixed formulations for the Stokes equations and Darcy equations earn different compatibility conditions, thus a straightforward application of the existing solvers for the Stokes equations and Darcy equations may not be feasible. To this end, a great amount of effort has been devoted to developing accurate and efficient numerical schemes for the coupled Stokes-Darcy problem, and a non-exhaustive list of these approaches include Lagrange multiplier methods [21, 17, 32, 18], weak Galerkin method [9, 22], strongly conservative methods [20, 16], stabilized mixed finite element method [28, 24], discontinuous Galerkin (DG) methods [26, 34], virtual element method [23, 33], a lowest-order staggered DG method [37] and penalty methods [38]. The coupled Stokes-Darcy problem describes multiphysics phenomena, and involves a Stokes subproblem and a Darcy subproblem, it is thus natural to resort to domain decomposition methods, which allows one to solve the coupled system sequentially with a low computational cost. Various domain decomposition methods have been developed for the coupled Stokes-Darcy problem, see, for example [13, 6, 8, 7, 31, 19], most of which are based on velocity-pressure formulation of the Stokes equations.

The first purpose of this paper is to devise and analyze a novel mixed-type method for the coupled Stokes-Darcy problem. In the proposed formulation, we use a mixed-type DG method in conjunction with stress-velocity formulation for the discretization of the Stokes subproblem, and we use staggered DG method introduced in [10] for the discretization of the Darcy subproblem. Unlike the schemes proposed in [27, 34], we enforce the normal continuity of velocity directly into the formulation of the method without resorting to Lagrange multiplier. This is based on our observation that the degrees of freedom for the Darcy velocity space is defined by using dual edge degrees of freedom and interior degrees of freedom, and the interface can be treated as the union of the primal edges. Therefore, the normal continuity of velocity can be simply imposed into the discrete formulation by replacing the Darcy’s normal velocity by the Stokes’ normal velocity, which in conjunction with a suitable decomposition of the stress variable on the interface yields the resulting formulation. It is worth mentioning that the normal continuity of velocity is satisfied exactly at the discrete level. The advantages of the proposed formulation are multifold: First, the symmetry of stress is strongly imposed. The stress variable has a physical meaning and its computation is also important from a practical point of view. Second, pressure is eliminated from the equation, which reduces the size of the global system. Third, the continuity of normal velocity can be imposed directly without using Lagrange multiplier. A rigorous convergence analysis is carried out for all the variables. The proof for the optimal convergence of $L^{2}$-error of Darcy velocity is non-trivial, to overcome this issue, we exploit a non-standard trace theorem. We remark that the proposed scheme utilizes stress-velocity formulation for the Stokes equations, which is rarely explored for the coupled Stokes-Darcy problem in the existing literature. Our work will shed new insights into the devising and analysis of novel numerical schemes for the coupled Stokes-Darcy problem. We also address that the numerical results demonstrate that the proposed scheme also works well for small values of viscosity.

The proposed global formulation involves four variables: stress, fluid velocity, Darcy velocity and Darcy pressure, which may require high computational costs especially for large scale problems. To reduce the computational costs, we aim to devise a domain decomposition method based on the proposed spatial discretization, where the global formulation is decomposed into the Stokes subproblem and the Darcy subproblem by using newly constructed Robin-type interface conditions. The construction of the novel Robin-type interface condition is not a simple extension of the method introduced in [30], instead it takes advantage of the special features of staggered DG method. Moreover, the compatibility conditions are derived to ensure the equivalence of the modified discrete formulation and the original discrete formulation. The convergence of the Robin-type domain decomposition method is rigorously analyzed with the help of the compatibility conditions. Our convergence analysis indicates that the convergence of the Robin-type domain decomposition method is $1-\mathcal{O}(h)$ for $\delta_{p}=\delta_{f}$, where $\delta_{p}$ and $\delta_{f}$ are parameters introduced in Section 5. Moreover, when $\delta_{p}$ and $\delta_{f}$ satisfy certain conditions, then the convergence rate of the Robin-type domain decomposition method is $h$-independent, which is particularly inspiring. Several numerical experiments are carried out to demonstrate the performance of the proposed scheme. We can observe that the proposed domain decomposition method converges as reflected by the theories. In particular, it also converges for small values of viscosity.

The rest of the paper is organized as follows. In the next section, we describe the model problem and derive the stress-velocity formulation for the Stokes equations. In Section 3, we derive the discrete formulation and prove the unique solvability. The convergence error estimates for all the variables measured in $L^{2}$-error are given in Section 4. Then, the Robin-type domain decomposition method is constructed and analyzed in Section 5. Several numerical experiments are carried out in Section 6 to verify the proposed theories. Finally, a concluding remark is given in Section 7.

Let $\Omega:=\Omega_{S}\cup\Omega_{D}$ denote the polygonal domain in $\mathbb{R}^{2}$, where $\Omega_{S}$ and $\Omega_{D}$ represent the fluid domain and porous media domain, respectively. Let $\Gamma$ denote the interface between $\Omega_{S}$ and $\Omega_{D}$, and let $\Gamma_{i}=\partial\Omega_{i}\backslash\Gamma$ for $i\in\{S,D\}$, see Figure 1 for an illustration of the computational domain. We use $\bm{n}_{i}(i=S,D)$ to represent the unit normal vector to $\partial\Omega_{i}$. Let $\bm{t}_{S}$ be an orthonormal system of tangential vectors on $\Gamma$.

In $\Omega_{S}$, the fluid flow is governed by the Stokes equations

	$\displaystyle-\mbox{div}\,\underline{\sigma}$	$\displaystyle=\bm{f}_{S}$	$\displaystyle\quad\mbox{in}\;\Omega_{S},$		(2.1)
	$\displaystyle\underline{\sigma}$	$\displaystyle=2\mu\varepsilon(\bm{u}_{S})-pI$	$\displaystyle\quad\mbox{in}\;\Omega_{S},$		(2.2)
	$\displaystyle\text{div}\,\bm{u}_{S}$	$\displaystyle=0$	$\displaystyle\quad\mbox{in}\;\Omega_{S},$		(2.3)
	$\displaystyle\bm{u}_{S}$	$\displaystyle=\bm{0}$	$\displaystyle\quad\mbox{on}\;\Gamma_{S},$		(2.4)

where $\mu$ is the viscosity coefficient which is assumed to be a positive constant, $\bm{u}_{S}$ is the fluid velocity, $p$ is the fluid pressure, $\underline{\sigma}$ is the stress tensor, $I$ is the $2\times 2$ identity matrix and $\varepsilon(\bm{u}_{S})$ is the deformation tensor defined by $\varepsilon(\bm{u}_{S})=\frac{\nabla\bm{u}_{S}+\nabla\bm{u}_{S}^{\prime}}{2}$. Hereafter, we use $A^{\prime}$ to represent the transpose of $A$.

In $\Omega_{D}$, the governing equations are Darcy equations

	$\displaystyle\mathrm{div}\bm{u}_{D}$	$\displaystyle=f_{D}\;\mbox{in}\;\Omega_{D},$		(2.5)
	$\displaystyle\bm{u}_{D}+K\nabla p_{D}$	$\displaystyle=0\quad\mbox{in}\;\Omega_{D},$		(2.6)
	$\displaystyle p_{D}$	$\displaystyle=0\quad\mbox{on}\;\Gamma_{D}.$		(2.7)

On the interface $\Gamma$, we prescribe the following Beavers-Joseph-Saffman conditions (cf. [1, 29])

	$\displaystyle\bm{u}_{S}\cdot\bm{n}_{S}$	$\displaystyle=\bm{u}_{D}\cdot\bm{n}_{S}$	$\displaystyle\quad\mbox{on}\;\Gamma,$		(2.8)
	$\displaystyle-\underline{\sigma}\bm{n}_{S}\cdot\bm{n}_{S}$	$\displaystyle=p_{D}$	$\displaystyle\quad\mbox{on}\;\Gamma,$		(2.9)
	$\displaystyle\bm{u}_{S}\cdot\bm{t}_{S}$	$\displaystyle=-G\underline{\sigma}\bm{n}_{S}\cdot\bm{t}_{S}$	$\displaystyle\quad\mbox{on}\;\Gamma,$		(2.10)

where $G>0$ is the phenomenological friction coefficient.

Now we will derive a stress-velocity formulation based on (2.1)-(2.3) by eliminating $p$. First, we have

	$\displaystyle\mbox{tr}(\varepsilon(\bm{u}_{S}))=\mbox{div}\bm{u}_{S}=0,$
	$\displaystyle\mbox{tr}(\underline{\sigma})=2\mu\mbox{tr}(\varepsilon(\bm{u}_{S}))-2p=-2p,$

thus, we obtain

$\displaystyle p=-\frac{1}{2}\mbox{tr}(\underline{\sigma}).$

Consequently, we can recast (2.2) into the following equivalent formulation

$\displaystyle\mathcal{A}\underline{\sigma}=2\mu\varepsilon(\bm{u}_{S}),$

where

$\displaystyle\mathcal{A}\underline{\sigma}=\underline{\sigma}-\frac{1}{2}\mbox{tr}(\underline{\sigma})I.$

Note that $\mathcal{A}\underline{\sigma}$ is a trace-free tensor called deviatoric part and $\text{ker}(\mathcal{A})=\{qI\;|\;q\;\text{is a scalar function}\}$.

Then we can rewrite (2.1)-(2.3) as the following equivalent system:

	$\displaystyle-\mbox{div}\,\underline{\sigma}$	$\displaystyle=\bm{f}$	$\displaystyle\quad\mbox{in}\;\Omega_{S},$		(2.11)
	$\displaystyle\mathcal{A}\underline{\sigma}$	$\displaystyle=2\mu\varepsilon(\bm{u}_{S})$	$\displaystyle\quad\mbox{in}\;\Omega_{S},$		(2.12)
	$\displaystyle\bm{u}_{S}$	$\displaystyle=\bm{0}$	$\displaystyle\quad\mbox{on}\;\Gamma_{S}.$		(2.13)

We introduce some notations that will be used throughout the paper. Let $D\subset\mathbb{R}^{d},d=1,2$. By $(\cdot,\cdot)_{D}$, we denote the scalar product in $L^{2}(D):(p,q)_{D}:=\int_{D}p\,q\,dx$. We use the same notation for the scalar product in $L^{2}(D)^{2}$ and in $L^{2}(D)^{2\times 2}$. More precisely, $(\bm{\xi},\bm{w})_{D}:=\sum_{i=1}^{2}(\xi^{i},w^{i})$ for $\bm{\xi},\bm{w}\in L^{2}(D)^{2}$ and $(\underline{\psi},\underline{\zeta})_{D}:=\sum_{i=1}^{2}\sum_{j=1}^{2}(\psi^{i,j},\zeta^{i,j})_{D}$ for $\underline{\psi},\underline{\zeta}\in L^{2}(D)^{2\times 2}$. The associated norm is denoted by $\|\cdot\|_{0,D}$. Given an integer $m>0$, we denote the scalar-valued Sobolev spaces by $H^{m}(D)=W^{m,2}(D)$ with the norm $\|\cdot\|_{m,D}$ and seminorm $|\cdot|_{m,D}$. In addition, we use $H^{m}(D)^{d}$ and $H^{m}(D)^{d\times d}$ to denote the vector-valued and tensor-valued Sobolev spaces, respectively. In the following, we use $C$ to denote a generic constant independent of the meshsize which may have different values at different occurrences.

In this section, we are devoted to the derivation of a novel mixed-type DG method for the coupled Stokes-Darcy system (2.5)-(2.10) and (2.11)-(2.12). To this end, we first introduce the following meshes. Following [10, 36, 35], we first let $\mathcal{T}_{u,D}$ be the initial partition of the domain $\Omega_{D}$ into non-overlapping triangular meshes as shown in the left panel of Figure 2. We use $\mathcal{F}_{h,\Gamma}$ to denote the set of edges lying on the interface $\Gamma$. We also let $\mathcal{F}_{pr,D}$ be the set of all edges excluding the interface edges in the initial partition $\mathcal{T}_{u,D}$ and $\mathcal{F}_{pr,D}^{0}\subset\mathcal{F}_{pr,D}$ be the subset of all interior edges of $\Omega_{D}$. For each triangular mesh $E$ in the initial partition $\mathcal{T}_{u,D}$, we construct the sub-triangulation by connecting an interior point $\nu$ to all the vertices of $E$. For simplicity, we select $\nu$ as the center point. We rename the union of these triangles sharing the common point $\nu$ by $S(\nu)$. Moreover, we will use $\mathcal{F}_{dl,D}$ to denote the set of all the new edges generated by this subdivision process and use $\mathcal{T}_{h,D}$ to denote the resulting quasi-uniform triangulation, on which our basis functions are defined. Here the sub-triangulation $\mathcal{T}_{h,D}$ satisfies standard mesh regularity assumptions (cf. [12]). This construction is illustrated in Figure 2, where the black solid lines are edges in $\mathcal{F}_{pr,D}^{0}$ and the red dotted lines are edges in $\mathcal{F}_{dl,D}$. For each interior edge $e\in\mathcal{F}_{pr,D}^{0}$, we use $D(e)$ to denote the union of the two triangles in $\mathcal{T}_{h,D}$ sharing the edge $e$, and for each boundary edge $e\in(\mathcal{F}_{pr,D}\cup\mathcal{F}_{h,\Gamma})\backslash\mathcal{F}_{pr,D}^{0}$, we use $D(e)$ to denote the triangle in $\mathcal{T}_{h,D}$ having the edge $e$, see Figure 2.

On the other hand, we let $\{\mathcal{T}_{h,S}\}$ be a family of shape-regular triangulations of $\bar{\Omega}_{S}$. For simplicity, we assume that the meshes between the two domains $\Omega_{S}$ and $\Omega_{D}$ match at the interface. We use $\mathcal{F}_{h,S}$ to denote the set of all edges of $\mathcal{T}_{h,S}$ excluding the edges lying on the interface. We use $\mathcal{F}_{h,S}^{0}$ to represent the subset of $\mathcal{F}_{h,S}$, i.e., $\mathcal{F}_{h,S}^{0}$ is the union of interior edges of $\bar{\Omega}_{S}$. In addition, we let $\mathcal{T}_{h}=\mathcal{T}_{h,S}\cup\mathcal{T}_{h,D}$ and $\mathcal{F}_{h}=\mathcal{F}_{h,S}\cup\mathcal{F}_{h,D}$. In what follows, $h_{e}$ stands for the length of edge $e\in\mathcal{F}_{h}$. For each triangle $T\in\mathcal{T}_{h}$, we let $h_{T}$ be the diameter of $T$ and $h=\max_{T\in\mathcal{T}_{h}}h_{T}$. For each interior edge $e$, we then fix $\bm{n}_{e}$ as one of the two possible unit normal vectors on $e$. When there is no ambiguity, we use $\bm{n}$ instead of $\bm{n}_{e}$ to simplify the notation. In addition, we use $\bm{t}$ to represent the corresponding unit tangent vector. To simplify the presentation, we only consider triangular meshes in this paper, and the extension to polygonal meshes will be investigated in the future paper.

For $k\geq 1$, $T\in\mathcal{T}_{h}$ and $e\in\mathcal{F}_{h}$, we define $P^{k}(T)$ and $P^{k}(e)$ as the spaces of polynomials of degree up to order $k$ on $T$ and $e$, respectively. For a scalar or vector function $v$ belonging to the broken Sobolev space, its jump and average on an interior edge $e$ are defined as

$$\left\llbracket v\right\rrbracket_{e}=v_{1}-v_{2},\quad\left\{\!\!\left\{v\right\}\!\!\right\}_{e}=\frac{v_{1}+v_{2}}{2},$$

where $v_{j}=v_{T_{j}},j=1,2$ and $T_{1}$, $T_{2}$ are the two triangles in $\mathcal{T}_{h}$ having the edge $e$. For the boundary edges, we simply define $\left\llbracket v\right\rrbracket_{e}=v_{1}$ and $\left\{\!\!\left\{v\right\}\!\!\right\}_{e}=v_{1}$. We can omit the subscript $e$ when it is clear which edge we are referring to.

We define the following spaces for the numerical approximations.

	$\displaystyle\Sigma_{h}^{S}:$	$\displaystyle=\{\underline{w}:\underline{w}=\underline{w}^{\prime},\underline{w}|_{T}\in P_{k-1}(T)^{2\times 2},\;\forall T\in\mathcal{T}_{h,S}\},$
	$\displaystyle U_{h}^{S}:$	$\displaystyle=\{\bm{v}_{S}:\bm{v}_{S}|_{T}\in P_{k}(T)^{2},\;\forall T\in\mathcal{T}_{h,S};\bm{v}_{S}|_{\Gamma_{S}}=\bm{0}\}$

and

	$\displaystyle U_{h}^{D}:$	$\displaystyle=\{\bm{v}_{D}:\bm{v}_{D}|_{T}\in P_{k}(T)^{2},\;\forall T\in\mathcal{T}_{h,D},\left\llbracket\bm{v}_{D}\cdot\bm{n}\right\rrbracket_{e}=0\;\forall e\in\mathcal{F}_{dl,D}\},$
	$\displaystyle P_{h}^{D}:$	$\displaystyle=\{q:q|_{T}\in P_{k}(T),\;\forall T\in\mathcal{T}_{h,D},\left\llbracket q\right\rrbracket_{e}=0\;\forall e\in\mathcal{F}_{pr,D}^{0};q|_{\Gamma_{D}}=0\}.$

In the following, we use $\varepsilon_{h}(\bm{v})$ to represent the element-wise defined deformation tensor, i.e., $\varepsilon_{h}(\bm{v})\mid_{T}=\frac{\nabla\bm{v}\mid_{T}+\nabla\bm{v}^{\prime}\mid_{T}}{2}$ for $\bm{v}\in U_{h}^{S}$. Hereafter, we use $\nabla_{h}$ and $\mathrm{div}_{h}$ to represent the gradient and divergence operators defined element by element, respectively.

For later use, we specify the degrees of freedom for $P_{h}^{D}$ as follows.

• (SD1)

  For $e\in\mathcal{F}_{pr,D}^{0}\cup\mathcal{F}_{h,\Gamma}$, we have

  $$\phi_{e}(q):=(q,p_{k})_{e}\quad\forall p_{k}\in P_{k}(e).$$

  (3.1)

• (SD2)

  For each $T\in\mathcal{T}_{h,D}$, we define

  $$\phi_{T}(q):=(q,p_{k-1})_{T}\quad\forall p_{k-1}\in P_{k-1}(T).$$

  (3.2)

Now we are ready to derive the discrete formulation based on the first order system (2.5)-(2.6) and (2.11)-(2.12). Multiplying (2.1) by a test function $\bm{v}_{S}\in U_{h}^{S}$ and applying integration by parts yield

$\displaystyle-(\mathrm{div}\bm{\sigma},\bm{v}_{S})_{\Omega_{S}}$

$\displaystyle=-\sum_{e\in\mathcal{F}_{h,S}^{0}}(\left\{\!\!\left\{\bm{\sigma}\bm{n}\right\}\!\!\right\},\left\llbracket\bm{v}_{S}\right\rrbracket)_{e}-\sum_{e\in\mathcal{F}_{h,\Gamma}}(\bm{\sigma}\bm{n}_{S},\bm{v}_{S})_{e}+(\bm{\sigma},\varepsilon_{h}(\bm{v}_{S}))_{\Omega_{S}}.$

where the second term can be recast into the following form by using (2.9) and (2.10)

	$\displaystyle-\sum_{e\in\mathcal{F}_{h,\Gamma}}(\bm{\sigma}\bm{n}_{S},\bm{v}_{S})_{e}$	$\displaystyle=-\sum_{e\in\mathcal{F}_{h,\Gamma}}(\bm{\sigma}\bm{n}_{S},(\bm{v}_{S}\cdot\bm{n}_{S})\bm{n}_{S}+(\bm{v}_{S}\cdot\bm{t}_{S})\bm{t}_{S})_{e}$
		$\displaystyle=\sum_{e\in\mathcal{F}_{h,\Gamma}}(p_{D},\bm{v}_{S}\cdot\bm{n}_{S})_{e}+\frac{1}{G}\sum_{e\in\mathcal{F}_{h,\Gamma}}(\bm{u}_{S}\cdot\bm{t}_{S},\bm{v}_{S}\cdot\bm{t}_{S})_{e}.$

Multiplying (2.6) by a test function $\bm{v}_{D}\in U_{h}^{D}$ and performing integration by parts lead to

$\displaystyle(\nabla p_{D},\bm{v}_{D})_{\Omega_{D}}=\sum_{e\in\mathcal{F}_{h,\Gamma}}(\bm{v}_{D}\cdot\bm{n}_{D},p_{D})_{e}+\sum_{e\in\mathcal{F}_{pr,D}^{0}}(p_{D},\left\llbracket\bm{v}_{D}\cdot\bm{n}\right\rrbracket)_{e}-(p_{D},\mathrm{div}_{h}\bm{v}_{D})_{\Omega_{D}}.$

Finally we multiply (2.5) by $q_{D}\in P_{h}^{D}$, then using the interface condition (2.8) implies that

	$\displaystyle(\mathrm{div}\bm{u}_{D},q_{D})_{\Omega_{D}}$	$\displaystyle=\sum_{e\in\mathcal{F}_{dl,D}}(\bm{u}_{D}\cdot\bm{n},\left\llbracket q_{D}\right\rrbracket)_{e}+\sum_{e\in\mathcal{F}_{h,\Gamma}}(\bm{u}_{D}\cdot\bm{n}_{D},q_{D})_{e}-(\bm{u}_{D},\nabla_{h}q_{D})_{\Omega_{D}}$
		$\displaystyle=\sum_{e\in\mathcal{F}_{dl,D}}(\bm{u}_{D}\cdot\bm{n},\left\llbracket q_{D}\right\rrbracket)_{e}-(\bm{u}_{D},\nabla_{h}q_{D})_{\Omega_{D}}-\sum_{e\in\mathcal{F}_{h,\Gamma}}(\bm{u}_{S}\cdot\bm{n}_{S},q_{D})_{e}.$

Based on the above derivations, we are now in position to define the following bilinear forms, which is instrumental for later use.

$\displaystyle a_{S}(\bm{\sigma}_{h},\bm{v}_{S})$

$\displaystyle=-\sum_{e\in\mathcal{F}_{h,S}^{0}}(\left\{\!\!\left\{\bm{\sigma}_{h}\bm{n}\right\}\!\!\right\},\left\llbracket\bm{v}_{S}\right\rrbracket)_{e}+(\bm{\sigma}_{h},\varepsilon_{h}(\bm{v}_{S}))_{\Omega_{S}}$

and

	$\displaystyle b_{D}^{*}(p_{D},\bm{v}_{D})$	$\displaystyle=-\sum_{e\in\mathcal{F}_{pr,D}^{0}}(p_{D},\left\llbracket\bm{v}_{D}\cdot\bm{n}\right\rrbracket)_{e}+(p_{D},\mathrm{div}_{h}\bm{v}_{D})_{\Omega_{D}}-\sum_{e\in\mathcal{F}_{h,\Gamma}}(\bm{v}_{D}\cdot\bm{n}_{D},p_{D})_{e},$
	$\displaystyle b_{D}(\bm{u}_{D},q_{D})$	$\displaystyle=\sum_{e\in\mathcal{F}_{dl,D}}(\bm{u}_{D}\cdot\bm{n},\left\llbracket q_{D}\right\rrbracket)_{e}-(\bm{u}_{D},\nabla_{h}q_{D})_{\Omega_{D}}.$

It follows from integration by parts that

$\displaystyle a_{S}(\bm{\sigma}_{h},\bm{v}_{S})=\sum_{e\in\mathcal{F}_{h,S}^{0}}(\left\llbracket\bm{\sigma}_{h}\bm{n}\right\rrbracket,\left\{\!\!\left\{\bm{v}_{S}\right\}\!\!\right\})_{e}-(\mathrm{div}_{h}\underline{\sigma}_{h},\bm{v}_{S})_{\Omega_{S}}\quad\forall\bm{v}_{S}\in U_{h}^{S}.$

The discrete formulation reads as follows: Find $(\underline{\sigma}_{h},\bm{u}_{S,h})\in\Sigma_{h}^{S}\times U_{h}^{S}$ and $(\bm{u}_{D,h},p_{D,h})\in U_{h}^{D}\times P_{h}^{D}$ such that

	$\displaystyle((2\mu)^{-1}\mathcal{A}\underline{\sigma}_{h},\underline{w})_{\Omega_{S}}=a_{S}(\underline{w},\bm{u}_{S,h})\quad\forall\underline{w}\in\Sigma_{h}^{S},$		(3.3)
	$\displaystyle a_{S}(\underline{\sigma}_{h},\bm{v}_{S})+(K^{-1}\bm{u}_{D,h},\bm{v}_{D})_{\Omega_{D}}-b_{D}^{*}(p_{D,h},\bm{v}_{D})+\sum_{e\in\mathcal{F}_{h,\Gamma}}(\bm{v}_{S}\cdot\bm{n}_{S},p_{D,h})_{e}$
	$\displaystyle\;+\sum_{e\in\mathcal{F}_{h,S}^{0}}\gamma h_{e}^{-1}(\left\llbracket\bm{u}_{S,h}\right\rrbracket,\left\llbracket\bm{v}_{S}\right\rrbracket)_{e}+\frac{1}{G}\sum_{e\in\mathcal{F}_{h,\Gamma}}(\bm{u}_{S,h}\cdot\bm{t}_{S},\bm{v}_{S}\cdot\bm{t}_{S})_{e}=(\bm{f}_{S},\bm{v}_{S})_{\Omega_{S}}\quad\forall\bm{v}_{S}\in U_{h}^{S},$		(3.4)
	$\displaystyle b_{D}(\bm{u}_{D,h},q_{D})-\sum_{e\in\mathcal{F}_{h,\Gamma}}(\bm{u}_{S,h}\cdot\bm{n}_{S},q_{D})_{e}=(f_{D},q_{D})_{\Omega_{D}}\quad\forall q_{D}\in P_{h}^{D},$		(3.5)

where $\gamma>0$ is a constant over each edge $e\in\mathcal{F}_{h,S}^{0}$.

Integration by parts implies the following discrete adjoint property

$\displaystyle b_{D}(\bm{v}_{D},q_{D})$

$\displaystyle=b_{D}^{*}(q_{D},\bm{v}_{D})\quad\forall(\bm{v}_{D},q_{D})\in U_{h}^{D}\times P_{h}^{D}.$

(3.6)

To facilitate later analysis, we define the following mesh-dependent norm/semi-norm

	$\displaystyle\|\bm{v}_{D}\|_{0,h}^{2}$	$\displaystyle=\|\bm{v}_{D}\|_{0,\Omega_{D}}^{2}+\sum_{e\in\mathcal{F}_{dl,D}}h_{e}\|\bm{v}_{D}\cdot\bm{n}\|_{0,e}^{2},$
	$\displaystyle\|q_{D}\|_{Z}^{2}$	$\displaystyle=\|\nabla_{h}q_{D}\|_{0,\Omega_{D}}^{2}+\sum_{e\in\mathcal{F}_{dl,D}}h_{e}^{-1}\|\left\llbracket q_{D}\right\rrbracket\|_{0,e}^{2},$
	$\displaystyle\|\bm{v}_{S}\|_{1,h}^{2}$	$\displaystyle=\|\nabla_{h}\bm{v}_{S}\|_{0,\Omega_{S}}^{2}+\sum_{e\in\mathcal{F}_{h,S}^{0}}h_{e}^{-1}\|\left\llbracket\bm{v}_{S}\right\rrbracket\|_{0,e}^{2}$

for any $(\bm{v}_{D},q_{D},\bm{v}_{S})\in U_{h}^{D}\times P_{h}^{D}\times U_{h}^{S}$.

Following [10, 11], we have the following inf-sup condition

$\displaystyle\|q_{D}\|_{Z}\leq C\sup_{\bm{v}_{D}\in U_{h}^{D}}\frac{b_{D}(\bm{v}_{D},q_{D})}{\|\bm{v}_{D}\|_{0,\Omega_{D}}}\quad\forall q_{D}\in P_{h}^{D}.$

(3.7)

For later use, we also introduce the following space

$\displaystyle U_{h}^{\text{RT}}:$

$\displaystyle=\{\bm{v}\in H(\text{div};\Omega_{S}),\bm{v}|_{T}\in RT_{k-1},\forall T\in\mathcal{T}_{h,S};\bm{v}\cdot\bm{n}=0\;\mbox{on}\;\Gamma_{S}\},$

where $H(\text{div};\Omega_{S})=\{\bm{v}:\bm{v}\in L^{2}(\Omega_{S})^{2},\mathrm{div}\bm{v}\in L^{2}(\Omega_{S})\}$ and $RT_{k}(T)$ is the Raviart-Thomas element of index $k$ introduced in [25].

In this section we will present the convergence error estimates for all the variables. To begin, we define the following interpolation operators, which play an important role for the subsequent analysis. We let $\Pi_{h}^{\text{BDM}}$ denote the BDM projection operator (cf. [4]), which satisfies the following error estimates for $0\leq\alpha\leq k+1$ (see, e.g., [4, 14])

$\displaystyle\|\bm{v}-\Pi_{h}^{\text{BDM}}\bm{v}\|_{\alpha,\Omega_{S}}$

$\displaystyle\leq Ch^{k+1-\alpha}\|\bm{v}\|_{k+1,\Omega_{S}}\quad\forall\bm{v}\in H^{k+1}(\Omega_{S})^{2}.$

Let $\mathcal{I}_{h}:H^{1}(\Omega_{D})\rightarrow P_{h}^{D}$ be defined by

$$\begin{split}(\mathcal{I}_{h}q-q,\phi)_{e}&=0\quad\forall\phi\in P^{k}(e),\forall e\in\mathcal{F}_{pr,D},\\ (\mathcal{I}_{h}q-q,\phi)_{T}&=0\quad\forall\phi\in P^{k-1}(T),\forall T\in\mathcal{T}_{h,D}\end{split}$$

and $\mathcal{J}_{h}:L^{2}(\Omega_{D})^{2}\cap H^{1/2+\delta}(\Omega_{D})^{2}\rightarrow U_{h}^{D}$, $\delta>0$ be defined by

$$\begin{split}((\mathcal{J}_{h}\bm{v}-\bm{v})\cdot\bm{n},\varphi)_{e}&=0\quad\forall\varphi\in P^{k}(e),\forall e\in\mathcal{F}_{dl,D},\\ (\mathcal{J}_{h}\bm{v}-\bm{v},\bm{\phi})_{T}&=0\quad\forall\bm{\phi}\in P^{k-1}(T)^{2},\forall T\in\mathcal{T}_{h,D}.\end{split}$$

It is easy to see that $\mathcal{I}_{h}$ and $\mathcal{J}_{h}$ are well defined polynomial preserving operators. In addition, the following approximation properties hold for $q\in H^{k+1}(\Omega_{D})$ and $\bm{v}\in H^{k+1}(\Omega_{D})^{2}$ (cf. [12, 10])

	$\displaystyle\|q-\mathcal{I}_{h}q\|_{\alpha,\Omega_{D}}$	$\displaystyle\leq Ch^{k+1-\alpha}|q|_{k+1,\Omega_{D}}\quad 0\leq\alpha\leq k+1,$		(4.1)
	$\displaystyle\|\bm{v}-\mathcal{J}_{h}\bm{v}\|_{\alpha,\Omega_{D}}$	$\displaystyle\leq Ch^{k+1-\alpha}|\bm{v}|_{k+1,\Omega_{D}}\quad 0\leq\alpha\leq k+1.$		(4.2)

Finally, we let $\Pi_{h}$ denote the $L^{2}$ orthogonal projection onto $\Sigma_{h}^{S}$, then it holds

$\displaystyle\|\underline{w}-\Pi_{h}\underline{w}\|_{\alpha,\Omega_{S}}\leq Ch^{k-\alpha}\|\underline{w}\|_{k,\Omega_{S}}\quad\forall\underline{w}\in H^{k}(\Omega_{S})^{2\times 2},0\leq\alpha\leq k.$

(4.3)

Performing integration by parts on the discrete formulation (3.3)-(3.5) and using the interface conditions (2.8)-(2.10), we can get the following error equations:

	$\displaystyle((2\mu)^{-1}\mathcal{A}(\underline{\sigma}-\underline{\sigma}_{h}),\underline{w})_{\Omega_{S}}=a_{S}(\underline{w},\bm{u}_{S}-\bm{u}_{S,h})\quad\underline{w}\in\Sigma_{h}^{S},$		(4.4)
	$\displaystyle a_{S}(\underline{\sigma}-\underline{\sigma}_{h},\bm{v}_{S})+(K^{-1}(\bm{u}_{D}-\bm{u}_{D,h}),\bm{v}_{D})_{\Omega_{D}}-b_{D}^{*}(p_{D}-p_{D,h},\bm{v}_{D})+\sum_{e\in\mathcal{F}_{h,\Gamma}}(\bm{v}_{S}\cdot\bm{n}_{S},p_{D}-p_{D,h})_{e}$
	$\displaystyle\;+\sum_{e\in\mathcal{F}_{h,S}^{0}}\gamma h_{e}^{-1}(\left\llbracket\bm{u}_{S}-\bm{u}_{S,h}\right\rrbracket,\left\llbracket\bm{v}_{S}\right\rrbracket)_{e}+\frac{1}{G}\sum_{e\in\mathcal{F}_{h,\Gamma}}((\bm{u}_{S}-\bm{u}_{S,h})\cdot\bm{t}_{S},\bm{v}_{S}\cdot\bm{t}_{S})_{e}=0\quad\forall\bm{v}_{S}\in U_{h}^{S},$		(4.5)
	$\displaystyle b_{D}(\bm{u}_{D}-\bm{u}_{D,h},q_{D})-\sum_{e\in\mathcal{F}_{h,\Gamma}}((\bm{u}_{S}-\bm{u}_{S,h})\cdot\bm{n}_{S},q_{D})_{e}=0\quad\forall q_{D}\in P_{h}^{D},$		(4.6)

which indicates that our discrete formulation is consistent.