text:
\section{Introduction}
In this paper, we consider a two-dimensional Stokes interface problem. For simplicity, we describe this problem using the following model. In the model, we consider a bounded domain $\Omega \subset \mathbb{R}^2$ that is divided into two subdomains $\Omega_1$ and $\Omega_2$. In each subdomain $\Omega_i,i=1,2$, the fluid flow is governed by the Stokes equation, i.e.,

\begin{eqnarray}
-\nabla\cdot(A_i \nabla {\bf u}_i)+\nabla p_i &=&{\bf f}_i ,\,\, \qquad\text{in}\ \Omega_i,\label{model 1} \\
\nabla\cdot {\bf u}_i &=&0 , \,\,\,\qquad\text{in}\ \Omega_i,\\
{\bf u}_i&=&{\bf g}_i , \qquad\text{on}\ \partial \Omega_i \backslash \Gamma ,
\end{eqnarray}
where the viscosity coefficient $A_i$ is piecewise constant. And  $\Gamma=\Omega_1 \cap \Omega_2$ denotes the interface of two subdomains. The interface conditions on $\Gamma$ are described by the following equations:
\begin{eqnarray}
{\bf u}_1-{\bf u}_2&=&\bm{\psi}, \qquad\text{on}\ \Gamma \label{interface 1}\\
(A_1 \nabla {\bf u}_1-p_1 I)  \mathbf{n_1} +(A_2 \nabla {\bf u}_2-p_2 I)  \mathbf{n_2}&=&\bm{\phi}, \qquad\text{on}\ \Gamma \label{interface 2}
\end{eqnarray}
where $\mathbf{n_1}$ and $\mathbf{n_2}$ denote unit outward normals of $\Omega_1$ and $\Omega_2$. The problem is usually used in geology and geophysics. For example, the mantle flow problem is Stokes equations with discontinuous viscosity $A$. This type of problem is investigated by many geophysicists and mathematicians \cite{Stokesinterfacebackground4,Stokesinterfacebackground1,Stokesinterfacebackground6,Stokesinterfacebackground2,Stokesinterfacebackground5,Stokesinterfacebackground3}.

In practical problems, the interface is usually a complex curve or surface. In two-dimensional
problems, for such interfaces, we usually use straight edges to approximate curved edges. Then the method on straight edge element is used to solve the problem numerically. However, when we use higher degree polynomials to approximate the solution, the geometric error of the straight edge approximating the curved edge will reduce the convergence order of the error between the numerical solution and the exact solution \cite{reduceerror3,reduceerror1,reduceerror2}. Therefore, for the domain with curved edges, we need to reduce the geometric error as much as possible. Many numerical methods solve the problem on the curved domain directly. For example, p-FEM method \cite{pFEM} establishes the relationship between local and Cartesian coordinates, isogeometric methods \cite{isogeometricmethods} uses the NURBS (non-uniform rational B-splines) to approximate the whole computational domain and the basis of the approximation functions. NURBS-enhanced finite element method (NEFEM) \cite{NEFEM,reduceerror3} has the similar idea. It employs NURBS to approximate the curved domain, and uses the finite element method to solve the equations. The standard finite element method is used in the interior elements (for elements not intersecting the curved boundary), and the NURBS is used for approximation of the solution in the boundary elements (for elements intersecting the curved boundary).

In this paper, we use the weak Galerkin (WG) finite element method to solve the interface problem on the domain with curved boundary. The WG method was first proposed in \cite{wang2013weak} to deal with second-order elliptic equations of the primal scheme. The WG method is characterized by replacing continuous functions with weak functions and differential operators with weak operators. Later, the WG method was also used to solve the Stokes equation \cite{WGStokes1,WGStokes2}, Brinkman equation \cite{WGBrinkman1,WGBrinkman2}, linear elasticity equation \cite{elasticityequation}. Due to the discontinuity of weak functions, the WG is more flexible in dealing with interface problems \cite{Stokesinterfacebackground6,ellipticinterface,stokesDarcy}. For example, when dealing with elliptic interface problems, it can be the same value on the interface \cite{ellipticinterface}, and the interface condition is reflected in the variational problem, or it can be two values, and the interface condition is transformed into Dirichlet boundary condition. In this paper, we take the idea of the NEFEM method. In the interior element, we use the WG to treat, that is, we define the basis function on the reference element and obtain the basis function on the actual element by affine transformation \cite{nonaffineisoparametric1,nonaffineisoparametric2}. For interface elements, we directly define basis functions on the actual elements. For numerical integration on the actual element, we establish a non-affine transformation on the reference element, and then perform numerical integration on the reference element.

An outline of this paper is as follows: In Section 2, some notations used in this paper are presented. We also give the definitions of weak finite element space and weak operator, and the numerical scheme of weak Galerkin finite element. The Section 3 is devoted to the proof of the existence, uniqueness and stability of the solution of the weak Galerkin finite element numerical scheme. In Section 4 and 5, the error estimation under the energy norm and the $L^2$ norm are proved seperately. Finally, we give some numerical examples to verify our proposed theory in Section 6.

\section{The numerical scheme}In the section, we first give the introduction of notations used in the paper. Next, we will define the weak Galerkin finite element space in the partiton of the domain $\Omega$ and the corresponding weak operators. Finally, the weak Galerkin finite element numerical scheme of Stokes interface problem is given.
\subsection{The notations}
Assume that $\mathcal{T}_h$ is the triangulation of the domain $\Omega$ containing both straight and curved triangular cell. For $T \in \mathcal{T}_h$, when all three edges of $T$ are straight, it is called a triangular cell. The set of such triangles is denoted by $\mathcal{T}_h^S$. When $T$ has a curved edge which intersects with the interface edge, it is called a curved cell. The set of curved elements is denoted by $\mathcal{T}_h^I$. The set of the curved elements in $\Omega_1$ is denoted by $\mathcal{T}_{1h}^I$. The set of the curved elements in $\Omega_2$ is denoted by $\mathcal{T}_{2h}^I$. The triangular elements and the straight elements formed by three vertices of curved elements should satisfy the regularity condition \cite{WGStokes1}.
For  $T \in \mathcal{T}_h$, let $|T|$ and $h_T$ be the area and the diameter of the element $T$, respectively. Denote $h=\max_{T \in \mathcal{T}_h} h_T$ is the mesh size. Let $\mathcal{E}_h$ be the set of all edges in the partition $\mathcal{T}_h$ and $\mathcal{E}_h^0$ be the set of all interior edges. Denote $\mathcal{E}_h^S$ is the set of all straight edges and $\mathcal{E}_h^I$ is the set of all interface edges.
\subsection{Weak finite element space}
We define the weak function space on the partition $\mathcal{T}_h$. For $T \in \mathcal{T}_h$, we define the weak function ${\mathbf v}= \{\ {\mathbf v}_0, {\mathbf v}_b \}$ on the element $T$, where ${\mathbf v}_0$ represents the value in the interior of T and ${\mathbf v}_b$ represents the value on the boundary $\partial T$. Note ${\mathbf v}_b$ has one unique value on an edge $e \in \mathcal{E}_h$. There is no relationship between the value of ${\mathbf v}_0$ and ${\mathbf v}_b$. Next, we give the weak Galerkin finite element space with respect to the velocity ${\bf u}$ in the domain. For the pressure $p$, we define piecewise polynomial space for approximation.

For a given integer $k \geqslant 1$,
\begin{align*}
V_h=&\{{\mathbf v}=\{ \mathbf{v_0}, \mathbf{v_b}\},\mathbf{v_0}|_{T} \in [P_k(T)]^2, \,T \in \mathcal{T}_h,\\
&\mathbf{v_b}|_e \in [P_{k-1}(e)]^2, e \in \mathcal{E}_h^s,
\tilde{{\mathbf v}}_b|_e \in [P_k(e)]^{2 \times 2}, e \in \mathcal{E}_h^I
 \}\\
 V_h^0=&\{{\mathbf v}=\{ \mathbf{v_0}, \mathbf{v_b}\} \in V_h, \mathbf{v_b}=\mathbf{0} ~on ~ \partial \Omega \}\\
 W_h=&\{q: q \in L_0^2(\Omega), q|_T \in P_{k-1}(T),\, T \in \mathcal{T}_h\}
\end{align*}

\subsection{Discrete weak operators}
 In this subsection, we introduce some weak differential operators.

\begin{definition}\label{weak gradient}
For each ${\mathbf v} \in V_h$, its discrete weak gradient is denoted by $\nabla_w {\mathbf v} \in [P_{k-1}(T)]^{2\times2}$ and satisfies the following equations:

(1) for T $\in \mathcal{T}_h^S$ and $\forall q \in [P_{k-1}(T)]^{2\times2}$,
\begin{equation}\label{weak gradient 1}
(\nabla_w {\mathbf v},q)_T=-({\mathbf v}_0,\nabla \cdot q)_T + \langle {\mathbf v}_b,q {\bf n} \rangle_{\partial T},
\end{equation}
\qquad (2) for T $\in \mathcal{T}_{1h}^I$ and $\forall q \in [P_{k-1}(T)]^{2\times2}$,
\begin{eqnarray}\label{weak gradient 2}
(\nabla_w {\mathbf v},q)_T=-({\mathbf v}_0,\nabla \cdot q)_T +\langle {\mathbf v}_{1b},q {\bf n} \rangle_{\partial T \backslash \Gamma}
+\langle \tilde{{\mathbf v}}_{1b},q \rangle_{\partial T \cap \Gamma},
\end{eqnarray}
\qquad (3) for T $\in \mathcal{T}_{2h}^I$ and $\forall q \in [P_{k-1}(T)]^{2\times2}$,
\begin{eqnarray}\label{weak gradient 3}
(\nabla_w {\mathbf v},q)_T=-({\mathbf v}_{20},\nabla \cdot q)_T +\langle {\mathbf v}_{2b},q {\bf n} \rangle_{\partial T \backslash \Gamma}
-\langle \tilde{{\mathbf v}}_{1b},q \rangle_{\partial T \cap \Gamma}.
\end{eqnarray}
\end{definition}
\qquad Similarly, we can give the definition of discrete weak divergence.
\begin{definition}\label{weak divergence}
For each ${\mathbf v} \in V_h$, its discrete weak divergence is denoted by $\nabla_w \cdot {\mathbf v} \in P_{k-1}(T)$ and satisfies the following equations:

(1) for T $\in \mathcal{T}_h^S$ and $\forall q \in P_{k-1}(T)$,
\begin{eqnarray}\label{weak divergence 1}
(\nabla_w \cdot {\mathbf v},q)_T=-({\mathbf v}_0, \nabla q)_T + \langle{\mathbf v}_b \cdot {\bf n}, q \rangle_{\partial T},
\end{eqnarray}
\qquad (2) for T $\in \mathcal{T}_{1h}^I$ and $\forall q \in P_{k-1}(T)$,
\begin{eqnarray}\label{weak divergence 2}
(\nabla_w \cdot {\mathbf v},q)_T=-( {\mathbf v}_0,\nabla q)_T + \langle {\mathbf v}_{1b} \cdot \mathbf{n},q \rangle_{\partial T \backslash \Gamma}
+\langle \tilde{{\mathbf v}}_{1b},qI \rangle_{\partial T \cap \Gamma},
\end{eqnarray}
\qquad (3) for T $\in \mathcal{T}_{2h}^I$ and $\forall q \in P_{k-1}(T)$,
\begin{eqnarray}\label{weak divergence 3}
(\nabla_w \cdot {\mathbf v},q)_T=-({\mathbf v}_{20},\nabla q)_T + \langle {\mathbf v}_{2b} \cdot \mathbf{n},q \rangle_{\partial T \backslash \Gamma}
-\langle \tilde{{\mathbf v}}_{1b},qI \rangle_{\partial T \cap \Gamma}.
\end{eqnarray}
where I is the identity matrix.
\end{definition}

\subsection{The numerical scheme}
Let $Q_0$ be the $L^2$ projection operator from $[L^2(T)]^2$ onto $[P_k(T)]^2$, $T \in \mathcal{T}_h$. For each edge $e \in \mathcal{E}_h^S$, $Q_b$ is defined as the $L^2$ projection operator from $L^2(e)$ onto $P_{k-1}(e)$. Denote by $\hat{Q}_b$ the $L^2$ projection operator from $[L^2(e)]^{2 \times 2}$ onto $[P_{k}(e)]^{2 \times 2}$, $e \in \mathcal{E}_h^I$. Next, we introduce two bilinear forms as follows:
\begin{align*}
s({\bf u},{\mathbf v})=&\sum_{T \in \mathcal{T}_h^S }h_T^{-1} \langle Q_b {\bf u}_0 -{\bf u}_b,Q_b {\mathbf v}_0 -{\mathbf v}_b \rangle_{\partial T}\\
+&\sum_{T \in \mathcal{T}_{1h}^I}h_T^{-1} \Big( \langle Q_b {\bf u}_0 -{\bf u}_b,Q_b {\mathbf v}_0 -{\mathbf v}_b \rangle_{\partial T \backslash \Gamma}\\
+&\langle \hat{Q}_b({\bf u}_0 {\bf n}_1^{\mathrm{T}})-\tilde{{\bf u}}_b,\hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})-\tilde{{\mathbf v}}_b \rangle_{\partial T \cap \Gamma} \Big)\\
+&\sum_{T \in \mathcal{T}_{2h}^I}h_T^{-1} \Big(\langle Q_b {\bf u}_0 -{\bf u}_b,Q_b {\mathbf v}_0 -{\mathbf v}_b \rangle_{\partial T \backslash \Gamma}\\
+&\langle \hat{Q}_b({\bf u}_0 {\bf n}_2^{\mathrm{T}})+\tilde{{\bf u}}_b,\hat{Q}_b({\mathbf v}_0 {\bf n}_2^{\mathrm{T}})+\tilde{{\mathbf v}}_b \rangle_{\partial T \cap \Gamma}\Big),\\
a({\bf u},{\mathbf v})=&\sum_{T \in \mathcal{T}_h}(A\nabla_w {\bf u},\nabla_w {\mathbf v})_T+s({\bf u},{\mathbf v}),\\
b({\mathbf v},p)=&\sum_{T \in \mathcal{T}_h}(\nabla_w \cdot {\mathbf v},p)_T.
\end{align*}

\begin{algorithm}
For the Stokes interface problem (\ref{model 1})-(\ref{interface 2}), seeking $\mathbf{u}_h \in V_h, {\bf u}_b=Q_b {\bf g} ~\text{on}\ \partial \Omega$ and $p_h \in W_h$ such that
\begin{eqnarray}\label{algorithm 1}
\begin{split}
a({\bf u}_h,{\mathbf v}_h)-&b({\mathbf v}_h,p_h)=({\bf f},{\mathbf v}_0)-\langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}}),A_2 \nabla_w {\mathbf v}_h \rangle_{\Gamma}\\
+& \sum_{T \in \mathcal{T}_{2h}^I}h_T^{-1} \langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}}),\hat{Q}_b({\mathbf v}_0 {\bf n}_2^{\mathrm{T}})+\tilde{{\mathbf v}}_b \rangle_{\partial T \cap \Gamma}+\langle \bm{\phi} {\bf n}_1^{\mathrm{T}},\tilde{{\mathbf v}}_b \rangle_{\Gamma},
\end{split}
\end{eqnarray}
\begin{eqnarray}\label{algorithm 2}
b({\bf u}_h,q_h)=-\langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}}),q_hI \rangle_{\Gamma}.
\end{eqnarray}
for all ${\mathbf v}_h \in V_h^0$ and $q_h \in W_h$, ${\bf n}_1$ and ${\bf n}_2$ are the unit outward noramal in $\Omega_1$ and $\Omega_2$, respectively.
\end{algorithm}

\section{Stability}
For $T \in \mathcal{T}_h$, let $Q_h$ be the $L^2$ projection operator onto $V_h$. When $T \in \mathcal{T}_h^S$, $Q_h {\bf u}=\{Q_0 {\bf u}, Q_b {\bf u} \}$. When $T \in \mathcal{T}_h^I$, $Q_h= \{Q_0,\tilde{Q}_b\}$, where $e \in \mathcal{E}_h^S$, $\tilde{Q}_b {\bf u} = Q_b {\bf u}$; $e \in \mathcal{E}_h^I$, $\tilde{Q}_b {\bf u} = \hat{Q}_b ({\bf u} {\bf n}^{\mathrm{T}})$. To make sure there's only one value on the interface, let $\hat{Q}_b ({\bf u} {\bf n}^{\mathrm{T}})=\hat{Q}_b ({\bf u}_1 {\bf n}_1^{\mathrm{T}})$. Therefore, when $T \in \mathcal{T}_{2h}^I$, $\hat{Q}_b ({\bf u}_2 {\bf n}_2^{\mathrm{T}})=\hat{Q}_b (\bm{\psi} {\bf n}_1^{\mathrm{T}})-\hat{Q}_b ({\bf u}_1 {\bf n}_1^{\mathrm{T}})$. $\mathcal{Q}_h$ and $\mathbb{Q}_h$ are defined as the $L^2$ projection operator onto $P_{k-1}(T)$ and $[P_{k-1}(T)]^{2 \times 2}$, respectively.
\begin{lemma}\label{weak gradient exchange}
On each element $T \in \mathcal{T}_h$, for the discrete weak gradient operator, we have the following properties:
$\forall \tau \in [P_{k-1}(T)]^{2 \times 2}$,

(1)for $T\nsubseteq \Omega_2$ or $ \partial T \cap \Gamma = \emptyset$, we have that
\begin{eqnarray}\label{weak gradient exchange 1}
(\nabla_w(Q_h {\bf u}),\tau)_T=(\mathbb{Q}_h(\nabla {\bf u}),\tau)_T,
\end{eqnarray}

(2)for $T\subseteq \Omega_2 $ and $\partial T \cap \Gamma \neq \emptyset$, we drive
\begin{eqnarray}\label{weak gradient exchange 2}
(\nabla_w(Q_h {\bf u}),\tau)_T=(\mathbb{Q}_h(\nabla {\bf u}),\tau)_T-\langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}} ),\tau \rangle_{\partial T \cap \Gamma}.
\end{eqnarray}
\end{lemma}

\begin{proof}According to Eq.(\ref{weak gradient 1}), definition of $Q_0$, $Q_b$ and $\mathbb{Q}_h$ and the integration by parts, for T $\in \mathcal{T}_h^S $ and $\forall \tau \in [P_{k-1}(T)]^{2 \times 2}$, we have
\begin{align*}
(\nabla_w(Q_h {\bf u}), \tau)_T=&-(Q_0 {\bf u}, \nabla \cdot \tau)_T+\langle Q_b {\bf u}, \tau \cdot {\bf n} \rangle_{\partial T}\\
=&-({\bf u},\nabla \cdot \tau)_T+\langle {\bf u},\tau \cdot {\bf n} \rangle_{\partial T}\\
=&(\nabla {\bf u},\tau)_T\\
=&(\mathbb{Q}_h(\nabla {\bf u}),\tau)_T.
\end{align*}
Next, using Eq.(\ref{weak gradient 2}) can lead to, for T $\in \mathcal{T}_{1h}^I$
\begin{align*}
(\nabla_w(Q_h {\bf u}), \tau)_T=&-(Q_0 {\bf u}, \nabla \cdot \tau)_T+\langle Q_b {\bf u}, \tau {\bf n}_1 \rangle_{\partial T \backslash \Gamma}+
\langle \hat{Q}_b({\bf u} {\bf n}_1^{\mathrm{T}}),\tau \rangle_{\partial T \cap \Gamma}\\
=&-({\bf u},\nabla \cdot \tau)_T+\langle {\bf u} ,\tau {\bf n}_1 \rangle_{\partial T \backslash \Gamma}+\langle {\bf u} {\bf n}_1^{\mathrm{T}},\tau \rangle_{\partial T \cap \Gamma}\\
=&-({\bf u},\nabla \cdot \tau)_T+\langle {\bf u},\tau {\bf n}_1 \rangle_{\partial T}\\
=&(\nabla {\bf u},\tau)_T\\
=&(\mathbb{Q}_h(\nabla {\bf u}),\tau)_T.
\end{align*}
Combining the above two equations, Eq.(\ref{weak gradient exchange 1}) is completed. Similarly, using (\ref{weak gradient 3}) and interface condition, for the T $\in \mathcal{T}_{2h}^I$ and $\forall \tau \in [P_{k-1}(T)]^{2 \times 2}$, we have
\begin{align*}
(\nabla_w(Q_h {\bf u}), \tau)_T=&-(Q_0 {\bf u}_2, \nabla \cdot \tau)_T+\langle Q_b {\bf u}_2, \tau {\bf n}_2 \rangle_{\partial T \backslash \Gamma}-
\langle \hat{Q}_b({\bf u}_1 {\bf n}_1^{\mathrm{T}}),\tau \rangle_{\partial T \cap \Gamma}\\
=&-({\bf u}_2,\nabla \cdot \tau)_T+\langle {\bf u}_2 ,\tau {\bf n}_2 \rangle_{\partial T \backslash \Gamma}-\langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}}),\tau \rangle_{\partial T \cap \Gamma}\\
&+\langle \hat{Q}_b({\bf u}_2 {\bf n}_2^{\mathrm{T}}),\tau \rangle_{\partial T \cap \Gamma}\\
=&-({\bf u}_2,\nabla \cdot \tau)_T+\langle {\bf u}_2,\tau {\bf n}_2 \rangle_{\partial T}-\langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}}),\tau \rangle _{\partial T \cap \Gamma}\\
=&(\nabla {\bf u}_2,\tau)_T-\langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}}),\tau \rangle _{\partial T \cap \Gamma}\\
=&(\mathbb{Q}_h(\nabla {\bf u}_2),\tau)_T-\langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}}),\tau \rangle _{\partial T \cap \Gamma}.
\end{align*}
\end{proof}

\begin{lemma}\label{weak divergence exchange}
On each element $T \in \mathcal{T}_h$, for the weak divergence operator, we have the following properties:
$\forall \tau \in P_{k-1}(T)$,

(1)for $T\nsubseteq \Omega_2$ or $\partial T \cap \Gamma = \emptyset$, we obtain that
\begin{eqnarray}\label{weak divergence exchange 1}
(\nabla_w \cdot (Q_h {\bf u}),\tau)_T=(\mathcal{Q}_h(\nabla \cdot {\bf u}),\tau)_T,
\end{eqnarray}

(2)for $T\subseteq \Omega_2$ and $\partial T \cap \Gamma \neq \emptyset$, we drive that
\begin{eqnarray}\label{weak divergence exchange 2}
(\nabla_w \cdot (Q_h {\bf u}),\tau)_T=(\mathcal{Q}_h(\nabla \cdot {\bf u}),\tau)_T-\langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}} ),\tau I\rangle_{\partial T \cap \Gamma}.
\end{eqnarray}
\end{lemma}

\begin{proof}According to the (\ref{weak divergence 1}), definition of $Q_0$, $Q_b$ and $\mathcal{Q}_h$ and the integration by parts, for the T $\in \mathcal{T}_h^S$ we have
\begin{align*}
(\nabla_w \cdot (Q_h {\bf u}),\tau)_T=&-(Q_0 {\bf u}, \nabla \tau)_T+\langle Q_b {\bf u} ,\tau {\bf n} \rangle_{\partial T}\\
=&-({\bf u},\nabla \tau)_T+\langle {\bf u} ,\tau {\bf n} \rangle_{\partial T}\\
=&(\nabla \cdot {\bf u},\tau)_T\\
=&(\mathcal{Q}_h(\nabla \cdot {\bf u}),\tau)_T.
\end{align*}
Next, using (\ref{weak divergence 2}), for the T $\in \mathcal{T}_{1h}^I$,
\begin{align*}
(\nabla_w \cdot (Q_h {\bf u}),\tau)_T=&-(Q_0 {\bf u},\nabla \tau)_T+\langle Q_b {\bf u}, \tau {\bf n}_1
\rangle_{\partial T \backslash \Gamma}
+\langle \hat{Q}_b({\bf u} {\bf n}_1^{\mathrm{T}}),\tau I \rangle _{\partial T \cap \Gamma}\\
=&-({\bf u},\nabla \tau)_T+\langle {\bf u},\tau {\bf n}_1 \rangle_{\partial T \backslash \Gamma}+\langle {\bf u} {\bf n}_1^{\mathrm{T}},\tau I \rangle_{\partial T \cap \Gamma}\\
=&-({\bf u},\nabla \tau)_T+\langle {\bf u},\tau {\bf n}_1 \rangle_{\partial T}\\
=&(\nabla \cdot {\bf u},\tau)_T\\
=&(\mathcal{Q}_h(\nabla \cdot {\bf u}),\tau )_T
\end{align*}
Using (\ref{weak divergence 3}) and the interface conditions, for the T $\subseteq \mathcal{T}_{2h}$ and $\partial T \cap \Gamma \neq \emptyset$,
\begin{align*}
(\nabla_w \cdot (Q_h {\bf u}),\tau)_T=&-(Q_0 {\bf u}_2,\nabla \tau)_T+\langle Q_b {\bf u}_2, \tau {\bf n}_2
\rangle_{\partial T \backslash \Gamma}
-\langle \hat{Q}_b({\bf u}_1 {\bf n}_1^{\mathrm{T}}),\tau I \rangle _{\partial T \cap \Gamma}\\
=&-({\bf u}_2,\nabla \tau)_T+\langle {\bf u}_2,\tau {\bf n}_2 \rangle_{\partial T \backslash \Gamma}\\
&+\langle \hat{Q}_b({\bf u}_2 {\bf n}_2^{\mathrm{T}}),\tau I \rangle_{\partial T \cap \Gamma}-\langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}}),\tau I \rangle_{\partial T \cap \Gamma}\\
=&-({\bf u}_2,\nabla \tau)_T+\langle {\bf u}_2,\tau {\bf n}_2 \rangle_{\partial T}
-\langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}}),\tau I \rangle_{\partial T \cap \Gamma}\\
=&(\nabla \cdot {\bf u}_2,\tau)_T-\langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}}),\tau I \rangle_{\partial T \cap \Gamma}\\
=&(\mathcal{Q}_h(\nabla \cdot {\bf u}_2),\tau )_T-\langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}}),\tau I \rangle_{\partial T \cap \Gamma}
\end{align*}
\end{proof}

We now define a semi-norm in the weak Galerkin finite element $V_h$ as follows:
\begin{eqnarray}
{|\hspace{-.02in}|\hspace{-.02in}|}{\mathbf v}{|\hspace{-.02in}|\hspace{-.02in}|}^2=a({\mathbf v},{\mathbf v}).
\end{eqnarray}
Then, we have the following property.

\begin{lemma}
${|\hspace{-.02in}|\hspace{-.02in}|} \cdot {|\hspace{-.02in}|\hspace{-.02in}|}$ provides a norm in $V_h^0$.
\end{lemma}
\begin{proof}
It's easy to see that the positive property and linear property for ${|\hspace{-.02in}|\hspace{-.02in}|} \cdot {|\hspace{-.02in}|\hspace{-.02in}|}$. Next, assume that ${|\hspace{-.02in}|\hspace{-.02in}|} \cdot {|\hspace{-.02in}|\hspace{-.02in}|} =0$ for some ${\mathbf v} \in V_h^0$, then it follows that
\begin{align*}
0=&\sum_{T \in \mathcal{T}_h}(A\nabla_w {\mathbf v},\nabla_w {\mathbf v})_T+\sum_{T \in \mathcal{T}_h^S}h_T^{-1}\langle Q_b {\mathbf v}_0-{\mathbf v}_b,Q_b {\mathbf v}_0-{\mathbf v}_b \rangle_{\partial T}\\
+&\sum_{T \in \mathcal{T}_{1h}^I} h_T^{-1} \Big( \langle Q_b {\mathbf v}_0 -{\mathbf v}_b,Q_b {\mathbf v}_0 -{\mathbf v}_b \rangle_{\partial T \backslash \Gamma}\\
+&\langle \hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})-\tilde{{\mathbf v}}_b,\hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})-\tilde{{\mathbf v}}_b) \rangle_{\partial T \cap \Gamma} \Big)\\
+&\sum_{T \in \mathcal{T}_{2h}^I }h_T^{-1} \Big( \langle Q_b {\mathbf v}_0 -{\mathbf v}_b,Q_b {\mathbf v}_0 -{\mathbf v}_b \rangle_{\partial T \backslash \Gamma}\\
+&\langle \hat{Q}_b({\mathbf v}_0 {\bf n}_2^{\mathrm{T}})+\tilde{{\mathbf v}}_b,\hat{Q}_b({\mathbf v}_0 {\bf n}_2^{\mathrm{T}})+\tilde{{\mathbf v}}_b) \rangle_{\partial T \cap \Gamma} \Big).\\
\end{align*}
The above equation implies the following conclusions:
\begin{enumerate}
  \item $\nabla_w {\mathbf v} =0$ on all $T \in \mathcal{T}_h $;\\
  \item $Q_b {\mathbf v}_0={\mathbf v}_b $ for $ e \in \mathcal{E}_h^S$;\\
  \item $\hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})=\tilde{{\mathbf v}}_b $ for $ e \in \mathcal{E}_h^I $;\\
  \item $\hat{Q}_b({\mathbf v}_0 {\bf n}_2^{\mathrm{T}})=-\tilde{{\mathbf v}}_b $ for $ e \in \mathcal{E}_h^I $.
\end{enumerate}
Therefore, by the above fact, we have for any $q \in [P_{k-1}(T)]^{2 \times 2}$,

(1)For the T $\in \mathcal{T}_{h}^S$,
\begin{align*}
0=&(\nabla_w {\mathbf v} ,q)_T\\
=&-({\mathbf v}_0,\nabla \cdot q)_T+\langle {\mathbf v}_b ,q {\bf n} \rangle_{\partial T}\\
=&(\nabla {\mathbf v}_0,q)_T-\langle {\mathbf v}_0 ,q {\bf n} \rangle_{\partial T}+\langle {\mathbf v}_b ,q {\bf n} \rangle_{\partial T}\\
=&(\nabla {\mathbf v}_0,q)_T-\langle Q_b {\mathbf v}_0 -{\mathbf v}_b ,q {\bf n} \rangle_{\partial T}\\
=&(\nabla {\mathbf v}_0,q)_T
\end{align*}
In the above equation letting $q=\nabla {\mathbf v}_0$ can drive $\nabla {\mathbf v}_0=0, \ T \in \mathcal{T}_{h}^I$. \\

(2)For the T $\in \mathcal{T}_{1h}^I$,
\begin{align*}
0=&(\nabla_w {\mathbf v} ,q)_T\\
=&-({\mathbf v}_0,\nabla \cdot q)_T+\langle {\mathbf v}_b ,q {\bf n}_1 \rangle_{\partial T \backslash \Gamma}+\langle \tilde{{\mathbf v}}_b , q \rangle_{\partial T \cap \Gamma}\\
=&(\nabla {\mathbf v}_0,q)_T-\langle {\mathbf v}_0 ,q {\bf n}_1 \rangle_{\partial T}+\langle {\mathbf v}_b ,q {\bf n}_1 \rangle_{\partial T \backslash \Gamma}
+\langle \tilde{{\mathbf v}}_b, q \rangle_{\partial T \cap \Gamma}\\
=&(\nabla {\mathbf v}_0,q)_T-\langle Q_b {\mathbf v}_0 -{\mathbf v}_b ,q {\bf n}_1 \rangle_{\partial T \backslash \Gamma}-\langle
\hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})-\tilde{{\mathbf v}}_b,q \rangle_{\partial T \cap \Gamma}\\
=&(\nabla {\mathbf v}_0,q)_T
\end{align*}
In the above equation letting $q=\nabla {\mathbf v}_0$ can drive $\nabla {\mathbf v}_0=0, T \in \mathcal{T}_{1h}^I$. \\

(3)For the T $\in \mathcal{T}_{2h}^I$,
\begin{align*}
0=&(\nabla_w {\mathbf v} ,q)_T\\
=&-({\mathbf v}_0,\nabla \cdot q)_T+\langle {\mathbf v}_b ,q {\bf n}_2 \rangle_{\partial T \backslash \Gamma}-\langle \tilde{{\mathbf v}}_b , q \rangle\\
=&(\nabla {\mathbf v}_0,q)_T-\langle {\mathbf v}_0 ,q {\bf n}_2 \rangle_{\partial T}+\langle {\mathbf v}_b ,q {\bf n}_2 \rangle_{\partial T \backslash \Gamma}
-\langle \tilde{{\mathbf v}}_b , q \rangle_{\partial T \cap \Gamma}\\
=&(\nabla {\mathbf v}_0,q)_T-\langle Q_b {\mathbf v}_0 -{\mathbf v}_b ,q {\bf n}_2 \rangle_{\partial T \backslash \Gamma}-\langle \hat{Q}_b({\mathbf v}_0 {\bf n}_2^{\mathrm{T}})+\tilde{{\mathbf v}}_b,q\rangle_{\partial T \cap \Gamma}\\
=&(\nabla {\mathbf v}_0,q)_T
\end{align*}
In the above equation letting $q=\nabla {\mathbf v}_0$ can drive $\nabla {\mathbf v}_0=0, T \in \mathcal{T}_{2h}^I$. \\

It follow that ${\mathbf v}_0=constant$, on $T \in \mathcal{T}_{h} $. This together with the fact that ${\mathbf v}_b=0 $ on $\partial \Omega$ implies that ${\mathbf v}_0=0$ and ${\mathbf v}_b=0$.
\end{proof}

\begin{lemma}\label{Abounded}
For $\forall~ {\mathbf v},{\bf w} \in V_h^0$,
\begin{eqnarray}
|a({\mathbf v},{\bf w})|\leqslant {|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|} \cdot {|\hspace{-.02in}|\hspace{-.02in}|} {\bf w} {|\hspace{-.02in}|\hspace{-.02in}|};\\
a({\mathbf v},{\mathbf v})= {|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|}^2.
\end{eqnarray}
\end{lemma}
\begin{lemma}\cite{WGStokes1}
	Let $\mathcal{T}_h$ be a finite element partition of $\Omega$ satisfying the shape regularity assumption and ${\bf w} \in [H^{r+1}(\Omega)]^d$ and $\rho \in H^r(\Omega)$ with
	$1 \leqslant r \leqslant k$. Then, for $1 \leqslant s \leqslant 1$ we have
	\begin{eqnarray}
		\sum_{T \in \mathcal{T}_h} h_T^{2s} \| {\bf w}- Q_0 {\bf w}\|_{T,s}^2 \leqslant Ch^{2(r+1)}\|{\bf w}\|_{r+1}^2,\label{projectorestimate1}\\
		\sum_{T \in \mathcal{T}_h} h_T^{2s} \|\nabla {\bf w}- \mathbb{Q}_h(\nabla {\bf w}) \|_{T,s}^2 \leqslant C h^{2r}\|{\bf w}\|_{r+1}^2,\label{projectorestimate2}\\
		\sum_{T \in \mathcal{T}_h} h_T^{2s} \|\rho- Q_h^p \rho \|_{T,s}^2 \leqslant Ch^{2r}\|\rho\|_{r+1}^2 \label{projectorestimate3}.
	\end{eqnarray}
Here C denotes a generic constant independent of the mesh size h and the functions in the estimates.
\end{lemma}

\begin{lemma}(inf-sup condition)\label{InfSupCondition}
There exists a positive constant C such that\\
\begin{eqnarray}\label{InfSup}
\sup_{{\mathbf v} \in V_h^0}\frac{b({\mathbf v},\rho)}{{|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|}} \geqslant C \| \rho\|
\end{eqnarray}
for all $\rho \in W_h $, where C is independent of mesh size h.
\end{lemma}
\begin{proof}
According to \cite{bookfiniteelementmethod1,bookfiniteelementmethod2,bookfiniteelementmethod3,bookfiniteelementmethod5,bookfiniteelementmethod4,WGStokes1},for any given $\rho \in W_h\subset L_0^2(\Omega)$, there exists a vector-value function $\tilde{{\mathbf v}} \in [H_0^1(\Omega)]^d$ such that
\begin{align*}
\frac{(\nabla \cdot \tilde{{\mathbf v}},\rho) }{\| \tilde{{\mathbf v}} \|_1} \geqslant C \| \rho\|,
\end{align*}
where $C \geqslant 0$ is a constant depending only on the domain $\Omega$. Let ${\mathbf v} = Q_h \tilde{{\mathbf v}} \in V_h$, then we need to verify the following inequality holds true:
\begin{eqnarray}\label{InfSup1}
{|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|} \leqslant C_0 \| \tilde{{\mathbf v}} \|_1.
\end{eqnarray}
Firstly, it follows from the definition of ${|\hspace{-.02in}|\hspace{-.02in}|} \cdot {|\hspace{-.02in}|\hspace{-.02in}|}$ and Lemma \ref{weak gradient exchange} that for $T \in \mathcal{T}_h$,
\begin{eqnarray}\label{InfSup1.1}
A\| \nabla_w {\mathbf v} \|_T^2 =A\| \nabla_w (Q_h \tilde{{\mathbf v}}) \|_T^2 = A\|\mathbb{Q}_h(\nabla \tilde{{\mathbf v}}) \|_T^2 \leqslant A\| \nabla \tilde{{\mathbf v}} \|_T^2 \leqslant C
\| \tilde{{\mathbf v}} \|_{1,T}^2.
\end{eqnarray}
Next, by using the definition of $Q_h$, the trace inequality and the estimate (\ref{projectorestimate1}), we can derive that for $T \in \mathcal{T}_h^S $,
\begin{eqnarray}\label{SvvS}
\begin{split}
s({\bf u},{\mathbf v})|_T=&h_T^{-1}\|Q_b {\mathbf v}_0 -{\mathbf v}_b \|_{\partial T}^{2}\\
=&h_T^{-1}\|Q_b (Q_0 \tilde{{\mathbf v}}) -Q_b{\tilde{{\mathbf v}}} \|_{\partial T}^{2}\\
\leqslant &h_T^{-1}\|Q_0 \tilde{{\mathbf v}}-\tilde{{\mathbf v}}\|_{\partial T}^{2}\\
\leqslant &C\Big( h_T^{-2}\|Q_0 \tilde{{\mathbf v}}-\tilde{{\mathbf v}}\|_T^2+\|\nabla(Q_0 \tilde{{\mathbf v}}-\tilde{{\mathbf v}}) \|_T^2 \Big)\\
\leqslant &C\|\tilde{{\mathbf v}}\|_{1,T}^2.
\end{split}
\end{eqnarray}
When the element T $\in \mathcal{T}_{1h}^I$, the stablizer term has the following definition
\begin{align*}
s({\bf u},{\mathbf v})|_T=h_T^{-1}\|Q_b{{\mathbf v}_0}-{\mathbf v}_b \|_{\partial T \backslash \Gamma}+h_T^{-1}\|\hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})-\tilde{{\mathbf v}}_b \|_{\partial T \cap \Gamma}.
\end{align*}
The first term is similar to (\ref{SvvS}), so we consider the second term. Proceeding as in the proof of (\ref{SvvS}), we have
\begin{eqnarray}
\begin{split}
&h_T^{-1}\|\hat{Q_b}({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})-\tilde{{\mathbf v}}_b \|_{\partial T \cap \Gamma}\\
=&h_T^{-1}\|\hat{Q_b}((Q_0 \tilde{{\mathbf v}}) {\bf n}_1^{\mathrm{T}})-\hat{Q}_b(\tilde{{\mathbf v}} {\bf n}_1^{\mathrm{T}})\|_{\partial T \cap \Gamma}\\
\leqslant & h_T^{-1}\|(Q_0 \tilde{{\mathbf v}}) {\bf n}_1^{\mathrm{T}}-\tilde{{\mathbf v}} {\bf n}_1^{\mathrm{T}} \|_{\partial T \cap \Gamma}\\
\leqslant & h_T^{-1}\|Q_0 \tilde{{\mathbf v}}- \tilde{{\mathbf v}}\|_{\partial T \cap \Gamma}\\
\leqslant & h_T^{-1}\|Q_0 \tilde{{\mathbf v}}- \tilde{{\mathbf v}}\|_{\partial T}\\
\leqslant &C\| \tilde{{\mathbf v}}\|_{1,T}^2.
\end{split}
\end{eqnarray}
Similarly, for $T \in \mathcal{T}_{2h}^I$, we have
\begin{align*}
s({\bf u},{\mathbf v})|_T=&h_T^{-1}\| Q_b {\mathbf v}_0- {\mathbf v}_b|_{\partial T \backslash \Gamma}+h_T^{-1} \| \hat{Q}_b({\mathbf v}_0 {\bf n}_2^{\mathrm{T}})+\tilde{{\mathbf v}}_b\|_{\partial T \cap \Gamma},
\end{align*}
then using the same estimate method as (\ref{SvvS}), we obtain that
\begin{eqnarray}\label{InfSup1.2}
\begin{split}
&h_T^{-1}\|\hat{Q_b}({\mathbf v}_0 {\bf n}_2^{\mathrm{T}})+\tilde{{\mathbf v}}_b \|_{\partial T \cap \Gamma}\\
=&h_T^{-1}\|\hat{Q_b}((Q_0 \tilde{{\mathbf v}}) {\bf n}_2^{\mathrm{T}})+\hat{Q}_b(\tilde{{\mathbf v}} {\bf n}_1^{\mathrm{T}})\|_{\partial T \cap \Gamma}\\
\leqslant & h_T^{-1}\|(Q_0 \tilde{{\mathbf v}}) {\bf n}_2^{\mathrm{T}}+\tilde{{\mathbf v}} {\bf n}_1^{\mathrm{T}} \|_{\partial T \cap \Gamma}\\
\leqslant & h_T^{-1}\|Q_0 \tilde{{\mathbf v}}- \tilde{{\mathbf v}}\|_{\partial T \cap \Gamma}\\
\leqslant & h_T^{-1}\|Q_0 \tilde{{\mathbf v}}- \tilde{{\mathbf v}}\|_{\partial T}\\
\leqslant &C\| \tilde{{\mathbf v}}\|_{1,T}^2
\end{split}
\end{eqnarray}
Consequently, Eq.(\ref{InfSup1}) follows by using (\ref{InfSup1.1})-(\ref{InfSup1.2}).
 Another step in the proof is to show $b({\mathbf v},\rho)$, on the basis of Lemma \ref{weak divergence exchange} and the definition of $\mathcal{Q}_h$, it's easy to see that
\begin{eqnarray}\label{InfSup2}
(\nabla \cdot (Q_h \tilde{{\mathbf v}}),\rho)_T=(\mathcal{Q}_h(\nabla \cdot \tilde{{\mathbf v}}),\rho)_T=(\nabla \cdot \tilde{{\mathbf v}},\rho)_T, ~ \forall \rho \in W_h
\end{eqnarray}
for all $T \in \mathcal{T}_h$. Combining (\ref{InfSup1}) with (\ref{InfSup2}), we have
$$
\frac{b({\mathbf v},p)}{{|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|}} \geqslant \frac{\sum_{T \in \mathcal{T}_h}(\nabla \cdot \tilde{{\mathbf v}},\rho)_T}{C_0  \|\tilde{{\mathbf v}} \|_1}
\geqslant \frac{\beta}{C_0} \| \rho \| =C \| \rho \| .
$$
This indicates that (\ref{InfSup}) holds true.
\end{proof}

According to the Lemma \ref{Abounded} and Lemma \ref{InfSupCondition}, it follows that the weak Galerkin finite element scheme (\ref{algorithm 1})-(\ref{algorithm 2}) exist a unique solution and the solution is stable.

\section{Error Equation}In the section, we will give the equation of the error of ${\bf u}$ and $p$. Firstly, we assume $A$ is a piecewise constant in the domain $\Omega$. We use $\{{\bf u}_h; p_h\}$ to denote the solution of (\ref{algorithm 1})-(\ref{algorithm 2}). Let $\{{\bf u},p\}$ be the exact solution of (\ref{model 1})-(\ref{interface 2}). Denote by $Q_h {\bf u} \in V_h$ the $L^2$ projection of ${\bf u}$. Let $Q_h^p$ be the $L^2$ projection of $p$ onto $W_h$. The error of ${\bf u} ~\text{and}\ p$ is given by
\begin{eqnarray}
e_h=Q_h {\bf u} - {\bf u}_h, ~ \varepsilon_h=Q_h^p p-p_h.
\end{eqnarray}

\begin{lemma}\label{EE}
For the $({\bf u}_i; p_i)\in [H^1(\Omega_i)]^2 \times H^1(\Omega_i)$ with $i=1,2$ satisfying the equations (\ref{model 1})-(\ref{interface 2}), we can derive the following equation
\begin{eqnarray}\label{EEE}
\begin{split}
&\sum_{T \in \mathcal{T}_h}(\nabla_w(Q_h {\bf u}), A\nabla_w {\mathbf v})_T-\sum_{T \in \mathcal{T}_h}(\nabla_w \cdot {\mathbf v}, Q_h^p p)_T \\
=&({\bf f},{\mathbf v}_0)+\ell_1({\bf u},{\mathbf v})+\tilde{\ell}_1({\bf u},{\mathbf v})-\ell_2(p,{\mathbf v})-\tilde{\ell}_2(p,{\mathbf v})\\
+&\sum_{e \in \mathcal{E}_h^I}
\langle \tilde{{\mathbf v}}_b,\bm{\phi}{\bf n}_1^{\mathrm{T}} \rangle_e - \sum_{e \in \mathcal{E}_h^I} \langle \hat{Q}_b(\bm{\psi}{\bf n}_1^{\mathrm{T}},
A_2 \nabla_w {\mathbf v} \rangle_e, \\
\end{split}
\end{eqnarray}
where
\begin{eqnarray}
\ell_1({\bf u},{\mathbf v})=\sum_{T \in \mathcal{T}_h^S} \langle {\mathbf v}_0 -{\mathbf v}_b, A \nabla {\bf u}-A \mathbb{Q}_h(\nabla {\bf u}) \rangle_{\partial T},
\end{eqnarray}
\begin{eqnarray}
\begin{split}
\tilde{\ell}_1({\bf u},{\mathbf v})=\sum_{T \in \mathcal{T}_{1h}^I}
\Big( \langle {\mathbf v}_0 -{\mathbf v}_b, A_1 \nabla {\bf u}-A_1 \mathbb{Q}_h(\nabla {\bf u}_1)   \rangle_{\partial T \backslash \Gamma}\\
+\langle {\mathbf v}_0 {\bf n}_1^{\mathrm{T}} -\tilde{{\mathbf v}}_b, A_1 \nabla {\bf u}_1-A_1 \mathbb{Q}_h(\nabla {\bf u}_1) \rangle_{\partial T \cap \Gamma} \Big)\\
+\sum_{T \in \mathcal{T}_{2h}^I}
\Big( \langle {\mathbf v}_0 -{\mathbf v}_b, A_2 \nabla {\bf u}-A_2 \mathbb{Q}_h(\nabla {\bf u}_2)  \rangle_{\partial T \backslash \Gamma}\\
+\langle {\mathbf v}_0 {\bf n}_2^{\mathrm{T}}+\tilde{{\mathbf v}}_b, A_2 \nabla {\bf u}_2-A_2 \mathbb{Q}_h(\nabla {\bf u}_2) \rangle_{\partial T \cap \Gamma} \Big),
\end{split}
\end{eqnarray}
\begin{eqnarray}
\ell_2(p,{\mathbf v})=\sum_{T \in \mathcal{T}_h^S} \langle {\mathbf v}_0 -{\mathbf v}_b, p{\bf n}-(Q_h^p p) {\bf n} \rangle_e,
\end{eqnarray}
\begin{eqnarray}
\begin{split}
\tilde{\ell}_2(p,{\mathbf v})=\sum_{T \in \mathcal{T}_{1h}^I}
\Big( \langle {\mathbf v}_0 -{\mathbf v}_b,p_1{\bf n}_1-(Q_h^p p_1) {\bf n}_1 \rangle_{\partial T \backslash \Gamma}
+ \langle {\mathbf v}_0 {\bf n}_1^{\mathrm{T}} -\tilde{{\mathbf v}}_b,p_1 I-(Q_h^p p_1) I \rangle_{\partial T \cap \Gamma} \Big)\\
+\sum_{T \in \mathcal{T}_{2h}^I} \Big(
\langle {\mathbf v}_0 -{\mathbf v}_b,p_2{\bf n}_2-(Q_h^p p_2) {\bf n}_2 \rangle_{\partial T \backslash \Gamma}
+\langle {\mathbf v}_0 {\bf n}_2^{\mathrm{T}}+\tilde{{\mathbf v}}_b,p_2I- (Q_h^p p_2) I \rangle_{\partial T \cap \Gamma} \Big).
\end{split}
\end{eqnarray}

\end{lemma}

\begin{proof}
According to the (\ref{weak gradient 1}), the definition of the weak gradient and the integration by parts, it follows that for $T \in \mathcal{T}_h^S$,
\begin{eqnarray}\label{EEP1}
\begin{split}
&(\nabla_w (Q_h {\bf u}),A \nabla_w {\mathbf v})_T\\
=&(\mathbb{Q}_h (\nabla {\bf u}),A \nabla_w {\mathbf v})_T\\
=&-({\mathbf v}_0,\nabla \cdot (A \mathbb{Q}_h (\nabla {\bf u})))_T+\langle {\mathbf v}_b,A \mathbb{Q}_h (\nabla {\bf u}) {\bf n} \rangle_{\partial T}\\
=&(A \mathbb{Q}_h (\nabla {\bf u}),\nabla {\mathbf v}_0)_T-\langle{\mathbf v}_0-{\mathbf v}_b,A \mathbb{Q}_h (\nabla {\bf u}) {\bf n}\rangle_{\partial T}\\
=&A(\nabla {\bf u},\nabla {\mathbf v}_0)_T-\langle{\mathbf v}_0-{\mathbf v}_b,A \mathbb{Q}_h(\nabla {\bf u}) {\bf n}\rangle_{\partial T}.\\
\end{split}
\end{eqnarray}
Similarly, by using the definition of the weak divergence and the integration by parts, we have
\begin{eqnarray}\label{EEP2}
\begin{split}
&(\nabla_w \cdot {\mathbf v}, Q_h^p p )_T\\
=&-({\mathbf v}_0,\nabla(Q_h^p p ))_T+\langle {\mathbf v}_b \cdot {\bf n},Q_h^p p \rangle_{\partial T}\\
=&(\nabla \cdot {\mathbf v}_0,Q_h^p p)_T-\langle {\mathbf v}_0 -{\mathbf v}_b, (Q_h^p p){\bf n}\rangle_{\partial T}\\
=&(\nabla \cdot {\mathbf v}_0,p)_T-\langle {\mathbf v}_0 -{\mathbf v}_b, (Q_h^p p){\bf n}\rangle_{\partial T}.
\end{split}
\end{eqnarray}
For $T \in \mathcal{T}_{1h}^I$, it follows from (\ref{weak gradient 2}), (\ref{weak gradient exchange 1}) that
\begin{eqnarray}\label{EEP3}
\begin{split}
&(\nabla_w (Q_h {\bf u}),A \nabla_w {\mathbf v})_T\\
=&(\mathbb{Q}_h (\nabla {\bf u}),A \nabla_w {\mathbf v})_T\\
=&-({\mathbf v}_0,\nabla \cdot (A \mathbb{Q}_h (\nabla {\bf u})))_T+\langle {\mathbf v}_b,A \mathbb{Q}_h (\nabla {\bf u}) {\bf n}_1 \rangle_{\partial T \setminus \Gamma}+\langle \tilde{{\mathbf v}}_b,A \mathbb{Q}_h(\nabla {\bf u}) \rangle_{\partial T \cap \Gamma}\\
=&(A \mathbb{Q}_h(\nabla {\bf u}), \nabla {\mathbf v}_0)_T-\langle {\mathbf v}_0,A \mathbb{Q}_h(\nabla {\bf u}) {\bf n}_1 \rangle_{\partial T}\\
+&\langle {\mathbf v}_b,A \mathbb{Q}_h(\nabla {\bf u}) {\bf n}_1 \rangle_{\partial T \setminus \Gamma}+\langle \tilde{{\mathbf v}}_b,A \mathbb{Q}_h(\nabla {\bf u})\rangle_{\partial T \cap \Gamma}\\
=&A( \nabla {\bf u},\nabla {\mathbf v}_0)_T-\langle {\mathbf v}_0 - {\mathbf v}_b,A \mathbb{Q}_h(\nabla {\bf u}) {\bf n}_1 \rangle_{\partial T \setminus \Gamma}
-\langle {\mathbf v}_0 {\bf n}_1^{\mathrm{T}}-\tilde{{\mathbf v}}_b,A \mathbb{Q}_h(\nabla {\bf u}) \rangle_{\partial T \cap \Gamma}.
\end{split}
\end{eqnarray}
From (\ref{weak divergence 2}), (\ref{weak divergence exchange 1}), we have
\begin{eqnarray}\label{EEP4}
\begin{split}
(\nabla_w \cdot {\mathbf v}, Q_h^p p_1 )_T=&-({\mathbf v}_0,\nabla(Q_h^p p_1 ))_T+\langle {\mathbf v}_b \cdot {\bf n}_1,Q_h^p p_1 \rangle_{\partial T \setminus \Gamma}\\
+&\langle \tilde{{\mathbf v}}_b,(Q_h^p p_1)I \rangle_{\partial T \cap \Gamma}\\
=&(\nabla \cdot {\mathbf v}_0,Q_h^p p_1)_T-\langle {\mathbf v}_0 -{\mathbf v}_b,
(Q_h^p p_1) {\bf n}_1 \rangle_{\partial T \setminus \Gamma}\\
-&\langle {\mathbf v}_0 {\bf n}_1^{\mathrm{T}}-\tilde{{\mathbf v}}_b,(Q_h^p p_1)I \rangle_{\partial T \cap \Gamma}\\
=&(\nabla \cdot {\mathbf v}_0,p_1)_T-\langle {\mathbf v}_0 -{\mathbf v}_b, (Q_h^p p_1) {\bf n}_1 \rangle_{\partial T \setminus \Gamma}\\
-&\langle {\mathbf v}_0 {\bf n}_1^{\mathrm{T}}-\tilde{{\mathbf v}}_b,(Q_h^p p_1)I \rangle_{\partial T \cap \Gamma}.
\end{split}
\end{eqnarray}
Similarly, by using Eq.(\ref{weak gradient 3}) and (\ref{weak gradient exchange 2}), we have for
$T \in \mathcal{T}_{2h}^I$,
\begin{eqnarray}\label{EEP5}
\begin{split}
&(\nabla_w (Q_h {\bf u}),A \nabla_w {\mathbf v})_T\\
=&(\mathbb{Q}_h (\nabla {\bf u}),A_2 \nabla_w {\mathbf v})_T+\langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}}),A_2 \nabla_w {\mathbf v} \rangle_{\partial T \cap \Gamma}\\
=&-({\mathbf v}_0,\nabla \cdot (A_2 \mathbb{Q}_h (\nabla {\bf u})))_T+\langle {\mathbf v}_b,A_2 \mathbb{Q}_h (\nabla {\bf u}) {\bf n} \rangle_{\partial T \setminus \Gamma}\\
-&\langle \tilde{{\mathbf v}}_b,A_2 \mathbb{Q}_h(\nabla {\bf u}) \rangle_{\partial T \cap \Gamma}+\langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}}),A_2 \nabla_w {\mathbf v} \rangle_{\partial T \cap \Gamma}\\
=&A_2 ( \nabla {\bf u},\nabla {\mathbf v}_0)_T-\langle {\mathbf v}_0-{\mathbf v}_b,A_2 \mathbb{Q}_h(\nabla {\bf u}) {\bf n} \rangle_{\partial T \setminus \Gamma}\\
-&\langle {\mathbf v}_0 {\bf n}_2^{\mathrm{T}}+\tilde{{\mathbf v}}_b, A_2 \mathbb{Q}_h (\nabla {\bf u})\rangle_{\partial T \cap \Gamma}+\langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}}),A_2 \nabla_w {\mathbf v} \rangle_{\partial T \cap \Gamma},
\end{split}
\end{eqnarray}
\begin{eqnarray}\label{EEP6}
\begin{split}
(\nabla_w \cdot {\mathbf v}, Q_h^p p_2 )_T=&-({\mathbf v}_0,\nabla(Q_h^p p_2 ))_T+\langle {\mathbf v}_b \cdot {\bf n}_2,Q_h^p p_2 \rangle_{\partial T \setminus \Gamma}\\
-&\langle \tilde{{\mathbf v}}_b,(Q_h^p p_2)I \rangle_{\partial T \cap \Gamma}\\
=&(\nabla \cdot {\mathbf v}_0,Q_h^p p_2)_T-\langle {\mathbf v}_0 -{\mathbf v}_b, (Q_h^p p_2) {\bf n}_2 \rangle_{\partial T \setminus \Gamma}\\
-&\langle {\mathbf v}_0 {\bf n}_2^{\mathrm{T}}+ \tilde{{\mathbf v}}_b,(Q_h^p p_2)I \rangle_{\partial T \cap \Gamma}\\
=&(\nabla \cdot {\mathbf v}_0,p_2)_T-\langle {\mathbf v}_0 -{\mathbf v}_b, (Q_h^p p_2) {\bf n}_2 \rangle_{\partial T \setminus \Gamma}\\
-&\langle {\mathbf v}_0 {\bf n}_2^{\mathrm{T}}+\tilde{{\mathbf v}}_b,(Q_h^p p_2)I \rangle_{\partial T \cap \Gamma}.
\end{split}
\end{eqnarray}
Combining (\ref{EEP1}), (\ref{EEP3}) with (\ref{EEP5}), we have
\begin{eqnarray}\label{EEP7}
\begin{split}
&\sum_{T \in \mathcal{T}_h}(\nabla_w (Q_h {\bf u}),A \nabla_w {\mathbf v})_T\\
=&\sum_{T \in \mathcal{T}_h}(A \nabla {\bf u}, \nabla {\mathbf v})_T -\sum_{T \in \mathcal{T}_{h}^S} \langle {\mathbf v}_0 -{\mathbf v}_b, A \mathbb{Q}_h(\nabla {\bf u})  {\bf n} \rangle_{\partial T}\\
-&\sum_{T \in \mathcal{T}_{1h}^I} \Big( \langle {\mathbf v}_0 -{\mathbf v}_b,A_1 \mathbb{Q}_h(\nabla {\bf u}) {\bf n}_1 \rangle_{\partial T \setminus \Gamma}+\langle {\mathbf v}_0 {\bf n}_1^{\mathrm{T}}-\tilde{{\mathbf v}}_b,A_1 \mathbb{Q}_h(\nabla {\bf u})\rangle_{\partial T \cap \Gamma} \Big)\\
-&\sum_{T \in \mathcal{T}_{2h}^I} \Big( \langle {\mathbf v}_0 -{\mathbf v}_b,A_2\mathbb{Q}_h(\nabla {\bf u}) {\bf n}_2 \rangle_{\partial T \setminus \Gamma}+\langle {\mathbf v}_0 {\bf n}_2^{\mathrm{T}}+\tilde{{\mathbf v}}_b,A_2\mathbb{Q}_h(\nabla {\bf u})\rangle_{\partial T \cap \Gamma} \\
+&\langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}}),A_2 \nabla_w {\mathbf v} \rangle_{\partial T \cap \Gamma}\Big),
\end{split}
\end{eqnarray}
It follows from Eqs.(\ref{EEP2}),(\ref{EEP4}) and (\ref{EEP6}) that
\begin{eqnarray}\label{EEP8}
\begin{split}
&(\nabla_w \cdot {\mathbf v},Q_h^p p)\\
=&\sum_{T \in \mathcal{T}_h}\langle \nabla \cdot {\mathbf v}_0, p \rangle_T - \sum_{T \in \mathcal{T}_{h}^S} \langle {\mathbf v}_0-{\mathbf v}_b, Q_h^p p {\bf n} \rangle_{\partial T}\\
-&\sum_{T \in \mathcal{T}_{1h}^I}\Big( \langle {\mathbf v}_0 -{\mathbf v}_b, (Q_h^p p_1)  {\bf n}_1 \rangle_{\partial T \setminus \Gamma} +\langle {\mathbf v}_0 {\bf n}_1^{\mathrm{T}}-\tilde{{\mathbf v}}_b, (Q_h^p p_1) I\rangle_{\partial T \cap \Gamma} \Big)\\
-&\sum_{T \in \mathcal{T}_{2h}^I}\Big( \langle {\mathbf v}_0 -{\mathbf v}_b, (Q_h^p p_2)  {\bf n}_2 \rangle_{\partial T \setminus \Gamma} +\langle {\mathbf v}_0 {\bf n}_2^{\mathrm{T}}+\tilde{{\mathbf v}}_b, (Q_h^p p_2) I\rangle_{\partial T \cap \Gamma} \Big),\\
\end{split}
\end{eqnarray}
The next thing in the proof to test (\ref{model 1}) by using ${\mathbf v}_0$ in ${\mathbf v}=\{ {\mathbf v}_0, {\mathbf v}_b\} \in V_h^0$ to obtain
\begin{eqnarray}\label{EEP9}
-(\nabla \cdot (A \nabla {\bf u}),{\mathbf v}_0)+(\nabla p, {\mathbf v}_0)=({\bf f}, {\mathbf v}_0).
\end{eqnarray}
According to the integration by parts and the interface condition (\ref{interface 2}), it follows that
\begin{eqnarray}\label{EEP10}
\begin{split}
-(\nabla \cdot (A \nabla {\bf u}),{\mathbf v}_0)=&\sum_{T \in \mathcal{T}_h} (-\nabla \cdot (A \nabla {\bf u}),{\mathbf v}_0)_T\\
=&\sum_{T \in \mathcal{T}_h}(A\nabla {\bf u}, \nabla {\mathbf v}_0)-\sum_{T \in \mathcal{T}_h} \langle A\nabla {\bf u} \cdot {\bf n} ,{\mathbf v}_0 \rangle_{\partial T}\\
=&\sum_{T \in \mathcal{T}_h}(A\nabla {\bf u}, \nabla {\mathbf v}_0)_T-\sum_{e \in \mathcal{E}_h^S} \langle {\mathbf v}_0 - {\mathbf v}_b ,A \nabla {\bf u} \cdot {\bf n} \rangle_e\\
-&\sum_{e \in \mathcal{E}_h^I} \Big( \langle {\mathbf v}_0 {\bf n}_1^{\mathrm{T}}-\tilde{{\mathbf v}}_b,A_1 \nabla {\bf u}_1  \rangle_e
+\langle {\mathbf v}_0 {\bf n}_2^{\mathrm{T}}+\tilde{{\mathbf v}}_b,A_2 \nabla {\bf u}_2 \rangle_e \Big)\\
-&\sum_{e \in \Gamma} \langle \tilde{{\mathbf v}}_b,\bm{\phi}{\bf n}_1^{\mathrm{T}}+p_1I-p_2I \rangle_e,\\
\end{split}
\end{eqnarray}
where we have used the fact that $\sum_{e \in \mathcal{E}_h^S} \langle {\mathbf v}_b,\nabla {\bf u} \cdot {\bf n}\rangle_e=0$. Then by the fact that $\sum_{e \in \mathcal{E}_h^S} \langle {\mathbf v}_b, p {\bf n} \rangle_e =0$, it follows that
\begin{eqnarray}\label{EEP11}
\begin{split}
(\nabla p, {\mathbf v}_0)=&\sum_{T \in \mathcal{T}_h}(\nabla p, {\mathbf v}_0)_T\\
=&\sum_{T \in \mathcal{T}_h}-(\nabla \cdot {\mathbf v}_0,p)_T+\sum_{T \in \mathcal{T}_h} \langle p, {\mathbf v}_0 \cdot {\bf n} \rangle_{\partial T}\\
=&\sum_{T \in \mathcal{T}_h}-(\nabla \cdot {\mathbf v}_0,p)_T+\sum_{e \in \mathcal{E}_h^S} \langle {\mathbf v}_0-{\mathbf v}_b, p{\bf n} \rangle_e\\
+&\sum_{T \in \mathcal{E}_h^I} \Big( \langle {\mathbf v}_0 {\bf n}_1^{\mathrm{T}}-\tilde{{\mathbf v}}_b,p_1 I \rangle_e+\langle {\mathbf v}_0 {\bf n}_2^{\mathrm{T}}+\tilde{{\mathbf v}}_b,p_2 I \rangle_e \Big)\\
+&\sum_{e \in \mathcal{E}_h^I}\langle \tilde{{\mathbf v}}_b,p_1 I- p_2 I \rangle_e.
\end{split}
\end{eqnarray}
By substituting (\ref{EEP10}) into (\ref{EEP7}), substituting (\ref{EEP11}) into (\ref{EEP8}), and adding (\ref{EEP8}) to (\ref{EEP7}), the proof of Eq.(\ref{EEE}) is completed.
\end{proof}

\begin{lemma}\label{error equation}
The errors ${\bf e}_h$ and $\varepsilon_h$ of the weak Galerkin finite element solution  satisfy the following equations:
\begin{eqnarray}\label{error equation 1}
a({\bf e}_h,{\mathbf v})-b({\mathbf v},\varepsilon_h)&=&\ell_1({\bf u},{\mathbf v})+\tilde{\ell}_1({\bf u},{\mathbf v})-\ell_2(p,{\mathbf v})-\tilde{\ell}_2(p,{\mathbf v})+\tilde{s}(Q_h {\bf u} ,{\mathbf v})\\
b({\bf e}_h,q)&=&0 \label{error equation 2}
\end{eqnarray}
for all ${\mathbf v} \in V_h^0$ and $q \in W_h$.
\end{lemma}
\begin{proof}
Since the solution $({\bf u}; p)$ satisfies the Eqs.(\ref{model 1})-(\ref{interface 2}), according to Lemma \ref{EE} we have
\begin{align*}
&\sum_{T \in \mathcal{T}_h}(\nabla_w (Q_h {\bf u}),A \nabla_w {\mathbf v})_T-\sum_{T \in \mathcal{T}_h}(\nabla_w \cdot {\mathbf v}, Q_h^p p)_T\\
=&({\bf f},{\mathbf v}_0)+\ell_1({\bf u},{\mathbf v})+\tilde{\ell}_1({\bf u},{\mathbf v})-\ell_2(p,{\mathbf v})-\tilde{\ell}_2(p,{\mathbf v})\\
+&\sum_{e \in \mathcal{E}_h^I}\langle \tilde{{\mathbf v}}_b,\bm{\phi}{\bf n}_1^{\mathrm{T}} \rangle_e -\sum_{e \in \mathcal{E}_h^I}\langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}}), A_2 \nabla_w {\mathbf v} \rangle_e
\end{align*}
Adding $S(Q_h {\bf u},{\mathbf v})$ to two side of the above equation yields
\begin{eqnarray}\label{eep1}
\begin{split}
&a(Q_h{\bf u},{\mathbf v})-b({\mathbf v},Q_h^p p)\\
=&({\bf f},{\mathbf v}_0)+\ell_1({\bf u},{\mathbf v})+\tilde{\ell}_1({\bf u},{\mathbf v})-\ell_2(p,{\mathbf v})-\tilde{\ell}_2(p,{\mathbf v})+s(Q_h {\bf u},{\mathbf v})\\
+&\sum_{e \in \mathcal{E}_h^I}\langle \tilde{{\mathbf v}}_b,\bm{\phi}{\bf n}_1^{\mathrm{T}} \rangle_e -\sum_{e \in \mathcal{E}_h^I}\langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}}), A_2 \nabla_w {\mathbf v} \rangle_e
\end{split}
\end{eqnarray}
For $s(Q_h {\bf u},{\mathbf v})$, we need to pay special attention to the elements $T \in \mathcal{T}_{2h}^I$. According to the definition of $Q_h$, we have
\begin{eqnarray}\label{eep2}
\begin{split}
s(Q_h {\bf u},{\mathbf v})|_T=&h_T^{-1}\langle Q_b(Q_0 {\bf u})-Q_b {\bf u},Q_b {\mathbf v}_0 -{\mathbf v}_b \rangle_{\partial T \setminus \Gamma}\\
+&h_T^{-1}\langle \hat{Q}_b((Q_0 {\bf u}_2){\bf n}_2^{\mathrm{T}})+\hat{Q}_b({\bf u}_1 {\bf n}_1^{\mathrm{T}}),
\hat{Q}_b({\mathbf v}_0 {\bf n}_2^{\mathrm{T}})+\tilde{{\mathbf v}}_b \rangle_{\partial T \cap \Gamma}\\
=&h_T^{-1}\langle Q_b(Q_0 {\bf u})-Q_b {\bf u},Q_b {\mathbf v}_0 -{\mathbf v}_b \rangle_{\partial T \setminus \Gamma}\\
+&h_T^{-1}\langle \hat{Q}_b((Q_0 {\bf u}_2){\bf n}_2^{\mathrm{T}})-\hat{Q}_b({\bf u}_2 {\bf n}_2^{\mathrm{T}}),
\hat{Q}_b({\mathbf v}_0 {\bf n}_2^{\mathrm{T}})+\tilde{{\mathbf v}}_b \rangle_{\partial T \cap \Gamma}\\
+&h_T^{-1}\langle \tilde{Q}_b(\bm{\psi}{\bf n}_1^{\mathrm{T}}),\hat{Q}_b({\mathbf v}_0 {\bf n}_2^{\mathrm{T}})+\tilde{{\mathbf v}}_b
\rangle_{\partial T \cap \Gamma}
\end{split}
\end{eqnarray}
Therefore, denote
\begin{eqnarray}\label{eep3}
\begin{split}
\tilde{s}(Q_h {\bf u},{\mathbf v})=&\sum_{T \in \mathcal{T}_h^S}h_T^{-1} \langle Q_b(Q_0 {\bf u})-Q_b {\bf u},
Q_b {\mathbf v}_0 -{\mathbf v}_b \rangle_{\partial T}\\
+&\sum_{ T \in \mathcal{T}_{1h}^I} h_T^{-1}\Big( \langle Q_b(Q_0 {\bf u})-Q_b {\bf u}, Q_b {\mathbf v}_0 -{\mathbf v}_b\rangle_{\partial T \setminus \Gamma} \\
+ &\langle \hat{Q}_b((Q_0 {\bf u}_1){\bf n}_1)-\hat{Q}_b({\bf u}_1 {\bf n}_1^{\mathrm{T}}),
\hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})-\tilde{{\mathbf v}}_b \rangle_{\partial T \cap \Gamma} \Big)\\
+&\sum_{T \in \mathcal{T}_{2h}^I} h_T^{-1}\Big( \langle Q_b(Q_0 {\bf u})-Q_b {\bf u}, Q_b {\mathbf v}_0 -{\mathbf v}_b\rangle_{\partial T \setminus \Gamma} \\
+ &\langle \hat{Q}_b((Q_0 {\bf u}_2){\bf n}_2)-\hat{Q}_b({\bf u}_2 {\bf n}_2^{\mathrm{T}}),
\hat{Q}_b({\mathbf v}_0 {\bf n}_2^{\mathrm{T}})+\tilde{{\mathbf v}}_b \rangle_{\partial T \cap \Gamma} \Big).\\
\end{split}
\end{eqnarray}
Then, by plugging Eq.(\ref{eep3}) back into Eq.(\ref{eep2}), we obtain
\begin{eqnarray}\label{eep4}
\begin{split}
a(Q_h {\bf u}, {\mathbf v})-b({\mathbf v},Q_h^p p)=&({\bf f},{\mathbf v}_0)+\ell_1({\bf u},{\mathbf v})+\tilde{\ell}_1({\bf u},{\mathbf v})-\ell_2(p,{\mathbf v})-\tilde{\ell}_2(p,{\mathbf v})+\tilde{s}(Q_h {\bf u},{\mathbf v})\\
+&\sum_{e \in \mathcal{E}_h^I} \langle \tilde{{\mathbf v}}_b, \bm{\phi}{\bf n}_1^{\mathrm{T}} \rangle_e-\sum_{e \in \mathcal{E}_h^I} \langle \hat{Q}_b(\bm{\psi} {\bf n}_1^{\mathrm{T}}),A_2 \nabla_w {\mathbf v} \rangle_e\\
+&\sum_{T \in \mathcal{T}_{2h}^I}
\langle \hat{Q}_b(\bm{\psi}{\bf n}_1^{\mathrm{T}}),\hat{Q}_b({\mathbf v}_0 {\bf n}_2^{\mathrm{T}})+\tilde{{\mathbf v}}_b \rangle_{\partial T \cap \Gamma}.
\end{split}
\end{eqnarray}
Subtracting Eq.(\ref{eep4}) from Eq.(\ref{algorithm 1}) can lead to
\begin{eqnarray}\label{eep5}
\begin{split}
&a({\bf e}_h, {\mathbf v})-b({\mathbf v},\varepsilon_h)\\
=&\ell_1({\bf u},{\mathbf v})+\tilde{\ell}_1({\bf u},{\mathbf v})-\ell_2(p,{\mathbf v})-\tilde{\ell}_2(p,{\mathbf v})+\tilde{s}(Q_h {\bf u}, {\mathbf v}).
\end{split}
\end{eqnarray}
Another step in the proof is to verify the Eq.(\ref{eep1}) by incorporating with
Eq.(\ref{weak divergence 1}). Therefore, we have the following equations: for $T \in \mathcal{T}_h^S$,
\begin{eqnarray}\label{eep6}
\begin{split}
(\nabla_w \cdot (Q_h {\bf u}), q)_T=&-(Q_0 {\bf u}, \nabla q)_T + \langle Q_b {\bf u} \cdot {\bf n}, q \rangle_{\partial T}\\
=&-({\bf u},\nabla q)_T + \langle {\bf u}, q{\bf n} \rangle_{\partial T}\\
=&(\nabla \cdot {\bf u} , q)_T,
\end{split}
\end{eqnarray}
and for $T \in \mathcal{T}_{1h}^I$,
\begin{eqnarray}\label{eep7}
\begin{split}
&(\nabla_w \cdot (Q_h {\bf u}), q)_T\\
=&-(Q_0 {\bf u}_1, \nabla q)_T+\langle Q_b {\bf u}_1 ,q{\bf n}\rangle_{\partial T \setminus \Gamma}+\langle \hat{Q}_b({\bf u}_1 {\bf n}_1^{\mathrm{T}}),
qI \rangle_{\partial T \cap \Gamma}\\
=&-({\bf u}_1, \nabla q)_T + \langle {\bf u}_1, q{\bf n}_1 \rangle_{\partial T \setminus \Gamma}+\langle {\bf u}_1 {\bf n}_1^{\mathrm{T}},qI \rangle_{\partial T \cap \Gamma}\\
=&-({\bf u}_1, \nabla q)_T+\langle {\bf u}_1, q{\bf n}_1 \rangle_{\partial T}\\
=&(\nabla \cdot {\bf u}_1, q)_T,
\end{split}
\end{eqnarray}
and for $T \in \mathcal{T}_{2h}^I$, using the interface condition (\ref{interface 1}) leads to
\begin{eqnarray}\label{eep8}
\begin{split}
&(\nabla_w \cdot (Q_h {\bf u}), q)_T\\
=&-(Q_0 {\bf u}_2, \nabla q)_T+\langle Q_b {\bf u}_2 ,q{\bf n}_2\rangle_{\partial T \setminus \Gamma}-\langle \hat{Q}_b({\bf u}_1 {\bf n}_1^{\mathrm{T}}),
qI \rangle_{\partial T \cap \Gamma}\\
=&-({\bf u}_2, \nabla q)_T + \langle {\bf u}_2, q{\bf n}_2 \rangle_{\partial T \setminus \Gamma}+\langle {\bf u}_2 {\bf n}_2^{\mathrm{T}},qI \rangle_{\partial T \cap \Gamma}-\langle \hat{Q}_b(\bm{\psi}{\bf n}_2^{\mathrm{T}}),qI \rangle_{\partial T \cap \Gamma}\\
=&-({\bf u}_2, \nabla q)_T+\langle {\bf u}_2, q{\bf n}_2 \rangle_{\partial T}-\langle \hat{Q}_b(\bm{\psi}{\bf n}_2^{\mathrm{T}}),qI \rangle_{\partial T \cap \Gamma}\\
=&(\nabla \cdot {\bf u}_2, q)_T-\langle \hat{Q}_b(\bm{\psi}{\bf n}_2^{\mathrm{T}}),qI \rangle_{\partial T \cap \Gamma},
\end{split}
\end{eqnarray}
Combing Eqs.(\ref{eep6}), (\ref{eep7}) and (\ref{eep8}), we see
\begin{eqnarray}\label{eep9}
\begin{split}
(\nabla_w \cdot (Q_h {\bf u}), q)=&\sum_{T \in \mathcal{T}_h}(\nabla_w \cdot (Q_h {\bf u}), q)_T\\
=&\sum_{T \in \mathcal{T}_h} (\nabla \cdot {\bf u}, q)_T-\sum_{T \in \mathcal{T}_{2h}^I} \langle \hat{Q}_b(\bm{\psi}{\bf n}_1^{\mathrm{T}}), qI\rangle_{\partial T \cap \Gamma},
\end{split}
\end{eqnarray}
Then, subtracting Eq.(\ref{eep9}) from Eq.(\ref{algorithm 2}), we drive that
\begin{eqnarray}
b({\bf e}_h,q)=0
\end{eqnarray}
Hence the error equations are proved.
\end{proof}

\section{Error Estimates in a Energy Norm}
In the section, we will obtain optimal order estimates for velocity error ${\bf e}_h$ in a energy norm and for pressure error $\varepsilon_h$ in a $L^2$ norm.
\begin{lemma}
For any ${\mathbf v}=\{ {\mathbf v}_0, {\mathbf v}_b\} \in V_h$, we obtain
\begin{eqnarray}\label{estimate 1}
\sum_{T \in \mathcal{T}_h} \| \nabla {\mathbf v}_0 \|^2_T \leqslant C {|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|}^2.
\end{eqnarray}
\end{lemma}
\begin{proof}
	For the proof of $T \in \mathcal{T}_{h}^S$, the reader is referred to \cite{WGStokes1}.For $T \in \mathcal{T}_{1h}^I$,
	\begin{align*}
		(\nabla {\mathbf v}_0 ,\nabla {\mathbf v}_0)_T &= -({\mathbf v}_0, \nabla \cdot (\nabla {\mathbf v}_0))_T +
		 \langle {\mathbf v}_0, \nabla {\mathbf v}_0 \cdot {\bf n} \rangle_{\partial T}\\
		 &=-({\mathbf v}_0, \nabla \cdot (\nabla {\mathbf v}_0))_T+
		 \langle {\mathbf v}_b, \nabla {\mathbf v}_0 \cdot {\bf n} \rangle_{\partial T \setminus \Gamma}
		 +\langle \tilde{{\mathbf v}}_b, \nabla {\mathbf v}_0 \rangle_{\partial T \cap \Gamma}\\
		 &+\langle Q_b {\mathbf v}_0 -{\mathbf v}_b, \nabla {\mathbf v}_0 \cdot {\bf n} \rangle_{\partial T \setminus \Gamma}
		 +\langle \hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})- \tilde{{\mathbf v}}_b ,\nabla {\mathbf v}_0 \cdot {\bf n} \rangle_{\partial T \cap \Gamma}\\
		 &=(\nabla_w {\mathbf v}, \nabla {\mathbf v}_0)_T+\langle Q_b {\mathbf v}_0 -{\mathbf v}_b, \nabla {\mathbf v}_0 \cdot {\bf n} \rangle_{\partial T \setminus \Gamma}+\langle \hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})- \tilde{{\mathbf v}}_b ,\nabla {\mathbf v}_0 \rangle_{\partial T \cap \Gamma}\\.
	\end{align*}
Then it follows from the Cauchy-Schwarz inequality and the trace inequality that
\begin{align*}
	\| \nabla {\mathbf v}_0\|_T^2 &\leqslant \|\nabla_w {\mathbf v} \|_T \| \nabla {\mathbf v}_0 \|_T +\| Q_b {\mathbf v}_0 -{\mathbf v}_b\|_{\partial T \setminus \Gamma} \| \nabla {\mathbf v}_0 \|_{\partial T \setminus \Gamma}
	+\| \hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})- \tilde{{\mathbf v}}_b\|_{\partial T \cap \Gamma}
	\|\nabla {\mathbf v}_0 \|_{\partial T \cap \Gamma}\\
	&\leqslant \|\nabla_w {\mathbf v} \|_T \| \nabla {\mathbf v}_0 \|_T+ h_T^{-\frac{1}{2}}\| Q_b {\mathbf v}_0 -{\mathbf v}_b\|_{\partial T \setminus \Gamma} \| \nabla {\mathbf v}_0 \|_T
	+h_T^{-\frac{1}{2}}\| \hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})- \tilde{{\mathbf v}}_b \|_{\partial T \setminus \Gamma} \| \nabla {\mathbf v}_0 \|_T.
\end{align*}
Therefore,
\begin{align*}
	\| \nabla {\mathbf v}_0\|_T \leqslant \|\nabla_w {\mathbf v} \|_T+h_T^{-\frac{1}{2}}\| Q_b {\mathbf v}_0 -{\mathbf v}_b\|_{\partial T \setminus \Gamma}+h_T^{-\frac{1}{2}}\| \hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})- \tilde{{\mathbf v}}_b \|_{\partial T \setminus \Gamma},
\end{align*}
it's easy to see that $\| \nabla {\mathbf v}_0\|_T^2 \leqslant C {|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|}^2$. For $T \in \mathcal{T}_{2h}^I$, we use the same method to prove $\| \nabla {\mathbf v}_0\|_T^2 \leqslant C {|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|}^2$. Finally, the proof is completed.
\end{proof}

\begin{lemma}\label{normal vector}
Suppose $T \in \mathcal{T}_h^I$ is the curved element which doesn't contain the origin. The equation for the interface can be written as $r=r(\theta), a \leqslant \theta \leqslant b$. Then the definition of the normal vector ${\bf n}$ on the interface extends to the whole element T, and ${\bf n}$ is twice continuously differentiable and bounded in the element T.
\end{lemma}
\begin{proof}
According to the polar equation of the interface, the coordinates in Cartesian coordinate system can be expressed as
\begin{eqnarray}\label{nvp1}
x(\theta)=r(\theta) \cos \theta,~ y(\theta)=r(\theta) \sin \theta,
\end{eqnarray}
So the normal vector ${\bf n}$ at the interface can be written as
\begin{eqnarray}\label{nvp2}
{\bf n}=\left(\begin{array}{ccc}
-r^{'}(\theta) \sin \theta-r(\theta) \cos \theta\\
r^{'}(\theta) \cos \theta-r(\theta) \sin \theta
\end{array}\right).
\end{eqnarray}
Since the interface does not pass through the origin, the following inverse function group can be determined from the function group (\ref{nvp1}),
\begin{eqnarray}\label{nvp3}
\left \{\begin{array}{l}
r=\sqrt{x^2+y^2} \\
\theta={\left \{\begin{array}{ll}
\arctan \frac{y}{x} & x > 0  \\
\frac{\pi}{2} & x = 0  \\
\pi+\arctan \frac{y}{x} & x < 0
\end{array}\right.}
\end{array}\right..
\end{eqnarray}
By plugging (\ref{nvp3}) back into the normal vector (\ref{nvp2}), we can write the normal vector as an expression for $x$ and $y$. Since the function $r(x,y)$ and $\theta(x,y)$ can define inside the element, the normal vector ${\bf n}$ can be extended from the boundary to the whole element.

The next step in the proof is to verify ${\bf n}$ is teice continuously differentiabe. Because $ r(x, y), \theta(x, y) $ are smooth, it's obvious that ${\bf n}$ is smooth. When $(x, y) \in T$, we have $\theta(x, y) \in [a,b]$, so $r^{'}(\theta(x,y))$ and ${\bf n}$  is bounded in T. For $\nabla {\bf n}$,
\begin{align*}
\nabla {\bf n} =\left(\begin{array}{ccc}
-r^{''}(\theta)\frac{\partial \theta}{\partial x} \sin \theta -r^{'}(\theta) \frac{\partial \theta}{\partial x} \cos \theta -1 & -r^{''}(\theta)\frac{\partial \theta}{\partial y} \sin \theta -r^{'}(\theta)\frac{\partial \theta}{\partial y} \cos \theta \\
 r^{''}(\theta)\frac{\partial \theta}{\partial x} \cos \theta -r^{'}(\theta) \frac{\partial \theta}{\partial x} \sin \theta  & r^{''}(\theta)\frac{\partial \theta}{\partial y} \cos \theta -r^{'}(\theta)\frac{\partial \theta}{\partial y} \sin \theta-1
\end{array}\right).
\end{align*}
Since $r(\theta)$ and $\theta (x, y)$ is smooth, $\nabla {\bf n}$ is continuous in element T. Moreover, we have
\begin{align*}
|\frac{\partial \theta}{\partial x}|=|\frac{-y}{x^2+y^2}|\leqslant M_1,\\
|\frac{\partial \theta}{\partial y}|=|\frac{x}{x^2+y^2}|\leqslant M_2.
\end{align*}
Therefore, $\nabla {\bf n}$ is bounded in the element T. The proof of the second derivative of n is bounded can be completed by the method analogous to that used above.
\end{proof}

\begin{lemma}\label{projector inequality}
Assume that e is the curved edge of the element T. Then there exists a positive constant C such that
\begin{eqnarray}
\| {\mathbf v}_0 {\bf n}_1^{\mathrm{T}}-\hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}}) \|_e^2 \leqslant C h_T^{2k-1} \| \nabla {\mathbf v}_{10} \|^2_{k,T},~\forall \, T \in \mathcal{T}_{1h}^I,\label{projector inequality 1}\\
\| {\mathbf v}_0 {\bf n}_2^{\mathrm{T}}-\hat{Q}_b({\mathbf v}_0 {\bf n}_2^{\mathrm{T}}) \|_e^2 \leqslant C h_T^{2k-1} \| \nabla {\mathbf v}_{20} \|_{k,T}^2,~\forall \, T \in \mathcal{T}_{2h}^I,\label{projector inequality 2}
\end{eqnarray}
where k=1,2.
\end{lemma}

\begin{proof}
First, let $\bar{Q}_b$ and $\bar{Q}_0$ be the $L^2$ projector onto $P_k(e)$ and $P_k(T)$. we give some notations used in the proof:
\begin{align*}
{\mathbf v}_0 =\left(\begin{array}{ccc}
v_{10}(x,y)\\
v_{20}(x,y)
\end{array}\right),
~{\bf n}_1=\left(\begin{array}{ccc}
n_{11}\\
n_{12}
\end{array}\right),
~\hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})=\left(\begin{array}{ccc}
\bar{Q}_b(v_{10}n_{11}) & \bar{Q}_b(v_{10}n_{12})\\
\bar{Q}_b(v_{20}n_{11}) & \bar{Q}_b(v_{20}n_{12})
\end{array}\right),
\end{align*}
According to the above notations, we have
\begin{align*}
\| {\mathbf v}_0 {\bf n}_1^{\mathrm{T}}-\hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}}) \|_e^2&=\| v_{10}n_{11}-\bar{Q}_b(v_{10}n_{11})\|_e^2+\| v_{10}n_{12}-\bar{Q}_b(v_{10}n_{12})\|_e^2\\
&+\| v_{20}n_{11}-\bar{Q}_b(v_{20}n_{11})\|_e^2+
\| v_{20}n_{12}-\bar{Q}_b(v_{20}n_{12})\|_e^2.
\end{align*}
By using the Lemma \ref{normal vector}, we have
\begin{align*}
&\| v_{10}n_{11}-\bar{Q}_b(v_{10}n_{11})\|_e^2
\leqslant \| v_{10}n_{11}-\bar{Q}_b(v_{10}n_{11})\|_{\partial T}^2\\
& \leqslant \| v_{10}n_{11}-\bar{Q}_b(v_{10}n_{11}) \|_{\partial T}^2\\
& \leqslant C_1\Big(h_T^{-1} \| v_{10}n_{11}-\bar{Q}_b(v_{10}n_{11})\|_T^2
+h_T \|\nabla (v_{10}n_{11}-\bar{Q}_b(v_{10}n_{11})) \|_{T}^2 \Big)\\
& \leqslant C_1 h_T^{2k-1} | v_{10}n_{11}|_{k,T}^2\\
\end{align*}
For k=1, we get
\begin{align*}
	&\| v_{10}n_{11}-\bar{Q}_b(v_{10}n_{11})\|_e^2 \\
	&\leqslant 2 C_1 h_T \Big( \| n_{11} \nabla v_{10}\|_T^2 +\| v_{10} \nabla n_{11} \|_T^2 \Big)\\
	& \leqslant 2 C_1 h_T \Big( \| n_{11} \|_T^2 \| \nabla v_{10} \|_T^2
	+ \| \nabla n_{11} \|_T^2 \| v_{10} \|_T^2 \Big)\\
	& \leqslant C h_T \| \nabla v_{10}\|_T^2
\end{align*}
Similarly, for k=2, we obtain
\begin{align*}
	&\| v_{10}n_{11}-\bar{Q}_b(v_{10}n_{11})\|_e^2 \\
	&\leqslant C h_T^3 | v_{10}n_{11}|_{2,T}^2\\
	&\leqslant C h_T^3 (\| v_{10}\|_T^2 \| \vartriangle n_{11}\|_T^2+\| \vartriangle v_{10}\|_T^2 \|n_{11}\|_T^2+\| \nabla v_{10}\|_T^2\| \nabla n_{11}\|_T^2)\\
	&\leqslant C_1(\| v_{10}\|_T^2+\| \vartriangle v_{10}\|_T^2+\| \nabla v_{10}\|_T^2)\\
	&\leqslant C|v_{10}|^2_{2,T}
\end{align*}
Similarly, the proof of (\ref{projector inequality 2}) can be completed.
\end{proof}

\begin{lemma}\label{H1 estimates}
Suppose that the partiton $\mathcal{T}_h$ of the domain $\Omega$ is shape regular. Then for ${\bf u}_i \in[H^{k+1}(\Omega_i)]^2 , p_i \in H^r(\Omega_i)$ with $i=1,2$, we have
\begin{eqnarray}
|\ell_1({\bf u},{\mathbf v})| \leqslant C h^k(\|{\bf u}_1\|_{k+1, \Omega_1}+\| {\bf u}_2 \|_{k+1, \Omega_2}) {|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|},\label{H1 estimates 1}\\
|\tilde{\ell}_1({\bf u},{\mathbf v})| \leqslant C h^k(\|{\bf u}\|_{k+1, \Omega_1}+\|{\bf u}\|_{k+1, \Omega_2}) {|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|},\label{H1 estimates 2}\\
|\ell_2(p,{\mathbf v})| \leqslant C h^r(\| p_1 \|_{r, \Omega_1}+\| p_2 \|_{r, \Omega_2}) {|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|},\label{H1 estimates 3}\\
|\tilde{\ell}_2(p,{\mathbf v})| \leqslant C h^r(\| p_1 \|_{r, \Omega_1}+\| p_2 \|_{r, \Omega_2}) {|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|},\label{H1 estimates 4}\\
|\tilde{s}(Q_h {\bf u}, {\mathbf v})| \leqslant C h^k(\|{\bf u}_1\|_{k+1, \Omega_1}+\| {\bf u}_2 \|_{k+1, \Omega_2}) {|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|},\label{H1 estimates 5}
\end{eqnarray}
where ${\mathbf v}=\{{\mathbf v}_0, {\mathbf v}_b\} \in  V_h$.
\end{lemma}

\begin{proof}
As to (\ref{H1 estimates 1}), according to Cauchy-Schwarz inequality and the Cauchy inequality, we have
\begin{eqnarray}\label{HEP1}
\begin{split}
|\ell_1({\bf u},{\mathbf v})|=&\sum_{T \in \mathcal{T}_h^S} \langle {\mathbf v}_0 -{\mathbf v}_b, A (\nabla {\bf u}){\bf n} -A \mathbb{Q}_h(\nabla {\bf u}){\bf n} \rangle_{\partial T}\\
&\leqslant A \sum_{T \in \mathcal{T}_h^S} \| {\mathbf v}_0 -{\mathbf v}_b  \|_{\partial T} \| (\nabla {\bf u}){\bf n} -\mathbb{Q}_h(\nabla {\bf u}){\bf n} \|_{\partial T}\\
&\leqslant A \Big( \sum_{T \in \mathcal{T}_h^S}  \| {\mathbf v}_0 -{\mathbf v}_b  \|^2_{\partial T} \Big)^{\frac{1}{2}} \Big( \sum_{T \in \mathcal{T}_h^S} \| \nabla {\bf u} - \mathbb{Q}_h(\nabla {\bf u}) \|^2_{\partial T}  \Big)^{\frac{1}{2}},
\end{split}
\end{eqnarray}
It follows from the trace inequality and the estimate (\ref{projectorestimate2}) that
\begin{eqnarray}\label{HEP2}
\begin{split}
&  \sum_{T \in \mathcal{T}_h^S} \| \nabla {\bf u} - \mathbb{Q}_h(\nabla u) \|^2_{\partial T}\\
& \leqslant \sum_{T \in \mathcal{T}_h^S} C \Big( h_T^{-1} \| \nabla {\bf u} - \mathbb{Q}_h(\nabla u) \|^2_T
+ h_T \| \nabla(\nabla {\bf u} - \mathbb{Q}_h(\nabla u))\|^2_T \Big)\\
& \leqslant C h^{2k-1}\| {\bf u} \|^2_{k+1}.
\end{split}
\end{eqnarray}
Using the triangle inequality, the trace inequality, Cauchy-Schwarz inequality and approximation properties of $L^2$ projection operators, as well as the estimate (\ref{estimate 1}), we obtain
\begin{eqnarray}\label{HEP3}
\begin{split}
& \sum_{T \in \mathcal{T}_h^S} \| {\mathbf v}_0 -{\mathbf v}_b  \|^2_{\partial T}\\
&\leqslant \sum_{T \in \mathcal{T}_h^S} 2 \Big( \|{\mathbf v}_0 - Q_b{\mathbf v}_0 \|^2_{\partial T}+ \|Q_b{\mathbf v}_0-{\mathbf v}_b \|^2_{\partial T} \Big)\\
&\leqslant C_1\Big(\sum_{T \in \mathcal{T}_h^S} h_T \| \nabla {\mathbf v}_0 \|_T^2 \Big)
+C_2 h \Big(\sum_{T \in \mathcal{T}_h^S} h_T^{-1}\|Q_b{\mathbf v}_0-{\mathbf v}_b \|^2_{\partial T} \Big)\\
&\leqslant C h {|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|}^2.
\end{split}
\end{eqnarray}
Using the Eqs.(\ref{HEP2}) and (\ref{HEP3}), we complete the proof of (\ref{H1 estimates 1}). \\

Next, as to (\ref{H1 estimates 2}), we only need to verify the integration on the interface, i.e.
\begin{eqnarray}\label{HEP4}
\begin{split}
&\sum_{T \in \mathcal{T}_{1h}^I} \langle {\mathbf v}_0{\bf n}_1^{\mathrm{T}}-\tilde{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}}), A_1 \nabla {\bf u}_1 - A_1 \mathbb{Q}_h (\nabla {\bf u}_1) \rangle_{\partial T \cap \Gamma}\\
&\leqslant \sum_{T \in \mathcal{T}_{1h}^I} \| {\mathbf v}_0{\bf n}_1^{\mathrm{T}}-\tilde{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}}) \|_{\partial T \cap \Gamma} \| \nabla {\bf u}_1 -  \mathbb{Q}_h (\nabla {\bf u}_1) \|_{\partial T \cap \Gamma}\\
&\leqslant \Big(\sum_{T \in \mathcal{T}_{1h}^I} \|{\mathbf v}_0{\bf n}_1^{\mathrm{T}}-\tilde{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})\|^2_{\partial T \cap \Gamma} \Big)^{\frac{1}{2}}
\Big(\sum_{T \in \mathcal{T}_{1h}^I} \|\nabla {\bf u}_1 -  \mathbb{Q}_h (\nabla {\bf u}_1)\|^2_{\partial T \cap \Gamma} \Big)^{\frac{1}{2}}\\
&\leqslant C h^{k} \|\nabla {\mathbf v}_0\| \|{\bf u}_1\|_{k+1,\Omega_1} \\
&\leqslant C h^{k} {|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|} \|{\bf u}_1\|_{k+1,\Omega_1}.
\end{split}
\end{eqnarray}
Similarly, we obtain
\begin{eqnarray}\label{HEP5}
\begin{split}
&\sum_{T \in \mathcal{T}_{1h}^I} \langle \tilde{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})-\tilde{{\mathbf v}}_b, A_1 \nabla {\bf u}_1 - A_1 \mathbb{Q}_h (\nabla {\bf u}_1) \rangle_{\partial T \cap \Gamma}\\
&\leqslant \Big(\sum_{T \in \mathcal{T}_{1h}^I}  h_T^{-1} \| \tilde{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})-\tilde{{\mathbf v}}_b \|^2_{\partial T \cap \Gamma} \Big)^{\frac{1}{2}}\\
&\Big(\sum_{T \in \mathcal{T}_{1h}^I} \| A_1 \nabla {\bf u}_1 - A_1 \mathbb{Q}_h (\nabla {\bf u}_1) \|^2_{\partial T \cap \Gamma}\Big)^{\frac{1}{2}}\\
&\leqslant C h^k {|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|} \|{\bf u}_1\|_{k+1,\Omega_1}.
\end{split}
\end{eqnarray}
Combining (\ref{HEP4}) and (\ref{HEP5}), we have
\begin{eqnarray}\label{HEP6}
\begin{split}
&\sum_{T \in \mathcal{T}_{1h}^I} \langle {\mathbf v}_0{\bf n}_1^{\mathrm{T}}-\tilde{{\mathbf v}}_b, A_1 \nabla {\bf u}_1 - A_1 \mathbb{Q}_h (\nabla {\bf u}_1) \rangle_{\partial T \cap \Gamma}\\
&= \sum_{T \in \mathcal{T}_{1h}^I} \langle {\mathbf v}_0{\bf n}_1^{\mathrm{T}}-\hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}}), A_1 \nabla {\bf u}_1 - A_1 \mathbb{Q}_h (\nabla {\bf u}_1)
\rangle_{\partial T \cap \Gamma}\\
&+\langle \hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}}) -\tilde{{\mathbf v}}_b,A_1 \nabla {\bf u}_1 - A_1 \mathbb{Q}_h (\nabla {\bf u}_1)  \rangle_{\partial T \cap \Gamma}\\
& \leqslant C h^k {|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|} \|{\bf u}_1\|_{k+1,\Omega_1}.
\end{split}
\end{eqnarray}
Using the same method for $T \in \mathcal{T}_{2h}^I$ leads to
\begin{eqnarray}\label{HEP7}
\begin{split}
&\sum_{T \in \mathcal{T}_{2h}^I} \langle {\mathbf v}_0{\bf n}_2^{\mathrm{T}}+\tilde{{\mathbf v}}_b, A_2 \nabla {\bf u}_2 - A_2 \mathbb{Q}_h (\nabla {\bf u}_2) \rangle_{\partial T \cap \Gamma}\\
&= \sum_{T \in \mathcal{T}_{2h}^I} \langle {\mathbf v}_0{\bf n}_2^{\mathrm{T}}-\hat{Q}_b({\mathbf v}_0 {\bf n}_2^{\mathrm{T}}), A_2 \nabla {\bf u}_2 - A_2 \mathbb{Q}_h (\nabla {\bf u}_2)
\rangle_{\partial T \cap \Gamma}\\
&+\langle \hat{Q}_b({\mathbf v}_0 {\bf n}_2^{\mathrm{T}}) + \tilde{{\mathbf v}}_b,A_2 \nabla {\bf u}_2 - A_2 \mathbb{Q}_h (\nabla {\bf u}_2)  \rangle_{\partial T \cap \Gamma}\\
& \leqslant C h^k {|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|} \|{\bf u}_2\|_{k+1,\Omega_2}.
\end{split}
\end{eqnarray}
Therefore, the proof of Eq.(\ref{H1 estimates 2}) is completed.

For the (\ref{H1 estimates 3}), the same techniques for proving (\ref{H1 estimates 1}) can be applied to obtain the following estimate:
\begin{eqnarray}\label{HEP8}
\begin{split}
|\ell_2(p,{\mathbf v})|=&\sum_{T \in \mathcal{T}_h^S} \langle {\mathbf v}_0 -{\mathbf v}_b,p{\bf n}-(Q_h^p p ){\bf n} \rangle_{\partial T}\\
&\leqslant C h^r (\|p_1\|_{r,\Omega_1}+\|p_2\|_{r,\Omega_2}){|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|}
\end{split}
\end{eqnarray}
Meanwhile, the proof of the (\ref{H1 estimates 4}) is same as the Eq.(\ref{H1 estimates 2}) and so is omitted. For $\tilde{S}(Q_h {\bf u},{\mathbf v})$, we divide our proof into four steps.
Firstly, we consider $T \in \mathcal{T}_h^S$,
\begin{eqnarray}\label{HEP9}
\begin{split}
&\sum_{T \in \mathcal{T}_h^S} h_T^{-1} \langle Q_b(Q_0 {\bf u})-Q_b {\bf u}, Q_b{\mathbf v}_0-{\mathbf v}_b \rangle_{\partial T}\\
=&\sum_{T \in \mathcal{T}_h^S}h_T^{-1} \langle Q_0 {\bf u}-{\bf u}, Q_b{\mathbf v}_0-{\mathbf v}_b \rangle_{\partial T}\\
& \leqslant \sum_{T \in \mathcal{T}_h^S} h_T^{-1} \| Q_0 {\bf u}-{\bf u} \|_{\partial T} \| Q_b{\mathbf v}_0-{\mathbf v}_b \|_{\partial T}\\
& \leqslant \Big( \sum_{T \in \mathcal{T}_h^S} h_T^{-1} \| Q_0 {\bf u}-{\bf u} \|_{\partial T}^2  \Big)^{\frac{1}{2}} \Big( \sum_{T \in \mathcal{T}_h^S} h_T^{-1} \| Q_b{\mathbf v}_0-{\mathbf v}_b \|_{\partial T}^2 \Big)^{\frac{1}{2}}.
\end{split}
\end{eqnarray}
For the estimate of the first term, we drive by using the trace inequality and the estimate (\ref{projectorestimate1}) that
\begin{eqnarray}\label{HEP10}
\begin{split}
&\sum_{T \in \mathcal{T}_h^S} h_T^{-1} \| Q_0 {\bf u}-{\bf u} \|_{\partial T}^2\\
&\leqslant \sum_{T \in \mathcal{T}_h^S} C h_T^{-1} (h_T^{-1}\|Q_0 {\bf u} -{\bf u}\|_T^2
+h_T \|\nabla (Q_0 {\bf u} -{\bf u})\|_T^2)\\
& \leqslant C h^{2k} (\|{\bf u}_1\|_{k+1,\Omega_1}^2+\|{\bf u}_2\|_{k+1,\Omega_2}^2).
\end{split}
\end{eqnarray}
Next, for $T \in \mathcal{T}_{1h}^I$, we see that
\begin{eqnarray}\label{HEP11}
\begin{split}
\tilde{s}(Q_h {\bf u},{\mathbf v})=&\sum_{T \in \mathcal{T}_{1h}^I } h_T^{-1}\langle Q_b(Q_0 {\bf u})-Q_b {\bf u}, Q_b{\mathbf v}_0-{\mathbf v}_b
\rangle_{\partial T \backslash \Gamma}\\
+&\sum_{T \in \mathcal{T}_{1h}^I } h_T^{-1} \langle \hat{Q}_b((Q_0{\bf u}_1){\bf n}_1^{\mathrm{T}})-\hat{Q}_b({\bf u}_1 {\bf n}_1^{\mathrm{T}}), \hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})-\tilde{{\mathbf v}}_b \rangle_{\partial T \cap \Gamma},
\end{split}
\end{eqnarray}
the first term is same as the triangle element, so we see on the curved edge
\begin{eqnarray}\label{HEP12}
\begin{split}
&\sum_{T \in \mathcal{T}_{1h}^I } h_T^{-1} \langle \hat{Q}_b((Q_0{\bf u}_1){\bf n}_1^{\mathrm{T}})-\hat{Q}_b({\bf u}_1 {\bf n}_1^{\mathrm{T}}), \hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})-\tilde{{\mathbf v}}_b \rangle_{\partial T \cap \Gamma}\\
=&\sum_{T \in \mathcal{T}_{1h}^I } h_T^{-1} \langle (Q_0{\bf u}_1){\bf n}_1^{\mathrm{T}}- {\bf u}_1 {\bf n}_1^{\mathrm{T}}, \hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})-\tilde{{\mathbf v}}_b  \rangle_{\partial T \cap \Gamma}\\
\leqslant &\sum_{T \in \mathcal{T}_{1h}^I }  h_T^{-1} \| Q_0{\bf u}_1- {\bf u}_1\|_{\partial T \cap \Gamma} \| \hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})-\tilde{{\mathbf v}}_b \|_{\partial T \cap \Gamma}\\
\leqslant &\Big( \sum_{T \in \mathcal{T}_{1h}^I} h_T^{-1} \| Q_0{\bf u}_1- {\bf u}_1 \|^2_{\partial T \cap \Gamma} \Big)^{\frac{1}{2}}
\Big( \sum_{T \in \mathcal{T}_{1h}^I} h_T^{-1} \|\hat{Q}_b({\mathbf v}_0 {\bf n}_1^{\mathrm{T}})-\tilde{{\mathbf v}}_b \|^2_{\partial T \cap \Gamma} \Big)^{\frac{1}{2}}\\
\leqslant &  C h^k \| {\bf u}_1 \|_{k+1,\Omega_1} {|\hspace{-.02in}|\hspace{-.02in}|} {\mathbf v} {|\hspace{-.02in}|\hspace{-.02in}|}
\end{split}
\end{eqnarray}
Finally, for the $T \in \mathcal{T}_{2h}^I$, the proof is similar to the above equation. By incorporating from Eq.(\ref{HEP9}) to (\ref{HEP12}), Eq.(\ref{H1 estimates 5}) is completed.
\end{proof}

\begin{theorem}\label{H1ERROR}
Assume the exact solution  $({\bf u}_i ,p_i) \in [H^{k+1}(\Omega_i)]^2 \times H^k (\Omega_i)$ with $i=1,2$  of the model (\ref{model 1})-(\ref{interface 2}) and let $({\bf u}_h,p_h) \in V_h \times W_h$ be the weak Galerkin finite element solutions, we have the following inequality holds true,
\begin{eqnarray}\label{H1error}
\begin{split}
&{|\hspace{-.02in}|\hspace{-.02in}|} Q_h {\bf u} - {\bf u}_h {|\hspace{-.02in}|\hspace{-.02in}|} +\| Q_h^p p -p_h\| \\
\leqslant & C h^k (\|{\bf u}_1\|_{k+1,\Omega_1}+\|{\bf u}_2\|_{k+1,\Omega_2}+\|p_1\|_{k,\Omega_1} +\|p_2\|_{k,\Omega_2})
\end{split}
\end{eqnarray}
\end{theorem}
\begin{proof}
Let ${\mathbf v}= {\bf e}_h$ in the Eq.(\ref{error equation 1}) and $q = \varepsilon_h$ in the Eq.(\ref{error equation 2}) and add the two equations. Then we have
\begin{eqnarray}\label{H1errorP1}
a({\bf e}_h, {\bf e}_h)=\ell_1({\bf u},{\bf e}_h)+\tilde{\ell}_1 ({\bf u},{\bf e}_h)-\ell_2(p,{\bf e}_h)-\tilde{\ell}_2 (p,{\bf e}_h)+\tilde{s}(Q_h {\bf u},{\bf e}_h).
\end{eqnarray}
According to Eqs.(\ref{H1 estimates 1})-(\ref{H1 estimates 5}), we drive the following estimate
\begin{eqnarray}\label{H1errorP2}
{|\hspace{-.02in}|\hspace{-.02in}|} {\bf e}_h {|\hspace{-.02in}|\hspace{-.02in}|}^2 \leqslant C h^k (\|{\bf u}_1\|_{k+1,\Omega_1}+\|{\bf u}_2\|_{k+1,\Omega_2}+\|p_1\|_{k,\Omega_1} +\|p_2\|_{k,\Omega_2}).
\end{eqnarray}
The proof of the first term ${|\hspace{-.02in}|\hspace{-.02in}|} {\bf e}_h {|\hspace{-.02in}|\hspace{-.02in}|}$ in Eq.(\ref{H1error}) is completed.

To estimate the second term  $ \| \varepsilon_h \|$, we drive from Eq.(\ref{error equation 1}) that
\begin{eqnarray}\label{H1errorP3}
b({\mathbf v},\varepsilon_h)=a({\bf e}_h,{\mathbf v})-\ell_1({\bf u},{\mathbf v})-\tilde{\ell}_1({\bf u},{\mathbf v})+\ell_2(p,{\mathbf v})+\tilde{\ell}_2(p,{\mathbf v})-\tilde{s}(Q_h {\bf u} ,{\mathbf v})
\end{eqnarray}
Using the equation above, (\ref{H1 estimates 1})-(\ref{H1 estimates 5}) and (\ref{H1errorP2}), we have
\begin{eqnarray}\label{H1errorP4}
|b({\mathbf v},\varepsilon_h)| \leqslant C h^k (\|{\bf u}_1\|_{k+1,\Omega_1}+\|{\bf u}_2\|_{k+1,\Omega_2}+\|p_1\|_{k,\Omega_1} +\|p_2\|_{k,\Omega_2}).
\end{eqnarray}
Incorporated with the inf-sup condition (\ref{InfSupCondition}), the estimate (\ref{H1errorP4}) gives
\begin{eqnarray}\label{H1errorP5}
\|\varepsilon_h\| \leqslant C h^k (\|{\bf u}_1\|_{k+1,\Omega_1}+\|{\bf u}_2\|_{k+1,\Omega_2}+\|p_1\|_{k,\Omega_1} +\|p_2\|_{k,\Omega_2}).
\end{eqnarray}
Adding (\ref{H1errorP2}) to (\ref{H1errorP5}) leads to the estimate (\ref{H1error}).
\end{proof}

\section{ Error Estimate in $L^2$ Norm} In general, we use the duality argument to arrive at an $L^2$-error estimate for ${\bf e}_0=Q_0 {\bf u} - {\bf u}_0$. Now we consider the following problem: seeking $({\bf w}; \rho)$ to satisfy
\begin{eqnarray}
- \nabla \cdot (A \nabla {\bf w})+\nabla \rho = {\bf e}_0 \quad in ~\Omega, \label{dual problem 1} \\
\nabla \cdot {\bf w} =0 \quad  in ~ \Omega \label{dual problem 2} \\
{\bf w} = 0 \quad on ~\partial \Omega \label{dual problem 3}
\end{eqnarray}
Assume the solution $({\bf w}; \rho) \in [H^2(\Omega)]^2 \times [H^1(\Omega)]$ of dual problem (\ref{dual problem 1})-(\ref{dual problem 3}) has $[H^2(\Omega)]^2 \times [H^1(\Omega)]$-regularity property, i.e.
\begin{eqnarray}
\| {\bf w} \|_2 +\| \rho \|_1 \leqslant C \|{\bf e}_0 \|.
\end{eqnarray}

\begin{theorem}
Under the assumption in theorem \ref{H1ERROR}, then we have the following optimal order error estimate holds true
\begin{eqnarray}
\| Q_0 {\bf u} -{\bf u}_0\| \leqslant C h^{k+1}(\| {\bf u}_1 \|_{k+1,\Omega_1}+\| {\bf u}_2\|_{k+1,\Omega_2}+\|p_1\|_{k,\Omega_1}+\|p_2\|_{k,\Omega_2}).
\end{eqnarray}
\end{theorem}
\begin{proof}
Due to $({\bf w};\rho)$ satisfies the Eq.(\ref{dual problem 1}) with ${\bf f} = {\bf e}_0=Q_0 {\bf u} -{\bf u}_0$, then letting ${\mathbf v} = {\bf e}_h $ in Eq.(\ref{error equation 1}) leads to
\begin{align*}
\begin{split}
&(\nabla_w(Q_h {\bf w}), A\nabla_w {\bf e}_h)-(\nabla_w \cdot {\bf e}_h, Q_h^p p)\\
=&({\bf e}_0,{\bf e}_0)+\ell_1({\bf w},{\bf e}_h)+\tilde{\ell}_1({\bf w},{\bf e}_h)-\ell_2(\rho,{\bf e}_h)-\tilde{\ell}_2(\rho,{\bf e}_h)\\
\end{split}
\end{align*}
Therefore, from the above equation we obtain
\begin{align*}
\begin{split}
\| {\bf e}_0\|^2=&(\nabla_w(Q_h {\bf w}), A\nabla_w {\bf e}_h)-(\nabla_w \cdot {\bf e}_h, Q_h^p p)\\
-&\ell_1({\bf w},{\bf e}_h)-\tilde{\ell}_1({\bf w},{\bf e}_h)+\ell_2(\rho,{\bf e}_h)+\tilde{\ell}_2(\rho,{\bf e}_h)\\
=&a(Q_h {\bf w},{\bf e}_h)-b({\bf e}_h,Q_h^p \rho)-\ell_1({\bf w},{\bf e}_h)-\tilde{\ell}_1({\bf w},{\bf e}_h)\\
+&\ell_2(\rho,{\bf e}_h)+\tilde{\ell}_2(\rho,{\bf e}_h)-\tilde{s}(Q_h {\bf w},{\bf e}_h),
\end{split}
\end{align*}
where we have used the fact that $ S(Q_h {\bf w},{\bf e}_h)=\tilde{S}(Q_h {\bf w}, {\bf e}_h)$.
Next, combining (\ref{error equation 2}) and the definition of $Q_h^p$, we have
$$
b({\bf e}_h,Q_h^p \rho)=0 \quad b(Q_h {\bf w}, \varepsilon_h)=0
$$
The above equation can change into
\begin{align*}
\begin{split}
\| {\bf e}_0\|^2=&a(Q_h {\bf w},{\bf e}_h)-b(Q_h {\bf w}, \varepsilon_h)-\ell_1({\bf w},{\bf e}_h)-\tilde{\ell}_1({\bf w},{\bf e}_h)\\
+&\ell_2(\rho,{\bf e}_h)+\tilde{\ell}_2(\rho,{\bf e}_h)-\tilde{s}(Q_h {\bf w},{\bf e}_h).
\end{split}
\end{align*}
By substituting (\ref{error equation 1}) into the above equation, we obtain
\begin{align*}
\| {\bf e}_0\|^2=&\ell_1({\bf u},Q_h {\bf w})+\tilde{\ell}_1({\bf u},Q_h {\bf w})-\ell_2(p,Q_h {\bf w})-\tilde{\ell}_2(p,Q_h {\bf w})+\tilde{s}(Q_h {\bf u},Q_h {\bf w})\\
-&\ell_1({\bf w},{\bf e}_h)-\tilde{\ell}_1({\bf w},{\bf e}_h)+\ell_2(\rho,{\bf e}_h)+\tilde{\ell}_2(\rho,{\bf e}_h)-\tilde{s}(Q_h {\bf w},{\bf e}_h).
\end{align*}
According to Lemma \ref{H1 estimates}, we have
\begin{align*}
&|\ell_1({\bf u},Q_h {\bf w})+\tilde{\ell}_1({\bf u},Q_h {\bf w})-\ell_2(p,Q_h {\bf w})-\tilde{\ell}_2(p,Q_h {\bf w})+\tilde{S}(Q_h {\bf u},Q_h {\bf w})|\\
\leqslant & Ch(\| {\bf w}\|_2+\| \rho \|_1){|\hspace{-.02in}|\hspace{-.02in}|} {\bf e}_h {|\hspace{-.02in}|\hspace{-.02in}|} \\
\leqslant &Ch{|\hspace{-.02in}|\hspace{-.02in}|} {\bf e}_h {|\hspace{-.02in}|\hspace{-.02in}|} \|{\bf e}_0\|
\end{align*}
Each of the remaining terms is handled as follows.\\

(1) For $\ell_1({\bf u},Q_h {\bf w})$, we first use the Cauchy-Schwarz inequality, trace inequality and the approximation property of the $L^2$ projection operator to derive
\begin{align*}
&\sum_{T \in \mathcal{T}_h} \langle Q_0 {\bf w} -{\bf w}, A (\nabla {\bf u}){\bf n} -A \mathbb{Q}_h(\nabla {\bf u}) {\bf n} \rangle_{\partial T}\\
\leqslant  & A \sum_{T \in \mathcal{T}_h} \|Q_0 {\bf w} -{\bf w}\|_{\partial T} \| \nabla {\bf u} - \mathbb{Q}_h(\nabla {\bf u}) \|_{\partial T}\\
\leqslant  & A\Big(\sum_{T \in \mathcal{T}_h} \| Q_0 {\bf w} -{\bf w}\|^2_{\partial T} \Big)^{\frac{1}{2}}
\Big(\sum_{T \in \mathcal{T}_h} \| \nabla {\bf u} - \mathbb{Q}_h(\nabla {\bf u}) \|^2_{\partial T} \Big)^{\frac{1}{2}}\\
\leqslant & \Big( C \sum_{T \in \mathcal{T}_h} (h_T^{-1} \| Q_0 {\bf w} -{\bf w} \|_T^2
+ h_T \| \nabla(Q_0 {\bf w} -{\bf w})\|_T^2 )\Big)^{\frac{1}{2}}\\
&\Big( C \sum_{T \in \mathcal{T}_h} (h_T^{-1} \| \nabla {\bf u} - \mathbb{Q}_h(\nabla {\bf u}) \|_T^2
+h_T \| \nabla(\nabla {\bf u} - \mathbb{Q}_h(\nabla {\bf u}))\|_T^2 )\Big)^{\frac{1}{2}}\\
\leqslant & \Big( \sum_{T \in \mathcal{T}_h} C h_T^3 \|{\bf w}\|_2^2 \Big)^{\frac{1}{2}}
\Big( \sum_{T \in \mathcal{T}_h} C h_T^{2k-1} \| {\bf u} \|_{k+1}^2 \Big)^{\frac{1}{2}}\\
\leqslant & Ch^{k+1} \|{\bf w}\|_2 \| {\bf u} \|_{k+1}
\end{align*}
Next, we use same techniques and the fact that $\sum_{e \in \mathcal{E}_h^S} \langle {\bf w}-Q_b {\bf w},A (\nabla {\bf u} ){\bf n}-A \mathbb{Q}_h(\nabla {\bf u}){\bf n} \rangle_e=0$  to drive
\begin{align*}
&\sum_{T \in \mathcal{T}_h} \langle {\bf w}-Q_b {\bf w}, A(\nabla {\bf u} ) {\bf n} -A \mathbb{Q}_h(\nabla {\bf u}){\bf n} \rangle_{\partial T}\\
=&\sum_{e \in \mathcal{E}_h^I} \langle {\bf w} -Q_b {\bf w}, A (\nabla {\bf u} ) {\bf n} -A \mathbb{Q}_h(\nabla {\bf u}) {\bf n} \rangle_e\\
\leqslant & \sum_{e \in \mathcal{E}_h^I }\| {\bf w} -Q_b {\bf w}\|_e \| A \nabla {\bf u}  -A \mathbb{Q}_h(\nabla {\bf u}) \|_e\\
\leqslant & A \Big( \sum_{e \in \mathcal{E}_h^I } \| {\bf w} -Q_b {\bf w}\|_e^2 \Big)^{\frac{1}{2}}
\Big( \sum_{e \in \mathcal{E}_h^I} \| \nabla {\bf u}  -\mathbb{Q}_h(\nabla {\bf u}) \|_e^2  \Big)^{\frac{1}{2}}\\
\leqslant & \Big( \sum_{e \in \mathcal{E}_h^I} \| {\bf w} -Q_b {\bf w}\|_e^2 \Big)^{\frac{1}{2}} \Big( C h^{2k-1} \| {\bf u} \|_{k+1}^2 \Big)^{\frac{1}{2}},
\end{align*}
where
\begin{align*}
\| {\bf w} -Q_b {\bf w}\|_e^2 \leqslant \| {\bf w} - Q_0 {\bf w}\|_e^2 \leqslant \|{\bf w} - Q_0 {\bf w} \|_{\partial T}^2 \leqslant C h_T^3 \|{\bf w}\|_{2}^2.
\end{align*}
Using the above inequality, we have
\begin{align*}
&\sum_{T \in \mathcal{T}_h} \langle {\bf w}-Q_b {\bf w}, A (\nabla {\bf u}) {\bf n} -A \tilde{Q}_h(\nabla {\bf u}) {\bf n} \rangle_{\partial T}\\
\leqslant& Ch^{k+1}\|{\bf w}\|_2 \| {\bf u} \|_{k+1}.
\end{align*}
Therefore,
\begin{align*}
\ell_1({\bf u},Q_h {\bf w})=&\sum_{T \in \mathcal{T}_h} \langle Q_0 {\bf w} -Q_b {\bf w}, A (\nabla {\bf u}) {\bf n} -A Q_h(\nabla {\bf u}) {\bf n} \rangle_{\partial T}\\
=&\sum_{T \in \mathcal{T}_h} \langle Q_0 {\bf w} -{\bf w}, A (\nabla {\bf u}) {\bf n} -A Q_h(\nabla {\bf u}) {\bf n} \rangle_{\partial T}\\
+&\sum_{T \in \mathcal{T}_h} \langle {\bf w} -Q_b {\bf w} , A (\nabla {\bf u}) {\bf n} -A Q_h(\nabla {\bf u}) {\bf n} \rangle_{\partial T}\\
\leqslant & Ch^{k+1} \|{\bf w}\|_2 \| {\bf u} \|_{k+1}
\end{align*}
(2)For $\tilde{\ell}_1 ({\bf u},Q_h {\bf w})$,
\begin{align*}
	\tilde{\ell}_1({\bf u},Q_h {\bf w})=&\sum_{T \in \mathcal{T}_{1h}^I}\Big(
	\langle Q_0 {\bf w} -Q_b {\bf w}, A_1 (\nabla {\bf u}_1){\bf n}_1 - A_1 \mathbb{Q}_h(\nabla {\bf u}_1){\bf n}_1 \rangle_{\partial T \backslash \Gamma} \\
	+ &\langle (Q_0{\bf w}) {\bf n}_1^{\mathrm{T}}-\hat{Q}_b({\bf w} {\bf n}_1^{\mathrm{T}}),A_1 \nabla {\bf u}_1 -A_1 \mathbb{Q}_h(\nabla {\bf u}_1) \rangle_{\partial T \cap \Gamma}
	\Big)\\
	+&\sum_{T \in \mathcal{T}_{2h}^I}\Big(
	\langle Q_0 {\bf w} -Q_b {\bf w}, A_2 (\nabla {\bf u}_2){\bf n}-A_2 \mathbb{Q}_h(\nabla {\bf u}_2){\bf n}_2 \rangle_{\partial T \backslash \Gamma} \\
	+ &\langle (Q_0{\bf w}) {\bf n}_2^{\mathrm{T}}-\hat{Q}_b({\bf w} {\bf n}_2^{\mathrm{T}}),A_1 \nabla {\bf u}_2 -A_2 \mathbb{Q}_h(\nabla {\bf u}_2) \rangle_{\partial T \cap \Gamma}
	\Big)\\
\end{align*}
Since the terms on the straight edges are similar to $\ell_1({\bf u},Q_h {\bf w})$, we only need to verify the terms on the interface. According to the trace inequality, the Cauchy-Schwarz inequality and the Lemma \ref{projector inequality 1}, we have
\begin{align*}
	&\sum_{T \in \mathcal{T}_{1h}^I} \langle (Q_0{\bf w}) {\bf n}_1^{\mathrm{T}}-{\bf w} {\bf n}_1^{\mathrm{T}}, A_1 (\nabla {\bf u}_1){\bf n}_1 - A_1 \mathbb{Q}_h(\nabla {\bf u}_1){\bf n}_1   \rangle_{\partial T \cap \Gamma}\\
	\leqslant & \| Q_0 {\bf w} -{\bf w} \|_{\partial T \cap \Gamma}
	\|  A_1 (\nabla {\bf u}_1)- A_1 \mathbb{Q}_h(\nabla {\bf u}_1)\|_{\partial T \cap \Gamma}\\
	\leqslant &\Big(\sum_{T \in \mathcal{T}_{1h}^I} \|Q_0 {\bf w} -{\bf w} \|^2_{\partial T \cap \Gamma} \Big)^{\frac{1}{2}} \Big(\sum_{T \in \mathcal{T}_{1h}^I}
	\|A_1 (\nabla {\bf u}_1)- A_1 \mathbb{Q}_h(\nabla {\bf u}_1) \|^2_{\partial T \cap \Gamma} \Big)^{\frac{1}{2}}\\
	\leqslant & C h^{k+1} \| {\bf u}_1\|_{k+1} \| {\bf w} \|_2.
\end{align*}
Similarly, we have
\begin{align*}
	&\sum_{T \in \mathcal{T}_{1h}^I} \langle {\bf w} {\bf n}_1^{\mathrm{T}}-\hat{Q}_b({\bf w} {\bf n}_1^{\mathrm{T}}), A_1 \nabla {\bf u}_1 -A_1 \mathbb{Q}_h(\nabla {\bf u}_1)  \rangle_{\partial T \cap \Gamma}\\
	\leqslant &\Big(\sum_{T \in \mathcal{T}_{1h}^I} \|{\bf w} {\bf n}_1^{\mathrm{T}}-\hat{Q}_b({\bf w} {\bf n}_1^{\mathrm{T}}) \|^2_{\partial T \cap \Gamma} \Big)^{\frac{1}{2}} \Big(\sum_{T \in \mathcal{T}_{1h}^I}
	\|A_1 (\nabla {\bf u}_1)- A_1 \mathbb{Q}_h(\nabla {\bf u}_1) \|^2_{\partial T \cap \Gamma} \Big)^{\frac{1}{2}}\\
	\leqslant & C h^{k+1} \| {\bf u}_1\|_{k+1} \| {\bf w} \|_2.
\end{align*}
Therefore, we have
\begin{align*}
	&\sum_{T \in \mathcal{T}_{1h}^I} \langle (Q_0{\bf w}) {\bf n}_1^{\mathrm{T}}-\hat{Q}_b({\bf w} {\bf n}_1^{\mathrm{T}}),A_1 \nabla {\bf u}_1 -A_1 \mathbb{Q}_h(\nabla {\bf u}_1) \rangle_{\partial T \cap \Gamma}\\
	=&\sum_{T \in \mathcal{T}_{1h}^I} \langle (Q_0{\bf w}) {\bf n}_1^{\mathrm{T}}-{\bf w} {\bf n}_1^{\mathrm{T}}, A_1 (\nabla {\bf u}_1){\bf n}_1 - A_1 \mathbb{Q}_h(\nabla {\bf u}_1){\bf n}_1   \rangle_{\partial T \cap \Gamma}\\
	+&\sum_{T \in \mathcal{T}_{1h}^I} \langle {\bf w} {\bf n}_1^{\mathrm{T}}-\hat{Q}_b({\bf w} {\bf n}_1^{\mathrm{T}}), A_1 \nabla {\bf u}_1 -A_1 \mathbb{Q}_h(\nabla {\bf u}_1)  \rangle_{\partial T \cap \Gamma}\\
	\leqslant & C h^{k+1} \| {\bf u}_1\|_{k+1} \| {\bf w} \|_2.
\end{align*}
Using the same method leads to
\begin{align*}
	&\sum_{T \in \mathcal{T}_{2h}^I} \langle (Q_0{\bf w}) {\bf n}_2^{\mathrm{T}}-\hat{Q}_b({\bf w} {\bf n}_2^{\mathrm{T}}),A_2 \nabla {\bf u}_2 -A_2 \mathbb{Q}_h(\nabla {\bf u}_2) \rangle_{\partial T \cap \Gamma}\\
	=&\sum_{T \in \mathcal{T}_{2h}^I} \langle (Q_0{\bf w}) {\bf n}_2^{\mathrm{T}}-{\bf w} {\bf n}_2^{\mathrm{T}}, A_2 (\nabla {\bf u}_2){\bf n}_2 - A_2 \mathbb{Q}_h(\nabla {\bf u}_2){\bf n}_2   \rangle_{\partial T \cap \Gamma}\\
	+&\sum_{T \in \mathcal{T}_{2h}^I} \langle {\bf w} {\bf n}_2^{\mathrm{T}}-\hat{Q}_b({\bf w} {\bf n}_2^{\mathrm{T}}), A_2 \nabla {\bf u}_2 -A_2 \mathbb{Q}_h(\nabla {\bf u}_2)  \rangle_{\partial T \cap \Gamma}\\
	\leqslant & C h^{k+1} \| {\bf u}_2\|_{k+1} \| {\bf w} \|_2.
\end{align*}
In conclusion, we obtain
\begin{align*}
	\tilde{\ell}_1({\bf u},Q_h {\bf w}) \leqslant Ch^{k+1}(\| {\bf u}_1 \|_{k+1}+\| {\bf u}_2 \|_{k+1}) \| {\bf w} \|_2.
\end{align*}
(3)for $\ell_2(p,Q_h {\bf w})$, we use the same method as the proof of $\ell_1({\bf u},Q_h {\bf w})$ to obtain
\begin{align*}
&\sum_{T \in \mathcal{T}_h} \langle Q_0 {\bf w} -{\bf w}, (p-Q_h^p p){\bf n} \rangle_{\partial T}\\
\leqslant & \sum_{T \in \mathcal{T}_h}  \| Q_0 {\bf w} -{\bf w} \|_{\partial T} \| p-Q_h^p p \|_{\partial T}\\
\leqslant & \Big( \sum_{T \in \mathcal{T}_h} \| Q_0 {\bf w} -{\bf w} \|_{\partial T}^2 \Big)^{\frac{1}{2}}
\Big( \sum_{T \in \mathcal{T}_h} \| p-Q_h^p p \|_{\partial T}^2\Big)^{\frac{1}{2}}\\
\leqslant & \Big( C \sum_{T \in \mathcal{T}_h} (h_T^{-1} \| Q_0 {\bf w} -{\bf w} \|_T^2
+ h_T \| \nabla(Q_0 {\bf w} -{\bf w})\|_T^2 )\Big)^{\frac{1}{2}}\\
&\Big( C \sum_{T \in \mathcal{T}_h} (h_T^{-1} \| p-Q_h^p p \|_T^2
+ h_T \| \nabla(p-Q_h^p p)\|_T^2 ) \Big)^{\frac{1}{2}}\\
\leqslant &\Big( \sum_{T \in \mathcal{T}_h} C h_T^3 \| {\bf w} \|_2^2 \Big)^{\frac{1}{2}} \Big( \sum_{T \in \mathcal{T}_h} Ch^{2k-1}
\| p\|_k^2  \Big)^{\frac{1}{2}}\\
\leqslant & C h^{k+1} \| {\bf w} \|_2 \| p \|_k.
\end{align*}
Similarly, we can arrive at
\begin{align*}
&\sum_{T \in \mathcal{T}_h} \langle {\bf w} - Q_b {\bf w}, (p-Q_h^p p){\bf n} \rangle_{\partial T}
\leqslant C h^{k+1} \| {\bf w} \|_2  \|p\|_k.
\end{align*}
Therefore, for $\ell_2(p,Q_h {\bf w})$ we have the following estimate
\begin{align*}
|\ell_2(p,Q_h {\bf w})| \leqslant C h^{k+1} \| {\bf w} \|_2 \| p \|_k.
\end{align*}
(4)For $\tilde{\ell}_2(p,Q_h {\bf w})$, we have
\begin{align*}
	\tilde{\ell}_2(p,Q_h {\bf w}) \leqslant Ch^{k+1}(\| p_1 \|_k+\| p_2 \|_k) \| {\bf w} \|_2.
\end{align*}
(5)for $\tilde{S}(Q_h {\bf u}, Q_h {\bf w})$,we first consider the element $T \in \mathcal{T}_h^S$ ,
\begin{align*}
&\sum_{T \in \mathcal{T}_h^S} h_T^{-1} \langle Q_b(Q_0 {\bf u})- Q_b {\bf u} ,Q_b(Q_0 {\bf w})-Q_b {\bf w} \rangle_{\partial T}\\
=&\sum_{T \in \mathcal{T}_h^S} h_T^{-1} \langle Q_0 {\bf u} -{\bf u} ,Q_b(Q_0 {\bf w})-Q_b {\bf w} \rangle_{\partial T}\\
\leqslant &  \sum_{T \in \mathcal{T}_h^S} h_T^{-1} \| Q_0 {\bf u} -{\bf u} \|_{\partial T}
\| Q_0 {\bf w} -{\bf w}\|_{\partial T}\\
\leqslant & \Big( \sum_{T \in \mathcal{T}_h^S} h_T^{-1} \| Q_0 {\bf u} -{\bf u} \|_{\partial T}^2 \Big)^{\frac{1}{2}}
\Big(\sum_{T \in \mathcal{T}_h^S}  h_T^{-1} \| Q_0 {\bf w} -{\bf w} \|_{\partial T}^2\Big)^{\frac{1}{2}}\\
\leqslant & \Big( C \sum_{T \in \mathcal{T}_h^S} (h_T^{-1} \| Q_0 {\bf u} -{\bf u} \|_T^2
+ h_T\| \nabla(Q_0 {\bf u} -{\bf u})\|_T^2 ) \Big)^{\frac{1}{2}} \\
&\Big( C \sum_{T \in \mathcal{T}_h^S} (h_T^{-1} \| Q_0 {\bf w} -{\bf w} \|_T^2
+ h_T \| \nabla(Q_0 {\bf w} -{\bf w})\|_T^2 ) \Big)^{\frac{1}{2}}\\
\leqslant & C h^{k+1} \| {\bf u} \|_{k+1} \| {\bf w} \|_2.
\end{align*}
Next, for the element $T \in \mathcal{T}_{1h}^I$, we have
\begin{align*}
&\tilde{s}(Q_h {\bf u},Q_h {\bf w})|_T \\
= &h_T^{-1} \langle Q_b(Q_0 {\bf u})- Q_b {\bf u} ,Q_b(Q_0 {\bf w})-Q_b {\bf w} \rangle_{\partial T \backslash \Gamma}\\
+&h_T^{-1} \langle \hat{Q}_b ((Q_0 {\bf u}_1){\bf n}_1^{\mathrm{T}}) - \hat{Q}_b({\bf u}_1 {\bf n}_1^{\mathrm{T}}),
 \hat{Q}_b ((Q_0 {\bf w}_1){\bf n}_1^{\mathrm{T}}) - \hat{Q}_b({\bf w}_1 {\bf n}_1^{\mathrm{T}}) \rangle_{\partial T \cap \Gamma}\\
\end{align*}
The first term is similar to stability term on triangular element, so we deal with the second term
\begin{align*}
&h_T^{-1} \langle \hat{Q}_b((Q_0 {\bf u}_1){\bf n}_1^{\mathrm{T}})-\hat{Q}_b({\bf u}_1 {\bf n}_1^{\mathrm{T}}), \hat{Q}_b ((Q_0 {\bf w}_1){\bf n}_1^{\mathrm{T}}) - \hat{Q}_b({\bf w}_1 {\bf n}_1^{\mathrm{T}}) \rangle_{\partial T \cap \Gamma}\\
=&h_T^{-1} \langle (Q_0 {\bf u}_1){\bf n}_1^{\mathrm{T}}- {\bf u}_1 {\bf n}_1^{\mathrm{T}},
\hat{Q}_b ((Q_0 {\bf w}_1){\bf n}_1^{\mathrm{T}}) - \hat{Q}_b({\bf w}_1 {\bf n}_1^{\mathrm{T}}) \rangle_{\partial T \cap \Gamma}\\
\leqslant & h_T^{-1} \| Q_0 {\bf u}_1 - {\bf u}_1\|_{\partial T \cap \Gamma}
\| \hat{Q}_b ((Q_0 {\bf w}_1){\bf n}_1^{\mathrm{T}}-{\bf w}_1 {\bf n}_1^{\mathrm{T}}) \|_{\partial T \cap \Gamma}\\
\leqslant & h_T^{-1} \| Q_0 {\bf u}_1 - {\bf u}_1\|_{\partial T \cap \Gamma}
\| (Q_0 {\bf w}_1){\bf n}_1^{\mathrm{T}}-{\bf w}_1 {\bf n}_1^{\mathrm{T}} \|_{\partial T \cap \Gamma}\\
\leqslant & Ch^{k+1}\| {\bf u}_1\|_{k+1} \| {\bf w} \|_2
\end{align*}
In addition, for the element $T \in \mathcal{T}_{2h}^I$, we can obtain
\begin{align*}
&\tilde{s}(Q_h {\bf u},Q_h {\bf w})|_T \\
= &h_T^{-1} \langle Q_b(Q_0 {\bf u})- Q_b {\bf u} ,Q_b(Q_0 {\bf w})-Q_b {\bf w} \rangle_{\partial T \backslash \Gamma}\\
+&h_T^{-1} \langle \hat{Q}_b ((Q_0 {\bf u}_2){\bf n}_2^{\mathrm{T}}) -\hat{Q}_b({\bf u}_2 {\bf n}_2^{\mathrm{T}}),
 \hat{Q}_b ((Q_0 {\bf w}_1){\bf n}_1^{\mathrm{T}}) + \hat{Q}_b({\bf w}_1 {\bf n}_1^{\mathrm{T}}) \rangle_{\partial T \cap \Gamma}\\
\leqslant & Ch^{k+1} \|u_2\|_{k+1,\Omega_2}\| {\bf w}\|_2
\end{align*}
Therefore,
$$
\tilde{s}(Q_h {\bf u}, Q_h {\bf w}) \leqslant Ch^{k+1} (\|{\bf u}_1\|_{k+1,\Omega_1}+\|{\bf u}_2\|_{k+1,\Omega_2})\|{\bf w}\|_2
$$
Combining the five estimates and the regularity condition, we can obtain
\begin{align*}
&\|{\bf e}_0\|^2\\
\leqslant &C h^{k+1}(\| {\bf u}_1 \|_{k+1,\Omega_1}+\|{\bf u}_2\|_{k+1,\Omega_2}+\|p_1\|_{k,\Omega_1}+\|p_2\|_{k,\Omega_2}) \| {\bf w} \|_2
+Ch{|\hspace{-.02in}|\hspace{-.02in}|} {\bf e}_h {|\hspace{-.02in}|\hspace{-.02in}|} \|{\bf e}_0\|\\
\leqslant &C h^{k+1}(\|{\bf u}_1\|_{k+1,\Omega_1}+\|{\bf u}_2\|_{k+1,\Omega_2}+\|p_1\|_{k,\Omega_1}+\|p_2\|_{k,\Omega_2})\| {\bf e}_0\|
+Ch{|\hspace{-.02in}|\hspace{-.02in}|} {\bf e}_h {|\hspace{-.02in}|\hspace{-.02in}|} \|{\bf e}_0\|
\end{align*}
Finally, according to Theorem \ref{H1ERROR}, it follows that
\begin{align*}
\| {\bf e}_0 \|
\leqslant & C h^{k+1}(\|{\bf u}_1\|_{k+1,\Omega_1}+\|{\bf u}_2\|_{k+1,\Omega_2}+\|p_1\|_{k,\Omega_1}+\|p_2\|_{k,\Omega_2})
+Ch{|\hspace{-.02in}|\hspace{-.02in}|} {\bf e}_h {|\hspace{-.02in}|\hspace{-.02in}|} \\
\leqslant & C h^{k+1}(\|{\bf u}_1\|_{k+1,\Omega_1}+\|{\bf u}_2\|_{k+1,\Omega_2}+\|p_1\|_{k,\Omega_1}+\|p_2\|_{k,\Omega_2})
\end{align*}
The proof of theorem is completed.
\end{proof}

\section{Numerical Results} In this section, we will give some numerical examples to validate the previous theory. In the numerical tests, we give the results of examples with discontinuous coefficients A, pressure p and velocity ${\bf u}$ at the interface for different interface shapes.

\subsection{Test Problem 1}In the test problem 1, we consider the following problems with discontinuous velocity and discontinuous pressure on the domain $\Omega=[-1,1] \times [-1,1]$ and the interface is described as
\begin{align*}
x^2+y^2=\frac{1}{4}.
\end{align*}

The exact solutions are
\begin{eqnarray*}
{\bf u}_1=\left(\begin{array}{ccc}
2 \sin y \cos y \cos x \\
(\sin^2 y-2)\sin x
\end{array}\right),
\end{eqnarray*}
\begin{eqnarray*}
	{\bf u}_2=\left(\begin{array}{ccc}
	- \cos(\pi x) \sin(\pi y)\\
      \sin(\pi x) \cos(\pi y)
	\end{array}\right),
\end{eqnarray*}

\begin{eqnarray}
	p=\left \{\begin{array}{rcl}
		1 & \mbox{in} & \Omega_1 \\
		\frac{\pi}{16-\pi} & \mbox{in} & \Omega_2
	\end{array}\right.,
   A=\left \{\begin{array}{rcl}
        1 & \mbox{in} & \Omega_2 \\
	    1 & \mbox{in} & \Omega_2
   \end{array}\right.
\end{eqnarray}

\begin{table}[H]
	\caption{Test Problem 1: curved triangular mesh}
	\centering
	\begin{tabular}{ccccccccc}
		\hline
		n&${|\hspace{-.02in}|\hspace{-.02in}|} Q_h {\bf u} -{\bf u}_h {|\hspace{-.02in}|\hspace{-.02in}|}$&order&$\|Q_0 {\bf u} -{\bf u}_0\|$&order&$\|Q_h^p p- p_h\|$&order&$\|{\bf u}-{\bf u}_h \|_{L^\infty(\Omega)}$&order\\
		\hline
		&&&&$k=1$&&&&\\
		\hline
		 1 & 3.9611e+00 &-- & 6.5572e-01 & --& 8.6957e-01&--&7.4688e-01&-- \\
		 2 & 1.9447e+00 &  1.026 & 1.5798e-01&  2.053  & 4.1120e-01&1.080&1.3056e-01 & 2.516 \\
		 3 & 9.6473e-01& 1.041 &3.9173e-02&2.072 & 1.9511e-01& 1.108&3.5175e-02& 1.948  \\
		 4 & 4.8011e-01 &1.055 &9.7663e-03 &2.101 & 9.5708e-02&1.077&8.8538e-03 &2.086\\
		\hline
		&&&$k=2$&&&\\
		\hline
        1 & 9.9607e-01 &--& 1.2002e-01 &--& 1.3309e-01 &-- & 5.7368e-02 & --\\
        2 & 2.5843e-01 & 1.946 & 1.6144e-02 & 2.894 & 2.9759e-02 &2.161  & 9.1152e-03& 2.654\\
        3 & 6.4942e-02 & 2.052  & 2.0504e-03 & 3.066 & 6.8283e-03 &2.187& 1.2027e-03 & 3.009 \\
        4 & 1.6237e-02 & 2.097 & 2.5745e-04 &3.138& 1.6139e-03 &2.182 & 1.5031e-04 & 3.145\\
		\hline
		&&&$k=3$&&&\\
		\hline
		1 & 1.9557e-01 &--& 2.2154e-02 &--& 1.8260e-02 &-- & 5.7239e-03&--\\
		2 & 2.4560e-02 & 2.993 & 1.3796e-03 &4.005& 2.1947e-03 &3.057 &5.5005e-04 & 3.379\\
		3 & 3.0687e-03 & 3.089 & 8.6147e-05 &4.120& 2.6353e-04 &3.149 & 3.8595e-05& 3.947\\
		4 & 3.8338e-04 & 3.146 & 5.3835e-06 &4.197& 3.2198e-05 &3.179 & 2.4991e-06&4.139\\
		\hline
	\end{tabular}
\end{table}

\begin{figure*}[t!]
	\centering
	\begin{subfigure}[t]{0.5\linewidth}
		\centering
	    \includegraphics[width=0.8\linewidth]{pics/C432-1.jpg}
		\caption{}
	\end{subfigure}%
	\begin{subfigure}[t]{0.5\linewidth}
		\centering
	    \includegraphics[width=0.8\linewidth]{pics/C432-2.jpg}
		\caption{}
	\end{subfigure}
	\begin{subfigure}[t]{0.5\linewidth}
			\centering \includegraphics[width=0.8\linewidth]{pics/C432-3.jpg}
			\caption{}
	\end{subfigure}
\caption{Test problem 1: (a) x-component of ${\bf u}_h$, (b) y-component of ${\bf u}_h$, (c) pressure $p_h$}
\end{figure*}
\begin{table}[H]
	\caption{Test Problem 1: straight triangular mesh}
	\centering
	\begin{tabular}{ccccccccc}
		\hline
		n&${|\hspace{-.02in}|\hspace{-.02in}|} Q_h {\bf u} -{\bf u}_h {|\hspace{-.02in}|\hspace{-.02in}|}$&order&$\|Q_0 {\bf u} -{\bf u}_0\|$&order&$\|Q_h^p p- p_h\|$&order&$\|{\bf u}-{\bf u}_h \|_{L^\infty(\Omega)}$&order\\
		\hline
		&&&&$k=1$&&&&\\
		\hline
        1 &  3.8965e+00 &--& 6.9519e-01 &--&  8.9801e-01 &-- &8.7996e-01&--\\
        2 &  1.9151e+00 & 1.025 & 1.6636e-01 & 2.063& 4.1484e-01&1.114&2.3388e-01 &1.912\\
        3 &9.5547e-01 &1.033& 4.1203e-02 &2.073 &  &1.118&5.9402e-02 &2.036\\
        4 & 4.7758e-01  & 1.049 & 1.0277e-02&2.100& 9.5744e-02 &1.079&1.4933e-02&2.088\\
		\hline
		&&&&$k=2$&&&&\\
		\hline
        1 &1.0515 &--& 1.4146e-01 &--& 2.1288e-01 &--&1.7603e-01 &--\\
        2 &2.7967e-01 &1.911& 2.2582e-02  &2.647&6.4754e-02 & 1.7170&3.1879e-02 &2.465\\
        3 & 7.4373e-02 & 1.968  & 4.2713e-03 &2.474 & 2.0841e-02 &1.684&6.1469e-03&2.445\\
        4 & 2.0444e-02  &1.953 & 9.5717e-04 &2.262 & 7.0023e-03 &1.649&1.7066e-03&1.938\\
		\hline
		&&&&$k=3$&&&&\\
		\hline
        1 & 3.7183e-01 &--& 6.1249e-02 &--& 1.6822e-01 &--&6.7480e-02 &--\\
        2 & 9.7064e-02 & 1.938 &1.4642e-02   &2.065& 5.7458e-02 &1.549&1.7359e-02& 1.959\\
        3 &  2.9548e-02  &1.767 & 3.6609e-03  &2.059&  1.9566e-02 &1.600&4.5858e-03& 1.978
        \\
        4 & 9.6547e-03 &1.692 & 9.1541e-04  &2.096&  6.7761e-03 &1.604  &1.16013e-03& 2.079
        \\
		\hline
	\end{tabular}
\end{table}
\qquad We perform the weak Galerkin approximation on the curved triangular mesh in table 7.1 and the approximation performance on the straight triangular mesh in table 7.2. The convergency rates on the curved triangular mesh for linear and high order weak Galerkin finite elements arrive at the optimal convergency order, but the convergency rates on the straight triangular mesh are less than or equal to second order. This results show the advantages for using curved element.

\subsection{Test Problem 2}In the test problem 2, we consider the following problems with discontinuous velocity and the following interface:
\begin{align*}
r=\frac{1}{7}+\frac{sin(5 \theta)}{7}
\end{align*}
The exact solutions are
\begin{eqnarray*}
	{\bf u}_1=\left(\begin{array}{ccc}
		2 \pi \sin(\pi x) \sin(\pi x) \cos(\pi y) \sin(\pi y) \\
		-2 \pi \sin(\pi x) \sin(\pi y) \cos(\pi x) \sin(\pi y)
	\end{array}\right)
\end{eqnarray*}
\begin{eqnarray*}
	{\bf u}_2=\left(\begin{array}{ccc}
		x^2 y^2+e^{-y}\\
		-\frac{2}{3} xy^3+2-\pi \sin(\pi x)
	\end{array}\right)
\end{eqnarray*}
\begin{eqnarray}
	p=\left \{\begin{array}{rcl}
		0 & \mbox{in} & \Omega_1 \\
		0 & \mbox{in} & \Omega_2
	\end{array}\right.
	A=\left \{\begin{array}{rcl}
		1 & \mbox{in} & \Omega_2 \\
		1 & \mbox{in} & \Omega_2
	\end{array}\right.
\end{eqnarray}
\begin{table}[H]
	\caption{Test Problem 2}
	\centering
	\begin{tabular}{ccccccccc}
		\hline
		n&${|\hspace{-.02in}|\hspace{-.02in}|} Q_h {\bf u} -{\bf u}_h {|\hspace{-.02in}|\hspace{-.02in}|}$&order&$\|Q_0 {\bf u} -{\bf u}_0\|$&order&$\|Q_h^p p- p_h\|$&order&$\|{\bf u}-{\bf u}_h \|_{L^\infty(\Omega)}$&order\\
		\hline
		&&&$k=1$&&&\\
		\hline
		1 & 8.3523e+00 &--& 6.8303e-01 &--& 1.3395e+00 &--&1.6144e+00&-- \\
		2 & 4.2334e+00 & 0.973 & 1.8409e-01 &1.876& 7.0321e-01 & 0.922&5.7339e-01&-- \\
		3 & 2.0738e+00 & 1.141 & 4.6527e-02 &2.198& 3.5766e-01 &1.081&1.5624e-01 &-- \\
		4 & 1.0517e+00 & 0.988 & 1.2162e-02 & 1.952 &1.4259e-01 & 1.338&4.5628e-02 &-- \\

	    \hline
	    &&&$k=2$&&&\\
	    \hline
        1 & 1.8413e+00 &--& 1.0910e-01 &--& 2.0817e-01 &-- &3.9152e-01&\\
        2 & 4.9217e-01 &1.888& 1.5530e-02 &2.789& 4.5704e-02 &2.169&9.2244e-02&2.069\\
        3 & 1.2483e-01 & 2.193 & 2.0196e-03 &3.260& 1.1769e-02 &2.168&1.0086e-02 &3.538 \\
        4 & 3.1483e-02 & 2.004 & 2.5750e-04 &2.997& 2.6264e-03 &2.182 &1.3225e-03&2.956\\

    	\hline
	    &&&$k=3$&&&\\
    	\hline
        1 & 2.8573e-01 &--& 1.4530e-02 &--& 2.7239e-02 &-- &9.0719e-02&\\
        2 & 3.8777e-02 & 2.858 & 1.0229e-03 &3.798& 3.2468e-03 &3.044 &7.0040e-03&3.665\\
        3 & 5.1489e-03 & 3.227 & 6.9235e-05 &4.304& 4.0691e-04 &3.319 &6.6249e-04&3.769\\
        4 & 6.6696e-04 & 2.974  & 4.5695e-06 &3.955& 4.6397e-05 &3.159 &4.6409e-05& 3.868\\

	    \hline
\end{tabular}
\end{table}

\begin{figure*}[t!]
	\centering
	\begin{subfigure}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=0.8\linewidth]{pics/F432-1.jpg}
			\caption{}
		\end{subfigure}%
	\begin{subfigure}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=0.8\linewidth]{pics/F432-2.jpg}
			\caption{}
		\end{subfigure}
	\begin{subfigure}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=0.8\linewidth]{pics/F432-3.jpg}
			\caption{}
		\end{subfigure}
	\caption{Test problem 2: (a) x-component of ${\bf u}_h$, (b) y-component of ${\bf u}_h$, (c) pressure $p_h$}
\end{figure*}
\subsection{Test Problem 3}In the test problem 3, we consider the following problems with discontinuous viscosity coefficient and  the following interface:
\begin{align*}
	r=\frac{1}{2}+\frac{sin(2 \theta)}{4}
\end{align*}
The exact solutions are
\begin{eqnarray*}
	{\bf u}_1=\left(\begin{array}{ccc}
		2 \pi \sin(\pi x) \sin(\pi x) \cos(\pi y) \sin(\pi y) \\
		-2 \pi \sin(\pi x) \sin(\pi y) \cos(\pi x) \sin(\pi y)
	\end{array}\right)
\end{eqnarray*}
\begin{eqnarray*}
	{\bf u}_2=\left(\begin{array}{ccc}
		x^2 y^2+e^{-y}\\
		-\frac{2}{3} xy^3+2-\pi \sin(\pi x)
	\end{array}\right)
\end{eqnarray*}
\begin{eqnarray}
	p=\left \{\begin{array}{rcl}
		0 & \mbox{in} & \Omega_1 \\
		0 & \mbox{in} & \Omega_2
	\end{array}\right.
	A=\left \{\begin{array}{rcl}
		1 & \mbox{in} & \Omega_2 \\
		10 & \mbox{in} & \Omega_2
	\end{array}\right.
\end{eqnarray}
\begin{table}[H]
	\caption{Test Problem 3}
	\centering
	\begin{tabular}{ccccccccc}
		\hline
		n&${|\hspace{-.02in}|\hspace{-.02in}|} Q_h {\bf u} -{\bf u}_h {|\hspace{-.02in}|\hspace{-.02in}|}$&order&$\|Q_0 {\bf u} -{\bf u}_0\|$&order&$\|Q_h^p p- p_h\|$&order&$\|{\bf u}-{\bf u}_h \|_{L^\infty(\Omega)}$&order\\
		\hline
		&&&&$k=1$&&&&\\
		\hline
        1 & 3.7339e+01 &-- & 3.4134e+00 &--& 7.7369e+00 &--&3.0523e+01& \\
        2 & 1.9557e+01 & 0.8879 & 9.3621e-01& 1.7762  & 3.3129e+00 &1.1646 & 1.2689e+00& 4.367\\
        3 & 9.9253e+00 &1.0701& 2.4258e-01 &2.1307& 1.3139e+00 &1.4591 & 2.3145e-01 &2.685  \\
        4 & 4.9796e+00 & 1.0953 & 6.1255e-02 &2.1855 & 5.3379e-01 &1.4303 & 6.0317e-02& 2.135 \\

		\hline
		&&&&$k=2$&&&&\\
		\hline
        1 & 4.7365e+00 &--& 3.3909e-01 &--& 6.5770e-01 &--&   7.4812e-01&\\
        2 & 1.2581e+00 & 1.8202 & 4.5822e-02 & 2.7482& 1.2309e-01 &2.3010 &1.3168e-01& 2.3852 \\
        3 & 3.1669e-01 & 2.1764 & 5.8042e-03 &3.2598& 2.6955e-02 &2.3961 &1.5557e-02& 3.3698 \\
        4 & 7.9984e-02 & 2.1852 & 7.3506e-04 &3.2813& 5.7394e-03 &2.4563 &2.5042e-03&2.901 \\

		\hline
	    &&&&$k=3$&&&&\\
		\hline
        1 & 4.7047e-01 &--& 2.7703e-02 &--& 7.6521e-02 &--&1.2870e-01&-- \\
        2 & 6.6885e-02 & 2.6785 & 1.9608e-03 &3.6361& 8.8654e-03 &2.9595& 9.8594e-03&3.527 \\
        3 & 8.6199e-03 & 3.2325 & 1.2957e-04 &4.2864& 7.0084e-04 &4.0036&6.8631e-04&4.204 \\
        4 & 1.0955e-03 & 3.2757 & 8.2702e-06 &4.3694& 7.6611e-05 &3.5150 &4.8863e-05&4.196\\

		\hline
	\end{tabular}
\end{table}
\begin{figure*}[t!]
	\centering
	\begin{subfigure}[t]{0.5\linewidth}
		\centering
		\includegraphics[width=0.8\linewidth]{pics/P432-1.jpg}
		\caption{}
	\end{subfigure}%
	\begin{subfigure}[t]{0.5\linewidth}
		\centering
		\includegraphics[width=0.8\linewidth]{pics/P432-2.jpg}
		\caption{}
	\end{subfigure}
	\begin{subfigure}[t]{0.5\linewidth}
		\centering
		\includegraphics[width=0.8\linewidth]{pics/P432-3.jpg}
		\caption{}
	\end{subfigure}
    \caption{Test problem 3: (a) x-component of ${\bf u}_h$, (b) y-component of ${\bf u}_h$, (c) pressure $p_h$}
\end{figure*}

\section{Conclusion}In this paper, we mainly use the weak Galerkin finite element method to deal with Stokes interface problems with complex interface. We present a weak finite element numerical scheme with only one value at the interface. It is proved that the optimal convergence order can be reached under the energy norm and the standard $L^2$ norm. At the same time, the results of numerical examples also show that the error convergence order under the energy norm and L2 norm is optimal, which is consistent with the theory.
\bibliographystyle{siam}